text:
using Plots
gr()
import AoC2019
input_data = AoC2019.read_numbers("input15.txt")

mutable struct BoardState
  dims
  labyrinth::Array{Int8, 2}
  current_position::Tuple{Int16, Int16}
end

function draw_state(state)
  color_map = Dict(0 => ' ', 1 => '.', 2 => 'S', 3 => 'w', 4 => 'D')
  println(AoC2019.visualizeArray(state.labyrinth, color_map))
end

function move_robot!(state, command, result)
  println("moving robot")
  if command == 1
    position_aimed = state.current_position .+ (-1, 0)
  elseif command == 2
    position_aimed = state.current_position .+ (1, 0)
  elseif command == 3
    position_aimed = state.current_position .+ (0, -1)
  elseif command == 4
    position_aimed = state.current_position .+ (0, 1)
  else
    error("Not implemented")
  end
  if result == 0
    state.labyrinth[position_aimed...] = 3
    # say that this way is bad
  elseif result == 1 || result == 2
    # take over state of previous way for previous direction
    state.labyrinth[state.current_position...] = 1
    state.current_position = position_aimed
    if result == 2
      state.labyrinth[position_aimed...] = 2
    end
  end
end

function generate_directions!(state, stack, steps_done)
  reverse_direction = Dict(1 => 2, 2 => 1, 3 => 4, 4 => 3)
  for i in 1:4
    if i == 1
      position_aimed = state.current_position .+ (-1, 0)
    elseif i == 2
      position_aimed = state.current_position .+ (1, 0)
    elseif i == 3
      position_aimed = state.current_position .+ (0, -1)
    elseif i == 4
      position_aimed = state.current_position .+ (0, 1)
    end
    if state.labyrinth[position_aimed...] == 1
      continue
    end
    push!(stack, (reverse_direction[i], steps_done))
    push!(stack, (i, steps_done + 1))
  end
  state.labyrinth[state.current_position...] = 1

end

function robot_searcher(input_data)

  robot_input = Channel{Int64}()
  robot_output = Channel{Int64}()
  #dims = (24, 42)
  dims = (50, 50)
  labyrinth = zeros(Int8, dims...)
  current_score = 0
  current_position = dims .÷ 2
  labyrinth[current_position...] = 4

  state = BoardState(dims, labyrinth, dims .÷ 2)

  #commands = [(1,0), (4,1), (1,0), (2,0), (4,0), (3,1), (3,0), (2,1), (2,0), (3,2)]
  #  for command in commands
  #    move_robot(state, command...)
  #  end
  #  @sync begin
  #    draw_state(state)
  #  end

  stack::Array{Tuple{Int8, Int32}} = []
  generate_directions!(state, stack, 0)
  println("Let's start")
  position_goal = (0,0)
  steps_needed = -1
  @sync begin
    @async AoC2019.run_program_async(input_data, robot_input, robot_output)

    @async while true
      try
        # inputting random direction
        if length(stack) == 0
          state.labyrinth[position_goal...] = 2
          close(robot_input)
          break
        end
        command = pop!(stack)
        direction = command[1]
        steps_done = command[2]
        put!(robot_input, direction)
        # get result
        result = take!(robot_output)
        if result == 0
          #didn't move, remove backtrace as well
          pop!(stack)
        end
        move_robot!(state, direction, result)
        if state.labyrinth[state.current_position...] == 2
          #close(robot_input)
          position_goal = state.current_position
          steps_needed = steps_done
        end
        if state.labyrinth[state.current_position...] == 0
          generate_directions!(state, stack, steps_done)
        end

        state.labyrinth[state.current_position...] = 4

        println("Position of goal: $position_goal, steps_needed = $steps_needed")
        plt = heatmap(state.labyrinth, yflip=true)
        gui(plt)
        #draw_state(state)
      catch err
        close(robot_input)
        close(robot_output)
        rethrow(err)
      end
    end
  end
  return state.labyrinth
end

labyrinth = robot_searcher(input_data)
#oxygen_generator = findfirst(x -> x == 2, labyrinth)
#status = AoC2019.bfs_depth(labyrinth, oxygen_generator, x->!isequal(x, 3))

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
#-----------------------------------------------------------------------
#   CHECK MINKOWSKI METRIC FUNCTIONS USED HERE
#-----------------------------------------------------------------------

tv  = [1.;0.;0.;0.]
xv  = [0.;1.;0.;0.]
yv  = [0.;0.;1.;0.]
zv  = [0.;0.;0.;1.]

@test ηdot(tv,tv) ≈ -1
@test ηdot(xv,xv) ≈ 1
@test ηdot(yv,yv) ≈ 1
@test ηdot(zv,zv) ≈ 1

@test ηdot(tv,xv) ≈ 0
@test ηdot(tv,yv) ≈ 0
@test ηdot(tv,zv) ≈ 0
@test ηdot(xv,yv) ≈ 0
@test ηdot(xv,zv) ≈ 0
@test ηdot(yv,zv) ≈ 0

k1  = π^2*[1.;1.;0.;0.]
k2  = π^2*[1.;0.;1.;0.]
k3  = π^2*[1.;0.;0.;1.]

kn1 = π^2*[1.;-1.;0.;0.]
kn2 = π^2*[1.;0.;-1.;0.]
kn3 = π^2*[1.;0.;0.;-1.]

@test ηdot(k1,k1) ≈ 0
@test ηdot(k2,k2) ≈ 0
@test ηdot(k3,k3) ≈ 0

@test ηdot(kn1,kn1) ≈ 0
@test ηdot(kn2,kn2) ≈ 0
@test ηdot(kn3,kn3) ≈ 0

@test ημν()[1,1] ≈ -1
@test ημν()[2,2] ≈ 1
@test ημν()[3,3] ≈ 1
@test ημν()[4,4] ≈ 1

@test ημν()[1,2] ≈ 0
@test ημν()[1,3] ≈ 0
@test ημν()[1,4] ≈ 0
@test ημν()[2,3] ≈ 0
@test ημν()[2,4] ≈ 0
@test ημν()[3,4] ≈ 0

ext: jl


text:
using Serialization
using FreeTypeAbstraction: iter_or_array

mutable struct TextureAtlas
    rectangle_packer::RectanglePacker
    mapping::Dict{Tuple{Char, String}, Int} # styled glyph to index in sprite_attributes
    index::Int
    data::Matrix{Float16}
    # rectangles we rendered our glyphs into in normalized uv coordinates
    uv_rectangles::Vector{Vec4f}
end

Base.size(atlas::TextureAtlas) = size(atlas.data)

@enum GlyphResolution High Low

const HIGH_PIXELSIZE = 64
const LOW_PIXELSIZE = 32

const CACHE_RESOLUTION_PREFIX = Ref("high")
const TEXTURE_RESOLUTION = Ref((2048, 2048))
const PIXELSIZE_IN_ATLAS = Ref(HIGH_PIXELSIZE)
# Small padding results in artifacts during downsampling. It seems like roughly
# 1.5px padding is required for a scaled glyph to be displayed without artifacts.
# E.g. for textsize = 8px we need 1.5/8 * 64 = 12px+ padding (in each direction)
# for things to look clear with a 64px glyph size.
const GLYPH_PADDING = Ref(12)

function set_glyph_resolution!(res::GlyphResolution)
    if res == High
        TEXTURE_RESOLUTION[] = (2048, 2048)
        CACHE_RESOLUTION_PREFIX[] = "high"
        PIXELSIZE_IN_ATLAS[] = HIGH_PIXELSIZE
        GLYPH_PADDING[] = 12
    else
        TEXTURE_RESOLUTION[] = (1024, 1024)
        CACHE_RESOLUTION_PREFIX[] = "low"
        PIXELSIZE_IN_ATLAS[] = LOW_PIXELSIZE
        GLYPH_PADDING[] = 6
    end
end

function TextureAtlas(initial_size = TEXTURE_RESOLUTION[])
    return TextureAtlas(
        RectanglePacker(Rect2(0, 0, initial_size...)),
        Dict{Tuple{Char, String}, Int}(),
        1,
        zeros(Float16, initial_size...),
        Vec4f[],
    )
end

begin
    # basically a singleton for the textureatlas

    function get_cache_path()
        return abspath(
            first(Base.DEPOT_PATH), "makiegallery", ".cache",
            "texture_atlas_$(CACHE_RESOLUTION_PREFIX[])_$(VERSION).jls"
        )
    end

    const _default_font = NativeFont[]
    const _alternative_fonts = NativeFont[]

    function defaultfont()
        if isempty(_default_font)
            push!(_default_font, to_font("Dejavu Sans"))
        end
        _default_font[]
    end

    function alternativefonts()
        if isempty(_alternative_fonts)
            alternatives = [
                "DejaVuSans.ttf",
                "NotoSansCJKkr-Regular.otf",
                "NotoSansCuneiform-Regular.ttf",
                "NotoSansSymbols-Regular.ttf",
                "FiraMono-Medium.ttf"
            ]
            for font in alternatives
                push!(_alternative_fonts, NativeFont(assetpath("fonts", font)))
            end
        end
        return _alternative_fonts
    end

    function load_ascii_chars!(atlas)
        # for c in '\u0000':'\u00ff' #make sure all ascii is mapped linearly
        #     insert_glyph!(atlas, c, defaultfont())
        # end
        for c in '0':'9' #make sure all ascii is mapped linearly
            insert_glyph!(atlas, c, defaultfont())
        end
        for c in 'A':'Z' #make sure all ascii is mapped linearly
            insert_glyph!(atlas, c, defaultfont())
        end
    end

    function cached_load()
        if isfile(get_cache_path())
            try
                return open(get_cache_path()) do io
                    dict = Serialization.deserialize(io)
                    fields = map(fieldnames(TextureAtlas)) do n
                        v = dict[n]
                        isa(v, Vector) ? copy(v) : v # otherwise there seems to be a problem with resizing
                    end
                    TextureAtlas(fields...)
                end
            catch e
                @info("You can likely ignore the following warning, if you just switched Julia versions for Makie")
                @warn(e)
                rm(get_cache_path())
            end
        end
        atlas = TextureAtlas()
        @info("Makie/Makie is caching fonts, this may take a while. Needed only on first run!")
        load_ascii_chars!(atlas)
        to_cache(atlas) # cache it
        return atlas
    end

    function to_cache(atlas)
        if !ispath(dirname(get_cache_path()))
            mkpath(dirname(get_cache_path()))
        end
        open(get_cache_path(), "w") do io
            dict = Dict(map(fieldnames(typeof(atlas))) do name
                name => getfield(atlas, name)
            end)
            Serialization.serialize(io, dict)
        end
    end

    const global_texture_atlas = Dict{Int, TextureAtlas}()

    function get_texture_atlas()
        return get!(global_texture_atlas, PIXELSIZE_IN_ATLAS[]) do
            return cached_load() # initialize only on demand
        end
    end
end

"""
    find_font_for_char(c::Char, font::NativeFont)
Finds the best font for a character from a list of fallback fonts, that get chosen
if `font` can't represent char `c`
"""
function find_font_for_char(c::Char, font::NativeFont)
    FT_Get_Char_Index(font, c) != 0 && return font
    # it seems that linebreaks are not found which messes up font metrics
    # if another font is selected just for those chars
    c in ('\n', '\r', '\t') && return font
    for afont in alternativefonts()
        if FT_Get_Char_Index(afont, c) != 0
            return afont
        end
    end
    error("Can't represent character $(c) with any fallback font nor $(font.family_name)!")
end

function glyph_index!(atlas::TextureAtlas, c::Char, font::NativeFont)
    if FT_Get_Char_Index(font, c) == 0
        for afont in alternativefonts()
            if FT_Get_Char_Index(afont, c) != 0
                font = afont
            end
        end
    end
    return insert_glyph!(atlas, c, font)
end

function glyph_uv_width!(atlas::TextureAtlas, c::Char, font::NativeFont)
    return atlas.uv_rectangles[glyph_index!(atlas, c, font)]
end

function glyph_uv_width!(c::Char)
    return glyph_uv_width!(get_texture_atlas(), c, defaultfont())
end

# function glyph_boundingbox(c::Char, font::NativeFont, pixelsize)
#     if FT_Get_Char_Index(font, c) == 0
#         for afont in alternativefonts()
#             if FT_Get_Char_Index(afont, c) != 0
#                 font = afont
#                 break
#             end
#         end
#     end
#     bb, ext = FreeTypeAbstraction.metrics_bb(c, font, pixelsize)
#     return bb
# end

function insert_glyph!(atlas::TextureAtlas, glyph::Char, font::NativeFont)
    return get!(atlas.mapping, (glyph, FreeTypeAbstraction.fontname(font))) do
        # We save glyphs as signed distance fields, i.e. we save the distance
        # a pixel is away from the edge of a symbol (continuous at the edge).
        # To get accurate distances we want to draw the symbol at high
        # resolution and then downsample to the PIXELSIZE_IN_ATLAS.
        downsample = 5
        # To draw a symbol from a sdf we essentially do `color * (sdf > 0)`. For
        # antialiasing we smooth out the step function `sdf > 0`. That means we
        # need a few values outside the symbol. To guarantee that we have those
        # at all relevant scales we add padding to the rendered bitmap and the
        # resulting sdf.
        pad = GLYPH_PADDING[]

        uv_pixel = render(atlas, glyph, font, downsample, pad)
        tex_size = Vec2f(size(atlas.data) .- 1) # starts at 1

        # 0 based
        idx_left_bottom = minimum(uv_pixel)
        idx_right_top = maximum(uv_pixel)

        # transform to normalized texture coordinates
        # -1 for indexing offset
        uv_left_bottom_pad = (idx_left_bottom) ./ tex_size
        uv_right_top_pad = (idx_right_top .- 1) ./ tex_size

        uv_offset_rect = Vec4f(uv_left_bottom_pad..., uv_right_top_pad...)
        i = atlas.index
        push!(atlas.uv_rectangles, uv_offset_rect)
        atlas.index = i + 1
        return i
    end
end

"""
    sdistancefield(img, downsample, pad)
Calculates a distance fields, that is downsampled `downsample` time,
with a padding applied of `pad`.
The padding is in units after downscaling!
"""
function sdistancefield(img, downsample, pad)
    # we pad before downsampling, so we need to have `downsample` as much padding
    pad = downsample * pad
    # padd the image
    padded_size = size(img) .+ 2pad

    # for the downsampling, we need to make sure that
    # we can divide the image size by `downsample` without reminder
    dividable_size = ceil.(Int, padded_size ./ downsample) .* downsample

    in_or_out = fill(false, dividable_size)
    # the size we fill the image up to
    wend, hend = size(img) .+ pad
    in_or_out[pad+1:wend, pad+1:hend] .= img .> (0.5 * 255)

    yres, xres = dividable_size .÷ downsample
    # divide by downsample to normalize distances!
    return Float16.(sdf(in_or_out, xres, yres) ./ downsample)
end

const font_render_callbacks = Dict{Int, Vector{Function}}()

function font_render_callback!(f)
    funcs = get!(font_render_callbacks, PIXELSIZE_IN_ATLAS[], Function[])
    push!(funcs, f)
end

function remove_font_render_callback!(f)
    for (s, callbacks) in font_render_callbacks
        filter!(f2-> f2 != f, callbacks)
    end
end

function render(atlas::TextureAtlas, glyph::Char, font, downsample=5, pad=6)
    #select_font_face(cc, font)
    if glyph == '\n' # don't render  newline
        glyph = ' '
    end
    # the target pixel size of our distance field
    pixelsize = PIXELSIZE_IN_ATLAS[]
    # we render the font `downsample` sizes times bigger
    # Make sure the font doesn't have a mutated font matrix from e.g. Cairo
    FreeTypeAbstraction.FreeType.FT_Set_Transform(font, C_NULL, C_NULL)
    bitmap, extent = renderface(font, glyph, pixelsize * downsample)
    # Our downsampeld & padded distancefield
    sd = sdistancefield(bitmap, downsample, pad)
    rect = Rect2(0, 0, size(sd)...)
    uv = push!(atlas.rectangle_packer, rect) # find out where to place the rectangle
    uv == nothing && error("texture atlas is too small. Resizing not implemented yet. Please file an issue at Makie if you encounter this") #TODO resize surface
    # write distancefield into texture
    atlas.data[uv.area] = sd
    for f in get(font_render_callbacks, pixelsize, ())
        # update everyone who uses the atlas image directly (e.g. in GLMakie)
        f(sd, uv.area)
    end
    # return the area we rendered into!
    return uv.area
end


text:
# Julia wrapper for header: include/declarations.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0 and then modified as some places

export read_prob

function triu(A::blockmatrix)
    ccall((:triu,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function store_packed(A::blockmatrix,B::blockmatrix)
    ccall((:store_packed,CSDP.csdp),Nothing,(blockmatrix,blockmatrix),A,B)
end

function store_unpacked(A::blockmatrix,B::blockmatrix)
    ccall((:store_unpacked,CSDP.csdp),Nothing,(blockmatrix,blockmatrix),A,B)
end

function alloc_mat_packed(A::blockmatrix,pB::Ref{blockmatrix})
    ccall((:alloc_mat_packed,CSDP.csdp),Nothing,(blockmatrix,Ref{blockmatrix}),A,pB)
end

function free_mat_packed(A::blockmatrix)
    ccall((:free_mat_packed,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function structnnz(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,constraints::Ptr{constraintmatrix})
    ccall((:structnnz,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{constraintmatrix}),n,k,C,constraints)
end

function actnnz(n::CSDP_INT,lda::CSDP_INT,A::Ptr{Cdouble})
    ccall((:actnnz,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,Ptr{Cdouble}),n,lda,A)
end

function bandwidth(n::CSDP_INT,lda::CSDP_INT,A::Ptr{Cdouble})
    ccall((:bandwidth,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,Ptr{Cdouble}),n,lda,A)
end

function qreig(n::CSDP_INT,maindiag::Ptr{Cdouble},offdiag::Ptr{Cdouble})
    ccall((:qreig,CSDP.csdp),Nothing,(CSDP_INT,Ptr{Cdouble},Ptr{Cdouble}),n,maindiag,offdiag)
end

function sort_entries(k::CSDP_INT,C::blockmatrix,constraints::Ptr{constraintmatrix})
    ccall((:sort_entries,CSDP.csdp),Nothing,(CSDP_INT,blockmatrix,Ptr{constraintmatrix}),k,C,constraints)
end

function norm2(n::CSDP_INT,x::Ptr{Cdouble})
    ccall((:norm2,CSDP.csdp),Cdouble,(CSDP_INT,Ptr{Cdouble}),n,x)
end

function norm1(n::CSDP_INT,x::Ptr{Cdouble})
    ccall((:norm1,CSDP.csdp),Cdouble,(CSDP_INT,Ptr{Cdouble}),n,x)
end

function norminf(n::CSDP_INT,x::Ptr{Cdouble})
    ccall((:norminf,CSDP.csdp),Cdouble,(CSDP_INT,Ptr{Cdouble}),n,x)
end

function Fnorm(A::blockmatrix)
    ccall((:Fnorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function Knorm(A::blockmatrix)
    ccall((:Knorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function mat1norm(A::blockmatrix)
    ccall((:mat1norm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function matinfnorm(A::blockmatrix)
    ccall((:matinfnorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function calc_pobj(C::blockmatrix,X::blockmatrix,constant_offset::Cdouble)
    ccall((:calc_pobj,CSDP.csdp),Cdouble,(blockmatrix,blockmatrix,Cdouble),C,X,constant_offset)
end

function calc_dobj(k::CSDP_INT,a::Ptr{Cdouble},y::Ptr{Cdouble},constant_offset::Cdouble)
    ccall((:calc_dobj,CSDP.csdp),Cdouble,(CSDP_INT,Ptr{Cdouble},Ptr{Cdouble},Cdouble),k,a,y,constant_offset)
end

function trace_prod(A::blockmatrix,B::blockmatrix)
    ccall((:trace_prod,CSDP.csdp),Cdouble,(blockmatrix,blockmatrix),A,B)
end

function linesearch(n::CSDP_INT,dX::blockmatrix,work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,cholinv::blockmatrix,q::Ptr{Cdouble},z::Ptr{Cdouble},workvec::Ptr{Cdouble},stepfrac::Cdouble,start::Cdouble,printlevel::CSDP_INT)
    ccall((:linesearch,CSDP.csdp),Cdouble,(CSDP_INT,blockmatrix,blockmatrix,blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,CSDP_INT),n,dX,work1,work2,work3,cholinv,q,z,workvec,stepfrac,start,printlevel)
end

function pinfeas(k::CSDP_INT,constraints::Ptr{constraintmatrix},X::blockmatrix,a::Ptr{Cdouble},workvec::Ptr{Cdouble})
    ccall((:pinfeas,CSDP.csdp),Cdouble,(CSDP_INT,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},Ptr{Cdouble}),k,constraints,X,a,workvec)
end

function dinfeas(k::CSDP_INT,C::blockmatrix,constraints::Ptr{constraintmatrix},y::Ptr{Cdouble},Z::blockmatrix,work1::blockmatrix)
    ccall((:dinfeas,CSDP.csdp),Cdouble,(CSDP_INT,blockmatrix,Ptr{constraintmatrix},Ptr{Cdouble},blockmatrix,blockmatrix),k,C,constraints,y,Z,work1)
end

function dimacserr3(k::CSDP_INT,C::blockmatrix,constraints::Ptr{constraintmatrix},y::Ptr{Cdouble},Z::blockmatrix,work1::blockmatrix)
    ccall((:dimacserr3,CSDP.csdp),Cdouble,(CSDP_INT,blockmatrix,Ptr{constraintmatrix},Ptr{Cdouble},blockmatrix,blockmatrix),k,C,constraints,y,Z,work1)
end

function op_a(k::CSDP_INT,constraints::Ptr{constraintmatrix},X::blockmatrix,result::Ptr{Cdouble})
    ccall((:op_a,CSDP.csdp),Nothing,(CSDP_INT,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble}),k,constraints,X,result)
end

function op_at(k::CSDP_INT,y::Ptr{Cdouble},constraints::Ptr{constraintmatrix},result::blockmatrix)
    ccall((:op_at,CSDP.csdp),Nothing,(CSDP_INT,Ptr{Cdouble},Ptr{constraintmatrix},blockmatrix),k,y,constraints,result)
end

function makefill(k::CSDP_INT,C::blockmatrix,constraints::Ptr{constraintmatrix},pfill::Ref{constraintmatrix},work1::blockmatrix,printlevel::CSDP_INT)
    ccall((:makefill,CSDP.csdp),Nothing,(CSDP_INT,blockmatrix,Ptr{constraintmatrix},Ref{constraintmatrix},blockmatrix,CSDP_INT),k,C,constraints,pfill,work1,printlevel)
end

function op_o(k::CSDP_INT,constraints::Ptr{constraintmatrix},byblocks::Ptr{Ptr{sparseblock}},Zi::blockmatrix,X::blockmatrix,O::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix)
    ccall((:op_o,CSDP.csdp),Nothing,(CSDP_INT,Ptr{constraintmatrix},Ptr{Ptr{sparseblock}},blockmatrix,blockmatrix,Ptr{Cdouble},blockmatrix,blockmatrix),k,constraints,byblocks,Zi,X,O,work1,work2)
end

function addscaledmat(A::blockmatrix,scale::Cdouble,B::blockmatrix,C::blockmatrix)
    ccall((:addscaledmat,CSDP.csdp),Nothing,(blockmatrix,Cdouble,blockmatrix,blockmatrix),A,scale,B,C)
end

function zero_mat(A::blockmatrix)
    ccall((:zero_mat,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function add_mat(A::blockmatrix,B::blockmatrix)
    ccall((:add_mat,CSDP.csdp),Nothing,(blockmatrix,blockmatrix),A,B)
end

function sym_mat(A::blockmatrix)
    ccall((:sym_mat,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function make_i(A::blockmatrix)
    ccall((:make_i,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function copy_mat(A::blockmatrix,B::blockmatrix)
    ccall((:copy_mat,CSDP.csdp),Nothing,(blockmatrix,blockmatrix),A,B)
end

function mat_mult(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix)
    ccall((:mat_mult,CSDP.csdp),Nothing,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix),scale1,scale2,A,B,C)
end

function mat_multspa(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspa,CSDP.csdp),Nothing,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_multspb(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspb,CSDP.csdp),Nothing,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_multspc(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspc,CSDP.csdp),Nothing,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_mult_raw(n::CSDP_INT,scale1::Cdouble,scale2::Cdouble,ap::Ptr{Cdouble},bp::Ptr{Cdouble},cp::Ptr{Cdouble})
    ccall((:mat_mult_raw,CSDP.csdp),Nothing,(CSDP_INT,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),n,scale1,scale2,ap,bp,cp)
end

function mat_mult_rawatlas(n::CSDP_INT,scale1::Cdouble,scale2::Cdouble,ap::Ptr{Cdouble},bp::Ptr{Cdouble},cp::Ptr{Cdouble})
    ccall((:mat_mult_rawatlas,CSDP.csdp),Nothing,(CSDP_INT,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),n,scale1,scale2,ap,bp,cp)
end

function matvec(A::blockmatrix,x::Ptr{Cdouble},y::Ptr{Cdouble})
    ccall((:matvec,CSDP.csdp),Nothing,(blockmatrix,Ptr{Cdouble},Ptr{Cdouble}),A,x,y)
end

function alloc_mat(A::blockmatrix,pB::Ref{blockmatrix})
    ccall((:alloc_mat,CSDP.csdp),Nothing,(blockmatrix,Ref{blockmatrix}),A,pB)
end

function free_mat(A::blockmatrix)
    ccall((:free_mat,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function initparams(params::Ptr{paramstruc},pprintlevel::Ptr{CSDP_INT})
    ccall((:initparams,CSDP.csdp),Nothing,(Ptr{paramstruc},Ptr{CSDP_INT}),params,pprintlevel)
end

function loaded_initsoln(problem::Ptr{Cvoid}, X::Ref{blockmatrix}, Z::Ref{blockmatrix})
    y = Ref{Ptr{Cdouble}}(C_NULL)
    ccall((:loaded_initsoln,CSDP.csdp),Nothing,(Ptr{Cvoid},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix}),problem,X,y,Z)
    return y[]
end
function initsoln(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix})
    X = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    y = Ref{Ptr{Cdouble}}(C_NULL)
    Z = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    ccall((:initsoln,CSDP.csdp),Nothing,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix}),n,k,C,a,constraints,X,y,Z)
    X[], y[], Z[]
end

function trans(A::blockmatrix)
    ccall((:trans,CSDP.csdp),Nothing,(blockmatrix,),A)
end

function chol_inv(A::blockmatrix,B::blockmatrix)
    ccall((:chol_inv,CSDP.csdp),Nothing,(blockmatrix,blockmatrix),A,B)
end

function chol(A::blockmatrix)
    ccall((:chol,CSDP.csdp),CSDP_INT,(blockmatrix,),A)
end

function solvesys(m::CSDP_INT,ldam::CSDP_INT,A::Ptr{Cdouble},rhs::Ptr{Cdouble})
    ccall((:solvesys,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,Ptr{Cdouble},Ptr{Cdouble}),m,ldam,A,rhs)
end

function user_exit(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,a::Ptr{Cdouble},dobj::Cdouble,pobj::Cdouble,constant_offset::Cdouble,constraints::Ptr{constraintmatrix},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix,params::paramstruc)
    ccall((:user_exit,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Cdouble,Cdouble,Cdouble,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},blockmatrix,paramstruc),n,k,C,a,dobj,pobj,constant_offset,constraints,X,y,Z,params)
end

function read_sol(fname::Ptr{UInt8},n::CSDP_INT,k::CSDP_INT,C::blockmatrix,pX::Ptr{blockmatrix},py::Ptr{Ptr{Cdouble}},pZ::Ptr{blockmatrix})
    ccall((:read_sol,CSDP.csdp),CSDP_INT,(Ptr{UInt8},CSDP_INT,CSDP_INT,blockmatrix,Ptr{blockmatrix},Ptr{Ptr{Cdouble}},Ptr{blockmatrix}),fname,n,k,C,pX,py,pZ)
end

function load_prob_from_file(fname::String,C::Ref{blockmatrix},printlevel::Integer=1)
    problem = Ref{Ptr{Cvoid}}(C_NULL)
    ret = ccall((:load_prob_from_file,CSDP.csdp),CSDP_INT,(Ptr{UInt8},Ref{blockmatrix},Ref{Ptr{Cvoid}},CSDP_INT),fname,C,problem,printlevel)
    if !iszero(ret)
        error("`CSDP.load_prob_from_file` failed.")
    end
    return LoadingProblem(problem[])
end
function read_prob(fname::String,printlevel::Integer=1)
    n = Ref{CSDP_INT}(0)
    k = Ref{CSDP_INT}(0)
    C = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    a = Ref{Ptr{Cdouble}}(C_NULL)
    constraints = Ref{Ptr{constraintmatrix}}(C_NULL)
    ccall((:read_prob,CSDP.csdp),CSDP_INT,(Ptr{UInt8},Ref{CSDP_INT},Ref{CSDP_INT},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{Ptr{constraintmatrix}},CSDP_INT),fname,n,k,C,a,constraints,printlevel)
    C = mywrap(C[])
    a = mywrap(a[], k[])
    constraints = mywrap(constraints[], k[])
    # Nothing to free as the constraintmatrix is immutable
    # the array is allocated as an array of bitstype in CSDP
    # the array of blocks, .. inside should be free'd though
    #constraints = map(c->ConstraintMatrix(c, k[]), constraints)
    @assert n[] == size(C, 1)
    C, a, constraints
end

function write_prob(fname::String,
                    n::CSDP_INT,
                    k::CSDP_INT,
                    C::blockmatrix,
                    a::Ptr{Cdouble},
                    constraints::Ptr{constraintmatrix})
    ccall((:write_prob,CSDP.csdp),CSDP_INT,(Ptr{UInt8},CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix}),
          fname,n,k,C,a,constraints)
end

function write_sol(fname::String,
                   n::CSDP_INT,
                   k::CSDP_INT,
                   X::blockmatrix,
                   y::Ptr{Cdouble},
                   Z::blockmatrix)
    ccall((:write_sol,CSDP.csdp),
          CSDP_INT,
          (Ptr{UInt8},
           CSDP_INT,
           CSDP_INT,
           blockmatrix,
           Ptr{Cdouble},
           blockmatrix),
          fname,n,k,X,y,Z)
end

function free_prob(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix)
    ccall((:free_prob,CSDP.csdp),Nothing,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},blockmatrix),n,k,C,a,constraints,X,y,Z)
end

#function new_blockmatrix(nblocks::Integer)
#    return ccall((:new_blockmatrix,CSDP.csdp),Ptr{Cvoid},(CSDP_INT,),nblocks)
#end
#function free_blockmatrix(C::Ptr{Cvoid})
#    ccall((:free_blockmatrix,CSDP.csdp),Nothing,(Ptr{Cvoid},),C)
#end

function Base.getindex(A::blockrec, i::Integer, j::Integer)
    return ccall((:getindex,CSDP.csdp), Cdouble, (blockrec,CSDP_INT,CSDP_INT), A, i, j)
end
function getblockrec(A::blockmatrix, i::Integer)
    return ccall((:getblockrec,CSDP.csdp), blockrec, (blockmatrix,CSDP_INT), A, i)
end

function allocate_loading_prob(pC::Ref{blockmatrix}, block_dims::Ptr{CSDP_INT}, num_constraints::Integer, num_entries::Ptr{CSDP_INT}, printlevel::Integer)
    return ccall((:allocate_loading_prob,CSDP.csdp),Ptr{Cvoid},(Ref{blockmatrix},Ptr{CSDP_INT},CSDP_INT,Ptr{CSDP_INT},CSDP_INT),pC,block_dims,num_constraints,num_entries,printlevel)
end
function free_loading_prob(problem::Ptr{Cvoid})
    ccall((:free_loading_prob,CSDP.csdp),Nothing,(Ptr{Cvoid},),problem)
end
function free_loaded_prob(problem::Ptr{Cvoid},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix)
    ccall((:free_loaded_prob,CSDP.csdp),Nothing,(Ptr{Cvoid},blockmatrix,Ptr{Cdouble},blockmatrix),problem,X,y,Z)
end

function setconstant(problem::Ptr{Cvoid}, mat::Integer, ent::Cdouble)
    ccall((:setconstant,CSDP.csdp),Nothing,(Ptr{Cvoid},CSDP_INT,Cdouble),problem,mat,ent)
end
function addentry(problem::Ptr{Cvoid}, mat::Integer, blk::Integer, indexi::Integer, indexj::Integer, ent::Cdouble, allow_duplicates::Integer)
    ccall((:addentry,CSDP.csdp),CSDP_INT,(Ptr{Cvoid},CSDP_INT,CSDP_INT,CSDP_INT,CSDP_INT,Cdouble,CSDP_INT),problem,mat,blk,indexi,indexj,ent,allow_duplicates)
end

function sdp(n::CSDP_INT, k::CSDP_INT,
             C::blockmatrix, a::Ptr{Cdouble}, constant_offset::Cdouble, constraints::Ptr{constraintmatrix},
             byblocks::Ptr{Ptr{sparseblock}},
             X::blockmatrix, y::Ptr{Cdouble}, Z::blockmatrix,
             workvec1::Ptr{Cdouble},workvec2::Ptr{Cdouble},workvec3::Ptr{Cdouble},workvec4::Ptr{Cdouble},workvec5::Ptr{Cdouble},workvec6::Ptr{Cdouble},workvec7::Ptr{Cdouble},workvec8::Ptr{Cdouble},
             diagO::Ptr{Cdouble},
             besty::Ptr{Cdouble},
             O::Ptr{Cdouble}, rhs::Ptr{Cdouble}, dy::Ptr{Cdouble}, dy1::Ptr{Cdouble}, Fp::Ptr{Cdouble},
             printlevel::CSDP_INT, parameters::paramstruc)

    pobj = Ref{Cdouble}(0.0)
    dobj = Ref{Cdouble}(0.0)

    work1    = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(work1));
    work2    = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(work2));
    work3    = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(work3));
    bestx    = blockmatrix(); alloc_mat_packed(C, Ref{blockmatrix}(bestx));
    bestz    = blockmatrix(); alloc_mat_packed(C, Ref{blockmatrix}(bestz));
    cholxinv = blockmatrix(); alloc_mat_packed(C, Ref{blockmatrix}(cholxinv));
    cholzinv = blockmatrix(); alloc_mat_packed(C, Ref{blockmatrix}(cholzinv));
    Zi       = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(Zi));
    dZ       = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(dZ));
    dX       = blockmatrix(); alloc_mat(C,        Ref{blockmatrix}(dX));

    fill = Ref{constraintmatrix}(constraintmatrix(C_NULL))
    makefill(k, C, constraints, fill, work1, printlevel)

    structnnz(n, k, C, constraints)

    sort_entries(k, C, constraints)

    status = ccall((:sdp, CSDP.csdp), CSDP_INT,
                   (CSDP_INT, CSDP_INT, # n, k
                    blockmatrix, Ptr{Cdouble}, Cdouble, Ptr{constraintmatrix},
                    Ptr{Ptr{sparseblock}}, constraintmatrix, # byblocks, fill
                    blockmatrix, Ptr{Cdouble}, blockmatrix,
                    blockmatrix, blockmatrix,            # cholxinv, cholzinv
                    Ref{Cdouble}, Ref{Cdouble},          # pobj, dobj
                    blockmatrix,blockmatrix,blockmatrix, # work1, work2, work3
                    Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},
                    Ptr{Cdouble},            # diagO
                    blockmatrix,             # bestx
                    Ptr{Cdouble},            # besty
                    blockmatrix,blockmatrix, # bestz, Zi
                    Ptr{Cdouble},Ptr{Cdouble},blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},
                    CSDP_INT,paramstruc),
                   n, k,
                   C, a, constant_offset, constraints,
                   byblocks, fill[],
                   X, y, Z,
                   cholxinv,cholzinv,
                   pobj, dobj,
                   work1, work2, work3,
                   workvec1, workvec2, workvec3, workvec4, workvec5, workvec6, workvec7, workvec8,
                   diagO,
                   bestx,
                   besty,
                   bestz, Zi,
                   O, rhs, dZ, dX, dy, dy1, Fp,
                   printlevel, parameters)

    free_mat(work1);
    free_mat(work2);
    free_mat(work3);
    free_mat_packed(bestx);
    free_mat_packed(bestz);
    free_mat_packed(cholxinv);
    free_mat_packed(cholzinv);

    free_mat(Zi);
    free_mat(dZ);
    free_mat(dX);

    return status, pobj[], dobj[]
end

function loaded_sdp(problem::Ptr{Cvoid},constant_offset::Cdouble,pX::Ref{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ref{blockmatrix},printlevel::CSDP_INT,parameters::paramstruc)
    pobj = Ref{Cdouble}(0.0)
    dobj = Ref{Cdouble}(0.0)
    status = ccall((:loaded_sdp,CSDP.csdp),CSDP_INT,
        (Ptr{Cvoid},Cdouble,        Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix},Ref{Cdouble},Ref{Cdouble},CSDP_INT,paramstruc),
         problem,   constant_offset,pX,              py,               pZ,pobj,dobj,printlevel,parameters)
    return status, pobj[], dobj[]
end
function parametrized_sdp(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},constant_offset::Cdouble,pX::Ptr{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ptr{blockmatrix},printlevel::CSDP_INT,parameters::paramstruc)
    pobj = Ref{Cdouble}(0.0)
    dobj = Ref{Cdouble}(0.0)
    status = ccall((:parametrized_sdp,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,Ptr{blockmatrix},Ref{Ptr{Cdouble}},Ptr{blockmatrix},Ref{Cdouble},Ref{Cdouble},CSDP_INT,paramstruc),n,k,C,a,constraints,constant_offset,pX,py,pZ,pobj,dobj,printlevel,parameters)
    return status, pobj[], dobj[]
end
function easy_sdp(n::CSDP_INT,k::CSDP_INT,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},constant_offset::Cdouble,pX::Ptr{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ptr{blockmatrix})
    pobj = Ref{Cdouble}(0.0)
    dobj = Ref{Cdouble}(0.0)
    status = ccall((:easy_sdp,CSDP.csdp),CSDP_INT,(CSDP_INT,CSDP_INT,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,Ptr{blockmatrix},Ref{Ptr{Cdouble}},Ptr{blockmatrix},Ref{Cdouble},Ref{Cdouble}),n,k,C,a,constraints,constant_offset,pX,py,pZ,pobj,dobj)
    return status, pobj[], dobj[]
end

function tweakgap(n::CSDP_INT,k::CSDP_INT,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},gap::Cdouble,Z::blockmatrix,dZ::blockmatrix,y::Ptr{Cdouble},dy::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,work4::blockmatrix,workvec1::Ptr{Cdouble},workvec2::Ptr{Cdouble},workvec3::Ptr{Cdouble},workvec4::Ptr{Cdouble},printlevel::CSDP_INT)
    ccall((:tweakgap,CSDP.csdp),Nothing,(CSDP_INT,CSDP_INT,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},blockmatrix,blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},CSDP_INT),n,k,a,constraints,gap,Z,dZ,y,dy,work1,work2,work3,work4,workvec1,workvec2,workvec3,workvec4,printlevel)
end

function bisect_(n::Ptr{CSDP_INT},eps1::Ptr{Cdouble},d::Ptr{Cdouble},e::Ptr{Cdouble},e2::Ptr{Cdouble},lb::Ptr{Cdouble},ub::Ptr{Cdouble},mm::Ptr{CSDP_INT},m::Ptr{CSDP_INT},w::Ptr{Cdouble},ind::Ptr{CSDP_INT},ierr::Ptr{CSDP_INT},rv4::Ptr{Cdouble},rv5::Ptr{Cdouble})
    ccall((:bisect_,CSDP.csdp),CSDP_INT,(Ptr{CSDP_INT},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{CSDP_INT},Ptr{CSDP_INT},Ptr{Cdouble},Ptr{CSDP_INT},Ptr{CSDP_INT},Ptr{Cdouble},Ptr{Cdouble}),n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5)
end


text:
struct GJF{TGV, TF<:AbstractFloat, TM} <: SecondOrderNonMetropolisSampler
    ∇V!::TGV
    β::TF
    γ::TF
    M::TM
    Δt::TF
    sqrtM::TM
    a::TF
    b::TF
    σ::TF
end

"""
    GJF(∇V!, β, γ, M, Δt)

Set up the G-JF integrator for inertial Langevin.

### Fields

* ∇V!   - In place gradient of the potential
* β     - Inverse temperature
* γ     - Damping Coefficient
* M     - Mass (either scalar or vector)
* Δt    - Time step
"""
function GJF(∇V!::TGV, β::TF, γ::TF, M::TM, Δt::TF) where {TGV, TF<:AbstractFloat,TM}
    a = (1 - 0.5 * γ * Δt)/(1 + 0.5 * γ * Δt);
    b = 1 / (1 + 0.5 * γ * Δt);
    σ = sqrt(2 * γ * Δt / β);
    sqrtM = sqrt.(M);
    return GJF(∇V!, β, γ, M, Δt, sqrtM, a, b, σ)
end

mutable struct GJFState{Tq, Tx<:AbstractVector{Tq}} <:SecondOrderNonMetropolisSamplerState
    x::Tx
    ∇V::Tq
    ∇V_new::Tq
    ξ::Tq
end

function InitState!(x₀, sampler::GJF)
    ∇V = similar(x₀[1]);
    sampler.∇V!(∇V, x₀[1]);
    return GJFState(x₀, copy(∇V), similar(x₀[1]), similar(x₀[1]));
end

function InitState(x₀, sampler::GJF)

    ∇V = similar(x₀[1]);
    sampler.∇V!(∇V, x₀[1]);
    return GJFState(deepcopy(x₀), copy(∇V), similar(x₀[1]), similar(x₀[1]));
end

function UpdateState!(state::GJFState, sampler::GJF)

    @. state.ξ = randn();

    @. state.x[1] = state.x[1] + sampler.b * sampler.Δt / sampler.M * state.x[2] - 0.5 * sampler.b * sampler.Δt^2 / sampler.M * state.∇V + 0.5 * sampler.b * sampler.Δt / sampler.sqrtM * sampler.σ * state.ξ;

    sampler.∇V!(state.∇V_new, state.x[1]);

    @. state.x[2] = sampler.a * state.x[2] - 0.5 * sampler.Δt * (sampler.a * state.∇V + state.∇V_new) + sampler.b * sampler.sqrtM * sampler.σ * state.ξ;

    @. state.∇V = state.∇V_new;

    state
end

max_stars_repo_name: cortner/JuBasicMD


text:
# #######################################
# Example 10.21 Probit Model of
# Robert, C. P., & Casella, G. (2010). Monte Carlo statistical methods (2. ed., softcover reprint of the hardcover 2. ed. 2004). New York, NY: Springer.
# #######################################

using Distributions
# generate data
function genD(n = 2000, sigma = sqrt(2.0), beta = 5.0)
    # sample r from Bernoulli(1/2)
    r = rand(Bernoulli(), n)
    d = ones(Int, n)
    for i = 1:n
        # sample z from N(r*beta, sigma2)
        zi = rand(Normal(r[i]*beta, sigma))
        if zi >= 0
            d[i] = 1
        else
            d[i] = 0
        end
    end
    return r, d
end

# Gibbs sampler
function gibbs(R, D, M = 100)
    n = length(D)
    # initial
    beta = 25
    sigma = 5
    BETA = ones(M)
    SIGMA = ones(M)
    for m = 1:M
        # sample zi
        z = ones(n)
        for i = 1:n
            z[i] = rand(TruncatedNormal(0, sigma / R[i], -Inf, beta)) * (2D[i] - 1)
        end
        # sample sigma
        rz = sum(R.*z)
        println(rz)
        invs2 = rand(Gamma(n/2+1.5, 1/( rz + 1/1.5 )) )
        sigma = 1 / sqrt(invs2)
        beta = rand(Normal(0, 1/( 2*(rz/sigma^2 + 1/200) )))
        SIGMA[m] = sigma
        BETA[m] = beta
    end
    return BETA, SIGMA
end

# run
R, D = genD()
gibbs(R, D)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# index functions
Base.getindex(cube::AbstractSpectralCube, i::Union{Colon, Integer}, j::Union{Colon, Integer}, k::Union{Colon, Integer}) = cube.data[i, j, k]
function Base.getindex(cube::AbstractSpectralCube, i::AbstractFloat, j::AbstractFloat, k::Colon)
    ii = closest_index(cube.xaxis, i)
    jj = closest_index(cube.yaxis, j)
    return cube.data[ii, jj, :]
end


text:
function CategoricalMarkovChain(time_period::Int)
    # Initial.
    init_x = rand(:x => 1, Categorical([0.5, 0.5]))
    observations = Vector{Int}(undef, time_period)
    observations[1] = init_x
    x = init_x
    for i in 2:time_period
        if x == 1
            trans = rand(:x => i, Categorical([0.5, 0.5]))
        else
            trans = rand(:x => i, Categorical([0.8, 0.2]))
        end
        observations[i] = trans
    end
    return observations
end

# ------------ CHMM ------------ #

# Observation model.
obs = z -> begin
    if z == 1
        rand(:x, Categorical([0.5, 0.5]))
    else
        rand(:x, Categorical([0.2, 0.8]))
    end
end

# Transition model.
trans = prev_z -> begin
    if prev_z == 1
        rand(:z, Categorical([0.5, 0.5]))
    else
        rand(:z, Categorical([0.2, 0.8]))
    end
end

function CategoricalHiddenMarkovModel(time_period::Int)
    # Initial.
    init_z = rand(:z => 1, Categorical([0.5, 0.5]))
    latents = Vector{Int}(undef, time_period)
    latents[1] = init_z
    observations = Vector{Int}(undef, time_period)

    observations[1] = rand(:obs => 1, obs, init_z)
    z = init_z

    for i in 2:time_period
        z = rand(:trans => i, trans, z)
        observation = rand(:obs => i, obs, z)
        latents[i] = z
        observations[i] = observation
    end
    return observations
end

@testset "Particle filtering" begin

    @testset "Categorical hidden Markov model" begin
        tol = 0.1
        checks = [-1.05, -2.18, -2.57, -2.907, -4.19]
        xs = [1, 1, 2, 2, 1]
        lmles = []

        # Testing.
        init_obs = Jaynes.target([(:obs => 1, :x) => xs[1]])
        ps = Jaynes.initialize_filter(init_obs, 50000, CategoricalHiddenMarkovModel, (1, ))
        push!(lmles, get_lmle(ps))
        for t=2:5
            obs = Jaynes.target([(:obs => t, :x) =>  xs[t]])
            Jaynes.filter_step!(obs, ps, (Δ(t, IntDiff(1)), ))
            push!(lmles, get_lmle(ps))
        end
        map(enumerate(checks)) do (k, v)
            @test v ≈ lmles[k] atol = tol
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ 263aba70-d11a-11eb-2ecb-05c55f92a63c
using GMT

# ╔═╡ 2b498c40-d119-11eb-10b0-bd14161ec72e
md"
# Geodesic buffers

Buffers polygons delimite areas that are within some distance of input features. We can have them arround points or lines. We will show here examples of both.
"

# ╔═╡ 14b3b990-d125-11eb-0f85-4166a82a1131
md"
We start by computing a great circle, also known as an orthodrome between two points and interpolated at 1000 km increments.
"

# ╔═╡ af387600-d11a-11eb-277b-1126ed8360ac
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k);

# ╔═╡ e37bd330-d11a-11eb-058d-019c85fceedc
coast(region=:global, proj=(name=:ortho, center=(0,45)), land=:peru, frame=:g)
plot!(ortho, lw=0.5, marker=:circ, ms=0.1, fill=:black, show=true)

# ╔═╡ df865940-d11e-11eb-2911-f17e7dd9501f
md"
Next we will draw geodesic circles with 500 km radius with center on the orthodrome vertices. And to do it we need ofc to compute those circles. We compute the circles with the `cirgeo` function.
"

# ╔═╡ 45d991d0-d11f-11eb-37d1-8f0b5c3a1354
c = circgeo(ortho, radius=500, unit=:k);

# ╔═╡ 5fe9c7c0-d11f-11eb-080b-9713769ab827
coast(region=:global, proj=(name=:ortho, center=(0,45)), land=:peru, frame=:g)
plot!(c, lw=0.1, fill=:gray)
plot!(ortho, lw=0.5, marker=:circ, ms=0.1, fill=:black, show=true)

# ╔═╡ d6b946f0-d11f-11eb-03e4-c135f729e105
md"
Now imagine that we plot many close circles and compute the union of them all. That's how we get a geodesic buffer.
"

# ╔═╡ 390a9460-d122-11eb-3c0c-138acdd74cfe
#
line = [-37. 1; -28 26; -45 35; -19 42; -9 55; 4 64; 32 72; 85 73; 135 73; 172 73; -144 73; -78 77; -27 72; -8 65; 8 54; 18 	39; 28 24; 34 1];
# Compute the buffer line
D = buffergeo(line, width=500000);
# and plot it
coast(region=:global, proj=(name=:ortho, center=(0,45)), land=:peru, frame=:g, plot=(data=D, fill=:green), show=true)

# ╔═╡ Cell order:
# ╠═2b498c40-d119-11eb-10b0-bd14161ec72e
# ╟─14b3b990-d125-11eb-0f85-4166a82a1131
# ╠═263aba70-d11a-11eb-2ecb-05c55f92a63c
# ╠═af387600-d11a-11eb-277b-1126ed8360ac
# ╠═e37bd330-d11a-11eb-058d-019c85fceedc
# ╠═df865940-d11e-11eb-2911-f17e7dd9501f
# ╠═45d991d0-d11f-11eb-37d1-8f0b5c3a1354
# ╠═5fe9c7c0-d11f-11eb-080b-9713769ab827
# ╠═d6b946f0-d11f-11eb-03e4-c135f729e105
# ╠═390a9460-d122-11eb-3c0c-138acdd74cfe


text:
#Image
function loadImage(str::String)
	a = pointer(str)
	return AFArray{Float32}(icxx"loadImage($a);")
end

function saveImage(str::String, a::AFArray{Float32})
	b = pointer(str)
	icxx"saveImage($b, $a);"
end

function rotate(a::AFArray{Float32}, theta)
	icxx"rotate($a, $theta);"
end

function scale(a::AFArray{Float32}, scale1, scale2)
	icxx"scale($a, $scale1, $scale2);"
end

function transform(a::AFArray{Float32}, b::AFArray)
	icxx"transform($a, $b);"
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_opt_rho2_7!(A_copy)
end

@inline function exp_mono_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=-0.8990318926117189
    coeff2=1.0529942676288555
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.8990318926117189
    coeff2=1.0529942676288555
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-3.857235106310176
    coeff2=2.559191063273799
    coeff3=-0.2675261675885479
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=2.6127319091880588
    coeff2=-3.5281474819181633
    coeff3=-4.966605691968948
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-34.51931508412978
    coeff2=-32.102642943446725
    coeff3=3.7030375734689573
    coeff4=4.5107167490732225
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=22.563420733949638
    coeff2=41.22336752504774
    coeff3=8.606161607327088
    coeff4=18.485870531498684
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-1.6333763880809702
    coeff2=-12.425428631305296
    coeff3=-7.020656636475915
    coeff4=15.340351060111168
    coeff5=4.961642608409374
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=2.576203242554731
    coeff2=10.289018418102025
    coeff3=0.34028419625059303
    coeff4=26.83688350147083
    coeff5=-12.951270741270744
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.027531215857697383
    coeff2=1.4691934416553218
    coeff3=-1.327421400455049
    coeff4=-1.5545952756429124
    coeff5=-0.47775953255938497
    coeff6=-4.727574599810485
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-1.0983539425730369
    coeff2=4.121404800758164
    coeff3=-0.10671574515634609
    coeff4=0.5170030863801119
    coeff5=-2.484455273478521
    coeff6=-5.578338573600123
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.20156975628485335
    coeff2=-0.11857297824908863
    coeff3=-0.13933557327081972
    coeff4=-0.008988745081115584
    coeff5=-0.6798469796207142
    coeff6=1.5814827606241575
    coeff7=1.5449585814160118
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=1.0897581376826212e-7
    coeff2=-1.555574375395291
    coeff3=2.8159169191071185
    coeff4=0.13432676772501073
    coeff5=2.5651632231929153
    coeff6=-0.3062199878200116
    coeff7=-0.022514410581526905
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=9.22780659770276
    coeff2=-7.818752581614379
    coeff3=-4.7517129085184555
    coeff4=29.297986883148763
    coeff5=-0.0069810802300877255
    coeff6=0.0011903433805721666
    coeff7=-1.1317145971276296e-9
    coeff8=-1.3350038182653414e-45
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
__precompile__()

module moeum
    using Logging; Logging.configure(level = Logging.OFF)

    export MOEUM, io

    #includes
    include("io/input.jl")
    include("structs/Core.jl")

    #exporting to root
    MOEUM = Core.MOEUM
    from_dict = input.from_dict
    from_dataframe = input.from_dataframe
    from_csv = input.from_csv

end


text:
using KitBase

np = 10000 # number of particles
ptc = Array{Particle1D}(undef, np)

for i = 1:np
    rd1 = rand(3)
    rd2 = rand(3)
    rd = rand()

    m = 1e-5
    v = @. sqrt(-log(rd1) / 1.0) * sin(2.0 * π * rd2)
    x = 0.0 + (rd - 0.5) * 1.0
    idx = i

    if v[1] < 0
        tb = (0.0 - 0.5 * 1.0 - x) / v[1]
    elseif v[1] > 0
        tb = (0.0 + 0.5 * 1.0 - x) / v[1]
    else
        tb = 1.0
    end

    e = 0.0
    flag = 1

    ptc[i] = Particle1D(m, x, v, e, idx, flag, tb)
end

v = zeros(np, 3)
for i = 1:np
    v[i, :] .= ptc[i].v
end

using Plots
histogram(v[:, 1], label = "u")
histogram!(v[:, 2], label = "v")
histogram!(v[:, 3], label = "w")


text:
using CUDA
using FRGRealTime
using Plots
vx1=collect(1:0.1:800)
vx2=collect(0.0:0.001:1.0)
cu(vx2)
using JuliaZH
sign(0)

FRGRealTime.dkVImintqs_gpu.(
    10.0,
    0.0001,
    Epi.(vx1, m2fun.(vx1))|>cu,
    vx1|>cu,
    vx2 .* (vx1 .- 800.0) .+ 800.0|>cu,
    m2fun.(vx2 .* (vx1 .- 800.0) .+ 800.0)|>cu,
    145.0,
    4.0,
    lamfun.(vx2 .* (vx1 .- 800.0) .+ 800.0)|>cu,
)

using FastGaussQuadrature

gausslegendre(1000)
cuv1,cuweight1=cu(gausslegendre(10000))
cuv2,cuweight2=cu.(collect.(transpose.(gausslegendre(10000))))

function domaintrans!(v1,v2,a,b)
    v1.=((b-a)/2).*v1.+((b-a)/2);
    v2.=((b-a)/2).*v2;
end

domaintrans!(cuv1,cuweight1,1.0,800.0)
domaintrans!(cuv2,cuweight2,0.0,1.0)

v2trans=cuv2 .* (cuv1 .- 800.0) .+ 800.0
lamvec=lamfun.(Array(v2trans))|>cu
mvec2=m2fun.(Array(v2trans))|>cu
mvec1=m2fun.(Array(cuv1))|>cu

function myfun5(w1,w2,v1,v2,mv1,mv2,lamv,p0,ps,Npi,T,UVScale)
    w1*w2*(v1 - UVScale) *2 *FRGRealTime.dkVImintqs_gpu(
        p0,
        ps,
        Epi(v1, mv1),
        v1,
        v2,
        mv2,
        T,
        Npi,
        lamv,
    )
end

function sumfun2(p0,ps)::Float32
    tempv=myfun5.(
        cuweight1,
        cuweight2,
        cuv1,
        v2trans,
        mvec1,
        mvec2,
        lamvec,
        p0,
        ps,
        4.0f0,
        145.0f0,
        800.0f0)
    return sum(tempv)
end

cup0vec=collect(1.0f0:0.1f0:4.9f0)
sumfun2.(cup0vec,0.001f0)

CUDA.memory_status()
CUDA.reclaim()

cpuweight1=cuweight1|>Array
cpuweight2=cuweight2|>Array
cpuv1=cuv1|>Array
cpuv2trans=v2trans|>Array
cpumvec1=mvec1|>Array
cpumvec2=mvec2|>Array
cpulamvec=lamvec|>Array

@btime(myfun5.(
    cpuweight1,
    cpuweight2,
    cpuv1,
    cpuv2trans,
    cpumvec1,
    cpumvec2,
    cpulamvec,
    10.0f0,
    0.001f0,
    4.0f0,
    145.0f0,
    800.0f0,
))

@time(1+1)

vx2 .* (vx1 .- 800.0) .+ 800.0|>cu

FRGRealTime.dkF1Allintqs_gpu.(
    100.0,
    10.001,
    10*vx2 |> cu,
    vx1 |> cu,
    m2fun.(vx1) |> cu,
    145.0,
)|>sum

function testfun2(x,y)::Float32
    sin(x)-y
end
cuvx1=vx1|>cu
cuvx2=vx2|>cu

function int4()
    mapreduce(x->testfun.(x,cuvx2),+,cuvx1)
end

mapreduce(x->testfun2.(x,cuvx2),+,cuvx1)

mapreduce(x->int3(x),+,cuvx2)

int2(1.0f0)

ps_temp=7.0
plot!(
    p0 -> FRGRealTime.loopfunpm(
        sqrt(p0^2+ps_temp^2-5^2),
        ps_temp,
        10.0,
        -80.0,
        10.0,
    ),
    0.0,
    100.0,
)


text:
module lipm

using StaticArrays
using MeshCat
using GeometryBasics
using Colors
using CoordinateTransformations, Rotations
using LinearAlgebra

include("dynamics/walking_dynamics.jl")
include("dynamics/interface.jl")
include("visualization/animate.jl")

export
    move_to_position,
    animate_walking_trajectory

end # module


text:
export DeltaPulse, SqExpPulse, BSplinePulse

"""
Pulse
Pixel response function for a radio image model. This makes
a discrete sampling continuous by picking a certain *smoothing*
kernel for the image.

# Notes
To see the implemented Pulses please use the subtypes function i.e.
`subtypes(Pulse)`
"""
abstract type Pulse <: AbstractModel end

visanalytic(::Type{<:Pulse}) = IsAnalytic()
imanalytic(::Type{<:Pulse}) = IsAnalytic()
isprimitive(::Type{<:Pulse}) = IsPrimitive()

@inline intensity_point(p::Pulse, x,y) = κ(p::Pulse, x)*κ(p::Pulse, y)
@inline visibility_point(p::Pulse, u,v) = ω(p::Pulse, u)*ω(p::Pulse, v)

flux(p::Pulse) = κflux(p)^2

"""
    $(TYPEDEF)
A dirac comb pulse function. This means the image is just the
dicrete Fourier transform
"""
struct DeltaPulse{T} <: Pulse end
DeltaPulse() = DeltaPulse{Float64}()
# This should really be a delta function but whatever
κflux(::DeltaPulse{T}) where {T} = one(T)
@inline κ(::DeltaPulse{T}, x) where {T} = one(T)
@inline ω(::DeltaPulse{T}, u) where {T} = one(T)
@inline radialextent(::Pulse) = 1.0

"""
    $(TYPEDEF)
Normalized square exponential kernel, i.e. a Gaussian. Note the
smoothness is modfied with `ϵ` which is the inverse variance in units of
1/pixels².
"""
struct SqExpPulse{T} <: Pulse
    ϵ::T
end
@inline @fastmath κ(b::SqExpPulse, x) = exp(-0.5*b.ϵ^2*x^2)/sqrt(2*π/b.ϵ^2)
@inline κflux(::SqExpPulse{T}) where {T} = one(T)
@inline @fastmath ω(b::SqExpPulse, u) = exp(-2*(π*u/b.ϵ)^2)
@inline radialextent(p::SqExpPulse) = 5/p.ϵ

@doc raw"""
    $(TYPEDEF)
Uses the basis spline (BSpline) kernel of order `N`. These are the kernel that come
from recursively convolving the tophat kernel
```math
    B_0(x) = \begin{cases} 1 & |x| < 1 \\ 0 & otherwise \end{cases}
```
`N` times.

## Notes

BSpline kernels have a number of nice properties:
1. Simple frequency response $\sinc(u/2)^N$
2. preserve total intensity

For `N`>1 these kernels aren't actually interpolation kernels however, this doesn't matter
for us.

Currently only the 0,1,3 order kernels are implemented.
"""
struct BSplinePulse{N} <: Pulse end
@inline ω(::BSplinePulse{N}, u) where {N} = sinc(u)^(N+1)
@inline κflux(::BSplinePulse) = 1.0

@inline κ(::BSplinePulse{0}, x::T) where {T} = abs(x) < 0.5 ? one(T) : zero(T)

@inline function κ(::BSplinePulse{1}, x::T) where {T}
    mag = abs(x)
    return mag < 1 ? 1-mag : zero(T)
end

@inline function κ(::BSplinePulse{3}, x::T) where {T}
    mag = abs(x)
    if mag < 1
        return evalpoly(mag, (4, 0, -6, 3))/6
    elseif 1 ≤ mag < 2
        return evalpoly(mag, (8, -12, 6, -1))/6
    else
        return zero(T)
    end
end


text:
module BasicTextRender

using ImageCore, ImageTransformations, ImageFiltering, Serialization

charStack = nothing
loadedFont = nothing
kern = ImageFiltering.KernelFactors.IIRGaussian((2,2))
loadedScale = nothing
warned = false

"""
    rendertext(str::String;
        height::Int=0,
        font::String="Courier",
        color::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,1.0),
        backgroundColor::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,0.0))

Render monospaced text from character lookup tables.
Not recommended for high quality output.
Loads slowly first time, or whenever the font is changed. Much faster next time.

font: "Courier", "Courier New", "Monaco"
"""
function rendertext(str::String;
        height::Int=0,
        font::String="Courier",
        color::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,1.0),
        backgroundColor::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,0.0))

    global charStack, loadedFont, kern, loadedScale, warned

    lookupdir = joinpath(dirname(@__DIR__),"gen","fonts",font)
    !isdir(lookupdir) && error("""Font "$font" not supported. Try "Courier", "Courier New" or "Monaco" """)

    if loadedFont != font
        charStack = Serialization.deserialize(joinpath(lookupdir,"charStack"))
        loadedFont = font
    end

    chars = collect(str)
    codes = Int.(codepoint.(chars)) .+ 1
    raw_chars = view(charStack, codes)
    raw_img = hcat(raw_chars...) ./ 255

    colormap = ImageCore.colorsigned(RGBA(0.0,0.0,0.0,0.0), backgroundColor, color)

    scale = height / size(raw_img, 1)
    if scale == 1 || (height <= 0)
        return colormap.(raw_img)
    elseif scale < 1
        if loadedScale != scale
            σ = 0.25 * (1/scale)
            if σ < 1
                kern = ImageFiltering.KernelFactors.gaussian((σ,σ))
            else
                kern = ImageFiltering.KernelFactors.IIRGaussian((σ,σ))
            end
        end

        return colormap.(ImageTransformations.imresize(ImageFiltering.imfilter(raw_img, kern, NA()), ratio=scale))
    else
        !warned && (@warn "Height exceeds font lookup table height. Sharpness will be suboptimal") && (warned = true)
        return colormap.(ImageTransformations.imresize(raw_img, ratio=scale))
    end
end

"""
    overlaytext!(img, str::String, height::Int, bottomleft::Tuple{Int,Int};
            font::String="Courier",
            color::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,1.0))

Overlay text on `img` in-place.
"""
function overlaytext!(img, str::String, height::Int, bottomleft::Tuple{Int,Int};
            font::String="Courier",
            color::ColorTypes.RGBA{Float64}=ColorTypes.RGBA(0.0,0.0,0.0,1.0))

    img_text = rendertext(str, font=font, height=height, color=color)
    width = size(img_text,2)
    x, y = bottomleft
    img[x:x+height-1, y:y+width-1] .= ((1 .- alpha.(img_text)) .* img[x:x+height-1, y:y+width-1]) .+ (alpha.(img_text) .* img_text)
end

export rendertext, overlaytext!
end #module


text:
# Process-local dictionary mapping from cluster names to instances of `Cluster`
# Might contain stale information
global clusters = Dict()

function create_cluster(;
    name::Union{String,Nothing} = nothing,
    instance_type::Union{String,Nothing} = "m4.4xlarge",
    max_num_workers::Union{Int,Nothing} = 2048,
    initial_num_workers::Union{Int,Nothing} = 16,
    min_num_workers::Union{Int,Nothing} = 0,
    iam_policy_arn::Union{String,Nothing} = nothing,
    s3_bucket_arn::Union{String,Nothing} = nothing,
    s3_bucket_name::Union{String,Nothing} = nothing,
    scaledown_time = 25,
    region = nothing,
    vpc_id = nothing,
    subnet_id = nothing,
    nowait=false,
    kwargs...,
)

    # Configure using parameters
    c = configure(; kwargs...)

    clusters = get_clusters(; kwargs...)
    if isnothing(name)
        name = "Cluster " * string(length(clusters) + 1)
    end
    if isnothing(region)
        aws_region = get_aws_config_region()
    end

    # Check if the configuration for this cluster name already exists
    # If it does, then recreate cluster
    if haskey(clusters, name)
        if clusters[name].status == :terminated
            @info "Started re-creating cluster named $name"
            send_request_get_response(
                :create_cluster,
                Dict("cluster_name" => name, "recreate" => true),
            )
            if !nowait
                wait_for_cluster(name)
            end
            return get_cluster(name)
        else
            error("Cluster with name $name already exists and its current status is $(string(clusters[name].status))")
        end
    end

    # Construct arguments
    if !isnothing(s3_bucket_name)
        s3_bucket_arn = "arn:aws:s3:::$s3_bucket_name*"
    elseif !isnothing(s3_bucket_arn)
        s3_bucket_name = last(split(s3_bucket_arn, ":"))
    end
    if isnothing(s3_bucket_arn)
        s3_bucket_arn = ""
    elseif !(s3_bucket_name in s3_list_buckets(get_aws_config()))
        error("Bucket $s3_bucket_name does not exist in the connected AWS account")
    end

    # Construct cluster creation
    cluster_config = Dict(
        "cluster_name" => name,
        "instance_type" => instance_type,
        "max_num_workers" => max_num_workers,
        "initial_num_workers" => initial_num_workers,
        "min_num_workers" => min_num_workers,
        "aws_region" => region,
        "s3_read_write_resource" => s3_bucket_arn,
        "scaledown_time" => scaledown_time,
        "recreate" => false,
    )
    if haskey(c["aws"], "ec2_key_pair_name")
        cluster_config["ec2_key_pair"] = c["aws"]["ec2_key_pair_name"]
    end
    if !isnothing(iam_policy_arn)
        cluster_config["additional_policy"] = iam_policy_arn
    end
    if !isnothing(vpc_id)
        cluster_config["vpc_id"] = vpc_id
    end
    if !isnothing(subnet_id)
        cluster_config["subnet_id"] = subnet_id
    end

    @info "Started creating cluster named $name"

    # Send request to create cluster
    send_request_get_response(:create_cluster, cluster_config)

    if !nowait
        wait_for_cluster(name)
    end

    # Cache info
    global clusters
    cl = Cluster(name, get_cluster_status(name), "", 0, s3_bucket_arn)
    clusters[name] = cl

    return cl
end

function destroy_cluster(name::String; kwargs...)
    configure(; kwargs...)
    @info "Destroying cluster named $name"
    send_request_get_response(:destroy_cluster, Dict{String,Any}("cluster_name" => name))
end

function delete_cluster(name::String; kwargs...)
    configure(; kwargs...)
    @info "Deleting cluster named $name"
    send_request_get_response(
        :destroy_cluster,
        Dict{String,Any}("cluster_name" => name, "permanently_delete" => true),
    )
end

function update_cluster(name::String; kwargs...)
    configure(; kwargs...)
    @info "Updating cluster named $name"
    send_request_get_response(
        :update_cluster,
        Dict{String, Any}("cluster_name" => name)
    )
end

function assert_cluster_is_ready(name::String; kwargs...)
    @info "Setting status of cluster named $name to running"

    # Configure
    configure(; kwargs...)

    send_request_get_response(:set_cluster_ready, Dict{String,Any}("cluster_name" => name))
end

struct Cluster
    name::String
    status::Symbol
    status_explanation::String
    num_jobs_running::Int32
    s3_bucket_arn::String
end

parsestatus(status) =
    if status == "creating"
        :creating
    elseif status == "destroying"
        :destroying
    elseif status == "updating"
        :updating
    elseif status == "failed"
        :failed
    elseif status == "starting"
        :starting
    elseif status == "stopped"
        :stopped
    elseif status == "running"
        :running
    elseif status == "terminated"
        :terminated
    elseif status == "unknown"
        :unknown
    else
        error("Unexpected status ", status)
    end

function get_clusters(cluster_name=nothing; kwargs...)
    @debug "Downloading description of clusters"
    configure(; kwargs...)
    filters = Dict()
    if !isnothing(cluster_name)
        filters["cluster_name"] = cluster_name
    end
    response = send_request_get_response(:describe_clusters, Dict{String,Any}("filters"=>filters))
    clusters_dict = Dict(
        name => Cluster(
            name,
            parsestatus(c["status"]),
            haskey(c, "status_explanation") ? c["status_explanation"] : "",
            c["num_jobs"],
            c["s3_read_write_resource"],
        ) for (name, c) in response["clusters"]
    )

    # Cache info
    global clusters
    for (name, c) in clusters_dict
        clusters[name] = c
    end

    return clusters_dict
end

function get_cluster_s3_bucket_arn(cluster_name=get_cluster_name(); kwargs...)
    configure(; kwargs...)
    global clusters
    # Check if cached, sine this property is immutable
    if !haskey(clusters, cluster_name)
        get_cluster(cluster_name)
    end
    return clusters[cluster_name].s3_bucket_arn
end

function get_cluster_s3_bucket_name(cluster_name=get_cluster_name(); kwargs...)
    configure(; kwargs...)
    return s3_bucket_arn_to_name(get_cluster_s3_bucket_arn(cluster_name))
end

get_cluster(name::String=get_cluster_name(), kwargs...) = get_clusters(name; kwargs...)[name]

get_running_clusters(args...; kwargs...) = filter(entry -> entry[2].status == :running, get_clusters(args...; kwargs...))

function get_cluster_status(name::String=get_cluster_name(), kwargs...)
    global clusters
    if haskey(clusters, name)
        if clusters[name].status == :failed
            @error c.status_explanation
        end
    end
    c = get_cluster(name; kwargs...)
    if c.status == :failed
        @error c.status_explanation
    end
    c.status
end

function wait_for_cluster(name::String=get_cluster_name(), kwargs...)
    t = 5
    cluster_status = get_cluster_status(name; kwargs...)
    p = nothing
    while (cluster_status == :creating || cluster_status == :updating)
        if isnothing(p)
            if cluster_status == :creating
                p = ProgressUnknown("Setting up cluster $(name)", spinner=true)
            else
                p = ProgressUnknown("Updating cluster $(name)", spinner=true)
            end
        end
        sleep(t)
        next!(p)
        if t < 80
            t *= 2
        end
        cluster_status = get_cluster_status(name; kwargs...)
    end
    if !isnothing(p)
        finish!(p, spinner = (cluster_status == :running ? '✓' : '✗'))
    end
    if cluster_status == :running
        # @info "Cluster $name is ready"
    elseif cluster_status == :terminated
        error("Cluster $name no longer exists")
    elseif cluster_status != :creating && cluster_status != :updating
        error("Failed to set up cluster named $name")
    else
        error("Cluster $name has unexpected status: $cluster_status")
    end
end

function upload_to_s3(src_path; dst_name=basename(src_path), cluster_name=get_cluster_name(), kwargs...)
    configure(; kwargs...)
    bucket_name = get_cluster_s3_bucket_name(cluster_name)
    s3_dst_path = S3Path("s3://$bucket_name/$dst_name", config=get_aws_config())
    if startswith(src_path, "http://") || startswith(src_path, "https://")
        Base.download(
            src_path,
            s3_dst_path
        )
    elseif startswith(src_path, "s3://")
        cp(
            S3Path(src_path),
            s3_dst_path
        )
    else
        if startswith(src_path, "file://")
            src_path = src_path[8:end]
        end
        if isfile(src_path)
            cp(
                Path(src_path),
                s3_dst_path
            )
        else # isdir
            for f_name in readdir(src_path)
                cp(
                    Path("$src_path/$f_name"),
                    S3Path(
                        "s3://$bucket_name/$(basename(src_path))/$(f_name)",
                        config=get_aws_config()
                    )
                )
            end
        end
    end
    return dst_name
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
function SolarRadiationPressure(arg0::ExtendedPVCoordinatesProvider, arg1::jdouble, arg2::RadiationSensitive)
    return SolarRadiationPressure((ExtendedPVCoordinatesProvider, jdouble, RadiationSensitive), arg0, arg1, arg2)
end

function SolarRadiationPressure(arg0::jdouble, arg1::jdouble, arg2::ExtendedPVCoordinatesProvider, arg3::jdouble, arg4::RadiationSensitive)
    return SolarRadiationPressure((jdouble, jdouble, ExtendedPVCoordinatesProvider, jdouble, RadiationSensitive), arg0, arg1, arg2, arg3, arg4)
end

function acceleration(obj::SolarRadiationPressure, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::SolarRadiationPressure, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_lighting_ratio(obj::SolarRadiationPressure, arg0::FieldVector3D, arg1::Frame, arg2::FieldAbsoluteDate)
    return jcall(obj, "getLightingRatio", RealFieldElement, (FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2)
end

function get_lighting_ratio(obj::SolarRadiationPressure, arg0::Vector3D, arg1::Frame, arg2::AbsoluteDate)
    return jcall(obj, "getLightingRatio", jdouble, (Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::SolarRadiationPressure)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end


text:
## I try to implement again a simple VAE using Julia
## There is a tutorial with a working VAE in the Flux documentation
## (https://github.com/jessebett/model-zoo/blob/master/vision/mnist/vae.jl)
## I try to do mine in a different way

## I am using the original Binarized MNIST
## you can download it from http://www.dmi.usherb.ca/~larocheh/mlpython/_modules/datasets/binarized_mnist.html
## I assume the files are into a folder called "Original_MNIST_binarized"
## I am not passing a directory to the loading function (possibly a TODO)

using Flux, Statistics
using Flux: Tracker, throttle, params, binarycrossentropy, gradient, @epochs, testmode!, glorot_uniform
#flux function for update parameters
using Flux.Tracker: update!
using Base.Iterators: partition
using Juno: @progress
using MLDataUtils
using Images
include("utils.jl")

img(x) = Gray.(permutedims(reshape(x, 28, 28)))

function log_standard_gaussian(x)
    #return torch.sum(-0.5 * math.log(2 * math.pi) - x ** 2 / 2, dim=-1)
    #log_pdf =  -0.5f0 * log(2 * pi) .- x.^2 / 2
    log_pdf =  @. -0.5 * log(2. * pi) - x^2 / 2.
    #@show(A)
    # alternatively, we can do
    #Norm = Normal(0)
    #A = log(pdf(Norm, x))
    #@show(A)
    return log_pdf
end

# function log_gaussian(x, mu, log_var)
#     std = exp.(log_var * 0.5f0)
#     #@show log_var
#     #@show std
#     Norm = Normal(mu, std)
#     A = log(pdf(Norm, x))
#     return A
# end

function log_gaussian(x, mu, log_var)
    log_pdf = - 0.5 * log(2. * pi) .- log_var / 2. - (x .- mu)^2 / (2. * exp.(log_var))
    #log_pdf = @. - 0.5 * log(2 * pi) - (log_var + eps(Float32))1 / 2 - (x - mu)^2 / (2 * exp(log_var))
    # print('Size log_pdf:', log_pdf.shape)
    return log_pdf
end

# function binary_cross_entropy(recon, x)
#     # we can also use "using Flux: binarycrossentropy"
#     # sum(binarycrossentropy.(recon, x))
#     s = @. x * log(recon + Float32(1e-8)) + (1 - x) * log(1 - recon + Float32(1e-8))
#     return -sum(s)
# end  THIS FUNCTION IS SUPER SLOW, MY GOODNESS! USE THE ONE PROVIDE BY FLUX

function reparametrization_trick(mu, log_var)
    z = mu + rand(Float64) * exp(log_var * 0.5)
    return z
end

function analytic_kl(mu, log_var)
    kl = (log_var - mu.^2 .+ 1.0 - exp.(log_var)) * 0.5
    return kl
end

## todo: I have to understand why the approaximate KL does not work
function _kl_divergence(z, mu_q, log_var_q, mu_p = nothing, log_var_p = nothing)
    # instead of computing it analytically, we can estimate it
    # using monte carlo methods
    qz = log_gaussian.(z, mu_q, log_var_q)
    ## we should do the same with p
    if isnothing(mu_p) & isnothing(log_var_p)
        pz = log_standard_gaussian.(z)
    else
        pz = log_gaussian.(z, mu_p, log_var_p)
    end
    kl = qz - pz
    return kl
end

const batch_size = 64
const latent_dim = 3
const hidden_dim = 400
const n_epochs = 200

# we can load the MNIST dataset
x_train, x_valid, x_test = load_binarized_mnist(batch_size, true)
example = x_train[:,1]
save("original2.png", img(example))

N = size(x_train, 2)
N_test = size(x_test, 2)
x_train_batch = [x_train[:,i] for i in Iterators.partition(1:N,batch_size)]
x_test_batch = [x_test[:,i] for i in Iterators.partition(1:N_test,batch_size)]
# we have to create the encoder
h, mu, log_var = Dense(28*28, hidden_dim, relu, initW = glorot_uniform, initb = zeros), Dense(hidden_dim, latent_dim, initW = glorot_uniform, initb = zeros), Dense(hidden_dim, latent_dim, initW = glorot_uniform, initb = zeros)
encoder(X) = (hidden_activation = h(X); (mu(hidden_activation), log_var(hidden_activation)))

# the decoder
decoder = Chain(Dense(latent_dim, hidden_dim, relu, initW = glorot_uniform, initb = zeros),
                Dense(hidden_dim, 28*28, sigmoid, initW = glorot_uniform, initb = zeros))

# we define the
# we have to define a method to samples
get_sample() = decoder(randn(Float64,latent_dim))
#function get_sample()
#    testmode!(decoder)
#    return decoder(randn(Float32,latent_dim))
#end
# now we have to define the callback --> used to observe the training process
function compute_loss(_data::Matrix)
    # we should pass our data into the encoder
    (_mu, _log_var) = encoder(_data)
    # then we use the reparametrization trick to sample z
    _z = reparametrization_trick.(_mu, _log_var)
    # we pass the z through the decoder to get the mean for each pixel
    mu_pixels = decoder(_z)

    # now we should compute the compute_loss
    likelihood = - sum(binarycrossentropy.(mu_pixels, _data))
    #@show likelihood
    kl = - sum(analytic_kl(_mu, _log_var))
    #kl = _kl_divergence(_z, _mu, _log_var)
    #©@show sum(kl)
    elbo = likelihood - kl
    loss = -elbo
    #@show loss

    #@show loss

    return loss, -likelihood, kl
end

opt = ADAM(3e-4)
ps = params(h, mu, log_var, decoder)

#
function train!()
    for ep in 1:n_epochs
        @info "Training Epoch $ep..."
        tot_recon = 0
        tot_kl = 0
        tot_loss = 0
        j = 0
        for data in x_train_batch
            #testmode!((h, mu, log_var, decoder), false)
            #@show j
            j += 1
            #we have to compute the loss of this
            (loss, _recon_err, _kl) = compute_loss(data)
            grad = Flux.Tracker.gradient(()-> loss/size(data,2), ps)
            update!(opt, ps, grad)
            tot_recon += _recon_err
            tot_kl += sum(_kl)
            tot_loss += loss
            #Flux.back!(l)<
            #Flux.Optimise._update_params!(opt, ps)
        end
        # we are at the end of the epochs
        println("Epoch: $ep, -ELBO: $(tot_loss/size(x_train,2)), recon err: $(tot_recon/size(x_train,2)), KL: $(tot_kl/size(x_train,2))")

        # ## TODO: understand how to sample during training, becuase at the moment it gets problems
        # (h, mu, log_var, decoder) = mapleaves(Flux.Tracker.data, ps)
        testmode!((h, mu, log_var, decoder), true)

        sample = img.([get_sample() for i = 1:64])
        @show(size(sample))
        grid_samples = vcat([hcat(imgs...) for imgs in partition(sample, 8)]...)
        save("samples/samples_from_random_epoch_$ep.png", grid_samples)

        example = [x_test[:,rand(1:N_test)] for i = 1:64]
        example2 = hcat(example...)
        grid_original = vcat([hcat(imgs...) for imgs in partition(img.(example), 8)]...)
        save("reconstructions/original_epoch_$ep.png", grid_original)
        m, lv = encoder(example2)
        zed2 = reparametrization_trick.(m,lv)
        recon = decoder(zed2)

        recon_to_save = img.([recon[:,i] for i = 1:64])
        grid_recon = vcat([hcat(imgs...) for imgs in partition(recon_to_save, 8)]...)
        save("reconstructions/reconstruction_epoch_$ep.png", grid_recon)
        testmode!((h, mu, log_var, decoder), false)
    end
end

train!()
#
#@epochs n_epochs train!()
# evalcb = throttle(() -> @show(ELBO(x_train[:, rand(1:N, batch_size)])/batch_size) , 30)
# @progress for i = 1:10
#   @info "Epoch $i"
#   Flux.train!(loss, ps, zip(x_train_batch), opt, cb=evalcb)
# end
(h, mu, log_var, decoder) = mapleaves(Flux.Tracker.data, (h, mu, log_var, decoder))
testmode!((h, mu, log_var, decoder), true)

sample = img.([get_sample() for i = 1:64])
@show(size(sample))
grid_samples = vcat([hcat(imgs...) for imgs in partition(sample, 8)]...)
save("samples/samples_from_random_epoch_$ep.png", grid_samples)

example = [x_test[:,rand(1:N_test)] for i = 1:64]
example2 = hcat(example...)
grid_original = vcat([hcat(imgs...) for imgs in partition(img.(example), 8)]...)
save("reconstructions/original_epoch_$ep.png", grid_original)
m, lv = encoder(example2)
zed2 = reparametrization_trick.(m,lv)
recon = decoder(zed2)

recon_to_save = img.([recon[:,i] for i = 1:64])
grid_recon = vcat([hcat(imgs...) for imgs in partition(recon_to_save, 8)]...)
save("reconstructions/reconstruction_epoch_$ep.png", grid_recon)


text:
import MPI

function main()
    MPI.Init()

    comm = MPI.COMM_WORLD
    println("Hello world, I am $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))")
    MPI.Barrier(comm)

    MPI.Finalize()
end

main()


text:
module Outlier
export cubeAnomalies, simpleAnomalies
using ..DAT
using ..Cubes
import ..Proc
using MultivariateAnomalies

function getDetectParameters(methodlist,trainarray,ntime)
  P = getParameters(methodlist,trainarray);
  init_detectAnomalies(zeros(Float64,ntime,size(trainarray,2)), P);
  P
end

"""
    cubeAnomalies(cube, methods, trainArray)

A simple wrapper around the function `detectAnomalies!` from the [MultivariateAnomalies](https://github.com/milanflach/MultivariateAnomalies.jl)
package.

* `cube` data cube with a axes: `TimeAxis`, `VariableAxis`
* `methods` vector of methods to be applied, choose from: `KDE`,`T2`,`REC`,`KNN-Gamma`,`KNN-Delta`,`SVDD`,`KNFST`
* `trainArray` a matrix of `nsample` x `nvar`, to estimate the training parameters for the model. Ideally does not contain any extreme values

**Input Axes** `TimeAxis`, `Variable`axis

**Output Axes** `TimeAxis`, `Method`axis
"""
function cubeAnomalies(c::AbstractCubeData,methods,trainArray;kwargs...)
  indims = InDims(TimeAxis,VariableAxis)
  outdims = OutDims(TimeAxis,CategoricalAxis("Method",methods))
  P = getDetectParameters(methods,trainArray,length(getAxis(TimeAxis,c)))
  mapCube(cubeAnomalies,c,P;indims=indims,outdims=outdims,kwargs...)
end

function cubeAnomalies(xout::AbstractArray, xin::AbstractArray, P::MultivariateAnomalies.PARAMS)
 detectAnomalies!(Float64.(xin), P)
 for i = 1:length(P.algorithms)
   copy!(view(xout, :, i), MultivariateAnomalies.return_scores(i, P))
 end
 return(xout)
end

function simpleAnomalies(c::AbstractCubeData,methods;kwargs...)
  indims = InDims(TimeAxis,VariableAxis)
  outdims = OutDims(TimeAxis,CategoricalAxis("Method",methods))
  mapCube(simpleAnomalies,c,methods;indims=indims,outdims=outdims,kwargs...)
end

function simpleAnomalies(xout::AbstractArray, xin::AbstractArray,methods)
  if !any(isnan,xin)
    P=getParameters(methods,xin)
    res=detectAnomalies(xin,P)
    for i=1:length(res)
      xout[:,i]=res[i]
    end
  else
    xout[:]=NaN
    end
end

end


text:
using OpenQuantumSystems
using DelimitedFiles

Nvib = 3

HR = 0.01
shift = (2.0 * HR)^0.5
omega = 300.0
mols = [
    Molecule([Mode(omega, shift)], Nvib, [12500.0, 12750.0]),
    Molecule([Mode(omega, shift)], Nvib, [12500.0, 12800.0]),
]

agg = Aggregate(mols)
for mol_i = 2:length(agg.molecules)
    agg.coupling[mol_i, mol_i+1] = 100
    agg.coupling[mol_i+1, mol_i] = 100
end

aggInds = getIndices(agg)
aggIndLen = length(aggInds)
base = GenericBasis([aggIndLen])
FCFact = getFranckCondonFactors(agg, aggInds)
Ham = getAggHamiltonian(agg, aggInds, FCFact; groundEnergy = true)
println(size(Ham.data))
writedlm("test/hamiltonian.csv", Ham.data, ',')

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems


text:
# Pressurized cylinder clamped around the ends.
#
# Example introduced in
# @article{Lee2004,
#    author = {Lee, P. S. and Bathe, K. J.},
#    title = {Development of MITC isotropic triangular shell finite elements},
#    journal = {Computers & Structures},
#    volume = {82},
#    number = {11-12},
#    pages = {945-962},
#    ISSN = {0045-7949},
#    DOI = {10.1016/j.compstruc.2004.02.004},
#    year = {2004},
#    type = {Journal Article}
# }

module cos_2t_press_cylinder_fixed_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

# Parameters:
E = 2.0e5
nu = 1/3;
pressure = 1.0;
R = 1.0;
L = 2.0;

# The cylinder axis is parallel to Y

function cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    n = cross(tangents[:, 1], tangents[:, 2])
    n = n/norm(n)
    # r = vec(XYZ); r[2] = 0.0
    csmatout[:, 3] .= n
    csmatout[:, 2] .= (0.0, 1.0, 0.0)
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    r = vec(XYZ); r[2] = 0.0
    r .= vec(r)/norm(vec(r))
    theta = atan(r[3], r[1])
    forceout[1:3] = r*pressure*cos(2*theta)
    forceout[4:6] .= 0.0
    # @show dot(n, forceout[1:3])
    return forceout
end

function _execute_dsg_model(formul, n = 8, thickness = R/100, visualize = true)
    tolerance = R/n/1000
    fens, fes = T3block(90/360*2*pi,L/2,n,n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a), y, R*cos(a))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

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

    # Apply EBC's
    # plane of symmetry perpendicular to Z
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0 0], inflate = tolerance)
    for i in [3,4,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # clamped edge perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the fixed edge
    # nl = selectnode(fens; box = Float64[R R L/2 L/2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))

    fi = ForceIntensity(FFlt, 6, computetrac!);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    strainenergy = 1/2 * U' * K * U
    @info "Strain Energy: $(round(strainenergy, digits = 9))"

    # Generate a graphical display of resultants
    ocsys = CSys(3, 3, cylindrical!)
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
    end
    vtkwrite("cos_2t_press_cylinder_fixed-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("n$nc", fld.values))
    end
    vtkwrite("cos_2t_press_cylinder_fixed-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
    end
    vtkwrite("cos_2t_press_cylinder_fixed-$(n)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if visualize
        scattersysvec!(dchi, (L/8)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
            #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end

    return strainenergy
end

function test_convergence(formul = FEMMShellT3FFModule, thickness = R/100)
    @info "Pressurized Cylindrical shell, fixed ends, formulation=$(formul)"
    for n in [8, 16, 32, 64, 128]
        _execute_dsg_model(formul, n, thickness, false)
    end
    return true
end

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

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
export documenter_keygen

using Base64
import LibGit2: GITHUB_REGEX

"""
    documenter_keygen()

Generates the SSH keys that are required for the automatic deployment of documentation with Documenter from a builder to GitHub Pages.

By default the links in the instructions need to be modified to correspond to actual URLs.

    documenter_keygen(; user="USER", repo="REPO")

The optional `user` and `repo` keyword arguments can be specified so that the URLs in the printed instructions could be copied directly. They should be the name of the GitHub user or organization where the repository is hosted and the full name of the repository,
respectively.

# Examples
```julia-repl
julia> using OpenSSH
julia> documenter_keygen()
[ Info: add the public key below to https://github.com/USER/REPO/settings/keys with read/write access:
ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname
[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.com/USER/REPO/settings (if you deploy using Travis CI) or https://github.com/USER/REPO/settings/secrets (if you deploy using GitHub Actions) with value:
LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==
julia> documenter_keygen(user="JuliaDocs", repo="OpenSSH.jl")
[Info: add the public key below to https://github.com/JuliaDocs/OpenSSH.jl/settings/keys with read/write access:
ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname
[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.com/JuliaDocs/OpenSSH.jl/settings (if you deploy using Travis CI) or https://github.com/JuliaDocs/OpenSSH.jl/settings/secrets (if you deploy using GitHub Actions) with value:
LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==
```
"""
function documenter_keygen(; user="USER", repo="REPO")
    # Error checking. Do the required programs exist?
    if Sys.iswindows()
        success(`where where`)      || error("'where' not found.")
        success(`where ssh-keygen`) || error("'ssh-keygen' not found.")
    else
        success(`which which`)      || error("'which' not found.")
        success(`which ssh-keygen`) || error("'ssh-keygen' not found.")
    end

    directory = pwd()
    filename  = "documenter-private-key"

    isfile(filename) && error("temporary file '$(filename)' already exists in working directory")
    isfile("$(filename).pub") && error("temporary file '$(filename).pub' already exists in working directory")

    # Generate the ssh key pair.
    success(`ssh-keygen -N "" -f $filename`) || error("failed to generate a SSH key pair.")

    # Prompt user to add public key to github then remove the public key.
    let url = "https://github.com/$user/$repo/settings/keys"
        @info("add the public key below to $url with read/write access:")
        println("\n", read("$filename.pub", String))
        rm("$filename.pub")
    end

    # Base64 encode the private key and prompt user to add it to travis. The key is
    # *not* encoded for the sake of security, but instead to make it easier to
    # copy/paste it over to travis without having to worry about whitespace.
    let travis_url = "https://travis-ci.com/$user/$repo/settings",
        github_url = "https://github.com/$user/$repo/settings/secrets"
        @info("add a secure environment variable named 'DOCUMENTER_KEY' to " *
              "$(travis_url) (if you deploy using Travis CI) or " *
              "$(github_url) (if you deploy using GitHub Actions) with value:")
        println("\n", base64encode(read(filename, String)), "\n")
        rm(filename)
    end
end

"""
    documenter_keygen(package::Module; remote="origin")

This method attempts to guess the package URLs from the Git remote.

`package` needs to be the top level module of the package. The `remote` keyword argument can be used to specify which Git remote is used for guessing the repository's GitHub URL.

This method requires `git` to be available from the command line.

!!! note
    the package must be in development mode. Make sure you run `pkg> develop pkg` from the Pkg REPL, or `Pkg.develop(\"pkg\")` before generating the SSH keys.

# Examples
```julia-repl
julia> using OpenSSH
julia> documenter_keygen(Keygen)
[Info: add the public key below to https://github.com/JuliaDocs/OpenSSH.jl/settings/keys with read/write access:
ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname
[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.com/JuliaDocs/OpenSSH.jl/settings (if you deploy using Travis CI) or https://github.com/JuliaDocs/OpenSSH.jl/settings/secrets (if you deploy using GitHub Actions) with value:
LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==
```
"""
function documenter_keygen(package::Module; remote="origin")
    # Error checking. Do the required programs exist?
    if Sys.iswindows()
        success(`where where`)      || error("'where' not found.")
        success(`where ssh-keygen`) || error("'ssh-keygen' not found.")
        success(`where git`)        || error("'git' not found.")
    else
        success(`which which`)      || error("'which' not found.")
        success(`which ssh-keygen`) || error("'ssh-keygen' not found.")
        success(`which git`)        || error("'git' not found.")
    end

    path = package_devpath(package)

    # Are we in a git repo?
    user, repo = cd(path) do
        success(`git status`) || error("Failed to run `git status` in $(path). 'Keygen.documenter_keygen' only works with Git repositories.")

        let r = readchomp(`git config --get remote.$remote.url`)
            m = match(GITHUB_REGEX, r)
            m === nothing && error("no remote repo named '$remote' found.")
            m[2], m[3]
        end
    end

    # Generate the ssh key pair.
    documenter_keygen(; user=user, repo=repo)
end


text:
module CompatHelper

using Base64: Base64
using Dates: Dates
using GitForge: GitForge, Forge, GitHub, GitLab
using MultilineStrings
using Mocking: Mocking, @mock
using Pkg: Pkg
using TimeZones: TimeZones
using TOML: TOML
using UUIDs: UUIDs, UUID

export DropEntry, KeepEntry, NewEntry
export CIService, GitHubActions, GitLabCI

@static if Base.VERSION >= v"1.7-"
    const VersionSpec = Pkg.Versions.VersionSpec
    const semver_spec = Pkg.Versions.semver_spec
else
    const VersionSpec = Pkg.Types.VersionSpec
    const semver_spec = Pkg.Types.semver_spec
end

const PRIVATE_SSH_ENVVAR = "COMPATHELPER_PRIV"

include(joinpath("utilities", "types.jl"))
include(joinpath("utilities", "utilities.jl"))
include(joinpath("utilities", "ci.jl"))
include(joinpath("utilities", "git.jl"))
include(joinpath("utilities", "ssh.jl"))
include(joinpath("utilities", "timestamps.jl"))
include(joinpath("utilities", "new_versions.jl"))

include("main.jl")
include("exceptions.jl")
include("dependencies.jl")
include("pull_requests.jl")

end # module


text:
# Common types and functions

type Rule
    P::Array{Int64} # Antecedent
    Q::Array{Int64} # Consequent
end

# Support Count: σ(x) = | {tᵢ|x ⊆ tᵢ,tᵢ∈ T}|
function σ(x, T)
    ret = 0
    for t in T
        ⊆(x,t) && (ret += 1)
    end
    ret
end

# Support of itemset x -> y, which x does not intersect y.
supp(x,y,T) = σ(∪(x,y),T)/length(T)

# Confidence of itemset x-> y, which x does not intersect y.
conf(x,y,T) = σ(∪(x,y),T)/σ(x,T)


text:
export readplist_string, readplist

function parse_obj(p::Parser)
    token = next_token(p)

    if     NUMBER == token.kind
        Meta.parse(token.lexeme)
    elseif STRING == token.kind || IDENT == token.kind
        token.lexeme
    elseif HEXBINARY == token.kind
        parse_hexbinary(token.lexeme)
    elseif LPAREN == token.kind
        backup_token(p)
        parse_array(p)
    elseif LBRACE == token.kind
        backup_token(p)
        parse_dict(p)
    end
end

parse_hexbinary(s::AbstractString) = hex2bytes(filter(isxdigit, s))

function parse_array(p::Parser)
    array = []

    expect(p, LPAREN)
    token = peek_token(p)
    if token.kind == RPAREN
        return array
    end
    push!(array, parse_obj(p))
    token = peek_token(p)
    while token.kind != RPAREN
        expect(p, COMMA)
        push!(array, parse_obj(p))
        token = peek_token(p)
    end
    expect(p, RPAREN)
    array
end

function parse_dict(p::Parser)
    dict =  Dict{Any, Any}()

    expect(p, LBRACE)
    token = peek_token(p)
    while token.kind != RBRACE
        key = parse_obj(p)::String
        expect(p, EQUAL)
        dict[key] = parse_obj(p)
        expect(p, SEMICOLON)
        token = peek_token(p)
    end
    expect(p, RBRACE)
    dict
end

function readplist_string(text::AbstractString)
    l = lex_plist(text)
    p = Parser(l)
    parse_obj(p)
end

function readplist(stream::IO)
    text = read(stream, String)
    readplist_string(text)
end

function readplist(filename::AbstractString)
    open(readplist, filename)
end


text:
using AdventOfCode2020
using Test

inputs_dir = "../Inputs"
dates = readdir(inputs_dir)

dates = ["Dec-1","Dec-2","Dec-3","Dec-4","Dec-5","Dec-6","Dec-7","Dec-8","Dec-9"]
results = Dict("Dec-1" => [989824, 66432240], "Dec-2" => [546,275], "Dec-3" => [280,4355551200], "Dec-4" => [192,101], "Dec-5" => [896,659], "Dec-6" => [6542,3299], "Dec-7" => [205,80902], "Dec-8" => [2058,1000], "Dec-9" => [144381670,20532569])
@testset "AdventOfCode" begin
    @testset "Date - $(date)" for date in dates
        @test AdventOfCode2020.run(1, date, joinpath(inputs_dir, date, "input.txt")) == results[date][1]
        @test AdventOfCode2020.run(2, date, joinpath(inputs_dir, date, "input.txt")) == results[date][2]
    end
end