text:
function f(n)
    s = n |> string |> collect |> sort
    g1 = parse(Int, foldl(*, reverse(s)))
    g2 = parse(Int, foldl(*, s))
    g1 - g2
end

function solve()
    n, k = [parse(Int, x) for x in split(readline())]
    for i=1:k
        n = f(n)
    end
    n
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# This file is a part of Pukeko.jl.
# License is MIT. https://github.com/IainNZ/Pukeko.jl

__precompile__(false)
module Pukeko

import Printf: @sprintf

export @test, @test_throws, @parametric

"""
    TEST_PREFIX

Functions with this string at the the start of their name will be treated as
self-contained sets of tests.
"""
const TEST_PREFIX = "test_"

"""
    TestException

The `Exception`` thrown when a Pukeko test fails. Used by `run_tests` to
distinguish between test errors and unexpected errors.
"""
struct TestException <: Exception
    message::String
    linfo::Union{Nothing, LineNumberNode}
end
function Base.show(io::IO, ex::TestException)
    if ex.linfo !== nothing
        print(io, basename(string(ex.linfo.file)), ":", ex.linfo.line, " ")
    end
    print(io, ex.message)
end

"""
    test_true(value, [linfo])

Throws iff `value` is not `true`. Calls to this are generated by `@test`.
"""
function test_true(value, linfo::Union{Nothing, LineNumberNode}=nothing)
    if value != true
        throw(TestException("Expression did not evaluate to `true`: " *
                            string(value), linfo))
    end
    return nothing
end

"""
    test_equal(value_left, value_right, [linfo])

Test that `value_left` is equal to `value_right`. Calls to this are
generated by `@test` for the case of `@test expr_left == expr_right`.
"""
function test_equal(value_left, value_right,
                    linfo::Union{Nothing, LineNumberNode}=nothing)
    if value_left != value_right
        throw(TestException("Expression did not evaluate to `true`: " *
                            string(value_left) * " != " *
                            string(value_right), linfo))
    end
    return nothing
end

"""
    @test expression

Test that `expression` is `true`.
"""
macro test(expression)
    # If `expression` is of form `expr_left == expr_right` -> `test_equal`.
    # Otherwise, use `test_true`.
    source = QuoteNode(__source__)
    if expression isa Symbol
        return quote
            test_true($(esc(expression)), $(source))
        end
    elseif (expression.head == :call && expression.args[1] == :(==) &&
            length(expression.args) == 3)
        return quote
            test_equal($(esc(expression.args[2])), $(esc(expression.args[3])),
                       $(source))
        end
    end
    return quote
        test_true($(esc(expression)), $(source))
    end
end

"""
    @test_throws exception_type expression

Test that `expression` throws an exception of type `exception_type`.
"""
macro test_throws(exception_type, expression)
    source = QuoteNode(__source__)
    return quote
        exception_thrown = true
        try
            $(esc(expression))
            exception_thrown = false
        catch exception
            expected_type = $(esc(exception_type))
            success = false
            if isa(expected_type, Type)
                success = exception isa expected_type
            else
                success = exception == expected_type
            end
            if !success
                throw(TestException("Expression threw exception of " *
                                    "type $(typeof(exception)), but " *
                                    "expected $(expected_type)",
                                    $(source)))
            end
        end
        if !exception_thrown
            throw(TestException("Expression did not throw an exception, " *
                                "expected $(expected_type) exception",
                                $(source)))
        end
    end
end

"""
    run_tests(module_to_test; fail_fast=false, timing=false,
              match_name=nothing, exclude_name=nothing)

Runs all the sets of tests in module `module_to_test`. Test sets are defined
as functions with names that begin with `TEST_PREFIX`. A summary is printed
after all test sets have been run and if there were any failures an
exception is thrown.

Configuration options:

  * If `fail_fast==false` (default), if any one test function fails, the
    others will still run. If `true`, testing will stop on the first
    failure. The commandline argument `--PUKEKO_FAIL_FAST` will override
    `fail_fast` to `true` for all `run_tests` calls.
  * If `timing==true` (default is `false`), print elapsed time and memory
    allocation statistics for every test function. The commandline
    argument `--PUKEKO_TIMING` will override `timing` to `true` for all
    `run_tests` calls.
  * If `match_name!=nothing` (default is `nothing`), only run tests that
    contain `match_name` in their names. The commandline argument
    `--PUKEKO_MATH=str` will override `match_name` to `str` for all
    `run_tests` calls.
  * If `exclude_name!=nothing` (default is `nothing`), do not tests that
    contain `exclude_name` in their names.
"""
function run_tests(module_to_test::Module;
    fail_fast::Bool=false, timing::Bool=false,
    match_name::Union{AbstractString, Nothing}=nothing,
    exclude_name::Union{AbstractString, Nothing}=nothing)
    # Parse commandline arguments.
    for arg in ARGS
        if arg == "--PUKEKO_FAIL_FAST"
            fail_fast = true
        elseif arg == "--PUKEKO_TIMING"
            timing = true
        elseif startswith(arg, "--PUKEKO_MATCH")
            match_name = split(arg, "=")[2]
        end
    end
    # Get a clean version of module name for logging messages.
    module_name = string(module_to_test)
    if startswith(module_name, "Main.")
        module_name = module_name[6:end]
    end
    # Keep track of failures to summarize at end.
    test_failures = Dict{String, TestException}()
    test_elapsed_time = Dict{String, UInt64}()
    test_start_mem = Dict{String, Base.GC_Num}()
    test_end_mem = Dict{String, Base.GC_Num}()
    test_functions = 0
    for maybe_function in names(module_to_test, all=true)
        maybe_function_name = string(maybe_function)
        # If not a function, skip to next
        if !(typeof(getfield(module_to_test, maybe_function)) <: Function)
             continue
        end
        # If not a test function, skip to next function.
        if !startswith(maybe_function_name, TEST_PREFIX)
            continue
        end
        # If it doesn't match, skip to next function.
        if match_name != nothing
            if !occursin(match_name, maybe_function_name)
                continue
            end
        end
        # If it does match, skip to next function.
        if exclude_name != nothing
            if occursin(exclude_name, maybe_function_name)
                continue
            end
        end
        # Track statistics for each test.
        test_functions += 1
        start_time = time_ns()
        if timing
            test_start_mem[maybe_function_name] = Base.gc_num()
        end
        # If we don't need to catch errors, don't even try.
        if fail_fast
            @eval module_to_test ($maybe_function)()
        else
            # Try to run the function. If it fails, figure out why.
            try
                @eval module_to_test ($maybe_function)()
            catch test_exception
                if isa(test_exception, TestException)
                    test_failures[maybe_function_name] = test_exception
                else
                    println("Unexpected exception occurred in test ",
                            "function `$(maybe_function_name)` ",
                            "in module `$(module_name)`")
                    throw(test_exception)
                end
            end
        end
        # Track statistics.
        if timing
            test_end_mem[maybe_function_name] = Base.gc_num()
        end
        test_elapsed_time[maybe_function_name] = time_ns() - start_time
    end
    # At least one test function failed, print out the exceptions.
    if length(test_failures) > 0
        println("Test failures occurred in module $(module_name)")
        println("Functions with failed tests:")
        for (function_name, test_exception) in test_failures
            println("    $(function_name): ", test_exception)
        end
        error("Some tests failed!")
    end
    # All passed, output statistics.
    total_time = sum(values(test_elapsed_time))
    println("$(test_functions) test function(s) ran successfully ",
            "in module $(module_name) ",
            @sprintf("(%.2f seconds)", total_time / 1e9))
    if timing
        # Sort by run time, descending.
        time_names = [(elapsed, function_name)
                      for (function_name, elapsed) in test_elapsed_time]
        sort!(time_names, rev=true)
        for (elapsed, function_name) in time_names
            gc_diff = Base.GC_Diff(test_end_mem[function_name],
                                   test_start_mem[function_name])
            print(function_name, ": ")
            Base.time_print(elapsed, gc_diff.allocd, gc_diff.total_time,
                            Base.gc_alloc_count(gc_diff))
            println()
        end
    end
end

"""
    parametric(module_to_test, func, iterable)

Create a version of `func` that is prefixed with `TEST_PREFIX` in
`module_to_test` for each value in `iterable`. If a value in `iterable` is
a tuple, it is splatted into the function arguments.
"""
function parametric(module_to_test, func, iterable)
    for (index, value) in enumerate(iterable)
        func_name = Symbol(string(TEST_PREFIX, func, index))
        if value isa Tuple
            @eval module_to_test $func_name() = $func($(value)...)
        else
            @eval module_to_test $func_name() = $func($(value))
        end
    end
end

"""
    @parametric func iterable

Create a version of `func` that is prefixed with `TEST_PREFIX` in the module
that this macro is called for each value in `iterable`. If a value in
`iterable` is a tuple, it is splatted into the function arguments.
"""
macro parametric(func, iterable)
    module_ = __module__
    return quote
        parametric($(module_), $(esc(func)), $(esc(iterable)))
    end
end

end

ext: jl


text:
# prior factor that passes a density belief straight through to inference without resampling

export PartialPriorPassThrough, PackedPartialPriorPassThrough

struct PartialPriorPassThrough{B <: Union{<:HeatmapGridDensity,<:LevelSetGridNormal}, T <:Tuple} <: AbstractPrior
  Z::B
  partial::T
end

getManifold(pppt::PartialPriorPassThrough) = getManifold(pppt.Z)

# this step is skipped during main inference process
getSample(cf::CalcFactor{<:PartialPriorPassThrough}) = sampleTangent(getManifold(cf.factor), cf.factor.Z)

## ====================================================================================================
## Serialize PartialPriorPassThrough
## ====================================================================================================

"""
    $TYPEDEF

Required internal density to store its type
"""
Base.@kwdef mutable struct PackedPartialPriorPassThrough <: AbstractPackedFactor
  Z::PackedSamplableBelief # PackedHeatmapGridDensity
  partial::Vector{Int}
end

function convert( ::Union{Type{<:AbstractPackedFactor}, Type{<:PackedPartialPriorPassThrough}},
                  obj::PartialPriorPassThrough )
  #

  po = convert(PackedSamplableBelief, obj.Z)
  PackedPartialPriorPassThrough(po, Int[obj.partial...])
end

function convert( ::Union{Type{<:AbstractFactor}, Type{<:PartialPriorPassThrough}},
                  obj::PackedPartialPriorPassThrough )
  #

  dens = convert(SamplableBelief, obj.Z)
  PartialPriorPassThrough(dens, tuple(obj.partial...))
end

#


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function ps_m5_opt(A)
    return ps_m5_opt!(copy(A))
end

