text:
export AbstractPreprocessor,
       Chain,
       Preprocessor,
       FourierPreprocessor,
       PolynomialPreprocessor,
       TilingPreprocessor,
       RadialBasisFunctions,
       RandomProjection,
       SparseRandomProjection,
       ImageResize,
       ImageResizeNearestNeighbour,
       ImageResizeBilinear,
       ImageCrop,
       StackFrames

using .Utils: Tiling, encode
using LinearAlgebra: norm
using Flux: Chain
using ImageTransformations: imresize!
using ShiftedArrays:CircShiftedVector

"""
Preprocess an [`Observation`](@ref) and return a new observation.
By default, the preprocessor is only applied to the state field
of the [`Observation`](@ref) and other fields remain unchanged.

For customized preprocessors that inherit `AbstractPreprocessor`,
you can change this behavior by rewriting
`(p::AbstractPreprocessor)(obs::Observation)` method.
"""
abstract type AbstractPreprocessor end

(p::AbstractPreprocessor)(obs::Observation) =
    Observation(
        obs.reward,
        obs.terminal,
        p(obs.state),
        merge(obs.meta, (; Symbol(:state_before_, typeof(p).name) => obs.state)),
    )

"""
    FourierPreprocessor(order::Int)

Transform a scalar to a vector of `order+1` Fourier bases.
"""
struct FourierPreprocessor <: AbstractPreprocessor
    order::Int
end

(p::FourierPreprocessor)(s::Number) = [cos(i * π * s) for i = 0:p.order]

"""
    PolynomialPreprocessor(order::Int)

Transform a scalar to vector of maximum `order` polynomial.
"""
struct PolynomialPreprocessor <: AbstractPreprocessor
    order::Int
end

(p::PolynomialPreprocessor)(s::Number) = [s^i for i = 0:p.order]

"""
    TilingPreprocessor(tilings::Vector{<:Tiling})

Use each `tilings` to encode the state and return a vector.
"""
struct TilingPreprocessor{Tt<:Tiling} <: AbstractPreprocessor
    tilings::Vector{Tt}
end

(p::TilingPreprocessor)(s::Union{<:Number,<:Array}) = [encode(t, s) for t in p.tilings]

"""
    ImageResize(img::Array{T, N})
    ImageResize(dims::Int...) -> ImageResize(Float32, dims...)
    ImageResize(T::Type{<:Number}, dims::Int...)

Using BSpline method to resize the `state` field of an [`Observation`](@ref) to size of `img` (or `dims`).
"""
struct ImageResize{T, N} <: AbstractPreprocessor
    img::Array{T, N}
end

ImageResize(dims::Int...) = ImageResize(Float32, dims...)
ImageResize(T::Type{<:Number}, dims::Int...) = ImageResize(Array{T}(undef, dims))

function (p::ImageResize)(obs::Observation)
    imresize!(p.img, obs.state)
    Observation(obs.reward, obs.terminal, p.img, obs.meta)
end

"""
    struct ImageCrop
        xidx::UnitRange{Int64}
        yidx::UnitRange{Int64}
    end

Select indices `xidx` and `yidx` from a 2 or 3 dimensional array.
TODO: Inefficient!

# Example:

```
c = ImageCrop(2:5, 3:2:9)
c([10i + j for i in 1:10, j in 1:10])
```
"""
struct ImageCrop{Tx,Ty} <: AbstractPreprocessor
    xidx::Tx
    yidx::Ty
end

(c::ImageCrop)(x::Array{T,2}) where {T} = x[c.xidx, c.yidx]
(c::ImageCrop)(x::Array{T,3}) where {T} = x[:, c.xidx, c.yidx]

"""
    struct ImageResizeNearestNeighbour
        outdim::Tuple{Int64, Int64}
    end

Resize any image to `outdim = (width, height)` by nearest-neighbour
interpolation (i.e. subsampling).

# Example:

```
r = ImageResizeNearestNeighbour((50, 50))
r(rand(200, 200))
r(rand(UInt8, 3, 100, 100))
```
"""
struct ImageResizeNearestNeighbour <: AbstractPreprocessor
    outdim::Tuple{Int64,Int64}
end
function (r::ImageResizeNearestNeighbour)(x)
    indim = size(x)
    xidx = round.(Int64, collect(1:r.outdim[1]) .* indim[end-1] / r.outdim[1])
    yidx = round.(Int64, collect(1:r.outdim[2]) .* indim[end] / r.outdim[2])
    length(indim) > 2 ? x[:, xidx, yidx] : x[xidx, yidx]
end

struct SparseRandomProjection <: AbstractPreprocessor
    w::Array{Float64,2}
    b::Array{Float64,1}
end

(p::SparseRandomProjection)(s) = clamp.(p.w * s + p.b, 0, Inf)

struct RandomProjection <: AbstractPreprocessor
    w::Array{Float64,2}
end
(p::RandomProjection)(s) = p.w * s

struct RadialBasisFunctions <: AbstractPreprocessor
    means::Array{Array{Float64,1},1}
    sigmas::Array{Float64,1}
    state::Array{Float64,1}
end

struct Box{T}
    low::Array{T,1}
    high::Array{T,1}
end

function RadialBasisFunctions(box::Box, n, sigma)
    dim = length(box.low)
    means = [rand(dim) .* (box.high - box.low) .+ box.low for _ = 1:n]
    RadialBasisFunctions(means, typeof(sigma) <: Number ? fill(sigma, n) : sigma, zeros(n))
end

function (p::RadialBasisFunctions)(s)
    @inbounds for i = 1:length(p.state)
        p.state[i] = exp(-norm(s - p.means[i]) / p.sigmas[i])
    end
    p.state
end

"""
    StackFrames(::Type{T}=Float32, d::Int...)

Use a pre-initialized [`CircularArrayBuffer`](@ref) to store the latest several states specified by `d`.
Before processing any observation, the buffer is filled with `zero{T}`.
"""
mutable struct StackFrames{T, N} <: AbstractPreprocessor
    buffer::CircularArrayBuffer{T, N}
    StackFrames(d::Int...) = StackFrames(Float32, d...)
    function StackFrames(::Type{T}, d::Vararg{Int, N}) where {T, N}
        p = new{T, N}(CircularArrayBuffer{T}(d...))
        for _ in 1:capacity(p.buffer)
            push!(p.buffer, zeros(T, size(p.buffer)[1:N-1]))
        end
        p
    end
end

function (p::StackFrames{T, N})(obs::Observation) where {T, N}
    push!(p.buffer, obs.state)

    Observation(
        obs.reward,
        obs.terminal,
        p.buffer,
        obs.meta
    )
end


text:
# ---
# title: 1123. Lowest Common Ancestor of Deepest Leaves
# id: problem1123
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/lowest-common-ancestor-of-deepest-leaves/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the lowest common ancestor of its
# deepest leaves_.
#
# Recall that:
#
#   * The node of a binary tree is a leaf if and only if it has no children
#   * The depth of the root of the tree is `0`. if the depth of a node is `d`, the depth of each of its children is `d + 1`.
#   * The lowest common ancestor of a set `S` of nodes, is the node `A` with the largest depth such that every node in `S` is in the subtree with root `A`.
#
# **Note:** This question is the same as 865:
# <https://leetcode.comhttps://leetcode.com/problems/smallest-subtree-with-all-the-deepest-nodes/>
#
#
#
# **Example 1:**
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/07/01/sketch1.png)
#
#
#
#     Input: root = [3,5,1,6,2,0,8,null,null,7,4]
#     Output: [2,7,4]
#     Explanation: We return the node with value 2, colored in yellow in the diagram.
#     The nodes coloured in blue are the deepest leaf-nodes of the tree.
#     Note that nodes 6, 0, and 8 are also leaf nodes, but the depth of them is 2, but the depth of nodes 7 and 4 is 3.
#
# **Example 2:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#     Explanation: The root is the deepest node in the tree, and it's the lca of itself.
#
#
# **Example 3:**
#
#
#
#     Input: root = [0,1,3,null,2]
#     Output: [2]
#     Explanation: The deepest leaf node in the tree is 2, the lca of one node is itself.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree will be in the range `[1, 1000]`.
#   * `0 <= Node.val <= 1000`
#   * The values of the nodes in the tree are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
module Gradient

using ..Omega: RandVar, Ω, SimpleΩ, linearize, unlinearize, apl

export gradient, lineargradient, back!, GradAlg, value, grad, gradarray

# Gradient interfaces

abstract type GradAlg end

"`lineargradient(::RandVar, ω::Ω, ::Alg)` Returns as vector gradient of ω components"
function lineargradient end

"`back!(::RandVar, ω::Ω, ::TrackerGradAlg)` update ω s.t ω[id].grad is gradient"
function back! end

"`grad(rv, ω)` returns ω where gradients in components"
function grad end

"Value associated with dual/tracked/etc number/array"
function value end

# Default behaviour is identiy
value(x) = x

# AD implementations

include("old.jl") # Delete me!

include("tracker.jl")
export TrackerGrad

include("forwarddiff.jl")
export ForwardDiffGrad

# include("zygote.jl")
# export ZygoteGrad

end


text:
module ErrorfreeArithmetic

export SysFloat,
       # error-free transformations
       two_sum, two_diff, two_prod, two_square,
       three_sum, three_diff, three_prod,
       four_sum, four_diff,
       ad_minus_bc, two_fma, three_fma, two_muladd,
       # error-free transformations with magnitude sorted arguments
       two_hilo_sum, two_lohi_sum, two_hilo_diff, two_lohi_diff,
       three_hilo_sum, three_lohi_sum, three_hilo_diff, three_lohi_diff,
       four_hilo_sum, four_lohi_sum, four_hilo_diff, four_lohi_diff,
       # least-error transformations, as close to error-free as possible
       two_div,
       two_inv, two_sqrt,
       # error-free remainders
       div_rem, sqrt_rem

using Base: IEEEFloat

"""
    SysFloat

SysFloats are floating point types with processor fma support.
"""
const SysFloat = Union{Float64, Float32}

"""
    FloatWithFMA

Floats with FMA support fused multiply-add, fma(x,y,z)
""" FloatWithFMA

if isdefined(Main, :DoubleFloats)
    const FloatWithFMA = Union{Float64, Float32, Float16, Double64, Double32, Double16}
else
    const FloatWithFMA = Union{Float64, Float32, Float16}
end

include("errorfree_all.jl")
include("errorfree_tuples.jl")
include("leasterror.jl")
include("remainder.jl")
include("complex.jl")

end # ErrorfreeArithmetic


text:
struct Settings
    horizon::Base.RefValue{Int}
    use_forecast_data::Bool
    use_parameters::Base.RefValue{Bool}
    time_series_cache_size::Int
    warm_start::Base.RefValue{Bool}
    balance_slack_variables::Bool
    services_slack_variables::Bool
    initial_time::Base.RefValue{Dates.DateTime}
    PTDF::Union{Nothing, PSY.PTDF}
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes}
    optimizer_log_print::Bool
    constraint_duals::Vector{Symbol}
    system_to_file::Bool
    export_pwl_vars::Bool
    allow_fails::Bool
    ext::Dict{String, Any}
end

function Settings(
    sys;
    initial_time::Dates.DateTime = UNSET_INI_TIME,
    use_parameters::Bool = false,
    use_forecast_data::Bool = true,
    time_series_cache_size::Int = IS.TIME_SERIES_CACHE_SIZE_BYTES,
    warm_start::Bool = true,
    balance_slack_variables::Bool = false,
    services_slack_variables::Bool = false,
    horizon::Int = UNSET_HORIZON,
    PTDF::Union{Nothing, PSY.PTDF} = nothing,
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes} = nothing,
    optimizer_log_print::Bool = false,
    constraint_duals::Vector{Symbol} = Vector{Symbol}(),
    system_to_file = true,
    export_pwl_vars = false,
    allow_fails = false,
    ext = Dict{String, Any}(),
)
    if time_series_cache_size > 0 &&
       sys.data.time_series_storage isa IS.InMemoryTimeSeriesStorage
        @info "Overriding time_series_cache_size because time series is stored in memory"
        time_series_cache_size = 0
    end

    return Settings(
        Ref(horizon),
        use_forecast_data,
        Ref(use_parameters),
        time_series_cache_size,
        Ref(warm_start),
        balance_slack_variables,
        services_slack_variables,
        Ref(initial_time),
        PTDF,
        optimizer,
        optimizer_log_print,
        constraint_duals,
        system_to_file,
        export_pwl_vars,
        allow_fails,
        ext,
    )
end

function log_values(settings::Settings)
    text = Vector{String}()
    for (name, type) in zip(fieldnames(Settings), fieldtypes(Settings))
        val = getfield(settings, name)
        if type <: Base.RefValue
            val = val[]
        end
        push!(text, "$name = $val")
    end

    @debug "Settings: $(join(text, ", "))"
end

function copy_for_serialization(settings::Settings)
    vals = []
    for name in fieldnames(Settings)
        if name == :optimizer
            # Cannot guarantee that the optimizer can be serialized.
            val = nothing
        else
            val = getfield(settings, name)
        end

        push!(vals, val)
    end

    return deepcopy(Settings(vals...))
end

function restore_from_copy(
    settings::Settings;
    optimizer::Union{Nothing, JuMP.MOI.OptimizerWithAttributes},
)
    vals = Dict{Symbol, Any}()
    for name in fieldnames(Settings)
        if name == :optimizer
            vals[name] = optimizer
        elseif name == :ext
            continue
        else
            val = getfield(settings, name)
            vals[name] = isa(val, Base.RefValue) ? val[] : val
        end
    end

    return vals
end

get_horizon(settings::Settings) = settings.horizon[]
get_use_forecast_data(settings::Settings) = settings.use_forecast_data
get_use_parameters(settings::Settings) = settings.use_parameters[]
get_initial_time(settings::Settings)::Dates.DateTime = settings.initial_time[]
get_PTDF(settings::Settings) = settings.PTDF
get_optimizer(settings::Settings) = settings.optimizer
get_ext(settings::Settings) = settings.ext
get_warm_start(settings::Settings) = settings.warm_start[]
get_constraint_duals(settings::Settings) = settings.constraint_duals
get_balance_slack_variables(settings::Settings) = settings.balance_slack_variables
get_services_slack_variables(settings::Settings) = settings.services_slack_variables
get_system_to_file(settings::Settings) = settings.system_to_file
get_export_pwl_vars(settings::Settings) = settings.export_pwl_vars
get_allow_fails(settings::Settings) = settings.allow_fails
get_optimizer_log_print(settings::Settings) = settings.optimizer_log_print
use_time_series_cache(settings::Settings) = settings.time_series_cache_size > 0

function set_horizon!(settings::Settings, horizon::Int)
    settings.horizon[] = horizon
    return
end

function set_initial_time!(settings::Settings, initial_time::Dates.DateTime)
    settings.initial_time[] = initial_time
    return
end

function set_use_parameters!(settings::Settings, val::Bool)
    settings.use_parameters[] = val
    return
end

function set_warm_start!(settings::Settings, warm_start::Bool)
    settings.warm_start[] = warm_start
    return
end


text:
using PreSumm
using CUDA
using GPUArrays:allowscalar
using Flux
using Transformers
using Transformers.Basic
using Transformers.Pretrain
using BSON: @load, @save

if CUDA.functional(true)
    # CUDA.device!(1)
    allowscalar(false)
    Flux.use_cuda[] = true
    enable_gpu(true)
else
    @warn "You're training the model without GPU support."
    Flux.use_cuda[] = false
    enable_gpu(false)
end

snapshots = readdir(out_dir(), join=true)

model_snapshots = filter(
    file -> isfile(file) && occursin("model.bson", file),
    snapshots
)

weights_snapshots = filter(
    file -> isfile(file) && occursin("weights.bson", file),
    snapshots
)

new_weights_snapshots = filter(
    file -> !(file ∈ weights_snapshots),
    map(
        file -> replace(file, "model.bson" => "weights.bson"),
        model_snapshots
    )
)

for weights_snapshot ∈ reverse(new_weights_snapshots)
    model_snapshot = replace(weights_snapshot, "weights.bson" => "model.bson")
    @assert isfile(model_snapshot)
    @load model_snapshot model
    @show model
    model = model |> cpu
    @show model
    @show params(model)
end


text:
"""
    Logger(logdir; [overwrite=false, time=time(), purge_step::Int])

Creates a TensorBoardLogger in the folder `logdir`. If `overwrite` is true it
deletes previously created events.
If `purge_step::Int` is passed, every step before `purge_step` will be ignored
by tensorboard (usefull in the case of restarting a crashed computation).
"""
function Logger(logdir; overwrite=false, time=time(), purge_step::Union{Int,Nothing}=nothing)
    if overwrite
        rm(logdir; force=true, recursive=true)
    end
    mkpath(logdir)

    hostname = gethostname()
    fname    = "events.out.tfevents.$time.$hostname"
    fpath    = joinpath(logdir, fname)
    file     = open(fpath, "w")

    all_files = Dict(fpath => file)

    # Create the initial log
    if purge_step != nothing
        start_step = purge_step
        ev_0 = Event(wall_time=time, step=start_step, file_version="brain.Event:2")
        write_event(file, ev_0)
        sess_log = TensorBoardLogger.SessionLog(status=TensorBoardLogger.SessionLog_SessionStatus.START)
        ev_0 = Event(wall_time=time, step=start_step, session_log=sess_log)
        write_event(file, ev_0)
    else
        start_step = 0
        ev_0 = Event(wall_time=time, step=start_step, file_version="brain.Event:2")
        write_event(file, ev_0)
    end

    Logger(realpath(logdir), file, all_files, start_step)
end

# normally the logs don't overwrite, but if you've not given a path, you clearly don't care.
Logger() = Logger("tensorboard_logs", overwrite=true)

# Accessors
"""
    logdir(lg::Logger) -> String

Returns the directory to which Logger `lg` is writing data.
"""
logdir(lg::Logger)   = lg.logdir

"""
    get_file(lg::Logger) -> IOS

Returns the main `file` IOStream object of Logger `lg`.
"""
get_file(lg::Logger) = lg.file
function add_log_file(lg::Logger, path::String)
    file = open(path, "w")
    lg.all_files[path] = file
    return file
end

"""
    get_file(lg, tags::String...) -> IOS

Returns the `file` IOStream object of Logger `lg` writing to the tag
`tags1/tags2.../tagsN`.
"""
function get_file(lg::Logger, tags::String...)
    key = joinpath(logdir(lg), tags...)
    if key ∈ lg.all_files
        return lg.all_files[key]
    else
        return add_log_file(lg, key)
    end
end

