text:
using Flux
using NNlib

struct Embedding
    W
end
Embedding(vocabsize, embedsize) = Embedding(param(randn(embedsize, vocabsize)))
(embedding::Embedding)(x) = embedding.W[:, x]
Flux.treelike(Embedding)

FeedForward(d_model, d_hidden) = Chain(Dense(d_model, d_hidden, relu), Dense(d_hidden, d_model))

struct ResidualBlock
    l
    norm
    # TODO dropout
end
ResidualBlock(l, d_model::Integer) = ResidualBlock(l, LayerNorm(d_model))
(l::ResidualBlock)(x...) = x[1] .+ l.norm(l.l(x...))
Flux.treelike(ResidualBlock)

struct Attention
    scale
    # TODO dropout
    # TODO causal
end
Attention(d_key::Integer, causal=false) = Attention(sqrt(d_key))
function (l::Attention)(q, k, v)
    alpha = k'*q
    # TODO causal mask
    return v*softmax(alpha, dim=1)
end
Flux.mapchildren(f, l::Attention) = Attention(l.scale)

# TODO multihead

struct EncoderLayer
    selfattn
    feedforward
end
EncoderLayer(d_model::Integer, d_hidden::Integer) = EncoderLayer(
    ResidualBlock(Attention(d_model), d_model),
    ResidualBlock(FeedForward(d_model, d_hidden), d_model))
(l::EncoderLayer)(x) = l.feedforward(l.selfattn(x, x, x))
Flux.treelike(EncoderLayer)

struct Encoder
    embed
    layers
end
Encoder(d_model, d_hidden, n_layers, vocabsize) = Encoder(
    Embedding(vocabsize, d_model),
    Chain((EncoderLayer(d_model, d_hidden) for _ in 1:n_layers)...))
function (l::Encoder)(x)
    x = posenc!(l.embed(x))
    encoding = [x = l(x) for l in l.layers.layers]
    return encoding
end
Flux.treelike(Encoder)


text:
# ---
# title: 798. Smallest Rotation with Highest Score
# id: problem798
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories:
# link: <https://leetcode.com/problems/smallest-rotation-with-highest-score/description/>
# hidden: true
# ---
#
#  Given an array `A`, we may rotate it by a non-negative integer `K` so that
# the array becomes `A[K], A[K+1], A{K+2], ... A[A.length - 1], A[0], A[1], ...,
# A[K-1]`.  Afterward, any entries that are less than or equal to their index
# are worth 1 point.
#
# For example, if we have `[2, 4, 1, 3, 0]`, and we rotate by `K = 2`, it
# becomes `[1, 3, 0, 2, 4]`.  This is worth 3 points because 1 > 0 [no points],
# 3 > 1 [no points], 0 <= 2 [one point], 2 <= 3 [one point], 4 <= 4 [one point].
#
# Over all possible rotations, return the rotation index K that corresponds to
# the highest score we could receive.  If there are multiple answers, return the
# smallest such index K.
#
#
#
#     **Example 1:**
#     Input: [2, 3, 1, 4, 0]
#     Output: 3
#     Explanation:
#     Scores for each K are listed below:
#     K = 0,  A = [2,3,1,4,0],    score 2
#     K = 1,  A = [3,1,4,0,2],    score 3
#     K = 2,  A = [1,4,0,2,3],    score 3
#     K = 3,  A = [4,0,2,3,1],    score 4
#     K = 4,  A = [0,2,3,1,4],    score 3
#
#
# So we should choose K = 3, which has the highest score.
#
#
#
#
#
#     **Example 2:**
#     Input: [1, 3, 0, 2, 4]
#     Output: 0
#     Explanation: A will always have 3 points no matter how it shifts.
#     So we will choose the smallest K, which is 0.
#
#
# **Note:**
#
#   * `A` will have length at most `20000`.
#   * `A[i]` will be in the range `[0, A.length]`.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
using Feather, DataFrames, Base.Test, NullableArrays

if VERSION < v"0.5.0-dev+4267"
    @eval is_windows() = $(OS_NAME == :Windows)
end

testdir = joinpath(dirname(@__FILE__), "data")
# testdir = joinpath(Pkg.dir("Feather"), "test/data")
files = map(x -> joinpath(testdir, x), readdir(testdir))

testnull{T}(v1::T, v2::T) = v1 == v2
testnull{T}(v1::Nullable{T}, v2::Nullable{T}) = (isnull(v1) && isnull(v2)) || (!isnull(v1) && !isnull(v2) && get(v1) == get(v2))
testnull{T}(v1::T, v2::Nullable{T}) = !isnull(v2) && get(v2) == v1
testnull{T}(v1::Nullable{T}, v2::T) = !isnull(v1) && get(v1) == v2

for f in files
    source = Feather.Source(f)
    df = Data.stream!(source, DataFrame)
    temp = tempname()
    sink = Feather.Sink(temp)
    Feather.Data.stream!(df, sink)
    Data.close!(sink)
    df2 = Feather.read(temp)

    for (c1,c2) in zip(df.columns,df2.columns)
        for i = 1:length(c1)
            @test testnull(c1[i], c2[i])
        end
    end
    @test Data.header(source) == Data.header(sink) && Data.types(source) == Data.types(sink)
    @test source.ctable.description == sink.ctable.description
    @test source.ctable.num_rows == sink.ctable.num_rows
    @test source.ctable.metadata == sink.ctable.metadata
    for (col1,col2) in zip(source.ctable.columns,sink.ctable.columns)
        @test col1.name == col2.name
        @test col1.metadata_type == col2.metadata_type
        @test typeof(col1.metadata) == typeof(col2.metadata)
        @test col1.user_metadata == col2.user_metadata

        v1 = col1.values; v2 = col2.values
        @test v1.type_ == v2.type_
        @test v1.encoding == v2.encoding
        # @test v1.offset == v2.offset # currently not python/R compatible due to wesm/feather#182
        @test v1.length == v2.length
        @test v1.null_count == v2.null_count
        # @test v1.total_bytes == v2.total_bytes
    end
    rm(temp)
end

# DataStreams interface
source_file = joinpath(testdir, "test_utf8.feather")
sink_file = joinpath(testdir, "test_utf8_new.feather")

ds = Feather.read(source_file)
@test size(ds) == (3,3)
Feather.read(source_file, Feather.Sink, sink_file)
@test isequal(ds, Feather.read(sink_file))

sink = Feather.Sink(joinpath(testdir, "test_utf8_new.feather"))
Feather.read(source_file, sink)
@test isequal(ds, Feather.read(sink_file))

source = Feather.Source(source_file)
ds = Feather.read(source)
@test size(ds) == (3,3)
source = Feather.Source(source_file)
Feather.read(source, Feather.Sink, sink_file)
@test isequal(ds, Feather.read(sink_file))

sink = Feather.Sink(joinpath(testdir, "test_utf8_new.feather"))
source = Feather.Source(source_file)
Feather.read(source, sink)
@test isequal(ds, Feather.read(sink_file))

si = Feather.write(sink_file, Feather.Source, source_file)
@test isequal(ds, Feather.read(sink_file))

source = Feather.Source(source_file)
Feather.write(sink_file, source)
@test isequal(ds, Feather.read(sink_file))

sink = Feather.Sink(joinpath(testdir, "test_utf8_new.feather"))
Feather.write(sink, Feather.Source, source_file)
@test isequal(ds, Feather.read(sink_file))

source = Feather.Source(source_file)
sink = Feather.Sink(joinpath(testdir, "test_utf8_new.feather"))
Feather.write(sink, source)
@test isequal(ds, Feather.read(sink_file))
rm(sink_file)

# check if valid, non-sudo docker is available
dockercheck = false
try
    dockercheck = success(`docker images`)
catch
    println("It seems that `docker` is not installed or has to be run with sudo, skipping python roundtrip tests")
end

if dockercheck

DOCKERTAG = "0.1"

println("Running python round-trip tests...")

println("Pulling feather docker image...")
run(`docker pull quinnj/feather:$DOCKERTAG`)

println("Create docker container from feather image...")
run(`docker run -it -d --name feathertest quinnj/feather:$DOCKERTAG /bin/sh`)

println("Generate a test.feather file from python...")
run(`docker cp runtests.py feathertest:/home/runtests.py`)
run(`docker exec feathertest python /home/runtests.py`)

println("Read test.feather into julia...")
run(`docker cp feathertest:/home/test.feather test.feather`)
df = Feather.read("test.feather")

@test df[:Autf8] == ["hey","there","sailor"]
@test df[:Abool] == [true, true, false]
@test df[:Acat] == CategoricalArrays.CategoricalArray(["a","b","c"])
@test df[:Acatordered] == CategoricalArrays.CategoricalArray(["d","e","f"])
@test df[:Adatetime] == [DateTime(2016,1,1), DateTime(2016,1,2), DateTime(2016,1,3)]
@test isequal(df[:Afloat32], NullableArray(Float32[1.0, 0.0, 0.0], [false, true, false]))
@test df[:Afloat64] == [Inf,1.0,0.0]

println("Writing test2.feather from julia...")
Feather.write("test2.feather", df)
df2 = Feather.read("test2.feather")

@test df2[:Autf8] == ["hey","there","sailor"]
@test df2[:Abool] == [true, true, false]
@test df2[:Acat] == CategoricalArrays.CategoricalArray(["a","b","c"])
@test df2[:Acatordered] == CategoricalArrays.CategoricalArray(["d","e","f"])
@test df2[:Adatetime] == [DateTime(2016,1,1), DateTime(2016,1,2), DateTime(2016,1,3)]
@test isequal(df2[:Afloat32], NullableArray(Float32[1.0, 0.0, 0.0], [false, true, false]))
@test df2[:Afloat64] == [Inf,1.0,0.0]

println("Read test2.feather into python...")
run(`docker cp test2.feather feathertest:/home/test2.feather`)
run(`docker cp runtests2.py feathertest:/home/runtests2.py`)
run(`docker exec feathertest python /home/runtests2.py`)

run(`docker stop feathertest`)
run(`docker rm feathertest`)
rm("test.feather")
rm("test2.feather")

end


text:
l=length(ARGS)

buildscript="BuildExecutable/src/build_executable.jl"

default_version=string(VERSION)

i=1
version=l<i?default_version:ARGS[i];i+=1
name=l<i?"app":ARGS[i];i+=1
tdir=l<i?"build":ARGS[i];i+=1
sdir=l<i?"source":ARGS[i];i+=1
script=l<i?"main.jl":ARGS[i];i+=1

debug = false
julia = replace(abspath(joinpath(JULIA_HOME, debug ? "julia-debug" : "julia")), default_version,version)

path = replace(dirname(Base.source_path()),"\\","/")

targetdir_build=joinpath(path,tdir)
targetdir_version=joinpath(targetdir_build,version)
targetdir_project=joinpath(targetdir_version,name)
targetdir = targetdir_project

sourcedir=replace(version,r"[.][0-9]+$","")
sourcescript=joinpath(path,"$sdir/$sourcedir/",script)
buildscriptpath=buildscript

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

if !isfile(julia * (is_windows() ? ".exe" : "")) error("'$julia' not found.") end

# find buildscript...
if !isfile(buildscriptpath)
  buildscriptpath=joinpath(path,buildscript)
  if !isfile(buildscriptpath)
    current = path
    prev = current
    found=false
    for i=1:10
      current=abspath(joinpath(current,"../"))
      buildscriptpath=joinpath(current,buildscript)
      if isfile(buildscriptpath) found=true; break end
      if current == prev break end #dublicate, sp break
      prev = current
    end
  end
  if !found error("'$buildscript' not found.") end
end

if !isfile(sourcescript) error("'$sourcescript' not found.") end

if !isdir(targetdir_build) mkdir(targetdir_build) end
if !isdir(targetdir_version) mkdir(targetdir_version) end
if !isdir(targetdir_project) mkdir(targetdir_project) end

if !isdir(targetdir_build) error("'$targetdir_build' failed to create.") end
if !isdir(targetdir_version) error("'$targetdir_version' failed to create.") end
if !isdir(targetdir_project) error("'$targetdir_project' failed to create.") end

cmd=`$julia $buildscriptpath --force "$name" $sourcescript $targetdir`

info(cmd)
run(cmd)

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
int_rules_4_7_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.7.6*(c+d*x)^m*trig(a+b*x)^n*trig(a+b*x)^p =#
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) ^ ~(n') * cos(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ ~m * sin(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) - ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * sin(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) * cos(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => (-((~c + ~d * ~x) ^ ~m) * cos(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) + ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * cos(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) ^ ~(n') * cos(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Antiderivative(ExpandTrigReduce((~c + ~d * ~x) ^ ~m, sin(~a + ~b * ~x) ^ ~n * cos(~a + ~b * ~x) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) ^ ~(n') * tan(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * sin(~a + ~b * ~x) ^ ~n * tan(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * sin(~a + ~b * ~x) ^ (~n - 2) * tan(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) ^ ~(n') * cot(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * cos(~a + ~b * ~x) ^ ~n * cot(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * cos(~a + ~b * ~x) ^ (~n - 2) * cot(~a + ~b * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sec(~(a') + ~(b') * ~x) ^ ~(n') * tan(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => ((~c + ~d * ~x) ^ ~m * sec(~a + ~b * ~x) ^ ~n) / (~b * ~n) - ((~d * ~m) / (~b * ~n)) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * sec(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~p, 1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ ~(n') * cot(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => (-((~c + ~d * ~x) ^ ~m) * csc(~a + ~b * ~x) ^ ~n) / (~b * ~n) + ((~d * ~m) / (~b * ~n)) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * csc(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (EqQ(~p, 1) && GtQ(~m, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sec(~(a') + ~(b') * ~x) ^ 2 * tan(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => ((~c + ~d * ~x) ^ ~m * tan(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) - ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * tan(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ 2 * cot(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => (-((~c + ~d * ~x) ^ ~m) * cot(~a + ~b * ~x) ^ (~n + 1)) / (~b * (~n + 1)) + ((~d * ~m) / (~b * (~n + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * cot(~a + ~b * ~x) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sec(~(a') + ~(b') * ~x) * tan(~(a') + ~(b') * ~x) ^ ~p, ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * sec(~a + ~b * ~x) * tan(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * sec(~a + ~b * ~x) ^ 3 * tan(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sec(~(a') + ~(b') * ~x) ^ ~(n') * tan(~(a') + ~(b') * ~x) ^ ~p, ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * sec(~a + ~b * ~x) ^ ~n * tan(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * sec(~a + ~b * ~x) ^ (~n + 2) * tan(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) * cot(~(a') + ~(b') * ~x) ^ ~p, ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * csc(~a + ~b * ~x) * cot(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * csc(~a + ~b * ~x) ^ 3 * cot(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ ~(n') * cot(~(a') + ~(b') * ~x) ^ ~p, ~x) => -(Antiderivative((~c + ~d * ~x) ^ ~m * csc(~a + ~b * ~x) ^ ~n * cot(~a + ~b * ~x) ^ (~p - 2), ~x)) + Antiderivative((~c + ~d * ~x) ^ ~m * csc(~a + ~b * ~x) ^ (~n + 2) * cot(~a + ~b * ~x) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && IGtQ(~p / 2, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * sec(~(a') + ~(b') * ~x) ^ ~(n') * tan(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Module([u = IntHide(sec(~a + ~b * ~x) ^ ~n * tan(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (IGtQ(~m, 0) && (IntegerQ(~n / 2) || IntegerQ((~p - 1) / 2)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ ~(n') * cot(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Module([u = IntHide(csc(~a + ~b * ~x) ^ ~n * cot(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (IGtQ(~m, 0) && (IntegerQ(~n / 2) || IntegerQ((~p - 1) / 2)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ ~(n') * sec(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => 2 ^ ~n * Antiderivative((~c + ~d * ~x) ^ ~m * csc(2 * ~a + 2 * ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IntegerQ(~n) && RationalQ(~m))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * csc(~(a') + ~(b') * ~x) ^ ~(n') * sec(~(a') + ~(b') * ~x) ^ ~(p'), ~x) => Module([u = IntHide(csc(~a + ~b * ~x) ^ ~n * sec(~a + ~b * ~x) ^ ~p, ~x)], Dist((~c + ~d * ~x) ^ ~m, u, ~x) - ~d * ~m * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * u, ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IntegersQ(~n, ~p) && (GtQ(~m, 0) && NeQ(~n, ~p)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~F)(~v) ^ ~(n') * (~G)(~w) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * (~F)(ExpandToSum(~v, ~x)) ^ ~n * (~G)(ExpandToSum(~v, ~x)) ^ ~p, ~x) <-- FreeQ([~m, ~n, ~p], ~x) && (TrigQ(~F) && (TrigQ(~G) && (EqQ(~v, ~w) && (LinearQ([~u, ~v, ~w], ~x) && Not(LinearMatchQ([~u, ~v, ~w], ~x))))))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(c') + ~(d') * ~x) * (~a + ~(b') * sin(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * sin(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * sin(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sin(~(c') + ~(d') * ~x) * (~a + ~(b') * cos(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * cos(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * cos(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sec(~(c') + ~(d') * ~x) ^ 2 * (~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * csc(~(c') + ~(d') * ~x) ^ 2 * (~a + ~(b') * cot(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * cot(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * cot(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sec(~(c') + ~(d') * ~x) * tan(~(c') + ~(d') * ~x) * (~a + ~(b') * sec(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => ((~e + ~f * ~x) ^ ~m * (~a + ~b * sec(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) - ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * sec(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * csc(~(c') + ~(d') * ~x) * cot(~(c') + ~(d') * ~x) * (~a + ~(b') * csc(~(c') + ~(d') * ~x)) ^ ~(n'), ~x) => (-((~e + ~f * ~x) ^ ~m) * (~a + ~b * csc(~c + ~d * ~x)) ^ (~n + 1)) / (~b * ~d * (~n + 1)) + ((~f * ~m) / (~b * ~d * (~n + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * (~a + ~b * csc(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (IGtQ(~m, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) ^ ~(p') * sin(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, sin(~a + ~b * ~x) ^ ~p * sin(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) ^ ~(p') * cos(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, cos(~a + ~b * ~x) ^ ~p * cos(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) ^ ~(p') * cos(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigReduce((~e + ~f * ~x) ^ ~m, sin(~a + ~b * ~x) ^ ~p * cos(~c + ~d * ~x) ^ ~q, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (IGtQ(~p, 0) && IGtQ(~q, 0))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * (~F)(~(a') + ~(b') * ~x) ^ ~(p') * (~G)(~(c') + ~(d') * ~x) ^ ~(q'), ~x) => Antiderivative(ExpandTrigExpand((~e + ~f * ~x) ^ ~m * (~G)(~c + ~d * ~x) ^ ~q, ~F, ~c + ~d * ~x, ~p, ~b / ~d, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (MemberQ([sin, cos], ~F) && (MemberQ([sec, csc], ~G) && (IGtQ(~p, 0) && (IGtQ(~q, 0) && (EqQ(~b * ~c - ~a * ~d, 0) && IGtQ(~b / ~d, 1))))))
end

max_stars_repo_path: src/rules/4 Trig functions/4.7 Miscellaneous/4.7.6 (c+d x)^m trig(a+b x)^n trig(a+b x)^p.jl


text:
function InverseLU(A::Array)

  #A = the matrix of coefficients 'A' must be squared

    m,n =size(A)

    if m != n
        error("Matrix must be squared")
    end

    In=zeros(m,m)
    B=Matrix{Float64}(I,m,m)
    D=zeros(m)

    L,U,P, p=DolittleLU(A)

    for j=1:m

        #Fordward Substitution
        D=zeros(m)
        D[1]=B[1,j]
        for i=2:m
           D[i]=(B[i,j] - (L[i,1:i-1]'*D[1:i-1]))
        end

        #Back substitution

          In[m,j]=D[end]/U[end,end]

    for i=m-1:-1:1
       In[i,j]=(D[i] - (U[i,i+1:m]'*In[i+1:m,j]))/U[i,i]
    end

    end
    return In, P
end


text:
using Firefly
using Test
using MCMCChains

@testset "Firefly.jl" begin
    @testset "Distributions" begin
        include("distributions.jl")
    end
    @testset "findpeak" begin
        X = ones(100)
        @test findpeak(X) ≈ 1 rtol = 1e-2
        @test findpeaks(X, 1) ≈ [1] rtol = 1e-2

        X = append!(ones(100), 2 .* ones(80))
        @test findpeak(X) ≈ 1 rtol = 1e-2
        @test findpeaks(X, 1) ≈ [1] rtol = 1e-2
        @test findpeaks(X, 2) ≈ [1, 2] rtol = 1e-2

        # interface
        c = Chains(reshape(X, (180, 1, 1)), ["x"])
        @test findpeak(c) == findpeak(X)
        @test findpeaks(c) == findpeaks(X)
        @test findpeak(c) == findpeak(c, :x)
        @test findpeaks(c) == findpeaks(c, :x)
        @test findpeaks(c, 1) == findpeaks(c, :x, 1)
        @test findpeaks(c, 1) == findpeaks(c, :x)[1:1]

        c2 = Chains(ones(100, 2, 1), ["x", "y"])
        @test_throws MethodError findpeak(c2)
        @test_throws MethodError findpeaks(c2)
        @test_throws MethodError findpeaks(ones(100, 100))
        @test_throws MethodError findpeaks(ones(100), 2.4)
        @test_throws MethodError findpeaks(c, 2.4)
    end
end


text:
__precompile__(true)
"""
StrBase package

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based partly on code in LegacyStrings that used to be part of Julia
"""
module StrBase

using ModuleInterfaceTools

@api extend! ChrBase

# Public Types

@api public Str, UniStr

# Convenience functions

@api public! unsafe_str, to_ascii, utf8, utf16, utf32, is_mutable, index

@api develop! check_string, unsafe_check_string, fast_check_string, skipascii, skipbmp,
              countmask, count_chars, _count_mask_al, _count_mask_ul, count_latin,
              _copysub, _cvtsize, _repeat, empty_str, _data, _pnt64, _str,
              ValidatedStyle, MutableStyle, EqualsStyle, CanContain

@api develop LineCounts, CharTypes, CharStat, maxbit, calcstats, check_continuation,
             UTF_LONG, UTF_LATIN1, UTF_UNICODE2, UTF_UNICODE3, UTF_UNICODE4, UTF_SURROGATE,
             UTF_INVALID, CHUNKSZ, CHUNKMSK,
             _memcmp, _memcpy, _memset, _fwd_memchr, _rev_memchr,
             empty_string, _calcpnt, _mask_bytes, _allocate,
             MS_UTF8, MS_UTF16, MS_UTF32, MS_SubUTF32, MS_Latin, MS_ByteStr, MS_RawUTF8,
             _wrap_substr, _empty_sub,
             AccessType, UInt16_U, UInt32_U, UInt16_S, UInt32_S, UInt16_US, UInt32_US,
             alignedtype, swappedtype, checkkeep, splitarr, __split, __rsplit, __replace

using Base: @_inline_meta, @propagate_inbounds, @_propagate_inbounds_meta, RefValue

# Convenience aliases
const SingleCU = SingleCodeUnitEncoding
const MultiCU  = MultiCodeUnitEncoding
const SCU = SingleCU()
const MCU = MultiCU()

include("types.jl")
include("fixparse.jl")
include("chars.jl")
include("access.jl")
include("traits.jl")
include("utf8proc.jl")
include("unicode.jl")
include("casefold.jl")
include("core.jl")
include("support.jl")
include("compare.jl")
include("ascii.jl")
include("latin.jl")
include("utf8.jl")
include("utf16.jl")
include("utf32.jl")
include("search.jl")
include("utf8search.jl")
include("utf16search.jl")
include("encode.jl")
include("stats.jl")
include("legacy.jl")
include("utf8case.jl")
include("utf16case.jl")
include("util.jl")
include("io.jl")
include("hash.jl")

@api freeze

end # module StrBase


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -48.987102771054815 <= q <= 45.974621434893855)

                     add_NL_constraint(m, :(log(1 + exp(-0.368688824745 + 0.553104331767496*log(1 + exp(-0.34679271547764046 + -0.03387429483407178*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.41639023650850415*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.7737560883943391*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.09930846098451918*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.9337943134926245*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.17925952214788987*log(1 + exp(0.6682158984675199 + 0.6233864552377666*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.35483759096347534*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4760912088664373*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + -0.23820709887370928*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9242931300611921*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.4698400936211189*log(1 + exp(0.8534714115378583 + 0.6703763846356168*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.08490767686166478*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.5656859981927576*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.4980009846607585*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.014109452266772848*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.7276938776175426*log(1 + exp(-0.09327714623589056 + -0.24359337663872616*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.2759577392076742*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4166615549845427*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.17410413013320136*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9439406543948436*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.4202042464647322*log(1 + exp(-0.09609361556421225 + -0.8712983668552394*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.8013603269787453*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.3524955117104245*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.28183131704908826*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.6019873452940772*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))))) + log(1 + exp(-0.8761821794177092 + 0.4562914248397121*log(1 + exp(-0.34679271547764046 + -0.03387429483407178*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.41639023650850415*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.7737560883943391*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.09930846098451918*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.9337943134926245*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.3827259206071312*log(1 + exp(0.6682158984675199 + 0.6233864552377666*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.35483759096347534*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4760912088664373*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + -0.23820709887370928*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9242931300611921*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.20062384412841983*log(1 + exp(0.8534714115378583 + 0.6703763846356168*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.08490767686166478*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.5656859981927576*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.4980009846607585*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.014109452266772848*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.28064116456943733*log(1 + exp(-0.09327714623589056 + -0.24359337663872616*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.2759577392076742*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4166615549845427*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.17410413013320136*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9439406543948436*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.4949171063694977*log(1 + exp(-0.09609361556421225 + -0.8712983668552394*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.8013603269787453*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.3524955117104245*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.28183131704908826*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.6019873452940772*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))))) + log(1 + exp(0.7716665307334409 + 0.20199447224764766*log(1 + exp(-0.34679271547764046 + -0.03387429483407178*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.41639023650850415*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.7737560883943391*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.09930846098451918*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.9337943134926245*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.3172295658101083*log(1 + exp(0.6682158984675199 + 0.6233864552377666*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.35483759096347534*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4760912088664373*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + -0.23820709887370928*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9242931300611921*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.10629001518557368*log(1 + exp(0.8534714115378583 + 0.6703763846356168*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.08490767686166478*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.5656859981927576*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.4980009846607585*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.014109452266772848*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.4930771030063288*log(1 + exp(-0.09327714623589056 + -0.24359337663872616*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.2759577392076742*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4166615549845427*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.17410413013320136*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9439406543948436*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.5476037221052898*log(1 + exp(-0.09609361556421225 + -0.8712983668552394*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.8013603269787453*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.3524955117104245*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.28183131704908826*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.6019873452940772*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))))) + log(1 + exp(0.6820199557891846 + 0.4693680154187745*log(1 + exp(-0.34679271547764046 + -0.03387429483407178*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.41639023650850415*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.7737560883943391*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.09930846098451918*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.9337943134926245*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.6552396777765495*log(1 + exp(0.6682158984675199 + 0.6233864552377666*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.35483759096347534*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4760912088664373*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + -0.23820709887370928*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9242931300611921*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.535415419107649*log(1 + exp(0.8534714115378583 + 0.6703763846356168*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.08490767686166478*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.5656859981927576*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.4980009846607585*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.014109452266772848*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.022291893207070146*log(1 + exp(-0.09327714623589056 + -0.24359337663872616*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.2759577392076742*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4166615549845427*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.17410413013320136*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9439406543948436*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.04141935716867051*log(1 + exp(-0.09609361556421225 + -0.8712983668552394*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.8013603269787453*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.3524955117104245*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.28183131704908826*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.6019873452940772*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))))) + log(1 + exp(-0.779255880286021 + 0.8816016628889711*log(1 + exp(-0.34679271547764046 + -0.03387429483407178*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.41639023650850415*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.7737560883943391*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.09930846098451918*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.9337943134926245*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.3409151924818685*log(1 + exp(0.6682158984675199 + 0.6233864552377666*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.35483759096347534*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4760912088664373*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + -0.23820709887370928*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9242931300611921*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.5165976261581497*log(1 + exp(0.8534714115378583 + 0.6703763846356168*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.08490767686166478*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.5656859981927576*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.4980009846607585*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.014109452266772848*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + -0.06878582640645536*log(1 + exp(-0.09327714623589056 + -0.24359337663872616*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + 0.2759577392076742*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + -0.4166615549845427*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.17410413013320136*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + 0.9439406543948436*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))) + 0.8469129987678201*log(1 + exp(-0.09609361556421225 + -0.8712983668552394*log(1 + exp(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3]))) + -0.8013603269787453*log(1 + exp(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3]))) + 0.3524955117104245*log(1 + exp(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3]))) + 0.28183131704908826*log(1 + exp(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3]))) + -0.6019873452940772*log(1 + exp(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function draw_synapse_panel(gui_data::GuiData, app_data::Model.AppData, sock::Comm.SocClient)
    if CImGui.TreeNode("Synapse")
        CImGui.PushItemWidth(80)

        gui_data.buffer = Model.prep_field(Model.alpha(app_data.model), 16)
        returned = CImGui.InputText("Alpha", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_alphas!(app_data.model, gui_data.buffer)
            else
                Model.set_alpha!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(180)

        gui_data.buffer = Model.prep_field(Model.ama(app_data.model), 16)
        returned = CImGui.InputText("ama", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_amas!(app_data.model, gui_data.buffer)
            else
                Model.set_ama!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(330)

        gui_data.buffer = Model.prep_field(Model.amb(app_data.model), 16)
        returned = CImGui.InputText("amb", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_ambs!(app_data.model, gui_data.buffer)
            else
                Model.set_amb!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(480)

        gui_data.buffer = Model.prep_field(Model.lambda(app_data.model), 16)
        returned = CImGui.InputText("lambda", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_lambdas!(app_data.model, gui_data.buffer)
            else
                Model.set_lambda!(app_data.model, gui_data.buffer)
            end
        end

        # new row
        gui_data.buffer = Model.prep_field(Model.learning_rate_fast(app_data.model), 16)
        returned = CImGui.InputText("Fast Learn Rate", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_learning_rate_fasts!(app_data.model, gui_data.buffer)
            else
                Model.set_learning_rate_fast!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(280)

        gui_data.buffer = Model.prep_field(Model.learning_rate_slow(app_data.model), 16)
        returned = CImGui.InputText("Slow Learn Rate", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_learning_rate_slows!(app_data.model, gui_data.buffer)
            else
                Model.set_learning_rate_slow!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(530)

        gui_data.buffer = Model.prep_field(Model.mu(app_data.model), 16)
        returned = CImGui.InputText("mu", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_mus!(app_data.model, gui_data.buffer)
            else
                Model.set_mu!(app_data.model, gui_data.buffer)
            end
        end

        # new row
        gui_data.buffer = Model.prep_field(Model.taoI(app_data.model), 16)
        returned = CImGui.InputText("taoI", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_taoIs!(app_data.model, gui_data.buffer)
            else
                Model.set_taoI!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(150)

        gui_data.buffer = Model.prep_field(Model.taoN(app_data.model), 16)
        returned = CImGui.InputText("taoN", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_taoNs!(app_data.model, gui_data.buffer)
            else
                Model.set_taoN!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.SameLine(300)

        gui_data.buffer = Model.prep_field(Model.taoP(app_data.model), 16)
        returned = CImGui.InputText("taoP", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_taoPs!(app_data.model, gui_data.buffer)
            else
                Model.set_taoP!(app_data.model, gui_data.buffer)
            end
        end

        # new row
        gui_data.buffer = Model.prep_field(Model.distance(app_data.model), 16)
        returned = CImGui.InputText("distance", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_distances!(app_data.model, gui_data.buffer)
            else
                Model.set_distance!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.PopItemWidth()
        CImGui.SameLine(200)

        CImGui.PushItemWidth(200)
        gui_data.buffer = Model.prep_field(Model.weight(app_data.model), 32)
        returned = CImGui.InputText("weight", gui_data.buffer, 32, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            if app_data.model.apply_to_all
                Model.set_weights!(app_data.model, gui_data.buffer)
            else
                Model.set_weight!(app_data.model, gui_data.buffer)
            end
        end
        CImGui.PopItemWidth()

        CImGui.PushItemWidth(500)
        active = Cint(Model.active_synapse(app_data.model))
        synapses = Cint(Model.synapses(app_data.model))

        @c CImGui.SliderInt("Synapse##7", &active, 1, synapses)
        Model.set_active_synapse!(app_data.model, Int64(active))
        CImGui.PopItemWidth()

        CImGui.SameLine(650)
        @c CImGui.Checkbox("Apply To All", &app_data.model.apply_to_all)

        CImGui.TreePop()
    end
end


text:
abstract type AbstractBandwidth end
struct Smoothed{B,BB} <: AbstractBandwidth
    base_bandwidth::BB
    γ::B
    α::Float64
end
function update!(s::Smoothed, x)
    b = update!(base_bandwidth, x)
    copyto!(b.γ, s.α * b.γ .+ (1 - b.α) * b)
end
struct Scheduled{B,S} <: AbstractBandwidth
    bandwidth::B
    schedule::S
end
bandwidth(b::Scheduled) = b.schedule(bandwidth(b.bandwidth))
function update!(s::Scheduled, x)
    update!(s.bandwidth, x)
    update!(s.schedule)
    s
end
mutable struct RosenblattSchedule
    n::Int
    exponent::Float64
end
RosenblattSchedule(d) = RosenblattSchedule(0, -1/(d + 4))
(s::RosenblattSchedule)(b) = b * s.n^s.exponent
update!(s::RosenblattSchedule) = s.n += 1
struct Bandwidth{H,B} <: AbstractBandwidth
    heuristic::H
    γ::B
    function Bandwidth(heuristic::H, gamma::T) where {H, T}
        if T <: Number && H <: AbstractHeuristic
            gamma = Ref(gamma)
        end
        new{H, typeof(gamma)}(heuristic, gamma)
    end
end
Bandwidth(; heuristic = MedianHeuristic(), gamma = NaN) = Bandwidth(heuristic, gamma)
Base.copyto!(x::Base.RefValue, y) = x[] = y
bandwidth(b::AbstractBandwidth) = b.γ
bandwidth(b::Bandwidth{<:Any, <:Base.RefValue}) = b.γ[]
abstract type AbstractHeuristic end
struct MedianHeuristic <: AbstractHeuristic
    factor::Float64
end
MedianHeuristic() = MedianHeuristic(1.)
struct ScottsHeuristic <: AbstractHeuristic end
update!(b::Bandwidth{MedianHeuristic}, x) = copyto!(b.γ, b.heuristic.factor * median(x))
update!(b::Bandwidth{ScottsHeuristic}, x) = copyto!(b.γ, std(x))
update!(b::AbstractBandwidth, ::Any) = b

struct Kernel{T,B,F}
    bandwidth::B
    distance::F
end
"""
    Kernel(; gamma = NaN,
             type = Gaussian,
             distance = euclidean,
             heuristic = MedianHeuristic(),
             bandwidth = Bandwidth(heuristic, gamma))

Creates a `Kernel`. See `subtypes(LikelihoodfreeInference.AbstractKernel)`,
`subtypes(LikelihoodfreeInference.AbstractBandwidth)`,
`subtypes(LikelihoodfreeInference.AbstractHeuristic)` for options.
The bandwidth parameter `gamma` can be a number or a vector.

# Example 1
```julia
k = Kernel()
k(2.) == exp(-2^2/2)
using LikelihoodfreeInference: update!, bandwidth  # some internals.
bandwidth(k) === NaN                               # default bandwidth.
update!(k, [2., 2.])                               # median heuristic bandwidth.
bandwidth(k) == 2.
a = [1, 2, 3]; b = [4, 5, 6]
k(a, b) == k(k.distance(a, b, bandwidth(k))) == exp(-sum(abs2, a - b)/(2*bandwidth(k)^2))
```
# Example 2
```julia
k = Kernel(gamma = zeros(2), heuristic = ScottsHeuristic())
data = [2*randn(2) for _ in 1:5]
update!(k, data)
bandwidth(k) == std(data)
a = [1, 2]; b = [4, 5]
k(a, b) ≈ exp(-sum(abs2, (a - b) ./ bandwidth(k))/2)
```
"""
function Kernel(; gamma = NaN,
                  type = Gaussian,
                  distance = euclidean,
                  heuristic = MedianHeuristic(),
                  bandwidth = Bandwidth(heuristic, gamma)
                 )
    Kernel{type, typeof(bandwidth), typeof(distance)}(bandwidth, distance)
end
bandwidth(k::Kernel) = bandwidth(k.bandwidth)
update!(k::Kernel, x) = update!(k.bandwidth, x)
function update!(k::Kernel, simulated::Vector{T}, data::T) where T
    if T <: AbstractVector{<:Number}
        d = [k.distance(x, data) for x in simulated]
    elseif T <: AbstractVector{<:AbstractVector}
        d = vcat([[k.distance(xi, d) for xi in x] for x in simulated, d in data]...)
    end
    update!(k, d)
end
update!(k::Kernel{<:Any,<:Bandwidth{Nothing}}, ::Vector{T}, ::T) where T = k
(k::Kernel)(a, b) = k(k.distance(a, b, bandwidth(k)))

abstract type AbstractKernel end
struct Gaussian <: AbstractKernel end
struct ModifiedGaussian <: AbstractKernel end
(k::Kernel{Gaussian})(x) = exp(-x^2/2)
(k::Kernel{ModifiedGaussian})(x) = exp(-1/2 * (x < 1 ? x^2 : x))


text:
using Weber
using Base.Test
include("find_timing.jl")

when_events,_,_ = find_timing() do
  @addtrials let test = true
    @addtrials if test
      addtrial(moment(() -> (record(:a); test = false)))
    end

    @addtrials if test
      addtrial(moment(() -> (record(:b))))
    end
  end
end

elseif_events,_,_ = find_timing() do
  @addtrials let test = true, test2 = false
    @addtrials if test
      addtrial(moment(() -> (record(:a); test = false)))
    else
      addtrial(moment(() -> (record(:bad_a); test = true)))
    end

    @addtrials if test
      addtrial(moment(() -> (record(:bad_b1); test = false)))
    elseif !test2
      addtrial(moment(() -> (record(:b); test2 = true)))
    else
      addtrial(moment(() -> (record(:bad_b2); test = false)))
    end

    @addtrials if !test2
      addtrial(moment(() -> record(:bad_c1)))
    elseif test
      addtrial(moment(() -> record(:bad_c2)))
    else
      addtrial(moment(() -> record(:c)))
    end
  end
end

@testset "Conditional Moments" begin
  @test when_events == [:a]
  @test elseif_events == [:a,:b,:c]
end


text:
## test objective function works in parallel (distributed and multithreaded)
# using Distributed, SharedArrays
# addprocs(4, exeflags = "--project")
# @everywhere include("parallel_setup.jl")

# vector of random parameters
xs = [rand(60) for _ in 1:100]

# dictionaries to hold results
res_ser = Dict{Int, Float64}()
res_thr = copy(res_ser)
res_dis = copy(res_ser)

# serial
[res_ser[i] = objective(x) for (i,x) in enumerate(xs)]

# multi-threaded - NOT WORKING
Threads.@threads for (i,x) in collect(enumerate(xs))
    res_thr[i] = objective(x)
end

# distributed - WORKING
@everywhere begin
    call_fun(dic, xs) = [dic[i] = objective(x) for (i, x) in enumerate(xs)]
    xs = $(copy(xs))
    res_dis = $(copy(res_dis))
    call_fun(res_dis, xs)
end

# alternative distributed - WORKING
res_dis2 = SharedArray{Float64}(100)
@sync @distributed for (i,x) in collect(enumerate(xs))
    res_dis2[i] = objective(x)
end

@show res_ser[2];
@show res_thr[2];
@everywhere @show res_dis[2]
@show res_dis2[2];

## test if better to run two simulations or one big one
function lorenz(u,p,t)
    σ, β, ρ = p
    x,y,z = u

    dx = σ * (y - x)
    dy = x * (ρ  - z) - y
    dz = x * y - β * z

    return SA[dx,dy,dz]
end

p = (10, 8/3, 28)
u₀ = SA[0.0, 1.0, 10.0]
tspan = (0,100)
prob = ODEProblem(lorenz, u₀, tspan, p)

simulate1(prob) = solve(prob, Tsit5(), callback=cb2, abstol=1e-7, reltol=1e-7)

function simulate2(prob)
    sol1 = solve(prob, Tsit5(), callback=cb1, abstol=1e-7, reltol=1e-7)

    u0new = sol1.u[end]
    tnew = (sol1.t[end], 100)

    newprob = remake(prob, u0=u0new, tspan=tnew)
    sol2 = solve(newprob, Tsit5(), callback = cb2, abstol=1e-7, reltol=1e-7)

    return sol1, sol2
end

condition1(u,t,integrator) = u[1]
affect1!(integrator) = terminate!(integrator)
cb1 = ContinuousCallback(condition1, affect1!, save_positions=(false,true))

condition2(u,t,integrator) = u[3] > 45.0
affect2!(int) = terminate!(int)
cb2 = DiscreteCallback(condition2, affect2!, save_positions=(false,true))

## not giving same answer


text:
# Tests automatic differentiation on s2.jl:
include("../src/transit_structure.jl")
include("../src/s2.jl")
if VERSION >= v"0.7"
  using Test
else
  using Base.Test
end

using ForwardDiff
using DiffResults

# Randomizer seed
using Random
Random.seed!(42)

function s2_grad_func(r::T,b::T) where {T <: Real}
  # Computes the derivative of s_n(r,b) with respect to r, b.
  # Create a vector for use with ForwardDiff
  x=[r,b]
  # Now, define a wrapper of s2 for use with ForwardDiff:
  function diff_s2(x::Array{T,1}) where {T <: Real}
  # x should be a two-element vector with values [r,b]
  rc,bc = x
  s2c= s2(rc,bc)
  return s2c
  end

  # Set up a type to store s_n and it's Jacobian with respect to x:
  out = DiffResults.GradientResult(x)
  # Compute the Jacobian (and value):
  out = ForwardDiff.gradient!(out,diff_s2,x)
  # Place the value in the s_2 vector:
  s_2 = DiffResults.value(out)
  # And, place the Jacobian in an array:
  s2_gradient= DiffResults.gradient(out)
return s_2,s2_gradient
end

function s2_grad_num(r::T,b::T) where {T <: Real}
# Now, carry out finite-difference:
dq = big(1e-18)

# Allocate an array for s_2:
# Make BigFloat versions of r & b:
r_big = big(r); b_big = big(b)
# Compute s_n to BigFloat precision:
s_2_big= s2(r_big,b_big)
# Now, compute finite differences:
s2_grad_big= zeros(BigFloat,2)
s2_plus= s2(r_big+dq,b_big)
s2_minus= s2(r_big-dq,b_big)
s2_grad_big[1] = (s2_plus-s2_minus)*.5/dq
s2_plus=s2(r_big,b_big+dq)
s2_minus=s2(r_big,b_big-dq)
s2_grad_big[2] = (s2_plus-s2_minus)*.5/dq
return convert(Float64,s_2_big),convert(Array{Float64,1},s2_grad_big)
end

# In the following lines I test the derivatives for the few special
# cases needed for computing s_2.  Note that in the paper I report
# \Lambda, while this gets multiplied by -\pi to produce s_2.

function runtest_s2(r,b,test_name)
@testset "$test_name" begin
s_2,s2_gradient= s2_grad_func(r,b)
s2_big,s2_grad_numeric = s2_grad_num(r,b)
s2_grad_ana = zeros(2)
diff1 = s2_grad_numeric-s2_gradient
s_2,Eofk,Em1mKdm=s2!(r,b,s2_grad_ana)
diff2 = s2_grad_ana-s2_gradient
println("b : ",b," r: ",r," flux diff: ",s_2-s2_big," grad diff(num-auto): ",diff1," diff(ana-auto): ",diff2 )
#println("gradient: ",s2_gradient," expected: ",-[2,-2./3.])
@test isapprox(s2_grad_numeric[1],s2_gradient[1],atol=1e-8)
@test isapprox(s2_grad_numeric[2],s2_gradient[2],atol=1e-8)
@test isapprox(s2_grad_ana[1],s2_gradient[1],atol=1e-8)
@test isapprox(s2_grad_ana[2],s2_gradient[2],atol=1e-8)
end
return
end

# Try r=b=1/2 special case:
runtest_s2(0.5,0.5,"r=0.5;b=0.5")

# Try b=0 special case:
runtest_s2(0.3,0.0,"r=0.3;b=0.0")

# Try r+b=1 special case:
runtest_s2(0.2,0.8,"r=0.2;b=0.8")

# Try r=b <1/2 special case:
runtest_s2(0.3,0.3,"r=0.3;b=0.3")

# Try r=b > 1/2 special case:
runtest_s2(3.0,3.0,"r=3.0;b=3.0")

# Now, try a random case with b+r < 1:
b=r=2.0
while b+r > 1
  global r = 2rand(); global b= 2rand()
end
runtest_s2(r,b,"r+b < 1")

# Now, try a random case with b+r > 1:
b=r=0.
while b+r < 1
  global r = 2rand(); global b= 2rand()
end
runtest_s2(r,b,"r+b > 1")

@testset "s2 grid" begin
# Now, try out hard cases - ingress/egress of small planets:
r0 = [0.01,0.01,1.0,10.0,100.0]
nb = 200
#l_max = 20
l_max = 10
n_max = l_max^2+2*l_max

using PyPlot
fig,axes = subplots(1,length(r0))
get_cmap("plasma")
epsilon = 1e-12; delta = 1e-3
for i=1:length(r0)
  global r=r0[i]
  if r < 1.0
    bgrid = [linearspace(1e-15,epsilon,nb); linearspace(epsilon,delta,nb); linearspace(delta,r-delta,nb);
     -logarithmspace(log10(delta),log10(epsilon),nb) .+ r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+ r;
     linearspace(r+delta,1-r-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+ (1-r); linearspace(1-r-epsilon,1-r+epsilon,nb);
     logarithmspace(log10(epsilon),log10(delta),nb) .+ (1-r); linearspace(1-r+delta,1+r-delta,nb); - logarithmspace(log10(delta),log10(epsilon),nb) .+ (1+r);
     linearspace(1+r-epsilon,1+r-1e-15,nb)]
  else
    bgrid = [linearspace(r-1+1e-15,r-1+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+ (r-1); linearspace(r-1+delta,r-delta,nb);
     -logarithmspace(log10(delta),log10(epsilon),nb) .+ r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+ r;
     linearspace(r+delta,r+1-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+ (r+1); linearspace(r+1-epsilon,r+1-1e-15,nb)]
  end
  nbgrid = length(bgrid)
  s2_jac_grid = zeros(nbgrid,2)
  s2_grid = zeros(nbgrid)
  s2_grid_big = zeros(nbgrid)
  s2_jac_grid_num = zeros(nbgrid,2)
  s2_grad_ana = zeros(2)
  for j=1:nbgrid
    s_2,Eofk,Em1mKdm= s2!(r,bgrid[j],s2_grad_ana)
    s2_grid[j]=s_2
    s2_jac_grid[j,:]=s2_grad_ana
    s2_big,s2_grad_numeric = s2_grad_num(r,bgrid[j])
    s2_grid_big[j]=s2_big
    s2_jac_grid_num[j,:]= s2_grad_numeric
    @test isapprox(s_2,s2_big,atol=1e-8)
    @test isapprox(s2_grad_ana[1],s2_grad_numeric[1],atol=1e-8)
    @test isapprox(s2_grad_ana[2],s2_grad_numeric[2],atol=1e-8)
  end
# Now, make plots:
  ax = axes[i]
  ax[:semilogy](bgrid,abs.(s2_grid[:,1]-s2_grid_big[:,1]) .+1e-18,lw=1,label="s 2 ")
  ax[:semilogy](bgrid,abs.(s2_jac_grid[:,1]-s2_jac_grid_num[:,1]) .+1e-18,lw=1,label="ds2/dr")
  ax[:semilogy](bgrid,abs.(s2_jac_grid[:,2]-s2_jac_grid_num[:,2]) .+1e-18,lw=1,label="ds2/db")
  ax[:legend](loc="upper right",fontsize=6)
  ax[:set_xlabel]("b values")
  ax[:set_ylabel]("Derivative Error")
end

savefig("test_s2_gradient.pdf", bbox_inches="tight")

end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Divergences
using Test

function testfun(𝒟, t₀, s)
	println("Testing "*string(𝒟))
	for (f, v) ∈ t₀
		str2 = "    "*string(f)
		print(str2)
		d = map(a -> f(𝒟, a), s)
		@test maximum(d[2:end] .- v[2:end]) <= 1e-04
		@test d[1] ≈ v[1]
		s₀ = f != Divergences.eval ? d[2:end] : sum(d[2:end])
		@test f(𝒟, s[2:end]) ≈ s₀
		printstyled(" "*repeat(".", 40-length(str2))*" [✓]"*"\n", color = :green)
	end
end

#=
Check that all Divergences satisfy the normalization
1.  γ(1) ≖ 0
2. γ'(1) ≖ 0
3.  γ(x) ⩾ 0
=#
#region

#=
Check that all Divergences satisfy the normalization
1.  γ(1) ≖ 0
2. γ'(1) ≖ 0
3.  γ(x) ⩾ 0
=#
seq = 0:.1:3
divs = (KullbackLeibler(), ReverseKullbackLeibler(), Hellinger(), [CressieRead(p) for p ∈ (-0.5, 0.5, 2. )]...)
for d ∈ divs
	str = "Testing normalization: "*string(d)
	println(str)
	str2 = "     γ(1) == 0"
	print(str2)
	@test Divergences.eval(d, [1.]) == 0
	printstyled(" "*repeat(".", 50-length(str2))*" [✓]"*"\n", color = :green)

	str2 = "    γ'(1) == 0"
	print(str2)
	@test Divergences.gradient(d, [1.]) == [0.0]
	printstyled(" "*repeat(".", 50-length(str2))*" [✓]"*"\n", color = :green)

	str2 = "     γ(x) ⩾  0"
	print(str2)
	@test map(u -> Divergences.eval(d, [u])[1], seq) > [0.0]
	printstyled(" "*repeat(".", 50-length(str2))*" [✓]"*"\n", color = :green)
end
#endregion

## ---- CressieRead ----
#=
Test Divergence.eval
=#
#region
seq = 0:.1:3
t₀ = Dict(
	2    => [0.3333333333, 0.2835, 0.234667, 0.187833, 0.144, 0.104167, 0.0693333, 0.0405, 0.0186667, 0.00483333, 0., 0.00516667, 0.0213333, 0.0495, 0.0906667, 0.145833, 0.216, 0.302167, 0.405333, 0.5265, 0.666667, 0.826833, 1.008, 1.21117, 1.43733, 1.6875, 1.96267, 2.26383, 2.592, 2.94817, 3.33333],
	0.5  => [0.6666666666, 0.50883, 0.385924, 0.285756, 0.203976, 0.138071, 0.086344, 0.0475494, 0.0207223, 0.00508662, 0., 0.00491964, 0.0193789, 0.0429707, 0.0753365, 0.116156, 0.165144, 0.222038, 0.286605, 0.358626, 0.437903, 0.524252, 0.617503, 0.717497, 0.824085, 0.937129, 1.0565, 1.18207, 1.31373, 1.45136, 1.59487],
	-0.5 => [2., 0.935089, 0.611146, 0.40911, 0.270178, 0.171573, 0.101613, 0.0533599, 0.0222912, 0.00526681, 0., 0.00476461, 0.0182195, 0.0392983, 0.0671362, 0.101021, 0.140356, 0.184638, 0.233437, 0.28638, 0.343146, 0.403449, 0.467041, 0.5337, 0.603227, 0.675445, 0.750194, 0.827329, 0.90672, 0.988245, 1.0718]
)
for (kv, val) ∈ t₀
	cr = CressieRead(kv)
	str = "Testing "*string(cr)
	print(str)
	d = map(a -> Divergences.eval(cr, a), seq)
	@test maximum(d[2:end] .- val[2:end]) <= 1e-05
	@test d[1] ≈ val[1]
	printstyled(" "*repeat(".", 40-length(str))*" [✓]"*"\n", color = :green)
	@test Divergences.eval(cr, seq[2:end]) ≈ sum(d[2:end])
end

#=
Test gradient
=#
t₀ = Dict(
	2    => [-0.5, -0.495, -0.48, -0.455, -0.42, -0.375, -0.32, -0.255, -0.18, -0.095, 0., 0.105, 0.22, 0.345, 0.48, 0.625, 0.78, 0.945, 1.12, 1.305, 1.5, 1.705, 1.92, 2.145, 2.38, 2.625, 2.88, 3.145, 3.42, 3.705, 4.],
	0.5  => [-2., -1.36754, -1.10557, -0.904555, -0.735089, -0.585786, -0.450807, -0.32668, -0.211146, -0.102633, 0., 0.0976177, 0.19089, 0.280351, 0.366432, 0.44949, 0.529822, 0.607681, 0.683282, 0.75681, 0.828427, 0.898275, 0.966479, 1.03315, 1.09839, 1.16228, 1.2249, 1.28634,  1.34664, 1.40588, 1.4641],
	-0.5 => [-Inf, -4.32456, -2.47214, -1.65148, -1.16228, -0.828427, -0.581989, -0.390457, -0.236068, -0.108185, 0., 0.0930748, 0.174258, 0.245884, 0.309691, 0.367007, 0.418861, 0.46607, 0.509288, 0.549047, 0.585786, 0.619869, 0.6516, 0.681239, 0.709006, 0.735089, 0.759653, 0.782839, 0.804771, 0.82556, 0.84529]
)
for (kv, val) ∈ t₀
	cr = CressieRead(kv)
	str = "Testing "*string(cr)
	print(str)
	d = map(a -> Divergences.gradient(cr, a), seq)
	@test maximum(d[2:end] .- val[2:end]) <= 1e-05
	@test d[1] ≈ val[1]
	printstyled(" "*repeat(".", 40-length(str))*" [✓]"*"\n", color = :green)
	@test Divergences.gradient(cr, seq[2:end]) ≈ d[2:end]
end

#=
Test Divergence.hessian
=#
seq = 0:0.1:2
t₀ = Dict(
	2    => [+Inf, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.],
	0.5  => [Inf,3.16228,2.23607,1.82574,1.58114,1.41421,1.29099,1.19523,1.11803,1.05409,1.,0.953463,0.912871,0.877058,0.845154,0.816497,0.790569,0.766965,0.745356,0.725476,0.707107]
)
for (kv, val) ∈ t₀
	cr = CressieRead(kv)
	str = "Testing "*string(cr)
	print(str)
	d = map(a -> Divergences.hessian(cr, a), seq)
	@test maximum(d[2:end] .- val[2:end]) <= 1e-05
	@test d[1] ≈ val[1]
	printstyled(" "*repeat(".", 40-length(str))*" [✓]"*"\n", color = :green)
	@test Divergences.hessian(cr, seq[2:end]) ≈ d[2:end]
end
#endregion

## ---- KullbackLeibler ----
#region
t₀ = Dict(
	Divergences.eval     => [Inf, 0.669741, 0.478112, 0.338808, 0.233484, 0.153426, .0935046, .0503275, .0214852, 0.00517554, 0., 0.0048412, 0.0187859, 0.0410735, 0.0710611, 0.108198, 							0.152006, 0.202068, 0.258016, 0.319522, 0.386294, 0.458068, 0.534606, 0.615691, 0.701125, 0.790727, 0.88433, 0.98178, 1.08293, 1.18766, 1.29584],
	Divergences.gradient => [-Inf, -2.30259, -1.60944, -1.20397, -0.916291, -0.693147, -0.510826, -0.356675, -0.223144, -0.105361, 0., 0.0953102, 0.182322, 0.262364, 0.336472, 0.405465, 0.470004, 0.530628, 0.587787, 0.641854, 0.693147, 0.741937, 0.788457, 0.832909, 0.875469, 0.916291, 0.955511, 0.993252, 1.02962, 1.06471, 1.09861],
	Divergences.hessian  => [+Inf, 10., 5., 3.33333, 2.5, 2., 1.66667, 1.42857, 1.25, 1.11111, 1., 0.909091, 0.833333, 0.769231, 0.714286, 0.666667, 0.625, 0.588235, 0.555556, 0.526316, 0.5, 0.47619, 0.454545, 0.434783, 0.416667, 0.4, 0.384615, 0.37037, 0.357143, 0.344828, 0.333333]
)
𝒟 = 𝒦ℒ()
testfun(𝒦ℒ(), t₀, 0:0.1:3)
#endregion

## ---- ReverseKullbackLeibler ----
#region
t₀ = Dict(
	Divergences.eval     => [Inf ,1.40259,0.809438,0.503973,0.316291,0.193147,0.110826,0.0566749,0.0231436,0.00536052,0.,0.00468982,0.0176784,0.0376357,0.0635278,0.0945349,0.129996,0.169372,0.212213,0.258146,0.306853,0.358063,0.411543,0.467091,0.524531,0.583709,0.644489,0.706748,0.770381,0.835289,0.901388],
	Divergences.gradient => [-Inf,-9.,-4.,-2.33333,-1.5,-1.,-0.666667,-0.428571,-0.25,-0.111111,0.,0.0909091,0.166667,0.230769,0.285714,0.333333,0.375,0.411765,0.444444,0.473684,0.5,0.52381,0.545455,0.565217,0.583333,0.6,0.615385,0.62963,0.642857,0.655172,0.666667],
	Divergences.hessian  => [Inf, 100., 25., 11.1111, 6.25, 4., 2.77778, 2.04082, 1.5625, 1.23457, 1., 0.826446, 0.694444, 0.591716, 0.510204, 0.444444, 0.390625, 0.346021, 0.308642, 0.277008, 0.25, 0.226757, 0.206612, 0.189036, 0.173611, 0.16, 0.147929, 0.137174, 0.127551, 0.118906, 0.111111]
)

testfun(ℬ𝓊𝓇ℊ(), t₀, 0:0.1:3)
#endregion

## ---- Hellinger ----
#region
t₀ = Dict(
	Divergences.eval     => [2., 0.935089, 0.611146, 0.40911, 0.270178, 0.171573, 0.101613,  0.0533599, 0.0222912, 0.00526681, 0., 0.00476461, 0.0182195, 0.0392983, 0.0671362, 0.101021, 0.140356, 0.184638, 0.233437,  0.28638, 0.343146],
	Divergences.gradient => [-Inf, -4.32456, -2.47214, -1.65148, -1.16228, -0.828427, -0.581989, -0.390457, -0.236068, -0.108185, 0., 0.0930748, 0.174258, 0.245884, 0.309691, 0.367007, 0.418861, 0.46607, 0.509288, 0.549047, 0.585786],
	Divergences.hessian  => [Inf, 31.6228, 11.1803, 6.08581, 3.95285, 2.82843, 2.15166, 1.70747, 1.39754, 1.17121, 1., 0.866784, 0.760726, 0.67466, 0.603682, 0.544331, 0.494106, 0.451156, 0.414087, 0.38183, 0.353553]
)
testfun(ℋ𝒟(), t₀, 0:0.1:2)
#endregion

## ---- Chi Squared ----
#region
seq = 0:0.1:2
t₀ = Dict(
	Divergences.eval     => (seq .- 1).^2/2,
	Divergences.gradient => (seq .- 1),
	Divergences.hessian  => [1, seq[2:end]./seq[2:end]...]
)
testfun(χ²(), t₀, 0:0.1:2)
#endregion

## ---- Modified Divergence ----
#region
#=
Given a divergence γ(x), the modified divergence is
γᵤ(x) if x > ρ
γ(x) if x <= ρ
where ρ > 1
=#

t₀ = Dict(
	Divergences.eval => [Inf, 0.669741, 0.478112, 0.338808, 0.233484, 0.153426, 0.0935046, 0.0503275, 0.0214852, 0.00517554, 0., 0.0048412, 0.0187859, 0.0411847, 0.0719168, 0.110982, 0.158381, 0.214113, 0.278179, 0.350578, 0.43131],
	Divergences.gradient => [-Inf,-2.30259,-1.60944,-1.20397,-0.916291,-0.693147,-0.510826,-0.356675,-0.223144,-0.105361,0.,0.0953102,0.18232155679395456,0.265655,0.348988,0.432322,0.515655,0.598988,0.682322,0.765655,0.848988],
	Divergences.hessian => [Inf, 10., 5., 3.33333, 2.5, 2., 1.66667, 1.42857, 1.25, 1.11111, 1., 0.909091, 0.833333, 0.833333, 0.833333, 0.833333, 0.833333, 0.833333, 0.833333, 0.833333, 0.833333]
)
testfun(ModifiedDivergence(𝒦ℒ(), 1.2), t₀, 0:0.1:2)
#endregion

## ---- FullyModified Divergence ----
#region
#=
Given a divergence γ(x), the modified divergence is
γᵤ(x) if x >= ρ
γₗ(x) if x <= φ
γ(x) if x ∈ (φ, ρ)
where ρ > 1 && φ <1
=#
t₀ = Dict(
	Divergences.eval => [1.7039728053,1.14786,0.726195,0.438663,0.260981,0.147837,0.0755154,0.0313648,0.00783337,0.0000503359,0.00468982,0.0193798,0.0426784,0.0743798,0.114484,0.162991,0.219901,0.285213,0.358928,0.441046,0.531567,0.630491,0.737817,0.853546,0.977678,1.11021,1.25115,1.40049],
	Divergences.gradient => [-5.666666666,-4.44444,-3.22222,-2.0303,-1.27273,-0.818182,-0.515152,-0.298701,-0.136364,-0.010101,0.0909091,0.173611,0.25,0.326389,0.402778,0.479167,0.555556,0.631944,0.708333,0.784722,0.861111,0.9375,1.01389,1.09028,1.16667,1.24306,1.31944,1.3958],
	Divergences.hessian => [11.11111111,11.11111111,11.11111111,9.18274,5.16529,3.30579,2.29568,1.68663,1.29132,1.0203,0.826446,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444]
)
𝒟 = FullyModifiedDivergence(ℬ𝓊𝓇ℊ(), 0.3, 1.2)
testfun(𝒟, t₀, 0:0.11:3)
#endregion

## Additional tests
@test_throws(DimensionMismatch, Divergences.eval(𝒦ℒ(), rand(10), rand(11)))
@test_throws(DimensionMismatch, Divergences.eval(ℬ𝓊𝓇ℊ(), rand(10), rand(11)))
@test_throws(DimensionMismatch, Divergences.eval(𝒞ℛ(1), rand(10), rand(11)))
@test_throws(DimensionMismatch, Divergences.eval(χ²(), rand(10), rand(11)))

divs = (𝒦ℒ(), ℬ𝓊𝓇ℊ(), 𝒞ℛ(1), ℋ𝒟(), χ²())

for d in divs
	@show d
	r = rand(10)
	@test Divergences.eval(d, r) == Divergences.eval(d, r, ones(length(r)))
	@test Divergences.gradient(d, r) == Divergences.gradient(d, r, ones(length(r)))
	@test Divergences.hessian(d, r) == Divergences.hessian(d, r, ones(length(r)))
end

ℳ𝒟 = ModifiedDivergence(𝒦ℒ(), 1.2)
ℱℳ𝒟 = FullyModifiedDivergence(𝒦ℒ(), 0.9, 1.2)

s = 1:0.1:2
t₀ = [0.0176784,0.00428749,0.,0.00405552,0.015811,0.0347323,0.0606034,0.0934189,0.133179,0.179883,0.233532]
d₀ = map(u -> Divergences.eval(ℳ𝒟, u, 1.2), s)
@test maximum(d₀ .- t₀) < 1e-06
@test Divergences.eval(ℳ𝒟, s, repeat([1.2], length(s))) ≈ sum(d₀)
t₀ = [-0.182322, -0.0870114, 0., 0.0800427, 0.154151, 0.223988, 0.293433, 0.362877, 0.432322, 0.501766, 0.5712]
d₀ = map(u -> Divergences.gradient(ℳ𝒟, u, 1.2), s)
@test maximum(d₀ - t₀) < 2e-05
@test Divergences.gradient(ℳ𝒟, s, repeat([1.2], length(s))) ≈ d₀
t₀ =[1.,0.909091,0.833333,0.769231,0.714286,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444]
d₀ = map(u -> Divergences.hessian(ℳ𝒟, u, 1.2), s)
@test maximum(d₀ - t₀) < 2e-05
@test Divergences.hessian(ℳ𝒟, s, repeat([1.2], length(s))) ≈ d₀

s = 0.1:0.1:2
t₀ = [0.554094,0.457446,0.370059,0.29193,0.22306,0.16345,0.113099,0.0720079,0.0401755,0.0176024,0.00428749,
	  0.,0.00405552,0.015811,0.0347323,0.0606034,0.0934189,0.133179,0.179883,0.233532]
d₀ = map(u -> Divergences.eval(ℱℳ𝒟, u, 1.2), s)
@test maximum(d₀ .- t₀) < 1e-06
@test Divergences.eval(ℱℳ𝒟, s, repeat([1.2], length(s))) ≈ sum(d₀)
t₀ = [-1.01277,-0.920175,-0.827583,-0.73499,-0.642398,-0.549805,-0.457212,-0.36462,-0.272027,-0.179435,-0.0870114,
	  2.22045*10^-16,0.0800427,0.154151,0.223988,0.293433,0.362877,0.432322,0.501766,0.57121]
d₀ = map(u -> Divergences.gradient(ℱℳ𝒟, u, 1.2), s)
@test maximum(d₀ - t₀) < 2e-05
@test Divergences.gradient(ℱℳ𝒟, s, repeat([1.2], length(s))) ≈ d₀
t₀ = [0.925926,0.925926,0.925926,0.925926,0.925926,0.925926,0.925926,0.925926,0.925926,0.925926,0.909091,0.833333,
	  0.769231,0.714286,0.694444,0.694444,0.694444,0.694444,0.694444,0.694444]
d₀ = map(u -> Divergences.hessian(ℱℳ𝒟, u, 1.2), s)
@test maximum(d₀ - t₀) < 2e-05
@test Divergences.hessian(ℱℳ𝒟, s, repeat([1.2], length(s))) ≈ d₀

x = rand(10)
@test sum(Divergences.gradient(ℱℳ𝒟, x)) ≈ Divergences.gradient_sum(ℱℳ𝒟, x)
@test sum(Divergences.hessian(ℱℳ𝒟, x)) ≈ Divergences.hessian_sum(ℱℳ𝒟, x)

@test Divergences.eval(ℱℳ𝒟, 3.2) ≈ Divergences.eval(ℱℳ𝒟, 3.2, 1.0)
@test Divergences.gradient(ℱℳ𝒟, 3.2) ≈ Divergences.gradient(ℱℳ𝒟, 3.2, 1.0)
@test Divergences.hessian(ℱℳ𝒟, 3.2) ≈ Divergences.hessian(ℱℳ𝒟, 3.2, 1.0)

# Divergences.hessian(ℱℳ𝒟, 3.2, 1)

# using ForwardDiff

# f(x) = Divergences.eval(ℱℳ𝒟, x)
# ForwardDiff.gradient(f, rand(10))

# f(x) = Divergences.eval(ℱℳ𝒟, x, rand(10))
# ForwardDiff.gradient(f, rand(10))

# ξ = rand(1_000_000);
# using BenchmarkTools
# @btime Divergences.eval(ℱℳ𝒟, ξ)

# @btime Divergences.eval(𝒦ℒ(), ξ)
# @btime Divergences.gradient(𝒦ℒ(), ξ);
# @btime Divergences.hessian(𝒦ℒ(), ξ);


text:
abstract type TreeNode end

# Leaves

mutable struct Leaf <: TreeNode
    pred :: Int
end

Leaf() = Leaf(1)

loss!(l :: Leaf, x :: Array{Float32, 2}, y :: Vector{Int}) = 0

function precompute!(l :: Leaf, x :: Array{Float32, 2}, y :: Vector{Int})
    if size(x, 2) != 0
        l.pred = mode(y)
    end
end

function predict(l :: Leaf, x :: Array{Float32, 2})
    fill(l.pred, size(x, 2))
end

# Nodes

mutable struct Node <: TreeNode
    splitter :: Flux.Chain
    subset :: Vector{Int}
    best :: Pair{Int, Int}
    depth :: Int
    impurity :: Float32
    left_split :: BitVector
    right_split :: BitVector
    left :: TreeNode
    right :: TreeNode
end

Flux.@functor Node

function Node(feat_sub :: Dict{Int, Vector{Int}}, hidden :: Int, depth :: Int, id :: Int)
    splitter = Chain(
        Dense(length(feat_sub[id]), hidden, tanh),
        Dense(hidden, hidden, tanh),
        Dense(hidden, 2),
        softmax
    )

    subset = feat_sub[id]

    if depth != 1
        id += 1
        left = Node(feat_sub, hidden, depth - 1, id)
        id += 1
        right = Node(feat_sub, hidden, depth - 1, id)
    else
        left = Leaf()
        right = Leaf()
    end

    Node(splitter, subset, Pair(1, 1), depth, Float32(0.0), BitVector(), BitVector(), left, right)
end

Flux.trainable(l :: Leaf) = ()
Flux.trainable(n :: Node) = (n.splitter, Flux.trainable(n.left)..., Flux.trainable(n.right)...,)

function precompute!(n :: Node, x :: Array{Float32, 2}, y :: Vector{Int})
    if size(x, 2) != 0
        decision = n.splitter(x[n.subset, :])
        n.left_split = decision[1, :] .>= 0.5
        n.right_split = decision[1, :] .< 0.5

        if length(y[n.left_split]) == 0
            left_best = 1
        else
            left_best = mode(y[n.left_split])
        end

        if length(y[n.right_split]) == 0
            right_best = 1
        else
            right_best = mode(y[n.right_split])
        end

        n.best = Pair(left_best, right_best)
        precompute!(n.left, x[:, n.left_split], y[n.left_split])
        precompute!(n.right, x[:, n.right_split], y[n.right_split])
    end
end

function predict(n :: Node, x :: Array{Float32, 2})
    decision = n.splitter(x[n.subset, :])
    left_split = decision[1, :] .>= 0.5
    right_split = decision[1, :] .< 0.5
    left_data = x[:, left_split]
    right_data = x[:, right_split]

    preds = Int.(zeros(size(x, 2)))
    preds[left_split] = predict(n.left, left_data)
    preds[right_split] = predict(n.right, right_data)

    return preds
end

function loss!(n :: Node, x :: Array{Float32, 2}, y :: Vector{Int})
    loss = 0
    if size(x, 2) != 0
        decision = n.splitter(x[n.subset, :])
        left = decision[1, :]
        right = decision[2, :]

        classes = 1:maximum(y)
        ȳ = Flux.onehotbatch(y, classes)

        left_weight = left' .* ȳ
        right_weight = right' .* ȳ

        left_best = Flux.onehotbatch(fill(n.best.first, length(y)), classes)
        right_best = Flux.onehotbatch(fill(n.best.second, length(y)), classes)

        n.impurity = Flux.crossentropy(left_weight, left_best)
        n.impurity += Flux.crossentropy(right_weight, right_best)

        loss += n.impurity
        loss += loss!(n.left, x, y)
        loss += loss!(n.right, x, y)

        return loss
    else
        return loss
    end
end

function splitter_importance(n :: Node, x :: Array{Float32, 2})
    val = min(200, size(x, 2))
    function predict_wrapper(n :: Node, data :: DataFrame)
        pred = DataFrame(y_pred = predict(n, Float32.(Array(data)')))
    end

    data_shap = ShapML.shap(explain = DataFrame(x[:, 1:val]'),
        reference = DataFrame(x[:, end-val+1:end]'),
        model = n,
        predict_function = predict_wrapper,
        sample_size = 150,
        parallel = :both,
    )

    s̄ = DataFrames.by(data_shap, [:feature_name],
        mean_effect = [:shap_effect] => x -> mean(abs.(x.shap_effect)))

    imp = Dict{Int, Float32}()
    for i ∈ 1:length(n.subset)
        @inbounds imp[n.subset[i]] = s̄.mean_effect[i]
    end
    return imp
end

function importance(n :: Node, x :: Array{Float32, 2})
    try
        scores = splitter_importance(n, x)
        for (k, v) ∈ scores
            scores[k] = v * 1/(n.impurity + 0.000001)
        end

        if n.depth > 1
            left_scores = importance(n.left, x[:, n.left_split])
            right_scores = importance(n.right, x[:, n.right_split])

            for d ∈ [left_scores, right_scores]
                for (k, v) ∈ d
                    if k ∈ keys(scores)
                        scores[k] += v
                    else
                        scores[k] = v
                    end
                end
            end
        end
        return scores
    catch e
        return Dict()
    end
end

function plot_importance(imp :: Dict{Int, Float32}, x :: Array{Float32, 2})
    bar(string.(1:size(x, 1)), [imp[x] for x ∈ 1:size(x, 1)], legend=:false, color=1:size(x, 1),
        title="Relative Feature Importance", xlabel="Feature", ylabel="Importance")
end

function tree_train!(epochs :: Int, n :: Node, x :: Array{Float32, 2}, y :: Vector{Int})
    ps = params(n)
    opt = NADAM(0.00025)

    pbar = Progress(epochs)
    for epoch ∈ 1:epochs
        local training_loss

        DeepForest.precompute!(n, x, y)
        gs = gradient(ps) do
            training_loss = loss!(n, x, y)
            return training_loss
        end
        next!(pbar; showvalues = [(:loss, training_loss), (:acc, mean(predict(n, x) .== y))])
        Flux.update!(opt, ps, gs)
    end
end


text:
##### Get prognostic variable list

import ..BalanceLaws: prognostic_vars

prognostic_vars(::DryModel) = ()
prognostic_vars(::EquilMoist) = (TotalMoisture(),)
prognostic_vars(::NonEquilMoist) =
    (TotalMoisture(), LiquidMoisture(), IceMoisture())

prognostic_vars(::NoPrecipitation) = ()
prognostic_vars(::RainModel) = (Rain(),)
prognostic_vars(::RainSnowModel) = (Rain(), Snow())

prognostic_vars(::NoTracers) = ()
prognostic_vars(::NTracers{N}) where {N} = (Tracers{N}(),)

prognostic_vars(m::AtmosModel) = (
    Mass(),
    Momentum(),
    Energy(),
    prognostic_vars(m.moisture)...,
    prognostic_vars(m.precipitation)...,
    prognostic_vars(m.tracers)...,
    prognostic_vars(m.turbconv)...,
)

max_forks_repo_name: jm-c/CLIMA


text:
import POMDPs.initialstate_distribution
const IVec8 = SVector{8, Int}

@with_kw struct AODiscreteVDPTagPOMDP <: POMDP{TagState, Int, IVec8}
    cpomdp::VDPTagPOMDP = VDPTagPOMDP()
    n_angles::Int       = 10
    binsize::Float64    = 0.5
end

@with_kw struct ADiscreteVDPTagPOMDP <: POMDP{TagState, Int, Vec8}
    cpomdp::VDPTagPOMDP = VDPTagPOMDP()
    n_angles::Int       = 10
end

const DiscreteVDPTagProblem = Union{AODiscreteVDPTagPOMDP, ADiscreteVDPTagPOMDP}

cproblem(p::AODiscreteVDPTagPOMDP) = p.cpomdp
cproblem(p::ADiscreteVDPTagPOMDP) = p.cpomdp

convert_s(::Type{T}, x::T, p) where T = x
convert_a(::Type{T}, x::T, p) where T = x
convert_o(::Type{T}, x::T, p) where T = x

# state
function convert_s(::Type{Int}, s::TagState, p::DiscreteVDPTagProblem)
    n = p.n_bins
    factor = n/(2*p.grid_lim)
    ai = clamp(ceil(Int, (s.agent[1]+p.grid_lim)*factor), 1, n)
    aj = clamp(ceil(Int, (s.agent[2]+p.grid_lim)*factor), 1, n)
    ti = clamp(ceil(Int, (s.target[1]+p.grid_lim)*factor), 1, n)
    tj = clamp(ceil(Int, (s.target[2]+p.grid_lim)*factor), 1, n)
    return sub2ind((n,n,n,n), ai, aj, ti, tj)
end
function convert_s(::Type{TagState}, s::Int, p::DiscreteVDPTagProblem)
    n = p.n_bins
    factor = 2*p.grid_lim/n
    ai, aj, ti, tj = ind2sub((n,n,n,n), s)
    return TagState((Vec2(ai, aj)-0.5)*factor-p.grid_lim, (Vec2(ti, tj)-0.5)*factor-p.grid_lim)
end

# action
function convert_a(::Type{Int}, a::Float64, p::DiscreteVDPTagProblem)
    i = ceil(Int, a*p.n_angles/(2*pi))
    while i > p.n_angles
        i -= p.n_angles
    end
    while i < 1
        i += p.n_angles
    end
    return i
end
convert_a(::Type{Float64}, a::Int, p::DiscreteVDPTagProblem) = (a-0.5)*2*pi/p.n_angles

function convert_a(T::Type{Int}, a::TagAction, p::DiscreteVDPTagProblem)
    i = convert_a(T, a.angle, p)
    if a.look
        return i + p.n_angles
    else
        return i
    end
end
function convert_a(::Type{TagAction}, a::Int, p::DiscreteVDPTagProblem)
    return TagAction(a > p.n_angles, convert_a(Float64, a % p.n_angles, p))
end

# observation
function convert_o(::Type{IVec8}, o::Vec8, p::AODiscreteVDPTagPOMDP)
    return floor.(Int, (o./p.binsize)::Vec8)::IVec8
end
# convert_o(::Type{Vec8}, o::Int, p::DiscreteVDPTagProblem) = (o-0.5)*2*pi/p.n_obs_angles

n_states(p::AODiscreteVDPTagPOMDP) = Inf
n_actions(p::DiscreteVDPTagProblem) = 2*p.n_angles
POMDPs.discount(p::DiscreteVDPTagProblem) = discount(cproblem(p))
isterminal(p::DiscreteVDPTagProblem, s) = isterminal(cproblem(p), convert_s(TagState, s, p))

POMDPs.actions(p::DiscreteVDPTagProblem) = 1:n_actions(p)

function POMDPs.gen(p::DiscreteVDPTagProblem, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return gen(cproblem(p), s, ca, rng)
end

function POMDPs.gen(::DDNOut{(:sp,:o,:r)}, p::ADiscreteVDPTagPOMDP, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return gen(DDNOut(:sp,:o,:r), cproblem(p), s, ca, rng)
end

function POMDPs.gen(n::DDNNode{:o}, p::ADiscreteVDPTagPOMDP, s::TagState, a::Int, sp::TagState, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    return POMDPs.gen(n, cproblem(p), s, ca, sp, rng)
end

function POMDPs.gen(::DDNOut{(:sp,:o,:r)}, p::AODiscreteVDPTagPOMDP, s::TagState, a::Int, rng::AbstractRNG)
    ca = convert_a(actiontype(cproblem(p)), a, p)
    csor = gen(DDNOut(:sp,:o,:r), cproblem(p), s, ca, rng)
    return (csor[1], convert_o(IVec8, csor[2], p), csor[3])
end

function POMDPs.gen(n::DDNNode{:o}, p::AODiscreteVDPTagPOMDP, s::TagState, a::Int, sp::TagState, rng::AbstractRNG)
    ca = convert_a(action_type(cproblem(p)), a, p)
    co = POMDPs.gen(n, cproblem(p), s, ca, sp, rng)
    return convert_o(IVec8, co, p)
end

POMDPs.initialstate_distribution(p::AODiscreteVDPTagPOMDP) = VDPInitDist()
POMDPs.initialstate_distribution(p::ADiscreteVDPTagPOMDP) = VDPInitDist()

#=
gauss_cdf(mean, std, x) = 0.5*(1.0+erf((x-mean)/(std*sqrt(2))))
function obs_weight(p::AODiscreteVDPTagPOMDP, a::Int, sp::TagState, o::Int)
    cp = cproblem(p)
    @assert cp.bearing_std <= 2*pi/6.0 "obs_weight assumes σ <= $(2*pi/6.0)"
    ca = convert_a(action_type(cp), a, p)
    co = convert_o(obs_type(cp), o, p) # float between 0 and 2pi
    upper = co + 0.5*2*pi/p.n_angles
    lower = co - 0.5*2*pi/p.n_angles
    if ca.look
        diff = sp.target - sp.agent
        bearing = atan(diff[2], diff[1])
        # three cases: o is in bin, below, or above
        if bearing <= upper && bearing > lower
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            prob = cdf_up - cdf_low
        elseif bearing <= lower
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            below_cdf_up = gauss_cdf(bearing, cp.bearing_std, upper-2*pi)
            below_cdf_low = gauss_cdf(bearing, cp.bearing_std, lower-2*pi)
            prob = cdf_up - cdf_low + below_cdf_up - below_cdf_low
        else # bearing > upper
            cdf_up = gauss_cdf(bearing, cp.bearing_std, upper)
            cdf_low = gauss_cdf(bearing, cp.bearing_std, lower)
            above_cdf_up = gauss_cdf(bearing, cp.bearing_std, upper+2*pi)
            above_cdf_low = gauss_cdf(bearing, cp.bearing_std, lower+2*pi)
            prob = cdf_up - cdf_low + above_cdf_up - above_cdf_low
        end
        return prob
    else
        return 1.0
    end
end

function obs_weight(p::ADiscreteVDPTagPOMDP, a::Int, sp::TagState, o::Float64)
    ca = convert_a(TagAction, a, p)
    return obs_weight(cproblem(p), ca, sp, o)
end
=#


text:
const S_N_NG86_LOOKUP = SingleCodonLookup{Tuple{Float64, Float64}}
const DS_DN_NG86_LOOKUP = PairwiseCodonLookup{Tuple{Float64, Float64}}

include("computation.jl")
include("lookups.jl")

"""
    dNdS_NG86(x, y, addone::Bool = true, code::Int = 1)

Compute dN and dS, using the [Nei and Gojobori 1986](https://www.ncbi.nlm.nih.gov/pubmed/3444411) method.

The genetic code that is used, is defined according to the numbering of
`ncbi_trans_table`. Code 1 is the standard genetic code.

This function requires two iterables `x` and `y`.
If these iterables yield `Codon{DNA}` or `Codon{RNA}` type variables. Then it is
assumed that `x` and `y` are iterables that yield a sequence of aligned codons.
If the iterables produce `DNA` or `RNA` type variables, then it is assumed `x`
and `y` iterables that conform to the behaviour of DNA or RNA sequences as
defined in the BioSequences package. In this case, a new `x` and `y` that do
have an element type of `Codon{DNA}` or `Codon{RNA}`.

NG86 is a counting method of computing dN/dS and is typically safer to use on
sequence data where codon usage, (esp. at 3rd position), is uniform,
the sequences are not very divergent, and transition/transversion rates, are similar.
"""
function dNdS_NG86(x, y; addone::Bool = false, code::Int = 1)
    snlookup = S_N_NG86_LOOKUPS[code]
    dsdnlookup = DS_DN_NG86_LOOKUPS[code]
    return _dNdS_NG86(x, y, addone, snlookup, dsdnlookup, eltype(x), eltype(y))
end

function pairwise_dNdS_NG86(x, addone::Bool = false, code::Int = 1)
    n = length(x)
    @assert n >= 2 "At least two sequences are required."
    snlookup = S_N_NG86_LOOKUPS[code]
    dsdnlookup = DS_DN_NG86_LOOKUPS[code]
    results = Matrix{Tuple{Float64, Float64}}(n, n)
    for i in 1:n
        results[i,i] = 0.0, 0.0
        for j in (i + 1):n
            results[i,j] = results[j,i] = _dNdS_NG86(x[i], x[j], addone, snlookup, dsdnlookup, eltype(x[i]), eltype(x[j]))
        end
    end
    return results
end


text:
"""
Code for convenience functions for calculating CCFs on RvSpectMLBase's types
At some point could consider moving this to RvSpectML.

Author: Eric Ford
Created: August 2020
"""

using Statistics

using NamedTupleTools
include("util.jl")

include("ccf_chunk.jl")

include("ccf_chunklist.jl")

using ThreadedIterables
include("ccf_chunklist_timeseries.jl")

include("order_ccf_chunklist.jl")

#using ThreadTools
include("order_ccf_chunklist_timeseries.jl")

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
struct DeterministicSamples{T<:Real} <: ContinuousMultivariateSampleable
    samples::Array{T,2}
end

Base.length(s::DeterministicSamples) = size(s.samples, 1)

function Distributions.rand(
    rng::AbstractRNG, s::DeterministicSamples{T}, num_requested_samples::Int
) where {T<:Real}
    samples = s.samples
    x = Array{T, 2}(undef, length(s), num_requested_samples)
    num_samples = size(samples, 2)
    num_first_samples = min(num_requested_samples, num_samples)
    x[:, 1:num_first_samples] = samples[:, end-num_first_samples+1:end]
    num_last_samples = num_requested_samples - num_first_samples
    if num_last_samples > 0
        x[:, num_first_samples+1:end] = samples[:, rand(rng, 1:num_samples, num_last_samples)]
    end
    return x
end

function Distributions.rand(
    ::AbstractRNG, s::DeterministicSamples{T}, ::Nothing
) where {T<:Real}
    return s.samples
end

mean(s::DeterministicSamples) = Statistics.mean(s.samples, dims=2)[:,1]
cov(s::DeterministicSamples) = Statistics.cov(s.samples, dims=2)

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
module TestMeasures
    using Test
    using BioEnergeticFoodWebs
    using LinearAlgebra

    # Internal functions
    i = ones(10)
    @test BioEnergeticFoodWebs.shannon(i) ≈ 1.0

    i = zeros(100).+0.001
    i[1] = 1.0
    @test BioEnergeticFoodWebs.shannon(i) ≈ 0.0 atol=0.2

    @test isnan(BioEnergeticFoodWebs.shannon(vec([1.0])))
    @test isnan(BioEnergeticFoodWebs.shannon(vec([-1.0])))

    i = ones(5)
    @test BioEnergeticFoodWebs.coefficient_of_variation(i) == 0.0

    i = collect(range(0.0, stop = 1.0, length = 3))
    @test BioEnergeticFoodWebs.coefficient_of_variation(i) ≈ 1+1/(4*length(i))

    # Test the total biomass thing
    B = Matrix{Float64}(I, 10, 10)
    A = Matrix{Float64}(I, 10, 10)
    fake_params = Dict{Symbol, Any}(:A => A)
    s = Dict{Symbol, Any}(:B => B, :p => fake_params)
    @test total_biomass(s, last=10) == 1.0
    @test_throws AssertionError total_biomass(s, last=1000)
    @test population_biomass(s, last=10)[1] == 0.1

    # Population stability
    @test isnan(population_stability(s, last=1))
    @test population_stability(s, last=2) ≈ -1.59099  atol=0.01
    @test species_richness(s, last=1) == 1.0
    @test species_persistence(s, last=1) == 0.1

    # Foodweb evenness
    B = repeat([0 .5 0 .5], 5)
    s = Dict{Symbol, Any}(:B => B)
    @test foodweb_evenness(s, last = 2) == 1.00

    # Test when the total biomass is 0
    empty_p = Dict{Symbol, Any}(:B => zeros(10, 10))
    @test isnan(population_biomass(empty_p, last=2))
    @test isnan(population_stability(empty_p, last=2))
    @test isnan(total_biomass(empty_p, last=2))
    @test isnan(foodweb_evenness(empty_p, last=2))
    @test isnan(species_richness(empty_p, last = 2))

end

# module TestSave
#     using Test
#     using BioEnergeticFoodWebs
#     using LinearAlgebra
#     using JLD
#     using JSON
#
#     A = [0 0 0 ; 0 0 0 ; 0 0 0]
#     b = rand(3)
#     p = model_parameters(A)
#     #default variable name and extension
#     def_vname = "befwm_simul"
#     def_ext = :jld
#     #test default arguments
#     cd(tempdir())
#     s = simulate(p,b)
#     # default file name
#     def_fname = "befwm_" * string(hash(s)) * ".jld"
#     BioEnergeticFoodWebs.save(s, as = def_ext)
#
#     # Test if the file is saved (under the default name)
#     @test isfile(def_fname)
#     # Test if the content is the same
#     sbis = load(def_fname, def_vname)
#     @test sbis == s
#
#     rm(def_fname)
#     # Test is it works with as = :JLD
#     ext = :JLD
#     BioEnergeticFoodWebs.save(s, as = ext)
#     @test isfile(def_fname)
#     sbis = load(def_fname, def_vname)
#     @test sbis == s
#     rm(def_fname)
#
#     # Test with .json
#     ext = :json
#     BioEnergeticFoodWebs.save(s, as = ext)
#     fname = "befwm_" * string(hash(s)) * ".json"
#     @test isfile(fname)
#     rm(fname)
#
#     # Test with .JSON
#     ext = :JSON
#     BioEnergeticFoodWebs.save(s, as = ext)
#     fname = "befwm_" * string(hash(s)) * ".json"
#     @test isfile(fname)
#     rm(fname)
#
# end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
module nav_geo

export cosineSimilarity,deltaF1,deltaF2,destination,forwardAzmuth,haversine,speed

cosineSimilarity(a,b) = (a' * b)/(norm(a)*norm(b))

function deltaF1(args::Array{Float64})
   gradient = Array(Float64,2);

   r = args[1];
   theta = args[2];
   s = args[3];
   gradient[1] = -2*r*(s/(s + 1)*e^-(2*r/(s + 1))*(2 - 0.5*cos(theta + pi))  + 2*e^-(0.25*r^2)) + 1;
   gradient[2] = 0.5*s*e^-(2*r/(s + 1))*sin(theta);
   return gradient;
end

function deltaF2(args::Array{Float64})
   gradient = Array(Float64,2);

   r = args[1];
   theta = args[2];
   s = args[3];
   denom = abs(r -s*(2*cos(theta) + 5))
   dtheta = (2*s*sin(theta)*(r - s*(2*cos(theta) + 5)))/denom + r/4*sin(theta);
   dr = (r - 2*s*cos(theta) - 5*s)/denom - cos(theta)/4;
   gradient[1] = cos(theta)*dr - 1/r*sin(theta)*dtheta;
   gradient[2] = sin(theta)*dr + 1/r*cos(theta)*dtheta;
   return gradient;
end

function destination(args::Array{Float64})
   res = Array(Float64,2);
   srcLat = args[1]*pi/180;
   srcLong = args[2]*pi/180;
   azmuth = args[3]
   distanceInMeters = args[4];
   distanceInRadians = distanceInMeters/earthRadius(args[1]);
   destLat = asin(sin(srcLat)*cos(distanceInRadians) + cos(srcLat)*sin(distanceInRadians)*cos(azmuth));
   dlong = srcLong + atan2(sin(azmuth)*sin(distanceInRadians)*cos(srcLat),cos(distanceInRadians) - sin(srcLat)*sin(destLat));
   destLong = (dlong + 3*pi) % (2*pi) - pi;
   res[1] = destLat*180/pi;
   res[2] = destLong*180/pi;
   return res;
end

function earthRadius(lat)
   a =  6378137.0;
   b =  6356752.3;
   phi = lat*pi/180;
   return sqrt(((a^2*cos(phi))^2 + (b^2*sin(phi))^2)/((a*cos(phi))^2 + (b*sin(phi))^2));
end

function forwardAzmuth(args::Array{Float64})
   lat1 = args[1];
   long1 = args[2];
   lat2 = args[3];
   long2 = args[4];

   phi1 = lat1*pi/180;
   phi2 = lat2*pi/180;
   lam1 = long1*pi/180;
   lam2 = long2*pi/180;

   x = sin(lam2 - lam1)*cos(phi2);
   y = cos(phi1)*sin(phi2) - sin(phi1)*cos(phi2)*cos(lam2-lam1)
   return atan2(x,y) % (2*pi);
end

function haversine(args::Array{Float64})
   lat1 = args[1];
   long1 = args[2];
   lat2 = args[3];
   long2 = args[4];

   phi1 = lat1*pi/180;
   phi2 = lat2*pi/180;
   lam1 = long1*pi/180;
   lam2 = long2*pi/180;

   return 2 * earthRadius(lat1) * asin(sqrt(sin((phi2 - phi1)/2)^2 + cos(phi1) * cos(phi2) * sin((lam2 - lam1)/2)^2))
end

speed(x) = 7*(2/(1 + e^-(abs(x/10)^(5/4))) - 1)

end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey