text:
using JpegTurbo_jll
libjpeg_wrap_dir = joinpath(@__DIR__, "..", "..", "gen", "libturbojpeg")
using CEnum
include(joinpath(libjpeg_wrap_dir, "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t
include(joinpath(libjpeg_wrap_dir, "libturbojpeg_common.jl"))
include(joinpath(libjpeg_wrap_dir, "libturbojpeg_api.jl"))

using ColorTypes, ImageCore, FixedPointNumbers

#taken from https://stackoverflow.com/questions/9094691/examples-or-tutorials-of-using-libjpeg-turbos-turbojpeg
function saveimage(filename::String, img::AbstractArray{T,2}; preallocate_max::Bool=false, jpegQual::Int = 93, chromacomp=TJSAMP_420) where T<:Colorant
    #93 chosen to match imagemagick (identified through https://www.imgonline.com.ua/eng/determine-jpeg-quality-result.php)
    #420 chosen to match imagemagick http://www.imagemagick.org/Usage/formats/#jpg_write

    rawimg = rawview(channelview(permuteddimsview(img, (2,1))))
    srcBuf = collect(vec(rawimg))

    nbands, width, height = size(rawimg)
    pitch = Int(width * nbands)
    pixelFormat = nbands == 3 ? TJPF_RGB : TJPF_GRAY
    jpegSubsamp = nbands == 3 ? chromacomp : TJSAMP_GRAY
    flags = 0

    handle = tjInitCompress()
    handle == C_NULL && error("TurboJPEG Error: Unable to initialize compressor object: $(unsafe_string(tjGetErrorStr()))")

    jpegSize = preallocate_max ? Ref(tjBufSize(width, height, jpegSubsamp)) : Ref(Culong(C_NULL))

    jpegBuf = Ref(tjAlloc(jpegSize[]))

    tj_stat = tjCompress2(handle, srcBuf, width, pitch, height, pixelFormat,
                        jpegBuf, jpegSize, jpegSubsamp, jpegQual, flags)

    if tj_stat != 0
        err = tjGetErrorStr()
        tjDestroy(handle)
        handle = C_NULL
        error("TurboJPEG Error: $(unsafe_string(err))")
    end

    jpegBuf_j = unsafe_wrap(Array, jpegBuf[], jpegSize[]; own = false)
    open(filename, "w") do io
       write(io, jpegBuf_j)
    end

    tjstat = tjDestroy(handle) #should deallocate data buffer
    handle = nothing
    tjFree(jpegBuf[])
    return nothing
end

using BenchmarkTools
img = rand(RGB{N0f8}, 1000, 2000);
#img = reshape(range(RGB(1,0,0), RGB(0,1,1), length=200), 10, 20)
@btime saveimage(joinpath(@__DIR__,"out_fileio_imageio.jpg"), $img)

using FileIO
@btime save(joinpath(@__DIR__,"out_fileio_quartimageio.jpg"), $img)

using ImageMagick
@btime save(joinpath(@__DIR__,"out_fileio_imagemagick.jpg"), $img)


text:
# Useful functions to create piecewise stationary bandit instances, stored in matrix of means:  Problem[a,t] is the mean of arm a at time t

using PyPlot

function NoChange(T,mu)
    # a stationary bandit problem with means mu
    K = length(mu)
    Problem = zeros(K,T)
    for k in 1:K
        Problem[k,:]=mu[k]*ones(T)
    end
    return Problem,zeros
end

function ProblemUnif(T,MeansMatrix)
    # generate means rewards with uniformly spaced breakpoints
    # MeansMatrix[a,k] gives the value of \mu_a^{(k+1)}
    (K,Episodes)=size(MeansMatrix)
    nbBreaks = Episodes - 1
    Problem = zeros(K,T)
    # size of each episode
    part = round(Int,T/(nbBreaks+1))
    # filling the matrix
    for c in 0:nbBreaks
	    for arm in 1:K
            Problem[arm,(1+c*part):(c+1)*part]=Problem[arm,(1+c*part):(c+1)*part].+MeansMatrix[arm,c+1]
        end
    end
    if (Episodes*part < T)
        for arm in 1:K
            Problem[arm,(Episodes*part+1):T]=Problem[arm,(Episodes*part+1):T].+MeansMatrix[arm,Episodes]
        end
    end
    BreakPoints = [i*part for i in 1:nbBreaks]
    return Problem,BreakPoints
end

function ProblemGeneral(MeansMatrix,Landmarks)
    # between Landmarks[k] and Landmarks[k+1] the mean of arm a is MeansMatrix[a,k]
    # one must have Landmarks[1]=1 and Landmarks[L]=T
    L = length(Landmarks)
    (K,l) = size(MeansMatrix)
    if L>(l+1)
        return "wrong dimensions"
    else
        T= Landmarks[L]
        Problem = zeros(Float32,K,T)
        for k in 1:(L-1)
            for arm in 1:K
                Problem[arm,Landmarks[k]:Landmarks[k+1]]=MeansMatrix[arm,k]*ones(Float32,Landmarks[k+1]-Landmarks[k]+1)
            end
        end
        return Problem,Landmarks[2:(L-1)]
    end
end

function RandomProblem(K,T,nbreaks,p,Deltamin = 0.05,Deltamax = 0.3,Spacemin=500)
    # generates a random instance with K arms up to horizon T with at most nbreaks breakpoint
    # p: probability that a change occurs for each arm in each breakpoint
    # Deltamin - Deltamax : minimal - maximal magnitude of a change-point
    # Spacemin : minimal distance between two breakpoints
    cont = true
    # choose the breakpoint (resample until they match the Spacemin condition)
    Landmarks = zeros(nbreaks+2)
    while (cont)
        cont = false
        Landmarks = sort(floor.(Int,T*rand(nbreaks+2)))
        Landmarks[1]=1
        Landmarks[nbreaks+2]=T
        for i in 1:(nbreaks+1)
            if ((Landmarks[i+1]-Landmarks[i])<Spacemin)
                cont = true
            end
        end
    end
    # matrix of means
    MeansMatrix = zeros(K,nbreaks+1)
    # random initial means
    for a in 1:K
        MeansMatrix[a,1]=rand()
    end
    # in each breakpoint, there is a probability p that and arm keeps its previous value
    for k in 2:(nbreaks+1)
        for a in 1:K
            prevmean = MeansMatrix[a,k-1]
            if (rand()>p)
                MeansMatrix[a,k]=prevmean
            else
                gap = Deltamin + (Deltamax - Deltamin)*rand()
                if (rand()<0.5)
                    gap = -gap
                end
                if ((prevmean + gap)>1)||((prevmean + gap)<0)
                    MeansMatrix[a,k] = prevmean-gap
                else
                    MeansMatrix[a,k] = prevmean + gap
                end
            end
        end
    end
    return ProblemGeneral(MeansMatrix,Landmarks)
end

function ProblemToMatrix(Problem,nbreaks)
    # return the matrix of means and landmarks that allow to generate the problem
    (K,T)=size(Problem)
    Landmarks = zeros(Int,nbreaks+2)
    MeansMatrix = zeros(K,nbreaks+1)
    ndetect = 0
    MeansMatrix[:,1]=MeansMatrix[:,1].+Problem[:,1]
    Landmarks[1]=1
    for t = 2:T
        if (Problem[:,t]!=Problem[:,t-1])
            ndetect = ndetect+1
            Landmarks[ndetect+1]=t
            MeansMatrix[:,ndetect+1]=MeansMatrix[:,ndetect+1].+Problem[:,t]
        end
    end
    Landmarks[ndetect+2]=T
    return MeansMatrix[:,1:(ndetect+1)],Landmarks[1:(ndetect+2)]
end

function PlotProblem(Pb)
    # visualize the evolution of means in a bandit instance
    (K,T)=size(Pb)
    for k in 1:K
        plot(1:T,Pb[k,:],label="arm $(k)")
    end
    legend()
    return
end

# matrix of means corresponding to Problem1 in the paper
MeansMatrix1 = reshape([0.9,0.5,0.3,0.9,0.2,0.3,0.1,0.2,0.3,0.1,0.7,0.3,0.1,0.7,0.5],3,5)
# matrix of means corresponding to Problem 2 in the paper
MeansMatrix2 = reshape([0.9,0.5,0.4,0.7,0.4,0.5,0.5,0.3,0.6,0.3,0.2,0.7,0.1,0.1,0.8],3,5)

# a few other saved Benchmarks (before trying any expes on them)
# WorseBecomeGood : changes not evenly spaced, worse arm eventually becomes good (PB3 in the paper)
# OneOptChange: only one meaningful change (PB4 in the paper)

ext: jl


text:
using BritishNationalGrid
using Test

const BNG = BritishNationalGrid

bngp_type(::BNGPoint{T}) where {T} = T

function check_bng_to_lonlat(easting, northing, lon, lat, tol)
    p = BNGPoint(easting, northing)
    plon, plat = lonlat(p)
    isapprox(plon, lon, atol=tol) && isapprox(plat, lat, atol=tol)
end
function check_lonlat_to_bng(lon, lat, easting, northing, tol)
    p = BNGPoint(lon=lon, lat=lat)
    isapprox(easting, p.e, atol=tol) && isapprox(northing, p.n, atol=tol)
end

## Construction
@testset "Construction" begin
    @test BNGPoint(500_000., 1_000_000).e == 500_000.
    @test BNGPoint(500_000., 1_000_000).n == 1_000_000.
    @test bngp_type(BNGPoint(0, 0)) == Int
    @test bngp_type(BNGPoint(0, Float64(0))) == Float64
    # Outside grid
    @test_throws ArgumentError BNGPoint(0, 1_300_001)
    @test_throws ArgumentError BNGPoint(700_001, 0)
    @test_throws ArgumentError BNGPoint(0, -1)
    @test_throws ArgumentError BNGPoint(-1, 0.)
    # Incorrect type
    @test_throws MethodError BNGPoint(1im, 1)
    # With squares
    @test_throws ArgumentError BNGPoint(0, 100_001, "SV")
    @test_throws ArgumentError BNGPoint(-1, 0, "SV")
    @test_throws ArgumentError BNGPoint(0, 0, "ZZ")
    @test BNGPoint(0, 0, "SV").e == 0
    @test BNGPoint(0, 0, "SV").n == 0
    @test BNGPoint(0, 0, "TV") == BNGPoint(500_000, 0)
    @test BNGPoint(1000, 1000, "TV") == BNGPoint(501_000, 1000)
end

## Conversion lonlat ↔ BNG
const lonlattol = 1e-6
const gridtol = 1 # Because the online calculator rounds of course

# Points taken from BGS converter at:
#   http://www.bgs.ac.uk/data/webservices/convertForm.cfm
@testset "Conversion" begin
    @test check_bng_to_lonlat(429157, 623009, -1.54000791003, 55.4999996103, lonlattol)
    @test check_bng_to_lonlat(0, 0, -7.55716018087, 49.766807227, lonlattol)
    @test check_lonlat_to_bng(-1.54000791003, 55.4999996103, 429157, 623009, gridtol)
    @test check_lonlat_to_bng(0, 60, 511648, 1125592, gridtol)
end

## Square identification
@testset "Square ID" begin
    @test BNG.SQUARE_NAMES[1,1] == "SV"
    @test BNG.SQUARE_NAMES[end,end] == "JM"
    @test square(BNGPoint(429157, 623009)) == "NU"
end

## Grid reference formatting
@testset "Gridref format" begin
    @test_throws ArgumentError gridref(BNGPoint(0, 0), 1)
    @test_throws ArgumentError gridref(BNGPoint(0, 0), 0)
    @test_throws ArgumentError gridref(BNGPoint(0, 0), 12, true)
    @test_throws ArgumentError gridref(BNGPoint(0, 0), 14)
    @test gridref(BNGPoint(429157, 623009), 8, true) == "NU 2915 2300"
    @test gridref(BNGPoint(429157, 623009)) == "4291 6230"
    @test gridref(BNGPoint(429157, 623009), 2) == "4 6"
    @test gridref(BNGPoint(429157, 623009), 2, true) == "NU 2 2"
    @test gridref(BNGPoint(429157, 623009), 4) == "42 62"
    @test gridref(BNGPoint(429157, 623009), 4, true) == "NU 29 23"
    @test gridref(BNGPoint(429157, 623009), 10, true, "") == "NU2915723009"
    @test gridref(BNGPoint(429157, 623009), 10, true, "_∘") == "NU_∘29157_∘23009"
end


text:
function enter(ex)
  # !isdefined(Atom, :Debugger) && @eval include(joinpath(dirname($@__FILE__), "debugger.jl"))
  # Base.invokelatest(Atom.Debugger.enter, :($ex))
end

isdebugging() = false # isdefined(Atom, :Debugger) && Debugger.isdebugging()


text:
# TODO: Need to finish
#       Need to test
function fit_fs_tp!(init::StateFS,
                    cfs::ConstantsFS,
                    dfs::DataFS;
                    nmcmc::Int, nburn::Int,
                    # These args are for iMCMC:
                    batchprop::Float64=0.1,
                    batchsizes=nothing,
                    prior_thin::Int=2,
                    temper::Float64=1.0, mb_update_burn_prop=0.6, anneal=false,
                    min_temper::Float64=1.0,
                    # End of iMCMC args.
                    tfs::Union{Nothing, Vector{TunersFS}}=nothing,
                    monitors=[monitor1, monitor2],
                    fix::Vector{Symbol}=Symbol[],
                    thins::Vector{Int}=[2, nsamps_to_thin(10, nmcmc)],
                    thin_dden::Int=1,
                    printFreq::Int=0,
                    checkpoint=0,
                    computeDIC::Bool=false, computeLPML::Bool=false,
                    computedden::Bool=false,
                    sb_ibp::Bool=false,
                    use_repulsive::Bool=true,
                    Z_marg_lamgam::Float64=1.0,
                    Z_marg_lamgam_decay_rate::Float64=100.0,
                    Z_marg_lamgam_min::Float64=0.05,
                    verbose::Int=1,
                    time_updates=false,
                    seed::Int=-1)

  @assert 0 <= mb_update_burn_prop <= 1

  printMsg(iter::Int, msg::String) = if printFreq > 0 && iter % printFreq == 0
    print(msg)
  end

  # Set random seed if needed
  if seed >= 0
    Random.seed!(seed)
  end

  # Print temperature
  println("temperature: $(temper)")

  # Create Tuners if needed
  if tfs == nothing
    tfs = TunersFS(Tuners(dfs.data.y, cfs.constants.K),
                   init.theta,
                   dfs.X)
  end

  if verbose >= 1
    fixed_vars_str = join(fix, ", ")
    if fixed_vars_str == ""
      fixed_vars_str = "nothing"
    end
    println("fixing: $fixed_vars_str")
    println("Use stick-breaking IBP: $(sb_ibp)")
    println("Z_marg_lamgam: $(Z_marg_lamgam)")
    println("Z_marg_lamgam_decay_rate: $(Z_marg_lamgam_decay_rate)")
    println("Z_marg_lamgam_min: $(Z_marg_lamgam_min)")
    println("use_repulsive: $(use_repulsive)")
    flush(stdout)
  end

  @assert printFreq >= -1
  if printFreq == 0
    numPrints = 10
    printFreq = Int(ceil((nburn + nmcmc) / numPrints))
  end

  # Instantiate (but not initialize) CPO stream
  if computeLPML
    cpoStream = MCMC.CPOstream{Float64}()
  end

  # DIC
  if computeDIC
    dicStream = initDicStream(dfs.data)
    loglikeDIC(param::DICparam) = computeLoglikeDIC(dfs.data, param)

    convertStateToDicParam(s::State)::DICparam = let
      _convertStateToDicParam(s, cfs.constants, dfs.data)
    end
  end

  # Cache for data density
  dden = Matrix{Vector{Float64}}[]

  # Cache for loglikelihood
  ll = Float64[]

  # Update function
  function update!(state, iter::Int, out)
    # Whether or not to marginalize over lambda and gamma.
    # We want to do this more often at the beginning, and less at the end.
    zmarg = ((Z_marg_lamgam - Z_marg_lamgam_min) *
             exp(-iter/Z_marg_lamgam_decay_rate) +
             Z_marg_lamgam_min) > rand()

    # Iterations to do minibatch update
    iters_mb_update_all_params = round(Int, nburn * mb_update_burn_prop)

    # Calculate current temperature
    curr_temper = if anneal
      # Gradually cool the temperature so that it is 1 when
      # we reach iteration `iters_mb_update_all_params`, which
      # is less than `nburn`.
      itmax = iters_mb_update_all_params
      power = (itmax - iter + 1) / itmax
      clamp(temper^power, min_temper, Inf)
    else
      temper
    end
    if verbose > 0
      print(" -- Current temper: $(curr_temper)")
    end

    # Update state using trained prior
    minibatch_update_all_params = iter < iters_mb_update_all_params || !anneal
    update_via_trained_prior!(state, dfs, cfs, tfs, batchprop, prior_thin,
                              batchsizes=batchsizes,
                              fix=fix, use_repulsive=use_repulsive,
                              Z_marg_lamgam=zmarg, sb_ibp=sb_ibp,
                              time_updates=time_updates, temper=curr_temper,
                              minibatch_update_all_params=minibatch_update_all_params,
                              verbose=verbose-2)

    # Pull out inner componenets for convenience
    s = state.theta
    c = cfs.constants
    d = dfs.data

    # Append loglike
    append!(ll, compute_marg_loglike(s, c, d, 1.0) / curr_temper)

    if computedden && iter > nburn && (iter - nburn) % thin_dden == 0
      # NOTE: `datadensity(s, c, d)` returns an (I x J) matrix of vectors of
      # length g.
      append!(dden, [datadensity(s, c, d)])
    end

    if computeLPML && iter > nburn
      # Inverse likelihood for each data point
      like = [[compute_like(i, n, s, c, d)
               for n in 1:d.N[i]] for i in 1:d.I]

      # Update (or initialize) CPO
      MCMC.updateCPO(cpoStream, vcat(like...))

      # Add to printMsg
      printMsg(iter, " -- LPML: $(MCMC.computeLPML(cpoStream))")
    end

    if computeDIC && iter > nburn
      # Update DIC
      MCMC.updateDIC(dicStream, s, updateParams,
                     loglikeDIC, convertStateToDicParam)

      # Add to printMsg
      printMsg(iter, " -- DIC: $(MCMC.computeDIC(dicStream, loglikeDIC,
                                                 paramMeanCompute))")

      DICg = MCMC.DIC_gelman(ll[(nburn+1):end])
      printMsg(iter, " -- DIC_gelman: $(DICg)")
    end

    printMsg(iter, "\n")
    flush(stdout)
  end

  println("Running Gibbs sampler ...")
  samples, state = MCMC.gibbs(deepcopy(init), update!, monitors=monitors,
                              thins=thins, nmcmc=nmcmc, nburn=nburn,
                              printFreq=printFreq, loglike=ll,
                              printlnAfterMsg=false)

  out = Dict(:samples => samples,
             :lastState => state,
             :init => init,
             :ll => ll,
             :c => cfs,
             :d => dfs)

  if computeDIC || computeLPML
    LPML = computeLPML ? MCMC.computeLPML(cpoStream) : NaN
    Dmean, pD = computeDIC ? MCMC.computeDIC(dicStream, loglikeDIC,
                                             paramMeanCompute,
                                             return_Dmean_pD=true) : (NaN, NaN)

    DICg = MCMC.DIC_gelman(ll[(nburn+1):end])

    metrics = Dict(:LPML => LPML,
                   :DIC => Dmean + pD,
                   :Dmean => Dmean,
                   :pD => pD,
                   :DICg => DICg)
    println()
    println("metrics:")
    for (k, v) in metrics
      out[k] = v
      println("$k => $v")
    end
    flush(stdout)
  end

  if computedden
    out[:dden] = dden
  end

  out[:nburn] = nburn
  out[:nmcmc] = nmcmc
  out[:batchprop] = batchprop
  out[:batchsizes] = batchsizes
  out[:prior_thin] = prior_thin
  out[:temper] = temper
  out[:anneal] = anneal

  return out
end


text:
"""
Struct for European Option

		euOption=EuropeanOptionND(T::num1,K::num2,isCall::Bool=true) where {num1 <: Number,num2 <: Number}

Where:\n
		T	=	Time to maturity of the Option.
		K	=	Strike Price of the Option.
		isCall  = true for CALL, false for PUT.
"""
mutable struct EuropeanOptionND{num1 <: Number, num2 <: Number, numtype <: Number} <: EuropeanBasketPayoff{numtype}
    T::num1
    K::num2
    isCall::Bool
    function EuropeanOptionND(T::num1, K::num2, isCall::Bool = true) where {num1 <: Number, num2 <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        elseif K <= 0.0
            error("Strike Price must be positive")
        else
            zero_typed = zero(num1) + zero(num2)
            return new{num1, num2, typeof(zero_typed)}(T, K, isCall)
        end
    end
end

export EuropeanOptionND;

function payoff(S::Array{abstractMatrix}, euPayoff::EuropeanOptionND, rfCurve::abstractZeroRateCurve, T1::num2 = maturity(euPayoff)) where {abstractMatrix <: AbstractMatrix{num}, num2 <: Number} where {abstractZeroRateCurve <: AbstractZeroRateCurve, num <: Number}
    r = rfCurve.r
    T = euPayoff.T
    iscall = euPayoff.isCall ? 1 : -1
    (Nsim, NStep) = size(S[1])
    NStep -= 1
    index1 = round(Int, T / T1 * NStep) + 1
    K = euPayoff.K
    #ST_all=[ sum(x_i[j,index1] for x_i in S) for j in 1:Nsim]
    ST_all = S[1][:, index1]
    for i = 2:length(S)
        ST_all += S[i][:, index1]
    end
    payoff2 = max.(iscall * (ST_all .- K), 0.0)

    return payoff2 * exp(-integral(r, T))
end


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

name = "librdkafka"
version = v"1.5.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/edenhill/librdkafka/archive/v1.5.0.tar.gz", "f7fee59fdbf1286ec23ef0b35b2dfb41031c8727c90ced6435b8cf576f23a656")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd librdkafka*
if [[ "${target}" != *-freebsd* ]]; then
rm -f /opt/${target}/${target}/sys-root/usr/lib/libcrypto.*
rm -f /opt/${target}/${target}/sys-root/usr/lib/libssl.*
rm -f /opt/${target}/${target}/sys-root/usr/lib/libsasl2.*
fi
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DRDKAFKA_BUILD_EXAMPLES=OFF -DRDKAFKA_BUILD_TESTS=OFF ..
make -j${nproc}
make install
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Lz4_jll", uuid="5ced341a-0733-55b8-9ab6-a4889d929147"))
    Dependency(PackageSpec(name="Zstd_jll", uuid="3161d3a3-bdf6-5164-811a-617609db77b4"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function find_tsp_tour(Ct::Matrix{Float64})
    scale_factor = 1000
    dist_mtx = round.(Int, Ct .* scale_factor)

    tsp_tour, tsp_tour_len = Concorde.solve_tsp(dist_mtx)

    @assert tsp_tour[1] == 1

    return tsp_tour
end

function find_tsp_tour(x::Vector{Float64}, y::Vector{Float64})
    scale_factor = 1000

    tsp_tour, tsp_tour_len = Concorde.solve_tsp(x * scale_factor, y * scale_factor; dist="EUC_2D")

    @assert tsp_tour[1] == 1

    return tsp_tour
end

function exact_partitioning(initial_tour, Ct, Cd; flying_range=MAX_DRONE_RANGE)
    n, _ = size(Ct)

    r = initial_tour
    T = fill(Inf, n, n)
    M = fill(-99, n, n)
    for i in 1:n-1
        for j in i+1:n
            if j == i + 1
                T[i, j] = Ct[r[i], r[j]]
                M[r[i], r[j]] = -1
            else
                for k in i+1:j-1
                    Tk1 = Cd[r[i], r[k]] + Cd[r[k], r[j]]
                    if Tk1 <= flying_range
                        Tk2 = sum([Ct[r[l], r[l+1]] for l in i:k-2]) +
                                Ct[r[k-1], r[k+1]] +
                                sum([Ct[r[l], r[l+1]] for l in k+1:j-1])
                        Tk = max(Tk1, Tk2)
                        if Tk < T[i, j]
                            T[i, j] = Tk
                            M[r[i], r[j]] = r[k]
                        end
                    end
                end

            end
        end
    end

    V = zeros(n)
    P = fill(-1, n)

    V[1] = 0
    for i in 2:n
        VV = [V[k] + T[k, i] for k in 1:i-1]
        V[i] = minimum(VV)
        P[i] = r[argmin(VV)]
    end

    # Retrieving solutions.
    combined_nodes = Int[]
    current_idx = n
    current = r[current_idx]
    while current != -1
        pushfirst!(combined_nodes, current)
        current_idx = findfirst(x -> x == current, r)
        current = P[current_idx]
    end
    drone_only_nodes = Int[]
    drone_route = Int[]
    push!(drone_route, combined_nodes[1])
    @assert combined_nodes[1] == r[1]
    for i in 1:length(combined_nodes)-1
        j1 = combined_nodes[i]
        j2 = combined_nodes[i+1]
        if M[j1, j2] != -1
            push!(drone_only_nodes, M[j1, j2])
            push!(drone_route, M[j1, j2])
        end
        push!(drone_route, j2)
    end
    truck_route = setdiff(initial_tour, drone_only_nodes)

    obj_val = objective_value(truck_route, drone_route, Ct, Cd)
    final_time = V[end]

    @assert isapprox(obj_val, final_time)
    return final_time, truck_route, drone_route
end

function two_point_move(tour, i, j)
    if i >= j
        return tour, false
    end

    tmp = copy(tour)
    tmp[j] = tour[i]
    tmp[i] = tour[j]
    return tmp, true
end

function one_point_move(tour, i, j)
    tmp = copy(tour)
    deleteat!(tmp, i)
    insert!(tmp, j, tour[i])
    return tmp, true
end

function two_opt_move(tour, i, j)
    if i >= j
        return tour, false
    end

    tmp = copy(tour)
    tmp[i:j] = tour[j:-1:i]
    return tmp, true
end

# Main function to call
function tsp_ep_all(
    x_coordinates,
    y_coordinates,
    truck_cost_factor,
    drone_cost_factor;
    local_search_methods=[two_point_move, one_point_move, two_opt_move],
    flying_range=MAX_DRONE_RANGE,
    time_limit=MAX_TIME_LIMIT
)
    """
    Runs `TSP-ep-all` heuristic algorithm of Agatz et al.

    `x_coordinates[1]`, `y_coordinates[1]`: the coordinates of the depot, then followed by all customer location coordinates
    `truck_cost_factor`: as defined in Agatz et al. instances
    `drone_cost_factor`: as defined in Agatz et al. instances
    """

    Ct, Cd = cost_matrices_with_dummy(x_coordinates, y_coordinates, truck_cost_factor, drone_cost_factor)
    n_nodes = length(x_coordinates)
    n1, n2 = size(Ct)
    @assert n_nodes + 1 == n1
    @assert size(Ct) == size(Cd)

    tsp_tour = find_tsp_tour(x_coordinates, y_coordinates)
    push!(tsp_tour, n_nodes+1) # adding a dummy node for the returning depot

    return tsp_ep_all(Ct, Cd, tsp_tour, flying_range=flying_range, local_search_methods=local_search_methods, time_limit=time_limit)
end

function tsp_ep_all(
    Ct,
    Cd,
    init_tour;
    local_search_methods=[two_point_move, one_point_move, two_opt_move],
    flying_range=MAX_DRONE_RANGE,
    time_limit=MAX_TIME_LIMIT
)
    n, _ = size(Ct)

    improved = true
    best_obj, t_route, d_route = exact_partitioning(init_tour, Ct, Cd, flying_range=flying_range)

    best_tour = copy(init_tour)
    best_t_route = copy(t_route)
    best_d_route = copy(d_route)

    if isempty(local_search_methods)
        return best_obj, best_t_route, best_d_route
    end

    time0 = time()

    is_time_over = false

    while improved && !is_time_over
        improved = false
        cur_best_obj = best_obj
        cur_best_tour = copy(best_tour)
        cur_best_t_route = copy(best_t_route)
        cur_best_d_route = copy(best_d_route)

        for i in 2:n-1
            if is_time_over
                break
            end

            for local_search in local_search_methods
                for j in 2:n-1
                    if is_time_over
                        break
                    end

                    new_tour, is_valid = local_search(best_tour, i, j)
                    is_time_over = time() - time0 > time_limit

                    if is_valid
                        ep_time, t_route, d_route = exact_partitioning(new_tour, Ct, Cd, flying_range=flying_range)
                        if ep_time < cur_best_obj
                            cur_best_tour = copy(new_tour)
                            cur_best_t_route = copy(t_route)
                            cur_best_d_route = copy(d_route)
                            cur_best_obj = ep_time
                            improved = true
                        end
                    end
                end
            end
        end

        if improved
            best_obj = cur_best_obj
            best_tour = cur_best_tour
            best_t_route = cur_best_t_route
            best_d_route = cur_best_d_route
        end
    end

    return best_obj, best_t_route, best_d_route
end


text:
struct Scordis
    dis
    fm
    Sinv::Matrix{Float64}
    cutoff::Real
    nlv::Int64
end

struct Odis
    dis
    fm
    cutoff::Real
    nlv::Int64
end

"""
    scordis(object::Union{Pca, Plsr}; nlv = nothing, rob = true, alpha = .01)
Compute the score distances (SDs) from a PCA or PLS model

* `object` : The fitted model.
* `nlv` : Nb. components (PCs or LVs) to consider. If nothing,
    it is the maximum nb. of components.
* `rob` : If true, the moment estimation of the distance cutoff is robustified.
    This may be relevant after robust PCA or PLS on small data sets
        containing extreme values.
* `alpha` : Risk-I level used for computing the cutoff detecting extreme values.

SDs are the Mahalanobis distances of the projections of row observations on the
score plan to the center of the score space.

They are computed for the training `Xtrain` and the eventual `X`.

A cutoff is computed from the training, using a moment estimation of the parameters of
a Chi-squared distrbution for SD^2 (see e.g. Pomerantzev 2008).
In the output, column dstand is a standardized distance, defined as SD / cutoff.
A value dstand > 1 may be considered as extreme.

The Winisi "GH" is also provided (usually considered as extreme if GH > 3).

## References

M. Hubert, P. J. Rousseeuw, K. Vanden Branden (2005). ROBPCA: a new approach to robust
principal components analysis. Technometrics, 47, 64-79.

Pomerantsev, A.L., 2008. Acceptance areas for multivariate classification derived by
projection methods. Journal of Chemometrics 22, 601-609. https://doi.org/10.1002/cem.1147
"""
function scordis(object::Union{Pca, Plsr};
    nlv = nothing, rob = true, alpha = .01)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    T = @view(object.T[:, 1:nlv])
    S = Statistics.cov(T, corrected = false)
    LinearAlgebra.inv!(cholesky!(S))   # ==> S := Sinv
    d2 = mahsq(T, zeros(nlv)', S)
    d = sqrt.(d2)
    if !rob
        mu = mean(d2)
        s2 = var(d2)
    else
        mu = median(d2)
        s2 = mad(d2)^2
    end
    nu = 2 * mu^2 / s2
    cutoff = sqrt(mu / nu * quantile.(Distributions.Chisq(nu), 1 - alpha))
    dstand = d / cutoff
    dis = DataFrame((d = d, dstand = dstand, gh = d2 / nlv))
    Scordis(dis, object, S, cutoff, nlv)
end

function predict(object::Scordis, X)
    nlv = object.nlv
    T = transform(object.fm, X; nlv = nlv)
    d2 = mahsq(T, zeros(nlv)', object.Sinv)
    d = sqrt.(d2)
    dstand = d / object.cutoff
    dis = DataFrame((d = d, dstand = dstand, gh = d2 / nlv))
    (dis = dis,)
end

"""
    odis(object::Union{Pca, Plsr}, X; nlv = nothing, rob = true, alpha = .01)
Compute the orthogonal distances (ODs) from a PCA or PLS model

* `object` : The fitted model.
* `X` : X-data that were used to compute the model (training).
* `nlv` : Nb. components (PCs or LVs) to consider. If nothing, it is the maximum
    nb. of components.
* `rob` : If true, the moment estimation of the distance cutoff is robustified.
    This may be relevant after robust PCA or PLS on small data sets
        containing extreme values.
* `alpha` : Risk-I level used for computing the cutoff detecting extreme values.

ODs ("X-residuals") are the Euclidean distances of row observations to their projections
to the score plan (see e.g. Hubert et al. 2005, Van Branden & Hubert 2005, p. 66;
Varmuza & Filzmoser, 2009, p. 79).

They are computed for the training `Xtrain` and the eventual `X`.

A cutoff is computed from the training, using a moment estimation of the parameters of a
Chi-squared distrbution for OD^2 (see Nomikos & MacGregor 1995, and Pomerantzev 2008).
In the output, column dstand is a standardized distance, defined as OD / cutoff.
A value dstand > 1 may be considered as extreme.

The cutoff for detecting extreme OD values is computed using a moment estimation of
a Chi-squared distrbution for the squared distance.

## References

M. Hubert, P. J. Rousseeuw, K. Vanden Branden (2005). ROBPCA: a new approach to robust principal
components analysis. Technometrics, 47, 64-79.

Nomikos, P., MacGregor, J.F., 1995. Multivariate SPC Charts for Monitoring Batch Processes.
null 37, 41-59. https://doi.org/10.1080/00401706.1995.10485888

Pomerantsev, A.L., 2008. Acceptance areas for multivariate classification derived by projection methods.
Journal of Chemometrics 22, 601â-609. https://doi.org/10.1002/cem.1147

K. Vanden Branden, M. Hubert (2005). Robuts classification in high dimension based on the SIMCA method.
Chem. Lab. Int. Syst, 79, 10-21.

K. Varmuza, P. Filzmoser (2009). Introduction to multivariate statistical analysis in chemometrics.
CRC Press, Boca Raton.
"""
function odis(object::Union{Pca, Plsr}, X;
        nlv = nothing, rob = true, alpha = .01)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    E = xresid(object, X; nlv = nlv)
    d = sqrt.(sum(E .* E, dims = 2))
    d2 = d.^2
    if(!rob)
        mu = mean(d2)
        s2 = var(d2)
    else
        mu = median(d2)
        s2 = mad(d2)^2
    end
    nu = 2 * mu^2 / s2
    cutoff = sqrt(mu / nu * quantile.(Distributions.Chisq(nu), 1 - alpha))
    dstand = d / cutoff
    dis = DataFrame((d = d, dstand = dstand))
    Odis(dis, object, cutoff, nlv)
end

function predict(object::Odis, X)
    nlv = object.nlv
    E = xresid(object.fm, X; nlv = nlv)
    d = sqrt.(sum(E .* E, dims = 2))
    dstand = d / object.cutoff
    dis = DataFrame((d = d, dstand = dstand))
    (dis = dis,)
end


text:
using LinearAlgebra, ForwardDiff, NLPModels, Random
Random.seed!(0)

export LBFGS_StrongWolfe

include("StrongWolfe.jl")

function LBFGS_StrongWolfe(
  nlp::AbstractNLPModel;
  m = 5,  # Memory reserved
  a1 = 1.0, # Initial Step
  η1::Float64 = 1e-4, # Armijo Condition
  η2::Float64 = 0.8, # StrongWolfe
  p::Float64 = 0.9, # Line Search Step
  atol::Real = 1e-6,
  rtol::Real = 1e-6,
  max_eval::Int = 1000,
  max_iter::Int = 0,
  max_time::Float64 = 10.0,
  )

  if !unconstrained(nlp)
    error("Problem is not unconstrained")
  end

  x = copy(nlp.meta.x0)
  n = nlp.meta.nvar

  f(x) = obj(nlp, x)
  ∇f(x) = grad(nlp, x)
  fx = f(x)
  ∇fx = ∇f(x)

  ϵ = atol + rtol * norm(∇fx)
  t₀ = time()

  iter = 0
  Δt = time() - t₀
  solved = norm(∇fx) < ϵ # First order stationary
  tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

  status = :unknown

  @info log_header(
  [:iter, :fx, :ngx, :nf, :Δt],
  [Int, Float64, Float64, Int, Float64],
  hdr_override=Dict(:fx => "f(x)", :ngx => "‖∇f(x)‖", :nf => "#f")
  )
  @info log_row(
  Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
  )

  Id = Matrix(1.0I, n, n)# Defino uma H0
  γ = 1
  r = - ∇fx
  y = zeros(n ,m)
  s = zeros(n, m)
  a = zeros(m)

  α, iteracao = StrongWolfe_LineSearch(nlp, x, r, a1, η1, η2, p) #Calcular um passo com α de Wolfe
  xt = x + α * r
  s[:, 1] = xt - x
  y[:, 1] = ∇f(xt) - ∇fx
  γ = dot(s[:, 1], y[:, 1]) / dot(y[:, 1], y[:, 1])
  H0 = γ * Id
  x = xt

  while !(solved || tired)

    r = - H0 * ∇fx
    α, iteracao = StrongWolfe_LineSearch(nlp, x, r, a1, η1, η2, p) #Calcular um passo com α de Wolfe
    xt = x + α * r #Atualiza o passo xₖ₊₁ = xₖ - αₖ r

    #calcular sₖ e yₖ
    s[:, iter%m+1] = xt - x
    y[:, iter%m+1] = ∇f(xt) - ∇fx
    x = xt

    H0 += ((s[:, iter%m+1] - H0 * y[:, iter%m+1]) * (s[:, iter%m+1] - H0 * y[:, iter%m+1])') / ((s[:, iter%m+1] - H0 * y[:, iter%m+1])' * y[:, iter%m+1])

    if iter > m
      q = ∇f(x)
      for i = ((iter - 1) % m + 1):-1:((iter - m) % m + 1)
        ρ = 1 / dot(y[:,i],s[:,i])
        a[i] = ρ * dot(s[:,i], q)
        q -= a[i] * y[:,i]
      end
      r = γ * H0 * q
      for i = ((iter - m) % m + 1):((iter - 1) % m + 1)
        ρ = 1 / dot(y[:,i],s[:,i])
        β = ρ * dot(y[:,i],r)
        r += s[:,i] * (a[i] - β) # retorna a direção r = Hₖ∇f
        r = -r
      end
    end

    fx = f(x)
    ∇fx = ∇f(x)

    iter += 1

    Δt = time() - t₀
    solved = norm(∇fx) < ϵ # First order stationary
    tired = neval_obj(nlp) ≥ max_eval > 0 || iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

    @info log_row(
    Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
    )
  end

  if solved
    status = :first_order
  elseif tired
    if neval_obj(nlp) ≥ max_eval > 0
      status = :max_eval
    elseif iter ≥ max_iter > 0
      status = :max_iter
    elseif Δt ≥ max_time > 0
      status = :max_time
    end
  end

  return GenericExecutionStats(
  status,
  nlp,
  solution=x,
  objective=f(x),
  dual_feas=norm(∇fx),
  elapsed_time=Δt,
  iter=iter
  )
end


text:
# implementation of MLJ measure interface for LossFunctions.jl

# Supervised Loss -- measure traits

is_measure_type(::Type{<:SupervisedLoss})          = true
orientation(::Type{<:SupervisedLoss})              = :loss
reports_each_observation(::Type{<:SupervisedLoss}) = true
is_feature_dependent(::Type{<:SupervisedLoss})     = false

MMI.supports_weights(::Type{<:SupervisedLoss}) = true
MMI.name(M::Type{<:SupervisedLoss})            = split(string(M), '.')[end]*"()"
MMI.docstring(M::Type{<:SupervisedLoss})       = name(M)

## DISTANCE BASED LOSS FUNCTION

MMI.prediction_type(::Type{<:DistanceLoss}) = :deterministic
MMI.target_scitype(::Type{<:DistanceLoss})  = AbstractArray{<:Continuous}

function value(measure::DistanceLoss, yhat, X, y, ::Nothing,
                ::Val{false}, ::Val{true})
    return measure(yhat, y)
end

function value(measure::DistanceLoss, yhat, X, y, w,
                ::Val{false}, ::Val{true})
    return w .* measure(yhat, y) ./ (sum(w)/length(y))
end

## MARGIN BASED LOSS FUNCTIONS

MMI.prediction_type(::Type{<:MarginLoss}) = :probabilistic
MMI.target_scitype(::Type{<:MarginLoss})  = AbstractArray{<:Binary}

# convert a Binary vector into vector of +1 or -1 values
# (for testing only):
pm1(y) = Int8(2) .* (Int8.(int(y))) .- Int8(3)

# rescale [0, 1] -> [-1, 1]
_scale(p) = 2p - 1

function value(measure::MarginLoss, yhat, X, y, ::Nothing,
                ::Val{false}, ::Val{true})
    check_pools(yhat, y)
    probs_of_observed = broadcast(pdf, yhat, y)
    return broadcast(measure, _scale.(probs_of_observed), 1)
end

function value(measure::MarginLoss, yhat, X, y, w,
                ::Val{false}, ::Val{true})
    return w .* value(measure, yhat, X, y, nothing) ./ (sum(w)/length(y))
end


text:
function CachedOperator(::Type{BandedBlockBandedMatrix}, op::Operator)
    l,u = blockbandwidths(op)
    λ,μ = subblockbandwidths(op)
    data = BandedBlockBandedMatrix{eltype(op)}(undef,
        blocklengths(rangespace(op))[1:0],blocklengths(domainspace(op))[1:0],
        (l,u), (λ,μ))

    CachedOperator(op,data,size(data),domainspace(op),rangespace(op),(-l,u),false)
end

# Grow cached operator
#
function resizedata!(B::CachedOperator{T,<:BandedBlockBandedMatrix{T}}, ::Colon, col::Integer) where {T<:Number}
    if col > size(B,2)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    if col > B.datasize[2]
        l,u,λ,μ = B.data.l,B.data.u,B.data.λ,B.data.μ
        J = Int(block(domainspace(B),col))

        rows = blocklengths(rangespace(B.op))[1:J+l]
        cols = blocklengths(domainspace(B.op))[1:J]

        # B.data = _BandedBlockBandedMatrix(PseudoBlockArray

        blocks = PseudoBlockArray(pad(B.data.data.blocks,:,(l+u+1)*sum(cols)), (axes(B.data.data,1), blockedrange(cols)))
        B.data = _BandedBlockBandedMatrix(blocks, rows, cols, (l, u), (λ, μ))

        jr=B.datasize[2]+1:col
        kr=colstart(B.data,jr[1]):colstop(B.data,jr[end])

        isempty(kr) || BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))

        B.datasize = (last(kr),col)
    end

    B
end

function resizedata!(B::CachedOperator{T,<:BandedBlockBandedMatrix{T}},n::Integer,m::Integer) where {T<:Number}
    resizedata!(B, :, m)
    if n < B.datasize[1]
        return B
    end

    l,u,λ,μ = B.data.l,B.data.u,B.data.λ,B.data.μ

    l,u,λ,μ = B.data.l,B.data.u,B.data.λ,B.data.μ

    # make sure we have enough rows
    K = Int(block(rangespace(B),n))
    rows = blocklengths(rangespace(B.op))[1:K]

    B.data = _BandedBlockBandedMatrix(B.data.data, blockedrange(rows), (l,u), (λ,μ))
    B.datasize = (n,B.datasize[2])

    B
end


text:
function recordarena(picsfolder, cardinals, elevations, background)
    p = Progress(length(cardinals)*length(elevations), 1, "Capturing...")
    for c in eachindex(cardinals), e in eachindex(elevations)
        i = index2led(c, elevations[e])
        filename = joinpath(picsfolder, string(cardinals[c], elevations[e], background, ".png"))
        takeapicture(filename, led = i, options = `--shutter 10000000`)
        next!(p)
    end
    dots = getdots()
    lightsoff(dots)
end
function record_arena_calibration()
    picsfolder = joinpath(tempdir(), constants["calibration"]["arena"]["picsfolder"])
    if !isdir(picsfolder)
        mkdir(picsfolder)
    else
        map(rm, readdir(picsfolder, join = true))
    end

    cardinals = ["NE", "SW", "SE", "NW"]
    elevations = 1:((constants["ledstrip"]["liveleds"] - 1)÷2 + 1)
    i = round.(Int, range(23, 48, length = 3))
    elevations = elevations[i]
    # cardinals = ["NE", "SW"]
    # elevations = [36]
    recordarena(picsfolder, cardinals, elevations, "bkg")

    n = 4
    for i in 1:n
        println("background done, press when ready with sticks, take $i of $n")
        readline()
        recordarena(picsfolder, cardinals, elevations, string(i))
    end
end


text:
function emulate(op::HloRng, a, b)
    dynamic_not_implemented(op)
end

function emulate(op::HloBroadcast, arg)
    if ndims(arg) == 0
        XRTArray(fill(convert(Array, arg)[], op.result_shape...))
    else
        dynamic_not_implemented(op)
    end
end

function emulate(op::HloReshape, arg)
    XRTArray(reshape(convert(Array, arg), op.result_shape))
end

function emulate(op::HloTranspose, arg)
    aarg = convert(Array, arg)
    if ndims(aarg) != length(op.permutation)
        aarg = reshape(aarg, (size(aarg)...,
            (1 for _ in 1:(length(op.permutation)-ndims(aarg)))...))
    end
    perm = map(x->x+1, op.permutation)
    @show perm
    XRTArray(collect(permutedims(aarg, )))
end

# TODO: Fix these
function emulate(op::HloDot, a, b)
    dynamic_not_implemented(op)
end

function emulate(op, args...)
    dynamic_not_implemented(op)
end


text:
mutable struct Domain
    dt::Float64
    present_year::Float64
    initial_year::Float64
    final_year::Float64
end

Domain(dt::Float64, initial_year::Float64, final_year::Float64) = (
    Domain(dt, initial_year, initial_year, final_year)
)
Domain(dt::Float64) = Domain(dt, initial_year, final_year)


text:
@testset "domain.jl" begin
    # Closed interval domain
    a, b = 0.0, 1.0
    domain = ClosedIntervalDomain(a, b)
    @test domain.left == a
    @test domain.right == b

    N = 10
    x = discretize(domain, N)
    @test length(x) == N + 1
    @test minimum(x) == a
    @test maximum(x) == b

    # Periodic interval domain
    a, b = -1.0, 1π
    domain = PeriodicIntervalDomain(a, b)
    @test domain.left == a
    @test domain.right == b

    N = 15
    x = discretize(domain, N)
    @test length(x) == N
    @test minimum(x) == a + (b - a) / N
    @test maximum(x) == b

    # Concrete unknown domain type
    domain = UnknownDomain()
    @test_throws Exception discretize(domain, 10)
end


text:
const EXAMPLES = filter(ex -> endswith(ex, ".jl") && ex != "run_examples.jl" && ex != "all_examples.jl" && ex != "goldsteinprice.jl" && ex != "sparse_polynomials.jl" && ex != "scaled_perm.jl" && ex != "symmetry.jl" && ex != "discourse_65377.jl" && ex != "discourse_65377_1.jl" && ex != "prajna.jl",
                        readdir(@__DIR__))

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# A hack to make it easier to pass in an array payload.
struct ArrayBody{T}
    xs::T
end

JSON3.write(kw::Pairs{Symbol, ArrayBody{T}, Tuple{Symbol}, NamedTuple{(:array,), Tuple{ArrayBody{T}}}}) where T =
    JSON3.write(kw.data.array.xs)

function api_call(c, method, path, Into=Nothing, params=Dict();  kwargs...)
    @debug "$method $path"

    headers = [
        "Authorization" => auth_header(c),
        "User-Agent" => USER_AGENT,
        "X-RateLimit-Precision" => "millisecond",
    ]

    body, query = if method in (:PATCH, :POST, :PUT)
        if haskey(kwargs, :file)
            kw_dict = Dict(kwargs)
            file = pop!(kw_dict, :file)
            # This is just a hack to allow for easier testing. No user should use this.
            boundary = NamedTuple()
            if haskey(kw_dict, :__boundary__)
                boundary = (; boundary=pop!(kw_dict, :__boundary__))
            end
            form = Form(
                Dict(:file => file, :payload_json => JSON3.write(kw_dict));
                boundary...,
            )
            form, params
        else
            push!(headers, "Content-Type" => "application/json")
            JSON3.write(kwargs), params
        end
    else
        "", kwargs
    end

    url = "$API_BASE/v$API_VERSION$path"
    if check_rate_limits(rate_limiter(c), path) === RATE_LIMIT_SENTINEL
        return RATE_LIMIT_SENTINEL
    end

    resp = request(method, url, headers, body; query=query, status_exception=false)
    rl_result = apply_rate_limits!(rate_limiter(c), resp)
    if rl_result === RATE_LIMIT_RETRY
        return api_call(c, method, path, Into, params; kwargs...)
    elseif rl_result === RATE_LIMIT_SENTINEL
        return RATE_LIMIT_SENTINEL
    end

    if 200 <= resp.status < 300
        return parse_response(resp, Into)
    else
        throw(StatusError(resp.status, resp))
    end
end

parse_response(resp::Response, ::Type{Nothing}) = nothing
function parse_response(resp::Response, Into)
    return if resp.status == 204
        nothing
    elseif header(resp, "Content-Type") == "application/json"
        x = JSON3.read(resp.body, Into)
        fix_datetimes!(x)
        x
    else
        resp.body
    end
end

# I hate this.
parse_datetime(s) = DateTime(replace(s, "+" => ".000+")[1:23], ISODateTimeFormat)
fix_datetimes!(xs::Vector) = foreach(fix_datetimes!, xs)
function fix_datetimes!(x::T) where T
    for name in fieldnames(T)
        if DateTime <: fieldtype(T, name)
            f = getfield(x, name)
            if f isa String
                setfield!(x, name, parse_datetime(f))
            elseif f isa Int
                d = unix2datetime(f)
                if year(d) == 1970
                    d = unix2datetime(1000f)
                end
                setfield!(x, name, d)
            end
        end
    end
end

path_argname(s) = s
path_argname(s::Symbol) = string("<", s === :emoji_id ? :emoji : s, ">")
path_argname(ex::Expr) = path_argname(ex.args[1])

# TODO: Refactor this to be not huge and ugly.
macro route(name, method, path, kwargs...)
    if path isa String
        path = Expr(:string, path)
    end

    # Copy all the interpolations as arguments.
    fun_args = [x isa Expr ? Expr(:kw, x.args...) : x for x in path.args if !(x isa String)]
    # We'll use this later, but construct it now while `path.args` still makes sense.
    doc_path = join(map(path_argname, path.args))

    # Update the string interpolation arguments to URI-encode them, among some other hacks.
    for (i, arg) in enumerate(path.args)
        done = false
        if !(arg isa String)
            arg = if arg isa Expr && arg.head === :(=)
                k, v = arg.args
                if v === :nothing && i == lastindex(path.args)
                    # When the argument is omitted, remove the trailing slash.
                    done = true
                    path.args[i-1] = path.args[i-1][1:end-1]
                    insert!(path.args, i, :($k === nothing ? "" : "/"))
                    i += 1
                    :($k === nothing ? "" : $k)
                else
                    # Remove the default value from the interpolation.
                    k
                end
            elseif arg === :emoji_id
                # Most endpoints that accept emojis use the name,
                # so escapeuri is implemented to encode the name.
                # For the few endpoints that require the ID, this special case exists.
                fun_args[findfirst(==(:emoji_id), fun_args)] = :emoji
                :(emoji isa Emoji ? emoji.id : emoji)
            else
                arg
            end
            path.args[i] = :(escapeuri($arg))
        end
        done && break
    end

    fun = :($name(client, $(fun_args...)) = api_call(client, $(QuoteNode(method)), $path))
    pushfirst!(fun.args[2].args, __source__)

    sig_args = fun.args[1].args
    call_args = fun.args[2].args[end].args
    insert!(sig_args, 2, Expr(:parameters))
    insert!(call_args, 2, Expr(:parameters))
    sig_kws = sig_args[2].args
    call_kws = call_args[2].args
    Into = :Nothing
    for ex in kwargs
        if ex === :kwargs
            # Add trailing keywords that go into the query parameters or body.
            push!(sig_kws, :(kwargs...))
            push!(call_kws, :(kwargs...))
        elseif ex isa Expr && ex.head === :(=)
            k, v = ex.args
            if k === :array
                # Pass one specific keyword as an array body.
                pushfirst!(sig_kws, v)
                pushfirst!(call_kws, Expr(:kw, :array, :(ArrayBody($v))))
            elseif k === :query
                # Add query parameters to a request that uses its keywords for the body.
                tuples = [Expr(:tuple, QuoteNode(t.args[1]), t.args[1]) for t in v.args]
                push!(call_args, Expr(:tuple, tuples...))
                for t in v.args
                    pushfirst!(sig_kws, Expr(:kw, t.args...))
                end
            end
        else
            Into = ex
        end
    end
    insert!(call_args, 6, Into)
    isempty(sig_kws) && deleteat!(sig_args, 2)
    isempty(call_kws) && deleteat!(call_args, 2)

    doc_call = replace(string(Expr(:call, sig_args...)), " = " => "=")
    doc = """
        $doc_call -> $Into

    Make a $method request to `$doc_path`.
    See [the Discord API documentation](https://discord.com/developers/docs/resources/$(RESOURCE[]))
    for more information.
    """

    block = quote
        export $name
        $fun
        @doc $doc $name
    end
    return esc(block)
end

const HasID = Union{Guild, DiscordChannel, User, Message, Overwrite, Role, Webhook}
HTTP.escapeuri(x::HasID) = string(x.id)
HTTP.escapeuri(e::Emoji) = escapeuri(e.name)
HTTP.escapeuri(i::Invite) = escapeuri(i.code)

const RESOURCE = Ref{String}()

RESOURCE[] = "audit-log"
@route get_guild_audit_log GET "/guilds/$guild/audit-logs" AuditLog kwargs

RESOURCE[] = "channel"
@route get_channel GET "/channels/$channel" DiscordChannel kwargs
@route update_channel PATCH "/channels/$channel" DiscordChannel kwargs
@route delete_channel DELETE "/channels/$channel" DiscordChannel
@route get_channel_messages GET "/channels/$channel/messages" Vector{Message} kwargs
@route get_channel_message GET "/channels/$channel/messages/$message" Message
@route create_message POST "/channels/$channel/messages" Message kwargs
@route create_reaction PUT "/channels/$channel/messages/$message/reactions/$emoji/@me"
@route delete_reaction DELETE "/channels/$channel/messages/$message/reactions/$emoji/$(user="@me")"
@route get_reactions GET "/channels/$channel/messages/$message/reactions/$emoji" Vector{User} kwargs
@route delete_all_reactions DELETE "/channels/$channel/messages/$message/reactions/$(emoji=nothing)"
@route update_message PATCH "/channels/$channel/messages/$message" Message kwargs
@route delete_message DELETE "/channels/$channel/messages/$message"
@route delete_messages POST "/channels/$channel/messages/bulk-delete" kwargs
@route update_channel_permissions PUT "/channels/$channel/permissions/$overwrite" kwargs
@route get_channel_invites GET "/channels/$channel/invites" Vector{Invite} kwargs
@route create_channel_invite POST "/channels/$channel/invites" Invite kwargs
@route delete_channel_permission DELETE "/channels/$channel/permissions/$overwrite"
@route create_typing_indicator POST "/channels/$channel/typing"
@route get_pinned_messages GET "/channels/$channel/pins" Vector{Message}
@route create_pinned_channel_message PUT "/channels/$channel/pins/$message"
@route delete_pinned_channel_message DELETE "/channels/$channel/pins/$message"
@route create_group_dm_recipient PUT "/channels/$channel/recipients/$user" kwargs
@route delete_group_dm_recipient DELETE "/channels/$channel/recipients/$user"

RESOURCE[] = "emoji"
@route get_guild_emojis GET "/guilds/$guild/emojis" Vector{Emoji}
@route get_guild_emoji GET "/guilds/$guild/emojis/$emoji_id" Emoji
@route create_guild_emoji POST "/guilds/$guild/emojis" Emoji kwargs
@route update_guild_emoji PATCH "/guilds/$guild/emojis" Emoji kwargs
@route delete_guild_emoji DELETE "/guilds/$guild/emojis/$emoji_id"

RESOURCE[] = "guild"
@route create_guild POST "/guilds" Guild kwargs
@route get_guild GET "/guilds/$guild" Guild kwargs
@route get_guild_preview GET "/guilds/$guild/preview" Guild
@route update_guild PATCH "/guilds/$guild" Guild kwargs
@route delete_guild DELETE "/guilds/$guild"
@route get_guild_channels GET "/guilds/$guild/channels" Vector{DiscordChannel}
@route create_guild_channel POST "/guilds/$guild/channels" DiscordChannel kwargs
@route update_guild_channel_positions PATCH "/guilds/$guild/channels" array=positions
@route get_guild_member GET "/guilds/$guild/members/$user" GuildMember
@route get_guild_members GET "/guilds/$guild/members" Vector{GuildMember} kwargs
@route create_guild_member PUT "/guilds/$guild/members/$user" GuildMember kwargs
@route update_guild_member PATCH "/guilds/$guild/members/$user" kwargs
@route modify_user_nick PATCH "/guilds/$guild/members/@me/nick" UserNickChange kwargs
@route create_guild_member_role PUT "/guilds/$guild/members/$user/roles/$role"
@route delete_guild_member_role DELETE "/guilds/$guild/members/$user/roles/$role"
@route delete_guild_member DELETE "/guilds/$guild/members/$user"
@route get_guild_bans GET "/guilds/$guild/bans" Vector{Guild}
@route get_guild_ban GET "/guilds/$guild/bans/$user" Ban
@route create_guild_ban PUT "/guilds/$guild/bans/$user" kwargs
@route delete_guild_ban DELETE "/guilds/$guild/bans/$user"
@route get_guild_roles GET "/guilds/$guild/roles" Vector{Role}
@route create_guild_role POST "/guilds/$guild/roles" Role kwargs
@route update_guild_role_positions PATCH "/guilds/$guild/roles" Vector{Role} array=positions
@route update_guild_role PATCH "/guilds/$guild/roles/$role" Role kwargs
@route delete_guild_role DELETE "/guilds/$guild/roles/$role"
@route get_guild_prune_count GET "/guilds/$guild/prune" PruneCount kwargs
@route create_guild_prunt POST "/guilds/$guild/prune" PruneCount kwargs
@route get_guild_voice_regions GET "/guilds/$guild/regions" Vector{VoiceRegion}
@route get_guild_invites GET "/guilds/$guild/invites" Vector{Invite}
@route get_guild_integrations GET "/guilds/$guild/integrations" Vector{Integration}
@route create_guild_integration POST "/guilds/$guild/integrations" kwargs
@route update_guild_integration PATCH "/guilds/$guild/integrations/$integration" kwargs
@route delete_guild_integration PATCH "/guilds/$guild/integrations/$integration"
@route sync_guild_integration POST "/guilds/$guild/integrations/$integration/sync"
@route get_guild_widget GET "/guilds/$guild/widget" GuildWidget
@route update_guild_widget PATCH "/guilds/$guild/widget" GuildWidget kwargs
@route get_guild_vanity_url GET "/guilds/$guild/vanity-url" Invite
@route get_guild_widget_image GET "/guilds/$guild/widget.png" String

RESOURCE[] = "invite"
@route get_invite GET "/invites/$invite" Invite kwargs
@route delete_invite DELETE "/invites/$invite" Invite

RESOURCE[] = "user"
@route get_user GET "/users/$(user="@me")" User
@route update_user PATCH "/users/@me" User kwargs
@route get_user_guilds GET "/users/@me/guilds" Vector{Guild} kwargs
@route leave_guild DELETE "/users/@me/guilds/$guild"
@route get_user_dms GET "/users/@me/channels" Vector{DiscordChannel}
@route create_dm POST "/users/@me/channels" DiscordChannel kwargs
@route create_group_dm POST "/users/@me/channels" DiscordChannel kwargs
@route get_user_connections GET "/users/@me/connections" Vector{Connection}

RESOURCE[] = "voice"
@route get_voice_regions GET "/voice/regions" Vector{VoiceRegion}

RESOURCE[] = "webhook"
@route create_webhook POST "/channels/$channel/webhooks" Webhook kwargs
@route get_channel_webhooks GET "/channels/$channel/webhooks" Vector{Webhook}
@route get_guild_webhooks GET "/guilds/$guild/webhooks" Vector{Webhook}
@route get_webhook GET "/webhooks/$webhook/$(token=nothing)" Webhook
@route update_webhook PATCH "/webhooks/$webhook/$(token=nothing)" Webhook kwargs
@route delete_webhook DELETE "/webhooks/$webhook/$(token=nothing)" Webhook kwargs
@route execute_webhook POST "/webhooks/$webhook/$token" Message query=(wait=true,) kwargs
@route execute_webhook_github POST "/webhooks/$webhook/$token/github" Message query=(wait=true,) kwargs
@route execute_webhook_slack POST "/webhooks/$webhook/$token/slack" Message query=(wait=true,) kwargs


text:
const image = let raw = read("input")
    last(raw) === 0x0a && pop!(raw)  # remove trailing newline
    raw
end

function part1(image)
    layers = collect(Iterators.partition(image, 25 * 6))
    i = argmin(map(x->count(==(0x30), x), layers))
    println("Part 1 answer: ", count(==(0x31), layers[i]) * count(==(0x32), layers[i]))
end

function firstcolor(slice::AbstractVector{UInt8})
    for i = 1:length(slice)
        x = slice[i]
        if x === 0x30 || x === 0x31
            return x
        end
    end
    0x32
end

function part2(image)
    # Do this little dance to get a row-major `reshape`
    X = permutedims(reshape(image, 25, 6, :), (2, 1, 3))
    println("Part 2 answer: ")
    for i = 1:6
        for j = 1:25
            c = firstcolor(view(X, i, j, :))
            print(c === 0x30 ? ' ' : '#')
        end
        println()
    end
end

part1(image)
part2(image)

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
"""
Wraps an optimizer in a hierarchy which enables one `bridge`.
"""
function lazy_relax(optimizer, bridge::Type{<:MOIB.Constraint.AbstractBridge})
        function construct()
                model = Model(optimizer)
                JuMP.add_bridge(model, bridge)
                backend(model).optimizer
        end
        construct
end


text:
# Compare Jain-Neal sampler to Reversible Jump on the galaxy dataset.
module GalaxyCompareJNtoRJ

using BayesianMixtures
B = BayesianMixtures
can_plot = (Pkg.installed("PyPlot")!=nothing)
can_plot ? using PyPlot : warn("Skipping plots since PyPlot is not installed.")

# ==================== Run Jain-Neal sampler ====================

# Read the data from file
x = readdlm("datasets/galaxy.dat",' ',Float64)[:]

# Specify model, data, and MCMC options
n_total = 1000000  # total number of MCMC sweeps to run
log_pk = "k -> log(k in (1:30)? 1/30 : 0)"   # log prior on k is Uniform{1,...,30}
options = B.options("Normal","MFM",x,n_total; n_keep=1000,log_pk=log_pk,t_max=30)

# Run MCMC sampler
result = B.run_sampler(options)

# ==================== Run RJMCMC sampler using Nmix ====================

# Remove old output files
if (!isdir("galx")); mkdir("galx"); end
outfiles = ["1.bdlog","1.bk","1.ent","1.k","1.log","1.out","1.pe"]
for outfile in outfiles; try rm(joinpath("galx",outfile)); end; end

# Run Peter Green's Nmix program
if is_apple() # if running Mac OSX
    run(`chmod +x Nmix`)
    tic()
    run(`./Nmix -n$n_total -nb0 -ns$n_total galx`)
    elapsed_time_RJ = toq()
elseif is_windows() # if running Windows
    tic()
    run(`Nmix.exe -n$n_total -nb0 -ns$n_total galx`)
    elapsed_time_RJ = toq()
else
    error("Only Mac and Windows are currently supported for running Green's Nmix program.")
end

# Read and delete the output
output = readdlm(joinpath("galx","1.out"))
for outfile in outfiles; rm(joinpath("galx",outfile)); end

# Parse the output
k_RJ = zeros(Int64,n_total)
t_RJ = zeros(Int64,n_total)
i = 1
for j = 1:n_total
    k,t = output[i,find(output[i,:].!="")]
    k_RJ[j] = k
    t_RJ[j] = t
    i += k+1
end

# ==================== Compare posterior on k ====================

# Compute the posterior on k from Jain-Neal sampler
k_posterior_JN = B.k_posterior(result)

# Compute posterior on k from RJMCMC sampler
use = options.n_burn+1:n_total
k_posterior_RJ = B.histogram(k_RJ[use], 0:options.t_max)[1]/length(use)

# Print posteriors on k
println("Posteriors on k:")
println("  Jain-Neal  RJMCMC")
for k=1:15
    @printf "k=%d: %.4f %.4f\n" k k_posterior_JN[k] k_posterior_RJ[k]
end

if can_plot

    # ==================== Compare traceplots of t ====================
    # Traceplot of t for Jain-Neal
    B.open_figure(1)
    B.traceplot(result.t[1:min(5000,n_total)])
    B.labels("Traceplot for t (Jain-Neal)","sweep number","t")

    # Traceplot of t for RJMCMC
    B.open_figure(2)
    B.traceplot(t_RJ[1:min(5000,n_total)])
    B.labels("Traceplot for t (RJMCMC)","sweep number","t")

    # ==================== Compare autocorrelation of t ====================
    # Plot autocorrelation for Jain-Neal
    B.open_figure(3)
    B.plot_autocorrelation(result,0.01)

    # Plot autocorrelation for RJMCMC
    B.plot_autocorrelation(t_RJ[use],0.01,elapsed_time_RJ/n_total; color="r")
    legend(["Jain-Neal","RJMCMC"],loc="upper right",fontsize=12)

    # ==================== Plot density estimate ====================
    # Plot density estimate using Jain-Neal output
    B.open_figure(4)
    B.plot_histogram(x; color="w",edgecolor="g",linewidth=0.5)
    B.rug_plot(x)
    B.plot_density_estimate(result; resolution=500)
end

end # module GalaxyCompareJNtoRJ


text:
rootpath = "/home/user/music"
pattern  = r".mp3$"

for (root, dirs, files) in walkdir(rootpath)
    for file in files
        if occursin(pattern, file) println(file) end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
push!(LOAD_PATH,"../src/")
using Documenter, AlphaStructures

makedocs(
	format = Documenter.HTML(
		prettyurls = get(ENV, "CI", nothing) == "true"
	),
	sitename = "AlphaStructures.jl",
	pages = [
		"1 - Home" => "index.md",
		"2 - Getting Started" => "gettingStarted.md",
		"Theory..." => [
			"3.0 - Theory Index" => "theory-index.md",
			"3.1 - Delaunay Triangulation" => "delaunay.md",
			"3.2 - Voronoi Diagrams" => "voronoi.md",
			"3.3 - Alpha Structures" => "alpha-structures.md",
		],
		"... and Practice" => [
			"4.0 - Module Introduction" => "this-module.md",
			"4.1 - DeWall Algorithm" => "delaunay-impl.md",
			"4.2 - Alpha Structures" => "alpha-structures-impl.md",
		],
		"A - About the Authors" => "authors.md",
		"B - Bibliography" => "bibliography.md"
	]
)

deploydocs(
    repo = "github.com/eOnofri04/AlphaStructures.jl.git",
)

# makedocs(
# 	format = :html,
# 	sitename = "AlphaStructures.jl",
# 	assets = ["assets/AlphaStructures.css", "assets/logo.png"],
# 	pages = [
# 		"1 - Home" => "index.md",
# 		"2 - Getting Started" => "gettingStarted.md",
# 		"Theory..." => [
# 			"3.0 - Theory Index" => "theory-index.md",
# 			"3.1 - Delaunay Triangulation" => "delaunay.md",
# 			"3.2 - Voronoi Diagrams" => "voronoi.md",
# 			"3.3 - Alpha Structures" => "alpha-structures.md",
# 			"3.4 - Persistent Homology" => "persistent-homology.md",
# 		],
# 		"... and Practice" => [
# 			"4.0 - Module Introduction" => "this-module.md",
# 			"4.1 - DeWall Algorithm" => "delaunay-impl.md",
# 			"4.2 - Alpha Structures" => "alpha-structures-impl.md",
# 			"4.3 - Persistent Homology" => "persistent-homology-impl.md",
# 		],
# 		"A - About the Authors" => "authors.md",
# 		"B - Bibliography" => "bibliography.md"
# 	]
# )
#

#
# makedocs(
#     modules = [AlphaStructures],
#     format = Documenter.HTML(
#         prettyurls = "deploy" in ARGS,
#     ),
#     sitename = "AlphaStructures.jl",
# 	pages = [
# 		"1 - Home" => "index.md",
# 		"2 - Getting Started" => "gettingStarted.md",
# 		"Theory..." => [
# 			"3.0 - Theory Index" => "theory-index.md",
# 			"3.1 - Delaunay Triangulation" => "delaunay.md",
# 			"3.2 - Voronoi Diagrams" => "voronoi.md",
# 			"3.3 - Alpha Structures" => "alpha-structures.md",
# 			"3.4 - Persistent Homology" => "persistent-homology.md",
# 		],
# 		"... and Practice" => [
# 			"4.0 - Module Introduction" => "this-module.md",
# 			"4.1 - DeWall Algorithm" => "delaunay-impl.md",
# 			"4.2 - Alpha Structures" => "alpha-structures-impl.md",
# 			"4.3 - Persistent Homology" => "persistent-homology-impl.md",
# 		],
# 		"A - About the Authors" => "authors.md",
# 		"B - Bibliography" => "bibliography.md"
# 	]
# )
#
# deploydocs(
#     repo = "github.com/eOnofri04/AlphaStructures.jl.git",
#     push_preview = true,
# )


text:
using QuantumOptics, CollectiveSpins

# System parameters
const γ = 1.
const e_dipole = [0,0,1.]
const T = [0:0.05:5]

const system = SpinCollection(geometry.chain(0.5, 8), e_dipole, γ)
const N = length(system.spins)

# Initial state (Bloch state)
const phi = 0.
const theta = pi/2.

const sx0 = cos(phi)*sin(theta)
const sy0 = sin(phi)*sin(theta)
const sz0 = cos(theta)

# Quantum: master equation
Ψ₀ = CollectiveSpins.quantum.blochstate(phi,theta,N)
ρ₀ = Ψ₀⊗dagger(Ψ₀)
@time tout, ρ_t = CollectiveSpins.quantum.timeevolution(T, system, ρ₀)

# Meanfield
state0 = CollectiveSpins.meanfield.blochstate(phi,theta,N)
@time tout, state_mf_t = CollectiveSpins.meanfield.timeevolution(T, system, state0)
@time tout, state_mf_t = CollectiveSpins.meanfield.timeevolution(T, system, state0)

# Meanfield + Correlations
state0 = CollectiveSpins.mpc.blochstate(phi,theta,N)
@time tout, state_cor_t = CollectiveSpins.mpc.timeevolution(T, system, state0)
@time tout, state_cor_t = CollectiveSpins.mpc.timeevolution(T, system, state0)

# Expectation values
sx_indep = exp(-0.5*γ*T)*sx0
sy_indep = exp(-0.5*γ*T)*sy0
sz_indep = 1 - exp(-γ*T)*(1-sz0)

sx_mf = map(s->(CollectiveSpins.meanfield.sx(s)[1]), state_mf_t)
sy_mf = map(s->(CollectiveSpins.meanfield.sy(s)[1]), state_mf_t)
sz_mf = map(s->(CollectiveSpins.meanfield.sz(s)[1]), state_mf_t)

sx_cor = map(s->(CollectiveSpins.mpc.sx(s)[1]), state_cor_t)
sy_cor = map(s->(CollectiveSpins.mpc.sy(s)[1]), state_cor_t)
sz_cor = map(s->(CollectiveSpins.mpc.sz(s)[1]), state_cor_t)

# Trace distances
ρmf_t = map(CollectiveSpins.meanfield.densityoperator, state_mf_t)
ρcor_t = map(CollectiveSpins.mpc.densityoperator, state_cor_t)
td_mf = [tracedistance(ρ_t[i],ρmf_t[i]) for i=1:length(T)]
td_cor = [tracedistance(ρ_t[i],ρcor_t[i]) for i=1:length(T)]

# rhop_t = ρ_t
# for i=N:-1:2
#     rhop_t = map(x->ptrace(x, [i]), rhop_t)
# end

embed(op::Operator) = QuantumOptics.embed(CollectiveSpins.quantum.basis(system), [1], [op])
sx_master = map(ρ->expect(embed(sigmax), ρ), ρ_t)
sy_master = map(ρ->expect(embed(sigmay), ρ), ρ_t)
sz_master = map(ρ->expect(embed(sigmaz), ρ), ρ_t)

# Visualization
using PyCall
@pyimport matplotlib.pyplot as plt

#name = "cube_atoms=8_d=22925"
#name = "cube_atoms=8_d=0627"
name = "cube_atoms=8_d=015"
# name = "chain_atoms=8_d=07"

plt.figure(figsize=(6,4))
plt.plot(T, -sz_master, lw=3.0, color="darkorange", label="Collective")
plt.plot(T, -sz_indep, lw=3.0, color="gray", label="Independent")
#plt.plot(T_smf, -sz_smf, label="symmetric meanfield")
plt.plot(T, -sz_mf, "-", lw=1.5, color="navy", label="Mean-field")
plt.plot(T, -sz_cor, "-", lw=1.5, color="green", label="Correlations")
plt.ylabel("\$-\\langle\\sigma_z\\rangle\$")
plt.ylim(-1,1)
plt.legend()
#plt.savefig("images/$(name)_sigmaz_timeevolution.pdf")

plt.figure(figsize=(6,4))
plt.plot(T, sx_master, lw=3.0, color="darkorange", label="Collective")
plt.plot(T, sx_indep, lw=3.0, color="gray", label="Independent")
#plt.plot(T_smf, sx_smf, label="symmetric meanfield")
plt.plot(T, sx_mf, "-", lw=1.5, color="navy", label="Mean-field")
plt.plot(T, sx_cor, "-", lw=1.5, color="green", label="Correlations")
plt.ylabel("\$\\langle\\sigma_x\\rangle\$")
plt.ylim(-1,1)
plt.legend()
#plt.savefig("images/$(name)_sigmax_timeevolution.pdf")

plt.figure(figsize=(6,4))
plt.plot(T, sy_master, lw=3.0, color="darkorange", label="Collective")
plt.plot(T, sy_indep, lw=3.0, color="gray", label="Independent")
#plt.plot(T_smf, sy_smf, label="symmetric meanfield")
plt.plot(T, sy_mf, "-", lw=1.5, color="navy", label="Mean-field")
plt.plot(T, sy_cor, "-", lw=1.5, color="green", label="Correlations")
plt.ylabel("\$\\langle\\sigma_y\\rangle\$")
plt.ylim(-1,1)
plt.legend()
#plt.savefig("images/$(name)_sigmay_timeevolution.pdf")

plt.figure(figsize=(6,4))
plt.plot(T, td_mf, lw=1.2, color="navy", label="Mean-field")
plt.plot(T, td_cor, lw=1.2, color="green", label="Correlation")
plt.xlabel("Time [\$\\gamma^{-1}\$]")
plt.ylabel("Error (Trace distance)")
plt.ylim(0,1)
plt.legend()
# plt.savefig("images/$(name)_error_timeevolution.pdf")

plt.show()
#println("⦑H⟩ = ", trunc(real(exp), 3))
#println("⦑H⟩ = ", trunc(real(exp_mean), 3))


text:
#Parameters of Neutral Diabatic Matrix Element
global const A_n = 457000.052788*kjpermol #Au--O: exponential repulsion: A
global const α_n = 3.75257871821/Å #Au--O: exponential repulsion: alpha
global const Au_O_cutoff_n= 10.0 * Å #Au--O: cutoff
global const B_n = 30788.8486039*kjpermol #Au--N: exponential repulsion: B
global const β_n = 2.9728905911/Å #Au--N: exponential repulsion: beta
global const Au_N_cutoff_n = 10.0  * Å #Au--N: cutoff
global const F_n = 638.5000000000*kjpermol #N--O: Morse: F
global const δ_n=2.7430000000/Å #N--O: Morse: delta
global const r_0_NO = 1.1507700000*Å #N--O: Morse: cutoff

#Parameters of Neutral Diabatic Matrix Element

global const C_i = 1.25581276843 * Å  # Image potential: C
global const D_i = 347.2225355 * kjpermol * Å  #Image potential: D
global const z0_i = 1.153606314 * Å  #Image potential: z0
global const A_i = 457000.052788 * kjpermol  #Au--O: exponential repulsion: A
global const α_i = 3.75257871821 / Å  #Au--O: exponential repulsion: alpha
global const Au_O_cutoff_i = 10 / Å #Au--O: cutoff
global const B_i = 23.8597594272 * kjpermol  #Au--N: Morse: B
global const β_i=1.91014033785 / Å  #Au--N: Morse: beta
global const rN_sur_e_i = 2.38958832878 * Å  #Au--N: Morse: rN_sur_e
global const Au_N_cutoff_i = 10 * Å  #Au--N: cutoff
global const F_i = 495.9809807 * kjpermol   #N--O: Morse: F
global const δ_i = 2.47093477934 / Å  #N--O: Morse: delta
global const rNO_e_i = 1.29289837288 / Å  #N--O: Morse: cutoff
global const KI = 512.06425722 * kjpermol  #N--O: Morse: cutoff

# DEFINE PARAMETERS FOR COUPLING FUNCTION

global const coup_a_N=-70.5259924491 * kjpermol  # Au--N: Exponential decay: a
global const coup_b_N=0.00470023958504 # Au--N: Exponential decay: b
global const coup_β_N=1.95982478112 / Å # Au--N: Exponential decay: beta
global const Au_N_coupling_cutoff=10 * Å # Au--N: Exponential decay: cutoff
global const coup_a_O=-16.7488672932 * kjpermol # Au--O: Exponential decay: a
global const coup_b_O=0.00617151653727 # Au--O: Exponential decay: b
global const coup_β_O=1.35353579356 / Å  # Au--O: Exponential decay: beta
global const Au_O_coupling_cutoff=10 * Å # Au--O: Exponential decay: cutoff

#Define parameters for Au-Au interaction potential
global const α = -4.94
global const β = 17.15
global const γ = 19.40

function def_d_matrices(α, β, γ)
    #see paper from 1947 for definitions
    d17 = [[α 0.0 0.0]; [0.0 β γ]; [0.0 γ β]]
    d28 = [[α 0.0 0.0]; [0.0 β -γ]; [0.0 -γ β]]
    d39 = [[β 0.0 γ]; [0.0 α 0.0]; [γ 0.0 β]]
    d410 = [[β 0.0 -γ]; [0.0 α 0.0]; [-γ 0.0 β]]
    d511 = [[β γ 0.0]; [γ β 0.0]; [0.0 0.0 α]]
    d612 = [[β -γ 0.0]; [-γ β 0.0]; [0.0 0.0 α]]

    return d17, d28, d39, d410, d511, d612
end
const global d17, d28, d39, d410, d511, d612 = def_d_matrices(α, β, γ)

function compute_d_new_basis(d17, d28, d39, d410, d511, d612)
    U = permutedims([[-1.0 0.0 1.0]/sqrt(2.0);[1.0 -2.0 1.0]/sqrt(6.0);[-1.0 -1.0 -1.0]/sqrt(3.0)])
    d1_new = U' * d17
    d2_new = U' * d28
    d3_new = U' * d39
    d4_new = U' * d410
    d5_new = U' * d511
    d6_new = U' * d612
    return d1_new, d2_new, d3_new, d4_new, d5_new, d6_new
end
const global d1_new, d2_new, d3_new, d4_new, d5_new, d6_new = compute_d_new_basis(d17, d28, d39, d410, d511, d612)

#set masses
const global mass_arr0 =  [14.00307440, 15.99491502, 196.966548]
const global mass_arr =  [14.00307440, 15.99491502, 196.966548]*amu


text:
"""
Operations on probability distributions.

Wouter Kouw
28-03-2029
"""

export Delta

import Base: *
import Statistics: mean, var
import Distributions: params, pdf, logpdf
# import Distributions: ContinuousUnivariateDistribution

struct Delta <: ContinuousUnivariateDistribution
    """
    Delta distribution.

    This is not a proper distribution, but this structure allows for a general
    interface for observed variables and unobserved variables (i.e. it is
    possible to call mean/var on a fixed variable).
    """

    # Attributes
    value::Float64

    function Delta(value::Float64)
        return new(value);
    end
end

function mean(d::Delta)
    return d.value
end

function std(d::Delta)
    return 0.0
end

function var(d::Delta)
    return 0.0
end

function meanvar(d::Delta)
    return mean(d), var(d)
end

struct Flat <: ContinuousUnivariateDistribution
    """
    Non-informative (flat) prior distribution.

    This is an improper prior that treats each element of its sample space as
    equally likely. It is used as an identity function for multiplication with
    other probability distrbutions. In other words:
        Flat(x) ⋅ Normal(x | 0,1) = Normal(x | 0,1).

    # TODO: Flat can be replaced by Uniform(-Inf, Inf)
    """
    Flat() = new()
end

function pdf(p::Flat, x::Float64)
    "Assuming sample space of Flat is [-Inf, Inf]"
    return 0.0
end

function logpdf(p::Flat, x::Float64)
    return log(pdf(p, x))
end

function mean(d::Flat)
    return 0.0
end

function var(d::Flat)
    return Inf
end

function meanvar(p::Flat)
    return mean(p), var(p)
end

function meanvar(p::Normal)
    "Mean and variance of a Normal distribution."
    return mean(p), var(p)
end

# function pdf(p::Normal, x::Float64)
#     "Probability of sample under Normal distribution."
#     return 1/sqrt(2*π*var(p))*exp(-(x-mean(p))^2/(2*var(p)))
# end
#
# function logpdf(p::Normal, x::Float64)
#     "Probability of sample under Normal distribution."
#     return -log(2*π)/2 - log(std(p)) -(x-mean(p))^2/(2*var(p))
# end

function meanvar(p::Gamma)
    "Mean and variance of a Gamma distribution."
    return mean(p), var(p)
end

function meanvar(graph::MetaGraph, x::Float64)
    "If variable is clamped to a value, then return that value with variance 0."
    return x, 0
end

function meanvar(graph::MetaGraph, x::Symbol)
    "Extract marginal from variable in graph and return mean and variance."
    return meanvar(graph[graph[x, :id], :node].marginal)
end

function *(a::Symbol, b::Integer)
    "Concatenate symbol and integer"
    return Symbol(string(a)*string(b))
end

function *(a::Symbol, b::String)
    "Concatenate symbol and integer"
    return Symbol(string(a)*b)
end

function *(px::Normal{Float64}, qx::Normal{Float64})
    "Multiplication of two normal distributions of the same variable."

    # Extract parameters
    μ_p, σ_p = params(px)
    μ_q, σ_q = params(qx)

    # Compute precisions
    τ_p = inv(σ_p^2)
    τ_q = inv(σ_q^2)

    # Add precisions
    τ = τ_p + τ_q

    # Add precision-weighted means
    τμ = τ_p*μ_p + τ_q*μ_q

    # Compute new variance
    σ2 = inv(τ)

    # Return normal distribution in mean-stddev parameterization
    return Normal(σ2*τμ, sqrt(σ2))
end

function *(px::Gamma{Float64}, qx::Gamma{Float64})
    "Multiplication of two Gamma distributions of the same variable."

    # Extract parameters
    shape_p, scale_p = params(px)
    shape_q, scale_q = params(qx)

    # Add shapes
    new_shape = shape_p + shape_q - 1.

    # Add inverse scales
    new_scale = inv(inv(scale_p) + inv(scale_q))

    # Return gamma distribution in shape-scale parameterization
    return Gamma(new_shape, new_scale)
end

function *(px::Flat, qx::Flat)
    return Flat()
end

function *(px::Union{Gamma{Float64},Normal{Float64}}, qx::Flat)
    "Multiplication of a distribution with a flat distribution."
    return px
end

function *(qx::Flat, px::Union{Gamma{Float64},Normal{Float64}})
    "Multiplication of a distribution with a flat distribution."
    return px
end

function expectation(q::Normal{Float64}, ν::Normal{Float64})
    "Compute the expectation E_q(x) [log ν(x)]"

    # Parameters
    μ_q, σ_q = params(q)
    μ_ν, σ_ν = params(ν)

    return -log(σ_ν) - 1/2*log(2*π) -1/(2*σ_ν^2)*(σ_q^2 + μ_q^2 -2*μ_q*μ_ν + μ_ν^2)
end

function expectation(q::Gamma{Float64}, ν::Gamma{Float64})
    "Compute the expectation E_q(x) [log ν(x)]"

    # Parameters
    α_q, θ_q = params(q)
    α_ν, θ_ν = params(ν)

    return -log(gamma(α_ν)) - α_ν*log(θ_ν) + (α_ν-1)*(digamma(α_q) + log(θ_q)) - (α_q*θ_q)/θ_ν
end

function expectation(q::Union{Gamma{Float64},Normal{Float64}}, ν::Flat)
    "If the message is a flat distribution, the integral E_q(x) [log ν(x)] diverges."
    return Inf
end

function expectation(ν::Flat, q::Union{Gamma{Float64},Normal{Float64}})
    "If the message is a flat distribution, the integral E_q(x) [log ν(x)] diverges."
    return Inf
end

function grad_expectation(q::Normal{Float64}, ν::Normal{Float64})
    "Compute the gradient of the expectation E_q(x) [log ν(x)] w.r.t. q(x)'s params"

    # Parameters
    μ_q, σ_q = params(q)
    μ_ν, σ_ν = params(ν)

    # Partial derivative w.r.t μ_q
    partial_μ = -(μ_q - μ_ν)/(σ_ν^2)

    # Partial derivative w.r.t σ_q
    partial_σ = - σ_q / σ_ν

    return (partial_μ, partial_σ)
end

function grad_expectation(q::Gamma{Float64}, ν::Gamma{Float64})
    "Compute the gradient of the expectation E_q(x) [log ν(x)] w.r.t. q(x)'s params"

    # Parameters
    α_q, θ_q = params(q)
    α_ν, θ_ν = params(ν)

    # Partial derivative w.r.t α_q
    partial_α = (α_ν-1)*polygamma(1, α_q) - θ_q / θ_ν

    # Partial derivative w.r.t σ_q
    partial_θ = (α_ν-1)/θ_q - α_q / θ_ν

    return (partial_α, partial_θ)
end

function grad_expectation(q::Union{Gamma{Float64},Normal{Float64}}, ν::Flat)
    "The gradient of a diverging integral, diverges."
    return (Inf, Inf)
end

function grad_expectation(ν::Flat, q::Union{Gamma{Float64},Normal{Float64}})
    "The gradient of a diverging integral, diverges."
    return (Inf, Inf)
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
const  WindThreshold=800
#meters
const windSpaceMin = 500
#meters
const windIdStartPos = 10
const aboveThresholdPercent= 50

#seconds
const WindingChangeTimeMaxInterval = 1
const windAlarmSwitchONMinTIME = 10
const windAlarmSwitchOFFMinTIME = 2*120

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
"""
    BaB(optimizer, ϵ::Float64)

BaB uses branch and bound to estimate the range of the output node.

# Problem requirement
1. Network: any depth, ReLU activation, single output
2. Input: hyperrectangle
3. Output: hyperrectangle (1d interval)

# Return
`CounterExampleResult` or `ReachabilityResult`

# Method
Branch and bound.
For branch, it uses iterative interval refinement.
For bound, it computes concrete bounds by sampling, approximated bound by optimization.
- `optimizer` default `GLPKSolverMIP()`
- `ϵ` is the desired accurancy for termination, default `0.1`.

# Property
Sound and complete.

# Reference
[R. Bunel, I. Turkaslan, P. H. Torr, P. Kohli, and M. P. Kumar,
"A Unified View of Piecewise Linear Neural Network Verification,"
*ArXiv Preprint ArXiv:1711.00455*, 2017.](https://arxiv.org/abs/1711.00455)
"""
@with_kw struct BaB
    optimizer::AbstractMathProgSolver = GLPKSolverMIP()
    ϵ::Float64                        = 0.1
end

function solve(solver::BaB, problem::Problem)
    (u_approx, u, x_u) = output_bound(solver, problem, :max)
    (l_approx, l, x_l) = output_bound(solver, problem, :min)
    bound = Hyperrectangle(low = [l], high = [u])
    reach = Hyperrectangle(low = [l_approx], high = [u_approx])
    return interpret_result(reach, bound, problem.output, x_l, x_u)
end

function interpret_result(reach, bound, output, x_l, x_u)
    if high(reach) < high(output) && low(reach) > low(output)
        return ReachabilityResult(:SAT, reach)
    end
    high(bound) > high(output)    && return CounterExampleResult(:UNSAT, x_u)
    low(bound) < low(output)      && return CounterExampleResult(:UNSAT, x_l)
    return ReachabilityResult(:Unknown, reach)
end

function output_bound(solver::BaB, problem::Problem, type::Symbol)
    nnet = problem.network
    global_concrete, x_star = concrete_bound(nnet, problem.input, type)
    global_approx = approx_bound(nnet, problem.input, solver.optimizer, type)
    doms = Tuple{Float64, Hyperrectangle}[(global_approx, problem.input)]
    index = ifelse(type == :max, 1, -1)
    while index * (global_approx - global_concrete) > solver.ϵ
        dom = pick_out(doms) # pick_out implements the search strategy
        subdoms = split_dom(dom[2]) # split implements the branching rule
        for i in 1:length(subdoms)
            dom_concrete, x = concrete_bound(nnet, subdoms[i], type) # Sample
            dom_approx = approx_bound(nnet, subdoms[i], solver.optimizer, type)
            if index * (dom_concrete - global_concrete) > 0
                (global_concrete, x_star) = (dom_concrete, x)
            end
            if index * (dom_approx - global_concrete) > 0
                add_domain!(doms, (dom_approx, subdoms[i]), type)
            end
        end
        global_approx = doms[1][1]
    end
    return (global_approx, global_concrete, x_star)
end

# Always pick the first dom
function pick_out(doms)
    return doms[1]
end

function add_domain!(doms::Vector{Tuple{Float64, Hyperrectangle}}, new::Tuple{Float64, Hyperrectangle}, type::Symbol)
    rank = length(doms) + 1
    index = ifelse(type == :max, 1, -1)
    for i in 1:length(doms)
        if index * (new[1] - doms[i][1]) >= 0
            rank = i
            break
        end
    end
    insert!(doms, rank, new)
end

# Always split the longest input dimention
function split_dom(dom::Hyperrectangle)
    max_value, index_to_split = findmax(dom.radius)
    return split_interval(dom, index_to_split)
end

# For simplicity
function concrete_bound(nnet::Network, subdom::Hyperrectangle, type::Symbol)
    points = [subdom.center, low(subdom), high(subdom)]
    values = Vector{Float64}(undef, 0)
    for p in points
        push!(values, sum(compute_output(nnet, p)))
    end
    value, index = ifelse(type == :min, findmin(values), findmax(values))
    return (value, points[index])
end

function approx_bound(nnet::Network, dom::Hyperrectangle, optimizer::AbstractMathProgSolver, type::Symbol)
    bounds = get_bounds(nnet, dom)
    model = Model(solver = optimizer)
    neurons = init_neurons(model, nnet)
    add_set_constraint!(model, dom, first(neurons))
    encode_Δ_lp!(model, nnet, bounds, neurons)
    index = ifelse(type == :max, 1, -1)
    o = sum(last(neurons))
    @objective(model, Max, index * o)
    status = solve(model, suppress_warnings = true)
    status == :Optimal && return getvalue(o)
    error("Could not find bound for dom: ", dom)
end


text:
using TextSegmentation
using Documenter

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

makedocs(;
    modules=[TextSegmentation],
    authors="Kento Kawasaki",
    repo="https://github.com/kawasaki-kento/TextSegmentation.jl/blob/{commit}{path}#{line}",
    sitename="TextSegmentation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kawasaki-kento.github.io/TextSegmentation.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    doctest = false,
)

deploydocs(;
    repo="github.com/kawasaki-kento/TextSegmentation.jl",
    devbranch="main",
)


text:
# Autogenerated wrapper script for Unicorn_jll for x86_64-apple-darwin
export libunicorn

JLLWrappers.@generate_wrapper_header("Unicorn")
JLLWrappers.@declare_library_product(libunicorn, "@rpath/libunicorn.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libunicorn,
        "lib/libunicorn.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestEvaluationMetrics

using Test
using LightGBM

@testset "merge_metrics" begin

    a = Dict(
        "dataset_a" => Dict(
            "some_metric" => [1., 2.],
            "another_metric" => [2.2, 3.3],
        ),
        "dataset_b" => Dict(
            "dodo_metric" => [0., 0.],
            "great_metric" => [0.1, 0.2],
        ),
    )

    b = Dict(
        "dataset_a" => Dict(
            "some_metric" => [3., 4.],
            "another_metric" => [4.4, 5.5],
        ),
        "dataset_b" => Dict(
            "dodo_metric" => [1., 1.],
            "great_metric" => [0.3, 0.4],
        ),
    )

    a_copy = deepcopy(a)
    b_copy = deepcopy(b)

    c = LightGBM.merge_metrics(a, b)

    # check no mutate
    @test a == a_copy
    @test b == b_copy

    # check no weird key munging from top/lower layer
    @test sort(collect(keys(c))) == ["dataset_a", "dataset_b"]
    @test sort(collect(keys(c["dataset_a"]))) == ["another_metric", "some_metric"]
    @test sort(collect(keys(c["dataset_b"]))) == ["dodo_metric", "great_metric"]

    # check the metrics concated properly
    @test c["dataset_a"]["some_metric"] == [1., 2., 3., 4.]
    @test c["dataset_a"]["another_metric"] == [2.2, 3.3, 4.4, 5.5]
    @test c["dataset_b"]["dodo_metric"] == [0., 0., 1., 1.]
    @test c["dataset_b"]["great_metric"] == [0.1, 0.2, 0.3, 0.4]

    d1 = Dict("breakme" => Dict{String, Vector{Float64}}())
    d2 = Dict("breakme2" => Dict{String, Vector{Float64}}())

    @test_throws ErrorException LightGBM.merge_metrics(d1, d2)

    e1 = Dict("ametric" => [1., 1., 1.])
    e2 = Dict("bmetric" => [1., 1., 1.])

    @test_throws ErrorException LightGBM.merge_metrics(e1, e2)

end

@testset "empty merge_metrics" begin

    # basically this test is just to make sure it works properly even when metrics is empty

    a = Dict{String, Vector{Float64}}()
    b = Dict{String, Vector{Float64}}()

    c = LightGBM.merge_metrics(a, b)

    @test typeof(c) == typeof(a)
    @test isempty(c)

    d1 = Dict{String, typeof(a)}()
    d2 = Dict{String, typeof(a)}()

    d3 = LightGBM.merge_metrics(d1, d2)

    @test typeof(d3) == typeof(d1)
    @test isempty(d3)

end

end # module


text:
# les bornes de for ne sont évaluées qu'une seule fois

function f(x)
    print(x, ".")
    x
end

for n = f(1) : f(10)
    print(n)
end
println()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
begin
    using LinearAlgebra, DataFrames,Optim, ForwardDiff, BenchmarkTools,Distributions,
        Expectations, QuantEcon, Statistics, GLM
    using Distributions: invsqrt2π, log2π, sqrt2, invsqrt2
    using DCDC
    using Test
    using Distributed, Suppressor
    using Plots
end

begin "Dynamic Decision Process"
    σ₀ = 1;
    ρ = 0.15;
    β = 1 / ( 1 + ρ);
    r = 0.2;
    α = 0.9;
    A = 1 + r;
    ddc = DynamicDecisionProcess(σ₀,β,α,A);
    # The log function has a solution
    value_function(s) = (log(1 - β * α) + (β * α) * log(β * α) /(1-β * α) ) / (1 - β) + α * log(s) /(1-β * α);
    policy_function(s) = (1 - β * α) * A * s^α;
    ddc.ValueFn.y = value_function.(ddc.ValueFn.xdata);
    ddc.PolicyFn.y = policy_function.(ddc.PolicyFn.xdata[:,1]);
    # computeEquilibrium(ddc);
    # ee = EulerEquation(σ₀,β,α);
end

# "Check Approximation of Policy Function"
begin
    scatter(ddc.PolicyFn.xdata[:,1],ddc.PolicyFn.y);
    scatter!(ddc.PolicyFn.xdata[:,1],ddc.PolicyFn(ddc.PolicyFn.xdata));
    scatter!(ddc.PolicyFn.xdata[:,1],policy_function.(ddc.PolicyFn.xdata[:,1]));
    # scatter!(ddc.PolicyFn.xdata[:,1],A * ddc.PolicyFn.xdata[:,1].^α)
end

# "Check Approximation of Value Function"
begin
    scatter(ddc.ValueFn.xdata[:,1] ,ddc.ValueFn.y);
    scatter!(ddc.ValueFn.xdata[:,1],ddc.ValueFn(ddc.ValueFn.xdata[:,1]));
end

begin "Check whether the optimal choice is correct"
    c_opt = deepcopy(ddc.PolicyFn.y);
    s = ddc.ValueFn.xdata;
    η = ddc.PolicyFn.xdata[:,2];
    ϵ = minimum(s)*0.9;
    s = ddc.PolicyFn.xdata[:,1];
    lb = ϵ;
    for n = 1:ddc.nSolve
        ub =  ddc.trans(0,s[n]) - ϵ;
        @show c_opt[n] = find_optim(s[n],η[n],lb,ub,ddc.β,ddc.trans,ddc.util,ddc.ValueFn);
    end
    s̃ = ddc.trans(c_opt,s);
    y=ddc.util(c_opt) + ddc.β * ddc.ValueFn(s̃);
    scatter(ddc.PolicyFn.xdata[:,1],c_opt);
    scatter!(ddc.PolicyFn.xdata[:,1],ddc.PolicyFn.y);
end

begin
    scatter(ddc.ValueFn.xdata[:,1],ddc.ValueFn.y);
    scatter!(ddc.ValueFn.xdata[:,1],value_function.(ddc.ValueFn.xdata[:,1]));
    scatter!(ddc.ValueFn.xdata[1:end-1,1],y);
end

begin
    c_opt = deepcopy(ddc.PolicyFn.y);
    s = ddc.PolicyFn.xdata[:,1];
    s̃ = ddc.trans(c_opt,s);
    y=ddc.util(c_opt) + ddc.β * ddc.ValueFn(s̃);
    scatter(ddc.ValueFn.xdata[:,1],ddc.ValueFn.y);
    scatter!(ddc.ValueFn.xdata[:,1],value_function.(ddc.ValueFn.xdata[:,1]));
    scatter!(ddc.ValueFn.xdata[1:end-1,1],y)
end

begin
    ee_scatter=[check_ee(ddc) for i = 1:100];
    @show mean(ee_scatter)
    scatter(ee_scatter);
end

#_____________________________________________________________________________
# linearize estimation method
#_____________________________________________________________________________

begin
    β_MC=zeros(10);
    for i = 1:10

        nM = 100;
        nT = 2;
        data = simulate_ddc(nM,nT,ddc);
        # data =simulate_ee(nM,nT,ee);
        # The moment condition is
        # u'(c_t) = β d_Trans(c_t,x_t) u'(c_t+1)
        a_t = []; a_t1 = [];
        s_t = []; s_t1 = [];
        for t = 1:(nT-1)
             a_t  = vcat(a_t,data.action[:,t]);
             s_t  = vcat(s_t, [s.s[1] for s in data.state[:,t]] );
             a_t1 = vcat(a_t1,data.action[:,t+1])
            s_t1  = vcat(s_t1, [s.s[1] for s in data.state[:,t+1]] );
        end
        a_t=convert(Array{Float64,1},a_t);
        s_t=convert(Array{Float64,1},s_t);
        a_t1=convert(Array{Float64,1},a_t1);
        s_t1=convert(Array{Float64,1},s_t1);
        dtrans=ddc.dtrans;
        dtrans_s = (a,s) -> dtrans(a,s)[2].data;
        R = dtrans_s.(a_t1,s_t1);
        y=log.(a_t1) - log.(a_t)
        X=hcat(ones(length(y)),log.(R));
        b=inv(X'*X)*X'*y;
        @show β_MC[i]=exp(b[1]/b[2]);
    end

    @show b=inv(X'*X)*X'*y;
    @show 1/b[2];
    @show "True beta"
    @show β
    @show exp(b[1]/b[2]);
    "End";
end

begin "The estimation"
    policy_approx = ApproxFn(s_t1,a_t1,:epan,2);
    scatter(s_t1,a_t1);
    scatter!(s_t1,policy_approx(s_t1));

end
# plot!(ee.PolicyFn.xdata[:,1],ee.PolicyFn.y);
# scatter(ddc.PolicyFn.xdata[:,2]);
# checked = [check_ee(ddc) for i = 1:100];
# scatter(checked[checked .< 100])
#

#_____________________________________________________________________________
# Nonlinearize estimation method
#_____________________________________________________________________________

using Flux #Get derivatives
using Flux.Tracker
using Flux.Tracker: update!
# Solve the Euler equations

lb = [0,0];
ub = [1,100];

σ̂=2;β̂=0.8;

begin
    nM = 100;
    nT = 3;
    data = simulate_ddc(nM,nT,ddc);
    # data =simulate_ee(nM,nT,ee);
    # The moment condition is
    # u'(c_t) = β d_Trans(c_t,x_t) u'(c_t+1)
    a_t = [];
    s_t = [];
    for t = 1:(nT)
         global a_t  = vcat(a_t,data.action[:,t]);
         global s_t  = vcat(s_t, [s.s[1] for s in data.state[:,t]] );
    end
    a_t=convert(Array{Float64,1},a_t);
    s_t=convert(Array{Float64,1},s_t);
end

trans=ddc.trans;
function _moment(θ,trans,a_t,s_t,policy_approx)
    σ̂,β̂=θ; s_t1_simul=trans(a_t,s_t);
    a_t1_simul=policy_approx(s_t1_simul);
    util=Utility(σ̂); dutil=(c,η) -> Tracker.gradient(util,c,η)[1].data;
    # η=randn(100)./100;
    dc_t=dutil.(a_t,zeros(length(a_t)));
    dc_t1=dutil.(a_t1_simul,zeros(length(a_t)));
    dtrans_s = (a,s) -> dtrans(a,s)[2].data;
    R = dtrans_s.(a_t1_simul,s_t1_simul);
    ϵ= (β̂ * dc_t1 .* R ) ./ dc_t .-1;
    ϵ.* s_t;
end

function f_obj(θ)
    CF =_moment(θ,ddc.trans,a_t,s_t,policy_approx);
    sum(CF)^2;
end

result = optimize(f_obj,[0.1,0.01])
result.minimizer
f_obj(result.minimizer)
f_obj([β,1]+randn(2))

# Results of Optimization Algorithm
#  * Algorithm: Nelder-Mead
#  * Starting Point: [0.1,0.01]
#  * Minimizer: [0.47126793573588627,1.0128364461302402]
#  * Minimum: 2.226862e-11
#  * Iterations: 49
#  * Convergence: true
#    *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: true
#    * Reached Maximum Number of Iterations: false
#  * Objective Calls: 98


text:
"""
    Scene

A Scene is a specific instance of the Frame type defined in Records. It represents a collection of vehicles at a given time.

# Constructors
    - `Scene(n::Int=100)`
    - `Scene(arr::Vector{Vehicle})`
"""
const Scene = Frame{Vehicle}
Scene(n::Int=100) = Frame(Vehicle, n)
Scene(arr::Vector{Vehicle}) = Frame{Vehicle}(arr, length(arr))

Base.show(io::IO, scene::Scene) = print(io, "Scene(with $(length(scene)) cars)")

"""
    Base.convert(::Type{Vehicle}, veh::Entity{VehicleState, D, Int64}) where D<:AbstractAgentDefinition

Converts an entity in Vehicle (it is converting the agent definition only)
"""
function Base.convert(::Type{Entity{VehicleState, VehicleDef, I}}, veh::Entity{VehicleState, D, I}) where {D<:AbstractAgentDefinition, I}
    vehdef = VehicleDef(class(veh.def), length(veh.def), width(veh.def))
    return Entity{VehicleState, VehicleDef, I}(veh.state, vehdef, veh.id)
end

"""
    SceneRecord
A SceneRecord is a specific instance of the QueueRecord type defined in Records.jl. It represents a collection of Scenes.

# constructor
    SceneRecord(capacity::Int, timestep::Float64, frame_capacity::Int=100)
"""
const SceneRecord = QueueRecord{Vehicle}
SceneRecord(capacity::Int, timestep::Float64, frame_capacity::Int=100) = QueueRecord(Vehicle, capacity, timestep, frame_capacity)
Base.show(io::IO, rec::SceneRecord) = print(io, "SceneRecord(nscenes=", nframes(rec), ")")


text:
using BinaryBuilder, Pkg.Artifacts

include("../common.jl")

name = "MegaRust"
version = v"1.18.3"

sources = [
    # TODO: Switch to musl once https://github.com/rust-lang/rustup.rs/pull/1882 is released
    "https://static.rust-lang.org/rustup/archive/$(version)/x86_64-unknown-linux-gnu/rustup-init" =>
    "a46fe67199b7bcbbde2dcbc23ae08db6f29883e260e23899a88b9073effc9076",
]

# The first thing we're going to do is to install Rust for all targets into a single prefix
script = raw"""
cd ${WORKSPACE}/srcdir

# We reinstall RustBase because that's easier than re-using it
export CARGO_HOME=${prefix}
export RUSTUP_HOME=${prefix}
mv *-rustup-init rustup-init
chmod +x rustup-init
./rustup-init -y --no-modify-path --default-host=${rust_host}

# Collection of all rust targets we will download toolchains for:
RUST_TARGETS=(
    aarch64-unknown-linux-gnu
    aarch64-unknown-linux-musl
    armv7-unknown-linux-gnueabihf
    armv7-unknown-linux-musleabihf
    i686-unknown-linux-gnu
    i686-unknown-linux-musl
    x86_64-unknown-linux-gnu
    x86_64-unknown-linux-musl
    powerpc64le-unknown-linux-gnu
    i686-pc-windows-gnu
    x86_64-pc-windows-gnu
    x86_64-apple-darwin
    x86_64-unknown-freebsd
)

for rust_target in "${RUST_TARGETS[@]}"; do
    # Install our target-specific stuffs
    ${CARGO_HOME}/bin/rustup target add ${rust_target}
done

# We're going to bundle cargo-edit since it's a useful dep
export OPENSSL_STATIC=yes
${CARGO_HOME}/bin/cargo install cargo-edit
"""

# We assemble this giant tarball, then will split it up immediately after this:
platforms = [Linux(:x86_64; libc=:glibc)]
products = [
    ExecutableProduct("cargo", :cargo),
]
dependencies = [
    "OpenSSL_jll",
]
ndARGS = filter(a -> !occursin("--deploy", a), ARGS)
build_info = build_tarballs(ndARGS, "MegaRust", version, sources, script, platforms, products, dependencies; skip_audit=true)

# We don't actually need the .tar.gz it creates, so delete that to save space
rm(joinpath("products", first(values(build_info))[1]))

# Take the hash of the unpacked MegaRust artifact, then split it into a bunch of smaller ones
mega_rust_path = artifact_path(first(values(build_info))[3])
rust_host = Linux(:x86_64; libc=:glibc)
rust_host_triplet = BinaryBuilder.map_rust_target(rust_host)

for target_platform in supported_platforms()
    rust_target_triplet = BinaryBuilder.map_rust_target(target_platform)
    @info("Generating artifacts for $(rust_target_triplet)...")
    unpacked_hash = create_artifact() do dir
        srcpath = joinpath(mega_rust_path, "toolchains", "stable-$(rust_host_triplet)", "lib", "rustlib", rust_target_triplet)
        dstpath = joinpath(dir, "toolchains", "stable-$(rust_host_triplet)", "lib", "rustlib")
        mkpath(dstpath)
        cp(srcpath, joinpath(dstpath, rust_target_triplet))

        # Our mingw rust shards need to have their crt2.o updated
        # https://github.com/rust-lang/rust/issues/48272#issuecomment-429596397
        if Sys.iswindows(target_platform)
            # Find the corresponding mingw toolchain within a GCC shard
            all_cs = BinaryBuilder.all_compiler_shards()
            cs = first(filter(cs -> cs.name == "GCCBootstrap" && cs.target == target_platform && cs.archive_type == :unpacked, all_cs))
            cs_path = BinaryBuilder.mount(cs, "")

            # Locate our GCC's mingw's crt2.o
            crt_src = joinpath(cs_path, triplet(target_platform), "sys-root", "lib", "crt2.o")
            # Overwrite the one that rust ships with
            crt_dst = joinpath(dstpath, rust_target_triplet, "lib", "crt2.o")
            cp(crt_src, crt_dst; force=true)
        end
    end
    squashfs_hash = unpacked_to_squashfs(unpacked_hash, "RustToolchain", version; platform=rust_host, target=target_platform)

    # Upload them both to GH releases on Yggdrasil
    upload_and_insert_shards("JuliaPackaging/Yggdrasil", "RustToolchain", version, unpacked_hash, squashfs_hash, rust_host; target=target_platform)
end

# Finally, we do RustBase:
unpacked_hash = create_artifact() do dir
    cp(joinpath(mega_rust_path, "bin"), joinpath(dir, "bin"))
    cp(joinpath(mega_rust_path, "toolchains"), joinpath(dir, "toolchains"))
    rm(joinpath(dir, "toolchains", "stable-$(rust_host_triplet)", "share"); recursive=true)
    rm(joinpath(dir, "toolchains", "stable-$(rust_host_triplet)", "etc"); recursive=true)
    for rust_target_triplet in BinaryBuilder.map_rust_target.(supported_platforms())
        rm(joinpath(dir, "toolchains", "stable-$(rust_host_triplet)", "lib", "rustlib", rust_target_triplet); recursive=true)
    end

    # Also generate "config" file for Cargo where we give it the linkers for all our targets
    open(joinpath(dir, "config"), "w") do io
        write(io, """
        # Configuration file for `cargo`
        """)
        for platform in supported_platforms()
            write(io, """
            [target.$(BinaryBuilder.map_rust_target(platform))]
            linker = "$(triplet(platform))-gcc"
            """)
        end
    end
end

squashfs_hash = unpacked_to_squashfs(unpacked_hash, "RustBase", version; platform=rust_host)
upload_and_insert_shards("JuliaPackaging/Yggdrasil", "RustBase", version, unpacked_hash, squashfs_hash, rust_host)

max_stars_repo_name: devmotion/Yggdrasil

max_issues_repo_name: devmotion/Yggdrasil

max_forks_repo_name: devmotion/Yggdrasil


text:
export outDims

function outDims(cLayer::PL, Ai::AbstractArray{T,N}) where {PL <: PaddableLayer, N, T}
    return outDims(cLayer, size(Ai))
end

function outDims(cLayer::PL, AiS::Tuple) where {PL <: PaddableLayer}

    f = cLayer.f
    s = cLayer.s
    inputS = cLayer.inputS[1:end-2]
    paddedS = paddedSize(cLayer, AiS)[1:end-2]
    ci, m = cLayer.inputS[end-1:end]
    co = cLayer.channels
    outputS = []
    for i=1:length(f)
        n = (paddedS[i] - f[i]) ÷ s[i] + 1
        push!(outputS, n)
    end

    return (outputS..., co, m)

end