text:
module twisting_circle_examples
# Twisting of a circle into shape with two turns

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexBeams.FEMMCorotBeamModule
using FinEtoolsFlexBeams.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
distribloads_global = FEMMCorotBeamModule.distribloads_global
restoringforce = FEMMCorotBeamModule.restoringforce
gyroscopic = FEMMCorotBeamModule.gyroscopic
using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield, linear_update_rotation_field!, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexBeams.VisUtilModule: plot_space_box, plot_midline, plot_solid, render, react!, default_layout_3d, save_to_json
using PlotlyJS
using JSON

function twisting_circle()
    # Parameters:
    # The  parameters given in the paper  are in terms of stiffness and mass coefficients.
    E = 200000.0;
    nu = 0.3;# Poisson ratio
    rho = 7.8e-9;
    b = 0.6; h = 6.0;  # cross-sectional dimensions
    radius = 120 # radius of the circle
    utol = 1e-3;
    maxit = 30
    Tmag = 700;
    dt = 0.0002;
    tend = 30.0;
    ng = 1/2; nb = 1/4*(1/2+ng)^2;
    # Select the number of elements around the circumference.
    n=32;
    # Choose the mass formulation:
    mass_type=1;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> h, s -> b, s -> [0.0, 0.0, 1.0])

    @show tolerance = radius/n/1000;
    fens, fes = frame_member([0 0 0; 2*pi 0 0], n, cs)
    for i in 1:count(fens)
        a = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
    end
    fens, fes = mergenodes(fens, fes, tolerance, [1, n+1])

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Locate some nodes
    box = fill(0.0, 6)
    # The pinned end is
    initbox!(box, [0.0,0.0,0.0])
    clampedn = selectnode(fens; box = box, inflate = tolerance)
    # The force is applied at the elbow
    initbox!(box, [2*radius,0,0])
    torquen = selectnode(fens; box = box, inflate = tolerance)

    # # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    for i in [2, 3, 5, 6]
        setebc!(dchi, torquen, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    stepdchi = deepcopy(dchi)
    u1 = deepcopy(u0)
    v0 = deepcopy(dchi) # need the numbering of the degrees of freedom
    v1 = deepcopy(dchi)
    a1 = deepcopy(dchi) # zero out the Acceleration
    a0 = deepcopy(dchi) # zero out the Acceleration
    Rfield1 = deepcopy(Rfield0)
    v0v = gathersysvec(v0)
    a0v = gathersysvec(a0)
    vpv = gathersysvec(v0);
    dchipv = gathersysvec(dchi);
    stepdchiv = gathersysvec(dchi);
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    tbox = plot_space_box([[0 -radius -radius]; [2*radius radius radius]])
    tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    plots = cat(tbox,  tshape0; dims = 1)
    pl = render(plots)
    sleep(0.5)

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    loadbdry = FESetP1(reshape(torquen, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))

    t = 0.0; #
    step = 0;
    while (t <= tend)
        t = t + dt;
        # println("Time $(t)"); # pause
        # Initialization
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        stepdchi.values[:] .= 0.0;# Total increment in current step
        a1.values[:] = -(1/nb/dt)*v0.values[:] -(1/2-nb)/nb*a0.values[:];
        v1.values[:] = v0.values[:] + dt*((1-ng)*a0.values[:] + ng*a1.values[:]);
        gathersysvec!(v0, v0v)
        gathersysvec!(a0, a0v)
        dchipv = dt*v0v + (dt^2/2*(1-2*nb))*a0v
        vpv = v0v +(dt*(1-ng))*a0v;
        fi = ForceIntensity(FFlt[0, 0, 0, (t < 0.08)*Tmag, 0, 0]);

        iter = 1;
        while true
            F = distribloads(lfemm, geom0, dchi, fi, 3);
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            @. rhs = F + Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi);
            M = mass(femm, geom0, u1, Rfield1, dchi);
            G = gyroscopic(femm, geom0, u1, Rfield1, v1, dchi);
            gathersysvec!(stepdchi, stepdchiv)
            @. TMPv = ((-1/(nb*dt^2))*stepdchiv+(1/(nb*dt^2))*dchipv)
            rhs .+= M*TMPv
            @. TMPv = ((-ng/nb/dt)*stepdchiv+(ng/nb/dt)*dchipv - vpv)
            rhs .+= G*TMPv;
            dchi = scattersysvec!(dchi, (K+(ng/nb/dt)*G+(1/(nb*dt^2))*M)\rhs); # Disp. incr
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            stepdchi.values[:] += dchi.values[:]
            v1.values[:] += (ng/nb/dt)*dchi.values[:];
            a1.values[:] += (1/nb/dt^2)*dchi.values[:];
            update_rotation_field!(Rfield1, dchi)
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations
        v0.values[:] = v1.values[:];       # update the velocities
        a0.values[:] = a1.values[:];       # update the accelerations

        if mod(step, 4) == 0
            @show step
            tshape1 = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
            plots = cat(tbox,  tshape0,  tshape1; dims = 1)
            pl.plot.layout[:title] = "t=$t"
            react!(pl, plots, pl.plot.layout)
            sleep(0.5)
        end

        step=step+1;
    end

    return true
end # twisting_circle

function allrun()
    println("#####################################################")
    println("# twisting_circle ")
    twisting_circle()
    return true
end # function allrun

end # module


text:
using Ladder
using Losses
using DataFrames
using Gadfly

function boostingattack(n,k)

    y = rand(0:1.0,n)

    l = ladder(y)

    U = rand(0:1.0,n,k)
    ks = Float64[]
    ls = Float64[]
    goodis = Int[]

    for i in 1:k
        # kaggle score with 5 bits of precision
        kaggle_score = round(empirical_loss(y,U[:,i],loss01),5)
        push!(ks,kaggle_score)
        push!(ls,score!(l,Submission(U[:,i]),loss01))
    end

    # select all variables that lowered the score in the ladder mechanism
    if ls[1]< 0.5
        push!(goodis,1)
    end
    for i in 2:k
        if ls[i] < ls[i-1]
            push!(goodis,i)
        end
    end

    # select all positions with loss < 0.5
    bk = float(ks .<= 0.5)
    # select all positions that lowered the score
    bl = zeros(k)
    bl[goodis] = 1.0

    # compute boosted vectors
    ku = float(U * bk .>= countnz(bk)/2)
    lu = float(U * bl  .>= countnz(bl)/2)

    return empirical_loss(y,ku,loss01), empirical_loss(y,lu,loss01)

end

function plotboostingattack(n,k,r)

    x = [i for i in 10:10:k]
    kaggleavg = zeros(int(k/10))
    ladderavg = zeros(int(k/10))
    testavg = zeros(int(k/10))

    for _ in 1:r
        for i in 10:10:k
            ks,ls = boostingattack(n,i)
            kaggleavg[div(i,10)] += ks/r
            ladderavg[div(i,10)] += ls/r
            testavg[div(i,10)] += empirical_loss(rand(0:1.0,2*n),rand(0:1.0,2*n),loss01)/r
        end
    end

    df1 = DataFrame(x=x,y=kaggleavg,label="Kaggle")
    df2 = DataFrame(x=x,y=ladderavg,label="Ladder")
    df3 = DataFrame(x=x,y=testavg,label="Test")
    df = vcat(df1,df2,df3)

    kaggleplot = plot(df, x="x", y="y",color="label",Geom.point,Geom.line,
                       Scale.color_discrete_manual("red","blue","green"),
                       Coord.Cartesian(ymin=0.42,ymax=0.51),
                       Guide.yticks(ticks=[0.42,0.44,0.46,0.48,0.5]),
                       Guide.xlabel("Number of queries"), Guide.ylabel("Reported loss"),
                       Guide.title("Ladder vs Kaggle (normal precision)"),
                       Theme(minor_label_font_size=15pt,major_label_font_size=15pt,
                             key_label_font_size=14pt,key_title_font_size=15pt,key_position = :none))

    draw(PDF(string("kaggleboosting-",n,"-",k,"-",r,".pdf"), 6inch, 4.5inch),kaggleplot)
    writetable(string("kaggleboosting-",n,"-",k,"-",r,".csv"),df)

    return df
end

ext: jl


text:
using Pkg
Pkg.add("COSMO")
Pkg.add("PyCall")
using COSMO
using PyCall

max_stars_repo_name: oxfordcontrol/cosmo-python

max_issues_repo_name: oxfordcontrol/cosmo-python

max_forks_repo_name: oxfordcontrol/cosmo-python


text:
module EmptyModule

end


text:
#using FileIO
using Images
using Statistics

"""
Version Comment
v0.1	hf, use yen-threshold select threshold for 3D stack at each time point
v0.2    hf, add multithreads
"""

function create3dmask(zstack)
    #zstack = imfilter(zstack, Kernel.gaussian((2,2,2)))
	#z_depth = 20
    z_depth = size(zstack)[3]
    local mask = zeros(size(zstack))
    #thresholds_z = [real(yen_threshold(zstack[:, :, i])) for i in 1:20]
	local thresholds_z = zeros(Float64, z_depth)
	#@inbounds Threads.@threads for z in 1:z_depth
	@inbounds for z in 1:z_depth
		zstack[:,:,z] = imfilter(zstack[:, :, z], Kernel.gaussian(2))
		local select_pixels = zstack[81:end-81, 81:end-81,  z] .> 0.001
		thresholds_z[z]
		if sum(select_pixels) > 2e2 #TODO
			thresholds_z[z] = real(otsu_threshold(
				zstack[81:end-81,81:end-81,z][select_pixels]))
		end
	end
    #threshold_3d = real(yen_threshold(imfilter(zstack, Kernel.gaussian((2,2,1)))))
    #threshold_3d = real(yen_threshold(zstack))
    threshold_3d = maximum(thresholds_z) # choose clearest one
    #mask = opening(zstack .> threshold_3d)
    mask = zstack .> threshold_3d
    #remove_small_area!(mask)
    mask, mask_size = remove_small_area(mask)
    return mask, mask_size, threshold_3d
end

function extract3dnucleus(stack)
    #z_depth = 20
    #t_len = size(stack)[3] ÷z_depth
    z_depth,t_len = size(stack)[3:4]
    local nucleus = zeros(Gray{Normed{UInt16,16}}, size(stack))
	local nucleus_len = length(stack)
    local thresholds = zeros(Float64, t_len)
    local mask_size = zeros(Float64, t_len)
    #nucleus_3dmask = zeros(size(stack)[1], size(stack)[2], z_depth)
	println("Extracting nucleus")
    @inbounds Threads.@threads for t in 1:t_len
		#local z = (t-1)*20+1 : 20*t
		nucleus_3dmask, mask_size[t], thresholds[t] = create3dmask(stack[:, :, :, t])
		nucleus_t = view(nucleus, :, :, :, t)
		stack_t = view(stack, :, :, :, t)
        #nucleus[:, :, z] = nucleus_3dmask .* stack[:, :, z]
		for i in eachindex(nucleus_3dmask)
			if nucleus_3dmask[i]
				nucleus_t[i] = stack_t[i]
			end
		end
    end
    nucleus, mask_size, thresholds
end

"""
Just remove small regions in binary images
"""
#function remove_small_area!(mask)
function remove_small_area(mask)
    mask_con = label_components(mask);
    #mask_res = BitArray(undef, size(mask));
    #mask_res .= false;
    con_size = component_lengths(mask_con)
    selected_con = (0:length(con_size)-1)[con_size .> 2e4]
    #@inbounds Threads.@threads for i in 1:length(mask)
	local mask_size = 0
    @inbounds for i in 1:length(mask)
        if mask_con[i] ∉ selected_con
            mask[i] = false
		else
			mask_size += 1
        end
    end
    mask, mask_size
end

"""
Computes threshold for grayscale image using Otsu's method. (modify from Imagas.jl)
Parameters:
-    img         = Grayscale input image
-    bins        = Number of bins used to compute the histogram. Needed for floating-point images.
"""
function otsu_threshold2(img::AbstractArray{T, N}, min, bins::Int = 256) where {T<:Union{Gray,Real}, N}

    #min, max = extrema(img)
    max = maximum(img)
	if max == 0 # avoid image without cell
		return T(0)
	end
    edges, counts = imhist(img, range(gray(min), stop=gray(max), length=bins))
	if sum(counts)<2e2 # don't count small region
		return T(0)
	end
    histogram = counts./sum(counts)

    ω0 = 0
    μ0 = 0
    μt = 0
    μT = sum((1:(bins+1)).*histogram)
    max_σb=0.0
    thres=1

    for t in 1:bins
        ω0 += histogram[t]
        ω1 = 1 - ω0
        μt += t*histogram[t]

        σb = (μT*ω0-μt)^2/(ω0*ω1)

        if(σb > max_σb)
            max_σb = σb
            thres = t
        end
    end

    return T((edges[thres-1]+edges[thres])/2)
end

#s3c2 = load("../mRNA_confocal_hamamatsu-60X-TIRF/20200316_result/s3_c2.tiff");
#@time nucleus_all, threshold_all = extract3dnucleus(s3c2);
#@time save(File(format"TIFF", "s5-c2_clear.ome.tiff"), N0f16.(nucleus_all))
#x, y, z_all = size(nucleus_all)
#@time embedxml(x, y, 20, z_all÷20, "s5-c2_clear.ome.tiff")


text:
#!/usr/bin/env julia
# Python and Advanced Data Structures
# https://teaching.healthtech.dtu.dk/22110/index.php/Python_and_Advanced_Data_Structures

# Ex 1
function read_experiment!(outdict, f)
    for line = eachline(f)
        line_split = split(line, "\t")
        key = string(line_split[1])
        nums = [parse(Float64, x) for x in line_split[2:end]]
        outdict[key] = append!(get(outdict, key, Float64[]), nums)
    end
end

function read_experiments(out, files...)
    expdict = Dict{String,Array{Float64,1}}()
    for file in files
        open(f->read_experiment!(expdict, f), file)
    end
    for (key, value) in sort(collect(expdict))
        print(out, key, "\t")
        println(out, join(value, "\t"))
    end
end

open("data/test.dat", "w") do f
    read_experiments(f, "data/test1.dat", "data/test2.dat", "data/test3.dat")
end

# Ex 2

matrix = readdlm("data/matrix.dat", '\t', Int, '\n')
#matrix1 = readdlm("data/matrix1", '\t', Int, '\n')
#matrix2 = readdlm("data/matrix2", '\t', Int, '\n')
#dot = matrix1*matrix2
transpose(matrix)


text:
"""
    e = fiedler(A;  nev = 1, tol = 0.0 )

Compute the Fiedler vector and value of the graph with adjacency matrix `A`.
Computes them using `eigs` from Arpack.jl, and returns `e` in the same format.
So, `e[1]` is the fiedler value, and `e[2]` is the vector.
If `nev` is more than 1, then it returns `nev` values and vectors.

Note that for small matrices, or for well-conditioned ones, it will be faster to simply use
~~~
    e = eigs(lap(a), which = :SM, nev = 1)
~~~
"""
function fiedler(a::AbstractArray; nev = 1, tol = 0.0)
    f = approxchol_lap(a)
    op = Laplacians.SqLinOp(true,1.0,size(a,1),f)
    e = eigs(op, which=:LM, nev=nev, tol=tol)
    e[1] .= 1 ./ e[1]
    return e
end


text:
function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in kwargs if k in keys)
end


text:
module PRM

    using Plots

    export run_sim, SimParam, SimState, SimOutput

    Base.@kwdef struct SimParam
        Tm::Float32 # Time constant of the PT1 Block [s]
        epsilon::Float32 # Local tolerance of update step size function
        t0::Float32 # Simulation start time [s]
        tf::Float32 # Simulation end time [s]
        ts::Float32 # Time at which input goes high
        x0::Float32 # Initial state of PT1-Block
        h0::Float32 # Intial step size [s]
        u0::Float32 # Rising edge size
    end

    mutable struct SimState
        x::Float32 # Current state of PT1 Block
        u1::Float32 # Current input
        u2::Float32
        u3::Float32
        y::Array{Float32} # Current output
        d::Float32 # Current local discretization error
        h::Float32 # Current step size
        t::Float32 # Current time
        i::Int32 # Current loop iteration
        hys_last_run::Int32 # Hysteresis memory
        hys_last_loop::Int32
    end
    SimState() = SimState(0.0, 0.0, 0.0, 0.0, [0.0, 0.0, 0.0], 0.0, 0.0, 0.0, 0, 0, 0)

    struct SimOutput{T}
        x_values::Vector{T} # State of PT1-Block
        u_values::Vector{T} # Input
        y_values::Vector{T} # Outputs of all blocks
        d_values::Vector{T} # Local Discretization Error
        t_values::Vector{T} # Time
        h_values::Vector{T} # Step size
    end
    SimOutput(T) = SimOutput{T}(T[], T[], T[], T[], T[], T[])
    SimOutput() = SimOutput(Float32)

    # PT1 Block
    function pt1(x, u)
        Tm = 10 # Time constant
        xdot = -1/Tm * x + 1/Tm * u # State equation
        return (; xdot, x)
    end

    # Hysteresis block
    function hysteresis(u, hys_last)
        h_e = 0.085 # Rising edge limit
        h_a = 0.065 # Falling edge limit

        if u >= h_e
            y = 1
        elseif u <= -h_e
            y = -1
        elseif -h_a <= u <= h_a
            y = 0
        end

        if h_a < u < h_e
            if hys_last == 1
                y = 1
            else
                y = 0
            end
        elseif -h_e < u < -h_a
            if hys_last == -1
                y = 1
            else
                y = 0
            end
        end
        return y
    end

    function step(u, t, t_step)
        if t < t_step
            return 0
        else
            return u
        end
    end

    function sim_prm!(x, u, t, state::SimState, output::Bool)
        # Block 3, PT1, only output
        _, y3 = pt1(x, 0)
        # Block 1, Subtractor
        y1 = -(u, y3)
        # Block 2, Hysteresis
        y2 = hysteresis(y1, state.hys_last_run)
        state.hys_last_run = y2
        # Block 3, PT1, only state
        xdot, _ = pt1(x, y2)

        if output
            state.y[1] = y1
            state.y[2] = y2
            state.y[3] = y3
        end

        return xdot
    end

      # Integration method
      function midpoint_method!(sim_topology::Function, state::SimState)
        x = state.x
        t = state.t
        h = state.h

        k1 = sim_topology(x, state.u1, t, state, true)
        k2 = sim_topology(x + h/2*k1, state.u2, t + h/2, state, false)
        k3 = sim_topology(x - h*k1 + 2*h*k2, state.u3, t + h, state, false)
        x = x + h*k2
        d = h/6*(k1 - 2*k2 + k3)

        state.x = x
        state.d = d
    end

    function log_output!(state::SimState, output::SimOutput)
        append!(output.u_values, state.u1)
        append!(output.x_values, state.x)
        append!(output.y_values, state.y)
        append!(output.t_values, state.t)
        append!(output.d_values, state.d)
        append!(output.h_values, state.h)
    end

    function render_output(output::SimOutput)
        y_values = reshape(output.y_values, (3, :))

        p1 = plot(output.t_values, output.u_values, label="Input")
        p2 = plot(output.t_values, y_values[1, :], label="Output Subtractor")
        p3 = plot(output.t_values, y_values[2, :], label="Output Hysteresis")
        p4 = plot(output.t_values, y_values[3, :], label="Output PT1")
        p5 = plot(output.t_values, output.h_values, label="Step Size")
        p6 = plot(output.t_values, output.d_values, label="LDE")

        plot(p1, p2, p3, p4, p5, p6, layout=(3, 2), xlabel="Time [s]")
    end

    function run_sim(simp::SimParam)
        output = SimOutput();
        state = SimState();

        state.x = simp.x0
        state.h = simp.h0
        state.i = 1

        state.t = simp.t0

        h_min = 12*simp.epsilon
        h_max = simp.Tm*2

        while state.t < simp.tf

            state.u1 = step(simp.u0, state.t, simp.ts)
            state.u2 = step(simp.u0, state.t + state.h/2, simp.ts)
            state.u3 = step(simp.u0, state.t + state.h, simp.ts)

            state.hys_last_run = state.hys_last_loop
            midpoint_method!(sim_prm!, state)

            if abs(state.d) > 0
                h_new = state.h*(simp.epsilon/abs(state.d))^(1/3)

                h_new = min(max(h_new, h_min), 0.99*h_max)

                if h_new > 2*state.h
                    state.h = h_new
                elseif h_new <= state.h
                    state.h = 0.75*h_new
                    continue
                end
            end

            state.hys_last_loop = state.hys_last_run
            log_output!(state, output)

            state.t = state.t + state.h
            state.i = state.i + 1
        end

        render_output(output)

    end

end # module
simp = PRM.SimParam(Tm=10, epsilon=1e-10, x0=0, h0=0.01, t0=0, tf=20, ts=1, u0=0.17);

PRM.run_sim(simp)

max_stars_repo_name: VeroAsghar/PRM

max_issues_repo_name: VeroAsghar/PRM

max_forks_repo_name: VeroAsghar/PRM


text:
@testset "basic" begin
  m = Chain(Conv((3, 3), 3 => 16), Conv((3, 3), 16 => 32))
  @test outputsize(m, (10, 10, 3, 1)) == (6, 6, 32, 1)

  m = Dense(10, 5)
  @test_throws DimensionMismatch outputsize(m, (5, 2)) == (5, 1)
  @test outputsize(m, (10,); padbatch=true) == (5, 1)

  m = Chain(Dense(10, 8, σ), Dense(8, 5), Dense(5, 2))
  @test outputsize(m, (10,); padbatch=true) == (2, 1)
  @test outputsize(m, (10, 30)) == (2, 30)

  m = Chain(Dense(10, 8, σ), Dense(8, 4), Dense(5, 2))
  @test_throws DimensionMismatch outputsize(m, (10,))

  m = Flux.Diagonal(10)
  @test outputsize(m, (10, 1)) == (10, 1)

  m = Maxout(() -> Conv((3, 3), 3 => 16), 2)
  @test outputsize(m, (10, 10, 3, 1)) == (8, 8, 16, 1)

  m = flatten
  @test outputsize(m, (5, 5, 3, 10)) == (75, 10)

  m = Flux.unsqueeze(3)
  @test outputsize(m, (5, 7, 13)) == (5, 7, 1, 13)

  m = Flux.Bilinear(10, 10, 7)
  @test outputsize(m, (10,)) == (7,)
  @test outputsize(m, (10, 32)) == (7, 32)

  m = Chain(Conv((3, 3), 3 => 16), BatchNorm(16), flatten, Dense(1024, 10))
  @test outputsize(m, (10, 10, 3, 50)) == (10, 50)
  @test outputsize(m, (10, 10, 3, 2)) == (10, 2)

  m = SkipConnection(Conv((3, 3), 3 => 16; pad = 1), (mx, x) -> cat(mx, x; dims = 3))
  @test outputsize(m, (10, 10, 3, 1)) == (10, 10, 19, 1)

  m = Parallel((mx, x) -> cat(mx, x; dims = 3), Conv((3, 3), 3 => 16; pad = 1), identity)
  @test outputsize(m, (10, 10, 3, 1)) == (10, 10, 19, 1)
end

@testset "multiple inputs" begin
  m = Parallel(vcat, Dense(2, 4, relu), Dense(3, 6, relu))
  @test outputsize(m, (2,), (3,)) == (10,)
  @test outputsize(m, ((2,), (3,))) == (10,)
  @test outputsize(m, (2,), (3,); padbatch=true) == (10, 1)
  @test outputsize(m, (2,7), (3,7)) == (10, 7)

  m = Chain(m, Dense(10, 13, tanh), softmax)
  @test outputsize(m, (2,), (3,)) == (13,)
  @test outputsize(m, ((2,), (3,))) == (13,)
  @test outputsize(m, (2,), (3,); padbatch=true) == (13, 1)
  @test outputsize(m, (2,7), (3,7)) == (13, 7)
end

@testset "activations" begin
  @testset for f in [celu, elu, gelu, hardsigmoid, hardtanh,
                     leakyrelu, lisht, logcosh, logσ, mish,
                     relu, relu6, rrelu, selu, σ, softplus,
                     softshrink, softsign, swish, tanhshrink, trelu]
    @test outputsize(Dense(10, 5, f), (10, 1)) == (5, 1)
  end
end

@testset "conv" begin
  m = Conv((3, 3), 3 => 16)
  @test outputsize(m, (10, 10, 3, 1)) == (8, 8, 16, 1)
  m = Conv((3, 3), 3 => 16; stride = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (2, 2, 16, 1)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 3, 1)) == (5, 5, 16, 1)
  m = Conv((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (4, 4, 16, 1)
  @test_throws DimensionMismatch outputsize(m, (5, 5, 2))
  @test outputsize(m, (5, 5, 3, 100)) == (4, 4, 16, 100)

  m = ConvTranspose((3, 3), 3 => 16)
  @test outputsize(m, (8, 8, 3, 1)) == (10, 10, 16, 1)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2)
  @test outputsize(m, (2, 2, 3, 1)) == (5, 5, 16, 1)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 3, 1)) == (5, 5, 16, 1)
  m = ConvTranspose((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test outputsize(m, (4, 4, 3, 1)) == (5, 5, 16, 1)

  m = DepthwiseConv((3, 3), 3 => 6)
  @test outputsize(m, (10, 10, 3, 1)) == (8, 8, 6, 1)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (2, 2, 6, 1)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 3, 1)) == (5, 5, 6, 1)
  m = DepthwiseConv((3, 3), 3 => 6; stride = 2, pad = 3, dilation = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (4, 4, 6, 1)

  m = CrossCor((3, 3), 3 => 16)
  @test outputsize(m, (10, 10, 3, 1)) == (8, 8, 16, 1)
  m = CrossCor((3, 3), 3 => 16; stride = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (2, 2, 16, 1)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 3, 1)) == (5, 5, 16, 1)
  m = CrossCor((3, 3), 3 => 16; stride = 2, pad = 3, dilation = 2)
  @test outputsize(m, (5, 5, 3, 1)) == (4, 4, 16, 1)

  m = AdaptiveMaxPool((2, 2))
  @test outputsize(m, (10, 10, 3, 1)) == (2, 2, 3, 1)

  m = AdaptiveMeanPool((2, 2))
  @test outputsize(m, (10, 10, 3, 1)) == (2, 2, 3, 1)

  m = GlobalMaxPool()
  @test outputsize(m, (10, 10, 3, 1)) == (1, 1, 3, 1)

  m = GlobalMeanPool()
  @test outputsize(m, (10, 10, 3, 1)) == (1, 1, 3, 1)

  m = MaxPool((2, 2))
  @test outputsize(m, (10, 10, 3, 1)) == (5, 5, 3, 1)
  m = MaxPool((2, 2); stride = 1)
  @test outputsize(m, (5, 5, 4, 1)) == (4, 4, 4, 1)
  m = MaxPool((2, 2); stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 2, 1)) == (5, 5, 2, 1)

  m = MeanPool((2, 2))
  @test outputsize(m, (10, 10, 3, 1)) == (5, 5, 3, 1)
  m = MeanPool((2, 2); stride = 1)
  @test outputsize(m, (5, 5, 4, 1)) == (4, 4, 4, 1)
  m = MeanPool((2, 2); stride = 2, pad = 3)
  @test outputsize(m, (5, 5, 2, 1)) == (5, 5, 2, 1)
end

@testset "normalisation" begin
  m = Dropout(0.1)
  @test outputsize(m, (10, 10)) == (10, 10)
  @test outputsize(m, (10,); padbatch=true) == (10, 1)

  m = AlphaDropout(0.1)
  @test outputsize(m, (10, 10)) == (10, 10)
  @test outputsize(m, (10,); padbatch=true) == (10, 1)

  m = LayerNorm(32)
  @test outputsize(m, (32, 32, 3, 16)) == (32, 32, 3, 16)
  @test outputsize(m, (32, 32, 3); padbatch=true) == (32, 32, 3, 1)

  m = BatchNorm(3)
  @test outputsize(m, (32, 32, 3, 16)) == (32, 32, 3, 16)
  @test outputsize(m, (32, 32, 3); padbatch=true) == (32, 32, 3, 1)

  m = InstanceNorm(3)
  @test outputsize(m, (32, 32, 3, 16)) == (32, 32, 3, 16)
  @test outputsize(m, (32, 32, 3); padbatch=true) == (32, 32, 3, 1)

  m = GroupNorm(16, 4)
  @test outputsize(m, (32, 32, 16, 16)) == (32, 32, 16, 16)
  @test outputsize(m, (32, 32, 16); padbatch=true) == (32, 32, 16, 1)
end


text:
include("utils/curves.jl")
include("utils/dict_utils.jl")
include("utils/rng_number.jl")
include("utils/zero_rates.jl")

include("utils/mc_mode.jl")

include("utils/mc_config.jl")
include("utils/sim_result.jl")
include("utils/underlying.jl")


text:
using MRIReco: Profile

function rawSignalToISMRMRD(signal,seq;
                                phantom=Phantom(name="Phantom",x=[0]),
                                sys=Scanner(),
                                simParams=Dict{String,Any}())
    version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"]))
    Nz = seq.DEF["Nz"]
    #Number of samples and FOV
    _, ktraj = get_kspace(seq) #kspace information
    mink = minimum(ktraj, dims=1)
    maxk = maximum(ktraj, dims=1)
    Wk = maxk .- mink
    Δx = 1 ./ Wk[1:2] #[m] Only x-y
    if haskey(seq.DEF, "FOV")
        FOVx, FOVy, _ = seq.DEF["FOV"] #[m]
        if FOVx > 1 FOVx *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
        if FOVy > 1 FOVy *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
        Nx = round(Int64, FOVx / Δx[1])
        Ny = round(Int64, FOVy / Δx[2])
    else
        Nx = seq.DEF["Nx"]
        Ny = seq.DEF["Ny"]
        FOVx = Nx * Δx[1]
        FOVy = Ny * Δx[2]
    end
    #It needs to be transposed for the raw data
    ktraj = transpose(ktraj) ./ maximum(2*abs.(ktraj[:]))

    #First we define the ISMRMRD data XML header
    params = Dict(
        #Simulator
        "systemVendor"                   => "Koma.jl", #String
        "systemModel"                    => "v"*version, #String
        "institutionName"                => "Pontificia Universidad Catolica de Chile", #String
        #Phantom
        "patientName"                    => phantom.name,
        #Sys
        "systemFieldStrength_T"          => sys.B0, #Float
        "H1resonanceFrequency_Hz"        => floor(Int64, γ * sys.B0), #Long (Int)
        #Seq
        "trajectory"                     => "other", #Must be: cartesian, epi, radial, goldenangle, spiral, and other
        "protocolName"                   => haskey(seq.DEF,"Name") ? seq.DEF["Name"] : "NoName", #String
        # "trajectoryDescription"          => Dict{String, Any}("comment"=>""), #You can put wathever you want here: comment, bandwidth, MaxGradient_G_per_cm, MaxSlewRate_G_per_cm_per_s, interleaves, etc
        "userParameters"                 => simParams, #Dict with parameters
        #encoding>
        "encodedFOV"                     => [FOVx, FOVy, 1.0],    #encodedSpace>fieldOfView_mm
        "reconFOV"                       => [FOVx, FOVy, 1.0],    #reconSpace>fieldOfView_mm
        "encodedSize"                    => [Nx, Ny, 1],            #encodedSpace>matrixSize
        "reconSize"                      => [Nx+Nx%2, Ny+Ny%2, 1],  #reconSpace>matrixSize
        #encodingLimits>
        "enc_lim_kspace_encoding_step_1" => Limit(0, Nx-1, ceil(Int, Nx / 2)),   #min, max, center, e.g. phase encoding line number
        "enc_lim_kspace_encoding_step_2" => Limit(0, 0, 0),     #min, max, center, e.g. partition encoding number
        "enc_lim_average"                => Limit(0, 0, 0),     #min, max, center, e.g. signal average number
        "enc_lim_slice"                  => Limit(0, 0, 0),     #min, max, center, e.g. imaging slice number
        "enc_lim_contrast"               => Limit(0, 0, 0),     #min, max, center, e.g. echo number in multi-echo
        "enc_lim_phase"                  => Limit(0, 0, 0),     #min, max, center, e.g. cardiac phase number
        "enc_lim_repetition"             => Limit(0, 0, 0),     #min, max, center, e.g. dynamic number for dynamic scanning
        "enc_lim_set"                    => Limit(0, 0, 0),     #min, max, center, e.g. flow encoding set
        "enc_lim_segment"                => Limit(0, 0, 0),     #min, max, center, #segment: e.g. segment number for segmented acquisition
        #sequenceParameters
        # "TR"                             => 0,
        # "TE"                             => 0,
        # "TI"                             => 0,
        # "flipAngle_deg"                  => 0,
        # "echo_spacing"                   => 0,
    )

    #Then, we define the Profiles
    profiles = Profile[]
    t_acq = get_sample_times(seq)
    Nadcs = sum(is_ADC_on.(seq))
    NadcsPerImage = floor(Int, Nadcs / Nz)
    scan_counter = 0
    nz = 0
    current = 1
    for s = seq #Iterate over sequence blocks
        if is_ADC_on(s)
            Nsamples = s.ADC.N[1]
            Δt_us = float32( s.ADC.T[1] / (Nsamples - 1) * 1e6 )
            t0_us = floor(Int32, t_acq[current]*1e6 )
            #Header of profile data, head::AcquisitionHeader
            head = AcquisitionHeader(
                0, #version: First unsigned int indicates the version
                0, #flags: bit field with flags, noise measure, calibration, coil sens, measure, etc
                0, #measurement_uid: Unique ID for the measurement
                scan_counter, #scan_counter: Current acquisition number in the measurement
                t0_us, #acquisition_time_stamp: Acquisition clock, I am "miss"-using this variable to store t0 in us
                (0, 0, 0), #physiology_time_stamp: Physiology time stamps, e.g. ecg, breating, etc.
                Nsamples, #number_of_samples
                1, #available_channels: Available coils
                1, #active_channels: Active coils on current acquisiton
                Tuple(0 for i=1:16), #channel_mask: Active coils on current acquisiton
                0, #discard_pre: Samples to be discarded at the beginning of acquisition
                0, #discard_post: Samples to be discarded at the end of acquisition
                0, #center_sample: Sample at the center of k-space
                0, #encoding_space_ref: Reference to an encoding space, typically only one per acquisition
                2, #trajectory_dimensions: Indicates the dimensionality of the trajectory vector (0 means no trajectory)
                Δt_us, #sample_time_us: Time between samples in micro seconds, sampling BW
                (0.0f0, 0.0f0, 0.0f0), #position: Three-dimensional spatial offsets from isocenter
                (1.0f0, 0.0f0, 0.0f0), #read_dir: Directional cosines of the readout/frequency encoding
                (0.0f0, 1.0f0, 0.0f0), #phase_dir: Directional cosines of the phase
                (0.0f0, 0.0f0, 1.0f0), #slice_dir: Directional cosines of the slice direction
                (0.0f0, 0.0f0, 0.0f0), #patient_table_position: Patient table off-center
                EncodingCounters( #idx: Encoding loop counters
                    scan_counter, #kspace_encode_step_1: e.g. phase encoding line number
                    0, #kspace_encode_step_2: e.g. partition encoding number
                    0, #average: e.g. signal average number
                    nz, #slice: e.g. imaging slice number
                    0, #contrast: e.g. echo number in multi-echo
                    0, #phase: e.g. cardiac phase number
                    0, #repetition: e.g. dynamic number for dynamic scanning
                    0, #set: e.g. flow encoding set
                    0, #segment: e.g. segment number for segmented acquisition
                    (0, 0, 0, 0, 0, 0, 0, 0) #user: Free user parameters
                ),
                (0, 0, 0, 0, 0, 0, 0, 0), #user_int: Free user parameters
                (0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0) #user_float: Free user parameters
            )
            #Trajectory information, traj::Array{Float32,2}, 1dim=DIM, 2dim=numsaples
            traj = ktraj[1:2, current:current+Nsamples-1]
            #Acquired data, data::Array{Complex{Float32},2}, 1dim=numsamples, 2dim=coils
            dat =  signal[current:current+Nsamples-1, :]
            #Saving profile
            push!(profiles, MRIReco.Profile(head, traj, dat))
            #Update counters
            scan_counter += 1
            current += Nsamples
            if scan_counter % NadcsPerImage == 0 #For now only Nz is considered
                nz += 1 #another image
                scan_counter = 0 #reset counter
            end
        end
    end

    return RawAcquisitionData(params, profiles)
end
# enum ISMRMRD_AcquisitionFlags {
#     ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP1               =  1,
#     ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1                =  2,
#     ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP2               =  3,
#     ISMRMRD_ACQ_LAST_IN_ENCODE_STEP2                =  4,
#     ISMRMRD_ACQ_FIRST_IN_AVERAGE                    =  5,
#     ISMRMRD_ACQ_LAST_IN_AVERAGE                     =  6,
#     ISMRMRD_ACQ_FIRST_IN_SLICE                      =  7,
#     ISMRMRD_ACQ_LAST_IN_SLICE                       =  8,
#     ISMRMRD_ACQ_FIRST_IN_CONTRAST                   =  9,
#     ISMRMRD_ACQ_LAST_IN_CONTRAST                    = 10,
#     ISMRMRD_ACQ_FIRST_IN_PHASE                      = 11,
#     ISMRMRD_ACQ_LAST_IN_PHASE                       = 12,
#     ISMRMRD_ACQ_FIRST_IN_REPETITION                 = 13,
#     ISMRMRD_ACQ_LAST_IN_REPETITION                  = 14,
#     ISMRMRD_ACQ_FIRST_IN_SET                        = 15,
#     ISMRMRD_ACQ_LAST_IN_SET                         = 16,
#     ISMRMRD_ACQ_FIRST_IN_SEGMENT                    = 17,
#     ISMRMRD_ACQ_LAST_IN_SEGMENT                     = 18,
#     ISMRMRD_ACQ_IS_NOISE_MEASUREMENT                = 19,
#     ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION             = 20,
#     ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING = 21,
#     ISMRMRD_ACQ_IS_REVERSE                          = 22,
#     ISMRMRD_ACQ_IS_NAVIGATION_DATA                  = 23,
#     ISMRMRD_ACQ_IS_PHASECORR_DATA                   = 24,
#     ISMRMRD_ACQ_LAST_IN_MEASUREMENT                 = 25,
#     ISMRMRD_ACQ_IS_HPFEEDBACK_DATA                  = 26,
#     ISMRMRD_ACQ_IS_DUMMYSCAN_DATA                   = 27,
#     ISMRMRD_ACQ_IS_RTFEEDBACK_DATA                  = 28,
#     ISMRMRD_ACQ_IS_SURFACECOILCORRECTIONSCAN_DATA   = 29,
#     ISMRMRD_ACQ_COMPRESSION1                        = 53,
#     ISMRMRD_ACQ_COMPRESSION2                        = 54,
#     ISMRMRD_ACQ_COMPRESSION3                        = 55,
#     ISMRMRD_ACQ_COMPRESSION4                        = 56,
#     ISMRMRD_ACQ_USER1                               = 57,
#     ISMRMRD_ACQ_USER2                               = 58,
#     ISMRMRD_ACQ_USER3                               = 59,
#     ISMRMRD_ACQ_USER4                               = 60,
#     ISMRMRD_ACQ_USER5                               = 61,
#     ISMRMRD_ACQ_USER6                               = 62,
#     ISMRMRD_ACQ_USER7                               = 63,
#     ISMRMRD_ACQ_USER8                               = 64
# };


text:
#= using Genie =#

module tmp

const API_VERSION = "0.0.1" # Change this to invalidate HTTP cache

import Genie
using Genie.Router: route, @params
using Genie.Requests: getpayload

# Generates HTML responses
using Genie.Renderer: html

using Memoize: @memoize
using ProgressMeter: @showprogress

using Base.Iterators: flatten, product
using Statistics

import HTTP
import GeoJSON
import Turf
import VegaLite

VegaLite.actionlinks(false) # Global setting - disable action button on all plots

# This converts its argument to json and sets the appropriate headers for content type
# We're customising it to set the CORS header
json(data; status::Int = 200) =
    Genie.Renderer.json(data; status = status, headers = Dict(
        "Access-Control-Allow-Origin" => "*",
        "Cache-Control" => "public, max-age=$(365 * 24 * 60 * 60)", # cache for a year (max recommended). Change API_VERSION to invalidate
    ))

Genie.config.session_auto_start = false
# Default headers are supposed to go here, but they don't seem to work.
#= Genie.config.cors_headers["Access-Control-Allow-Origin"] = "*" =#

# Test route
route("/") do
    json(Dict("Answer" => 42))
end

# Used solely to invalidate caches - increment API_VERSION if you need to do so
route("/version") do
    Genie.Renderer.json(Dict("version" => API_VERSION); status = 200, headers = Dict(
        "Access-Control-Allow-Origin" => "*",
        "Cache-Control" => "max-age=0", # Don't cache
    ))
end

### APP ###

## State (all constant)

include("$(@__DIR__)/scenarios.jl")

const zones = GeoJSON.parsefile("$(@__DIR__)/../data/geometry/zones.geojson")

const zone_centroids = Array{Array{Float64,1},1}(undef,length(zones.features))
for f in zones.features
    zone_centroids[f.properties["fid"]] = Turf.centroid(f.geometry).coordinates
end

const NUM_ZONES = zone_centroids |> length

# This should probably be reloaded periodically so server doesn't need to be restarted?
const links, mats, metadata = load_scenarios(packdir)

const NUM_LINKS = (links |> first)[2] |> size |> first

## Funcs

# Get list of scenarios (scenario name, id, years active)
list_scenarios() = metadata["scenarios"]

# Get list of variable names and ids
function list_variables(domain)
    if domain == "links"
        metadata["links"]["columns"]
    elseif domain == "od_matrices"
        metadata["od_matrices"]["columns"]
    else
        throw(DomainError(domain, "is invalid"))
    end
end

function mat_data(scenario, year, variable)
    try
        mats[(scenario, year)][variable]
    catch e
        if e isa KeyError
            @warn e
            return fill(missing, NUM_ZONES, NUM_ZONES)
        else
            rethrow()
        end
    end
end

function link_data(scenario, year, variable)
    try
        links[(scenario, year)][!, Symbol(variable)]
    catch e
        if e isa KeyError
            @warn e
            return fill(missing, NUM_LINKS)
        else
            rethrow()
        end
    end
end

# Get colour to draw each shape for (scenario, year, variable, comparison scenario, comparison year) -> (colours)
mat_comp(args...; kwargs...) = comp(mat_data, args...; kwargs...)
link_comp(args...; kwargs...) = comp(link_data, args...; kwargs...)
function comp(data, scenario, year, variable, comparison_scenario, comparison_year;percent=true)

    main = data(scenario, year, variable)
    comparator = data(comparison_scenario, comparison_year, variable)

    # dims = 2 sums rows; dims = 1 sums cols
    # Returns Array{Float,2} but we want Array{Float,1} so we flatten it
    if percent
        # Avoid NaN from dividing by zero.
        # Skip missing values here because we don't want eps to be `missing`
        # and poison every value.
        eps = mean(sum(skipmissing(vec(main)))) / 10000

        result = (sum(main, dims = 2) .+ eps) ./ (sum(comparator, dims = 2) .+ eps) .- 1
    else
        result = sum(main .- comparator, dims = 2)
    end

    return collect(flatten(result))
end

"""
    var_stats(domain, variable, quantiles, percent=false)

Return quantiles to use when comparing `variable` between scenarios.

We approximate the quantiles for all pair-wise differences with this algorithm:

1. Compute the element-wise difference between each pair of scenarios for the
given `domain` and `variable`
2. sample from these differences to create a vector of differences
3. return the quantiles of the sampled vector
"""
# This is expensive, so it is memoized and the cache is warmed up before the
# webserver starts.
@memoize function var_stats(domain, variable, quantiles, percent=false)
    vars = []
    if domain == "od_matrices"
        vars = [scen[variable] for scen in values(mats)]
    elseif domain == "links"
        vars = [df[!, Symbol(variable)] for df in values(links)]
    end

    # Choose what type of diff to calculate.
    function percent_diff(l, r)
        # dims = 2 sums rows; dims = 1 sums cols
        eps = mean(sum(vars[l], dims = 2)) / 10000
        result = (sum(vars[l], dims = 2) .+ eps) ./ (sum(vars[r], dims = 2) .+ eps) .- 1
    end
    function absolute_diff(l, r)
        vars[l] .- vars[r]
    end
    diff = percent ? percent_diff : absolute_diff

    # Sample the differences between the values of this variable for each pair
    # of scenarios.
    pairs = product(1:length(vars), 1:length(vars))
    diffs = [rand(diff(a, b), 10000) for (a, b) in pairs]

    return quantile(flatten(vcat(diffs...)), quantiles)
end

"""
    var_stats_1d(domain, variable, quantiles)

Collect all values for `variable` from all scenarios in `domain` into a vector
and return the quantiles of that vector.

This is a bit expensive, so it is memoized.
"""
# memoized functions can't have docstrings - don't delete this comment
@memoize function var_stats_1d(domain, variable, quantiles)
    # dims = 2 sums rows; dims = 1 sums cols
    vars = []
    if domain == "od_matrices"
        vars = [get(v,variable,[]) for (k,v) in mats]
    elseif domain == "links"
        vars = [get(v,Symbol(variable),[]) for (k,v) in links]
    end
    vcat(vars...) |> flatten |> x -> quantile(x,quantiles)
end

"Given a pair of variable_name => meta, return true if it should be used"
is_used((name, meta)) = get(meta, "use", true)

println("Warming up the cache: links")
@showprogress for variable in keys(filter(is_used, list_variables("links")))
    # get these quantiles from colourMap in index.js
    var_stats("links", variable, (0.1, 0.9))
    var_stats("links", variable, (0.05, 0.95), true)
end

println("Warming up the cache: matrices")
@showprogress for variable in keys(filter(is_used, list_variables("od_matrices")))
    # get these quantiles from colourMap in index.js
    var_stats("od_matrices", variable, (0.0001, 0.9999))
    var_stats("od_matrices", variable, (0.05, 0.95), true)
end

## Routes

route("/scenarios") do
    list_scenarios() |> json
end

route("/variables/:domain") do
    list_variables(@params(:domain)) |> json
end

# Return a vector: [low_quantile, high_quantile]
route("/stats") do
    # Was removed this in ac81797 but is 'needed' below
    defaults = Dict(
        :domain => "od_matrices",
        :comparewith => "DoMin",
        :variable => "Total_GHG",
        :percent => "true",
        :quantiles => "0.0001,0.9999",  # These default percentiles seem v. generous
                                        # but we look at all possible differences
                                        # so overwhelming majority of differences are
                                        # tiny. Seems to work OK in practice.
    )
    d = merge(defaults, getpayload())
    quantiles = parse.(Float64,split(d[:quantiles],","))
    if d[:domain] in ["od_matrices", "links"]
        d[:comparewith] != "none" && return var_stats(d[:domain],d[:variable],Tuple(quantiles),d[:percent]=="true") |> json
        # TODO: If metadata holds the quantiles return them here below
        return var_stats_1d(d[:domain],d[:variable],Tuple(quantiles)) |> json
    else
        throw(DomainError(:domain))
    end
end

# Return a vector of data, usually floats, for the requested scenario,
# comparison parameters, and so on.
route("/data") do
    defaults = Dict(
        :domain => "od_matrices",
        :scenario => "Rail",
        :year => "2030",
        :comparewith => "DoMin",
        :compareyear => "auto",
        :variable => "Total_GHG",
        :percent => "true",
        :row => "false",
    )
    d = merge(defaults, getpayload())
    scenario = d[:scenario]
    year = parse(Int, d[:year])
    variable = d[:variable]
    compareyear = d[:compareyear] == "auto" ? year : parse(Int, d[:compareyear])
    percent = d[:percent] == "true"
    comparewith = d[:comparewith]

    if d[:domain] == "od_matrices"
        if d[:row] != "false"
            rows = parse.(Int,split(d[:row],','))
            return sum([mat_data(scenario, year, variable)[row,:] for row in rows]) |> json
        end
        if d[:comparewith] == "none"
            reshape(sum(mat_data(scenario, year, variable), dims = 2), :) |> json
        else
            mat_comp(scenario, year, variable,
                     comparewith, compareyear, percent=percent) |> json
        end
    elseif d[:domain] == "links"
        if d[:comparewith] == "none"
            link_data(scenario, year, variable) |> json
        else
            link_comp(scenario, year, variable,
                      comparewith, compareyear, percent=percent) |> json
        end
    else
        throw(DomainError(d[:domain]))
    end
end

route("/centroids") do
    zone_centroids |> json
end

# Want to provide:
# Zone name
# Absolute change
# Relative change
#= route("/data/popup/:zone") do =#
#=     d = getpayload() =#
#=     json(d[:zone]) =#
#= end =#

route("/oembed") do
    defaults = Dict(
        :maxwidth => 100,
        :maxheight => 100,
        :format => "json"
       )
    params = merge(defaults, getpayload())
    @assert haskey(params, :url)

    json(Dict(
        :type => "rich",
        :version => 1.0,
        # TODO: manipulate the src line to add parameters if required.
        # HTTP.URI has some code for this.
        # TODO: sanitise arguments?
        # Should check that domain and port are right
        # Shouldn't allow insertion of arbitrary HTTP
        :html => "",
        :width => params[:maxwidth],
        :height => params[:maxheight],
        :provider_name => "Integrated Transport Planning",
        :provider_url => "https://www.itpworld.net/",
       ))
end

# Adapted from VegaLite.jl
function vegalite_to_html(vl;title="Greener Connectivity Plot",width=200,height=200)
    spec = VegaLite.convert_vl_to_vg(vl)
    """

          <title>$title</title>
          <meta charset="UTF-8">

    """
end

# TODO:
# display units
# sensible ticks for years
#
# Usage:
# E.g. http://localhost:2016/api/charts?scenarios=Rail,Fleet,DoNothing&width=300&height=300
# Probably embed in an iframe
route("/charts") do
    defaults = Dict(
        :domain => "od_matrices",
        :scenarios => "Rail,DoMin",
        :variable => "Total_GHG",
        :rows => "all", # Unused
        :width => "200",
        :height => "200",
    )
    d = merge(defaults, getpayload())
    d[:rows] = d[:rows] == "all" ? Colon() : parse.(Int, split(d[:rows], ","))

    scenarios = split(d[:scenarios], ",")
    # All years that any selected scenario has data for
    years = sort(unique(vcat((metadata["scenarios"][s]["at"] for s in scenarios)...)))

    datafnc = d[:domain] == "links" ? link_data : mat_data
    df = DataFrame(year=String[], val=[], scenario=String[])
    for year in years
        for scenario in scenarios
            value = sum(datafnc(scenario, year, d[:variable])[d[:rows], :])
            scenario_name = get(metadata["scenarios"][scenario], "name", scenario)
            push!(df, (year = string(year), val = value, scenario = scenario_name))
        end
    end

    width = parse(Int, d[:width])
    height = parse(Int, d[:height])
    unit = get(metadata[d[:domain]]["columns"][d[:variable]], "unit", d[:variable])
    vl = df |> VegaLite.@vlplot(
        # Awful heuristic - these control size of plot excluding legend, labels etc
        width = width * .8,
        height = height * .5,

        mark = {
            :line,
            point = { filled = false,fill = :white },
        },
        color = {
            :scenario,
            legend = {
                title = nothing,
                orient = width < 500 ? :bottom : :right,
            },
        },
        x = {
            :year,
            title = "Year",
            type = "temporal"
        },
        y = {
            :val,
            title = unit,
            type = "quantitative",
            axis = {
                formatType = "number",
                format = ".3~s",
            },
        },
    )
    vegalite_to_html(vl; width=width, height=height)
end

Genie.AppServer.startup(parse(Int,get(ENV,"GENIE_PORT", "8000")),"0.0.0.0", async = false)

end

max_stars_repo_name: PolyCement/gcvt

max_issues_repo_name: PolyCement/gcvt

max_forks_repo_name: PolyCement/gcvt


text:
using Logging
using Glob
using ImageCore
using ImageMagick
using Test
using PNGFiles
import PNGFiles._inspect_png_read

const DEBUG_FILE = joinpath(@__DIR__, "_debugfile.jmd")

logger = ConsoleLogger(stdout, Logging.Info)
global_logger(logger)

PNG_TEST_PATH = joinpath(@__DIR__, "temp")
isdir(PNG_TEST_PATH) && rm(PNG_TEST_PATH, recursive = true)
mkdir(PNG_TEST_PATH)

include("utils.jl")

open(DEBUG_FILE, "w") do f
    write(f, "```julia\n")
    write(f, """include(\"$(joinpath(@__DIR__, "utils.jl"))\")\n """)
    write(f, "```\n")
end
function _add_debugging_entry(fpath, case, imdiff_val=missing)
    apath = joinpath(@__DIR__, abspath(fpath))
    open(DEBUG_FILE, "a") do f
        write(f, "\n\n")
        write(f, "# $(case) ############################################################\n")
        write(f, "```julia\n")
        write(f, "p = _standardize_grayness(load(\"$(apath)\"))\n")
        write(f, "i = _standardize_grayness(ImageMagick.load(\"$(apath)\"))\n")
        write(f, "plotdiffs(p, i)\n")
        write(f, "```\n")
        write(f, "```julia\n")
        write(f, "imdiff(i, p) # $(imdiff_val)\n")
        write(f, "_inspect_png_read(\"$(apath)\");\n")
        write(f, "```\n")
    end
end

@testset "PNGFiles" begin
    include("test_invalid_inputs.jl")
    include("test_pngsuite.jl")
    include("test_synthetic_images.jl")
    include("test_testimages.jl")
end

# Cleanup
isdir(PNG_TEST_PATH) && rm(PNG_TEST_PATH, recursive = true)


text:
# only should agree for pure MassAction or pure ConstantRateJumps
# mixed systems will be reordered when generating jump problems!

# pure mass action
dnanetwork = @reaction_network begin
    c1, G --> G + M
    c2, M --> M + P
    c3, M --> 0
    c4, P --> 0
    c5, 2P --> P2
    c6, P2 --> 2P
    c7, P2 + G --> P2G
    c8, P2G --> P2 + G
end c1 c2 c3 c4 c5 c6 c7 c8
rnpar = [.09, .05, .001, .0009, .00001, .0005, .005, .9]
varlabels = ["G", "M", "P", "P2","P2G"]
u0 = [1000, 0, 0, 0,0]
tf = 4000.
prob = DiscreteProblem(dnanetwork, u0, (0.0, tf), rnpar)
jprob = JumpProblem(prob, RSSA(), dnanetwork)
@test all(rxtospecies_depgraph(dnanetwork) .== jprob.discrete_jump_aggregation.jumptovars_map)
@test all(speciestorx_depgraph(dnanetwork) .== jprob.discrete_jump_aggregation.vartojumps_map)

jprob = JumpProblem(prob, NRM(), dnanetwork)
@test all(DiffEqBiological.rxtorx_depgraph(dnanetwork) .== jprob.discrete_jump_aggregation.dep_gr)

# pure ConstantRateJump
hillnetwork = @reaction_network begin
    hillr(m₃,α,K,n), ∅ --> m₁
    hillr(m₁,α,K,n), ∅ --> m₂
    hill(m₂,α,K,n), ∅ --> m₃
end α K η
p = (1.,1.,1.)
u0 = [10,10,10]
tf = 10.
prob = DiscreteProblem(hillnetwork, u0, (0.0, tf), rnpar)
jprob = JumpProblem(prob, RSSA(), hillnetwork)
@test all(rxtospecies_depgraph(hillnetwork) .== jprob.discrete_jump_aggregation.jumptovars_map)
@test all(speciestorx_depgraph(hillnetwork) .== jprob.discrete_jump_aggregation.vartojumps_map)

jprob = JumpProblem(prob, NRM(), hillnetwork)
@test all(DiffEqBiological.rxtorx_depgraph(hillnetwork) .== jprob.discrete_jump_aggregation.dep_gr)


text:
#created in including file
#comm = MPIComm{UInt64, UInt16, UInt32}()

pid = myPid(comm)

macro MPIMapTests()
    @test isa(map, BlockMap{UInt64, UInt16, UInt32})

    @test uniqueGIDs(map)

    for i = 1:5
        @test myLID(map, i)
        @test pid*5 + i - 5 == gid(map, i)
    end

    for i = 1:20
        if cld(i, 5) == pid
            @test myGID(map, i)
            @test (i-1)%5+1 == lid(map, i)
        else
            @test !myGID(map, i)
            @test 0 == lid(map, i)
        end
    end

    @test !myLID(map, -1)
    @test !myLID(map, 0)
    @test !myLID(map, 6)
    @test !myLID(map, 46)
    @test !myGID(map, -1)
    @test !myGID(map, 0)
    @test !myGID(map, 21)
    @test !myGID(map, 46)

    @test 0 == lid(map, -1)
    @test 0 == lid(map, 0)
    @test 0 == lid(map, 21)
    @test 0 == lid(map, 46)
    @test 0 == gid(map, -1)
    @test 0 == gid(map, 0)
    @test 0 == gid(map, 6)
    @test 0 == gid(map, 46)

    @test distributedGlobal(map)

    @test linearMap(map)

    @test 20 == numGlobalElements(map)
    @test 5 == numMyElements(map)

    @test pid*5 - 4 == minMyGID(map)
    @test pid*5 == maxMyGID(map)
    @test 1 == minAllGID(map)
    @test 20 == maxAllGID(map)
    @test 1 == minLID(map)
    @test 5 == maxLID(map)

    @test ([1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4],
            [1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5]) == remoteIDList(map, collect(1:20))

    @test collect((1:5) .+ 5*(pid - 1)) == myGlobalElementIDs(map)
end

map = BlockMap(20, comm)
@MPIMapTests

map = BlockMap(20, 5, comm)
@MPIMapTests

map = BlockMap(5*numProc(comm), collect((1:5) .+ 5*(pid - 1)), comm)
@MPIMapTests


text:
include("common.jl")
include("test_brngs.jl")
include("test_continuous_distributions.jl")
include("test_discrete_distributions.jl")


text:
# This file is part of Kpax3. License is MIT.

function kpax3mcmc!(data::Matrix{UInt8},
                    priorR::PriorRowPartition,
                    priorC::PriorColPartition,
                    settings::KSettings,
                    support::MCMCSupport,
                    state::AminoAcidState)
  fpS = open(string(settings.ofile, "_settings.bin"), "w")
  fpR = open(string(settings.ofile, "_row_partition.bin"), "w")
  fpC = open(string(settings.ofile, "_col_partition.bin"), "w")

  # total number of states that will be saved
  N = fld(settings.T, settings.tstep)

  # effective number of simulations
  T = settings.tstep * N

  # indices of units i and j
  ij = zeros(Int, 2)

  # neighbour indices
  neighbours = zeros(Int, support.n)

  try
    write(fpS, support.n)
    write(fpS, support.m)
    write(fpS, N)
    write(fpS, settings.α)
    write(fpS, settings.θ)
    write(fpS, settings.γ[1])
    write(fpS, settings.γ[2])
    write(fpS, settings.γ[3])
    write(fpS, settings.r)
  finally
    close(fpS)
  end

  try
    write(fpR, support.n)
    write(fpR, support.m)
    write(fpR, N)

    write(fpC, support.n)
    write(fpC, support.m)
    write(fpC, N)

    if settings.burnin > 0
      if settings.verbose
        println("Starting burnin phase...")
      end

      # sample which operators we are going to use
      operator = StatsBase.sample(UInt8[1; 2; 3], settings.op, settings.burnin)

      for t in 1:settings.burnin
        if operator[t] == 0x01
          splitmerge!(ij, neighbours, data, priorR, priorC, settings, support,
                      state)
        elseif operator[t] == 0x02
          gibbs!(data, priorR, priorC, settings, support, state)
        else
          biased_random_walk!(data, priorR, priorC, settings, support, state)
        end

        sampleC!(priorC, state)

        if settings.verbose && (t % settings.verbosestep == 0)
          println("Burnin: step ", t, " done.")
        end
      end

      if settings.verbose
        println("Burnin phase completed.")
      end
    end

    if settings.verbose
      println("Starting collecting samples...")
    end

    operator = StatsBase.sample(UInt8[1; 2; 3], settings.op, T)

    for t in 1:T
      if operator[t] == 0x01
        splitmerge!(ij, neighbours, data, priorR, priorC, settings, support,
                    state)
      elseif operator[t] == 0x02
        gibbs!(data, priorR, priorC, settings, support, state)
      else
        biased_random_walk!(data, priorR, priorC, settings, support, state)
      end

      sampleC!(priorC, state)

      if t % settings.tstep == 0
        savestate!(fpR, fpC, state)
      end

      if settings.verbose && (t % settings.verbosestep == 0)
        flush(fpR)
        flush(fpC)
        println("Step ", t, " done.")
      end
    end

    if settings.verbose
      println("Markov Chain simulation complete.")
    end
  finally
    close(fpR)
    close(fpC)
  end

  nothing
end

function kpax3mcmc(settings::KSettings)
  if settings.verbose
    Printf.@printf("Computing pairwise distances... ")
  end

  tmp = zeros(UInt8, length(settings.miss))
  idx = 0
  for c in 1:length(settings.miss)
    if settings.miss[c] != UInt8('-')
      idx += 1
      tmp[idx] = settings.miss[c]
    end
  end

  miss = if idx > 0
           copyto!(zeros(UInt8, idx), 1, tmp, 1, idx)
         else
           zeros(UInt8, 1)
         end

  (data, id, ref) = readfasta(settings.ifile, settings.protein, miss,
                              settings.l, false, 0)

  n = size(data, 2)

  d = if settings.protein
        distaamtn84(data, ref)
      else
        distntmtn93(data, ref)
      end

  D = zeros(Float64, n, n)
  idx = 1
  for j in 1:(n - 1), i in (j + 1):n
    D[i, j] = D[j, i] = d[idx]
    idx += 1
  end

  if settings.verbose
    Printf.@printf("done.\n")
  end

  # expected number of cluster approximately between cbrt(n) and sqrt(n)
  g = ceil(Int, n^(2 / 5))
  kset = max(1, g - 20):min(n, g + 20)

  x = AminoAcidData(settings)
  (m, n) = size(x.data)

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = initializestate(x.data, D, kset, priorR, priorC, settings)

  support = MCMCSupport(state, priorC)

  kpax3mcmc!(x.data, priorR, priorC, settings, support, state)

  if settings.verbose
    Printf.@printf("Processing Markov Chain output... ")
  end

  processchain(x, settings.ofile)

  if settings.verbose
    Printf.@printf("done.\n")
  end

  nothing
end

function kpax3mcmc(x::AminoAcidData,
                   partition::Vector{Int},
                   settings::KSettings)
  (m, n) = size(x.data)

  R = normalizepartition(partition, n)

  v = zeros(Int, n)
  g = 0
  k = 0
  u = 0
  for a in 1:n
    g = R[a]

    if v[g] == 0
      k += 1
    end

    v[g] += 1

    if v[g] > u
      u = v[g]
    end
  end

  maxclust = min(n, max(k, settings.maxclust))
  maxunit = min(n, max(u, settings.maxunit))

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = AminoAcidState(x.data, R, priorR, priorC, settings)

  support = MCMCSupport(state, priorC)

  kpax3mcmc!(x.data, priorR, priorC, settings, support, state)

  if settings.verbose
    Printf.@printf("Processing Markov Chain output... ")
  end

  processchain(x, settings.ofile)

  if settings.verbose
    Printf.@printf("done.\n")
  end

  nothing
end

function kpax3mcmc(x::AminoAcidData,
                   partition::String,
                   settings::KSettings)
  (m, n) = size(x.data)

  R = normalizepartition(partition, x.id)

  v = zeros(Int, n)
  g = 0
  k = 0
  u = 0
  for a in 1:n
    g = R[a]

    if v[g] == 0
      k += 1
    end

    v[g] += 1

    if v[g] > u
      u = v[g]
    end
  end

  maxclust = min(n, max(k, settings.maxclust))
  maxunit = min(n, max(u, settings.maxunit))

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = AminoAcidState(x.data, R, priorR, priorC, settings)

  support = MCMCSupport(state, priorC)

  kpax3mcmc!(x.data, priorR, priorC, settings, support, state)

  if settings.verbose
    Printf.@printf("Processing Markov Chain output... ")
  end

  processchain(x, settings.ofile)

  if settings.verbose
    Printf.@printf("done.\n")
  end

  nothing
end

function splitmerge!(ij::Vector{Int},
                     neighbours::Vector{Int},
                     data::Matrix{UInt8},
                     priorR::PriorRowPartition,
                     priorC::PriorColPartition,
                     settings::KSettings,
                     support::MCMCSupport,
                     state::AminoAcidState)
  # cluster founders (units i and j)
  StatsBase.sample!(1:support.n, ij, replace=false, ordered=false)

  # clusters of i and j respectively
  gi = state.R[ij[1]]
  gj = state.R[ij[2]]

  # total number of neighbours
  S = 0

  # generic unit
  u = 0

  if gi == gj
    for l in 1:state.v[gi]
      u = state.unit[gi][l]
      if (u != ij[1]) && (u != ij[2])
        S += 1
        neighbours[S] = u
      end
    end

    shuffleunits!(neighbours, S)

    split!(ij, neighbours, S, data, priorR, priorC, settings, support, state)
  else
    for l in 1:state.v[gi]
      u = state.unit[gi][l]
      if u != ij[1]
        S += 1
        neighbours[S] = u
      end
    end

    for l in 1:state.v[gj]
      u = state.unit[gj][l]
      if u != ij[2]
        S += 1
        neighbours[S] = u
      end
    end

    shuffleunits!(neighbours, S)

    merge!(ij, neighbours, S, data, priorR, priorC, settings, support, state)
  end

  nothing
end


text:
###
### Mer specific specializations of src/biosequence/transformations.jl
###

"""
    complement(x::T) where {T <: Skipmer}

Return the complement of a short sequence type `x`.
"""
function BioSymbols.complement(x::T) where {T<:AbstractMer}
    return T(complement_bitpar(encoded_data(x), Alphabet(x)))
end

"""
    reverse(x::T) where {T <: Skipmer}

Return the reverse of short sequence type variable `x`.
"""
function Base.reverse(x::T) where {T<:AbstractMer}
    bits = encoded_data(x)
    rbits = reversebits(bits, BitsPerSymbol{2}())
    return T(rbits >> (sizeof(bits) * 8 - 2 * length(x)))
end

"""
    reverse_complement(x::Skipmer)

Return the reverse complement of `x`.
"""
reverse_complement(x::AbstractMer) = complement(reverse(x))

"""
    canonical(kmer::Skipmer)

Return the canonical sequence of `x`.

A canonical sequence is the numerical lesser of a k-mer and its reverse complement.
This is useful in hashing/counting sequences in data that is not strand specific,
and thus observing the short sequence is equivalent to observing its reverse complement.
"""
@inline canonical(x::AbstractMer) = min(x, reverse_complement(x))

function swap(x::T, i, j) where {T<:AbstractMer}
    i = 2 * length(x) - 2i
    j = 2 * length(x) - 2j
    b = encoded_data(x)
    x = ((b >> i) ⊻ (b >> j)) & encoded_data_type(x)(0x03)
    return T(b ⊻ ((x << i) | (x << j)))
end

function Random.shuffle(x::T) where {T<:AbstractMer}
    # Fisher-Yates shuffle for mers.
    j = lastindex(x)
    for i in firstindex(x):(j - 1)
        j′ = rand(i:j)
        x = swap(x, i, j′)
    end
    return x
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestAttributes

using Test
using MathOptInterface
const MOI = MathOptInterface

include("dummy.jl")

function test_attributes_is_set_by_optimize()
    @test MOI.is_set_by_optimize(MOI.TerminationStatus())
    @test !MOI.is_set_by_optimize(MOI.ConstraintSet())
    @test !MOI.is_set_by_optimize(MOI.ObjectiveSense())
    @test MOI.is_set_by_optimize(MOI.CallbackNodeStatus(1))
end

function test_attributes_is_copyable()
    @test !MOI.is_copyable(MOI.TerminationStatus())
    @test !MOI.is_copyable(MOI.ConstraintSet())
    @test MOI.is_copyable(MOI.ObjectiveSense())
end

function test_attributes_supports()
    model = DummyModel()
    @test_throws ArgumentError MOI.supports(model, MOI.TerminationStatus())
    @test_throws ArgumentError begin
        MOI.supports(
            model,
            MOI.ConstraintSet(),
            MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
        )
    end
    @test MOI.supports(model, MOI.ObjectiveSense())
end

function test_attributes_set_vector()
    attr = MOI.VariablePrimalStart()
    err = DimensionMismatch(
        "Number of indices (1) does not match the " *
        "number of values (2) set to `$attr`.",
    )
    model = DummyModel()
    x = MOI.VariableIndex(1)
    @test_throws err MOI.set(model, MOI.VariablePrimalStart(), [x], ones(2))
end

function test_attributes_integration_compute_conflict_1()
    optimizer = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.Model{Float64}(),
        MOI.Bridges.full_bridge_optimizer(optimizer, Float64),
    )
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, x, MOI.LessThan(0.0))
    c2 = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.optimize!(model)
    @test MOI.get(optimizer, MOI.ConflictStatus()) ==
          MOI.COMPUTE_CONFLICT_NOT_CALLED
    MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.set(
        optimizer,
        MOI.ConstraintConflictStatus(),
        MOI.get(
            optimizer,
            MOI.ListOfConstraintIndices{
                MOI.VariableIndex,
                MOI.LessThan{Float64},
            }(),
        )[1],
        MOI.NOT_IN_CONFLICT,
    )
    MOI.set(
        optimizer,
        MOI.ConstraintConflictStatus(),
        MOI.get(
            optimizer,
            MOI.ListOfConstraintIndices{
                MOI.VariableIndex,
                MOI.GreaterThan{Float64},
            }(),
        )[1],
        MOI.IN_CONFLICT,
    )
    MOI.compute_conflict!(model)
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) ==
          MOI.NOT_IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.IN_CONFLICT
end

MOI.Utilities.@model(
    OnlyScalarConstraints,
    (),
    (MOI.GreaterThan, MOI.LessThan),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)

function test_attributes_integration_compute_conflict_2()
    optimizer = MOI.Utilities.MockOptimizer(OnlyScalarConstraints{Float64}())
    model = MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.Interval(0.0, 1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(optimizer, MOI.ConflictStatus()) ==
          MOI.COMPUTE_CONFLICT_NOT_CALLED
    MOI.set(optimizer, MOI.ConflictStatus(), MOI.CONFLICT_FOUND)
    MOI.compute_conflict!(model)
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test_throws ArgumentError MOI.get(model, MOI.ConstraintConflictStatus(), c)
end

struct _NoConstraintName <: MOI.AbstractOptimizer end

function test_no_constraint_name()
    model = _NoConstraintName()
    @test_throws(
        MOI.VariableIndexConstraintNameError(),
        MOI.supports(
            model,
            MOI.ConstraintName(),
            MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne},
        ),
    )
    @test_throws(
        MOI.VariableIndexConstraintNameError(),
        MOI.set(
            model,
            MOI.ConstraintName(),
            MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(1),
            "name",
        ),
    )
end

function test_get_fallback()
    model = DummyModelWithAdd()
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.EqualTo(0.0))
    err = ArgumentError(
        "$(typeof(model)) does not support getting the attribute " *
        "$(MOI.SolveTimeSec()).",
    )
    @test_throws(err, MOI.get(model, MOI.SolveTimeSec()))
    output = Ref{Cdouble}()
    @test_throws(err, MOI.get!(output, model, MOI.SolveTimeSec()))
    errv = ArgumentError(
        "$(typeof(model)) does not support getting the attribute " *
        "$(MOI.VariablePrimal()).",
    )
    @test_throws(errv, MOI.get(model, MOI.VariablePrimal(), x))
    errc = ArgumentError(
        "$(typeof(model)) does not support getting the attribute " *
        "$(MOI.ConstraintPrimal()).",
    )
    @test_throws(errc, MOI.get(model, MOI.ConstraintPrimal(), c))
    @test_throws(
        ArgumentError(
            "Unable to get attribute $(MOI.VariablePrimal()): invalid " *
            "arguments $((c,)).",
        ),
        MOI.get(model, MOI.VariablePrimal(), c),
    )
end

function test_ConstraintBasisStatus_fallback()
    model = DummyModelWithAdd()
    c = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    @test_throws(
        ErrorException(
            "Querying the basis status of a `VariableIndex` constraint is " *
            "not supported. Use [`VariableBasisStatus`](@ref) instead.",
        ),
        MOI.get(model, MOI.ConstraintBasisStatus(), c),
    )
end

function test_UnsupportedSubmittable()
    model = DummyModelWithAdd()
    sub = MOI.LazyConstraint{Int}(1)
    @test_throws(
        MOI.UnsupportedSubmittable(
            sub,
            "submit(::$(typeof(model)), ::$(typeof(sub))) is not supported.",
        ),
        MOI.submit(model, sub, 1),
    )
end

function test_attribute_value_type()
    @test MOI.attribute_value_type(MOI.CallbackNodeStatus(1)) ==
          MOI.CallbackNodeStatusCode
    @test MOI.attribute_value_type(MOI.LazyConstraintCallback()) == Function
    @test MOI.attribute_value_type(MOI.RelativeGap()) == Float64
    @test MOI.attribute_value_type(MOI.SimplexIterations()) == Int64
    @test MOI.attribute_value_type(MOI.BarrierIterations()) == Int64
    @test MOI.attribute_value_type(MOI.NodeCount()) == Int64
    @test MOI.attribute_value_type(
        MOI.ConstraintBridgingCost{MOI.VariableIndex,MOI.ZeroOne}(),
    ) == Float64
    @test MOI.attribute_value_type(MOI.VariableBridgingCost{MOI.ZeroOne}()) ==
          Float64
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestAttributes.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
"""
estimate_condeig(V)

Estimate the condition number of an upper triangular matrix V; the matrix
V is expected to be either a `Matrix{ComplexF64}` or `Matrix{Complex{BigFloat}}`.

The condition number is first estimated converting the matrix in standard floating
point; if the estimate is not accurate enough, it is replaced with the condition number
of |V|; if the latter is too loose, it is computed in higher precision.
"""
function estimate_condeig(V::UpperTriangular)

    # Compute a floating point version of V
    VF = convert(Matrix{ComplexF64}, V)
    p = 1

    nF = opnorm(VF, p)
    kV = nF * opnorm(inv(VF), p)

    if kV <= 1e-2 / eps()
        return kV
    end

    VV = -abs.(VF)
    VV = VV - diagm(2 * diag(VV))
    kV2 = opnorm(inv(VV), p) * nF

    if kV2 <= kV * 1.0e4
        return kV2
    end

    kV3 = opnorm(VF, p) * opnorm(convert(Matrix{ComplexF64}, inv(V)), p)

    return kV3
end

"""
build_eigendecomp_tree!(T::BivMatTree, maxkV::Float64, A::Union{Matrix{ComplexF64}, Matrix{Float64}}, meth::Algorithm, use_mp = true)

This recursive routine is used to construct the eigenvector matrices related to the atomic blocks of the
matrix A (i.e., the diagonal blocks corresponding to the selected partitioning).
"""
function build_eigendecomp_tree!(T::BivMatTree, maxkV::Float64, A::Union{Matrix{ComplexF64}, Matrix{Float64}}, meth::Algorithm, use_mp = true)
    if isa(T.A11, BivMatTreeTail)

        uh = eps() / T.kV / maxkV;
        d_uh = convert(Int64, max(32, ceil(log10(1 / uh))));

        if use_mp
            ok = false
            while !ok
                with_digits(d_uh) do
                    ok = true
                    T1 = mp(T.Ttilde, d_uh)

                    # FIXME: Need to handle the 2x2 blocks in the real Schur form
                    (T.D, T.V) = eigentri(T1);
                    T.kV = estimate_condeig(T.V);

                    uh = eps() / T.kV / maxkV;
                    d_uh_new = convert(Int64, max(32, ceil(log10(1 / uh))));

                    if d_uh_new > d_uh
                        d_uh = d_uh_new
                        ok = false
                    end
                end
            end
        else
            (T.D, T.V) = eigentri(T.Ttilde);
            T.kV = cond(T.V)
        end
    else
        build_eigendecomp_tree!(T.A11, maxkV, A, meth, use_mp);
        build_eigendecomp_tree!(T.A22, maxkV, A, meth, use_mp);
    end
end