"""
    set_step(lg, iter::Int)

Sets the iteration counter in the logger to `iter`. This counter is used by the
logger when no value is passed by the user.
"""
set_step(lg::Logger, iter::Int) = lg.global_step = iter

"""
    iteration(lg)

Returns the internal iteration counter of the logger. When no step keyword
is provided to the loggers, it will use this value.
"""
step(lg::Logger) = lg.global_step

# Additional things

#const default_logging_session = Ref(Logger())

"""
    set_tb_logdir(logdir, overwrite=false)
Start a new log in the given directory
"""
function set_tb_logdir(logdir, overwrite=false)
    default_logging_session[] = Logger(logdir, overwrite=overwrite)
end

"""
    reset_tb_logs()
Reset the current log, deleteing all information
"""
function reset_tb_logs()
    logdir = default_logging_session[].logdir
    default_logging_session[] = Logger(logdir, overwrite=true)
end


text:
using CSV; using DataFrames; using DecisionTree; using Statistics
# load data
pwd()
readdir()
submission = CSV.read("house_prices/sample_submission.csv", DataFrame)
test = CSV.read("house_prices/test.csv", DataFrame)
train = CSV.read("house_prices/train.csv", DataFrame)

# y_train variable
y_train = train.SalePrice

# removing ID and train.SalePrice
select!(train, Not(:"Id"))
train = select!(train, Not(:"SalePrice"))
test = select!(test, Not(:"Id"))

# concatenating train and test data
all_data = vcat(train, test)

# casting variables
for (i, col) in enumerate(eachcol(all_data))
    if eltype(col) == String7
        all_data[!, i] = String.(col)
    elseif eltype(col) == String3
        all_data[!, i] = String.(col)
    elseif eltype(col) == String1
        all_data[!, i] = String.(col)
    end
end

# checking types:
eltype.(eachcol(all_data))

# separate string cols and int cols:
int_cols = select(all_data, findall(col -> eltype(col) <: Int64, eachcol(all_data)))
string_cols = select(all_data, findall(col -> eltype(col) <: String, eachcol(all_data)))

# missing values
for col in eachcol(int_cols)
    replace!(col, missing => mean(col))
end

for col in eachcol(string_cols)
    replace!(col, missing => "missing_value")
end

# col bind:
all_data = hcat(string_cols, int_cols)

# run model
model = RandomForestRegressor(n_trees = 500)

# separating columns in X_train and X_test
X_train = all_data[1:nrow(train), :]
X_test = all_data[nrow(train)+1:end, :]

# conversions to fit the model
y_train = Int.(y_train)
y_train = convert(AbstractVector, y_train)
X_train = Matrix(X_train)

fit!(model, X_train, y_train)

# make predictions
X_test = Matrix(X_test)
y_hat = predict(model, X_test)

#submission file
submission.SalePrice = y_hat
CSV.write("./house_prices/submission_julia.csv", submission)


text:
#####
##### `map`
#####

function rrule(::typeof(map), f, xs...)
    y = map(f, xs...)
    function map_pullback(ȳ)
        ntuple(length(xs)+2) do full_i
            full_i == 1 && return NO_FIELDS
            full_i == 2 && return DoesNotExist()
            i = full_i-2
            @thunk map(ȳ, xs...) do ȳi, xis...
                _, pullback = _checked_rrule(f, xis...)
                ∂xis = pullback(ȳi)
                extern(∂xis[i+1])  #+1 to skp ∂self
            end
        end
    end
    return y, map_pullback
end

#####
##### `mapreduce`, `mapfoldl`, `mapfoldr`
#####

for mf in (:mapreduce, :mapfoldl, :mapfoldr)
    sig = :(rrule(::typeof($mf), f, op, x::AbstractArray{<:Real}))
    call = :($mf(f, op, x))
    if mf === :mapreduce
        insert!(sig.args, 2, Expr(:parameters, Expr(:kw, :dims, :(:))))
        insert!(call.args, 2, Expr(:parameters, Expr(:kw, :dims, :dims)))
    end
    pullback_name = Symbol(mf, :_pullback)
    body = quote
        y = $call
        function $pullback_name(ȳ)
            ∂x = @thunk broadcast(x, ȳ) do xi, ȳi
                _, pullback_f = _checked_rrule(f, xi)
                _, ∂xi = pullback_f(ȳi)
                extern(∂xi)
            end
            (NO_FIELDS, DoesNotExist(), DoesNotExist(), ∂x)
        end
        return y, $pullback_name
    end
    eval(Expr(:function, sig, body))
end

#####
##### `sum`
#####

function frule((_, ẋ), ::typeof(sum), x)
    return sum(x), sum(ẋ)
end

function rrule(::typeof(sum), f, x::AbstractArray{<:Real}; dims=:)
    y, mr_pullback = rrule(mapreduce, f, Base.add_sum, x; dims=dims)
    function sum_pullback(ȳ)
        return NO_FIELDS, DoesNotExist(), last(mr_pullback(ȳ))
    end
    return y, sum_pullback
end

function rrule(::typeof(sum), x::AbstractArray{<:Real}; dims=:)
    y, inner_pullback = rrule(sum, identity, x; dims=dims)
    function sum_pullback(ȳ)
        return NO_FIELDS, last(inner_pullback(ȳ))
    end
    return y, sum_pullback
end

function rrule(::typeof(sum), ::typeof(abs2), x::AbstractArray{<:Real}; dims=:)
    y = sum(abs2, x; dims=dims)
    function sum_abs2_pullback(ȳ)
        return (NO_FIELDS, DoesNotExist(), @thunk(2ȳ .* x))
    end
    return y, sum_abs2_pullback
end


text:
include("MLP.jl")
include("MLPData.jl")
include("FFMLP.jl")
include("RecMLP.jl")

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using TypeClasses
using Documenter

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

makedocs(;
    modules=[TypeClasses],
    authors="Stephan Sahm <stephan.sahm@gmx.de> and contributors",
    repo="https://github.com/JuliaFunctional/TypeClasses.jl/blob/{commit}{path}#{line}",
    sitename="TypeClasses.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaFunctional.github.io/TypeClasses.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
            "Introduction" => "manual.md",
            "TypeClasses" => "manual-TypeClasses.md",
            "DataTypes" => "manual-DataTypes.md",
        ],
        "Library" => "library.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaFunctional/TypeClasses.jl",
    devbranch="main",
)


text:
function phase_space_summary_plot(sim, dir; species="electron", save=false)
    r_mean, p_mean,
    r_mean_plus, p_mean_plus,
    r_mean_minus, p_mean_minus = @withprogress name="Computing phase space mean" begin
        r_mean, p_mean = phase_space_mean(sim, dir; species)
        @logprogress 1/3
        r_mean_plus, p_mean_plus = phase_space_mean(sim, dir; species, rcond=(i,r)->r>r_mean[i])
        @logprogress 2/3
        r_mean_minus, p_mean_minus = phase_space_mean(sim, dir; species, rcond=(i,r)->r<r_mean[i])
        @logprogress 1

        r_mean, p_mean, r_mean_plus, p_mean_plus, r_mean_minus, p_mean_minus
    end

    ts = ustrip.(unit_t, get_time.(sim))
    r = uconvert.(unit_l, r_mean)
    r_plus = uconvert.(unit_l, r_mean_plus)
    r_minus = uconvert.(unit_l, r_mean_minus)
    rm, rmp, rmm = mean(r), mean(r_plus), mean(r_minus)
    p = uconvert.(pₑ, p_mean)
    p_plus = uconvert.(pₑ, p_mean_plus)
    p_minus = uconvert.(pₑ, p_mean_minus)
    pm, pmp, pmm = mean(p), mean(p_plus), mean(p_minus)

    r_label = "$dir - $(dir)ₘ"
    plt1 = Plots.plot(ts, r .- rm, xlabel="t", ylabel=r_label, label="all")
    Plots.plot!(plt1, ts, r_plus .- rmp, xlabel="t", ylabel=r_label, label="r > r̅")
    Plots.plot!(plt1, ts, r_minus .- rmm, xlabel="t", ylabel=r_label, label="r < r̅")

    plt2 = Plots.plot(ts, p, xlabel="t", ylabel="p$dir", label="all")
    Plots.plot!(plt2, ts, p_plus, xlabel="t", ylabel="p$dir", label="r > r̅")
    Plots.plot!(plt2, ts, p_minus, xlabel="t", ylabel="p$dir", label="r < r̅")

    plt3 = Plots.scatter(r, p, zcolor=ts, xlabel=string(dir), ylabel="p$dir", legend=false)

    plt4 = Plots.scatter(r .- rm, p, xlabel=r_label, ylabel="p$dir", label="all")
    Plots.scatter!(r_plus .- rmp, p_plus, xlabel=r_label, ylabel="p$dir", label="r > r̅")
    Plots.scatter!(r_minus .- rmm, p_minus, xlabel=r_label, ylabel="p$dir", label="r < r̅")

    plt = Plots.plot(plt1, plt2, plt3, plt4, layout=4, size=(900,700))

    save && Plots.savefig(plt, joinpath(sim.dir, "phase_space_summary.png"))

    return plt
end


text:
#' # Optimizing water networks with Knitro
#'
#' This notebook explains how to solve effectively NLP and MINLP
#' problem with the Knitro solver.
#'
#' The problem studied hereafter was originally designed by
#' Pierre Carpentier, for educational purposes.
#' The original code is available [here](http://perso.ensta-paristech.fr/~pcarpent/TP_Reseau/ENSMP/).
#' The author credits him for the following tutorial.
#'
#' # Description of the problem
#'
#' The tutorial studies the optimization of water flows in a water
#' distribution network (at steady state). The optimizer aims at
#' minimizing the network's energy while statisfying some linear constraints
#' corresponding to Kirchhoff's first law.
#'
#' Let $G = (N, E)$ be an oriented graph. We note $n = |E|$ the number of
#' arcs and $m = |N|$ the number of nodes in this graph.
#'
#' We suppose that the network has $r$ tanks storing some waters to fulfill
#' the demands in $d$ nodes (distinct from the tanks).
#' We split the set of nodes $N$ accordingly:
#' $$
#' N = N_r \cup N_d.
#' $$
#'
#' We suppose further that the graph is *connected*, implying
#' $$
#' n \geq m - 1 .
#' $$
#' We note $A \in R^{m \times n}$ the incidence matrix of the graph.
#'
#' Lets introduce physical variables to describe further the water network.
#' - the vector of resulting flows at nodes, denoted $f = (f_d, f_r)$. $f_d$ is given for demands' nodes.
#' - the vector of pressures at nodes, denoted $p = (p_d, p_r)$. Pressures $p_r$ are given for the reservoirs.
#' - the vector of resistances in all arcs, denoted $r$ (parameters of the problem)
#' - the vector of flows across arcs, denoted $q$.
#'
#' The decision variable is the vector of flows $q$.
#'
#' ### Constraints
#'
#' The first Kirchhoff law states that:
#' $$
#' A q -f = 0,
#' $$
#' as we suppose that no accumulation occurs in the nodes.
#' The second Kirchhoff law takes into account the losses in the pipes,
#' which is given on each arc by a function $\phi_{alpha}$ (corresponding
#' to the Colebrooks law):
#' $$
#' \phi_{\alpha}(q_\alpha) = r_\alpha q_\alpha | q_\alpha |
#' $$
#' The second Kirchhoff law writes, in a vectorial manner,
#' $$
#' A^\top p + r \circ q \circ | q | = 0
#' $$
#'
#' ### Objective
#' On each arc $\alpha$, we define the energy function $\Phi_\alpha$
#' as
#' $$
#' \Phi_\alpha(q_\alpha) = \dfrac{1}{3} r_\alpha q_\alpha^2 | q_\alpha |
#' $$
#' On the graph, the overall energy equates
#' $$
#' J(q, f_r) = \dfrac{1}{3} q^\top (r \circ q \circ | q |) + p_r^\top f_r
#' $$
#'
#' The global problem writes
#' $$
#' \min_{q} \dfrac{1}{3} q^\top (r \circ q \circ | q |) + p_r^\top f_r \qquad
#' s.t. \quad Aq  - f = 0
#' $$
#'
#' ### Reformulation
#' By applying some mathematical tricks, we are able to reformulate
#' the problem in the following manner.
#'
#' By considering some properties of the overall graph, we split the vector
#' of flows in two $q = (q_T, q_C)$, where $q_T$ depends linearly on $q_C$.
#' Then, we are able to prove the existence of a matrix $B$ and a fixed
#' vector $q_0$ such that the vector of flows on arcs writes
#' $$
#' q = q_0 + B \; q_C
#' $$
#' By using this formulation, we can reduce the dimension of the search
#' space by selecting as decision variable the subvector $q_C$ instead of $q$.
#'
#' Note that the problem becomes also constraint-free, as the vector $q$
#' statisfying the previous equation sastifies also the first Kirchhoff law.
#'
#' We reformulate the optimization problem as
#' $$
#' \min_{q_c} \dfrac{1}{3} q^\top (r \circ q \circ | q |) + p_r^\top f_r \qquad
#' s.t. \quad q = q_0 + B q_C
#' $$

#' ---
#' # Numerical resolution
#' We start by importing JuMP and Knitro.

using JuMP, KNITRO

#' Import data.
include("data.jl");

#' Plotting utilities
using Pkg

# Turned off by default.
PLOT_GRAPH = false
if haskey(Pkg.installed(), "PyPlot")
    PLOT_GRAPH = true
    include("utils.jl");
end

#' ---
#' ## First case study: the non-linear problem
#'
#' Define Knitro Optimizer with JuMP.
model = Model(with_optimizer(KNITRO.Optimizer, outlev=3))

#' Write non-linear optimization problem first.
# Parameters.
α1 = Ar' * pr
# Dimension of the problem.
nx = n - md
@variable(model, qc[1:nx])
# Add dummy variable in the model...
@variable(model, q[1:n])
# ... as we restrict the problem inside a 9-dimensional manifold:
@constraint(model, q .== q0 + B*qc)
@NLobjective(model, Min,
             sum(r[i] * abs(q[i]) * q[i]^2 / 3 + α1[i] * q[i] for i in 1:n))
optimize!(model)

#' Display results
if PLOT_GRAPH
    optimal_flow = JuMP.value.(q)
    fig = figure()
    plot_network(flow=optimal_flow)
    display(fig)
end

#' Observe that some arcs seem almost useless here...

#' ---
#' ## Extension to Mixed-integer non-linear programming

#' We were able to quickly formule the NLP problem with JuMP, and Knitro
#' finds a solution in few iterations.
#' We now modify slightly the problem. Consider the $n$ arcs inside the graph.
#' The questions are:
# - *how many arcs in $E$ can we remove before the problem becomes infeasible?*
# - *how much is the cost impacted by arc removals?*

#' We define hereafter the MINLP version of the problem.
# Add a maximum flows through the pipes
const QMAX = 10.

function load_mip_model!(model::JuMP.Model; nremovals=3)
    # Reconsider the NLP problem introduced previously.
    @variable(model, qc[1:nx])
    @variable(model, q[1:n])
    @constraint(model, q .== q0 + B*qc)

    # Introduce switch: $z$ is set to 0 if the corresponding arc is removed
    # from the graph.
    @variable(model, z[1:n], Bin)

    # Bounds the abs flows through the arcs by QMAX
    # The constraints write:
    # | q | <= q_{max} * z
    # Note that z_a = 0 implies that q_a = 0.
    @constraint(model,  q .<= QMAX * z)
    @constraint(model, -q .<= QMAX * z)
    # Ensure that we remove exactly nremovals arcs.
    @constraint(model, sum(z) == n - nremovals)

    # Same cost as previously.
    @NLobjective(model, Min,
                 sum(r[i] * abs(q[i]) * q[i]^2 / 3 + α1[i] * q[i] for i in 1:n))
    return
end

#' **Optimization related remark:** Note that the constraint
#' $$
#' | q | \leq q_{max} \; z
#' $$
#' could be reformulated as a complementarity constraint (also supported by Knitro
#' but not natively by JuMP).

#' **Physics related remark:** We are a bit loosy with the physics of the
#' problem if we remove the arcs in the previous manner. Imagine that for
#' a given arc $a \in E$, $z_a =0$ thus implying $q_a =0$. Then, if we
#' note $i$ and $j$ the two adjacent nodes, the physics tells us that
#' $$
#' p_i = 0 \, \quad p_j = 0 ,
#' $$
#' which is not the case if we solve the previous optimization problem.
#' However, the goal of this tutorial is purely pedagogical and as a
#' consequence we allow us to play a bit with the physics.

#' ### Default Knitro

#' Solve with default Knitro.
# Build non-linear solver.
model = Model(with_optimizer(KNITRO.Optimizer, outlev=3))
load_mip_model!(model)
@time JuMP.optimize!(model)

#' We plot the solution by displaying in red the arcs removed from the graph.
# Plot!
if PLOT_GRAPH
    optimal_flow = JuMP.value.(model[:q])
    fig = figure()
    plot_network(flow=optimal_flow)
    display(fig)
end

#' In the resolving, Knitro uses the Branch & Bound algorithm to find
#' the optimal solution corresponding to the (convex) MINLP problem.
#' We refer to the [documentation](https://www.artelys.com/docs/knitro/2_userGuide/minlp.html)
#' for further details.
#' The approach used here is different than in BONMIN (outer approximation).

#' Observe that Knitro with default options computes 165 nodes before
#' finding the solution.
#' Is it possible to find a better tuning for Knitro?

#' ### MIP-Tuner

#' Since Knitro 12.0, a MINLP tuner was added to compute the optimal
#' parameterization to solve a given MINLP problem with Knitro. The MINLP
#' tuner uses an exhaustive search to find the optimal setting.
#'
#' By default, the tuner tests 36 combinations.
model = Model(with_optimizer(KNITRO.Optimizer, outlev=3, tuner=1))
load_mip_model!(model)
@time JuMP.optimize!(model)

#' Setting `mip_branchrule=2` and `mip_selectrule=3` seems to give
#' better results, according to the tuner.
#' Note that setting `outlev` to 0 switch off Knitro's output.
mip_solver = with_optimizer(KNITRO.Optimizer, mip_knapsack=0, mip_branchrule=2,
                            mip_selectrule=3, outlev=0)
model = Model(mip_solver)
load_mip_model!(model)
@time JuMP.optimize!(model)

#' ### Quantify impact of arcs' removals

#' The fewer arcs, the harder the problem.
#' If `nremovals >= 10`, Knitro takes too long to solve the problem.
#' However, the problem remains feasible till n_removals = 13.
#' Compute for instance the solution of the following problem:
feas_model = Model(with_optimizer(KNITRO.Optimizer))
@variable(feas_model, qc[1:nx])
@variable(feas_model, q[1:n])
@constraint(feas_model, q .== q0 + B*qc)
@variable(feas_model, z[1:n], Bin)
@constraint(feas_model,  q .<= QMAX * z)
@constraint(feas_model, -q .<= QMAX * z)
@objective(feas_model, Min, sum(z))

#' We know study the evolution of the objective cost w.r.t. the
#' number of removals.
max_removals = 9

# Save results in some arrays.
cost_values = Float64[]

for nremove in 1:max_removals
    println("Remove ", nremove, " arcs from graph.")
    model = Model(mip_solver)
    load_mip_model!(model, nremovals=nremove)
    @time JuMP.optimize!(model)

    push!(cost_values, JuMP.objective_value(model))

    # Plot!
    if PLOT_GRAPH
        optimal_flow = JuMP.value.(model[:q])
        fig = figure()
        plot_network(flow=optimal_flow)
        title("Optimal solution with $nremove removals")
        display(fig)
    end
end

#' Plot evolution of costs w.r.t. number of removals.
if PLOT_GRAPH
    fig = figure()
    plot(1:max_removals, cost_values, lw=3, c="k")
    xlabel("#removals")
    ylabel("Objective value")
    display(fig)
end

#' We observe that removing up to 6 arcs does not impact significantly
#' the cost.

max_stars_repo_name: frapac/Knitro-Zoo

max_issues_repo_name: frapac/Knitro-Zoo

max_forks_repo_name: frapac/Knitro-Zoo


text:
# includes -
include("Kinetics.jl")
include("Control.jl")
include("Inputs.jl")
include("DataDictionary.jl")
include("SolveBalances.jl")
include("Balances.jl")
using ODE

# Required by the tf_network_lib -
using JSON
using POETs
include("Error.jl")
include("Simulation.jl")
include("tf_network_lib.jl")

# Required by the error function -
using PyCall
@pyimport numpy as np

# Plotting -
using PyPlot
include("PlotLib.jl")

# Required by simple ensemble -
using Sobol


text:
function readdata(path::AbstractString)
	a = readlines(path)
	if a[1][1] == '>'
		fasta = []
		dna = []
		for i in a
			if i[1] == '>'
				push!(fasta,i)
				push!(dna,"")
			else
				dna[end] *= i
			end
		end
		return fasta,dna
	else
		entity = ""
		for i in a
			entity *= i
		end
		return entity

	end
end

path = "C:\\Users\\vsger\\Downloads\\rosalind_splc.txt"
fasta,dna = readdata(path)

dna_string = dna[1]
introns = dna[2:end]

DNACodonTable = Dict("TTT" => "F", "CTT" => "L", "ATT" => "I",
                       "GTT" => "V", "TTC" => "F", "CTC" => "L",
                       "ATC" => "I", "GTC" => "V", "TTA" => "L",
                       "CTA" => "L", "ATA" => "I", "GTA" => "V",
                       "TTG" => "L", "CTG" => "L", "ATG" => "M",
                       "GTG" => "V", "TCT" => "S", "CCT" => "P",
                       "ACT" => "T", "GCT" => "A", "TCC" => "S",
                       "CCC" => "P", "ACC" => "T", "GCC" => "A",
                       "TCA" => "S", "CCA" => "P", "ACA" => "T",
                       "GCA" => "A", "TCG" => "S", "CCG" => "P",
                       "ACG" => "T", "GCG" => "A", "TAT" => "Y",
                       "CAT" => "H", "AAT" => "N", "GAT" => "D",
                       "TAC" => "Y", "CAC" => "H", "AAC" => "N",
                       "GAC" => "D", "TAA" => "stop", "CAA" => "Q",
                       "AAA" => "K", "GAA" => "E", "TAG" => "stop",
                       "CAG" => "Q", "AAG" => "K", "GAG" => "E",
                       "TGT" => "C", "CGT" => "R", "AGT" => "S",
                       "GGT" => "G", "TGC" => "C", "CGC" => "R",
                       "AGC" => "S", "GGC" => "G", "TGA" => "stop",
                       "CGA" => "R", "AGA" => "R", "GGA" => "G",
                       "TGG" => "W", "CGG" => "R", "AGG" => "R",
                       "GGG" => "G")

function translate(dna)
	prot = ""
	for i in 1:3:(length(dna)-2)
		if DNACodonTable[dna[i:(i+2)]] != "stop"
			prot *= DNACodonTable[dna[i:(i+2)]]
		else
			break
		end
	end
	return prot
end

function splice(d,intron)
	l = length(intron)
	d = replace(d,intron => "")
	return d
end

function main()
	DNA = dna_string
	for i in introns
		DNA = splice(DNA,i)
	end
	protein = translate(DNA)
	return protein
end
main()

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
using JuMP
using ClimateMARGO
using ClimateMARGO.Models, ClimateMARGO.Optimization, ClimateMARGO.Diagnostics
using Test

function temperature_optimization_works(name::String, temp_goal::Float64)
    config_path = "../configurations"
    model = ClimateModel(ClimateMARGO.IO.included_configurations[name])
    status = optimize_controls!(model, temp_goal=temp_goal)
    return (
        (raw_status(status) == "Solve_Succeeded") &
        isapprox(
            maximum(T(model, M=true, R=true, G=true, A=true)),
            temp_goal,
            rtol=1.e-5
        )
    )
end

function tests()
    @testset "Subset of tests" begin
        for temp_goal in 1.5:0.5:4.0
            @test temperature_optimization_works("default", temp_goal)
        end
    end
end

tests()


text:
@testset "AbstractDemoCard" begin
    include("markdown.jl")
end


text:
@inline _select(::StaticInt{M}, ::StaticInt{M}) where {M} = StaticInt{M}()
@noinline _select(::StaticInt{M}, ::StaticInt{N}) where {M,N} = throw("$M ≠ $N")
@inline _select(::StaticInt{M}, _) where {M} = StaticInt{M}()
@inline _select(_, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline _select(x, _) = x

"""
Checks sizes for compatibility, and preserves the static size information if
given a mix of static and dynamic sizes.
"""
@inline function matmul_sizes(C,A,B)
    MC, NC = size(C)
    MA, KA = size(A)
    KB, NB = size(B)
    @assert ((MC == MA) & (KA == KB) & (NC == NB)) "Size mismatch."
    (_select(MA, MC), _select(KA, KB), _select(NB, NC))
end

function unsafe_copyto_avx!(pB, pA, M, N)
    LoopVectorization.@avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
        pB[m,n] = pA[m,n]
    end
end

function default_stridedpointer_quote(::Type{T}, N, Ot) where {T}
    C = 1
    B = 0
    R = Expr(:tuple)
    o = Expr(:tuple)
    xt = Expr(:tuple)
    st = Expr(:call, Expr(:curly, :StaticInt, sizeof(T)))
    for n ∈ 1:N
        push!(xt.args, Expr(:call, :*, :st, Expr(:ref, :x, n)))
        push!(R.args, n)
        push!(o.args, Expr(:call, Ot))
    end
    quote
        $(Expr(:meta,:inline))
        st = $st
        StridedPointer{$T,$N,$C,$B,$R}(ptr, $xt, $o)
    end
end

@generated function default_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :One)
end
@generated function default_zerobased_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :Zero)
end


text:
module SorbetHelperKillZone

using ..Ahorn, Maple

@mapdef Entity "SorbetHelper/KillZone" KillZone(
    x::Integer,
    y::Integer,
    width::Integer=8,
    height::Integer=8,
    flag::String="",
    inverted::Bool=false,
    fastKill::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Kill Zone (Sorbet Helper)" => Ahorn.EntityPlacement(
        KillZone,
        "rectangle",
    )
)

Ahorn.minimumSize(entity::KillZone) = 8, 8
Ahorn.resizable(entity::KillZone) = true, true
Ahorn.editingOrder(entity::KillZone) = String["x", "y", "width", "height", "flag", "inverted", "fastKill"]

function Ahorn.selection(entity::KillZone)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::KillZone, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (176 / 255, 99 / 255, 100 / 255, 0.45), (145 / 255, 59 / 255, 95 / 255, 0.82))
end

end

max_stars_repo_name: Kitty-Cats/SorbetHelper

max_issues_repo_name: Kitty-Cats/SorbetHelper

max_forks_repo_name: Kitty-Cats/SorbetHelper


text:
module RT

include("vec3.jl")
import .V3: Vec3, at_len, x, y, z

using LinearAlgebra
using Images

struct Ray
    origin::Vec3
    arrow::Vec3
    Ray(o::Vec3, a::Vec3) = new(o, normalize(a))
end

struct Camera
    scr_w::UInt16
    scr_h::UInt16
    view_w::Float64
    view_h::Float64
    eye_pos::Vec3
    eye_up::Vec3
    eye_right::Vec3
    view_tl::Vec3
end

function Camera(scr_w, scr_h)
    ASPECT_R = scr_w / scr_h
    VIEW_HH = 1.0
    VIEW_HW = ASPECT_R * VIEW_HH
    view_h = VIEW_HH * 2
    view_w = VIEW_HW * 2

    eye_pos = Vec3([0, 0, 0])
    EYE_DIR = Vec3([0, 0, -1])
    eye_up = Vec3([0, 1, 0])
    eye_right = EYE_DIR × eye_up

    FOCAL_LEN = 1.0
    VIEW_CENTER = eye_pos + at_len(EYE_DIR, FOCAL_LEN)
    view_tl = VIEW_CENTER + at_len(eye_up, VIEW_HH) - at_len(eye_right, VIEW_HW)

    Camera(scr_w, scr_h, view_w, view_h, eye_pos, eye_up, eye_right, view_tl)
end

get_ray(self::Camera, u, v) =
    Ray(
        self.eye_pos,
        self.view_tl
        + at_len(self.eye_right, u * self.view_w)
        - at_len(self.eye_up, v * self.view_h)
        - self.eye_pos
    )

function render(self::Camera, ray_color_fn, print_progress=true)
    img = zeros(RGB, self.scr_h, self.scr_w)
    for j in 1:self.scr_h
        if print_progress && j % 100 == 0
            println("$j of $(self.scr_h)")
        end
        v = (j - .5) / self.scr_h
        for i in 1:self.scr_w
            u = (i - .5) / self.scr_w
            ray = get_ray(self, u, v)
            img[j, i] = ray_color_fn(ray)
        end
    end
    img
end

abstract type Hittable end

struct HitRecordTest
    t::Float64
    ray::Ray
    obj::Hittable
end

struct HitRecordFull
    front::Bool
    normal::Vec3
    point::Vec3
end

struct Sphere <: Hittable
    center::Vec3
    r::Float64
end

function hit_test(self::Sphere, ray::Ray, tmin=0, tmax=Inf)::Union{HitRecordTest, Nothing}
    oc = ray.origin - self.center
    half_b = oc ⋅ ray.arrow
    c = oc ⋅ oc - self.r^2
    Δ = half_b^2 - c
    if Δ < 0
        return nothing
    else
        sΔ = sqrt(Δ)
        root = -half_b - sΔ
        if root < tmin || root > tmax
            root = -half_b + sΔ
            if root < tmin || root > tmax
                return nothing
            end
        end
        return HitRecordTest(root, ray, self)
    end
end

function hit_full(self::Sphere, hit::HitRecordTest)::HitRecordFull
    p = at_len(hit.ray.arrow, hit.t)
    normal = normalize(p - self.center)
    front = true
    if hit.ray.arrow ⋅ normal > 0
        front = false
        normal = -normal
    end
    HitRecordFull(front, normal, p)
end

function bg_color(ray::Ray)::RGB
    t = (y(ray.arrow) + 1) / 2
    blue = RGB(.5, .7, 1)
    white = RGB(1, 1, 1)
    t * blue + (1 - t) * white
end

spheres = [
    Sphere(Vec3([0, 0, -1]), .5),
    Sphere(Vec3([0, -500.5, -1]), 500),
]

function ray_color(ray::Ray)::RGB
    hit = nothing
    for sphere in spheres
        new_hit = hit_test(sphere, ray)
        if !isnothing(new_hit) && (isnothing(hit) || new_hit.t < hit.t)
            hit = new_hit
        end
    end
    if !isnothing(hit)
        full = hit_full(hit.obj, hit)
        factor = full.front ? 1 : .5
        col = (full.normal + Vec3([1, 1, 1])) * factor / 2
        RGB(x(col), y(col), z(col))
    else
        bg_color(ray)
    end
end

end # module RT

max_stars_repo_name: hdtx/raytracer-comparison

max_issues_repo_name: hdtx/raytracer-comparison

max_forks_repo_name: hdtx/raytracer-comparison