text:
using AbstractPlotting: get_texture_atlas, glyph_bearing!, glyph_uv_width!, glyph_scale!, calc_position, calc_offset

gpuvec(x) = GPUVector(GLBuffer(x))

to_range(x, y) = to_range.((x, y))
to_range(x::ClosedInterval) = (minimum(x), maximum(x))
to_range(x::VecTypes{2}) = x
to_range(x::AbstractRange) = (minimum(x), maximum(x))
to_range(x::AbstractVector) = (minimum(x), maximum(x))

function to_range(x::AbstractArray)
    if length(x) in size(x) # assert that just one dim != 1
        to_range(vec(x))
    else
        error("Can't convert to a range. Please supply a range/vector/interval or a tuple (min, max)")
    end
end

function to_glvisualize_key(k)
    k == :rotations && return :rotation
    k == :markersize && return :scale
    k == :glowwidth && return :glow_width
    k == :glowcolor && return :glow_color
    k == :strokewidth && return :stroke_width
    k == :strokecolor && return :stroke_color
    k == :positions && return :position
    k == :linewidth && return :thickness
    k == :marker_offset && return :offset
    k == :colormap && return :color_map
    k == :colorrange && return :color_norm
    k == :transform_marker && return :scale_primitive
    return k
end

make_context_current(screen::Screen) = GLFW.MakeContextCurrent(to_native(screen))

function cached_robj!(robj_func, screen, scene, x::AbstractPlot)
    # poll inside functions to make wait on compile less prominent
    GLFW.PollEvents()
    robj = get!(screen.cache, objectid(x)) do
        filtered = filter(x.attributes) do (k, v)
            !(k in (:transformation, :tickranges, :ticklabels, :raw))
        end
        gl_attributes = Dict{Symbol, Any}(map(filtered) do key_value
            key, value = key_value
            gl_key = to_glvisualize_key(key)
            gl_value = lift_convert(key, value, x)
            gl_key => gl_value
        end)
        if haskey(gl_attributes, :scale)
            gl_attributes[:use_pixel_marker] = lift(x-> x isa Vec{2, <:AbstractPlotting.Pixel}, gl_attributes[:scale])
            gl_attributes[:scale] = lift(x-> AbstractPlotting.number.(x), gl_attributes[:scale])
            if haskey(gl_attributes, :offset)
                gl_attributes[:offset] = lift(x-> AbstractPlotting.number.(x), gl_attributes[:offset])
            end
        end
        if haskey(gl_attributes, :lightposition)
            eyepos = scene.camera.eyeposition
            gl_attributes[:lightposition] = lift(gl_attributes[:lightposition], eyepos) do pos, eyepos
                return pos == :eyeposition ? eyepos : pos
            end
        end
        robj = robj_func(gl_attributes)
        for key in (:pixel_space, :view, :projection, :resolution, :eyeposition, :projectionview)
            robj[key] = getfield(scene.camera, key)
        end
        screen.cache2plot[robj.id] = x
        robj
    end
    push!(screen, scene, robj)
    robj
end

function remove_automatic!(attributes)
    filter!(attributes) do (k, v)
        to_value(v) != automatic
    end
end

index1D(x::SubArray) = parentindices(x)[1]

handle_view(array::AbstractVector, attributes) = array
handle_view(array::Node, attributes) = array

function handle_view(array::SubArray, attributes)
    A = parent(array)
    indices = index1D(array)
    attributes[:indices] = indices
    return A
end

function handle_view(array::Node{T}, attributes) where T <: SubArray
    A = lift(parent, array)
    indices = lift(index1D, array)
    attributes[:indices] = indices
    A
end

function lift_convert(key, value, plot)
    lift(value) do value
        return convert_attribute(value, Key{key}(), Key{AbstractPlotting.plotkey(plot)}())
    end
end

pixel2world(scene, msize::Number) = pixel2world(scene, Point2f0(msize))[1]

function pixel2world(scene, msize::StaticVector{2})
    # TODO figure out why Vec(x, y) doesn't work correctly
    p0 = AbstractPlotting.to_world(scene, Point2f0(0.0))
    p1 = AbstractPlotting.to_world(scene, Point2f0(msize))
    diff = p1 - p0
    return diff
end

pixel2world(scene, msize::AbstractVector) = pixel2world.(scene, msize)

function handle_intensities!(attributes)
    if haskey(attributes, :color) && attributes[:color][] isa AbstractVector{<: Number}
        c = pop!(attributes, :color)
        attributes[:intensity] = lift(x-> convert(Vector{Float32}, x), c)
    else
        delete!(attributes, :intensity)
        delete!(attributes, :color_map)
        delete!(attributes, :color_norm)
    end
end

function Base.insert!(screen::GLScreen, scene::Scene, x::Combined)
    # poll inside functions to make wait on compile less prominent
    GLFW.PollEvents()
    if isempty(x.plots) # if no plots inserted, this truely is an atomic
        draw_atomic(screen, scene, x)
    else
        foreach(x.plots) do x
            # poll inside functions to make wait on compile less prominent
            GLFW.PollEvents()
            insert!(screen, scene, x)
        end
    end
end

struct FastPixel end
AbstractPlotting.to_spritemarker(x::FastPixel) = x

function draw_atomic(screen::GLScreen, scene::Scene, x::Union{Scatter, MeshScatter})
    robj = cached_robj!(screen, scene, x) do gl_attributes
        # signals not supported for shading yet
        gl_attributes[:shading] = to_value(get(gl_attributes, :shading, true))
        marker = lift_convert(:marker, pop!(gl_attributes, :marker), x)
        if isa(x, Scatter)
            gl_attributes[:billboard] = map(rot-> isa(rot, Billboard), x.rotations)
            gl_attributes[:distancefield][] == nothing && delete!(gl_attributes, :distancefield)
            gl_attributes[:uv_offset_width][] == Vec4f0(0) && delete!(gl_attributes, :uv_offset_width)
        end
        positions = handle_view(x[1], gl_attributes)
        if marker[] isa FastPixel
            filter!(gl_attributes) do (k, v,)
                k in (:color_map, :color, :color_norm, :scale, :fxaa, :model)
            end
            if !(gl_attributes[:color][] isa AbstractVector{<: Number})
                delete!(gl_attributes, :color_norm)
                delete!(gl_attributes, :color_map)
            end
            visualize(positions, Style(:speed), Dict{Symbol, Any}(gl_attributes))
        else
            handle_intensities!(gl_attributes)
            visualize((marker, positions), Style(:default), Dict{Symbol, Any}(gl_attributes))
        end
    end
end

function draw_atomic(screen::GLScreen, scene::Scene, x::Lines)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        linestyle = pop!(gl_attributes, :linestyle)
        data = Dict{Symbol, Any}(gl_attributes)
        data[:pattern] = to_value(linestyle)
        positions = handle_view(x[1], data)
        handle_intensities!(data)
        visualize(positions, Style(:lines), data)
    end
end

function draw_atomic(screen::GLScreen, scene::Scene, x::LineSegments)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        linestyle = pop!(gl_attributes, :linestyle)
        data = Dict{Symbol, Any}(gl_attributes)
        data[:pattern] = to_value(linestyle)
        positions = handle_view(x.converted[1], data)
        if haskey(data, :color) && data[:color][] isa AbstractVector{<: Number}
            c = pop!(data, :color)
            data[:color] = lift(AbstractPlotting.el32convert, c)
        else
            delete!(data, :color_map)
            delete!(data, :color_norm)
        end
        visualize(positions, Style(:linesegment), data)
    end
end

function to_gl_text(string, startpos::AbstractVector{T}, textsize, font, align, rot, model) where T <: VecTypes
    atlas = get_texture_atlas()
    N = length(T)
    positions, uv_offset_width, scale = Point{3, Float32}[], Vec4f0[], Vec2f0[]
    # toffset = calc_offset(string, textsize, font, atlas)
    char_str_idx = iterate(string)
    broadcast_foreach(1:length(string), startpos, textsize, (font,), align) do idx, pos, tsize, font, align
        char, str_idx = char_str_idx
        _font = isa(font[1], NativeFont) ? font[1] : font[1][idx]
        mpos = model * Vec4f0(to_ndim(Vec3f0, pos, 0f0)..., 1f0)
        push!(positions, to_ndim(Point{3, Float32}, mpos, 0))
        push!(uv_offset_width, glyph_uv_width!(atlas, char, _font))
        if isa(tsize, Vec2f0) # this needs better unit support
            push!(scale, tsize) # Vec2f0, we assume it's already in absolute size
        else
            push!(scale, glyph_scale!(atlas, char,_font, tsize))
        end
        char_str_idx = iterate(string, str_idx)
    end
    positions, Vec2f0(0), uv_offset_width, scale
end

function to_gl_text(string, startpos::VecTypes{N, T}, textsize, font, aoffsetvec, rot, model) where {N, T}
    atlas = get_texture_atlas()
    mpos = model * Vec4f0(to_ndim(Vec3f0, startpos, 0f0)..., 1f0)
    pos = to_ndim(Point{3, Float32}, mpos, 0f0)
    rscale = Float32(textsize)
    chars = Vector{Char}(string)
    scale = glyph_scale!.(Ref(atlas), chars, (font,), rscale)
    positions2d = calc_position(string, Point2f0(0), rscale, font, atlas)
    # font is Vector{FreeType.NativeFont} so we need to protec
    aoffset = AbstractPlotting.align_offset(Point2f0(0), positions2d[end], atlas, rscale, font, aoffsetvec)
    aoffsetn = to_ndim(Point{3, Float32}, aoffset, 0f0)
    uv_offset_width = glyph_uv_width!.(Ref(atlas), chars, (font,))
    positions = map(positions2d) do p
        pn = rot * (to_ndim(Point{3, Float32}, p, 0f0) .+ aoffsetn)
        pn .+ pos
    end
    positions, Vec2f0(0), uv_offset_width, scale
end

function draw_atomic(screen::GLScreen, scene::Scene, x::Text)
    robj = cached_robj!(screen, scene, x) do gl_attributes

        liftkeys = (:position, :textsize, :font, :align, :rotation, :model)

        gl_text = lift(to_gl_text, x[1], getindex.(Ref(gl_attributes), liftkeys)...)
        # unpack values from the one signal:
        positions, offset, uv_offset_width, scale = map((1, 2, 3, 4)) do i
            lift(getindex, gl_text, i)
        end

        atlas = get_texture_atlas()
        keys = (:color, :stroke_color, :stroke_width, :rotation)
        signals = getindex.(Ref(gl_attributes), keys)
        visualize(
            (DISTANCEFIELD, positions),
            color = signals[1],
            stroke_color = signals[2],
            stroke_width = signals[3],
            rotation = signals[4],
            scale = scale,
            offset = offset,
            uv_offset_width = uv_offset_width,
            distancefield = get_texture!(atlas),
            visible = gl_attributes[:visible]
        )
    end
end

function draw_atomic(screen::GLScreen, scene::Scene, x::Heatmap)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        gl_attributes[:ranges] = lift(to_range, x[1], x[2])
        interp = to_value(pop!(gl_attributes, :interpolate))
        interp = interp ? :linear : :nearest
        tex = Texture(x[3], minfilter = interp)
        pop!(gl_attributes, :color)
        if haskey(gl_attributes, :nan_color)
            gl_attributes[:nan_color] = lift(to_color, gl_attributes[:nan_color])
        end
        gl_attributes[:stroke_width] = pop!(gl_attributes, :thickness)
        GLVisualize.assemble_shader(GLVisualize.gl_heatmap(tex, gl_attributes))
    end
end

function vec2color(colors, cmap, crange)
    AbstractPlotting.interpolated_getindex.((to_colormap(cmap),), colors, (crange,))
end

function get_image(plot)
    if isa(plot[:color][], AbstractMatrix{<: Number})
        lift(vec2color, pop!.(Ref(plot), (:color, :color_map, :color_norm))...)
    else
        delete!(plot, :color_norm)
        delete!(plot, :color_map)
        return pop!(plot, :color)
    end
end

function draw_atomic(screen::GLScreen, scene::Scene, x::Image)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        gl_attributes[:ranges] = lift(to_range, x[1], x[2])
        img = get_image(gl_attributes)
        # remove_automatic!(gl_attributes)
        visualize(img, Style(:default), gl_attributes)
    end
end

convert_mesh_color(c::AbstractArray{<: Number}, cmap, crange) = vec2color(c, cmap, crange)
convert_mesh_color(c, cmap, crange) = c

function draw_atomic(screen::GLScreen, scene::Scene, x::Mesh)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        # signals not supported for shading yet
        gl_attributes[:shading] = to_value(pop!(gl_attributes, :shading))
        color = pop!(gl_attributes, :color)
        cmap = get(gl_attributes, :color_map, Node(nothing)); delete!(gl_attributes, :color_map)
        crange = get(gl_attributes, :color_norm, Node(nothing)); delete!(gl_attributes, :color_norm)
        if to_value(color) isa Colorant
            gl_attributes[:vertex_color] = color
        end
        if to_value(color) isa AbstractMatrix{<:Colorant}
            gl_attributes[:image] = color
        end
        mesh = x[1]
        if to_value(color) isa AbstractVector{<: Number}
            mesh = lift(x[1], color, cmap, crange) do mesh, color, cmap, crange
                color_sampler = AbstractPlotting.sampler(cmap, color, crange)
                GeometryBasics.pointmeta(mesh, color=color_sampler)
            end
        end

        if to_value(color) isa AbstractMatrix{<: Number}
            mesh = lift(x[1], color, cmap, crange) do mesh, color, cmap, crange
                color_sampler = convert_mesh_color(color, cmap, crange)
                mesh, uv = GeometryBasics.pop_pointmeta(mesh, :uv)
                uv_sampler = AbstractPlotting.sampler(color_sampler, uv)
                GeometryBasics.pointmeta(mesh, color=uv_sampler)
            end
        end

        if to_value(color) isa AbstractVector{<: Colorant}
            mesh = lift(x[1], color, cmap, crange) do mesh, color, cmap, crange
                GeometryBasics.pointmeta(mesh, color=color)
            end
        end

        visualize(mesh, Style(:default), gl_attributes)
    end
end

function draw_atomic(screen::GLScreen, scene::Scene, x::Surface)
    robj = cached_robj!(screen, scene, x) do gl_attributes
        color = pop!(gl_attributes, :color)
        img = nothing
        # signals not supported for shading yet
        # We automatically insert x[3] into the color channel, so if it's equal we don't need to do anything
        if isa(to_value(color), AbstractMatrix{<: Number}) && to_value(color) !== to_value(x[3])
            crange = pop!(gl_attributes, :color_norm)
            cmap = pop!(gl_attributes, :color_map)
            img = lift(color, cmap, crange) do img, cmap, norm
                AbstractPlotting.interpolated_getindex.((cmap,), img, (norm,))
            end
        elseif isa(to_value(color), AbstractMatrix{<: Colorant})
            img = color
            gl_attributes[:color_map] = nothing
            gl_attributes[:color] = nothing
            gl_attributes[:color_norm] = nothing
        end
        gl_attributes[:image] = img
        args = x[1:3]
        gl_attributes[:shading] = to_value(get(gl_attributes, :shading, true))
        if all(v-> to_value(v) isa AbstractMatrix, args)
            return visualize(args, Style(:surface), gl_attributes)
        else
            gl_attributes[:ranges] = to_range.(to_value.(args[1:2]))
            return visualize(args[3], Style(:surface), gl_attributes)
        end
    end
    return robj
end

function draw_atomic(screen::GLScreen, scene::Scene, vol::Volume)
    robj = cached_robj!(screen, scene, vol) do gl_attributes
        model = vol[:model]
        x, y, z = vol[1], vol[2], vol[3]
        gl_attributes[:model] = lift(model, x, y, z) do m, xyz...
            mi = minimum.(xyz)
            maxi = maximum.(xyz)
            w = maxi .- mi
            m2 = Mat4f0(
                w[1], 0, 0, 0,
                0, w[2], 0, 0,
                0, 0, w[3], 0,
                mi[1], mi[2], mi[3], 1
            )
            return convert(Mat4f0, m) * m2
        end
        return visualize(vol[4], Style(:default), gl_attributes)
    end
end


text:
export run_keycode_helper, run_calibrate, create_new_project
using DataStructures
using Lazy: @>
"""
    run_keycode_helper(;extensions=[])

Runs a program that will display the keycode for each key that you press.
"""
function run_keycode_helper(;extensions=Extension[])
  result = @spawn run_keycode_helper__(extensions=extensions)
  fetch(result)
end

function run_keycode_helper__(;extensions=Extension[])
  exp = Experiment(hide_output=true,extensions=extensions)
  setup(exp) do
    addbreak(instruct("Press keys to see their codes."))
    addtrial(show_cross(),await_response(iskeydown(key":escape:"))) do event
      if iskeydown(event)
        display(visual(string(keycode(event))))
      end
    end
  end

  run(exp,await_input=false)
end

"""
   run_calibrate()

Runs a program that will allow you to play pure tones and adjust their level.

This program provides one means of calibrating the levels of sound in
your experiment. Using a sound-level meter you can determine the dB SPL of
each tone, and adjust the attenuation to achieve a desired sound level.
"""
function run_calibrate()
  result = @spawn run_calibrate__()
  fetch(result)
end

