text:
export iauC2ibpn
"""
Form the celestial-to-intermediate matrix for a given date given
the bias-precession-nutation matrix.  IAU 2000.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2 double       TT as a 2-part Julian Date (Note 1)
   rbpn        double[3][3] celestial-to-true matrix (Note 2)

Returned:
   rc2i        double[3][3] celestial-to-intermediate matrix (Note 3)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix rbpn transforms vectors from GCRS to true equator (and
   CIO or equinox) of date.  Only the CIP (bottom row) is used.

   3. The matrix rc2i is the first stage in the transformation from
   celestial to terrestrial coordinates:

      [TRS] = RPOM * R_3(ERA) * rc2i * [CRS]

            = RC2T * [CRS]

   where [CRS] is a vector in the Geocentric Celestial Reference
   System and [TRS] is a vector in the International Terrestrial
   Reference System (see IERS Conventions 2003), ERA is the Earth
   Rotation Angle and RPOM is the polar motion matrix.

   4. Although its name does not include "00", This function is in fact
   specific to the IAU 2000 models.

Called:
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauC2ixy     celestial-to-intermediate matrix, given X,Y

References:
   "Expressions for the Celestial Intermediate Pole and Celestial
   Ephemeris Origin consistent with the IAU 2000A precession-
   nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

   n.b. The celestial ephemeris origin (CEO) was renamed "celestial
         intermediate origin" (CIO) by IAU 2006 Resolution 2.

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauC2ibpn(double date1, double date2, double rbpn[3][3],
#                double rc2i[3][3])

function iauC2ibpn(date1::Real, date2::Real, rbpn::AbstractMatrix{<:Real})

   # Allocate return value
   rc2i = zeros(Float64, 3, 3)

   ccall((:iauC2ibpn, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Ptr{Cdouble},
         Ptr{Cdouble}),
         convert(Float64, date1), convert(Float64, date2), convert(Array{Float64, 2}, rbpn'),
         rc2i)

   return SMatrix{3,3}(rc2i')
end


text:
# Set ARCH dependant types / constants
# Directly copied or inspired from https://github.com/JuliaParallel/PETSc.jl.git

"""
    Retrieve a PETSc datatype from a String
"""
function DataTypeFromString(name::AbstractString)
    dtype_ref = Ref{PetscDataType}()
    found_ref = Ref{PetscBool}()
    ccall((:PetscDataTypeFromString, libpetsc), PetscErrorCode,
             (Cstring, Ptr{PetscDataType}, Ptr{PetscBool}),
             name, dtype_ref, found_ref)
    @assert found_ref[] == PETSC_TRUE
    return dtype_ref[]
end

"""
    Retrieve a PETSc datatype from a PETSc datatype
"""
function PetscDataTypeGetSize(dtype::PetscDataType)
    datasize_ref = Ref{Csize_t}()
    ccall((:PetscDataTypeGetSize, libpetsc), PetscErrorCode,
             (PetscDataType, Ptr{Csize_t}),
             dtype, datasize_ref)
    return datasize_ref[]
end

"""
    Find the Julia type for `PetscReal`
"""
function PetscReal2Type()
    # Workaround for RegistryCI
    (libpetsc == "JULIA_REGISTRYCI_AUTOMERGE") && (return Cdouble)

    PETSC_REAL = DataTypeFromString("Real")
    result =
        PETSC_REAL == PETSC_DOUBLE ? Cdouble :
        PETSC_REAL == PETSC_FLOAT ? Cfloat :
        error("PETSC_REAL = $PETSC_REAL not supported.")
    return result
end

"""
    Find the Julia type for PetscScalar. This function must be called after
    the line `const PetscReal = PetscReal2Type` since it uses `PetscReal`.
"""
function PetscScalar2Type()
    # Workaround for RegistryCI
    (libpetsc == "JULIA_REGISTRYCI_AUTOMERGE") && (return Cdouble)

    PETSC_REAL = DataTypeFromString("Real")
    PETSC_SCALAR = DataTypeFromString("Scalar")
    result =
        PETSC_SCALAR == PETSC_REAL ? PetscReal :
        PETSC_SCALAR == PETSC_COMPLEX ? Complex{PetscReal} :
        error("PETSC_SCALAR = $PETSC_SCALAR not supported.")
    return result
end

"""
    Find the Julia type for `PetscInt`
"""
function PetscInt2Type()
    # Workaround for RegistryCI
    (libpetsc == "JULIA_REGISTRYCI_AUTOMERGE") && (return Int32)

    PETSC_INT_SIZE = PetscDataTypeGetSize(PETSC_INT)
    result =
        PETSC_INT_SIZE == 4 ? Int32 :
        PETSC_INT_SIZE == 8 ? Int64 :
        error("PETSC_INT_SIZE = $PETSC_INT_SIZE not supported.")
    return result
end

const PetscReal = PetscReal2Type()
const PetscScalar = PetscScalar2Type()
const PetscInt = PetscInt2Type()

const PetscIntOne = PetscInt(1) # Integer `1` with the type of PetscInt, usefull to go back and forth with julia/petsc indexing

ext: jl


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# usage: plot_lp_results.py [results directory]

using CSV
using DataFrames
using Plots

# If set to false, all restart frequencies are plotted.
PLOT_ONLY_BEST_FIXED_RESTART = true

gr()

fntsm = Plots.font("serif", pointsize=12)
fntlg = Plots.font("serif", pointsize=18)
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)

pdhg = DataFrame(CSV.File(joinpath(ARGS[1],"pdhg.csv")))
pdhg_adaptive = DataFrame(CSV.File(joinpath(ARGS[1],"pdhg_adaptive.csv")))
pdhg_restart = Dict()

last_active_set_update = findall(pdhg_adaptive[!, "num_active_bounds_changed"] .> 0)[end]
@show last_active_set_update

# This list should match up with lp_experiment.py.
restart_frequencies = [64, 256, 1024, 4096, 16384, 65536]
for restart_frequency in restart_frequencies
  filename = "pdhg_restart_$(restart_frequency).csv"
  pdhg_restart[restart_frequency] = DataFrame(CSV.File(joinpath(ARGS[1], filename)))
end

iter_num = pdhg[!, "iteration_num"]
feature_names = [("current_kkt_err_l2", "L2 norm of KKT residual (current iterate)"),
                 ("current_kkt_err_linf", "LInf norm of KKT residual (current iterate)"),
                 ("avg_kkt_err_l2", "L2 norm of KKT residual (averaged iterate)"),
                 ("avg_kkt_err_linf", "LInf norm of KKT residual (averaged iterate)")]

# Find the fixed restart frequency for which current_kkt_err_l2 feature first
# drops below 1e-6.
best_restart_frequency = nothing
earliest_time = typemax(Int32)
for restart_frequency in restart_frequencies
  indices = findall(pdhg_restart[restart_frequency][!, "current_kkt_err_l2"] .<= 1e-6)
  if length(indices) > 0 && indices[1] < earliest_time
    global earliest_time = indices[1]
    global best_restart_frequency = restart_frequency
  end
end

for (feature, name) in feature_names
  plot(iter_num, pdhg[!, feature], label="pdhg", line=(1,:solid))
  plot!(iter_num, pdhg_adaptive[!, feature], label="pdhg adaptive",
    line = (1,:solid))

  if PLOT_ONLY_BEST_FIXED_RESTART
    @assert best_restart_frequency !== nothing "No best fixed restart frequency"
    plot!(iter_num, pdhg_restart[best_restart_frequency][!, feature],
          label="best fixed restart ($(best_restart_frequency))",
          line=(1,:solid))
  else
    for restart_frequency in restart_frequencies
      plot!(iter_num, pdhg_restart[restart_frequency][!, feature],
            label="pdhg restart $(restart_frequency)",
            line=(1,:solid))
    end
  end

  scatter!([last_active_set_update], [pdhg_adaptive[last_active_set_update, feature]], label="last active set update", markershape=:star6, markersize=6)

  yaxis!(name, :log10)
  xaxis!("Iteration")
  y_pow_10_range = -9:2
  ylims!((10.0^y_pow_10_range[1], 10.0^y_pow_10_range[end]))
  yticks!(10.0.^y_pow_10_range[1:2:end], ["10^$i" for i in y_pow_10_range[1:2:end]])
  savefig(joinpath(ARGS[1], "$(feature).pdf"))
end

# A custom summary plot with choice of average/current for different methods.
plot(iter_num, pdhg[!, "current_kkt_err_l2"], label="No restarts", line=(1,:solid))
plot!(iter_num, pdhg_adaptive[!, "avg_kkt_err_l2"], label="Our adaptive",
    line = (1,:solid))
plot!(iter_num, pdhg_restart[best_restart_frequency][!, "avg_kkt_err_l2"],
          label="Best fixed ($(best_restart_frequency))",
          line=(1,:solid))
did_restart_indices = findall(pdhg_adaptive[!, "did_restart"] .== 1)
scatter!(did_restart_indices, pdhg_adaptive[did_restart_indices, "avg_kkt_err_l2"], label="", color=2)
did_restart_indices = findall(pdhg_restart[best_restart_frequency][!, "did_restart"] .== 1)
scatter!(did_restart_indices, pdhg_restart[best_restart_frequency][did_restart_indices, "avg_kkt_err_l2"], label="", color=3)
scatter!([last_active_set_update], [pdhg_adaptive[last_active_set_update, "avg_kkt_err_l2"]], label="Last active set change", markershape=:star6, markersize=6)
yaxis!("Residual", :log10)
xaxis!("Iteration")
y_pow_10_range = -9:2
ylims!((10.0^y_pow_10_range[1], 10.0^y_pow_10_range[end]))
yticks!(10.0.^y_pow_10_range[1:2:end], ["10^$i" for i in y_pow_10_range[1:2:end]])
savefig(joinpath(ARGS[1], "summary.pdf"))


text:
mutable struct MainWindow <: GtkWindow

    handle::Ptr{Gtk.GObject}
    style_and_language_manager::StyleAndLanguageManager
    editor#TODO type this ? (circular def)
    console_manager
    pathCBox
    statusBar
    project
    menubar
    sidepanel_ntbook

    function MainWindow()

        w = GtkWindow("GtkIDE.jl - v$(VERSION)", 1800, 1200)
        signal_connect(main_window_key_press_cb, w, "key-press-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)
        signal_connect(main_window_quit_cb, w, "delete-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)

        sl_mng = StyleAndLanguageManager()
        n = new(w.handle, sl_mng)
        Gtk.gobject_move_ref(n, w)
    end
end

function init!(main_window::MainWindow, editor, c_mng, pathCBox, statusBar, project,
    menubar, sidepanel_ntbook)

    main_window.editor = editor
    main_window.console_manager = c_mng
    main_window.pathCBox = pathCBox
    main_window.statusBar = statusBar
    main_window.project = project
    main_window.menubar = menubar
    main_window.sidepanel_ntbook = sidepanel_ntbook
    load(main_window) #load last session
end

## exiting
@guarded (PROPAGATE) function main_window_quit_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    main_window = convert(GtkWindow, widgetptr)

    if typeof(main_window.project) == Project
        save(main_window.project)
    end
    save(main_window)

    global is_running = false
    #REDIRECT_STDOUT && stop_console_redirect(main_window)

    return PROPAGATE
end

function toggle_editor()#use visible ?
    if Gtk.G_.position(mainPan) > 0
        mainPanPos = Gtk.G_.position(mainPan)
        Gtk.G_.position(mainPan, 0)
    else
        Gtk.G_.position(mainPan, 650) #FIXME need a layout type to save all these things
    end
end

@guarded (PROPAGATE) function main_window_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    event = convert(Gtk.GdkEvent, eventptr)
    main_window = convert(GtkWindow, widgetptr)

    mod = get_default_mod_mask()

    if doing(Action("r", PrimaryModifier), event)
        #@async begin
            #crashes if we are still in the callback
            #sleep(0.2)
#            eval(Main, :(restart()))
            #restart(main_window)
        #end
    end
    if event.keyval == GdkKeySyms.F1
        toggle_sidepanel()
    end
    if event.keyval == GdkKeySyms.F2
        toggle_editor()
    end

    if doing(Actions["console_editor_switch"], event)
        c = current_console(main_window.console_manager).view
        if !get_gtk_property(c, :has_focus, Bool)
            grab_focus(c)
        else
            grab_focus(current_tab(main_window.editor).view)
        end
    end

    return PROPAGATE
end

function sidePanelButton_clicked_cb(widgetptr::Ptr, user_data)
    toggle_sidepanel()
    return nothing
end

function editorButtonclicked_cb(widgetptr::Ptr, user_data)
    toggle_editor()
    return nothing
end

function toggle_sidepanel()
    visibility = !visible(sidepanel_ntbook)
    visible(sidepanel_ntbook, visibility)
    for panel in sidepanel_ntbook.panels
        visible(panel, visibility)
        visible(panel) && update!(panel)#Panels can opt out of updating while hidden, so update them when revealed
    end
end

# Not ideal, it always refresh when using the pathdisplay
@guarded nothing function on_path_change(main_window::MainWindow,
    doUpdate=false, console=current_console(main_window))

    c_path = unsafe_string(Gtk.G_.active_text(main_window.pathCBox))
    path = pwd(console)

    update_pathEntry(main_window, path)

    if path != c_path || doUpdate
        push!(main_window.pathCBox, path)
        for panel in main_window.sidepanel_ntbook.panels
            on_path_change(panel, path)
        end
        save(main_window.project)
    end
    nothing
end

function reload()

    GtkREPL.reload()
    Core.eval(GtkIDE, quote
    include(joinpath(HOMEDIR, "gtk_utils.jl"))
    include(joinpath(HOMEDIR, "PlotWindow.jl"))
    include(joinpath(HOMEDIR, "StyleAndLanguageManager.jl"))
    include(joinpath(HOMEDIR, "MainWindow.jl"))
    include(joinpath(HOMEDIR, "Project.jl"))
    include(joinpath(HOMEDIR, "Console.jl"))
    include(joinpath(HOMEDIR, "ConsoleCommands.jl"))
    include(joinpath(HOMEDIR, "Editor.jl"))
    include(joinpath(HOMEDIR, "NtbookUtils.jl"))
    include(joinpath(HOMEDIR, "PathDisplay.jl"))
    include(joinpath(HOMEDIR, "MainMenu.jl"))
    include(joinpath(HOMEDIR, "SidePanels.jl"))
    include(joinpath(HOMEDIR, "Logo.jl"))
    init_console_commands()
    nothing
    end)

end

##
function restart(main_window::MainWindow, new_workspace=false)

        println("restarting...")
        sleep(0.1)
        is_running = false

        update!(main_window.project)
        save(main_window.project)

        REDIRECT_STDOUT && stop_console_redirect(main_window)

        new_workspace && workspace()
        destroy(main_window)
#        gtkide()

        #include( joinpath(HOMEDIR, "GtkIDE.jl") )

        reload()

        #include("init.jl")
        __init__()

end

function run_tests()
    include( joinpath(HOMEDIR, "..", "test", "runtests.jl") )
end

#this allows to save some info about the session
JSON.lower(w::MainWindow) = Dict(
    "project.name" => w.project.name,
 )

function save(w::MainWindow)
    !isdir( joinpath(HOMEDIR, "config") ) && mkdir( joinpath(HOMEDIR, "config") )
    open( joinpath(HOMEDIR, "config", "main_window.json") , "w") do io
        JSON.print(io, w)
    end
end

function load(w::MainWindow)
    !isdir( joinpath(HOMEDIR, "config") ) && mkdir( joinpath(HOMEDIR, "config") )
    pth = joinpath(HOMEDIR, "config", "main_window.json")
    if !isfile(pth)
        w.project.name = "default"
        return
    end

    j = JSON.parsefile(pth)

    w.project.name = j["project.name"]
end


text:
"""
    new_variables_to_attach(formulater::DisjunctiveFormulater)::Vector

Returns a vector of the kinds (from `_V_INT_SETS`) of the variables added
_every time_ that `formulater` is applied to a model. The contract is: these
variables will always be included in a model, regardless of which node it is
created at. Variable creation will be handled by `populate_lp_model!`; it is
not the role of `formulater` to create them.

This function returns the CVIs used to reference these variables in the model.

Note that the new variables need not all be integer; extended formulations can
register additional continuous variables (though this is not necessary--
formulaters can add additional continuous variables not registered with the
formulation in this manner, provided they handle basis getting/setting
appropriately). But any integer variables that you wish to branch on _must_ be
 registered in this manner.
"""
function new_variables_to_attach end

"""
    compute_disjunction_activity(
        form::DMIPFormulation,
        formulater::DisjunctiveFormulater{DisjunctiveConstraints.NaiveBigM},
        node::Node,
        config::AlgorithmConfig,
    )::Vector{Bool}

Computes the activity of a disjunctive constraint: That is, which alternatives
are still feasible for the current `node`.

Returns a vector `activity` of Bools. The entry `activity[i]` is `true` if the
`i`-th disjunction in `formulater.disjunction` is feasible, and false
otherwise.

Notes:
* If `activity` is all falses, the caller can infer that the entire problem is
    infeasible.
* Here, an alternative being "feasible" means that the associated values for
    the integer variables are feasible. Therefore, an alternative may be
    "infeasible" and yet still contained in the feasible region for the problem
    (if, for example, the alternative is completely contained in another).
"""
function compute_disjunction_activity end

"""
    delete_all_constraints!(model::Gurobi.Optimizer, disjunction_state)