"""
build_partition_tree!(T, A, ind, strategy, method, delta = 0.1, bs = 1, mergeblocks = true)

Recursive function constructing the partition tree for a matrix A, according to the given parameters for
    separating the eigenvalues.
    """
    function build_partition_tree!(T, A, ind, strategy, method, delta = 0.1, bs = 1, mergeblocks = true)

        l = length(ind);

        maxkV = 1.0;

        if l == 1
            T.ind = ind[1];
            nblocks = 1;

            if method == Diag
                (T.Ttilde, T.kV) = condeigvec(A[T.ind, T.ind]);
                maxkV = T.kV;
            else
                T.Ttilde = A[T.ind, T.ind];
            end

            T.A11 = BivMatTreeTail();
            T.A22 = BivMatTreeTail();

            return nblocks, maxkV;
        elseif sum(map(length, ind)) <= bs
            T.ind = reduce(vcat, ind);

            nblocks = 1;

            if method == Diag
                (T.Ttilde, T.kV) = condeigvec(A[T.ind, T.ind]);
                maxkV = T.kV;
            else
                T.Ttilde = A[T.ind, T.ind];
            end

            T.A11 = BivMatTreeTail();
            T.A22 = BivMatTreeTail();

            return nblocks, maxkV;
        else
            T.ind = reduce(vcat, ind);

            if strategy == "balanced"
                v = cumsum(map(length, ind));
                _, ii = findmin(map(abs, v .- v[end]/2));
                ind1 = ind[1:ii];
                ind2 = ind[ii + 1: end];
            elseif strategy == "single"
                v = cumsum(map(length, ind[end:-1:1]));
                i = max(1, sum(v .< bs));
                ind1 = ind[1:end - i];
                ind2 = ind[end-i+1:end];
            else
                error("Unsupported strategy argument")
            end

            i1 = ind1[1][1]; i2 = ind1[end][end];
            j1 = ind2[1][1]; j2 = ind2[end][end];

            T.A11 = newBivMatTree();
            T.A22 = newBivMatTree();

            nblocks1, kV1 = build_partition_tree!(T.A11, A, ind1, strategy, method, delta, bs, mergeblocks);
            nblocks2, kV2 = build_partition_tree!(T.A22, A, ind2, strategy, method, delta, bs, mergeblocks);

            maxkV = max(kV1, kV2);

            nblocks = nblocks1 + nblocks2;

            T.sylv = sylvester(A[i1:i2, i1:i2], -A[j1:j2, j1:j2], -A[i1:i2, j1:j2])

            ninv = norm(T.sylv) / norm(A[i1:i2, j1:j2]);

            # If the Sylvester equation is
            # too ill-conditioned we merge the corresponding blocks
            if ninv > 10 / delta
                if mergeblocks
                    T.ind = reduce(vcat, ind);
                    T.sylv = nothing;
                    T.A11 = BivMatTreeTail();
                    T.A22 = BivMatTreeTail();
                    nblocks = 1;
                    if method == Diag
                        T.Ttilde, T.kV = condeigvec(A[T.ind, T.ind]);
                        maxkV = T.kV;
                    else
                        T.Ttilde = A[T.ind, T.ind]
                    end
                else
                    warn("Ill-conditioned Sylvester equation: results may be inaccurate");
                end
            end

            return nblocks, maxkV
        end
    end

    function fun2m_preprocessing(A, deltaA, method, bs, strategy)

        if istriu(A)
            TA = A; UA = one(A);
            diagTA = diag(TA);
            FA = Schur(TA, UA, diag(TA))
        else
            FA = schur(A);
            TA = FA.T; UA = FA.Z;
            diagTA = diag(TA);
        end

        # Determine reordering of the Schur forms into block forms.
        ordA = blocking(TA, deltaA);
        ordA, indA = swapping(ordA);  # Gives the blocking.
        ordA::Array{LinearAlgebra.BlasInt,1} = maximum(ordA) .- ordA .+ 1;    # Since ORDSCHUR puts highest index top left.
        ordschur!(FA, ordA)
        TA = FA.T; UA = FA.Z;

        mergeblocks = method != Taylor;

        treeA = newBivMatTree();

        nblocksA, maxkVA = build_partition_tree!(treeA, TA, indA, strategy, method, deltaA, bs, mergeblocks);

        if method == Diag || method == DiagNoHp
            use_mp = method == Diag
            build_eigendecomp_tree!(treeA, maxkVA, A, method, use_mp);
        end

        return UA, TA, treeA, nblocksA
    end

    """
    fun2m(f, A, B, C; <keyword arguments>)

    Evaluate the bivariate matrix functions X = f{A, B^T}(C); the function f should have
    the prototype:

    f = (x,y,i,j) -> f^(i,j)(x,y),

    which denotes the i-th derivative of f with respect to x, the j-th derivative with respect to y, evaluate at a point (x, y).

    Note that, unless `BivMatFun.Taylor` is used (see below), this function is always called with `i = j = 0`.

    # Arguments
    Several optional arguments can be specified to control the behavior of `fun2m`:
    - `delta::Float64 = 0.05`: Minimum separation between the eigenvalues, used to determined the partitioning.
    - `method::Algorithm = BivMatFun.Diag`: Algorithm used to evaluate the function. The default is `BivMatFun.Diag`, which implements a mixed precision perturb and diagonalize approach. Other options are `BivMatFun.Taylor` (using Taylor expansions at the atomic blocks), or `BivMatFun.DiagNoHp`, which uses perturb-and-diagonalize avoiding multiprecision. The last option is provided for benchmarks only, and is not expected to be accurate. The only option that needs to evaluate the derivatives is `BivMatFun.Taylor`. For all other cases the function `f` is always called with `i = j = 0`.
    - `bs::Union{Nothing, Int64} = nothing`: Minimum block-size to use in the partitioning. Specify `nothing` to automatically choose an optimal value based on the chosen `method`.
    - `min_digits::Int64 = 0`: Minimum number of digits to use in multiprecision computations. Use `0` (default) to not pose any constraint.
    - `strategy::String = "balanced"`: strategy to use for splitting the blocks. Valid options are "balanced" or "single".
    - `user_function = nothing`: if different from nothing, this function is used to evaluate the atomic blocks. It needs to have the prototype (A, B, C) -> X.
    - `parallel::Bool = true`: Use the multithreaded implementation.
    """
    function fun2m(f, A, B, C; delta::Float64 = 0.05, method::Algorithm = Diag, bs::Union{Nothing, Int64} = nothing,
        min_digits::Int64 = 0, strategy::String = "balanced", user_function = nothing, parallel::Bool = true)

        info = Fun2MInfo(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing);

        # Default parameters.
        if bs === nothing
            if method == Taylor
                bs = 1;
            elseif method == Diag
                bs = 4;
            elseif method == DiagNoHp
                bs = 4;
            else
                bs = 1;
            end
        end

        # Checking input dimensions
        m = size(A, 1);
        n = size(B, 1);

        if m != size(A, 2)
            error("A must be square")
        end

        if n != size(B, 2)
            error("B must be square")
        end

        if m != size(C, 1) || n != size(C, 2)
            error("Incompatible dimensions");
        end

        # When at least one among A and B is a scalar
        if m == 1 || n == 1
            F = scalar_fun(f, A, B, C);
            return;
        end

        if parallel
            t1 = Threads.@spawn fun2m_preprocessing(A, delta, method, bs, strategy)
            t2 = Threads.@spawn fun2m_preprocessing(B, delta, method, bs, strategy)

            UA, TA, treeA, info.nblocksA = fetch(t1)
            UB, TB, treeB, info.nblocksB = fetch(t2)
        else
            UA, TA, treeA, info.nblocksA = fun2m_preprocessing(A, delta, method, bs, strategy)
            UB, TB, treeB, info.nblocksB = fun2m_preprocessing(B, delta, method, bs, strategy)
        end

        info.time_fun2mric = @timed begin

            C = UA' * C * UB;
            F = similar(C)
            info.digits = fun2m_ric!(f, F, TA, TB, C, treeA, treeB, method, min_digits, 0, user_function, parallel);
            F = UA * F * UB';

        end # End of time fun2mric

        if isreal(A) && isreal(B) && isreal(C) && norm(imag(F),1) <= 10*n*eps()*norm(F,1)
            F = real(F);
        end

        return F, info;
    end

    # ----------------------------------Recursive part---------------------------------------------------
    function fun2m_ric!(f, F, A, B, C, treeA::BivMatTree, treeB::BivMatTree, method::Algorithm, min_digits::Int64, level::Int64, user_function, parallel::Bool)
        # Computes the bivariate matrix functions on triangularized and blocked
        # coefficients via recursion and Sylvester solution
        # Base cases

        digits = 0;

        (m, n) = size(C);

        if isa(treeA, BivMatTreeTail) || isa(treeB, BivMatTreeTail)
            F .= 0;
            return digits;
        end

        iA = treeA.ind; iB = treeB.ind;

        if isa(treeA.A11, BivMatTreeTail) && isa(treeB.A11, BivMatTreeTail) # atomic block
            F[:], digits = fun2m_atom(f, A[iA, iA], B[iB, iB], C, treeA, treeB, method, min_digits, user_function)
            return digits;
        else # if at least one of the two has children

            treeA_A11 = treeA.A11; treeA_A22 = treeA.A22;
            treeB_A11 = treeB.A11; treeB_A22 = treeB.A22;

            if ! isa(treeA.A11, BivMatTreeTail)
                m1 = length(treeA.A11.ind);
            else
                m1 = length(iA);
                treeA_A11 = treeA; # copy the node into the first child
            end
            if ! isa(treeB.A11, BivMatTreeTail)
                n1 = length(treeB.A11.ind);
            else
                n1 = length(iB);
                treeB_A11 = treeB; # copy the node into the first child
            end

            # Split the problem by halving the blocks of A and B
            V1 = (treeA.sylv !== nothing) ? treeA.sylv : zeros(m1, 0)
            W1 = (treeB.sylv !== nothing) ? treeB.sylv : zeros(n1, 0)

            C11 = C[1:m1, 1:n1] + V1 * C[m1+1:end,1:n1]
            C21 = C[m1+1:end,1:n1]
            C12 = -C[1:m1,1:n1] * W1 + C[1:m1,n1+1:end] - V1 * C[m1+1:end,1:n1] * W1 + V1 * C[m1+1:end,n1+1:end]
            C22 = -C[m1+1:end,1:n1] * W1 + C[m1+1:end,n1+1:end]

            F11 = similar(C11)
            F21 = similar(C21)
            F12 = similar(C12)
            F22 = similar(C22)

            if parallel && level <= 2
                t11 = Threads.@spawn fun2m_ric!(f, F11, A, B, C11, treeA_A11, treeB_A11, method, min_digits, level+1, user_function, parallel);
                t21 = Threads.@spawn fun2m_ric!(f, F21, A, B, C21, treeA_A22, treeB_A11, method, min_digits, level+1, user_function, parallel);
                t12 = Threads.@spawn fun2m_ric!(f, F12, A, B, C12, treeA_A11, treeB_A22, method, min_digits, level+1, user_function, parallel);
                t22 = Threads.@spawn fun2m_ric!(f, F22, A, B, C22, treeA_A22, treeB_A22, method, min_digits, level+1, user_function, parallel);
                digits = max(fetch(t11), fetch(t12), fetch(t21), fetch(t22))
            else
                digits11 = fun2m_ric!(f, F11, A, B, C11, treeA_A11, treeB_A11, method, min_digits, level+1, user_function, parallel);
                digits21 = fun2m_ric!(f, F21, A, B, C21, treeA_A22, treeB_A11, method, min_digits, level+1, user_function, parallel);
                digits12 = fun2m_ric!(f, F12, A, B, C12, treeA_A11, treeB_A22, method, min_digits, level+1, user_function, parallel);
                digits22 = fun2m_ric!(f, F22, A, B, C22, treeA_A22, treeB_A22, method, min_digits, level+1, user_function, parallel);
                digits = max(digits11, digits12, digits21, digits22);
            end

            F[1:m1,1:n1] = F11 - V1 * F21;
            F[1:m1,n1+1:end] = F11 * W1 - V1 * F21 * W1 + F12 - V1 * F22;
            F[m1+1:end,1:n1] = F21;
            F[m1+1:end,n1+1:end] = F21 * W1 + F22;

            return digits;
        end

    end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun