text:
"""
BVB Inference, general purpose.
* input
    * ll: log-likelihood, should return a float number for a vector, or a vector
          for a matrix (each column corresponds to a vector of parameter)
    * init_approx_d: a Distribution, initial approximation, usually set to be a flat
                     Normal() or MvNormal()
    * method: "inclusive" (default) or "exclusive"
    * max_iter: 100 (default)
    * n: number of MC samples used to search α, for approximating integral in Newton's method
    * verbose: >=0, more display with higher verbose
    * eval_interval: for every such # of rounds, evaluate cross_entropy and ELBO as trace information
    * eval_cross_entropy: true or false (default)
    * n_eval_cross_entropy: number of samples for evaluating cross_entropy (default: 100000)
    * eval_ELBO: true or false (default)
    * n_eval_ELBO: number of samples for evaluating cross_entropy (default: 100000)
    * show_plot: show unicode plots if true (default: false)
    * opt_grtol: gradient norm tol for optimization to terminate (default: 1e-8).
                 Can be very slow if set too small.
    * opt_max_iter: max iter allowed for optimization (default: 1000)
    * opt_n_init: pick Θ with maximum L/q from such number of samples ~ q (default: 1, just randomly initialize)
    * opt_show_trace: show trace of LBFGS if true (default: false)
    * ...: extra functions of approx_d that is called every `eval_interval`, traced and reported
* return
    * approx_d: as a MixtureModel
    * trace: a Dict storing trace information
"""
function infer(ll::Function, init_approx_d::Distribution; method="exclusive",
               max_iter=100, n=100, verbose=1, eval_interval=10,
               diagonal = false,
               eval_cross_entropy=false, eval_ELBO=false, show_plot=false,
               opt_max_iter=1000, sgd_max_iter=1000, opt_show_trace=false,
               n_eval_cross_entropy=100000, n_eval_ELBO=100000, kwargs...)
    iter = 0
    p = length(init_approx_d)
    @assert method == "inclusive" || method == "exclusive"
    println("\nInferring $p parameters with $method KL for $max_iter iterations.")
    println("Using $n MC samples for searching α.")
    trace_eval_iter = Vector{Int64}()
    trace_cross_entropy = Vector{Float64}()
    trace_ELBO = Vector{Float64}()
    trace_α = Vector{Float64}()
    trace_docall = Dict{Any, Vector}()
    if !isempty(kwargs)
        println("Running with extra reporting of $kwargs")
        for (name, do_call_f) in kwargs
            isa(do_call_f, Function) || error("Wrong keyword argument: $name is not a function.")
            trace_docall[name] = Vector{Float64}()
        end
    end

    approx_d = init_approx_d

    while iter < max_iter
        iter += 1
        println("... Iter $iter ...")
        # get new component
        h = try
                BVB.get_new_component_laplacian(ll, approx_d;
                                                verbose=(verbose>2),
                                                diagonal = diagonal,
                                                max_iter=opt_max_iter,
                                                show_trace=opt_show_trace)
            catch error_msg
                warn("Error in getting h, retrying: ", error_msg)
                continue
            end
        if show_plot && p>1
            fig = scatterplot(1:p, mean(h), color=:red, name="h", title="mean")
            if isa(approx_d, MixtureModel)
                scatterplot!(fig, 1:p, mean(approx_d.components[end]), color=:blue, name="last h")
            else
                scatterplot!(fig, 1:p, mean(approx_d), color=:blue, name="last h")
            end
            print(fig)
        end
        # determine α
        if method == "inclusive"
            α = BVB.search_α(approx_d, h, ll;
                            n=n, verbose=(verbose>0), show_plot=show_plot)
        else
            α = BVB.search_α_exclusive_sgd(approx_d, h, ll;
                                              n=n, verbose=(verbose>0), max_iter=sgd_max_iter)
        end
        # add new component
        approx_d = BVB.add_to_mixture(approx_d, h, α)
        # evaluate trace information
        push!(trace_α, α)
        if show_plot
            print(barplot(collect(1:length(trace_α)), trace_α, title="trace of α"))
        end
        if eval_interval == 1 || iter % eval_interval == 1
            push!(trace_eval_iter, iter)
            if eval_cross_entropy
                if verbose > 1
                    println("Computing cross entropy...")
                end
                push!(trace_cross_entropy, BVB.cross_entropy(approx_d, ll, n=n_eval_cross_entropy))
                if verbose > 0
                    println(@sprintf "iter %4d: cross entropy = %g" iter trace_cross_entropy[end])
                end
            end
            if eval_ELBO
                if verbose > 1
                    println("Computing ELBO...")
                end
                push!(trace_ELBO, BVB.ELBO(approx_d, ll, n=n_eval_ELBO))
                if verbose > 0
                    println(@sprintf "iter %4d: ELBO = %g" iter trace_ELBO[end])
                end
            end
            # extra do_call funtions
            for (name, do_call_f) in kwargs
                s_f = do_call_f(approx_d)
                println("Iter $iter: $name => $s_f")
                push!(trace_docall[name], s_f)
            end
            if show_plot
                _lineplot_trace(trace_docall)
            end
        end
    end

    trace = Dict{Any, Any}()
    trace["method"] = method
    trace["α"] = trace_α
    trace["iter"] = trace_eval_iter
    trace["cross_entropy"] = trace_cross_entropy
    trace["ELBO"] = trace_ELBO
    trace["max_iter"] = max_iter
    trace["n"] = n
    for name in keys(trace_docall)
        trace[name] = trace_docall[name]
    end

    return approx_d, trace
end

"""
BVB Inference, for mini-batch estimated log-likelihood.
* input
    * ll_factory: each call of ll_factory() will generate an ll(...) function, as
                  stochastic approximator of full-data ll.
    * init_approx_d: a Distribution, initial approximation, usually set to be a flat
                     Normal() or MvNormal()
    * learning_rate: a function of iter (default: t -> 1.0/t). α_t = learning_rate(t) x (α_t∗)
    * method: "inclusive" or "exclusive" (default)
    * max_iter: 100 (default)
    * n: number of MC samples used to search α, for approximating integral in Newton's method
    * first_n: number of MC samples used in the first round (default to n)
    * verbose: >=0, more display with higher verbose
    * eval_interval: for every such # of rounds, evaluate cross_entropy and ELBO as trace information
    * eval_cross_entropy: true or false (default)
    * n_eval_cross_entropy: number of samples for evaluating cross_entropy (default: 100000)
    * eval_ELBO: true or false (default)
    * n_eval_ELBO: number of samples for evaluating cross_entropy (default: 100000)
    * show_plot: show unicode plots if true (default: false)
    * opt_grtol: gradient norm tol for optimization to terminate (default: 1e-8).
                 Can be very slow if set too small.
    * opt_max_iter: max iter allowed for optimization (default: 1000)
    * opt_n_init: pick Θ with maximum L/q from such number of samples ~ q (default: 1, just randomly initialize)
    * opt_show_trace: show trace of LBFGS if true (default: false)
    * ...: extra functions of approx_d that is called every `eval_interval`, traced and reported
* return
    * approx_d: as a MixtureModel
    * trace: a Dict storing trace information
"""
# TODO(guo): make this work
function minibatch_infer(ll_factory::Function, init_approx_d::Distribution, learning_rate::Function=(t->1.0/t);
               method="exclusive",
               max_iter=100, n=50000, first_n=n, verbose=1, eval_interval=10,
               eval_cross_entropy=false, eval_ELBO=false,
               n_eval_cross_entropy=100000, n_eval_ELBO=100000,
               show_plot=false,
               opt_grtol=1e-8, opt_max_iter=1000, opt_n_init=1, opt_show_trace=false,
               kwargs...)
    iter = 0
    p = length(init_approx_d)
    @assert method == "inclusive" || method == "exclusive"
    println("\nMini-batched inferring $p parameters with $method KL for $max_iter iterations.")
    println("Using $n samples for searching α.")
    trace_eval_iter = Vector{Int64}()
    trace_cross_entropy = Vector{Float64}()
    trace_ELBO = Vector{Float64}()
    trace_α = Vector{Float64}()
    trace_α_effective = Vector{Float64}()
    trace_docall = Dict{Any, Vector}()
    if !isempty(kwargs)
        println("Running with extra reporting of $kwargs")
        for (name, do_call_f) in kwargs
            isa(do_call_f, Function) || error("Wrong keyword argument: $name is not a function.")
            trace_docall[name] = Vector{Float64}()
        end
    end

    approx_d = init_approx_d

    while iter < max_iter
        iter += 1
        println("... Iter $iter ...")
        # renew stochastic approximation of ll
        ll = ll_factory()
        # get new component
        h = try
                BVB.get_new_component_laplacian(ll, approx_d, verbose=(verbose>2),
                                                n_init=opt_n_init, grtol=opt_grtol,
                                                max_iter=opt_max_iter, show_trace=opt_show_trace)
            catch error_msg
                println("Retrying: get h", error_msg)
                iter -= 1
                continue
            end
        if show_plot && p>1
            fig = scatterplot(1:p, mean(h), color=:red, name="h", title="mean")
            if isa(approx_d, MixtureModel)
                scatterplot!(fig, 1:p, mean(approx_d.components[end]), color=:blue, name="last h")
            else
                scatterplot!(fig, 1:p, mean(approx_d), color=:blue, name="last h")
            end
            print(fig)
        end
        # determine α
        if method == "inclusive"
            α = BVB.search_α(approx_d, h, ll,
                             n=(isa(approx_d, MixtureModel) ? n : first_n),
                             verbose=(verbose>0), show_plot=show_plot)
        else
            α = BVB.search_α_exclusive_newton(approx_d, h, ll,
                                              n=(isa(approx_d, MixtureModel) ? n : first_n),
                                              verbose=(verbose>0), show_plot=show_plot)
        end
        # add new component
        lr = learning_rate(1.0 * iter)
        α_effective = α * lr
        if verbose > 0
            println("effective α <- (α = $α) x (lr = $lr) = $α_effective")
        end
        approx_d = BVB.add_to_mixture(approx_d, h, α_effective)
        # evaluate trace information
        push!(trace_α, α)
        push!(trace_α_effective, α_effective)
        if show_plot
            print(barplot(collect(1:length(trace_α)), trace_α, title="trace of α"))
            print(barplot(collect(1:length(trace_α_effective)), trace_α_effective, title="trace of effective α"))
        end
        if eval_interval==1 || iter % eval_interval == 1
            push!(trace_eval_iter, iter)
            if eval_cross_entropy
                if verbose > 1
                    println("Computing cross entropy...")
                end
                push!(trace_cross_entropy, BVB.cross_entropy(approx_d, ll, n=n_eval_cross_entropy))
                if verbose > 0
                    println(@sprintf "iter %4d: cross entropy = %g" iter trace_cross_entropy[end])
                end
            end
            if eval_ELBO
                if verbose > 1
                    println("Computing ELBO...")
                end
                push!(trace_ELBO, BVB.ELBO(approx_d, ll, n=n_eval_ELBO))
                if verbose > 0
                    println(@sprintf "iter %4d: ELBO = %g" iter trace_ELBO[end])
                end
            end
            # extra do_call funtions
            for (name, do_call_f) in kwargs
                s_f = do_call_f(approx_d)
                println("Iter $iter: $name => $s_f")
                push!(trace_docall[name], s_f)
            end
            if show_plot
                _lineplot_trace(trace_docall)
            end
        end
    end

    trace = Dict{Any, Any}()
    trace["method"] = method
    trace["α"] = trace_α
    trace["iter"] = trace_eval_iter
    trace["cross_entropy"] = trace_cross_entropy
    trace["ELBO"] = trace_ELBO
    trace["max_iter"] = max_iter
    trace["n"] = n
    for name in keys(trace_docall)
        trace[name] = trace_docall[name]
    end

    return approx_d, trace
end

function _lineplot_trace(trace::Dict)
    if !isempty(trace)
        for name in keys(trace)
            info = trace[name]
            fig = lineplot(1:length(info), info, title="extra trace", name=string(name))
            print(fig)
        end
    end
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
# test simple methods
#---------------------
withblas = TensorContractionsXD.use_blas() ? "with" : "without"
withcache = TensorContractionsXD.use_cache() ? "with" : "without"
@testset "Method syntax $withblas BLAS and $withcache cache" begin
    @testset "tensorcopy" begin
        A = randn(Float64, (3,5,4,6))
        p = randperm(4)
        C1 = permutedims(A,p)
        C2 = @inferred tensorcopy(A,(1:4...,),(p...,))
        @test C1 ≈ C2
        @test_throws TensorContractionsXD.IndexError tensorcopy(A,1:3,1:4)
        @test_throws TensorContractionsXD.IndexError tensorcopy(A,[1,2,2,4],1:4)
    end

    @testset "tensoradd" begin
        A = randn(Float64, (3,5,4,6))
        B = randn(Float64, (5,6,3,4))
        p = (3,1,4,2)
        C1 = @inferred tensoradd(A,p,B,(1:4...,))
        C2 = A+permutedims(B,p)
        @test C1 ≈ C2
        @test_throws DimensionMismatch tensoradd(A,1:4,B,1:4)
    end

    @testset "tensortrace" begin
        A = randn(Float64, (50,100,100))
        C1 = tensortrace(A,[:a,:b,:b])
        C2 = zeros(50)
        for i = 1:50
            for j = 1:100
                C2[i] += A[i,j,j]
            end
        end
        @test C1 ≈ C2
        A = randn(Float64, (3,20,5,3,20,4,5))
        C1 = @inferred tensortrace(A,(:a,:b,:c,:d,:b,:e,:c),(:e,:a,:d))
        C2 = zeros(4,3,3)
        for i1 = 1:4, i2 = 1:3, i3 = 1:3
            for j1 = 1:20,j2 = 1:5
                C2[i1,i2,i3] += A[i2,j1,j2,i3,j1,i1,j2]
            end
        end
        @test C1 ≈ C2
    end

    @testset "tensorcontract" begin
        A = randn(Float64, (3,20,5,3,4))
        B = randn(Float64, (5,6,20,3))
        C1 = @inferred tensorcontract(A,(:a,:b,:c,:d,:e),B,(:c,:f,:b,:g),(:a,:g,:e,:d,:f))
        C2 = @inferred tensorcontract(A,(:a,:b,:c,:d,:e),B,(:c,:f,:b,:g),(:a,:g,:e,:d,:f))
        C3 = zeros(3,3,4,3,6)
        for a = 1:3, b = 1:20, c = 1:5, d = 1:3, e = 1:4, f = 1:6, g = 1:3
            C3[a,g,e,d,f] += A[a,b,c,d,e]*B[c,f,b,g]
        end
        @test C1 ≈ C3
        @test C2 ≈ C3
        @test_throws TensorContractionsXD.IndexError tensorcontract(A,[:a,:b,:c,:d],B,[:c,:f,:b,:g])
        @test_throws TensorContractionsXD.IndexError tensorcontract(A,[:a,:b,:c,:a,:e],B,[:c,:f,:b,:g])
    end

    @testset "tensorproduct" begin
        A = randn(Float64, (5,5,5,5))
        B = rand(ComplexF64,(5,5,5,5))
        C1 = reshape((@inferred tensorproduct(A,(1,2,3,4),B,(5,6,7,8),(1,2,5,6,3,4,7,8))),(5*5*5*5,5*5*5*5))
        C2 = kron(reshape(B,(25,25)),reshape(A,(25,25)))
        @test C1 ≈ C2
        @test_throws TensorContractionsXD.IndexError tensorproduct(A,[:a,:b,:c,:d],B,[:d,:e,:f,:g])
        @test_throws TensorContractionsXD.IndexError tensorproduct(A,[:a,:b,:c,:d],B,[:e,:f,:g,:h],[:a,:b,:c,:d,:e,:f,:g,:i])
    end

    # test in-place methods
    #-----------------------
    # test different versions of in-place methods,
    # with changing element type and with nontrivial strides

    @testset "tensorcopy!" begin
        Abig = randn(Float64, (30,30,30,30))
        A = view(Abig,1 .+ 3*(0:9),2 .+ 2*(0:6),5 .+ 4*(0:6),4 .+ 3*(0:8))
        p = (3,1,4,2)
        Cbig = zeros(ComplexF64,(50,50,50,50))
        C = view(Cbig,13 .+ (0:6),11 .+ 4*(0:9),15 .+ 4*(0:8),4 .+ 3*(0:6))
        Acopy = tensorcopy(A,1:4,1:4)
        Ccopy = tensorcopy(C,1:4,1:4)
        TensorContractionsXD.tensorcopy!(A,1:4,C,p)
        TensorContractionsXD.tensorcopy!(Acopy,1:4,Ccopy,p)
        @test C ≈ Ccopy
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensorcopy!(A,1:3,C,p)
        @test_throws DimensionMismatch TensorContractionsXD.tensorcopy!(A,p,C,p)
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensorcopy!(A,1:4,C,[1,1,2,3])
    end

    @testset "tensoradd!" begin
        Abig = randn(Float64, (30,30,30,30))
        A = view(Abig,1 .+ 3*(0:9),2 .+ 2*(0:6),5 .+ 4*(0:6),4 .+ 3*(0:8))
        p = (3,1,4,2)
        Cbig = zeros(ComplexF64,(50,50,50,50))
        C = view(Cbig,13 .+ (0:6),11 .+ 4*(0:9),15 .+ 4*(0:8),4 .+ 3*(0:6))
        Acopy = tensorcopy(A,1:4,p)
        Ccopy = tensorcopy(C,1:4,1:4)
        α = randn(Float64)
        β = randn(Float64)
        TensorContractionsXD.tensoradd!(α,A,1:4,β,C,p)
        Ccopy = β*Ccopy+α*Acopy
        @test C ≈ Ccopy
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensoradd!(1.2,A,1:3,0.5,C,p)
        @test_throws DimensionMismatch TensorContractionsXD.tensoradd!(1.2,A,p,0.5,C,p)
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensoradd!(1.2,A,1:4,0.5,C,[1,1,2,3])
    end

    @testset "tensortrace!" begin
        Abig = rand(Float64, (30,30,30,30))
        A = view(Abig,1 .+ 3*(0:8),2 .+ 2*(0:14),5 .+ 4*(0:6),7 .+ 2*(0:8))
        Bbig = rand(ComplexF64,(50,50))
        B = view(Bbig,13 .+ (0:14),3 .+ 5*(0:6))
        Acopy = tensorcopy(A,1:4)
        Bcopy = tensorcopy(B,1:2)
        α = randn(Float64)
        β = randn(Float64)
        TensorContractionsXD.tensortrace!(α,A,[:a,:b,:c,:a],β,B,[:b,:c])
        Bcopy = β*Bcopy
        for i = 1 .+ (0:8)
            Bcopy += α*view(A,i,:,:,i)
        end
        @test B ≈ Bcopy
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensortrace!(α,A,[:a,:b,:c],β,B,[:b,:c])
        @test_throws DimensionMismatch TensorContractionsXD.tensortrace!(α,A,[:a,:b,:c,:a],β,B,[:c,:b])
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensortrace!(α,A,[:a,:b,:a,:a],β,B,[:c,:b])
        @test_throws DimensionMismatch TensorContractionsXD.tensortrace!(α,A,[:a,:b,:a,:c],β,B,[:c,:b])
    end

    @testset "tensorcontract!" begin
        Abig = rand(Float64, (30,30,30,30))
        A = view(Abig,1 .+ 3*(0:8),2 .+ 2*(0:14),5 .+ 4*(0:6),7 .+ 2*(0:8))
        Bbig = rand(ComplexF64,(50,50,50))
        B = view(Bbig,3 .+ 5*(0:6),7 .+ 2*(0:7),13 .+ (0:14))
        Cbig = rand(ComplexF32,(40,40,40))
        C = view(Cbig,3 .+ 2*(0:8),13 .+ (0:8),7 .+ 3*(0:7))
        Acopy = tensorcopy(A,1:4)
        Bcopy = tensorcopy(B,1:3)
        Ccopy = tensorcopy(C,1:3)
        α = randn(Float64)
        β = randn(Float64)
        Ccopy = β*Ccopy
        for d = 1 .+ (0:8),a = 1 .+ (0:8),e = 1 .+ (0:7)
            for b = 1 .+ (0:14),c = 1 .+ (0:6)
                Ccopy[d,a,e] += α*A[a,b,c,d]*conj(B[c,e,b])
            end
        end
        TensorContractionsXD.tensorcontract!(α,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'C',β,C,[:d,:a,:e])
        @test C ≈ Ccopy
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensorcontract!(α,A,[:a,:b,:c,:a],'N',B,[:c,:e,:b],'N',β,C,[:d,:a,:e])
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensorcontract!(α,A,[:a,:b,:c,:d],'N',B,[:c,:b],'N',β,C,[:d,:a,:e])
        @test_throws TensorContractionsXD.IndexError TensorContractionsXD.tensorcontract!(α,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'N',β,C,[:d,:e])
        @test_throws DimensionMismatch TensorContractionsXD.tensorcontract!(α,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'N',β,C,[:d,:e,:a])
    end
end

ext: jl


text:
#=move_user_resign.jl - the current user resigns, accepts 1 command line argument,<filename> => database
=#

module move_user_resign

using SQLite
function resign(database)

  #database = ARGS[1] #/path/to/database/file {string}
  db =  SQLite.DB(database) #Opens the database gamefile

  #= ---- Determines the move_number ---- =#

  res = SQLite.query(db,"""SELECT MAX("move_number") FROM moves;""") #Finds the last played move (maximum move_number)
  lastMoveIDNullable = res[1,1] #SQL query with max move_number (POSSIBLY "NULL" if no moves have been made)

  if (!isnull(lastMoveIDNullable)) #Checks that lastMoveID was not NULL
    lastMoveID = get(res[1,1]) #Parses the last move move_number as an Int

  else #lastMoveID is NULL
    lastMoveID = 0 #move_number "0" is unsused. This implies no moves have been made

  end
  move_number = lastMoveID+1 #Current move number

  #=-----UPDATE DATABASE W/MOVE-----=#
  #Option will be dropped pieces abv
    SQLite.query(db,"INSERT INTO moves (move_number, move_type)
    VALUES ($(move_number),'resign');")
end
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

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

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

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# convolution!(a, out, k);
# Array(out)
# outc = similar(img)
# copy!(outc, out)


text:
using Documenter, DataFrameTools

makedocs(
    modules = [DataFrameTools],
    sitename="DataFrameTools.jl",
    authors = "Chip Kent",
    format = Documenter.HTML(),
)

deploydocs(
    repo = "github.com/chipkent/DataFrameTools.jl.git",
    devbranch = "main",
    push_preview = true,
)


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

name = "Graphene"
version = v"1.10.0"

# Collection of sources required to build Graphene
sources = [
    ArchiveSource("https://github.com/ebassi/graphene/releases/download/$(version)/graphene-$(version).tar.xz",
                  "406d97f51dd4ca61e91f84666a00c3e976d3e667cd248b76d92fdb35ce876499")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/graphene-*/
mkdir build && cd build

meson .. \
    -Dgtk_doc=false \
    -Dgobject_types=true \
    -Dintrospection=false \
    -Dtests=false \
    -Dinstalled_tests=false \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libgraphene-1", "libgraphene-1.0"], :libgraphene),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Glib_jll",
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
    @docstring(macro_name::Symbol) -> Expr

Redefine the `docstring` variable in the innermost macro.

# Arguments
- `macro_name::Symbol`: name of the macro that generated the exception

# Returns
- `Expr`: `docstring` redefinition

# Example
```jldoctest; output = false
using PDFHighlights
using SyntaxTree

macro_name = QuoteNode(:name)

@macroexpand(PDFHighlights.Internal.Exceptions.@docstring(name)) |> linefilter! ==
quote
    macro_name = \$(macro_name)
    docstring = Meta.quot(
        \"\"\"
            \$(exception_name) <: Exception

        \$(docstring)

        \$(fields)

        See also: [`@\$(macro_name)`](@ref)
        \"\"\"
    )
end |> linefilter!

# output

true
```
"""
macro docstring(macro_name::Symbol)
    macro_name = QuoteNode(macro_name)
    return esc(
        quote
            macro_name = $(macro_name)
            docstring = Meta.quot(
                """
                    $(exception_name) <: Exception

                $(docstring)

                $(fields)

                See also: [`@$(macro_name)`](@ref)
                """
            )
        end
    )
end


text:
using DataFrames
df = DataFrame(
    :person=>["bob","phil","nick"],
    :london=>[1,1,0],
    :spain=>[1,0,0],
    Symbol("london,spain")=>[1,1,1])

function row_spread(row)
    dict = Dict{String, Int}()
    for (colname, val) in zip(keys(row), values(row))
        if colname == :person
            continue
        end
        for country in split(string(colname), ",")
            dict[country] = get(dict, country, 0) + val
        end
    end
    new_row = hcat(DataFrame(person = row.person), DataFrame(dict))
    new_row
end

reduce(vcat, row_spread(row) for row in eachrow(df))

b=a[1]

df1 = stack(df, names(df[:, Not(:person)]))

dict = Dict(n=>split(n) for n in names(df))

max_stars_repo_name: xiaodaigh/my-rl-gym

max_issues_repo_name: xiaodaigh/my-rl-gym

max_forks_repo_name: xiaodaigh/my-rl-gym


text:
#=
@reexport module BlazeCommon

import ..AbstractEchelleObservation, ..AbstractEchelleObservationSet
import ..AbstractBlazeType, ..AbstractBlaze1DType, ..AbstractBlaze2DType
import ..max_num_pixels_per_order
import ..default_pixel_buffer_lo, ..default_pixel_buffer_hi

# Util functions on data from EchelleTools
import ..num_pixels_per_order, ..num_orders, ..num_obs
=#

using Statistics

function calc_blaze(b::BT, pixel_range::AbstractRange{T2}) where {BT<:AbstractBlaze1DType, T2<:Real}
      map(p->calc_blaze(b,p), pixel_range)
end

function calc_blaze(b::BT) where {BT<:AbstractBlaze1DType}
      calc_blaze(b, b.pixel_range)
end

function calc_blaze(blaze_list::BlazeLT, data::EOT;
      output::AbstractArray{T2,2} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOT <: AbstractEchelleObservation  }
      @assert length(blaze_list.blaze_list) == num_orders(data)
      for blaze in blaze_list.blaze_list
            output[blaze.pixel_range,blaze.order] .= calc_blaze(blaze)
      end
      return output
end

function calc_blaze(blaze_list::BlazeLT, data::EOST;
      output::AbstractArray{T2,3} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet  }
      @assert length(blaze_list.blaze_list) == num_orders(data)
      for blaze in blaze_list.blaze_list
            output[blaze.pixel_range,blaze.order,:] .= calc_blaze(blaze)
      end
      return output
end
export calc_blaze

function calc_normalizations(blaze_list::BlazeLT, data::EOST;
      output::OutputT = zeros(eltype(data.flux),num_obs(data))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      # BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet, OutputT <: AbstractArray{T2,1}  }
   @assert length(blaze_list.blaze_list) == num_orders(data)
   @assert length(output) == num_obs(data)
   output .= map(obsid->
      mean(map(order->
            mean(data.flux[blaze_list.blaze_list[order].pixel_range,order,obsid] ./
                  calc_blaze(blaze_list.blaze_list[order]) )
            , 1:num_orders(data) ))
      , 1:num_obs(data))
end
export calc_normalizations

#end # module BlazeComm


text:
using Pkg
Pkg.activate(".")

# this is the code for the specific model
include("SVlib.jl")

# controls and characteristics of model
global const n = 1000
global const burnin = 100

# the following defines a function that returns a draw of the statistic, given the parameter
global const WileE_model = θ -> sqrt(n)*aux_stat(SVmodel(θ::Array{Float64,1}, n::Int64, burnin::Int64)[1])

global const θtrue = [exp(-0.736/2.0), 0.9, 0.363]
global const lb = [0.05, 0.0, 0.05]
global const ub = [2.0, 0.99, 1.0]
global const nParams = 3

# controls of raw draws from prior, for training net
global const TrainingTestingSize = Int64(1e4) # set fairly large, to limit monte carlo error
global const TrainingProportion = 0.5
global const Epochs = 1000 # passes through entire training set

println("Ready to run the SV project!");


text:
"""
    SqExponentialKernel()

The squared exponential kernel is a Mercer kernel given by the formula:
```
    κ(x, y) = exp(-‖x - y‖² / 2)
```
Can also be called via `RBFKernel`, `GaussianKernel` or `SEKernel`.
See also [`ExponentialKernel`](@ref) for a
related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
"""
struct SqExponentialKernel <: SimpleKernel end

kappa(κ::SqExponentialKernel, d²::Real) = exp(-d² / 2)

metric(::SqExponentialKernel) = SqEuclidean()

iskroncompatible(::SqExponentialKernel) = true

Base.show(io::IO,::SqExponentialKernel) = print(io,"Squared Exponential Kernel")

## Aliases ##
const RBFKernel = SqExponentialKernel
const GaussianKernel = SqExponentialKernel
const SEKernel = SqExponentialKernel

"""
    ExponentialKernel()

The exponential kernel is a Mercer kernel given by the formula:
```
    κ(x,y) = exp(-‖x-y‖)
```
"""
struct ExponentialKernel <: SimpleKernel end

kappa(κ::ExponentialKernel, d::Real) = exp(-d)

metric(::ExponentialKernel) = Euclidean()

iskroncompatible(::ExponentialKernel) = true

Base.show(io::IO, ::ExponentialKernel) = print(io, "Exponential Kernel")

## Alias ##
const LaplacianKernel = ExponentialKernel

"""
    GammaExponentialKernel(; γ = 2.0)

The γ-exponential kernel [1] is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-‖x-y‖^γ)
```
Where `γ > 0`, (the keyword `γ` can be replaced by `gamma`)
For `γ = 2`, see `SqExponentialKernel` and `γ = 1`, see `ExponentialKernel`.

[1] - Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I.
    Williams, MIT Press, 2006.
"""
struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
    γ::Vector{Tγ}
    function GammaExponentialKernel(; gamma::T=2.0, γ::T=gamma) where {T<:Real}
        @check_args(GammaExponentialKernel, γ, γ >= zero(T), "γ > 0")
        return new{T}([γ])
    end
end

@functor GammaExponentialKernel

kappa(κ::GammaExponentialKernel, d::Real) = exp(-d^first(κ.γ))

metric(::GammaExponentialKernel) = Euclidean()

iskroncompatible(::GammaExponentialKernel) = true

function Base.show(io::IO, κ::GammaExponentialKernel)
    print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
end


text:
module ArrayChannels

# Precompile imports
import Distributed
import Serialization
import Sockets
import Core

# Dependency unqualified imports
import Base: ReentrantLock
import Distributed: RRID, WorkerPool, @spawnat, @fetchfrom, @everywhere, nprocs, procs, myid, remotecall_wait, remotecall_fetch, copy!, axes
import Serialization: AbstractSerializer, serialize, deserialize, serialize_cycle_header, serialize_type, writetag, deserialize_fillarray!, AbstractSerializer, OBJECT_TAG
import Base: AbstractChannel, put!, take!, size, show, getindex, setindex!, fill!, length, broadcast, broadcast!, iterate

export
    ArrayChannel,
    ac_get_from,
    put!,
    take!,
    reduce!,
    getindex,
    setindex!,
    get_local_data

include("inplacearray.jl")
include("arraychannels.jl")
include("precompile.jl")

end


text:
"""
    parse_input_data(input_data::ImageHDU)
    parse_input_data(input_data::Tuple{AbstractArray, WCSTransform})
    parse_input_data(input_data::String, hdu_in)
    parse_input_data(input_data::FITS, hdu_in)

Parse input data and returns an Array and WCS object.

# Arguments
- `input_data`: image to reproject which can be name of a FITS file,
                an ImageHDU or a FITS file.
- `hdu_in`: used to set HDU to use when more than one HDU is present.
"""
function parse_input_data(input_data::ImageHDU)
    return read(input_data), WCS.from_header(read_header(input_data, String))[1]
end

function parse_input_data(input_data::Tuple{AbstractArray, WCSTransform})
    return input_data[1], input_data[2]
end

function parse_input_data(input_data::String, hdu_in)
    return parse_input_data(FITS(input_data), hdu_in)
end

function parse_input_data(input_data::FITS, hdu_in)
    return parse_input_data(input_data[hdu_in])
end

# TODO: extend support for passing FITSHeader when FITSHeader to WCSTransform support is possible.

"""
    parse_output_projection(output_projection::WCSTransform, shape_out)
    parse_output_projection(output_projection::ImageHDU; shape_out)
    parse_output_projection(output_projection::String, hdu_number)
    parse_output_projection(output_projection::FITS, hdu_number)

Parse output projection and returns a WCS object and shape of output.

# Arguments
- `output_projection`: WCS information about the image to be reprojected which can be
                       name of a FITS file, an ImageHDU or WCSTransform.
- `shape_out`: shape of the output image.
- `hdu_number`: specifies HDU number when file name is given as input.
"""
function parse_output_projection(output_projection::WCSTransform, shape_out)
    if length(shape_out) == 0
        throw(DomainError(shape_out, "The shape of the output image should not be an empty tuple"))
    end

    return output_projection, shape_out
end

function parse_output_projection(output_projection::ImageHDU, shape_out)
    wcs_out = WCS.from_header(read_header(output_projection, String))[1]
    if shape_out === nothing
        shape_out = size(output_projection)
    end
    if length(shape_out) == 0
        throw(DomainError(shape_out, "The shape of the output image should not be an empty tuple"))
    end
    return wcs_out, shape_out
end

function parse_output_projection(output_projection::String, hdu_number)
    parse_output_projection(FITS(output_projection), hdu_number)
end

function parse_output_projection(output_projection::FITS, hdu_number)
    wcs_out = WCS.from_header(read_header(output_projection[hdu_number], String))[1]

    if output_projection[hdu_number] isa ImageHDU
        shape_out = size(output_projection[hdu_number])
    else
        throw(ArgumentError("Given FITS file doesn't have ImageHDU"))
    end

    return wcs_out, shape_out
end


text:
"""
    plotresults(dp; <keyword arguments>)

Plot results from DPClustering object. Will plot histogram of raw data with density estimates from Gibbs sampling.
...
## Arguments
- `save = false`: Set to `true` if you want the plot to be saved
- `dir = ""`: Directory where the plot will be saved to. Default is the current working directory.
- `plotname = "DPclustering"`: Name to call plot when saving.
...
"""
function plotresults(dp; save = false, dir = "", plotname = "DPClustering")

  DFres = dp.DF
  DF = dp.data.DF

  histogram(DF[:VAF], nbins = 100, normed = true, linecolor = :white,
  color = RGBA(0.42, 0.5, 0.5, 0.6))
  plot!(DFres[:x], DFres[:mean], color = RGBA(0.325, 0.525, 0.608))
  plot!(DFres[:x], DFres[:uq], fillrange = DFres[:lq],
               fillalpha = 0.5,
               fillcolor = RGBA(0.325, 0.525, 0.608), linecolor = false,
               markerstrokecolor=:white, titlefont = font(12, "Calibri"), ytickfont = font(10, "Calibri"), xtickfont = font(10, "Calibri"), legend = false, grid = false,
               yaxis = ("Density"), xaxis = ("VAF"))

  if save == true
      Plots.savefig(joinpath(dir, "$(plotname).pdf"))
  end

end


text:
@testset "Tests on enforce_causality_and_invertibility!" begin

    # Setup
    A = [0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1];
    B = 1.25 * [0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1];
    A_companion = companion_form(A, extended=false);
    B_companion = companion_form(B, extended=false);
    A_view = @view A_companion[1:2, :];
    B_view = @view B_companion[1:2, :];
    benchmark = [0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1];

    # Test 1a: vector input does not change
    MessyTimeSeriesOptim.enforce_causality_and_invertibility!(A);
    @test A == benchmark;

    # Test 1b: vector input changes
    MessyTimeSeriesOptim.enforce_causality_and_invertibility!(B);
    @test round.(B, digits=16) == benchmark;

    # Test 2a: view input does not change
    MessyTimeSeriesOptim.enforce_causality_and_invertibility!(A_view);
    @test A_view == benchmark;

    # Test 2b: view input changes
    MessyTimeSeriesOptim.enforce_causality_and_invertibility!(B_view);
    @test round.(B_view, digits=16) == benchmark;

    # View consequences on the original companion forms
    @test round.(B_companion, digits=16) == A_companion;

    # Reset the companion forms to perform final tst
    A_companion = companion_form([0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1], extended=false);
    B_companion = companion_form(1.25 * [0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1], extended=false);
    new_benchmark = 0.7*A_companion + 0.3*B_companion;

    # Test on companion forms
    MessyTimeSeriesOptim.enforce_causality_and_invertibility!(B_companion, A_companion);
    @test B_companion == new_benchmark;
end


text:
using Adapt
using Statistics
using Oceananigans.Grids
using Oceananigans.BoundaryConditions: zero_halo_regions!

"""
    struct AveragedField{X, Y, Z, A, G, N, O} <: AbstractReducedField{X, Y, Z, A, G, N}

Type representing an average over a field-like object.
"""
struct AveragedField{X, Y, Z, S, A, G, N, O} <: AbstractReducedField{X, Y, Z, A, G, N}
       data :: A
       grid :: G
       dims :: NTuple{N, Int}
    operand :: O
     status :: S

    function AveragedField{X, Y, Z}(data, grid, dims, operand, recompute_safely) where {X, Y, Z}

        dims = validate_reduced_dims(dims)
        validate_reduced_locations(X, Y, Z, dims)
        validate_field_data(X, Y, Z, data, grid)

        status = recompute_safely ? FieldStatus(0.0) : nothing

        return new{X, Y, Z, typeof(status), typeof(data),
                   typeof(grid), length(dims), typeof(operand)}(data, grid, dims,
                                                                operand, status)
    end
end

"""
    AveragedField(operand::AbstractField; dims, data=nothing, recompute_safely=false)

Returns an AveragedField.
"""
function AveragedField(operand::AbstractField; dims, data=nothing, recompute_safely=false)

    arch = architecture(operand)
    loc = reduced_location(location(operand), dims=dims)
    grid = operand.grid

    if isnothing(data)
        data = new_data(arch, grid, loc)
        recompute_safely = true
    end

    return AveragedField{loc[1], loc[2], loc[3]}(data, grid, dims, operand, recompute_safely)
end

"""
    compute!(avg::AveragedField)

Compute the average of `avg.operand` and store the result in `avg.data`.
"""
function compute!(avg::AveragedField)
    compute!(avg.operand)

    zero_halo_regions!(avg.operand, dims=avg.dims)

    sum!(parent(avg.data), parent(avg.operand))

    sz = size(avg.grid)
    avg.data.parent ./= prod(sz[d] for d in avg.dims)

    return nothing
end

compute!(avg::AveragedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(avg, time)

#####
##### Very sugar
#####

Statistics.mean(ϕ::AbstractField; kwargs...) = AveragedField(ϕ; kwargs...)

#####
##### Adapt
#####

Adapt.adapt_structure(to, averaged_field::AveragedField{X, Y, Z}) where {X, Y, Z} =
    AveragedField{X, Y, Z}(Adapt.adapt(to, averaged_field.data),
                           Adapt.adapt(to, averaged_field.grid),
                           Adapt.adapt(to, averaged_field.dims),
                           Adapt.adapt(to, averaged_field.operand),
                           Adapt.adapt(to, averaged_field.status))


text:
using DelimitedFiles
using Printf

"""
    Image -> Array, Array -> Image
"""
function rgb_to_array(img_rgb)
    # Image -> 3D array
    # CHW -> WHC
    array = permutedims(Float32.(channelview(img_rgb)), (3, 2, 1))
    # [0, 1] -> [-1, 1]
    array = @. 2f0 * array - 1f0
    # 3D -> 4D
    return reshape(array, size(array)..., 1)
end

function save_array_as_image(path, array::AbstractArray{Float32,3})
    array = clamp.(array, -1f0, 1f0)
    img = permutedims(array, (3, 2, 1)) |> cpu
    img = @. (img + 1f0) / 2f0
    save(path, colorview(RGB, img))
end

"""
    Make directories
"""
function generate_dirs(max_step)
    dirs = ["./output/real",
            "./output/loss",
            "./output/weights",
            "./output/animation"]
    for s in 1:max_step
        push!(dirs, @sprintf("./output/fake/step%03d/adv", s))
        push!(dirs, @sprintf("./output/fake/step%03d/rec", s))
    end

    for d in dirs
        isdir(d) || mkpath(d)
    end
end

"""
    Load/Save paths
"""
scaled_real_savepath(n) = @sprintf("./output/real/step%03d.png", n)
fake_adv_savepath(n, epoch) = @sprintf("./output/fake/step%03d/adv/epoch%05d.png", n, epoch)
fake_rec_savepath(n, epoch) = @sprintf("./output/fake/step%03d/rec/epoch%05d.png", n, epoch)
discriminator_savepath(n) = @sprintf("./output/weights/dscr_step%03d.bson", n)
generator_savepath(n) = @sprintf("./output/weights/gen_step%03d.bson", n)
animation_savepath(n) = @sprintf("./output/animation/img%04d.png", n)

"""
    Save images
"""
function save_scaled_reals(real_img_p)
    # save real images
    for (n, img) in enumerate(real_img_p)
        save_array_as_image(scaled_real_savepath(n), view(img, :, :, :, 1))
    end
end

function save_generated_images(genp::GeneratorPyramid, noise_rec::AbstractVector{T},
            noise_adv::AbstractVector{T}, st::Integer, ep::Integer) where {T <: AbstractArray{Float32,4}}
    g_fake_rec = genp(noise_rec, st, false)
    save_array_as_image(fake_rec_savepath(st, ep), view(g_fake_rec, :, :, :, 1))

    g_fake_adv = genp(noise_adv, st, false)
    save_array_as_image(fake_adv_savepath(st, ep), view(g_fake_adv, :, :, :, 1))
end

"""
    Load/Save hyperparameters
"""
function save_hyperparams(path, hp::HyperParams)
    open(path, "w") do io
        JSON.print(io, show_dict(hp), 4)
    end
end

function load_hyperparams(path)
    hp = HyperParams()
    js = JSON.parsefile(path, inttype=Int64)
    foreach(js) do (k, v)
        sk = Symbol(k)
        setfield!(hp, sk, oftype(getfield(hp, sk), v))
    end
    return hp
end

"""
    Save params/result
"""
function save_training_loss(st, ep, loss_dscr, loss_gen_adv, loss_gen_rec)
    path = @sprintf("./output/loss/step%03d_loss.csv", st)
    if ep == 1
        open(path, "w") do io
            println(io, "epoch,loss_dscr,loss_gen_adv,loss_gen_rec")
        end
    end

    open(path, "a") do io
        println(io, ep, ",", loss_dscr, ",", loss_gen_adv, ",", loss_gen_rec)
    end
end

function save_noise_amplifiers(st, noise_amplifier)
    path = "./output/noise_amplifiers.csv"
    if st == 1
        open(path, "w") do io
            println(io, "noise_amplifiers")
        end
    end

    open(path, "a") do io
        println(io, noise_amplifier)
    end
end

function load_with_batchnorm!(chain::Chain, data)
    Flux.loadparams!(chain, params(data))
    for (lm, l) in zip(chain, data)
        if isa(l, BatchNorm)
            copyto!(lm.μ, l.μ)
            copyto!(lm.σ², l.σ²)
        end
    end
end

function load_model_params!(dscrp::DiscriminatorPyramid, genp::GeneratorPyramid, max_stage::Integer)
    for i in 1:max_stage
        @load discriminator_savepath(i) dscr
        load_with_batchnorm!(dscrp.chains[i], dscr |> gpu)

        @load generator_savepath(i) gen
        load_with_batchnorm!(genp.chains[i].layers, gen |> gpu)
    end
end

function save_model_params(dscrp::DiscriminatorPyramid, genp::GeneratorPyramid, st::Integer)
    dscr = dscrp.chains[st] |> cpu
    @save discriminator_savepath(st) dscr

    gen = genp.chains[st].layers |> cpu
    @save generator_savepath(st) gen
end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
__precompile__(true)
module Coord

using Reexport
@reexport using GridInterpolations
@reexport using ODE
@reexport using JLD

include("tables.jl")
include("state.jl")
include("dynamics.jl")
include("uav.jl")
include("encounter.jl")
include("strategies.jl")
include("environment.jl")
include("build_simulation.jl")
include("run_simulation.jl")

end # module


text:
# Measure the performance of the transverse Mercator projections

import Geodesy
import Proj4

function perf_geodesy(n)
    lla = Geodesy.LLA(1.0,1.0,0.0)
    trans1 = Geodesy.UTMfromLLA(Geodesy.utm_zone(lla)..., Geodesy.wgs84)
    trans2 = Geodesy.LLAfromUTM(Geodesy.utm_zone(lla)..., Geodesy.wgs84)
    for i = 1:n
        utm = Geodesy.transform(trans1, lla)
        lla2 = Geodesy.transform(trans2, utm)
        #if !(lla ≈ lla2)
    #        error("Not invertable")
#        end
    end
end

function perf_proj4(n)
    wgs84 = Proj4.Projection("+proj=longlat +datum=WGS84 +no_defs")
    utm56 = Proj4.Projection("+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs")
    lla = [150.0, -27.0, 0.0]
    for i = 1:n
        utm = Proj4.transform(wgs84, utm56, lla)
        lla2 = Proj4.transform(utm56, wgs84, utm)
#        if !(lla ≈ lla2)
#            error("Not invertable")
#        end
    end
end

println("Comparing LLA->UTM->LLA transformation speed for Geodesy and Proj4")
println("\nwarming up...")
perf_geodesy(1)
@time perf_geodesy(1)
perf_proj4(1)
@time perf_proj4(1)

n = 1_000_000
println("\nPerforming (2×) $n transformations with Geodesy")
perf_geodesy(1)
@time perf_geodesy(n)
println("\nPerforming (2×) $n transformations with Proj4")
perf_proj4(1)
@time perf_proj4(n)

# 5/5/2016 1:40pm (~14 × slower) (probably more like 25 × slower excluding ≈ test)
#
# Comparing LLA->UTM->LLA transformation speed for Geodesy and Proj4
#
# warming up...
#   0.000042 seconds (469 allocations: 19.448 KB)
#   0.000026 seconds (18 allocations: 736 bytes)
#
# Performing (2×) 1000000 transformations with Geodesy
#  16.281498 seconds (315.00 M allocations: 8.717 GB, 4.76% gc time)
#
# Performing (2×) 1000000 transformations with Proj4
#   1.126085 seconds (12.00 M allocations: 488.281 MB, 2.59% gc time)
#
# ------------------------------------------------------------------------------
# 12/5/2016 1:35pm (~2.4 × slower, ran excluding ≈ test)
#
# Comparing LLA->UTM->LLA transformation speed for Geodesy and Proj4
#
# warming up...
#   0.000067 seconds (856 allocations: 32.526 KB)
#   0.000022 seconds (9 allocations: 512 bytes)
#
# Performing (2×) 1000000 transformations with Geodesy
#   1.445283 seconds (712 allocations: 22.531 KB)
#
# Performing (2×) 1000000 transformations with Proj4
#   0.608355 seconds (3.00 M allocations: 274.658 MB, 1.71% gc time)


text:
using Test
using Dates
using ExtensibleScheduler

@testset "TriggerIterator" begin
    dt_start = DateTime(2010, 1, 1, 20, 30)

    @testset "InfiniteTrigger FiniteIterator" begin
        n_trig = -1
        n_itr = 4
        trigger = Trigger(Dates.Time(20, 30), n=n_trig)
        itr = iterate(trigger, dt_start; n=n_itr)

        @test_throws MethodError length(trigger)
        @test length(itr) == n_itr

        n_triggered = 0
        for dt in itr
            #println(dt)
            n_triggered += 1
        end
        @test n_triggered == 4
    end
    @testset "FiniteTrigger FiniteIterator" begin
        n_trig = 3
        n_itr = 4
        trigger = Trigger(Dates.Time(20, 30), n=n_trig)
        itr = iterate(trigger, dt_start; n=n_itr)

        @test length(trigger) == n_trig
        @test length(itr) == min(n_trig, n_itr)

        n_triggered = 0
        for dt in itr
            n_triggered += 1
        end
        @test n_triggered == min(n_trig, n_itr)
    end

    @testset "FiniteTrigger InfiniteIterator" begin
        n_trig = 3
        n_itr = -1
        trigger = Trigger(Dates.Time(20, 30), n=n_trig)
        dt_start = DateTime(2010, 1, 1, 20, 30)
        itr = iterate(trigger, dt_start; n=n_itr)

        @test length(trigger) == n_trig
        @test length(itr) == n_trig

        n_triggered = 0
        for dt in itr
            println(dt)
            n_triggered += 1
        end
        @test n_triggered == n_trig
    end

    @testset "InfiniteTrigger InfiniteIterator" begin
        n_trig = -1
        n_itr = -1
        trigger = Trigger(Dates.Time(20, 30), n=n_trig)
        dt_start = DateTime(2010, 1, 1, 20, 30)
        itr = iterate(trigger, dt_start; n=n_itr)

        @test_throws MethodError length(trigger)
        @test_throws MethodError length(itr)

        n_stop = 10
        n_triggered = 0
        for (i, dt) in enumerate(itr)
            n_triggered += 1
            if i == n_stop
                break
            end
        end
        @test n_triggered == n_stop
    end

end


text:
using StochasticProcesses

bm = BrownianMotion()
t = linspace(0, 1, 1000)

@time cumsim(bm, t, 1000)
Profile.clear()
@profile @time cumsim(bm, t, 1000)

f = open(perf/"cumsim_prof.prof","w")
withenv("COLUMNS"=>"10000") do
  Profile.print(f)
end
close(f)


text:
################################
## definition of Timenum type ##
################################

## constraints:
## - DataArray, where individual entries are subtype of Date
## - idx and vals must have same number of observations
## - values must be numeric

type Timenum{T} <: AbstractTimenum
    vals::DataFrame
    idx::Array{T, 1}

    function Timenum(vals::DataFrame, idx::Array{T, 1})
        chkIdx(idx)
        chkNumDf(vals)
        if(size(vals, 1) != length(idx))
            if (length(idx) == 0) | (size(vals, 1) == 0)
                df = convert(DataFrame, rand(2, 2))
                return new(df[1:2, [false, false]], Array{T, 1}[])
            end
            error(length(idx), " idx entries, but ", size(vals, 1), " rows of data")
        end
        return new(vals, idx)
    end
end

function Timenum{T}(vals::DataFrame, idx::Array{T, 1})
    return Timenum{T}(vals, idx)
end


text:
module Boltz

using Lux
using NNlib
using Random
using Statistics

# Loading Pretained Weights
using Artifacts, LazyArtifacts
using JLD2

# TODO: We want to have generic Lux implementaions for Metalhead models
# We can automatically convert several Metalhead.jl models to Lux
using Metalhead

# Mark certain parts of layers as non-differentiable
import ChainRulesCore

# Utility Functions
include("utils.jl")

# General Layer Implementations
include("layers.jl")

# Vision Models
## Vision Transformer Implementation
include("vision/vit.jl")
## VGG Nets
include("vision/vgg.jl")
## Metalhead to Lux
include("vision/metalhead.jl")

# Exports
export alexnet, convmixer, densenet, googlenet, mobilenet, resnet, resnext, vgg,
       vision_transformer

end


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables,Distributions
include("lgwt.jl")
#data
function PS4()
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code;

#Question 1
#from his solution:
function mlogit_with_Z(theta, X, Z, y)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

#values from PS3
startvals=[0.05570767876416688, 0.08342649976722213, -2.344887681361976, 0.04500076157943125, 0.7365771540890512, -3.153244238810631, 0.09264606406280998, -0.08417701777996893, -4.273280002738097, 0.023903455659102114, 0.7230648923377259, -3.749393470343111, 0.03608733246865346, -0.6437658344513095, -4.2796847340030375, 0.0853109465190059, -1.1714299392376775, -6.678677013966667, 0.086620198654063, -0.7978777029320784, -4.969132023685069, -0.0941942241795243];
#look at this to see why you might need this: http://julianlsolvers.github.io/Optim.jl/v0.9.3/algo/autodiff/
td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward);

theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))

#Hessian from ps3
theta_hat_mle_ad = theta_hat_optim_ad.minimizer
println(theta_hat_mle_ad)
H  = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([theta_hat_mle_ad theta_hat_mle_ad_se])

#question 2
# now it is positive which makes more sense because increase in log wage should increase utility

#question 3
#part a

#part b
#1
nodes, weights = lgwt(7,-10,10);
d=Normal(0,2)
sum(weights.*pdf.(d,nodes))
println(sum((nodes.^2).*weights.*pdf.(d,nodes)))

#2
nodes, weights = lgwt(10,-10,10);
sum(weights.*pdf.(d,nodes))
println(sum((nodes.^2).*weights.*pdf.(d,nodes)))

#3
# it is close, with more quadrature points it will converge to 4

#Part C
#1
d=Normal(0,2)
N=1000000
x=rand(Uniform(-10, 10 ),N)

println((20/N)*sum(x.^2 .*(pdf.(d,x))))

#2
println((20/N)*sum(x .*(pdf.(d,x))))

#3
println((20/N)*sum(pdf.(d,x)))

#4
#the difference is around 0.2 which is high in my opinion

#Question 4
function mxlogit_with_quadrature(theta, X, Z, y,rn)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        nodes, weights = lgwt(rn,-4,4);
        d=Normal(0,1)

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = sum(exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma).*weights.*pdf.(d,nodes))
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( log.(P.^bigY) )

        return loglike
    end

#Question 5
function mxlogit_with_simulation(theta, X, Z, y,rn)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        d=Normal(0,1)
        rx=rand(Uniform(0,1 ),rn)

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = sum(exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma).*(pdf.(d,x)))
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( log.(P.^bigY) )

        return loglike
    end

end

PS4()


text:
using Test, Documenter, PkgTPL
DocMeta.setdocmeta!(PkgTPL, :DocTestSetup, :(using PkgTPL), recursive = true)

@testset "PkgTPL" begin
    doctest(PkgTPL)
end


text:
@testset "Testing Example 1" begin
    local m
    makemap = true
    placement = true
    routing = true
    # Try making a map
    m = Example1.make_map()

    # Try placement - use the Search Move Generator to give it a go.
    place!(m; move_attempts = 5000)

    # Try routing
    route!(m)

    # Get statistics from the map
    MapperCore.report_routing_stats(m)

    expected_links = 12
    hist = MapperCore.global_link_histogram(m)
    # Get the number of global links from this
    found_links = 0
    for (k, v) in hist
        found_links += k * v
    end
    @test found_links == expected_links
    @test MapperCore.check_routing(m)
end

@testset "Testing Example 2" begin
    # Try making a map
    m = Example2.make_map()

    # Try placement
    place!(m)

    # Try routing
    route!(m)

    # Get statistics from the map
    MapperCore.report_routing_stats(m)
    # expect routing to fail because of inadequate resources.
    @test MapperCore.check_routing(m) == false
end


text:
ex = quote begin
    function foo(x,z)
        y = x+2
        return y,y/z
        end
    end
    println(foo(3, 1))
    println(foo(3.0, 1))
end
ex′ = postwalk(annotate, ex)
ex2 = wrap(ex′)

calls = Edges(eval(ex2))
@test length(calls) == 2

edgelist = @typegraph begin
    function foo(x)
        return x/2
    end
    function bar(x)
        return 3x+1
    end

    function g(y)
        if y < 2
            return 1
        end
        if mod(y,2) == 0
            return g(foo(y))
        else
            return g(bar(y))
        end
        return 0
    end
    #g(4)
    println("====")
    g(5)
end

#TODO: figure out why unique(edgelist doesn't work, I think it is == of arrays of types or something.
@test length(unique([(f.func, f.args, f.ret) for f in  Edges(edgelist)])) == 5


text:
## `fmpz` to GAP integer
function GAP.julia_to_gap(obj::fmpz)
  Nemo._fmpz_is_small(obj) && return GAP.julia_to_gap(Int(obj))
  GC.@preserve obj begin
    x = Nemo._as_bigint(obj)
    return ccall((:MakeObjInt, GAP.libgap), GAP.GapObj, (Ptr{UInt64}, Cint), x.d, x.size)
  end
end

# Basically the same as the usual image function but without a type check since
# we don't have elem_type(C) in this case
function image(M::MapFromFunc{D, C}, a) where {D, C <: GAP.GapObj}
  if isdefined(M, :header)
    if isdefined(M.header, :image)
      return M.header.image(a)
    else
      error("No image function known")
    end
  else
    return M(a)
  end
end

function preimage(M::MapFromFunc{D, C}, a) where {D, C <: GAP.GapObj}
  if isdefined(M, :header)
    if isdefined(M.header, :preimage)
      return M.header.preimage(a)
    else
      error("No preimage function known")
    end
  else
    return M\(a)
  end
end

################################################################################
#
#  Ring isomorphism
#
################################################################################

# computes the isomorphism between the Oscar field F and the corresponding GAP field
function _ring_iso_oscar_gap(F::T) where T <: Union{Nemo.GaloisField, Nemo.GaloisFmpzField}
   p = characteristic(F)
   G = GAP.Globals.GF(GAP.Obj(p))
   e = GAP.Globals.One(G)

   f(x::Union{Nemo.gfp_elem, Nemo.gfp_fmpz_elem}) = GAP.Obj(lift(x))*e
   finv(x) = F(fmpz(GAP.Globals.IntFFE(x)))

   return MapFromFunc(f, finv, F, G)
end

function _ring_iso_oscar_gap(F::T) where T <: Union{Nemo.FqNmodFiniteField, Nemo.FqFiniteField}
   p = characteristic(F)
   d = degree(F)
   p_gap = GAP.Obj(p)
   Fp_gap = GAP.Globals.GF(p_gap) # the prime field in GAP
   e = GAP.Globals.One(Fp_gap)

   # prime fields are easy and efficient to deal with, handle them separately
   if d == 1
      f1(x::Union{Nemo.fq_nmod, Nemo.fq}) = GAP.Obj(coeff(x, 0))*e
      finv1(x) = F(fmpz(GAP.Globals.IntFFE(x)))
      return MapFromFunc(f1, finv1, F, Fp_gap)
   end

   # non-prime fields: convert the defining polynomial to GAP...
   L = [ GAP.Obj(lift(coeff(defining_polynomial(F), i)))*e for i in 0:d - 1 ]
   push!(L, e)
   L_gap = GAP.GapObj(L)
   f_gap = GAP.Globals.UnivariatePolynomial(Fp_gap, L_gap)

   # ... and compute a GAP field G defined via this polynomial
   # (If the given polynomial is a Conway polynomial then we may call
   # GAP's `GF(p, d)`, which avoids GAP's `AlgebraicExtension`.)
   if GAP.Globals.IsCheapConwayPolynomial(p_gap, d) &&
      f_gap == GAP.Globals.ConwayPolynomial(p_gap, d)
      G = GAP.Globals.GF(p_gap, d)
   else
      G = GAP.Globals.GF(Fp_gap, f_gap)
   end

   # compute matching bases of both fields
   if GAP.Globals.IsAlgebraicExtension(G)
#FIXME:
# As soon as the problem from https://github.com/gap-system/gap/issues/4694
# is fixed, go back to `Basis_G = GAP.Globals.Basis(G)` also in this case.
# Note that the above `GF` call delegates to `FieldExtension`,
# and we want a basis in the filter `IsCanonicalBasisAlgebraicExtension`.
      Basis_G = GAP.Globals.Objectify(GAP.Globals.NewType(
                GAP.Globals.FamilyObj(G),
                GAP.Globals.IsCanonicalBasisAlgebraicExtension),
                GAP.NewPrecord(0))
      GAP.Globals.SetUnderlyingLeftModule(Basis_G, G)
   else
      Basis_G = GAP.Globals.Basis(G)
   end
   Basis_F = Vector{elem_type(F)}(undef, d)
   Basis_F[1] = F(1)
   for i = 2:d
      Basis_F[i] = Basis_F[i - 1]*gen(F)
   end

   function f(x::Union{fq_nmod, fq})
      v = [ GAP.Obj(coeff(x, i)) for i in 0:d - 1 ]
      return sum([ v[i]*Basis_G[i] for i in 1:d ])
   end

   # For "small" primes x should be an FFE, but for bigger ones it's GAP_jll.Mptr (?)
   function finv(x)
      v = GAP.Globals.Coefficients(Basis_G, x)
      v_int = [ fmpz(GAP.Globals.IntFFE(v[i])) for i = 1:length(v) ]
      return sum([ v_int[i]*Basis_F[i] for i = 1:d ])
   end

   return MapFromFunc(f, finv, F, G)
end

function __to_gap(h, x::Vector)
  return GAP.Globals.GModuleByMats(GAP.julia_to_gap([GAP.julia_to_gap(map(h, Matrix(y))) for y in x]), codomain(h))
end

function __gap_matrix_to_julia(h, g)
  return matrix(domain(h), [map(y -> preimage(h, y), gg) for gg in GAP.gap_to_julia(g)])
end

function __to_julia(h, C)
  return [ matrix(domain(h), [map(y -> preimage(h, y), gg) for gg in GAP.gap_to_julia(g)]) for g in GAP.Globals.MTX.Generators(C)]
end

function Hecke.stub_composition_factors(x::Vector{T}) where {T <: MatElem}
  F = base_ring(x[1])
  h = _ring_iso_oscar_gap(F)
  V = __to_gap(h, x)
  Vcf = GAP.Globals.MTX.CompositionFactors(V)
  res = Vector{T}[]
  for C in Vcf
    push!(res, _to_julia(h, C))
  end
  return res
end

function Hecke.stub_basis_hom_space(x::Vector{T}, y::Vector{T}) where {T <: MatElem}
  F = base_ring(x[1])
  h = _ring_iso_oscar_gap(F)
  @assert base_ring(x[1]) == base_ring(y[1])
  @assert length(x) == length(y)
  hb = GAP.Globals.MTX.BasisModuleHomomorphisms(__to_gap(h, x), __to_gap(h, y))
  hbb = [__gap_matrix_to_julia(h, g) for g in GAP.gap_to_julia(hb)]
  return hbb
end


text:
using Flux
using Flux: @epochs
using Flux: params
using Flux: throttle
using Flux: glorot_uniform
using Flux: @functor
using Plots
using Distributions
using Zygote
using Flux.Optimise: update!
using ForwardDiff
using ReverseDiff
using FiniteDiff
using ZygoteRules

```
Workflow:

# 1. Define network architecture:
    # 1. define 1 LSTM-like layer
    # 2. define 1 dense layer with activation for the start
    # 3. define 1 dense layer without activation for the end
    # 4. connect them all together -> model
# 2. Define the finite element solution function
# 3. Define sampling function
# 4. Define the loss function
# 5. Train the model

```

# define a DGM layer
mutable struct LSTMLayer
    Ug :: Array{Float64}
    Uz :: Array{Float64}
    Ur :: Array{Float64}
    Uh :: Array{Float64}
    Wg :: Array{Float64}
    Wz :: Array{Float64}
    Wr :: Array{Float64}
    Wh :: Array{Float64}
    bg :: Array{Float64}
    bz :: Array{Float64}
    br :: Array{Float64}
    bh :: Array{Float64}
    act1
    act2
end

LSTMLayer(input_dim :: Integer, output_dim :: Integer, activation_1, activation_2) = LSTMLayer(glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, input_dim),
    glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, input_dim), glorot_uniform(output_dim, output_dim), glorot_uniform(output_dim, output_dim),
    glorot_uniform(output_dim, output_dim), glorot_uniform(output_dim, output_dim), zeros(output_dim), zeros(output_dim), zeros(output_dim), zeros(output_dim),
    activation_1, activation_2)

function (a :: LSTMLayer)(S,x)

    # assign all opbject variables
    Ug, Uz, Ur, Uh, Wg, Wz, Wr, Wh, bg, bz, br, bh, act1, act2 = a.Ug, a.Uz, a.Ur, a.Uh, a.Wg, a.Wz, a.Wr, a.Wh, a.bg, a.bz, a.br, a.bh, a.act1, a.act2

    # define all operations
    G = act1.(Ug*x + Wg*S .+ bg)
    Z = act1.(Uz*x + Wz*S .+ bz)
    R = act1.(Ur*x + Wr*S .+ br)
    H = act2.(Uh*x + Wh*(S.*R) .+ bh)
    S_new = (1 .- G).*H + Z.*S

    return S_new
end

# define dense layer with activation
mutable struct Dense_start
    W1 :: Array{Float64}
    b1 :: Array{Float64}
    act
end

Dense_start(input_dim::Integer, output_dim::Integer, act) = Dense_start(glorot_uniform(output_dim, input_dim), zeros(output_dim), act)

function (a :: Dense_start)(x)

    # assign all opbject variables
    W1, b1, act = a.W1, a.b1, a.act

    # define all operations
    S_1 = act.(W1*x .+ b1)

    return S_1
end

# define dense layer without activation
mutable struct Dense_end
    W :: Array{Float64}
    b :: Array{Float64}
end

Dense_end(input_dim::Integer, output_dim::Integer) = Dense_end(glorot_uniform(output_dim, input_dim), zeros(output_dim))

function (a :: Dense_end)(x)

    # assign all opbject variables
    W, b = a.W, a.b

    # define all operations
    y = W*x .+ b

    return y
end

# connect them together into one model
mutable struct DGM
    D1 :: Dense_start
    LSTMLayers :: Array{LSTMLayer}
    D2 :: Dense_end
end

DGM(layer_width, n_layers, spatial_dim) = DGM(Dense_start(spatial_dim + 1, layer_width, tanh),
[LSTMLayer(spatial_dim + 1, layer_width, tanh, tanh) for i in 1:n_layers], Dense_end(layer_width, 1))

function (a :: DGM)(t, x)

    # assign object variables
    D1, LSTMLayers, D2 = a.D1, a.LSTMLayers, a.D2

    #concatenate time and space inputs
    X = vcat(t, x)

    # define all operations
    S = D1(X)

    for lstm_layer in LSTMLayers
        S = lstm_layer(S, X)
    end

    y = D2(S)

    return y
end

# convey trainable parameters to Flux
@functor Dense_start
@functor LSTMLayer
@functor Dense_end
@functor DGM

# check if Flux recognizes parameters
m1 = DGM(50, 3, 1)
params(m1)

# define analytical solution function
function european_put(S, K, r, sigma, t)

   d1 = (log.(S./K) .+ (r + sigma^2/2)*(t_term-t))/(sigma*sqrt(t_term-t))
   d2 = d1 .- (sigma*sqrt(t_term-t))
   put_price = -S.*cdf.(Normal(0,1), -d1) .+ K*exp(-r * (t_term-t))*cdf.(Normal(0,1), -d2)

   return put_price
end

# define the fintie element solution for American put
function finite_element_solution()
    dt = t_term/n_steps
    t = collect(0:dt:t_term + dt)
    dx = sigma*sqrt(3*dt)
    alpha = 0.5*(sigma^2)*dt/(dx^2)
    beta = (r - 0.5*sigma^2)*dt/(2*dx)

    # diagnostics
    # println("r: $(r)")
    # println("sigma: $(sigma)")
    #
    # println("dt: $(dt)")
    # # println("t: $(t)")
    # println("dx: $(dx)")
    # println("alpha: $(alpha)")
    # println("beta: $(beta)")

    # log space grid
    x_max = ceil(10*sigma*sqrt(t_term)/dx)*dx
    x_min = -x_max
    x = collect(x_min:dx:x_max + dx)
    Ndx = length(x) - 1
    x_int = collect(2:Ndx)

    # diagnostics
    # println("x_max: $(x_max)")
    # println("x_min: $(x_min)")
    # # println("x: $(x)")
    # println("len(x): $(length(x))")
    # println("Ndx: $(Ndx)")
    # # println("x_int: $(x_int)")
    # println("len(x_int): $(length(x_int))")

    # put option price grid
    h = ones(Ndx + 1, n_steps + 1)
    h[:, end] = max.(K .- S0*exp.(x), 0)
    hc = ones(Ndx + 1)
    S = S0.*exp.(x)

    # # println("h: $(x)")
    # println("len(h): $(length(h))")
    # # println("hc: $(x)")
    # println("len(hc): $(length(hc))")
    # # println("S: $(x)")
    # println("len(S): $(length(S))")

    # calculate optimal exercise boundary
    exer_bd = ones(n_steps + 1)
    exer_bd[end] = K

    # back-step in time
    for n in collect(n_steps + 1:-1:2)
        # calculate continuation values
        hc[x_int] =  h[x_int, n] .- r*dt*h[x_int, n] .+ beta.*(h[x_int .+ 1, n] .- h[x_int .- 1, n]) + alpha*(h[x_int .+ 1, n] -2*h[x_int, n] + h[x_int .- 1,n])
        hc[1] = 2*hc[2] - hc[3]
        hc[Ndx + 1] = 2*hc[Ndx] - hc[Ndx]

        # compare with intrinsic values
        exer_val = K .- S
        h[:, n-1] = max.(hc, exer_val)

        # calculate optimal exercise boundaries
        checkIdx = (hc > exer_val)*1
        idx = argmax(checkIdx) - 1
        if max(checkIdx) > 0
            exer_bd[n-1] = S[idx]
        end
    end
    putPrice = h
    return t, S, putPrice, exer_bd
end

function sampling_function(space_internal, space_terminal)

    # within-domain sample
    t_internal = rand(Uniform(t_initial, t_term), space_internal)
    S_internal = rand(Uniform(S_low, S_high*S_multiplier), space_internal)

    # terminal sample
    t_terminal = t_term .*ones(space_terminal)
    S_terminal = rand(Uniform(S_low, S_high*S_multiplier), space_terminal)

    return t_internal', S_internal', t_terminal', S_terminal'
end

# define the loss function
function loss_function(N)

    t_internal, S_internal, t_terminal, S_terminal = sampling_function(nSim_interior, nSim_terminal)

    # differential operator loss
    ϵ = 0.01
    model_output_internal = model(t_internal, S_internal)
    ∂g∂x = (model(t_internal, S_internal .+ ϵ) - model_output_internal)./ϵ
    ∂g∂t = (model(t_internal .+ ϵ, S_internal) - model_output_internal)./ϵ
    ∂g∂xx = (model(t_internal, S_internal .+ 2*ϵ) - 2*model(t_internal, S_internal .+ ϵ) + model_output_internal)./(ϵ^2)

    operator_loss_vec = ∂g∂t + r.*S_internal.*∂g∂x + (0.5*(sigma^2)).*(S_internal.^2).*∂g∂xx - r.*model_output_internal

    payoff = relu.(K .- S_internal)
    value = model(t_internal, S_internal)
    L1 = mean((operator_loss_vec.*(value-payoff)).^2)

    temp = relu.(operator_loss_vec)
    L2 = mean(temp.^2)

    V_ineq = relu.(-(payoff - value))
    L3 = mean(V_ineq.^2)

    target_payoff = relu.(K .- S_terminal)
    fitted_payoff = model(t_terminal, S_terminal)

    L4 = mean((fitted_payoff - target_payoff).^2)

    return L1 + L2 + L3 + L4
end

# initialize all problem parameters
r = 0.05           # Interest rate
sigma = 0.5       # Volatility
K = 50             # Strike
t_term = 1              # Terminal time
S0 = 50           # Initial price

# Solution parameters (domain on which to solve PDE)
t_initial = 0 + 1e-10    # time lower bound
S_low = 0.0 + 1e-10  # spot price lower bound
S_high = 2*K         # spot price upper bound

# neural network parameters
n_steps = 10000
num_layers = 3
nodes_per_layer = 50
learning_rate = 0.001

# Training parameters
sampling_stages  = 100   # number of times to resample new time-space domain points
steps_per_sample = 10    # number of SGD steps to take before re-sampling

# Sampling parameters
nSim_interior = 1000
nSim_terminal = 100
S_multiplier  = 1.5   # multiplier for oversampling i.e. draw S from [S_low, S_high * S_multiplier]

# define model
model = DGM(nodes_per_layer, num_layers, 1)

S_model = reshape(S_values, 1, length(S_values))
t_model = 0.6667.*ones(1, length(S_values))

function train_model()
    # set optimizer as ADAM
    evalcb = function ()
    @show(loss_function())
    end
    opt = ADAM(learning_rate)
    dataset = [(1) for i in 1:10]
    Flux.train!(loss_function, params(model), zip(dataset), opt, cb = throttle(evalcb, 25))
end

@epochs 1000 train_model()

params(model)

# compare with analytical solution
compare_times = [0, 0.3333, 0.6667, 1]
n_plots = 100
S_values = [S_low + i*(S_high - S_low)/n_plots for i in 0:n_plots]
S_model = reshape(S_values, 1, length(S_values))

t, S, price, exer_bd = finite_element_solution()
S_idx = [i for i in 1:length(S) if(S[i] < S_high)]
t_idx = [1, 3334, 6668, n_steps]

for (i, t) in enumerate(compare_times)
    # analytical_soln(S, K, r, sigma, t)
    true_values = price[S_idx, t_idx[i]]
    t_model = t.*ones(1, length(S_values))
    fitted_values = model(t_model, S_model)
    fitted_values = reshape(fitted_values, length(fitted_values))
    display(plot!(S_values, fitted_values, label = "Fitted values"))
    display(plot!(S[S_idx], true_values, label = "True values"))
end


text:
module AbstractGPsMakie

using AbstractGPs
using LinearAlgebra
using Makie

using AbstractGPs: AbstractGP, FiniteGP

export gpsample, gpsample!

include(joinpath("recipes", "symband.jl"))
include(joinpath("recipes", "gpsample.jl"))

include("conversions.jl")

# default plot type
Makie.plottype(::AbstractVector, ::FiniteGP) = SymBand
Makie.plottype(::FiniteGP) = SymBand
Makie.plottype(::AbstractVector, ::AbstractGP) = SymBand

end


text:
cd("lesson-plots")
using Pkg
Pkg.add("Plots")
# gr()
#%%
using Plots

x=1:0.01:10*π
y=sin.(x)

plot(x,y, label="sin(x)")
plot!(xlab="x", ylab="f(x)")

savefig("img1a.png")

y2=sin.(x).^2
plot!(x, y2, label="sin(x)^2", color=:red, line=:dash)

savefig("img1b.png")

xaxis!(:log10)
plot!(legend=:bottomleft)
savefig("img1c.png")

plot!(title="My first plot!")
#%%
plotly()
x=1:0.1:3*π
y=1:0.1:3*π

xx = reshape([xi for xi in x for yj in y],  length(y), length(x))
yy = reshape([yj for xi in x for yj in y],  length(y), length(x))
zz = sin.(xx).*cos.(yy)
plot3d(xx,yy,zz, label=:none, st = :surface)
plot!(xlab="x", ylab="y", zlab="sin(x)*cos(y)")
savefig("img2")

Pkg.add("ORCA")
using ORCA

savefig("img2.png")
#%%
pyplot()
x=0:0.1:2*π
y=sin.(x).^2

plot(x, y, label="sin(x)^2")
savefig("img3.png")
#%%

Pkg.add("LaTeXStrings")
using LaTeXStrings

plot(x, y, label=L"$\sin(x)^2$")
savefig("img3b.png")
#%%

PyPlot.plot(x,y)
PyPlot.xscale("log")

PyPlot.gcf()

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# Apply kinetic operator to wave function in reciprocal space

function op_K( pw::PWGrid, psi::Array{ComplexF64,2}, ik::Int64 )
    #
    out = zeros(ComplexF64,size(psi))
    Ncol = size(psi)[2]

    Ngwk = pw.gkvec.Ngw[ik]
    idx_gkk = pw.gkvec.idx_gk[ik]
    G = pw.gvec.G[:,idx_gkk]

    for is = 1:Ncol
        for ig = 1:Ngwk
            Gk = G[:,ig] + pw.gkvec.kpts[:,ik]
            out[ig,is] = psi[ig,is]*( Gk[1]^2 + Gk[2]^2 + Gk[3]^2 )
        end
    end
    return 0.5*out # two minus signs
end


text:
###############################################################################
#
# Initialize RTDimensions
#
###############################################################################
"""
    create_rt_dims(
                can::Canopy4RT{FT},
                wls::WaveLengths{FT}
    ) where {FT<:AbstractFloat}

Create [`RTDimensions`](@ref), given
- `can` [`Canopy4RT`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
"""
function create_rt_dims(
            can::Canopy4RT{FT},
            wls::WaveLengths{FT}
) where {FT<:AbstractFloat}
    @unpack nAzi, nIncl, nLayer = can;
    @unpack nPAR, nWL, nWLE, nWLF = wls;

    return RTDimensions(nAzi   = nAzi  ,
                        nIncl  = nIncl ,
                        nLayer = nLayer,
                        nPAR   = nPAR  ,
                        nWL    = nWL   ,
                        nWLE   = nWLE  ,
                        nWLF   = nWLF  )
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using YAML
using Statistics
using Cambrian
import Cambrian.mutate

# mutate must be overriden in the global scope (or use eval)
cfg = get_config("../cfg/oneplus.yaml")
mutate(i::Cambrian.FloatIndividual) = Cambrian.mutate(i, cfg.m_rate)

function test_oneplus_evo(fitness::Function, d_fitness::Int)
    cfg = get_config("../cfg/oneplus.yaml"; d_fitness=d_fitness, id="test")
    e = OnePlusEvo{Cambrian.FloatIndividual}(cfg, fitness)

    @test length(e.population) == cfg.population.size
    for i in e.population
        @test all(i.fitness .== -Inf)
    end

    evaluate!(e)

    for i in e.population
        @test i.fitness == fitness(i)
    end
    best = sort(e.population)[end]

    populate!(e)
    new_best = sort(e.population)[end]
    @test !(new_best < best)
    max_count = 0
    for i in e.population
        if i.fitness == new_best.fitness
            max_count += 1
        else
            for f in i.fitness
                @test f == -Inf
            end
        end
    end
    @test max_count == 1
    e.gen += 1

    step!(e)
    @test length(e.population) == cfg.population.size
    for i in e.population
        @test i.fitness == fitness(i)
    end
    new_best = sort(e.population)[end]
    @test !(new_best < best)
    @test e.gen == 2

    print("1+λ step ", string(fitness))
    @timev step!(e)

    run!(e)
    @test length(e.population) == cfg.population.size
    for i in e.population
        @test i.fitness == fitness(i)
    end
    new_best = sort(e.population)[end]
    @test !(new_best < best)
    @test e.gen == cfg.n_gen
end

@testset "1+λ" begin
    @testset "Sum of genes" begin
        function sum_eval(i::Individual)
            [sum(i.genes)]
        end
        test_oneplus_evo(sum_eval, 1)
    end

    @testset "Rosenbrock" begin
        function rosenbrock(i::Individual)
            x = i.genes
            y = -(sum([(1.0 - x[i])^2 + 100.0 * (x[i+1] - x[i]^2)^2
                      for i in 1:(length(x)-1)]))
            [y]
        end
        test_oneplus_evo(rosenbrock, 1)
    end

    @testset "Multi-objective" begin
        function moo(i::Individual)
            [mean(i.genes), std(i.genes), minimum(i.genes)]
        end
        test_oneplus_evo(moo, 3)
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment

basedir = @__DIR__

@testset "read_egret_solution" begin
    solution = UnitCommitment.read_egret_solution(
        "$basedir/../fixtures/egret_output.json.gz",
    )
    for attr in ["Is on", "Production (MW)", "Production cost (\$)"]
        @test attr in keys(solution)
        @test "115_STEAM_1" in keys(solution[attr])
        @test length(solution[attr]["115_STEAM_1"]) == 48
    end
    @test solution["Production cost (\$)"]["315_CT_6"][15:20] ==
          [0.0, 0.0, 884.44, 1470.71, 1470.71, 884.44]
    @test solution["Startup cost (\$)"]["315_CT_6"][15:20] ==
          [0.0, 0.0, 5665.23, 0.0, 0.0, 0.0]
    @test length(keys(solution["Is on"])) == 154
end


text:
import Test

@test 31 == 31
@test 32 == 32
@test 33 == 33


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

m5_1_2_1 = OrderedDict(
  :d => [:a, :m],
  :m => :a
);

# DMA_dag1:

dag5_1_2_1 = DAG("sr5_1_2_1", m5_1_2_1);
show(dag5_1_2_1)

fname = joinpath(mktempdir(), "m5.1.2.1.dot")
to_graphviz(dag5_1_2_1, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

bs = basis_set(dag5_1_2_1)
display(bs)

adjustmentsets = adjustment_sets(dag5_1_2_1, :a, :d)
println("Adjustment sets:")
adjustmentsets |> display

# DMA_dag2:

m5_1_2_2 = OrderedDict(
  [:d, :m] => :a
);

dag5_1_2_2 = DAG("sr5_1_2_12", m5_1_2_2);
show(dag5_1_2_1)

fname = joinpath(mktempdir(), "m5.1.2.2.dot")
to_graphviz(dag5_1_2_2, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

bs = basis_set(dag5_1_2_2)
display(bs)

adjustmentsets = adjustment_sets(dag5_1_2_2, :a, :d)
println("Adjustment sets:")
adjustmentsets |> display

#end