function run_calibrate__()
  exp = Experiment(hide_output=true)
  freqs = [100,250,500,1000,2000,4000,6000,8000]
  atten = collect(20.0 for f in freqs)
  old_tone = nothing

  playing = 1
  key_to_tone(e::KeyDownEvent) = key_to_tone(keycode(e))
  key_to_tone(e::KeyboardKey) = e.code - key"1".code + 1
  key_to_tone(e) = 0

  tone_text() = visual("$(freqs[playing]) Hz, $(-round(atten[playing],2))"*
                       " dB tone",y=-0.5)

  setup(exp) do
    instructions = moment() do
      @> tone(freqs[playing]) attenuate(atten[playing]) fadeto
      display(visual("Hit 1-8 to play a tone.",y=0.75,duration=Inf) +
              visual("Hit - or = to adjust dB by 10 dB\n"*
                     "Hit [ or ] to adjust dB by 1 dB\n"*
                     "Hit ; or ' to adjust dB by 0.1 dB\n"*
                     "Hit , or . to adjust dB by 0.01 dB\n",
                     clean_whitespace=false,y=0.25,duration=Inf) +
              visual("$(freqs[1]) Hz, $(-atten[1]) dB tone",y=-0.5))
    end
    atten_adjust = Dict(key"-" => -10,   key"=" => +10,
                        key"[" => -1,    key"]" => +1,
                        key";" => -0.1,  key"'" => +0.1,
                        key"," => -0.01, key"." => +0.01)

    addtrial(instructions,await_response(iskeydown(key":escape:"))) do e
      if iskeydown(e)
        if 1 <= (p = key_to_tone(e)) <= 8
          playing = p
          @> tone(freqs[playing]) attenuate(atten[playing]) fadeto
          display(tone_text())
        elseif keycode(e) ∈ keys(atten_adjust)
          atten[playing] -= atten_adjust[keycode(e)]
          @> tone(freqs[playing]) attenuate(atten[playing]) fadeto
          display(tone_text())
        end
      end
    end
  end

  run(exp,await_input=false)
  return OrderedDict(f => atten for (f,atten) in zip(freqs,atten))
end

"""
    create_new_project(name,dir=".")

Creates a set of files to help you get started on a new experiment.

This creates a file called run_[name].jl, and a README.md and setup.jl file for
your experiment. The readme provides useful information for running the
experiment that is common across all experiments. The run file provides some
guidelines to get you started creating an experiment and the setup file is a
script that can be used to install Weber and any additional dependencies for the
project, for anyone who wants to download and run your experiment.
"""

function create_new_project(name,dir=".")
  values = Dict(r"{{project}}" => name,
                r"{{weber_dir}}" => abspath(joinpath(dirname(@__FILE__),"..")),
                r"{{version}}" => "v\"$(Weber.version)\"")

  apply_template(joinpath(dirname(@__FILE__),"..","templates","README.md"),
                 joinpath(dir,"README.md"),values)
  apply_template(joinpath(dirname(@__FILE__),"..","templates","setup.jl"),
                 joinpath(dir,"setup.jl"),values)
  apply_template(joinpath(dirname(@__FILE__),"..","templates","run_project.jl"),
                 joinpath(dir,"run_$name.jl"),values)
  apply_template(joinpath(dirname(@__FILE__),"..","templates",".gitignore"),
                 joinpath(dir,".gitignore"),values)
end

function apply_template(sourcefile,destfile,values)
  open(sourcefile,"r") do source
    open(destfile,"w") do dest
      for line in readlines(source)
        for (pat,result) in values
          line = replace(line,pat,result)
        end
        write(dest,line)
      end
    end
  end
end

ext: jl


text:
using IndexFunArrays
using BenchmarkTools

function rr_test(s)
    x = ones(s)
    y = rr2(x) .+ sqrt.(1.2 .* rr2(x, scale=(1,1), offset=CtrCorner));
    @btime $y .= rr2($x) .+ sqrt.(1.2 .* rr2($x, scale=(1,1), offset=CtrCorner));
    @btime $y .= rr2($x) .+ sqrt.(1.2 .* rr2($x, scale=(1,1), offset=CtrCorner));
    @btime $y .= rr2($x) .+ sqrt.(1.2 .* rr2($x, scale=(1,1), offset=CtrCorner));
    return
end

function xx_test(s)
    x = ones(s)
    y = xx(x) .+ sqrt.(1.2 .* abs.(xx(x)));
    @btime $y .= xx($x) .+ sqrt.(1.2 .* abs.(xx($x)));
    @btime $y .= xx($x) .+ sqrt.(1.2 .* abs.(xx($x)));
    return
end

function compare_to_CartesianIndices()
    x = ones(1000,1000);
    qq2(index) = sum(Tuple(index));
    ww2 = qq2.(CartesianIndices(x));
    y = rr2(x) .+ sqrt.(1.2.*rr2(x).*rr2(x));

    @info "rr2 based"
    @btime $y .= rr2($x) .+ sqrt.(1.2.*rr2($x).*rr2($x));
    @btime $y .= rr2($x) .+ sqrt.(1.2.*rr2($x).*rr2($x));
    @info "CartesianIndices based"
    @btime $y .= ($qq2).(CartesianIndices($x)) .+ sqrt.(1.2.*($qq2).(CartesianIndices($x)).*($qq2).(CartesianIndices($x)));
    @btime $y .= ($qq2).(CartesianIndices($x)) .+ sqrt.(1.2.*($qq2).(CartesianIndices($x)).*($qq2).(CartesianIndices($x)));
    @info "CartesianIndices based with initialized function"
    @btime $y .= $ww2 .+ sqrt.(1.2.*$ww2.*$ww2);
    @btime $y .= $ww2 .+ sqrt.(1.2.*$ww2.*$ww2);
end


text:
using Flux.Optimise: update!
using Flux.Losses: mse, mae
using Statistics
using Zygote

include("generator.jl")
include("discriminator.jl")

export Discriminator
export Generator

struct Cyclegan
    gen_AB
    gen_BA
    disc_AB
    disc_BA
end

@functor Cyclegan

function calculate_loss_discr(real, fake)
    return (real + fake) / 2
end

function calculate_loss_gen(loss_G_A, loss_G_B, cycle_A_loss,
    cycle_B_loss)
    return loss_G_A
        + loss_G_B
        + cycle_A_loss * 10
        + cycle_B_loss * 10
end

function train_discr(cycleGAN, original_A, original_B,
                     fake_A, fake_B, opt_discr)
    ps = Flux.params(cycleGAN.disc_AB, cycleGAN.disc_BA)
    loss, back = Zygote.pullback(ps) do
                        #calculate A
                        D_A_real = cycleGAN.disc_AB(original_A)
                        D_A_fake = cycleGAN.disc_AB(fake_A)
                        D_A_real_loss = mse(D_A_real, ones(size(D_A_real)))
                        D_A_fake_loss = mse(D_A_fake, zeros(size(D_A_fake)))
                        D_A_loss = D_A_real_loss + D_A_fake_loss

                        #calculate A
                        D_B_real = cycleGAN.disc_BA(original_B)
                        D_B_fake = cycleGAN.disc_BA(fake_B)
                        D_B_real_loss = mse(D_B_real, ones(size(D_B_real)))
                        D_B_fake_loss = mse(D_B_fake, zeros(size(D_B_fake)))
                        D_B_loss = D_B_real_loss + D_B_fake_loss
                        calculate_loss_discr(D_A_loss, D_B_loss)
    end
    grads = back(1.f0)
    update!(opt_discr, ps, grads)

    return loss
end

Zygote.@nograd train_discr

function train_gan(cycleGAN, original_A, original_B, opt_gen, opt_discr)
    loss = Dict()
    ps = Flux.params(cycleGAN.gen_AB, cycleGAN.gen_BA)
    loss["G_loss"], back = Zygote.pullback(ps) do
                            fake_A = cycleGAN.gen_AB(original_B)
                            fake_B = cycleGAN.gen_BA(original_A)
                            loss["D_loss"]= train_discr(cycleGAN,
                                                        original_A, original_B,
                                                        fake_A, fake_B, opt_discr)
                            # adversarial loss for both generators
                            D_A_fake = cycleGAN.disc_AB(fake_A)
                            D_B_fake = cycleGAN.disc_BA(fake_B)
                            loss_G_A = mse(D_A_fake, ones(size(D_A_fake)))
                            loss_G_B = mse(D_B_fake, ones(size(D_B_fake)))

                            #cycle loss
                            cycle_B = cycleGAN.gen_BA(fake_A)
                            cycle_A = cycleGAN.gen_AB(fake_B)
                            cycle_B_loss = mae(original_B, cycle_B)
                            cycle_A_loss = mae(original_A, cycle_A)

                            # identity loss (remove these for efficiency if you set lambda_identity=0)
                            #identity_B = gen_B(original_B)
                            #identity_A = gen_A(original_A)
                            #identity_B_loss = mae(original_B, identity_B)
                            #identity_A_loss = mae(original_A, identity_A)

                            calculate_loss_gen(loss_G_A, loss_G_B, cycle_A_loss, cycle_B_loss)
    end
    grads = back(1.f0)
    update!(opt_gen, ps, grads)
    return loss
end


text:
## Lifted operators

import Base: promote_op, scalarmin, scalarmax
using Compat

if isdefined(Base, :fieldname) && Base.fieldname(Nullable, 1) == :hasvalue # Julia 0.6
    _Nullable(R, x, hasvalue::Bool) = Nullable{R}(x, hasvalue)
else
    _Nullable(R, x, hasvalue::Bool) = Nullable{R}(x, !hasvalue)
end

"""
    null_safe_op(f::Any, ::Type...)::Bool

Returns whether an operation `f` can safely be applied to any value of the passed type(s).
Returns `false` by default.

Custom types should implement methods for some or all operations `f` when applicable:
returning `true` means that the operation may be called on any bit pattern without
throwing an error (though returning invalid or nonsensical results is not a problem).
In particular, this means that the operation can be applied on the whole domain of the
type *and on uninitialized objects*. As a general rule, these proporties are only true for
safe operations on `isbits` types.

Types declared as safe can benefit from higher performance for operations on nullable: by
always computing the result even for null values, a branch is avoided, which helps
vectorization.
"""
null_safe_op(f::Any, ::Type...) = false

@compat const SafeSignedInts = Union{Int128,Int16,Int32,Int64,Int8}
@compat const SafeUnsignedInts = Union{Bool,UInt128,UInt16,UInt32,UInt64,UInt8}
@compat const SafeInts = Union{SafeSignedInts,SafeUnsignedInts}
@compat const SafeFloats = Union{Float16,Float32,Float64}

# Float types appear in both since they promote to themselves,
# and therefore can't fail due to conversion of negative numbers
@compat const SafeSigned = Union{SafeSignedInts,SafeFloats}
@compat const SafeUnsigned = Union{SafeUnsignedInts,SafeFloats}
@compat const SafeTypes = Union{SafeInts,SafeFloats}

# Unary operators

# Note this list does not include sqrt since it can raise an error,
# nor cbrt (for which there is no functor on Julia 0.4)
for op in (:+, :-, :abs, :abs2)
    @eval begin
        null_safe_op{T<:SafeTypes}(::typeof($op), ::Type{T}) = true
    end
end

null_safe_op{T<:SafeInts}(::typeof(~), ::Type{T}) = true
null_safe_op{T<:SafeTypes}(::typeof(cbrt), ::Type{T}) = true
null_safe_op(::typeof(!), ::Type{Bool}) = true

# Temporary workaround until JuliaLang/julia#18803
null_safe_op(::typeof(cbrt), ::Type{Float16}) = false

for op in (:+, :-, :!, :~, :abs, :abs2, :sqrt, :cbrt)
    @eval begin
        @inline function Base.$op{S}(x::Nullable{S})
            R = promote_op($op, S)
            if null_safe_op($op, S)
                _Nullable(R, $op(x.value), !isnull(x))
            else
                isnull(x) ? Nullable{R}() :
                            Nullable{R}($op(x.value))
            end
        end
        Base.$op(x::Nullable{Union{}}) = Nullable()
    end
end

# Binary operators

# Note this list does not include ^, ÷ and %
# Operations between signed and unsigned types are not safe: promotion to unsigned
# gives an InexactError for negative numbers
for op in (:+, :-, :*, :/, :&, :|, :<<, :>>, :(>>>),
           :(==), :<, :>, :<=, :>=,
           :scalarmin, :scalarmax,
           :isless)
    @eval begin
        # to fix ambiguities
        null_safe_op{S<:SafeFloats,
                     T<:SafeFloats}(::typeof($op), ::Type{S}, ::Type{T}) = true
        null_safe_op{S<:SafeSigned,
                     T<:SafeSigned}(::typeof($op), ::Type{S}, ::Type{T}) = true
        null_safe_op{S<:SafeUnsigned,
                     T<:SafeUnsigned}(::typeof($op), ::Type{S}, ::Type{T}) = true
    end
end

for op in (:+, :-, :*, :/, :%, :÷, :&, :|, :^, :<<, :>>, :(>>>),
           :(==), :<, :>, :<=, :>=,
           :scalarmin, :scalarmax)
    @eval begin
        @inline function Base.$op{S,T}(x::Nullable{S}, y::Nullable{T})
            R = promote_op(Base.$op, S, T)
            if null_safe_op(Base.$op, S, T)
                _Nullable(R, Base.$op(x.value, y.value), !(isnull(x) | isnull(y)))
            else
                (isnull(x) | isnull(y)) ? Nullable{R}() :
                                          Nullable{R}(Base.$op(x.value, y.value))
            end
        end
        Base.$op(x::Nullable{Union{}}, y::Nullable{Union{}}) = Nullable()
        Base.$op{S}(x::Nullable{Union{}}, y::Nullable{S}) = Nullable{S}()
        Base.$op{S}(x::Nullable{S}, y::Nullable{Union{}}) = Nullable{S}()
    end
end

if !method_exists(isless, Tuple{Nullable{Int}, Nullable{Int}})
    function Base.isless{S,T}(x::Nullable{S}, y::Nullable{T})
        # NULL values are sorted last
        if null_safe_op(isless, S, T)
            (!isnull(x) & isnull(y)) |
            (!isnull(x) & !isnull(y) & isless(x.value, y.value))
        else
            if isnull(x)
                return false
            elseif isnull(y)
                return true
            else
                return isless(x.value, y.value)
            end
        end
    end
    Base.isless(x::Nullable{Union{}}, y::Nullable{Union{}}) = false
    Base.isless(x::Nullable{Union{}}, y::Nullable) = false
    Base.isless(x::Nullable, y::Nullable{Union{}}) = !isnull(x)
end


text:
import Base: +, -, *, <
using HCubature, ForwardDiff

# Function composition
(+)(a::Function, b::Function) = (x, y) -> (@inline a)(x, y) + (@inline b)(x, y)
(*)(a::Function, b::Function) = (x, y) -> (@inline a)(x, y) * (@inline b)(x, y)

(+)(a::Function, b::Number) = (x, y) -> (@inline a)(x, y) + b
(+)(a::Number, b::Function) = (x, y) -> (@inline b)(x, y) + a
(-)(a::Function, b::Number) = (x, y) -> (@inline a)(x, y) - b
(-)(a::Number, b::Function) = (x, y) -> (@inline b)(x, y) - a

# Cartesian position comparison
(<)(a::Tuple, b::Tuple) = all(a .< b)

# Derivatives
∂x(f) = (x, y) -> ForwardDiff.derivative(x -> f(x, y), x)
∂y(f) = (x, y) -> ForwardDiff.derivative(y -> f(x, y), y)

# Integrals
∫, ∫∫ = hquadrature, hcubature
∑ = sum

# Notice valid only for linear functions
∫dlX(f, a, b, y, slope) = ∫(x -> f(x, y(x)), a, b)[1] * sqrt(1 + slope^2) # helper
∫dlY(f, a, b, x, slope) = ∫(y -> f(x(y), y), a, b)[1] * sqrt(1 + slope^2) # helper

∫dl(f, e::EdgeX) = ∫dlX(f, e.range..., e.y, e.a)
∫dl(f, e::EdgeY) = ∫dlY(f, e.range..., e.x, e.a)
∫dl(f, edges) = sum(∫dl.(f, edges))

initdiv = trunc(Int, ceil(max(1/dx,1/dy, 1)))
∫∫dS(f, A, B) = hcubature(v -> f(v...), A, B, initdiv=initdiv)[1]

max_stars_repo_name: Pan-Maciek/FEM

max_issues_repo_name: Pan-Maciek/FEM

max_forks_repo_name: Pan-Maciek/FEM


text:
function assignStructure(df, tmArr)
	nrows = size(df, 1)
    df.element = repeat([""], nrows)
    df.structMap = zeros(nrows)
    df.resNum = 1:nrows

    init = 1
    for region in tmArr
        df.element[collect(region)] .= "tm"
        df.structMap[collect(region)] .= init

        init += 1
    end

    init = 0.5
    for row in eachrow(df)
        if row.structMap != 0
            init = row.structMap + 0.5
        else
            row.structMap = init
            row.element = "loop"
        end
    end

    nt = @view df[df.structMap .== minimum(df.structMap), :]
    ct = @view df[df.structMap .== maximum(df.structMap), :]
    nt[:,:element] .= "nterm"
    ct[:,:element] .= "cterm"

    return df
end


text:
"""
    is_temporal(x) -> StaticBool

Returns `True` if `x` refers to time.
"""
is_temporal(x::Symbol) = is_temporal(static(x))
is_temporal(x) = is_temporal(typeof(x))
is_temporal(::Type{T}) where {T} = static(false)
is_temporal(::Type{<:Union{StaticSymbol{:time},StaticSymbol{:Time}}}) = static(true)
is_temporal(::Type{<:Dates.AbstractTime}) = static(true)

function ArrayInterface.to_dims(::Type{T}, ::typeof(is_temporal)) where {T}
    d = _find_timedim(static(ndims(T)), dimnames(T))
    if d === nothing
        throw(ArgumentError("$T does not have a temporal dimension."))
    else
        return d
    end
end

_find_timedim(d::StaticInt{0}, n::Tuple) where {D} = nothing
@inline function _find_timedim(d::StaticInt{D}, n::Tuple) where {D}
    _find_timedim(is_temporal(@inbounds(getfield(n, D))), d, n)
end
_find_timedim(::True, d::StaticInt{D}, n::Tuple) where {D} = d
_find_timedim(::False, d::StaticInt{D}, n::Tuple) where {D} = _find_timedim(d - static(1), n)

"""
    has_timedim(x) -> Bool

Returns 'true' if `x` has a time dimension.
"""
has_timedim(x) = _find_timedim(static(ndims(x)), dimnames(x)) !== nothing

"""
    timedim(::Type{T}) -> Integer

Returns the dimension associated with time.
"""
timedim(x) = timedim(typeof(x))
timedim(::Type{T}) where {T} = ArrayInterface.to_dims(T, is_temporal)

"""
    ntimes(x) -> Integer

Returns the number of elements along the time dimension.
"""
ntimes(x) = size(x, timedim(x))

"""
    times(x)

Returns a collection of times associated with `x`. If `x` has a specific dimension
associated with time, this returns the keys along that dimension.
"""
@inline times(x) = keys(ArrayInterface.axes(x, is_temporal))

"""
    time_last(x)

Last time point along the time axis.
"""
time_last(x) = last(times(x))

"""
    time_first(x)

First time point along the time axis.
"""
time_first(x) = first(times(x))

"""
    time_step(x)

The time step/interval between each element.
"""
time_step(x) = step(times(x))

"""
    duration(x)

Duration of the event along the time axis.
"""
duration(x) = time_last(x) - time_first(x) + time_step(x)

"""
    sampling_rate(x)

Number of samples per second.
"""
sampling_rate(x) = inv(time_step(x))

"""
    assert_timedim_last(x)

Throw an error if the `x` has a time dimension that is not the last dimension.
"""
@inline function assert_timedim_last(x)
    if _find_timedim(static(ndims(x)), dimnames(x)) != ndims(x)
        throw(ArgumentError("time dimension is not last"))
    else
        return nothing
    end
end


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
np = pyimport("numpy")

"""
1  -  2  -  3  -  4  -  5

local node id: left to right
"""

ndofs = 2
nx = 4
nnodes, neles = (nx + 1), nx
ngp=2

NT = 1000
Δt = 0.05
T = NT*Δt

function BoundaryCondition(tid::Int64)
# PARAMETER: distance between nodes 0-1 and 0-2
l = 1.0
nodes = zeros(Float64, nnodes, ndofs)
for i = 1:nx+1
    nodes[i, 1], nodes[i, 2] = (i - 1)*l, 0.0
end

EBC, g = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)
#fix node 1
EBC[1, :] .= -1
FBC, fext = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)

#tid = 2 # pull from right with force control
#tid = 4 # pull from right with displacement control

gt = nothing
ft = nothing

if tid == 2
    FBC[nx+1, 1]  = -1
    fext[nx+1, 1] = 1.0e3
elseif tid == 4
    EBC[nx+1, 1] = -2
    EBC[nx+1, 2] = -1

    function ggt(t)
         vel = 0.2*1e-3#m/ms
   	 return ones(sum(EBC.==-2))*vel*t, ones(sum(EBC.==-2))*vel, zeros(sum(EBC.==-2))
    end
    gt = ggt
end

npoints = nnodes
node_to_point = collect(1:npoints)

return nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point
end


text:
# restore tagless final encoding from Core
export re_tagless

@nospecialize
function re_tagless(pi :: PatternInfo, ln :: Union{Nothing, LineNumberNode} = nothing)
    config = (type=pi.typetag, ln=ln)
    re_tagless(config, pi.pattern)
end

function re_tagless(config :: NamedTuple, p :: And)
    UserSitgnature.and([re_tagless(e) for e in p.ps], config)
end

function re_tagless(config :: NamedTuple, p::Or)
    UserSitgnature.or([re_tagless(e) for e in p.ps], config)
end

function re_tagless(config :: NamedTuple, p :: Literal)
    UserSitgnature.literal(p.val, config)
end

function re_tagless(config :: NamedTuple, p :: Wildcard)
    UserSitgnature.wildcard(config)
end

function re_tagless(config :: NamedTuple, p :: Deconstrucution)
    UserSitgnature.decons(p.comp::PComp, p.extract, [re_tagless(p) for p in p.params], config)
end

function re_tagless(config :: NamedTuple, p :: Guard)
    UserSitgnature.guard(p.predicate, config)
end

function re_tagless(config :: NamedTuple, p :: Effect)
    UserSitgnature.effect(p.perform, config)
end
@specialize

max_stars_repo_name: KristofferC/AbstractPattern

max_issues_repo_name: KristofferC/AbstractPattern

max_forks_repo_name: KristofferC/AbstractPattern


text:
# ------------------------------------------------------------------------------------------
# # Multiple dispatch
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In this notebook we'll explore **multiple dispatch**, which is a key feature of Julia.
#
# Multiple dispatch makes software *generic* and *fast*!
#
# #### Starting with the familiar
#
# To understand multiple dispatch in Julia, let's start with what we've already seen.
#
# We can declare functions in Julia without giving Julia any information about the types of
# the input arguments that function will receive:
# ------------------------------------------------------------------------------------------

f(x) = x^2

# ------------------------------------------------------------------------------------------
# and then Julia will determine on its own which input argument types make sense and which
# do not:
# ------------------------------------------------------------------------------------------

f(10)

f([1, 2, 3])

# ------------------------------------------------------------------------------------------
# #### Specifying the types of our input arguments
#
# However, we also have the *option* to tell Julia explicitly what types our input arguments
# are allowed to have.
#
# For example, let's write a function `foo` that only takes strings as inputs.
# ------------------------------------------------------------------------------------------

foo(x::String, y::String) = println("My inputs x and y are both strings!")

# ------------------------------------------------------------------------------------------
# We see here that in order to restrict the type of `x` and `y` to `String`s, we just follow
# the input argument name by a double colon and the keyword `String`.
#
# Now we'll see that `foo` works on `String`s and doesn't work on other input argument
# types.
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# To get `foo` to work on integer (`Int`) inputs, let's tack `::Int` onto our input
# arguments when we declare `foo`.
# ------------------------------------------------------------------------------------------

foo(x::Int, y::Int) = println("My inputs x and y are both integers!")

foo(3, 4)

# ------------------------------------------------------------------------------------------
# Now `foo` works on integers! But look, `foo` also still works when `x` and `y` are
# strings!
# ------------------------------------------------------------------------------------------

foo("hello", "hi!")

# ------------------------------------------------------------------------------------------
# This is starting to get to the heart of multiple dispatch. When we declared
#
# ```julia
# foo(x::Int, y::Int) = println("My inputs x and y are both integers!")
# ```
# we didn't overwrite or replace
# ```julia
# foo(y::String, y::String)
# ```
# Instead, we just added an additional ***method*** to the ***generic function*** called
# `foo`.
#
# A ***generic function*** is the abstract concept associated with a particular operation.
#
# For example, the generic function `+` represents the concept of addition.
#
# A ***method*** is a specific implementation of a generic function for *particular argument
# types*.
#
# For example, `+` has methods that accept floating point numbers, integers, matrices, etc.
#
# We can use the `methods` to see how many methods there are for `foo`.
# ------------------------------------------------------------------------------------------

methods(foo)

# ------------------------------------------------------------------------------------------
# Aside: how many methods do you think there are for addition?
# ------------------------------------------------------------------------------------------

methods(+)

# ------------------------------------------------------------------------------------------
# So, we now can call `foo` on integers or strings. When you call `foo` on a particular set
# of arguments, Julia will infer the types of the inputs and dispatch the appropriate
# method. *This* is multiple dispatch.
#
# Multiple dispatch makes our code generic and fast. Our code can be generic and flexible
# because we can write code in terms of abstract operations such as addition and
# multiplication, rather than in terms of specific implementations. At the same time, our
# code runs quickly because Julia is able to call efficient methods for the relevant types.
#
# To see which method is being dispatched when we call a generic function, we can use the
# @which macro:
# ------------------------------------------------------------------------------------------

@which foo(3, 4)

# ------------------------------------------------------------------------------------------
# Let's see what happens when we use `@which` with the addition operator!
# ------------------------------------------------------------------------------------------

@which 3.0 + 3.0

# ------------------------------------------------------------------------------------------
# And we can continue to add other methods to our generic function `foo`. Let's add one that
# takes the ***abstract type*** `Number`, which includes subtypes such as `Int`, `Float64`,
# and other objects you would think of as numbers:
# ------------------------------------------------------------------------------------------

foo(x::Number, y::Number) = println("My inputs x and y are both numbers!")

# ------------------------------------------------------------------------------------------
# This method for `foo` will work on, for example, floating point numbers:
# ------------------------------------------------------------------------------------------

foo(3.0, 4.0)

# ------------------------------------------------------------------------------------------
# We can also add a fallback, duck-typed method for `foo` that takes inputs of any type:
# ------------------------------------------------------------------------------------------

foo(x, y) = println("I accept inputs of any type!")

# ------------------------------------------------------------------------------------------
# Given the methods we've already written for `foo` so far, this method will be called
# whenever we pass non-numbers to `foo`:
# ------------------------------------------------------------------------------------------

v = rand(3)
foo(v, v)

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 9.1
#
# Extend the function `foo`, adding a method that takes only one input argument, which is of
# type `Bool`, and returns "foo with one boolean!"
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### 9.2
#
# Check that the method being dispatched when you execute
# ```julia
# foo(true)
# ```
# is the one you wrote.
# ------------------------------------------------------------------------------------------

@assert foo(true) == "foo with one boolean!"

# ------------------------------------------------------------------------------------------
# Please click on `Validate` on the top, once you are done with the exercises.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using StanSample, MCMCChains, CSV

df = filter(row -> row[:age] >= 18,
  CSV.read(joinpath(@__DIR__, "..", "..", "data", "Howell1.csv"), DataFrame))

stan4_1 = "
// Inferring a Rate
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 20);
  sigma ~ uniform( 0 , 50 );

  // Observed heights
  h ~ normal(mu, sigma);
}
";

m4_1s = SampleModel("m4.1s", stan4_1);

m4_1_data = Dict("N" => length(df[!, :height]), "h" => df[!, :height]);

rc4_1s = stan_sample(m4_1s, data=m4_1_data);

if success(rc4_1s)
  #chn4_1s = read_samples(m4_1s; output_format=:mcmcchains)
  # Update parameter names
  #chn4_1s = replacenames(chn4_1s, Dict("mu" => "μ", "sigma" => "σ"))
  #chn4_1s

  part4_1s = read_samples(m4_1s; output_format=:particles)
  part4_1s |> display
end


text:
# scale the point set down such that the union of all circumcircles lies within the square frame.
# this does not only allow for arbitrary point distributions, but more importantly
# ensures convexity (and correctness) of the resulting Delaunay triangulation.

function scaleShiftPoints( points::Array{T,1} ) where T<:AbstractPoint2D
  # convex hull of point set to find the outer edges
  hull = quickHull( points )
  # union of the outer circumcircles
  ccu = circumcircleUnion( hull, points )
  # frame around the union of circumcircles
  ranges = frameRanges( ccu )
  # scale the point set down to the initial square
  scaledPoints = shrink( points, ranges )
  # do the tessellation on the scaled point set,
  # use the ranges to convert back to original scale afterwards
  return scaledPoints, ranges
end

# convert the edge points back to original scale after the tessellation

function expand( points::Array{Point2D,1}, ranges::NTuple{4,Float64} )
  xmin = ranges[1]
  ymin = ranges[3]
  scale = max( ranges[4] - ranges[3], ranges[2] - ranges[1] ) / 0.98
  offset = 1.01
  scaledPoints = [ Point2D( ( p._x - offset ) * scale + xmin, ( p._y - offset ) * scale + ymin ) for p in points ]
end

function shrink( points::Array{Point2D,1}, ranges::NTuple{4,Float64} )
  h = ranges[4] - ranges[3]
  b = ranges[2] - ranges[1]
  offset = 1.01
  scale = 0.98 / max( h, b )
  scaledPoints = [ Point2D( ( p._x - ranges[1] ) * scale + offset, ( p._y - ranges[3] ) * scale + offset ) for p in points ]
end

#=== auxiliary functions and structs ===#

# convex hull of point set

function quickHull(points::Array{T,1}) where T<:AbstractPoint2D
    return quickHull(Point2D.(points))
end

function quickHull( points::Array{Point2D,1} )
  ConvexHull = Array{Point2D,1}(undef,0)
  A = points[ argmin( getx.( points ) ) ]   # leftmost point
  B = points[ argmax( getx.( points ) ) ]   # rightmost point
  push!( ConvexHull, A )
  push!( ConvexHull, B )
  pr = Array{Point2D,1}(undef,0)            # points to the right of AB
  pl = Array{Point2D,1}(undef,0)            # points to the left of AB
  ( pl, pr ) = dividePointSet( Line2D( A, B ), setdiff( points, ConvexHull ) )
  findHull!( ConvexHull, pr, A, B )         # divide-and-conquer approach
  findHull!( ConvexHull, pl, B, A )
  return ConvexHull
end

function findHull!( ConvexHull::Array{T,1}, points::Array{T,1}, A::T, B::T ) where T<:AbstractPoint2D
  if isempty( points )
    return
  end
  C = findFarthestPoint( A, B, points )
  pos = findfirst( x -> x == A, ConvexHull ) + 1
  insert!( ConvexHull, pos, C )
  pAC = dividePointSet( Line2D(A,C), setdiff(points,[C]) )[2]
  pCB = dividePointSet( Line2D(C,B), setdiff(points,[C]) )[2]
  findHull!( ConvexHull, pAC, A, C )
  findHull!( ConvexHull, pCB, C, B )
end

function dividePointSet( l::Line2D, points::Array{T,1} ) where T <: AbstractPoint2D
  pr = Array{T,1}(undef,0)
  pl = Array{T,1}(undef,0)
  for i in 1:length(points)
    if orientation(l,points[i]) == -1
      push!(pr,points[i])
    else
      push!(pl,points[i])
    end
  end
  return (pl, pr)
end

function findFarthestPoint(a::T,b::T,points::Array{T,1}) where T<:AbstractPoint2D
  distances = Array{Float64,1}(undef,size(points,1))
  for i in 1:length(distances)
    distances[i] = pointLineDistance(a,b,points[i])
  end
  return points[argmax(distances)]
end

function pointLineDistance(a::T,b::T,c::T) where T<:AbstractPoint2D
  return abs( (gety(b)-gety(a))*getx(c) - (getx(b)-getx(a))*gety(c) + getx(b)*gety(a) - gety(b)*getx(a) ) / sqrt( (gety(b)-gety(a))^2 + (getx(b)-getx(a))^2 )
end

mutable struct Circle{T<:AbstractPoint2D}
  c::T
  r::Float64
end

function getcen(circ::Circle{T}) where T<:AbstractPoint2D
    return circ.c
end

function getcenx(circ::Circle{T}) where T<:AbstractPoint2D
    return getx( getcen( circ ) )
end
function getceny(circ::Circle{T}) where T<:AbstractPoint2D
    return gety( getcen( circ ) )
end
function getrad(circ::Circle{T}) where T<: AbstractPoint2D
    return circ.r
end

function circumcircleUnion( hull::Array{T,1}, points::Array{S,1} ) where {T<:AbstractPoint2D,S<:AbstractPoint2D}
  ccU = Array{Circle{Point2D},1}(undef,length(hull))
  hullcirc = copy(hull)
  push!( hullcirc, hull[1] )
  for i in 1:length( hull )
    cc2 = circumcircle( hullcirc[i], hullcirc[i+1] )
    ccU[i] = cc2
    for j in 1:length(points)
      if pointsDistance( points[j], cc2.c ) < cc2.r
        cc3 = circumcircle( hullcirc[i], hullcirc[i+1], points[j] )
        if getrad(cc3) > getrad(ccU[i])
          ccU[i] = cc3
        end
      end
    end
  end
  return ccU
end

function pointsDistance( p1::T, p2::S ) where {T<:AbstractPoint2D, S<:AbstractPoint2D}
  sqrt( ( getx(p1) - getx(p2) ) ^ 2 + ( gety(p1) - gety(p2) ) ^ 2 )
end

function circumcircle( a::T, b::S ) where {T<:AbstractPoint2D, S<:AbstractPoint2D}
  px = (getx(a) + getx(b)) / 2.0
  py = (gety(a) + gety(b)) / 2.0
  dx = (getx(b) - getx(a)) / 2.0
  dy = (gety(b) - gety(a)) / 2.0
  r = sqrt( dx^2 + dy^2 )
  return Circle{Point2D}( Point2D(px,py), r )
end

function circumcircle( a::T, b::S, c::U ) where {T<:AbstractPoint2D,S<:AbstractPoint2D, U<:AbstractPoint2D}
  la = getx(a)^2 + gety(a)^2
  lb = getx(b)^2 + gety(b)^2
  lc = getx(c)^2 + gety(c)^2
  xyy = getx(a) * (gety(b) - gety(c))
  yxx = gety(a) * (getx(b) - getx(c))
  xy = getx(b) * gety(c)
  yx = gety(b) * getx(c)
  z = 2 * ( xyy - yxx + xy - yx )
  px = ( la*(gety(b)-gety(c)) + lb*(gety(c)-gety(a)) + lc*(gety(a)-gety(b)) ) / z
  py = ( la*(getx(c)-getx(b)) + lb*(getx(a)-getx(c)) + lc*(getx(b)-getx(a)) ) / z
  r = sqrt( (px-getx(a))^2 + (py-gety(a))^2 )
  return Circle{Point2D}( Point2D(px,py), r )
end

function frameRanges( ccU::Array{Circle{T},1} ) where T<:AbstractPoint2D
  xmin = minimum( getcenx.(ccU) .- getrad.(ccU) )
  xmax = maximum( getcenx.(ccU) .+ getrad.(ccU) )
  ymin = minimum( getceny.(ccU) .- getrad.(ccU) )
  ymax = maximum( getceny.(ccU) .+ getrad.(ccU) )
  return xmin, xmax, ymin, ymax
end


text:
# # "Point exponential decay" time-stepper convergence test

using PyPlot
using Oceananigans

# Define a few utilities for running tests and unpacking and plotting results
include("ConvergenceTests/ConvergenceTests.jl")

using .ConvergenceTests
using .ConvergenceTests.PointExponentialDecay: run_test

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

""" Run advection-diffusion test for all Nx in resolutions. """
function run_convergence_test(stop_time, proposal_Δt, timestepper)

    # Adjust time-step
    stop_iterations = [round(Int, stop_time / dt) for dt in proposal_Δt]
                 Δt = [stop_time / stop_iter for stop_iter in stop_iterations]

    # Run the tests
    results = [run_test(timestepper=timestepper, Δt=dt, stop_iteration=stop_iter)
               for (dt, stop_iter) in zip(Δt, stop_iterations)]

    return results, Δt
end

unpack_errors(results) = map(r -> r.L₁, results)

stop_time = 3
Δt = 10 .^ range(-3, 0, length=30)  # Equally spaced in log space.
ab2_results, Δt = run_convergence_test(stop_time, Δt, :QuasiAdamsBashforth2)
rk3_results, Δt = run_convergence_test(stop_time, Δt, :RungeKutta3)

ab2_L₁ = unpack_errors(ab2_results)
rk3_L₁ = unpack_errors(rk3_results)

fig, axs = subplots()

loglog(Δt, ab2_L₁, "o--", alpha=0.8, linewidth=1, label="Measured error, Quasi-second-order Adams-Bashforth")
loglog(Δt, rk3_L₁, "^-.", alpha=0.8, linewidth=1, label="Measured error, Third-order Runge-Kutta")

# Guide line to confirm second-order scaling
loglog(Δt, ab2_L₁[1] .* Δt / Δt[1], "k-", alpha=0.4, linewidth=2, label=L"\sim \Delta t")
loglog(Δt, rk3_L₁[1] .* (Δt / Δt[1]).^3, "r-", alpha=0.4, linewidth=2, label=L"\sim \Delta t^3")

xlabel(L"\Delta t")
ylabel("\$L_1\$-norm of \$ | c_\\mathrm{simulated}(t) - \\mathrm{e}^{-t} | \$")
title("Oceananigans time-stepper convergence for \$ c(t) = \\mathrm{e}^{-t} \$")
removespines("top", "right")
legend()

filepath = joinpath(@__DIR__, "figs", "point_exponential_decay_time_stepper_convergence.png")
savefig(filepath, dpi=480)

test_rate_of_convergence(ab2_L₁, Δt, data_points=9, expected=1.0, atol=0.01)
test_rate_of_convergence(rk3_L₁, Δt, data_points=9, expected=3.0, atol=0.01)


text:
@testset "Function interpolation" begin
    rₘₐₓ = 20
    @testset "t₀ = $(t₀), ϕ = $(ϕ)" for (t₀,ϕ) in [(0.0,0.0)# ,
                                                   # (rₘₐₓ/2, π/3)
                                                   ]
        R = FEDVR(range(0,stop=rₘₐₓ,length=11), 10, t₀=t₀, ϕ=ϕ)
        R⁻¹ = pinv(R)
        r = axes(R,1)
        r̃ = range(0,stop=rₘₐₓ,length=1001)
        χ = R[r̃,:]

        fu = r -> r^2*exp(-r)
        fv = r -> r^6*exp(-8r)

        u = R*(R \ fu.(r))
        @test norm(χ * (R⁻¹*u) - fu.(r̃)) < 1e-6
        v = R*(R \ fv.(r))
        @test norm(χ * (R⁻¹*v) - fv.(r̃)) < 5e-5

        @testset "Restricted basis" begin
            R̃ = R[:,2:end-1]
            χ̃ = R̃[r̃,:]

            # # Actually, we want
            # ũ = R̃*(R̃ \ fu.(r))
            # # but R̃*... expands the coefficient vector by two elements
            # # (which are "hidden" by the restriction matrix), which
            # # causes a dimension mismatch below.
            ũ = R̃ \ fu.(r)
            # @test norm(χ̃ * ũ.args[2] - fu.(r̃)) < 5e-6
            @test norm(χ̃ * ũ - fu.(r̃)) < 2e-6
            # ṽ = R̃*(R̃ \ fv.(r))
            ṽ = R̃ \ fv.(r)
            # @test norm(χ̃ * ṽ.args[2] - fv.(r̃)) < 1e-4
            @test norm(χ̃ * ṽ - fv.(r̃)) < 5e-5

            h = r -> (r-rₘₐₓ/2)^2
            c = R̃ \ h.(r)
            # Vandermonde matrix to compare with; its generation is
            # costly, which is why we prefer to compute overlaps with
            # basis functions instead.
            rl = locs(R̃)
            V = R̃[rl,:]
            c̃ = V \ h.(rl)

            @test c ≈ c̃ atol=1e-13
        end
    end
    @testset "Restricted interpolation" begin
        L = 20.0
        n = 10
        t = range(-L/2, stop=L/2, length=n)
        B = FEDVR(t, [4,3,4,2,5,4,2,4,8])
        N = size(B,2)
        xx = axes(B,1)

        Bl = B[:,1:round(Int, 0.6N)+1]
        Bv = B[:,round(Int, 0.4N):end]

        f = sin
        c = B \ f.(xx)
        cl = Bl \ f.(xx)
        cv = Bv \ f.(xx)

        @test length(cl) == size(Bl,2)
        @test c[1:length(cl)-1] == cl[1:end-1]

        @test length(cv) == size(Bv,2)
        @test c[end-length(cv)+1:end] == cv
    end
end


text:
using QuartetNetworkGoodnessFit
using CSV
using PhyloNetworks
using Test
using Distributed

@testset "quartet network GoF tests" begin
  # each file should be its own testset, to run them all even if one has a failure
  include("test_ticr.jl")
  include("test_qnetGoF.jl")
end


text:
legalMoves = [2 3 4
    1 5 7
    1 6 8
    1 9 10
    2 6 11
    3 5 12
    2 9 13
    3 10 14
    4 7 15
    4 8 16
    5 13 17
    6 14 17
    7 11 18
    8 12 19
    9 16 18
    10 15 19
    11 12 20
    13 15 20
    14 16 20
    17 18 19]
# Each line of this matrix legalMoves represents a vertex and specifies three neighboring vertices.

#= Function move expects one number as a current vertex and a matrix of legal moves.
It returns one number as a new, randomly chosen, legal, current vertex. =#
function move(n, A)
    return A[n, rand(1:3)]
end

#= Function simulate expects one number as the number of simulations and a matrix of legal moves.
It returns the average number of moves the spider needed in to get to an ant. =#
function simulate(noSimulations, B)
    sum = 0
    for i in 1:noSimulations
        #= The state of the simulation is reset.
        Ant is positioned in a vertex numbered 1.
        Spider is positioned in a vertex numbered 20.
        Number of moves is reset to 0. =#
        ant = 1
        spider = 20
        noMoves = 0

        # Until the spider moves to the desired position it moves in random, legal moves that are counted.
        while (ant != spider)
            spider = move(spider, B)
            noMoves += 1
        end

        # The number of moves needed in this simulation is added to the sum of all moves needed in previous simulations.
        sum += noMoves
    end

    # The average number of moves is calculated by dividing the sum of all moves needed in all simulations by the number of sumulations.
    return(sum / noSimulations)
end

# Specify the number of simulations for every function simulate call.
noSimulations = 1_000_000

# Run the function 10 times and print their averages.
for j in 1:10
    println(simulate(noSimulations, legalMoves))
end

max_stars_repo_name: zznidar/Random-walk-edges-dodecahedron

max_issues_repo_name: zznidar/Random-walk-edges-dodecahedron

max_forks_repo_name: zznidar/Random-walk-edges-dodecahedron


text:
## ExoplanetsSysSim/src/stellar_table.jl
## (c) 2015 Eric B. Ford

module StellarTable
using ExoplanetsSysSim
using DataFrames
using JLD

if VERSION >= v"0.5-"
  import Compat: UTF8String, ASCIIString
end

export setup_star_table, star_table, num_usable_in_star_table, set_star_table

df = DataFrame()
usable = Array{Int64}(0)
#data = Array{Float64}(0,0)
#colid = Dict()

function setup(sim_param::SimParam; force_reread::Bool = false)
  global df
  if haskey(sim_param,"read_stellar_catalog") && !force_reread
     return df
     #return data
  end
  stellar_catalog_filename = convert(ASCIIString,joinpath(Pkg.dir("ExoplanetsSysSim"), "data", convert(ASCIIString,get(sim_param,"stellar_catalog","q1_q17_dr24_stellar.csv")) ) )
  df = setup(stellar_catalog_filename)
  add_param_fixed(sim_param,"read_stellar_catalog",true)
  return df
end

function setup(filename::ASCIIString; force_reread::Bool = false)
  global df, usable
  if ismatch(r".jld$",filename)
  try
    data = load(filename)
    df = data["stellar_catalog"]
    usable = data["stellar_catalog_usable"]
    Core.typeassert(df,DataFrame)
    Core.typeassert(usable,Array{Int64,1})
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in jld format."))
  end
  else
  try
    df = readtable(filename)
  catch
    error(string("# Failed to read stellar catalog >",filename,"< in ascii format."))
  end

  has_mass = ! (isna(df[:mass]) | isna(df[:mass_err1]) | isna(df[:mass_err2]))
  has_radius = ! (isna(df[:radius]) | isna(df[:radius_err1]) | isna(df[:radius_err2]))
  has_dens = ! (isna(df[:dens]) | isna(df[:dens_err1]) | isna(df[:dens_err2]))
  has_rest = ! (isna(df[:rrmscdpp04p5]) | isna(df[:dataspan]) | isna(df[:dutycycle]))
  is_usable = has_mass & has_radius & has_dens & has_rest
  # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
  # TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
  symbols_to_keep = [ :kepid, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :rrmscdpp04p5, :dataspan, :dutycycle ]
  delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
  usable = find(is_usable)
  df = df[usable, symbols_to_keep]
  end
  return df
  #global data = convert(Array{Float64,2}, df) # df[usable, symbols_to_keep] )
  #global colid = Dict(zip(names(df),[1:length(names(df))]))
  #return data
end

setup_star_table(sim_param::SimParam; force_reread::Bool = false) = setup(sim_param, force_reread=force_reread)
setup_star_table(filename::ASCIIString; force_reread::Bool = false) = setup(filename, force_reread=force_reread)

function num_usable()
  global usable
  @assert typeof(usable) == Array{Int64,1}
  length(usable)
end
num_usable_in_star_table() = num_usable()

function idx(i::Integer)
  global usable
  usable[i]
end

#function col( sym::Symbol )
#  global colid
#  colid[sym]
#end

function star_table(i::Integer, sym::Symbol)
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

#function star_table_data(i::Integer, sym::Symbol)
#  global data
#  return data[i,col(sym)]
#end

function star_table(i::Integer)
  global df, usable
  return df[i,:]
  #return df[usable[i],:]
end

function star_table(i::Integer, sym::Vector{Symbol})
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

function star_table(i::Vector{Integer}, sym::Symbol)
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

function star_table(i::Vector{Integer}, sym::Vector{Symbol})
  global df, usable
  return df[i,sym]
  #return df[usable[i],sym]
end

function set_star_table(df2::DataFrame)
  global df
  df = df2
end

function set_star_table(df2::DataFrame, usable2::Array{Int64,1})
  global df, usable
  df = df2
  usable = usable2
end

end # module StellarTable

using ExoplanetsSysSim.StellarTable


text:
# Discrete distributions
for D in [:Categorical, :Poisson]
    @eval struct $D{P<:$(_params[D].parameters[2])} <: Distribution{Integer,P}
        params :: P
        function $D(; kwargs...)
            params = Tuple(values(kwargs))
            new{typeof(params)}(params)
        end
    end
end

# Continuous distributions
for D in [:Uniform, :Normal, :Exponential, :Gamma]
    @eval struct $D{P<:$(_params[D].parameters[2])} <: Distribution{Real,P}
        params :: P
        function $D(; kwargs...)
            params = Tuple(values(kwargs))
            new{typeof(params)}(params)
        end
    end
end

##########################
# Sampling and densities
##########################

# Categorical
function randCategorical(p::Tuple{Vararg{T}}) where T<:Real
    draw = rand(T)
    cp = zero(T)
    i = zero(Int)
    while cp < draw && i < length(p)
        cp += p[i+=1]
    end
    max(i, 1)
end
randCategorical(p::T...) where T<:Integer = randCategorical(Float64.(p)...)

@cufunc function logpdfCategorical(i::I, p::Tuple{Vararg{Real}}) where I<:Integer
    # p = homogenize(p_) # a bit wasteful, only need to convert p[i]
    # T = eltype(eltype(p))
    T = promote_type(eltype.(p)...)
    ifelse(1 ≤ i ≤ length(p), @inbounds log(T(p[i])), -T(Inf))
end

# Poisson
# TODO eliminar dependencia de StatsFuns:
randPoisson(λ::Real) = convert(Int, poisinvcdf(λ, rand()))

@cufunc function logpdfPoisson(i::Integer, λ::T) where T<:Real
    x = convert(eltype(T), i)
    iszero(λ) && return ifelse(iszero(x), zero(T), -T(Inf))
    x * log(λ) - λ - lgamma(x + one(T))
end

# Uniform
randUniform(a::T, b::T) where T<:Real = a + (b - a)rand(T)
randUniform(a::T, b::T) where T<:Integer = Float64(a) + Float64(b - a)rand()

@cufunc function logpdfUniform(x_::Real, a_::Real, b_::Real)
    x, a, b = promote(x_, a_, b_)
    ifelse(a ≤ x ≤ b, -log(b - a), -typeof(x)(Inf))
end

# Normal
randNormal(μ::T, σ::T) where T<:Real = μ + σ * randn(T)
randNormal(μ::T, σ::T) where T<:Integer = Float64(μ) + Float64(σ) * randn()

@cufunc function logpdfNormal(x_::Real, μ_::Real, σ_::Real)
    x, m, s = promote(x_, μ_, σ_) # cannot use σ because the macro would transform that to CuArrays.cufunc(σ)
    T = eltype(x)
    iszero(s) && return ifelse(x == m, T(Inf), -T(Inf))
    -(((x - m) / s)^2 + T(log2π))/2 - log(s)
end

# Exponential
randExponential(λ::Real) = λ * randexp(eltype(λ))
randExponential(λ::Integer) = λ * randexp()

@cufunc function logpdfExponential(x_::Real, λ_::Real)
    x, λ = promote(x_, λ_)
    ifelse(x < zero(x), -eltype(x)(Inf), log(λ) - λ * x)
end

# Gamma
function _MarsagliaTsang2000(α::T) where T<:Real
    d = α - one(T)/3
    c = one(T) / sqrt(9d)
    while true
        x = randn(T)
        v = (one(T) + c*x)^3
        while v < zero(T)
            x = randn(T)
            v = (one(T) + c*x)^3
        end
        u = rand(T)
        # u < one(F) - F(0.0331)x^4 && return d*v
        log(u) < x^2/2 + d*(one(T) - v + log(v)) && return d*v
    end
end

function randGamma(α::T, θ::T) where T<:Real
    if α < one(T) # use the γ(1+α)*U^(1/α) trick from Marsaglia and Tsang (2000)
        x = θ * _MarsagliaTsang2000(α + one(T)) # Gamma(α + 1, θ)
        e = randexp(T)
        return x * exp(-e / α)
    elseif α == one(T)
        return θ * randexp(eltype(T))
    else
        return θ * _MarsagliaTsang2000(α)
    end
end
randGamma(α::T, θ::T) where T<:Integer = randGamma(Float64(α), Float64(θ))

@cufunc function logpdfGamma(x_::Real, α_::Real, θ_::Real)
    x, α, θ = promote(x_, α_, θ_)
    T = eltype(x)
    -lgamma(α) - α*log(θ) + (α-one(T))log(x) - x/θ
end