text:
module Day20

export get_inputs, get_solution1, get_solution2
export Tile, rotate!, flipy!, layout_image, concatenate_image

using UnPack: @unpack
using OffsetArrays: OffsetArray

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_input("test_input1.txt")
    test_output1 = 20899048083289
    test_output2 = 273
    data = read_input("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_input(io)
    str = read(joinpath(@__DIR__, io), String)
    tile_strs = split(str, r"\n\n|\r\n\r\n", keepempty=false)
    tiles = Tile.(tile_strs)
    return Dict(tile.id => tile for tile in tiles)
end

lines_to_image(lines) = BitArray(permutedims(reduce(hcat, collect.(lines)), [2,1]) .== '#')

## Struct definitions
const Maybe{X} = Union{X, Nothing}
const PlainIndex = Union{Int, <:AbstractVector}

mutable struct Connections{T}
    R::Maybe{T}
    U::Maybe{T}
    L::Maybe{T}
    D::Maybe{T}
end

Base.@inline Base.getproperty(c::Connections, ::Val{B}) where B = getfield(c, B)
Base.@inline Base.setproperty!(c::Connections, ::Val{B}, x) where B = setfield!(c, B, x)

Base.@kwdef mutable struct Tile <: AbstractMatrix{Bool}
    id::Int = 0
    image::Matrix{Bool}
    rotation::Int = 0
    flipy::Bool = false
    parent_connections::Connections = Connections()
    child_connections::Connections = Connections()
end

Tile(image::AbstractMatrix) = Tile(; image)

function Tile(str::AbstractString)
    lines = split(str, r"\n|\r\n", keepempty=false)
    id = parse(Int, only(match(r"Tile.([0-9]+)\:", lines[1]).captures))
    image = lines_to_image(lines[2:end])
    return Tile(; id, image)
end

Connections(; R=nothing, U=nothing, L=nothing, D=nothing) = Connections{Tile}(R, U, L, D)

# Base.show(io::IO, tile::Tile) = show(io, MIME("text/plain"), tile)
function Base.show(io::IO, mime::MIME"text/plain", tile::Tile)
    println("Tile $(tile.id):")
    for line in eachrow(tile)
        for char in line
            print(char ? '#' : '.')
        end
        print('\n')
    end
end

Base.size(tile::Tile) = size(tile.image)

Base.getindex(tile::Tile, ij::CartesianIndex) = getindex(tile, ij.I...)
Base.getindex(tile::Tile, ::Colon, j::PlainIndex) = getindex(tile, axes(tile.image, 1), j)
Base.getindex(tile::Tile, i::PlainIndex, ::Colon) = getindex(tile, i, axes(tile.image, 2))
Base.getindex(tile::Tile, ::Colon, ::Colon) = getindex(tile, axes(tile)...)
function Base.getindex(tile::Tile, i::PlainIndex, j::PlainIndex)
    # rot = tile.rotation
    # ax_i, ax_j = axes(tile.image)
    # if rot == 1
    #     i, j = reverse(ax_j)[j], i
    # elseif rot == 2
    #     i, j = reverse(ax_i)[i], reverse(ax_j)[j]
    # elseif rot == 3
    #     i, j = j, reverse(ax_i)[i]
    # end

    # flipy = tile.flipy
    # if flipy
    #     i = reverse(ax_i)[i]
    # end

    return tile.image[i, j]
end

# Base.@inline Base.getproperty(tile::Tile, ::Val{B}) where B = getfield(tile, B)
# Base.@inline Base.setproperty!(tile::Tile, ::Val{B}, x) where B = setfield!(tile, B, x)

get_border(tile, dir) = get_border(tile, Val(dir))
get_border(tile, ::Val{:R}) = @view tile[:, end]
get_border(tile, ::Val{:U}) = @view tile[begin, :]
get_border(tile, ::Val{:L}) = @view tile[:, begin]
get_border(tile, ::Val{:D}) = @view tile[end, :]

const BORDERS = Val.((:R, :U, :L, :D))

opposite(::Val{:R}) = Val(:L)
opposite(::Val{:U}) = Val(:D)
opposite(::Val{:L}) = Val(:R)
opposite(::Val{:D}) = Val(:U)

# Apply clockwise rotation to a tile and update all tiles attached to it
rotate!(::Nothing, rot=1) = nothing
function rotate!(tile::Tile, rot=1)
    children = tile.child_connections
    @unpack R, U, L, D = children

    for r in 1:rot
        tile_image = copy(tile.image)
        for i in 1:size(tile.image, 1), (j, j_inv) in enumerate(size(tile.image, 2):-1:1)
            tile.image[i, j] = tile_image[j_inv, i]
        end
    end

    tile.rotation = mod(tile.rotation + rot, 4)
    children.R = rotate!(U, rot)
    children.U = rotate!(L, rot)
    children.L = rotate!(D, rot)
    children.D = rotate!(R, rot)

    return tile
end

flipy!(::Nothing) = nothing
function flipy!(tile::Tile)
    children = tile.child_connections
    @unpack R, U, L, D = children

    tile_image = copy(tile.image)
    for (i, i_inv) in enumerate(size(tile.image, 1):-1:1), j in 1:size(tile.image, 2)
        tile.image[i, j] = tile_image[i_inv, j]
    end

    tile.flipy = !tile.flipy
    children.R = flipy!(L)
    children.L = flipy!(R)

    return tile
end

Base.@inline function connect!(tile, other_tile, B)
    setproperty!(tile.child_connections, B, other_tile)
    setproperty!(other_tile.parent_connections, opposite(B), tile)
    return tile
end

## Solution functions
function set_matches!(tile_set, unused_tiles=copy(tile_set))
    connected_tiles = Dict(pop!(unused_tiles))
    while !isempty(unused_tiles)
        for (this_id, this_tile) in connected_tiles
            for b in BORDERS
                getproperty(this_tile.child_connections, b) !== nothing && continue
                this_border = get_border(this_tile, b)
                opp_b = opposite(b)
                for (id, tile) in unused_tiles
                    connected = false
                    for rot in 1:4
                        rotate!(tile)
                        for flip in 1:2
                            flipy!(tile)
                            getproperty(tile.child_connections, opp_b) !== nothing && continue
                            if get_border(tile, opp_b) == this_border
                                # println("Connecting $id to $this_id")
                                connect!(this_tile, tile, b)
                                delete!(unused_tiles, id)
                                push!(connected_tiles, id=>tile)
                                connected = true
                                break
                            end
                        end
                        connected && break
                    end
                    connected && break
                end
            end
        end
    end
    return connected_tiles
end

const x̂ = (0, 1)
const ŷ = (1, 0)

_layout_image!(image, ::Nothing, idx) = nothing

function _layout_image!(image, tile, idx=(0, 0))
    tile in values(image) && return image
    push!(image, idx=>tile)

    @unpack L, R, U, D = tile.child_connections
    _layout_image!(image, R, idx .+ x̂)
    _layout_image!(image, L, idx .- x̂)
    _layout_image!(image, U, idx .- ŷ)
    _layout_image!(image, D, idx .+ ŷ)

    return image
end

layout_image(tiles) = _layout_image!(Dict{Tuple{Int, Int}, Tile}(), first(values(tiles)))

function concatenate_image(layout; keep_borders=false)
    sz_i, sz_j = size(last(first(layout)).image) .- 2
    inds = keys(layout)
    min_i, max_i = extrema(first.(inds))
    min_j, max_j = extrema(last.(inds))
    ax = axes(layout[1,1].image)
    ax = keep_borders ? ax : (ax[1][2:end-1], ax[2][2:end-1])

    block_image = [layout[i, j][ax...] for i in min_i:max_i, j in min_j:max_j]
    cols = reduce.(vcat, eachcol(block_image))
    image = reduce(hcat, cols)

    # image = falses(sz_i*(max_i-min_i+1), sz_j*(max_j-min_j+1))
    # for (i, j) in inds
    #     i_range = (i-min_i)*sz_i .+ (1:sz_i)
    #     j_range = (j-min_j)*sz_j .+ (1:sz_j)
    #     image[i_range, j_range] .= layout[i, j][2:end-1, 2:end-1]
    # end

    return Tile(; image)
end

# Part 1
function get_solution1(data)
    tiles = set_matches!(deepcopy(data))
    image = layout_image(tiles)
    inds = keys(image)
    return prod(image[i, j].id for i in extrema(first.(inds)), j in extrema(last.(inds)))
end

# Part 2
function sea_monster_image()
    str = """
                      #
    #    ##    ##    ###
     #  #  #  #  #  #
    """
    image = lines_to_image(split(str, r"\n|\r\n", keepempty=false)) |> collect
end

function get_solution2(data)
    tiles = set_matches!(deepcopy(data))
    image = concatenate_image(layout_image(tiles))
    sea_monster = sea_monster_image()
    sz_sm_i, sz_sm_j = size(sea_monster)
    sz_im_i, sz_im_j = size(image)
    count_sea_monster = count(sea_monster)
    last_i = sz_im_i - sz_sm_i + 1
    last_j = sz_im_j - sz_sm_j + 1
    num_waters = Int[]
    this_image = falses(size(image)...)
    num_sea_monsters = 0
    for rot in 1:4
        rotate!(image)
        for flip in 1:2
            flipy!(image)
            this_image .= image
            for j in 1:last_j, i in 1:last_i
                i_range = i:i+sz_sm_i-1
                j_range = j:j+sz_sm_j-1
                this_sub_image = @view this_image[i_range, j_range]
                this_num_waters = count(this_sub_image .& sea_monster)
                # println("($i_range, $j_range)")
                if (i, j) == (17, 2)
                    # println("rot=$(mod(rot, 4)), flip=$(mod(flip, 2))")
                    # show(stdout, MIME("text/plain"), Tile(; id=0, image=this_image[i_range, :]))
                end
                if this_num_waters == count_sea_monster
                    # println("($i_range, $j_range)")
                    # show(stdout, MIME("text/plain"), image)
                    num_sea_monsters += 1
                end
            end
        end
    end
    return count(image) - num_sea_monsters*count_sea_monster
end

# 303 too low

end


text:
using NCDatasets
using Statistics
using Interpolations
using LinearAlgebra
using Glob
using JLD

"""
generate_cm_params(cm_params::Array{Float64},
                            cm_param_names::Array{String})

Generate a CLIMAParameters file, setting the values (cp_params) of
a group of CLIMAParameters (cm_param_names).
"""
function generate_cm_params(cm_params::Array{Float64}, cm_param_names::Array{String})
    # Generate version
    version = rand(11111:99999)
    if length(cm_params) == length(cm_param_names)
        open("clima_param_defs_$(version).jl", "w") do io
            for i in 1:length(cm_params)
                write(io, "CLIMAParameters.Atmos.SubgridScale.
                    $(cm_param_names[i])(::EarthParameterSet) =
                    $(cm_params[i])\n")
            end
        end
    else
        throw(ArgumentError("Number of parameter names must be equal to number of values provided."))
    end
    return version
end

"""
get_clima_profile(sim_dir::String,
                     var_name::Array{String,1};
                     ti::Float64=0.0,
                     tf::Float64=0.0,
                     z_obs::Union{Array{Float64, 1}, Nothing} = nothing,
                     get_variance::Bool=false)

Get a time-averaged (from ti to tf) vertical profile from a ClimateMachine.jl diagnostics
file, interpolating to observation heights (z_obs) if necessary. If get_variance, return
the variance over the same averaging period.
"""
function get_clima_profile(sim_dir::String,
                     var_name::Array{String,1};
                     ti::Float64=0.0,
                     tf::Float64=0.0,
                     z_obs::Union{Array{Float64, 1}, Nothing} = nothing,
                     get_variance::Bool=false)

    ds_name = glob("$(sim_dir)/*.nc")[1]
    ds = NCDataset(ds_name)

    if length(var_name) == 1 && occursin("z", var_name[1])
        prof_vec = Array(ds[var_name[1]])
    else
        # Time in Datetime format
        t = Array(ds["time"])
        # Convert to float (seconds)
        t = [(time_ .- t[1]).value/1000.0 for time_ in t]
        dt = t[2]-t[1]
        ti_diff, ti_index = findmin( broadcast(abs, t.-ti) )
        tf_diff, tf_index = findmin( broadcast(abs, t.-tf) )

        prof_vec = zeros(0)
        ts_vec = zeros(0, length(ti_index:tf_index))
        # If simulation does not contain values for ti or tf, return high value
        if ti_diff > dt || tf_diff > dt
            for i in 1:length(var_name)
                var_ = Array(ds[var_name[i]])
                append!(prof_vec, 1.0e4*ones(length(var_[:, 1])))
            end
            println("Initial or final times are not contained in netcdf file.")
        else
            for i in 1:length(var_name)
                var_ = Array(ds[var_name[i]])
                if !isnothing(z_obs)
                    z_ = Array(ds["z"])
                    # Form interpolation
                    var_itp = interpolate( (z_, t), var_,
                        ( Gridded(Linear()), Gridded(Linear()) ))
                    # Evaluate interpolation
                    var_ = var_itp(z_obs, t)
                end
                # Average in time
                append!(prof_vec, mean(var_[:, ti_index:tf_index], dims=2))
                ts_vec = cat(ts_vec, var_[:, ti_index:tf_index], dims=1)
            end
        end
    end
    close(ds)
    if get_variance
        cov_mat = cov(ts_vec, dims=2)
        return prof_vec, cov_mat
    else
        return prof_vec
    end
end

"""
agg_clima_ekp(n_params::Integer, output_name::String="ekp_clima")

Aggregate all iterations of the parameter ensembles and write to file, given the
number of parameters in each parameter vector (p).
"""
function agg_clima_ekp(n_params::Integer, output_name::String="ekp_clima")
    # Get versions
    version_files = glob("versions_*.txt")
    # Recover parameters of last iteration
    last_up_versions = readlines(version_files[end])

    ens_all = Array{Float64, 2}[]
    for (it_num, file) in enumerate(version_files)
        versions = readlines(file)
        u = zeros(length(versions), n_params)
        for (ens_index, version_) in enumerate(versions)
            if it_num == length(version_files)
                open("../../../ClimateMachine.jl/test/Atmos/EDMF/$(version_)", "r") do io
                    u[ens_index, :] = [parse(Float64, line) for (index, line) in enumerate(eachline(io)) if index%3 == 0]
                end
            else
                open("$(version_).output/$(version_)", "r") do io
                    u[ens_index, :] = [parse(Float64, line) for (index, line) in enumerate(eachline(io)) if index%3 == 0]
                end
            end
        end
        push!(ens_all, u)
    end
    save( string(output_name,".jld"), "ekp_u", ens_all)
    return
end

ext: jl


text:
using CUDAnative

function kernel_vadd(a, b, c)
    i = blockIdx().x
    return i
end

Base.code_llvm(stdout, kernel_vadd, Tuple{CuDeviceArray{Float32,1,AS.Global},
                                          CuDeviceArray{Float32,1,AS.Global},
                                          CuDeviceArray{Float32,1,AS.Global}})

using TinyCassette
struct Ctx end

Base.code_llvm(stdout, TinyCassette.execute,
               Tuple{Ctx, typeof(kernel_vadd),
                     CuDeviceArray{Float32,1,AS.Global},
                     CuDeviceArray{Float32,1,AS.Global},
                     CuDeviceArray{Float32,1,AS.Global}})


text:
__precompile__()

module GeoInterface

    using RecipesBase

    export  AbstractPosition, Position,
            AbstractGeometry, AbstractGeometryCollection, GeometryCollection,
            AbstractPoint, Point,
            AbstractMultiPoint, MultiPoint,
            AbstractLineString, LineString,
            AbstractMultiLineString, MultiLineString,
            AbstractPolygon, Polygon,
            AbstractMultiPolygon, MultiPolygon,
            AbstractFeature, Feature,
            AbstractFeatureCollection, FeatureCollection,

            geotype, # methods
            xcoord, ycoord, zcoord, hasz,
            coordinates,
            geometries,
            geometry, bbox, crs, properties,
            features

    abstract type AbstractPosition{T <: Real} <: AbstractVector{T} end
    geotype(::AbstractPosition) = :Position
    xcoord(::AbstractPosition) = error("xcoord(::AbstractPosition) not defined.")
    ycoord(::AbstractPosition) = error("ycoord(::AbstractPosition) not defined.")
    # optional
    zcoord(::AbstractPosition) = error("zcoord(::AbstractPosition) not defined.")
    hasz(::AbstractPosition) = false
    coordinates(p::AbstractPosition) = hasz(p) ? Float64[xcoord(p),ycoord(p),zcoord(p)] : Float64[xcoord(p),ycoord(p)]
    # (Array-like indexing # http://julia.readthedocs.org/en/latest/manual/arrays/#arrays)
    Base.eltype(p::AbstractPosition{T}) where {T <: Real} = T
    Base.ndims(AbstractPosition) = 1
    Base.length(p::AbstractPosition) = hasz(p) ? 3 : 2
    Base.size(p::AbstractPosition) = (length(p),)
    Base.size(p::AbstractPosition, n::Int) = (n == 1) ? length(p) : 1
    Base.getindex(p::AbstractPosition, i::Int) = (i==1) ? xcoord(p) : (i==2) ? ycoord(p) : (i==3) ? zcoord(p) : nothing
    Base.convert(::Type{Vector{Float64}}, p::AbstractPosition) = coordinates(p)
    # Base.linearindexing{T <: AbstractPosition}(::Type{T}) = LinearFast()

    abstract type AbstractGeometry end
    coordinates(obj::AbstractGeometry) = error("coordinates(::AbstractGeometry) not defined.")

        abstract type AbstractPoint <: AbstractGeometry end
        geotype(::AbstractPoint) = :Point

        abstract type AbstractMultiPoint <: AbstractGeometry end
        geotype(::AbstractMultiPoint) = :MultiPoint

        abstract type AbstractLineString <: AbstractGeometry end
        geotype(::AbstractLineString) = :LineString

        abstract type AbstractMultiLineString <: AbstractGeometry end
        geotype(::AbstractMultiLineString) = :MultiLineString

        abstract type AbstractPolygon <: AbstractGeometry end
        geotype(::AbstractPolygon) = :Polygon

        abstract type AbstractMultiPolygon <: AbstractGeometry end
        geotype(::AbstractMultiPolygon) = :MultiPolygon

        abstract type AbstractGeometryCollection <: AbstractGeometry end
        geotype(::AbstractGeometryCollection) = :GeometryCollection
        geometries(obj::AbstractGeometryCollection) = error("geometries(::AbstractGeometryCollection) not defined.")

    abstract type AbstractFeature end
    geotype(::AbstractFeature) = :Feature
    geometry(obj::AbstractFeature) = error("geometry(::AbstractFeature) not defined.")
    # optional
    properties(obj::AbstractFeature) = Dict{String,Any}()
    bbox(obj::AbstractFeature) = nothing
    crs(obj::AbstractFeature) = nothing

    abstract type AbstractFeatureCollection end
    geotype(::AbstractFeatureCollection) = :FeatureCollection
    features(obj::AbstractFeatureCollection) = error("features(::AbstractFeatureCollection) not defined.")
    # optional
    bbox(obj::AbstractFeatureCollection) = nothing
    crs(obj::AbstractFeatureCollection) = nothing

    include("operations.jl")
    include("geotypes.jl")
    include("plotrecipes.jl")
end


text:
include("5.3.1 (a+b arctan(c x^n))^p.jl")
include("5.3.2 (d x)^m (a+b arctan(c x^n))^p.jl")
include("5.3.3 (d+e x)^m (a+b arctan(c x^n))^p.jl")
include("5.3.4 u (a+b arctan(c x))^p.jl")
include("5.3.5 u (a+b arctan(c+d x))^p.jl")
include("5.3.6 Exponentials of inverse tangent.jl")
include("5.3.7 Miscellaneous inverse tangent.jl")

max_stars_repo_path: src/rules/5 Inverse trig functions/5.3 Inverse tangent/.jl


text:
"""
    matex(input::Vector)
Takes a one-dimensional vector as input and formats it for LaTeX display.
"""
function matex(input::Vector)
	S = ""
	for i in input[1:end-1]
		S = string(S, " ", i, " \\\\")
	end
	S = string(
		"\\begin{bmatrix}",
		S, " ", input[end],
		"\\end{bmatrix}")
	S
end

"""
    matex(input::Matrix)
Takes a two-dimensional matrix as input and formats it for LaTeX display.
"""
function matex(input::Matrix)
	S = ""
	rows = collect(eachrow(input))
	for row in rows
		for i in row[1:end-1]
			S = string(S, i, " & ")
		end
		S = string(S, last(row))
		if row != last(rows)
			S = string(S, " \\\\ ")
		end
	end
	S = string("\\begin{bmatrix}", S, "\\end{bmatrix}")
	S
end


text:
# This file includes testset for Plugin

@testset "PluginTestSet" begin
    @info "Running PluginTestSet ..."

    # Construction of a new plugin
    Base.@kwdef struct MeanPlugin{PR}  <: AbstractPlugin
        process::PR = x -> mean(x)
    end

    # Try equip a writer in a model.
    model = Model()
    clock=Clock(0., 0.01, 10.)
    addnode!(model, SinewaveGenerator(), label=:gen)
    addnode!(model, Writer(buflen=50, plugin=MeanPlugin()), label=:writer)
    addbranch!(model, :gen => :writer)
    simulate!(model, clock)

    # Test the simulation data
    data = read(getnode(model, :writer).component, flatten=false)
    @test length(data) == 20
    for (t,x) in data
        @test isapprox(x, mean(sin.(2 * pi * t)))
    end

    @info "Done PluginTestSet ..."
end # testset


text:
if isfile("version_git.jl")
    include("version_git.jl")
end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
"Evaluate formula as fully as possible."
function partialeval(domain::Domain, state::GenericState, term::Term)
    funcs = merge(GLOBAL_FUNCTIONS, state.values, get_funcdefs(domain))
    return eval_term(term, Subst(), funcs)
end

"Get domain constant type declarations as a list of clauses."
function get_const_clauses(domain::Domain)
   return [@julog($ty(:o) <<= true) for (o, ty) in get_constypes(domain)]
end

"Get domain type hierarchy as a list of clauses."
function get_type_clauses(domain::Domain)
    clauses = [[Clause(@julog($ty(X)), Term[@julog($s(X))]) for s in subtys]
               for (ty, subtys) in get_typetree(domain) if length(subtys) > 0]
    return length(clauses) > 0 ? reduce(vcat, clauses) : Clause[]
end

"Get all proof-relevant Horn clauses for PDDL domain."
function get_clauses(domain::Domain)
   return [collect(values(get_axioms(domain)));
           get_const_clauses(domain); get_type_clauses(domain)]
end

"""
    find_matches(term, state, domain=nothing)

Returns a list of all matching substitutions of `term` with respect to
a given `state` and `domain`.
"""
function find_matches(domain::GenericDomain, state::GenericState, term::Term)
    if term.name in keys(state.values)
        clauses = Vector{Clause}(get_fluents(state))
        _, subst = resolve(term, clauses; mode=:all)
    else
        clauses = isnothing(domain) ? Clause[] : get_clauses(domain)
        clauses = Clause[clauses; collect(state.types); collect(state.facts)]
        funcs = state.values
        _, subst = resolve(term, clauses; funcs=funcs, mode=:all)
    end
    matches = Term[substitute(term, s) for s in subst]
    return matches
end


text:
module DIDebug

using Cxx; using JITTools; using ELF; using ObjFileBase; using DataStructures; using DWARF

import JITTools: buffer, datasize, datapointer
import ObjFileBase: readmeta, getSectionLoadAddress

export julia_passes

include(Pkg.dir("Cxx","test","llvmincludes.jl"))
if isdefined(Base, :active_repl)
   include(Pkg.dir("Cxx","src","CxxREPL","replpane.jl"))
end

addHeaderDir(joinpath(JULIA_HOME,"..","..","deps","llvm-svn","lib"))
cxx"""
#include "llvm/AsmParser/Parser.h"
#include "llvm/ExecutionEngine/JITEventListener.h"
#include "llvm/ExecutionEngine/ObjectMemoryBuffer.h"
#include "llvm/IR/LegacyPassManager.h"
#include "llvm/IR/LegacyPassManager.h"
#include "llvm/IR/DiagnosticPrinter.h"
#include "llvm/CodeGen/Passes.h"
#include "ExecutionEngine/MCJIT/MCJIT.h"
#include "llvm/Analysis/TargetTransformInfo.h"
#include "llvm/Analysis/Passes.h"
#include "llvm/Transforms/Scalar.h"
#include "llvm/Transforms/Vectorize.h"
#include "llvm/CodeGen/MachineFunctionAnalysis.h"
#include "llvm/CodeGen/MachineModuleInfo.h"
#include "llvm/Target/TargetLoweringObjectFile.h"
#include "llvm/CodeGen/MachineFunction.h"
#include "llvm/IR/DIBuilder.h"
#include "llvm/Support/TargetSelect.h"
#include "llvm/Support/TargetRegistry.h"
#include "llvm/MC/SubtargetFeature.h"
#include "llvm/MC/MCDisassembler.h"
#include "llvm/MC/MCRegisterInfo.h"
#include "llvm/MC/MCSubtargetInfo.h"
#include "llvm/MC/MCAsmInfo.h"
#include "llvm/MC/MCInstPrinter.h"
#include "llvm/MC/MCInstrInfo.h"
#include "llvm/MC/MCInst.h"
using namespace llvm;
"""

cxxparse(readall(Pkg.dir("DIDebug","src","FunctionMover.cpp")))

immutable LoadedObjectInfo
    LOI::pcpp"llvm::RuntimeDyld::LoadedObjectInfo"
end

immutable ObjectFile
    file::pcpp"llvm::object::ObjectFile"
end
Base.convert(::Type{pcpp"llvm::object::ObjectFile"}, x::ObjectFile) = x.file
Cxx.cppconvert(x::ObjectFile) = x.file

datasize(obj::ObjectFile) = icxx"$obj->getData().size();"
datapointer(obj::ObjectFile) = icxx"$obj->getData().data();"
readmeta(obj::ObjectFile) = readmeta(buffer(obj))

function parseIR(str)
    icxx"""
        llvm::SMDiagnostic Err;
        auto mod = llvm::parseAssembly(
            llvm::MemoryBufferRef(
            llvm::StringRef($(pointer(str)),$(sizeof(str))),
            llvm::StringRef("<in-memory>")),
            Err, jl_LLVMContext);
        if (mod == nullptr)
        {
            std::string message = "Failed to parse LLVM Assembly: \n";
            llvm::raw_string_ostream stream(message);
            Err.print("julia",stream,true);
            jl_error(stream.str().c_str());
        }
        return mod.release();
    """
end

function parseBitcode(str)
    icxx"""
        std::string Message;
        raw_string_ostream Stream(Message);
        DiagnosticPrinterRawOStream DP(Stream);
        ErrorOr<std::unique_ptr<Module>> ModuleOrErr = llvm::parseBitcodeFile(
            llvm::MemoryBufferRef(
            llvm::StringRef($(pointer(str)),$(sizeof(str))),
            llvm::StringRef("<in-memory>")),
            jl_LLVMContext, [&](const DiagnosticInfo &DI) { DI.print(DP); });

        if (ModuleOrErr.getError()) {
            Stream.flush();
            jl_error(Message.c_str());
        }

        return ModuleOrErr.get().release();
    """
end

function getSectionLoadAddress(L::LoadedObjectInfo, Name::ByteString)
    Addr = icxx"$(L.LOI)->getSectionLoadAddress($(pointer(Name)));"
    Addr == C_NULL && error("Section Not loaded or not found")
    Addr
end

getSectionLoadAddress(L::LoadedObjectInfo, sec) = getSectionLoadAddress(L,sectionname(sec))

function findObjForAddr(AddrMap, objects, addr)
    searchsortedlast(AddrMap, addr)
end

cxx"""
static llvm::object::OwningBinary<llvm::object::ObjectFile>
emitObject(llvm::TargetMachine *TM, llvm::Module *M) {
  llvm::legacy::PassManager PM;

  M->setDataLayout(TM->createDataLayout());

  // The RuntimeDyld will take ownership of this shortly
  llvm::SmallVector<char, 4096> ObjBufferSV;
  llvm::raw_svector_ostream ObjStream(ObjBufferSV);
  llvm::MCContext *Ctx;

  // Turn the machine code intermediate representation into bytes in memory
  // that may be executed.
  if (TM->addPassesToEmitMC(PM, Ctx, ObjStream, true))
    llvm::report_fatal_error("Target does not support MC emission!");

  // Initialize passes.
  PM.run(*M);

  std::unique_ptr<llvm::MemoryBuffer> CompiledObjBuffer(
                                new llvm::ObjectMemoryBuffer(std::move(ObjBufferSV)));

  llvm::ErrorOr<std::unique_ptr<llvm::object::ObjectFile>> LoadedObject =
    llvm::object::ObjectFile::createObjectFile(CompiledObjBuffer->getMemBufferRef());
  llvm::object::OwningBinary<llvm::object::ObjectFile> Obj(std::move(LoadedObject.get()),std::move(CompiledObjBuffer));

  return Obj;
}

static std::unique_ptr<llvm::RuntimeDyld::LoadedObjectInfo>
loadObject(llvm::ExecutionEngine *EE, llvm::object::OwningBinary<llvm::object::ObjectFile> Obj)
{
    return std::unique_ptr<llvm::RuntimeDyld::LoadedObjectInfo>{nullptr}; // ((llvm::MCJIT*)EE)->addObjectFile(std::move(Obj));
}

static std::unique_ptr<llvm::RuntimeDyld::LoadedObjectInfo>
loadModule(llvm::TargetMachine *TM, llvm::ExecutionEngine *EE, llvm::Module *M)
{
    return loadObject(EE,emitObject(TM,M));
}

void RunFPM(llvm::legacy::FunctionPassManager *FPM, llvm::Module *M)
{
    llvm::Module::FunctionListType &flist = M->getFunctionList();
    auto it = flist.begin();
    for(;it != flist.end(); ++it)
        FPM->run(*it);
}

/// addPassesToX helper drives creation and initialization of TargetPassConfig.
legacy::PassManager *createPassManager(LLVMTargetMachine *TM,
                                          MachineFunctionAnalysis *MFA,
                                          AnalysisID StopAfter) {
  legacy::PassManager *PM = new legacy::PassManager;

  // Add internal analysis passes from the target machine.
  PM->add(createTargetTransformInfoWrapperPass(TM->getTargetIRAnalysis()));

  // Targets may override createPassConfig to provide a target-specific
  // subclass.
  TargetPassConfig *PassConfig = TM->createPassConfig(*PM);
  PassConfig->setStartStopPasses(nullptr, nullptr, StopAfter);

  // Set PassConfig options provided by TargetMachine.
  PassConfig->setDisableVerify(false);

  PM->add(PassConfig);

  PassConfig->addIRPasses();

  PassConfig->addCodeGenPrepare();

  PassConfig->addPassesToHandleExceptions();

  PassConfig->addISelPrepare();

  // Install a MachineModuleInfo class, which is an immutable pass that holds
  // all the per-module stuff we're generating, including MCContext.
  MachineModuleInfo *MMI = new MachineModuleInfo(
      *TM->getMCAsmInfo(), *TM->getMCRegisterInfo(), TM->getObjFileLowering());
  PM->add(MMI);

  // Set up a MachineFunction for the rest of CodeGen to work on.
  PM->add(MFA);

  //TM->setFastISel(true);

  // Ask the target for an isel.
  if (PassConfig->addInstSelector())
    return nullptr;

  PassConfig->addMachinePasses();

  PassConfig->setInitialized();

  return PM;
}
"""

function setup_FPM(mod,passes)
    # Initialize the optimization passes
    FPM = @cxxnew llvm::legacy::FunctionPassManager(mod)
    icxx"$FPM->add(createTargetTransformInfoWrapperPass(jl_TargetMachine->getTargetIRAnalysis()));";
    for p in passes
        add(FPM, @eval @cxx $p())
    end
    @cxx FPM->doInitialization();
    FPM
end

# This is due to a bug in staged functions
# force add to specialize
add{T}(FPM,x::T) = @cxx FPM->add(x)

julia_passes = [
    :createTypeBasedAliasAnalysisPass,
    :createCFGSimplificationPass,
    :createPromoteMemoryToRegisterPass,
    :createInstructionCombiningPass,
    :createScalarReplAggregatesPass,
    :createInstructionCombiningPass,
    :createJumpThreadingPass,
    :createInstructionCombiningPass,
    :createReassociatePass,
    :createEarlyCSEPass,
    :createLoopIdiomPass,
    :createLoopRotatePass,
#    :createLowerSimdLoopPass,
    :createLICMPass,
    :createLoopUnswitchPass,
    :createInstructionCombiningPass,
    :createIndVarSimplifyPass,
    :createLoopDeletionPass,
    :createSimpleLoopUnrollPass,
    :createLoopVectorizePass,
    :createInstructionCombiningPass,
    :createGVNPass,
    :createSCCPPass,
    :createSinkingPass,
    :createInstructionSimplifierPass,
    :createInstructionCombiningPass,
    :createJumpThreadingPass,
    :createDeadStoreEliminationPass,
    :createAggressiveDCEPass,
];

counter = 0
function setup_test(f, passes = julia_passes, dependencies = true)
    mod = @cxxnew llvm::Module(pointer("test"),icxx"jl_LLVMContext;")
    CU = icxx"""
        Module *m = $mod;
        m->addModuleFlag(llvm::Module::Warning, "Dwarf Version",4);
        m->addModuleFlag(llvm::Module::Error, "Debug Info Version",
          llvm::DEBUG_METADATA_VERSION);
        DIBuilder DIB(*m);
        DICompileUnit *CU = DIB.createCompileUnit(0x01, "test.jl", ".", "julia", true, "", 0);
        DIB.finalize();
        CU;
    """
    FPM = setup_FPM(mod,passes)
    mover = @cxxnew FunctionMover2(mod, CU, dependencies)
    r = @cxx MapFunction(f,mover)
    global counter += 1
    name = string("foobar",counter)
    @cxx r->setName(pointer(name))
    @cxx pcpp"llvm::Function"(r.ptr)->setLinkage(@cxx llvm::GlobalValue::ExternalLinkage)
    #@cxx r->dump()
    @cxx RunFPM(FPM,mod)
    #icxx"jl_ExecutionEngine->addModule(std::unique_ptr<Module>{$mod});"
    name, mod
end

MF_passes = [
    (@cxx &llvm::ExpandISelPseudosID)
    (@cxx &llvm::PrologEpilogCodeInserterID)
    (@cxx &llvm::LocalStackSlotAllocationID)
]

function run_pm(f, passes = julia_passes, dependencies = true, StopAfter = (@cxx &llvm::VirtRegRewriterID))
    name, mod = setup_test(f,passes,dependencies)
    tm = @cxx (@cxx jl_ExecutionEngine)->getTargetMachine()
    ltm = pcpp"llvm::LLVMTargetMachine"(tm.ptr)
    mfa = @cxxnew MachineFunctionAnalysis(*(ltm))
    PM = @cxx createPassManager(ltm,mfa,StopAfter)
    @cxx PM->run(*(mod))
    mod, mfa
end

function FindFunctionNamed(str)
    icxx"jl_ExecutionEngine->FindFunctionNamed($(pointer(str)));"
end

function WriteBitcodeToFile(fname,mod)
    icxx"""
        std::error_code err;
        StringRef fname_ref = StringRef($(pointer(fname)));
        raw_fd_ostream OS(fname_ref, err, sys::fs::F_None);
        WriteBitcodeToFile($mod, OS);
    """
end

function getLam(llvmf)
  lam = icxx"""cast<ConstantInt>(cast<ConstantAsMetadata>(
    $llvmf->getMetadata("lam")->getOperand(0).get())->getValue())->getZExtValue();"""
  unsafe_pointer_to_objref(convert(Ptr{Void},lam))
end

include("disassembler.jl")
include("investigate.jl")

end # module


text:
import Literate
using Base.Filesystem: cp, basename

println("Start generating tests for Literate.jl examples")

EXAMPLEDIR = joinpath(@__DIR__, "..", "examples")
GENERATEDDIR = joinpath(@__DIR__, "examples")
mkpath(GENERATEDDIR)
for name in readdir(EXAMPLEDIR)
    path = abspath(joinpath(EXAMPLEDIR, name))
    if endswith(path, ".jl")
        example = path
        script = Literate.script(example, GENERATEDDIR)
    elseif isdir(path)
        folder = path
        target = joinpath(GENERATEDDIR, basename(folder))
        cp(folder, target, force=true)
        @info "cp $folder -> $target"
    end
end

println("Finished generating tests for Literate.jl examples")


text:
"""
    CoxIngersollRoss()

[Cox-Ingersoll-Ross model](https://en.wikipedia.org/wiki/Cox%E2%80%93Ingersoll%E2%80%93Ross_model).
"""
struct CoxIngersollRoss <: Model end

"""
    evaluate(IRD, CoxIngersollRoss(), n)

Evaluate interest rate derivative `IRD` using `CoxIngersollRoss` model.

# Arguments
- `IRD::InterestRateDerivative`: interest rate derivative
- `n`: number of paths to simulate
"""
function evaluate(IRD::InterestRateDerivative, m::CoxIngersollRoss, n::Int64 = 12)
    Δt = IRD.t / n
    rates = [IRD.r]
    for i = 0:n
        dr = IRD.k * (IRD.θ - rates[end]) * Δt + IRD.σ * √rates[end] * randn()
        append!(rates, rates[end] + dr)
    end
    return rates
end


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

using Random
export TIHCSearch

struct TIHCSearch <: LocalSearchAlgorithm
    ssize::Int32 # sample size
    restarts::Int32
    use_local_improvement::Bool
    TIHCSearch() = new(1, 1, false)
    TIHCSearch(ssize, restarts; use_local_improvement=false) = new(Int32(ssize), Int32(restarts), use_local_improvement)
 end

function randomsearch(index::SearchGraph{T}, dist::Function, q::T, res::KnnResult, best::KnnResult, navigation_state, ssize::Integer) where T
    n = length(index.db)
    ssize = min(n, max(maxlength(res), ssize))
    range = 1:n
    @inbounds for i in 1:ssize
        # while abs(prev_dist - last(res).dist) < rs.tol
        objID = rand(range)
        S = get(navigation_state, objID, UNKNOWN)
        if S == UNKNOWN
            navigation_state[objID] = VISITED
            d = convert(Float32, dist(q, index.db[objID]))
            push!(res, objID, d)
            push!(best, objID, d)
        end
    end

    res
end

function search(isearch::TIHCSearch, index::SearchGraph{T}, dist::Function, q::T, res::KnnResult, navigation_state, hints=EMPTY_INT_VECTOR) where T
    n = length(index.db)
    restarts = min(isearch.restarts, n)
    best = KnnResult(1)
    # if length(hints) == 0 hints = rand(1:n, bs.bsize) end

    @inbounds for i in 1:restarts
        empty!(best)
        randomsearch(index, dist, q, res, best, navigation_state, isearch.ssize)
        length(best) == 0 && break  # all items in the neighborhood were visited
        # @show i, restarts, isearch, best, length(best)
        start_point = first(best).objID
        navigation_state[start_point] = EXPLORED
        d = convert(Float32, dist(q, index.db[start_point]))
        push!(res, start_point, d)
        hill_climbing(index, dist, q, res, navigation_state, start_point, isearch.use_local_improvement)
    end

    res
end

function opt_expand_neighborhood(fun, algo::TIHCSearch, n::Int, iter::Int)
    f_(w) = ceil(Int, w * (rand() - 0.5))
    f(x, w) = max(1, x + f_(w))
    logn = ceil(Int, log(2, n+1))
    w = 1

    while w <= logn  ## log log n
        #TIHCSearch(min(200, n), 1, use_local_improvement=true) |> fun
        TIHCSearch(algo.ssize, f(algo.restarts, w)) |> fun
        TIHCSearch(f(algo.ssize, w), algo.restarts) |> fun
        w += w
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayPrivateEndpointConnectionProperties <: SwaggerModel
    privateEndpoint::Any # spec type: Union{ Nothing, PrivateEndpoint } # spec name: privateEndpoint
    privateLinkServiceConnectionState::Any # spec type: Union{ Nothing, PrivateLinkServiceConnectionState } # spec name: privateLinkServiceConnectionState
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    linkIdentifier::Any # spec type: Union{ Nothing, String } # spec name: linkIdentifier

    function ApplicationGatewayPrivateEndpointConnectionProperties(;privateEndpoint=nothing, privateLinkServiceConnectionState=nothing, provisioningState=nothing, linkIdentifier=nothing)
        o = new()
        validate_property(ApplicationGatewayPrivateEndpointConnectionProperties, Symbol("privateEndpoint"), privateEndpoint)
        setfield!(o, Symbol("privateEndpoint"), privateEndpoint)
        validate_property(ApplicationGatewayPrivateEndpointConnectionProperties, Symbol("privateLinkServiceConnectionState"), privateLinkServiceConnectionState)
        setfield!(o, Symbol("privateLinkServiceConnectionState"), privateLinkServiceConnectionState)
        validate_property(ApplicationGatewayPrivateEndpointConnectionProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ApplicationGatewayPrivateEndpointConnectionProperties, Symbol("linkIdentifier"), linkIdentifier)
        setfield!(o, Symbol("linkIdentifier"), linkIdentifier)
        o
    end
end # type ApplicationGatewayPrivateEndpointConnectionProperties

const _property_map_ApplicationGatewayPrivateEndpointConnectionProperties = Dict{Symbol,Symbol}(Symbol("privateEndpoint")=>Symbol("privateEndpoint"), Symbol("privateLinkServiceConnectionState")=>Symbol("privateLinkServiceConnectionState"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("linkIdentifier")=>Symbol("linkIdentifier"))
const _property_types_ApplicationGatewayPrivateEndpointConnectionProperties = Dict{Symbol,String}(Symbol("privateEndpoint")=>"PrivateEndpoint", Symbol("privateLinkServiceConnectionState")=>"PrivateLinkServiceConnectionState", Symbol("provisioningState")=>"ProvisioningState", Symbol("linkIdentifier")=>"String")
Base.propertynames(::Type{ ApplicationGatewayPrivateEndpointConnectionProperties }) = collect(keys(_property_map_ApplicationGatewayPrivateEndpointConnectionProperties))
Swagger.property_type(::Type{ ApplicationGatewayPrivateEndpointConnectionProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPrivateEndpointConnectionProperties[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPrivateEndpointConnectionProperties }, property_name::Symbol) =  _property_map_ApplicationGatewayPrivateEndpointConnectionProperties[property_name]

function check_required(o::ApplicationGatewayPrivateEndpointConnectionProperties)
    true
end

function validate_property(::Type{ ApplicationGatewayPrivateEndpointConnectionProperties }, name::Symbol, val)
end


text:
ENV["PYTHON"]=""
import Pkg
Pkg.build("PyCall")


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1968a1b0-9090-11eb-3424-ad2150889fd1
using PlutoUI

# ╔═╡ c09c0e20-d8be-11ea-11be-0966527e6f2d
using Plots

# ╔═╡ b711bf40-c912-11eb-120d-1bb53ca25d17
using ControlSystems

# ╔═╡ cfbeb720-90a8-11eb-1ac2-3bda60b2285c
using OrdinaryDiffEq

# ╔═╡ 889652b2-dbe8-11ea-0835-e966a39af471
md"""# Problema de estabilidad absoluta de Lur`e
## Partimos de una planta con función de transferencia
"""

# ╔═╡ bf0a5c82-d8bc-11ea-3b1c-b397e3b526fa
H1(s)=(s+3)/(s^2+7s+10)

# ╔═╡ 6b0b4680-9095-11eb-3afe-9b792355080b
md"lo primero que miramos es qué polos tiene, factorizando el denominador:

$H_1=\frac{s+3}{(s+5)(s+2)}$

No hay polos ni ceros en el semiplano derecho ni tampoco en el origen.
"

# ╔═╡ 9999f2c0-dbdd-11ea-04a5-dd56cf1eae5c
md"""# Conjeturas A-K
## Se trata de comprobar el rango para el cuál la realimentación negativa con ganancia K es estable:
+ Se calcula el número de polos $n$ en el semiplano derecho (inestables)
+ El criterio de Nyquist dice que la planta es estable si la respuesta en frecuecia rodea al punto $-1/K+0*i$ $n$ veces en sentido antihorario.
"""

# ╔═╡ 4628f590-dbe3-11ea-2a7a-83cabda49ee7
function conjeturas(H,K)
	#respuesta en frecuencia
	wp=10.0.^(-50:0.01:50) #rango de 10^-50 a 10ymin=^50
	wn=-10.0.^(50:-0.01:-50)
	w=vcat(wn,0,wp)#todas las frecuencias incluidas las negativas
	Hs=H.(w*im)

	#si hay polos en el eje imaginario habrá valores infinitos o NaN para pintar bien esos casos los "recorto", en este caso no se cumplen las condiciones de las conjeturas pero puede igualemnte ser interesante ver el rango de estabilidad.

	recortar(x)=min(max(x,-10^3),10^3)

	#pinta la flecha en la frecuencia w0
if isfinite(maximum(norm.(Hs)))
	#pinto normalmente
	plot(real.(Hs),imag.(Hs), label=false)
	#pinta la flecha
	quiver!([real(H(0))],[imag(H(0))],
			quiver=([real(H(0.01*im)-H(0))],
				[imag(H(0.01*im)-H(0))]),label=false)
else
	#en este caso es cuando hay que recortar
	plot(recortar.(real.(Hs)),recortar.(imag.(Hs)), label=false)
	ylims!((-10,10))
	xlims!((-10,10))
end

#para que el movimiento sea intuitivo el slider ha de ir desde la derecha

#pinta el punto
	scatter!([-1/K],[0],label=false)
	title!("""La respuesta en frecuencia ha de rodear al punto
	en sentido antihorario tantas veces como polos inestables
	hay, K=$(K)""")
end

# ╔═╡ 490c8852-909f-11eb-2522-ad23744ce9d8
	#CSS que hace los Sliders más grandes
    html"""<style type="text/css">
    input[type=range] {
    width: 75%;
    } </style>
    """

# ╔═╡ 203ea980-9090-11eb-25f8-8ff76ab5262a
begin
	K_max_conjeturas=10
	@bind sliderAK Slider((-1:0.001:1)*K_max_conjeturas,default=1.0)
end

# ╔═╡ 93bfc0d0-cb03-11eb-283b-6d01e65bed3c
conjeturas(H1,sliderAK)

# ╔═╡ 29e11c20-90ae-11eb-1348-0f5a981d3361
md"Esto sugiere que es estable en el sector $[-\frac{10}{3}, \infty)$
Podemos ver analíticamente el -10/3 observando que el corte es el $w_c=0$ y por lo tanto $H(w_c)=\frac{0+3}{(0+5)(0+2)}=\frac{3}{10}
 \to H(0)k=-1 \to k=-10/3$

Veamos qué dicen los criterios *de verdad* y no las conjeturas
"

# ╔═╡ 66f0b770-d8c5-11ea-0ee3-7ba9b6dfb585
md"""# Criterio del círculo
## Antes de empezar comprueba que:
+ Como es una funcion de transferencia ya es contolable/observable (OK)
+ No hay integradores ni realimentación directa $$d=0$$
¿Es estable? Usa los casos A y B. Si no lo es usa los casos C y D.

**Cada caso** de los que se puede usar da **un rango** posible de ganancias
"""

# ╔═╡ 00b900b0-d8c6-11ea-0e42-0f4bfaf17fd2
md"## Planta estable"

# ╔═╡ ed754f20-db56-11ea-02fa-3581c244dcbd
md"""### Caso A
$$0=k_1<k_2$$
La recta ha de quedar a la izquierda
"""

# ╔═╡ 9975e6e0-db57-11ea-30e1-3ffbd104839c
K_maxA=100

# ╔═╡ adee0ba0-d988-11ea-19c7-e324f1f86031
@bind k2a Slider((0:0.001:1)*K_maxA,default=1.0)

# ╔═╡ 839b4c72-9092-11eb-2237-1dc284e3041e
k2a

# ╔═╡ 2b1196f0-90dc-11eb-2e20-89a61fcb1498
md"Este caso ha sido más conservador que las conjeturas $[0, \infty)$"

# ╔═╡ b47374ee-db55-11ea-2b8f-5b38e3a2c5a0
md"""### Caso B
$$k_1<0<k_2$$
El círculo ha de rodear la respuesta en frecuencia
"""

# ╔═╡ bd8ef6f0-dbb8-11ea-3fe0-3b17ccdb1fcd
begin
	K_maxB1=10
	K_maxB2=100
end

# ╔═╡ fe58cf80-dbb8-11ea-1c6c-93fef94ed79c
@bind k1b Slider((-1:0.001:0)*K_maxB1,default=-1.0)

# ╔═╡ d2de7080-dbb8-11ea-1d7e-11b40cd2ce8d
@bind k2b Slider((0:0.001:1)*K_maxB2,default=1.0)

# ╔═╡ d3077b80-9092-11eb-0f8b-85c6edde81f5
(k1=k1b,k2=k2b)

# ╔═╡ 65071560-90dc-11eb-24f3-3d03b8461507
md"Este caso es igual que las conjeturas AK es decir $[-\frac{10}{3}, \infty)$ "

# ╔═╡ 1ae51110-d98c-11ea-3b02-1ddccacadb56
md"## Planta inestable
En este caso **la planta de antes no nos sirve** ya que era estable, usamos estra otra

$H_2=\frac{s+3}{(s+5)(s-2)}$

Que tiene **un** polo inestable

"

# ╔═╡ 8054e0c0-90a2-11eb-20fb-ad26bded4db5
H2(s)=(s+3)/(s^2+3s-10)

# ╔═╡ 9ca66610-d985-11ea-0426-13ef5f785091
function circulo(H,K1,K2,caso)
	"""hay cuatro casos se pueden determinar por los signos de k1 y k2
	------ Estables
	A) 0=k1<k2 recta a la izquierda
	B) k1<0<k2 dentro del circulo
	------ Inestables
	C) 0<k1<=k2 La respueta rodea al circulo en sentido antihorario n veces donde n es el número de polos inestables
	D) k1<k2<0 La respuesta de -H rodea en el sentido antihorario n veces el circulo D(-k1,-k2) es como C pero con todo negado"""

	#compruebo que las entradas son correctas
	if !isfinite(H(0))
		return "No puede haber polos en el origen"
	end
	if caso=='A'
		if (K1==0)&(0<K2)
			titulo="Recta a la izquierda"
		else
			return "Error: en el caso A 0=k1<k2"
		end
	elseif caso=='B'
		if (K1<0)&(0<K2)
			titulo="El circulo rodea la respuesta en frecuencia"
		else
			return "Error: en el caso B k1<0<k2"
		end
	elseif caso=='C'
		if (0<K1)&(K1<=K2)
			titulo="La respueta rodea al circulo en sentido antihorario n veces\n donde n es el número de polos inestables"
		else
			return "Error: en el caso C 0<k1<=k2 "
		end
	elseif caso=='D'
		if (K1<=K2)&(K2<0)
			H2(x)=-H(x)
            return circulo( H2,-K2,-K1,'C')
		else
			return "Error: en el caso D k1<k2<0"
		end
	else
		return "No es un caso valido, ha de ser 'A','B','C' o 'D'"
	end
	#por fin los cálculos...
	#respuesta en frecuencia
	    wp=10.0.^(-50:0.01:50) #rango de 10^-50 a 10ymin=^50
	    wn=-10.0.^(50:-0.01:-50)
		w=vcat(wn,wp)#todas las frecuencias incluidas las negativas

	    Hs=H.(w*im)
	    plot(Hs,line=:arrow,label=false)
	    Hs=collect(Hs)
	    plot(real.(Hs),imag.(Hs), line=:arrow,label=false)
	    #pinta la flecha
	    quiver!([real(H(0))],[imag(H(0))],
		        quiver=([real(H(0.01*im)-H(0))],[imag(H(0.01*im)-H(0))])
		        ,label=false)

	#pinta el "circulo"
	if caso=='A'
		#el círculo es una recta 1/K1=inf
		ymin=minimum(imag.(Hs))
	    ymax=maximum(imag.(Hs))
	    y=[ymin,ymax]
		x=[-1/K2,-1/K2]
	    plot!(x,y,label=false)
	    title!(titulo)
	else
		#el círuclo es un círculo de vedad
		p1=-1/K1
		p2=-1/K2 #puntos de corte
	    centro=(p1+p2)/2
	    radio=abs(p2-p1)/2
	    theta=collect(0:0.01:2*pi)
	    circunferencia=centro.+radio*exp.(im*theta)
		scatter!([-1/K1 -1/K2],[0 0],label=false)
	    plot!(circunferencia,label="circulo")
		title!(titulo)
	end
end

# ╔═╡ b642eb40-d988-11ea-26a4-374a74d5dc08
circulo(H1,0,k2a,'A')

# ╔═╡ 9d833cf0-db03-11ea-3b28-33343306f4dc
circulo(H1,k1b,k2b,'B')

# ╔═╡ 7697fac0-90a4-11eb-1687-5ff2a5cf415d
md"Primero veamos que dicen las conjeturas"

# ╔═╡ ab676b50-90a4-11eb-1632-f1f60fd49f07
begin
	K_max_conjeturas2=10
	@bind sliderAK2 Slider((-1:0.001:1)*K_max_conjeturas2,default=1.0)
end

# ╔═╡ 81144e40-90a4-11eb-2941-f53c46a9300c
conjeturas(H2,sliderAK2)

# ╔═╡ c2d18d10-90dc-11eb-1792-3779513d19cb
md"Las conjeturas dicen $[\frac{10}{3}, \infty)$"

# ╔═╡ f9938e60-d98b-11ea-1e6d-2fc1399f7bef
md"""### Caso C
$$0<k_1 \le k_2$$

Como es inestable habrá $n$ polos en el semiplano derecho.
El círculo ha de **estar rodeado** por la respuesta en frecuencia $n$ veces en sentido antihorario.
"""

# ╔═╡ 1be6ec20-dbba-11ea-2479-ed82409ebdd7
begin
	K_maxC1=10
	K_maxC2=100
end

# ╔═╡ 76a9f030-dbba-11ea-3326-95356d82a2c5
@bind k1c Slider((0:0.001:1)*K_maxC1,default=1.0)

# ╔═╡ 783969d0-dbba-11ea-10d0-b1cea6c3b727
@bind k2c Slider((0:0.001:1)*K_maxC2,default=2.0)

# ╔═╡ 5d374e70-9093-11eb-3126-bfa8080d6bc2
(k1=k1c,k2=k2c)

# ╔═╡ e3f469a0-dbb9-11ea-30bb-813e5688c466
circulo(H2,k1c,k2c,'C')

# ╔═╡ e748e0ce-90dc-11eb-14cf-5f0c6c4ae075
md"Lo mismo que las conjeturas   $[\frac{10}{3}, \infty)$"

# ╔═╡ c45e9020-d989-11ea-028f-5f36703d75da
md"""### Caso D
$$k_1<k_2<0$$

Es como el caso C pero con todos los signos cambiados y H cambiado por -H:
El círculo ha de **estar rodeado** por la respuesta en frecuencia (de -H) $n$ veces en sentido antihorario.
"""

# ╔═╡ 245ff590-dbba-11ea-1deb-1976fc8d165d
begin
	K_maxD1=100
	K_maxD2=100
end

# ╔═╡ 8d0b4eee-dbba-11ea-2fad-5d7ed6f8035b
@bind k1d Slider((-1:0.001:0)*K_maxD1,default=-2.0 )

# ╔═╡ 8f55db30-dbba-11ea-0370-ddadc3b11dfe
@bind k2d Slider((-1:0.001:0)*K_maxD1,default=-1.0)

# ╔═╡ 3f2dfc70-9094-11eb-340c-592edc764cd1
(k1d=k1d,k2d=k2d)

# ╔═╡ 11a30cf0-9094-11eb-2af0-6f424889df66
circulo(H2,k1d,k2d,'D')

# ╔═╡ 105992d0-90dd-11eb-1a3d-653111528c53
md"No es posible meter el círculo dentro de la respuesta en frecuencia, este caso **no nos dice nada**"

# ╔═╡ 8e2736c0-d8bb-11ea-2fde-43a428dd9281
md"""# Criterio de popov
## Antes de empezar comprueba que:
+ Como es una funcion de transferencia ya es contolable/observable
+ Hay a lo sumo un integrador y si lo hay la planta se puede poner como $$d/s+ H_{estable}$$ con d>0 es decir
    + La parte del ingregrador es positiva ($$H(s) \to \infty$$  si  $$s\to 0^+$$)
    + Quitando el integrador la planta es estable"""

# ╔═╡ 7c27fd1e-c90f-11eb-0edf-29b2a8b2867b
md"Como ejemplo podemos usar la planta H1 añadiendo un polo en el origen:

$H_3=\frac{s+3}{s(s+5)(s+2)}$"

# ╔═╡ b7044bc0-9880-11eb-31b6-b1578086897a
H3(s)=H1(s)/s

# ╔═╡ bc4b38a0-c913-11eb-2eb9-25b5e6468a61
md"""### Primero a ver que nos dicen las cojneturas AK
Siendo rigurosos no nos dicen **nada** ¿Por qué?

El problema es que no **son aplicables** si hay un polo con parte real igual a cero como es el caso...

Pero al fin y al cabo son conjeturas y lo que queremos es una idea *optimista* de dónde buscar y las conjeturas nos la dan ya que si una función cualquiera del sector es estable ha de serlo también el caso particular $f(x)=kx$ que es lo que calcularemos aquí.
"""

# ╔═╡ 31d44760-c90f-11eb-3a7a-ffb4d6c20b53
begin
	K_max_conjeturas3=10
	@bind sliderAK3 Slider((-1:0.001:1)*K_max_conjeturas3,default=1.0)
end

# ╔═╡ 439f2e10-c90f-11eb-1b76-eb23238ba418
begin
	conjeturas(H3,sliderAK3)
	ylims!(-1,1) #hay que hacer zoom para ver algo
	xlims!(()) #esto directamente quita el límite de x y lo pone automático
end

# ╔═╡ 710513f0-c910-11eb-2abb-3df94293cdb4
md"**OJO** Ahora no está claro cuando la respuesta *rodea* al punto ya que es abierta..

Lo que si es importante es ver cuando **cambiamos de lado** (la estabilidad cambia al pasar de un lado de la curva al otro lado).

En este caso se cambia de lado cuando $k=0$ y nos aproximamos a la curva cuando $k=\pm \infty$

Luego basta con tomar un ejemplo concreto de ganancia ($k=1$ por ejemplo) y ver si es estable para ver en qué lado es estable y en cual inestable)

Realimentando con $k=1$ tenemos

$H_{LazoCerrado}=\frac{H_3}{1+H_3}=\frac{\frac{s+3}{s(s+5)(s+2)}}{1+\frac{s+3}{s(s+5)(s+2)}}=\frac{s+3}{s(s+5)(s+2)+(s+3)}$

$H_{LazoCerrado}=\frac{s+3}{s^3+7s^2+11s+3}$

Que es estable por Routh–Hurwitz ya que

$7>0$
$3>0$
$7 \cdot 11>3$

Luego el lado que contine $k=1$ es estable y el otro inestable por tanto el intervalo estable es $[0, \infty)$.

Si no queremos hacerlo a meno podemos usar la toolbox de control:"

# ╔═╡ d8394bbe-c912-11eb-250f-1b7385c0a536
S=tf([1.0, 0.0],[1.0])

# ╔═╡ 76c049a0-90ac-11eb-1892-a31cb390383d
md"Evaluando H(s) con S nos da la función de transferencia (esto funciona sólo por que H es racional)"

# ╔═╡ 33a205b0-c913-11eb-22a9-d51437c6252c
H3s=H3(S)

# ╔═╡ 420138b0-c913-11eb-2509-558e30b8048b
#el lazo cerrado con ganancia 1
H_LC=minreal(H3s/(1+H3s)) #minreal es para cancelar los ceros y polos sobrantes

# ╔═╡ 5c92fab0-c913-11eb-0542-2bd8ac7560de
#y por fin los ceros y los polos
zpk(H_LC)

# ╔═╡ 2cfce6fe-c915-11eb-1f34-5be3aed8045a
md"## Ahora a ver qué dice Popov"

# ╔═╡ 497f9250-d8bc-11ea-194d-77979b94814d
function popov(H,K,alpha)
	if K<0
		return "Error, en Popov K>=0"
	end
	W(w)=real(H(w*im))+w*im*imag(H(w*im)) #funcion asociada
	wp=10.0.^(-50:0.01:50) #rango de 10^-50 a 10ymin=^50
	wn=-10.0.^(-50:0.01:50)
	w=hcat(wn,wp) #todas las frecuencias incluidas las negativas
	Ws=W.(w)
	p=plot(Ws, line = :arrow,label=false)
	#recta
	xmin=minimum(real.(Ws))
	xmin=min(xmin,-1/K) #me aseguro de pintar el corte en el eje
	xmax=maximum(real.(Ws))
	x=[xmin,xmax]
	y=(x.+1/K)./alpha
	scatter!([-1/K],[0],label=false)
	plot!(x.+y.*im,label=false)
	title!("La recta ha de estar por encima de la función asociada")
end

# ╔═╡ fee20250-d8bb-11ea-299c-6d093b0edd20
K_max=1000

# ╔═╡ cb075200-d8bb-11ea-234b-6d4b6fca9083
@bind K Slider(0.0:0.01:K_max,default=10.0)

# ╔═╡ f2567110-dbe0-11ea-24a0-4d1a72151516
@bind pendiente Slider(0.0:0.1:90.0,default=45.0)

# ╔═╡ f8ad6aa0-d8bb-11ea-0dea-85d5051066fe
md"K=$(K)"

# ╔═╡ e443a3e0-d8c0-11ea-362f-75314f8b2cbb
alpha=1/tan(pendiente*pi/180)

# ╔═╡ dbcceee0-d985-11ea-2a3e-f9c3ec08cb75
popov(H3,K,alpha)

# ╔═╡ 367c9ac0-90dd-11eb-011b-0f7c4fca2c89
md"Es decir  $[0, \infty)$"

# ╔═╡ 3b303230-90a5-11eb-148b-d9c93a8c3772
md"# Simulando el sistema
Hasta ahora hemos usado funciones de tranferencia, pero ¿cómo lo simulamos?
Lo más fácil es pasar al espacio de estados

$\dot x =Ax+Bu$

$y=Cx$

Asumiremos siempre en este tema que $D=0$
"

# ╔═╡ e7264d30-90a6-11eb-07fc-d1f2120385cd
H=H1(S)

# ╔═╡ 98a5d8b2-90dd-11eb-3061-1713c5934c8c
md"Obtenido esto podemos pasar al espacio de estados"

# ╔═╡ f2cdb330-90a6-11eb-3169-8fb43ab0c79c
estados=ss(H)

# ╔═╡ a4378a10-90ac-11eb-0041-a9145a07928c
md"Definimos ahora el sistema"

# ╔═╡ 17f33990-90ad-11eb-2fa5-b762b4f1138b
md"# Sector
Una función f(x) pertenece al sector $[k_1 k_2]$ si está entre las rectas $y=k_1x$ y $y=k_2x$, algebráicamente:

$k_1x^2 \le xf(x) \le k_2x^2$

u obviamdo la división por cero (suponiendo que $\frac{f(x)}{x}$ sea continua en $0$ definiendola como el $lim_{x\to 0}\frac{f(x)}{x}$)

$k_1\le \frac{f(x)}{x} \le k_2 \  \forall x\ne0$
"

# ╔═╡ b1a75e90-90ad-11eb-0f7d-1d2b98b9fce4
f(x)=x*abs(x)-10/3*sin(x)

# ╔═╡ a9dd2000-90df-11eb-0edb-e364fefafe32
@bind k1 Slider(-10:0.01:10,default=-1.0)

# ╔═╡ 3387bd60-90e0-11eb-1963-7dfc892de25c
@bind k2 Slider(0:0.01:100,default=1.0)

# ╔═╡ a134fad0-90e0-11eb-1145-df3b4d6da246
@bind maximo Slider(1:0.01:100,default=1.0)

# ╔═╡ 041d31e0-90b3-11eb-1d07-c3b6bb15e4ae
begin
	x=-maximo:0.01:maximo
	fx=f.(x)
	f_max=maximum(fx)
	f_min=minimum(fx)
	fig1=plot(x,fx)
	plot!(x,k1*x,label="$(k1)*x",color=:red)
	plot!(x,k2*x,label="$(k2)*x",color=:red)
	sector1 = Shape([(0.0, 0.0), (maximo, k1*maximo),
			   (maximo, k2*maximo), (0.0, 0.0)])
	sector2 = Shape([(0.0, 0.0), (-maximo, -k1*maximo),
			   (-maximo, -k2*maximo), (0.0, 0.0)])
    plot!(sector1, fillcolor = plot_color(:yellow, 0.3),label=false)
	plot!(sector2, fillcolor = plot_color(:yellow, 0.3),label=false)
	plot!(ylims = (f_min,f_max))
	title!("La función ha de estar dentro de la zona amarilla")

	fig2=plot(x,fx./x,label="f(x)/x")
	sector3 = Shape([(-maximo, k1), (maximo, k1),
			   (maximo, k2), (-maximo, k2),(-maximo, k1)])
	plot!(sector3, fillcolor = plot_color(:yellow, 0.3),label=false)

	l = @layout [a; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ 0eae14a0-c923-11eb-115a-616a0e34d031
max(1,NaN)

# ╔═╡ d6fd4a4e-90e0-11eb-34bf-ab3bb29cc905
md"Jugando un poco es fácil ver que es  $[-\frac{10}{3}, \infty)$
Para verlo analíticamente habría que hacer

$\frac{f(x)}{x}=|x|-\frac{10}{3}\frac{sin(x)}{x}$

Y ver que toma cualquier valor entre $-\frac{10}{3}$ (cuando $x \to 0$) e $\infty$ (cuando $x \to \infty$)
"

# ╔═╡ af179aae-90a7-11eb-2a7b-1b96a28c8ebc
function control(y)
	r=0
    u=r-f(y)
end

# ╔═╡ 21c1cfa0-90a7-11eb-3c54-23fa9b677080
function derivadas(x,p,t)
	A,B,C=p #extraemos los parámetros
	y=C*x
	y=y[1] #no quiero un vector sino un escalar [2.5] no es lo mismo que 2.5
	u=control(y) #calculamos el control
	dx=A*x+B*u
end

# ╔═╡ 1fc9b3b0-90a8-11eb-3b34-77d759a763d5
let
	x0 = zeros(size(estados.B))
	x0[1]=1
	parametros=(estados.A, estados.B, estados.C)
	tspan = (0.0, 30.0)
	prob = ODEProblem(derivadas,x0,tspan,parametros);
	sol = solve(prob,Tsit5());

	# Salida y referencia
	#la salida y el control no se guardan así que los reconstruyo
	ts=tspan[1]:0.01:tspan[end]
	X=sol.(ts)
	y=[(estados.C*x)[1] for x in X] #explicar paso a paso
	u=control.(y)

	plot(ts,y,xaxis="t",yaxis="y(t)",label="y(t)")
	#referencia
	fig1=plot!(ts,0.0*ts, linestyle=:dot, linewidth=2, label="referencia")

	#Señal de control, u no se guarda, hay que recalcularlo
	fig2=plot(ts,u,xaxis="t",yaxis="u(t)",legend=false)
	l = @layout [a ; b]
	plot(fig1,fig2, layout=l)
end

# ╔═╡ Cell order:
# ╟─889652b2-dbe8-11ea-0835-e966a39af471
# ╠═bf0a5c82-d8bc-11ea-3b1c-b397e3b526fa
# ╟─6b0b4680-9095-11eb-3afe-9b792355080b
# ╟─9999f2c0-dbdd-11ea-04a5-dd56cf1eae5c
# ╠═4628f590-dbe3-11ea-2a7a-83cabda49ee7
# ╠═1968a1b0-9090-11eb-3424-ad2150889fd1
# ╠═490c8852-909f-11eb-2522-ad23744ce9d8
# ╠═203ea980-9090-11eb-25f8-8ff76ab5262a
# ╠═c09c0e20-d8be-11ea-11be-0966527e6f2d
# ╠═93bfc0d0-cb03-11eb-283b-6d01e65bed3c
# ╟─29e11c20-90ae-11eb-1348-0f5a981d3361
# ╟─66f0b770-d8c5-11ea-0ee3-7ba9b6dfb585
# ╠═9ca66610-d985-11ea-0426-13ef5f785091
# ╟─00b900b0-d8c6-11ea-0e42-0f4bfaf17fd2
# ╟─ed754f20-db56-11ea-02fa-3581c244dcbd
# ╠═9975e6e0-db57-11ea-30e1-3ffbd104839c
# ╠═adee0ba0-d988-11ea-19c7-e324f1f86031
# ╠═839b4c72-9092-11eb-2237-1dc284e3041e
# ╠═b642eb40-d988-11ea-26a4-374a74d5dc08
# ╟─2b1196f0-90dc-11eb-2e20-89a61fcb1498
# ╟─b47374ee-db55-11ea-2b8f-5b38e3a2c5a0
# ╠═bd8ef6f0-dbb8-11ea-3fe0-3b17ccdb1fcd
# ╠═fe58cf80-dbb8-11ea-1c6c-93fef94ed79c
# ╠═d2de7080-dbb8-11ea-1d7e-11b40cd2ce8d
# ╟─d3077b80-9092-11eb-0f8b-85c6edde81f5
# ╠═9d833cf0-db03-11ea-3b28-33343306f4dc
# ╟─65071560-90dc-11eb-24f3-3d03b8461507
# ╟─1ae51110-d98c-11ea-3b02-1ddccacadb56
# ╠═8054e0c0-90a2-11eb-20fb-ad26bded4db5
# ╟─7697fac0-90a4-11eb-1687-5ff2a5cf415d
# ╟─ab676b50-90a4-11eb-1632-f1f60fd49f07
# ╠═81144e40-90a4-11eb-2941-f53c46a9300c
# ╟─c2d18d10-90dc-11eb-1792-3779513d19cb
# ╟─f9938e60-d98b-11ea-1e6d-2fc1399f7bef
# ╠═1be6ec20-dbba-11ea-2479-ed82409ebdd7
# ╠═76a9f030-dbba-11ea-3326-95356d82a2c5
# ╠═783969d0-dbba-11ea-10d0-b1cea6c3b727
# ╟─5d374e70-9093-11eb-3126-bfa8080d6bc2
# ╠═e3f469a0-dbb9-11ea-30bb-813e5688c466
# ╟─e748e0ce-90dc-11eb-14cf-5f0c6c4ae075
# ╟─c45e9020-d989-11ea-028f-5f36703d75da
# ╠═245ff590-dbba-11ea-1deb-1976fc8d165d
# ╠═8d0b4eee-dbba-11ea-2fad-5d7ed6f8035b
# ╠═8f55db30-dbba-11ea-0370-ddadc3b11dfe
# ╟─3f2dfc70-9094-11eb-340c-592edc764cd1
# ╠═11a30cf0-9094-11eb-2af0-6f424889df66
# ╟─105992d0-90dd-11eb-1a3d-653111528c53
# ╟─8e2736c0-d8bb-11ea-2fde-43a428dd9281
# ╟─7c27fd1e-c90f-11eb-0edf-29b2a8b2867b
# ╠═b7044bc0-9880-11eb-31b6-b1578086897a
# ╟─bc4b38a0-c913-11eb-2eb9-25b5e6468a61
# ╟─31d44760-c90f-11eb-3a7a-ffb4d6c20b53
# ╠═439f2e10-c90f-11eb-1b76-eb23238ba418
# ╟─710513f0-c910-11eb-2abb-3df94293cdb4
# ╠═b711bf40-c912-11eb-120d-1bb53ca25d17
# ╠═d8394bbe-c912-11eb-250f-1b7385c0a536
# ╟─76c049a0-90ac-11eb-1892-a31cb390383d
# ╠═33a205b0-c913-11eb-22a9-d51437c6252c
# ╠═420138b0-c913-11eb-2509-558e30b8048b
# ╠═5c92fab0-c913-11eb-0542-2bd8ac7560de
# ╟─2cfce6fe-c915-11eb-1f34-5be3aed8045a
# ╠═497f9250-d8bc-11ea-194d-77979b94814d
# ╠═fee20250-d8bb-11ea-299c-6d093b0edd20
# ╠═cb075200-d8bb-11ea-234b-6d4b6fca9083
# ╠═f2567110-dbe0-11ea-24a0-4d1a72151516
# ╟─f8ad6aa0-d8bb-11ea-0dea-85d5051066fe
# ╟─e443a3e0-d8c0-11ea-362f-75314f8b2cbb
# ╠═dbcceee0-d985-11ea-2a3e-f9c3ec08cb75
# ╟─367c9ac0-90dd-11eb-011b-0f7c4fca2c89
# ╟─3b303230-90a5-11eb-148b-d9c93a8c3772
# ╠═e7264d30-90a6-11eb-07fc-d1f2120385cd
# ╟─98a5d8b2-90dd-11eb-3061-1713c5934c8c
# ╠═f2cdb330-90a6-11eb-3169-8fb43ab0c79c
# ╠═a4378a10-90ac-11eb-0041-a9145a07928c
# ╠═21c1cfa0-90a7-11eb-3c54-23fa9b677080
# ╟─17f33990-90ad-11eb-2fa5-b762b4f1138b
# ╠═b1a75e90-90ad-11eb-0f7d-1d2b98b9fce4
# ╠═a9dd2000-90df-11eb-0edb-e364fefafe32
# ╠═3387bd60-90e0-11eb-1963-7dfc892de25c
# ╠═a134fad0-90e0-11eb-1145-df3b4d6da246
# ╟─041d31e0-90b3-11eb-1d07-c3b6bb15e4ae
# ╠═0eae14a0-c923-11eb-115a-616a0e34d031
# ╟─d6fd4a4e-90e0-11eb-34bf-ab3bb29cc905
# ╠═af179aae-90a7-11eb-2a7b-1b96a28c8ebc
# ╠═cfbeb720-90a8-11eb-1ac2-3bda60b2285c
# ╠═1fc9b3b0-90a8-11eb-3b34-77d759a763d5

max_stars_repo_path: Tema 3.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 3.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 3.jl

max_forks_repo_name: Dictino/CNL


text:
using Luxor, Test, Random

# test empty drawings

Drawing()
[Drawing() for x in 1:10]
[Drawing(x, x) for x in 1:10]
Drawing(NaN, NaN, :rec)
[Drawing(x, x, :rec) for x in 1:10]
[Drawing(x, x, :image) for x in 1:10]

p = Point(1, 2)
x, y = p

@test x == 1.0
@test y == 2.0

bb = BoundingBox(p, 2p)

rp = rand(bb)

@test rp in bb

@testset "center3pts" begin
    @test center3pts(Point(1,1), Point(3,1), Point(2,0)) == (Point(2.0, 1.0), 1.0)
    Random.seed!(42)
    p1 = Point(rand(), rand())
    p2 = Point(rand(), rand())
    p3 = Point(rand(), rand())
    c, r = center3pts(p1, p2, p3)
    @test distance(p1, c) ≈ distance(p2, c) ≈ distance(p2, c)
end

@testset "concatenate" begin
    d1 = Drawing(200,100,:svg)
    origin()
    circle(O,60,:fill)
    finish()

    d2 = Drawing(200,200,:svg)
    rect(O,200,200,:fill)
    finish()
    @test typeof(hcat(d1,d2; hpad=10, valign=:top, clip = true)) <: Drawing
    @test typeof(hcat(d1,d2; hpad=10, valign=:middle, clip = false)) <: Drawing
    @test typeof(hcat(d1,d2; hpad=10, valign=:bottom, clip = true)) <: Drawing

    @test typeof(vcat(d1,d2; vpad=10, halign=:left, clip = true)) <: Drawing
    @test typeof(vcat(d1,d2; vpad=10, halign=:center, clip = false)) <: Drawing
    @test typeof(vcat(d1,d2; vpad=10, halign=:right, clip = true)) <: Drawing

end

println("...finished basic-test")


text:
# simple package

include("../packages/SimplePackage/src/SimplePackage.jl")

@test isdefined(SimplePackage, :hello)

@test SimplePackage.hello() == "Hello, World!"


text:
using JetPack, Jets, Test

n1,n2,n3 = 10,11,12

@testset "JopCircShift" begin
    @testset "CircShift linearop is linear, size=$(N), T=$(T)" for T in (Float64,Float32), N in ((n1,n2), (n1,n2,n3))
        op = JopCircShift(JetSpace(T,N), (1,2,3)[1:length(N)])
        m1 = -1 .+ 2 .* rand(domain(op))
        m2 = -1 .+ 2 .* rand(domain(op))
        @test (op*m1 .+ op*m2) ≈ op*(m1 .+ m2)
        d1 = -1 .+ 2 .* rand(range(op))
        d2 = -1 .+ 2 .* rand(range(op))
        @test (op'*d1 .+ op'*d2) ≈ op'*(d1 .+ d2)
    end

    @testset "CircShift dot product test, size=$(N), T=$(T)" for T in (Float64,Float32), N in ((n1,n2), (n1,n2,n3))
        op = JopCircShift(JetSpace(T,N), (1,2,3)[1:length(N)])
        lhs,rhs = dot_product_test(op, rand(domain(op)), rand(range(op)))
        @test isapprox(lhs, rhs, rtol=1e-6)
    end
end


text:
include("sampler_test_utils.jl")

describe("Constrain Parameters Sampler") do

	describe("construction") do

		subsA = DataGenerators.CategoricalSampler(2)

		test("numparams, paramranges, getparams") do
			cprs = [(0.2, 0.7), (0.3, 0.6)]
			s = DataGenerators.ConstrainParametersSampler(subsA, cprs)
			@check DataGenerators.paramranges(s) == cprs
			@check DataGenerators.numparams(s) == DataGenerators.numparams(subsA)
			@check DataGenerators.getparams(s) == DataGenerators.getparams(subsA)
		end

		test("wrong number of parameter ranges") do
			cprs = [(0.2, 0.7)]
			@check_throws DataGenerators.ConstrainParametersSampler(subsA, cprs)
			cprs = [(0.2, 0.7), (0.3, 0.6), (0.4, 0.8)]
			@check_throws DataGenerators.ConstrainParametersSampler(subsA, cprs)
		end

		test("invalid order of parameter range bounds") do
			cprs = [(0.2, 0.7), (0.6, 0.3)]
			@check_throws DataGenerators.ConstrainParametersSampler(subsA, cprs)
		end

		test("constrained parameter ranges are outside subsampler parameter") do
			cprs = [(1.2, 1.7), (0.3, 0.6)]
			@check_throws DataGenerators.ConstrainParametersSampler(subsA, cprs)
		end

		test("set params") do
			cprs = [(0.2, 0.7), (0.3, 0.6)]
			s = DataGenerators.ConstrainParametersSampler(subsA, cprs)
			DataGenerators.setparams!(s, [0.5, 0.5])
			@check DataGenerators.getparams(subsA) == [0.5, 0.5]
			DataGenerators.setparams!(s, [0.7, 0.3])
			@check DataGenerators.getparams(subsA) == [0.7, 0.3]
			DataGenerators.setparams!(s, [0.2, 0.6])
			ps = DataGenerators.getparams(subsA)
			@check abs(ps[1]-0.25) <= 1e-5
			@check abs(ps[2]-0.75) <= 1e-5
			@check_throws DataGenerators.setparams!(s, [0.2])
			@check_throws DataGenerators.setparams!(s, [0.2, 0.6, 0.2])
			@check_throws DataGenerators.setparams!(s, [prevfloat(0.2), 0.5])
			@check_throws DataGenerators.setparams!(s, [nextfloat(0.7), 0.5])
			@check_throws DataGenerators.setparams!(s, [0.5, prevfloat(0.3)])
			@check_throws DataGenerators.setparams!(s, [0.5, nextfloat(0.6)])
		end

	end

	describe("sampling") do

		subsA = DataGenerators.NormalSampler([2.0, 1.0])
		s = DataGenerators.ConstrainParametersSampler(subsA, [(1.0,10.0), (1.0,5.0)])

		test("sampling is of subsampler") do
			@check isconsistentnormal(s, DataGenerators.getparams(subsA))
		end

	end

	describe("estimate parameters") do

		test("estimates parameters of subsampler") do

			cprs = [(0.2, 0.8)]

			subs1A = DataGenerators.GeometricSampler()
			s1 = DataGenerators.ConstrainParametersSampler(subs1A, cprs)
			params = [0.25]
			DataGenerators.setparams!(s1, params)

			subs2A = DataGenerators.GeometricSampler()
			s2 = DataGenerators.ConstrainParametersSampler(subs2A, cprs)
			otherparams = [0.6]
			DataGenerators.setparams!(s2, otherparams)

			traces = map(1:100) do i
				x, trace = DataGenerators.sample(s1, (0, 1))
				trace
			end

			estimateparams!(s2, traces)

			@check isconsistentgeometric(subs2A, params[1:1])

		end

	end

end


text:
module TensorDecompositions
    using TensorOperations
    using Distributions
    using ProgressMeter
    using Base.Cartesian
    using StatsBase
    using LinearAlgebra

    export

    # types
    SparseArray,
    TensorDecompositions, PARAFAC2, CANDECOMP, CUR, Tucker,

    # TensorDecomposition methods
    rel_residue,
    compose, compose!,

    # additional tensor algebra methods
    tensorcontractmatrix, tensorcontractmatrices!, tensorcontractmatrices,

    # tensor decomposition algorithms
    hosvd,                   # High-order SVD
    candecomp,               # Canonical polyadic decomposition (CANDECOMP/PARAFAC)
    sshopm,                  # Shifted symmetric higher-order power method
    nncp,                    # Non-negative CANDECOMP
    spnntucker,              # sparse (semi) non-negative Tucker
    tensorcur3,              # Tensor-CUR for 3-mode tensors
    parafac2                 # PARAFAC2 model

    include("array_pool.jl")
    include("utils.jl")
    include("abstract_decomposition.jl")
    include("tucker.jl")
    include("hosvd.jl")
    include("candecomp.jl")
    include("sshopm.jl")
    include("nncp.jl")
    include("tensorcur.jl")
    include("parafac2.jl")
    include("spnntucker.jl")

end


text:
using ITensors
using Test
using Suppressor

include("util.jl")

@testset "SVD Algorithms" begin

  @testset "Matrix With Zero Sing Val" begin
    M = [1.0 2.0 5.0 4.0;
         1.0 1.0 1.0 1.0;
         0.0 0.5 0.5 1.0;
         0.0 1.0 1.0 2.0]
    U,S,V = NDTensors.svd_recursive(M)
    @test norm(U*LinearAlgebra.Diagonal(S)*V'-M) < 1E-13
  end

  @testset "Real Matrix" begin
    M = rand(10,20)
    U,S,V = NDTensors.svd_recursive(M)
    @test norm(U*LinearAlgebra.Diagonal(S)*V'-M) < 1E-12

    M = rand(20,10)
    U,S,V = NDTensors.svd_recursive(M)
    @test norm(U*LinearAlgebra.Diagonal(S)*V'-M) < 1E-12
  end

  @testset "Cplx Matrix" begin
    M = rand(ComplexF64,10,15)
    U,S,V = NDTensors.svd_recursive(M)
    @test norm(U*LinearAlgebra.Diagonal(S)*V'-M) < 1E-13

    M = rand(ComplexF64,15,10)
    U,S,V = NDTensors.svd_recursive(M)
    @test norm(U*LinearAlgebra.Diagonal(S)*V'-M) < 1E-13
  end

  @testset "Regression Test 1" begin
    # Implementation of the SVD was giving low
    # accuracy for this case
    M = rand(2,2,2,2)

    M[:, :, 1, 1] =
    [7.713134067177845   -0.16367628720441685;
     -1.5253996568409225   1.3577749944302373]

    M[:, :, 2, 1] =
     [0.0                -2.1219889218225276;
      -8.320068013774126   0.43565608213298096]

    M[:, :, 1, 2] =
     [0.0  -8.662721825820825;
      0.0  -0.46817091771736885]

    M[:, :, 2, 2] =
     [0.0   0.0;
      0.0  -8.159570989998151]

    t1 = Index(2,"t1")
    t2 = Index(2,"t2")
    u1 = Index(2,"u1")
    u2 = Index(2,"u2")

    T = itensor(M,t1,t2,u1,u2)

    U,S,V = svd(T,(u1,t1))
    @test norm(U*S*V-T)/norm(T) < 1E-10
  end

  @testset "Ill-conditioned matrix" begin
    d = 5000
    i = Index(d, "i")
    T = itensor(make_illconditioned_matrix(dim(i)), i', i)

    @suppress begin
      F = svd(T, i'; alg = "divide_and_conquer")
    end
    # Depending on the LAPACK implementation,
    # this sometimes works so don't test it
    #@test isnothing(F)

    F = svd(T, i'; alg = "qr_iteration")
    @test !isnothing(F)
    @test F.U * F.S * F.V ≈ T
  end

end

nothing


text:
module TestNHModel
  using Test
  using EcologicalNetworks

  @test typeof(nestedhierarchymodel(10, 20)) <: UnipartiteNetwork

  #@test richness(nestedhierarchymodel(10, 20)) == 10
  #@test richness(nestedhierarchymodel(10, 49)) == 10
  #@test richness(nestedhierarchymodel(20, 20)) == 20

  # Need to design test, and nh model doesnt have any argument error yet.
  # @test_throws ArgumentError nestedhierarchymodel(10, 99)
  # @test_throws ArgumentError nestedhierarchymodel(10, 51)
  # @test_throws ArgumentError nestedhierarchymodel(10, 50)
  # @test_throws ArgumentError nestedhierarchymodel(10, 200)
  # @test_throws ArgumentError nestedhierarchymodel(10, 0)

  #empirical_foodweb = EcologicalNetworks.nz_stream_foodweb()[1]
  #@test typeof(nichemodel(empirical_foodweb)) <: UnipartiteNetwork

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Performance testing

import Base: unsafe_getindex
# @inline unsafe_getindex(xs...) = Base.getindex(xs...)

function sumelt(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    for k = 1:n
        for a in A
            s += a
        end
    end
    s
end

function sumeach(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    for k = 1:n
        for I in eachindex(A)
            val = unsafe_getindex(A, I)
            s += val
        end
    end
    s
end

function sumlinear(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    for k = 1:n
        for I in 1:length(A)
            val = unsafe_getindex(A, I)
            s += val
        end
    end
    s
end
function sumcartesian(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    for k = 1:n
        for I in CartesianIndices(size(A))
            val = unsafe_getindex(A, I)
            s += val
        end
    end
    s
end

function sumcolon(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    c = Colon()
    for k = 1:n
        @simd for i = 1:ncols
            val = unsafe_getindex(A, c, i)
            s += first(val)
        end
    end
    s
end

function sumrange(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = 1:nrows
    for k = 1:n
        @simd for i = 1:ncols
            val = unsafe_getindex(A, r, i)
            s += first(val)
        end
    end
    s
end

function sumlogical(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = falses(nrows)
    r[1:4:end] = true
    for k = 1:n
        @simd for i = 1:ncols
            val = unsafe_getindex(A, r, i)
            s += first(val)
        end
    end
    s
end

function sumvector(A, n)
    s = zero(eltype(A)) + zero(eltype(A))
    nrows = size(A, 1)
    ncols = size(A, 2)
    r = rand(1:nrows, 5)
    for k = 1:n
        @simd for i = 1:ncols
            val = unsafe_getindex(A, r, i)
            s += first(val)
        end
    end
    s
end

abstract type MyArray{T,N} <: AbstractArray{T,N} end

struct ArrayLS{T,N} <: MyArray{T,N}  # IndexCartesian
    data::Array{T,N}
end
struct ArrayLSLS{T,N} <: MyArray{T,N}  # IndexCartesian with IndexCartesian similar
    data::Array{T,N}
end
Base.similar(A::ArrayLSLS, ::Type{T}, dims::Tuple{Vararg{Int}}) where {T} = ArrayLSLS(similar(A.data, T, dims))
@inline Base.setindex!(A::ArrayLSLS, v, I::Int...) = A.data[I...] = v
@inline Base.unsafe_setindex!(A::ArrayLSLS, v, I::Int...) = Base.unsafe_setindex!(A.data, v, I...)
Base.first(A::ArrayLSLS) = first(A.data)

struct ArrayLF{T,N} <: MyArray{T,N}  # IndexLinear
    data::Array{T,N}
end
struct ArrayStrides{T,N} <: MyArray{T,N}
    data::Array{T,N}
    strides::NTuple{N,Int}
end
ArrayStrides(A::Array) = ArrayStrides(A, strides(A))

struct ArrayStrides1{T} <: MyArray{T,2}
    data::Matrix{T}
    stride1::Int
end
ArrayStrides1(A::Array) = ArrayStrides1(A, size(A,1))

Base.size(A::MyArray) = size(A.data)

@inline Base.getindex(A::ArrayLF, i::Int) = getindex(A.data, i)
@inline Base.getindex(A::ArrayLF, i::Int, i2::Int) = getindex(A.data, i, i2)
@inline Base.getindex(A::Union{ArrayLS, ArrayLSLS}, i::Int, j::Int) = getindex(A.data, i, j)
@inline Base.unsafe_getindex(A::ArrayLF, indx::Int) = unsafe_getindex(A.data, indx)
@inline Base.unsafe_getindex(A::Union{ArrayLS, ArrayLSLS}, i::Int, j::Int) = unsafe_getindex(A.data, i, j)

@inline Base.getindex(A::ArrayStrides{T,2}, i::Real, j::Real) where {T} = getindex(A.data, 1+A.strides[1]*(i-1)+A.strides[2]*(j-1))
@inline Base.getindex(A::ArrayStrides1, i::Real, j::Real) = getindex(A.data, i + A.stride1*(j-1))
@inline Base.unsafe_getindex(A::ArrayStrides{T,2}, i::Real, j::Real) where {T} = unsafe_getindex(A.data, 1+A.strides[1]*(i-1)+A.strides[2]*(j-1))
@inline Base.unsafe_getindex(A::ArrayStrides1, i::Real, j::Real) = unsafe_getindex(A.data, i + A.stride1*(j-1))

# Using the qualified Base.IndexLinear() in the IndexStyle definition
# requires looking up the symbol in the module on each call.
import Base: IndexLinear
Base.IndexStyle(::Type{T}) where {T<:ArrayLF} = IndexLinear()

if !applicable(unsafe_getindex, [1 2], 1:1, 2)
    @inline Base.unsafe_getindex(A::Array, I...) = @inbounds return A[I...]
    @inline Base.unsafe_getindex(A::MyArray, I...) = @inbounds return A[I...]
    @inline Base.unsafe_getindex(A::SubArray, I...) = @inbounds return A[I...]
    @inline Base.unsafe_getindex(A::BitArray, I1::BitArray, I2::Int) = unsafe_getindex(A, Base.to_index(I1), I2)
end

function makearrays(::Type{T}, sz) where T
    L = prod(sz)
    A = reshape(convert(Vector{T}, [1:L;]), sz)
    AS = ArrayLS(A)
    ASS = ArrayLSLS(A)
    AF = ArrayLF(A)
    Astrd = ArrayStrides(A)
    Astrd1 = ArrayStrides1(A)
    outersz = (sz[1]+1,sz[2]+2)
    B = reshape(convert(Vector{T}, [1:prod(outersz);]), outersz)
    Asub = view(B, 1:sz[1], 2:sz[2]+1)
    Bit = trues(sz)
    (A, AF, AS, ASS, Asub, Bit,)
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
# import RoMEPlotting: plotFactor, reportFactors, plotFactorMeasurements, getTimeEasy

"""
    $SIGNATURES

Plot Pose2Pose2 measurement and predicted values in factor centric manner -- i.e. used with inverse solve.
"""
function plotFactorValues(asMeasured::AbstractMatrix{<:Real},
                          asPredicted::AbstractMatrix{<:Real},
                          fct::Type{T}=Pose2Pose2;
                          title="inv. solve",
                          fctsym="",
                          hdl=[],
                          dist::Vector{Float64}=zeros(1) ) where {T <: FunctorInferenceType}
  #
  PPg = manikde!(asMeasured[1:2,:], Point2)
  PP  = manikde!(asPredicted[1:2,:], Point2)
  # dist[1] = minimum(abs.([kld(PPg, PP)[1]; kld(PP, PPg)[1]]))
  mmd!(dist, asPredicted, asMeasured, SE2_Manifold)
  pt = plotKDE([PPg;PP], c=["blue";"red"], legend=["meas";"pred";], levels=3, title="$title $fctsym,\nmmd(..)=$(round(dist[1],digits=6))")

  pc = plotKDECircular([manikde!(asMeasured[3:3,:], Sphere1);manikde!(asPredicted[3:3,:], Sphere1);], c=["blue";"red";], legend=["meas";"pred";], title="$title $fctsym,\n$(T)")

  push!(hdl, pt)
  push!(hdl, pc)

  hstack(pt, pc)
end

"""
    $SIGNATURES

Calculate the "inverse" SLAM solution to compare measured and predicted noise model samples.
"""
function plotFactorMeasurements(dfg::AbstractDFG,
                                fctsym::Symbol,
                                fct::FunctorInferenceType;
                                hdl=[],
                                dist::Vector{Float64}=zeros(1) )
  #
  @error "plotFactorMeasurements not implemented yet for $(typeof(fct))."
end

function plotFactorMeasurements(dfg::AbstractDFG,
                                fctsym::Symbol,
                                fct::Pose2Point2BearingRange;
                                hdl=[],
                                dist::Vector{Float64}=[0.0;] )
  #
  asMeas, asPred = solveFactorMeasurements(dfg, fctsym)

  pc = plotKDECircular([manikde!(asPred[1:1,:], Sphere1);manikde!(asMeas[1:1,:], Sphere1)], c=["blue";"red";], legend=["meas";"pred";], title="inv. solve, $fctsym,\nPose2Point2BearingRange")
  pcl = plotKDE([manikde!(asPred[1:1,:], ContinuousScalar);manikde!(asMeas[1:1,:], ContinuousScalar)], c=["blue";"red";], legend=["meas";"pred";], title="unwrapped rotation, should be [-pi,pi)")
  pl = plotKDE([manikde!(asPred[2:2,:], ContinuousScalar);manikde!(asMeas[2:2,:], ContinuousScalar)], c=["blue";"red";], legend=["meas";"pred";])

  push!(hdl, pc)
  push!(hdl, pcl)
  push!(hdl, pl)

  hstack(vstack(pc,pcl), pl)
end

function plotFactorMeasurements(dfg::AbstractDFG,
                                fctsym::Symbol;
                                hdl=[],
                                dist::Vector{Float64}=[0.0;]  )
  #
  fct = getFactorType(dfg, fctsym)
  plotFactorMeasurements(dfg, fctsym, fct, hdl=hdl, dist=dist)
end

"""
    $SIGNATURES

Return plot of each specific factor.
"""
function plotFactor(dfg::AbstractDFG,
                    fctsym::Symbol,
                    fct::Pose2Point2Range;
                    hdl=[],
                    dist::Vector{Float64}=Float64[0.0;]  )
  # variables
  vars = ls(dfg, fctsym)

  # the pose
  pose = intersect(vars, ls(dfg, Pose2))[1]
  poin = intersect(vars, ls(dfg, Point2))[1]

  # calculate predicted range (inverse solve)
  ptss = map(x->getPoints(getKDE(dfg, x)), vars)
  dpts = (ptss[1][1:2,:] - ptss[2][1:2,:]).^2
  pred = sqrt.(sum(dpts, dims=1))
  plpr = Gadfly.plot(x=pred, Geom.histogram(density=true), Theme(default_color=colorant"deepskyblue"), Guide.title("predicted"))

  # measured
  # meas = getSamples(fct, length(pred))
  smps = rand(fct.Z, length(pred))
  plme = Gadfly.plot(x=smps, Geom.histogram(density=true), Theme(default_color=colorant"magenta"), Guide.title("measured"))

  # measure kl divergence
  PP = manikde!(smps, ContinuousScalar)
  PPg = manikde!(pred, ContinuousScalar)
  dist[1] = minimum(abs.([kld(PPg, PP)[1]; kld(PP, PPg)[1]]))

  modl = Gadfly.plot(
  Gadfly.layer(x=smps, Geom.histogram(density=true), Theme(default_color=colorant"magenta")),
    Gadfly.layer(x=pred, Geom.histogram(density=true), Theme(default_color=colorant"deepskyblue")),
    Guide.manual_color_key("Legend", ["samples";"predicted"], ["magenta";"deepskyblue"]),
    Guide.title("Pose2Point2Range $fctsym\nmin(|kld...|)=$(round(dist[1],digits=3))"),
    Guide.xlabel("range")
  )

  plxy = plotKDE(map(x->getKDE(dfg, x), vars), dims=[1;2], legend=string.(vars), title="Pose2Point2Range: XY plane")
  # plx = plotKDE([getKDE(dfg, vars[1]);], dims=[1;2], legend=[string(vars[1]);], title="Pose2Point2Range: XY plane")
  # ply = plotKDE([getKDE(dfg, vars[2]);], dims=[1;2], legend=[string(vars[2]);], title="Pose2Point2Range: XY plane")

  # mock projection
  tfg = initfg()
  addVariable!(tfg, pose, Pose2)
  addVariable!(tfg, poin, Point2)
  addFactor!(tfg, [pose;poin], fct, autoinit=false)
  manualinit!(tfg, pose, getKDE(dfg,pose))
  manualinit!(tfg, poin, getKDE(dfg,poin))
  apts=Array{Float64,2}(undef, 2, 0)
  loop =true
  while loop
    try
      apts = approxConv(tfg, ls(tfg,pose)[1], poin)
      loop = false
    catch
      @warn "plotFactor Pose2Point2BearingRange non-linear solve fail on approxConv, retrying"
    end
  end
  plhist2 = Gadfly.plot(x=apts[1,:], y=apts[2,:], Geom.histogram2d)
  spc = mean(smps)
  plt = drawPosesLandms(tfg, point_size=5pt, spscale=0.2*spc)
  union!(plhist2.layers, plt.layers)
  pla = plotKDE(manikde!(apts, Point2().manifolds),levels=3,c=["gray50"])

  # plot pose by itself
  # poseplot = plotPose(dfg, pose)
  posepl1 = plotKDE(dfg, pose, dims=[1;2], c=["green"])
  posepl2 = plotKDECircular(marginal(getKDE(dfg, pose), [3]))
  landmpl = plotKDE(dfg, poin, c=["red"])

  # plot handles
  push!(hdl, landmpl)
  push!(hdl, posepl1)
  push!(hdl, posepl2)
  push!(hdl, modl)
  push!(hdl, plme)
  push!(hdl, plpr)
  push!(hdl, plhist2)

  # put plot together
  botr = vstack(modl, plme, plpr) # 4,5,6
  # bot = hstack(plhist2, botr)     # 7,(4,5,6)
  farl = vstack(landmpl, posepl1, posepl2) # 1,2,3

  hstack(farl, botr, plhist2)     # (1,2,3),(4,5,6),7
end

function plotFactor(dfg::AbstractDFG,
                    fctsym::Symbol,
                    fct::Pose2Point2Bearing;
                    hdl=[],
                    dist::Vector{Float64}=Float64[0.0;]  )
  #

  # variables
  vars = ls(dfg, fctsym)

  # the pose
  pose = intersect(vars, ls(dfg, Pose2))[1]
  poin = intersect(vars, ls(dfg, Point2))[1]

  # basic max estimates
  xp = getKDEMax(getKDE(dfg,pose))
  lp = getKDEMax(getKDE(dfg,poin))

  # convolve the yaw angle with bearing rotation model
  pX = marginal(getBelief(dfg, pose), [3])
  pts = approxConvCircular(pX, fct.bearing) # TODO fix by using approxConvBelief instead

  # draw plots
  measest = manikde!(pts, Sphere1)

  # inverse solve for predicted bearing
  dx = getPoints(getKDE(dfg, poin))[1,:] - getPoints(getKDE(dfg, pose))[1,:]
  dy = getPoints(getKDE(dfg, poin))[2,:] - getPoints(getKDE(dfg, pose))[2,:]
  pred = reshape(atan.(dy,dx), 1,:)

  ppX = manikde!(pred, Sphere1)

  scal = 0.2*norm(xp[1:2]-lp)
  plcl = plotKDECircular( [measest; ppX], logpdf=true, legend=["Meas. Est.";"Predicted"], radix=scal, scale=0.2*scal, rVo=[xp[1:2];0.0] )

  # plot pose and point by itself
  posepl1 = plotKDE(dfg, pose, dims=[1;2], c=["green"])
  posepl2 = plotKDECircular(marginal(getKDE(dfg, pose), [3]))
  landmpl = plotKDE(dfg, poin, c=["red"])

  tfg = initfg()
  addVariable!(tfg, pose, Pose2)
  addVariable!(tfg, poin, Point2)
  addFactor!(tfg, [pose;poin], fct, autoinit=false)
  manualinit!(tfg, pose, getKDE(dfg,pose))
  manualinit!(tfg, poin, getKDE(dfg,poin))
  plt = drawPosesLandms(tfg, point_size=5pt)

  union!(plt.layers, plcl.layers)
  plt.coord = Coord.cartesian(aspect_ratio=1.0)

  # draw line from pose to point
  push!(plt.layers, Gadfly.layer(x=[xp[1];lp[1]], y=[xp[2];lp[2]], Geom.line)[1])

  # plot handles
  push!(hdl, landmpl)
  push!(hdl, posepl1)
  push!(hdl, posepl2)
  push!(hdl, plcl)
  push!(hdl, plt)

  hstack(vstack(landmpl,posepl1, posepl2), vstack(plcl, plt))
end

function plotFactor(dfg::AbstractDFG,
                    fctsym::Symbol,
                    fct::Pose2Point2BearingRange;
                    hdl=[],
                    dist::Vector{Float64}=Float64[0.0;]  )
    #
    hdlb = []

    # variables
    vars = ls(dfg, fctsym)

    # the pose
    pose = intersect(vars, ls(dfg, Pose2))[1]
    poin = intersect(vars, ls(dfg, Point2))[1]

    # plot current pose & point
    # pl_poin = plotKDE(dfg, poin, levels=5, c=["red";])
    # pl_pose = plotPose(dfg, pose, c=["black"])

    # project landmark

    # do range model separately too
    pr = Pose2Point2Range(fct.range)
    plotFactor(dfg, fctsym, pr, hdl=hdl)

    br = Pose2Point2Bearing(fct.bearing)
    plotFactor(dfg, fctsym, br, hdl=hdlb)

    newguid = Gadfly.GuideElement[]
    for guid in hdl[4].guides
      if typeof(guid) == Gadfly.Guide.Title
        push!(newguid, Gadfly.Guide.Title("Pose2Point2BearingRange $fctsym"))
      else
        push!(newguid, guid)
      end
    end
    hdl[4].guides = newguid

    union!(hdl, hdlb)

    # plot a prediction of landmark
    predpoints = approxConv(dfg, fctsym, poin)
    PP = manikde!(predpoints, Point2)
    PPg= getKDE(dfg, poin)
    dist[1] = minimum(abs.([kld(PPg, PP)[1]; kld(PP, PPg)[1]]))
    predLandm = plotKDE([PPg; PP], levels=2, c=["red"; "deepskyblue"], legend=["landmark";"predicted"], title="min(|kld...|)=$(round(dist[1],digits=3))")

    push!(hdl, predLandm)

    # AMP.mmd!(dist, PPg, predpoints)

    # measurement solution
    plotFactorMeasurements(dfg, fctsym, hdl=hdl)

    return hstack(vstack(hdl[1],hdl[2],hdl[3]),vstack(hdl[4],hdl[5],hdl[6]),vstack(predLandm, hdlb[5],hdl[7]), vstack(hdl[14],hdl[15],hdl[16]))
end

function plotFactor(dfg::AbstractDFG,
                    fctsym::Symbol;
                    hdl=[],
                    dist::Vector{Float64}=Float64[0.0;]  )
  #
  plotFactor(dfg, fctsym, getFactorType(dfg, fctsym), hdl=hdl, dist=dist)
end

## Reports

getTimeEasy() = split(split("$(now())", 'T')[end],'.')[1]

# import RoMEPlotting: getTimeEasy, reportFactors
# import ApproxManifoldProducts: mmd!

## Also in IIF >v0.7.9
import IncrementalInference: isMultihypo
isMultihypo(fct) = isa(solverData(fct).fnc.hypotheses, Distribution)

"""
$SIGNATURES

Methods to generate plots of all requested factor user supplied vs. current predicted measurement -- i.e. inverse solve.
"""
reportFactors() = @info "Empty function to support both automated documentation and conditional features (@require)."

#


text:
module Equality

using Test

using Darwin
import Darwin: crossover!, mutate!, selection

using LearningStrategies: learn!, MaxIter, strategy, Verbose
using Distributions

import Base: copy
import Random: AbstractRNG, rand, SamplerType

# This example is taken from Westley Argentum:
# https://github.com/WestleyArgentum/GeneticAlgorithms.jl

struct EqualityMonster
    # a + 2b + 3c + 4d + 5e = 42
    abcde::Vector{Int}
end

EqualityMonster() = EqualityMonster(Vector{Int}(undef, 5))

copy(m::EqualityMonster) = EqualityMonster(copy(m.abcde))
rand(rng::AbstractRNG, ::SamplerType{EqualityMonster}) = EqualityMonster(rand(rng, 0:42, 5))

@fitness function fitness(ent::EqualityMonster)
    # we want the expression `a + 2b + 3c + 4d + 5e - 42`
    # to be as close to 0 as possible
    score = sum(ent.abcde .* (1:5))
    1 / abs(score - 42)
end

const EMCrossover = LiftedCrossover{EqualityMonster, UniformCrossover{Vector{Int}, 2, 1}}

function crossover!((p₁, p₂)::NTuple{2, EqualityMonster}, operator::EMCrossover, generation::Integer)
    EqualityMonster.(crossover!((p₁.abcde, p₂.abcde), operator.inner, generation))
end

const EMMutation = LiftedMutation{EqualityMonster, PointwiseMutation{Vector{Int}}}
const EMMutation2 = LiftedMutation{EqualityMonster, AdditiveMutation{Vector{Int}, DiscreteUniform}}

function mutate!(child::EqualityMonster,
                 operator::Union{EMMutation, EMMutation2},
                 generation::Integer)
    mutate!(child.abcde, operator.inner, generation)
    child
end

function run_with(selection, crossover, mutation, generations)
    initial_population = rand(Individual{EqualityMonster}, 64)
    model = PopulationModel(initial_population, fitness)
    strat = strategy(Verbose(GAStrategy(selection, crossover, mutation)),
                     MaxIter(generations))
    learn!(model, strat)
end

@testset "EqualityMonster" begin
    # @test isinf(fitness(result.population[fittest]))

    # let's go crazy and mutate 20% of the time
    r1 = run_with(PairWithBestSelection{EqualityMonster, 1}(),
                  EMCrossover(),
                  EMMutation(0.2, DiscreteUniform(0, 42)),
                  300)
    @test sum(r1.fittest.genome.abcde .* (1:5)) == 42

    r2 = run_with(TournamentSelection{EqualityMonster, 5, 1, 2}(),
                  EMCrossover(),
                  EMMutation2(0.2, 20, 0, 42),
                  300)
    @test sum(r2.fittest.genome.abcde .* (1:5)) == 42
end

end # module


text:
################################################################################
# back propagation.
hasbackpropagation{F, i, o, tf}(::Type{FullyConnectedLayer{F, i, o, tf}}) = true

#TODO - make all of the i/o value parameters to force recompilation for each
#value.

function reversematrixfma{F}(input_array::AbstractVector{F}, matrix::AbstractMatrix{F}, output_array::AbstractVector{F}, i, o)
  for idx = 1:i
      input_array[idx] = zero(F)
    for jdx = 1:o
      input_array[idx] += matrix[jdx, idx] * output_array[jdx]
    end
  end
end

function scaledouterproductfma{F}(matrix::AbstractMatrix{F}, output_deltas::AbstractVector{F}, input_array::AbstractVector, alpha::F, i, o)
  for idx = 1:o
    for jdx = 1:i
      matrix[idx,jdx] -=  alpha * (output_deltas[idx] * input_array[jdx]) # + lambda * mlp.transitions[l_idx - 1])
    end
  end
end

function scaledsubtract{F}(target_vector::AbstractVector{F}, value_vector::AbstractVector{F}, alpha::F, o)
  for idx = 1:o
    target_vector[idx] -= alpha * value_vector[idx]
  end
end

function dxasychainrule{F}(outer_differential::AbstractVector{F}, inner_value::AbstractVector{F}, f::Function, o)
  for idx = 1:o
    outer_differential[idx] = outer_differential[idx] .* (dxasy(f))(inner_value[idx])
  end
end

@generated function backpropagate!{F, i, o, tf}(fcl::FullyConnectedLayer{F, i, o, tf},
#==#                                            input::AbstractArray,
#==#                                            output::AbstractArray{F},
#==#                                            output_deltas::AbstractArray{F},
#==#                                            input_deltas::VoidableDeltas{F} = nothing)

  code = quote
    #for now.
    const alpha = F(0.1)

    #overwrite the output delta values with the adjusted values taking the
    dxasychainrule(output_deltas, output, tf, o)

    scaledsubtract(fcl.bias, output_deltas, alpha, o)

    scaledouterproductfma(fcl.transition, output_deltas, input, alpha, i, o)
  end

  (input_deltas == Void) && return code

  #handle the backpropagating delta data.
  quote
    $code
    reversematrixfma(input_deltas, fcl.transition, output_deltas, i, o)
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
struct ColorScheme
    cs::Dict{Symbol, Union{Int, Symbol}}
end
ColorScheme() = ColorScheme(Dict())
ColorScheme(ps::Pair...) = ColorScheme(Dict(ps...))

Base.getindex(cs::ColorScheme, key) = get(cs.cs, key, key)

function Base.print(io::IO, tokens::Vector{Token})
    for token in tokens
        printstyled(io, token.s; color = token.c.c)
    end
end

Base.print(io::IO, cs::ColorScheme, xs...) = Base.print(IOContext(io, :colorscheme => cs), xs...)

function Base.print(ioc::IOContext, tokens::Vector{Token})
    cs = get(ioc, :colorscheme, ColorScheme())
    for token in tokens
        printstyled(ioc, token.s; color = cs[token.c.c])
    end
end


text:
import JSON

include("ClientAI/ClientAI.jl")
using ClientAI

include("messages.jl")
include("websocket.jl")

type TykitingClient
  host::ASCIIString
  port::Int
  name::ASCIIString
  ai_source::AbstractString
  ai::Nullable{AbstractAI}
  ai_module::Module
  team_id::Int
  socket
  function TykitingClient(host="localhost", port=3000, name="bot", ai="ngcbot")
    # load the ai script
    ai_source = "tyckiting_client/ai/"*ai*".jl"
    # check already here, that the file exists
    @assert isfile(ai_source)
    ai_module = nothing
    @checked ai_module = include(ai_source)

    # ensure compilation of the bot functions
    dummy = ai_module.create(0, Config())
    info("running bot in dummy computation to compile")
    @time begin
      # TODO make events believable! This sight event is probably not possible, so it might cause errors for bots
      # that rely on that.
      events = AbstractEvent[SightEvent(1, 0, Position(0,0)), NoActionEvent(0)]
      ai_module.init_round(dummy, AbstractBot[ClientAI.WebSockOwnBot(0, "dummy", 0, true, Position(0,0), 10)], events, 0)
      event_dispatch(dummy, events)
      ai_module.decide(dummy)
    end

    new(host, port, name, ai_source, nothing, ai_module, 0, nothing)
  end
end

function on_connected(client::TykitingClient, message::ConnectedMsg)
  info("Connected to server with id $(message.team_id)")
  client.team_id = message.team_id

  @checked client.ai = client.ai_module.create(message.team_id, message.config)
  send(client, Dict("type" => "join", "teamName" => client.name))
end

function on_start(client::TykitingClient, message::StartMsg)
#   Handles game start event from server.
#    Args:  message: Message from server containing team compositions
  info("Game started")
  if isdefined(client.ai_module, :on_start)
    enemies = vcat([t.bots for t in message.other_teams]...)
    client.ai_module.on_start(get(client.ai), message.you.bots, enemies)
  end
end

function on_events(client::TykitingClient, message::EventsMsg)
  info("Round $(message.round_id)")
  responses = Any[]
  timing = @elapsed try
    # event dispatcher loop
    client.ai_module.init_round(get(client.ai), message.you.bots, message.events, message.round_id)
    event_dispatch(get(client.ai), message.events)
    responses = client.ai_module.decide(get(client.ai))
  catch e
    println(e)
    Base.show_backtrace(STDOUT, catch_backtrace())
  end
  if timing > 0.25
    warn("Computation took longer than 250ms ($(round(Int, 1000*timing))). Limit: 300ms!")
  end

  actions = Dict("type" => "actions", "roundId" => message.round_id, "actions" => map(to_dict, responses))
  send(client, actions)
end

function on_end(client::TykitingClient, message::EndMsg)
  # Handles game end event from server
  #    Args:  message: Message from server containing the winning team id
  if message.winner_team_id == -1
      result = "It's a draw."
  elseif message.winner_team_id == client.team_id
      result = "You win!"
  else
    result = "You loose."
  end

  info("Game ended. $result")
  close(client)
end

function on_message(client::TykitingClient, message::AbstractString)
  json = JSON.parse(message)
  msgtype = json["type"]
  if msgtype == "connected"
    on_connected(client, ConnectedMsg(json))
  elseif msgtype == "start"
    on_start(client, StartMsg(json))
  elseif msgtype == "events"
    on_events(client, EventsMsg(json))
  elseif msgtype == "end"
    on_end(client, EndMsg(json))
  elseif msgtype == "error"
    error(json)
  else
    warn("Unkown message type $msgtype")
  end
end

function run(client::TykitingClient, debug)
  ws = WebSock.run(client.host, client.port, x->on_message(client, x), debug)
  client.socket = ws
  return
end

function send(client::TykitingClient, data::Dict)
  js = JSON.json(data)
  WebSock.send(client.socket, js)
end
close(client::TykitingClient) = WebSock.close(client.socket)

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
using Redpitaya

ip = "rp-f07083.local"
rp = RedPitaya(ip)

# enable LED2
state(rp, "LED2", true)

# ask if LED2 is on
b = state(rp, "LED2")

println("LED2 is $b")


text:
# Set up input data
test_case = Dict(
    "gen_pmin" => [20, 40, 60],
    "gen_pmax" => [50, 100, 150],
    "gencost" => [
        [2, 0, 0, 3, 1, 4, 7]'
        [2, 0, 0, 3, 2, 5, 8]'
        [2, 0, 0, 3, 3, 6, 9]'
    ],
)

# Build expected returns
expected_one_segment = [
    [1, 0, 0, 2, 20, (1 * 20^2 + 4 * 20 + 7), 50, (1 * 50^2 + 4 * 50 + 7)]'
    [1, 0, 0, 2, 40, (2 * 40^2 + 5 * 40 + 8), 100, (2 * 100^2 + 5 * 100 + 8)]'
    [1, 0, 0, 2, 60, (3 * 60^2 + 6 * 60 + 9), 150, (3 * 150^2 + 6 * 150 + 9)]'
]

expected_two_segment = zeros(3, 10)
expected_two_segment[:, 1] .= 1
expected_two_segment[:, 4] .= 3
expected_two_segment[:, 5:6] = expected_one_segment[:, 5:6]
expected_two_segment[1, 7:8] = [35, (1 * 35^2 + 4 * 35 + 7)]
expected_two_segment[2, 7:8] = [70, (2 * 70^2 + 5 * 70 + 8)]
expected_two_segment[3, 7:8] = [105, (3 * 105^2 + 6 * 105 + 9)]
expected_two_segment[:, 9:10] = expected_one_segment[:, 7:8]

expected_all_equal = copy(test_case["gencost"])
expected_all_equal[:, 5:6] .= 0
expected_all_equal[:, 7] = [2707, 20508, 68409]

expected_some_equal = zeros(3, 8)
expected_some_equal[1, 1:7] = expected_all_equal[1, :]
expected_some_equal[2:3, :] = expected_one_segment[2:3, :]

# Actually run tests
@testset "test default" begin
    this_case = copy(test_case)
    gencost_new = REISE._linearize_gencost(this_case)
    @test gencost_new == expected_one_segment
end

@testset "test one segment" begin
    this_case = copy(test_case)
    gencost_new = REISE._linearize_gencost(this_case; num_segments=1)
    @test gencost_new == expected_one_segment
end

@testset "test two segments" begin
    this_case = copy(test_case)
    gencost_new = REISE._linearize_gencost(this_case; num_segments=2)
    @test gencost_new == expected_two_segment
end

@testset "test all Pmin = Pmax" begin
    this_case = copy(test_case)
    this_case["gen_pmin"] = this_case["gen_pmax"]
    gencost_new = REISE._linearize_gencost(this_case; num_segments=3)
    @test gencost_new == expected_all_equal
end

@testset "test some Pmin = Pmax" begin
    this_case = copy(test_case)
    this_case["gen_pmin"][1] = this_case["gen_pmax"][1]
    gencost_new = REISE._linearize_gencost(this_case; num_segments=1)
    @test gencost_new == expected_some_equal
end