@inline function ps_m5_opt!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    max_memslots=8
    memslots=Vector{Matrix{T}}(undef,max_memslots)
    n=size(A,1)
    for j=1:max_memslots
        memslots[j]=Matrix{T}(undef,n,n)
    end
    # The first slots are precomputed nodes [:A]
    memslots[1]=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=-0.15614734226906415
    coeff2=1.0138950449120168
    memslots[2] .= coeff2.*memslots[1]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.1561459961705319
    coeff2=1.0138962847063822
    memslots[3] .= coeff2.*memslots[1]
    mul!(memslots[3],true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots[4],memslots[3],memslots[2])
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-1.1763384175098646
    coeff2=0.14819920360945626
    coeff3=0.7004279880959638
    memslots[2] .= coeff2.*memslots[1] .+ coeff3.*memslots[4]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=-0.15896658336808092
    coeff2=1.03474409787952
    coeff3=-0.6899181523788153
    memslots[3] .= coeff2.*memslots[1] .+ coeff3.*memslots[4]
    mul!(memslots[3],true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots[5],memslots[2],memslots[3])
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.32868543911603565
    coeff2=0.04948779256327628
    coeff3=1.1645163087525163
    coeff4=0.33379574828605685
    memslots[2] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.2715485387807492
    coeff2=0.4515173495475029
    coeff3=-0.45819098373171513
    coeff4=0.5179573824690393
    memslots[3] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5]
    mul!(memslots[3],true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots[6],memslots[2],memslots[3])
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.1784350608747074
    coeff2=-0.9109875948763699
    coeff3=0.20063142586799526
    coeff4=-1.6463200998393868
    coeff5=0.8066435555717062
    memslots[2] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5] .+ coeff5.*memslots[6]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.9977699469177798
    coeff2=0.7154077419527756
    coeff3=-1.0589444472962166
    coeff4=0.32574547503223705
    coeff5=0.04488529740138668
    memslots[3] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5] .+ coeff5.*memslots[6]
    mul!(memslots[3],true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots[7],memslots[2],memslots[3])
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.7570005286093737
    coeff2=0.007598464345057781
    coeff3=-0.8669783274817463
    coeff4=-0.14407160251709974
    coeff5=-0.6983909370796325
    coeff6=1.3011465317844473
    memslots[2] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5] .+ coeff5.*memslots[6] .+ coeff6.*memslots[7]
    mul!(memslots[2],true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.6207474625776077
    coeff2=-0.8881294075403886
    coeff3=0.24425285144993336
    coeff4=-0.9312258295819543
    coeff5=-1.0814610150998485
    coeff6=1.0407649022964955
    memslots[3] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5] .+ coeff5.*memslots[6] .+ coeff6.*memslots[7]
    mul!(memslots[3],true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots[8],memslots[2],memslots[3])
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.2989395453295784
    coeff2=2.080346787066744
    coeff3=-1.5251262809722197
    coeff4=0.004244973897574398
    coeff5=-0.0326540011916379
    coeff6=0.8226732267898158
    coeff7=0.7171476120487754
    # Smart lincomb recycle A
    memslots[1] .= coeff2.*memslots[1] .+ coeff3.*memslots[4] .+ coeff4.*memslots[5] .+ coeff5.*memslots[6] .+ coeff6.*memslots[7] .+ coeff7.*memslots[8]
    mul!(memslots[1],true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots[1] # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using ITensors
using ITensors.NDTensors
using LinearAlgebra
using Test

@testset "diagITensor" begin
  d = 3
  i = Index(d, "i")
  j = Index(d, "j")
  k = Index(d, "k")
  l = Index(d, "l")
  m = Index(d, "m")
  n = Index(d, "n")
  o = Index(d, "o")
  p = Index(d, "p")
  q = Index(d, "q")

  v = collect(1:d)
  vr = randn(d)

  @testset "non-uniform diagonal values" begin
    @testset "diagITensor constructor (no vector, order 2)" begin
      D = diagITensor(i, j)

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d
        if ii == jj
          @test D[i => ii, j => jj] == 0.0
        else
          @test D[i => ii, j => jj] == 0.0
        end
      end
    end

    @testset "diagITensor constructor (no vector, order 3)" begin
      D = diagITensor(i, j, k)

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test D[i => ii, j => jj, k => kk] == 0.0
        else
          @test D[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "diagITensor constructor (no vector, complex)" begin
      D = diagITensor(ComplexF64, i, j)

      @test eltype(D) == ComplexF64
      for ii in 1:d, jj in 1:d
        if ii == jj
          @test D[i => ii, j => jj] == complex(0.0)
        else
          @test D[i => ii, j => jj] == complex(0.0)
        end
      end
    end

    @testset "diag" for ElType in (Float64, ComplexF64)
      A = diagITensor(randn(ElType, d), i, j)
      dA = diag(A)
      @test dA isa DenseTensor{ElType,1}
      @test dA[1] == A[1, 1]
      @test dA[2] == A[2, 2]
    end

    @testset "diagITensor constructor (vector, order 2)" begin
      D = diagITensor(v, i, j)

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d
        if ii == jj
          @test D[i => ii, j => jj] == v[ii]
        else
          @test D[i => ii, j => jj] == 0.0
        end
      end
    end

    @testset "diagITensor constructor (vector, order 3)" begin
      D = diagITensor(v, i, j, k)

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test D[i => ii, j => jj, k => kk] == v[ii]
        else
          @test D[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "diagITensor constructor (complex)" begin
      vc = v + im * v
      D = diagITensor(vc, i, j, k)

      @test eltype(D) == ComplexF64
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test D[i => ii, j => jj, k => kk] == vc[ii]
        else
          @test D[i => ii, j => jj, k => kk] == complex(0.0)
        end
      end
    end

    @testset "Complex operations" begin
      xr = randn(d)
      xi = randn(d)
      D = diagITensor(xr + im * xi, i, j, k)
      @test eltype(D) == ComplexF64
      rD = real(D)
      iD = imag(D)
      @test eltype(rD) == Float64
      @test eltype(iD) == Float64
      @test typeof(storage(rD)) <: NDTensors.Diag
      @test norm(rD + im * iD - D) < 1E-8
    end

    @testset "Constructor AllowAlias/NeverAlias" begin
      vv = ones(d)
      D = diagITensor(vv, i, j)
      @test eltype(D) === Float64
      D[1, 1] = 5.0
      @test vv[1] == 1.0
      @test vv[1] != D[1, 1]

      vv = ones(Int, d)
      D = diagITensor(vv, i, j)
      @test eltype(D) === Float64
      D[1, 1] = 5.0
      @test vv[1] == 1.0
      @test vv[1] != D[1, 1]

      vv = ones(Int, d)
      D = diagITensor(Int, vv, i, j)
      @test eltype(D) === Int
      D[1, 1] = 5
      @test vv[1] == 1
      @test vv[1] != D[1, 1]

      vv = ones(d)
      D = diagitensor(vv, i, j)
      @test eltype(D) === Float64
      D[1, 1] = 5.0
      @test vv[1] == 5.0
      @test vv[1] == D[1, 1]

      vv = ones(Int, d)
      D = diagitensor(vv, i, j)
      @test eltype(D) === Float64
      D[1, 1] = 5.0
      @test vv[1] == 1.0
      @test vv[1] != D[1, 1]

      vv = ones(Int, d)
      D = diagitensor(Int, vv, i, j)
      @test eltype(D) === Int
      D[1, 1] = 5
      @test vv[1] == 5
      @test vv[1] == D[1, 1]

      D = diagITensor(1, i, j)
      @test eltype(D) === Float64
      D[1, 1] = 5
      @test D[1, 1] == 5

      D = diagITensor(Int, 1, i, j)
      @test eltype(D) === Int
      D[1, 1] = 5
      @test D[1, 1] == 5
    end

    @testset "fill!" begin
      D = diagITensor(ones(d), i, j, k)
      D = fill!(D, 2.0)
      for ii in 1:d
        @test D[i => ii, j => ii, k => ii] == 2.0
      end

      @test eltype(D) == Float64
    end

    @testset "Set elements" begin
      D = diagITensor(i, j, k)

      for ii in 1:d
        D[i => ii, j => ii, k => ii] = ii
      end

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test D[i => ii, j => jj, k => kk] == ii
        else
          @test D[i => ii, j => jj, k => kk] == 0.0
        end
      end

      # Can't set off-diagonal elements
      @test_throws ErrorException D[i => 2, j => 1, k => 1] = 0.0
      @test_throws ErrorException D[i => 1, j => 2, k => 1] = 0.0
    end

    @testset "Convert diag to dense" begin
      D = diagITensor(v, i, j, k)
      T = dense(D)

      @test storage(T) isa NDTensors.Dense{Float64}
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test T[ii, ii, ii] == ii
        else
          @test T[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "Convert diag to dense with denseblocks" begin
      D = diagITensor(v, i, j, k)
      T = denseblocks(D)

      @test storage(T) isa NDTensors.Dense{Float64}
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test T[ii, ii, ii] == ii
        else
          @test T[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "Add (Diag + Diag)" begin
      v1 = randn(d)
      v2 = randn(d)
      D1 = diagITensor(v1, i, j, k)
      D2 = diagITensor(v2, k, i, j)

      v3 = v1 + v2
      D3 = D1 + D2

      @test D3 ≈ dense(D1) + dense(D2)
      for ii in 1:d
        @test D3[ii, ii, ii] == v3[ii]
      end
    end

    @testset "Add ( number * Diag + Diag)" begin
      v1 = randn(d)
      v2 = randn(d)
      D1 = Float32(2.0) * diagITensor(v1, i, j, k)
      D2 = diagITensor(v2, k, i, j)

      v3 = 2 * v1 + v2
      D3 = D1 + D2

      @test D3 ≈ dense(D1) + dense(D2)
      for ii in 1:d
        @test D3[ii, ii, ii] == v3[ii]
      end
    end

    @testset "Add (Diag uniform + Diag uniform)" begin
      D1 = δ(i, j, k)
      D2 = δ(k, i, j)

      D3 = D1 + D2

      @test D3 ≈ dense(D1) + dense(D2)
    end

    @testset "Add (Diag + Dense)" begin
      D = diagITensor(vr, i, j, k)
      A = randomITensor(k, j, i)

      R = D + A

      @test R ≈ dense(D) + A
      for ii in 1:d
        @test R[ii, ii, ii] ≈ D[ii, ii, ii] + A[ii, ii, ii]
      end
    end

    @testset "Add (Dense + Diag)" begin
      D = diagITensor(vr, i, j, k)
      A = randomITensor(i, k, j)

      R = A + D

      @test R ≈ dense(D) + A
      for ii in 1:d
        @test R[ii, ii, ii] ≈ D[ii, ii, ii] + A[ii, ii, ii]
      end
    end

    @testset "Contraction (all contracted)" begin
      D = diagITensor(v, i, j, k)
      A = randomITensor(j, k, i)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (all contracted) with different types" begin
      D = diagITensor(v, i, j, k)
      A = randomITensor(Float32, j, k, i)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A

      D = diagITensor(v, i, j, k)
      A = randomITensor(ComplexF32, j, k, i)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (all dense contracted)" begin
      D = diagITensor(v, j, k, i)
      A = randomITensor(i, j)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction Diag*Dense (general)" begin
      D = diagITensor(v, l, i, k, j)
      A = randomITensor(m, k, n, l)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction Diag*Dense (outer)" begin
      D = diagITensor(v, l, i, k, j)
      A = randomITensor(m, n)

      @test order(D * A) == 6
      @test D * A ≈ dense(D) * A
    end

    @testset "Contraction Diag*Diag (outer)" begin
      D1 = diagITensor(v, l, i)
      D2 = diagITensor(v, m, n)

      @test order(D1 * D2) == 4
      @test D1 * D2 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction Diag*Diag (all contracted)" begin
      D1 = diagITensor(v, l, i, k, j)
      D2 = diagITensor(vr, j, l, i, k)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction Diag*Diag (general)" begin
      D1 = diagITensor(v, l, i, k, j)
      D2 = diagITensor(vr, m, k, n, l)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction Diag*Diag (no contracted)" begin
      D1 = diagITensor(v, i, j)
      D2 = diagITensor(vr, k, l)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction Diag*Scalar" begin
      D = diagITensor(v, i, j)
      x = 2.0

      @test x * D ≈ x * dense(D)
      @test D * x ≈ x * dense(D)

      xc = 2 + 3im

      @test xc * D ≈ xc * dense(D)
      @test D * xc ≈ xc * dense(D)
    end
  end

  @testset "Uniform diagonal ITensor" begin
    @testset "delta constructor (order 2)" begin
      D = δ(i, j)

      @test eltype(D) == Float64
      for ii in 1:d, jj in 1:d
        if ii == jj
          @test D[i => ii, j => jj] == 1.0
        else
          @test D[i => ii, j => jj] == 0.0
        end
      end
    end

    @testset "delta constructor (order 3)" begin
      D = δ(i, j, k)

      @test eltype(D) == Float64
      for ii in 1:dim(i), jj in 1:dim(j), kk in 1:dim(k)
        if ii == jj == kk
          @test D[i => ii, j => jj, k => kk] == 1.0
        else
          @test D[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "Set elements" begin
      D = δ(i, j, k)

      @test eltype(D) == Float64

      # Can't set elements of uniform diag tensor
      # TODO: should we make a function that converts
      # to a version that can?
      @test_throws ErrorException D[i => 1, j => 1, k => 1] = 2.0
      @test_throws ErrorException D[i => 2, j => 1, k => 1] = 4.3
      @test_throws ErrorException D[i => 1, j => 2, k => 1] = 2
    end

    @testset "Convert diag uniform to dense" begin
      D = δ(i, j, k)
      T = dense(D)

      @test storage(T) isa NDTensors.Dense{Float64}
      for ii in 1:d, jj in 1:d, kk in 1:d
        if ii == jj == kk
          @test T[ii, ii, ii] == 1.0
        else
          @test T[i => ii, j => jj, k => kk] == 0.0
        end
      end
    end

    @testset "Add (Diag uniform + Dense)" begin
      D = δ(i, j, k)
      A = randomITensor(k, j, i)

      R = D + A

      @test R ≈ dense(D) + A
      for ii in 1:d
        @test R[ii, ii, ii] ≈ D[ii, ii, ii] + A[ii, ii, ii]
      end
    end

    @testset "Contraction (Diag uniform * Dense, all contracted)" begin
      D = δ(i, j, k)
      A = randomITensor(j, k, i)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (Diag uniform * Dense, all dense contracted)" begin
      D = δ(j, k, i)
      A = randomITensor(i, j)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (Diag uniform * Dense, general)" begin
      D = δ(l, i, k, j)
      A = randomITensor(m, k, n, l)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction with different bond dimensions" begin
      idim = 3
      mdim = 2

      i = Index(idim, "i")
      m = Index(mdim, "m")

      A = randomITensor(i, i', m)
      D = δ(i, i')

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (Diag uniform * Dense, replace index)" begin
      D = δ(i, k)
      A = randomITensor(m, k, n, l)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (Diag uniform * Dense, replace index 2)" begin
      D = δ(k, i)
      A = randomITensor(m, n, k, l)

      @test D * A ≈ dense(D) * A
      @test A * D ≈ dense(D) * A
    end

    @testset "Contraction (Diag uniform * Diag uniform, all contracted)" begin
      D1 = δ(l, i, k, j)
      D2 = δ(j, l, i, k)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction (Diag uniform * Diag uniform, general)" begin
      D1 = δ(l, i, k, j)
      D2 = δ(m, k, n, l)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end

    @testset "Contraction (Diag uniform * Diag uniform, no contracted)" begin
      D1 = δ(i, j)
      D2 = δ(k, l)

      @test D1 * D2 ≈ dense(D1) * dense(D2)
      @test D2 * D1 ≈ dense(D1) * dense(D2)
    end
  end
end

nothing


text:
#################################################################
#
#    1st order derivation testing (getindex-setindex)
#
#################################################################

##  ref  testing
@compare x[2]                v1ref
@compare sum(x[2:3])         v1ref
@compare sum(x[2:4])         v2ref

@compare sum(x[:,2])         v2ref
@compare sum(x[1,:])         v2ref
@compare sum(x[2:end,:])     v2ref
@compare sum(x[:,2:end])     v2ref
@compare sum(x[1:end-1,:])   v2ref
_idx1=2
@compare sum(x[_idx1:end,:]) v2ref

@compare getindex(x,2)          v1ref
@compare sum(getindex(x, 2:3))  v1ref
@compare sum(getindex(x, 2:4))  v2ref
@compare sum(getindex(x, :,2))  v2ref
@compare sum(getindex(x, 1,:))  v2ref

@compare x[2]+x[1]           v2ref
@compare log(x[2]^2+x[1]^2)  v2ref

@compare (a = zeros(5) ; a[1]   = x[3]   ; a[3]) v1ref
@compare (a = zeros(5) ; a[1]   = x[3]   ; a[1]) v1ref
@compare (a = zeros(5) ; a[1:2] = x[3]   ; a[1]) v1ref
@compare (a = zeros(5) ; a[1:2] = x[3:4] ; a[1]) v1ref
@compare (a = zeros(5) ; a[1:2] = x[3:4] ; a[4]) v1ref

@compare (a = zeros(5) ; setindex!(a, x[3]  , 1)   ; a[3]) v1ref
@compare (a = zeros(5) ; setindex!(a, x[3]  , 1)   ; a[1]) v1ref
@compare (a = zeros(5) ; setindex!(a, x[3]  , 1:2) ; a[1]) v1ref
@compare (a = zeros(5) ; setindex!(a, x[3:4], 1:2) ; a[1]) v1ref
@compare (a = zeros(5) ; setindex!(a, x[3:4], 1:2) ; a[4]) v1ref

_idx2 = 1

ex = quote
    a = zeros(2,2)
    b = fill(0.0, size(x))
    b[2] = b[2] + x[2]
    b[1] = b[1] + x[1]
    a[1:2] = b
    a[_idx2]
end
compare(ex, ones(2))

ex = quote
    a = zeros(2,2)
    a[1:2] = x
    a[_idx2]
end
compare(ex, ones(2))
compare(ex, 1.)

ex = quote
    a = zeros(2,2)
    a[1:2] = x[1]
    a[_idx2]
end
compare(ex, ones(2))

ex = quote
    a = ones(5)
    b = sum(a)*x
    a[2] += x
    c = sum(a)
    b + c
end
compare(ex, 2.)

ex = quote
    a=zeros(3)
    b=zeros(3)
    b[2]=x
    a[1]=x
    sum(a)+sum(b)
end
compare(ex, 2.)

ex = quote
    a=zeros(5)
    a[2:3] = x
    a[3:4] = x
    sum(a)
end
compare(ex, -1.)

ex = quote
    a = zeros(2,2)
    b = zeros(2)

    b[2] = b[2] + x[2]
    b[1] = b[1] + x[1]
    a[1:2] = b

    b[2] = b[2] + x[2]
    b[1] = b[1] + x[1]
    a[3:4] = b

    a[_idx2]
end
_idx2 = 1
compare(ex, [1., 1.])
compare(ex, [0., 1.])
_idx2 = 3
compare(ex, [1., 1.])
compare(ex, [0., 1.])

ex = quote
    a = zeros(2,2)
    b = zeros(2)

    b[1] = b[1] + x[1]
    b[2] = b[2] + x[2]
    a[1] = b[1]
    a[2] = b[2]

    b[2] = b[2] + x[2]
    b[1] = b[1] + x[1]
    a[3] = b[1]
    a[4] = b[2]

    a[_idx2]
end
_idx2 = 1
compare(ex, [1., 1.])   #   pb : super slow
compare(ex, [0., 1.])
_idx2 = 3
compare(ex, [1., 1.])
compare(ex, [0., 1.])

ex = quote
    a = zeros(4)
    b = zeros(2)

    b += x
    a[1:2] = b
    b[1] = b[1] + x[1]
    a[3:4] = b

    a[_idx2]
end
compare(ex, [1., 1.])
compare(ex, [1., 0.])

ex = quote
    a = zeros(2,2)
    a[1] = x[1]
    a[3] = x[1]
    a[4] = x[1]
    a[1] = x[1]
    a[_idx2]
end
compare(ex, [1., 1.])
compare(ex, ones(10))

ex = quote
    a = zeros(5,5)
    a[1:4,2] = x
    a[3,1:5] = x
    sum(a)
end
compare(ex, 1.)


text:
using Distributions

f(x) = x^2 - 5 * x + 6

function get_best_solution(soln)
    best = soln[1]
    for s in soln
        if f(s) < f(best)
            best = s
        end
    end
    return best
end

function flower_pollination_algorithm(f;total_flowers=25,λ = 1.5,γ = 0.1,p = 0.8, MaxGenerations = 10)
    flowers = rand(Float64, total_flowers)
    L = Levy()

    global_best = get_best_solution(flowers)
    print("Start Global best:", global_best)
    println("Out of answers:", flowers)
    println("Best answer so far:", global_best)

    t = 1
    while t < MaxGenerations
        for i = 1:total_flowers
            new_value = nothing
            if rand() < p
                new_value = flowers[i] + γ * rand(L) * (global_best - flowers[i])
            else
                xj = rand(1:total_flowers)
                xk = rand(1:total_flowers)
                new_value = flowers[i] + rand() * (xj - xk)
            end
            if f(new_value) < f(flowers[i])
                flowers[i] = new_value
            end
        end
        global_best = get_best_solution(flowers)
        println("Global Best so far: ", global_best)
        t = t + 1
    end
    return global_best
end

global_minimum=flower_pollination_algorithm(f,total_flowers=25,λ = 1.5,γ = 0.1,p = 0.8)

println("Final Global Minimum: ", global_minimum)


text:
# using PyCall
# using JLD
using PyCallJLD
using MetaICVI
using ClusterValidityIndices
using Test
using Logging

# Include some test utilities (data loading, etc.)
include("test_utils.jl")

# MetaICVI module testing
@testset "MetaICVI.jl" begin
    # Create the module
    opts = MetaICVIOpts(
        fail_on_missing = true
    )
    metaicvi = MetaICVIModule(opts)

    # Display some aspects of the module
    @info fieldnames(typeof(metaicvi))
    @info metaicvi
    @info metaicvi.classifier

    # Point to the correct data directories
    data_dir(args...) = joinpath("../data/testing", args...)
    results_dir(args...) = joinpath("../data/results", args...)

    # Load the data and test across all supervised modules
    data = load_iris(data_dir("Iris.csv"))
    data.train_y = relabel_cvi_data(data.train_y)

    # Iterate over the data
    n_data = length(data.train_y)
    performances = zeros(n_data)
    for i = 1:n_data
        sample = data.train_x[:, i]
        label = data.train_y[i]
        performances[i] = get_metaicvi(metaicvi, sample, label)
    end

    # Perform some simple tests
    @test all(performances .>= 0)
    @test all(performances .<= 1)
end

# Rocket testing
@testset "Rocket.jl" begin
    # Test saving and loading
    filepath = "my_rocket"

    # Create the rocket module
    my_rocket = MetaICVI.Rocket.RocketModule()

    # Save the module to the filepath
    MetaICVI.Rocket.save_rocket(my_rocket, filepath)

    # Test that the file exists
    @test isfile(filepath)

    # Load the rocket module
    my_new_rocket = MetaICVI.Rocket.load_rocket(filepath)

    # Delete the saved file
    rm(filepath)
end


text:
## Posterior distribution of the flexible links model (applied to all food webs archived in Mangal)

# counterfactuals and posterior distribution of flexible links model
fl_post = CSV.read(joinpath("data", "sims", "beta_binomial_posterior.csv"))

# parameters mu and phi from the entire joint posterior distribution
mu_post = fl_post[:,:mu]
phi_post = exp.(fl_post[:,:phi])

max_stars_repo_name: FrancisBanville/Trophic-METE

max_issues_repo_name: FrancisBanville/Trophic-METE

max_forks_repo_name: FrancisBanville/Trophic-METE


text:
#=
struct ExecuteFieldParallel
	#enter::Function
	#leave::Function
	result::Dict

	function ExecuteFieldParallel(resolvers)
      result=Dict()
		function strands_execution(strand)
			#=

			 query{
	      event{
	        attends{
	          name
	        }
	      }
        }

			la hebra seria:  query - event - attends - name
			se ejecuta el resolve   resolve["event"]["attends"] -> se alamacena el resultado en una cache para evitar
			que se ejecuten muchas veces la misma consulta
			si ya tienes los datos solo pasalos

			y se pasa a la hoja "name" en lo campo root para la ejecucion de la hoja
			se evita la ejecucion del resolve  resolve["query"]["event"]

            soluciona el Multi-layered data fetching

			aun pienso en como resolver el n+1 con este mismo enfoque

	    query{
	      event{ -> []
	      	name
	        attends{ -> []
	          name
	        }
	      }
        }

			tal vez :

			query - event - name =>  toda la tabla event

			 query - event - attends - name => toda la tabla attends

			 y despues que se fusionen los datos, pero quien sabe
			=#
		end

		new(result)

    end

end
=#

struct ExecuteField
    visitante
    datos::Dict

    function ExecuteField(symbol_table::Dict, resolvers::Dict,ctx; Variables=nothing)
    	global datos=Dict("data"=>Dict())

    	function resolvefieldValue(nombrecampo::String,args::Dict,type_padre::String,root::Union{Nothing,Dict},path::String,tipoactual::String)
				if (haskey(resolvers, type_padre)) && (haskey(resolvers[type_padre], nombrecampo))
						return resolvers[type_padre][nombrecampo](root,args,ctx,Dict("fieldName"=>nombrecampo,"parentType"=>type_padre,"path"=>path,"dato=Type"=> tipoactual))
				elseif (typeof(root)<:Dict) # si root es un diccionario
			        	return root[nombrecampo]
			    else
					return "null"
				end
		end

    	function mivisitante(x::Array{Node,1},path::String,type_padre::String,root::Union{Nothing,Dict})
    		for c in x
		     mivisitante(c,path,type_padre,root)
			end
		end

		function mivisitante(node::Node,path::String,type_padre::String,root::Union{Nothing,Dict})
			t = typeof(node)

            if (node.kind == "Field")

            	nombre_nodo = node.name.value
            	tipoactual= ""

            	    if (haskey(symbol_table[type_padre], nombre_nodo))
			        	tipoactual= symbol_table[type_padre][nombre_nodo]["tipo"]
			        else
			        	throw(GraphQLError("{\"data\": null,\"errors\": [{\"message\": \"The $(nombre_nodo) field does not exists.\"}]}"))
					end

                args=Dict()
			        			if length(node.arguments)>0
			        			   for data in node.arguments
			        			   	if typeof(data.value[2])<:Diana.Object_
			        			   		push!(args, data.name.value => Dict() )
			        			   		for i in data.value[2].fields
			        			   		    if typeof(i.value[2])<:Diana.Variable
			        			   		    	push!(args[data.name.value], i.name.value => Variables[i.value[2].name.value] )
			        			   			else
			        			   		    	push!(args[data.name.value], i.name.value => i.value[2].value )
			        			   			end
			        			   		end
			        			   	else
			        			   			if typeof(data.value[2])<:Diana.Variable
			        			   		    	push!(args, data.name.value => Variables[data.value[2].name.value] )
			        			   			else
			        			   		    	push!(args, data.name.value => data.value[2].value )
			        			   			end
			        			   	end
			        			   end
			        		    end

            		if ((tipoactual == "String") | (tipoactual == "ID"))
            			eval(Meta.parse("push!(datos$(path), \"$(nombre_nodo)\" => \"$(resolvefieldValue(nombre_nodo,args,type_padre,root,path,tipoactual))\")"))
            		elseif ( (tipoactual == "Float") | (tipoactual == "Int") | (tipoactual == "Boolean"))
                             eval(Meta.parse("push!(datos$(path), \"$(nombre_nodo)\" => $(resolvefieldValue(nombre_nodo,args,type_padre,root,path,tipoactual)))"))
			        else # el tipo lista [ ] apenas lo pondre

			        	if (haskey(resolvers, type_padre)) && (haskey(resolvers[type_padre], nombre_nodo))

                                root= resolvers[type_padre][nombre_nodo](nothing,args,ctx,Dict("fieldName"=>nombre_nodo,"parentType"=>type_padre,"path"=>path,"returnType"=> tipoactual))
			  	        end

						if !(nombre_nodo in collect(keys(eval(Meta.parse("datos$(path)")))))
            				eval(Meta.parse("push!(datos$(path), \"$(nombre_nodo)\"=>Dict())"))
                        end

            			path= path*"[\"$(nombre_nodo)\"]"
			  	        type_padre= tipoactual

	                end
            end

			for f in fieldnames(t)
				subarbol= getfield(node, f)
				if( (typeof(subarbol )<: Node) & !(isa(subarbol, Diana.Name) ) )
					mivisitante(subarbol,path,type_padre,root)
				elseif (typeof(subarbol) <: Array{Node,1})
				    mivisitante(subarbol,path,type_padre,root)
				end
			end
			#rules.leave(x)
		end

	    function visitante(ast::Node, raiz::String)
			mivisitante(ast,"[\"data\"]",raiz,nothing)
	    end
      new(visitante,datos)
    end
end

function ExecuteQuery(operation::Node, symbol_table::Dict, resolvers::Dict,context; Variables=nothing, initialValue=nothing)
    resultexec = ExecuteField(symbol_table,resolvers,context,Variables=Variables)
    resultexec.visitante(operation,symbol_table["query"])
    return resultexec.datos
end

function ExecuteMutation(operation::Node, symbol_table::Dict, resolvers::Dict,context; Variables=nothing, initialValue=nothing)
	resultexec = ExecuteField(symbol_table,resolvers,context,Variables=Variables)
    resultexec.visitante(operation,symbol_table["mutation"])
    return resultexec.datos
end

#function Subscribe(operation::Node, symbol_table::Dict; Variables=nothing, initialValue=nothing)

#end

#function ExecuteQueryParallel(operation::Node, symbol_table::Dict, resolvers; Variables=nothing, initialValue=nothing)

#end


text:
"""
    function koala

From the canonical koala collection

# Example

```jldoctest
julia> cowsay("Welcome to Down Under", cow=Cowsay.koala)
 _______________________
< Welcome to Down Under >
 -----------------------
  \\
   \\
       ___
     {~._.~}
      ( Y )
     ()~*~()
     (_)-(_)
```
"""
function koala(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts
   $thoughts
       ___
     {~._.~}
      ( Y )
     ()~*~()
     (_)-(_)
"""
return the_cow
end


text:
include("filtering.jl");
## Protocols Decay

Df = Prew(1:20)
Protocol_colors!(Df)
@df Df plot(:Poke, :Prew,
    group = :env, linecolor = :color)
@df Df scatter!(:Poke, :Prew,
    group = :env, color = :color)
savefig(joinpath(figs_loc,"Fig2/A_ProtocolsDecay.pdf"))

## Rewards per trial

Df = combine(groupby(streaks,:Treatment)) do dd
    summarize(dd,:Protocol,:Num_Rewards)
end
rename!(Df, :Xaxis=>:Protocol)
Protocol_colors!(Df)
gd = groupby(Df,:Treatment)
tp = [@df subdf bar(:Protocol, :Mean,
        xlabel = :Treatment[1],
        yerror = :SEM,
        color = :color,
        legend = false,
        yaxis = "Rewards per trial") for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/B_RewardsPerProtocol.pdf"))

## Marginal value theorem scatter plot
MVT(pokes)
savefig(joinpath(figs_loc,"Fig2/D1_reward_rate.pdf"))
MVT(pokes; group =:Phase)
savefig(joinpath(figs_loc,"Fig2/Phase_reward_rate.pdf"))

## Marginal value theorem bar plot
gd = groupby(pokes,[:Treatment,:MouseID,:Day,:Trial])
Rrate = combine([:InstRewRate, :AverageRewRate, :Reward] => (i,a,r) -> (Leaving = i[end],
    Average = a[end],
    Reward = r[end]),gd)
filter!(r -> r.Reward, Rrate)
df1 = stack(Rrate,[:Leaving,:Average])
rename!(df1, [:variable => :Rate_on, :value => :InstRewRate])
df2 = combine(groupby(df1,[:Treatment])) do dd
    summarize(dd,:Rate_on,:InstRewRate)
end
Drug_colors!(df2)
gd = groupby(df2,:Treatment)
tp = [@df subdf bar(string.(:Xaxis), :Mean,
        xlabel = :Treatment[1],
        yerror = :SEM,
        # ylims = (0,0.7),
        color = :color,
        label = false,
        yaxis = "Rewards rate") for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/D2_reward_rate.pdf"))

## Marginal value theorem SignedRankTest
df1 = combine(groupby(pokes,:Treatment)) do dd
    wilcoxon(dd,:InstRewRate,:AverageRewRate)
end
Drug_colors!(df1)
plot_wilcoxon(df1)
savefig(joinpath(figs_loc,"Fig2/D3_SignedRank_rewardrate.pdf"))

## Cumulative pokes before leaving
Df = combine(groupby(streaks,[:Treatment,:Protocol])) do dd
    ecdf(dd,:Num_pokes)
end
Protocol_colors!(Df)
gd = groupby(Df,:Treatment)
tp = [@df subdf plot(:Xaxis, :Mean,
    xlabel = :Treatment[1],
    group = :Protocol,
    linecolor = :color,
    color= :color,
    legend = false,
    ribbon = :ERR) for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/E_pokes_per_trial.pdf"))
## Number of pokes before leaving

Df = combine(groupby(streaks,:Treatment)) do dd
    summarize(dd,:Protocol,:Num_pokes)
end
rename!(Df, :Xaxis=>:Protocol)
Protocol_colors!(Df)
gd = groupby(Df,:Treatment)
tp = [@df subdf bar(:Protocol, :Mean,
        xlabel = :Treatment[1],
        yticks = 0:2:20,
        ylims = (0,18),
        yerror = :SEM,
        color = :color,
        legend = :topleft,
        label = false,
        yaxis = "Pokes per trial") for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/F_pokes_per_trial.pdf"))

## P reward at leaving cumulative
gd = groupby(pokes,[:Treatment,:MouseID,:Day,:Protocol,:Trial])
Rrate = combine(:Pnextrew => i -> (RewRateLeaving = i[end],),gd)
gd = groupby(Rrate,[:Treatment,:Protocol])
Df = combine(gd) do dd
    ecdf(dd,:RewRateLeaving)
end
Protocol_colors!(Df)
gd = groupby(Df,:Treatment)
tp = [@df subdf plot(:Xaxis, :Mean,
    # xlabel = :Treatment[1],
    ribbon = :ERR,
    group = :Protocol,
    legend = false,
    color = :color,
    linecolor = :color,
    ylabel = "Cumulative",
    xlabel = "P reward at leaving \n $(:Treatment[1])"
    ) for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/G_pokes_per_trial.pdf"))

## P reward at leaving bar plot
gd = groupby(pokes,[:Treatment,:MouseID,:Day,:Protocol,:Trial])
Rrate = combine(:Pnextrew => i -> (RewRateLeaving = i[end],),gd)
gd = groupby(Rrate,[:Treatment])
Df = combine(gd) do dd
    summarize(dd,:Protocol,:RewRateLeaving)
end
rename!(Df, :Xaxis=>:Protocol)
Protocol_colors!(Df)
gd = groupby(Df,:Treatment)
tp = [@df subdf bar(:Protocol, :Mean,
        xlabel = :Treatment[1],
        yerror = :SEM,
        ylims = (0,0.25),
        color = :color,
        legend = :topleft,
        label = false,
        yaxis = "P reward at leaving") for subdf in gd]
ord = [1,5,4,2,7,9,8,6,3]
plot(tp[ord]...)
savefig(joinpath(figs_loc,"Fig2/H_PrewardatLeaving.pdf"))
##


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using FactCheck

facts("Testing embed_id layer") do
    net = Net(8)
    data, data_value = MemoryDataLayer(net, :data, (1,))
    embed = EmbedIDLayer(:embed, net, data, 100, 128)
    for i in 1:8
        data_value[1,i] = i
    end
    init(net)

    params = SolverParameters(
        lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
        mom_policy   = MomPolicy.Fixed(0.9),
        max_epoch    = 300,
        regu_coef    = .0005)
    sgd = SGD(params)
    forward(net; solver=sgd)
    expected = get_buffer(net, :embedweights)
    actual = get_buffer(net, :embedvalue)
    for i in 1:8
        @fact actual[:, i][:] --> roughly(expected[i, :][:])
    end
    rand!(get_buffer(net, :embed∇))
    backward(net)
    expected = get_buffer(net, :embed∇)
    actual = get_buffer(net, :embed∇weights)
    for i in 1:8
        @fact sum(actual[i, :, :], 3)[:] --> roughly(expected[:, i][:])
    end
end

FactCheck.exitstatus()


text:
for (func, typ) in [(:clblasSsymm, cl.CL_float),
                    (:clblasDsymm, cl.CL_double),
                    (:clblasCsymm, CL_float2),
                    (:clblasZsymm, CL_double2)]

    @eval @blasfun $func(order::clblasOrder, side::clblasSide,
                              uplo::clblasUplo,
                              M::Csize_t, N::Csize_t,
                              alpha::$typ,
                              A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                              B::cl.CL_mem, offb::Csize_t, ldb::Csize_t,
                              beta::$typ,
                              C::cl.CL_mem, offc::Csize_t, ldc::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, side::clblasSide,
                               uplo::clblasUplo,
                               M::Csize_t, N::Csize_t,
                               alpha::$typ,
                               A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                               B::cl.CL_mem, offb::Csize_t, ldb::Csize_t,
                               beta::$typ,
                               C::cl.CL_mem, offc::Csize_t, ldc::Csize_t)

end


text:
# affine_map.jl

"""
An affine map has the general form `y = a*x + b`, with types for `a`, `b`, `x`
and `y` such that the expression is valid.

We use matrix and vector to denote `a` and `b` respectively.
"""
abstract type AbstractAffineMap{S,T} <: AbstractMap{S,T}
end

jacobian(m::AbstractAffineMap, x) = matrix(m)

function jacobian(m::AbstractAffineMap{S,T}) where {S,T}
    A = matrix(m) * one(jac_type(S,T))
    ConstantMap{S,typeof(A)}(A)
end

islinear(map::AbstractMap) = false
islinear(map::AbstractAffineMap) = true

update_eltype(map::AbstractAffineMap{S,T}, ::Type{T}) where {S,T} = map

update_eltype(map::AbstractAffineMap{S,T}, ::Type{U}) where {S,T,U} =
    map_update_eltype(map, U)

"""
A `LinearMap` is an affine map that represents `y = a*x`, where `a` can have any
type such that `a*x` maps type `S` to type `T.`
"""
struct LinearMap{S,T,A} <: AbstractAffineMap{S,T}
    a   ::  A
end

LinearMap{T}(a) where {T} = LinearMap{T,T}(a)

LinearMap{S,T}(a) where {S,T} = LinearMap{S,T,typeof(a)}(a)

LinearMap(a::SMatrix{M,N,T}) where {M,N,T} = LinearMap{SVector{N,T},SVector{M,T},typeof(a)}(a)

LinearMap(a::Matrix{T}) where {T} = LinearMap{Vector{T},Vector{T},typeof(a)}(a)

LinearMap(a::T) where {T <: Number} = LinearMap{T,T,T}(a)

matrix(m::LinearMap) = m.a

vector(m::LinearMap) = zero(codomaintype(m))

function map_update_eltype(m::LinearMap, T)
    a = map(T, matrix(m))
    LinearMap(a)
end

(m::LinearMap)(x) = applymap(m, x)

applymap(m::LinearMap, x) = matrix(m) * x

inv(m::LinearMap{S,T}) where {S,T} = LinearMap{T,S}(inv(matrix(m)))

# Because StaticArrays does not currently support `pinv` we include a workaround:
pinv(m::SMatrix{M,N}) where {M,N}  = SMatrix{N,M}(pinv(convert(Array,m)))

left_inverse(m::LinearMap{S,T}) where {S,T} =  LinearMap{T,S}(pinv(matrix(m)))
right_inverse(m::LinearMap{S,T}) where {S,T} = LinearMap{T,S}(pinv(matrix(m)))

"""
Translation represents `y = x + v`, where `v` is a vector in the same space as
`x` and `y`.
"""
struct Translation{T} <: AbstractAffineMap{T,T}
    vector  ::  T
end

translation_map(vector::T) where T = Translation{T}(vector)

matrix(m::Translation{T}) where T = one(jac_type(T,T))

vector(m::Translation) = m.vector

function map_update_eltype(m::Translation, T)
    b = map(T, vector(m))
    Translation(b)
end

(m::Translation)(x) = applymap(m, x)

applymap(m::Translation, x) = x + vector(m)

apply_inverse(m::Translation, y) = y - vector(m)

inv(m::Translation) = Translation(-vector(m))

"""
`AffineMap` represents `y = a*x + b`, i.e. it combines a `LinearMap` and a
`Translation`.
"""
struct AffineMap{S,T,A} <: AbstractAffineMap{S,T}
    a   ::  A
    b   ::  T
end

AffineMap{S,T}(a::A, b) where {S,T,A} = AffineMap{S,T,A}(a, b)

AffineMap(a::T, b::T) where {T} = AffineMap{T,T,T}(a, b)

AffineMap(a::SMatrix{M,N,T}, b::SVector{M,T}) where {M,N,T} =
    AffineMap{SVector{N,T},SVector{M,T},typeof(a)}(a, b)

AffineMap(a::Number, b::SVector{N,T}) where {N,T} =
     AffineMap(SMatrix{N,N,T}(1.0I), b)

matrix(m::AffineMap) = m.a

vector(m::AffineMap) = m.b

function map_update_eltype(m::AffineMap, T)
    a = map(T, matrix(m))
    b = map(T, vector(m))
    AffineMap(a, b)
end

(m::AffineMap)(x) = applymap(m, x)

applymap(m::AffineMap, x) = m.a * x .+ m.b

# If y = a*x+b, then x = inv(a)*(y-b).
inv(m::AffineMap{S,T}) where {S,T} = AffineMap{T,S}(inv(m.a), -inv(m.a)*m.b)

# Todo: implement `full` for affine maps, which would result in `a` always being
# a dense matrix.

########################
# Some useful functions
########################

"Make the linear map y = a*x + b."
linear_map(a, b) = AffineMap(a, b)

"Map the interval [a,b] to the interval [c,d]."
interval_map(a, b, c, d) = linear_map((d-c)/(b-a), c - a*(d-c)/(b-a))

## Simple scaling maps

"Scale all variables by a."
scaling_map(a) = LinearMap(a)

"Scale the variables by a and b."
scaling_map(a, b) = LinearMap(SMatrix{2,2}(a, 0, 0, b))

"Scale the variables by a, b and c."
scaling_map(a, b, c) = LinearMap(SMatrix{3,3}(a, 0, 0,  0, b, 0,  0, 0,c))

"Scale the variables by a, b, c and d."
scaling_map(a, b, c, d) = LinearMap(SMatrix{4,4}(a,0,0,0, 0,b,0,0, 0,0,c,0, 0,0,0,d))

#############################
# Rotations around the origin
#############################

# Rotation in positive (counterclockwise) direction
# (Note: the SMatrix constructor expects the arguments column-first)
rotationmatrix(theta) = SMatrix{2,2}(cos(theta), sin(theta), -sin(theta), cos(theta))

# Rotation about X-axis (phi), Y-axis (theta) and Z-axis (psi)
# As above, the matrix is given column-by-column
rotationmatrix(phi,theta,psi) =
    SMatrix{3,3}(cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta),
        -cos(phi)*sin(psi)+sin(phi)*sin(theta)*cos(psi), cos(phi)*cos(psi)+sin(phi)*sin(theta)*sin(psi), sin(phi)*cos(theta),
        sin(phi)*sin(psi)+cos(phi)*sin(theta)*cos(psi), -sin(phi)*cos(psi)+cos(phi)*sin(theta)*sin(psi), cos(phi)*cos(theta))

rotation_map(theta) = LinearMap(rotationmatrix(theta))

rotation_map(phi, theta, psi) = LinearMap(rotationmatrix(phi,theta,psi))

##############
# Arithmetic
##############

∘(m2::AbstractAffineMap, m1::AbstractAffineMap) = affine_composition(m2, m1)

"""
Compute the affine map that represents map2 after map1, that is:
y = a2*(a1*x+b1)+b2 = a2*a1*x + a2*b1 + b2.
"""
affine_composition(map2::AbstractAffineMap, map1::AbstractAffineMap) =
    AffineMap(matrix(map2) * matrix(map1), matrix(map2)*vector(map1) + vector(map2))

affine_composition(map2::LinearMap{U,T}, map1::LinearMap{T,S}) where {S,T,U} =
    LinearMap{S,U}(matrix(map2) * matrix(map1))

affine_composition(map2::Translation{T}, map1::Translation{T}) where {T} =
    Translation{T}(translation_vector(map2) + translation_vector(map1))


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using TakingBroadcastSeriously: Broadcasted, unwrap

for f in :[%,
           tan, asin, acos, atan,
           sinh, cosh, tanh, asinh, acosh, atanh,
           min, max,
           hypot].args
  # copy from TakingBroadcastSeriously
  @eval Base.$f(a::Broadcasted...) = Broadcasted(broadcast_($f, unwrap.(a)...))
  @eval Base.$f(a::Broadcasted, b) = Broadcasted(broadcast_($f, unwrap(a), b))
  @eval Base.$f(b, a::Broadcasted) = Broadcasted(broadcast_($f, b, unwrap(a)))
end

for f in :[σ, sigmoid, relu, softmax, log_softmax].args
  # copy from TakingBroadcastSeriously
  @eval $f(a::Broadcasted...) = Broadcasted(broadcast_($f, unwrap.(a)...))
  @eval $f(a::Broadcasted, b) = Broadcasted(broadcast_($f, unwrap(a), b))
  @eval $f(b, a::Broadcasted) = Broadcasted(broadcast_($f, b, unwrap(a)))
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("NeidSolarScripts")   end
 using Pkg
 Pkg.activate(".")
 using HDF5, JLD2, FITSIO, FileIO, DataFrames, Query, RvSpectMLBase, Statistics, MultivariateStats, Plots

if occursin(r"NeidSolarScripts$", pwd())   cd("output56")   end
files = DataFrame(:filename=>readdir())
jld_files = files |> @filter(contains(_.filename,"new.jld2")) |> DataFrame

plt = plot()
 for day in eachrow(jld_files)
    #f = h5open("solar_20210104_new.jld2")
    f = h5open(day.filename)
    ccfs = read(f,"ccfs_espresso")
    M = fit(PCA,ccfs[401:end-401,60:end-60],maxoutdim=10,pratio=1.0)
    plot!(plt,M.mean./maximum(M.mean), label="Mean")
    map(i->plot!(plt,sign(M.proj[200,i]).*M.proj[:,i]./maximum(abs.(extrema(M.proj[:,i]))), label=string(i)),1:4);
    break
 end
 display(plt)
 #savefig("ccf_pca_20120104.png")

v_grid =read(h5open(jld_files.filename[1]),"v_grid")
v_grid = 100.0 * (-616:616)
jld_files.mean_ccf = Vector{Vector{Float64}}(undef,size(jld_files,1))
jld_files.smooth_ccf = Vector{Vector{Float64}}(undef,size(jld_files,1))
jld_files.ccf_resid = Vector{Array{Float64,2}}(undef,size(jld_files,1))
#plt = plot()
for (i,day) in enumerate(eachrow(jld_files))
     #f = h5open("solar_20210104_new.jld2")
     f = h5open(day.filename)
     #jld_files.mean_ccf[i] = read(f,"mean_ccf")
     #jld_files.smooth_ccf[i] = read(f,"ccf_template_smooth")
     jld_files.ccf_resid[i] = read(f,"ccf_resid_minus_rv_proj")
  end
  #display(plt)
  #savefig("ccf_pca_20120104.png")
mean_ccf_matrix = reduce(hcat,jld_files.mean_ccf)
smooth_ccf_matrix = reduce(hcat,jld_files.smooth_ccf)

norm_mean = vec((mean(mean_ccf_matrix[1:300,:],dims=1).+mean(mean_ccf_matrix[end-300:end,:],dims=1))/2)
norm_smooth = vec((mean(smooth_ccf_matrix[1:300,:],dims=1).+mean(smooth_ccf_matrix[end-300:end,:],dims=1))/2)
norm_mean = vec((mean(mean_ccf_matrix[600:632,:],dims=1))/2)
norm_smooth = vec((mean(smooth_ccf_matrix[600:632,:],dims=1))/2)

plot(v_grid,mean_ccf_matrix./norm_mean'.-smooth_ccf_matrix./norm_smooth')
plot(v_grid,(mean_ccf_matrix.-smooth_ccf_matrix)./norm_smooth')
xlims!(-12e3,12e3)

plot(v_grid,(mean_ccf_matrix-smooth_ccf_matrix)[:,1])

cols_use = 301:(size(mean_ccf_matrix,1)-301)
M = fit(PCA,mean_ccf_matrix[cols_use,:],maxoutdim=10,pratio=1.0)
principalvars(M)
scatter(log10.(principalvars(M)./tprincipalvar(M)));
 xlabel!("Number of PCs");
 ylabel!("log frac variance remaining")

plot(v_grid[cols_use,:],mean(mean_ccf_matrix[cols_use,:],dims=2),label="Mean");
 plot!(v_grid[cols_use,:],M.proj[:,1].+0.5,label="1");
 plot!(v_grid[cols_use,:],M.proj[:,2].+0.4,label="2");
 plot!(v_grid[cols_use,:],M.proj[:,3].+0.25,label="3");
 plot!(v_grid[cols_use,:],M.proj[:,4].+0.15,label="4");
 plot!(v_grid[cols_use,:],M.proj[:,5].+0.,label="5");
 xlims!(-1e4,1e4)
#savefig("daily_ccf_PCA.png")

f8 = h5open(jld_files.filename[8])
ccfs =read(f8,"ccfs_espresso")

plot(v_grid,ccfs[:,100:120].-mean(ccfs[:,100:120],dims=2))

plot(v_grid,ccfs[:,100:120].-mean(ccfs[:,100:120],dims=2))
xlims!(-1.2e3,1.2e3)


text:
module TestJulia15

using CompositionsBase
using Documenter: doctest
using Test

@testset "⨟" begin
    @test tuple ∘ inv === inv ⨟ tuple
    @test ⨟(tuple) === opcompose(tuple)
end

@testset "doctest" begin
    doctest(CompositionsBase; manual = false)
end

@testset "all derived from ∘" begin
    @test compose   === (∘)
    @test opcompose === (⨟)
    struct FreeMagma
        word
    end
    FM = FreeMagma
    Base.:(∘)(a::FM, b::FM) = FM((a.word, b.word))
    @test_throws MethodError opcompose()
    @test_throws MethodError compose()
    @test opcompose(FM(1))                      === FM(1)
    @test opcompose(FM(1), FM(2))               === compose(FM(2), FM(1))               === FM((2,1))
    @test opcompose(FM(1), FM(2), FM(3))        === compose(FM(3), FM(2), FM(1))        === FM(((3,2),1))
    @test opcompose(FM(1), FM(2), FM(3), FM(4)) === compose(FM(4), FM(3), FM(2), FM(1)) === FM((((4,3),2),1))

    # test that opcompose(::Vararg{<:Any, N})
    # only depends on compose(::Vararg{<:Any, N})
    # for fixed N
    struct S end
    Base.:(∘)(::S) = 1
    Base.:(∘)(::S, ::S) = 2
    Base.:(∘)(::S, ::S, ::S) = 3
    Base.:(∘)(::S, ::S, ::S, ::S) = 4
    s = S()
    @test opcompose(s,         ) === 1
    @test opcompose(s, s,      ) === 2
    @test opcompose(s, s, s,   ) === 3
    @test opcompose(s, s, s, s,) === 4
end

@testset "inference" begin
    @inferred opcompose(sin)
    @inferred opcompose(sin, cos)
    @inferred opcompose(sin, cos, tan)
    @inferred opcompose(sin, cos, tan, cot)
    @inferred opcompose(sin, cos, tan, cot, exp)
end

end  # module


text:
function value_of(::Type{TideSystem}, arg0::JString)
    return jcall(TideSystem, "valueOf", TideSystem, (JString,), arg0)
end

function values(::Type{TideSystem})
    return jcall(TideSystem, "values", Vector{TideSystem}, ())
end


text:
#=
Created on Saturday 26 December 2020
Last update: Monday 27 September 2021

@author: Bram De Jaegher
bram.de.jaegher@gmail.com
=#

# --- Types --- #

abstract type AbstractQuestion end
abstract type AbstractDifficulty end
abstract type AbstractQuestionBlock end
struct NoDiff <: AbstractDifficulty end
struct Easy <: AbstractDifficulty end
struct Intermediate <: AbstractDifficulty end
struct Hard <: AbstractDifficulty end

"""
Defines a question.

# Arguments

- `description` (md_str): a markdown string to be displayed above the validators
- `validators` (Array{Bool}): an array of booleans with the tests the answer to the question should solve.
- `status` (md_str): a markdown string used to change/display the question state (correct, missing, incorrect).
	the default value is missing and should probably not be changed.

Difficulty types
------
- `T ∈ [NoDiff, Easy, Intermediate, Hard]`
- At this point purely easthetic since the difficulty is only used as display
"""
mutable struct Question{T<:AbstractDifficulty}  <: AbstractQuestion
	description::Markdown.MD
	validators::Any
	status::Markdown.MD
	difficulty::T

	Question{T}(;description=md"",
						validators=[],
						status=md"") where {T<:AbstractDifficulty} = new{T}(description, validators, status, T())
end

"""
Defines a Question block.

# Arguments
- `title` (md_str): a markdown string to be displayed as title
- `description` (md_str): the general descriptions, this is diplayed directly below the title.
- `hints` (Array{md_str}): an array of markdown strings for hint admonitions but this can me any markdown and will be displayed in the \"hints\" section.
- `questions` (Array{AbstractQuestion}): an array of questions. Currently exactly one mandatory question is expected but 0-∞ optional questions can be defined.

# Examples
```julia
q₁ = Question(;
	description=md\"\"\"
	Complete the function `myclamp(x)` that clamps a number `x` between 0 and 1.

	Open assignments always return `missing`.
	\"\"\",
	validators= @safe[myclamp(-1)==0, myclamp(0.3)==0.3, myclamp(1.1)==1.0]
)

q₂ = Question{Easy}(;
	description=md\"\"\"
	Try to make the clamping also work for arrays.
	\"\"\",
	validators= @safe[myclamp([2.0, 0.3])==[1.0, 0.3]]
)

q₃ = Question{Intermediate}(;
	description=md\"\"\"
	This is an intermediate question. Surely you can complete this
	\"\"\",
	validators= @safe[true]
)

q₄ = Question{Hard}(;
	description=md\"\"\"
	I admit, this one is definitely harder
	\"\"\",
	validators= @safe[false]
)

qb = QuestionBlock(;
title=md"### Question 1.0: What a crazy exercise",
description=md\"\"\"
	Some additional general kind off description and all.
	Anything markdowny. Just make sure to use the triple accolades.

	\"\"\",
questions = [q₁, q₂, q₃, q₄],
hints=[	hint(md"Have you tried this?"),
		hint(md"Have you tried switching it on and off again?")]
);
```
"""
mutable struct QuestionBlock <: AbstractQuestionBlock
	title::Markdown.MD
	description::Markdown.MD
	hints::Array{Markdown.MD}
	questions::Array{T} where {T<:AbstractQuestion}

	QuestionBlock(;title=md"",
									description=md"",
									hints = Markdown.MD[],
									questions = [Question()]) = new(title, description, hints, questions)
end

# --- Rendering --- #
Base.show(io::IO, ::MIME"text/html", q::AbstractQuestionBlock) = print(io::IO, tohtml(q))

function tohtml(q::QuestionBlock)

	hint_string = ""
	if length(q.hints) > 0
		hint_string = "

Hints:
"
		for hint in q.hints
			hint_string *= "
" * html(hint) * "
"
		end
	end

	N_mandatory = sum(isa.(q.questions,Question))
	state_string = ""
	for index in 1:N_mandatory
		if q.questions[index].description !== ""
			state_string *= "
 $(html(q.questions[index].description))
"
		end
		state_string *= "
 $(html(q.questions[index].status))
"
	end

	out = """

			$(html(q.title))

 $(html(q.description))

			$state_string

			$hint_string

		$question_css
	"""
	return out
end

# --- Validation --- #
function check_answer(q::AbstractQuestion)
	validators = q.validators
	all_valid = all(validators)
	some_valid = any(validators)

	if length(validators) < 1 && return md""
	elseif ismissing(all_valid)
		status = still_missing()
	elseif some_valid && !all_valid
		status = keep_working()
	elseif !all_valid
		status = keep_working()
	elseif all_valid
		status = correct()
	end
	return status
end

function validate(q::QuestionBlock)
	for (index, question) in enumerate(q.questions)
		q.questions[index].status = check_answer(question)
	end
	return q
end

max_stars_repo_name: Beramos/PlutoTeaching

max_issues_repo_name: Beramos/PlutoTeaching

max_forks_repo_name: Beramos/PlutoTeaching


text:
module MLGoodies
using Reexport

include("utils/Utils.jl")
@reexport using .Utils

include("transforms/Transforms.jl")
@reexport using .Transforms

include("datasets/Datasets.jl")

include("batches/Batches.jl")
@reexport using .Batches

include("evaluation/Evaluation.jl")
@reexport using .Evaluation

end # module

max_stars_repo_name: dhonza/MLGoodies


text:
using Documenter, VegaDatasets

makedocs(
	modules=[VegaDatasets],
	sitename="VegaDatasets.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo="github.com/queryverse/VegaDatasets.jl.git"
)


text:
module Bridge
export ContinuousTimeProcess, SamplePath, VSamplePath, GSamplePath
export stack

export LinPro, Wiener, WienerBridge, CSpline
export sample, sample!, .., quvar, ito, bracket, lp, llikelihood, transitionprob, girsanov

export BridgeProp, DHBridgeProp, FilterProp, PBridgeProp, GuidedProp, GuidedBridge, innovations, innovations!, lptilde
export ubridge!, ubridge

export ullikelihood, ullikelihoodtrapez, uinnovations!, ubridge

# LinearProcess
export LinearProcess

# Levy
export LevyProcess, GammaProcess, GammaBridge, VarianceGammaProcess, LocalGammaProcess

# Poisson
export ThinningAlg, InhomogPoisson

# mclog
export mcstart, mcnext, mcbandmean, mcband

# euler
export SDESolver, Euler, BridgePre, EulerMaruyama, EulerMaruyama!, StochasticHeun,
    StochasticRungeKutta, bridge!, bridge, solve, solve!

# ode
export solve, solve!, solvebackward!, ODESolver, BS3

# guip
export LeftRule

using Random
using LinearAlgebra
using Statistics

include("chol.jl")

using StaticArrays
using Distributions
using Colors

import Base: rand
import Random: rand!

import Distributions: sample, sample!

function Btilde
end
function βtilde
end
function mcsvd3
end
function visualize_uncertainty
end
function B!
end
function a!
end

"""
    b!(t, y, tmp1, P)

Compute drift ``b`` in `y` (without factor ``Δt``, modifying `tmp1`.
"""
function b!
end

"""
    σ!(t, y, Δw, tmp2, P)

Compute stochastic increment at `y`, ``σ Δw``, modifying `tmp2`.
"""
function σ!
end

P(Po) = Po.Target
Pt(Po) = Po.Pt

target(Po) = Po.Target
auxiliary(Po) = Po.Pt

include("expint.jl")
#include("setcol.jl")

include("fsa.jl")
include("gaussian.jl")

include("types.jl")
include("sizedtype.jl")
include("misc.jl")

#_b((i,s)::IndexedTime, x, P) = b(s, x, P)
#_b!((i,s)::IndexedTime, x, tmp, P) = b!(s, x, tmp, P)
_b((i,s), x, P) = b(s, x, P)
_b!((i,s), x, tmp, P) = b!(s, x, tmp, P)

#btilde((i,s)::IndexedTime, x, P) = btilde(s, x, P)
#a((i,s)::IndexedTime, x, P) = a(s, x, P)
r((i,s)::IndexedTime, x, P) = r(s, x, P)
H((i,s)::IndexedTime, x, P) = H(s, x, P)

include("cspline.jl")
include("wiener.jl")
include("ellipse.jl")
include("ode.jl")
include("gode.jl")
include("lyap.jl")

include("ode!.jl")
include("diffusion.jl")
include("linpro.jl")
include("guip.jl")
include("guip!.jl")

include("partialbridge.jl")
include("partialbridgenuH.jl")
include("partialbridgen!.jl")

include("euler.jl")
include("sde.jl")

include("sde!.jl")
include("levy.jl")
include("poisson.jl")

include("timechange.jl")
include("uncertainty.jl")
include("mclog.jl")

include("bessel.jl")

include("Models.jl")

include("recipes.jl")

include("deprecated.jl")

end


text:
using Gadfly, Distributions

set_default_plot_size(6inch, 3inch)

plot(y=rand(Normal(), 100), Geom.boxplot)


text:
using Base.Test
using Meganet

npixel = 20
nex = 8
nchannel = 3
for TYPE=[Float64,Float32]
    @testset "normLayer $TYPE" begin
    for k=1:3
     L     = getNormLayer(TYPE,[npixel,nchannel,nex],k)
     testAbstractMeganetElement(L)
    end
    end

    @testset "getTVNormLayer (no scaling) $TYPE" begin
    Lb = getTVNormLayer(TYPE,[npixel,nchannel,nex],isTrainable=false)
     testAbstractMeganetElement(Lb)
    end

    @testset "getBatchNormLayer (no scaling) $TYPE" begin
    Lb = getBatchNormLayer(TYPE,[npixel,nchannel,nex],isTrainable=false)
     testAbstractMeganetElement(Lb)
    end

    @testset "getTVNormLayer (with scaling) $TYPE" begin
    Lb = getTVNormLayer(TYPE,[npixel,nchannel,nex],isTrainable=true)
     testAbstractMeganetElement(Lb,out=false)
    end

    @testset "getBatchNormLayer (with scaling) $TYPE" begin
    Lb = getBatchNormLayer(TYPE,[npixel,nchannel,nex],isTrainable=true)
     testAbstractMeganetElement(Lb)
    end
end


text:
include("init.jl")

@testset "ValueFederate Creation" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()
    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate State" begin

    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    state = h.helicsFederateGetState(vFed)
    @test state == 0

    h.helicsFederateEnterExecutingMode(vFed)

    state = h.helicsFederateGetState(vFed)
    @test state == 2

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate publication registration" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    pubid1 = h.helicsFederateRegisterTypePublication(vFed, "pub1", "string", "")
    pubid2 = h.helicsFederateRegisterGlobalTypePublication(vFed, "pub2", "int", "")
    pubid3 = h.helicsFederateRegisterTypePublication(vFed, "pub3", "double", "V")

    h.helicsFederateEnterExecutingMode(vFed)

    @test h.helicsPublicationGetKey(pubid1) == "TestA Federate/pub1"
    @test h.helicsPublicationGetKey(pubid2) == "pub2"

    @test h.helicsPublicationGetKey(pubid3) == "TestA Federate/pub3"
    @test h.helicsPublicationGetType(pubid3) == "double"
    @test h.helicsPublicationGetUnits(pubid3) == "V"

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate named point" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = "start of a longer string in place of the shorter one and now this should be very long"
    defVal = 5.3
    #testValue1 = "inside of the functional relationship of helics"
    testValue1 = "short string"
    testVal1 = 45.7823
    testValue2 = "I am a string"
    testVal2 = 0.0

    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_NAMED_POINT, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")

    h.helicsInputSetDefaultNamedPoint(subid, defaultValue, defVal)

    h.helicsFederateEnterExecutingMode(vFed)

    # publish string1 at time=0.0
    h.helicsPublicationPublishNamedPoint(pubid, testValue1, testVal1)

    @test h.helicsInputGetNamedPoint(subid) == (defaultValue, defVal)

    @test h.helicsFederateRequestTime(vFed, 1.0) == 0.01

    # # get the value
    @test h.helicsInputGetNamedPoint(subid) == (testValue1, testVal1)

    # publish a second string
    h.helicsPublicationPublishNamedPoint(pubid, testValue2, testVal2)

    # # make sure the value is still what we expect
    @test h.helicsInputGetNamedPoint(subid) == (testValue1, testVal1)
    # # make sure the string is what we expect
    # # @test value3 == testValue1
    # @test val3 == [testValue1, testVal1]

    # # advance time
    @test h.helicsFederateRequestTime(vFed, 2.0) == 0.02

    # # make sure the value was updated
    @test h.helicsInputGetNamedPoint(subid) == (testValue2, testVal2)

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)

end

@testset "ValueFederate Test Bool" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = true
    testValue1 = true
    testValue2 = false

    # register the publications
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_BOOLEAN, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")

    h.helicsInputSetDefaultBoolean(subid, defaultValue)

    h.helicsFederateEnterExecutingMode(vFed)

    # publish string1 at time=0.0
    h.helicsPublicationPublishBoolean(pubid, testValue1)
    val = h.helicsInputGetBoolean(subid)

    @test val == defaultValue

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    # get the value
    val = h.helicsInputGetBoolean(subid)

    # make sure the string is what we expect
    @test val == testValue1

    # publish a second string
    h.helicsPublicationPublishBoolean(pubid, testValue2)

    # make sure the value is still what we expect
    val = h.helicsInputGetBoolean(subid)
    @test val == testValue1
    # advance time
    grantedtime = h.helicsFederateRequestTime(vFed, 2.0)
    # make sure the value was updated
    @test grantedtime == 0.02

    val = h.helicsInputGetBoolean(subid)
    @test val == testValue2

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)

end

@testset "ValueFederate publisher registration" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    pubid1 = h.helicsFederateRegisterPublication(vFed, "pub1", h.HELICS_DATA_TYPE_STRING, "")
    pubid2 = h.helicsFederateRegisterGlobalPublication(vFed, "pub2", h.HELICS_DATA_TYPE_INT, "")
    pubid3 = h.helicsFederateRegisterPublication(vFed, "pub3", h.HELICS_DATA_TYPE_DOUBLE, "V")
    h.helicsFederateEnterExecutingMode(vFed)

    publication_key = h.helicsPublicationGetKey(pubid1)
    @test publication_key == "TestA Federate/pub1"
    publication_type = h.helicsPublicationGetType(pubid1)
    @test publication_type == "string"
    publication_key = h.helicsPublicationGetKey(pubid2)
    @test publication_key == "pub2"
    publication_key = h.helicsPublicationGetKey(pubid3)
    @test publication_key == "TestA Federate/pub3"
    publication_type = h.helicsPublicationGetType(pubid3)
    @test publication_type == "double"
    publication_units = h.helicsPublicationGetUnits(pubid3)
    @test publication_units == "V"
    publication_type = h.helicsPublicationGetType(pubid2)
    @test publication_type == "int64"

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate subscription and publication registration" begin

    broker = createBroker()
    vFed, fedinfo = createValueFederate(1, "fed0")

    pubid = h.helicsFederateRegisterPublication(vFed, "pub1", h.HELICS_DATA_TYPE_STRING, "")
    pubid2 = h.helicsFederateRegisterGlobalPublication(vFed, "pub2", h.HELICS_DATA_TYPE_INT, "volts")
    pubid3 = h.helicsFederateRegisterTypePublication(vFed, "pub3", "double", "V")

    subid1 = h.helicsFederateRegisterSubscription(vFed, "sub1", "")
    subid2 = h.helicsFederateRegisterSubscription(vFed, "sub2", "")
    subid3 = h.helicsFederateRegisterSubscription(vFed, "sub3", "V")

    h.helicsFederateEnterExecutingMode(vFed)

    publication_type = h.helicsPublicationGetType(pubid3)
    @test publication_type == "double"

    sub_key = h.helicsSubscriptionGetKey(subid1)
    @test sub_key == "sub1"
    sub_type = h.helicsInputGetType(subid1)
    @test sub_type == ""
    sub_key = h.helicsSubscriptionGetKey(subid2)
    @test sub_key == "sub2"
    sub_key = h.helicsSubscriptionGetKey(subid3)
    @test sub_key == "sub3"
    sub_type = h.helicsInputGetType(subid3)
    @test sub_type == ""
    sub_units = h.helicsInputGetUnits(subid3)
    @test sub_units == "V"
    sub_type = h.helicsInputGetType(subid2)
    @test sub_type == ""

    subid_b = h.helicsFederateGetSubscription(vFed, "sub1")
    tmp = h.helicsSubscriptionGetKey(subid_b)
    @test tmp == "sub1"
    # check the getSubscriptionByIndex function
    subid_c = h.helicsFederateGetInputByIndex(vFed, 2)
    tmp = h.helicsInputGetUnits(subid_c)
    @test tmp == "V"
    # check publications

    sv = h.helicsPublicationGetKey(pubid)
    sv2 = h.helicsPublicationGetKey(pubid2)
    @test sv == "Testfed0/pub1"
    @test sv2 == "pub2"
    pub3name = h.helicsPublicationGetKey(pubid3)
    @test pub3name == "Testfed0/pub3"

    type = h.helicsPublicationGetType(pubid3)
    @test type == "double"
    units = h.helicsPublicationGetUnits(pubid3)
    @test units == "V"

    # check the getSubscription function

    pubid_b = h.helicsFederateGetPublication(vFed, "Testfed0/pub1")
    tmp = h.helicsPublicationGetType(pubid_b)
    @test tmp == "string"
    # check the getSubscriptionByIndex function
    pubid_c = h.helicsFederateGetPublicationByIndex(vFed, 1)
    tmp = h.helicsPublicationGetUnits(pubid_c)
    @test tmp == "volts"

    #this one should be invalid
    @test_throws h.HELICSErrorInvalidArgument pubid_d = h.helicsFederateGetPublicationByIndex(vFed, 5)

    h.helicsFederateFinalize(vFed)

    state = h.helicsFederateGetState(vFed)
    @test state == h.HELICS_STATE_FINALIZE

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)

end

@testset "ValueFederate single transfer" begin

    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_STRING, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")

    h.helicsFederateEnterExecutingMode(vFed)

    h.helicsPublicationPublishString(pubid, "string1")

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    s = h.helicsInputGetString(subid)
    @test s == "string1"

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test double" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = 1.0
    testValue = 2.0
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_DOUBLE, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")
    h.helicsInputSetDefaultDouble(subid, defaultValue)

    h.helicsFederateEnterExecutingMode(vFed)

    # publish string1 at time=0.0
    h.helicsPublicationPublishDouble(pubid, testValue)

    value = h.helicsInputGetDouble(subid)
    @test value == defaultValue

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    value = h.helicsInputGetDouble(subid)
    @test value == testValue

    # publish string1 at time=0.0
    h.helicsPublicationPublishDouble(pubid, testValue + 1)

    grantedtime = h.helicsFederateRequestTime(vFed, 2.0)
    @test grantedtime == 0.02

    value = h.helicsInputGetDouble(subid)
    @test value == testValue + 1

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test complex" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    rDefaultValue = 1.0
    iDefaultValue = 1.0
    rTestValue = 2.0
    iTestValue = 2.0
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_COMPLEX, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")
    h.helicsInputSetDefaultComplex(subid, (rDefaultValue + im * iDefaultValue))

    h.helicsFederateEnterExecutingMode(vFed)

    # publish string1 at time=0.0
    h.helicsPublicationPublishComplex(pubid, (rTestValue + im * iTestValue))

    @test rDefaultValue + im * iDefaultValue == h.helicsInputGetComplex(subid)

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    @test (rTestValue + im * iTestValue) == h.helicsInputGetComplex(subid)

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test integer" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = 1
    testValue = 2
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_INT, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")
    h.helicsInputSetDefaultInteger(subid, defaultValue)

    h.helicsFederateEnterExecutingMode(vFed)

    h.helicsPublicationPublishInteger(pubid, testValue)

    value = h.helicsInputGetInteger(subid)
    @test value == defaultValue

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    value = h.helicsInputGetInteger(subid)
    @test value == testValue

    h.helicsPublicationPublishInteger(pubid, testValue + 1)

    grantedtime = h.helicsFederateRequestTime(vFed, 2.0)
    @test grantedtime == 0.02

    value = h.helicsInputGetInteger(subid)
    @test value == testValue + 1

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test string" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = "String1"
    testValue = "String2"
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_STRING, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")
    h.helicsInputSetDefaultString(subid, defaultValue)

    h.helicsFederateEnterExecutingMode(vFed)

    h.helicsPublicationPublishString(pubid, testValue)

    value = h.helicsInputGetString(subid)
    @test value == defaultValue

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    value = h.helicsInputGetString(subid)
    @test value == testValue

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test vectorD" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    defaultValue = [0.0, 1.0, 2.0]
    testValue = [3.0, 4.0, 5.0]
    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_VECTOR, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")
    h.helicsInputSetDefaultVector(subid, defaultValue)

    h.helicsFederateEnterExecutingMode(vFed)

    h.helicsPublicationPublishVector(pubid, testValue)

    @test h.helicsInputGetVector(subid) == defaultValue

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    value = h.helicsInputGetVector(subid)

    @test value == testValue

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test single transfer" begin
    broker = createBroker()
    vFed, fedinfo = createValueFederate()

    s = "n2"

    pubid = h.helicsFederateRegisterGlobalPublication(vFed, "pub1", h.HELICS_DATA_TYPE_STRING, "")
    subid = h.helicsFederateRegisterSubscription(vFed, "pub1", "")

    h.helicsFederateEnterExecutingMode(vFed)

    h.helicsPublicationPublishString(pubid, "string1")

    grantedtime = h.helicsFederateRequestTime(vFed, 1.0)
    @test grantedtime == 0.01

    s = h.helicsInputGetString(subid)

    @test s == "string1"

    time = h.helicsInputLastUpdateTime(subid)
    @test time == 0.01

    h.helicsPublicationPublishString(pubid, "string2")

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)

end

@testset "ValueFederate default value tests" begin
    broker = createBroker()
    vFed1, fedinfo = createValueFederate(1, "fed0")

    inp_raw1 = h.helicsFederateRegisterInput(vFed1, "key1", h.HELICS_DATA_TYPE_RAW, "raw")
    inp_raw2 = h.helicsFederateRegisterInput(vFed1, "key2", h.HELICS_DATA_TYPE_RAW, "raw")

    inp_bool = h.helicsFederateRegisterInput(vFed1, "key3", h.HELICS_DATA_TYPE_BOOLEAN, "")

    inp_time = h.helicsFederateRegisterInput(vFed1, "key4", h.HELICS_DATA_TYPE_TIME, "")

    inp_char = h.helicsFederateRegisterInput(vFed1, "key5", h.HELICS_DATA_TYPE_STRING, "")

    inp_vect = h.helicsFederateRegisterInput(vFed1, "key6", h.HELICS_DATA_TYPE_VECTOR, "V")

    inp_double = h.helicsFederateRegisterInput(vFed1, "key7", h.HELICS_DATA_TYPE_DOUBLE, "kW")

    inp_double2 = h.helicsFederateRegisterInput(vFed1, "key8", h.HELICS_DATA_TYPE_DOUBLE, "")

    inp_np = h.helicsFederateRegisterInput(vFed1, "key9", h.HELICS_DATA_TYPE_NAMED_POINT, "")

    h.helicsInputSetMinimumChange(inp_double, 1100.0)
    h.helicsInputSetDefaultDouble(inp_double, 10000.0)

    h.helicsInputSetOption(inp_double2, h.HELICS_HANDLE_OPTION_CONNECTION_REQUIRED, true)

    pub = h.helicsFederateRegisterPublication(vFed1, "", h.HELICS_DATA_TYPE_INT, "MW")
    h.helicsPublicationSetOption(pub, h.HELICS_HANDLE_OPTION_CONNECTION_REQUIRED, true)
    h.helicsPublicationAddTarget(pub, "Testfed0/key7")
    h.helicsPublicationAddTarget(pub, "Testfed0/key8")

    h.helicsInputSetDefaultRaw(inp_raw1, "")
    data = "this is a string"
    h.helicsInputSetDefaultRaw(inp_raw2, data)

    h.helicsInputSetDefaultBoolean(inp_bool, true)

    h.helicsInputSetDefaultTime(inp_time, 12.3)
    h.helicsInputSetDefaultChar(inp_char, 'q')
    h.helicsInputSetDefaultVector(inp_vect, Float64[])
    h.helicsInputSetDefaultNamedPoint(inp_np, data, 15.7)

    h.helicsFederateEnterExecutingMode(vFed1)
    @test h.helicsInputGetInjectionUnits(inp_double) == "MW"
    @test h.helicsInputGetInjectionUnits(inp_double2) == "MW"
    @test h.helicsInputGetType(inp_double) == "double"
    @test h.helicsInputGetPublicationType(inp_double) == "int64"

    c2 = h.helicsInputGetChar(inp_char)
    @test c2 == 'q'
    h.helicsInputGetVector(inp_vect)

    optset = h.helicsInputGetOption(inp_double2, h.HELICS_HANDLE_OPTION_CONNECTION_REQUIRED)
    @test optset == true

    optset = h.helicsPublicationGetOption(pub, h.HELICS_HANDLE_OPTION_CONNECTION_REQUIRED)
    @test optset == true
    h.helicsPublicationPublishInteger(pub, 12)

    h.helicsFederateRequestNextStep(vFed1)
    @test h.helicsInputGetDouble(inp_double) == 12000.0
    @test h.helicsInputGetDouble(inp_double2) == 12.0

    h.helicsPublicationPublishInteger(pub, 13)

    h.helicsFederateRequestNextStep(vFed1)
    @test h.helicsInputIsUpdated(inp_double) == false
    @test h.helicsInputIsUpdated(inp_double2) == true

    @test h.helicsInputGetDouble(inp_double) == 12000.0
    @test h.helicsInputGetDouble(inp_double2) == 13.0

    h.helicsPublicationPublishInteger(pub, 15)

    h.helicsFederateRequestNextStep(vFed1)

    @test h.helicsInputIsUpdated(inp_double) == true
    @test h.helicsInputIsUpdated(inp_double2) == true

    h.helicsInputClearUpdate(inp_double)
    h.helicsInputClearUpdate(inp_double2)

    @test h.helicsInputIsUpdated(inp_double) == false
    @test h.helicsInputIsUpdated(inp_double2) == false

    _, rval = h.helicsInputGetNamedPoint(inp_np)
    @test rval == 15.7

    out, rval = h.helicsInputGetNamedPoint(inp_np)
    @test out == "this is a string"
    @test rval == 15.7

    h.helicsFederateFinalize(vFed1)

    destroyFederate(vFed1, fedinfo)
    destroyBroker(broker)

end

@testset "ValueFederate test info filed" begin

    broker = createBroker()
    vFed, fedinfo = createValueFederate(1, "fed0")

    h.helicsFederateSetFlagOption(vFed, h.HELICS_HANDLE_OPTION_CONNECTION_OPTIONAL, true)
    # register the publications/subscriptions

    subid1 = h.helicsFederateRegisterSubscription(vFed, "sub1", "")
    pubid1 = h.helicsFederateRegisterTypePublication(vFed, "pub1", "string", "")
    pubid2 = h.helicsFederateRegisterGlobalTypePublication(vFed, "pub2", "string", "")

    # Set info fields
    h.helicsInputSetInfo(subid1, "sub1_test")
    h.helicsPublicationSetInfo(pubid1, "pub1_test")
    h.helicsPublicationSetInfo(pubid2, "pub2_test")
    h.helicsFederateEnterExecutingMode(vFed)

    @test h.helicsInputGetInfo(subid1) == "sub1_test"
    @test h.helicsPublicationGetInfo(pubid1) == "pub1_test"
    @test h.helicsPublicationGetInfo(pubid2) == "pub2_test"

    cr = h.helicsFederateGetCoreObject(vFed)
    h.helicsFederateFinalize(vFed)

    wait = h.helicsCoreWaitForDisconnect(cr, 70)
    if (wait == false)
        wait = h.helicsCoreWaitForDisconnect(cr, 500)
    end
    @test wait == true

    destroyFederate(vFed, fedinfo)
    destroyBroker(broker)
end

@testset "ValueFederate test file load" begin

    filename = joinpath(@__DIR__, "valuefederate.json")
    vFed = h.helicsCreateValueFederateFromConfig(filename)

    name = h.helicsFederateGetName(vFed)
    @test name == "valueFed"

    @test h.helicsFederateGetInputCount(vFed) == 3
    @test h.helicsFederateGetPublicationCount(vFed) == 2

    h.helicsFederateFinalize(vFed)
    h.helicsFederateFree(vFed)

end


text:
struct QMSHeuristic
    k::Int
    p
    η
end

function QMSHeuristic(k; p=0.0, η = 0.0)
    QMSHeuristic(k, p, η)
end

function QMSHeuristic(;p = 0.05, η = 0.0)
    QMSHeuristic(0, p=p, η = η)
end

nu_to_C(nu, n) = 1/(n*nu)

function find_parameters(method::QMSHeuristic, data, labels, oracle, solver, qm, seed)

    if method.k == 0
        k = trunc(Int, method.p * size(data, 2)) + 1 # cf. p. 5069: K = ceil(p * n)
    else
        k = method.k
    end

    kdtree = KDTree(data; leafsize = k + 1)
    nn = knn(kdtree, data, k + 1, true)

    S_K = [1/k * sum(nn[2][i][2:end]) for i in 1:size(data, 2)]

    order = sortperm(S_K)

    # Detect change point, cf. Z. Ghafoori, S. Rajasegarar, S. M. Erfani, S. Karunasekera, and C. A. Leckie, “Unsupervised Parameter Estimation for One-Class Support Vector Machines,” in Advances in Knowledge Discovery and Data Mining, 2016, pp. 183–195.
    S_K_delta = S_K[order][2:end] - S_K[order][1:end-1]
    S_K_delta2 = S_K_delta[2:end] - S_K_delta[1:end-1];
    curvature = S_K_delta2 ./ (1 .+ S_K_delta[1:end-1].^2).^0.5;

    m = findmax(curvature)[2]

    # H1.nu
    nu = (size(data, 2) - (1 - method.η) * m) / size(data, 2)
    C = nu_to_C(nu, size(data, 2))

    # H2.gamma
    δ_min, q = findmin([nn[2][i][2] for i in 1:size(data,2)])
    δ_max = 1/(size(data,2) - 1) * sum(euclidean(data[:, q], data[:, i]) for i in 1:size(data, 2) if i != q)
    γ = -log(δ_min/δ_max) / (δ_max^2 - δ_min^2)

    return (γ, C)
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
using SafeTestsets

@time begin
	@time @safetestset "Interface Tests" begin include("interface.jl") end
	@time @safetestset "Block Encryption/Decryption tests" begin include("blocktest.jl") end
	@time @safetestset "CBC Mode tests" begin include("cbc.jl") end
	@time @safetestset "CTR Mode tests" begin include("ctr.jl") end
	@time @safetestset "ECB Mode tests" begin include("ecb.jl") end
end