text:
lines = readlines("input.txt")
numbers = parse.(Int, split(lines[1], ','))
fields_s = reshape(lines[2:end], 6, :)[2:end, :]
fields = parse.(Int, reshape(reduce(vcat, split.(fields_s)), 5,5,:))
marks = falses(size(fields))

for n = numbers
    marks[findall(x -> x == n, fields)] .= true

    rowhits = mapslices(prod, marks, dims=[1])
    if any(rowhits) # rows
        board = findfirst(isone, rowhits)[3]
        println(sum(@. fields[:,:,board] * !marks[:,:,board]) * n)
        break
    end

    colhits = mapslices(prod, marks, dims=[2])
    if any(colhits)
        board = findfirst(isone, colhits)[3]
        println(sum(@. fields[:,:,board] * !marks[:,:,board]) * n)
        break
    end
end


text:
using Test, Onda, Dates, Random, UUIDs

@testset "pretty printing" begin
    @test repr(TimeSpan(6149872364198, 123412345678910)) ==
          "TimeSpan(01:42:29.872364198, 34:16:52.345678910)"

    signal = Signal([:a, :b, Symbol("c-d")], Nanosecond(3), Nanosecond(Second(12345)),
                    :unit, 0.25, -0.5, Int16, 50.2, Symbol("lpcm.zst"), nothing)
    @test sprint(show, signal, context=(:compact => true)) ==
          "Signal([:a, :b, Symbol(\"c-d\")])"
    @test sprint(show, signal) == """
                                  Signal:
                                    channel_names: [:a, :b, Symbol(\"c-d\")]
                                    start_nanosecond: 3 nanoseconds (00:00:00.000000003)
                                    stop_nanosecond: 12345000000000 nanoseconds (03:25:45.000000000)
                                    sample_unit: :unit
                                    sample_resolution_in_unit: 0.25
                                    sample_offset_in_unit: -0.5
                                    sample_type: Int16
                                    sample_rate: 50.2 Hz
                                    file_extension: Symbol(\"lpcm.zst\")
                                    file_options: nothing"""

    samples = Samples(signal, true,
                      rand(Random.MersenneTwister(0), signal.sample_type, 3, 5))
    @test sprint(show, samples, context=(:compact => true)) == "Samples(3×5 Array{Int16,2})"
    @test sprint(show, samples) == """
                                   Samples (00:00:00.099601594):
                                     signal.channel_names: [:a, :b, Symbol(\"c-d\")]
                                     signal.start_nanosecond: 3 nanoseconds (00:00:00.000000003)
                                     signal.stop_nanosecond: 12345000000000 nanoseconds (03:25:45.000000000)
                                     signal.sample_unit: :unit
                                     signal.sample_resolution_in_unit: 0.25
                                     signal.sample_offset_in_unit: -0.5
                                     signal.sample_type: Int16
                                     signal.sample_rate: 50.2 Hz
                                     signal.file_extension: Symbol(\"lpcm.zst\")
                                     signal.file_options: nothing
                                     encoded: true
                                     data:
                                   3×5 Array{Int16,2}:
                                    20032  4760  27427  -20758   24287
                                    14240  5037   5598   -5888   21784
                                    16885   600  20880  -32493  -19305"""
    annotations = Set(Annotation("val", TimeSpan(i, i + 1)) for i in 1:10)
    recording = Recording(Dict(:test => signal), annotations)
    @test sprint(show, recording) == """
                                     Recording (03:25:45.000000000; 12345.0 seconds)
                                       signals:
                                         :test => Signal([:a, :b, Symbol(\"c-d\")])
                                       annotations: (10 total)"""
    mktempdir() do root
        dataset = Dataset(joinpath(root, "test.onda"); create=true)
        return @test sprint(show, dataset) == "Dataset($(dataset.path), 0 recordings)"
    end
end

ext: jl


text:
using DelftBlue
using Documenter

DocMeta.setdocmeta!(DelftBlue, :DocTestSetup, :(using DelftBlue); recursive=true)

makedocs(;
    modules=[DelftBlue],
    authors="Oriol Colomes <oriol.colomes@gmail.com>",
    repo="https://github.com/oriolcg/DelftBlue.jl/blob/{commit}{path}#{line}",
    sitename="DelftBlue.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://oriolcg.github.io/DelftBlue.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/oriolcg/DelftBlue.jl",
    devbranch="main",
)


text:
using Base.Test
using RealInterface
using SpecialFunctions
using DiffBase

srand(1)

function finitediff(f, x)
    ϵ = cbrt(eps(typeof(x))) * max(one(typeof(x)), abs(x))
    return (f(x + ϵ) - f(x - ϵ)) / (ϵ + ϵ)
end

for f in vcat(RealInterface.UNARY_ARITHMETIC, RealInterface.UNARY_MATH, RealInterface.UNARY_SPECIAL_MATH)
    if DiffBase.hasdiffrule(f, 1)
        deriv = DiffBase.diffrule(f, :x)
        modifier = in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)) ? 1 : 0
        @eval begin
            x = rand() + $modifier
            @test isapprox($deriv, finitediff($f, x), rtol=0.05)
        end
    end
end

for f in vcat(RealInterface.BINARY_MATH, RealInterface.BINARY_ARITHMETIC, RealInterface.BINARY_SPECIAL_MATH)
    if DiffBase.hasdiffrule(f, 2)
        derivs = DiffBase.diffrule(f, :x, :y)
        @eval begin
            x, y = rand(1:10), rand()
            dx, dy = $(derivs[1]), $(derivs[2])

            if !(isnan(dx))
                @test isapprox(dx, finitediff(z -> $f(z, y), float(x)), rtol=0.05)
            end

            if !(isnan(dy))
                @test isapprox(dy, finitediff(z -> $f(x, z), y), rtol=0.05)
            end
        end
    end
end


text: Dict{String,Any}("title" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "Client: \"XXXX\", Job: XXXX","type" => "string"), Dict{String,Any}("value" => "Code: XXXX","type" => "string")],"type" => "array"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
module OperationsTests

using Test
using TensorValues

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

# Comparison

@test a==a
@test a ≈ a
@test a!=b
@test [a,a] == [a,a]
@test [a,a] ≈ [a,a]

# Addition / subtraction

c = +a
r = VectorValue(1,2,3)
@test c == r

c = -a
r = VectorValue(-1,-2,-3)
@test c == r

c = a + b
r = VectorValue(3,3,9)
@test c == r

c = a - b
r = VectorValue(-1,1,-3)
@test c == r

# Matrix Division

t = one(TensorValue{3,Int,9})

c = t\a

@test c == a

# Operations by a scalar

a = VectorValue(1,2,3)
r = VectorValue(2,4,6)

c = 2 * a
@test isa(c,VectorValue{3,Int})
@test c == r

c = a * 2
@test isa(c,VectorValue{3,Int})
@test c == r

r = VectorValue(3,4,5)

c = 2 + a
@test isa(c,VectorValue{3,Int})
@test c == r

c = a + 2
@test isa(c,VectorValue{3,Int})
@test c == r

r = VectorValue(1/2,1.0,3/2)

c = a / 2
@test isa(c,VectorValue{3,Float64})
@test c == r

# Dot product (simple contraction)

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

t = TensorValue(1,2,3,4,5,6,7,8,9)
s = TensorValue(9,8,3,4,5,6,7,2,1)

c = a * b
@test isa(c,Int)
@test c == 2+2+18

c = t * a
@test isa(c,VectorValue{3,Int})
r = VectorValue(30,36,42)
@test c == r

c = s * t
@test isa(c,TensorValue{3,Int})
r = TensorValue(38,24,18,98,69,48,158,114,78)
@test c == r

# Inner product (full contraction)

c = inner(2,3)
@test c == 6

c = inner(a,b)
@test isa(c,Int)
@test c == 2+2+18

c = inner(t,s)
@test isa(c,Int)
@test c == 185

# Reductions

a = VectorValue(1,2,3)

c = sum(a)
@test isa(c,Int)
@test c == 1+2+3

c = maximum(a)
@test isa(c,Int)
@test c == 3

c = minimum(a)
@test isa(c,Int)
@test c == 1

# Outer product (aka dyadic product)

a = VectorValue(1,2,3)
e = VectorValue(2,5)

c = outer(2,3)
@test c == 6

r = VectorValue(2,4,6)
c = outer(2,a)
@test isa(c,VectorValue{3,Int})
@test c == r

c = outer(a,2)
@test isa(c,VectorValue{3,Int})
@test c == r

c = outer(a,e)
@test isa(c,MultiValue{Tuple{3,2},Int})
r = MultiValue{Tuple{3,2},Int}(2,4,6,5,10,15)
@test c == r

# Linear Algebra

t = TensorValue(10,2,30,4,5,6,70,8,9)

c = det(t)
@test c ≈ -8802.0

c = inv(t)
@test isa(c,TensorValue{3})

# Measure

a = VectorValue(1,2,3)
c = meas(a)
@test c ≈ 3.7416573867739413

t = TensorValue(10,2,30,4,5,6,70,8,9)
c = meas(t)
@test c ≈ 8802.0

# Broadcasted operations

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

A = Array{VectorValue{3,Int},2}(undef,(4,5))
A .= a
@test A == fill(a,(4,5))

C = .- A

r = VectorValue(-1,-2,-3)
R = fill(r,(4,5))
@test C == R

B = fill(b,(4,1))

C = A .+ B

r = VectorValue(3,3,9)
R = fill(r,(4,5))
@test C == R

C = A .- B

r = VectorValue(-1,1,-3)
R = fill(r,(4,5))
@test C == R

# conj

v = VectorValue(1,0)
@test v == v'
@test conj(v) == v'

t = TensorValue(1,2,3,4)
@test trace(t) == 5
@test tr(t) == 5

t = TensorValue(1,2,3,4,5,6,7,8,9)
@test trace(t) == 15
@test tr(t) == 15

@test symmetic_part(t) == TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0)

a = TensorValue(1,2,3,4)
b = a'
@test adjoint(a) == b
@test b == TensorValue(1,3,2,4)
@test a*b == TensorValue(10,14,14,20)

u = VectorValue(1.0,2.0)
v = VectorValue(2.0,3.0)
@test dot(u,v) ≈ inner(u,v)
@test norm(u) ≈ sqrt(inner(u,u))

end # module OperationsTests


text:
using PackageCompiler
PackageCompiler.create_sysimage(:OhMyREPL; precompile_statements_file="ohmyrepl_precompile.jl", replace_default=true)


text:
# Model m4.2s.jl

using Pkg, DrWatson

begin
    using StanQuap
    using StatisticalRethinking
end

begin
    df = CSV.read(sr_datadir("Howell1.csv"), DataFrame; delim=';')
    df = filter(row -> row[:age] >= 18, df);
end;

stan4_2 = "
data {
  int N;
  real<lower=0> height[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 0.1);
  sigma ~ uniform(0 , 50);

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

data = (N = length(df.height), height = df.height)
init = (mu = 180.0, sigma = 10.0)
q4_2s, m4_2s, o4_2s = stan_quap("m4.2s", stan4_2; data, init);

if q4_2s.converged
  quap4_2s_df = sample(q4_2s)          # DataFrame with samples
  precis(quap4_2s_df)
  post4_2s = read_samples(m4_2s)
  post4_2s |> display
end

# End of m4.2s.jl


text:
struct Simulation{DT <: Number}
    Δt::DT
    x::OffsetArray{DT,3,Array{DT,3}}

    function Simulation(x₀::AbstractArray{DT,2}, Δt::DT, nt::Int) where {DT}
        x = OffsetArray(zeros(DT, size(x₀)..., nt+1), axes(x₀)..., 0:nt)
        x[:,:,0] .= x₀
        new{DT}(Δt, x)
    end
end

ntimesteps(sim::Simulation) = lastindex(sim.x,3)
nsamples(sim::Simulation) = lastindex(sim.x,2)

eachtimestep(sim::Simulation) = axes(sim.x,3)[1:end]
eachsample(sim::Simulation) = axes(sim.x,2)

function run!(sim::Simulation)
    for n in eachtimestep(sim)
        for i in eachsample(sim)
            sim.x[1,i,n] = sim.x[1,i,n-1] + sim.Δt * sim.x[2,i,n-1]
            sim.x[2,i,n] = sim.x[2,i,n-1] - sim.Δt * sin(sim.x[1,i,n-1])
        end
    end
    return sim.x
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
#############################################################################
# lp_over_assoccpoly.jl
# solves linear programming over associated polyhedra.
#############################################################################

export solveLPOA!

function solveLPOA!(objective::AbstractExpr, p::AssocPoly)
  if objective.head == :*
    if isa(objective.children[2], Variable)
      if isa(objective.children[1], Constant)
        objective.children[2].value = fenchel(p, objective.children[1].value[1, :])
      else
        objective.children[2].value = fenchel(p, evaluate(objective.children[1]))
      end
      return objective.children[2].value
    elseif isa(objective.children[1], Variable)
      if isa(objective.children[2], Constant)
        objective.children[1].value = fenchel(p, objective.children[2].value[1, :])
      else
        objective.children[1].value = fenchel(p, evaluate(objective.children[2]))
      end
      return objective.children[1].value
    end
  end
end


text:
# [try, catch](https://docs.julialang.org/en/v1/manual/control-flow/#The-try/catch-statement-1) can be used to deal with errors as with many other languages
try
    push!(a,1)
catch err
    showerror(stdout, err, backtrace());println()
end
#> UndefVarError: a not defined
#> Stacktrace:
#>  [1] top-level scope at C:\JuliaByExample\src\error_handling.jl:5
#>  [2] include at .\boot.jl:317 [inlined]
#>  [3] include_relative(::Module, ::String) at .\loading.jl:1038
#>  [4] include(::Module, ::String) at .\sysimg.jl:29
#>  [5] exec_options(::Base.JLOptions) at .\client.jl:229
#>  [6] _start() at .\client.jl:421
println("Continuing after error")
#> Continuing after error

max_stars_repo_name: talgatomarov/JuliaByExample

max_issues_repo_name: talgatomarov/JuliaByExample

max_forks_repo_name: talgatomarov/JuliaByExample


text:
for i in 1:1000
    println("Hello World")
end


text:
mutable struct StepDate{T,F} <: AbstractStep
    selections::T
    params::Union{Nothing,Vector{Symbol}}
    features::F
    trained::Bool
    skip::Bool
    keep::Bool
    roles::Union{Nothing,Symbol,Vector{Symbol}}
end
function fit!(step::StepDate, df)
    step.params = getselectionkeys(df, step.selections)
    step.trained = true
end

function transform!(s::StepDate, df)
    for k in s.params
        for f in s.features
            df[Symbol(string(k) * "_" * funname(f))] = f.(df[k])
        end
        if !s.keep
            deletecols!(df, k)
        end
    end
    df
end
function transform(s::StepDate, df)
    transform!(s, copy(df))
end

function update!(roles::Dict, step::StepDate)
    for k in step.params
        kroles = whichkeys(roles, k)
        newroles = step.roles === nothing ? kroles : asarray(step.roles)
        if newroles !== nothing
            for newrole in newroles
                !(newrole in keys(roles)) && get!(roles, newrole, Symbol[])
                for f in step.features
                    push!(roles[newrole], Symbol(string(k) * "_" * funname(f)))
                end
            end
        end
        if !step.keep
            if kroles !== nothing
                for role in kroles
                    roles[role] = filter(x -> !(k in x), roles[role])
                end
            end
        end
    end
    roles
end

function step_date!(r::Recipe, s...; features=[year, month, dayofweek], keep::Bool=false, roles=nothing, skip::Bool=false)
    push!(r.steps,
          StepDate([s...],
                   nothing,
                   asarray(features),
                   false,
                   skip,
                   keep,
                   roles))
end
function step_date!(r::Recipe; features=[year, month, dayofweek], keep::Bool=false, roles=nothing, skip=false)
    push!(r.steps,
          StepDate([if_eltype(Date)],
                   nothing,
                   asarray(features),
                   false,
                   skip,
                   keep,
                   roles))
end

function step_datetime!(r::Recipe, s...; features=[year, month, dayofweek, hour, minute, second], keep::Bool=false, roles=nothing, skip::Bool=false)
    push!(r.steps,
          StepDate([s...],
                   nothing,
                   asarray(features),
                   false,
                   skip,
                   keep,
                   roles))
end
function step_datetime!(r::Recipe; features=[year, month, dayofweek, hour, minute, second], keep::Bool=false, roles=nothing, skip::Bool=false)
    push!(r.steps,
          StepDate([if_eltype(DateTime)],
                   nothing,
                   asarray(features),
                   false,
                   skip,
                   keep,
                   roles))
end


text:
!isempty(ARGS) || error("No experiment config supplied.")
isfile(ARGS[1]) || error("Cannot read '$(ARGS[1])'")
isabspath(ARGS[1]) || error("Please use an absolute path for the experiment config.")
@info "Experiment config supplied: '$(ARGS[1])'"
config_file = ARGS[1]

worker_config = "$(@__DIR__)/../config/config.jl"
@info "Loading worker config: $worker_config"
include(worker_config)

using JuMP, Gurobi
using SVDD, OneClassActiveLearning
using Serialization, Logging, Random

push!(LOAD_PATH, joinpath(@__DIR__, "../src/KernelLearning"))
using KernelLearning

function init_result_dir(exp_dir)
    if isdir(exp_dir)
        print("Type 'yes' or 'y' to delete and overwrite experiment $(exp_dir): ")
        argin = readline()
        if argin == "yes" || argin == "y"
            rm(exp_dir, recursive=true)
        else
            error("Overwriting anyways... Just kidding, nothing happened.")
        end
    else
        mkpath(exp_dir)
        mkpath(joinpath(exp_dir, "results"))
        mkpath(joinpath(exp_dir, "log"))
    end
end

function generate_experiments(methods, data_files, solver, oracle_param, quality_metrics)
    exp_dir = joinpath(data_output_root, exp_name)
    init_result_dir(exp_dir)

    experiments = []
    for data_file in data_files
        out_dir = split(data_file, '/')[end-1]
        output_path = joinpath(exp_dir, "results", out_dir)
        isdir(output_path) || mkdir(output_path)

        data, labels = load_data(data_file)

        oracle = OneClassActiveLearning.initialize_oracle(oracle_param[:type], data, labels, oracle_param[:param])
        for seed in seeds
            for method in methods
                experiment = Dict{Symbol, Any}(
                    :data_file => data_file,
                    :data_set_name => out_dir,
                    :oracle => oracle,
                    :method => method,
                    :seed => seed,
                    :solver => solver,
                    :quality_metrics => quality_metrics,
                    :selection_metric => selection_metric)
                exp_hash = hash(sprint(print, experiment))
                @assert exp_hash == hash(sprint(print, deepcopy(experiment)))
                experiment[:hash] = "$exp_hash"

                out_name = splitext(splitdir(data_file)[2])[1]
                out_name = joinpath(output_path, "$(out_name)_$(string(method))_$(exp_hash).json")
                experiment[:output_file] = out_name
                experiment[:log_dir] = joinpath(exp_dir, "log")
                push!(experiments, deepcopy(experiment))
            end
        end
    end

    cp(@__FILE__, joinpath(exp_dir, splitdir(@__FILE__)[2]), follow_symlinks=true)
    @info "Created $exp_dir with $(length(experiments)) instances."

    open(joinpath(exp_dir, "experiment_hashes"), "a") do f
        for e in experiments
            write(f, "$(e[:hash])\n")
        end
    end
    serialize(joinpath(exp_dir, "experiments.jser"), experiments)
    return nothing
end

include(config_file)

all(isdir.(joinpath.(data_input_root, data_dirs))) || error("Not all data dirs are valid.")
data_files = vcat(([joinpath.(data_input_root, x, readdir(joinpath(data_input_root, x))) for x in data_dirs])...)
@info("Found $(length(data_files)) data files.")

generate_experiments(learning_methods, data_files, solver, oracle_param, quality_metrics)

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
using Plots
#plotlyjs()
pyplot()

x_min, x_max = -5, 5
y_min, y_max = -5, 5

alpha, beta = 0.2, 0.1

# Axes
gs = 3
z = linspace(x_min, x_max, gs)
x = zeros(gs)
y = zeros(gs)
plot(x, y, z, color=:black, linewidth=2, alpha=0.5, label="")
plot!(z, x, y, color=:black, linewidth=2, alpha=0.5, label="")
plot!(y, z, x, color=:black, linewidth=2, alpha=0.5, label="")

# Fixed linear function, to generate a plane
f(x, y) = alpha .* x + beta .* y

# Vector locations, by coordinate
x_coords = [3, 3]
y_coords = [4, -4]
z = f(x_coords, y_coords)

# Lines to vectors
n = 2
x_vec = zeros(n, n)
y_vec = zeros(n, n)
z_vec = zeros(n, n)
labels = []
for i=1:n
  x_vec[:, i] = [0; x_coords[i]]
  y_vec[:, i] = [0; y_coords[i]]
  z_vec[:, i] = [0; f(x_coords[i], y_coords[i])]
  lab = string("a", i)
  push!(labels, lab)
end
plot!(x_vec, y_vec, z_vec, color=[:blue :red], linewidth=1.5, alpha=0.6, label=reshape(labels,1,length(labels)))

# Draw the plane
grid_size = 20
xr2 = linspace(x_min, x_max, grid_size)
yr2 = linspace(y_min, y_max, grid_size)
z2 = Array{Float64}(grid_size, grid_size)
for i in 1:grid_size
    for j in 1:grid_size
        z2[j, i] = f(xr2[i], yr2[j])
    end
end
surface!(xr2, yr2, z2, cbar=false, alpha=0.2, fill=:blues, xlims=(x_min, x_max), ylims=(x_min, x_max), zlims=(x_min, x_max), xticks=[0], yticks=[0], zticks=[0])


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to process and print the table title.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _print_title(buf::IO, title_tokens::Vector{String}, has_color::Bool, title_crayon::Crayon)

Print the table title to the buffer `buf`.
"""
function _print_title(
    buf::IO,
    title_tokens::Vector{String},
    # Configurations
    has_color::Bool,
    title_crayon::Crayon
)
    num_tokens = length(title_tokens)

    num_tokens == 0 && return nothing

    has_color && print(buf, title_crayon)

    @inbounds for i = 1:num_tokens
        print(buf, rstrip(title_tokens[i]))

        # In the last line we must not add the new line character
        # because we need to reset the crayon first if the display
        # supports colors.
        i != num_tokens && println(buf)
    end

    has_color && print(buf, _reset_crayon)
    println(buf)

    return nothing
end

"""
    _tokenize_title(title::AbstractString, display_width::Int, table_width::Int, title_alignment::Symbol, title_autowrap::Bool, title_same_width_as_table::Bool)

Split the table title into tokens considering the line break character.
"""
function _tokenize_title(
    title::AbstractString,
    display_width::Int,
    table_width::Int,
    # Configurations
    title_alignment::Symbol,
    title_autowrap::Bool,
    title_same_width_as_table::Bool
)
    # Process the title separating the tokens.
    title_tokens = String[]

    if length(title) > 0
        # Compute the title width.
        title_width = title_same_width_as_table ? table_width : display_width

        # If the title width is not higher than 0, then we should only print the
        # title.
        if title_width ≤ 0
            push!(title_tokens, title)

            # Otherwise, we must check for the alignments.
        else
            title_tokens_raw = string.(split(title, '\n'))
            title_autowrap && (title_tokens_raw = _str_autowrap(title_tokens_raw, title_width))
            num_tokens = length(title_tokens_raw)

            @inbounds for i = 1:num_tokens
                token = title_tokens_raw[i]
                token_str = _str_aligned(token, title_alignment, title_width)
                push!(title_tokens, token_str)
            end
        end
    end

    return title_tokens
end


text:
# Conditional HINT network from Kruse et al. (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random
Random.seed!(11)

# Define network
nx = 64
ny = 64
n_in = 2
n_hidden = 4
batchsize = 2
L = 2
K = 2

# Multi-scale and single scale network
CH0 = NetworkMultiScaleConditionalHINT(nx, ny, n_in, batchsize, n_hidden, L, K; split_scales=false, k1=3, k2=1, p1=1, p2=0)
CH1 = NetworkConditionalHINT(nx, ny, n_in, batchsize, n_hidden, L*K; k1=3, k2=1, p1=1, p2=0)

nets = [CH0, CH1, reverse(CH1)]

function test_inv(CH, nx, ny, n_in)
    print("\nInvertibility test HINT network\n")
    # Test layers
    test_size = 10
    X = randn(Float32, nx, ny, n_in, test_size)
    Y = X + .1f0*randn(Float32, nx, ny, n_in, test_size)

    # Forward-backward
    Zx, Zy, logdet = CH.forward(X, Y)
    X_, Y_ = CH.backward(0f0.*Zx, 0f0.*Zy, Zx, Zy)[3:4]
    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

    # Forward-inverse
    Zx, Zy, logdet = CH.forward(X, Y)
    X_, Y_ = CH.inverse(Zx, Zy)
    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-3)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-3)

    # Y-lane only
    Zyy = CH.forward_Y(Y)
    Yy = CH.inverse_Y(Zyy)
    @test isapprox(norm(Y - Yy)/norm(Y), 0f0; atol=1f-3)
end

# Loss
function loss(CH, X, Y)
    Zx, Zy, logdet = CH.forward(X, Y)
    f = -log_likelihood(tensor_cat(Zx, Zy)) - logdet
    ΔZ = -∇log_likelihood(tensor_cat(Zx, Zy))
    ΔZx, ΔZy = tensor_split(ΔZ)
    ΔX, ΔY = CH.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    return f, ΔX, ΔY
end

function test_grad(CH, nx, ny, n_in)
    test_size = 10
    X = randn(Float32, nx, ny, n_in, test_size)
    Y = X + .1f0*randn(Float32, nx, ny, n_in, test_size)
    X0 = randn(Float32, nx, ny, n_in, test_size)
    Y0 = X0 + .1f0*randn(Float32, nx, ny, n_in, test_size)
    dX = X - X0
    dY = Y - Y0

    f0, ΔX, ΔY = loss(CH, X0, Y0)
    h = 0.1f0
    maxiter = 6
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)

    print("\nGradient test cond. HINT net: input\n")
    for j=1:maxiter
        f = loss(CH, X0 + h*dX, Y0 + h*dY)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX, ΔX) - h*dot(dY, ΔY))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

# Loop over networks and reversed counterpart
for CH in nets
    # Invertibility
    test_inv(CH, nx, ny, n_in)
    test_grad(CH, nx, ny, n_in)
end


text:
# Solves for the temporary equilibrium given factor price guess (x in eq 5 on p 750 → xdot in eq 11 on p 754)
function solve_temp_equilibrium(
    ω, ω0, phat, labor, VARjn0, VALjn0, λ, λ_0, λ_1, Sn_prime,
    κ_hat, Z_hat, α, io, θinv, θinv_stack, Ψ, G, G_vec, γ,
    mp, ap
)

    ωmax = 1
    itw = 1

    λ_prime = similar(λ)
    X_prime = similar(VARjn0)

    while (itw <= ap.maxit) & (ωmax > ap.tol_in)

        # calculate good prices and input bundle consistent with factor prices [eq (11) and (12) on page 754]
        phat, xhat = solve_prices(phat, ω, ω0, κ_hat, Z_hat, θinv, θinv_stack, G, G_vec, γ, λ, λ_0, λ_1, mp, ap)

        # calculate bilateral trade shares [eq (13) on page 754]
        λ_prime = solve_λprime(λ, λ_0, λ_1, κ_hat, Z_hat, xhat, phat, θinv, θinv_stack, γ, mp)

        # calculate total expenditures [eq (14) on page 754]
        X_prime = solve_expenditures(α, Ψ, G, λ_prime, ω, labor, Sn_prime, VARjn0, VALjn0, io, mp)

        # calculate factot prices [eq (15) on page 754]
        ω_new = solve_factor_prices(X_prime, λ_prime, Ψ, γ, labor, VARjn0, VALjn0, mp)

        ωmax = maximum(abs.((ω_new .- ω))./ω) # checking tolerance

        # new factor prices
        ω1 = ap.damp*ω_new .+ (1 - ap.damp)*ω

        ## CDP error criterion
        # ωmax_usa = sum(abs.((ω1[:,1:mp.R] .- ω[:,1:mp.R])).^2)
        # ωmax_row = sum(abs.((ω1[1,mp.R+1:mp.N] .- ω[1,mp.R+1:mp.N])).^2)
        # ωmax = ωmax_usa + ωmax_row

        if mod(itw, 5) == 0
            println("Factor price error: $(ωmax).")
        end

        ω = ω1
        itw += 1

    end

    # recovering wages and rental rates frω matrix of factor prices ω
    wages = similar(ω)
    rental_rate = similar(ω)
    wages[:,1:mp.R] = ω[:,1:mp.R].*(labor[:,1:mp.R].^(-Ψ[:,1:mp.R]))
    wages[:,mp.R+1:mp.N] = ω[:,mp.R+1:mp.N]  # outside the US, wage dot = rent dot frω footnote 63 on p.811, we assume no reallocation of labor outside the US
    rental_rate[:,1:mp.R] = wages[:,1:mp.R].*labor[:,1:mp.R]
    rental_rate[:,mp.R+1:mp.N] = ω[:,mp.R+1:mp.N]   # outside the US, wage dot = rent dot

    # recovering deficits
    VARjn_prime = (VARjn0.*ω .*(labor.^(1 .- Ψ)))
    VAR_prime = sum(VARjn_prime, dims = 1)'
    chi_prime = sum(VAR_prime)
    Bn_prime = Sn_prime .- io.*chi_prime + VAR_prime

    # recovering X
    X_prime = reshape(X_prime', 1, mp.J*mp.N)'
    xbilat_prime = repeat(X_prime, 1, mp.N).*λ_prime

    # recovering other variables
    VALjn_prime = wages .* labor .* VALjn0
    VAjn_prime = VALjn_prime .+ VARjn_prime
    VAjn_prime = VAjn_prime[:,1:mp.R]

    Phat = prod(phat.^(α), dims = 1)

    return ω, wages, Phat, rental_rate, phat, VARjn_prime, VALjn_prime, X_prime, λ_prime, xbilat_prime

end

function solve_factor_prices(X_prime, λ_prime, Ψ, γ, labor, VARjn0, VALjn0, mp)

    X = reshape(X_prime', 1, mp.J*mp.N)'

    DP = λ_prime.*X

    Exjnp = similar(γ)
    for j=1:mp.J
        for n=1:mp.N
            Exjnp[j,n] = sum(DP[1+mp.N*(j-1):mp.N*j,n]) # rightmost term in eqn (15)
        end
    end

    aux = γ.*Exjnp
    ω0 = ones(mp.J,mp.N)
    ω0[:,1:mp.R] = aux[:,1:mp.R]./((labor[:,1:mp.R].^(1 .- Ψ[:,1:mp.R])).*(VARjn0[:,1:mp.R] + VALjn0[:,1:mp.R]))
    VAR = sum(VARjn0, dims = 1)
    VAL = sum(VALjn0, dims = 1)
    aux = sum(aux, dims = 1)
    ω0[:,mp.R+1:mp.N] .= (aux[:,mp.R+1:mp.N]./(VAR[:,mp.R+1:mp.N] + VAL[:,mp.R+1:mp.N]))

    return ω0

end

function solve_expenditures(α, Ψ, G, λ_prime, ω, labor, Sn_prime, VARjn0, VALjn0, io, mp)

    # counterfactual equilibrium structures
    VARjn_prime = VARjn0 .* ω .* (labor.^(1 .- Ψ))
    VAR_prime = sum(VARjn_prime, dims = 1)'

    # total revenues in the global portfolio
    chi_prime = sum(VAR_prime)

    # new trade balance
    Bn_prime = Sn_prime .- io.*chi_prime + VAR_prime

    # Calculate the Ω matrix as in CP 2015.
    NBP = similar(λ_prime')

    for n_out in 1:mp.N
        for n_in in 1:mp.N
            @inbounds NBP[n_out, 1 + (n_in-1)*mp.J : n_in*mp.J] = λ_prime[n_in:mp.N:mp.J*mp.N, n_out]
        end
    end

    NNBP = kron(NBP, ones(mp.J, 1))
    GG = kron(ones(1,mp.N), G')
    GP = GG.*NNBP
    Ω = I - GP

    # calculate total expenditures
    aux = sum(ω.*(labor.^(1 .- Ψ)).*(VARjn0 + VALjn0), dims = 1)' .- Bn_prime
    aux = kron(aux, ones(mp.J,1))
    X = Ω \ (reshape(α, mp.N * mp.J, 1) .* aux)
    X_prime = reshape(X, mp.J, mp.N)

    return X_prime

end

function solve_λprime(λ, λ_0, λ_1, κ_hat, Z_hat, xhat, phat, θinv, θinv_stack, γ, mp)

    # Calculating bilateral trade shares using eqn (13)
    xhat_prime = xhat.^(-1 ./ θinv)
    phat_prime = phat.^(-1 ./ θinv)

    λ_k = ( λ_1 ./ λ_0 ) .* λ .* (κ_hat.^(-1 ./ θinv_stack))
    λ_k[isnan.(λ_k)] .= 0

    Dλ = similar(κ_hat)
    λ_prime = similar(Dλ)
    for n in 1:mp.N
        index = n:mp.N:size(θinv_stack,1)-(mp.N-n)
        Dλ[index,:] = λ_k[index,:] .* (xhat_prime.*(Z_hat.^(γ./(repeat(θinv, 1, mp.N)))))
        λ_prime[index,:] = Dλ[index,:]./repeat(phat_prime[:,n], 1, mp.N);
    end

    return λ_prime

end

function solve_prices(phat_guess, ω, ω0, κ_hat, Z_hat, θinv, θinv_stack, G, G_vec, γ, λ, λ_0, λ_1, mp, ap)

    log_p_guess = [phat_guess[:,i] for i = 1:mp.N]
    xhat = similar(ω, Float64)
    phat = similar(ω, Float64)
    pfmax = 1.
    it = 1
    GP = [ones(mp.J) for i = 1:mp.N]

    # create vector of vectors
    log_ω = log.(ω./ω0)
    # preallocate outside loop
    log_xhat = similar(ω, Float64)
    # first term in log(xhat), can be precomputed
    log_xhat_1 = similar(ω, Float64)

    for n in 1:mp.N
        log_xhat_1[:,n] = γ[:,n].*log_ω[:,n]
    end

    while (it <= ap.maxit) & (pfmax > ap.tol_in)

        # G*log(p), use mul! to make faster
        GP_array = fastmultiply(GP, G_vec, log_p_guess)

        # calculate input bundle costs [log(xhat) frω eq (11) on p 754]
        log_xhat = log_xhat_1 .+ GP_array

        xhat = exp.(log_xhat)

        # calculate change in prices [eq (12) on p 754]
        λ_k = ( λ_1 ./ λ_0 ) .* λ .* (κ_hat.^(-1 ./ θinv_stack))
        λ_k[isnan.(λ_k)] .= 0

        for j in 1:mp.J
            for n in 1:mp.N
                phat[j,n] = λ_k[n+(j-1)*mp.N,:]' *
                    (Z_hat[j,:].^(gamma[j,:]./(θinv[j])).*(xhat[j,:].^(-1/θinv[j])))
                phat[j,n] = phat[j,n]^(-θinv[j])
            end
        end

        ## CDP error criterion
        pfmax = maximum(abs.(phat - phat_guess))

        phat_guess = deepcopy(phat)
        log_p_guess = [log.(phat_guess[:,n]) for n = 1:size(phat_guess,2)]
        it += 1

    end

    return phat_guess, xhat

end

max_stars_repo_name: rafserqui/class-repo

max_issues_repo_name: rafserqui/class-repo

max_forks_repo_name: rafserqui/class-repo


text:
using Luxor

const julia_blue    = (0.251, 0.388, 0.847)
const julia_purple  = (0.584, 0.345, 0.698)
const julia_green   = (0.22, 0.596, 0.149)
const julia_red     = (0.796, 0.235, 0.2)

function draw_full(type)
    currentwidth = 900 # pts
    currentheight = 250 # pts
    Drawing(currentwidth, currentheight, "logo.$type")

    sethue(julia_red)
    circle(Point(40, 30), 25, :fill)

    sethue(julia_purple)
    circle(Point(115, 30), 25, :fill)

    sethue(julia_green)
    circle(Point(830, 30), 25, :fill)

    # origin()

    fontsize(300)
    fontface("Literata SemiBold")

    str = "actScı" # Uses "dotless i"

    translate(0, 225)

    textpath(str)                             # make text into a path
    setline(3)
    setcolor("black")
    fillpreserve()                            # fill but keep

    finish()
end

draw_full("png")
draw_full("svg")
preview()

function draw_square(type)
    currentwidth = 250 # pts
    currentheight = 250 # pts
    Drawing(currentwidth, currentheight, "logo_square.$type")

    sethue(julia_red)
    circle(Point(85, 30), 25, :fill)

    sethue(julia_purple)
    circle(Point(160, 30), 25, :fill)

    # origin()

    fontsize(300)
    fontface("Literata SemiBold")

    str = "a" # Uses "dotless i"

    translate(45, 225)

    textpath(str)                             # make text into a path
    setline(3)
    setcolor("black")
    fillpreserve()                            # fill but keep

    finish()
end

draw_square("png")
draw_square("svg")
preview()

max_stars_repo_name: alecloudenback/ActuarialScienceLogo

max_issues_repo_name: alecloudenback/ActuarialScienceLogo

max_forks_repo_name: alecloudenback/ActuarialScienceLogo


text:
# some utility functions
const Γ = gamma

function generalized_binomial(α::T, n::S) where {T,S <: Integer}
    R = Base.promote_op(/,T,S)
    out = one(R)
    for i in n:-1:1
        out *= α/i
        α -= 1
    end
    out
end

## memoized version of the above
@memoize function generalized_binomial′(a, k::Int)
    k < 0 && return NaN * one(a)
    k == 0 && return one(a)/1
    generalized_binomial′(a,k-1) * (a-k+1)/k
end

## x⁽ⁱ⁾ is falling (x)⋅(x-1)⋯(x-n+1)
function Pochhammer(::Val{:falling}, z::T,n::Int) where {T}
    iszero(n) && return one(T)
    out = z
    for i in 1:n-1
        z -= 1
        out *= z
    end
    out
end

## (x)ᵢ is rising (x)⋅(x+1)⋯(x+n-1) = Γ(z+n)/Γ(z)
function Pochhammer(::Val{:rising}, z::T,n::Int) where {T}
    iszero(n) && return one(T)
    out = z
    for i in 1:n-1
        z += 1
        out *= z
    end
    out
end
# default is (x)ᵢ
Pochhammer(z,n::Int) = Pochhammer(Val(:rising),  z, n)

# (x)_n/n!
Pochhammer_factorial(z, n) = Pochhammer_factorial(Val(:rising), z, n)

function Pochhammer_factorial(::Val{:rising}, z, n)
    iszero(n) && return one(z)/1
    prod((z+i)/(n-i) for i in 0:n-1)
end

## Alternative to HypergeometricFunctions.jl so that
## variables can be inserted into the numerator
"""
    pFq(as, bs, z; [maxevals=1000])

Compute the generalized [hypergeometric](https://en.wikipedia.org/wiki/Generalized_hypergeometric_function) function. The `HypergeometricFunction.jl` package would normally be used for such calculations, but this one accepts polynomial values for the `as` and `z` values.

# Example

From [mathworld](https://mathworld.wolfram.com/HypergeometricFunction.html)

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> import SpecialPolynomials: pFq, Pochhammer

julia> pFq((1/3,2/3), 5/6, 27/32) ≈ 8/5
true

julia> pFq([1/4, 1/2], [3/4], 80/81; maxevals=2000) ≈ 9/5
true

julia> x = variable()
Polynomial(x)

julia> n = 5
5

julia> pFq((-n,n+1), 1, (1-x)/2) ≈ basis(Legendre,n)(x)
true

julia> α, β, n = 1/2, 1/2, 5;

julia> Pochhammer(α+1,n)/factorial(n) * pFq((-n,n+α+β+1), α+1, (1-x)/2) ≈ basis(Jacobi{α,β}, n)(x)
true
```
"""
function pFq(as, bs, z; maxevals=1000)
    n = 1
    acc = trm = one(z)

    p,q = length(as), length(bs)

    # element in `as` is negative integer or 0 && converges
    # p ≤ q && converges
    # p = q + 1 |x| < 1 && converge
    # p > q + 1 && diverges
    while n < maxevals
        a = isempty(as) ? 1 : prod(as)
        b = isempty(bs) ? 1 : prod(bs)

        iszero(a) && return acc
        iszero(b) && return Inf  # check on b

        trm *= a /b * z  /n
        acc += trm

        as = plus_1(as)

        bs = plus_1(bs)
        n += 1
    end

    if (p > q + 1 || (p == q+1 && abs(z) < 1))
        return acc
    else
        NaN*acc
    end
end
## tuples
@inline plus_1(as) = map(x->x+1, as)

# use HypergeometricFunction.mFn if possible
function pFq(as::Tuple{A,B}, bs::Tuple{C}, z::AbstractFloat) where {A,B,C}
    HypergeometricFunctions._₂F₁(as[1],as[2],bs[1],z)
end
function pFq(as::Tuple{A,B,C}, bs::Tuple{D,E}, z::AbstractFloat) where {A,B,C,D,E}
    HypergeometricFunctions._₃F₂(as[1],as[2],as[3],bs[1],bs[2],z)
end

## Try to speed up  quadgk by not specializing on F
mutable struct Wrapper
    F
end
(F::Wrapper)(x) = F.F(x)

_quadgk(f, a, b) = quadgk(Wrapper(f), a, b)[1]
const ∫ = _quadgk

checked_div(a, b) = (iszero(a) && iszero(b)) ? a : a/b


text:
mutable struct Prop3DAcoTTIDenQ_DEO2_FDTD
    p::Ptr{Cvoid}
end

function Prop3DAcoTTIDenQ_DEO2_FDTD(;
        nz=0,
        ny=0,
        nx=0,
        nsponge=60,
        nbz=512,
        nby=8,
        nbx=8,
        dz=5.0,
        dy=8.0,
        dx=5.0,
        dt=0.001,
        nthreads=Sys.CPU_THREADS,
        freesurface=true,
        freqQ=5.0,
        qMin=0.1,
        qInterior=100.0)
    nz,ny,nx,nsponge,nbz,nby,nbx,nthreads = map(x->round(Int,x), (nz,ny,nx,nsponge,nbz,nby,nbx,nthreads))
    dz,dy,dx,dt = map(x->Float32(x), (dz,dy,dx,dt))

    @assert nx > 0
    @assert ny > 0
    @assert nz > 0
    @assert nsponge > 0
    @assert nthreads > 0
    @assert nbx > 0
    @assert nby > 0
    @assert nbz > 0

    fs = freesurface ? 1 : 0

    p = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_alloc, libprop3DAcoTTIDenQ_DEO2_FDTD),
        Ptr{Cvoid},
        (Clong, Clong,    Clong, Clong, Clong, Clong,   Cfloat, Cfloat, Cfloat, Cfloat, Clong, Clong, Clong),
         fs,    nthreads, nx,    ny,    nz,    nsponge, dx,     dy,     dz,     dt,     nbx,   nby,   nbz)

    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_SetupDtOmegaInvQ, libprop3DAcoTTIDenQ_DEO2_FDTD),
        Cvoid,
        (Ptr{Cvoid}, Cfloat, Cfloat, Cfloat),
         p,          freqQ,  qMin,   qInterior)

    Prop3DAcoTTIDenQ_DEO2_FDTD(p)
end

free(prop::Prop3DAcoTTIDenQ_DEO2_FDTD) = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_free, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

function size(prop::Prop3DAcoTTIDenQ_DEO2_FDTD)
    nz = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNz, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    ny = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNy, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nx = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNx, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    (nz,ny,nx)
end

for _f in (:V, :Eps, :Eta, :B, :F, :PSpace, :MSpace, :PCur, :POld, :MCur, :MOld, :SinTheta, :CosTheta, :SinPhi, :CosPhi, :DtOmegaInvQ)
        symf = "Prop3DAcoTTIDenQ_DEO2_FDTD_get" * string(_f)
    @eval $(_f)(prop::Prop3DAcoTTIDenQ_DEO2_FDTD) = unsafe_wrap(Array, ccall(($symf, libprop3DAcoTTIDenQ_DEO2_FDTD), Ptr{Float32}, (Ptr{Cvoid},), prop.p), size(prop), own=false)
end

propagateforward!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD) = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_TimeStep, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)
propagateadjoint!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD) = propagateforward!(prop) # self-adjoint

scale_spatial_derivatives!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD) =
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_ScaleSpatialDerivatives, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

abstract type Prop3DAcoTTIDenQ_DEO2_FDTD_Model end
struct Prop3DAcoTTIDenQ_DEO2_FDTD_Model_VEA <: Prop3DAcoTTIDenQ_DEO2_FDTD_Model end
struct Prop3DAcoTTIDenQ_DEO2_FDTD_Model_V <: Prop3DAcoTTIDenQ_DEO2_FDTD_Model end

# v,ϵ and η in model-space
function forwardBornInjection!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_VEA, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_ForwardBornInjection_VEA, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat},        Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["ϵ"], dmodel["η"], wavefields["pold"], wavefields["mold"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_VEA,
        imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_AdjointBornAccumulation_VEA, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat},        Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["ϵ"], dmodel["η"], wavefields["pold"], wavefields["mold"], wavefields["pspace"], wavefields["mspace"])
end

# v in model-space
function forwardBornInjection!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_V, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_ForwardBornInjection_V, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_AdjointBornAccumulation_V, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"])
end

function adjointBornAccumulation!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionWaveFieldSeparationFWI, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep_V, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat},          Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"], 1)
end

function adjointBornAccumulation!(prop::Prop3DAcoTTIDenQ_DEO2_FDTD, modeltype::Prop3DAcoTTIDenQ_DEO2_FDTD_Model_V,
        imagingcondition::ImagingConditionWaveFieldSeparationRTM, dmodel, wavefields)
    ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep_V, libprop3DAcoTTIDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Ptr{Cfloat},          Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], wavefields["mspace"], 0)
end

function show(io::IO, prop::Prop3DAcoTTIDenQ_DEO2_FDTD)
    nx = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNx, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    ny = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNy, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nz = ccall((:Prop3DAcoTTIDenQ_DEO2_FDTD_getNz, libprop3DAcoTTIDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    write(io, "Prop3DAcoTTIDenQ_DEO2_FDTD -- nx,ny,nz; $nx,$ny,$nz")
end

function tilesize(::Type{Prop3DAcoTTIDenQ_DEO2_FDTD}; nz, ny, nx, rngy=4:2:64, rngx=4:2:64, N=10, nthreads=Sys.CPU_THREADS-1)
    speeds = Matrix{Float64}(undef, length(rngy), length(rngx))
    for (iby,by) in enumerate(rngy), (ibx,bx) in enumerate(rngx)
        @info "by=$by, bx=$bx"
        prop = Prop3DAcoTTIDenQ_DEO2_FDTD(nz=nz, ny=ny, nx=nx, nbz=nz, nby=by, nbx=bx, nthreads=nthreads)

        v = V(prop)
        ϵ = Eps(prop)
        η = Eta(prop)
        sinθ = SinTheta(prop)
        cosθ = CosTheta(prop)
        sinϕ = SinPhi(prop)
        cosϕ = CosPhi(prop)
        b = B(prop)
        f = F(prop)
        pcur = PCur(prop)
        pold = POld(prop)
        mcur = MCur(prop)
        mold = MOld(prop)

        v .= 1500
        ϵ .= 1
        η .= 1
        sinθ .= 0
        cosθ .= 1
        sinϕ .= 0
        cosϕ .= 1
        b .= 1
        f .= 1
        rand!(pcur)
        rand!(pold)
        rand!(mcur)
        rand!(mold)

        t = 0.0
        set_zero_subnormals(true)
        for i = 1:N
            t += @elapsed propagateforward!(prop)
        end
        set_zero_subnormals(false)
        speeds[iby,ibx] = N*prod(size(prop))/1000/1000/t
        free(prop)
    end
    rngx,rngy,speeds
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
module LocalCoverage

using Coverage
using DocStringExtensions
import Pkg

export generate_coverage, open_coverage, clean_coverage

"Directory for coverage results."
const COVDIR = "coverage"

"Coverage tracefile."
const LCOVINFO = "lcov.info"

"""
$(SIGNATURES)

Open the HTML coverage results in a browser for `pkg` if they exist.

See [`generate_coverage`](@ref).
"""
function open_coverage(pkg)
    htmlfile = Pkg.dir(pkg, COVDIR, "index.html")
    if !isfile(htmlfile)
        @warn("Not found, run generate_coverage(pkg) first.")
        return nothing
    end
    try
        if Sys.isapple()
            run(`open $htmlfile`)
        elseif Sys.islinux() || Sys.isbsd()
            run(`xdg-open $htmlfile`)
        elseif Sys.iswindows()
            run(`start $htmlfile`)
        end
    catch e
        error("Failed to open the generated $(htmlfile)\n",
              "Error: ", sprint(Base.showerror, e))
    end
    nothing
end

"""
$(SIGNATURES)

Generate a coverage report for package `pkg`.

When `genhtml`, the corresponding external command will be called to generate a
HTML report. This can be found in eg the package `lcov` on Debian/Ubuntu.

`*.cov` files are near the source files as generated by Julia, everything else
is placed in `Pkg.dir(pkg, \"$(COVDIR)\")`. The summary is in
`Pkg.dir(pkg, \"$(COVDIR)\", \"$(LCOVINFO)\")`.

Use [`clean_coverage`](@ref) for cleaning.
"""
function generate_coverage(pkg; genhtml = true)
    Pkg.test(pkg; coverage = true)
    cd(Pkg.dir(pkg)) do
        coverage = Coverage.process_folder()
        isdir(COVDIR) || mkdir(COVDIR)
        tracefile = "$(COVDIR)/lcov.info"
        Coverage.LCOV.writefile(tracefile, coverage)
        if genhtml
            branch = strip(read(`git rev-parse --abbrev-ref HEAD`, String))
            title = "on branch $(branch)"
            run(`genhtml -t $(title) -o $(COVDIR) $(tracefile)`)
        end
    end
end

"""
$(SIGNATURES)

Clean up after [`generate_coverage`](@ref).
"""
function clean_coverage(pkg)
    Coverage.clean_folder(Pkg.dir(pkg))
    rm(Pkg.dir(pkg, COVDIR); force = true, recursive = true)
end

end # module


text:
# PS2 CHEME 5440/7770 Cornell University
#
# instructions:
# To execute this script execute the command (from the Julia REPL):
# julia> include("execute_me_PS2.jl")

# load the include file -
include("Include.jl")

# main method -
function main(parameters_file_path::String, inducer_array::Array{Float64,1})::Array{Float64,2}

    # initialize -
    number_of_inducer_levels = length(inducer_array)
    mRNA_steady_state_array = Array{Float64,2}(undef,number_of_inducer_levels,2)

    try

        # load the parameter file, and create the problem dictionary -
        problem_dictionary = generate_problem_dictionary(parameters_file_path)

        # For PS2, we are calculating the steady state level of mRNA as a function of inducer -
        initial_condition_array = problem_dictionary["initial_condition_array"]
        for (index, inducer_level) in enumerate(inducer_array)

            # set the inducer concentration -
            initial_condition_array[3] = inducer_level
            problem_dictionary["initial_condition_array"] = initial_condition_array

            # estimate the steady state -
            XSS = estimate_steady_state(problem_dictionary)

            # capture the steady state mRNA -
            mRNA_steady_state_array[index,1] = inducer_level
            mRNA_steady_state_array[index,2] = XSS[1]
        end

        # return -
        return mRNA_steady_state_array
    catch error

        # grab the error message, and post -
        @error "Oooops! Hmmmmm. Something happend" exception=(error, catch_backtrace())
    end
end

# setup the path to the parameters file -
path_to_parameters_file = joinpath(_PATH_TO_CONFIG, "Parameters.toml")

# setup the inducer array: units of nM
inducer_array = 10 .^(range(-3,stop=3, length=100))

# execute -
mRNA_steady_state_array = main(path_to_parameters_file, inducer_array)

# to plot - uncomment me -
plot(mRNA_steady_state_array[:,1],mRNA_steady_state_array[:,2], xscale=:log10, legend=false)
xlabel!("Inducer concentration (nM)")
ylabel!("Steady state mRNA concentration (nM)")

path_to_plot = joinpath(_PATH_TO_OUTPUT,"PS2_Plot.png")
savefig(path_to_plot)


text:
module CliffsideFlag

using ..Ahorn, Maple

validFlagIndicdes = collect(0:10)

const placements = Ahorn.PlacementDict(
    "Cliffside Big Flag" => Ahorn.EntityPlacement(
        Maple.CliffsideFlag,
    )
)

Ahorn.editingOptions(entity::Maple.CliffsideFlag) = Dict{String, Any}(
    "index" => validFlagIndicdes
)

function flagSprite(entity::Maple.CliffsideFlag)
    index = Int(get(entity.data, "index", 0))
    lookup = lpad(string(index), 2, "0")

    return "scenery/cliffside/flag$(lookup)"
end

function Ahorn.selection(entity::Maple.CliffsideFlag)
    x, y = Ahorn.position(entity)
    sprite = flagSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.0, jy=0.0)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CliffsideFlag, room::Maple.Room)
    sprite = flagSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.0, jy=0.0)
end

end


text:
module CrumbleWallOnRumble

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Crumble Wall on Rumble" => Ahorn.EntityPlacement(
        Maple.CrumbleWallOnRumble,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    )
)

Ahorn.editingOptions(entity::Maple.CrumbleWallOnRumble) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::Maple.CrumbleWallOnRumble) = 8, 8
Ahorn.resizable(entity::Maple.CrumbleWallOnRumble) = true, true

Ahorn.selection(entity::Maple.CrumbleWallOnRumble) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CrumbleWallOnRumble, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end


text:
import PredictMD
import Test

tmpdir = PredictMD.maketempdir()
tmpfile = joinpath(tmpdir, "Project.toml")
Test.@test_throws(ErrorException, PredictMD.TomlFile(tmpfile))


text:
using Pkg
Pkg.activate("./metrics")
Pkg.instantiate()

using CSV
using DataFrames
using Distributed
using Dates

metric_target_raw = get(ENV, "METRIC_TARGET", "")
if isempty(metric_target_raw)
    error("No metric target set.")
end

run(`git fetch`)

metric_targets = if startswith(metric_target_raw, "regex ")
    # a workflow_dispatch input of "regex some_regex" matches tags against some_regex
    tags = strip.(split(read(`git tag`, String)))
    regex = match(r"regex (.*)", metric_target_raw)[1] |> strip
    r = Regex(regex)
    ts = filter(x -> !isnothing(match(r, x)), tags)
    @info "available tags: $tags"
    if isempty(ts)
        error("No metric targets found with regex: $r")
    end
    ts
else
    [metric_target_raw]
end

@info "metric targets: $metric_targets"

results = DataFrame()

# one vector of code parts for every file in metrics folder
parts_vectors = map(readdir(joinpath("metrics", "metrics"), join = true)) do file
    code = read(file, String)
    parts = split(code, r"^(?=## )"m, keepempty = false)
end

for metric_target in metric_targets
    @info "checking out metric target $metric_target"
    run(`git checkout $metric_target --`)

    makieversion = match(r"version = \"(.*?)\"", read("Project.toml", String))[1]
    glmakieversion = match(r"version = \"(.*?)\"", read("GLMakie/Project.toml", String))[1]
    cairomakieversion = match(r"version = \"(.*?)\"", read("CairoMakie/Project.toml", String))[1]
    commit_date = DateTime(
        strip(String(read(`git show -s --format=%ci`)))[1:end-6],
        "yyyy-mm-dd HH:MM:SS")

    df = DataFrame()
    date = now()

    for parts in parts_vectors

        local i_proc
        try
            i_proc = addprocs(1)[1]

            @everywhere i_proc begin
                using Pkg
            end

            @everywhere i_proc begin
                pkg"activate --temp"
                pkg"dev . MakieCore GLMakie CairoMakie"
                Pkg.precompile()
                @timed begin end
            end

            for part in parts
                partname = match(r"^## (.*)", part)[1] |> strip
                @info "executing part: $partname"
                partcode = """
                    @timed begin
                        $part
                        nothing
                    end
                """
                timing = remotecall_fetch(i_proc, partcode) do p
                    include_string(Main, p)
                end

                push!(df, (
                    date = date,
                    commit_date = commit_date,
                    metric_target = metric_target,
                    juliaversion = string(Sys.VERSION),
                    makie = makieversion,
                    glmakie = glmakieversion,
                    cairomakie = cairomakieversion,
                    name = partname,
                    time = timing.time,
                    allocations = timing.bytes,
                    gc_time = timing.gctime,
                ))
            end

        finally
            rmprocs(i_proc)
        end
    end
    append!(results, df, cols = :union)
end

branch_name = "metrics"
remote_branch = "origin/metrics"

run(`git config --global user.email "41898282+github-actions[bot]@users.noreply.github.com"`)
run(`git config --global user.name "github-actions[bot]"`)

@info "Checking out $branch_name."
if !success(`git checkout -b $branch_name $remote_branch`)
    @info "branch $branch_name doesn't exist, creating new orphan branch"
    run(`git checkout --orphan $branch_name --`)
    run(`git rm -rf .`)
end

println("\nFiles in working directory:")
run(`ls`)

filename = "compilation_latencies.csv"

df = if !isfile(filename)
    @info "$filename doesn't exist, creating empty DataFrame."
    DataFrame()
else
    @info "Loading DataFrame from $filename."
    CSV.read(filename, DataFrame)
end

append!(df, results, cols = :union)
sort!(df, :commit_date)

@info "Writing out DataFrame to $filename."
CSV.write(filename, df)

run(`git add $filename`)
run(`git commit -m "update metrics"`)
run(`git push origin $branch_name`)


text:
using SparsityOperators
using Test
using Random
using LinearAlgebra
using FFTW

include("testOperators.jl")


text:
"""
    abstract type AbstractMesh{N,T}

An abstract mesh structure in dimension `N` with primite data of type `T`. By
default `T=Float64`, and `N=3`.

The `AbstractMesh` interface expects the following methods to be implemented:

- `etypes(msh)` : return a `Vector{AbstractElement}` representing the element
  types composing the mesh.
- `elements(msh,etype)` : return an iterator over the mesh elements of type `etype::AbstractElement`
- `msh[E][i]` : return the *i-th* element of type `E`.
"""
abstract type AbstractMesh{N,T} end

ambient_dimension(M::AbstractMesh{N}) where {N} = N

# we define the geometric dimension of the mesh to be the largest of the geometric
# dimension of its entities.
geometric_dimension(M::AbstractMesh) = maximum(x -> geometric_dimension(x), etypes(M))

Base.eltype(M::AbstractMesh{N,T}) where {N,T} = T

Base.length(mesh::AbstractMesh) = length(nodes(mesh))

"""
    struct GenericMesh{N,T} <: AbstractMesh{N,T}

Data structure representing a generic mesh in an ambient space of dimension `N`, with data of type `T`.
"""
Base.@kwdef struct GenericMesh{N,T} <: AbstractMesh{N,T}
    nodes::Vector{SVector{N,T}} = Vector{SVector{N,T}}()
    # element types
    etypes::Vector{DataType} = Vector{DataType}()
    # for each element type (key), get the data required to reconstruct the
    # elements (value)
    elements::OrderedDict{DataType,Any} = OrderedDict{DataType,Any}()
    # mapping from elementary entity to (etype,tags)
    ent2tags::OrderedDict{AbstractEntity,OrderedDict{DataType,Vector{Int}}} = OrderedDict{AbstractEntity,OrderedDict{DataType,Vector{Int}}}()
end

# convert a mesh to 2d by ignoring third component. Note that this also requires
# converting various element types to their 2d counterpart.
function convert_to_2d(mesh::GenericMesh{3})
    @assert all(x -> geometric_dimension(x) < 3, etypes(mesh))
    T = eltype(mesh)
    # create new dictionaries for elements and ent2tags with 2d elements as keys
    elements  = empty(mesh.elements)
    ent2tags  = empty(mesh.ent2tags)
    for (E, tags) in mesh.elements
        E2d = convert_to_2d(E)
        elements[E2d] = tags
    end
    for (ent, dict) in mesh.ent2tags
        new_dict = empty(dict)
        for (E, tags) in dict
            E2d = convert_to_2d(E)
            new_dict[E2d] = tags
        end
        ent2tags[ent] = new_dict
    end
    # construct new 2d mesh
    GenericMesh{2,T}(;
        nodes=[x[1:2] for x in nodes(mesh)],
        etypes=convert_to_2d.(etypes(mesh)),
        elements=elements,
        ent2tags=ent2tags
    )
end

nodes(m::GenericMesh)    = m.nodes
elements(m::GenericMesh) = m.elements
ent2tags(m::GenericMesh) = m.ent2tags

function convert_to_2d(::Type{LagrangeElement{R,N,SVector{3,T}}}) where {R,N,T}
    LagrangeElement{R,N,SVector{2,T}}
end
convert_to_2d(::Type{SVector{3,T}}) where {T} = SVector{2,T}

"""
    etypes(M::GenericMesh)

Return the element types contained in the mesh.
"""
etypes(mesh::GenericMesh) = mesh.etypes

"""
    dom2elt(m::GenericMesh,Ω,E)

Compute the element indices `idxs` of the elements of type `E` composing `Ω`, so that `m[E][idxs]` gives all
the elements of type `E` meshing `Ω`.
"""
function dom2elt(m::GenericMesh,Ω,E::DataType)
    idxs = Int[]
    for ent in entities(Ω)
        tags = get(m.ent2tags[ent],E,Int[])
        append!(idxs,tags)
    end
    return idxs
end

"""
    dom2elt(m::GenericMesh,Ω)

Return a `OrderedDict` with keys being the `etypes` of `m`, and values being the
indices of the elements in `Ω` of type `E`.
"""
function dom2elt(m::GenericMesh,Ω)
    dict = OrderedDict{DataType,Vector{Int}}()
    for E in etypes(m)
        tags = dom2elt(m,Ω,E)
        if !isempty(tags)
            dict[E] = tags
        end
    end
    return dict
end

"""
    _qrule_for_mesh(m,p)

Given a mesh `m`, create a dictionary mapping each element type of `m` to an
appropriate quadrature rule of order `p` over that element type.

See also [`Integration._qrule_for_reference_shape`](@ref)
"""
function _qrule_for_mesh(m,order)
    OrderedDict(E=>Integration._qrule_for_element(E,order) for E in etypes(m))
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# ----------------------------------------------------
# --- Import modules
# ----------------------------------------------------
using DigitalComm
using Test
# ----------------------------------------------------
# --- Tests
# ----------------------------------------------------
# Note --> The mapping system is described in bitMapping.jl

println("Tests for Max-Log soft  symbol demapper");
for alphabet = [4;16;64;256]
	@testset  "$alphabet-QAM" begin
		# Minimal Tx
		mcs	  = alphabet;		  # Constellation size
		n	  = Int(log2(mcs));		  # Bit per symbol
		N	  = 4096;	  # Numbver of test symbols
		bitSeq	= genBitSequence(N*n);
		qamSeq	= bitMappingQAM(mcs,bitSeq);
		# Assuming unitary channel
		channel = ones(Complex{Float64},length(qamSeq));
		# Calculating LLR
		llr		= symbolDemappingQAM(mcs,qamSeq,channel);
		hardB	= llrToHardBits(llr);
		# Checking size of decoding vector
		@test length(llr)	== N*n;
		@test length(hardB) == N*n;
		# Checking that bang and classic methods are the same
		buff  = zeros(Float64,N*n);
		symbolDemappingQAM!(buff,mcs,qamSeq,channel);
		buffB = zeros(UInt8,n*N);
		llrToHardBits!(buffB,buff);
		@test all(hardB .== buffB);
		# Checking that decoded method is initial method
		@test all( buffB .≈ bitSeq)
	end
end


text:
module ConvexHull

struct Point
    x::Float64
    y::Float64
end

function Base.isless(p::Point, q::Point)
    p.x < q.x || (p.x == q.x && p.y < q.y)
end

function isrightturn(p::Point, q::Point, r::Point)
    (q.x - p.x) * (r.y - p.y) - (q.y - p.y) * (r.x - p.x) < 0
end

function grahamscan(points::Vector{Point})
    sorted = sort(points)
    upperhull = halfhull(sorted)
    lowerhull = halfhull(reverse(sorted))
    [upperhull..., lowerhull[2:end-1]...]
end

function halfhull(points::Vector{Point})
    halfhull = points[1:2]
    for p in points[3:end]
        push!(halfhull, p)
        while length(halfhull) > 2 && !isrightturn(halfhull[end-2:end]...)
            deleteat!(halfhull, length(halfhull) - 1)
        end
    end
    halfhull
end

end # module


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

struct DensityWithShape{D<:AbstractDensity,S<:AbstractValueShape} <: AbstractDensity
    density::D
    shape::S
end

Base.parent(density::DensityWithShape) = density.density

var_bounds(density::DensityWithShape) = var_bounds(parent(density))

ValueShapes.varshape(density::DensityWithShape) = density.shape

DensityInterface.logdensityof(density::DensityWithShape, v::Any) = logdensityof(parent(density), v)

struct ReshapedDensity{D<:AbstractDensity,S<:AbstractValueShape} <: AbstractDensity
    density::D
    shape::S
end

Base.parent(density::ReshapedDensity) = density.density

var_bounds(density::ReshapedDensity) = var_bounds(density.density)

ValueShapes.varshape(density::ReshapedDensity) = density.shape

function _reshapeddensity_variate(density::ReshapedDensity, v::Any)
    vs = varshape(density)
    orig_vs = varshape(parent(density))
    reshape_variate(orig_vs, vs, v)
end

function DensityInterface.logdensityof(density::ReshapedDensity, v::Any)
    v_reshaped = _reshapeddensity_variate(density, v)
    logdensityof(parent(density), v_reshaped)
end

function checked_logdensityof(density::ReshapedDensity, v::Any)
    v_reshaped = _reshapeddensity_variate(density, v)
    checked_logdensityof(parent(density), v_reshaped)
end

ValueShapes.unshaped(density::AbstractDensity) = _unshaped_density(density, varshape(density))

_unshaped_density(density::AbstractDensity, ::ArrayShape{<:Real,1}) = density
_unshaped_density(density::AbstractDensity, ::AbstractValueShape) = ReshapedDensity(density, ArrayShape{Real}(totalndof(density)))

(shape::AbstractValueShape)(density::AbstractDensity) = _reshaped_density(density, shape, varshape(density))

function _reshaped_density(density::AbstractDensity, new_shape::VS, orig_shape::VS) where {VS<:AbstractValueShape}
    if orig_shape == new_shape
        density
    else
        throw(ArgumentError("Value shapes are incompatible"))
    end
end

function _reshaped_density(density::AbstractDensity, new_shape::AbstractValueShape, orig_shape::AbstractValueShape)
    if totalndof(orig_shape) == totalndof(new_shape)
        ReshapedDensity(density, new_shape)
    else
        throw(ArgumentError("Value shapes are incompatible"))
    end
end


text:
using DSP
using Test
using DelimitedFiles: readdlm

@testset "rational ratio" begin
    # AM Modulator
    # sig(t) = [(1 + sin(2π*0.005*t)) * sin(2π*.05*t) for t in t]
    x_ml   = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_x.txt"),'\t'))

    #
    # [y1,b1] = resample(x, 1, 2)
    #
    rate   = 1//2
    h1_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_1_2.txt"),'\t'))
    y1_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_1_2.txt"),'\t'))
    y1_jl  = resample(x_ml, rate, h1_ml)
    @test y1_jl ≈ y1_ml

    #
    # [y2,b2] = resample(x, 2, 1)
    #
    rate   = 2//1
    h2_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_2_1.txt"),'\t'))
    y2_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_2_1.txt"),'\t'))
    y2_jl  = resample(x_ml, rate, h2_ml)
    @test y2_jl ≈ y2_ml

    #
    # [y3,b3] = resample(x, 3, 2)
    #
    rate   = 3//2
    h3_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_3_2.txt"),'\t'))
    y3_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_3_2.txt"),'\t'))
    y3_jl  = resample(x_ml, rate, h3_ml)
    @test y3_jl ≈ y3_ml

    #
    # [y4,b4] = resample(x, 2, 3)
    #
    rate  = 2//3
    h4_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_2_3.txt"),'\t'))
    y4_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_2_3.txt"),'\t'))
    y4_jl = resample(x_ml, rate, h4_ml)
    @test y4_jl ≈ y4_ml
end

@testset "array signal" begin
    rate   = 1//2
    x_ml   = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_x.txt"),'\t'))
    h1_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_1_2.txt"),'\t'))
    y1_ml  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_1_2.txt"),'\t'))

    expected_result = [y1_ml 2y1_ml]
    X = [x_ml 2x_ml]
    y1_jl  = resample(X, rate, h1_ml, dims=1)
    @test y1_jl ≈ expected_result

    y1_jl  = resample(X', rate, h1_ml, dims=2)
    @test y1_jl ≈ expected_result'

    expected_result_3d = permutedims(reshape(expected_result, (size(expected_result, 1), size(expected_result, 2), 1)), (3, 1, 2))
    X_3d = permutedims(reshape(X, (size(X, 1), size(X, 2), 1)), (3, 1, 2))
    y1_jl  = resample(X_3d, rate, h1_ml, dims=2)
    @test y1_jl ≈ expected_result_3d

    expected_result_3d = permutedims(expected_result_3d, (1, 3, 2))
    X_3d = permutedims(X_3d, (1, 3, 2))
    y1_jl  = resample(X_3d, rate, h1_ml, dims=3)
    @test y1_jl ≈ expected_result_3d
end

@testset "irrational ratio" begin
    ratio    = 3.141592653589793
    cycles   = 2
    tx       = range(0, stop=cycles, length=1000)
    x        = sinpi.(2*tx)
    y        = resample(x, ratio)
    yLen     = length(y)
    ty       = range(0, stop=cycles, length=yLen)
    yy       = sinpi.(2*ty)
    idxLower = round(Int, yLen/3)
    idxUpper = idxLower*2
    yDelta   = abs.(y[idxLower:idxUpper].-yy[idxLower:idxUpper])
    @test all(map(delta -> abs(delta) < 0.005, yDelta))

    # Test Float32 ratio (#302)
    f32_ratio = convert(Float32, ratio)
    f32_y     = resample(x, f32_ratio)
    ty        = range(0, stop =cycles, length =yLen)
    yy        = sinpi.(2*ty)
    idxLower  = round(Int, yLen/3)
    idxUpper  = idxLower*2
    yDelta    = abs.(f32_y[idxLower:idxUpper].-yy[idxLower:idxUpper])
    @test all(map(delta -> abs(delta) < 0.005, yDelta))
end

@testset "resample_filter" begin
    @testset "decimation" begin
        ratio = 1//2
        h     = resample_filter(ratio)
        r0    = abs.(freqresp(PolynomialRatio(h, [1]), 0))
        rc    = abs.(freqresp(PolynomialRatio(h, [1]), ratio*π))
        @test isapprox(r0, 1.0)
        @test isapprox(rc, numerator(ratio)/2, rtol=0.001)

        ratio = 1//32
        h     = resample_filter(ratio)
        r0    = abs.(freqresp(PolynomialRatio(h, [1]), 0))
        rc    = abs.(freqresp(PolynomialRatio(h, [1]), ratio*π))
        @test isapprox(r0, numerator(ratio))
        @test isapprox(rc, numerator(ratio)/2, rtol=0.001)
    end

    @testset "interpolation" begin
        ratio = 2//1
        h     = resample_filter(ratio)
        r0    = abs.(freqresp(PolynomialRatio(h, [1]), 0))
        rc    = abs.(freqresp(PolynomialRatio(h, [1]), 1/ratio*π))
        @test isapprox(r0, numerator(ratio))
        @test isapprox(rc, numerator(ratio)/2, rtol=0.001)

        ratio = 32//1
        h     = resample_filter(ratio)
        r0    = abs.(freqresp(PolynomialRatio(h, [1]), 0))
        rc    = abs.(freqresp(PolynomialRatio(h, [1]), 1/ratio*π))
        @test isapprox(r0, numerator(ratio))
        @test isapprox(rc, numerator(ratio)/2, rtol=0.001)
    end

    @testset "arbitrary rate" begin
        ratio = 3.141592653589793
        Nϕ    = 32
        fc    = 1/Nϕ
        h     = resample_filter(ratio, Nϕ)
        r0    = abs.(freqresp(PolynomialRatio(h, [1]), 0))
        rc    = abs.(freqresp(PolynomialRatio(h, [1]), fc*π))
        @test isapprox(r0, Nϕ)
        @test isapprox(rc, Nϕ/2, rtol=0.001)
    end
end


text:
type PiecewiseDefaultCurve{B <: BootstrapHelper, DC <: DayCount, P <: Interpolation, T <: BootstrapTrait, BT <: Bootstrap} <: InterpolatedDefaultProbabilityCurve{P}
  lazyMixin::LazyMixin
  settlementDays::Int
  referenceDate::Date
  instruments::Vector{B}
  dc::DC
  interp::P
  trait::T
  accuracy::Float64
  boot::BT
  times::Vector{Float64}
  dates::Vector{Date}
  data::Vector{Float64}
  errors::Vector{Function}
  validCurve::Bool
end

function PiecewiseDefaultCurve{B <: BootstrapHelper, DC <: DayCount, P <: Interpolation, T <: BootstrapTrait, BT <: Bootstrap}(referenceDate::Date, instruments::Vector{B}, dc::DC, interp::P, trait::T,
                                accuracy::Float64, boot::BT = IterativeBootstrap())
  # get the initial length of instruments
  n = length(instruments)
  # create an initial state of the curve
  pyc = PiecewiseDefaultCurve{B, DC, P, T, BT}(LazyMixin(),
                            0,
                            referenceDate,
                            instruments,
                            dc,
                            interp,
                            trait,
                            accuracy,
                            boot,
                            Vector{Float64}(n + 1),
                            Vector{Date}(n + 1),
                            Vector{Float64}(n + 1),
                            Vector{Function}(n + 1),
                            false)

  # initialize the bootstrapping
  initialize(pyc.boot, pyc)

  return pyc
end

function perform_calculations!(curve::InterpolatedDefaultProbabilityCurve)
  _calculate!(curve.boot, curve)
  return curve
end

survival_probability(ts::AbstractDefaultProbabilityTermStructure, d::Date) = survival_probability(ts, time_from_reference(ts, d))

function survival_probability(ts::AbstractDefaultProbabilityTermStructure, t::Float64)
  # TODO handle jumps
  return survival_probability_impl(ts, t)
end

function survival_probability_impl(ts::PiecewiseDefaultCurve, t::Float64)
  calculate!(ts)
  if t == 0.0
    return 1.0
  end

  if t <= ts.times[end]
    integral = get_primitive(ts.interp, t, true)
  else
    # flat hazard rate extrapolation
    integral = get_primitive(ts.interp, ts.times[end], true) + ts.data[end] * (t - ts.times[end])
  end

  return exp(-integral)
end

default_probability(ts::AbstractDefaultProbabilityTermStructure, d::Date) = 1.0 - survival_probability(ts, d)
default_probability(ts::AbstractDefaultProbabilityTermStructure, t::Float64) = 1.0 - survival_probability(ts, t)

function default_probability(ts::AbstractDefaultProbabilityTermStructure, d1::Date, d2::Date)
  p1 = d1 < reference_date(ts) ? 0.0 : default_probability(ts, d1)
  p2 = default_probability(ts, d2)

  return p2 - p1
end

function default_probability(ts::AbstractDefaultProbabilityTermStructure, t1::Float64, t2::Float64)
  t1 <= t2 || error("time mismatch")

  p1 = t1 < 0.0 ? 0.0 : default_probability(ts, t1)
  p2 = default_probability(ts, t2)
  return p2 - p1
end

default_density(ts::AbstractDefaultProbabilityTermStructure, d::Date) = default_density(ts, time_from_reference(ts, d))
function default_density(ts::AbstractDefaultProbabilityTermStructure, t::Float64)
  # TODO check range
  return default_density_impl(ts, t)
end

default_density_impl(ts::AbstractDefaultProbabilityTermStructure, t::Float64) = hazard_rate_impl(ts, t) * survival_probability_impl(ts, t)
default_density_impl(ts::PiecewiseDefaultCurve, t::Float64) = _default_density_impl(ts, t) * survival_probability_impl(ts, t)

# sub method
function _default_density_impl(ts::PiecewiseDefaultCurve, t::Float64)
  calculate!(ts)
  if t <= ts.times[end]
    return QuantLib.Math.value(ts.interp, t)
  end

  return ts.data[end]
end

hazard_rate(ts::AbstractDefaultProbabilityTermStructure, d::Date) = hazard_rate(ts, time_from_reference(ts, d))

function hazard_rate(ts::AbstractDefaultProbabilityTermStructure, t::Float64)
  S = survival_probability(ts, t)
  return S == 0.0 ? 0.0 : default_density(ts, t) / S
end

function nodes(ts::PiecewiseDefaultCurve)
  calculate!(ts)
  results = Vector{Tuple}(length(ts.dates))
  for i in eachindex(ts.dates)
    results[i] = (ts.dates[i], ts.data[i])
  end

  return results
end


text:
using Printf
import ModiaResult

# the sol is different for modia than ModelingToolkit
function reportMinMax(name::String,modelName::String,result)

    (time, signal, signalType) = ModiaResult.rawSignal(result, name)

    y_min, i_min = findmin(signal[1])
    y_max, i_max = findmax(signal[1])
    t_min = time[1][i_min]
    t_max = time[1][i_max]
    print("$(name): $(modelName)\n")
    print("var\ttime\tvalue\n")
    print("min\t$(t_min)\t$(y_min)\n")
    print("max\t$(t_max)\t$(y_max)\n\n")
    return y_min, t_min, y_max, t_max
end

function errorCheck(measured, predicted)
    error_abs = measured-predicted
    error_rel = (measured-predicted)/measured

    println("Measured\tPredicted\tM-P\t\t(M-P)/M")
    @printf("%.2e\t%.2e\t%.2e\t%.2e\n\n",ustrip(measured),ustrip(predicted),ustrip(error_abs),ustrip(error_rel))
end

max_stars_repo_path: TRANSFORM - Modia/Utilities/summary.jl

max_issues_repo_path: TRANSFORM - Modia/Utilities/summary.jl

max_forks_repo_path: TRANSFORM - Modia/Utilities/summary.jl


text:
# quick script to draw Pose2D what is in the DB

using Caesar, RoME, Gadfly

## Uncomment out for command line operation
cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
session = addrdict["session"]
Nparticles = parse(Int,addrdict["num particles"])

# interactive operation
# session = "SESSROX"
# Nparticles = 100
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

# fieldnames(fg.cg.neo4j)

fg = initfg(sessionname=session, cloudgraph=cloudGraph)

# setBackendWorkingSet!(conn, session)
fullLocalGraphCopy!(fg)

println("Plotting...")
pl = drawPosesLandms(fg)

println("Rendering plotpose2_$(session).pdf...")
draw(PDF("plotpose2_$(session).pdf",40cm,40cm),pl)

#


text:
"""
    get_prop_from_block(block::Array, index::Integer)
Get the full property from `block` at `index`. Returns an array or scalar depending on the block type.
"""
function get_prop_from_block(block::Array, index::Integer)
    return ndims(block) == 1 ? block[index] : block[:, index]
end

"""
    read_halo_prop(filebase::String, haloid::HaloID, blockname::String; verbose::Bool=true)
Get halo property from block `blockname` by halo id. Returns an array or scalar depending on the block type.
"""
function read_halo_prop(filebase::String, haloid::HaloID, blockname::String; verbose::Bool=true)

    if verbose
        @info "Reading property $blockname of halo $haloid"
    end

    # get full specified block for all halos in file
    sub_input = select_file(filebase, haloid.file)
    block = read_subfind(sub_input, blockname)

    # return array or scalar for given halo
    return get_prop_from_block(block, haloid.id)
end

"""
    read_halo_prop_and_id(filebase::String, i_global::Integer, blockname::String, nfiles::Integer=1; verbose::Bool=true)
Get halo property from block `blockname` by global halo index `i_global` (zero-based index).
`nfiles` should generally be set to `h.num_files`, obtained from `read_header`.
"""
function read_halo_prop_and_id(filebase::String, i_global::Integer, blockname::String, nfiles::Integer=1; verbose::Bool=true)
    if verbose
        @info "Reading property $blockname of halo at index $i_global"
    end

    # store here for error handling
    i_global0 = i_global

    # run through files until halo at index i_global is reached
    for i in 0:nfiles - 1
        sub_input = select_file(filebase, i)

        if verbose
            @info "Reading file $i of $(nfiles - 1)"
        end

        block = read_subfind(sub_input, blockname)
        len = size(block, ndims(block))
        if len ≤ i_global # halo is not in current file
            i_global -= len
        else
            return get_prop_from_block(block, i_global + 1), HaloID(i, i_global + 1)
        end
    end
    error("Halo at index $i_global0 does not exist!")
end


text:
using Test, MyPackage

@test 42 == (40 + 2)

aShape = Box()
@test aShape.position == (0.0, 0.0)
@test getPosition(aShape) == (0.0, 0.0)
@test getPosition(aShape) == (0, 0)

setPosition!(aShape, (2.0, 3.0))
@test getPosition(aShape) == (2.0, 3.0)
@test getPosition(aShape) == (2, 3)

aCircle = Circle()
@test getSize(aCircle) == 5
setSize!(aCircle, 15.0)
@test getSize(aCircle) == 15
setSize!(aCircle, 24)
@test getSize(aCircle) == 24

composite = Composite()
@test numberOfShapes(composite) == 0
addShape!(composite, aCircle)
@test numberOfShapes(composite) == 1
addShape!(composite, aShape)
@test numberOfShapes(composite) == 2


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# In this notebook, we'll demonstrate an extremely simple approach for computing basic inverse kinematics (IK) and controlling the position of some point on our robot using the Jacobian transpose.
#
# For a brief technical introduction, see <https://groups.csail.mit.edu/drl/journal_club/papers/033005/buss-2004.pdf> or <https://homes.cs.washington.edu/~todorov/courses/cseP590/06_JacobianMethods.pdf>

# ## Setup

using Pkg # hide
Pkg.activate(@__DIR__) # hide
Pkg.instantiate() # hide
using RigidBodyDynamics
using StaticArrays
using MeshCat, MeshCatMechanisms, Blink

# Fix the random seed, so we get repeatable results
using Random
Random.seed!(42);

# First we'll load our double pendulum robot from URDF:

srcdir = dirname(pathof(RigidBodyDynamics))
urdf = joinpath(srcdir, "..", "test", "urdf", "Acrobot.urdf")
mechanism = parse_urdf(urdf)
state = MechanismState(mechanism)
mechanism

# Now we choose a point on the robot to control. We'll pick the end of the second link, which is located 2m from the origin of the `lower_link` body:

body = findbody(mechanism, "lower_link")
point = Point3D(default_frame(body), 0., 0, -2)

# Let's visualize the mechanism and its attached point. For visualization, we'll use [MeshCatMechanisms.jl](https://github.com/JuliaRobotics/MeshCatMechanisms.jl) with [Blink.jl](https://github.com/JunoLab/Blink.jl).

## Create the visualizer
vis = MechanismVisualizer(mechanism, URDFVisuals(urdf))

## Render our target point attached to the robot as a sphere with radius 0.07
setelement!(vis, point, 0.07)

# Open the visualizer in a new Blink window:
#nb ##NBSKIP
#nb open(vis, Window())
#md ## open(vis, Window());

# Create a `MechanismVisualizer`:

mvis = MechanismVisualizer(mechanism, URDFVisuals(urdf));

# ## Inverse Kinematics

# First, let's use the point jacobian to solve a simple inverse kinematics problem. Given a target location `desired` expressed in world frame, we want to find the joint angles `q` such that the `point` attached to the robot is at the desired location.
#
# To do that, we'll iteratively update `q` by applying:
#
# \begin{align}
# \Delta q = \alpha \, J_p^\top \, \Delta p
# \end{align}
#
# where $\alpha$ is our step size (equivalent to a learning rate in gradient descent) and $\Delta p$ is the error in the position of our target point.

function jacobian_transpose_ik!(state::MechanismState,
                               body::RigidBody,
                               point::Point3D,
                               desired::Point3D;
                               α=0.1,
                               iterations=100)
    mechanism = state.mechanism
    world = root_frame(mechanism)

    ## Compute the joint path from world to our target body
    p = path(mechanism, root_body(mechanism), body)
    ## Allocate the point jacobian (we'll update this in-place later)
    Jp = point_jacobian(state, p, transform(state, point, world))

    q = copy(configuration(state))

    for i in 1:iterations
        ## Update the position of the point
        point_in_world = transform(state, point, world)
        ## Update the point's jacobian
        point_jacobian!(Jp, state, p, point_in_world)
        ## Compute an update in joint coordinates using the jacobian transpose
        Δq = α * Array(Jp)' * (transform(state, desired, world) - point_in_world).v
        ## Apply the update
        q .= configuration(state) .+ Δq
        set_configuration!(state, q)
    end
    state
end

# To use our IK method, we just have to set our current state and choose a desired location for the tip of the robot's arm:

rand!(state)
set_configuration!(vis, configuration(state))

# Choose a desired location. We'll move the tip of the arm to
# [0.5, 0, 2]
desired_tip_location = Point3D(root_frame(mechanism), 0.5, 0, 2)
# Run the IK, updating `state` in place
jacobian_transpose_ik!(state, body, point, desired_tip_location)
set_configuration!(vis, configuration(state))

# We asked for our point to be close to [0.5, 0, 2],
# but since the arm cannot move in the y direction at all
# we end up near [0.5, 0.25, 2] instead
transform(state, point, root_frame(mechanism))

# We can try varying the target and watching the IK solution change:

qs = typeof(configuration(state))[]

# Vary the desired x position from -1 to 1
for x in range(-1, stop=1, length=100)
    desired = Point3D(root_frame(mechanism), x, 0, 2)
    jacobian_transpose_ik!(state, body, point, desired)
    push!(qs, copy(configuration(state)))
end
ts = collect(range(0, stop=1, length=length(qs)))
setanimation!(vis, Animation(vis, ts, qs))

# ## Control

# Now let's use the same principle to generate torques and actually control the robot. To make things more interesting, let's get the end of the robot's arm to trace out a circle.

circle_origin = SVector(0., 0.25, 2)
radius = 0.5
ω = 1.0  # radians per second at which the point should move in its circle

using GeometryTypes: Point

## Draw the circle in the viewer
θ = repeat(range(0, stop=2π, length=100), inner=(2,))[2:end]
cx, cy, cz = circle_origin
geometry = PointCloud(Point.(cx .+ radius .* sin.(θ), cy, cz .+ 0.5 .* cos.(θ)))
setobject!(vis[:circle], LineSegments(geometry, LineBasicMaterial()))

# This function will take in the parameters of the circle
# and the target point and return a function we can use
# as the controller. By wrapping the construction of the
# controller in this way, we avoid any issues with accessing
# non-const global variables.
function make_circle_controller(state::MechanismState,
                                body::RigidBody,
                                point::Point3D,
                                circle_origin::AbstractVector,
                                radius,
                                ω)
    mechanism = state.mechanism
    world = root_frame(mechanism)
    joint_path = path(mechanism, root_body(mechanism), body)
    point_world = transform(state, point, root_frame(mechanism))
    Jp = point_jacobian(state, joint_path, point_world)
    v̇ = similar(velocity(state))

    function controller!(τ, t, state)
        desired = Point3D(world, circle_origin .+ radius .* SVector(sin(t / ω), 0, cos(t / ω)))
        point_in_world = transform_to_root(state, body) * point
        point_jacobian!(Jp, state, joint_path, point_in_world)
        Kp = 200
        Kd = 20
        Δp = desired - point_in_world
        v̇ .= Kp * Array(Jp)' * Δp.v .- 20 .* velocity(state)
        τ .= inverse_dynamics(state, v̇)
    end
end
controller! = make_circle_controller(state, body, point, circle_origin, radius, ω)
ts, qs, vs = simulate(state, 10, controller!);

# Animate the resulting trajectory:

setanimation!(vis, Animation(vis, ts, qs))

# Now we can plot the behavior of the controller. The initial state is quite far from the target, so there's some significant overshoot early in the trajectory, but the controller eventually settles into tracking the desired circular path. This controller isn't very well-tuned, and we could certainly do better with a more advanced approach, but this is still a nice demonstration of a very simple control policy.

using Plots: gr, plot
gr()

xs = Float64[]
zs = Float64[]

## Downsample by 100 just so the plot doesn't become a huge file:
for q in qs[1:100:end]
    set_configuration!(state, q)
    p = transform(state, point, root_frame(mechanism))
    push!(xs, p.v[1])
    push!(zs, p.v[3])
end

plot(xs, zs, xlim=(-1, 1), ylim=(1, 3), aspect_ratio=:equal, legend=nothing)

max_stars_repo_path: examples/4. Jacobian IK and Control/4. Jacobian IK and Control.jl


text:
using JuMP, ConditionalJuMP, Cbc
using Test

const viscosity = 100
const stiffness = 100
const penetration = 0.1
const normal_force = stiffness * penetration
const μ = 0.5
const Δt = 0.05

# A very simple stick-slip dynamics. Frictional force is equal to
# -(viscosity * velocity), clipped to lie within the friction cone
# at ±(μ * normal_force).
function dynamics(x)
    q = x[1]
    v = x[2]
    viscous_force = -viscosity * v
    frictional_force = @switch(
        (viscous_force ≤ -μ * normal_force)                                        => -μ * normal_force,
        ((viscous_force ≥ -μ * normal_force) & (viscous_force ≤ μ * normal_force)) => viscous_force,
        (viscous_force ≥ μ * normal_force)                                         => μ * normal_force
    )
    ẋ = [v, frictional_force]
    ẋ
end

# Naive forward-euler discretization of the continuous dynamics
function update(x)
    ẋ = dynamics(x)
    x .+ ẋ .* Δt
end

# Simulate forward using the given dynamics to get a trajectory of states
N = 20
x0 = [0, 3.0]
xsim = [x0]
for i in 2:N
    push!(xsim, update(xsim[end]))
end

# Now set up a mixed-integer program to solve for all N states simultaneously
m = Model(solver=CbcSolver(logLevel=0))
xvars = [@variable(m, [1:2], basename="x", lowerbound=-4, upperbound=4) for i in 1:N]
@constraint(m, xvars[1] .== x0)
for i in 2:N
    @constraint(m, xvars[i] .== update(xvars[i - 1]))
end
warmstart!(m)
solve(m)
xvals = getvalue.(xvars)

# Verify that the solution to the mixed-integer program is exactly the same as
# the simulation
@test all([xvals[i] ≈ xsim[i] for i in 1:N])

expected = Array{Float64,1}[[0.0, 3.0], [0.15, 2.75], [0.2875, 2.5], [0.4125, 2.25], [0.525, 2.0], [0.625, 1.75], [0.7125, 1.5], [0.7875, 1.25], [0.85, 1.0], [0.9, 0.75], [0.9375, 0.5], [0.9625, 0.25], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0], [0.975, 0.0]]
@test all([xvals[i] ≈ expected[i] for i in 1:N])
@test all([xsim[i] ≈ expected[i] for i in 1:N])


text:
module Node_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

B = BnBModel([Interval(1.0,2.0),Interval(3.0,4.0)])
A1,A2,A3,A4,A5 = EAGOBranchBound.NS_depth_breadth(B)
@test A1 == [Interval(1.0,2.0),Interval(3.0,4.0)]
@test A2 == -Inf
@test A3 == Inf
@test A4 == 1
@test A5 == 1

end