text:
function FESpaces.FESpace(model::DistributedDiscreteModel,
                          reffe::Tuple{RaviartThomas,Any,Any};
                          conformity=nothing,kwargs...)

  cell_reffes = map_parts(model.models) do m
    basis,reffe_args,reffe_kwargs = reffe
    cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
  end
  _common_fe_space_constructor(model,cell_reffes;conformity,kwargs...)
end

function FESpace(model::DistributedDiscreteModel,
                 reffe::GenericRefFE{RaviartThomas};
                 conformity=nothing, kwargs...)
  cell_reffes = map_parts(model.models) do m
    Fill(reffe,num_cells(m))
  end
  _common_fe_space_constructor(model,cell_reffes;conformity,kwargs...)
end

function _common_fe_space_constructor(model,cell_reffes;conformity,kwargs...)
  sign_flips=_generate_sign_flips(model,cell_reffes)
  spaces = map_parts(model.models,sign_flips,cell_reffes) do m,sign_flip,cell_reffe
     conf = Conformity(testitem(cell_reffe),conformity)
     cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
     FESpace(m, cell_fe; kwargs...)
  end
  gids =  generate_gids(model,spaces)
  vector_type = _find_vector_type(spaces,gids)
  DistributedSingleFieldFESpace(spaces,gids,vector_type)
end

function _generate_sign_flips(model,cell_reffes)
  sign_flips=map_parts(model.models,model.gids.partition,cell_reffes) do m, p, cell_reffe

    D = num_cell_dims(m)
    gtopo = get_grid_topology(m)

    # Extract composition among cells and facets
    cell_wise_facets_ids = get_faces(gtopo, D, D - 1)
    cache_cell_wise_facets_ids = array_cache(cell_wise_facets_ids)

    # Extract cells around facets
    cells_around_facets = get_faces(gtopo, D - 1, D)
    cache_cells_around_facets = array_cache(cells_around_facets)

    ncells = num_cells(m)
    ptrs = Vector{Int32}(undef,ncells+1)
    for cell in 1:ncells
      reffe=cell_reffe[cell]
      ptrs[cell+1] = num_dofs(reffe)
    end
    PArrays.length_to_ptrs!(ptrs)
    ndata = ptrs[end]-1
    data  = Vector{Bool}(undef,ndata)
    data .= false

    for cell in p.oid_to_lid
      sign_flip = view(data,ptrs[cell]:ptrs[cell+1]-1)
      reffe=cell_reffe[cell]
      D = num_dims(reffe)
      face_own_dofs = get_face_own_dofs(reffe)
      facet_lid = get_offsets(get_polytope(reffe))[D] + 1
      cell_facets_ids = getindex!(cache_cell_wise_facets_ids,
                                  cell_wise_facets_ids,
                                  cell)
      for facet_gid in cell_facets_ids
          facet_cells_around = getindex!(cache_cells_around_facets,
                cells_around_facets,
                facet_gid)
          is_slave=false
          if (length(facet_cells_around)==1)
            is_slave == false
          else
            mx=maximum(p.lid_to_gid[facet_cells_around])
            is_slave = (p.lid_to_gid[cell] == mx)
          end
          if is_slave
              for dof in face_own_dofs[facet_lid]
                  sign_flip[dof] = true
              end
          end
          facet_lid = facet_lid + 1
      end
    end
    PArrays.Table(data,ptrs)
  end
  exchange!(sign_flips,model.gids.exchanger)
  sign_flips
end


text:
@testset "Solids tests" begin
    @testset "Barite" begin
        barite=Barite()

        @test barite.ρ≈rho(101325.0,298.15,barite)
        @test barite.cₚ≈Cp(101325.0,298.15,barite)
        @test barite.k≈therm_cond(101325.0,298.15,barite)
    end

    @testset "Hematite" begin
        hematite=Hematite()

        @test hematite.ρ≈rho(101325.0,298.15,hematite)
        @test hematite.cₚ≈Cp(101325.0,298.15,hematite)
        @test hematite.k≈therm_cond(101325.0,298.15,hematite)
    end

    @testset "CaCO3" begin
        carbonite=CaCO3()

        @test carbonite.ρ≈rho(101325.0,298.15,carbonite)
        @test carbonite.cₚ≈Cp(101325.0,298.15,carbonite)
        @test carbonite.k≈therm_cond(101325.0,298.15,carbonite)
    end

end

ext: jl


text:
@testset "readcards Tests" begin
    using CircuitSimulator: readnoncommentline, cardpart, readcards
    @test readnoncommentline(IOBuffer("* comment line\n not a comment")) == " not a comment"
    @test readnoncommentline(IOBuffer("* comment line\n")) == ""
    testbuffer = IOBuffer("""
    * comment line
    * also should remove blank lines

    * line that
    this is not a comment
    * but this is
    2nd non comment line
    """)
    @test readnoncommentline(testbuffer) == "this is not a comment"
    @test readnoncommentline(testbuffer) == "2nd non comment line"
    @test cardpart("   abcd \$ comment") == (false, "abcd ")
    @test cardpart("   abcd // comment") == (false, "abcd ")
    @test cardpart("  + abcd // comment") == (true, " abcd ")
    @test cardpart("  + a\"\$\"bcd // comment") == (true, " a\"\$\"bcd ")
    @test cardpart("  + a\'\$\'bcd // comment") == (true, " a\'\$\'bcd ")
    @test cardpart("   abcd μ\$μ comment") == (false, "abcd μ")
    testbuffer = IOBuffer("""
    * Title Line
    *V1 N001 0 PULSE(0 1 0 .01 .01 {.5-.01} 1 10)
    V1 N001 0 1V
    R1 N003 N001 1k
    R2 0 N003 67.9
    R3 N003 N002 234\$ comment
    + a = 2 b = 3// comment
    + c = 6
    * comment
    + d = 7
    R4 N002 0 1Meg
    C1 N002 N001 1µ// comment
    .tran 20
    .backanno\$ comment
    .end
    """)
    verified = IOBuffer("""
    V1 N001 0 1V
    R1 N003 N001 1k
    R2 0 N003 67.9
    R3 N003 N002 234 a = 2 b = 3 c = 6 d = 7
    R4 N002 0 1Meg
    C1 N002 N001 1µ
    .tran 20
    .backanno
    .end
    """)
    testbuffercollected = collect(readcards(testbuffer))
    @test (length(testbuffercollected)) == 9
    for card in testbuffercollected
        @test card == readline(verified)
    end
end;


text:
@testset "315.count-of-smaller-numbers-after-self.jl" begin
    f1 = count_smaller_method1
    f2 = count_smaller_method2
    @test f1([5, 2, 6, 1]) == f2([5, 2, 6, 1]) == [2, 1, 1, 0]
    @test f1([2, 0, 1]) == f2([2, 0, 1]) == [2, 0, 0]
    @test f1([-1]) == f2([-1]) == [0]
    @test f1([-1, -1]) == f2([-1, -1]) == [0, 0]
    @test f1([7, 6, 5, 4, 3, 2, 1]) == f2([7, 6, 5, 4, 3, 2, 1]) == [6, 5, 4, 3, 2, 1, 0]
end


text:
# An inflexible stationary demand.

# Import packages.

using Dates

# Demand without flexibility.

"""
A stationary demand that has no flexibility.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location

# Fields
- `location`: network location of the demand
- `demands` : timeseries of the demands

# Example
```
example = StationaryInflexibleDemand(
    "a fixed location",
    TimeArray(
        [Time(0), Time(8), Time(9), Time(17), Time(18)],
        [     5.,     10.,      8.,      11.,       3.]
    )
) :: StationaryInflexibleDemand{Time,String}
```
"""
struct StationaryInflexibleDemand{T,L} <: InflexibleDemand{T,L}
    location :: L
    demands  :: TemporalDemand{T}
end

"""
Time-varying demands and their locations.

# Arguments
- `demand :: StationaryInflexibleDemand{T,L}`: the demand
"""
function demands(demand :: StationaryInflexibleDemand{T,L}) :: LocatedDemand{T,L} where L where T <: TimeType
    TimeArray(
        timestamp(demand.demands),
        map(x -> (demand.location, x), values(demand.demands))
    )
end


text:
"""
Data specifying a bilevel linear-linear problem
"""
struct BilevelProblem
    cx::Vector{Float64}
    cy::Vector{Float64}
    G::Matrix{Float64}
    H::Matrix{Float64}
    q::Vector{Float64}
    d::Vector{Float64}
    A::Matrix{Float64}
    B::Matrix{Float64}
    b::Vector{Float64}
    nu::Int
    nl::Int
    mu::Int
    ml::Int
    function BilevelProblem(cx, cy, G, H, q, d, A, B, b)
        nu = length(cx)
        nl = length(cy)
        nl == length(d) || DimensionMismatch("Objectives")
        mu = length(q)
        ml = length(b)
        size(A) == (ml,nu) && size(B) == (ml,nl) || DimensionMismatch("Lower constraints")
        size(G) == (mu,nu) && size(H) == (mu,nl) || DimensionMismatch("Higher constraints")
        new(cx,cy,G,H,q,d,A,B,b,nu,nl,mu,ml)
    end
end

abstract type NearOptimalMethod end

"""
    ExtendedMethod

Solution method for the near-optimal problem using the extended formulation
"""
struct ExtendedMethod <: NearOptimalMethod end

"""
    LazyExtended

Solve the model and lazily verify which sub-problems must be added
"""
struct LazyExtended <: NearOptimalMethod end

"""
    BilinearTerms

Solution approach for the near-optimal problem using the bilinear formulation
"""
struct BilinearTerms <: NearOptimalMethod end

"""
    LazyBatched(batch_size::Int)

Solve the model with the lazy subproblem expansion and maximum batch size `batch_size`.
"""
struct LazyBatched <: NearOptimalMethod
    batch_size::Int
end


text:
using CUDA
import MPI
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using BenchmarkTools
using Printf
using ProgressBars

include("part1_utils.jl")

macro qx(ix, iy, iz)
    esc(:(-D_dx * (Hτ[$ix, $iy, $iz] - Hτ[$ix-1, $iy, $iz])))
end
macro qy(ix, iy, iz)
    esc(:(-D_dy * (Hτ[$ix, $iy, $iz] - Hτ[$ix, $iy-1, $iz])))
end
macro qz(ix, iy, iz)
    esc(:(-D_dz * (Hτ[$ix, $iy, $iz] - Hτ[$ix, $iy, $iz-1])))
end

struct BenchResults
    Δt :: Float64
    Work :: Float64
    Performance :: Float64
    Memory :: Float64
    Intensity :: Float64
    Throughput :: Float64
end

""" 3D diffusion step using kernel programming

Cells         = (nx-2) * (ny-2) * (nz-2)

Work per cell = 3 * (3 * sub + 3 * mul) + 1 * (1 * sub + 1 * mul) + 3 * add + 1 * (1 * sub + 1 * mul)
              = 25 FLOP
Work          = 25 FLOP * Cells

Memory moved per cell  = 9 * read Hτ + 1 * read Ht + 1 * read/write dHdτ + 1 * read/write Hτ2
                       = 14 Float64 read or write
Memory moved           = 14 Float64 * Cells

Note that we can probably reduce the memory reads by a lot using shared memory
"""
@parallel_indices (ix, iy, iz) function diffusion_3D_step_τ(Ht, Hτ, Hτ2, dHdτ, dτ, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz)
    if (1 < ix < size(Hτ, 1) && 1 < iy < size(Hτ, 2) && 1 < iz < size(Hτ, 3))
        dHdτ[ix, iy, iz] = (  # We read full stencil from Hτ ( 9 elements in 3D ) and write once
            (@qx(ix + 1, iy, iz) - @qx(ix, iy, iz)) * _dx +  # Work = 3 * sub + 3 * mul
            (@qy(ix, iy + 1, iz) - @qy(ix, iy, iz)) * _dy +
            (@qz(ix, iy, iz + 1) - @qz(ix, iy, iz)) * _dz +
            (Hτ[ix, iy, iz] - Ht[ix, iy, iz]) * _dt  # Work = 1 * sub + 1 * mul
        )
        Hτ2[ix, iy, iz] =
          Hτ[ix, iy, iz] - dτ * dHdτ[ix, iy, iz]  # Work = 1 * sub + 1 * mul (or 1 * fma)
    end
    return nothing
end

""" 3D diffusion step using kernel programming, with shared memory

Similar to `diffusion_3D_step_τ`, but massively reducing memory loads by using shared memory (on GPU).

Cells         = (nx-2) * (ny-2) * (nz-2)

Work per cell = 3 * (3 * sub + 3 * mul) + 1 * (1 * sub + 1 * mul) + 3 * add + 1 * (1 * sub + 1 * mul)
              = 25 FLOP
Work          = 25 FLOP * Cells

Memory moved per cell  = 1 * read Hτ + 1 * read Ht + 1 * read/write dHdτ + 1 * read/write Hτ2
                       = 6 Float64 read or write
Memory moved           = 6 Float64 * Cells

"""
@parallel_indices (ix, iy, iz) function diffusion_3D_step_τ_shared_memory(Ht, Hτ_, Hτ2, dHdτ, dτ, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz)
    tx, ty, tz = @threadIdx().x + 1, @threadIdx().y + 1, @threadIdx().z + 1
    Hτ = @sharedMem(eltype(Ht), (@blockDim().x+2, @blockDim().y+2, @blockDim().z+2))
    Hτ[tx, ty, tz] = Hτ_[ix, iy, iz]
    if (1 < ix < size(Hτ_, 1) && 1 < iy < size(Hτ_, 2) && 1 < iz < size(Hτ_, 3))
        if (@threadIdx().x == 1)             Hτ[tx-1,ty,tz] = Hτ_[ix-1,iy,iz] end
        if (@threadIdx().x == @blockDim().x) Hτ[tx+1,ty,tz] = Hτ_[ix+1,iy,iz] end
        if (@threadIdx().y == 1)             Hτ[tx,ty-1,tz] = Hτ_[ix,iy-1,iz] end
        if (@threadIdx().y == @blockDim().y) Hτ[tx,ty+1,tz] = Hτ_[ix,iy+1,iz] end
        if (@threadIdx().z == 1)             Hτ[tx,ty,tz-1] = Hτ_[ix,iy,iz-1] end
        if (@threadIdx().z == @blockDim().z) Hτ[tx,ty,tz+1] = Hτ_[ix,iy,iz+1] end
        @sync_threads()
        dHdτ[ix, iy, iz] = (  # We read full stencil from Hτ ( 9 elements in 3D ) and write once
            (@qx(tx + 1, ty, tz) - @qx(tx, ty, tz)) * _dx +  # Work = 3 * sub + 3 * mul
            (@qy(tx, ty + 1, tz) - @qy(tx, ty, tz)) * _dy +
            (@qz(tx, ty, tz + 1) - @qz(tx, ty, tz)) * _dz +
            (Hτ_[ix, iy, iz] - Ht[ix, iy, iz]) * _dt  # Work = 1 * sub + 1 * mul
        )
        Hτ2[ix, iy, iz] =
          Hτ_[ix, iy, iz] - dτ * dHdτ[ix, iy, iz]  # Work = 1 * sub + 1 * mul (or 1 * fma)
    end
    return nothing
end

@views function diffusion_3D_kernel_programming(; nx, ny, nz, ttot = 1., tol=1e-8, use_shared_memory=true, do_vis = false, verbose=true, init_and_finalize_MPI=!isinteractive(), scale_physical_size=false)
    me, dims, _nprocs, coords, comm_cart =
        init_global_grid(nx, ny, nz; init_MPI = init_and_finalize_MPI, quiet=!verbose)

    # physics
    D = 1.0

    # Allow scaling the physical problem dimensions with the number of grid points.
    # This is useful for weak scaling experiments where we want to keep the work per worker fixed.
    # In particular, note that usually increasing nx -> 2*nx => dx -> dx/2 => dτ = dτ/4  => !nonlinear!
    # Therefore, by scaling the physical size together with (nx, ny, nz) we can keep (dx, dy, dz) constant and therefore keep the work per worker constant.
    if scale_physical_size
        lx, ly, lz = dims .* (10., 10., 10.)
    else
        lx, ly, lz = 10.0, 10.0, 10.0
    end

    # derived numerics
    dx, dy, dz = lx / nx_g(), ly / ny_g(), lz / nz_g()

    # bind MPI ranks to GPUs
    if USE_GPU
        select_device()
    end

    total_N = prod(dims) * nx * ny * nz

    # numerics
    dt = 0.2
    dτ = min(dx, dy, dz)^2 ./ D / 8.1
    nt = cld(ttot, dt)
    iter_max = 1e5

    # array allocation
    center = [lx / 2, ly / 2, lz / 2]
    qx = @zeros(nx - 1, ny - 2, nz - 2)
    qy = @zeros(nx - 2, ny - 1, nz - 2)
    qz = @zeros(nx - 2, ny - 2, nz - 1)
    Ht = @zeros(nx, ny, nz)
    Ht .= init_local_gaussian(center, dx, dy, dz, Ht)
    apply_boundary_conditions!(Ht, coords, dims)
    Hτ = copy(Ht)
    Hτ2 = @zeros(nx, ny, nz)
    residual_H = @zeros(nx, ny, nz)

    H_g = zeros(nx * dims[1], ny * dims[2], nz * dims[3])

    _dt = 1.0 / dt
    _dx = 1.0 / dx
    _dy = 1.0 / dy
    _dz = 1.0 / dz
    D_dx = D / dx
    D_dy = D / dy
    D_dz = D / dz

    iter_outer = 0
    timed_iter_total = 0  # without warmup
    tic = 0

    # GPU shared_memory setup
    if use_shared_memory
        threads = (32, 4, 4)  # these could be fined tuned still
        blocks  = ((nx, ny, nz) .- 1) .÷ threads .+ 1
        shmem = prod(threads.+2)*sizeof(Float64)
    end

    # for t ∈ ProgressBar(0:dt:ttot-dt)
    for t ∈ (0:dt:ttot-dt)
        if verbose && me == 0
            println("Iter: $(iter_outer)")
        end
        if iter_outer == 3  # manual warmup
            if verbose && me == 0
                println("Starting to measure")
            end
            tic = time()
            timed_iter_total = 0
        end
        iter_inner = 0
        err = 2 * tol
        while err > tol && iter_inner < iter_max
            if use_shared_memory
                @parallel blocks threads shmem=shmem diffusion_3D_step_τ_shared_memory(Ht, Hτ, Hτ2, residual_H, dτ, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz)
                update_halo!(Hτ)
            else
                # currently `@hide_communication` doesn't work with shared memory
                @hide_communication (8, 8, 8) begin
                    @parallel diffusion_3D_step_τ(Ht, Hτ, Hτ2, residual_H, dτ, _dt, _dx, _dy, _dz, D_dx, D_dy, D_dz)
                    update_halo!(Hτ)
                end
            end
            Hτ, Hτ2 = Hτ2, Hτ  # pointer swap, i.e. no memory movement
            err = dist_norm_L2(residual_H * dt, comm_cart) / sqrt(total_N)  # ≈ (1 * read & 1 * pow + 1 * add) * (nx-2) * (ny-2) * (nz-2)
            iter_inner += 1
        end
        if verbose && me == 0
            if err <= tol
                println("Converged after $iter_inner iterations.")
            else
                println("Couldn't converge within $iter_max iterations.")
            end
        end
        timed_iter_total += iter_inner
        iter_outer += 1
        Ht .= Hτ  # ≈ 1 * read/write + 1 * read * (nx-2) * (ny-2) * (nz-2), but these are "few" compared to the inner iter
    end
    toc = time()
    Δt = toc - tic  # seconds

    # /local/ benchmark results * number of MPI ranks
    nranks = MPI.Comm_size(MPI.COMM_WORLD)
    Work = nranks * timed_iter_total * (25+2) * (nx-2) * (ny-2) * (nz-2)  # see kernel docstring
    Performance = Work / Δt
    Memory = nranks * (use_shared_memory ?
                  timed_iter_total * ( 6+1) * sizeof(Float64) * (nx-2) * (ny-2) * (nz-2)
                : timed_iter_total * (14+1) * sizeof(Float64) * (nx-2) * (ny-2) * (nz-2))

    Intensity = Work / Memory
    Throughput = Memory / Δt

    println("Finished after $(iter_outer-2) outer iterations in $(@sprintf("%3.3f", Δt)) seconds of compute!")

    X_g = LinRange(0 + dx / 2, lx - dx / 2, nx * dims[1])

    gather!(Array(Ht), H_g)

    finalize_global_grid(; finalize_MPI = init_and_finalize_MPI)

    return X_g, H_g, BenchResults(Δt, Work, Performance, Memory, Intensity, Throughput)
end


text:
return Dict{String, String}(
    "motor" => "models/SLICOT/motor/motor.jl",
    "building" => "models/SLICOT/building/building.jl",
    "pde" => "models/SLICOT/pde/pde.jl",
    "heat" => "models/SLICOT/heat/heat.jl",
    "iss" => "models/SLICOT/iss/iss.jl",
    "beam" => "models/SLICOT/beam/beam.jl",
    "mna1" => "models/SLICOT/mna1/mna1.jl",
    "fom" => "models/SLICOT/fom/fom.jl",
    "mna5" => "models/SLICOT/mna5/mna5.jl"
    )

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module M2_0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    -3.53674108537383e-11*l1.^2 + 7.19493382490823e-12*l1.*l2 + 6.38787038845961e-11*l1.*l3 - 5.19058379413529e-14*l1 + 8.86717956244728e-12*l2.^2 - 2.4932305151384e-11*l2.*l3 + 7.55074773517067e-16*l2 - 1.96418867892903e-11*l3.^2 + 5.13914262466604e-14*l3 - 1.85095554911322e-17
end
end


text: println("Hello world with Julia!")


text:
using JOLI # remove when not needed any more to avoid reloading JOLI

# Replace all instances of joOPER string with name of your operator.
# Once you define helper function in module below
# you wil refer to them later as joOPER_etc.function_name
# where joOPER_etc is the name of module with helper functions.

## helper module
module joOPER_etc
    using JOLI: jo_convert, LocalVector
    # if your functions here can handle multi-vectors (matricies)
    # use LocalMatrix above and below instead of LocalVector

    # put your helper functions here and do not export them
    # refer to them outside of the module as joOPER_etc.function_name
    # in each function use jo_convert to convert your final output to
    # range of the operator; jo_convert(RDT,your_output_array)
    # .... e.g.
    function frwrd(rdt::DataType,x::LocalVector,m::Integer,otherargs...)
        y=Vector{rdt}(undef,m) # if needed
        # do what is needed here
        jo_convert(rdt,y,false)
    end
    function trnsp(rdt::DataType,x::LocalVector,m::Integer,otherargs...)
        y=Vector{rdt}(undef,m) # if needed
        # do what is needed here
        jo_convert(rdt,y,false)
    end
    function adjnt(rdt::DataType,x::LocalVector,m::Integer,otherargs...)
        y=Vector{rdt}(undef,m) # if needed
        # do what is needed here
        jo_convert(rdt,y,false)
    end
    function cnjgt(rdt::DataType,x::LocalVector,m::Integer,otherargs...)
        y=Vector{rdt}(undef,m) # if needed
        # do what is needed here
        jo_convert(rdt,y,false)
    end
end
using .joOPER_etc

export joOPER
# replace anything following % in doc with appropriate information
"""
    julia> op = joOPER(% typeless arguments)

% Description

# Signature

    % full function signature

# Arguments

- % argument list -> `arg`: meaning
- optional
    - % argument list -> `arg`: meaning
- keywords
    - % argument list -> `arg`: meaning (default ???)
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- % extra information

# Examples

% description

    % example

examples with DDT/RDT

    % joOPER(...; DDT=...)
    % joOPER(...; DDT=...,RDT=...)

"""
function joOPER(n::Integer,your_otherargs...;
                DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joOPER",your_otherkwargs...)
    # add/modify signature above as needed, but keep DDT,RDT,name
    # do what is needed here
    # like e.g. calculate m if depending on n
    m=n
    fWVok=false # or true if your functions can handle multi-vectors
    return joLinearFunctionFwd(m,n,
        v1->joOPER_etc.frwrd(RDT,v1,m),
        v2->joOPER_etc.trnsp(DDT,v2,m),
        v3->joOPER_etc.adjnt(DDT,v3,m),
        v4->joOPER_etc.cnjgt(RDT,v4,m),
        DDT,RDT;
        name=name,
        fMVok=fWVok
        )
end


text:
# EFTfitter.jl - Advanced Tutorial
# This tutorial introduces some of the advanced functionalities of EFTfitter.jl
# using a generic example. Please see the [tutorial] for basic usage of EFTfitter.

# ============= Parameters =============================================#
# We use the same parameters & priors as in the basic tutorial:
parameters = BAT.NamedTupleDist(
    C1 = -3..3, # short for: Uniform(-3, 3)
    C2 = Normal(0, 0.5) # Normal distribution
)

# ============= Observables =============================================#
# We use the same observables as in the basic tutorial. However, we use a different
# way for creating the vector of functions for the MeasurementDistribution.

function xsec1(params)
    coeffs = [20.12, 5.56, 325.556]
    return myfunc(params, coeffs)
end

function xsec2(params)
    coeffs = [2.12, 4.3, 12.6]
    return myfunc(params, coeffs)
end

function myfunc(params, c)
    return c[1] * params.C1 + c[2] * params.C1 * params.C2+ c[3] * params.C2
end

# When using distributions of measurements, a vector of functions with the predictions
# for the observable needs to be passed containing a function for each of the bins which
# have only the model parameters as their argument. Defining a separate function for each
# bin can, however, become tedious for a large number of bins, especially since typically
# the bins of a distribution have a similar functional dependence on the model parameters
# and only differ in some coefficients. In such cases, it is possible to use Julia's
# [metaprogramming](https://docs.julialang.org/en/v1/manual/metaprogramming/) features to
# create the vector of functions. The distribution in our [basic tutorial](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/tutorial/)
# has been defined by implementing three functions that all call the same function `myfunc`
# but with different values for the coefficients
# The same result can also be achieved like this:

function get_coeffs(i) # return the coefficients for bin i
    coeffs = [[2.2, 5.5, 6.6], [2.2, 5.5, 6.6], [2.2, 5.5, 6.6]]
    return coeffs[i]
end

function my_dist_func(params, i)
    coeffs = get_coeffs(i)
    return coeffs[1] * params.C1 + coeffs[2] * params.C1 * params.C2+ coeffs[3] * params.C2
end

# create an array of Functions with names `diff_xsec_binX`:
diff_xsec=Function[]
for i in 1:3
    @eval begin
        function $(Symbol("diff_xsec_bin$i"))(params)
            return my_dist_func(params, $i)
        end
        push!(diff_xsec, $(Symbol("diff_xsec_bin$i")))
    end
end

# ============= Measurements =============================================#
# Information about the uncertainties of measurements need to be provided to EFTfitter.jl
# in terms of the uncertainty values and corresponding correlation matrices.
# If you have these information in terms of covariance matrices, you need to convert it
# to correlation matrices and uncertainty values before. The function
# [`cov_to_cor`](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/api/#EFTfitter.cov_to_cor-Tuple{Array{var%22#s58%22,2}%20where%20var%22#s58%22%3C:Real})
# can be used for this:

cov_syst = [3.24   0.81   0.378  0.324  0.468;
            0.81   0.81   0.126  0.162  0.234;
            0.378  0.126  0.49   0.126  0.182;
            0.324  0.162  0.126  0.81   0.234;
            0.468  0.234  0.182  0.234  1.69]

cor_syst, unc_syst = cov_to_cor(cov_syst)

measurements = (

    Meas1 = Measurement(xsec1, 21.6,
            uncertainties = (stat=0.8, syst=unc_syst[1], another_unc=2.3)),

    Meas2 = Measurement(Observable(xsec2, min=0), 1.9,
            uncertainties = (stat=0.6, syst=unc_syst[2], another_unc=1.1), active=true),

    MeasDist = MeasurementDistribution(diff_xsec, [1.9, 2.93, 4.4],
                uncertainties = (stat = [0.7, 1.1, 1.2], syst= unc_syst[3:5], another_unc = [1.0, 1.2, 1.9]),
                active=[true, false, true]),
)

# ============= Correlations =============================================#
# The correlations are again the same as in the basic tutorial.
dist_corr = [1.0 0.5 0.0;
             0.5 1.0 0.0;
             0.0 0.0 1.0]

another_corr_matrix = to_correlation_matrix(measurements,
    (:Meas1, :Meas2, 0.4),
    (:Meas1, :MeasDist, 0.1),
    (:MeasDist, :MeasDist, dist_corr),
    (:MeasDist_bin2, :MeasDist_bin3, 0.3),
)

correlations = (
    stat = NoCorrelation(active=true), # will use the identity matrix of the correct size

    syst = Correlation([1.0 0.5 0.3 0.2 0.2;
                        0.5 1.0 0.2 0.2 0.2;
                        0.3 0.2 1.0 0.2 0.2;
                        0.2 0.2 0.2 1.0 0.2;
                        0.2 0.2 0.2 0.2 1.0], active=true), # `active = false`: ignore all uncertainty values and correlations for this type of uncertainty

    another_unc = Correlation(another_corr_matrix, active=true)
)

# ============= Nuisance Correlations =============================================#
# When performing an analysis with unknown correlation coefficients, it is possible
# to treat them as nuisance parameters in the fit.
# For this, we define a further `NamedTuple` consisting of `NuisanceCorrelation` objects:

nuisance_correlations = (
    ρ1  = NuisanceCorrelation(:syst, :Meas1, :Meas2, -1..1),
    ρ2  = NuisanceCorrelation(:syst, :MeasDist_bin1, :MeasDist_bin3, truncated(Normal(0.5, 0.1), -1, 1)),
)

# In the `NuisanceCorrelation` object we specify the name of the uncertainty type, the
# names of the two measurements we want to correlate using the nuisance correlations and
# a prior for the nuisance parameter. Note that the nuisance parameters should only be
# varied in the interval (-1, 1) as they represent correlation coefficients.
# For `ρ1` we choose a flat prior between -1 and 1. For `ρ2` we have some expectations and
# formulate them using a Gaussian prior with μ=0.5 and σ=0.1. However, to ensure that `ρ2`
# is only varied in the allowed region of (-1, 1), we [truncate](https://juliastats.org/Distributions.jl/stable/truncate/#Distributions.truncated)
# the normal distribution accordingly.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
######## CONSTRAINTS ############

# Generic fallback functions
function get_startup_shutdown(
    device,
    ::Type{<:VariableType},
    ::Type{<:AbstractDeviceFormulation},
) #  -> Union{Nothing, NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}}
    nothing
end

function get_min_max_limits(
    device,
    ::Type{<:VariableType},
    ::Type{<:AbstractDeviceFormulation},
) #  -> Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
    return (min=0.0, max=0.0)
end

@doc raw"""
Constructs min/max range constraint from device variable.

If min and max within an epsilon width:

``` variable[name, t] == limits.max ```

Otherwise:

``` limits.min <= variable[name, t] <= limits.max ```

where limits in constraint_infos.

# LaTeX

`` x = limits^{max}, \text{ for } |limits^{max} - limits^{min}| < \varepsilon ``

`` limits^{min} \leq x \leq limits^{max}, \text{ otherwise } ``
"""
function add_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: VariableType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function add_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: RangeConstraintLBExpressions,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function add_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: RangeConstraintUBExpressions,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function add_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ComponentActivePowerVariableLimitsConstraint,
    U <: VariableType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_parameterized_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
    )
    return
end

function _add_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation}
    time_steps = get_time_steps(container)
    device_names = [PSY.get_name(d) for d in devices]

    con_lb =
        add_constraints_container!(container, T(), V, device_names, time_steps, meta="lb")

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type
        con_lb[ci_name, t] =
            JuMP.@constraint(container.JuMPmodel, array[ci_name, t] >= limits.min)
    end
    return
end

function _add_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation}
    time_steps = get_time_steps(container)
    device_names = [PSY.get_name(d) for d in devices]

    con_ub =
        add_constraints_container!(container, T(), V, device_names, time_steps, meta="ub")

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type
        con_ub[ci_name, t] =
            JuMP.@constraint(container.JuMPmodel, array[ci_name, t] <= limits.max)
    end
    return
end

@doc raw"""
Constructs min/max range constraint from device variable and on/off decision variable.

If device min = 0:

``` varcts[name, t] <= limits.max*varbin[name, t]) ```

``` varcts[name, t] >= 0.0 ```

Otherwise:

``` varcts[name, t] <= limits.max*varbin[name, t] ```

``` varcts[name, t] >= limits.min*varbin[name, t] ```

where limits in constraint_infos.

# LaTeX

`` 0 \leq x^{cts} \leq limits^{max} x^{bin}, \text{ for } limits^{min} = 0 ``

`` limits^{min} x^{bin} \leq x^{cts} \leq limits^{max} x^{bin}, \text{ otherwise } ``
"""
function add_semicontinuous_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: VariableType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_semicontinuous_lower_bound_range_constraints_impl!(
        container,
        T,
        array,
        devices,
        model,
    )
    _add_semicontinuous_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        devices,
        model,
    )
    return
end

function add_semicontinuous_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: RangeConstraintLBExpressions,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_semicontinuous_lower_bound_range_constraints_impl!(
        container,
        T,
        array,
        devices,
        model,
    )
    return
end

function add_semicontinuous_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: RangeConstraintUBExpressions,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_semicontinuous_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        devices,
        model,
    )
    return
end

function _add_semicontinuous_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [OnVariable()]

    con_lb = add_constraints_container!(container, T(), V, names, time_steps, meta="lb")

    @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1"
    varbin = get_variable(container, only(binary_variables), V)

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type
        con_lb[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] >= limits.min * varbin[ci_name, t]
        )
    end
    return
end

function _add_semicontinuous_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {T <: ConstraintType, V <: PSY.Component, W <: AbstractDeviceFormulation}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [OnVariable()]

    con_ub = add_constraints_container!(container, T(), V, names, time_steps, meta="ub")

    @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1"
    varbin = get_variable(container, only(binary_variables), V)

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W) # depends on constraint type and formulation type
        con_ub[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] <= limits.max * varbin[ci_name, t]
        )
    end
    return
end

@doc raw"""
Constructs min/max range constraint from device variable and reservation decision variable.

``` varcts[name, t] <= limits.max * (1 - varbin[name, t]) ```

``` varcts[name, t] >= limits.min * (1 - varbin[name, t]) ```

where limits in constraint_infos.

# LaTeX

`` 0 \leq x^{cts} \leq limits^{max} (1 - x^{bin}), \text{ for } limits^{min} = 0 ``

`` limits^{min} (1 - x^{bin}) \leq x^{cts} \leq limits^{max} (1 - x^{bin}), \text{ otherwise } ``
"""
function add_reserve_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: InputActivePowerVariableLimitsConstraint,
    U <: VariableType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_reserve_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_reserve_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function add_reserve_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: InputActivePowerVariableLimitsConstraint,
    U <: ExpressionType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), W)
    _add_reserve_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_reserve_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function _add_reserve_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
) where {
    T <: InputActivePowerVariableLimitsConstraint,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [ReservationVariable()]

    IS.@assert_op length(binary_variables) == 1
    varbin = get_variable(container, only(binary_variables), V)

    names = [PSY.get_name(x) for x in devices]
    # MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    # In the future this can be updated
    con_lb = add_constraints_container!(container, T(), V, names, time_steps, meta="lb")

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W)
        con_lb[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] >= limits.min * (1 - varbin[ci_name, t])
        )
    end
    return
end

function _add_reserve_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {
    T <: InputActivePowerVariableLimitsConstraint,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [ReservationVariable()]

    IS.@assert_op length(binary_variables) == 1
    varbin = get_variable(container, only(binary_variables), V)

    names = [PSY.get_name(x) for x in devices]
    # MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    # In the future this can be updated
    con_ub = add_constraints_container!(container, T(), V, names, time_steps, meta="ub")

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, W)
        con_ub[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] <= limits.max * (1 - varbin[ci_name, t])
        )
    end
    return
end

@doc raw"""
Constructs min/max range constraint from device variable and reservation decision variable.

``` varcts[name, t] <= limits.max * varbin[name, t] ```

``` varcts[name, t] >= limits.min * varbin[name, t] ```

where limits in constraint_infos.

# LaTeX

`` limits^{min} x^{bin} \leq x^{cts} \leq limits^{max} x^{bin},``
"""
function add_reserve_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{W},
    model::DeviceModel{W, X},
    ::Type{Y},
) where {
    T <:
    Union{ReactivePowerVariableLimitsConstraint, OutputActivePowerVariableLimitsConstraint},
    U <: VariableType,
    W <: PSY.Component,
    X <: AbstractDeviceFormulation,
    Y <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), W)
    _add_reserve_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_reserve_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

@doc raw"""
Constructs min/max range constraint from device variable and reservation decision variable.

``` varcts[name, t] <= limits.max * varbin[name, t] ```

``` varcts[name, t] >= limits.min * varbin[name, t] ```

where limits in constraint_infos.

# LaTeX

`` limits^{min} x^{bin} \leq x^{cts} \leq limits^{max} x^{bin},``
"""
function add_reserve_range_constraints!(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{W},
    model::DeviceModel{W, X},
    ::Type{Y},
) where {
    T <:
    Union{ReactivePowerVariableLimitsConstraint, OutputActivePowerVariableLimitsConstraint},
    U <: ExpressionType,
    W <: PSY.Component,
    X <: AbstractDeviceFormulation,
    Y <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), W)
    _add_reserve_upper_bound_range_constraints_impl!(container, T, array, devices, model)
    _add_reserve_lower_bound_range_constraints_impl!(container, T, array, devices, model)
    return
end

function _add_reserve_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{W},
    ::DeviceModel{W, X},
) where {
    T <:
    Union{ReactivePowerVariableLimitsConstraint, OutputActivePowerVariableLimitsConstraint},
    W <: PSY.Component,
    X <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [ReservationVariable()]

    con_lb = add_constraints_container!(container, T(), W, names, time_steps, meta="lb")

    @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1"
    varbin = get_variable(container, only(binary_variables), W)

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, X) # depends on constraint type and formulation type
        con_lb[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] >= limits.min * varbin[ci_name, t]
        )
    end
    return
end

function _add_reserve_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    devices::IS.FlattenIteratorWrapper{W},
    ::DeviceModel{W, X},
) where {
    T <:
    Union{ReactivePowerVariableLimitsConstraint, OutputActivePowerVariableLimitsConstraint},
    W <: PSY.Component,
    X <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    binary_variables = [ReservationVariable()]

    con_ub = add_constraints_container!(container, T(), W, names, time_steps, meta="ub")

    @assert length(binary_variables) == 1 "Expected $(binary_variables) for $U $V $T $W to be length 1"
    varbin = get_variable(container, only(binary_variables), W)

    for device in devices, t in time_steps
        ci_name = PSY.get_name(device)
        limits = get_min_max_limits(device, T, X) # depends on constraint type and formulation type
        con_ub[ci_name, t] = JuMP.@constraint(
            container.JuMPmodel,
            array[ci_name, t] <= limits.max * varbin[ci_name, t]
        )
    end
    return
end

function add_parameterized_lower_bound_range_constraints(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    ::Type{P},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: ExpressionType,
    P <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_parameterized_lower_bound_range_constraints_impl!(
        container,
        T,
        array,
        P,
        devices,
        model,
    )
    return
end

function add_parameterized_lower_bound_range_constraints(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    ::Type{P},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: VariableType,
    P <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_parameterized_lower_bound_range_constraints_impl!(
        container,
        T,
        array,
        P,
        devices,
        model,
    )
    return
end

# This function is re-used in SemiContinuousFeedforward
function lower_bound_range_with_parameter!(
    jump_model::JuMP.Model,
    constraint_container::JuMPConstraintArray,
    lhs_array,
    param_multiplier::JuMPFloatArray,
    param_array::Union{JuMPParamArray, JuMPFloatArray},
    devices::IS.FlattenIteratorWrapper{V},
) where {V <: PSY.Component}
    time_steps = axes(constraint_container)[2]
    for device in devices, t in time_steps
        name = PSY.get_name(device)
        constraint_container[name, t] = JuMP.@constraint(
            jump_model,
            lhs_array[name, t] >= param_multiplier[name, t] * param_array[name, t]
        )
    end
    return
end

function _add_parameterized_lower_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    ::Type{U},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {
    T <: ConstraintType,
    U <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]

    constraint = add_constraints_container!(container, T(), V, names, time_steps, meta="lb")

    parameter = get_parameter_array(container, U(), V)
    multiplier = get_parameter_multiplier_array(container, U(), V)
    jump_model = get_jump_model(container)
    lower_bound_range_with_parameter!(
        jump_model,
        constraint,
        array,
        multiplier,
        parameter,
        devices,
    )
    return
end

function add_parameterized_upper_bound_range_constraints(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    ::Type{P},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: ExpressionType,
    P <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_expression(container, U(), V)
    _add_parameterized_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        P,
        devices,
        model,
    )
    return
end

function add_parameterized_upper_bound_range_constraints(
    container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
    ::Type{P},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {
    T <: ConstraintType,
    U <: VariableType,
    P <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
    X <: PM.AbstractPowerModel,
}
    array = get_variable(container, U(), V)
    _add_parameterized_upper_bound_range_constraints_impl!(
        container,
        T,
        array,
        P,
        devices,
        model,
    )
    return
end

# This function is re-used in SemiContinuousFeedforward
function upper_bound_range_with_parameter!(
    jump_model::JuMP.Model,
    constraint_container::JuMPConstraintArray,
    lhs_array,
    param_multiplier::JuMPFloatArray,
    param_array::Union{JuMPParamArray, JuMPFloatArray},
    devices::IS.FlattenIteratorWrapper{V},
) where {V <: PSY.Component}
    time_steps = axes(constraint_container)[2]
    for device in devices, t in time_steps
        name = PSY.get_name(device)
        constraint_container[name, t] = JuMP.@constraint(
            jump_model,
            lhs_array[name, t] <= param_multiplier[name, t] * param_array[name, t]
        )
    end
    return
end

function _add_parameterized_upper_bound_range_constraints_impl!(
    container::OptimizationContainer,
    ::Type{T},
    array,
    ::Type{P},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
) where {
    T <: ConstraintType,
    P <: ParameterType,
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]

    constraint = add_constraints_container!(container, T(), V, names, time_steps, meta="ub")

    parameter = get_parameter_array(container, P(), V)
    multiplier = get_parameter_multiplier_array(container, P(), V)
    jump_model = get_jump_model(container)
    upper_bound_range_with_parameter!(
        jump_model,
        constraint,
        array,
        multiplier,
        parameter,
        devices,
    )
    return
end


text:
module qccdSimulator
export readJSON, readTimeJSON

include("./types/description.jl")
include("./types/control.jl")
include("./utils/qccdevctrlUtils.jl")
include("./utils/qccdevctrlChecks.jl")
include("./qccdevcontrol.jl")

using .QCCDevDes_Types
using JSON3

"""
Creates an object QCCDevDescription from JSON.
Throws ArgumentError an error if input is not a valid file.
"""
function readJSON(path::String)::QCCDevDescription
    if !isfile(path)
        throw(ArgumentError("Input is not a file"))
    end
    # Parsing JSON
    return topology::QCCDevDescription  = try
        JSON3.read(read(path, String), QCCDevDescription)
    catch err
        throw(ArgumentError(err.msg))
    end
end

"""
Fills OperationTimes global variable from JSON file.
Throws ArgumentError if input is not a valid file.
Throws ArgumentError if there is a negative time value.
Throws ArgumentError if file contents are not of the format Dict{String, Int}.
"""
function readTimeJSON(path ::String)
    if !isfile(path)
        throw(ArgumentError("Input is not a file"))
    end
    # Parsing JSON
    times = try
        JSON3.read(read(path, String), Dict{Symbol, Int64})
    catch err
        throw(ArgumentError(err.msg))
    end
    isempty(filter(v -> v < 0, collect(values(times)))) ||
        throw(ArgumentError("Time values can't be negative"))

    setOperationTimes(times)
end

end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
using Base.Test
using Compat

using ArbFloats


text:
#!/usr/bin/env julia
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "align_short_reads"
    help = "Align short reads with minimap2 & generate consensus fasta."
    action = :command

    "align_assembly"
    help = "Align assembly with minimap2 & generate re-aligned fasta."
    action = :command

    "merge_alignments"
    help = "Merge several samples into a combined fasta and compute distance matrix."
    action = :command
end

@add_arg_table s["align_short_reads"] begin
    "out_prefix"
    help = "Prefix for output files"
    arg_type = String
    required = true

    "ref"
    help = "Reference fasta. Must be gzipped and indexed by samtools faidx."
    arg_type = String
    required = true

    "fastqs"
    help = "Fastq files to align to reference."
    nargs = '+'
    arg_type = String
    required = true

    "--min_ac"
    help = "Minimum read count of consensus allele."
    arg_type = Int
    default = 10

    "--min_af"
    help = "Minimum frequency of consensus allele."
    arg_type = Float64
    default = 0.9

    "--max_dp"
    help = "Maximum read depth; positions with higher depth are masked."
    arg_type = Int
    default = typemax(Int)

    "--threads"
    help = "Threads. Note: currently only used in alignment stage, but not in the pileup stage."
    arg_type = Int
    default = 3
end

@add_arg_table s["align_assembly"] begin
    "out_prefix"
    help = "Prefix for output files"
    arg_type = String
    required = true

    "ref"
    help = "Reference fasta. Must be gzipped and indexed by samtools faidx."
    arg_type = String
    required = true

    "query"
    help = "Query assembly."
    arg_type = String
    required = true

    "--preset"
    help = "Preset for minimap2"
    arg_type = String
    default = "asm5"

    "--min_ac"
    help = "Minimum read count of consensus allele."
    arg_type = Int
    default = 1

    "--min_af"
    help = "Minimum frequency of consensus allele."
    arg_type = Float64
    default = 1.0

    "--max_dp"
    help = "Maximum read depth; positions with higher depth are masked."
    arg_type = Int
    default = 1

    "--threads"
    help = "Threads. Note: currently only used in alignment stage, but not in the pileup stage."
    arg_type = Int
    default = 3
end

@add_arg_table s["merge_alignments"] begin
    "out_prefix"
    help = "Prefix for output files"
    arg_type = String
    required = true

    "fastas"
    help = "Fasta files to merge."
    nargs = '+'
    arg_type = String
    required = true

    "--fasta2sample_regex"
    help = "Regex for fasta path. Sample name is the first capturing group."
    default = "(?:.*/)?([^/]*)\\.fa"
end

parsed_args = parse_args(s)
cmd = parsed_args["%COMMAND%"]

# import Spid at the end, in case we're just calling --help
using Spid

if cmd == "align_short_reads"
    parsed_args = parsed_args[cmd]

    align_short_reads(
        parsed_args["ref"], parsed_args["fastqs"],
        string(parsed_args["out_prefix"], ".bam"),
        string(parsed_args["out_prefix"], ".fa"),
        min_ac=parsed_args["min_ac"], min_af=parsed_args["min_af"],
        max_dp=parsed_args["max_dp"],
        threads=parsed_args["threads"]
    )
elseif cmd == "align_assembly"
    parsed_args = parsed_args[cmd]

    align_assembly(
        parsed_args["ref"], parsed_args["query"],
        string(parsed_args["out_prefix"], ".bam"),
        string(parsed_args["out_prefix"], ".fa"),
        preset=parsed_args["preset"], min_ac=parsed_args["min_ac"],
        min_af=parsed_args["min_af"], max_dp=parsed_args["max_dp"],
        threads=parsed_args["threads"]
    )
elseif cmd == "merge_alignments"
    parsed_args = parsed_args[cmd]

    merge_alignments(
        parsed_args["out_prefix"],
        parsed_args["fastas"],
        Regex(parsed_args["fasta2sample_regex"])
    )
else
    println("Unrecognized command $cmd")
end


text:
Base.@kwdef mutable struct SSPStruct
  alpha = nothing
end
Base.@kwdef mutable struct RosenbrockStruct
  transformed = nothing
  nStage = nothing
  alpha = nothing
  Gamma = nothing
  b = nothing
  a = nothing
  c = nothing
  m = nothing
  d = nothing
  b2 = nothing
  beta0 = nothing
  beta = nothing
  SSP = SSPStruct(;)
end

function RosenbrockMethod(Method)
  str = Method
  ROS = RosenbrockStruct(;)
if str == "SSP-Knoth"
    ROS.transformed=true;
    ROS.nStage=3;
    ROS.alpha=[0 0 0
      1 0 0
      1/4 1/4 0];
    ROS.d=ROS.alpha*ones(ROS.nStage,1);
    ROS.b=[1/6 1/6 2/3];
    ROS.Gamma=[1 0 0
      0 1 0
      -3/4 -3/4 1];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.SSP.alpha=[1 0 0
                   3/4 1/4 0
                   1/3 0 2/3];
elseif str == "RK3_H"
    ROS.transformed=true;
    ROS.nStage=3;
    ROS.alpha=[0 0 0
      1/3 0 0
      0 1/2 0];
    g=(3.0+sqrt(3.0))/6.0;
    ROS.Gamma=[g 0 0
               (1-12*g^2)/(-9+36*g) g 0
               -1/4+2*g 1/4-3*g g];
    ROS.b=[0,0,1];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.d=ROS.Gamma[1,1];
elseif str == "RODAS_N"
    ROS.transformed=false;
    ROS.nStage=4;
    ROS.alpha=[0 0 0 0
      0 0 0 0
      1. 0 0 0
      0.75 -0.25 0.5 0];
    ROS.Gamma=[0.5 0 0 0
      1. 0.5 0 0
      -0.25 -0.25 0.5 0
      1.0 / 12 1.0 / 12 -2. / 3 0.5];
    ROS.b  = [5. / 6 -1.0 / 6 -1.0 / 6 0.5];
elseif str == "RODAS"
    ROS.transformed=true;
    ROS.nStage=4;
    ROS.alpha=[0 0 0 0
      0 0 0 0
      1. 0 0 0
      0.75 -0.25 0.5 0];
    ROS.Gamma=[0.5 0 0 0
      1. 0.5 0 0
      -0.25 -0.25 0.5 0
      1.0 / 12 1.0 / 12 -2. / 3 0.5];
    ROS.b  = [5. / 6 -1.0 / 6 -1.0 / 6 0.5];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.d=ROS.Gamma[1,1];
elseif str == "TSROSWSANDU3_N"
    ROS.transformed=false;
    ROS.nStage=3;
    ROS.alpha=[0 0 0
      0.43586652150845899941601945119356 0 0
      0.43586652150845899941601945119356 0 0];
    ROS.Gamma=[0.43586652150845899941601945119356 0 0
      -0.19294655696029095575009695436041 0.43586652150845899941601945119356 0
      0 1.74927148125794685173529749738960 0.43586652150845899941601945119356];
    ROS.b=[-0.75457412385404315829818998646589,1.94100407061964420292840123379419
      ,-0.18642994676560104463021124732829];
elseif str == "TSROSWSANDU3"
    ROS.transformed=true;
    ROS.nStage=3;
    ROS.alpha=[0 0 0
      0.43586652150845899941601945119356 0 0
      0.43586652150845899941601945119356 0 0];
    ROS.Gamma=[0.43586652150845899941601945119356 0 0
      -0.19294655696029095575009695436041 0.43586652150845899941601945119356 0
      0 1.74927148125794685173529749738960 0.43586652150845899941601945119356];
    ROS.b=[-0.75457412385404315829818998646589,1.94100407061964420292840123379419
      ,-0.18642994676560104463021124732829];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.d=ROS.Gamma[1,1];
elseif str == "TROSWLASSP3P4S2C"
    r=5;
    ROS.transformed=false;
    ROS.nStage=4;
    ROS.alpha = [0 0 0 0
      1.0 / 2. 0 0 0
      1.0 / 2. 1.0 / 2. 0 0
      1.0 / 6. 1.0 / 6. 1.0 / 6. 0 ];
    ROS.Gamma= [1.0 / 2. 0 0 0
      0.0 3. / 4. 0 0
      -2. / 3. -23. / 9. 2. / 9. 0
      1.0 / 18. 65. / 108. -2. / 27 0];
    ROS.b = [1.0 / 6.,1.0 / 6.,1.0 / 6.,1.0 / 2.];
    ROS.b2= [3. / 16.,10. / 16.,3. / 16.,0];

    #ROS.m
elseif str == "ROSEul"
    ROS.nStage=1;

elseif str == "ROS3Pw"
    ROS.nStage=3;
    ROS.beta0=7.88675134594812865529e-01;
    ROS.alpha[1,1]=2.0e0;
    ROS.alpha[2,1]=6.33974596215561403412e-01;
    ROS.alpha[3,1]=1.63397459621556140341e+00;
    ROS.alpha[3,2]=2.94228634059947813384e-01;
    ROS.alpha[3,3]=1.07179676972449078320e+00;

    ROS.beta[2,1]=-2.53589838486224561365e+00;
    ROS.beta[3,1]=-1.62740473580835520728e+00;
    ROS.beta[3,2]=-2.74519052838329002952e-01;
    ROS.alpha=ROS.alpha*ROS.beta0;
    ROS.beta=ROS.beta*ROS.beta0;
    ROS.d=ROS.beta0;
elseif str == "ROS-AMF"
    ROS.nStage=2;
    ROS.transformed=true;
    d=1/3+1/6*sqrt(3);

    ROS.alpha=[0 0
               2/3 0];
    ROS.Gamma=[d 0
               -4/3*d d];
    ROS.b=[1/4,3/4];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.d=ROS.Gamma[1,1];
elseif str == "ROS-AMF_N"
    ROS.nStage=2;
    ROS.transformed=false;
    d=1/3+1/6*sqrt(3);

    ROS.alpha=[0 0
               2/3 0];
    ROS.Gamma=[d 0
               -4/3*d d];
    ROS.b=[1/4,3/4];
    ROS.a=ROS.alpha/ROS.Gamma;
    ROS.c=-inv(ROS.Gamma);
    ROS.m=ROS.b/ROS.Gamma;
    ROS.d=ROS.Gamma[1,1];
elseif str == "ROS2"
    ROS.nStage=2;
    # Matrix form aftrer Wolfbrandt
    #     d=.5;
    #     a(2,1)=2. / 3.;
    #     d(2,1)=.5;
    #     b(1)=.25;
    #     b(2)=.75
    alpha=zeros(ROS.nStage,ROS.nStage);
    gamma=zeros(ROS.nStage,ROS.nStage);
    ROS.d=0.5+sqrt(3)/6;
    d=ROS.d;
    alpha[2,1]=2/3;
    b[1,1]=0.25;
    b[1,2]=0.75;
    gamma[1,1]=d;
    gamma[2,1]=-d/b(1,2);
    gamma[2,2]=d;
    ROS.a=alpha/gamma;
    ROS.m=b/gamma;
    ROS.c=inv(diag(diag(gamma)))-inv(gamma);
elseif str == "ROSRK3"
    ROS.transformed=true;
    ROS.nStage=3;
    ROS.a=zeros(ROS.nStage,ROS.nStage);
    ROS.c=zeros(ROS.nStage,ROS.nStage);
    beta0=1;
    #beta0=1.0/2.0
    #beta0=(3.0+sqrt(3.0))/6.0;
    alpha = zeros(3,3) # TODO: Check with Oswald
    beta = zeros(3,2) # TODO: Check with Oswald
    alpha[1,1]=1.0/(3.0*beta0);
    alpha[2,1]=(-1.0 + 12.0*beta0^2.0)/(18.0*beta0^2.0*(-1.0 + 4.0*beta0));
    alpha[2,2]=1.0/(2.0*beta0);
    alpha[3,1]=   (1.0 - 3.0*beta0*(7.0 + 16.0*beta0*(-2.0 + 3.0*beta0)))/
      (36.0*beta0^3.0*(-1.0 + 4.0*beta0));
    alpha[3,2]=   (-1.0 + 12.0*beta0)/(4.0*beta0^2.0);
    alpha[3,3]=   1.0/beta0;
    beta[2,1]= (1.0 - 12.0*beta0^2.0)/(beta0^2.0*(-9.0 + 36.0*beta0));
    beta[3,1]= (-1.0 + 3.0*beta0*(7.0 + 16.0*beta0*(-2.0 + 3.0*beta0)))/
      (36. * beta0^3.0*(-1.0 + 4.0*beta0));
    beta[3,2]= (0.250 - 3.0*beta0)/beta0^2.0;
    alpha=alpha*beta0;
    beta=beta*beta0;
    ROS.a[2:3,1:2]=alpha[1:2,1:2];
    ROS.m=alpha[3,:];
    ROS.c=beta;
    ROS.d=1;
end
return ROS
end


text:
@testset "Construction" begin

    # basic constructor
    xs = rand(5)
    ys = rand(5)
    p = Lagrange(xs, ys)
    @test all(p.(xs) .≈  ys)

    # test weights
    @test all(p.ws .≈  SP.lagrange_barycentric_weights(xs))
    zs = rand(5)
    q = Lagrange(xs, p.ws, zs)
    @test all(q.(xs) .≈  zs)

    # fit is alternate
    q = fit(Lagrange, xs, ys)
    us = rand(5)
    @test all(p.(us) .≈ q.(us))

    # test one, zero, variable
    x = variable()
    q = fit(Lagrange, [1,2,3],[2,3,1])
    @test iszero(zero(q))
    @test convert(Polynomial, one(q)) == one(x)
    @test convert(Polynomial, variable(q)) == x

end

@testset "Conversion" begin

    xs = [1,2,3]
    ys = xs.^2
    p = fit(Lagrange, xs, ys)
    q = convert(Polynomial, p)
    @test coeffs(q)  == [0,0,1]

    pp = fit(Lagrange, xs, q)
    @test coeffs(pp) == coeffs(p)

end

@testset "Evaluation" begin

    xs = rand(5)
    ys = sin.(xs)
    p = Lagrange(xs, ys)
    @test all(p.(xs) == ys)   # shortcuts  for xs
    us = rand(5)
    @test all(abs.(p.(us) - sin.(us)) .<= 1e-1)

    xs, ws = SP.lagrange_barycentric_nodes_weights(Chebyshev, 64)
    ys = sin.(xs)
    p = Lagrange(xs, ws, ys)
    us = 0.1 .+ 0.8*rand(5) # avoid edges...
    @test all(abs.(sin.(us) - p.(us)) .<=  1e-3)

end

@testset "Arithmetic" begin

    p = fit(Lagrange, [1,2,3], x -> x^2)
    q = fit(Lagrange, [2,3,4,5], x -> x^3)

    us = rand(5)

    @test maximum(abs, p.(us) + q.(us) - (p + q).(us)) <= 1e-5
    @test maximum(abs, p.(us) .* q.(us) - (p * q).(us)) <= 1e-2  # this  can  be way  off for even moderate sized polynomials.

end

@testset "Elementwise operations" begin

    p = fit(Lagrange, [1,2,3], x -> x^2)
    p + 1
    -p
    2p

end

@testset "roots" begin

    q = fit(Lagrange, [2,3,4,5], x -> (x-1)*(x-1/2)*(x+1/2))
    @test all(roots(q) .≈ [-1/2, 1/2, 1])

end


text:
# # Tutorial
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/guide/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/guide/@__NAME__.ipynb)

#nb # DFTK is a Julia package for playing with plane-wave
#nb # density-functional theory algorithms. In its basic formulation it
#nb # solves periodic Kohn-Sham equations.
#
# This document provides an overview of the structure of the code
# and how to access basic information about calculations.
# Basic familiarity with the concepts of plane-wave density functional theory
# is assumed throughout. Feel free to take a look at the
#md # [Periodic problems](@ref periodic-problems)
#nb # [Periodic problems](https://juliamolsim.github.io/DFTK.jl/stable/guide/periodic_problems/)
# or the
#md # [density-functional theory](@ref density-functional-theory)
#nb # [density-functional theory](https://juliamolsim.github.io/DFTK.jl/stable/guide/density_functional_theory/)
# chapters for some introductory material on the topic.
#
# !!! note "Convergence parameters in the documentation"
#     We use rough parameters in order to be able
#     to automatically generate this documentation very quickly.
#     Therefore results are far from converged.
#     Tighter thresholds and larger grids should be used for more realistic results.
#
# For our discussion we will use the classic example of
# computing the LDA ground state of the
# [silicon crystal](https://www.materialsproject.org/materials/mp-149).
# Performing such a calculation roughly proceeds in three steps.

using DFTK
using Plots
using Unitful
using UnitfulAtomic

## 1. Define lattice and atomic positions
a = 5.431u"angstrom"          # Silicon lattice constant
lattice = a / 2 * [[0 1 1.];  # Silicon lattice vectors
                   [1 0 1.];  # specified column by column
                   [1 1 0.]]

# By default, all numbers passed as arguments are assumed to be in atomic
# units.  Quantities such as temperature, energy cutoffs, lattice vectors, and
# the k-point grid spacing can optionally be annotated with Unitful units,
# which are automatically converted to the atomic units used internally. For
# more details, see the [Unitful package
# documentation](https://painterqubits.github.io/Unitful.jl/stable/) and the
# [UnitfulAtomic.jl package](https://github.com/sostock/UnitfulAtomic.jl).

## Load HGH pseudopotential for Silicon
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))

## Specify type and positions of atoms
atoms = [Si => [ones(3)/8, -ones(3)/8]]

## 2. Select model and basis
model = model_LDA(lattice, atoms)
kgrid = [4, 4, 4]     # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 7              # kinetic energy cutoff
## Ecut = 190.5u"eV"  # Could also use eV or other energy-compatible units
basis = PlaneWaveBasis(model; Ecut, kgrid)
## Note the implicit passing of keyword arguments here:
## this is equivalent to PlaneWaveBasis(model; Ecut=Ecut, kgrid=kgrid)

## 3. Run the SCF procedure to obtain the ground state
scfres = self_consistent_field(basis, tol=1e-8);

# That's it! Now you can get various quantities from the result of the SCF.
# For instance, the different components of the energy:
scfres.energies

# Eigenvalues:
hcat(scfres.eigenvalues...)
# `eigenvalues` is an array (indexed by k-points) of arrays (indexed by
# eigenvalue number). The "splatting" operation `...` calls `hcat`
# with all the inner arrays as arguments, which collects them into a
# matrix.
#
# The resulting matrix is 7 (number of computed eigenvalues) by 8
# (number of k-points). There are 7 eigenvalues per k-point because
# there are 4 occupied states in the system (4 valence electrons per
# silicon atom, two atoms per unit cell, and paired spins), and the
# eigensolver gives itself some breathing room by computing some extra
# states (see `n_ep_extra` argument to `self_consistent_field`).
#
# We can check the occupations ...
hcat(scfres.occupation...)
# ... and density, where we use that the density objects in DFTK are
# indexed as ρ[iσ, ix, iy, iz], i.e. first in the spin component and then
# in the 3-dimensional real-space grid.
rvecs = collect(r_vectors(basis))[:, 1, 1]  # slice along the x axis
x = [r[1] for r in rvecs]                   # only keep the x coordinate
plot(x, scfres.ρ[1, :, 1, 1], label="", xlabel="x", ylabel="ρ", marker=2)

# We can also perform various postprocessing steps:
# for instance compute a band structure
plot_bandstructure(scfres; kline_density=10)
# or get the cartesian forces (in Hartree / Bohr)
compute_forces_cart(scfres)[1]  # Select silicon forces
# The `[1]` extracts the forces for the first kind of atoms,
# i.e. `Si` (silicon) in the setup of the `atoms` list of step 1 above.
# As expected, they are almost zero in this highly symmetric configuration.

# ## Where to go from here
# Take a look at the
#md # [example index](@ref example-index)
#nb # [example index](https://juliamolsim.github.io/DFTK.jl/dev/#example-index-1)
# to continue exploring DFTK.


text:
using Test
using QCMaterial

@testset "computation.apply_qubit_op" begin
    n_qubit = 2
    # Prepare |Psi> = |00>
    state = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, 0b00)

    # Apply c^dagger_1 to |Psi>, which will yield |10>.
    op = jordan_wigner(FermionOperator("1^"))

    # Prepare |phi> = |10>
    state0_bra = QulacsQuantumState(n_qubit)
    set_computational_basis!(state0_bra, 0b10)
    circuit_bra = uccgsd(n_qubit, orbital_rot=true, conserv_Sz_singles=false)

    squared_norm = apply_qubit_op(op, state, circuit_bra, state0_bra)

    # Verify the result
    @test squared_norm ≈ 1.
end

@testset "computation.apply_qubit_op2" begin
    n_qubit = 2
    # Prepare |Psi> = (|01> + |10>)/sqrt(2)
    state = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, 0b01)

    state2 = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, 0b10)

    #state + state2

end


text:
atol = rtol = 1e-4

m = 2
n = 2
row = Int32[1,2,2]
col = Int32[1,1,2]
val = Float64[1.,.1,2.]

coo = MadNLP.SparseMatrixCOO(m,n,row,col,val)
csc = sparse(row,col,val,m,n)
dense = Array(csc)

b = [1.0,3.0]
x = similar(b)

@testset "SparseMatrixCOO" begin
    @test MadNLP.size(coo) == (2,2)
    @test coo[1,1] == 1.
end

@testset "LAPACK" begin
    sol= [0.8542713567839195, 1.4572864321608041]
    M = MadNLPLapackCPU.Solver(dense)
    MadNLP.introduce(M)
    MadNLP.improve!(M)
    MadNLP.factorize!(M)
    MadNLP.is_inertia(M) && (MadNLP.inertia(M) = (2,0,0))
    x = MadNLP.solve!(M,copy(b))
    @test solcmp(x,sol)
end

macro test_linear_solver(name)
    str=string(name)
    quote
        if isdefined(MadNLP,Symbol($str))
            @testset $str begin
                sol= [0.8542713567839195, 1.4572864321608041]
                M = MadNLP.$name.Solver(csc)
                MadNLP.introduce(M)
                MadNLP.improve!(M)
                MadNLP.factorize!(M)
                MadNLP.is_inertia(M) && (MadNLP.inertia(M) = (2,0,0))
                x = MadNLP.solve!(M,copy(b))
                @test solcmp(x,sol)
            end
        end
    end
end

@test_linear_solver ma27
@test_linear_solver ma57
@test_linear_solver ma77
@test_linear_solver ma86
@test_linear_solver ma97
@test_linear_solver pardiso
@test_linear_solver pardisomkl
@test_linear_solver umfpack
@test_linear_solver mumps


text:
#---------------------------------------------------------------
#---------------------------------------------------------------
# Differential Tags
type DTag
    tag::Array
end

DTag(x) = DTag([x])

DTag(x...) = DTag([i for i in x] |> sort)

Base.length(t::DTag) = length(t.tag)
Base.:(==)(x::DTag,y::DTag) = x.tag == y.tag
Base.getindex(t::DTag, i::Int) = (t.tag)[i]

function Base.isless(t1::DTag,t2::DTag)
    if length(t1) == length(t2)
        for i in 1:length(t1)
            if t1[i] >= t2[i]
                return false
            end
        end
        t1 == t2 ? (return false ) : (return true)
    elseif length(t1) < length(t2)
        true
    else
        false
    end
end

function hasDuplicates(arr)
    for i in 1:(length(arr)-1)
        for j in (i+1):length(arr)
            if arr[i] == arr[j]
                return true
            end
        end
    end
    false
end

Base.setdiff(t1::DTag, t2::DTag) = DTag(setdiff(t1.tag, t2.tag) |> sort)
Base.intersect(t1::DTag, t2::DTag) = DTag(intersect(t1.tag, t2.tag) |> sort)

tagRemove(t1::DTag, t2::DTag) = intersect(t1,t2) != DTag() ? setdiff(t1,t2) : DTag(-1)
t1 = DTag([1,2])
t2 = DTag([1])

#---------------------------------------------------------------
#---------------------------------------------------------------
# Differentials
type Differential
    terms::SortedDict
end

function printEpsilons(t::DTag)
    str = ""
    for i in t.tag
        str = str*"ϵ$i"
    end
    str
end

function Base.show(io::IO, diff::Differential)
    str = ""
    for (tag, term) in diff.terms
        if tag == DTag()
            str = str*"$(term)$(printEpsilons(tag))"
        elseif str == ""
            str = str*"($(term))$(printEpsilons(tag))"
        else
            str = str*" + ($(term))$(printEpsilons(tag))"
        end
    end
    print(io, str)
end

function Differential(iterable)
    Differential(try delete!(SortedDict(iterable),DTag(-1)) catch SortedDict(iterable) end)
end

function Differential(keys::Union{Array,Tuple}, values::Union{Array,Tuple})
    Differential(try delete!(SortedDict(zip(keys,values)), DTag(-1)) catch SortedDict(zip(keys,values)) end)
end

Base.length(t::Differential) = length(t.terms)
Base.:(==)(x::Differential,y::Differential) = x.terms == y.terms
Base.getindex(t::Differential, i::DTag) = (t.terms)[i]
getTagList(Dx::Differential) = [key for (key, _) in Dx.terms]

function Base.getindex(Dx::Differential, i::Int)
    key = getTagList(Dx)[i]
    key == DTag() ? Dx.terms[key] : Differential(key => Dx.terms[key])
end

function Base.getindex(Dx::Differential, i::UnitRange)
    keys = getTagList(Dx)[i]
    if length(i) == 1
        out = (Dx.terms)[keys...]
    else
        Differential(k => Dx.terms[k] for k in keys)
    end
end

Base.endof(Dx::Differential) = length(Dx |> getTagList)

lastTag(Dx::Differential) = getTagList(Dx)[end]

function tagUnion(x::Differential, y::Differential)
    vcat(x |> getTagList, y |> getTagList) |> sort |> fastuniq
end

function fastuniq(v)
  v1 = Vector{eltype(v)}()
  if length(v)>0
    laste = v[1]
    push!(v1,laste)
    for e in v
      if e != laste
        laste = e
        push!(v1,laste)
      end
    end
  end
  return v1
end

tagCount = 0

#---------------------------------------------------------------
#---------------------------------------------------------------
# Dtype
struct Dtype <: Function end
const D = Dtype()
const ∂ = Dtype()

Base.:^(D::Dtype, i::Int) = t -> foldl((x,y)->D(x),t,1:i)

max_stars_repo_name: MasonProtter/ModellingToolkitDeriv

max_issues_repo_name: MasonProtter/ModellingToolkitDeriv

max_forks_repo_name: MasonProtter/ModellingToolkitDeriv


text:
@testset "tasks" begin
    c = Condition()
    ch = Channel{String}(1)
    f() = "original"
    function background()
        # Wait until notified allowing us to control when this async code is executed
        wait(c)
        put!(ch, @mock f())
        return nothing
    end

    p = @patch f() = "mocked"

    @sync begin
        # Task started outside patched context should not call patched functions.
        @async background()
        yield()

        apply(p) do
            @test (@mock f()) == "mocked"

            notify(c)
            @test_broken take!(ch) == "original"

            # Task started inside patched context should call patched functions.
            @async background()
            yield()
            notify(c)
            @test take!(ch) == "mocked"

            # Task started inside patched context should call patched functions even when
            # execution finishes outside of patched context.
            @async background()
            yield()
        end

        notify(c)
        @test_broken take!(ch) == "mocked"
    end
end


text:
module ParallelGLM
    using Distributions, Docile, StatsBase
    export GLM,                         # types
           IdentityLink,
           InverseLink,
           Link,
           LogLink,
           LogitLink,
           PGLM,
           SGLM,

           canonical,                   # functions
           devresid2,
           initμη!,
           invlink,
           link,
           updateXtW!,
           varfunc,
           μη

    include("link.jl")
    include("dist.jl")
    include("pglm.jl")
    include("irls.jl")

end # module


text:
using Xiao2021
using Test

@testset "solve_eqm" begin
    # Write your tests here.
end


text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

using DSP, MultiThreadingTools

export preprocess_transform
function preprocess_transform(env::DLEnv, lib::EventLibrary, steps::Vector{String}; copyf=deepcopy)
  initialize(lib)
  result = copyf(lib)
  setname!(result, lib[:name]*"_preprocessed")
  put_label!(result, :FailedPreprocessing, zeros(Float32, eventcount(lib)))
  result.prop[:preprocessing] = steps

  # perform the steps
  for (i,step) in enumerate(steps)
    info(env, 3, "Preprocesing $step on $(eventcount(result)) events...")
    pfunction = eval(parse(step))
    result = pfunction(result)
  end

  failed_count = length(find(f -> f != 0, result[:FailedPreprocessing]))
  if failed_count > 0
    info(env, 1, "Preprocesing failed for $failed_count events in $(lib[:name]). These have been tagged with the label :FailedPreprocessing = 1")
  end

  return result
end

export charge_pulses
function charge_pulses(events::EventLibrary; create_new=true)
  if events[:waveform_type] == "charge"
    return events
  elseif events[:waveform_type] == "current"
    if create_new
      events = deepcopy(events)
    end
    return integrate(events)
  elseif events[:waveform_type] == "raw"
    if create_new
      events = deepcopy(events)
    end
    return raw_to_charge(events)
  else
    throw(ArgumentError("EventLibrary must be of type current or charge."))
  end
end

export current_pulses
function current_pulses(events::EventLibrary; create_new=true)
  if events[:waveform_type] == "current"
    return events
  elseif events[:waveform_type] == "charge"
    if create_new
      events = deepcopy(events)
    end
    return differentiate(events)
  elseif events[:waveform_type] == "raw"
    if create_new
      events = deepcopy(events)
    end
    return current_pulses(raw_to_charge(events); create_new=false)
  else
    throw(ArgumentError("EventLibrary must be of type current or charge."))
  end
end

function current_pulses(libs::Dict{Symbol,EventLibrary}; create_new=true)
  return mapvalues(libs, current_pulses, create_new)
end

function raw_to_charge(events::EventLibrary)
  if events[:waveform_type] != "raw"
    throw(ArgumentError("EventLibrary must contain pulses of type raw"))
  end
  events.waveforms = - events.waveforms
  events.prop[:waveform_type] = "charge"
  return events
end

export baseline
function baseline(events::EventLibrary)
  events = charge_pulses(events; create_new=false)

  bl_size = Int64(round(sample_size(events) / 5))
  weights = hamming(bl_size)
  weights /= sum(weights)

    weights2 = hamming(Int64(bl_size/2))
    weights2 /= sum(weights2)

  bl_levels = zeros(eventcount(events))
  bl_std = zeros(eventcount(events))
  bl_slope = zeros(eventcount(events))

  @everythread for i in threadpartition(1:eventcount(events))
    bl_level = dot(events.waveforms[1:bl_size, i], weights)
    bl_levels[i] = bl_level
    bl_slope[i] = (bl_level - dot(events.waveforms[Int64(bl_size/2)+1:bl_size, i], weights2)) / (bl_size/2)
    events.waveforms[:, i] -= bl_level
    bl_std[i] = std(events.waveforms[1:bl_size, i])
  end

  put_label!(events, :baseline_level, bl_levels)
  put_label!(events, :baseline_std, bl_std)
  put_label!(events, :baseline_slope, bl_slope)
  return events
end

export HE
function HE(events::EventLibrary; cut=1000)
  return filter(events, :E, E -> E>cut)
end

export align_peaks
function align_peaks(events::EventLibrary; target_length=256)
  currents = current_pulses(events; create_new=true)

  s = sample_size(events)
  half = Int64(target_length/2)
  rwf = zeros(Float32, target_length, eventcount(events))

  @everythread for i in threadpartition(1:eventcount(events))
    max_index = findmax(currents.waveforms[:,i])[2]
    if (max_index < half) || (max_index > s - half)
      events[:FailedPreprocessing][i] = 1
    else
      rwf[:,i] = events.waveforms[(max_index-half+1) : (max_index+half) , i]
    end
  end

  events.waveforms = rwf

  return events
end

export align_midpoints
function align_midpoints(events::EventLibrary; center_y=0.5, target_length=256)
  charges = charge_pulses(events; create_new=true)

  s = sample_size(events)
  half = Int64(target_length/2)
  rwf = zeros(Float32, target_length, eventcount(events))

  @everythread for i in threadpartition(1:eventcount(events))
    index = findmin(abs.(charges.waveforms[:,i] - center_y))[2]
    if (index < half) || (index > s - half)
      events[:FailedPreprocessing][i] = 1
    else
      rwf[:,i] = events.waveforms[(index-half+1) : (index+half) , i]
    end
  end

  events.waveforms = rwf

  return events
end

export extract_noise
function extract_noise(events::EventLibrary; target_length=256)
  s = sample_size(events)
  rwf = zeros(Float32, target_length, eventcount(events))

  @everythread for i in threadpartition(1:eventcount(events))
    rwf[:,i] = events.waveforms[(s-target_length+1):s , i]
    rwf[:,i] -= mean(rwf[:,i])
  end

  events.waveforms = rwf

  return events
end

export normalize_energy
function normalize_energy(events::EventLibrary; value=1)
  charges = charge_pulses(events; create_new=true)

  top_size = Int64(round(sample_size(events) / 5))
  weights = hamming(top_size)
  weights /= sum(weights)

  top_levels = zeros(Float32, eventcount(events))

  @everythread for i in threadpartition(1:eventcount(events))
    top_level = dot(charges.waveforms[(end-top_size+1) : end, i], weights)
    top_levels[i] = top_level
    events.waveforms[:,i] *= value / top_level
  end

  put_label!(events, :top_level, top_levels)
  return events
end

export denoise_waveforms!
function denoise_waveforms!(events::EventLibrary; designmethod=Butterworth(5), lowpass_cutoff=0.2)
  waveforms = events[:wf]
  ff = digitalfilter(Lowpass(lowpass_cutoff),designmethod)
  for i in 1:size(waveforms,2)
        waveforms[:,i] = filtfilt(ff, waveforms[:,i])
  end
end

function denoise_waveforms!(data::DLData; designmethod=Butterworth(5), lowpass_cutoff=0.2)
  for lib in data
    denoise_waveforms!(lib)
  end
end

export integrate
function integrate(events::EventLibrary)
    events = initialize(events)
  @everythread for i in threadpartition(1:eventcount(events))
    events.waveforms[:,i] = integrate_array(events.waveforms[:,i])
  end
  if events[:waveform_type] == "charge"
    events.prop[:waveform_type] = "other"
  elseif events[:waveform_type] == "current"
    events.prop[:waveform_type] = "charge"
  end
  return events
end

export integrate_array
function integrate_array(a)
  b = copy(a)
  for i in 2:length(a)
    b[i] += b[i-1]
  end
  return b
end

export differentiate
function differentiate(events::EventLibrary)
    events = initialize(events)
  @everythread for i in threadpartition(1:eventcount(events))
      events.waveforms[:,i] = vcat(0, diff(events.waveforms[:,i]))
  end
  if events[:waveform_type] == "current"
    events.prop[:waveform_type] = "other"
  elseif events[:waveform_type] == "charge"
    events.prop[:waveform_type] = "current"
  end
  return events
end

export scale_waveforms
function scale_waveforms(lib::EventLibrary, val::Number)
  lib = copy(initialize(lib))
  lib.waveforms = lib.waveforms * val
  return lib
end
function scale_waveforms(data::DLData, val::Number)
  return DLData([scale_waveforms(lib, val) for lib in data])
end

type NaNError <: Exception
    msg::AbstractString
end

export check_nan
function check_nan(lib::EventLibrary)
  wf = waveforms(lib)
  nan_indices = find(x->isnan(x), wf)
  if length(nan_indices) > 0
    throw(NaNError("At indices $nan_indices"))
  end
end
check_nan(data::DLData) = for lib in data check_nan(lib) end

export correct_tau_decay
function correct_tau_decay(lib::EventLibrary; tau=52.098e-6/4e-9) # tau in units of sampling interval
  alpha = 1 - exp(- 1/tau)

  events = charge_pulses(lib)
  waveforms = events[:wf]

  @everythread for pulse in threadpartition(1:eventcount(events))
      acc = 0
      for i in 1:size(waveforms,1)
          x = waveforms[i, pulse]
          res = x + acc
          acc += x * alpha
          waveforms[i, pulse] = res
      end
  end

  return events
end

export calculate_deviation
function calculate_deviation(pulses::EventCollection, reconst::EventCollection, noise::EventCollection)
    std_reconst = zeros(Float32, eventcount(noise))
    std_noise = zeros(Float32, eventcount(noise))

    wf_pulses = waveforms(pulses)
    wf_reconst = waveforms(reconst)
    wf_noise = waveforms(noise)

    for i in 1:eventcount(pulses)
        std_reconst[i] = std(wf_reconst[:,i]-wf_pulses[:,i])
        std_noise[i] = std(wf_noise[:,i])
    end

    return std_reconst, std_noise
end


text:
include("iPredict.jl/src/MixtureModelling.jl")
using DelimitedFiles
using Distributions
using StatsPlots
using DataFrames
using Plots, Plots.Measures

scnorm_counts = ["H1_24_scnorm_counts.txt",
"H1_96_scnorm_counts.txt",
"H9_24_scnorm_counts.txt",
"H9_96_scnorm_counts.txt"]
isoChoiceModels = ["Weibull", "Random", "UniformObserved", "CellVariable"]

function make_k_plot(p1,p2,p3,p4,pReal)

    #initialise mu, sigma and relationships
    mu_arr = Float64[p1[2], p2[2], p3[2], p4[2]]
    sigma_arr = Float64[p1[3], p2[3], p3[3], p4[3]]
    k_arr = Float64[0.25, 0.25, 0.25, 0.25]
    relationships = Array{Float64}(undef, length(pReal[4]), 4)

    df = DataFrame(iteration = Int64[0, 0, 0, 0], NumIsos = Int64[1,2,3,4],
    k = Float64[0.25,0.25,0.25,0.25])

    #Do expectation maximisation for k.
    for i in 1:100
        expectation!(relationships, length(pReal[4]), 4, mu_arr, sigma_arr, k_arr, pReal[4])
        maximisation_k!(relationships, k_arr)

        #save results
        newRow = DataFrame(iteration = repeat([i], 4), NumIsos = [1,2,3,4],
        k = k_arr)
        append!(df, newRow)
    end

    l = @layout [a{0.001h}; b c{0.13w}]

    p1 = Plots.plot(df.iteration, df.k, group = df.NumIsos,
    xlabel = "Iteration", ylabel = "Mixing Fraction", legend = :none,
    axis = (font(20)))

    p2 = Plots.plot(df.iteration, df.k, group = df.NumIsos, grid=false,
    xlims=(-1,-0.5), showaxis=false, legendfontsize = 20)

    p0=plot(title ="", grid=false, showaxis=false)

    #make a plot of k vs iterations to see whether we reached convergence
    out_plot = plot(p0,p1,p2,layout=l, size = (250 * 8, 200 * 8),
    linestyle=:solid, linealpha=0.5, linewidth=2*8, axis = (font(60)),
    legendfontsize=60, margin = 20mm)

    return out_plot
end

for countsMatrix in scnorm_counts

    for model in 1:length(isoChoiceModels)

        #Get the data
        data_store = "figures/figure5_data/" * isoChoiceModels[model] * "_" * countsMatrix
        out_plot_path = "figures/figure5/" * isoChoiceModels[model] * "_" * countsMatrix

        #Extract a plot, mu, sigma and the data
        pReal = GetLogNormal(data_store * "_real_isoforms_array.txt", "Real")
        p1 = GetLogNormal(data_store * "_1_isoforms_array.txt", "1 Isoform")
        p2 = GetLogNormal(data_store * "_2_isoforms_array.txt", "2 Isoforms")
        p3 = GetLogNormal(data_store * "_3_isoforms_array.txt", "3 Isoforms")
        p4 = GetLogNormal(data_store * "_4_isoforms_array.txt", "4 Isoforms")

        #make a plot of k vs iterations to see whether we reached convergence
        out_plot = make_k_plot(p1,p2,p3,p4,pReal)
        savefig(out_plot, out_plot_path)

        #See how good the log normal fit is
        ln_plot = plot([p1[1], p2[1], p3[1], p4[1], pReal[1]]...,
        legend = :none, layout=(5,1),  size = (250*8, 600*8),
        axis = (font(180)), titlefontsize=180, margin = 30mm,
 linestyle=:solid, linealpha=0.5, linewidth=2*10,)
        savefig(ln_plot, out_plot_path * "logNormalFits.png")

    end
