text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using FMIExport
using FMICore

originalGetReal = nothing # function pointer to the original fmi2GetReal c-function

# custom function for fmi2GetReal!(fmi2Component, Union{Array{fmi2ValueReference}, Ptr{fmi2ValueReference}}, Csize_t, value::Union{Array{fmi2Real}, Ptr{fmi2Real}}::fmi2Status
# for information on how the FMI2-functions are structured, have a look inside FMICore.jl/src/FMI2_c.jl or the FMI2.0.3-specification on fmi-standard.org
function myGetReal!(c::fmi2Component, vr::Union{Array{fmi2ValueReference}, Ptr{fmi2ValueReference}}, nvr::Csize_t, value::Union{Array{fmi2Real}, Ptr{fmi2Real}})
    global originalGetReal

    # first, we do what the original function does
    status = fmi2GetReal!(originalGetReal, c, vr, nvr, value)

    # if we have a pointer to an array, we must interprete it as array to access elements
    if isa(value, Ptr{fmi2Real})
        value = unsafe_wrap(Array{fmi2Real}, value, nvr, own=false)
    end
    if isa(vr, Ptr{fmi2Real})
        vr = unsafe_wrap(Array{fmi2Real}, vr, nvr, own=false)
    end

    # now, we multiply the position sensor output by two (just for fun!)
    for i in 1:nvr
        if vr[i] == 335544320 # value reference for "positionSensor.s"
            value[i] *= 2.0
        end
    end

    # ... and we return the original status
    return status
end

# this function is called, as soon as the DLL is loaded and Julia is initialized
# must return a FMU2-instance to work with
FMIBUILD_CONSTRUCTOR = function(resPath)
    global originalGetReal

    # loads an existing FMU inside the FMU
    fmu = fmiLoad(joinpath(resPath, "SpringDamperPendulum1D.fmu"))

    # save, where the original `fmi2GetReal` function was stored, so we can access it in our new function
    originalGetReal = fmu.cGetReal

    # now we overwrite the original function
    fmiSetFctGetReal(fmu, myGetReal!)

    return fmu
end

### FMIBUILD_NO_EXPORT_BEGIN ###
# The line above is a start-marker for excluded code for the FMU compilation process!

import FMIZoo

tmpDir = mktempdir(; prefix="fmibuildjl_test_", cleanup=false)
@info "Saving example files at: $(tmpDir)"
fmu_save_path = joinpath(tmpDir, "Manipulation.fmu")

sourceFMU = FMIZoo.get_model_filename("SpringDamperPendulum1D", "Dymola", "2022x")
fmu = FMIBUILD_CONSTRUCTOR(dirname(sourceFMU))
import FMIBuild:fmi2Save        # <= this must be excluded during export, because FMIBuild cannot execute itself (but it is able to build)
fmi2Save(fmu, fmu_save_path; resources=Dict(sourceFMU=>"SpringDamperPendulum1D.fmu"))    # <= this must be excluded during export, because fmi2Save would start an infinte build loop with itself

# The following line is a end-marker for excluded code for the FMU compilation process!
### FMIBUILD_NO_EXPORT_END ###


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Pango"
version = v"1.42.4"

# Collection of sources required to build Pango
sources = [
    "http://ftp.gnome.org/pub/GNOME/sources/pango/$(version.major).$(version.minor)/pango-$(version).tar.xz" =>
    "1d2b74cd63e8bd41961f2f8d952355aa0f9be6002b52c8aa7699d9f5da597c9d"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pango-*/
./configure --prefix=$prefix --host=$target --enable-introspection=no
# The generated Makefile tries to build some examples in the "tests" directory,
# but this would fail.  Let's skip it.
sed -i 's/^\(SUBDIRS = .*\) tests/\1/' Makefile
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpango", :libpango)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Pango-only dependencies
    "https://github.com/giordano/FriBidiBuilder/releases/download/v1.0.5/build_FriBidi.v1.0.5.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libffi-v3.2.1-0/build_Libffi.v3.2.1.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Gettext-v0.19.8-0/build_Gettext.v0.19.8.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libiconv-v1.15-0/build_Libiconv.v1.15.0.jl",
    # Freetype2-related dependencies
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-4/build_FreeType2.v2.10.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Bzip2-v1.0.6-2/build_Bzip2.v1.0.6.jl",
    # Glib-related dependencies
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Glib-v2.59.0%2B0/build_Glib.v2.59.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/PCRE-v8.42-2/build_PCRE.v8.42.0.jl",
    # Cairo-related dependencies
    "https://github.com/giordano/Yggdrasil/releases/download/Cairo-v1.14.12/build_Cairo.v1.14.12.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Pixman-v0.36.0-0/build_Pixman.v0.36.0.jl",
    "https://github.com/JuliaIO/LibpngBuilder/releases/download/v1.0.3/build_libpng.v1.6.37.jl",
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: giordano/PangoBuilder

max_issues_repo_name: giordano/PangoBuilder

max_forks_repo_name: giordano/PangoBuilder


text:
"""
    xor_swap(a, b)

This is an old trick to exchange the values of variables a and b without using
an additional temporary variable.
"""
macro xor_swap(a, b)
    esc(:($a ⊻= $b; $b ⊻= $a; $a ⊻= $b))
end

"""
    swapbits{T<:Unsigned}(x::T, i::Integer, j::Integer)

Swap the i'th and j'th bits in an unsigned integer.
Note this uses zero based indexes for `i` and `j`.

E.g. to swap the LSB and MSB of a byte: 1001 1000 (0x98) -> 0001 1001 (0x19)

```@example
swapbits(0x98, 0, 7)
```
"""
@inline function swapbits(x::T, i::Integer, j::Integer) where {T <: Unsigned}
    ibit = (x >> i) & T(1)
    jbit = (x >> j) & T(1)
    ixj = ibit ⊻ jbit
    ixj = (ixj << i) | (ixj << j)
    return x ⊻ ixj
end

@inline function _swap_bits_chunks_kernel(x, mask, offset)
    msk = repeatpattern(typeof(x), mask)
    return ((x >> offset) & msk) | ((x & msk) << offset)
end

@inline swap_odd_even_bits(x::Unsigned)        = _swap_bits_chunks_kernel(x, 0x55, 1)
@inline swap_consecutive_bitpairs(x::Unsigned) = _swap_bits_chunks_kernel(x, 0x33, 2)
@inline swap_nibbles(x::Unsigned)              = _swap_bits_chunks_kernel(x, 0x0F, 4)
@inline swap_bytes(x::Unsigned)                = _swap_bits_chunks_kernel(x, 0x00FF, 8)
@inline swap_uint16s(x::Unsigned)              = _swap_bits_chunks_kernel(x, 0x0000FFFF, 16)
@inline swap_uint32s(x::Unsigned)              = _swap_bits_chunks_kernel(x, 0x00000000FFFFFFFF, 32)
@inline swap_uint64s(x::Unsigned)              = _swap_bits_chunks_kernel(x, 0x0000000000000000FFFFFFFFFFFFFFFF, 64)


text:
using SurfaceWaves
using Test

@testset "SurfaceWaves.jl" begin
    # Write your tests here.
end


text:
using Random
abstract type AbstractDataStream end

"""
Information representing a data stream.
"""
mutable struct Stream
    n_samples::Int
    n_targets::Int
    n_features::Int
    n_num_features::Int
    n_cat_features::Int
    n_classes::Int
    n_cat_features_idx::Vector{Int}
    current_sample_x::Vector{Float64}
    current_sample_y::Vector{Int}
    sample_idx::Int
    feature_names::Vector{String}
    target_names::Vector{String}
    target_values::Vector{Int}
    name::String
    random_state_seed::Int
    random_state::AbstractRNG
end

function Stream(;n_samples::Int,
    n_targets::Int,
    n_features::Int,
    n_num_features::Int,
    n_cat_features::Int,
    n_classes::Int,
    n_cat_features_idx::Vector{Int},
    current_sample_x::Vector{Float64},
    current_sample_y::Vector{Int},
    sample_idx::Int,
    feature_names::Vector{String},
    target_names::Vector{String},
    target_values::Vector{Int},
    name::String,
    random_state_seed::Int,
    random_state::AbstractRNG)

    Stream(n_samples, n_targets, n_features, n_num_features, n_cat_features,
        n_classes, n_cat_features_idx, current_sample_x, current_sample_y, sample_idx,
        feature_names, target_names, target_values, name, random_state_seed, random_state)
end

function Stream(name::String, n_features::Int64)
    seed = 42
    rng = MersenneTwister(seed)
    return Stream(-1, 1, n_features, n_features, 0, 2, [], [], [], 0, ["f$(i)" for i in 1:n_features], ["label"], [0, 1], name, seed, rng)
end

function classify_instance(s::AbstractDataStream, X)
    return 0
end

function get_stream_name(s::AbstractDataStream)
    return s.stream.name
end

"""
Generate the next X,y sample in the stream.
X is a 2 dimensional vector, uniformly sampled between 0-1.
y is determined by the classification_function of the SineGenerator.
"""
function next_sample!(s::AbstractDataStream, batch_size)
    data_X = Array{Float64}(undef, s.stream.n_features, batch_size)
    data_Y = Vector{Int64}(undef, batch_size)
    for j in 1:batch_size
        s.stream.sample_idx += 1
        for f in 1:s.stream.n_num_features
            feature = next_num_feature(s, f)
            data_X[f, j] = feature
        end
        for f in s.stream.n_num_features+1:s.stream.n_features
            feature = next_cat_feature(s, f)
            data_X[f, j] = feature
        end
        data_Y[j] = classify_instance(s, data_X[:, j])
        s.stream.current_sample_x = view(data_X, :, j)
        s.stream.current_sample_y = view(data_Y, j:j)
    end
    return data_X, data_Y
end

"""
Reset the RNG state sample counter to starting state.
"""
function stream_reset!(s::AbstractDataStream)
    rng = MersenneTwister(s.stream.random_state_seed)
    s.stream.random_state = rng
    s.stream.sample_idx = 0
    return
end

"""
Each feature is uniform random in range [0, 1]
"""
function next_num_feature(s::AbstractDataStream, f_id::Int)
    rand(s.stream.random_state)
end

"""
Each feature is an int between 1 and 10
"""
function next_cat_feature(s::AbstractDataStream, f_id::Int)
    rand(s.stream.random_state, 1:10)
end


text:
@eval fun1(x) = $a
Base.@eval fun2(x) = $a

@eval function fun3(x)
    $a
end
Base.@eval function fun4(x)
    $a
end

@eval struct Test1
    a::$T
end
Base.@eval struct Test2
    a::$T
end

@eval macro macro1(x)
    $a
end
Base.@eval macro macro2(x)
    $a
end


text:
"""
    DiracDelta{T,TM<:Manifolds.Manifold} <: Manifolds.MPointDistribution{TM}

Dirac distribution around a point ``p``, that is, a distribution whose support is a single
point ``p`` on a manifold ``M``.

# Constructor

    DiracDelta(M::Manifold, p)
"""
struct DiracDelta{P,M<:Manifolds.Manifold} <: Manifolds.MPointDistribution{M}
    manifold::M
    point::P
end

Base.eltype(::Type{<:DiracDelta{T}}) where {T} = T

function Base.rand(::AbstractRNG, s::Random.SamplerTrivial{<:DiracDelta})
    d = s.self
    return copy(d.point)
end

Distributions._rand!(::AbstractRNG, d::DiracDelta, q) = copyto!(q, d.point)

function Distributions.convolve(
    d1::D1, d2::D1
) where {P1,P2,G<:Manifolds.AbstractGroupManifold,D1<:DiracDelta{P1,G},D2<:DiracDelta{P2,G}}
    M = d1.manifold
    p = Manifolds.compose(M, d1.point, d2.point)
    return DiracDelta(M, p)
end

function Distributions.insupport(d::DiracDelta, p)
    M = d.manifold
    return Manifolds.is_manifold_point(M, p) && isapprox(M, p, d.point)
end

Statistics.mean(d::DiracDelta) = d.point

Statistics.median(d::DiracDelta) = Statistics.mean(d)

Statistics.std(d::DiracDelta) = zero(real(Manifolds.number_eltype(d.point)))

Statistics.var(d::DiracDelta) = zero(real(Manifolds.number_eltype(d.point)))

Distributions.mode(d::DiracDelta) = Statistics.mean(d)

inversion(d::DiracDelta) = d


text:
# const ID = Int
const VectorID = Vector{Int}

# @inline combine(a::Vector{T}, b::Vector{T}) where T =
@inline append(a::Vector{T}, b::Vector{T}) where T = vcat(a, b)
# @inline append(a::Vector{T}, b::T) where T = @show vcat(a, T[b])
@inline singletonid(::Type{Vector{T}}, i::Vararg{T, N}) where {T, N} = T[i...]

# @inline firstelem(::Vector{T}) where T = zero(T)
# function increment(a::Vector)
#   b = copy(a)
#   b[end] += 1
#   b
# end`


text:
function check_boundary!(r::Robot, side)
    while !isborder(r, side)
        move!(r, side)
    end
end
function mark_side!(r::Robot, side)
    while !isborder(r, side)
        putmarker!(r)
        move!(r, side)
    end
end
function mark_boundaries!(r::Robot)
    check_boundary!(r, Sud)
    check_boundary!(r, West)
    mark_side!(r, Nord)
    mark_side!(r, Ost)
    mark_side!(r, Sud)
    mark_side!(r, West)
end


text:
module SimplicialComplexes

using DataStructures

export SimplicialComplex, alllabels

include("utils.jl")
include("simplicial_complex.jl")

end


text:
function lorenz(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 [dx,dy,dz]
end

using DifferentialEquations, BenchmarkTools
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,Tsit5(),save_everystep=false)

function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,Tsit5(),save_everystep=false)

tspan = (0.0,500.0) # 5x longer than before
prob = ODEProblem(lorenz!,u0,tspan)
@benchmark solve(prob,Tsit5(),save_everystep=false)

using StaticArrays
A = @SVector [2.0,3.0,5.0]

function lorenz_static(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 @SVector [dx,dy,dz]
end

u0 = @SVector [1.0,0.0,0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan)
@benchmark solve(prob,Tsit5())

@benchmark solve(prob,Tsit5(),save_everystep=false)

A = rand(1000,1000); B = rand(1000,1000); C = rand(1000,1000)
test(A,B,C) = A + B + C
@benchmark test(A,B,C)

test2(A,B,C) = map((a,b,c)->a+b+c,A,B,C)
@benchmark test2(A,B,C)

function test3(A,B,C)
    D = similar(A)
    @inbounds for i in eachindex(A)
        D[i] = A[i] + B[i] + C[i]
    end
    D
end
@benchmark test3(A,B,C)

test4(A,B,C) = A .+ B .+ C
@benchmark test4(A,B,C)

sin.(A) .+ sin.(B)

test5(A,B,C) = @. A + B + C #only one array allocated
@benchmark test5(A,B,C)

D = zeros(1000,1000);

test6!(D,A,B,C) = D .= A .+ B .+ C #only one array allocated
@benchmark test6!(D,A,B,C)

test7!(D,A,B,C) = @. D = A + B + C #only one array allocated
@benchmark test7!(D,A,B,C)

test8!(D,A,B,C) = map!((a,b,c)->a+b+c,D,A,B,C)
@benchmark test8!(D,A,B,C)

@benchmark A*B

using LinearAlgebra
@benchmark mul!(D,A,B) # same as D = A * B

# Generate the constants
p = (1.0,1.0,1.0,10.0,0.001,100.0) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
Ax[2,1] = 2.0
Ax[end-1,end] = 2.0
Ay[1,2] = 2.0
Ay[end,end-1] = 2.0

function basic_version!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = r[:,:,1]
  v = r[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  dr[:,:,1] = Du .+ a.*u.*u./v .+ ubar .- α*u
  dr[:,:,2] = Dv .+ a.*u.*u .- β*v
end

a,α,ubar,β,D1,D2 = p
uss = (ubar+β)/α
vss = (a/β)*uss^2
r0 = zeros(100,100,2)
r0[:,:,1] .= uss.+0.1.*rand.()
r0[:,:,2] .= vss

prob = ODEProblem(basic_version!,r0,(0.0,0.1),p)

@benchmark solve(prob,Tsit5())

A = rand(4)
@show A
B = @view A[1:3]
B[2] = 2
@show A

function gm2!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  Du = D1*(Ay*u + u*Ax)
  Dv = D2*(Ay*v + v*Ax)
  @. du = Du + a.*u.*u./v + ubar - α*u
  @. dv = Dv + a.*u.*u - β*v
end
prob = ODEProblem(gm2!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

Ayu = zeros(N,N)
uAx = zeros(N,N)
Du = zeros(N,N)
Ayv = zeros(N,N)
vAx = zeros(N,N)
Dv = zeros(N,N)
function gm3!(dr,r,p,t)
  a,α,ubar,β,D1,D2 = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm3!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2
function gm4!(dr,r,p,t)
  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm4!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

p = (1.0,1.0,1.0,10.0,0.001,100.0,N)
function fast_gm!(du,u,p,t)
  a,α,ubar,β,D1,D2,N = p

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end

  @inbounds for j in 2:N-1, i in 2:N-1
    du[i,j,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for j in 2:N-1
    i = 1
    du[1,j,1] = D1*(2u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = 1
    du[1,j,2] = D2*(2u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,1] = D1*(2u[i-1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
           a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for j in 2:N-1
    i = N
    du[end,j,2] = D2*(2u[i-1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
           a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds for i in 2:N-1
    j = 1
    du[i,1,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = 1
    du[i,1,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
  end
  @inbounds for i in 2:N-1
    j = N
    du[i,end,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
  end

  @inbounds begin
    i = 1; j = 1
    du[1,1,1] = D1*(2u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,1,2] = D2*(2u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
              a*u[i,j,1]^2 - β*u[i,j,2]

    i = 1; j = N
    du[1,N,1] = D1*(2u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[1,N,2] = D2*(2u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = 1
    du[N,1,1] = D1*(2u[i-1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[N,1,2] = D2*(2u[i-1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]

    i = N; j = N
    du[end,end,1] = D1*(2u[i-1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
             a*u[i,j,1]^2/u[i,j,2] + ubar - α*u[i,j,1]
    du[end,end,2] = D2*(2u[i-1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
             a*u[i,j,1]^2 - β*u[i,j,2]
   end
end
prob = ODEProblem(fast_gm!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

prob = ODEProblem(fast_gm!,r0,(0.0,10.0),p)
@benchmark solve(prob,Tsit5())

using Sundials
@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))

prob = ODEProblem(fast_gm!,r0,(0.0,100.0),p)
# Will go out of memory if we don't turn off `save_everystep`!
@benchmark solve(prob,Tsit5(),save_everystep=false)

@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES))

@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)

prob = ODEProblem(fast_gm!,r0,(0.0,500.0),p)
@benchmark solve(prob,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module _Math

using Brainstorm._Math
using Base.Test

function test_math_setbits()
    @test setbits(0, []) == 0
    @test setbits(0, [0]) == 2^0
    @test setbits(0, [1]) == 2^1
    @test setbits(0, [10]) == 2^10
    @test setbits(0, [1,3]) == 2^1 + 2^3

    @test setbits(12, []) == 12
    @test setbits(12, [2,3]) == 12
    @test setbits(12, [0,1,2,3]) == 15

    @test setbits(big(0), []) == big(0)
    @test setbits(big(0), [1000]) == big(2)^1000
end

include("Series/series.jl")
include("NumberTheory/numbertheory.jl")
include("Combinatorics/combinatorics.jl")
include("Algebra/algebra.jl")
include("EllipticCurves/ellipticcurves.jl")

function test_all()
    println("")
    print(rpad("Math...", 50, ' '))

    test_math_setbits()

    println("PASS")

    println("")
    _Series.test_all()
    println("")
    _NumberTheory.test_all()
    println("")
    _Combinatorics.test_all()
    println("")
    _Algebra.test_all()
    println("")
    _EllipticCurves.test_all()
end

end


text:
using ADCME
using AdFem

include("../chip_unstructured_solver.jl")
include("../chip_unstructured_geometry.jl")

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

function k_exact(x, y)
    k_mold + 1000 * k_chip_ref * (x-0.49)^2 / (1 + x^2)
end

function k_nn(x, y, θ)
    fc(x, [20,20,20,1], θ)
end

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 14    # number of iterations for Newton's method

xy = mesh.nodes
xy2 = zeros(mesh.nedge, 2)
for i = 1:mesh.nedge
    xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
end
xy = [xy;xy2]

x, y = xy[chip_fem_idx, 1], xy[chip_fem_idx, 2]
k_chip_exact = constant( @. k_exact(x, y) )

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

B = constant(compute_interaction_matrix(mesh))
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
      bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

S0 = constant(zeros(nelem+3*ndof))
S_data = matread("diagnose1/data.mat")["V"]

ntest = 100
max_iter = 20
dist_test = zeros(ntest)
loss_test = zeros(ntest)

θ_trained = matread("k_diagnose_theta.mat")["theta"]

for t in 1:ntest
    @info t
# begin of FOR LOOP
# generate one set of NN weights and biases

    n = 901
    θ = Variable(randn(n))
    d = rand() # for diagnose1+ and diagnose3, sample distance between 1 and 5
    θ = θ / norm(θ) * d
    k_chip = k_nn(x, y, θ .+ θ_trained)

    S = solve_navier_stokes(S0, NT, k_chip)
    S_computed = S[end, :]

    loss =  mean((S_computed .- S_data)^2)
    loss = loss * 1e10

    # ---------------------------------------------------
    _loss = Float64[]
    cb = (vs, iter, loss)->begin
        _loss
        push!(_loss, loss)
        printstyled("[#t $t] [#iter $iter] loss=$loss\n", color=:green)
        if mod(iter, 5)==1
            matwrite(string("diagnose2/test$t","_loss$iter.mat"), Dict("iter"=>iter,"loss"=>_loss, "k_chip"=>vs[1], "d"=>d[1]))
        end
    end

    sess = Session(); init(sess)
    loss = BFGS!(sess, loss, max_iter; vars = [k_chip], callback = cb)

    # record distance d; loss after a fixed number of iterations (converge or not?)
    dist_test[t] = d
    loss_test[t] = loss[end, 1]

    matwrite("diagnose2/summary$t.mat", Dict("dist"=>dist_test,"loss"=>loss_test))
# end of FOR LOOP
end

# ---------------------------------------------------


text:
@testset ExtendedTestSet "dice" begin
    @testset ExtendedTestSet "dice" begin
        ŷ = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        @test dice(ŷ, y) == 0.0
    end

    @testset ExtendedTestSet "dice" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y2 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y2; dims=3)

        ŷ1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        ŷ2 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        ŷ = cat(ŷ1, ŷ2; dims=3)
        @test dice(ŷ, y) == 0.0
    end
end


text:
export smooth, iscontinuous, isdifferentiable, isanalytic
module smoothFunction

#--- Real Functions
function iscontinuous(x)
    """
    checks whether function x is continuous
    """
return true #skip for now
end

function isdifferentiable(x)
    """
    Checks whether function x isDifferentiable
        TODO: pick a Suitable differentiation Module
        TODO: Call a differentiation Module
    """
    return true # skip for now
end

function smooth(x=exp(x))
    """
    A smooth function basis function, that is:
        1.continuous #TODO
        2.differentiablle #TODO
    """
    smooth = exp(x)
    if isdifferentiable(smooth) &&  iscontinuous(smooth) #TODO: verification functions
        return smooth
    end
end

function smooth(x = exp(x),numerator=a,denominator=b)
    """
    A smooth function basis function, that is:
        1.continuous
        2.differentiablle
    """
    smooth = exp(x)
    if isdifferentiable(smooth) && iscontinuous(smooth)
        return denominator!=0 ? (numerator/denominator)*smooth : NaN
    end
end

#--- complex functions #testing

function analytic()
    """

    TODO: needs further review
    analytic function to be used as input for a laplacian
    """
    isanalytic = nothing
    if  isanalytic
        isanalytic = true

    else
        isanalytic = false

    end
 return isanalytic
end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
mutable struct Data3D
	mat::Matrix{Float32}
	bound::Shape{Float32}

	Data3D(bound = empty!(AABB()), mat = Matrix{Float32}(I, 4, 4)) = new(mat, bound)
end

function assure_bound_type(data::Data3D, bound::Shape)
	if typeof(data.bound) != typeof(bound)
		data.bound = similar(bound)
	end
	nothing
end

mutable struct Spatial <: LeafObj
	parent::NodeObj
	local3D::Data3D
	world3D::Data3D
	next3D::Data3D

	Spatial(matLocal, boundLocal) = new(NoNode(), Data3D(boundLocal, matLocal), Data3D(similar(boundLocal)), Data3D(similar(boundLocal)))
end

function init(engine::Engine, ::Type{Spatial}, def::Dict{Symbol, Any})::Spatial
	matLocal = get_transform(def)
	boundLocal = Empty{Float32}()
	if haskey(def, :boundLocal)
		boundLocal = load_def(engine, def[:boundLocal])
	end
	statial = Spatial(matLocal, boundLocal)
end

get_id(t::Spatial) = :spatial

get_local_transform(spatial::Spatial) = spatial.local3D.mat
function set_local_transform(spatial::Spatial, m::Matrix{Float32})
	spatial.next3D[:] = m
	world = top(spatial)
	!isa(world, BaseWorld) && return
	register_transform_update(world, spatial.parent)
	nothing
end

get_local_bound(spatial::Spatial) = spatial.local3D.bound
function set_local_bound(spatial::Spatial, bound::Shape{Float32})
	spatial.next3D.bound = bound
	world = top(spatial)
	!isa(world, BaseWorld) && return
	register_bound_update(world, spatial.parent)
	nothing
end

get_world_transform(spatial::Spatial) = spatial.world3D.mat
get_world_bound(spatial::Spatial) = spatial.world3D.bound

function next_transform(spatial::Spatial)
	spatial.local3D.mat[:] = spatial.next3D.mat
end

function next_bound(spatial::Spatial)
	assure_bound_type(spatial.local3D, spatial.next3D.bound)
	assign(spatial.local3D.bound, spatial.next3D.bound)
end

function update_world_transform(spatial::Spatial, worldMat::Matrix{Float32})
	spatial.world3D.mat = worldMat * spatial.local3D.mat
end

function update_world_bound(spatial::Spatial)
	assure_bound_type(spatial.world3D, spatial.local3D.bound)
	transform(spatial.world3D.bound, spatial.world3D.mat, spatial.local3D.bound)
end

mutable struct Visual <: LeafObj
	id::Symbol
	parent::NodeObj
	visual::GRU.Renderable
	matLocal::Matrix{Float32}

	Visual(id::Symbol, visual::GRU.Renderable, matLocal::Matrix{Float32} = Matrix{Float32}(I, 4, 4)) = new(id, NoNode(), visual, matLocal)
end

function init(engine::Engine, ::Type{Visual}, def::Dict{Symbol, Any})::Visual
	id = get_id!(def)
	visual = load_def(engine, def[:visual])
	matLocal = get_transform(def)
	vis = Visual(id, visual, matLocal)
	vis
end

get_id(v::Visual) = v.id
get_local_bound(v::Visual) = transform(v.matLocal, GRU.localbound(v.visual))

function update(v::Visual)
	parentWorld = get_world_transform(v.parent)
	GRU.settransform(v.visual, parentWorld * v.matLocal)
end

function render(v::Visual, engine::Engine)
	update(v)
	GRU.add(engine.renderer, v.visual)
end


text:
module RocketSafeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: safe()" begin

    println("Testing: operator safe()")

    run_proxyshowcheck("Safe", safe())

    run_testset([
        (
            source = from(1:5) |> safe(),
            values = @ts([ 1:5, c ])
        ),
        (
            source = from([ 0, 1, 2 ]) |> safe() |> map(Int, d -> d === 0 ? 0 : throw(d)),
            values = @ts([ 0, e(1) ])
        ),
        (
            source = completed() |> safe(),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> safe(),
            values = @ts(e("e"))
        ),
        (
            source = never() |> safe(),
            values = @ts()
        )
    ])

end

end


text:
#=
sino_plot.jl
sinogram geometry plots
2020-07-24, Jeff Fessler, University of Michigan
=#

export sino_geom_plot!, sino_geom_plot_grids

# using MIRT: jim, image_geom, ImageGeom
# using MIRT: SinoGeom, sino_geom,
# using MIRT: SinoPar, SinoFan, SinoMoj, SinoFanArc, SinoFanFlat

"""
    sino_geom_plot_grid(sg, Plots.scatter) <- trick to avoid Plots.dependency
Scatter plot of `(r,phi)` sampling locations from `sg.grid`
"""
function sino_geom_plot_grid(sg::SinoGeom, scatter::Function)
	(r, phi) = sg.grid
	ylim = [min(0, rad2deg(minimum(phi))), max(360, rad2deg(maximum(phi)))]
	rmax = ceil(maximum(abs.(r))/10, digits=0)*10
	scatter(r, rad2deg.(phi), label="", markersize=1, markerstrokecolor=:auto,
		markershape=:circle, linewidth=0, # emulate scatter
		ylim = ylim, xlim = [-1,1]*rmax, # ylabel="ϕ",
		title="$(typeof(sg))", xtick=(-1:1)*rmax, ytick=[0,360])
end

"""
    sino_geom_plot_grids()
scatter plot of (r,phi) sampling locations for all geometries
"""
function sino_geom_plot_grids(
	plot::Function ;
	orbit::Real = 360, down::Int = 30,
)
	geoms = (
		sino_geom(:par ; nb = 888, na = 984, down=down, d = 0.5, orbit=orbit,
			offset = 0.25),
		sino_geom(:fan ; nb = 888, na = 984, d = 1.0, orbit = orbit,
			offset = 0.75, dsd = 949, dod = 408, down=down),
		sino_geom(:fan ; nb = 888, na = 984, d = 1.0, orbit = orbit,
			offset = 0.75, dsd = 949, dod = 408, down=down,
			dfs = Inf, source_offset = 0.7), # flat fan
		sino_geom(:moj ; nb = 888, na = 984, down=down, d = 1.0, orbit=orbit,
			offset = 0.25),
	)

	ngeom = length(geoms)
	pl = Array{Any}(undef, ngeom)

	for ii=1:ngeom
		sg = geoms[ii]
		pl[ii] = sg.plot_grid(plot)
	end
	return pl
end

"""
    sino_geom_plot!(sg, plot! ; ig)
Picture of the source position / detector geometry
"""
function sino_geom_plot!(
	sg::SinoGeom,
	plot!::Function ; # trick to avoid Plot dependency in MIRT
	ig::Union{Nothing,ImageGeom} = nothing,
)
	plot!(aspect_ratio=1)

	scat!(args... ; kwargs...) =
		plot!(args... ; kwargs..., markerstrokecolor=:auto, linewidth=0)

	xmax = sg.rfov; xmin = -xmax; (ymin,ymax) = (xmin,xmax)
	if !isnothing(ig)
#		plot!(jim(ig.x, ig.y, ig.mask[:,:,1], clim=(0,1))) # todo: jim!
		xmin = minimum(ig.x); xmax = maximum(ig.x)
		ymin = minimum(ig.y); ymax = maximum(ig.y)
	end
	plot!([xmax, xmin, xmin, xmax, xmax],
		[ymax, ymax, ymin, ymin, ymax], color=:green, label="")
	plot!(xtick=round.([xmin, 0, xmax], digits=0))
	plot!(ytick=round.([ymin, 0, ymax], digits=2))

	θ = LinRange(0, 2*pi, 1001)
	rfov = sg.rfov
	scat!([0], [0], marker=:circle, label="")
	plot!(rfov * cos.(θ), rfov * sin.(θ), color=:magenta, label="") # rfov circle
	rfov = round(sg.rfov, digits=1)
	plot!(xlabel="x", ylabel="y", title = "$(typeof(sg))")

#=
	if sg isa SinoPar
	end
=#

	if sg isa SinoFan
		x0 = 0
		y0 = sg.dso
		t = LinRange(0, 2π, 100)
		rot = sg.ar[1]
		rot = [cos(rot) -sin(rot); sin(rot) cos(rot)]
		p0 = rot * [x0; y0]
		pd = rot * [sg.xds'; sg.yds'] # detector points

		tmp = sg.ar .+ π/2 # trick: angle beta defined ccw from y axis
		scat!([p0[1]], [p0[2]], color=:blue, marker=:square, label="") # source
		plot!(sg.dso * cos.(t), sg.dso * sin.(t), color=:cyan, label="") # source circle
		scat!(sg.dso * cos.(tmp), sg.dso * sin.(tmp),
			color=:blue, marker=:circle, markersize=2, label="") # source points
		scat!(vec(pd[1,:]), vec(pd[2,:]), marker=:circle,
			color=:orange, markersize=1, label="") # detectors

		plot!([pd[1,1], p0[1], pd[1,end]], [pd[2,1], p0[2], pd[2,end]],
			color=:red, label="")
		plot!(title="$(typeof(sg))")
	end

	if sg isa SinoMoj
		θ = LinRange(0, 2π, 100)
		rphi = sg.nb/2 * sg.d_moj.(θ)
		plot!(rphi .* cos.(θ), rphi .* sin.(θ), color=:blue, label="")
	#	rmax = maximum(sg.s)
	#	axis([-1 1 -1 1] * max([rmax ig.fov/2]) * 1.1)
	end

	plot!()
end


text:
module Strawberry

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Strawberry" => Ahorn.EntityPlacement(
        Maple.Strawberry
    ),
    "Golden Strawberry" => Ahorn.EntityPlacement(
        Maple.GoldenStrawberry
    ),
    "Space Berry" => Ahorn.EntityPlacement(
        Maple.Strawberry,
        "point",
        Dict{String, Any}(
            "moon" => true
        )
    ),

    "Strawberry (Winged)" => Ahorn.EntityPlacement(
        Maple.Strawberry,
        "point",
        Dict{String, Any}(
            "winged" => true
        )
    ),
    "Golden Strawberry (Winged)" => Ahorn.EntityPlacement(
        Maple.GoldenStrawberry,
        "point",
        Dict{String, Any}(
            "winged" => true
        )
    ),

    "Golden Strawberry (Winged, No Dash)" => Ahorn.EntityPlacement(
        Maple.GoldenStrawberryNoDash,
        "point",
    ),
)

const strawberryUnion = Union{Maple.Strawberry, Maple.GoldenStrawberry, Maple.GoldenStrawberryNoDash}

# name, winged, has pips, moon
const sprites = Dict{Tuple{String, Bool, Bool, Bool}, String}(
    ("strawberry", false, false, false) => "collectables/strawberry/normal00",
    ("strawberry", true, false, false) => "collectables/strawberry/wings01",
    ("strawberry", false, true, false) => "collectables/ghostberry/idle00",
    ("strawberry", true, true, false) => "collectables/ghostberry/wings01",

    ("strawberry", false, false, true) => "collectables/moonBerry/normal00",
    ("strawberry", true, false, true) => "collectables/moonBerry/ghost00",
    ("strawberry", false, true, true) => "collectables/moonBerry/ghost00",
    ("strawberry", true, true, true) => "collectables/moonBerry/ghost00",

    ("goldenBerry", false, false, false) => "collectables/goldberry/idle00",
    ("goldenBerry", true, false, false) => "collectables/goldberry/wings01",
    ("goldenBerry", false, true, false) => "collectables/ghostgoldberry/idle00",
    ("goldenBerry", true, true, false) => "collectables/ghostgoldberry/wings01",

    ("memorialTextController", true, false, false) => "collectables/goldberry/wings01",
    ("memorialTextController", true, true, false) => "collectables/goldberry/wings01",
)

const seeds = Dict{String, String}(
    "strawberry" => "collectables/strawberry/seed00",
    "goldenBerry" => "collectables/goldberry/seed00",
    "memorialTextController" => "collectables/goldberry/seed00",
)

const fallback = "collectables/strawberry/normal00"

Ahorn.nodeLimits(entity::strawberryUnion) = 0, -1

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

    nodes = get(entity.data, "nodes", ())
    moon = get(entity.data, "moon", false)
    winged = get(entity.data, "winged", false) || entity.name == "memorialTextController"
    hasPips = length(nodes) > 0

    sprite = sprites[(entity.name, winged, hasPips, moon)]
    seedSprite = seeds[entity.name]

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = node

        push!(res, Ahorn.getSpriteRectangle(seedSprite, nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::strawberryUnion)
    x, y = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = node

        Ahorn.drawLines(ctx, Tuple{Number, Number}[(x, y), (nx, ny)], Ahorn.colors.selection_selected_fc)
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::strawberryUnion, room::Maple.Room)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())
    moon = get(entity.data, "moon", false)
    winged = get(entity.data, "winged", false) || entity.name == "memorialTextController"
    hasPips = length(nodes) > 0

    sprite = sprites[(entity.name, winged, hasPips, moon)]
    seedSprite = seeds[entity.name]

    for node in nodes
        nx, ny = node

        Ahorn.drawSprite(ctx, seedSprite, nx, ny)
    end

    Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
"""
    logit10(x::Float64)
Function transforming interval [0,1] to [-Inf, Inf] using logit transformation.
"""
function logit10(x::Float64)
    log10(x / (1.0 - x))
end

"""
    logistic10(x::Float64)
Function transforming interval [-Inf, Inf] to [0,1] using logistic transformation.
Inverse function for [`logit10`](@ref).
"""
function logistic10(x::Float64)
    exp10(x) / (exp10(x) + 1.0)
end

"""
    scaling(x::Float64, scale::Symbol = :direct)
Transforms values from specific scale to range [-Inf, Inf] based on option.

## Return
Transformed value.

## Arguments
* `x`: input value.
* `scale`: transformation type: `:direct, :log, :logit`.
"""
function scaling(x::Float64, scale::Symbol = :direct)
    if scale == :direct
        return x
    elseif scale == :log
        return log10(x)
    elseif scale == :logit
        return logit10(x)
    else
        throw(DomainError(scale, "scale type is not supported"))
    end
end

"""
    unscaling(x::Float64, scale::Symbol = :direct)
Transforms values from [-Inf, Inf] to specific scale based on option. Inverse function
for [`scaling`](@ref).

## Return
Transformed value.

## Arguments
* `x`: input value.
* `scale`: transformation type: `:direct, :log, :logit`.
"""
function unscaling(x::Float64, scale::Symbol = :direct)
    if scale == :direct
        return x
    elseif scale == :log
        return exp10(x)
    elseif scale == :logit
        return logistic10(x)
    else
        throw(DomainError(scale, "scale type is not supported"))
    end
end

scaling(x_tup::Tuple{Float64,Float64}, scale::Symbol = :direct) = scaling.(x_tup, scale)

unscaling(x_tup::Tuple{Float64,Float64}, scale::Symbol = :direct) = unscaling.(x_tup, scale)

scaling(::Nothing, ::Symbol) = nothing

unscaling(::Nothing, ::Symbol) = nothing

"""
    function get_endpoint(
        theta_init::Vector{Float64},
        theta_num::Int,
        loss_func::Function,
        method::Symbol,
        direction::Symbol = :right;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bound::Float64 = unscaling(
            (direction==:left) ? -9.0 : 9.0,
            scale[theta_num]
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 1e-3,
        local_alg::Symbol = :LN_NELDERMEAD,
        kwargs...
        )

Calculates right or left endpoint of CI for parameter component. It is a wripper
of `get_right_endpoint` functions for selection of direction and using different
transformations for faster optimization.

## Return
[`EndPoint`](@ref) object storing confidence endpoint and profile points found on fly.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values is not necessary to be the optimum values for `loss_func` but it the value of `loss_func` must be lower than `loss_crit`.
- `theta_num`: number ``n`` of vector component to compute confidence interval ``\\theta^n``.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` the profile of which is analyzed. Usually we use log-likelihood for profile analysis in form ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to evaluate interval endpoint. Currently the following methods are implemented: `:CICO_ONE_PASS`, `:LIN_EXTRAPOL`, `:QUADR_EXTRAPOL`.
- `direction`: `:right` or `:left` endpoint to estimate.

## Keyword arguments
see [`get_interval`](@ref)

"""
function get_endpoint(
    theta_init::Vector{Float64},
    theta_num::Int,
    loss_func::Function,
    method::Symbol,
    direction::Symbol = :right;

    loss_crit::Float64 = 0.0,
    # :direct, :log, :logit
    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_bound::Float64 = unscaling(
        (direction==:left) ? -9.0 : 9.0,
        scale[theta_num]
        ),
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3,
    local_alg::Symbol = :LN_NELDERMEAD,
    kwargs... # other options for get_right_endpoint
    )
    isLeft = direction == :left

    # checking arguments
    # theta_bound[1] < theta_init < theta_bound[2]
    theta_init_outside_theta_bounds = .! [theta_bounds[i][1] < theta_init[i] < theta_bounds[i][2] for i in 1:length(theta_init)]
    if any(theta_init_outside_theta_bounds)
        throw(ArgumentError("theta_init is outside theta_bound: $(findall(theta_init_outside_theta_bounds))"))
    end
    # scan_bound should be within theta_bounds
    !(theta_bounds[theta_num][1] < scan_bound < theta_bounds[theta_num][2]) &&
        throw(ArgumentError("scan_bound are outside of the theta_bounds $(theta_bounds[theta_num])"))
    # theta_init should be within scan_bound
    if (theta_init[theta_num] >= scan_bound && !isLeft) || (theta_init[theta_num] <= scan_bound && isLeft)
        throw(ArgumentError("init values are outside of the scan_bound $scan_bound"))
    end
    # 0 <= theta_bound[1] for :log
    less_than_zero_theta_bounds = (scale .== :log) .& [theta_bounds[i][1] < 0 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":log scaled theta_bound min is negative: $(findall(less_than_zero_theta_bounds))"))
    end
    # 0 <= theta_bounds <= 1 for :logit
    less_than_zero_theta_bounds = (scale .== :logit) .& [theta_bounds[i][1] < 0 || theta_bounds[i][2] > 1 for i in 1:length(theta_init)]
    if any(less_than_zero_theta_bounds)
        throw(ArgumentError(":logit scaled theta_bound min is outside range [0,1]: $(findall(less_than_zero_theta_bounds))"))
    end
    # loss_func(theta_init) < loss_crit
    !(loss_func(theta_init) < loss_crit) &&
        throw(ArgumentError("Check theta_init and loss_crit: loss_func(theta_init) should be < loss_crit"))

    # set counter in the scope
    counter::Int = 0
    # set supreme, maximal or minimal value of scanned parameter inside critical
    supreme_gd = nothing

    # transforming
    theta_init_gd = scaling.(theta_init, scale)
    if isLeft theta_init_gd[theta_num] *= -1 end # change direction
    function loss_func_gd(theta_gd::Vector{Float64})
        theta_g = copy(theta_gd)
        if isLeft theta_g[theta_num] *= -1 end # change direction
        theta = unscaling.(theta_g, scale)
        # calculate function
        loss_norm = loss_func(theta) - loss_crit

        # update counter
        counter += 1
        # update supreme ?
        update_supreme = (loss_norm < 0.) &&
            (typeof(supreme_gd)==Nothing || (theta_gd[theta_num] > supreme_gd))
        if update_supreme
            supreme_gd = theta_gd[theta_num]
        end

        return loss_norm
    end
    theta_bounds_gd = scaling.(theta_bounds, scale)
    if isLeft theta_bounds_gd[theta_num] = (-1*theta_bounds_gd[theta_num][2], -1*theta_bounds_gd[theta_num][1]) end # change direction
    scan_bound_gd = scaling(scan_bound, scale[theta_num])
    if isLeft scan_bound_gd *= -1 end # change direction

    # calculate endpoint using base method
    (optf_gd, pp_gd, status) = get_right_endpoint(
        theta_init_gd,
        theta_num,
        loss_func_gd,
        Val(method);
        theta_bounds = theta_bounds_gd,
        scan_bound = scan_bound_gd,
        scan_tol = scan_tol,
        loss_tol = loss_tol,
        local_alg = local_alg,
        kwargs... # options for local fitter
    )

    # transforming back
    if (isLeft && typeof(optf_gd)!==Nothing) optf_gd *= -1 end # change direction
    optf = unscaling(optf_gd, scale[theta_num])
    temp_fun = (pp::ProfilePoint) -> begin
        if isLeft pp.params[theta_num] *= -1 end # change direction
        ProfilePoint(
            unscaling(pp.params[theta_num], scale[theta_num]),
            pp.loss + loss_crit,
            unscaling.(pp.params, scale),
            pp.ret,
            pp.counter
        )
    end
    pps = [ temp_fun(pp_gd[i]) for i in 1:length(pp_gd) ]
    # transforming supreme back
    if (isLeft && typeof(supreme_gd)!==Nothing) supreme_gd *= -1 end # change direction
    supreme = unscaling(supreme_gd, scale[theta_num])

    EndPoint(optf, pps, status, direction, counter, supreme)
end

"""
    function get_right_endpoint(
        theta_init::Vector{Float64},
        theta_num::Int,
        loss_func::Function,
        method::Val{:CICO_ONE_PASS};

        theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
            (-Inf, Inf), length(theta_init)
            ),
        scan_bound::Float64 = 9.0,
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 0.,
        local_alg::Symbol = :LN_NELDERMEAD,
        kwargs...
        )
Interface for current and future methods for endpoint estimation.

## Return
Tuple of three values:
- Right end point value: `::Float64`.
- Profile points estimated on fly: `::Array{ ProfilePoint, 1}`, see [`ProfilePoint`](@ref).
- Status of sulution: `::Symbol`. One of values: `:BORDER_FOUND_BY_SCAN_TOL`, `:SCAN_BOUND_REACHED`.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values is not necessary to be the optimum values for `loss_func` but it the value of `loss_func` must be lower than `loss_crit`.
- `theta_num`: number ``n`` of vector component to compute confidence interval ``\\theta^n``.
- `loss_func`: loss function the profile of which is analyzed, see [`get_interval`](@ref). In this `function loss` crit is always equal 0 for code simplification.
- `method`: this value is always fixed. Implemented methods are: `Val{:CICO_ONE_PASS}`. It is implemented for easy switching between different implemented and future methods.

## Keyword arguments
- `theta_bound`: vector of bounds for each component in format `(left_bound, right_bound)`. This bounds define the ranges for possible parameter values.
- `scan_bound`: right scan bound for `theta_num` component. It must be within the `theta_bounds` for the scanned component.
- `scan_tol`: Absolute tolerance of scanned component (stop criterion).
- `loss_tol`: Absolute tolerance of `loss_func` at `loss_crit` (stop criterion). *Restriction*. Currently is not effective for `:CICO_ONE_PASS` methods because of limitation in `LN_AUGLAG` interface.
- `local_alg`: algorithm of optimization. Currently the local derivation free algorithms form NLOPT pack were tested. The methods: `:LN_NELDERMEAD, :LN_COBYLA, :LN_PRAXIS` show good results. Methods: `:LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA` is not recommended.
- `kwargs...`: the additional keyword arguments passed to `get_right_endpoint` for specific `method`.
"""
function get_right_endpoint
end


text:
# First, we define a nice unicode summation function
∑ = sum

# then, make an ℓ2 norm function
ℓ2(θ::Vector) = √θ'θ
ℓ2([-0.5,0,1,2,3,4,5,6])

# Let's make an ℓp norm function with methods for ints, rationals and floats
ℓp(θ::Vector, p::Int64) = ∑(abs.(θ).^p)^(1//p)
ℓp(θ::Vector, p::Rational) = ∑(abs.(θ).^p)^(1//p)
ℓp(θ::Vector, p::Float64) = ∑(abs.(θ).^p)^(1/p)

# Now let's compare their precision
𝑥 = ℓp([-0.5,0,1,2,3,4,5,6], 122/3)
𝑦 = ℓp([-0.5,0,1,2,3,4,5,6], 122//3)

# 𝑥 and 𝑦 differ on the last digit. I trust 𝑦 more because it uses rationals!
𝑥 == 𝑦
𝑥 ≈ 𝑦

max_stars_repo_name: vankesteren/julia-tryout

max_issues_repo_name: vankesteren/julia-tryout

max_forks_repo_name: vankesteren/julia-tryout


text:
mutable struct Primitive
    attributes::Dict{String,Int}
    indices::Union{Nothing,Int}
    material::Union{Nothing,Int}
    mode::Union{Nothing,Int}
    targets::Union{Nothing,Vector{Dict}}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Primitive()
        obj = new()
        obj.indices = nothing
        obj.material = nothing
        obj.mode = nothing
        obj.targets = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::Primitive, sym::Symbol)
    x = getfield(obj, sym)
    sym === :mode && x === nothing && return 4
    return x
end

function Base.setproperty!(obj::Primitive, sym::Symbol, x)
    if sym === :indices && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the accessor that contains mesh indices should be ≥ 0"))
    elseif sym === :material && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the material to apply to this primitive when rendering should be ≥ 0"))
    elseif sym === :mode && x !== nothing
        x === POINTS || x === LINES || x === LINE_LOOP || x === LINE_STRIP ||
            x === TRIANGLES || x === TRIANGLE_STRIP || x === TRIANGLE_FAN ||
                throw(ArgumentError("the type of primitives to render should be one of: POINTS(0), LINES(1), LINE_LOOP(2), LINE_STRIP(3), TRIANGLES(4), TRIANGLE_STRIP(5), TRIANGLE_FAN(6)"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Primitive}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Primitive}) = (:indices, :material, :mode, :targets, :extensions, :extras)

mutable struct Mesh
    primitives::ZVector{Primitive}
    weights::Union{Nothing,Vector{Cfloat}}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Mesh()
        obj = new()
        obj.weights = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

JSON3.StructType(::Type{Mesh}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Mesh}) = (:weights, :name, :extensions, :extras)


text:
# Entry for script
function main_rbc_2_joint(args=ARGS)
    d = parse_commandline_rbc_2_joint(args)
    return estimate_rbc_2_joint((; d...)) # to named tuple
end

function estimate_rbc_2_joint(d)
    # Or move these into main package when loading?
    Turing.setadbackend(:zygote)
    HMCExamples.set_BLAS_threads()
    use_tensorboard = true # could add toggle later

    # load data relative to the current path
    data_path = joinpath(pkgdir(HMCExamples), d.data_path)
    df = Matrix(DataFrame(CSV.File(data_path)))
    z = [df[i, :] for i in 1:size(df, 1)]
    ϵ0 = Matrix(DataFrame(CSV.File(joinpath(pkgdir(HMCExamples), "data/epsilons_burnin_rbc_2.csv");header=false)))
    # Create the perturbation and the turing models
    m = SecondOrderPerturbationModel(rbc_2)
    turing_model = rbc_second(
        z, m, d.p_f, d.alpha_prior, d.beta_prior, d.rho_prior, allocate_cache(m)
    )

    # Sampler
    name = "rbc-second-joint-s$(d.num_samples)-seed$(d.seed)"
    include_vars = ["α", "β_draw", "ρ"]  # variables to log
    callback = TensorBoardCallback(d.results_path; name, include=include_vars)
    num_adapts = convert(Int64, floor(d.num_samples * d.adapts_burnin_prop))

    Random.seed!(d.seed)
    @info "Generating $(d.num_samples) samples with $(num_adapts) adapts"
    alg = NUTS(num_adapts, d.target_acceptance_rate)

    chain = sample(
        turing_model,
        NUTS(num_adapts, d.target_acceptance_rate; d.max_depth),
        d.num_samples;
        init_params=[d.p; ϵ0],
        progress=true,
        save_state=true,
        callback,
    )

    # Calculate and save results into the logdir
    calculate_experiment_results(chain, callback.logger, include_vars, d.full_results)

    # Store parameters in log directory
    parameter_save_path = joinpath(callback.logger.logdir, "parameters.json")

    @info "Storing Parameters at $(parameter_save_path) "
    open(parameter_save_path, "w") do f
        write(f, JSON.json(d))
    end
end

function parse_commandline_rbc_2_joint(args)
    s = ArgParseSettings(; fromfile_prefix_chars=['@'])

    # See the appropriate _defaults.txt file for the default vvalues.
    @add_arg_table! s begin
        "--data_path"
        help = "relative path to data from the root of the package"
        arg_type = String
        "--p"
        help = "Initialization of parameters"
        arg_type = Vector{Float64}
        "--p_f"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}
        "--alpha_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--beta_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--num_samples"
        help = "samples to draw in chain"
        arg_type = Int64
        "--adapts_burnin_prop"
        help = "Proportion of Adaptations burned in"
        arg_type = Float64
        "--target_acceptance_rate"
        help = "Target acceptance rate for dual averaging for NUTS"
        arg_type = Float64
        "--max_depth"
        help = "Maximum depth for NUTS"
        arg_type = Int64
        "--seed"
        help = "Random number seed"
        arg_type = Int64
        "--results_path"
        arg_type = String
        help = "Location to store results and logs"
        "--full_results"
        arg_type = Bool
        help = "Save the complete set of figures and results for the chain"
    end

    args_with_default = vcat("@$(pkgdir(HMCExamples))/src/rbc_2_joint_defaults.txt", args)
    return parse_args(args_with_default, s; as_symbols=true)
end


text:
"""
    simple_interest(principal::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Simple Interest.
Simple Interest = principal × interest_rate × number_period
returns a AbstractFloat
# Example
```
julia> simple_interest(100.00 , 0.1 ,2)
=20.0
```
"""
function simple_interest(principal::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return principal *  interest_rate * number_period
end

"""
    compound_interest(principal::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Compound Interest.
Compound Interest = (principal × (1+ interest_rate)^ number_period) - principal
returns a AbstractFloat
# Example
```
julia> compound_interest(100.00 , 0.1 ,2)
=20.0
```
"""
function compound_interest(principal::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return future_value(principal,interest_rate,number_period,1) - principal
end

"""
    future_value(present_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Future Value of Money.
future_value = present_value * (1 + interest_rate)^number_period
returns a AbstractFloat
# Example
```
julia> future_value(100.00 , 0.1 ,2)
=121.00000000000001
```
"""
function future_value(present_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return present_value * (1 + interest_rate)^number_period
end

"""
    present_value(present_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Present Value of Money.
present_value = future_value * (1 + interest_rate)^-number_period
returns a AbstractFloat
# Example
```
julia> present_value(121.00000000000001 , 0.1 ,2)
=100.00
```
"""
function present_value(future_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return future_value * (1 + interest_rate)^-number_period
end

"""
    present_value_interest_factor(interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Present Value Interest Factor.
It is a factor to determine PV of a given FV with interest rate and period
present_value_interest_factor = 1/((1 + interest_rate)^number_period)
returns a AbstractFloat
# Example
```
julia> present_value_interest_factor(0.05 ,5)
=0.7835261664684589
```
"""
function present_value_interest_factor(interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return 1/((1 + interest_rate)^number_period)
end

"""
    future_value(present_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64,compound_period::Int64 )::AbstractFloat
Future Value of Money with compounting period.
future_value = present_value * (1 + interest_rate / compound_period)^(number_period * compound_period)
returns a AbstractFloat
# Example
```
julia> future_value(100.00 , 0.1 ,2, 2)
=121.55062500000003
```
"""
function future_value(present_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64,compound_period::Int64 )::AbstractFloat
    return present_value * (1 + interest_rate / compound_period)^(number_period * compound_period)
end

"""
    future_value_interest_factor(interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
Future Value Interest Factor.
It is a factor to determine FV of a given PV with interest rate and period
future_value_interest_factor = (1 + interest_rate)^number_period)
returns a AbstractFloat
# Example
```
julia> future_value_interest_factor(0.05 ,5)
=1.2762815625000004
```
"""
function future_value_interest_factor(interest_rate::AbstractFloat, number_period::Int64 )::AbstractFloat
    return (1 + interest_rate)^number_period
end

"""
    present_value(future_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64, compound_period::Int64)::AbstractFloat
Present Value of Money with compounting period.
present_value = present_value * (1 + interest_rate / compound_period)^(number_period * compound_period)
returns a AbstractFloat
# Example
```
julia> present_value(121.00000000000001 , 0.1 ,2)
=100.00
```
"""
function present_value(future_value::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64, compound_period::Int64)::AbstractFloat
    return future_value * (1 + interest_rate / compound_period)^-(number_period * compound_period)
end

"""
    present_annuity(cash_flow_each_pay::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64)::AbstractFloat
An annuity is a series of payments made at equal intervals.
Examples of annuities are regular deposits to a savings account,
monthly home mortgage payments, etc.
returns a AbstractFloat
# Example
```
julia> present_annuity(50000.00 , 0.06 ,25)
=639167.8079134206
```
"""
function present_annuity(cash_flow_each_pay::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64)::AbstractFloat
    return cash_flow_each_pay * (1-(1/(1+interest_rate)^number_period))/interest_rate

end

"""
    present_annuity_due(cash_flow_each_pay::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64)::AbstractFloat
Annuity due: in which payments are made at the beginning of each period.
returns a AbstractFloat
# Example
```
julia> present_annuity_due(50000.00 , 0.06 ,25)
=677517.8763882258
```
"""
function present_annuity_due(cash_flow_each_pay::AbstractFloat, interest_rate::AbstractFloat, number_period::Int64)::AbstractFloat
    return cash_flow_each_pay * ((1-(1/(1+interest_rate)^number_period))/interest_rate)* (1 + interest_rate)
end


text:
using RunMetrics

data_dir = joinpath(@__DIR__, "..", "sample_data")

strava_zip = joinpath(data_dir, "stravaSmall.zip")
strava_dest = joinpath(data_dir, "stravaFiles")

garmin_zip = joinpath(data_dir, "garminSmall.zip")
garmin_dest = joinpath(data_dir, "garminFiles")

import_strava(strava_zip, strava_dest)
import_garmin(garmin_zip, garmin_dest)


text:
addprocs(2)

import ProgressMeter
import Base.Test.@test
import Base.Test.@test_throws

using Compat

srand(123)

include("test.jl")
include("test_showvalues.jl")
include("test_pmap.jl")

println("")
println("All tests complete")


text:
__precompile__()

module InfiniteArrays
   using Base, Statistics, LinearAlgebra, FillArrays, LazyArrays, DSP

import Base: *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
            fld, cld, div, min, max, minimum, maximum, mod,
            <, ≤, >, ≥, promote_rule, convert, copy,
            size, step, isempty, length, first, last,
            getindex, setindex!, intersect, @_inline_meta,
            sort, sort!, issorted, sortperm, sum, in, broadcast,
            eltype, parent, real, imag,
            conj, transpose,
            exp, log, sqrt,
                      cos, sin, tan, csc, sec, cot,
                      cosh, sinh, tanh, csch, sech, coth,
                      acos, asin, atan, acsc, asec, acot,
                      acosh, asinh, atanh, acsch, asech, acoth, (:),
            AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo,
           to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
         checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
         _in_range, _range, _rangestyle, Ordered,
         ArithmeticWraps, floatrange, reverse, unitrange_last,
         AbstractArray, AbstractVector, Array, Vector, Matrix,
         axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
         cumsum, show_delim_array, show_circular, Int

using Base.Broadcast
import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        @nexprs, @ncall, combine_eltypes, DefaultArrayStyle, instantiate

import LinearAlgebra: BlasInt, BlasFloat, norm, diag, diagm, ishermitian, issymmetric,
                             det, logdet, istriu, istril, adjoint, tr

import Statistics: mean, median

import FillArrays: AbstractFill, getindex_value
import LazyArrays: LazyArrayStyle, _materialize, ArrayMulArray, AbstractBandedLayout,
                    ZerosLayout, VcatLayout, MatMulVec

import DSP: conv

export ∞, Hcat, Vcat, Zeros, Ones, Fill, Eye, BroadcastArray, cache

include("Infinity.jl")
include("infrange.jl")
include("infarrays.jl")

##
# Temporary hacks for base support
##
OneTo(::Infinity) = OneToInf()
OneTo{T}(::Infinity) where T<:Integer = OneToInf{T}()
UnitRange(start::Integer, ::Infinity) = InfUnitRange(start)
UnitRange{T}(start::Integer, ::Infinity) where T<:Real = InfUnitRange{T}(start)

Int(::Infinity) = ∞

# stay lazy if infinite
_materialize(M::ArrayMulArray, ::Tuple{<:OneToInf}) = M
_materialize(M::ArrayMulArray, ::Tuple{<:OneToInf,<:OneToInf}) = M
_materialize(M::ArrayMulArray, ::Tuple{<:OneToInf,<:OneTo}) = M
_materialize(M::ArrayMulArray, ::Tuple{<:OneTo,<:OneToInf}) = M

_materialize(M::MatMulVec{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}}, ::Tuple{<:OneToInf}) =
    copyto!(similar(M), M)

end # module


text:
struct MetaProp
    data::Vector{ConTab}
    metric::Symbol
    y::Vector
    var::Vector
end

struct MetaPropResult{Symbol}
    data::MetaProp
    wts::Vector{Float64}
    est::Float64
    var::Float64
    chisq::Float64
    hetq::Float64
    heti::Float64
    hettau::Float64
end

"""
    metaprop(d, metric; adj = 0)

Meta-analysis for 2x2 tables.

`ds`:

`metric`:
- :rr
- :or
- :diff
"""
function metaprop(d::DataSet, metric; adj = 0)
    if metric == :diff
        cty = contabdiff.(d.ds; adj = adj)
    elseif metric == :or
        cty = contabor.(d.ds; adj = adj)
    elseif metric == :rr
        cty = contabrr.(d.ds; adj = adj)
    else
    end
    y   = getindex.(cty, 1)
    var = getindex.(cty, 2)
    MetaProp(d.ds, metric, y, var)
end

"""
    metapropfixed(mp; weights = :default)

Inverce Variance method used by default.

`weights`:

- `:iv` | `:default`
- `:mh`

"""
function metapropfixed(mp; weights = :default)
    varwts    = 1 ./ mp.var
    if weights == :default || weights == :iv
        wts    = varwts
    elseif weights == :mh
        if mp.metric == :diff
            wts    = mhwdiff(mp.data)
        elseif mp.metric == :or
            wts    = mhwor(mp.data)
        elseif mp.metric == :rr
            wts    = mhwrr(mp.data)
        end
    #elseif weights == :peto
    else
        error("weights keyword unknown!")
    end
    k       = length(mp.y)
    var     = 1 / sum(varwts)
    est     = sum(wts .* mp.y) / sum(wts)
    chisq   = sum(varwts .* (mp.y .^ 2))
    q       = sum(varwts .* ((mp.y .- est) .^ 2))
    i²      = max(0, (q - (k - 1))/q * 100)
    MetaPropResult{:fixed}(mp, wts, est, var,  chisq, q, i², NaN)
end

"""
    metaproprandom(mp; tau = :default)

tau - τ² calculation method:
- `:dl` DerSimonian-Laird
- `:ho`
- `:hm` Hartung and Makambi
- `:sj` Sidik and Jonkman

"""
function metaproprandom(mp; tau = :default)
    k       = length(mp.y)
    wts     = 1 ./ mp.var
    swts    = sum(wts)
    var     = 1 / swts
    est     = sum(wts .* mp.y) / swts
    chisq   = sum(wts .* (mp.y .^ 2))
    q       = sum(wts .* ((mp.y .- est) .^ 2))
    i²      = max(0, (q - (k - 1))/q * 100)
    s       = sum(wts) - sum(wts .^ 2) / sum(wts)
    if tau == :dl || tau == :default # DerSimonian-Laird
        τ² = max(0, (q - (k - 1)) / s)
    elseif tau == :ho
        τ² = max(0, 1 / (k - 1) * sum((mp.y .- mean(mp.y)) .^ 2) - 1 / k * sum(mp.var))
    elseif tau == :hm # Hartung and Makambi
        τ² = q ^ 2 / (2 * (k - 1) + q) / (sum(wts) - sum(wts .^ 2) / sum(wts))
    elseif tau == :sj # Sidik and Jonkman
        qi  =  1 ./ (mp.var ./ max(0.01, 1 / (k - 1) * sum((mp.y .- mean(mp.y)) .^ 2) - 1 / k * sum(mp.var)) .+ 1)
        rest = sum(qi .* mp.y) / sum(qi)
        τ²   = 1 / (k - 1) * sum(qi .* ((mp.y .- rest) .^ 2))
    else
        error("tau keyword unknown!")
    end
    rwts = 1 ./ (mp.var .+ τ²)
    est     = sum(rwts .* mp.y) / sum(rwts)
    var     = 1 / sum(rwts) #?
    i²      = (τ² / (τ² + (k-1)/s)) * 100
    MetaPropResult{:random}(mp, rwts, est, var, chisq, q, i², τ²)
end

"""
    StatsBase.confint(mpr::MetaPropResult; level = 0.95)

Confidence interval for pooled proportion.
"""
function StatsBase.confint(mpr::MetaPropResult; level = 0.95)
    alpha = 1 - level
    se = sqrt(mpr.var)
    d = quantile(Normal(), 1 - alpha / 2) * se
    mpr.est - d, mpr.est + d
end

function mhwdiff(data)
    wts = Vector{Float64}(undef, length(data))
    for i = 1:length(data)
        wts[i] = (data[i].tab[1,1] + data[i].tab[1,2]) * (data[i].tab[2,1] + data[i].tab[2,2] ) / (data[i].tab[1,1] + data[i].tab[1,2] + data[i].tab[2,1] + data[i].tab[2,2])
    end
    wts
end
function mhwor(data)
    wts = Vector{Float64}(undef, length(data))
    for i = 1:length(data)
        wts[i] = data[i].tab[1,2] * data[i].tab[2,1] / (data[i].tab[1,1] + data[i].tab[1,2] + data[i].tab[2,1] + data[i].tab[2,2])
    end
    wts
end
function mhwrr(data)
    wts = Vector{Float64}(undef, length(data))
    for i = 1:length(data)
        wts[i] = (data[i].tab[1,1] + data[i].tab[1,2]) * data[i].tab[2,1] / (data[i].tab[1,1] + data[i].tab[1,2] + data[i].tab[2,1] + data[i].tab[2,2])
    end
    wts
end

#function petow(data)
#end
function contabdiff(contab; adj = 0)
    a = contab.tab[1,1] + adj
    b = contab.tab[1,2] + adj
    c = contab.tab[2,1] + adj
    d = contab.tab[2,2] + adj
    pt = a / (a + b)
    pc = c / (c + d)
    return  pt - pc, pt * (1 - pt) / (a + b) + pc * (1 - pc) / (c + d)
end
function contabor(contab; adj = 0)
    a = contab.tab[1,1] + adj
    b = contab.tab[1,2] + adj
    c = contab.tab[2,1] + adj
    d = contab.tab[2,2] + adj
    pt = a / (a + b)
    pc = c / (c + d)
    return log((pt / (1 - pt)) / (pc / (1 - pc))), 1 / a + 1 / b + 1 / c + 1 / d
end
function contabrr(contab; adj = 0)
    a = contab.tab[1,1] + adj
    b = contab.tab[1,2] + adj
    c = contab.tab[2,1] + adj
    d = contab.tab[2,2] + adj
    pt = a / (a + b)
    pc = c / (c + d)
    return log(pt / pc), 1 / a - 1 / (a + b) + 1 / c - 1 /(c + d)
end


text:
# Before running this, please make sure to activate and instantiate the
# tutorial-specific package environment, using this
# [`Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-horse/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-horse/Manifest.toml), or by following
# [these](https://juliaai.github.io/DataScienceTutorials.jl/#learning_by_doing) detailed instructions.

# ## Initial data processing
#
# In this example, we consider the [UCI "horse colic" dataset](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)
#
# This is a reasonably messy classification problem with missing values etc and so some work should be expected in the feature processing.
#
# ### Getting the data
#
# The data is pre-split in training and testing and we will keep it as such

using MLJ
using HTTP
using CSV
import DataFrames: DataFrame, select!, Not
req1 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.data")
req2 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.test")
header = ["surgery", "age", "hospital_number",
    "rectal_temperature", "pulse",
    "respiratory_rate", "temperature_extremities",
    "peripheral_pulse", "mucous_membranes",
    "capillary_refill_time", "pain",
    "peristalsis", "abdominal_distension",
    "nasogastric_tube", "nasogastric_reflux",
    "nasogastric_reflux_ph", "feces", "abdomen",
    "packed_cell_volume", "total_protein",
    "abdomcentesis_appearance", "abdomcentesis_total_protein",
    "outcome", "surgical_lesion", "lesion_1", "lesion_2", "lesion_3",
    "cp_data"]
csv_opts = (header=header, delim=' ', missingstring="?",
            ignorerepeated=true)
data_train = CSV.read(req1.body, DataFrame; csv_opts...)
data_test  = CSV.read(req2.body, DataFrame; csv_opts...)
@show size(data_train)
@show size(data_test)

# ### Inspecting columns
#
# To simplify the analysis, we will drop the columns `Lesion *` as they would need specific re-encoding which would distract us a bit.

unwanted = [:lesion_1, :lesion_2, :lesion_3]
data = vcat(data_train, data_test)
select!(data, Not(unwanted));

# Let's also keep track of the initial train-test split

train = 1:nrows(data_train)
test = last(train) .+ (1:nrows(data_test));

# We know from reading the description that some of these features represent multiclass data; to facilitate the interpretation, we can use `autotype` from `ScientificTypes`.
# By default, `autotype` will check all columns and suggest a Finite type assuming there are relatively few distinct values in the column.
# More sophisticated rules can be passed, see [ScientificTypes.jl](https://alan-turing-institute.github.io/ScientificTypes.jl/dev/):

datac = coerce(data, autotype(data));

# Let's see column by column whether it looks ok now

sch = schema(datac)
for (name, scitype) in zip(sch.names, sch.scitypes)
    println(rpad("$name", 30), scitype)
end

# Most columns are now treated as either Multiclass or Ordered, this
# corresponds to the [description of the data](https://archive.ics.uci.edu/ml/datasets/Horse+Colic). For instance:
#
# - `Surgery` is described as `1=yes / 2=no`
# - `Age` is described as `1=adult / 2=young`
#
# Inspecting the rest of the descriptions and the current scientific type,
# there are a few more things that can be observed:
#
# - hospital number is still a count, this means that there are relatively many hospitals and so  that's  probably not very useful,
# - pulse and respiratory rate are still as count but the data description suggests that they can be considered as continuous

length(unique(datac.hospital_number))

# yeah let's drop that

datac = select!(datac, Not(:hospital_number));

# let's also coerce the pulse and respiratory rate, in fact we can do that with
# `autotype` specifying as rule the `discrete_to_continuous`

datac = coerce(datac, autotype(datac, rules=(:discrete_to_continuous,)));

# ### Dealing with missing values
#
# There's quite a lot of missing values, in this tutorial we'll be a bit rough in how we deal with them applying the following rules of thumb:
#
# - drop the rows where the outcome is unknown
# - drop columns with more than 20% missing values
# - simple imputation of whatever's left

missing_outcome = ismissing.(datac.outcome)
idx_missing_outcome = missing_outcome |> findall

# Ok there's only two row which is nice, let's remove them from the train/test indices and drop the rows

train = setdiff!(train |> collect, idx_missing_outcome)
test = setdiff!(test |> collect, idx_missing_outcome)
datac = datac[.!missing_outcome, :];

# Now let's look at how many missings there are per features

for name in names(datac)
    col = datac[:, name]
    ratio_missing = sum(ismissing.(col)) / nrows(datac) * 100
    println(rpad(name, 30), round(ratio_missing, sigdigits=3))
end

# Let's drop the ones with more than 20% (quite a few!)

unwanted = [:peripheral_pulse, :nasogastric_tube, :nasogastric_reflux,
        :nasogastric_reflux_ph, :feces, :abdomen, :abdomcentesis_appearance, :abdomcentesis_total_protein]
select!(datac, Not(unwanted));

# Note that we could have done this better and investigated the nature of the features for which there's a lot of missing values but don't forget that our goal is to showcase MLJ!
#
# Let's conclude by filling all missing values and separating the feature matrix from the  target

@load FillImputer
filler = machine(FillImputer(), datac)
fit!(filler)
datac = MLJ.transform(filler, datac)

y, X = unpack(datac, ==(:outcome), name->true);
X = coerce(X, autotype(X, :discrete_to_continuous));

# ## A baseline model
#
# Let's define a first sensible model and get a baseline, basic steps are:
# - one-hot-encode the categoricals
# - feed all this into a classifier

@load OneHotEncoder
MultinomialClassifier = @load MultinomialClassifier pkg="MLJLinearModels"

# Let's have convenient handles over the training data

Xtrain = X[train,:]
ytrain = y[train];

# And let's define a pipeline corresponding to the operations above

SimplePipe = Pipeline(
    OneHotEncoder(),
    MultinomialClassifier(),
    prediction_type=:probabilistic
)
mach = machine(SimplePipe, Xtrain, ytrain)
res = evaluate!(
    mach;
    resampling=Holdout(fraction_train=0.9),
    measure=cross_entropy
)
round(res.measurement[1], sigdigits=3)

# This is the cross entropy on some held-out 10% of the training set.
# We can also just for the sake of getting a baseline, see the misclassification on the whole training data:

mcr = misclassification_rate(predict_mode(mach, Xtrain), ytrain)
println(rpad("MNC mcr:", 10), round(mcr, sigdigits=3))

# That's not bad at all actually.
# Let's tune it a bit and see if we can get a bit better than that, not much point in going crazy, we might get a few percents but not much more.

model = SimplePipe
lambdas = range(model, :(multinomial_classifier.lambda), lower=1e-3, upper=100, scale=:log10)
tm = TunedModel(model=SimplePipe, ranges=lambdas, measure=cross_entropy)
mtm = machine(tm, Xtrain, ytrain)
fit!(mtm)
best_pipe = fitted_params(mtm).best_model

# So it looks like it's useful to regularise a fair bit to get a lower cross entropy

ŷ = MLJ.predict(mtm, Xtrain)
cross_entropy(ŷ, ytrain) |> mean

# Interestingly this does not improve our missclassification rate

mcr = misclassification_rate(mode.(ŷ), ytrain)
println(rpad("MNC mcr:", 10), round(mcr, sigdigits=3))

# We've probably reached the limit of a simple linear model.
#
# ## Trying another model
#
# There are lots of categoricals, so maybe  it's just better to use something that deals well with that like a tree-based classifier.

XGBC = @load XGBoostClassifier
dtc = machine(XGBC(), Xtrain, ytrain)
fit!(dtc)
ŷ = MLJ.predict(dtc, Xtrain)
cross_entropy(ŷ, ytrain) |> mean

# So we get a worse cross entropy but...

misclassification_rate(mode.(ŷ), ytrain)

# a significantly better misclassification rate.
#
# We could investigate more, do more tuning etc, but the key points of this tutorial was to show how to handle data with missing values.

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials