text:
# Autogenerated wrapper script for Ghostscript_jll for aarch64-linux-musl
export dvipdf, eps2eps, gs, gsbj, gsdj, gsdj500, gslj, gslp, gsnd, pdf2dsc, pdf2ps, pf2afm, pfbtopfa, pphs, printafm, ps2ascii, ps2epsi, ps2pdf, ps2pdf12, ps2pdf13, ps2pdf14, ps2pdfwr, ps2ps, ps2ps2

JLLWrappers.@generate_wrapper_header("Ghostscript")
JLLWrappers.@declare_file_product(dvipdf)
JLLWrappers.@declare_file_product(eps2eps)
JLLWrappers.@declare_executable_product(gs)
JLLWrappers.@declare_file_product(gsbj)
JLLWrappers.@declare_file_product(gsdj)
JLLWrappers.@declare_file_product(gsdj500)
JLLWrappers.@declare_file_product(gslj)
JLLWrappers.@declare_file_product(gslp)
JLLWrappers.@declare_file_product(gsnd)
JLLWrappers.@declare_file_product(pdf2dsc)
JLLWrappers.@declare_file_product(pdf2ps)
JLLWrappers.@declare_file_product(pf2afm)
JLLWrappers.@declare_file_product(pfbtopfa)
JLLWrappers.@declare_file_product(pphs)
JLLWrappers.@declare_file_product(printafm)
JLLWrappers.@declare_file_product(ps2ascii)
JLLWrappers.@declare_file_product(ps2epsi)
JLLWrappers.@declare_file_product(ps2pdf)
JLLWrappers.@declare_file_product(ps2pdf12)
JLLWrappers.@declare_file_product(ps2pdf13)
JLLWrappers.@declare_file_product(ps2pdf14)
JLLWrappers.@declare_file_product(ps2pdfwr)
JLLWrappers.@declare_file_product(ps2ps)
JLLWrappers.@declare_file_product(ps2ps2)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        dvipdf,
        "bin/dvipdf",
    )

    JLLWrappers.@init_file_product(
        eps2eps,
        "bin/eps2eps",
    )

    JLLWrappers.@init_executable_product(
        gs,
        "bin/gs",
    )

    JLLWrappers.@init_file_product(
        gsbj,
        "bin/gsbj",
    )

    JLLWrappers.@init_file_product(
        gsdj,
        "bin/gsdj",
    )

    JLLWrappers.@init_file_product(
        gsdj500,
        "bin/gsdj500",
    )

    JLLWrappers.@init_file_product(
        gslj,
        "bin/gslj",
    )

    JLLWrappers.@init_file_product(
        gslp,
        "bin/gslp",
    )

    JLLWrappers.@init_file_product(
        gsnd,
        "bin/gsnd",
    )

    JLLWrappers.@init_file_product(
        pdf2dsc,
        "bin/pdf2dsc",
    )

    JLLWrappers.@init_file_product(
        pdf2ps,
        "bin/pdf2ps",
    )

    JLLWrappers.@init_file_product(
        pf2afm,
        "bin/pf2afm",
    )

    JLLWrappers.@init_file_product(
        pfbtopfa,
        "bin/pfbtopfa",
    )

    JLLWrappers.@init_file_product(
        pphs,
        "bin/pphs",
    )

    JLLWrappers.@init_file_product(
        printafm,
        "bin/printafm",
    )

    JLLWrappers.@init_file_product(
        ps2ascii,
        "bin/ps2ascii",
    )

    JLLWrappers.@init_file_product(
        ps2epsi,
        "bin/ps2epsi",
    )

    JLLWrappers.@init_file_product(
        ps2pdf,
        "bin/ps2pdf",
    )

    JLLWrappers.@init_file_product(
        ps2pdf12,
        "bin/ps2pdf12",
    )

    JLLWrappers.@init_file_product(
        ps2pdf13,
        "bin/ps2pdf13",
    )

    JLLWrappers.@init_file_product(
        ps2pdf14,
        "bin/ps2pdf14",
    )

    JLLWrappers.@init_file_product(
        ps2pdfwr,
        "bin/ps2pdfwr",
    )

    JLLWrappers.@init_file_product(
        ps2ps,
        "bin/ps2ps",
    )

    JLLWrappers.@init_file_product(
        ps2ps2,
        "bin/ps2ps2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MotorControl

println("\nMotorControl: Demonstrating the ability to simulate hierarchical mixed domain models")

using Modia
@usingModiaPlot
using  ModiaLang.Measurements
export MotorControl2

include("$(Modia.modelsPath)/Blocks.jl")
include("$(Modia.modelsPath)/Electric.jl")
include("$(Modia.modelsPath)/Rotational.jl")

MotorControl1 = Model(
    step = Step | Map(height=4.7*u"A", offset=0u"A"),
    feedback = Feedback,
    PI = PI | Map(T=0.005u"s", k=30, x = Var(init=0.0u"A")),
    firstOrder = FirstOrder | Map(k=1u"V/A", T=0.001u"s", x = Var(init=0.0u"V")),

    signalVoltage = SignalVoltage,
    resistor = Resistor | Map(R=13.8u"Ω"),
    inductor = Inductor | Map(L=0.061u"H"),
    emf = EMF | Map(k=1.016u"N*m/A"),
    ground = Ground,
    currentSensor = CurrentSensor,

    motorInertia = Inertia | Map(J=0.0025u"kg*m^2"),
    idealGear = IdealGear | Map(ratio=105),
    springDamper = SpringDamper | Map(c=5.0e5u"N*m/rad", d=500u"N*m*s/rad"),
    loadInertia = Inertia | Map(J=100u"kg*m^2"),
    tload = Torque,

    equations = :[
      tload.tau = 0.0u"N*m",
    ],

    connect = :[
      (step.y, feedback.u1)
      (feedback.y, PI.u)
      (PI.y, firstOrder.u)
      (firstOrder.y, signalVoltage.v)

      (signalVoltage.p, resistor.p)
      (resistor.n, inductor.p)
      (inductor.n, emf.p)
      (emf.n, ground.p, currentSensor.p)
      (currentSensor.n, signalVoltage.n)
      (currentSensor.i, feedback.u2)

      (emf.flange, motorInertia.flange_a)
      (motorInertia.flange_b, idealGear.flange_a)
      (idealGear.flange_b, springDamper.flange_a)
      (springDamper.flange_b, loadInertia.flange_a)
      (tload.flange, loadInertia.flange_b) ]
)

model = @instantiateModel(MotorControl1, log=false)
println("Simulate")
@time simulate!(model, stopTime=0.1, tolerance=1e-6, log=false, requiredFinalStates =
    [3.487844601078223, 106.82874860310305, 4.616087152802267, 2.014120727821878, 41.98048886114646, 0.018332432934516536, 0.3725930536373392])
plot(model, [("currentSensor.i", "step.y"), "loadInertia.w"], figure=1)

# Hierarchical model

ControlledMotor = Model(
  refCurrent = input,
  feedback = Feedback,
  PI = PI | Map(T=0.005u"s", k=30, x = Var(init=0.0u"A")),
  firstOrder = FirstOrder | Map(k=1u"V/A", T=0.001u"s", x = Var(init=0.0u"V")),

  signalVoltage = SignalVoltage,
  resistor = Resistor | Map(R=13.8u"Ω"),
  inductor = Inductor | Map(L=0.061u"H"),
  emf = EMF | Map(k=1.016u"N*m/A"),
  ground = Ground,
  currentSensor = CurrentSensor,

  motorInertia = Inertia | Map(J=0.0025u"kg*m^2"),
  flange = Flange,

  connect = :[
    (refCurrent, feedback.u1)
    (feedback.y, PI.u)
    (PI.y, firstOrder.u)
    (firstOrder.y, signalVoltage.v)

    (signalVoltage.p, resistor.p)
    (resistor.n, inductor.p)
    (inductor.n, emf.p)
    (emf.n, ground.p, currentSensor.p)
    (currentSensor.n, signalVoltage.n)
    (currentSensor.i, feedback.u2)

    (emf.flange, motorInertia.flange_a)
    (motorInertia.flange_b, flange) ]
)

MotorControl2 = Model(
  step = Step | Map(height=4.7*u"A", offset=0u"A"),

  controlledMotor = ControlledMotor,

  idealGear = IdealGear | Map(ratio=105),
  springDamper = SpringDamper | Map(c=5.0e5u"N*m/rad", d=500u"N*m*s/rad"),
  loadInertia = Inertia | Map(J=100.0u"kg*m^2"),
  tload = Torque,

  equations = :[
    tload.tau = 0.0u"N*m",
  ],

  connect = :[
    (step.y, controlledMotor.refCurrent)

    (controlledMotor.flange, idealGear.flange_a)
    (idealGear.flange_b, springDamper.flange_a)
    (springDamper.flange_b, loadInertia.flange_a)
    (tload.flange, loadInertia.flange_b) ]
)

model = @instantiateModel(MotorControl2)
println("Simulate")
@time simulate!(model, stopTime=0.1, tolerance=1e-6, log=false, requiredFinalStates =
    [3.487844601078223, 106.82874860310305, 4.616087152802267, 2.014120727821878, 41.98048886114646, 0.018332432934516536, 0.3725930536373392])
plot(model, [("controlledMotor.currentSensor.i", "step.y"), "loadInertia.w"], figure=1)

# Model with uncertainties

MotorControlWithUncertainties = MotorControl2 | Map( loadInertia = Map(J=(100.0 ± 10)*u"kg*m^2"), controlledMotor = Map(PI = Map(k=30 ± 3) ) )

model = @instantiateModel(MotorControlWithUncertainties, FloatType = Measurement{Float64})
println("Simulate")
@time simulate!(model, stopTime=0.1, tolerance=1e-6, log=false)
plot(model, [("controlledMotor.currentSensor.i", "step.y"), "loadInertia.w"], figure=1)

end

#=
module MotorControlModuleMonteCarlo

using ModiaBase.Modia
using Unitful
using Main.MotorControlModule
include("../test/SimulateAndPlot.jl")

using MonteCarloMeasurements

# Model with Monte Carlo

MotorControlWithUncertainties = MotorControl2 | Map( loadInertia = Map(J=(100.0 ∓ 10)), controlledMotor = Map(PI = Map(k=30 ∓ 3) ) )

model = @instantiateModel(MotorControlWithUncertainties, FloatType = StaticParticles{Float64,100}, unitless=true)
println("Simulate")
@time simulate!(model, stopTime=0.1, tolerance=1e-6, log=false)
plot(model, [("controlledMotor.currentSensor.i", "step.y"), "loadInertia.w"], figure=1)

end
=#

ext: jl


text:
export eliminate, eliminate!, recover!

"""
    eliminate!(eg::EliminateGraph, vertices)

Eliminate vertices from a graph.
"""
function eliminate!(eg::EliminateGraph, vi::Int)
    N = nv0(eg)
    @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
    for j in N-eg.nv+1:N
        @inbounds vj = eg.vertices[j]
        if vj==vi
            iptr += 1
            eg.level += 1
            @inbounds eg.ptr[eg.level] = iptr
            unsafe_swap!(eg.vertices, j, iptr)
            break
        end
    end
    eg.nv -= 1
    return eg
end

function eliminate!(eg::EliminateGraph, vs)
    N = nv0(eg)
    @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
    for vi in vs
        for j in N-eg.nv+1:N
            @inbounds vj = eg.vertices[j]
            if vj==vi
                iptr += 1
                unsafe_swap!(eg.vertices, j, iptr)
                break
            end
        end
    end
    eg.level += 1
    eg.nv -= length(vs)
    @inbounds eg.ptr[eg.level] = iptr
    return eg
end

# the fallback
function eliminate!(eg::EliminateGraph, nc::AbstractVertexSet)
    eliminate!(eg, generate_set(eg, nc))
end

"""
    eliminate([func], eg::EliminateGraph, vertices)
    eg \\ vertices

Eliminate vertices from a graph, return the value of `func(eliminated_graph)` if `func` provided.
"""
eliminate(eg, vertices) = eliminate!(copy(eg), vertices)

@inline function eliminate(func, eg::EliminateGraph, vi)
    eliminate!(eg, vi)
    res = func(eg)
    recover!(eg)
    return res
end

Base.:\(eg::EliminateGraph, vertices) = eliminate(eg, vertices)

"""restore eliminated vertices for a level (one call of elimintion function)."""
function recover!(eg::EliminateGraph)
    @inbounds eg.nv += eg.ptr[eg.level] - (eg.level==1 ? 0 : eg.ptr[eg.level-1])
    eg.level -= 1
    eg
end

##################### Specialized ######################
@generated function eliminate!(eg::EliminateGraph, nc::Neighbors{SP,1}) where SP
    condition = SP == CLOSED ? :(vi == nc.i || isconnected(eg, vi, nc.i)) : :(isconnected(eg, vi, nc.i))
    quote
        @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
        for i in nv0(eg)-eg.nv+1:nv0(eg)
            @inbounds vi = eg.vertices[i]
            $condition || continue
            iptr += 1
            eg.nv -= 1
            unsafe_swap!(eg.vertices, i, iptr)
        end
        eg.level += 1
        @inbounds eg.ptr[eg.level] = iptr
        return eg
    end
end

function eliminate!(eg::EliminateGraph, set::Mirrors{CLOSED})
    nc1 = generate_set(eg, NeighborCover(set.i))

    @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
    for i in nv0(eg)-eg.nv+1:nv0(eg)
        @inbounds vi = eg.vertices[i]
        vi == set.i || ((!isconnected(eg,vi,set.i) && any(nb->isconnected(eg,vi,nb), nc1)) && isclique(eg, setdiff(nc1,neighbors(eg,vi)))) || continue
        iptr += 1
        eg.nv -= 1
        unsafe_swap!(eg.vertices, i, iptr)
    end
    eg.level += 1
    @inbounds eg.ptr[eg.level] = iptr
    return eg
end

function eliminate!(eg::EliminateGraph, set::Mirrors{OPEN}) # can not generated since any is not pure
    nc1 = generate_set(eg, NeighborCover(set.i))

    @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
    for i in nv0(eg)-eg.nv+1:nv0(eg)
        @inbounds vi = eg.vertices[i]
        (vi!=set.i && !isconnected(eg,vi,set.i) && any(nb->isconnected(eg,vi,nb), nc1)) && isclique(eg, setdiff(nc1,neighbors(eg,vi))) || continue
        iptr += 1
        eg.nv -= 1
        unsafe_swap!(eg.vertices, i, iptr)
    end
    eg.level += 1
    @inbounds eg.ptr[eg.level] = iptr
    return eg
end

@generated function eliminate!(eg::EliminateGraph, nc::UnionOf)
    if nc <: UnionOf{Neighbors{SP,1}, Neighbors{SP,1}} where SP
        SP = nc.parameters[1].parameters[1]
        condition = SP == CLOSED ? :(vi == nc.A.i || vi == nc.B.i || isconnected(eg, vi, nc.A.i) ||
            isconnected(eg, vi, nc.B.i)) : :(isconnected(eg, vi, nc.A.i) || isconnected(eg, vi, nc.B.i))
    elseif nc <: UnionOf{Neighbors{SP,1}, UnionOf{Neighbors{SP,1},Vertex}} where SP
        SP = nc.parameters[1].parameters[1]
        condition = SP == CLOSED ? :(vi == nc.A.i || vi == nc.B.A.i || vi == nc.B.B.i || isconnected(eg, vi, nc.A.i) ||
            isconnected(eg, vi, nc.B.A.i)) : :(vi == nc.B.B.i || isconnected(eg, vi, nc.A.i) || isconnected(eg, vi, nc.B.A.i))
    else
        return :(invoke(eliminate!, Tuple{EliminateGraph,AbstractVertexSet}, eg, nc))
    end

    quote
        @inbounds iptr = eg.level == 0 ? 0 : eg.ptr[eg.level]
        for i in nv0(eg)-eg.nv+1:nv0(eg)
            @inbounds vi = eg.vertices[i]
            $condition || continue
            iptr += 1
            eg.nv -= 1
            unsafe_swap!(eg.vertices, i, iptr)
        end
        eg.level += 1
        @inbounds eg.ptr[eg.level] = iptr
        return eg
    end
end


text:
result = 0
n = 1000_000
for i in 1:n
    result += sin(rand())
end
show(result)

# @time include("1_5_glob.jl")
# 459652.1136757382  0.351386 seconds (5.00 M allocations: 91.616 MB, 3.12% gc time)

max_stars_repo_path: Section 1/1_5_glob.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_5_glob.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_5_glob.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
####
#### Regenerate the training curves displayed on the tutorial and on the README
####

using AlphaZero
using JSON3

# Defining the benchmarks
Training = Examples.ConnectFour.Training
netonly = Benchmark.NetworkOnly(τ=0.5)
alphazero = Benchmark.Full(Training.arena.mcts)
minmax = Benchmark.MinMaxTS(depth=5, amplify_rewards=true, τ=0.2)
mcts = Training.mcts_baseline
alphazero_benchmark = [
  Benchmark.Duel(alphazero, mcts, Training.benchmark_sim),
  Benchmark.Duel(alphazero, minmax, Training.benchmark_sim)]
netonly_benchmark = [
  Benchmark.Duel(netonly, mcts, Training.benchmark_sim),
  Benchmark.Duel(netonly, minmax, Training.benchmark_sim)]
full_benchmark = [alphazero_benchmark; netonly_benchmark]

dir = UserInterface.default_session_dir("connect-four")
logger = UserInterface.Log.Logger()

function rerun_benchmarks()
  UserInterface.run_new_benchmark(dir, "alphazero", alphazero_benchmark; logger)
  UserInterface.run_new_benchmark(dir, "netonly", netonly_benchmark; logger)
  UserInterface.run_new_benchmark(dir, "full", full_benchmark; logger)
end

# Reuse the JSON benchmark data to make a different plot.
function replot_from_json(itcmax=nothing)
  repdir = joinpath(dir, "netonly")
  repfile = joinpath(repdir, UserInterface.BENCHMARK_FILE)
  reports = open(repfile, "r") do io
    JSON3.read(io, Vector{Report.Benchmark})
  end
  params = UserInterface.load_env(dir).params
  if !isnothing(itcmax)
    reports = reports[1:itcmax+1]
  end
  UserInterface.plot_benchmark(params, reports, repdir)
end


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: CUFFTAPI __attribute__ ( ( visibility ( "default" ) ) )

const CUFFT_VER_MAJOR = 10
const CUFFT_VER_MINOR = 2
const CUFFT_VER_PATCH = 1
const CUFFT_VER_BUILD = 245
const CUFFT_VERSION = 10201
const MAX_CUFFT_ERROR = 0x11
const CUFFT_FORWARD = -1
const CUFFT_INVERSE = 1

@cenum cufftCompatibility_t::UInt32 begin
    CUFFT_COMPATIBILITY_FFTW_PADDING = 1
end

const CUFFT_COMPATIBILITY_DEFAULT = CUFFT_COMPATIBILITY_FFTW_PADDING
const MAX_SHIM_RANK = 3

@cenum cufftResult_t::UInt32 begin
    CUFFT_SUCCESS = 0
    CUFFT_INVALID_PLAN = 1
    CUFFT_ALLOC_FAILED = 2
    CUFFT_INVALID_TYPE = 3
    CUFFT_INVALID_VALUE = 4
    CUFFT_INTERNAL_ERROR = 5
    CUFFT_EXEC_FAILED = 6
    CUFFT_SETUP_FAILED = 7
    CUFFT_INVALID_SIZE = 8
    CUFFT_UNALIGNED_DATA = 9
    CUFFT_INCOMPLETE_PARAMETER_LIST = 10
    CUFFT_INVALID_DEVICE = 11
    CUFFT_PARSE_ERROR = 12
    CUFFT_NO_WORKSPACE = 13
    CUFFT_NOT_IMPLEMENTED = 14
    CUFFT_LICENSE_ERROR = 15
    CUFFT_NOT_SUPPORTED = 16
end

const cufftResult = cufftResult_t
const cufftReal = Cfloat
const cufftDoubleReal = Cdouble
const cufftComplex = cuComplex
const cufftDoubleComplex = cuDoubleComplex

@cenum cufftType_t::UInt32 begin
    CUFFT_R2C = 42
    CUFFT_C2R = 44
    CUFFT_C2C = 41
    CUFFT_D2Z = 106
    CUFFT_Z2D = 108
    CUFFT_Z2Z = 105
end

const cufftType = cufftType_t
const cufftCompatibility = cufftCompatibility_t
const cufftHandle = Cint


text:
function start()
    sysimg = project_relative_path("precompile", "StatisticsWithJulia.so")
    root = project_relative_path()
    exe = joinpath(Sys.BINDIR, "julia")
    script = :(using Pkg;
               Pkg.activate($root);
               using Pluto;
               Pluto.run(notebook = $(joinpath(root, "index.jl"))))
    if isfile(sysimg)
        run(`$exe -J$sysimg -e $script`)
    else
        run(`$exe -e $script`)
    end
end

function update()
    @info "Performing an automatic update while keeping local changes.
    If this fails, please run manually `git pull` in the directory
    `$(project_relative_path())`."
    current_dir = pwd()
    cd(project_relative_path())
    if !isempty(readlines(`git diff --stat`))
        run(`git add -u`)
        run(`git commit -m "automatic commit of local changes"`)
    end
    run(`git pull -s recursive -X patience -X ours -X ignore-all-space --no-edit`)
    cd(current_dir)
    Pkg.activate(project_relative_path())
    Pkg.instantiate()
end

function create_sysimage()
    exe = joinpath(Sys.BINDIR, "julia")
    run(`$exe $(project_relative_path("precompile", "precompile.jl"))`)
end

if isfile(project_relative_path("precompile", "StatisticsWithJulia.so"))
    @warn "You may have to create a new system image with this update."
end


text:
@testset "UnicycleGame" begin
    # Test UnicycleGame
    p = 3
    model = UnicycleGame(p=p)
    @test model.n == 4*p
    @test model.m == 2p
    @test model.p == p
    @test model.ni == 4*ones(Int,p)
    @test model.mi == 2*ones(Int,p)
    @test model.pu[1] == SVector{2,Int}((1,4))
    @test model.pu[2] == SVector{2,Int}((2,5))
    @test model.pu[3] == SVector{2,Int}((3,6))

    @test model.px[1] == SVector{2,Int}((1,4))
    @test model.px[2] == SVector{2,Int}((2,5))
    @test model.px[3] == SVector{2,Int}((3,6))

    @test model.pz[1] == SVector{4,Int}((1,4,7,10))
    @test model.pz[2] == SVector{4,Int}((2,5,8,11))
    @test model.pz[3] == SVector{4,Int}((3,6,9,12))

    @test size(model) == (4*p, 2p, model.pu, p)

    T = Float64
    x = ones(SVector{model.n,T})
    u = 100*ones(SVector{model.m,T})
    @test (@ballocated dynamics($model, $x, $u)) == 0
    # @btime dynamics($model, $x, $u)
end


text:
using Base.Test
using ProximalOperators

srand(123)
x = randn(100)
y0 = randn(100)
y = copy(y0)

prox_col = [NormL1(0.1),IndBallL0(1)]
ind_col = [1:50,51:100]
M1 = ones(Bool,size(x))
M1[51:end] = false
M2 = M1.==false

f = SlicedSeparableSum(prox_col,ind_col)
#test second constructor and mask
ff = SlicedSeparableSum([(prox_col[1], M1), (prox_col[2], M2)])
fff = SlicedSeparableSum([prox_col[1] => M1, prox_col[2] => M2])
y,fy = prox(f,x,1.)
yy,fyy = prox(ff,x,1.)
yyy,fyyy = prox(fff,x,1.)

@test norm(fyy-fy)<1e-11
@test norm(y-yy)<1e-11

@test norm(fyyy-fy)<1e-11
@test norm(y-yyy)<1e-11

y1,fy1 = prox(prox_col[1],x[ind_col[1]],1.)
y2,fy2 = prox(prox_col[2],x[ind_col[2]],1.)

@test norm((fy1+fy2)-fy)<1e-11
@test norm(y-[y1;y2])<1e-11

X1,X2 = randn(10,10),randn(10,10)
X = [X1 X2]

f = SlicedSeparableSum([NormL1(1.), NormL21(0.1)], [1:10,11:20], 2)
#test second constructor
ff = SlicedSeparableSum([(NormL1(1.),1:10), (NormL21(0.1),11:20)], 2)

y,fy = prox(f,X,1.)
yy,fyy = prox(ff,X,1.)

@test norm(fyy-fy)<1e-11
@test norm(y-yy)<1e-11

y1,fy1 = prox(NormL1(1.),X1,1.)
y2,fy2 = prox(NormL21(0.1),X2,1.)

@test norm((fy1+fy2)-fy)<1e-11
@test norm(y-[y1 y2])<1e-11

println(ff)
# println(SlicedSeparableSum(repmat([NormL1(1.)],6),repmat([1:10],6),2 ))


text:
using PhysicsEngine
using Test

cd("../")

projectName = "PhysicsEngine"
testFiles = [
  "simulation_test.jl"
  "grid_test.jl"
  "vector_test.jl"
  "particle_test.jl"
  "renderer_test.jl"
  "text_test.jl"
]

for file in testFiles
  include(file)
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
# string of ints, describe a tree structure.
# header: number of children, number of "metadata" entries
# children: (length unknown)
# metadata:

struct 🎄
    children::Vector{🎄}
    metadata::Vector{Int}
end

function extract_tree(content)
    nchil, nmeta, content = content[1], content[2], @view(content[3:end])
    children = 🎄[]
    for _ = 1:nchil
        child, content = extract_tree(content)
        push!(children, child)
    end
    metadata, content = content[1:nmeta], @view(content[nmeta+1:end])
    🎄(children, metadata), content
end

star1(tree::🎄) = sum(tree.metadata) + mapreduce(star1, +, tree.children, init=0)

# I thougth this might be faster but I was way wrong
star11(tree::🎄) = sum(tree.metadata) +
    (isempty(tree.children) ? 0 : sum(star11(c) for c in tree.children))

test_input = parse.(Int, split("2 3 0 3 10 11 12 1 1 0 1 99 2 1 1 2"))
test_tree, _ = extract_tree(test_input);
star1(test_tree) == 138

input = parse.(Int, split(read("08.input", String)));
tree, _ = extract_tree(input);
star1(tree) # = 36891

star2(tree::🎄, child::Int) =
    child ∈ 1:length(tree.children) ? star2(tree.children[child]) : 0
star2(tree::🎄) =
    if length(tree.children) > 0
        mapreduce(i->star2(tree, i), +, tree.metadata, init=0)
    else
        sum(tree.metadata)
    end

star2(test_tree) == 66
star2(tree)

################################################################################
# some optimizations??

using BenchmarkTools

# variant that pulls metadata out into its own vector
function extract_tree_noview(content)
    nchil, nmeta, content = content[1], content[2], @view(content[3:end])

    children = 🎄[]
    for _ = 1:nchil
        child, content = extract_tree(content)
        push!(children, child)
    end

    metadata, content = content[1:nmeta], @view(content[nmeta+1:end])

    🎄(children, metadata), content
end

tree_noview, _ = extract_tree_noview(input);

# really does not seem to make any difference.  I'd have expect some difference
# in the allocations but I guess if the vectors are all small allocating the
# view isn't saving you anything over allocating the vector directly
#
# julia> @btime extract_tree($input);
#   99.348 μs (10210 allocations: 479.88 KiB)
# julia> @btime extract_tree_noview($input);
#   99.560 μs (10210 allocations: 480.00 KiB)

# julia> @btime star1($tree)
#   410.655 μs (6380 allocations: 147.94 KiB)
# 36891

# julia> @btime star1($tree_noview)
#   411.760 μs (6380 allocations: 147.94 KiB)
# 36891

# julia> @btime star2($tree)
#   95.911 μs (1291 allocations: 27.67 KiB)
# 20083

# julia> @btime star2($tree_noview)
#   95.952 μs (1287 allocations: 27.58 KiB)
# 20083

function extract_tree_noview_plus(content)
    nchil, nmeta, content = content[1], content[2], content[3:end]

    children = 🎄[]
    for _ = 1:nchil
        child, content = extract_tree(content)
        push!(children, child)
    end

    metadata, content = content[1:nmeta], content[nmeta+1:end]

    🎄(children, metadata), content
end

@btime extract_tree_noview_plus($input);

# take a SMALL hit for duplicating the content every time but still, not nearly
# as I would have thought.
#
# julia> @btime extract_tree_noview_plus($input);
#   112.935 μs (10211 allocations: 591.88 KiB)


text:
########
# Dual #
########

"""
    ForwardDiff.can_dual(V::Type)

Determines whether the type V is allowed as the scalar type in a
Dual. By default, only `<:Real` types are allowed.
"""
can_dual(::Type{<:Real}) = true
can_dual(::Type) = false

struct Dual{T,V,N} <: Real
    value::V
    partials::Partials{N,V}
    function Dual{T, V, N}(value::V, partials::Partials{N, V}) where {T, V, N}
        can_dual(V) || throw_cannot_dual(V)
        new{T, V, N}(value, partials)
    end
end

##############
# Exceptions #
##############

struct DualMismatchError{A,B} <: Exception
    a::A
    b::B
end

Base.showerror(io::IO, e::DualMismatchError{A,B}) where {A,B} =
    print(io, "Cannot determine ordering of Dual tags $(e.a) and $(e.b)")

@noinline function throw_cannot_dual(V::Type)
    throw(ArgumentError("Cannot create a dual over scalar type $V." *
        " If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{$V}) = true."))
end

"""
    ForwardDiff.≺(a, b)::Bool

Determines the order in which tagged `Dual` objects are composed. If true, then `Dual{b}`
objects will appear outside `Dual{a}` objects.

This is important when working with nested differentiation: currently, only the outermost
tag can be extracted, so it should be used in the _innermost_ function.
"""
≺(a,b) = throw(DualMismatchError(a,b))

################
# Constructors #
################

@inline Dual{T}(value::V, partials::Partials{N,V}) where {T,N,V} = Dual{T,V,N}(value, partials)

@inline function Dual{T}(value::A, partials::Partials{N,B}) where {T,N,A,B}
    C = promote_type(A, B)
    return Dual{T}(convert(C, value), convert(Partials{N,C}, partials))
end

@inline Dual{T}(value, partials::Tuple) where {T} = Dual{T}(value, Partials(partials))
@inline Dual{T}(value, partials::Tuple{}) where {T} = Dual{T}(value, Partials{0,typeof(value)}(partials))
@inline Dual{T}(value) where {T} = Dual{T}(value, ())
@inline Dual{T}(x::Dual{T}) where {T} = Dual{T}(x, ())
@inline Dual{T}(value, partial1, partials...) where {T} = Dual{T}(value, tuple(partial1, partials...))
@inline Dual{T}(value::V, ::Chunk{N}, p::Val{i}) where {T,V,N,i} = Dual{T}(value, single_seed(Partials{N,V}, p))
@inline Dual(args...) = Dual{Nothing}(args...)

# we define these special cases so that the "constructor <--> convert" pun holds for `Dual`
@inline Dual{T,V,N}(x::Dual{T,V,N}) where {T,V,N} = x
@inline Dual{T,V,N}(x) where {T,V,N} = convert(Dual{T,V,N}, x)
@inline Dual{T,V,N}(x::Number) where {T,V,N} = convert(Dual{T,V,N}, x)
@inline Dual{T,V}(x) where {T,V} = convert(Dual{T,V}, x)

##############################
# Utility/Accessor Functions #
##############################

@inline value(x) = x
@inline value(d::Dual) = d.value

@inline value(::Type{T}, x) where T = x
@inline value(::Type{T}, d::Dual{T}) where T = value(d)
@inline function value(::Type{T}, d::Dual{S}) where {T,S}
    if S ≺ T
        d
    else
        throw(DualMismatchError(T,S))
    end
end

@inline partials(x) = Partials{0,typeof(x)}(tuple())
@inline partials(d::Dual) = d.partials
@inline partials(x, i...) = zero(x)
@inline Base.@propagate_inbounds partials(d::Dual, i) = d.partials[i]
@inline Base.@propagate_inbounds partials(d::Dual, i, j) = partials(d, i).partials[j]
@inline Base.@propagate_inbounds partials(d::Dual, i, j, k...) = partials(partials(d, i, j), k...)

@inline Base.@propagate_inbounds partials(::Type{T}, x, i...) where T = partials(x, i...)
@inline Base.@propagate_inbounds partials(::Type{T}, d::Dual{T}, i...) where T = partials(d, i...)
@inline function partials(::Type{T}, d::Dual{S}, i...) where {T,S}
    if S ≺ T
        zero(d)
    else
        throw(DualMismatchError(T,S))
    end
end

@inline npartials(::Dual{T,V,N}) where {T,V,N} = N
@inline npartials(::Type{Dual{T,V,N}}) where {T,V,N} = N

@inline order(::Type{V}) where {V} = 0
@inline order(::Type{Dual{T,V,N}}) where {T,V,N} = 1 + order(V)

@inline valtype(::V) where {V} = V
@inline valtype(::Type{V}) where {V} = V
@inline valtype(::Dual{T,V,N}) where {T,V,N} = V
@inline valtype(::Type{Dual{T,V,N}}) where {T,V,N} = V

@inline tagtype(::V) where {V} = Nothing
@inline tagtype(::Type{V}) where {V} = Nothing
@inline tagtype(::Dual{T,V,N}) where {T,V,N} = T
@inline tagtype(::Type{Dual{T,V,N}}) where {T,V,N} = T

####################################
# N-ary Operation Definition Tools #
####################################

macro define_binary_dual_op(f, xy_body, x_body, y_body)
    FD = ForwardDiff
    defs = quote
        @inline $(f)(x::$FD.Dual{Txy}, y::$FD.Dual{Txy}) where {Txy} = $xy_body
        @inline $(f)(x::$FD.Dual{Tx}, y::$FD.Dual{Ty}) where {Tx,Ty} = Ty ≺ Tx ? $x_body : $y_body
    end
    for R in AMBIGUOUS_TYPES
        expr = quote
            @inline $(f)(x::$FD.Dual{Tx}, y::$R) where {Tx} = $x_body
            @inline $(f)(x::$R, y::$FD.Dual{Ty}) where {Ty} = $y_body
        end
        append!(defs.args, expr.args)
    end
    return esc(defs)
end

macro define_ternary_dual_op(f, xyz_body, xy_body, xz_body, yz_body, x_body, y_body, z_body)
    FD = ForwardDiff
    defs = quote
        @inline $(f)(x::$FD.Dual{Txyz}, y::$FD.Dual{Txyz}, z::$FD.Dual{Txyz}) where {Txyz} = $xyz_body
        @inline $(f)(x::$FD.Dual{Txy}, y::$FD.Dual{Txy}, z::$FD.Dual{Tz}) where {Txy,Tz} = Tz ≺ Txy ? $xy_body : $z_body
        @inline $(f)(x::$FD.Dual{Txz}, y::$FD.Dual{Ty}, z::$FD.Dual{Txz}) where {Txz,Ty} = Ty ≺ Txz ? $xz_body : $y_body
        @inline $(f)(x::$FD.Dual{Tx}, y::$FD.Dual{Tyz}, z::$FD.Dual{Tyz}) where {Tx,Tyz} = Tyz ≺ Tx ? $x_body  : $yz_body
        @inline function $(f)(x::$FD.Dual{Tx}, y::$FD.Dual{Ty}, z::$FD.Dual{Tz}) where {Tx,Ty,Tz}
            if Tz ≺ Tx && Ty ≺ Tx
                $x_body
            elseif Tz ≺ Ty
                $y_body
            else
                $z_body
            end
        end
    end
    for R in AMBIGUOUS_TYPES
        expr = quote
            @inline $(f)(x::$FD.Dual{Txy}, y::$FD.Dual{Txy}, z::$R) where {Txy} = $xy_body
            @inline $(f)(x::$FD.Dual{Tx}, y::$FD.Dual{Ty}, z::$R)  where {Tx, Ty} = Ty ≺ Tx ? $x_body : $y_body
            @inline $(f)(x::$FD.Dual{Txz}, y::$R, z::$FD.Dual{Txz}) where {Txz} = $xz_body
            @inline $(f)(x::$FD.Dual{Tx}, y::$R, z::$FD.Dual{Tz}) where {Tx,Tz} = Tz ≺ Tx ? $x_body : $z_body
            @inline $(f)(x::$R, y::$FD.Dual{Tyz}, z::$FD.Dual{Tyz}) where {Tyz} = $yz_body
            @inline $(f)(x::$R, y::$FD.Dual{Ty}, z::$FD.Dual{Tz}) where {Ty,Tz} = Tz ≺ Ty ? $y_body : $z_body
        end
        append!(defs.args, expr.args)
        for Q in AMBIGUOUS_TYPES
            Q === R && continue
            expr = quote
                @inline $(f)(x::$FD.Dual{Tx}, y::$R, z::$Q) where {Tx} = $x_body
                @inline $(f)(x::$R, y::$FD.Dual{Ty}, z::$Q) where {Ty} = $y_body
                @inline $(f)(x::$R, y::$Q, z::$FD.Dual{Tz}) where {Tz} = $z_body
            end
            append!(defs.args, expr.args)
        end
        expr = quote
            @inline $(f)(x::$FD.Dual{Tx}, y::$R, z::$R) where {Tx} = $x_body
            @inline $(f)(x::$R, y::$FD.Dual{Ty}, z::$R) where {Ty} = $y_body
            @inline $(f)(x::$R, y::$R, z::$FD.Dual{Tz}) where {Tz} = $z_body
        end
        append!(defs.args, expr.args)
    end
    return esc(defs)
end

function unary_dual_definition(M, f)
    FD = ForwardDiff
    Mf = M == :Base ? f : :($M.$f)
    work = qualified_cse!(quote
        val = $Mf(x)
        deriv = $(DiffRules.diffrule(M, f, :x))
    end)
    return quote
        @inline function $M.$f(d::$FD.Dual{T}) where T
            x = $FD.value(d)
            $work
            return $FD.Dual{T}(val, deriv * $FD.partials(d))
        end
    end
end

function binary_dual_definition(M, f)
    FD = ForwardDiff
    dvx, dvy = DiffRules.diffrule(M, f, :vx, :vy)
    Mf = M == :Base ? f : :($M.$f)
    xy_work = qualified_cse!(quote
        val = $Mf(vx, vy)
        dvx = $dvx
        dvy = $dvy
    end)
    dvx, _ = DiffRules.diffrule(M, f, :vx, :y)
    x_work = qualified_cse!(quote
        val = $Mf(vx, y)
        dvx = $dvx
    end)
    _, dvy = DiffRules.diffrule(M, f, :x, :vy)
    y_work = qualified_cse!(quote
        val = $Mf(x, vy)
        dvy = $dvy
    end)
    expr = quote
        $FD.@define_binary_dual_op(
            $M.$f,
            begin
                vx, vy = $FD.value(x), $FD.value(y)
                $xy_work
                return $FD.Dual{Txy}(val, $FD._mul_partials($FD.partials(x), $FD.partials(y), dvx, dvy))
            end,
            begin
                vx = $FD.value(x)
                $x_work
                return $FD.Dual{Tx}(val, dvx * $FD.partials(x))
            end,
            begin
                vy = $FD.value(y)
                $y_work
                return $FD.Dual{Ty}(val, dvy * $FD.partials(y))
            end
        )
    end
    return expr
end

#####################
# Generic Functions #
#####################

Base.copy(d::Dual) = d

Base.eps(d::Dual) = eps(value(d))
Base.eps(::Type{D}) where {D<:Dual} = eps(valtype(D))

function Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N}
    ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
end

function Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N}
    ForwardDiff.Dual{T}(prevfloat(d.value), d.partials)
end

Base.rtoldefault(::Type{D}) where {D<:Dual} = Base.rtoldefault(valtype(D))

Base.floor(::Type{R}, d::Dual) where {R<:Real} = floor(R, value(d))
Base.floor(d::Dual) = floor(value(d))

Base.ceil(::Type{R}, d::Dual) where {R<:Real} = ceil(R, value(d))
Base.ceil(d::Dual) = ceil(value(d))

Base.trunc(::Type{R}, d::Dual) where {R<:Real} = trunc(R, value(d))
Base.trunc(d::Dual) = trunc(value(d))

Base.round(::Type{R}, d::Dual) where {R<:Real} = round(R, value(d))
Base.round(d::Dual) = round(value(d))

if VERSION ≥ v"1.4"
    Base.div(x::Dual, y::Dual, r::RoundingMode) = div(value(x), value(y), r)
else
    Base.div(x::Dual, y::Dual) = div(value(x), value(y))
end

Base.hash(d::Dual) = hash(value(d))
Base.hash(d::Dual, hsh::UInt) = hash(value(d), hsh)

function Base.read(io::IO, ::Type{Dual{T,V,N}}) where {T,V,N}
    value = read(io, V)
    partials = read(io, Partials{N,V})
    return Dual{T,V,N}(value, partials)
end

function Base.write(io::IO, d::Dual)
    write(io, value(d))
    write(io, partials(d))
end

@inline Base.zero(d::Dual) = zero(typeof(d))
@inline Base.zero(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(zero(V), zero(Partials{N,V}))

@inline Base.one(d::Dual) = one(typeof(d))
@inline Base.one(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(one(V), zero(Partials{N,V}))

@inline Random.rand(rng::AbstractRNG, d::Dual) = rand(rng, value(d))
@inline Random.rand(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(rand(V), zero(Partials{N,V}))
@inline Random.rand(rng::AbstractRNG, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(rand(rng, V), zero(Partials{N,V}))
@inline Random.randn(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randn(V), zero(Partials{N,V}))
@inline Random.randn(rng::AbstractRNG, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randn(rng, V), zero(Partials{N,V}))
@inline Random.randexp(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randexp(V), zero(Partials{N,V}))
@inline Random.randexp(rng::AbstractRNG, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randexp(rng, V), zero(Partials{N,V}))

# Predicates #
#------------#

isconstant(d::Dual) = iszero(partials(d))

for pred in UNARY_PREDICATES
    @eval Base.$(pred)(d::Dual) = $(pred)(value(d))
end

for pred in BINARY_PREDICATES
    @eval begin
        @define_binary_dual_op(
            Base.$(pred),
            $(pred)(value(x), value(y)),
            $(pred)(value(x), y),
            $(pred)(x, value(y))
        )
    end
end

########################
# Promotion/Conversion #
########################

function Base.promote_rule(::Type{Dual{T1,V1,N1}},
                                      ::Type{Dual{T2,V2,N2}}) where {T1,V1,N1,T2,V2,N2}
    # V1 and V2 might themselves be Dual types
    if T2 ≺ T1
        Dual{T1,promote_type(V1,Dual{T2,V2,N2}),N1}
    else
        Dual{T2,promote_type(V2,Dual{T1,V1,N1}),N2}
    end
end

function Base.promote_rule(::Type{Dual{T,A,N}},
                           ::Type{Dual{T,B,N}}) where {T,A,B,N}
    return Dual{T,promote_type(A, B),N}
end

for R in (Irrational, Real, BigFloat, Bool)
    if isconcretetype(R) # issue #322
        @eval begin
            Base.promote_rule(::Type{$R}, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T,promote_type($R, V),N}
            Base.promote_rule(::Type{Dual{T,V,N}}, ::Type{$R}) where {T,V,N} = Dual{T,promote_type(V, $R),N}
        end
    else
        @eval begin
            Base.promote_rule(::Type{R}, ::Type{Dual{T,V,N}}) where {R<:$R,T,V,N} = Dual{T,promote_type(R, V),N}
            Base.promote_rule(::Type{Dual{T,V,N}}, ::Type{R}) where {T,V,N,R<:$R} = Dual{T,promote_type(V, R),N}
        end
    end
end

Base.convert(::Type{Dual{T,V,N}}, d::Dual{T}) where {T,V,N} = Dual{T}(convert(V, value(d)), convert(Partials{N,V}, partials(d)))
Base.convert(::Type{Dual{T,V,N}}, x) where {T,V,N} = Dual{T}(convert(V, x), zero(Partials{N,V}))
Base.convert(::Type{Dual{T,V,N}}, x::Number) where {T,V,N} = Dual{T}(convert(V, x), zero(Partials{N,V}))
Base.convert(::Type{D}, d::D) where {D<:Dual} = d

Base.float(d::Dual{T,V,N}) where {T,V,N} = convert(Dual{T,promote_type(V, Float16),N}, d)
Base.AbstractFloat(d::Dual{T,V,N}) where {T,V,N} = convert(Dual{T,promote_type(V, Float16),N}, d)

###################################
# General Mathematical Operations #
###################################

@inline Base.conj(d::Dual) = d

for (M, f, arity) in DiffRules.diffrules()
    in((M, f), ((:Base, :^), (:NaNMath, :pow), (:Base, :/), (:Base, :+), (:Base, :-))) && continue
    if arity == 1
        eval(unary_dual_definition(M, f))
    elseif arity == 2
        eval(binary_dual_definition(M, f))
    else
        error("ForwardDiff currently only knows how to autogenerate Dual definitions for unary and binary functions.")
    end
end

#################
# Special Cases #
#################

# +/- #
#-----#

@define_binary_dual_op(
    Base.:+,
    begin
        vx, vy = value(x), value(y)
        Dual{Txy}(vx + vy, partials(x) + partials(y))
    end,
    Dual{Tx}(value(x) + y, partials(x)),
    Dual{Ty}(x + value(y), partials(y))
)

@define_binary_dual_op(
    Base.:-,
    begin
        vx, vy = value(x), value(y)
        Dual{Txy}(vx - vy, partials(x) - partials(y))
    end,
    Dual{Tx}(value(x) - y, partials(x)),
    Dual{Ty}(x - value(y), -partials(y))
)

@inline Base.:-(d::Dual{T}) where {T} = Dual{T}(-value(d), -partials(d))

# * #
#---#

@inline Base.:*(d::Dual, x::Bool) = x ? d : (signbit(value(d))==0 ? zero(d) : -zero(d))
@inline Base.:*(x::Bool, d::Dual) = d * x

# / #
#---#

# We can't use the normal diffrule autogeneration for this because (x/y) === (x * (1/y))
# doesn't generally hold true for floating point; see issue #264
@define_binary_dual_op(
    Base.:/,
    begin
        vx, vy = value(x), value(y)
        Dual{Txy}(vx / vy, _div_partials(partials(x), partials(y), vx, vy))
    end,
    Dual{Tx}(value(x) / y, partials(x) / y),
    begin
        v = value(y)
        divv = x / v
        Dual{Ty}(divv, -(divv / v) * partials(y))
    end
)

# exponentiation #
#----------------#

for f in (:(Base.:^), :(NaNMath.pow))
    @eval begin
        @define_binary_dual_op(
            $f,
            begin
                vx, vy = value(x), value(y)
                expv = ($f)(vx, vy)
                powval = vy * ($f)(vx, vy - 1)
                if isconstant(y)
                    logval = one(expv)
                elseif iszero(vx) && vy > 0
                    logval = zero(vx)
                else
                    logval = expv * log(vx)
                end
                new_partials = _mul_partials(partials(x), partials(y), powval, logval)
                return Dual{Txy}(expv, new_partials)
            end,
            begin
                v = value(x)
                expv = ($f)(v, y)
                if y == zero(y)
                    new_partials = zero(partials(x))
                else
                    new_partials = partials(x) * y * ($f)(v, y - 1)
                end
                return Dual{Tx}(expv, new_partials)
            end,
            begin
                v = value(y)
                expv = ($f)(x, v)
                deriv = (iszero(x) && v > 0) ? zero(expv) : expv*log(x)
                return Dual{Ty}(expv, deriv * partials(y))
            end
        )
    end
end

@inline Base.literal_pow(::typeof(^), x::Dual{T}, ::Val{0}) where {T} =
    Dual{T}(one(value(x)), zero(partials(x)))

for y in 1:3
    @eval @inline function Base.literal_pow(::typeof(^), x::Dual{T}, ::Val{$y}) where {T}
        v = value(x)
        expv = v^$y
        deriv = $y * v^$(y - 1)
        return Dual{T}(expv, deriv * partials(x))
    end
end

# hypot #
#-------#

@inline function calc_hypot(x, y, z, ::Type{T}) where T
    vx = value(x)
    vy = value(y)
    vz = value(z)
    h = hypot(vx, vy, vz)
    p = (vx / h) * partials(x) + (vy / h) * partials(y) + (vz / h) * partials(z)
    return Dual{T}(h, p)
end

@define_ternary_dual_op(
    Base.hypot,
    calc_hypot(x, y, z, Txyz),
    calc_hypot(x, y, z, Txy),
    calc_hypot(x, y, z, Txz),
    calc_hypot(x, y, z, Tyz),
    calc_hypot(x, y, z, Tx),
    calc_hypot(x, y, z, Ty),
    calc_hypot(x, y, z, Tz),
)

# fma #
#-----#

@generated function calc_fma_xyz(x::Dual{T,<:Any,N},
                                 y::Dual{T,<:Any,N},
                                 z::Dual{T,<:Any,N}) where {T,N}
    ex = Expr(:tuple, [:(fma(value(x), partials(y)[$i], fma(value(y), partials(x)[$i], partials(z)[$i]))) for i in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        v = fma(value(x), value(y), value(z))
        return Dual{T}(v, $ex)
    end
end

@inline function calc_fma_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
    vx, vy = value(x), value(y)
    result = fma(vx, vy, z)
    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
end

@generated function calc_fma_xz(x::Dual{T,<:Any,N},
                                y::Real,
                                z::Dual{T,<:Any,N}) where {T,N}
    ex = Expr(:tuple, [:(fma(partials(x)[$i], y,  partials(z)[$i])) for i in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        v = fma(value(x), y, value(z))
        Dual{T}(v, $ex)
    end
end

@define_ternary_dual_op(
    Base.fma,
    calc_fma_xyz(x, y, z),                         # xyz_body
    calc_fma_xy(x, y, z),                          # xy_body
    calc_fma_xz(x, y, z),                          # xz_body
    Base.fma(y, x, z),                             # yz_body
    Dual{Tx}(fma(value(x), y, z), partials(x) * y), # x_body
    Base.fma(y, x, z),                              # y_body
    Dual{Tz}(fma(x, y, value(z)), partials(z))      # z_body
)

# muladd #
#--------#

@generated function calc_muladd_xyz(x::Dual{T,<:Any,N},
                                    y::Dual{T,<:Any,N},
                                    z::Dual{T,<:Any,N}) where {T,N}
    ex = Expr(:tuple, [:(muladd(value(x), partials(y)[$i], muladd(value(y), partials(x)[$i], partials(z)[$i]))) for i in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        v = muladd(value(x), value(y), value(z))
        return Dual{T}(v, $ex)
    end
end

@inline function calc_muladd_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
    vx, vy = value(x), value(y)
    result = muladd(vx, vy, z)
    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
end

@generated function calc_muladd_xz(x::Dual{T,<:Any,N},
                                   y::Real,
                                   z::Dual{T,<:Any,N}) where {T,N}
    ex = Expr(:tuple, [:(muladd(partials(x)[$i], y,  partials(z)[$i])) for i in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        v = muladd(value(x), y, value(z))
        Dual{T}(v, $ex)
    end
end

@define_ternary_dual_op(
    Base.muladd,
    calc_muladd_xyz(x, y, z),                         # xyz_body
    calc_muladd_xy(x, y, z),                          # xy_body
    calc_muladd_xz(x, y, z),                          # xz_body
    Base.muladd(y, x, z),                             # yz_body
    Dual{Tx}(muladd(value(x), y, z), partials(x) * y), # x_body
    Base.muladd(y, x, z),                             # y_body
    Dual{Tz}(muladd(x, y, value(z)), partials(z))      # z_body
)

# sincos #
#--------#

@inline sincos(x) = (sin(x), cos(x))

@inline function sincos(d::Dual{T}) where T
    sd, cd = sincos(value(d))
    return (Dual{T}(sd, cd * partials(d)), Dual{T}(cd, -sd * partials(d)))
end

# sincospi #
#----------#

if VERSION >= v"1.6.0-DEV.292"
    @inline function Base.sincospi(d::Dual{T}) where T
        sd, cd = sincospi(value(d))
        return (Dual{T}(sd, cd * π * partials(d)), Dual{T}(cd, -sd * π * partials(d)))
    end
end

# Symmetric eigvals #
#-------------------#

function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
    λ,Q = eigen(Symmetric(value.(parent(A))))
    parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N)
    Dual{Tg}.(λ, tuple.(parts...))
end

function LinearAlgebra.eigvals(A::Hermitian{<:Complex{<:Dual{Tg,T,N}}}) where {Tg,T<:Real,N}
    λ,Q = eigen(Hermitian(value.(real.(parent(A))) .+ im .* value.(imag.(parent(A)))))
    parts = ntuple(j -> diag(real.(Q' * (getindex.(partials.(real.(A)) .+ im .* partials.(imag.(A)), j)) * Q)), N)
    Dual{Tg}.(λ, tuple.(parts...))
end

function LinearAlgebra.eigvals(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
    λ,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev)))
    parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N)
    Dual{Tg}.(λ, tuple.(parts...))
end

# A ./ (λ - λ') but with diag special cased
function _lyap_div!(A, λ)
    for (j,μ) in enumerate(λ), (k,λ) in enumerate(λ)
        if k ≠ j
            A[k,j] /= μ - λ
        end
    end
    A
end

function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
    λ = eigvals(A)
    _,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev)))
    parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
    Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
end

function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
    λ = eigvals(A)
    _,Q = eigen(SymTridiagonal(value.(parent(A))))
    parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
    Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
end

###################
# Pretty Printing #
###################

function Base.show(io::IO, d::Dual{T,V,N}) where {T,V,N}
    print(io, "Dual{$(repr(T))}(", value(d))
    for i in 1:N
        print(io, ",", partials(d, i))
    end
    print(io, ")")
end

function Base.typemin(::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
    ForwardDiff.Dual{T,V,N}(typemin(V))
end

function Base.typemax(::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
    ForwardDiff.Dual{T,V,N}(typemax(V))
end

if VERSION >= v"1.6.0-rc1"
    Printf.tofloat(d::Dual) = Printf.tofloat(value(d))
end


text:
# [src/interpd.jl]

using Interpolations

"""
An interpolation function for time series which may have multiple values
at a single timestep (e.g., discrete updates). It works like MATLAB's
built-in interp1.

xi = interpd(t, x, ti);
xi = interpd(t, x, ti, mode);

Inputs:

    t     Times  (n-by-1)
    x     States (n-by-m)
    ti    Output times (p-by-1)
    mode  '-' for left-most values, '+' for right-most values (default)

    (NOTE that unlike the MATLAB version, no additional arguments can be passed in yet (e.g., 'nearest'))

Outputs:

    xi  States corresponding to output times (p-by-m)

Example:

    # Create some times and states. Note that there are two states at t=3.

    t = [2.76, 2.91, 3,   3,   3.12];
    x = [0.2,  0.3,  0.4, 1.1, 1.2];

    # Create the desired output times.

    ti = (2.8:0.1:3.1);

    # Interpolate (t, x) linearly at ti, keeping the right-most values.

    xi = interpd(t, x, ti)

    # Interpolate (t, x) linearly at ti, keeping the left-most values.

    xi = interpd(t, x, ti, '-')

See also: examples_odehybrid.

Online doc: http://www.anuncommonlab.com/doc/odehybrid/interpd.html

Copyright 2014 An Uncommon Lab
"""
function interpd(t, x, ti, mode)

    # Check that t and ti are columns.
    if (length(size(t)) > 1) && (~any(size(t) == 1))
        @error "INTERPD: Invalid inputs - input t should be a column!"
    else
        t = t[:];
    end

    if (length(size(ti)) > 1) && (~any(size(ti) == 1))
        @error "INTERPD: Invalid inputs - input ti should be a column!"
    else
        ti = ti[:];
    end

    # Check that the rows of x match the rows of t.
    if size(x, 1) != length(t) && size(x, 2) == length(t)
        x = x';
    end

    # Set the default mode.
    if isempty(mode)
        mode = '+';
    end

    # Create the output.
    xi = zeros(maximum(size(ti)), size(x, 2));

    # Find where there are doubled steps (t(doubled(k)) and t(doubled(k)+1) are the same.
    doubled = findall(t[1:end-1] .== t[2:end])

    # If nothing is doubled, just pass along to interp1.
    if isempty(doubled)
        xi = interp1((t,), x, ti);
        return xi
    end

    # We will break the interpolation into separate pieces -- those between the doubled
    #   (or tripled or otherwise duplicated) time steps. Then we'll interpolate that span.
    #   We'll keep the edges according to whether we're using - or + mode.

    # Make sure we get the first and last spans.
    if doubled[end] != maximum(size(t))
        doubled = [doubled; maximum(size(t))]
    end

    # Start with the initial point.
    if ((mode == "-") || (mode == '-')) && doubled[1] != 1
        xi[1, :] .= interp1((t[1:2],), x[1:2, :], ti[1])
    else
        # We may overwrite this for +
        xi[1, :] .= x[1, :]
    end

    # For each doubled index...
    k_last = 1;
    for k = doubled'

        # Get the span from the last index to this one.
        tk  = t[k_last:k];
        xk  = x[k_last:k, :];

        # Get the range for the requested outputs.
        out_indices = (ti .≥ t[k_last]) .& (ti .≤ t[k]);

        # If we're using - mode, drop the first index; we don't want to
        # overwrite the value from the end of the last span with the value
        # from the beginning of this one.

        if (mode == "-" || mode == '-')
            first_one = findall(out_indices .== 1)
            if length(first_one) > 0
                first_one = first_one[1]
                out_indices[first_one] = false
            end
        end

        # If there are any outputs in this span...
        if any(out_indices)

            # Select the subset of output times.
            tik = ti[out_indices];

            # If there's only one point, don't interp. Otherwise, use interp1.
            if k - k_last == 0
                xi[out_indices, :] .= xk;
            else
                xi[out_indices, :] .= interp1((tk,), xk, tik);
            end

        end

        # Get ready to move to the next span.
        k_last = k + 1;

    end

    if size(xi, 2) == 1
        xi = xi[:, 1]
    end

    return xi
end;

"""
Simplified Interp1 (Ported from MATLAB), which returns interpolated values of a 1-D function at query points using linear interpolation.

Inputs:
    x:  Vector containing sample points (in a tuple, pass in (x,) if needs be)
    v:  Vector containing values corresponding to x
    xq: Vector containing coordinates of the query points
    options: (Optional) Interpolation method

Outputs:
    interp_vals:  Interpolated values at query points
"""
function interp1(x, v, xq, options = Gridded(Linear()))

    squeeze = false

    if length(size(v)) == 2
        N, M = size(v)
    elseif length(size(v)) == 1
        N = size(v, 1)
        M = 1
        squeeze = true
    elseif length(size(v)) == 0
        N, M = 1, 1
    else
        @error "INTERP1: Invalid size for dependent variable v!"
    end

    if !isa(xq, Array)
        xq = [xq]
    end
    P = maximum(size(xq))
    interp_vals = zeros(P, M)

    # Allow v to be an [n × m] matrix
    if length(size(v)) > 1
        for i = 1:M
            interpolation = interpolate(x, v[:, i], options)
            interp_vals[:, i] .= interpolation[xq]
        end
    elseif size(v, 2) == 1
        interpolation = interpolate(x, v[:, 1], options)
        interp_vals .= interpolation[xq]
    else
        interpolation = interpolate(x, v, options)
        interp_vals .= interpolation[xq]
    end

    # Return correct array size (or scalar)
    if squeeze
        if size(interp_vals, 2) == 1
            interp_vals = interp_vals[:, 1]
        end
        if (size(interp_vals, 1) == 1)
            interp_vals = interp_vals[1, :]
        end
    end

    if maximum(size(interp_vals)) == 1
        interp_vals = interp_vals[1]
    end

    return interp_vals
end;

function interpd(t, x, ti)
    """
        Wrapper function for interpd that defaults "mode" to right-most ("+")
    """
    return interpd(t, x, ti, "+")
end;


text:
#=
 Methods and types associated with single and double layer potentials
=#

module ImmersedLayers

using Reexport
using DocStringExtensions
@reexport using CartesianGrids
@reexport using RigidBodyTools
#using CartesianGrids
#using RigidBodyTools

using LinearAlgebra

using UnPack

export DoubleLayer, SingleLayer, MaskType, Mask, ComplementaryMask,
        DoubleLayer!, SingleLayer!, Mask!,
        SurfaceCache,SurfaceScalarCache,SurfaceVectorCache,
        @ilmproblem,
        BasicScalarILMProblem,BasicVectorILMProblem,prob_cache,
        AbstractScalingType,GridScaling,IndexScaling,
        BasicILMCache,ILMSystem,
        similar_grid,similar_gridgrad,similar_gridcurl,similar_surface,
        zeros_grid,zeros_gridgrad,zeros_gridcurl,zeros_surface,
        ones_grid,ones_gridgrad,ones_gridcurl,ones_surface,
        x_grid,y_grid,x_gridcurl,y_gridcurl,
        AbstractExtraILMCache,AbstractScalarILMProblem,AbstractVectorILMProblem,
        ConvectiveDerivativeCache,convective_derivative, convective_derivative!,
        regularize!, interpolate!,
        regularize_normal!,normal_interpolate!,
        regularize_normal_cross!,normal_cross_interpolate!,
        surface_curl!,surface_divergence!,surface_grad!,inverse_laplacian!,
        laplacian!,
        surface_curl_cross!,surface_divergence_cross!,surface_grad_cross!,
        mask,mask!,complementary_mask,complementary_mask!,
        create_CLinvCT,create_CL2invCT,create_GLinvD,create_nRTRn,create_RTLinvR,
        create_GLinvD_cross,create_surface_filter,
        solve

abstract type LayerType{N} end

include("cartesian_extensions.jl")
include("tools.jl")
include("cache.jl")
include("problem.jl")
include("system.jl")
include("layers.jl")
include("grid_operators.jl")
include("surface_operators.jl")
include("matrix_operators.jl")

include("plot_recipes.jl")

@deprecate DoubleLayer surface_divergence!
@deprecate SingleLayer regularize!

end


text:
using DataFrames
showln(x) = (show(x); println())
# TODO: needs more links to docs.

# A DataFrame is an in-memory database
df = DataFrame(A = [1, 2], B = [e, pi], C = ["xx", "xy"])
showln(df)
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# The columns of a DataFrame can be indexed using numbers or names
showln(df[1])
#> [1,2]
showln(df[:A])
#> [1,2]

showln(df[2])
#> [2.718281828459045,3.141592653589793]
showln(df[:B])
#> [2.718281828459045,3.141592653589793]

showln(df[3])
#> ASCIIString["xx","xy"]
showln(df[:C])
#> ASCIIString["xx","xy"]

# The rows of a DataFrame can be indexed only by using numbers
showln(df[1, :])
#> 1x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
showln(df[1:2, :])
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# DataFrames can be loaded from CSV files using readtable()
iris = readtable("iris.csv")

# Check the names and element types of the columns of our new DataFrame
showln(names(iris))
#> [:SepalLength,:SepalWidth,:PetalLength,:PetalWidth,:Species]
showln(eltypes(iris))
#> Type[Float64,Float64,Float64,Float64,UTF8String]

# Subset the DataFrame to only include rows for one species
showln(iris[iris[:Species] .== "setosa", :])
#> 50x5 DataFrame
#> |-------|-------------|------------|-------------|------------|----------|
#> | Row # | SepalLength | SepalWidth | PetalLength | PetalWidth | Species  |
#> | 1     | 5.1         | 3.5        | 1.4         | 0.2        | "setosa" |
#> | 2     | 4.9         | 3.0        | 1.4         | 0.2        | "setosa" |
#> | 3     | 4.7         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 4     | 4.6         | 3.1        | 1.5         | 0.2        | "setosa" |
#> | 5     | 5.0         | 3.6        | 1.4         | 0.2        | "setosa" |
#> | 6     | 5.4         | 3.9        | 1.7         | 0.4        | "setosa" |
#> | 7     | 4.6         | 3.4        | 1.4         | 0.3        | "setosa" |
#> | 8     | 5.0         | 3.4        | 1.5         | 0.2        | "setosa" |
#> | 9     | 4.4         | 2.9        | 1.4         | 0.2        | "setosa" |
#> ⋮
#> | 41    | 5.0         | 3.5        | 1.3         | 0.3        | "setosa" |
#> | 42    | 4.5         | 2.3        | 1.3         | 0.3        | "setosa" |
#> | 43    | 4.4         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 44    | 5.0         | 3.5        | 1.6         | 0.6        | "setosa" |
#> | 45    | 5.1         | 3.8        | 1.9         | 0.4        | "setosa" |
#> | 46    | 4.8         | 3.0        | 1.4         | 0.3        | "setosa" |
#> | 47    | 5.1         | 3.8        | 1.6         | 0.2        | "setosa" |
#> | 48    | 4.6         | 3.2        | 1.4         | 0.2        | "setosa" |
#> | 49    | 5.3         | 3.7        | 1.5         | 0.2        | "setosa" |
#> | 50    | 5.0         | 3.3        | 1.4         | 0.2        | "setosa" |

# Count the number of rows for each species
showln(by(iris, :Species, df -> size(df, 1)))
#> 3x2 DataFrame
#> |-------|--------------|----|
#> | Row # | Species      | x1 |
#> | 1     | "setosa"     | 50 |
#> | 2     | "versicolor" | 50 |
#> | 3     | "virginica"  | 50 |

# Discretize entire columns at a time
iris[:SepalLength] = iround(iris[:SepalLength])
iris[:SepalWidth] = iround(iris[:SepalWidth])

# Tabulate data according to discretized columns to see "clusters"
tabulated = by(
    iris,
    [:Species, :SepalLength, :SepalWidth],
    df -> size(df, 1)
)
showln(tabulated)
#> 17x4 DataFrame
#> |-------|--------------|-------------|------------|----|
#> | Row # | Species      | SepalLength | SepalWidth | x1 |
#> | 1     | "setosa"     | 4           | 3          | 4  |
#> | 2     | "setosa"     | 5           | 2          | 1  |
#> | 3     | "setosa"     | 5           | 3          | 23 |
#> | 4     | "setosa"     | 5           | 4          | 17 |
#> | 5     | "setosa"     | 6           | 4          | 5  |
#> | 6     | "versicolor" | 5           | 2          | 3  |
#> | 7     | "versicolor" | 5           | 3          | 3  |
#> | 8     | "versicolor" | 6           | 2          | 6  |
#> | 9     | "versicolor" | 6           | 3          | 29 |
#> | 10    | "versicolor" | 7           | 3          | 9  |
#> | 11    | "virginica"  | 5           | 3          | 1  |
#> | 12    | "virginica"  | 6           | 2          | 1  |
#> | 13    | "virginica"  | 6           | 3          | 22 |
#> | 14    | "virginica"  | 7           | 3          | 19 |
#> | 15    | "virginica"  | 7           | 4          | 1  |
#> | 16    | "virginica"  | 8           | 3          | 4  |
#> | 17    | "virginica"  | 8           | 4          | 2  |

max_issues_repo_name: svaksha/JuliaByExample

max_forks_repo_name: svaksha/JuliaByExample


text:
using MultivariateSeries
X = @ring x0 x1 x2
d = 4
F = (x0+x1+x2)^d + 1.5*(x0+x1)^d -2.0*(x0-x2)^d
w, Xi = decompose(F)
F1 = tensor(w, Xi, X, d)
@assert norm(F-F1)<1.e-6


text:
module Util
export get_variable, get_param, build_constraint, get_previous_timestep,
get_loc_tech_carriers, get_loc_tech, get_timestep_weight, split_comma_list,
get_conversion_plus_io, loc_tech_is_in

using AxisArrays; using NCDatasets; using JuMP; using CalliopeJuMP; using Missings

function get_variable(dataset, var::String)
    """
    returns an AxisArray of the given variable, indexed over its relevant
    dimensions

    Parameters
    ----------
    var : string
        From the variables available in the loaded Dataset
    Returns
    -------
    AxisArray
    """
    var_dimnames = NCDatasets.dimnames(dataset[var])
    var_out = AxisArray(dataset[var][:],
        Tuple([AxisArrays.Axis{Symbol(i)}(dataset[i][:])
            for i in var_dimnames]))

    return var_out
end

function get_param(model_dict, param, non_time_indices, timestep=nothing)
    """
    returns an AxisArray of the given variable, indexed over all dimensions
    except timesteps, in order to remove timesteps from a possible timeseries
    variable (which may also be a static variable, we just don"t know!)

    Parameters
    ----------
    param : string
        From the variables available in the loaded Dataset
    non_time_indices : string or tuple
        List of all search indices that aren"t the timesteps. Should be ordered
        as reversed alphabetical (e.g. [loc_techs, costs])
    timestep : nothing or string, default = nothing
        one value from from sets["timesteps"]
    Returns
    -------
    AxisArray without the time axis
    """
    sets = model_dict["sets"]
    dataset = model_dict["dataset"]
    param_index = []

    if !haskey(dataset, param)
        return_val = model_dict["defaults"][param]
    else
        axes = NCDatasets.dimnames(dataset[param])
        for i in 1:length(non_time_indices)
            if "timesteps" in axes
                j = i + 1
            else
                j = i
            end
            if non_time_indices[i] in sets[axes[j]]
                push!(param_index, findin(sets[axes[j]], [non_time_indices[i]])[1])
            else
                return_val = model_dict["defaults"][param]
                break
            end
        end

        if length(param_index) == length(non_time_indices)
            if axes[1] == "timesteps"
                unshift!(param_index, findin(sets["timesteps"], [timestep])[1])
            end
            return_val = getindex(dataset[param], param_index...)
        end
    end

    if ismissing(return_val)
        return_val = model_dict["defaults"][param]
    end

    if return_val == nothing
        return 0
    else
        return return_val
    end
end

function build_constraint(constraint_sets, constraint_name, model_dict; tier=nothing)
    """
    create a constraint object, indexed over the given sets.

    Parameters
    ----------
    constraint_sets : array of strings
        names of the sets over which the constraint will be indexed
    constraint_name : string
        Name of constraint obect to build. Will also be used to refer to the
        correct constraint rule, which will be constraint_name + "_constraint_rule"
    tier : string
        For us in conversion_plus.jl `balance_conversion_plus_tiers_constraint_rule`
        whereby the same rule is used for multiple carrier tiers
    Return
    ------
    No direct return, the constraint Dictionary constraint_dict will contain an
    additional entry with the key `constraint_name`
    """

    sets = model_dict["sets"]
    backend_model = model_dict["backend_model"]

    # create a tuple of integer lists
    indexed_sets = [1:length(sets[i]) for i in constraint_sets]

    # Use the Julia @constraintref macro to build an empty constraint object
    temporary_constraint = Array{ConstraintRef}((length(i) for i in indexed_sets)...)
    # get the relevant function from this module
    temporary_constraint_rule = getfield(
        CalliopeJuMP.load_constraints,
        Symbol(string(constraint_name, "_constraint_rule"))
    )

    # Loop through to populate the previously built empty constraint object
    for i in CartesianRange(size(temporary_constraint))
        idx = [item for item in i.I]
        if tier != nothing
            temporary_constraint[i] = temporary_constraint_rule(backend_model,
                [sets[constraint_sets[constr]][idx[constr]]
                for constr in 1:length(constraint_sets)],
                model_dict, tier
            )
        else
            temporary_constraint[i] = temporary_constraint_rule(backend_model,
                [sets[constraint_sets[constr]][idx[constr]]
                for constr in 1:length(constraint_sets)],
                model_dict
            )
        end
    end

    return temporary_constraint
end

function get_previous_timestep(model_dict, timestep)
    timestep_idx = findfirst(model_dict["sets"]["timesteps"][DateTime(timestep)])

    return model_dict["sets"]["timesteps"][timestep_idx - 1]
end

function get_loc_tech_carriers(model_dict, loc_carrier)

    lookup = model_dict["parameters"]["lookup_loc_carriers"]
    sets = model_dict["sets"]
    loc_tech_carriers = split_comma_list(lookup[loc_carrier])

    loc_tech_carriers_prod = [
        i for i in loc_tech_carriers if i in sets["loc_tech_carriers_prod"]
    ]
    loc_tech_carriers_con = [
        i for i in loc_tech_carriers if i in sets["loc_tech_carriers_con"]
    ]

    if haskey(sets, "loc_tech_carriers_export")
        loc_tech_carriers_export = [
            i for i in loc_tech_carriers if i in sets["loc_tech_carriers_export"]
        ]
    else
        loc_tech_carriers_export = []
    end

    return (
        loc_tech_carriers_prod,
        loc_tech_carriers_con,
        loc_tech_carriers_export
    )
end

function get_loc_tech(loc_tech_carrier)
    return rsplit(loc_tech_carrier, "::", limit=2)[1]
end

function get_timestep_weight(model_dict)
    time_res = model_dict["parameters"]["timestep_resolution"].data
    weights = model_dict["parameters"]["timestep_weights"].data

    return sum(time_res .* weights) / 8760
end

function split_comma_list(comma_list)
    """
    Take a comma deliminated string and split it into a list of strings
    """
    return split(comma_list, ",")
end

function get_conversion_plus_io(model_dict, tier)
    """
    from a carrier_tier, return the primary tier (of `in`, `out`) and
    corresponding decision variable (`carrier_con` and `carrier_prod`, respectively)
    """
    if contains(tier, "out")
        return "out", model_dict["variables"]["carrier_prod"]
    elseif contains(tier, "in")
        return "in", model_dict["variables"]["carrier_con"]
    end
end

function loc_tech_is_in(sets, loc_tech, model_set)
    """
    Check if set exists and if loc_tech is in the set

    Parameters
    ----------
    loc_tech : string
    model_set : string
    """

    if haskey(sets, model_set) && loc_tech in sets[model_set]
        return true
    else
        return false
    end
end

end

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


text:
using NLopt
count = 0;
function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
         grad[1] = 0
         grad[2] = 0.5/sqrt(x[2])
    end
    global count;
    count::Int += 1;
    sqrt(x[2]);
end

function mycons(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 3*a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

opt = Opt(:LD_MMA, 2);
lower_bounds!(opt, [-Inf, 0.]);
xtol_rel!(opt,1e-4);
min_objective!(opt, myfunc);

inequality_constraint!(opt, (x,g) -> mycons(x,g,2,0), 1e-8);
inequality_constraint!(opt, (x,g) -> mycons(x,g,-1,1), 1e-8);

(minf,minx,ret) = optimize(opt, [1.2, 5.6]);
@printf "Ans => %6.3f at (%6.3f,%6.3f)" minf minx[1] minx[2];

max_stars_repo_path: Module 3/Chapter06/main-nlopt.jl

max_issues_repo_path: Module 3/Chapter06/main-nlopt.jl

max_forks_repo_path: Module 3/Chapter06/main-nlopt.jl


text:
include("tree.jl")
include("jats.jl")

path = ARGS[1]
tree = readjats(path)
open("sample-ljats.xml", "w") do f
    println(f, toxml(tree))
end

max_stars_repo_name: paperai/LightJATS

max_issues_repo_name: paperai/LightJATS

max_forks_repo_name: paperai/LightJATS


text:
@testset "Splines" begin
    @testset "Discontinuous knot set" begin
        t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3)
        B = BSpline(t)

        @test all(B[0.0, :] .== 0)
        @test all(B[0.0:0.5:1.0, :][1, :] .== 0)

        @test B.B == B[B.x,:]
        # Reference generated using AnalyticBSplines.jl
        @test B.S ≈ [ 3//5    2//9     2//45    0//1    0//1
                      2//9    7//15   83//270   1//270  0//1
                      2//45  83//270  26//27   83//270  2//45
                      0//1    1//270  83//270   7//15   2//9
                      0//1    0//1     2//45    2//9    2//5]
    end

    @testset "Cardinal splines" begin
        for k = 1:5
            t = ArbitraryKnotSet(k, 1.0:6.0, 1, 1)
            R = BSpline(t)
            @test size(R.B,2) == 6-k
            @test size(R.S) == (6-k,6-k)
            @test (bandwidth(R.S,1),bandwidth(R.S,2)) == (k-1,k-1)
        end
    end

    @testset "Evaluate B-splines" begin
        k = 3
        t = LinearKnotSet(k, 0, 1, 2)
        B = BSpline(t)
        @testset "Eval on subintervals" begin
            x₁ = range(-1,stop=-0.5,length=10)
            χ₁ = B[x₁, :]
            @test norm(χ₁) == 0

            function testbasis(x)
                χ = B[x, :]
                B̃ = spzeros(Float64, length(x), 4)
                B̃[:,1] = (x .>= 0) .* (x .< 0.5) .* ((2x .- 1).^2)
                B̃[:,2] = (x .>= 0) .* (x .< 0.5) .* (2/3*(1 .- (3x .- 1).^2)) +
                    (x .>= 0.5) .* (x .< 1)  .* (2*(x .- 1).^2)
                B̃[:,3] = (x .>= 0) .* (x .< 0.5) .* (2*x.^2) +
                    (x .>= 0.5) .* (x .< 1)  .* (2/3*(1 .- (3x .- 2).^2))
                B̃[:,4] = (x .>= 0.5) .* (x .<= 1) .* ((2x .- 1).^2)
                for j = 1:4
                    δ,δr = vecdist(χ[:,j],B̃[:,j])
                    @test δ < 1e-15
                    @test δr < 1e-15
                end
            end
            testbasis(range(0,stop=1,length=50))
            testbasis(range(-0.5,stop=0.6,length=40))
            testbasis(range(0.5,stop=1.6,length=40))
        end
    end

    @testset "Restricted basis" begin
        t = LinearKnotSet(3, 0.0, 1.0, 6)
        B = BSpline(t,3)
        B̃ = B[:,2:end-1]
        @test size(B̃) == (ℵ₁, 6)
        χ̃ = B̃[CompactBases.locs(B), :]
        @test size(χ̃,2) == 6
        @test χ̃ == B.B[:,2:end-1]
    end

    @testset "Evaluate spline" begin
        @testset "Linear knot set" begin
            k = 4
            t = LinearKnotSet(k, 0, 1, 5)
            B = BSpline(t)
            x = range(0, stop=1, length=301)
            χ = B[x, :]

            i = 1:size(B,2)
            c = sin.(i)

            # Whether we evaluate the spline or the basis functions should
            # not make a difference.
            s = (B*c)[x]
            s′ = χ*c
            @test all(isfinite, s)
            @test all(isfinite, s′)
            @test all(!isnan, s)
            @test all(!isnan, s′)
            @test s ≈ s′ rtol=1e-14
        end
        @testset "Discontinuous knot set" begin
            t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3)
            B = BSpline(t,3)

            χ = B[B.x, :]
            @test all(isfinite, χ)
            @test all(!isnan, χ)
        end
    end

    @testset "Mass matrices" begin
        k = 4
        t = LinearKnotSet(k, 0, 1, 5)
        B = BSpline(t)
        @test B'B == B.S
        @test_broken B.S isa Symmetric

        B̃ = B[:,2:end-1]
        @test B̃'B̃ == B.S[2:end-1,2:end-1]

        @test B̃'B == B.S[2:end-1,:]

        # TODO: Test mass matrices between different B-splines bases
        # but same knot set/quadrature points.
    end

    @testset "Function interpolation" begin
        B = BSpline(LinearKnotSet(7, 0, 1, 10))
        x = axes(B,1)
        c = B \ sin.(2π*x)
        x̃ = range(0, stop=1, length=301)
        χ = B[x̃,:]
        @test all(vecdist(χ*c, sin.(2π*x̃)) .< (3e-6, 2e-7))

        B̃ = B[:,2:end-1]
        c̃ = B̃ \ sin.(2π*x)
        χ̃ = B̃[x̃,:]
        @test all(vecdist(χ̃*c̃, sin.(2π*x̃)) .< (3e-6, 2e-7))
    end
end


text:
"""
   @relpath(name)

Return the absolute path to file `name` relative to the executing script.

### Input

- `name` -- file name

### Output

A string.

### Notes

This macro is equivalent to `joinpath(@__DIR__, name)`.

The `@relpath` macro is used in model scripts to load data files relative to the
location of the model, without having to change the directory of the Julia session.

For instance, suppose that the folder `/home/projects/models` contains the script
`my_model.jl`, and suppose that the data file `my_data.dat` located in the same
directory is required to be loaded by `my_model.jl`.
Then,

```julia
# suppose the working directory is /home/julia/ and so we ran the script as
# julia -e "include("../projects/models/my_model.jl")"

# in the model file /home/projects/models/my_model.jl we write:
d = open(@relpath "my_data.dat")
# do stuff with d
```
In this example, the macro `@relpath "my_data.dat"` evaluates to the string
`/home/projects/models/my_data.dat`. If the script `my_model.jl` only had
`d = open("my_data.dat")`, without `@relpath`, this command would fail as julia
would have looked for `my_data.dat` in the *working* directory, resulting in an
error that the file `/home/julia/my_data.dat` is not found.
"""
macro relpath(name::String)
    __source__.file === nothing && return nothing
    _dirname = dirname(String(__source__.file))
    dir = isempty(_dirname) ? pwd() : abspath(_dirname)
    return joinpath(dir, name)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


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

const AbstractSamples{T<:Real} = AbstractVector{T}


text:
function distance(s1::AbstractString, s2::AbstractString)

end


text:
function solve()
    x = parse(Int, readline())
    for a=-118:119, b=-118:119
        if a^5 - b^5 == x
            return a, b
        end
    end
end

a, b = solve()
println(a, " ", b)

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using LinearAlgebra: iszero
using PWS
using Distributions
using LinearAlgebra
using DataFrames
using CSVFiles

data = DataFrame(map(range(0,20,step=0.1)) do duration
    if iszero(duration)
        return (DiscreteTimes=duration, Value=0.0)
    end

    system = GaussianSystem(delta_t=0.05, duration=duration)

    n_dim = size(system.joint.Σ, 1) ÷ 2
    c_xx = Hermitian(system.joint.Σ[n_dim + 1:end, n_dim + 1:end])
    c_ss = Hermitian(system.joint.Σ[1:n_dim, 1:n_dim])
    c_z = Hermitian(system.joint.Σ)

    (DiscreteTimes=duration, Value=0.5 * (logdet(c_ss) + logdet(c_xx) - logdet(c_z)))
end)

save(joinpath(pwd(), "plots", "figure_gene_expr", "data", "gene-expression_2021-07-19_gaussian.csv"), data)

marginal = MvNormal(c_xx)
val1 = logpdf(marginal, initial[n_dim + 1:end]) # this is the analytically correct value of the log marginal


text:
function compute_crowding_distance(pop)
    crowding = zeros(length(pop))

    if length(pop[1].f) == 2
        sort!(pop, by = x -> x.f[1])
        pop[1].f[1] == pop[end].f[1] && return #Don't waste time if all indivs are the same
        crowding[1] = crowding[end] = Inf

        width_y1 = (pop[end].f[1] - pop[1].f[1])
        width_y2 = (pop[1].f[2] - pop[end].f[2])
        for i = 2:length(pop)-1
            crowding[i] = (pop[i+1].f[1] - pop[i-1].f[1]) / width_y1 + (pop[i-1].f[2] - pop[i+1].f[2]) / width_y2
        end

        return crowding
    end

    for j = 1:length(first(pop).f) # Foreach objective
        let j = j
            sort!(pop, by = x -> x.f[j]) #sort by the objective value
        end
        crowding[1] = crowding[end] = Inf #Assign infinite value to extremas
        if pop[1].f[j] != pop[end].f[j]
            for i = 2:length(pop)-1
                crowding[i] += (pop[i+1].f[j] - pop[i-1].f[j]) / (pop[end].f[j] - pop[1].f[j])
            end
        end
    end

    crowding
end

function update_crowding_distance!(pop)
    crowding = compute_crowding_distance(pop)
    for (i, sol) in enumerate(pop)
        sol.crowding = crowding[i]
    end

    pop
end

function truncate_population!(population, N)
    fast_non_dominated_sort!(population)

    let f::Int = 1
        ind = 0
        indnext = findlast(x -> x.rank == f, population)
        while 0 < indnext <= N
            ind = indnext
            f += 1
            indnext = findlast(x -> x.rank == f, population)
        end
        indnext == 0 && (indnext = length(population))
        update_crowding_distance!(view(population, ind+1:indnext))
        sort!(view(population, ind+1:indnext), by = x -> x.crowding, rev = true, alg = PartialQuickSort(N-ind))
    end
    deleteat!(population, N+1:length(population))
end


text:
using Spark

using Test

Spark.init()

@testset "Spark" begin

include("basic.jl")
include("map.jl")
include("map_partitions.jl")
include("attach.jl")
include("reduce.jl")
include("text_file.jl")
include("share_variable.jl")
include("flat_map.jl") # still fails
include("cartesian.jl")
include("group_by_key.jl")
include("reduce_by_key.jl")
include("collect_pair.jl")
include("map_pair.jl")
include("repartition_coalesce.jl")
include("filter.jl")
include("pipe.jl")
include("sql.jl")

end


text:
# This file contains code that was a part of Julia
# License is MIT: see LICENSE.md

size(cp::Chr) = ()
ndims(cp::Chr) = 0
ndims(::Type{<:Chr}) = 0
length(cp::Chr) = 1
lastindex(cp::Chr) = 1
getindex(cp::Chr) = cp
first(cp::Chr) = cp
last(cp::Chr) = cp

@static if !NEW_ITERATE
    start(cp::Chr) = false
    next(cp::Chr, state) = (cp, true)
    done(cp::Chr, state) = state
end

isempty(cp::Chr) = false
in(x::AbsChar, y::Chr) = x == y
in(x::Chr, y::AbsChar) = x == y
in(x::Chr, y::Chr) = x == y
-(x::AbsChar, y::Chr) = Int(x) - Int(y)
-(x::Chr, y::AbsChar) = Int(x) - Int(y)
-(x::Chr, y::Chr)     = Int(x) - Int(y)
-(x::T, y::Integer) where {T<:Chr} = T((Int32(x) - Int32(y))%UInt32)
+(x::T, y::Integer) where {T<:Chr} = T((Int32(x) + Int32(y))%UInt32)
+(x::Integer, y::Chr) = y + x
show(io::IO, cp::Chr)  = show(io, Char(cp))
print(io::IO, cp::Chr) = print(io, Char(cp))

Base.hex(chr::C)      where {C<:Chr} = hex(codepoint(chr))
Base.hex(chr::C, pad) where {C<:Chr} = hex(codepoint(chr), pad)


text:
module YLEFlow

using Symbolics

# Write your package code here.

@variables t k ρ a0 a1 a2 a3 a4 a5 a6 a7 a8
@variables U(..) af1(..) af2(..) af3(..) af4(..) af5(..) af6(..) af7(..) af8(..)

msigma2 = 2*ρ*Differential(ρ)(Differential(ρ)(U(ρ))) + Differential(ρ)(U(ρ))
#dUdk = exp(5*t)/(exp(2*t) + msigma2)
dUdk = k^5/(k^2 + msigma2)

S1 =  [U(ρ),af1(ρ),af2(ρ),af3(ρ),af4(ρ),af5(ρ),af6(ρ),af7(ρ),af8(ρ)]
S2 =  [af1(ρ),af2(ρ),af3(ρ),af4(ρ),af5(ρ),af6(ρ),af7(ρ),af8(ρ),0]
S3 = [a0, a1, a2, a3, a4, a5, a6, a7, a8]

R=[ρ,ρ,ρ,ρ,ρ,ρ,ρ,ρ,ρ]
R[1]=dUdk

Q=similar(R)
for (c,p) in zip(S1,S2)
    R[1]=simplify(substitute(R[1],Dict( Differential(ρ)(c)  => p)))
end

for i in 2:length(R)
    R[i] = Symbolics.derivative(R[i-1],ρ)

    for (c,p) in zip(S1,S2)
        R[i]=(substitute(R[i],Dict( Differential(ρ)(c)  => p)))
    end
end

Q = copy(R)

for (c,p) in zip(S1,S3)
    for i in 1:length(R)
        Q[i]=substitute(Q[i],Dict( c  => p))
    end
end

for i in 1:length(R)
    Q[i]=simplify(Q[i])
end

Symbolics.derivative(Q[2],a2)

dU1dk = Symbolics.build_function(Q[2],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU2dk = Symbolics.build_function(Q[3],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU3dk = Symbolics.build_function(Q[4],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU4dk = Symbolics.build_function(Q[5],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU5dk = Symbolics.build_function(Q[6],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU6dk = Symbolics.build_function(Q[7],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU7dk = Symbolics.build_function(Q[8],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})
dU8dk = Symbolics.build_function(Q[9],k,a1,a2,a3,a4,a5,a6,a7,a8,ρ,expression=Val{false})

function drodk(dU1dk,dU2dk,a2,a3,ρ)
    return  -(dU1dk + 2*ρ*dU2dk)/(3*a2+2*ρ*a3)
end

function flow!(da,a,m2,t)
 k = exp(t)
 ρ = a[1]
 a1 = m2 - 2ρ * a[2]

 dU1dk_s = dU1dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU2dk_s = dU2dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU3dk_s = dU3dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU4dk_s = dU4dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU5dk_s = dU5dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU6dk_s = dU6dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU7dk_s = dU7dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)
 dU8dk_s = dU8dk(t, a1, a[2], a[3], a[4], a[5], a[6], a[7], a[8], ρ)

 da[1] = drodk(dU1dk_s,dU2dk_s,a[2],a[3],ρ)

 da[2] = dU2dk_s + a[3] * da[1]
 da[3] = dU3dk_s + a[4] * da[1]
 da[4] = dU4dk_s + a[5] * da[1]
 da[5] = dU5dk_s + a[6] * da[1]
 da[6] = dU6dk_s + a[7] * da[1]
 da[7] = dU7dk_s + a[8] * da[1]
 da[8] = dU8dk_s
end

function g()

end

end


text:
using OffsetArrays

@testset "Associated Legendre Functions" begin
    nmax = 1_000
    min_length = (nmax * (nmax + 3)) ÷ 2 + 1

    for T in [BigFloat, Float64, Float32]
        ϵ = eps(T)
        P̄ = fill(T(NaN), min_length)
        P̄′ = Vector{T}(undef, min_length)
        recursion_coefficients = ALFRecursionCoefficients(nmax, T)

        for β in βrange(T)
            expiβ = cis(β)

            # Make sure we're checking the size of P̄ against nmax
            @test_throws ArgumentError ALFcompute!(P̄, expiβ, nmax+1, recursion_coefficients)

            ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)
            ALFcompute!(P̄′, expiβ, nmax)
            P̄″ = ALFcompute(expiβ, nmax)
            @test eltype(P̄″) === eltype(P̄)

            offset = -1
            for n in 0:nmax
                offset -= n
                P̄ₙ = OffsetVector(P̄, offset)
                # Test sum of squares
                total = sum(P̄ₙ[m]^2 for m in 0:n)
                @test total ≈ T(2n+1) rtol=2n*ϵ
                # Test in-place computation without recursion_coefficients
                P̄ₙ′ = OffsetVector(P̄′, offset)
                @test P̄ₙ′[0:n] ≈ P̄ₙ[0:n] rtol=2n*ϵ
                # Test allocating computation without recursion_coefficients
                P̄ₙ″ = OffsetVector(P̄″, offset)
                @test P̄ₙ″[0:n] ≈ P̄ₙ[0:n] rtol=2n*ϵ
            end
        end
    end
end


text:
using REPL.REPLCompletions
using JSON

const home = homedir()

function getcompletion(base::AbstractString, pos::Int, context_module=Main)
    if pos == 0
        return [], false
    end
    if pos >= 2 && base[1:2] in ("\"~", "`~")
        nbase = base[1] * home * base[3:end]
        base = nbase
        pos += (length(home)-1)
    end
    ret, range, should_complete = completions(base, pos, context_module)
    completionlist = unique(completion_text.(ret))
    maxsymbolpos = 0
    for i in ('`', '"', '/', '.', ' ')
        symbolpos = findprev(isequal(i), base, pos)
        if !isnothing(symbolpos)
            if symbolpos > maxsymbolpos
                if i == '.'
                    if maxsymbolpos == 0
                        maxsymbolpos = symbolpos
                    end
                else
                    maxsymbolpos = symbolpos
                end
            end
        end
    end
    if maxsymbolpos != 0
        completionlist =  base[1:maxsymbolpos] .* completionlist
    end
    return completionlist, should_complete
end

function evalstr(str::AbstractString, context_module=Main)
    result = nothing
    try
        result = context_module.eval(Meta.parse(str))
    catch error
        result = error
    end
    return result
end

function writecompletion(io, base::AbstractString, pos::Int, context_module=Main)
    completionlist, should_complete = getcompletion(base, pos, context_module)
    if should_complete
        JSON.Writer.print(io, completionlist)
    else
        JSON.Writer.print(io, [base])
    end
    nothing
end

function writecompletion(io, base::AbstractString, pos::AbstractString, context_module=Main)
    writecompletion(io, base, Meta.parse(pos), context_module)
end

function writeeval(io, str::AbstractString, context_module=Main)
    result = evalstr(str, context_module)
    if isnothing(result)
        JSON.Writer.print(io, ["Nothing"])
    elseif isa(result, Function)
        JSON.Writer.print(io, ["Function"])
    else
        JSON.Writer.print(io, [typeof(result), result])
    end
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

n = 10
plot(x=rand(n), y=rand(n), color=sinh.(rand(n)), Scale.color_asinh)


text:
function getNeighbors(arr, idx)
    i = idx[1]
    j = idx[2]
    maxrow = size(arr)[1]
    maxcol = size(arr)[2]
    neighboridxs = Tuple[]
    if i == 1 || i == maxrow || j == 1 || j == maxcol
        if (i,j) == (1,1)
            neighboridxs = push!(neighboridxs, (1,2))
            neighboridxs = push!(neighboridxs, (2,1))
            neighboridxs = push!(neighboridxs, (2,2))
        elseif (i,j) == (1, maxcol)
            neighboridxs = push!(neighboridxs, (1,maxcol-1))
            neighboridxs = push!(neighboridxs, (2,maxcol))
            neighboridxs = push!(neighboridxs, (2,maxcol-1))
        elseif (i,j) == (maxrow, 1)
            neighboridxs = push!(neighboridxs, (maxrow-1,1))
            neighboridxs = push!(neighboridxs, (maxrow,2))
            neighboridxs = push!(neighboridxs, (maxrow-1,2))
        elseif (i,j) == (maxrow, maxcol)
            neighboridxs = push!(neighboridxs, (maxrow-1,maxcol))
            neighboridxs = push!(neighboridxs, (maxrow,maxcol-1))
            neighboridxs = push!(neighboridxs, (maxrow-1, maxcol-1))
        elseif i == 1
            neighboridxs = push!(neighboridxs, (i, j-1))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
            neighboridxs = push!(neighboridxs, (i+1, j-1))
            neighboridxs = push!(neighboridxs, (i+1, j+1))
        elseif j == 1
            neighboridxs = push!(neighboridxs, (i-1, j))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
            neighboridxs = push!(neighboridxs, (i-1, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j+1))
        elseif j == maxcol
            neighboridxs = push!(neighboridxs, (i-1, maxcol))
            neighboridxs = push!(neighboridxs, (i+1, maxcol))
            neighboridxs = push!(neighboridxs, (i, maxcol - 1))
            neighboridxs = push!(neighboridxs, (i-1, j-1))
            neighboridxs = push!(neighboridxs, (i+1, j-1))
        elseif i == maxrow
            neighboridxs = push!(neighboridxs, (maxrow, j-1))
            neighboridxs = push!(neighboridxs, (maxrow, j+1))
            neighboridxs = push!(neighboridxs, (maxrow-1, j))
            neighboridxs = push!(neighboridxs, (maxrow-1, j-1))
            neighboridxs = push!(neighboridxs, (maxrow-1, j+1))
        end
    else
        neighboridxs = push!(neighboridxs, (i, j-1))
        neighboridxs = push!(neighboridxs, (i, j+1))
        neighboridxs = push!(neighboridxs, (i+1, j))
        neighboridxs = push!(neighboridxs, (i-1, j))
        neighboridxs = push!(neighboridxs, (i-1, j-1))
        neighboridxs = push!(neighboridxs, (i-1, j+1))
        neighboridxs = push!(neighboridxs, (i+1, j-1))
        neighboridxs = push!(neighboridxs, (i+1, j+1))
    end
    return neighboridxs
end

function step(arr)
    hasflashed = zeros(Int, size(arr))
    flashcandidates = Tuple[]
    for c in CartesianIndices(arr)
        arr[c] += 1
        if arr[c] > 9
            i = c[1]
            j = c[2]
            flashcandidates = push!(flashcandidates, (i,j))
        end
    end
    flashcounter = 0
    while (length(flashcandidates) > 0)
        for el in flashcandidates
            if arr[el[1], el[2]] > 9 && hasflashed[el[1], el[2]] == 0
                #println("el: "*string(el))
                flashcounter += 1
                #= if el == (10,10)
                    println(getNeighbors(arr, el))
                end =#
                flashcandidates = append!(flashcandidates, getNeighbors(arr, el))
                hasflashed[el[1], el[2]] = 1
                for idx in getNeighbors(arr, el)
                    #println(idx)
                    arr[idx[1], idx[2]] += 1
                end
                flashcandidates = deleteat!(flashcandidates, findfirst(x -> x == el, flashcandidates))
            elseif arr[el[1], el[2]] <= 9 || hasflashed[el[1], el[2]] != 0
                flashcandidates = deleteat!(flashcandidates, findfirst(x -> x == el, flashcandidates))
            end
        end
    end
    for c in CartesianIndices(arr)
        if arr[c] > 9
            arr[c] = 0
        end
    end
    return flashcounter
end

function solve(filename)
    values = readlines(filename*".txt")
    values = split.(values, "")
    #println(values)
    intvalues = Vector{Int}[]
    for line in values
        convline = parse.(Int, line)
        intvalues = push!(intvalues, convline)
    end
    valarr = zeros(Int, 10,10)
    for i in 1:length(intvalues)
        for j in 1:length(intvalues)
            valarr[i,j] = intvalues[i][j]
        end
    end
    #for i in 1:10
        #println(valarr[i, :])
    #end
    #println(" ")
    #println(getNeighbors(valarr, (1,3)))
    res = 0
    not_synced = true
    i = 1
    res_i = 0
    while not_synced
        stepres = step(valarr)
        res += stepres
        if stepres == 100
            not_synced = false
            res_i = i
        end
        i += 1
        #= if st == 10
            for i in 1:10
            println(valarr[i, :])
            end
            println(res)
        end  =#
    end
    println("Number of total flashes: "*string(res))
    println("synced after step: "*string(res_i))
end

solve("example")
solve("input")


text:
using Documenter
import BusinessDays

makedocs(
    sitename = "BusinessDays.jl",
    modules = [ BusinessDays ],
    pages = [ "Home" => "index.md",
              "API Reference" => "api.md" ]
)

deploydocs(
    repo = "github.com/JuliaFinance/BusinessDays.jl.git",
    target = "build",
)


text:
export Box

"""
Box set.
"""
struct Box{T <: Real} <: ConvexSet
    a::Vector{T}
    b::Vector{T}
    function Box(a::Vector{T}, b::Vector{T}) where T <: Real
        length(a) == length(b) || throw(ArgumentError("lengths of lower/upper bounds do not match"))
        new{T}(a, b)
    end
end

function project!(s::Box{T}, v::Vector{T}, y::Vector{T}) where T <: Real
    v .= clamp.(y, s.a, s.b)
end

function project(s::Box{T}, y::Vector{T}) where T <: Real
    project!(s, similar(y), y)
end


text:
using GeomDAE
using PyPlot

const Δt    = 0.1
const ntime = 100000
const neps  = 1E-14
const nmax  = 20

function f(t, x, fx)
    fx[1] = x[2]
    fx[2] = sin(x[1])
    nothing
end

function run_pendulum(tableau, filename)
    x0 = [acos(0.4), 0.0]
    ode = ODE(f, x0)
    int = Integrator(ode, tableau, Δt)
    sol = Solution(ode, Δt, ntime)

    integrate!(int, sol)
    set_initial_conditions!(sol, ode)

    print("Running ", tableau.name, "...")
    @time integrate!(int, sol)

    fig = figure(figsize=(6,6))
    plot(sol.x[1,:,1], sol.x[2,:,1])
    xlim(0, 6)
    ylim(-2,+2)
    savefig(filename)
end

run_pendulum(getTableauExplicitEuler(), "pendulum_explicit_euler.png")
run_pendulum(getTableauExplicitMidpoint(), "pendulum_explicit_midpoint.png")
run_pendulum(getTableauHeun(), "pendulum_heun.png")
run_pendulum(getTableauKutta(), "pendulum_kutta.png")
run_pendulum(getTableauERK4(), "pendulum_explicit_rk4.png")
run_pendulum(getTableauERK438(), "pendulum_explicit_rk438.png")

run_pendulum(getTableauImplicitEuler(), "pendulum_implicit_euler.png")
run_pendulum(getTableauGLRK1(), "pendulum_implicit_glrk1.png")
run_pendulum(getTableauGLRK2(), "pendulum_implicit_glrk2.png")
run_pendulum(getTableauGLRK3(), "pendulum_implicit_glrk3.png")

function qf(x, fx)
    fx[1] = x[1]
    nothing
end

function pf(t, x, fx)
    fx[1] = sin(x[1])
    nothing
end

function run_pendulum_partitioned(tableau, filename)
    q0 = [acos(0.4)]
    p0 = [0.0]

    ode = PODE(qf, pf, q0, p0)
    int = Integrator(ode, tableau, Δt)
    sol = Solution(ode, Δt, ntime)

    integrate!(int, sol)
    set_initial_conditions!(sol, ode)

    print("Running ", tableau.name, "...")
    @time integrate!(int, sol)

    fig = figure(figsize=(6,6))
    plot(sol.x[1,1,:,1], sol.x[2,1,:,1])
    # plot(sol[1], sol[2])
    xlim(0, 6)
    ylim(-2,+2)
    savefig(filename)
end

run_pendulum_partitioned(getTableauSymplecticEulerA(), "pendulum_symplectic_euler_a.png")
run_pendulum_partitioned(getTableauSymplecticEulerB(), "pendulum_symplectic_euler_b.png")
run_pendulum_partitioned(TableauEPRK(:PERK4, 4, getTableauERK4(), getTableauERK4()), "pendulum_explicit_prk4.png")


text:
# Autogenerated wrapper script for libLLVM_assert_jll for i686-linux-musl-cxx11
export libllvm, llvm_config

using Zlib_jll
JLLWrappers.@generate_wrapper_header("libLLVM_assert")
JLLWrappers.@declare_library_product(libllvm, "libLLVM-13jl.so")
JLLWrappers.@declare_executable_product(llvm_config)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libllvm,
        "lib/libLLVM-13jl.so",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        llvm_config,
        "tools/llvm-config",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# generate a random 2x4 matrix
# with values between 1 and 20
matrix = rand(1:20, 2, 4)

# determines the greates value within
# the previously generated matrix
max = maximum(matrix)

println("The randomly generated matrix is: ", matrix)
println("... and has the maximum: ", max)

max_stars_repo_name: DominicSchiller/osmi-module-datascience

max_issues_repo_name: DominicSchiller/osmi-module-datascience

max_forks_repo_name: DominicSchiller/osmi-module-datascience