Delete all constraints in `model` that were added to formulate a given
disjunctive constraint.
"""
function delete_all_constraints! end

function formulate!(
    state::CurrentState,
    form::DMIPFormulation,
    formulater::DisjunctiveFormulater,
    node::Node,
    node_result::NodeResult,
    config::AlgorithmConfig,
)
    cvis = form.disjunction_formulaters[formulater]
    # TODO: If activity has exactly one true, can model that alternative
    # directly. If activity is all falses, the problem is infeasible, and we
    # can "short circuit" and not worry about formulating or solving.
    activity = compute_disjunction_activity(form, formulater, node, config)
    _f = [instantiate(v, state) for v in formulater.disjunction.f]
    f = MOIU.vectorize(_f)
    # NOTE: Could be a bit more clever here, by dropping alternatives which are
    # proven not to be active. However, this may require us to slice into
    # `cvis` below in the call to `DisjunctiveConstraints.formulate!`.
    # TODO: Do what's described above.
    masked_lbs = copy(formulater.disjunction.s.lbs)
    masked_ubs = copy(formulater.disjunction.s.ubs)
    masked_lbs[:, (!).(activity)] .= -Inf
    masked_ubs[:, (!).(activity)] .= Inf
    s = DisjunctiveConstraints.DisjunctiveSet(masked_lbs, masked_ubs)

    disj_state = DisjunctiveConstraints.formulate!(
        state.gurobi_model,
        formulater.method,
        DisjunctiveConstraints.Disjunction(f, s),
        [instantiate(cvi, state) for cvi in cvis],
    )
    state.disjunction_state[formulater] = disj_state
    if disj_state.proven_infeasible
        node_result.status = INFEASIBLE_LP
        node_result.cost = Inf
    end
    return nothing
end


text:
function findBoundaryIndex(boundaries, name)
    for boundary in boundaries
        boundary.name == name && return boundary.index
    end
end

function assignBC(boundaryName, field, type, value=0)
    iBoundary = findBoundaryIndex(mesh.boundaries, boundaryName)
    boundary = mesh.boundaries[iBoundary]
    boundary.bcType = type

    iCellStart = mesh.nCells + boundary.startFace - mesh.nInternalFaces
    iCellEnd = iCellStart + boundary.nFaces - 1

    iFaceStart = boundary.startFace
    iFaceEnd = iFaceStart + boundary.nFaces - 1

    if type == "fixedValue"
        field.cellValues[iCellStart:iCellEnd] = [value for _ in 1:boundary.nFaces]
    elseif type == "zeroGradient"
        owners = [mesh.faces[iFace].iOwner for iFace in iFaceStart:iFaceEnd]
        field.cellValues[iCellStart:iCellEnd] = [field.cellValues[iOwner] for iOwner in owners]
    end
    field.faceValues[iFaceStart:iFaceEnd] = field.cellValues[iCellStart:iCellEnd]
end


text:
# double factorial
factorial2(n) = prod(n:-2:1)

dist2(dx,dy,dz) = dx*dx + dy*dy + dz*dz

triangle(i::Int64) = div(i*(i+1),2)
triangle(i::Int64,j::Int64) = i < j ? triangle(j-1)+i : triangle(i-1)+j

iindex(i::Int64,j::Int64,k::Int64,l::Int64) = triangle(triangle(i,j),triangle(k,l))

trace2(A,B) = sum(A.*B)


text:
#! /usr/bin/env julia

using ArgParse

function args_main(s)
  @add_arg_table s begin
    "apply"
    help = "apply a provided model on images"
    action = :command
    "train"
    help = "train a model with provided images with annotations"
    action = :command
    "lesson"
    help = "train a model according to a given lesson file (.dcct)"
    action = :command
    "test"
    help = "test the performance of a provided model on images with annotations"
    action = :command
    "synth"
    help = "use DCellC on data constructed by a synthetic module"
    action = :command
  end
end

include("dcellc-apply.jl")
include("dcellc-train.jl")
include("dcellc-lesson.jl")

function cmd_main(args)
  if args["%COMMAND%"] == "apply"
    cmd_apply(args["apply"])
  elseif args["%COMMAND%"] == "train"
    cmd_train(args["train"])
  elseif args["%COMMAND%"] == "lesson"
    cmd_lesson(args["lesson"])
  elseif args["%COMMAND%"] == "test"
    cmd_test(args["test"])
  elseif args["%COMMAND%"] == "synth"
    cmd_synth(args["synth"])
  end
end

s = ArgParseSettings()

args_main(s)

args_apply(s["apply"])
args_train(s["train"])
args_lesson(s["lesson"])
#args_test(s)
#args_synth(s)

# apply
# train
# test
# synth
#   - train
#   - test
#   - generate

using Knet
include("../src/DCellC.jl")
using DCellC

cmd_main(parse_args(s))

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
#
# Common functions for ReferenceSystem
#
function nndistance(ref::ReferenceSystem, g::Array{Function,2})::Array{Float64,2}
  N::Int = ncomp(ref)
  σ::Array{Float64,2} = hsdiameter(ref)
  d = Array{Float64,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    opt = Optim.optimize(r -> -g[i,j](r), 0.5σ[i,j], 1.5σ[i,j])
    d[i,j] = d[j,i] = Optim.minimizer(opt)
  end

  return d
end

function nndistance(ref::ReferenceSystem)::Array{Float64,2}
  g::Array{Function,2} = paircorrelation(ref)
  return nndistance(ref, g)
end

# pairwise contribution to internal energy
function internal(ref::ReferenceSystem, pert::Perturbation, u::Matrix{Pairpotential})::Float64
  N::Int = ncomp(ref)
  ρ::Float64 = totalnumberdensity(ref)
  c::Vector{Float64} = composition(ref)
  g::Matrix{Function} = paircorrelation(ref)
  rcore::Matrix{Float64} = emptyradius(ref)
  rcut::Array{Float64,2} = cutoffradius(pert)

  U::Float64 = 0

  for i in 1:N, j in 1:N
    i > j && continue
    n = i == j ? 1 : 2

    gs = spline(g[i,j], rcore[i,j], rcut[i,j], 32)
    us = spline(u[i,j], rcore[i,j], rcut[i,j], 32)
    U += 2π*ρ * n*c[i]*c[j] * ∫(r -> us(r)*gs(r)*r^2, rcore[i,j], rcut[i,j])
  end

  return U
end

function internal(ref::ReferenceSystem, pert::Perturbation)::Float64
  internal(ref, pert, pairpotential(pert))
end


text:
# TODO: add more target distributions and make them iteratable
# TODO: Integrate with https://github.com/xukai92/VecTargets.jl to achieve goal noted
#       above.

# Dimension of testing distribution
const D = 5
# Tolerance ratio
const TRATIO = Int == Int64 ? 1 : 2
# Deterministic tolerance
const DETATOL = 1e-3 * D * TRATIO
# Random tolerance
const RNDATOL = 5e-2 * D * TRATIO * 2

# Hand-coded multivariate Gaussian

const gaussian_m = zeros(D)
const gaussian_s = ones(D)

struct Gaussian{Tm, Ts}
    m::Tm
    s::Ts
end

function ℓπ_gaussian(g::AbstractVecOrMat{T}, s) where {T}
    return .-(log(2 * T(pi)) .+ 2 .* log.(s) .+ abs2.(g) ./ s.^2) ./ 2
end

ℓπ_gaussian(m, s, x) = ℓπ_gaussian(m .- x, s)

function ∇ℓπ_gaussianl(m, s, x)
    g = m .- x
    v = ℓπ_gaussian(g, s)
    return v, g
end

function get_ℓπ(g::Gaussian)
    ℓπ(x::AbstractVector) = sum(ℓπ_gaussian(g.m, g.s, x))
    ℓπ(x::AbstractMatrix) = dropdims(sum(ℓπ_gaussian(g.m, g.s, x); dims=1); dims=1)
    return ℓπ
end

function get_∇ℓπ(g::Gaussian)
    function ∇ℓπ(x::AbstractVector)
        val, grad = ∇ℓπ_gaussianl(g.m, g.s, x)
        return sum(val), grad
    end
    function ∇ℓπ(x::AbstractMatrix)
        val, grad = ∇ℓπ_gaussianl(g.m, g.s, x)
        return dropdims(sum(val; dims=1); dims=1), grad
    end
    return ∇ℓπ
end

ℓπ = get_ℓπ(Gaussian(gaussian_m, gaussian_s))
∂ℓπ∂θ = get_∇ℓπ(Gaussian(gaussian_m, gaussian_s))

# For the Turing model
# @model gdemo() = begin
#     s ~ InverseGamma(2, 3)
#     m ~ Normal(0, sqrt(s))
#     1.5 ~ Normal(m, sqrt(s))
#     2.0 ~ Normal(m, sqrt(s))
#     return s, m
# end

using Distributions: logpdf, InverseGamma, Normal
using Bijectors: invlink, logpdf_with_trans

function invlink_gdemo(θ)
    s = invlink(InverseGamma(2, 3), θ[1])
    m = θ[2]
    return [s, m]
end

function ℓπ_gdemo(θ)
    s, m = invlink_gdemo(θ)
    logprior = logpdf_with_trans(InverseGamma(2, 3), s, true) + logpdf(Normal(0, sqrt(s)), m)
    loglikelihood = logpdf(Normal(m, sqrt(s)), 1.5) + logpdf(Normal(m, sqrt(s)), 2.0)
    return logprior + loglikelihood
end

test_show(x) = test_show(s -> length(s) > 0, x)
function test_show(pred, x)
    io = IOBuffer(; append = true)
    show(io, x)
    s = read(io, String)
    @test pred(s)
end


text:
#!/usr/bin/env julia --project

using ArgParse
using LinearAlgebra
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.TurbulenceClosures
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics:
    air_density,
    air_temperature,
    total_energy,
    internal_energy,
    PhasePartition,
    air_pressure
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet
using CLIMAParameters.Atmos.SubgridScale

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# from Thatcher and Jablonowski (2016)
CLIMAParameters.Planet.press_triple(::EarthParameterSet) = 610.78

# driver-specific parameters added here
T_sfc_pole(::EarthParameterSet) = 271.0

function init_baroclinic_wave!(problem, bl, state, aux, coords, t)
    FT = eltype(state)

    # parameters
    _grav::FT = grav(bl.param_set)
    _R_d::FT = R_d(bl.param_set)
    _Ω::FT = Omega(bl.param_set)
    _a::FT = planet_radius(bl.param_set)
    _p_0::FT = MSLP(bl.param_set)

    k::FT = 3
    T_E::FT = 310
    T_P::FT = 240
    T_0::FT = 0.5 * (T_E + T_P)
    Γ::FT = 0.005
    A::FT = 1 / Γ
    B::FT = (T_0 - T_P) / T_0 / T_P
    C::FT = 0.5 * (k + 2) * (T_E - T_P) / T_E / T_P
    b::FT = 2
    H::FT = _R_d * T_0 / _grav
    z_t::FT = 15e3
    λ_c::FT = π / 9
    φ_c::FT = 2 * π / 9
    d_0::FT = _a / 6
    V_p::FT = 1
    M_v::FT = 0.608
    p_w::FT = 34e3             ## Pressure width parameter for specific humidity
    η_crit::FT = p_w / _p_0 ## Critical pressure coordinate
    q_0::FT = 1e-4                ## Maximum specific humidity (default: 0.018)
    q_t::FT = 1e-12            ## Specific humidity above artificial tropopause
    φ_w::FT = 2π / 9           ## Specific humidity latitude wind parameter

    # grid
    φ = latitude(bl.orientation, aux)
    λ = longitude(bl.orientation, aux)
    z = altitude(bl.orientation, bl.param_set, aux)
    r::FT = z + _a
    γ::FT = 1 # set to 0 for shallow-atmosphere case and to 1 for deep atmosphere case

    # convenience functions for temperature and pressure
    τ_z_1::FT = exp(Γ * z / T_0)
    τ_z_2::FT = 1 - 2 * (z / b / H)^2
    τ_z_3::FT = exp(-(z / b / H)^2)
    τ_1::FT = 1 / T_0 * τ_z_1 + B * τ_z_2 * τ_z_3
    τ_2::FT = C * τ_z_2 * τ_z_3
    τ_int_1::FT = A * (τ_z_1 - 1) + B * z * τ_z_3
    τ_int_2::FT = C * z * τ_z_3
    I_T::FT =
        (cos(φ) * (1 + γ * z / _a))^k -
        k / (k + 2) * (cos(φ) * (1 + γ * z / _a))^(k + 2)

    # base state virtual temperature, pressure, specific humidity, density
    T_v::FT = (τ_1 - τ_2 * I_T)^(-1)
    p::FT = _p_0 * exp(-_grav / _R_d * (τ_int_1 - τ_int_2 * I_T))

    # base state velocity
    U::FT =
        _grav * k / _a *
        τ_int_2 *
        T_v *
        (
            (cos(φ) * (1 + γ * z / _a))^(k - 1) -
            (cos(φ) * (1 + γ * z / _a))^(k + 1)
        )
    u_ref::FT =
        -_Ω * (_a + γ * z) * cos(φ) +
        sqrt((_Ω * (_a + γ * z) * cos(φ))^2 + (_a + γ * z) * cos(φ) * U)
    v_ref::FT = 0
    w_ref::FT = 0

    # velocity perturbations
    F_z::FT = 1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3
    if z > z_t
        F_z = FT(0)
    end
    d::FT = _a * acos(sin(φ) * sin(φ_c) + cos(φ) * cos(φ_c) * cos(λ - λ_c))
    c3::FT = cos(π * d / 2 / d_0)^3
    s1::FT = sin(π * d / 2 / d_0)
    if 0 < d < d_0 && d != FT(_a * π)
        u′::FT =
            -16 * V_p / 3 / sqrt(3) *
            F_z *
            c3 *
            s1 *
            (-sin(φ_c) * cos(φ) + cos(φ_c) * sin(φ) * cos(λ - λ_c)) /
            sin(d / _a)
        v′::FT =
            16 * V_p / 3 / sqrt(3) * F_z * c3 * s1 * cos(φ_c) * sin(λ - λ_c) /
            sin(d / _a)
    else
        u′ = FT(0)
        v′ = FT(0)
    end
    w′::FT = 0
    u_sphere = SVector{3, FT}(u_ref + u′, v_ref + v′, w_ref + w′)
    u_cart = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    ## Compute ure profile
    ## Pressure coordinate η
    ## η_crit = p_t / p_w ; p_t = 10000 hPa, p_w = 340 hPa
    η = p / _p_0
    if η > η_crit
        q_tot = q_0 * exp(-(φ / φ_w)^4) * exp(-((η - 1) * _p_0 / p_w)^2)
    else
        q_tot = q_t
    end
    phase_partition = PhasePartition(q_tot)

    ## temperature & density
    T::FT = T_v / (1 + M_v * q_tot)
    ρ::FT = air_density(bl.param_set, T, p, phase_partition)

    ## potential & kinetic energy
    e_pot::FT = gravitational_potential(bl.orientation, aux)
    e_kin::FT = 0.5 * u_cart' * u_cart
    e_tot::FT = total_energy(bl.param_set, e_kin, e_pot, T, phase_partition)

    ## Assign state variables
    state.ρ = ρ
    state.ρu = ρ * u_cart
    state.ρe = ρ * e_tot

    if bl.moisture isa EquilMoist
        state.moisture.ρq_tot = ρ * q_tot
    end

    nothing
end

"""
Defines analytical function for prescribed T_sfc and q_sfc, following
Thatcher and Jablonowski (2016), used to calculate bulk surface fluxes.

T_sfc_pole = SST at the poles (default: 271 K), specified above
"""
struct Varying_SST_TJ16{PS, O, MM}
    param_set::PS
    orientation::O
    moisture::MM
end
function (st::Varying_SST_TJ16)(state, aux, t)
    FT = eltype(state)
    φ = latitude(st.orientation, aux)

    _T_sfc_pole = T_sfc_pole(st.param_set)::FT

    Δφ = FT(26) * FT(π) / FT(180)       # latitudinal width of Gaussian function
    ΔSST = FT(29)                       # Eq-pole SST difference in K
    T_sfc = ΔSST * exp(-φ^2 / (2 * Δφ^2)) + _T_sfc_pole

    eps = FT(0.622)
    ρ = state.ρ

    q_tot = state.moisture.ρq_tot / ρ
    q = PhasePartition(q_tot)

    e_int = internal_energy(st.moisture, st.orientation, state, aux)
    T = air_temperature(st.param_set, e_int, q)
    p = air_pressure(st.param_set, T, ρ, q)

    _T_triple = T_triple(st.param_set)::FT          # triple point of water
    _press_triple = press_triple(st.param_set)::FT  # sat water pressure at T_triple
    _LH_v0 = LH_v0(st.param_set)::FT                # latent heat of vaporization at T_triple
    _R_v = R_v(st.param_set)::FT                    # gas constant for water vapor

    q_sfc =
        eps / p *
        _press_triple *
        exp(-_LH_v0 / _R_v * (FT(1) / T_sfc - FT(1) / _T_triple))

    return T_sfc, q_sfc
end

function config_baroclinic_wave(
    FT,
    param_set,
    poly_order,
    resolution,
    with_moisture,
)
    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))
    ref_state = HydrostaticState(temp_profile_ref)

    # Set up the atmosphere model
    exp_name = "BaroclinicWave"
    domain_height::FT = 30e3 # distance between surface and top of atmosphere (m)
    if with_moisture
        hyperdiffusion = EquilMoistBiharmonic(FT(8 * 3600))
        moisture = EquilMoist{FT}()
    else
        hyperdiffusion = DryBiharmonic(FT(8 * 3600))
        moisture = DryModel()
    end

    _C_drag = C_drag(param_set)::FT
    bulk_flux = Varying_SST_TJ16(param_set, SphericalOrientation(), moisture)
    problem = AtmosProblem(
        boundarycondition = (
            AtmosBC(
                energy = BulkFormulaEnergy(
                    (state, aux, t, normPu_int) -> _C_drag,
                    (state, aux, t) -> bulk_flux(state, aux, t),
                ),
                moisture = BulkFormulaMoisture(
                    (state, aux, t, normPu_int) -> _C_drag,
                    (state, aux, t) -> begin
                        _, q_tot = bulk_flux(state, aux, t)
                        q_tot
                    end,
                ),
            ),
            AtmosBC(),
        ),
        init_state_prognostic = init_baroclinic_wave!,
    )

    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        problem = problem,
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        hyperdiffusion = hyperdiffusion,
        moisture = moisture,
        source = (Gravity(), Coriolis()),
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_baroclinic_wave!;
        model = model,
    )

    return config
end

function main()
    # add a command line argument to specify whether to use a moist setup
    # TODO: this will move to the future namelist functionality
    bw_args = ArgParseSettings(autofix_names = true)
    add_arg_group!(bw_args, "BaroclinicWave")
    @add_arg_table! bw_args begin
        "--with-moisture"
        help = "use a moist setup"
        metavar = "yes|no"
        arg_type = String
        default = "yes"
    end

    cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = bw_args)
    moisture_arg = lowercase(cl_args["with_moisture"])
    if moisture_arg == "yes"
        with_moisture = true
    elseif moisture_arg == "no"
        with_moisture = false
    else
        error("invalid argument to --with-moisture: " * moisture_arg)
    end

    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 3                           # discontinuous Galerkin polynomial order
    n_horz = 12                              # horizontal element number
    n_vert = 6                               # vertical element number
    n_days::FT = 1
    timestart::FT = 0                        # start time (s)
    timeend::FT = n_days * day(param_set)    # end time (s)

    # Set up driver configuration
    driver_config = config_baroclinic_wave(
        FT,
        param_set,
        poly_order,
        (n_horz, n_vert),
        with_moisture,
    )

    # Set up experiment
    ode_solver_type = ClimateMachine.IMEXSolverType(
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
        split_explicit_implicit = true,
        discrete_splitting = false,
    )

    CFL = FT(0.1) # target acoustic CFL number

    # time step is computed such that the horizontal acoustic Courant number is CFL
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_config = config_diagnostics(FT, driver_config)

    # Set up user-defined callbacks
    filterorder = 20
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            ("moisture.ρq_tot",),
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        #user_callbacks = (cbtmarfilter, cbfilter),
        check_euclidean_distance = true,
    )
end

function config_diagnostics(FT, driver_config)
    interval = "40000steps" # chosen to allow a single diagnostics collection

    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(1), FT(1), FT(1000)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = setup_atmos_default_diagnostics(
        AtmosGCMConfigType(),
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

main()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
import OpenCL
using ModernGL, GLWindow, GLFW, GLAbstraction, Images
const cl = OpenCL

const width = 1024 # Also needs changing in qjulia_kernel.cl
const height = 1024 # Also needs changing in qjulia_kernel.cl

const window = createwindow("QJulia", width, height)

const glVersion = bytestring(glGetString(GL_VERSION))
const glVendor = bytestring(glGetString(GL_VENDOR))
const glRenderer = bytestring(glGetString(GL_RENDERER))

println("Using OpenGL $glVersion on $glRenderer provided by $glVendor")

const glLib = @windows? "opengl32" : @linux? "libGL" : ""

@windows? begin
    function getProperties()
        currentContext = ccall(("wglGetCurrentContext", glLib), Ptr{Void}, ())
        currentDC = ccall(("wglGetCurrentDC", glLib), Ptr{Void}, ())

        [(cl.CL_GL_CONTEXT_KHR, currentContext),
        (cl.CL_WGL_HDC_KHR, currentDC)]
    end

end : @linux? begin
    function getProperties()
        currentContext = ccall(("glXGetCurrentContext", glLib), Ptr{Void}, ())
        currentDC = ccall(("glXGetCurrentDisplay", glLib), Ptr{Void}, ())

        [(cl.CL_GL_CONTEXT_KHR, currentContext),
        (cl.CL_GLX_DISPLAY_KHR, currentDC)]
    end

end : error("Can't handle this")

function getDevices(platform)
    props = [getProperties(), (cl.CL_CONTEXT_PLATFORM, platform)]
    parsed_props = cl._parse_properties(props)

    nbytes = cl.Csize_t[0]
    err = cl.api.clGetGLContextInfoKHR(parsed_props, cl.CL_DEVICES_FOR_GL_CONTEXT_KHR, cl.C_NULL, cl.C_NULL, nbytes)
    if (err != cl.CL_SUCCESS)
        error("Failed to query devices! ", err)
    end

    ndevices = div(nbytes[1], sizeof(cl.CL_device_id))
    devices = Array(cl.CL_device_id, ndevices)

    err = cl.api.clGetGLContextInfoKHR(parsed_props, cl.CL_DEVICES_FOR_GL_CONTEXT_KHR, nbytes[1], devices, cl.C_NULL)
    if (err != cl.CL_SUCCESS)
        error("Failed to obtain devices! ", err)
    end

    return [cl.Device(id) for id in devices]
end

# Given a filter function
function getDevice(f :: Function)
    devices = cl.Device[]
    platforms = cl.platforms()

    isempty(platforms) && error("Could not find any OpenCL platforms.")

    for platform in platfroms
        append!(devices, getDevices(platform))
    end

    devices = filter(f, devices)

    isempty(devices) && error("Could not get any devices")

    gpu_devices = filter(devices) do dev
        dev[:device_type] == :gpu
    end

    isempty(gpu_devices) ? first(devices) : first(gpu_devices)
end

const device = getDevice(dev -> "cl_khr_gl_sharing" in dev[:extensions] && dev[:has_image_support])
const platform = device[:platform]
const clVersion = cl.opencl_version(platform)

println("Using OpenCL $clVersion on $(cl.info(device, :name)) driven by platform $(cl.info(platform, :name))")

const props = [getProperties(), (cl.CL_CONTEXT_PLATFORM, platform)]

# Setup OpenCL
const ctx = cl.Context(device, properties = props)
const queue = cl.CmdQueue(ctx)

# Setup program parameter

gl_texture = Texture(C_NULL, GL_TEXTURE_2D, GL_RGBA8, [width, height], GL_RGBA, GL_UNSIGNED_BYTE)

err_code = Array(cl.CL_int, 1)

if clVersion < v"1.2.0"
	const image = cl.api.clCreateFromGLTexture2D(ctx.id, cl.CL_MEM_READ_WRITE, GL_TEXTURE_2D, 0, gl_texture.id, err_code)
else
	const image = cl.api.clCreateFromGLTexture(ctx.id, cl.CL_MEM_READ_WRITE, GL_TEXTURE_2D, 0, gl_texture.id, err_code)
end
if err_code[1] != cl.CL_SUCCESS
    error(err_code[1])
else
    println("cl texture created from opengl texture successfully")
end

const ε = 0.0003f0

const colorT = 0.0f0
const colorA = [0.25f0, 0.45f0, 1.0f0, 1.0f0]
const colorB = [0.9f0, 0.45f0, 0.1f0, 1.0f0]
const colorC = [0.7f0, 0.7f0, 1.0f0, 1.0f0]

const μT = 0.0f0
const μA = [-0.278f0, -0.479f0, 0.2f0, 0.0f0]
const μB = [0.278f0, 0.479f0, 0.0f0, 0.4f0]
const μC = [-0.278f0, -0.03f0, -0.24f0, 0.235f0]
#const μC = [0.285f0, -0.001f0, -0.231f0, 0.235f0]
const kernelsource = open(readall, joinpath(dirname(Base.source_path()), "qjulia_kernel.cl"))
const qjulia_program = cl.Program(ctx, source=kernelsource) |> cl.build!
const qjulia_kernel = cl.Kernel(qjulia_program, "QJuliaKernel")

# Setup clImage and clBuffer
# Original Code CreateComputeResult

# Setup opencl image ref
#const image = cl.Image{cl.RGBA, Float32}(ctx, :w, shape = (width, height, 1, 1, 1, 1)) # OpenCL.jl excpects a shape of dim + nchannels
# Setup OCL buffer
const buffer = cl.Buffer(Float32, ctx, :w, sizeof(Float32) * cl.nchannels(cl.RGBA) * width * height)

function compute()
	# Blocking call to compute the context
    cl.call(queue, qjulia_kernel, (width, height), nothing, buffer, μC..., colorC..., ε)

    glFinish()
    ret_event = Array(cl.CL_event, 1)

    err = cl.api.clEnqueueAcquireGLObjects(queue.id, cl.cl_uint(1), [image], cl.cl_uint(0), cl.C_NULL, ret_event)
    if (err != cl.CL_SUCCESS)
        error("Failed to acquire GL object! ", err)
    end
    evt_1 = cl.Event(ret_event[1], retain=false)

    origin = Csize_t[ 0, 0, 0 ]
    region = Csize_t[width, height, 1 ]

    err = cl.api.clEnqueueCopyBufferToImage(queue.id, buffer.id, image,
                                     0, origin, region, cl.cl_uint(1), [evt_1.id], ret_event)

    if err != cl.CL_SUCCESS
        println("Failed to copy buffer to image! ", err)
    end
    evt_2 = cl.Event(ret_event[1], retain=false)

    err = cl.api.clEnqueueReleaseGLObjects(queue.id, 1, [image], 1, [evt_2.id], ret_event)
    if err != cl.CL_SUCCESS
        println("Failed to release GL object! ", err)
    end

    evt_3 = cl.Event(ret_event[1], retain=false)
    cl.wait(evt_3)
    cl.finish(queue)
end

# TODO workgroup size

function main()
	compute()

    const fullscreenQuad = RenderObject([
        :position       => GLBuffer(GLfloat[-1,-1, -1,1, 1,1, 1,-1], 2),
        :indexes        => GLBuffer(GLuint[0, 1, 2,  2, 3, 0], 1, bufferType = GL_ELEMENT_ARRAY_BUFFER),
        :uv             => GLBuffer(GLfloat[0,1,  0,0,  1,0, 1,1], 2),
        :fullscreenTex  => gl_texture
    ], GLProgram("simple"))
    glClearColor(1f0, 1f0, 1f0, 0f0)
    glViewport(0,0,width, height)
    runner = 0.0f0
    while !GLFW.WindowShouldClose(window)

        μC[1] = float32(sin(runner))
        μC[2] = float32(sin(runner /0.8f0))
        μC[3] = float32(sin(runner /0.6f0))
        μC[4] = float32(sin(runner /0.4f0))

        runner += 0.01f0
        compute()
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

        glDisable(GL_DEPTH_TEST)
        programID = fullscreenQuad.vertexArray.program.id
        if programID!= glGetIntegerv(GL_CURRENT_PROGRAM)
            glUseProgram(programID)
        end
        render(fullscreenQuad.uniforms)
        glBindVertexArray(fullscreenQuad.vertexArray.id)
        glDrawElements(GL_TRIANGLES, fullscreenQuad.vertexArray.indexLength, GL_UNSIGNED_INT, GL_NONE)
        # Swap front and back buffers
        GLFW.SwapBuffers(window)

        # Poll for and process events
        GLFW.PollEvents()
    end

    GLFW.Terminate()
end

main()


text:
using Documenter
using OrthographicProjection

# makedocs(
# 	format = Documenter.HTML(),
# 	sitename = "OrthographicProjection.jl",
# 	#assets = ["assets/OrthographicProjection.css", "assets/logo.jpg"],
# 	modules = [OrthographicProjection]
# )
#
#
# deploydocs(
# 	repo = "github.com/marteresagh/OrthographicProjection.jl.git"
# )

makedocs(
    modules = [OrthographicProjection],
    format = Documenter.HTML(
        prettyurls = "deploy" in ARGS,
    ),
    sitename = "OrthographicProjection.jl",
    pages = [
        "Home" => "index.md",
        "Description" => "description.md",
        "Results" => "results.md",
        "Script" => "script.md",
        "References" => "refs.md",
    ]
)

deploydocs(
    repo = "github.com/marteresagh/OrthographicProjection.jl.git",
    push_preview = true,
    devurl = "dev",
)


text:
module frame_

export
	frame

include("../src/libiup_h.jl")

using IUP

#------------------------example html/examples/C/frame.c ----------------------------------
# IupFrame Example
# Draws a frame around a button. Note that "FGCOLOR" is added to the frame but it is inherited by the button.

function frame()
	IupOpen()	#Initializes IUP

	# Creates frame with a label
	thisframe = IupFrame(
				IupHbox(IupFill(),
					IupLabel ("IupFrame Attributes:\nFGCOLOR = \"255 0 0\"\nSIZE = \"EIGHTHxEIGHTH\"\nTITLE = \"This is the frame\"\nMARGIN = \"10x10\""),
					IupFill(),
 				)
			)

	# Sets frame's attributes
	IupSetAttributes(thisframe, "FGCOLOR=\"255 0 0\", SIZE=EIGHTHxEIGHTH, TITLE=\"This is the frame\", MARGIN=10x10")

	# Creates dialog
	dialog = IupDialog(thisframe)

	# Sets dialog's title
	IupSetAttribute(dialog, "TITLE", "IupFrame")

	IupShow(dialog)                          # Shows dialog in the center of the screen */
	IupMainLoop()                            # Initializes IUP main loop
	IupClose()                               # And close it when ready
end

end 	# module


text:
## parallel elastic component
# passive joint torques from Reiner and Edrich (1999)

# convert model joint angles to joint angle definitions from paper
convertHip(θ) = rad2deg(π - θ)
convertKnee(θ) = rad2deg(π - θ)
convertAnkle(θ) = rad2deg(θ - (π / 2))

# passive torques
torqueHip(ϕh, ϕk) = ℯ^(1.4655 - .0034ϕk - .075ϕh)  - ℯ^(1.3403 - .0226ϕk + .0305ϕh) + 8.072
torqueKnee(ϕh, ϕk, ϕa) = ℯ^(1.8 - .046ϕa - .0352ϕk + .0217ϕh) - ℯ^(-3.971 - .0004ϕa + .0495ϕk - .0128ϕh) - 4.82 + ℯ^(2.22 - .15ϕk)
torqueAnkle(ϕk, ϕa) = ℯ^(2.1016 - .0843ϕa - .0176ϕk) - ℯ^(-7.9763 + .1949ϕa + .0008ϕk) - 1.792

function PEtorque(θh, θk, θa)
    # convert to match joint angle definitions
    ϕh = convertHip(θh)
    ϕk = convertKnee(θk)
    ϕa = convertAnkle(θa)

    # calculate passive torques
    τh = torqueHip(ϕh, ϕk)
    τk = torqueKnee(ϕh, ϕk, ϕa)
    τa = torqueAnkle(ϕk, ϕa)

    return τh, τk, τa
end


text:
"""
    get_affine(vol::Vector{DICOM.DICOMData})
    get_affine(slice::DICOM.DICOMData)

Extracts the affine matrix from the DICOM header of the input
volume or slice.

# Arguments

- `vol::Vector{DICOM.DICOMData}`: 3D DICOM volume or sequence of 2D DICOM
    slices
- `slice::DICOM.DICOMData`: 2D DICOM slice
"""
function get_affine(vol::Vector{DICOM.DICOMData})
	ImageOrientation = (0x0020, 0x0037)
	PixelSpacing = (0x0028, 0x0030)
	ImagePosition = (0x0020, 0x0032)

	img_ornt = vol[1][ImageOrientation]
	pix_space = vol[1][PixelSpacing]
	img_position = vol[1][ImagePosition]

	F11, F21, F31 = img_ornt[1:4]
	F12, F22, F32 = img_ornt[4:end]
	dr, dc = pix_space
	Sx, Sy, Sz = img_position

	a = [F11*dr F12*dc 0 Sx
		 F21*dr F22*dc 0 Sy
		 F31*dr F32*dc 0 Sz
		 0 0 0 1]
    return a
end

function get_affine(slice::DICOM.DICOMData)
	ImageOrientation = (0x0020, 0x0037)
	PixelSpacing = (0x0028, 0x0030)
	ImagePosition = (0x0020, 0x0032)

	img_ornt = vol[1][ImageOrientation]
	pix_space = vol[1][PixelSpacing]
	img_position = vol[1][ImagePosition]

	F11, F21, F31 = img_ornt[1:4]
	F12, F22, F32 = img_ornt[4:end]
	dr, dc = pix_space
	Sx, Sy, Sz = img_position

	a = [F11*dr F12*dc 0 Sx
		 F21*dr F22*dc 0 Sy
		 F31*dr F32*dc 0 Sz
		 0 0 0 1]
    return a
end

"""
    io_orientation(affine, tol=nothing)

Orientation of input axes in terms of output axes for `affine`

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L22)

Parameters
----------
affine : (q+1, p+1) Array
Transformation affine from ``p`` inputs to ``q`` outputs.  Usually this
will be a shape (4,4) matrix, transforming 3 inputs to 3 outputs, but
the code also handles the more general case
tol : {None, float}, optional
threshold below which SVD values of the affine are considered zero. If
`tol` is None, and `S` is an array with singular values for `affine`,
and `eps` is the epsilon value for datatype of `S`, then `tol` set
to `maximum(S) * maximum(size(RS)) * floatmin(typeof(S[1]))`

Returns
-------
orientations : (p, 2) Array
one row per input axis, where the first value in each row is the closest
corresponding output axis. The second value in each row is 1 if the
input axis is in the same direction as the corresponding output axis and
-1 if it is in the opposite direction.  If a row is [NaN, NaN],
which can happen when p > q, then this row should be considered dropped.
"""
function io_orientation(affine, tol=nothing)
    q, p = size(affine)[1] - 1, size(affine)[2] - 1
    RZS = affine[1:q, 1:p]
    A = RZS * RZS
    zooms = sqrt.(sum(A; dims=1))
    zooms[map(x -> x == 0.0, zooms)] .= 1
    RS = RZS ./ zooms
    P, S, Qs = svd(collect(RS))
    if (tol === nothing)
        tol = maximum(S) * maximum(size(RS)) * floatmin(Float64.(S[1]))
    end
    keep = S .> tol
    R = P[:, keep] * Qs[keep, :]
    ornt = Int8.(ones((p, 2))) .* NaN
    for in_ax in 1:p
        col = R[:, in_ax]
        if ((all(x -> x == 0.0, col)) == false)
            out_ax = argmax(abs.(col))
            ornt[in_ax, 1] = out_ax
            @assert col[out_ax] != 0
            if col[out_ax] < 0
                ornt[in_ax, 2] = -1
            else
                ornt[in_ax, 2] = 1
            end
            R[out_ax, :] .= 0
        end
    end
    return ornt
end

"""
	axcodes2ornt(axcodes, labels=nothing)

Convert axis codes `axcodes` to an orientation

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L309)

Parameters
----------
axcodes : (N,) tuple
	axis codes - see ornt2axcodes docstring
labels : optional, None or sequence of (2,) sequences
	(2,) sequences are labels for (beginning, end) of output axis.  That
	is, if the first element in `axcodes` is ``front``, and the second
	(2,) sequence in `labels` is ('back', 'front') then the first
	row of `ornt` will be ``[1, 1]``. If None, equivalent to
	``(('L','R'),('P','A'),('I','S'))`` - that is - RAS axes.
Returns
-------
ornt : (N,2) array-like
	orientation array - see io_orientation docstring

"""
function axcodes2ornt(axcodes, labels=nothing)
    if (labels === nothing)
        labels = zip(["L" "P" "I"], ["R", "A", "S"])
    else
        labels = labels
    end
    n_axes = length(axcodes)
    ornt = Int8.(ones(n_axes, 2)) .* NaN
    for (code_idx, code) in enumerate(axcodes)
        for (label_idx, codes) in enumerate(labels)
            # println(label_idx, codes)
            if (code == nothing)
                continue
            end
            if (code in codes)
                if (code == codes[1])
                    ornt[code_idx, :] = [label_idx, -1]
                else
                    ornt[code_idx, :] = [label_idx, 1]
                    break
                end
            end
        end
    end
    return ornt
end

"""
    ornt2axcodes(ornt, labels=nothing)

Convert orientation `ornt` to labels for axis directions

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L309)

Parameters
----------
ornt : (N,2) array-like
    orientation array - see io_orientation docstring
labels : optional, None or sequence of (2,) sequences
    (2,) sequences are labels for (beginning, end) of output axis.  That
    is, if the first row in `ornt` is ``[1, 1]``, and the second (2,)
    sequence in `labels` is ('back', 'front') then the first returned axis
    code will be ``'front'``.  If the first row in `ornt` had been
    ``[1, -1]`` then the first returned value would have been ``'back'``.
    If None, equivalent to ``(('L','R'),('P','A'),('I','S'))`` - that is -
    RAS axes.
Returns
-------
axcodes : (N,) tuple
    labels for positive end of voxel axes.  Dropped axes get a label of
    None.
"""
function ornt2axcodes(ornt, labels=nothing)
	if labels == nothing
		labels = collect(zip(["L" "P" "I"], ["R", "A", "S"]))
	end
    axcodes = []
	 for row in eachrow(ornt)
		axno, direction = row[1], row[2]
        if isnan(axno)
            push!(axcodes, nothing)
            continue
		end
		axint = Int(round(axno))
		if axint != axno
			error("Non integer axis number $(axno)")
		elseif direction == 1
			axcode = labels[axint][2]
		elseif direction == -1
			axcode = labels[axint][1]
		else
			error("Direction should be -1 or 1")
		end
		push!(axcodes, axcode)
	end
    return axcodes
end

"""
	ornt_transform(start_ornt, end_ornt)

Return the orientation that transforms from `start_ornt` to `end_ornt`

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L309)

Parameters
----------
start_ornt : (n,2) orientation array
	Initial orientation.
end_ornt : (n,2) orientation array
	Final orientation.

Returns
-------
orientations : (p, 2) array
	The orientation that will transform the `start_ornt` to the `end_ornt`.
"""
function ornt_transform(start_ornt, end_ornt)
    @assert size(start_ornt) == size(end_ornt)

    result = Array{Float64}(undef, size(start_ornt))
    for (end_in_idx, (end_out_idx, end_flip)) in enumerate(eachrow(end_ornt))
        for (start_in_idx, (start_out_idx, start_flip)) in enumerate(eachrow(start_ornt))
            if (end_out_idx == start_out_idx)
                if (start_flip == end_flip)
                    flip = 1
                else
                    flip = -1
                end
                result[start_in_idx, :] = [end_in_idx, flip]
                break
            end
        end
    end
    return result
end

"""
	apply_orientation(arr, ornt)

Apply transformations implied by `ornt` to the first n axes of the array `arr`

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L309)

Parameters
----------
arr : array-like of data with ndim >= n
ornt : (n,2) orientation array
   orientation transform. ``ornt[N,1]` is flip of axis N of the
   array implied by `shape`, where 1 means no flip and -1 means
   flip.  For example, if ``N==0`` and ``ornt[0,1] == -1``, and
   there's an array ``arr`` of shape `shape`, the flip would
   correspond to the effect of ``np.flipud(arr)``.  ``ornt[:,0]`` is
   the transpose that needs to be done to the implied array, as in
   ``arr.transpose(ornt[:,0])``

Returns
-------
t_arr : array
   data array `arr` transformed according to ornt
"""
function apply_orientation(arr, ornt)
    t_arr = arr
    n = size(ornt)[1]
    if (length(size(t_arr)) < n)
        error("Data array has fewer dimensions than `orientation`")
    end
    for (ax, flip) in enumerate(ornt[:, 2])
        if (flip == -1)
            t_arr = reverse(t_arr; dims=ax)
        end
    end
    full_transpose = collect(1:length(size(t_arr)))
    full_transpose[1:n] = sortperm(ornt[:, 1])
    return t_arr = permutedims(t_arr, full_transpose)
end

"""
	inv_ornt_aff(ornt, shape::AbstractArray)

Affine transform reversing transforms implied in `ornt`

Ported from nibabel
(https://github.com/nipy/nibabel/blob/e51bcb43d9c6f5ad329ffb230afda0c26b9e8617/nibabel/orientations.py#L309)

Imagine you have an array ``arr`` of shape `shape`, and you apply the
    transforms implied by `ornt` (more below), to get ``tarr``.
    ``tarr`` may have a different shape ``shape_prime``.  This routine
    returns the affine that will take a array coordinate for ``tarr``
    and give you the corresponding array coordinate in ``arr``.
    Parameters
    ----------
    ornt : (p, 2) ndarray
       orientation transform. ``ornt[P, 1]` is flip of axis N of the array
       implied by `shape`, where 1 means no flip and -1 means flip.  For
       example, if ``P==0`` and ``ornt[0, 1] == -1``, and there's an array
       ``arr`` of shape `shape`, the flip would correspond to the effect of
       ``np.flipud(arr)``.  ``ornt[:,0]`` gives us the (reverse of the)
       transpose that has been done to ``arr``.  If there are any NaNs in
       `ornt`, we raise an ``OrientationError`` (see notes)
    shape : length p sequence
       shape of array you may transform with `ornt`

    Returns
    -------
    transform_affine : (p + 1, p + 1) ndarray
       An array ``arr`` (shape `shape`) might be transformed according to
       `ornt`, resulting in a transformed array ``tarr``.  `transformed_affine`
       is the transform that takes you from array coordinates in ``tarr`` to
       array coordinates in ``arr``.

    Notes
    -----
    If a row in `ornt` contains NaN, this means that the input row does not
    influence the output space, and is thus effectively dropped from the output
    space.  In that case one ``tarr`` coordinate maps to many ``arr``
    coordinates, we can't invert the transform, and we raise an error

"""
function inv_ornt_aff(ornt, shape::AbstractArray)
    if any(x -> x === NaN, ornt)
        error("Cannot invert the orientation transform")
    end
    p = size(ornt)[1]
    shape = shape[1:p]
    axis_transpose = [Int(v) for v in ornt[:, 1]]
    undo_reorder = I(p + 1)[append!(axis_transpose, p + 1), :]
    undo_flip = Float64.(diagm(append!(ornt[:, 2], 1)))
    center_trans = -(shape .- 1) ./ 2.0
    undo_flip[1:p, p + 1] = (ornt[:, 2] .* center_trans) .- center_trans
    return undo_flip * undo_reorder
end

"""
	orientation(data_array, axcodes=nothing, affine=nothing, as_closest_canonical=false,
				labels=zip(["L" "P" "I"], ["R", "A", "S"]))

Change the input image's orientation into the specified based on `axcodes`

Ported from Monai
(https://docs.monai.io/en/latest/_modules/monai/transforms/spatial/array.html#Orientation)
"""
function orientation(
    data_array,
    axcodes=nothing,
    affine=nothing,
    as_closest_canonical=false,
    labels=zip(["L" "P" "I"], ["R" "A" "S"]),
)
    sr = length(size(data_array)) - 1
    if (sr ≤ 0)
        error("data_array must have at least one spatial dimension.")
    end
    if (affine === nothing)
        affine = Float64.(I(sr + 1))
        affine_ = Float64.(I(sr + 1))
    else
        affine_ = affine
    end
    src = io_orientation(affine_)

    if as_closest_canonical
        spatial_ornt = src
    else
        if (axcodes === nothing)
            error("error")
        end
        dst = axcodes2ornt(axcodes[1:sr], labels)
        if (length(dst) .< sr)
            error("axcodes must match data_array spatially")
        end
        spatial_ornt = ornt_transform(src, dst)
    end
    ornt = copy(spatial_ornt)
    ornt[:, 1] .+= 1
    ornt = cat([1 1], ornt; dims=1)
    shape = size(data_array)[2:end]
    data_array = apply_orientation(data_array, ornt)
    new_affine = inv_ornt_aff(spatial_ornt, collect(shape))
    return data_array, affine, new_affine
end


text:
# Commonly used constants. The factors convert from the respective unit
# to atomic units
module units
    const Ha = 1                   # Hartree -> Hartree
    const eV = 0.03674932248       # electron volt -> Hartree
    const Å  = 1.8897261254578281  # Ǎngström -> Bohr
    const K  = 3.166810508e-6      # Kelvin -> Hartree
    const Ry = 0.5                 # Rydberg -> Hartree
    const Ǎ  = Å                   # Ǎngström -> Bohr  (deprecated)
end

"""
    unit_to_ao(symbol)
Get the factor converting from the unit `symbol` to atomic units.
E.g. `unit_to_au(:eV)` returns the conversion factor from electron volts to Hartree.
"""
unit_to_au(symbol::Symbol) = getfield(units, symbol)


text:
@safetestset "Process New Patterns" begin
    using ParameterSpacePartitions, Test
    import ParameterSpacePartitions: process_new_patterns!, Chain, is_new

    all_patterns = [[1,2],[1,3]]
    patterns = [[1,2],[1,4]]

    chains = [Chain([.3,.3], [1,2], .2),Chain([.3,.3], [1,3], .2)]
    chains[1].acceptance[1] = false
    chains[2].acceptance[1] = false

    parms = [[.3,.4],[.3,.5]]
    options = Options(;
        radius = .2,
        bounds = [(0,1),(0,1)],
        n_iters = 100,
        init_parms = [[.2,.3]]
    )

    @test !is_new(all_patterns, patterns[1])
    @test is_new(all_patterns, patterns[2])

    process_new_patterns!(all_patterns, patterns, parms, chains, options)

    @test length(chains) == 3
    @test chains[3].pattern == [1,4]
end

@safetestset "in_bounds" begin
    using ParameterSpacePartitions, Test
    import ParameterSpacePartitions: in_bounds, Chain

    bounds = [(0,1),(0,1)]
    parms = [.5,.5]
    in_bounds(parms, bounds)
    @test in_bounds(parms, bounds)

    bounds = [(0,1),(0,1)]
    parms = [.5,-1]
    @test !in_bounds(parms, bounds)

    bounds = [(0,1),(0,1)]
    parms = [1.2,.5]
    @test !in_bounds(parms, bounds)
end

@safetestset "update_position!" begin
    using ParameterSpacePartitions
    using Test
    import ParameterSpacePartitions: update_position!, Chain

    parms = [.4]
    proposal = [.3]
    pattern = [4,3]
    bounds = [(0,1)]
    chain = Chain(parms, pattern, .3)
    update_position!(chain, proposal, pattern, bounds)
    @test chain.parms == [.3]
    @test chain.pattern == [4,3]
    @test chain.acceptance[2] == true

    parms = [.4]
    proposal = [.3]
    pattern = [4,3]
    chain = Chain(parms, pattern, .3)
    update_position!(chain, proposal, [4,5], bounds)
    @test chain.parms == [.4]
    @test chain.pattern == [4,3]
    @test chain.acceptance[2] == false

    parms = [.4]
    proposal = [-.3]
    pattern = [4,3]
    chain = Chain(parms, pattern, .3)
    update_position!(chain, proposal, pattern, bounds)
    @test chain.parms == [.4]
    @test chain.pattern == [4,3]
    @test chain.acceptance[2] == false

end

@safetestset "Cube" begin
    using ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Test, Random, Distributions
    using DataFrames

    Random.seed!(541)
    # dimensions of the hypbercue
    n_dims = 3
    # partitions per dimension
    n_part = 4
    # number of starting points
    n_start = 1

    # partition boundaries
    p_bounds = range(0, 1, n_part + 1)
    hypercube = HyperCube(p_bounds)
    bounds = fill((0, 1), n_dims)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .18,
        bounds,
        n_iters = 100,
        parallel = false,
        init_parms
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        hypercube
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1, :p2, :p3]
    )

    groups = groupby(df, :pattern)
    length(groups)
    combine(groups,
        :p1 => minimum, :p1 => maximum,
        :p2 => minimum, :p2 => maximum,
        :p3 => minimum, :p3 => maximum
    ) |> sort

    @test length(groups) == n_part^n_dims
end

@safetestset "5D Hypercube" begin
    using ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Test, Random, Distributions
    using DataFrames

    Random.seed!(103)
    # partitions per dimension
    n_part = 4
    # dimensions of hypbercue
    n_dims = 5
    # number of starting points
    n_start = 1

    # bounds of the hypercube
    bounds = fill((0, 1), n_dims)
    # partition boundaries
    p_bounds = range(0, 1, n_part + 1)
    hypercube = HyperCube(p_bounds)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .15,
        bounds,
        n_iters = 100,
        parallel = false,
        init_parms
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        hypercube
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1,:p2,:p3,:p4,:p5]
    )

    groups = groupby(df, :pattern)
    combine(groups,
        :p1 => minimum, :p1 => maximum,
        :p2 => minimum, :p2 => maximum,
        :p3 => minimum, :p3 => maximum,
        :p4 => minimum, :p4 => maximum,
        :p5 => minimum, :p5 => maximum
    ) |> sort

    @test length(groups) == n_part^n_dims
end

@safetestset "5D Polytope" begin
    using ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Test, Random, Distributions
    using DataFrames, LinearAlgebra
    Random.seed!(778)

    # dimensions of the hypbercue
    n_dims = 5
    # number of partitions
    n_part = 50
    # number of starting points
    n_start = 1

    # partition boundaries
    polytopes = [Polytope(rand(n_dims)) for i in 1:n_part]
    bounds = fill((0, 1), n_dims)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .1,
        bounds,
        n_iters = 500,
        parallel = false,
        init_parms
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        polytopes
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1, :p2, :p3, :p4, :p5]
    )

    groups = groupby(df, :pattern)
    @test length(groups) == n_part
end

@safetestset "HyperSpheres" begin
    using Test, ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using LinearAlgebra, Random, DataFrames, Distributions
    Random.seed!(383)

    # dimensions of the hypbercue
    n_dims = 7
    # number of hyperspheres
    n_obj = 50
    # distribution of radii
    r_dist = () -> rand(Uniform(.2, .25))
    # number of starting points
    n_start = 10

    # partition boundaries
    hyperspheres = set_locations(r_dist, n_dims, n_obj)
    bounds = fill((0, 1), n_dims)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .15,
        bounds,
        n_iters = 5000,
        init_parms,
        λ = .05
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        hyperspheres
    )

    df = DataFrame(results)
    groups = groupby(df, :pattern)

    @test length(groups) == (n_obj + 1)
end

@safetestset "Adapt Radius 5D HyperCube" begin
    using ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Test, Random, Distributions
    using DataFrames

    Random.seed!(587)
    # partitions per dimension
    n_part = 4
    # dimensions of hypbercue
    n_dims = 5
    # number of starting points
    n_start = 1

    # bounds of the hypercube
    bounds = fill((0, 1), n_dims)
    # partition boundaries
    p_bounds = range(0, 1, n_part + 1)
    hypercube = HyperCube(p_bounds)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    t_rate = .4

    options = Options(;
        radius = .40,
        bounds,
        n_iters = 500,
        parallel = false,
        init_parms,
        t_rate,
        λ = .0
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        hypercube
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1,:p2,:p3,:p4,:p5]
    )

    groups = groupby(df, :chain_id)
    mean_accept = combine(groups, :acceptance=>mean=>:mean)
    # show that initial radius results in poor acceptance
    @test mean(mean_accept.mean) < .10

    options = Options(;
        radius = .40,
        bounds,
        n_iters = 500,
        parallel = false,
        init_parms,
        t_rate,
        λ = .2
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        hypercube
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1,:p2,:p3,:p4,:p5]
    )

    groups = groupby(df, :chain_id)
    mean_accept = combine(groups, :acceptance=>mean=>:mean)
    @test mean(mean_accept.mean) ≈ t_rate atol = .03
    @test std(mean_accept.mean) < .05
end

@safetestset "Adapt Radius 5D Polytope" begin
    using ParameterSpacePartitions
    using ParameterSpacePartitions.TestModels
    using Test, Random, Distributions
    using DataFrames, LinearAlgebra
    Random.seed!(2002)

    # dimensions of the hypbercue
    n_dims = 5
    # number of partitions
    n_part = 50
    # number of starting points
    n_start = 1

    # partition boundaries
    polytopes = [Polytope(rand(n_dims)) for i in 1:n_part]
    bounds = fill((0, 1), n_dims)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .5,
        bounds,
        n_iters = 500,
        parallel = false,
        init_parms,
        adapt_radius! = no_adaption!
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        polytopes
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1, :p2, :p3, :p4, :p5]
    )

    groups = groupby(df, :pattern)
    @test length(groups) == n_part

    groups = groupby(df, :chain_id)
    mean_accept = combine(groups, :acceptance=>mean=>:mean)
    # show that initial radius results in poor acceptance
    @test mean(mean_accept.mean) < .10

    t_rate = .4

    options = Options(;
        radius = 1.0,
        bounds,
        n_iters = 1000,
        parallel = false,
        init_parms,
        λ = 0.2,
        t_rate
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        polytopes
    )

    df = DataFrame(results)
    transform!(
        df,
        :parms => identity => [:p1, :p2, :p3, :p4, :p5]
    )

    groups = groupby(df, :pattern)
    @test length(groups) == n_part

    groups = groupby(df, :chain_id)
    mean_accept = combine(groups, :acceptance=>mean=>:mean)
    # show that initial radius results in poor acceptance
    @test mean(mean_accept.mean) ≈ t_rate atol = .04
    @test std(mean_accept.mean) < .05
end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
using LazySets, MathematicalPredicates

function motor_specification()
    # initial set:
    # xᵢ ∈ [0.002, 0.0025] for i = 1
    # xᵢ ∈ [0.001, 0.0015] for i = 5
    # xᵢ = 0 otherwise
    X0 = Hyperrectangle(low=[0.002, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0],
                        high=[0.0025, 0.0, 0.0, 0.0, 0.0015, 0.0, 0.0, 0.0])
    X1 = Hyperrectangle([0.00225, 0.0, 0.0, 0.0, 0.00125, 0.0, 0.0, 0.0],
                        [0.00025, 0.0, 0.0, 0.0, 0.00025, 0.0, 0.0, 0.0])

    # safety property: x1 ≤ 0.35 ∨ x5 ≤ 0.45
    property = Disjunction([
        is_contained_in(HalfSpace([1.; zeros(7)], 0.35)),
        is_contained_in(HalfSpace([zeros(4); 1.; zeros(3)], 0.45))])

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# functions relating to counter examples of the bapat sunder conjecture
# and its physical realization

"""

    cholesky_semi_definite_positive(A)

cholesky decomposition (`A` = R' * R) for a sdp but not strictly positive
definite matrix
"""
function cholesky_semi_definite_positive(A)

	# method found in some forum
	X = sqrt(A)

	QR_dec = qr(X)

	R = QR_dec.R

	try
		@test R' * R ≈ A
		return R
	catch err
		throw(Exception(err))
		# println(err)
		# println(R)
		# pretty_table(R)
	end

end

"""

    incorporate_in_a_unitary(X)

incorporates the renormalized matrix X in a double sized unitary through the
proof of Lemma 29  of Aaronson Arkipov seminal [The Computational Complexity of Linear Optics](https://arxiv.org/abs/1011.3245)
"""
function incorporate_in_a_unitary(X)

	Y = X / norm(X)

	I_yy = Matrix(I, size(Y)) - Y'*Y

	@test is_hermitian(I_yy)
	@test is_positive_semidefinite(I_yy)

	Z = cholesky_semi_definite_positive(I_yy)

	@test Y'*Y + Z'*Z ≈ Matrix(I, size(Y))

	W = rand(ComplexF64, 2 .* size(Y))

	n = size(Y)[1]
	W[1:n, 1:n] = Y
	W[n+1:2n, 1:n] = Z

	W = modified_gram_schmidt(W)

	@test is_unitary(W)

	W

end

"""

	incorporate_in_a_unitary_non_square(X)

same as `incorporate_in_a_unitary` but for a matrix renormalized `X` of type (m,n) with m >= n
generates a minimally sized unitary (ex 9*9 interferometer for the 7*2 M' of the first counter example of drury)
"""
function incorporate_in_a_unitary_non_square(X)

	# pad X with zeros

	(m,n) = size(X)

	if m<n
		throw(ArgumentError("m<n not implemented"))
	end

	#padding with zeros
	X_extended = zeros(eltype(X), (m,m))
	X_extended[1:m, 1:n] = X

	#now the same as incorporate_in_a_unitary but adapted to take only the minimally necessary elements
	Y = X_extended / norm(X_extended)

	I_yy = Matrix(I, size(Y)) - Y'*Y

	@test is_hermitian(I_yy)
	@test is_positive_semidefinite(I_yy)

	Z = cholesky_semi_definite_positive(I_yy)

	@test Y'*Y + Z'*Z ≈ Matrix(I, size(Y))

	W = rand(ComplexF64, (m+n, m+n))

	W[1:m, 1:n] = Y[1:m, 1:n] #not taking the padding
	W[m+1:m+n, 1:n] = Z[1:n, 1:n]  #taking only the minimal elements necessary to make the first n columns an orthonormal basis

	W = modified_gram_schmidt(W)

	@test is_unitary(W)

	W

end

"""

	add_columns_to_make_square_unitary(M_dagger)

Makes a square matrix `U` of which the first columns are `M_dagger`, which has to be
unitary by itself by a random choice of vectors that are then orthonormalized.
"""
function add_columns_to_make_square_unitary(M_dagger)

	n,r = size(M_dagger)
	U = rand(ComplexF64, (n,n))

	U[1:n,1:r] = M_dagger

	U = modified_gram_schmidt(U)

	is_unitary(U)

	U

end


text:
abstract type HH3DTDBIO{T} <: RetardedPotential{T} end

struct HH3DSingleLayerTDBIO{T} <: HH3DTDBIO{T}
    "speed of light"
    speed_of_light::T
    "number of temporal differentiations"
    num_diffs::Int
end

HH3DSingleLayerTDBIO(c) = HH3DSingleLayerTDBIO(c, 0)

struct HH3DHyperSingularTDBIO{T} <: HH3DTDBIO{T}
    speed_of_light::T
    weight_of_weakly_singular_term::T
    weight_of_hyper_singular_term::T
    num_diffs_on_weakly_singular_term::Int
    num_diffs_on_hyper_singular_term::Int
end

function HH3DHyperSingularTDBIO(;speed_of_light, numdiffs)
    id = one(speed_of_light)
    HH3DHyperSingularTDBIO(speed_of_light, id, id, numdiffs+1, numdiffs-1)
end

struct HH3DDoubleLayerTDBIO{T} <: HH3DTDBIO{T}
    speed_of_light::T
    weight::T
    num_diffs::Int
end

HH3DDoubleLayerTDBIO(;speed_of_light) = HH3DDoubleLayerTDBIO(speed_of_light,one(speed_of_light),0)

# See: ?BEAST.quaddata for help
function quaddata(operator::HH3DTDBIO,
        test_local_space, trial_local_space, time_local_space,
        test_element, trial_element, time_element)

    dmax = numfunctions(time_local_space)-1
    bn = binomial.((0:dmax),(0:dmax)')

    V = eltype(test_element[1].vertices)
    ws = WiltonInts84.workspace(V)
    order = 4
    @show order
    quadpoints(test_local_space, test_element, (order,)), bn, ws

end

# See: ?BEAST.quadrule for help
function quadrule(operator::HH3DTDBIO,
        test_local_space, trial_local_space, time_local_space,
        p, test_element, q, trial_element, r, time_element,
        quad_data)

    # WiltonInts84Strat(quad_data[1,p])
    qd = quad_data
    WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])

end

function innerintegrals!(zlocal, operator::HH3DSingleLayerTDBIO,
        test_point,
        test_local_space, trial_local_space, time_local_space,
        test_element, trial_element, time_element,
        quad_rule, quad_weight)

    # error("Here!!!")

    dx = quad_weight
    x = cartesian(test_point)
    # n = normal(test_point)

    # a = trial_element[1]
    # ξ = x - dot(x -a, n) * n

    r = time_element[1]
    R = time_element[2]
    @assert r < R

    N = max(degree(time_local_space), 1)
    ∫G, ∫vG, ∫∇G = WiltonInts84.wiltonints(
        trial_element[1],
        trial_element[2],
        trial_element[3],
        x, r, R, Val{2}, quad_rule.workspace)

    a = dx / (4*pi)
    D = operator.num_diffs
    @assert D == 0
    @assert numfunctions(test_local_space)  == 1
    @assert numfunctions(trial_local_space) == 1

    @inline function tmRoR_sl(d, iG)
        sgn = isodd(d) ? -1 : 1
        r = sgn * iG[d+2]
    end

    # bns = quad_rule.binomials

    @assert D == 0
    for k in 1 : numfunctions(time_local_space)
        d = k - 1
        d < D && continue
        q = reduce(*, d-D+1:d ,init=1)
        zlocal[1,1,k] += a * q * tmRoR_sl(d-D, ∫G)
    end # k
end

function innerintegrals!(zlocal, operator::HH3DHyperSingularTDBIO,
        test_point,
        test_local_space, trial_local_space, time_local_space,
        test_element, trial_element, time_element,
        quad_rule, quad_weight)

    dx = quad_weight
    x = cartesian(test_point)

    nx = normal(test_point)
    ny = normal(trial_element)
    ndotn = dot(nx,ny)

    a = trial_element[1]
    ξ = x - dot(x -a, ny) * ny

    r = time_element[1]
    R = time_element[2]
    @assert r < R

    N = max(degree(time_local_space), 1)
    ∫G, ∫Gξy, = WiltonInts84.wiltonints(
        trial_element[1],
        trial_element[2],
        trial_element[3],
        x, r, R, Val{N-1},quad_rule.workspace)

    @assert numfunctions(test_local_space)  <= 3
    @assert numfunctions(trial_local_space) == 3

    @inline function tmRoR(d, iG)
        r = (isodd(d) ? -1 : 1) * iG[d+2]
    end

    @inline function tmRoRf(d, ∫G, ∫Gξy, bξ, h, m)
        (isodd(d) ? -1 : 1) * ((1 - h*dot(m,bξ)) * ∫G[d+2] - h*dot(m, ∫Gξy[d+2]))
    end

    test_values = test_local_space(test_point, Val{:withcurl})
    trial_values = trial_local_space(center(trial_element), Val{:withcurl})

    # weakly singular term
    α = dx / (4π) * operator.weight_of_weakly_singular_term
    Ds = operator.num_diffs_on_weakly_singular_term
    for i in 1 : numfunctions(test_local_space)
        g, curlg = test_values[i]
        for j in 1 : numfunctions(trial_local_space)
            b = trial_element[j]
            opp_edge = trial_element[mod1(j+2,3)] - trial_element[mod1(j+1,3)]
            h = norm(opp_edge)/2/volume(trial_element)
            m = normalize(cross(opp_edge, ny))
            for k in 1 : numfunctions(time_local_space)
                d = k-1
                d < Ds && continue
                q = reduce(*, d-Ds+1:d ,init=1)
                zlocal[i,j,k] += α * q * g * ndotn * tmRoRf(d-Ds, ∫G, ∫Gξy, ξ-b, h, m)
            end
        end
    end

    # Hyper-singular term
    β = dx / (4π) * operator.weight_of_hyper_singular_term
    Dh = operator.num_diffs_on_hyper_singular_term
    for i in 1 : numfunctions(test_local_space)
        g, curlg = test_values[i]
        for j in 1 : numfunctions(trial_local_space)
            _, curlf = trial_values[j]
            for k in 1 : numfunctions(time_local_space)
                d = k - 1
                d < Dh && continue
                q = reduce(*, d-Dh+1:d ,init=1)
                zlocal[i,j,k] += β * q * dot(curlg, curlf) * tmRoR(d-Dh, ∫G)
            end
        end
    end
end

function innerintegrals!(zlocal, operator::HH3DDoubleLayerTDBIO,
        test_point,
        test_local_space, trial_local_space, time_local_space,
        test_element, trial_element, time_element,
        quad_rule, quad_weight)

    @assert numfunctions(test_local_space)  <= 3
    @assert numfunctions(trial_local_space) == 1

    dx = quad_weight
    x = cartesian(test_point)
    nx = normal(test_point)

    # a = trial_element[1]
    # ξ = x - dot(x -a, ny) * ny

    r = time_element[1]
    R = time_element[2]
    @assert r < R

    N = max(degree(time_local_space), 1)
    ∫G, ∫Gξy, ∫∇G = WiltonInts84.wiltonints(
        trial_element[1],
        trial_element[2],
        trial_element[3],
        x, r, R, Val{N-1},quad_rule.workspace)

    @inline function grad_tmRoR(d,∇G)
        σ = isodd(d) ? -1 : 1
        return σ * ∇G[d+1]
    end

    test_values = test_local_space(test_point, Val{:withcurl})
    trial_values = trial_local_space(center(trial_element), Val{:withcurl})

    @assert all(getindex.(trial_values,1) .≈ [1])
    # @assert all(getindex.(trial_values,2) .≈ Ref([0,0,0]))

    # weakly singular term
    α = dx / (4π) * operator.weight
    D = operator.num_diffs
    for i in 1 : numfunctions(test_local_space)
        g, curlg = test_values[i]
        for j in 1 : numfunctions(trial_local_space)
            f, curlf = trial_values[j]
            for k in 1 : numfunctions(time_local_space)
                d = k-1
                d < D && continue
                q = reduce(*, d-D+1:d ,init=1)
                @assert q == 1
                zlocal[i,j,k] += α * q * g * dot(nx, grad_tmRoR(d-D, ∫∇G)) * f
            end
        end
    end

    # @assert all(zlocal[:,:,2] .≈ 0)
end


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using RigidBodyDynamics
using Random
using Profile
using BenchmarkTools

const ScalarType = Float64
# const ScalarType = Float32

function create_floating_atlas()
    url = "https://raw.githubusercontent.com/RobotLocomotion/drake/6e3ca768cbaabf15d0f2bed0fb5bd703fa022aa5/drake/examples/Atlas/urdf/atlas_minimal_contact.urdf"
    urdf = RigidBodyDynamics.cached_download(url, "atlas.urdf")
    atlas = parse_urdf(urdf, scalar_type=ScalarType, floating=true)
    atlas
end

function create_benchmark_suite()
    suite = BenchmarkGroup()
    mechanism = create_floating_atlas()
    remove_fixed_tree_joints!(mechanism)

    state = MechanismState{ScalarType}(mechanism)
    result = DynamicsResult{ScalarType}(mechanism)
    nv = num_velocities(state)
    mat = MomentumMatrix(root_frame(mechanism), Matrix{ScalarType}(undef, 3, nv), Matrix{ScalarType}(undef, 3, nv))
    torques = similar(velocity(state))
    rfoot = findbody(mechanism, "r_foot")
    lhand = findbody(mechanism, "l_hand")
    p = path(mechanism, rfoot, lhand)
    jac = GeometricJacobian(default_frame(lhand), default_frame(rfoot), root_frame(mechanism),
        Matrix{ScalarType}(undef, 3, nv), Matrix{ScalarType}(undef, 3, nv))

    suite["mass_matrix!"] = @benchmarkable(begin
        setdirty!($state)
        mass_matrix!($(result.massmatrix), $state)
    end, setup = rand!($state), evals = 10)

    suite["dynamics_bias!"] = @benchmarkable(begin
        setdirty!($state)
        dynamics_bias!($result, $state)
    end, setup = begin
        rand!($state)
    end, evals = 10)

    suite["inverse_dynamics!"] = @benchmarkable(begin
        setdirty!($state)
        inverse_dynamics!($torques, $(result.jointwrenches), $(result.accelerations), $state, v̇, externalwrenches)
    end, setup = begin
        v̇ = similar(velocity($state))
        rand!(v̇)
        externalwrenches = RigidBodyDynamics.BodyDict(convert(BodyID, body) => rand(Wrench{ScalarType}, root_frame($mechanism)) for body in bodies($mechanism))
        rand!($state)
    end, evals = 10)

    suite["dynamics!"] = @benchmarkable(begin
        setdirty!($state)
        dynamics!($result, $state, τ, externalwrenches)
    end, setup = begin
        rand!($state)
        τ = similar(velocity($state))
        rand!(τ)
        externalwrenches = RigidBodyDynamics.BodyDict(convert(BodyID, body) => rand(Wrench{ScalarType}, root_frame($mechanism)) for body in bodies($mechanism))
    end, evals = 10)

    suite["momentum_matrix!"] = @benchmarkable(begin
        setdirty!($state)
        momentum_matrix!($mat, $state)
    end, setup = rand!($state), evals = 10)

    suite["geometric_jacobian!"] = @benchmarkable(begin
        setdirty!($state)
        geometric_jacobian!($jac, $state, $p)
    end, setup = rand!($state), evals = 10)

    suite["mass_matrix! and geometric_jacobian!"] = @benchmarkable(begin
        setdirty!($state)
        mass_matrix!($(result.massmatrix), $state)
        geometric_jacobian!($jac, $state, $p)
    end, setup = begin
        rand!($state)
    end, evals = 10)

    suite["momentum"] = @benchmarkable(begin
        setdirty!($state)
        momentum($state)
    end, setup = rand!($state), evals = 10)

    suite["momentum_rate_bias"] = @benchmarkable(begin
        setdirty!($state)
        momentum_rate_bias($state)
    end, setup = rand!($state), evals = 10)

    suite["kinetic_energy"] = @benchmarkable(begin
        setdirty!($state)
        kinetic_energy($state)
    end, setup = rand!($state), evals = 10)

    suite["gravitational_potential_energy"] = @benchmarkable(begin
        setdirty!($state)
        gravitational_potential_energy($state)
    end, setup = rand!($state), evals = 10)

    suite["center_of_mass"] = @benchmarkable(begin
        setdirty!($state)
        center_of_mass($state)
    end, setup = rand!($state), evals = 10)

    mcmechanism = maximal_coordinates(mechanism)
    mcstate = MechanismState{ScalarType}(mcmechanism)
    mcresult = DynamicsResult{ScalarType}(mcmechanism)

    suite["constraint_jacobian!"] = @benchmarkable(begin
        setdirty!($mcstate)
        RigidBodyDynamics.constraint_jacobian!($(mcresult.constraintjacobian), $(mcresult.constraintrowranges), $mcstate)
    end, setup = rand!($mcstate), evals = 10)

    suite["constraint_bias!"] = @benchmarkable(begin
        setdirty!($mcstate)
        RigidBodyDynamics.constraint_bias!($(mcresult.constraintbias), $mcstate)
    end, setup = rand!($mcstate), evals = 10)

    suite
end

function runbenchmarks()
    suite = create_benchmark_suite()
    Profile.clear_malloc_data()
    overhead = BenchmarkTools.estimate_overhead()
    Random.seed!(1)
    results = run(suite, verbose=true, overhead=overhead, gctrial=false)
    for result in results
        println("$(first(result)):")
        display(last(result))
        println()
    end
end

runbenchmarks()


text:
push!(LOAD_PATH,".")
import cosmosis
block=cosmosis.make_datablock()
cosmosis.put_int(block, "my_section", "my_int", 17)
x = cosmosis.get_int(block, "my_section", "my_int")

cosmosis.put_double(block, "my_section", "my_double", 12.99)
y = cosmosis.get_double(block, "my_section", "my_double")

cosmosis.put_bool(block, "my_section", "my_bool", 1)
z = cosmosis.get_bool(block, "my_section", "my_bool")

cosmosis.put_string(block, "my_section", "my_str", "hat")
w = cosmosis.get_string(block, "my_section", "my_str")

a = fill(19.4, 2, 3)
cosmosis.put_double_array(block, "my_section", "myarr", a)
b= cosmosis.get_double_array(block, "my_section", "myarr")

c = fill(Int32(19), 2, 3)
cosmosis.put_int_array(block, "my_section", "myarr2", c)
d = cosmosis.get_int_array(block, "my_section", "myarr2")

println("x = ", x)
println("y = ", y)
println("z = ", z)
println("w = ", w)
println("b = ", b)
println("d = ", d)

max_stars_repo_name: annis/cosmosis

max_issues_repo_name: annis/cosmosis

max_forks_repo_name: annis/cosmosis


text:
"""
```
t = find_threshold(Unimodal(), histogram, edges)
```

Generates a threshold value for array `histogram` and interval `edges`
using Rosin's algorithm.

# Output
Returns `t`, a real number that specifies the image's threshold using `edges`,
an interval that specifies how the histogram is divided into bins, and `histogram`,
an array that contains the frequencies of each bin.

# Details
This algorithm first selects the bin in the histogram with the highest frequency.
The algorithm then searches from the location of the maximum bin to the last bin
of the histogram for the first bin with a frequency of 0 (known as the minimum bin.).
A line is then drawn that passes through both the maximum and minimum bins. The
bin with the greatest orthogonal distance to the line is chosen as the threshold
value.

## Assumptions
This algorithm assumes that:

* The histogram is unimodal
* There is always at least one bin that has a frequency of 0. If not, the
algorithm will use the last bin as the minimum bin.
* If the histogram includes multiple bins with a frequency of 0, the algorithm
will select the first zero bin as its minimum.
* If there are multiple bins with the greatest orthogonal distance, the leftmost
bin is selected as the threshold.

# Example

Compute the threshold for the "moonsurface" image in the TestImages package.
```julia

using TestImages, HistogramThresholding

img = testimage("moonsurface")
edges, counts = build_histogram(img,256)
t = find_threshold(UniThreshold(),counts, edges)
```

# Reference
[1] P. L. Rosin, “Unimodal thresholding,” Pattern Recognition, vol. 34, no. 11, pp. 2083–2096, Nov. 2001.

# See Also
"""
function find_threshold(algorithm::UniThreshold, histogram::AbstractArray, edges::AbstractRange)
    #= Calculates the orthogonal distance between the line (slope m,
    y-intercept c), and the point (co-ordinates x₀, y₀).
    =#
    function calculate_distance(m::Real, c::Real, x₀::Real, y₀::Real)
        abs(m*x₀ + -1*y₀ + c) / sqrt((m)^2 + 1)
    end

    # Finds the index & value of bin with highest value.
    max_val, max_index = findmax(histogram)

    # Initialise minimum bin index.
    min_index = 1
    found_min = false

    # Find the first zero bin.
    for j in max_index:last(first(axes(histogram)))
        if histogram[j] == 0.0
            min_index = j
            found_min = true
            break
        end
    end

    # If no zero bin is found, use the last bin.
    if found_min == false
        min_index = last(first(axes(histogram)))
    end

    # Find slope and y-intercept of the line.
    m = max_val / (max_index - min_index)
    c = -m*min_index

    # Initialise the maximum bin as having the inital greatest distance.
    t = max_index
    max_distance = calculate_distance(m, c, t, histogram[t])

    # Find bin with the greatest distance.
    for k in max_index:min_index
        if calculate_distance(m, c, k, histogram[k]) > max_distance
            max_distance = calculate_distance(m, c, k, histogram[k])
            t = k
        end
    end

    # Return the greatest distance.
    return edges[t]
end


text:
using Oceananigans.Utils: prettytime, ordered_dict_show
using Oceananigans: short_show

"""Show the innards of a `Model` in the REPL."""
function Base.show(io::IO, model::IncompressibleModel{TS, C, A}) where {TS, C, A}
    print(io, "IncompressibleModel{$A, $(eltype(model.grid))}",
        "(time = $(prettytime(model.clock.time)), iteration = $(model.clock.iteration)) \n",
        "├── grid: $(short_show(model.grid))\n",
        "├── tracers: $(tracernames(model.tracers))\n",
        "├── closure: $(typeof(model.closure))\n",
        "├── buoyancy: $(typeof(model.buoyancy))\n")

    if isnothing(model.particles)
        print(io, "└── coriolis: $(typeof(model.coriolis))")
    else
        particles = model.particles.properties
        properties = propertynames(particles)
        print(io, "├── coriolis: $(typeof(model.coriolis))\n")
        print(io, "└── particles: $(length(particles)) Lagrangian particles with $(length(properties)) properties: $properties")
    end
end


text:
import GitLab
using Base.Test

myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

myrepo = GitLab.repo_by_name("TestProject1"; headers=options)
## @show myrepo

statuses = GitLab.statuses(myrepo, "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8"; params=options)
@test get(first(statuses)[1].sha) == "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8"
## @show statuses

status_data = Dict{AbstractString, Any}()
status_data["state"] = "running"
status = GitLab.create_status(myrepo, "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8"; headers=options, params=status_data)
@test get(status.sha) == "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8"
## @show status

status_data["state"] = "invalid"
try
    status = GitLab.create_status(myrepo, "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8"; headers=options, params=status_data)
    @test 1 == 0
catch e
    println("Failed (Expected)  to update status : $(e)")
end
println("Statuses Tests Done !!!")