end

H1_96_1 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_1_isoforms_array.txt",
"1 Isoform")
H1_96_2 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_2_isoforms_array.txt",
"2 Isoforms")
H1_96_3 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_3_isoforms_array.txt",
"3 Isoforms")
H1_96_4 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_4_isoforms_array.txt",
"4 Isoforms")
H1_96_real =  GetLogNormal(
"figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_real_isoforms_array.txt",
"Real")

H1_96_kplot = make_k_plot( H1_96_1, H1_96_2, H1_96_3, H1_96_4, H1_96_real)

H1_24_1 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_1_isoforms_array.txt",
"1 Isoform")
H1_24_2 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_2_isoforms_array.txt",
"2 Isoforms")
H1_24_3 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt",
"3 Isoforms")
H1_24_4 =  GetLogNormal(
"figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_4_isoforms_array.txt",
"4 Isoforms")
H1_24_real =  GetLogNormal(
"figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_real_isoforms_array.txt",
"Real")

H1_24_kplot = make_k_plot( H1_24_1, H1_24_2, H1_24_3, H1_24_4, H1_24_real)

pAB = plot( H1_96_1[1], H1_24_1[1],
              H1_96_2[1], H1_24_2[1],
              H1_96_3[1], H1_24_3[1],
              H1_96_4[1], H1_24_4[1],
              H1_96_real[1], H1_24_real[1],
              #H1_96_kplot, H1_24_kplot,
              layout=(5,2), size = (500*8, 600*8),
              axis = (font(100)), titlefontsize=100, margin = 30mm,
       linestyle=:solid, linealpha=0.5, linewidth=2*10,)

savefig(pAB, "figures/figure5/figure5AB.png")
savefig(H1_96_kplot, "figures/figure5/figure5C.png")
savefig(H1_24_kplot, "figures/figure5/figure5D.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
# TODO: Fix filter to ignore directories.
# TODO: Add handling of errors to allow partial results.

export head, head!

@with_kw struct HeadOptions
    lines = 10
    quiet = false
    verbose = false
end

function parse_options(::Type{HeadOptions}, options)
    # TODO: Allow negative number of lines to read until n from the end.
    max_lines = get_option(
        options, :n, :lines; default = 10).second
    quiet = get_option(
        options, :q, :quiet, :silent, default = false).second
    verbose = get_option(
        options, :v, :verbose, default = false).second
    if verbose && quiet
        error(":verbose and :quiet both supplied. At most one can be.")
    end
    HeadOptions(
        lines = max_lines,
        quiet = quiet,
        verbose = verbose,
    )
end

"""
    head(options..., files...)

Return array of first 10 lines of each file.

With more than one file, return each as pair `file => lines` or print each
with a header giving the file name.
With no file, or when file is "-", read standard input.

Mandatory arguments to long options are mandatory for short options too.

    :n, :lines => [-]NUM    print the first NUM lines instead of the first 10;
                            with the leading '-', print all but the last
                            NUM lines of each file
    :q, :quiet, :silent     never print headers giving file names
    :v, :verbose            always print headers giving file names
      :help     display this help and exit

    FOLLOWING NOT IMPLEMENTED
    :c, :bytes => [-]NUM    print the first NUM bytes of each file;
                            with the leading '-', print all but the last
                            NUM bytes of each file
    :z, :zero-terminated    line delimiter is NUL, not newline
      :version  output version information and exit

NUM may have a multiplier suffix:
b 512, kB 1000, K 1024, MB 1000\\*1000, M 1024\\*1024,
GB 1000\\*1000\\*1000, G 1024\\*1024\\*1024, and so on for T, P, E, Z, Y.

"""
function head(args...)
    options_seq, pos_args = parse_args(args)
    do_print_usage(options_seq, @doc head) && return
    options = parse_options(HeadOptions, options_seq)
    _head(options, pos_args)
end

function _head(options, args)
    multiple_files = length(args) > 1
    max_lines = options.lines
    quiet = options.quiet
    verbose = options.verbose

    lines = []
    if length(args) == 0
        args = ["-"]
    end
    for fn in args
        file_lines = String[]
        n_lines = 0
        fn == "-" && fn == stdin
        for fline = eachline(fn)
            push!(file_lines, fline)
            n_lines += 1
            n_lines ≥ max_lines && break
        end
        if multiple_files
            if quiet
                push!(lines, file_lines)
            else
                push!(lines, fn => file_lines)
            end
        else
            if verbose
                push!(lines, fn => file_lines)
            else
                append!(lines, file_lines)
            end
        end
    end
    lines
end

function head!(args...)
    options, pos_args = parse_args(args)
    do_print_usage(options, @doc head) && return
    options = parse_options(HeadOptions, options)
    lines = _head(options, pos_args)
    lines === nothing && return

    multiple_files = length(pos_args) > 1
    if multiple_files
        for file_lines in lines
            _head_printlines(file_lines)
        end
    else
        if lines[1] isa Pair
            _head_printlines(lines[1])
        else
            _head_printlines(lines)
        end
    end
end

function _head_printlines(file_lines)
    for line in file_lines
        println(line)
    end
    println()
end

function _head_printlines(file_lines::Pair)
    println("==> ", file_lines.first, " <==")
    _head_printlines(file_lines.second)
end


text:
include("CaboCha.jl")

sentence = "クロールで泳いでいる少女を見た"

# 解析結果、Dict型、添え字は1から
@show cabocha_parser(sentence)

# 標準出力のみ(第二引数指定) #
println(cabocha_parser(sentence, "f0"))
println(cabocha_parser(sentence, "f1"))
println(cabocha_parser(sentence, "f3"))


text:
using CSV
using DataFrames
using Markdown
using Printf

julia_results = CSV.read("bench_results/GaussianProcesses_jl.csv", header=["label", "GaussianProcesses.jl"])
GPy_results = CSV.read("bench_results/GPy.csv", header=["label", "GPy"])
gpml_results = CSV.read("bench_results/gpml.csv", header=["label", "gpml"])

comparison = join(julia_results, GPy_results, kind=:outer, on=:label, makeunique=false)
comparison = join(comparison, gpml_results,   kind=:outer, on=:label, makeunique=false)
sort!(comparison, Symbol("GaussianProcesses.jl"))

CSV.write("bench_results/compare.csv", comparison);
comparison

""" Return the column index of the fastest execution for row i """
function rowindmin(i)
    row = [comparison[i, col] for col in 2:4]
    row_nomiss = [(ismissing(x) ? Inf : x) for x in row]
    imin = argmin(row_nomiss)
    return imin +1
end

imin = rowindmin.(1:nrow(comparison))
print(imin)

fmt(::Missing) = ""
fmt(x::Int64) = @sprintf("%d", x)
fmt(x::Float64) = @sprintf("%.0f", x)
fmt(x) = String(x)

rows = [String.(names(comparison))] # start with just header

for irow in 1:nrow(comparison)
    row = [fmt(comparison[irow, col]) for col in names(comparison)]
    row[imin[irow]] = @sprintf("**%s**", row[imin[irow]])
    push!(rows, row)
end

t=Markdown.Table(rows, [:r, :r, :r, :r, :r])

b = IOBuffer()
Markdown.plain(b, t)
print(String(take!(b)))

fmt(::Missing) = ""
fmt(x::Int64) = @sprintf("%d", x)
fmt(x::Float64) = @sprintf("%.0f", x)
fmt(x) = String(x)

rows = [String.(names(comparison))] # start with just header

for irow in 1:nrow(comparison)
    row = [fmt(comparison[irow, col]) for col in names(comparison)]
    row[imin[irow]] = @sprintf("**%s**", row[imin[irow]])
    push!(rows, row)
end

t=Markdown.Table(rows, [:r, :r, :r, :r, :r])

b = IOBuffer()
Markdown.plain(b, t)
print(String(take!(b)))


text:
β = [0.2, -1.0]
γ = [-0.4, 0.5]

θtrue = [β..., γ...]

function θstata(model::Symbol)
    model == :probit && return [.20154344,  -.95584277,  -.36636477,    .4843079]
    model == :logit  && return [.33638928, -1.6251676, -.62388885, .82423386]
    throw(error())
end

function vcovstata(model::Symbol)

    model == :probit && return [ .00159179  -.00028434   -.0000851   .00003368;
                                -.00028434   .00252983   .00032415  -.00026669;
                                 -.0000851   .00032415   .00210196   .00095648;
                                 .00003368  -.00026669   .00095648   .00218769]
    model == :logit && return [.00470505  -.00105662  -.00032196   .00017456;
                              -.00105662   .00868674   .00130306  -.00135786;
                              -.00032196   .00130306   .00606232   .00258262;
                               .00017456  -.00135786   .00258262   .00647028]
    throw(error())
end

function θpolr(model::Symbol)
    model == :logit  && return [0x1.58766ddc938dfp-2, -0x1.a00afbb34ab94p+0, -0x1.3f6e5c3a0abc6p-1, 0x1.a601fa3d8b336p-1]
    model == :probit && return [0x1.9cc2cf36aaba6p-3, -0x1.e9643962405bap-1, -0x1.772853f0930d9p-2, 0x1.efee6817f6849p-2 ]
    throw(error())
end

function vcovpolr(model::Symbol)
    model == :logit  && return reshape([0x1.3459b8009a617p-8, -0x1.14fcb4689d42cp-10, -0x1.5199a6ddeb182p-12, 0x1.6e13a5ae9d64cp-13, -0x1.14fcb4689d428p-10, 0x1.1ca59d7d22ea6p-7, 0x1.55971fae05968p-10, -0x1.63f420acf772ap-10, -0x1.5199a6ddeb18cp-12, 0x1.55971fae0598p-10, 0x1.8d4cd2d9c3a5cp-8, 0x1.52825fc23e733p-9, 0x1.6e13a5ae9d678p-13, -0x1.63f420acf771cp-10, 0x1.52825fc23e738p-9, 0x1.a8094dae89bfep-8], 4, 4)
    model == :probit && return reshape([0x1.a1478d41c4cc5p-10, -0x1.2a27562a3841ap-12, -0x1.64f3df22abaacp-14, 0x1.1a8ad6e51a082p-15, -0x1.2a27562a3841p-12, 0x1.4b970d34ca5eap-9, 0x1.53e4445a38accp-12, -0x1.17a49b5a7ce5cp-12, -0x1.64f3df22abaafp-14, 0x1.53e4445a38ad6p-12, 0x1.1381eb141b246p-9, 0x1.f579486d2078p-11, 0x1.1a8ad6e51a068p-15, -0x1.17a49b5a7ce5ep-12, 0x1.f579486d2078ep-11, 0x1.1ebe907559363p-9], 4, 4)
    throw(error())
end


text:
using PDDL, PDDL.Parser, Test

# Define equivalence shorthand for abstract interpreter testing
≃(a, b) = PDDL.equiv(a, b)

include("strips/test.jl")
include("typing/test.jl")
include("axioms/test.jl")
include("adl/test.jl")
include("numeric/test.jl")
include("constants/test.jl")
include("arrays/test.jl")
include("sets/test.jl")
include("functions/test.jl")


text:
# BFS demo using abhinavmehndiratta SuiteSparseGraphBLAS Julia wrapper
# (can be found at https://github.com/abhinavmehndiratta/SuiteSparseGraphBLAS.jl)

using sparse, GraphBLASInterface, SuiteSparseGraphBLAS, BenchmarkTools

# Begin GraphBLAS session
GrB_init(GrB_NONBLOCKING)

# GraphBLAS matrix and vector dimensions
dim = 100

# Build a random sparse GraphBLAS matrix
A = GrB_Matrix(Int64, dim, dim)
GraphBLAS_sparse_square_matrix!(A, dim)

# Create an empty GraphBLAS integer vector
v = GrB_Vector(Int64, dim)

bench = @benchmark GraphBLAS_bfs!(A, v, ZeroBasedIndex(0))

dropzeros!(v)
GrB_wait()

# Free workspace
GrB_Matrix_free(A)
GrB_Vector_free(v)

# End GraphBLAS session
GrB_finalize()

display(bench)