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 = "GDB"
version = v"10.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/gdb/gdb-10.1.tar.xz", "f82f1eceeec14a3afa2de8d9b0d3c91d5a3820e23e0a01bbb70ef9f0276b62c0")
]

# Bash recipe for building across all platforms
script = raw"""
apk add texinfo
cd $WORKSPACE/srcdir/gdb-10.1/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-expat
make -j${nproc} all
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "musl"),
    Platform("i686", "windows"; ),
    Platform("x86_64", "windows"; )
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("gdbserver", :gdbserver),
    ExecutableProduct("gdb", :gdb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d")),
    Dependency("Expat_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
Return functional and gradient of the LS objective
* `last_x::Vector{Float64}` : buffer is only updated when x!=last_x, and mediumified such that last_x=x
"""
function func(x::Vector{Float64}, last_x::Vector{Float64}, pa::PFWI)
	global fwi_to

	if(!isequal(x, last_x))
		copyto!(last_x, x)
		# do forward modelling, apply F x
		@timeit_debug fwi_to "F!" F!(pa, x)
	end

	# compute misfit
	f = func_grad!(pa.paTD)
	return f
end

function grad!(storage, x::Vector{Float64}, last_x::Vector{Float64}, pa::PFWI)
	global fwi_to
	reset_timer!(fwi_to)

	# (inactive when applied on same model)
	if(!isequal(x, last_x))
		copyto!(last_x, x)
		# do forward modelling, apply F x
		@timeit_debug fwi_to "F!" F!(pa, x)
	end

	# compute functional and get ∇_d J (adjoint sources)
	f = func_grad!(pa.paTD, :dJx);

	# update adjoint sources after time reversal
	update_adjsrc!(pa.adjsrc, pa.paTD.dJx, pa.adjageom)

	# do adjoint modelling here with adjoint sources Fᵀ F P x
	@timeit_debug fwi_to "Fadj!" Fadj!(pa)

	# adjoint of interpolation
	spray_gradient!(storage,  pa)

	return storage
end

function ζfunc(x, last_x, pa::PFWI{T1,T2,T3}) where {T1,T2,T3<:Union{LS,Migr,Migr_FD}}
	return func(x, last_x, pa)
end

function ζgrad!(storage, x, last_x, pa::PFWI{T1,T2,T3})  where {T1,T2,T3<:Union{LS,Migr,Migr_FD}}
	return grad!(storage, x, last_x, pa)
end

function ζfunc(x, last_x, pa::PFWI{T1,T2,LS_prior})  where {T1,T2}
	f1=func(x, last_x, pa)

	# calculate the generalized least-squares error
	# note: change the inverse model covariance matrix `pmgls.Q` accordingly
	# f2=func_grad!(nothing, x, pa.mx.prior, obj.pmgls)

	return f1*obj.pdgls+f2
end

function ζgrad!(storage, x, last_x, pa::PFWI{T1,T2,LS_prior}) where {T1,T2}
	g1=pa.mx.gm[1]
	grad!(g1, x, last_x, pa)

	g2=pa.mx.gm[2]
	# func_grad!(g2, x, pa.mx.prior, obj.pmgls)

	rmul!(g1, obj.pdgls)

	for i in eachindex(storage)
		@inbounds storage[i]=g1[i]+g2[i]
	end
	return storage
end

"""
Perform a forward simulation.
Update `pa.paTD.x`.
This simulation is common for both functional and gradient calculation.
During the computation of the gradient, we need an adjoint simulation.n
Update the buffer, which consists of the modelled data
and boundary values for adjoint calculation.

# Arguments

* `x::Vector{Float64}` : inversion variable
* `pa::PFWI` : parameters that are constant during the inversion
* if x is absent, using `pa.mediumm` for modeling
"""
function F!(pa::PFWI{FdtdAcou,T1,T2}, x) where {T1,T2}

	# switch off born scattering
	#pa.paf.c.born_flag=false

	# initialize boundary, as we will record them now
	initialize_boundary!(pa.paf)

	if(!(x===nothing))
		# project x, which lives in mediumi, on to model space (mediumm)
		x_to_mediumm!(pa, x)
	end

	# update model in the forward engine
	update!(pa.paf.c, pa.mediumm)

	pa.paf.c.activepw=[1,]
	pa.paf.c.illum_flag=false
	pa.paf.c.sflags=[2, 0]
	pa.paf.c.rflags=[1, 0] # record only after first scattering
	update!(pa.paf,[pa.srcwav,pa.adjsrc])
	pa.paf.c.backprop_flag=1
	pa.paf.c.gmodel_flag=false

	update!(pa.paf);

	# copy data to evaluate misfit
	dcal=pa.paf.c.data[1]
	copyto!(pa.paTD.x,dcal)
end

"""
Born modeling with `mediumm` as the perturbed model and `mediumm0` as the background model.
"""
function F!(pa::PFWI{FdtdAcouBorn,T1,T2}, x) where {T1,T2}

	# update background model in the forward engine
	update!(pa.paf.c, pa.mediumm0)
	if(!(x===nothing))
		# project x, which lives in mediumi, on to model space (mediumm)
		x_to_mediumm!(pa, x)
	end
	# update perturbed models in the forward engine
	update_δmods!(pa.paf.c, pa.mediumm)

	Fbornmod!(pa::PFWI)
end

function Fbornmod!(pa::PFWI{FdtdAcouBorn,T1,T2}) where {T1,T2}

	# switch on born scattering
	#pa.paf.c.born_flag=true

	pa.paf.c.activepw=[1,2] # two wavefields are active
	pa.paf.c.illum_flag=false
	pa.paf.c.sflags=[2, 0] # no sources on second wavefield
	pa.paf.c.rflags=[0, 1] # record only after first scattering

	# source wavelets (for second wavefield, they are dummy)
	update!(pa.paf,[pa.srcwav,pa.adjsrc])

	# actually, should record only when background field is changed
	pa.paf.c.backprop_flag=1 # store boundary values for gradient later

	pa.paf.c.gmodel_flag=false # no gradient

	update!(pa.paf);
	dcal=pa.paf.c.data[2]
	copyto!(pa.paTD.x,dcal)

	# switch off born scattering once done
	#pa.paf.c.born_flag=false
end

"""
Perform adjoint modelling in `paf` using adjoint sources `adjsrc`.
"""
function Fadj!(pa::PFWI)

	# need to explicitly turn off the born flag for adjoint modelling
	#pa.paf.c.born_flag=false

	# require gradient, switch on the flag
	pa.paf.c.gmodel_flag=true

	# both wavefields are active
	pa.paf.c.activepw=[1,2]

	# no need of illum during adjoint modeling
	pa.paf.c.illum_flag=false

	# force boundaries in first pw and back propagation for second pw
	pa.paf.c.sflags=[-2,2]
	pa.paf.c.backprop_flag=-1

	# update source wavelets in paf using adjoint sources
	update!(pa.paf,[pa.srcwav,pa.adjsrc])

	# no need to record data during adjoint propagation
	pa.paf.c.rflags=[0,0]

	# adjoint modelling
	update!(pa.paf);

	# put rflags back
	pa.paf.c.rflags=[1,1]

	return pa.paf.c.gradient
end

"""
```
F=LinearMap(pa)
```
If `pa` is an instance of `SeisInvExpt`, then
return the linearized forward modeling operator `F`, such that
`F*x` can be computed without explicitly storing the operator matrix (see `LinearMaps.jl`).
The imaging/migration operator is given by `transpose(F)`.
These operators are the building blocks of iterative optimization schemes.
"""
function LinearMaps.LinearMap(pa::PFWI{FdtdAcouBorn,T2,T3}) where {T2,T3}
	fw=(y,x)->Fborn_map!(y, x, pa)
	bk=(y,x)->Fadj_map!(y, x, pa)

	return LinearMap(fw, bk,
		sum(length.(pa.paTD.dJx)),  # length of output
		  xfwi_ninv(pa), # length of input
		  ismutating=true)
end

function Fborn_map!(δy, δx, pa)
	δx_to_δmods!(pa, δx)
	Fbornmod!(pa)
	copyto!(δy, pa.paTD.x)
end

function Fadj_map!(δy, δx, pa)
	copyto!(pa.paTD.dJx, δx)

	# adjoint sources
	update_adjsrc!(pa.adjsrc, pa.paTD.dJx, pa.adjageom)

	# adjoint simulation
	Fadj!(pa)

	# chain rule corresponding to reparameterization
	pert_chainrule!(pa.mxm.gx, pa.paf.c.gradient, pa.mediumm0, pa.parameterization)

	# finally, adjoint of interpolation
	Interpolation.interp_spray!(δy,
			     pa.mxm.gx, pa.paminterp, :spray,
			     count(pa.parameterization.≠:null))
end

ext: jl


text:
Base.eval(:(have_color = true))

using Documenter
using Literate
using Glob
using Random
using HMM2d

Random.seed!(2019)
ENV["GKSwstype"] = "nul"
#=
find_examples() = map(x -> joinpath("examlpes/", split(basename(x), ".")[1]), glob("*.jl", "examples/"))

if !("SKIP_EXAMPLES" in keys(ENV))
    for example in find_examples()
        Literate.markdown("$(example).jl", "docs/src/examples", documenter = true)
    end
end
=#

makedocs(sitename = "HMM2d", modules = [HMM2d], pages = ["index.md", "Manual" => ["placeholder.md"], "theory.md", "_index.md" ])

deploydocs( repo = "github.com/benjaminwsebastian/HMM2d.jl.git" )


text:
using Dates
using PowerSystems
using Setfield
using Test
using TimeSeries

const EVIPRO_DATA = abspath(joinpath(dirname(Base.find_package("PowerSystems")), "..", "data", "evi-pro", "FlexibleDemand_1000.mat"))

macro trytotest(f)
    if isfile(EVIPRO_DATA)
        f
    else
        @warn string("Demand-response tests not run because file '", EVIPRO_DATA, "' is missing.")
    end
end

@testset "Reading EVIpro dataset" begin
    @trytotest begin
        @test begin
           @info string("Testing using EVIpro dataset \"", EVIPRO_DATA, "\".")
           bevs = populate_BEV_demand(EVIPRO_DATA)
           length(bevs) == 1000
        end
    end
end

function augment(bev)
    bev = @set bev.power     = update(bev.power    , Time(23,59,59,999), NaN                       )
    bev = @set bev.locations = update(bev.locations, Time(23,59,59,999), ("", (ac = NaN, dc = NaN)))
    bev
end

function checkcharging(f)
    bevs = populate_BEV_demand(EVIPRO_DATA)
    i = 0
    for bev in bevs
        i += 1
        bev = augment(bev)
        @test begin
            charging = f(bev)
            (balance, rates, battery) = verify(bev, charging, message=string("BEV ", i, " in '", EVIPRO_DATA, "'"))
            balance && rates && battery
        end
    end
end

@testset "Earliest demands on EVIpro dataset" begin
    @trytotest begin
        checkcharging(earliestdemands)
    end
end

@testset "Latest demands on EVIpro dataset" begin
    @trytotest begin
        checkcharging(latestdemands)
    end
end

@testset "Greedy demands on EVIpro dataset" begin
    # FIXME: `PowerSystems.greedycharging` fails some tests.
    @error "`PowerSystems.greedycharging` fails some tests."
#   @trytotest begin
#       checkcharging(x -> greedydemands(x, map(v -> 1., x.power)) |> locateddemand)
#   end
end


text:
import TESOptimalFiltering: filter5lag, smoother, calculate_filter, autocorrelation
using Base.Test, HDF5
srand(1000)
# write your own tests here
pulse_model = 100*(exp(-(1:1000)/150)-exp(-(1:1000)/350));
noise_data = randn(1000000);
dt = 9.6e-6;
f_3db = 1e30;
noise_autocorr = autocorrelation(noise_data, length(pulse_model));
filter, variance = calculate_filter(pulse_model, noise_autocorr, f_3db, dt)
pulses = Vector{Float64}[pulse_model*i+randn(length(pulse_model)) for i=1:10] # create pulses with heights 1 thru 1-10
pulses2d = zeros(Float64, length(pulse_model), 10 )
for i in eachindex(pulses) pulses2d[:,i] = pulses[i] end
filt_value, filt_phase = filter5lag(filter, pulses[1]) # filter the first pulse
filt_values, filt_phases = filter5lag(filter, pulses) # filter all the pulses
filt_values2d, filt_phases2d = filter5lag(filter, pulses2d) # filter all the pulses

@test filt_value == filt_values[1]
@test filt_phase == filt_phases[1]
@test all(filt_values2d .== filt_values)
@test all(filt_phases2d .== filt_phases)

# results seem noiser than expected
#@test all(mean(abs(filt_values - collect(1:10))) .< sqrt(variance)*5) # check that pulse heights came out about right

pmm = pulse_model-mean(pulse_model);
approx_filter = pmm/dot(pmm,pmm);
@test maximum(approx_filter-filter)<1e-4

# compare to data from 20150730_50mM_irontris_100ps_xes_chan479.ljh caclulated by python mass on sept 14 2015
h5 = h5open(joinpath(Pkg.dir("Mass2/test"),"pythonmass_refdata.hdf5"))

filter, variance = calculate_filter(h5["average_pulse"][:], h5["noise_autocorr"][:], read(h5["f_3db"]), read(h5["dt"]))
#normalize so filter returns pulse height of average pulse
filter*=maximum(h5["average_pulse"][:])
vdv = 1/sqrt(variance)/sqrt(8*log(2))

@test abs(vdv/read(h5["vdv"])-1)<1e-3
# mass filter is shortened by 2 points on either end
# check that the filters are roughly similar by comparing their filtered values
@test abs(dot(filter[3:end-2],filter[3:end-2])/dot(filter[3:end-2],h5["filter"][:])-1) < 2/1000


text:
export StripTod, split_into_n, allocate_tod
export pwrq1, pwrq2, pwru1, pwru2
export demq1, demq2, demu1, demu2
export PWR_Q1_RANK, PWR_Q2_RANK, PWR_U1_RANK, PWR_U2_RANK
export DEM_Q1_RANK, DEM_Q2_RANK, DEM_U1_RANK, DEM_U2_RANK

import Printf: @printf
using AstroTime

"""
    split_into_n(length, num_of_segments)

Given the `length` of the array, the function returns a split of the
elements in `num_of_segments` subsets whose length is roughly the
same.

The result is an array containing the number of elements which must go
in each section; all the elements in the array have the property of
being similar and their sum being equal to `length`.

# Examples
```julia-repl
julia> split = split_into_n(20, 3)
3-element Array{Int64,1}:
6
7
7

julia> sum(split)
20
```
"""
function split_into_n(length, num_of_segments)
    @assert num_of_segments > 0
    @assert length >= num_of_segments
    start_pos = zeros(Int, num_of_segments+1)

    for i in 1:num_of_segments+1
        start_pos[i] =  floor(((i-1)*length/num_of_segments))
    end

    return start_pos[2:end]-start_pos[1:end-1]
end

@doc raw"""
A matrix of time-ordered data for a subset of polarimeters

The fields of this structure are the following:

- `time_range`: a range representing the time of each sample in the
  TOD.

- `samples`: a N×8×m matrix, where N is the number of samples in the
  timeline, and m is the number of polarimeters. The middle dimension
  spans the three Stokes parameters `I`, `Q`, and `U` (in this order).
  The value `N` is equal to `length(time_range)`

"""
struct StripTod{T <: Real, S}
    polarimeters::Any
    time_range::AbstractRange{S}
    samples::Array{T,3}
end

function Base.show(io::IO, tod::StripTod{T, S}) where {T, S}
    println(io, string(
        "TOD($(length(tod.polarimeters)) polarimeters, ",
        "$(length(tod.time_range)) samples/polarimeter, ",
        "time range from $(minimum(tod.time_range)) to ",
        "$(maximum(tod.time_range)))",
    ))
end

const PWR_Q1_RANK = 1
const PWR_Q2_RANK = 2
const PWR_U1_RANK = 3
const PWR_U2_RANK = 4
const DEM_Q1_RANK = 5
const DEM_Q2_RANK = 6
const DEM_U1_RANK = 7
const DEM_U2_RANK = 8

for (symb, ch_name, value) in [
    (:PWR_Q1_RANK, "PWR_Q1", PWR_Q1_RANK),
    (:PWR_Q2_RANK, "PWR_Q2", PWR_Q2_RANK),
    (:PWR_U1_RANK, "PWR_U1", PWR_U1_RANK),
    (:PWR_U2_RANK, "PWR_U2", PWR_U2_RANK),
    (:DEM_Q1_RANK, "DEM_Q1", DEM_Q1_RANK),
    (:DEM_Q2_RANK, "DEM_Q2", DEM_Q2_RANK),
    (:DEM_U1_RANK, "DEM_U1", DEM_U1_RANK),
    (:DEM_U2_RANK, "DEM_U2", DEM_U2_RANK),
]
    @eval begin
        @doc """
    $($ch_name)_RANK = $($value)

The rank of `$($ch_name)` in the `samples` field of a [`StripTod`](@ref) variable.

"""
        $symb
    end
end

for (index, ch_name, fn_name) in [
    (PWR_Q1_RANK, "PWR_Q1", :pwrq1),
    (PWR_Q2_RANK, "PWR_Q2", :pwrq2),
    (PWR_U1_RANK, "PWR_U1", :pwru1),
    (PWR_U2_RANK, "PWR_U2", :pwru2),
    (DEM_Q1_RANK, "DEM_Q1", :demq1),
    (DEM_Q2_RANK, "DEM_Q2", :demq2),
    (DEM_U1_RANK, "DEM_U1", :demu1),
    (DEM_U2_RANK, "DEM_U2", :demu2),
]
    @eval begin
        function ($fn_name)(t::StripTod, pol::T) where T <: Integer
            @view t.samples[:, $index, pol]
        end

        @doc """
    $($fn_name)(tod::Stripeline.StripTod, pol::T) where T <: Integer

Return a view of the timeline for channel `$($ch_name)` in the
time-ordered data for the polarimeter indexed by `pol` in `tod`. It is
a shortcut for

    tod.samples[:, $($ch_name)_RANK, pol]

"""
        $fn_name
    end
end

@doc raw"""
    allocate_tod(t, time_range, polarimeters; mpi_rank=0, mpi_size=1, zero_tod=true)

Allocate memory for a time-ordered data matrix appropriate for the set
of polarimeters in the variable `polarimeters` (an array or range,
only its length is used and it is copied as-is in the returned
[`StripTod`](@ref) value). The type of each sample is set to `t`. The
variable `time_range` must be a range of floating-point values or
epochs (as defined by the
[AstroTime](https://juliaastro.github.io/AstroTime.jl/stable/)
package) which represent the full set of samples acquired during an
observation.

The optional parameters `mpi_rank` and `mpi_size` can be used in MPI
scripts to allocate only the subset of samples that must be processed
by the current MPI process.

If `zero_tod` is true, each sample in the TOD is set to zero;
otherwise, its value is undefined. Using `zero_tod = false` is
slightly faster, but it should be used only if you are going to set by
hand every sample in the TOD once `allocate_tod` returns.

Return a [`StripTod`](@ref) value.

### Example ###

```jldoctest; setup = :(using Stripeline)
julia> tod = allocate_tod(Float32, 0.0:0.01:100.0, 1:5)
TOD(5 polarimeters, 10001 samples/polarimeter, time range from 0.0 to 100.0)
```

"""
function allocate_tod(
    t::Type{T},
    time_range::R,
    polarimeters;
    mpi_rank = 0,
    mpi_size = 1,
    zero_tod = true,
) where {S, V, R <: AbstractRange, T <: Real}

    rank_idx = mpi_rank + 1 # Same as rank, but it starts from 1

    nsamples_arr = split_into_n(length(time_range), mpi_size)
    nsamples = nsamples_arr[rank_idx]

    first_idx = sum(nsamples_arr[1:(rank_idx - 1)]) + 1
    last_idx = first_idx + nsamples - 1

    # Determine the time range for the current rank
    cur_range = time_range[first_idx]:step(time_range):time_range[last_idx]

    # Size of the field `samples`
    matrix_size = (length(cur_range), 8, length(polarimeters))

    StripTod(
        polarimeters,
        cur_range,
        zero_tod ? zeros(T, matrix_size) : Array{T}(undef, matrix_size),
    )
end


text:
using Colors, Compose
using Cairo, Fontconfig
using Plots
import FileIO
import Cairo, Fontconfig

using Compose: circle, rectangle

include("Assisting Files/Disk.jl")

set_default_graphic_size(10cm, 10cm)

function main_func()

    num_steps = 5

    expression = "a(b(), c())"

    Frame1 = Dict{String, Disk}()
    Frame1["a"] = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex(0.8 * cosd(20), 0.8 * sind(20)), complex(0.8 * cosd(-20), 0.8 * sind(-20))])
    Frame1["b"] = Disk("orange", [0.2], [])
    Frame1["c"] = Disk("red", [0.2], [])

    Frame2 = deepcopy(Frame1)
    Frame2["a"] = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex(0.8 * cosd(200), 0.8 * sind(200)), complex(0.8 * cosd(160), 0.8 * sind(160))])

    Transition = deepcopy(Frame1)

    mkpath("Saved Images/Disk Operad")

    for j in 1:60
        frame = compose(context(),
                        compose(context(), Compose.text(0.1, 0.1, string(floor(j/60 * 100), "%"), hcenter, vcenter), fontsize(5)),
                        compose(context(), circle(0.5, 0.5, Transition["a"].radius[1]*0.5),  fill("yellow")),
                        compose(context(), rectangle(), fill(colorant"white")))
        for i in 0:(60/(num_steps-1)):60
            #file_path = string("Saved Images/Disk Operad/", string(i), ".png")

            new1_angle1 = (((60 - i)/60) * angle(Frame1["a"].parameters[1])) + ((i/60) * angle(Frame2["a"].parameters[1]))
            new1_radius1 = (((60 - i)/60) * abs(Frame1["a"].parameters[1])) + ((i/60) * abs(Frame2["a"].parameters[1]))
            new1_angle2 = ((((60 - i)/60) * angle(Frame1["a"].parameters[2])) + ((i/60) * angle(Frame2["a"].parameters[2])))
            new1_angle2 = new1_angle2 - (2*pi)
            new1_radius2 = (((60 - i)/60) * abs(Frame1["a"].parameters[2])) + ((i/60) * abs(Frame2["a"].parameters[2]))
            polar1 = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex( new1_radius1 * cos(new1_angle1) + 0.2, new1_radius1 * sin(new1_angle1)), complex( new1_radius2 * cos(new1_angle2 - 2 *(new1_angle2 + deg2rad(20))) - 0.2, new1_radius2 * sin(new1_angle2 - 2 *(new1_angle2 + deg2rad(20))))])

            polar2 = Disk(colorant"rgba(0,0,0,0)", [0.7], [(((60 - i)/60) * Frame1["a"].parameters[1]) + ((i/60) * Frame2["a"].parameters[1]), (((60 - i)/60) * Frame1["a"].parameters[2]) + ((i/60) * Frame2["a"].parameters[2])])
            polar2 = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex( abs(polar2.parameters[1]) * cos(angle(polar2.parameters[1])) + 0.2, abs(polar2.parameters[1]) * sin(angle(polar2.parameters[1]))), complex( abs(polar2.parameters[2]) * cos(angle(polar2.parameters[2])) - 0.2, abs(polar2.parameters[2]) * sin(angle(polar2.parameters[2])))])

            Transition["a"] = Disk(colorant"rgba(0,0,0,0)", [0.7], [( ((60-j)/60) * polar1.parameters[1] ) + ( ((j)/60) * polar2.parameters[1] ), ( ((60-j)/60) * polar1.parameters[2] ) + ( ((j)/60) * polar2.parameters[2] )])

            frame = compose(context(), disk_compose_single_base(Meta.parse(expression), Transition, 0), frame)
            #draw(PNG(file_path, 10cm, 10cm, dpi=250), frame)
        end
        file_path = string("Saved Images/Disk Operad/", string(j), ".png")
        draw(PNG(file_path, 10cm, 10cm, dpi=250), frame)
    end

    for j in 1:60
        frame = compose(context(),
                        compose(context(), Compose.text(0.1, 0.1, string(floor(j/60 * 100), "%"), hcenter, vcenter), fontsize(5)),
                        compose(context(), circle(0.5, 0.5, Transition["a"].radius[1]*0.5),  fill("yellow")),
                        compose(context(), rectangle(), fill(colorant"white")))
        for i in 0:(60/(num_steps-1)):60
            #file_path = string("Saved Images/Disk Operad/", string(i), ".png")

            new1_angle1 = (((60 - i)/60) * angle(Frame1["a"].parameters[1])) + ((i/60) * angle(Frame2["a"].parameters[1]))
            new1_radius1 = (((60 - i)/60) * abs(Frame1["a"].parameters[1])) + ((i/60) * abs(Frame2["a"].parameters[1]))
            new1_angle2 = ((((60 - i)/60) * angle(Frame1["a"].parameters[2])) + ((i/60) * angle(Frame2["a"].parameters[2])))
            new1_angle2 = new1_angle2 - (2*pi)
            new1_radius2 = (((60 - i)/60) * abs(Frame1["a"].parameters[2])) + ((i/60) * abs(Frame2["a"].parameters[2]))
            polar1 = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex( new1_radius1 * cos(new1_angle1) + 0.2, new1_radius1 * sin(new1_angle1)), complex( new1_radius2 * cos(new1_angle2 - 2 *(new1_angle2 + deg2rad(20))) - 0.2, new1_radius2 * sin(new1_angle2 - 2 *(new1_angle2 + deg2rad(20))))])

            polar2 = Disk(colorant"rgba(0,0,0,0)", [0.7], [(((60 - i)/60) * Frame1["a"].parameters[1]) + ((i/60) * Frame2["a"].parameters[1]), (((60 - i)/60) * Frame1["a"].parameters[2]) + ((i/60) * Frame2["a"].parameters[2])])
            polar2 = Disk(colorant"rgba(0,0,0,0)", [0.7], [complex( abs(polar2.parameters[1]) * cos(angle(polar2.parameters[1])) + 0.2, abs(polar2.parameters[1]) * sin(angle(polar2.parameters[1]))), complex( abs(polar2.parameters[2]) * cos(angle(polar2.parameters[2])) - 0.2, abs(polar2.parameters[2]) * sin(angle(polar2.parameters[2])))])

            Transition["a"] = Disk(colorant"rgba(0,0,0,0)", [0.7], [( ((60-j)/60) * polar2.parameters[1] ) + ( ((j)/60) * polar1.parameters[1] ), ( ((60-j)/60) * polar2.parameters[2] ) + ( ((j)/60) * polar1.parameters[2] )])

            frame = compose(context(), disk_compose_single_base(Meta.parse(expression), Transition, 0), frame)
            #draw(PNG(file_path, 10cm, 10cm, dpi=250), frame)
        end
        file_path = string("Saved Images/Disk Operad/", string(j + 60), ".png")
        draw(PNG(file_path, 10cm, 10cm, dpi=250), frame)
    end

    anim = @animate for i in 1:120
        file_path = string("Saved Images/Disk Operad/", string(i), ".png")
        image = FileIO.load(file_path)
        plot(image, axis = nothing, background_color=:transparent)
    end

    gif(anim, "Saved Images/Rotating Cicles.gif", fps = 30)

end


text:
# Define the default algorithms for various pieces of X-ray physics data.
# If you want to change the default data that NeXL uses, make the change in here.

import Unitful: @u_str, ustrip
import PhysicalConstants.CODATA2018: PlanckConstant, SpeedOfLightInVacuum

const plancksConstant = ustrip(PlanckConstant |> u"eV*s")
const hc = ustrip((PlanckConstant * SpeedOfLightInVacuum) |> u"eV*cm") # (plancks⋅speed-of-light)
const speedOfLight = ustrip(SpeedOfLightInVacuum |> u"cm/s")

"""
    edgeenergy(z::Int, ss::Int, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64
    edgeenergy(cxr::CharXRay, ::Type{<:NeXLAlgorithm}=FFASTDB)

Return the minimum energy (in eV) necessary to ionize the specified sub-shell in the specified atom
or the ionized shell for the specified characteristic X-ray.
"""
edgeenergy(z::Int, ss::Int)::Float64 =
    edgeenergy(z, ss, FFASTDB)
edgeenergy(cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB) =
    edgeenergy(cxr.z, cxr.transition.innershell.index, alg)

hasedge(z::Int, ss::Int) = hasedge(z, ss, FFASTDB)

"""
     energy(ass::AtomicSubShell, ty::Type{<:NeXLAlgorithm}=FFASTDB)
     energy(elm::Element, ss::SubShell, ty::Type{<:NeXLAlgorithm}=FFASTDB)

 The edge energy in eV for the specified AtomicSubShell

Example:

    julia> energy(n"Fe L3")
    708.0999999999999
    julia> energy(n"Fe", n"L3")
    708.0999999999999

    energy(elm::Element, tr::Transition)::Float64
    energy(cxr::CharXRay)

The characteristic X-ray energy for the specified element / transition or characteristic X-ray.
"""
energy(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    edgeenergy(ass.z, ass.subshell.index, alg)
energy(elm::Element, ss::SubShell, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    edgeenergy(z(elm), ss.index, alg)
energy(elm::Element, tr::Transition, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    energy(z(elm), tr.innershell.index, tr.outershell.index, alg)
energy(cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB) = energy(
        cxr.z,
        cxr.transition.innershell.index,
        cxr.transition.outershell.index,
        alg
    )

"""
    eachelement()

Return the range of atomic numbers for which there is a complete set of energy, weight, MAC, ... data
"""
eachelement

let allelements = NTuple{length(FFAST.eachelement()), Element}( elements[FFAST.eachelement()] )
    global eachelement() = allelements
end

"""
    subshellindices(z::Int, alg::Type{<:NeXLAlgorithm} = FFASTDB)

Return the shells occupied in a neutral, ground state atom of the specified atomic number.
"""
subshellindices(z::Int) = subshellindices(z, FFASTDB)

"""
    jumpratio(z::Int, ss::Int) =

Compute the jump ratio.
"""
jumpratio(z::Int, ss::Int) = jumpratio(z::Int, ss::Int, FFASTDB)

"""
    mac(elm::Element, energy::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64
    mac(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64

The mass absorption coefficient for an X-ray of the specified energy (eV) or
characteristic X-ray line in the specified element.  In cm²/g.
"""
mac(elm::Element, energy::Float64)::Float64 =
    mac(elm, energy, FFASTDB)

const userMacs = Dict{Tuple{Element, CharXRay}, Float64}()
mac(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm}=FFASTDB)::Float64 =
    get(userMacs, (elm, cxr), mac(elm, energy(cxr), alg))
set_user_mac!(elm::Element, cxr::CharXRay, mac::Float64) =
    userMacs[(elm, cxr)]=mac
delete_user_mac!(elm::Element, cxr::CharXRay) =
    delete!(userMacs, (elm, cxr))
clear_user_macs!() =
    empty!(userMacs)

"""
    macU(elm::Element, energy::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)
    macU(elm::Element, cxr::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue
    macU(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue

The mass absorption coefficient (with uncertainty estimate) for an X-ray of the specified energy (eV)
or characteristix X-ray line in the specified element.
"""
macU(elm::Element, energy::Float64)::UncertainValue =
    macU(elm, energy, FFASTDB)
macU(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue =
    macU(elm, energy(cxr), alg)

"""
    characteristicXRayAvailable(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64

Is the weight associated with this transition greater than zero?
"""
charactericXRayAvailable(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Bool =
    isAvailable(z, inner, outer, alg)

"""
    ionizationcrosssection(ass::AtomicSubShell, energy::AbstractFloat, ty::Type{<:NeXLAlgorithm}=Bote2009)
    ionizationcrosssection(z::Int, ss::Int, energy::AbstractFloat)

Computes the absolute ionization crosssection (in cm²) for the specified AtomicSubShell and
electon energy (in eV) using the default algorithm.

Example:

    julia> (/)(map(e->NeXLCore.ionizationcrosssection(n"Fe K",e),[10.0e3,20.0e3])...)
    0.5672910174711278
"""
ionizationcrosssection(
    ass::AtomicSubShell,
    energy::AbstractFloat,
    ty::Type{<:NeXLAlgorithm} = Bote2009,
) = ionizationcrosssection(ass.z, ass.subshell.index, energy, ty)
ionizationcrosssection(z::Int, ss::Int, energy::AbstractFloat) =
    ionizationcrosssection(z, ss, energy, Bote2009)

"""
Represents the fractional number of X-rays emitted following the ionization of the sub-shell `ionized` via
the characteristic X-ray `z inner-outer`.  Due to cascades, `inner` does not necessarily equal `ionized`.
The `ionized` subshell may transition to a valency in `inner` via a combination of Auger, fluorescence or
Koster-Kronig transitions.  The various different forms make assumptions about the relationship between
`ionized` and `inner`, and about `outer`.

    fluorescenceyield(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL)::Float64

The fraction of relaxations from the specified shell that relax via any radiative transition. (`inner`==`ionized`)

    fluorescenceyield(cxr::CharXRay, alg::Type{<:NeXLAlgorithm}=CullenEADL)

The fraction of ionizations of `inner(cxr)` that relax via the one path `cxr`. `ionized==inner` && outer(cxr)

    fluorescenceyield(ash::AtomicSubShell, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64

The fractional number of `cxr` X-rays emitted (on average) for each ionization of `ash`.  This makes no
assumptions about `inner`, `outer` and `ionized`
"""
fluorescenceyield(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL) =
    fluorescenceyield(ass.z, ass.subshell.index, alg)
fluorescenceyield(cxr::CharXRay, alg::Type{<:NeXLAlgorithm}=CullenEADL) =
    fluorescenceyield(cxr.z, cxr.transition.innershell.index, cxr.transition.outershell.index, alg)
fluorescenceyield(ash::AtomicSubShell, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = CullenEADL) =
    ash.z == cxr.z ? fluorescenceyield(ash.z, ash.subshell.index, cxr.inner.index, cxr.outer.index, alg) : 0.0


text:
"""
    delete_const_columns!(table::CleanTable)

Removes in-place from a `CleanTable` each column filled with just a constant value.
"""
function delete_const_columns!(table::CleanTable)
    columns = cols(table)
    ndel = 0

    if Threads.nthreads() > 1 && length(columns) > 1 && length(columns[1]) >= 1_000_000
        col_state = replace!(Vector{Bool}(undef, length(columns)), true => false)

        Threads.@threads for i in 1:length(columns)
            if _is_const_col(columns[i])
                col_state[i] = true
            end
        end

        to_delete = findall(col_state)
        for j in to_delete
            deleteat!(columns, j - ndel)
            deleteat!(names(table), j - ndel)
            ndel += 1
        end
    else
        for i in 1:length(columns)
            if _is_const_col(columns[i - ndel])
                deleteat!(columns, i - ndel)
                deleteat!(names(table), i - ndel)
                ndel += 1
            end
        end
    end

    return table
end

"""
    delete_const_columns(table)

Creates a `CleanTable` with copied columns and removes each column filled with just a
constant value.
"""
function delete_const_columns(table)
    return delete_const_columns!(CleanTable(table))
end

function _is_const_col(col)
    known = []
    missing_found = false

    for el in col
        if !ismissing(el)
            if !in(el, known)
                push!(known, el)
                length(known) + missing_found >= 2 && return false
            end
        else
            length(known) >= 1 && return false

            if !missing_found
                missing_found = true
            end
        end
    end

    return true
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
module QuartetNetworkGoodnessFit

using DataFrames
using Distributed
using NLopt
using PhyloNetworks
using Random: seed!
using SharedArrays
using SpecialFunctions: loggamma
using Statistics: mean, median
using StatsFuns: normccdf, chisqccdf, betacdf, betaccdf

export
quarnetGoFtest!,
ticr,
ticr!

include("utils.jl")
include("ticr.jl")
include("quarnetGoF.jl")

# check hybrid-Lambda was downloaded (during build) and define variable hybridlambda
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
function sethybridlambdapath(path::String) global hybridlambda = path; end
warnmsg = """Functions using hybrid-Lambda for simulations will not work. Run:
    using Pkg; Pkg.build("QuartetNetworkGoodnessFit") then restart Julia and try again.
    If the error persists, install the Hybrid-Lambda dependency manually and run:
    QuartetNetworkGoodnessFit.sethybridlambdapath("path_to_your_hybrid-Lambda_executable")
    """
if !isfile(depsjl_path)
    @warn "QuartetNetworkGoodnessFit not installed properly: no deps.jl file.\n" * warnmsg
else
    include(depsjl_path)
    try # check dependencies. `check_deps` defined in `deps.jl`
        check_deps()
    catch e
        @warn "QuartetNetworkGoodnessFit not installed properly: hybrid-Lambda executable not found.\n" * warnmsg
    end
end
end # module


text:
# TODO: implement CUDA kernels for below

function leaky_relu(x; alpha=0.2)
    neg = min(0, x) * alpha
    pos = max(0, x)
    return neg + pos
end

# TODO: merge this with `special.jl`

# import StatsFuns: log1pexp, logexpm1

# # This function is adapted from
# # https://github.com/denizyuret/Knet.jl/blob/687dba214c3f3326272c7c2191e14e8fdbc35c6a/src/unary.jl#L7-L28
# function unary_op(f)
#     J = Symbol(f)
#     M = which(@__MODULE__, J)
#     @eval begin
#         function Base.Broadcast.broadcasted(::typeof($J), x::Knet.KnetArray{T}) where {T<:AbstractFloat}
#             y = similar(x)
#             @inbounds for i=1:length(y)
#                 xi = x[i]
#                 y[i] = $J(xi)
#             end
#             return y
#         end
#         # TODO: see if I can implement gradient as well; ref:
#         # https://github.com/denizyuret/Knet.jl/blob/687dba214c3f3326272c7c2191e14e8fdbc35c6a/src/unary.jl#L62-L70
#         # Bcasted methods
#         ($M).$J(x::Knet.Bcasted{<:Knet.KnetArray{T}}) where {T<:AbstractFloat} = Base.Broadcast.broadcasted($J, x.value) |> Knet.Bcasted
#         Base.Broadcast.broadcasted(::typeof($J), x::Knet.Bcasted{<:Knet.KnetArray{T}}) where {T<:AbstractFloat} = Base.Broadcast.broadcasted($J, x.value) |> Knet.Bcasted
#     end
#     @eval begin # so we do not trigger some default Base implementation
#         ($M).$J(x::Knet.Bcasted) = throw(MethodError($J,(x,)))
#         Base.Broadcast.broadcasted(::typeof($J), x::Knet.Bcasted) = throw(MethodError($J,(x,)))
#     end
# end

# for f in [:log1pexp, :logexpm1]
#     unary_op(f)
# end

# Knet.AutoGrad.@primitive log1pexp(x::Array),dy,y  dy.*exp.(x .- y)
# Knet.AutoGrad.@primitive log1pexp(x::Knet.KnetArray),dy,y  dy.*exp.(x .- y)

# const softplus = log1pexp
# const invsoftplus = logexpm1


text:
export real_requests

import Requests

type RequestsHttpResponse
    r::Requests.Response
end

statuscode(resp::RequestsHttpResponse) = Requests.statuscode(resp.r)
text(resp::RequestsHttpResponse) = Requests.text(resp.r)

type RequestsHttp <: AbstractHttp end

post(::RequestsHttp, uri::AbstractString; args...) =
    RequestsHttpResponse(Requests.post(Requests.URI(uri); args...))

real_requests = RequestsHttp()


text:
struct convGauss{T,S} <: AbstractPDF{1}
    source::T
    σ::S
end

pars(d::convGauss, isfree::Bool) = pars(d.source, isfree) + pars(d.σ, isfree)
lims(d::convGauss) = lims(d.source)
updatevalueorflag(d::convGauss, s::Symbol, isfree::Bool, v=getproperty(pars(d),s)) =
    convGauss(
        ispar(d.source,s) ? updatevalueorflag(d.source,s,isfree,v) : d.source,
        ispar(d.σ,     s) ? updatevalueorflag(d.σ,     s,isfree,v) : d.σ)

function func(d::convGauss, x::NumberOrTuple; p=pars(d))
    σ = func(d.σ, x; p)
    g(z) = AlgebraPDF.standardgauss(z, σ)
    f(z) = func(d.source, z; p)
    return quadgk(y->f(x-y) * g(y), -5*σ, +5*σ)[1]
end


text:
@testset "linear objective function tests" begin
    PR.silence!()
    f, partition = a -> a^3, collect(-1.0:0.25:1.0)
    ufd = UnivariateFunctionData(
        f,
        a -> 3 * a^2,
        partition,
        NaN64,
        1e-3,
        1e-3,
        0,
        length(partition),
    )

    # create LP relaxation of the MILP relaxation of the univariate function.
    milp = Model(cbc_optimizer)
    @variable(milp, -1.0 <= x <= 1.0)
    @variable(milp, y)
    construct_univariate_relaxation!(milp, f, x, y, partition, true)

    # create LP relaxation of the univariate function using the convex hull formulation.
    lp = Model(cbc_optimizer)
    @variable(lp, -1.0 <= x_lp <= 1.0)
    @variable(lp, y_lp)
    construct_univariate_relaxation!(lp, f, x_lp, y_lp, partition, false)

    # univariate function data

    # Collect possible solutions of relaxations.
    sec_verts, tan_verts = PR._collect_vertices(ufd)
    sln_verts = PR.Vertex2d[]
    push!(sln_verts, sec_verts[1])
    append!(sln_verts, tan_verts)
    push!(sln_verts, sec_verts[end])

    tol = 1e-5

    for α in collect(π/20:π/20:(2*π/5))
        set_objective_function(milp, y - x * tan(α))
        set_objective_function(lp, y_lp - x_lp * tan(α))

        for s in [MOI.MIN_SENSE, MOI.MAX_SENSE]
            # Solve MILP relaxation.
            set_objective_sense(milp, s)
            optimize!(milp)
            @test termination_status(milp) == MOI.OPTIMAL

            # Verify that the solution is one of the candidate solutions in `sln_verts`.
            milp_obj = objective_value(milp)
            sln_found = false
            x_sln = value.(x)
            y_sln = value.(y)
            for v in sln_verts
                if isapprox(v[1], x_sln, atol = tol) &&
                   isapprox(v[2], y_sln, atol = tol)
                    sln_found = true
                    break
                end
            end
            @test sln_found == true

            # Verify that the LP relaxation obtained by dropping binary restrictions gives the
            # same objective value as the MILP relaxation.
            JuMP.relax_integrality(milp)
            optimize!(milp)
            @test termination_status(milp) == MOI.OPTIMAL
            lp1_obj = objective_value(milp)
            @test milp_obj ≈ lp1_obj atol = tol

            # Solve LP relaxation formulated in the convex hull form.
            # Verify that the LP relaxation formulated in the convex hull form has the same
            # objective as the above 2 relaxations.
            set_objective_sense(lp, s)
            optimize!(lp)
            @test lp1_obj ≈ objective_value(lp) atol = tol
        end
    end
end


text:
using QuComputStates
using Base.Test

t = ComputState([1,2,3,4])
@test rawcoeffs(t) == [1,2,3,4]

# write your own tests here
@test 1 == 1


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

# author: Travis Askham (askhamwhat@gmail.com)
#
#TODO: there are a number of stumbling points for the
# FAST_KER=ON option when cross-compiling. it is completely
# disabled for now. Some issues:
# - on musl, darwin, and freebsd various standard libraries
#   are not found
# - on mingw there is an invalid regiser error on x86_64 arch
# - generally, these seem to require a later version of gcc
#   but we would like to compile for various libgfortran
#   major versions

using BinaryBuilder, Pkg

name = "FMM3D"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/flatironinstitute/FMM3D.git", "2fd129984e48abe437fddd12013dc434563a773a")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd FMM3D/

touch make.inc

# clang has some issue with openmp (lomp not found)
# force gcc
if [[ ${target} = *apple* || ${target} = *freebsd* ]]; then
  export CC="gcc"
fi

echo "CC=${CC}" >> make.inc;
echo "CXX=${CXX}" >> make.inc;
echo "FC=${FC}" >> make.inc;

if [[ ${target} = *mingw* ]]; then
    echo "FFLAGS= -fPIC -O3 -fno-asynchronous-unwind-tables -funroll-loops" >> make.inc;
else
    echo "FFLAGS= -fPIC -O3 -funroll-loops" >> make.inc;
fi

export OMPFLAGS="-fopenmp"
export OMPLIBS="-lgomp"
echo "OMPFLAGS= ${OMPFLAGS}" >> make.inc;
echo "OMPLIBS= ${OMPLIBS}" >> make.inc;

export LIBS="-lm -lgfortran"
echo "LIBS = ${LIBS}" >> make.inc;

make -j${nproc} lib OMP=ON

echo "${TARGET}"

# build the correct type of dynamic library for
# current target

export SHAREFLAGS="-shared -fPIC"

cd lib-static
ar x libfmm3d.a
${CC} ${SHAREFLAGS} ${OMPFLAGS} *.o -o "${libdir}/libfmm3d.${dlext}" ${LIBS} ${OMPLIBS}
cd ..

install_license ${WORKSPACE}/srcdir/FMM3D/LICENSE
rm make.inc
exit
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using SDLGamesLib
using Test

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

@testset "anim.jl" begin
    include("anim.jl")
end


text:
include(Pkg.dir("MethylUtils", "src","Features.jl"))
include(Pkg.dir("MethylUtils", "src","Functions.jl"))

using Distributions

function rand_CytosineCounts(num_cgs, max_cov, min_cov=0)
   coverage_distribution = [ x for x= 1:max_cov]
   cgs = Array(CytosineCount,num_cgs)
   k=1
   while k <= num_cgs
       (Cs, Ts) = rand_c_counts(0.8, coverage_distribution)
       if ( Cs + Ts ) >= min_cov
           cgs[k] = CytosineCount(SequencePos(1,2),true,int16(Cs) ,int16(Ts))
           k +=1
       end
   end
   return cgs
end

function rand_c_counts(prob, coverage_distribution)
    num_reads = sample(coverage_distribution)
    binom = Binomial(num_reads,prob)
    Cs = rand(binom)
    Ts = num_reads - Cs
    return (Cs,Ts)
end

SF_B = SequenceFeatures(Dict(),Dict())
for k=1:10
    println( "Making CytsonineCounts" )
    rc = rand_CytosineCounts(100,100,5)
    SF_B.features["chr$k"] = rc
end

SF_O = SequenceFeatures(Dict(),Dict())
for k=1:10
    println( "Making CytsonineCounts" )
    rc = rand_CytosineCounts(100,100,5)
    SF_O.features["chr$k"] = rc
end

meth_c = run_mlml(SF_B,SF_O)
println(meth_c)


text:
using Documenter, ExponentialUtilities

include("pages.jl")

makedocs(
    sitename="ExponentialUtilities.jl",
    authors="Chris Rackauckas",
    modules=[ExponentialUtilities],
    clean=true,doctest=false,
    format = Documenter.HTML(analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://exponentialutilities.sciml.ai/stable/"),
    pages=pages
)

deploydocs(
   repo = "github.com/SciML/ExponentialUtilities.jl.git";
   push_preview = true
)

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using Documenter
using AdaptiveResonance
using DemoCards
# using JSON

# Generate the demo files
# this is the relative path to docs/
demopage, postprocess_cb, demo_assets = makedemos("examples")

assets = [
    joinpath("assets", "favicon.ico")
]

# if there are generated css assets, pass it to Documenter.HTML
isnothing(demo_assets) || (push!(assets, demo_assets))

# Make the documentation
makedocs(
    modules=[AdaptiveResonance],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = assets,
    ),
    pages=[
        "Home" => "index.md",
        "Getting Started" => [
            "getting-started/whatisart.md",
            "getting-started/basic-example.md",
        ],
        "Tutorial" => [
            "Guide" => "man/guide.md",
            demopage,
            # "Examples" => "man/examples.md",
            "Modules" => "man/modules.md",
            "Contributing" => "man/contributing.md",
            "Index" => "man/full-index.md",
        ],
    ],
    repo="https://github.com/AP6YC/AdaptiveResonance.jl/blob/{commit}{path}#L{line}",
    sitename="AdaptiveResonance.jl",
    authors="Sasha Petrenko",
    # assets=String[],
)

# 3. postprocess after makedocs
postprocess_cb()

# a workdaround to github action that only push preview when PR has "push_preview" labels
# issue: https://github.com/JuliaDocs/Documenter.jl/issues/1225
# function should_push_preview(event_path = get(ENV, "GITHUB_EVENT_PATH", nothing))
#     event_path === nothing && return false
#     event = JSON.parsefile(event_path)
#     haskey(event, "pull_request") || return false
#     labels = [x["name"] for x in event["pull_request"]["labels"]]
#     return "push_preview" in labels
#  end

deploydocs(
    repo="github.com/AP6YC/AdaptiveResonance.jl.git",
    devbranch="develop",
    # push_preview = should_push_preview(),
)


text:
const sources = Dict{Symbol, Any}()
const actfuns = Dict{Symbol, Any}()
const rnnactfuns = Dict{Symbol, Any}() # Recurrent layers have activation functions as attributes and use different parameter names compared to their respective operations.
const actlayers = Dict{Symbol, Any}()
const fluxlayers = Dict{Symbol, Any}()
const fluxrecurrentlayers = Dict{Symbol, Any}()
const invariantops = Dict{Symbol, Any}()
const pseudotransparentops = Dict{Symbol, Any}()
const verts = Dict{Symbol, Any}()
const fluxlayertypes = Dict{Symbol, Any}()

# Rundown of the basic idea here:

# Aspect 1
# ONNX does not have activation functions as an attribute to its layers but rather represents them as a separate node
# This would indeed be workable, but...
# 1. It is a bit annoying that model -> serialize -> deserialize does not result in the exact same thing
# 2. If one wants to use the mutation functionality of NaiveNASflux it might not be desirable to have activation
#    functions as separate vertices in the graph as this invites for things like inserting something else between
#    the layer and its activation function.

# To be able to have activation functions back inside their layers when deserializing, whenever an op which is a key
# in actlayers is encountered there is a "lookahead" to see if the op of the next node is in actfuns. If it is, the
# two ops will be merged into one vertex containing the layer and its activation function.
# A very similar thing is done for global pooling operations followed by squeeze or reshape.

# Aspect 2
# The vertices of NaiveNASflux require a few inputs when creating them. One in particular is knowledge of the size
# trait which is obviously not possible to obtain from the ONNX data. In order to spare users from having to supply
#  this extra input with each operation there is one dict per "general type".

# As NaiveNASflux already has the knowledge what is needed for all layers in Flux, they have their own dict
#  (fluxlayers) which just outsources the vertex creation to NaiveNASflux. Note that all actlayers are inserted
# in this dict.

# Functions which always produce the same number of outputs as inputs and are not defined in Flux, e.g.
#  GlobalAveragePool end up in invariantops.

# Functions which have dedicated vertex construction methods, such as Concat and Add end up in verts.

sources[:Constant] = params -> constant(Val.(keys(params))..., values(params)...)
constant(::Val{:value}, val::ONNX.TensorProto) = val |> array
constant(::Val{:value}, val) = val

actfuns[:Relu] = params -> Flux.relu

actfuns[:Elu] = function(params)
    α = get(params, :alpha, 1)
    return x -> Flux.elu(x, oftype(x, α))
end
rnnactfuns[:Elu] = (ind, params) -> actfuns[:Elu](Dict(:alpha => get(params, :activation_alpha, ntuple(i -> 1, ind))[ind]))

actfuns[:Selu] = function(params)
    haskey(params, :alpha) || haskey(params, :gamma) && return Flux.selu
    γ = get(params, :gamma, Float32(1.05070102214813232421875))
    α = get(params, :alpha, Float32(1.67326319217681884765625))
    return x -> selu(x, oftype(x, γ), oftype(x, α))
end
Flux.selu(x, γ, α) = γ * ifelse(x > 0, x/1, α * (exp(x) - 1))

actfuns[:Tanh] = params -> tanh
rnnactfuns[:Tanh] = (ind, params) -> tanh

mrev(x) = x
mrev(x::AbstractVector) = reverse(x)
prev(x) = x
prev(x::AbstractVector) = reshape(permutedims(reverse(reshape(x, length(x) ÷ 2,:);dims=1)),:)

# mrev = maybe reverse. prev = rearrange padding, e.g. (1,2,1,2) => (2,2,1,1) or (1,2,3,1,2,3) => (3,3,2,2,1,1)
_akpsd(params) = get(params, :activation, identity), mrev(get(params, :kernel_shape, 1)), prev(get(params, :pads, 0)), mrev(get(params, :strides, 1)), mrev(get(params, :dilations, 1))
akpsd(params) = a2t.(_akpsd(params))
a2t(x) = x
a2t(a::AbstractArray) = Tuple(a)

actlayers[:Conv] = function(params, weight::AbstractArray{T, N}, bias=Flux.Zeros()) where {T, N}
    a,_,p,s,d = akpsd(params)
    @assert get(params, :group, 1) == 1 "Group size not supported!" #Or?
    return Conv(flipweights(FluxConv{N-2}(), weight), bias, a, pad=p, stride=s, dilation=d)
end
fluxlayertypes[:Conv] = (weight, bias=nothing) -> FluxConv{length(size(weight))-2}()

biasarray(b::Flux.Zeros, esize, β) = b
biasarray(b::AbstractArray, esize, β) = length(b) === 1 ? repeat(β .* vec(b), esize) : β .* reshape(b, :)
biasarray(b::Number, esize, β) = repeat([β * b], esize)

actlayers[:Gemm] = function(params, weight::AbstractArray{T, N}, bias=Flux.Zeros()) where {T,N}
    act = get(params, :activation, identity)
    wt = Bool(get(params, :transB, 0)) ? permutedims : identity
    α = get(params, :alpha, 1)
    β = get(params, :beta, 1)

    weight = α .* wt(weight)
    bias = biasarray(bias, size(weight, 1), β)

    return Dense(weight , bias, act)
end
fluxlayertypes[:Gemm] = (pars...) -> FluxDense()

actlayers[:BatchNormalization] = function(params, γ, β, μ, σ²)
    λ = get(params, :activation, identity)
    ϵ = get(params, :epsilon, 1f-5)
    momentum = get(params, :momentum, 9f-1)

    return BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, true, true, nothing, length(γ))
end
fluxlayertypes[:BatchNormalization] = (pars...) -> FluxBatchNorm()

default_Wb_Rb(Wh_WBh) = fill!(similar(Wh_WBh, (size(Wh_WBh, 2) * 2, size(Wh_WBh, 3))), 0)
default_init_h(Wb_Rb, sc) = fill!(similar(Wb_Rb, (size(Wb_Rb,1) ÷ sc, size(Wb_Rb,2))), 0)
# TODO when https://github.com/FluxML/Flux.jl/issues/1279 is resolved default_init_h(Wh_WBh, sc) = fill!(similar(Wh_WBh, (size(Wh_WBh, 2) ÷ sc, size(Wh_WBh, 3))), 0)

fluxrecurrentlayers[:RNN] = function(params, Wi_WBi, Wh_WBh, Wb_Rb=default_Wb_Rb(Wh_WBh), seqlen=[], h3d = default_init_h(Wb_Rb, 2))
    @assert size(Wi_WBi, 3) == 1 "Num directions must be 1! Bidirectional (num directions = 2) not supported!" # TODO: Add...

    Wi,Wh,b,h = recurrent_arrays(FluxRnn(), Wi_WBi, Wh_WBh, Wb_Rb, h3d)
    act = rnnactfuns[Symbol(get(params, :activations, ["Tanh"])[])](1, params)
    cell = Flux.RNNCell(act, Wi, Wh, b, fill!(similar(h), 0))
    return Flux.Recur(cell, h)
end
fluxlayertypes[:RNN] = (pars...) -> FluxRnn()

fluxrecurrentlayers[:LSTM] = function(params, Wi_WBi, Wh_WBh, Wb_Rb=default_Wb_Rb(Wh_WBh), seqlen=[1], h3d = default_init_h(Wb_Rb, 8), c3d=default_init_h(Wb_Rb,8), peep=nothing)
    @assert size(Wi_WBi, 3) == 1 "Num directions must be 1! Bidirectional (num directions = 2) not supported!" # TODO: Add...
    @assert isnothing(peep) "Peepholes not supported!" # Or?
    Wi,Wh,b,h,c = recurrent_arrays(FluxLstm(), Wi_WBi, Wh_WBh, Wb_Rb, h3d, c3d)
    # Flux only supports default activation functions
    # We can only check that given values doesn't deviate
    supported = [:Sigmoid, :Tanh, :Tanh]
    acts = get(params, :activations, supported)
    @assert all(zip(supported, acts)) do (e,a)
        e == a
    end "Got unsupported activation function: $acts"

    # b, h and c must all be of the same type when creating a cell, but
    # it is actually Recur which has the state
    cell = Flux.LSTMCell(Wi, Wh, b, (fill!(similar(h), 0), fill!(similar(c), 0)))
    return Flux.Recur(cell, (h, c))
end
fluxlayertypes[:LSTM] = (pars...) -> FluxLstm()

function recurrent_arrays(lt, Wi_WBi, Wh_WBh, Wb_Rb, h3ds...)
    # ONNX weights are on the form [num_directions, hidden_size, input_size] (where num_directions is 2 for bidirectional else 1)
    # Flux weights are of shape [hidden_size, input_size]
    # To spice things up a bit, all julia arrays are loaded in reverse order, i.e we get an array with the arrangement [input_size, hidden_size, num_directions].
    # First remove the num_directions dimension, then transpose into the correct shape
    hsize = size(Wh_WBh, 1)
    Wi = unflipweights(lt, permutedims(dropdims(Wi_WBi, dims=3)), hsize)
    Wh = unflipweights(lt, permutedims(dropdims(Wh_WBh, dims=3)), hsize)
    b = Wb_Rb isa Flux.Zeros ? Wb_Rb : dropdims(unflipweights(lt, sum(reshape(Wb_Rb, :, 2), dims=2), hsize),dims=2)
    return Wi, Wh, b, h3ds...
end

fluxlayers[:MaxPool] = function(params)
    _,k,p,s,_ = akpsd(params)
    return MaxPool(k, pad=p, stride=s)
end
fluxlayertypes[:MaxPool] = (pars...) -> FluxPoolLayer()

fluxlayers[:AveragePool] = function(params)
    _,k,p,s,_ = akpsd(params)
    return MeanPool(k, pad=p, stride=s)
end
fluxlayertypes[:AveragePool] = (pars...) -> FluxPoolLayer()

fluxlayers[:Dropout] = params -> Dropout(get(params, :ratio, 0.5))
fluxlayertypes[:Dropout] = (pars...) -> FluxDropOut()

invariantops[:GlobalAveragePool] = function(params)
    wrap = get(params, :wrap, identity)
    return x -> globalmeanpool(x, wrap)
end
fluxlayertypes[:GlobalAveragePool] = (pars...) -> FluxPoolLayer()

function globalmeanpool(x::AbstractArray{T,N}, wrap) where T where N
    wrap(MeanPool(size(x)[1:N-2])(x))
end

invariantops[:GlobalMaxPool] = function(params)
    wrap = get(params, :wrap, identity)
    return x -> globalmaxpool(x, wrap)
end
fluxlayertypes[:GlobalMaxPool] = (pars...) -> FluxPoolLayer()

function globalmaxpool(x::AbstractArray{T,N}, wrap) where T where N
    wrap(MaxPool(size(x)[1:N-2])(x))
end

invariantops[:Squeeze] = function(params)
    np_axes = get(params, :axes, missing)
    dimfun = ismissing(np_axes) ? x -> Tuple(findall(i -> i == 1, size(x))) : x -> Tuple(numpy2fluxdim.(np_axes, ndims(x)))
    return x -> dropdims(x, dims=dimfun(x))
end

invariantops[:ReduceMean] = function(params)
    np_axes = get(params, :axes, missing)
    keepdims = Bool(get(params, :keepdims, 1))

    dimexp =
    if keepdims && ismissing(np_axes)
        # As mean returns a scalar when no dimensions are provided
        expanddims
    elseif !keepdims
        (out, x, dims) -> dropdims(out, dims=dims)
    else
        (out, x, dims) -> out
    end

    ismissing(np_axes) && return x -> dimexp(mean(x), x, missing)

    return function(x)
        dims = Tuple(numpy2fluxdim.(np_axes, ndims(x)))
        out = mean(x, dims=dims)
        return dimexp(out, x, dims)
    end
end
expanddims(out, x, dims) = fill(out, ntuple(i -> 1, ndims(x)))

invariantops[:Softmax] = params -> x -> onnxsoftmax(x; np_axis = get(params, :axis, 1))

function onnxsoftmax(x::AbstractArray{T, 2}; np_axis=1) where T
    dim = numpy2fluxdim(np_axis, 2)
    Flux.softmax(x; dims=dim)
end
function onnxsoftmax(x::AbstractArray{T, N}; np_axis=1) where {T,N}
    dim = numpy2fluxdim(np_axis, N)
    sz = size(x)
    reshape(Flux.softmax(reshape(x, prod(sz[1:dim]), :)), sz...)
end

pseudotransparentops[:Reshape] = function(params, shape)
    shape_t = Tuple(reverse(replace(shape, -1 => Colon())))
    return MeasureNout(Reshape(shape_t))
end

pseudotransparentops[:Flatten] = function(params)
    dim = -get(params,:axis, 1)
    return MeasureNout(Flatten(dim))
end

verts[:Input] = function(name, inputs, params; kwargs...)
    inshape = params[:size]
    ltype = params[:ltype]
    indims = length(inshape)

    insize = indims > 0 ? inshape[max(1, actdim(ltype))] : 1 # assume scalar
    return inputvertex(name, insize, ltype)
end

verts[:Add] = (name, inputs, params; kwargs...) -> elemwisevertex(name, inputs, params, +, 0; kwargs...)
verts[:Mul] = (name, inputs, params; kwargs...) -> elemwisevertex(name, inputs, params, *, 1; kwargs...)

function elemwisevertex(name, inputs, params, op, id; traitdecoration=identity, layerfun=identity, kwargs...)
    c = reduce((c1,c2) -> op.(c1, c2), get(params, :Constant, id))
    c = length(c) == 1 ? c[] : c
    let cc = c
        opp, wrap = cc == id ? (op, layerfun) : (identity, f -> layerfun((x...) -> op.(cc, x...)))
        conf = VertexConf(traitdecoration = named(name) ∘ traitdecoration, outwrap = wrap, kwargs...)
        return NaiveNASlib.elemwise(opp, conf, inputs...)
    end
end

verts[:Concat] = function(name, inputs, params; traitdecoration=identity, layerfun=identity, kwargs...)
    dims = numpy2fluxdim(params[:axis], inputs[1])
    return conc(inputs..., dims=dims, traitdecoration = named(name) ∘ traitdecoration, outwrap=layerfun, kwargs...)
end

function refresh()
    for (s, f) in actlayers
        fluxlayers[s] = f
    end

    for (s, f) in fluxrecurrentlayers
        fluxlayers[s] = f
    end

    for (s, f) in actfuns
        invariantops[s] = function(args...;kwargs...)
                                actfun = f(args...; kwargs...)
                                return x -> actfun.(x)
                            end
    end

    for (s, f) in fluxlayers
        verts[s] = (name, inputs, args...;kwargs...) -> fluxvertex(name, f(args...), inputs...; kwargs...)
    end

    for (s, f) in invariantops
        verts[s] = (name, inputs, args...;traitdecoration=identity, layerfun=identity, kwargs...) -> invariantvertex(layerfun(f(args...)), inputs...; traitdecoration = named(name) ∘ traitdecoration, kwargs...)
    end

    for (s,f) in pseudotransparentops
        verts[s] = function(name, inputs, args...;traitdecoration=identity, layerfun=identity, kwargs...)
            comp = f(args...)
            return absorbvertex(layerfun(comp), inputs...; traitdecoration = named(name) ∘ traitdecoration ∘ SizePseudoTransparent, kwargs...)
        end
    end

    for (s,f) in sources
        verts[s] = function(name, inputs, args...;kwargs...)
            @assert isempty(inputs) "Source of type $s got inputs $(inputs)!"
            return sourcevertex_with_outputs(f(args...), name)
        end
    end

    for s in keys(verts)
        get!(fluxlayertypes, s, (args...) -> missing)
    end

end

refresh()

list_supported_ops(io::IO=stdout) = foreach(ot -> println(io, ot), filter(ot -> ot != :Input, sort(collect(keys(verts)))))


text:
@testset "activate" begin
    patch = @patch run(cmd::Base.AbstractCmd, args...) = println("Mocked: $cmd")

    Mocking.apply(patch) do
        activate(TEST_CONFIG, join(TEST_TMP_DIR, "test-playground"))
        activate(TEST_CONFIG, "myproject")
        activate(TEST_CONFIG)

        # Test activate when the SHELL environment isn't set.
        withenv("SHELL" => nothing) do
            activate(TEST_CONFIG)
        end
    end

    activate(Environment("myproject"); shell=false)
    deactivate()
end


text:
using DiffEqOperators

nx = 200
x = range(0,2pi,length=nx)
dx = x[2]-x[1]
y = sin.(x)
dy = [cos.(x),-sin.(x),-cos.(x),y,cos.(x),-sin.(x)]

for dor in 1:6, aor in 2:2:10

    D1 = CenteredDifference(dor,aor,dx,nx-2)
    #take derivative
    dyt = D1*y
    #test result
    @test dy[dor][2:end-1] ≈ dyt atol=1.0^(1-aor)

    #TODO: implement specific tests for the left and right boundary regions, waiting until after update
end


text:
@marginalrule GCV(:y_x) (m_y::UniNormalOrExpLinQuad, m_x::UniNormalOrExpLinQuad, q_z::Any, q_κ::Any, q_ω::Any) = begin

    y_mean, y_precision = mean_precision(m_y)
    x_mean, x_precision = mean_precision(m_x)

    z_mean, z_var = mean_var(q_z)
    κ_mean, κ_var = mean_var(q_κ)
    ω_mean, ω_var = mean_var(q_ω)

    ksi = κ_mean ^ 2 * z_var + z_mean ^ 2 * κ_var + z_var * κ_var
    A = exp(-ω_mean + ω_var / 2)
    B = exp(-κ_mean * z_mean + ksi / 2)
    W = [ y_precision + A * B -A * B; -A * B x_precision + A * B ]
    ξ = [ y_mean * y_precision; x_mean * x_precision ]

    return MvNormalWeightedMeanPrecision(ξ, W)
end


text:
__precompile__(true)

module SHA

# Export convenience functions, context types, update!() and digest!() functions
export sha1, SHA1_CTX, update!, digest!
export sha224, sha256, sha384, sha512
export sha2_224, sha2_256, sha2_384, sha2_512
export sha3_224, sha3_256, sha3_384, sha3_512
export SHA224_CTX, SHA256_CTX, SHA384_CTX, SHA512_CTX
export SHA2_224_CTX, SHA2_256_CTX, SHA2_384_CTX, SHA2_512_CTX
export SHA3_224_CTX, SHA3_256_CTX, SHA3_384_CTX, SHA3_512_CTX
export HMAC_CTX, hmac_sha1
export hmac_sha224, hmac_sha256, hmac_sha384, hmac_sha512
export hmac_sha2_224, hmac_sha2_256, hmac_sha2_384, hmac_sha2_512
export hmac_sha3_224, hmac_sha3_256, hmac_sha3_384, hmac_sha3_512

include("constants.jl")
include("types.jl")
include("base_functions.jl")
include("sha1.jl")
include("sha2.jl")
include("sha3.jl")
include("common.jl")
include("hmac.jl")

# Create data types and convenience functions for each hash implemented
for (f, ctx) in [(:sha1, :SHA1_CTX),
                 (:sha224, :SHA224_CTX),
                 (:sha256, :SHA256_CTX),
                 (:sha384, :SHA384_CTX),
                 (:sha512, :SHA512_CTX),
                 (:sha2_224, :SHA2_224_CTX),
                 (:sha2_256, :SHA2_256_CTX),
                 (:sha2_384, :SHA2_384_CTX),
                 (:sha2_512, :SHA2_512_CTX),
                 (:sha3_224, :SHA3_224_CTX),
                 (:sha3_256, :SHA3_256_CTX),
                 (:sha3_384, :SHA3_384_CTX),
                 (:sha3_512, :SHA3_512_CTX),]
    g = Symbol(:hmac_, f)

    @eval begin
        # Our basic function is to process arrays of bytes
        function $f(data::T) where T<:Union{Array{UInt8,1},NTuple{N,UInt8} where N}
            ctx = $ctx()
            update!(ctx, data)
            return digest!(ctx)
        end
        function $g(key::Vector{UInt8}, data::T) where T<:Union{Array{UInt8,1},NTuple{N,UInt8} where N}
            ctx = HMAC_CTX($ctx(), key)
            update!(ctx, data)
            return digest!(ctx)
        end

        # AbstractStrings are a pretty handy thing to be able to crunch through
        $f(str::AbstractString) = $f(Vector{UInt8}(codeunits(str)))
        $g(key::Vector{UInt8}, str::AbstractString) = $g(key, Vector{UInt8}(str))

        # Convenience function for IO devices, allows for things like:
        # open("test.txt") do f
        #     sha256(f)
        # done
        function $f(io::IO, chunk_size=4*1024)
            ctx = $ctx()
            buff = Vector{UInt8}(undef, chunk_size)
            while !eof(io)
                num_read = readbytes!(io, buff)
                update!(ctx, buff[1:num_read])
            end
            return digest!(ctx)
        end
        function $g(key::Vector{UInt8}, io::IO, chunk_size=4*1024)
            ctx = HMAC_CTX($ctx(), key)
            buff = Vector{UInt8}(chunk_size)
            while !eof(io)
                num_read = readbytes!(io, buff)
                update!(ctx, buff[1:num_read])
            end
            return digest!(ctx)
        end
    end
end

end #module SHA

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
primitive type Salsa512 512 end

Salsa512(x::Vector{UInt128}) = reinterpret(Salsa512, x) |> first

import SIMD.valloc
valloc(t::Type, n) = valloc(t, sizeof(t), n)

import Base.zero
zero(::Type{Salsa512}) = zeros(UInt128, 4) |> Salsa512


text:
const IVec8 = SVector{8, Int}

struct DSubHuntPOMDP <: POMDP{SubState, Int, IVec8}
    cp::SubHuntPOMDP
    binsize::Float64
end

struct DSubObsDist
    cd::SubObsDist
    binsize::Float64
end

rand(rng::AbstractRNG, d::DSubObsDist) = floor.(Int, (rand(rng, d.cd)./d.binsize)::Vec8)::IVec8

function pdf(d::DSubObsDist, o::IVec8)
    p = 1.0
    cd = d.cd
    for i in 1:length(o)
        lo = d.binsize * o[i]
        hi = lo + d.binsize
        if i == cd.abeam
            p *= cdf(cd.an, hi) - cdf(cd.an, lo)
        else
            p *= cdf(cd.n, hi) - cdf(cd.n, lo)
        end
    end
    return p
end

observation(p::DSubHuntPOMDP, a::Int, sp::SubState) = DSubObsDist(observation(p.cp, a, sp), p.binsize)

n_states(p::DSubHuntPOMDP) = n_states(p.cp)
state_index(p::DSubHuntPOMDP, s::SubState) = state_index(p.cp, s)
states(p::DSubHuntPOMDP) = states(p.cp)
initial_state_distribution(p::DSubHuntPOMDP) = initial_state_distribution(p.cp)
actions(p::DSubHuntPOMDP) = actions(p.cp)
n_actions(p::DSubHuntPOMDP) = n_actions(p.cp)
transition(p::DSubHuntPOMDP, s::SubState, a::Int) = transition(p.cp, s, a)
discount(p::DSubHuntPOMDP) = discount(p.cp)
isterminal(p::DSubHuntPOMDP, s::SubState) = isterminal(p.cp, s)
reward(p::DSubHuntPOMDP, s::SubState, a::Int, sp::SubState) = reward(p.cp, s, a, sp)


text:
abstract type GraphQLNode end

struct ListType <: GraphQLNode
    type::Any
    # true for trailing !, false otherwise
    non_null::Bool

    function ListType(type, non_null::Some{Nothing})
        new(type, false)
    end

    function ListType(type, non_null::Some{T}) where {T}
        new(type, true)
    end

    function ListType(type, non_null::Bool)
        new(type, non_null)
    end
end

struct NamedType <: GraphQLNode
    type::Any
    # true for !, false otherwise
    non_null::Bool

    function NamedType(type, non_null::Some{Nothing})
        new(type, false)
    end

    function NamedType(type, non_null::Some{T}) where {T}
        new(type, true)
    end

    function NamedType(type, non_null::Bool)
        new(type, non_null)
    end
end

struct Document <: GraphQLNode
    definitions::Vector{Any}
end
Base.getindex(d::Document, ind...) = d.definitions[ind...]
Base.length(d::Document) = length(d.definitions)
Base.firstindex(d::Document) = 1
Base.lastindex(d::Document) = length(d.definitions)
Base.iterate(d::Document, state = 1) = state > length(d) ? nothing : (d[state], state + 1)

struct SelectionSet <: GraphQLNode
    selections::Vector{Any}
end

struct TypeDefinition <: GraphQLNode
    description::Any
    type::Any
end

struct Argument <: GraphQLNode
    name::Any
    value::Any
end

struct Arguments <: GraphQLNode
    arguments::Vector{Argument}
end

struct Directive <: GraphQLNode
    name::Any
    arguments::Vector{Argument}
end

struct Directives
    directives::Vector{Directive}
end

struct Alias <: GraphQLNode
    name::Any
end

struct Field <: GraphQLNode
    alias::Union{Alias,Nothing}
    name::Any
    arguments::Vector{Argument}
    directives::Vector{Directive}
    selection_set::Vector{Any}
end

struct TypeCondition <: GraphQLNode
    named_type::NamedType
end

struct ListValue <: GraphQLNode
    values::Vector{Any}

    function ListValue(values)
        new(_force_any(values))
    end
end

struct ObjectField <: GraphQLNode
    name::Any
    value::Any
end

struct ObjectValue <: GraphQLNode
    object_fields::Vector{ObjectField}
end

struct Variable <: GraphQLNode
    name::Any
end

struct DefaultValue <: GraphQLNode
    value::Any
end

struct VariableDefinition <: GraphQLNode
    var::Variable
    type::Any
    value::DefaultValue
end

struct VariableDefinitions <: GraphQLNode
    variable_definitions::Vector{VariableDefinition}
end

struct OperationTypeDefinition <: GraphQLNode
    operation_type::String
    named_type::NamedType
end

struct SchemaDefinition <: GraphQLNode
    directives::Vector{Directive}
    operation_type_definitions::Vector{OperationTypeDefinition}
end

struct ScalarTypeDefinition <: GraphQLNode
    name::Any
    directives::Vector{Directive}
end

struct ScalarTypeExtension <: GraphQLNode
    name::Any
    directives::Vector{Directive}
end

struct InputValueDefinition <: GraphQLNode
    description::Any
    name::Any
    type::Any
    default_value::Union{DefaultValue,Nothing}
    directives::Vector{Directive}
end

struct ArgumentsDefinition <: GraphQLNode
    input_value_definitions::Vector{InputValueDefinition}
end

struct FieldDefinition <: GraphQLNode
    description::Any
    name::Any
    arguments_definition::Union{ArgumentsDefinition,Nothing}
    type::Any
    directives::Vector{Directive}
end

struct FieldsDefinition <: GraphQLNode
    field_definitions::Vector{FieldDefinition}
end

struct EnumValueDefinition <: GraphQLNode
    description::Any
    value::Any
    directives::Vector{Directive}
end

struct EnumValuesDefinition <: GraphQLNode
    values::Vector{EnumValueDefinition}
end

struct EnumTypeDefinition <: GraphQLNode
    name::Any
    directives::Vector{Directive}
    enums::Vector{EnumValueDefinition}
end

struct InputFieldsDefinition <: GraphQLNode
    input_value_definitions::Vector{InputValueDefinition}
end

struct InputObjectTypeDefinition <: GraphQLNode
    name::Any
    directives::Vector{Directive}
    fields::Vector{InputValueDefinition}
end

struct DirectiveDefinition <: GraphQLNode
    description::Any
    name::Any
    arguments_definition::Union{ArgumentsDefinition,Nothing}
    directive_locations::Any
end

struct UnionTypeDefinition <: GraphQLNode
    name::Any
    directives::Vector{Directive}
    types::Vector{NamedType}
end

struct InterfaceTypeDefinition <: GraphQLNode
    name::Any
    directives::Vector{Directive}
    fields_definition::Vector{FieldDefinition}
end

struct ObjectTypeDefinition <: GraphQLNode
    name::Any
    implements_interfaces::Any
    directives::Vector{Directive}
    fields_definition::Vector{FieldDefinition}
end

struct FragmentSpread <: GraphQLNode
    fragment_name::Any
    directives::Vector{Directive}
end

struct InlineFragment <: GraphQLNode
    type_condition::TypeCondition
    directives::Vector{Directive}
    selection_set::Vector{Any}
end

struct FragmentDefinition <: GraphQLNode
    fragment_name::Any
    type_condition::TypeCondition
    directives::Vector{Directive}
    selection_set::Vector{Any}
end

# Copied from OpenQASM.jl
#
# NOTE:
# In order to preserve some line number
# we usually don't annote types to AST

# work around JuliaLang/julia/issues/38091
function _force_any(x)
    if isnothing(x)
        return Any[]
    else
        return Vector{Any}(x)
    end
end

function Base.convert(
    ::Type{Vector{Directive}},
    x::Some{T},
) where {T<:Union{Directives,Nothing}}
    s = something(x)
    if isnothing(s)
        return Directive[]
    else
        Vector{Directive}(s.directives)
    end
end

function Base.convert(
    ::Type{Vector{Argument}},
    x::Some{T},
) where {T<:Union{Arguments,Nothing}}
    s = something(x)
    if isnothing(s)
        return Argument[]
    else
        Vector{Argument}(s.arguments)
    end
end

function Base.convert(
    ::Type{Vector{Any}},
    x::Some{T},
) where {T<:Union{SelectionSet,Nothing}}
    s = something(x)
    if isnothing(s)
        return Any[]
    else
        Vector{Any}(s.selections)
    end
end

function Base.convert(
    ::Type{Vector{FieldDefinition}},
    x::Some{T},
) where {T<:Union{FieldsDefinition,Nothing}}
    s = something(x)
    if isnothing(s)
        return FieldDefinition[]
    else
        Vector{FieldDefinition}(s.field_definitions)
    end
end

function Base.convert(
    ::Type{Vector{InputValueDefinition}},
    x::Some{T},
) where {T<:Union{ArgumentsDefinition,InputFieldsDefinition,Nothing}}
    s = something(x)
    if isnothing(s)
        return InputValueDefinition[]
    else
        Vector{InputValueDefinition}(s.input_value_definitions)
    end
end

function Base.convert(
    ::Type{Vector{EnumValueDefinition}},
    x::Some{T},
) where {T<:Union{EnumValuesDefinition,Nothing}}
    s = something(x)
    if isnothing(s)
        return EnumValueDefinition[]
    else
        Vector{EnumValueDefinition}(s.values)
    end
end

function Base.convert(
    ::Type{Vector{NamedType}},
    x::Some{T},
) where {T<:Union{Vector{NamedType},Nothing}}
    s = something(x)
    return isnothing(s) ? NamedType[] : s
end

function Base.convert(::Type{Union{T,Nothing}}, x::Some{Nothing}) where {T<:GraphQLNode}
    nothing
end

function Base.convert(::Type{Union{T,Nothing}}, x::Some{T}) where {T<:GraphQLNode}
    something(x)
end

Base.convert(::Type{String}, ::Some{Nothing}) = ""


text:
function othermaxplus{T}(dim::Int64,li::Vector{Int64},lj::Vector{Int64},
                        lw::Vector{T},m::Int64,n::Int64)
  # OTHERMAXPLUS Apply the other-max-plus operator to a sparse matrix
  #
  # The other-max-plus operator applies the max-plus aggegration function
  # over the rows or columns of the matrix, where, for
  # each non-zeros, the non-zero itself cannot be the maximum.  Consequently,
  # Hence, it's the max-plus of all the "other" elements in the row or column
  # (which is controlled by dim).
  #
  # omp = othermaxplus(dim,li,lj,lw,m,n) applies the other-max-plus operator
  # to the matrix sparse(li, lj, lw, m, n).  The return value omp gives
  # the other-max-plus value for each non-zero element, e.g. the matrix
  # is sparse(li, lj, omp, m, n).

  if dim == 1
    # max-plus over cols
    i1 = lj
    i2 = li
    N = n
  else
    # max-plus over rows
    i1 = li
    i2 = lj
    N = m
  end
  # the output of the other-max-plus is either the maximum element
  # in the row or column (if the element itself isn't the maximum) or
  # the second largest element in the row or column (if the element itself
  # IS the maximum).

  dimmax1 = zeros(Float64,N)      # largest value
  dimmax2 = zeros(Float64,N)     # second largest value
  #this is correct because of the definition of the max-plus function.
  dimmaxind = zeros(Float64,N)   # index of largest value
  nedges = length(li)
  for i = 1:nedges
    if lw[i] > dimmax2[i1[i]]
      if lw[i] > dimmax1[i1[i]]
        dimmax2[i1[i]] = dimmax1[i1[i]]
        dimmax1[i1[i]] = lw[i]
        dimmaxind[i1[i]] = i2[i]
      else
        dimmax2[i1[i]] = lw[i]
      end
    end
  end
  omp = zeros(Float64,size(lw))
  for i = 1:nedges
    if i2[i] == dimmaxind[i1[i]]
      omp[i] = dimmax2[i1[i]]
    else
      omp[i] = dimmax1[i1[i]]
    end
  end
  return omp
end

max_stars_repo_name: nassarhuda/NetworkAlign

max_issues_repo_name: nassarhuda/NetworkAlign

max_forks_repo_name: nassarhuda/NetworkAlign


text:
@testset "5.3 Inverse tangent" begin
    include("5.3.2 (d x)^m (a+b arctan(c x^n))^p.jl")
    include("5.3.3 (d+e x)^m (a+b arctan(c x^n))^p.jl")
    include("5.3.4 u (a+b arctan(c x))^p.jl")
    include("5.3.5 u (a+b arctan(c+d x))^p.jl")
    include("5.3.6 Exponentials of inverse tangent.jl")
    include("5.3.7 Inverse tangent functions.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.3 Inverse tangent/.jl


text:
@testset "006" begin
    @test tdpe.e006(10) == 2640
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
@debug "Loading DataFrames support into Gadfly"

using .DataFrames

function meltdata(U::AbstractDataFrame, colgroups::Vector{Col.GroupedColumn})
    um, un = size(U)

    # Figure out the size of the new melted matrix
    allcolumns = Set{Symbol}(propertynames(U))

    vm = um
    colidxs = [colgroup.columns===nothing ? collect(allcolumns) : colgroup.columns for colgroup in colgroups]
    vm *=  prod(length.(colidxs))
    grouped_columns = reduce(vcat, colidxs)

    ungrouped_columns = setdiff(allcolumns, grouped_columns)
    vn = length(colgroups) + length(ungrouped_columns)

    V = AbstractArray[]
    vnames = Symbol[]
    colmap = Dict{Any, Int}()

    eltypd = Dict{Symbol, DataType}(k=>v for (k,v) in zip(propertynames(U), eltype.(eachcol(U))))
    # allocate vectors for grouped columns
    for (j, (colgroup, colidx)) in enumerate(zip(colgroups, colidxs))
        eltyp = promote_type(getindex.([eltypd], colidx)...)
        push!(V, eltyp == Vector ? Array{eltyp}(undef, vm) : Array{Union{Nothing,eltyp}}(undef, vm))
        name = gensym()
        push!(vnames, name)
        colmap[colgroup] = j
    end

    # allocate vectors for ungrouped columns
    for (j, col) in enumerate(ungrouped_columns)
        push!(V, Array{eltypd[col]}(undef, vm))
        colmap[col] = j + length(colgroups)
        push!(vnames, col)
    end

    # Indicator columns for each colgroup
    col_indicators = Array{Symbol}(undef, vm, length(colgroups))
    row_indicators = Array{Int}(undef, vm, length(colgroups))

    vi = 1
    for ui in 1:um
        for colidx in product(colidxs...)
            # copy grouped columns
            for (vj, uj) in enumerate(colidx)
                V[vj][vi] = U[ui, uj]
                col_indicators[vi, vj] = uj
                row_indicators[vi, vj] = ui
            end

            # copy ungrouped columns
            for (vj, uj) in enumerate(ungrouped_columns)
                V[vj + length(colgroups)][vi] = U[ui, uj]
            end

            vi += 1
        end
    end

    df = DataFrame(; collect(zip(vnames, V))...)
    return MeltedData(U, df, row_indicators, col_indicators, colmap)
end

evalmapping(source::MeltedData{T}, arg::Col.GroupedColumnValue) where T<:AbstractDataFrame =
    source.melted_data[:,source.colmap[Col.GroupedColumn(arg.columns)]]

evalmapping(source::AbstractDataFrame, arg::Symbol) = source[:,arg]
evalmapping(source::AbstractDataFrame, arg::AbstractString) = evalmapping(source, Symbol(arg))
evalmapping(source::AbstractDataFrame, arg::Integer) = source[:,arg]
evalmapping(source::AbstractDataFrame, arg::Expr) = with(source, arg)


text:
using Random
import LinearAlgebra: norm, mul!, ldiv!

module mytestmodule

using FourierFlows, FFTW
import LinearAlgebra: norm, mul!, ldiv!

export updatevars!, set_u!

struct Params <: AbstractParams
    mu::Float64              # drag
end

function Equation(p, g)
  LC = 0.0*ones(size(g.kr))
  @. LC = -p.mu
  FourierFlows.Equation(LC, calcN!)
end

# Construct Vars type
physvars = [:u]
transvars = [:uh]
expr = FourierFlows.structvarsexpr(:Vars, physvars, transvars; vardims=1)
eval(expr)

"Returns the Vars object for Kuramoto-Sivashinsky."
function Vars(g)
  @createarrays Float64 (g.nx,) u
  @createarrays Complex{Float64} (g.nkr,) uh
  Vars(u, uh)
end

function calcN!(N, sol, t, s, v, p, g)
  nothing
end

function set_u!(prob::AbstractProblem, u)
  v, s, g = prob.vars, prob.state, prob.grid
  mul!(s.sol, g.rfftplan, u)
  updatevars!(v, s, g)
  nothing
end

function updatevars!(v, s, g)
  v.uh .= s.sol
  ldiv!(v.u, g.rfftplan, v.uh)
  nothing
end

updatevars!(prob::AbstractProblem) = updatevars!(prob.vars, prob.state, prob.grid)

@inline function pseudoenergy(prob)
  s = prob.state
  sum(abs.(s.sol).^2)
end

@inline function pseudoenstrophy(prob)
  s = prob.state
  0.5*sum(abs.(s.sol).^2)
end

end #module

function testsimplediagnostics()
    nx = 16
    Lx = 2π
    mu = 0.02
    dt = 0.01

    g = OneDGrid(nx, Lx)
    p = mytestmodule.Params(mu)
    v = mytestmodule.Vars(g)
    eq = mytestmodule.Equation(p, g)
    ts = FourierFlows.ETDRK4TimeStepper(dt, eq.LC)
    prob = FourierFlows.Problem(g, v, p, eq, ts)

    Random.seed!(1234)
     u0 = randn(size(g.x))

    mytestmodule.set_u!(prob, u0)

    nsteps = 200
    extrasteps = 20

    freqE = 2
    E = Diagnostic(mytestmodule.pseudoenergy, prob; nsteps=nsteps, freq=freqE)

    freqZ = 1
    Z = Diagnostic(mytestmodule.pseudoenstrophy, prob; nsteps=nsteps, freq=freqZ)

    diags = [E, Z]

    nstepstot = nsteps + extrasteps
    while prob.step < nstepstot
      stepforward!(prob, diags, 1)
    end

    stepforward!(prob, 1)
    increment!(Z)
    stepforward!(prob, 1)
    increment!([E,Z])

     resize!(Z, Z.count+2)
     Z.count += 1
     resize!(E, E.count+1)
     E.count += 1

    stepforward!(prob, 1)
    update!(Z)
    Z.count += 1
    stepforward!(prob, 1)
    update!([E, Z])
    E.count += 1
    Z.count += 1

    (
     isapprox(length(E.data), 2+Int(round(nstepstot/freqE+1))) &&
     isapprox(length(Z.data), 4+Int(round(nstepstot/freqZ+1))) &&
     isapprox(norm(E.data), norm(E.data[1]*exp.(-2*mu*E.time)), rtol=1e-13) &&
     isapprox(norm(Z.data), norm(Z.data[1]*exp.(-2*mu*Z.time)), rtol=1e-13)
    )
end

@test testsimplediagnostics()


text:
include("Kalman.jl")
# using Main.Kalman
include("SimWorld.jl")

struct Agent{KMs<:Array{KalmanModel,1}, FN<:Function, I<:Int}#, SWs<:Array{SimWorld,1}}
    models::KMs
    choice::FN
    dims::I

    # worlds::SWs
end

function getMSE(MusLL::Array{Array{Float64,2},1}, dims::Int, Vars::Array{Array{Float64,2},1})
    MVar = reshape(vcat([mean(reshape(diag(data), dims, :), dims=1) for data in Vars]...), :, size(Vars,1))
    MSEs = [mean(reshape(data[1:end-1,:], dims, :, size(data,2)), dims=1) for data in MusLL]
    MSE = dropdims(mean(vcat(MSEs...), dims=1), dims=1)

    return MSE, MVar
end
#
# function plotMSE(MSE::Array{Float64,2}, mVars::Array{Float64,2}, allT::Array{Float64,1}, opts::PlotOpts)
#
#     mVars = hcat([diag(M) for M in kworld.Vars]...)
#
#     subplot(2, 2, 1)
#     plot(allT, MSE[1,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
#     plot(allT, mVars[1,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
#     # plotshade(RMSE[1,:], eVars[1,:], allT, opts)
#     # fill_between(allT, RMSE[1,:]+mVars[1,:], RMSE[1,:]-mVars[1,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
#     xlabel("time")
#     ylabel("mean square error")
#     title("Position MSE")
#
#     if size(MSE,1) > 3
#         subplot(2, 2, 2)
#         plot(allT, MSE[2,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
#         plot(allT, mVars[2,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
#         # plotshade(RMSE[2,:], eVars[2,:], allT, opts)
#         # fill_between(allT, RMSE[2,:]+mVars[2,:], RMSE[2,:]-mVars[2,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
#         xlabel("time")
#         ylabel("mean square error")
#         title("Object Distance MSE")
#     end
#
#     subplot(2, 2, 3)
#     plot(allT, MSE[end-1,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
#     plot(allT, mVars[end-1,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
#     # plotshade(RMSE[end,:], eVars[end,:], allT, opts)
#     # fill_between(allT, RMSE[end,:]+mVars[end,:], RMSE[end,:]-mVars[end,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
#     xlabel("time")
#     ylabel("mean square error")
#     title("Velocity MSE")
#
#     subplot(2, 2, 4)
#     plot(allT, MSE[end,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
#     # plotshade(RMSE[end,:], eVars[end,:], allT, opts)
#     # fill_between(allT, RMSE[end,:]+mVars[end,:], RMSE[end,:]-mVars[end,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
#     xlabel("time")
#     ylabel("Log Likelihood")
#     title("Overall Model Performance")
#
#     return
# end

max_stars_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_issues_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_forks_repo_name: rhoadesScholar/Drugowitsch-Rotation


text:
module MyPackage1
using ForwardDiff
# Write your package code here.
greet = println("hello world")
include("extra_file.jl")

export my_f,my_probability,my_sigmoid
println("hello world")
println("My name is sandesh")
println("My father name is ramesh")
println("My mother name is padma")
println("My brother name is rishikesh")

end


text:
###########################################################################
#  Hamiltonian Monte Carlo (HMC)
#
#  Parameters :
#    - nLeaps : number of intermediate jumps wihtin each step
#    - leapStep : inner steps scaling
#    - tuner: for tuning the HMC parameters
#
###########################################################################

export HMC

println("Loading HMC(nLeaps, leapStep, tuner) sampler")

###########################################################################
#                  HMC specific 'tuners'
###########################################################################
abstract HMCTuner <: MCMCTuner

type EmpiricalHMCTune
  nLeaps::Int
  leapStep::Float64
  accepted::Int
  proposed::Int
  rate::Float64

  function EmpiricalHMCTune(nLeaps::Int, leapStep::Float64, accepted::Int, proposed::Int, rate::Float64)
    @assert nLeaps > 0 "Number of leapfrog steps ($nLeaps) should be > 0"
    @assert leapStep > 0 "leapfrog step size ($leapStep) should be > 0"
    @assert 0 <= accepted "Number of accepted Monte Carlo steps ($accepted) should be non negative"
    @assert 0 <= proposed "Number of proposed Monte Carlo steps ($proposed) should be non negative"
    new(nLeaps, leapStep, accepted, proposed)
  end
end

EmpiricalHMCTune(nLeaps::Int, leapStep::Float64, accepted::Int, proposed::Int) =
  EmpiricalHMCTune(nLeaps::Int, leapStep::Float64, accepted::Int, proposed::Int, NaN)

function adapt!(tune::EmpiricalHMCTune, tuner::EmpMCTuner)
  tune.rate = tune.accepted/tune.proposed
  tune.leapStep *= (1/(1+exp(-11*(tune.rate-tuner.targetRate)))+0.5)
  tune.nLeaps = min(tuner.maxStep, ceil(tuner.targetPath/tune.leapStep))
end

count!(tune::EmpiricalHMCTune) = (tune.accepted += 1)

reset!(tune::EmpiricalHMCTune) = ((tune.accepted, tune.proposed) = (0, 0))

###########################################################################
#                  HMC type
###########################################################################

immutable HMC <: MCMCSampler
  nLeaps::Int
  leapStep::Float64
  tuner::Union(Nothing, MCMCTuner)

  function HMC(nLeaps::Int, leapStep::Real, tuner::Union(Nothing, MCMCTuner))
    @assert nLeaps>0 "inner steps should be > 0"
    @assert leapStep>0 "inner steps scaling should be > 0"
    new(nLeaps, leapStep, tuner)
  end
end
HMC(tuner::Union(Nothing, MCMCTuner)=nothing) = HMC(10, 0.1, tuner)
HMC(nLeaps::Int, tuner::Union(Nothing, MCMCTuner)=nothing) = HMC(nLeaps, 0.1, tuner)
HMC(nLeaps::Int, leapStep::Float64) = HMC(nLeaps, leapStep, nothing)
HMC(leapStep::Float64, tuner::Union(Nothing, MCMCTuner)=nothing) = HMC(10, leapStep, tuner)
HMC(;init::Int=10, scale::Float64=0.1, tuner::Union(Nothing, MCMCTuner)=nothing) =
  HMC(init, scale, tuner)

###########################################################################
#                  HMC task
###########################################################################

####  Helper functions and types for HMC sampling task
type HMCSample
  pars::Vector{Float64} # sample position
  grad::Vector{Float64} # gradient
  m::Vector{Float64}    # momentum
  logTarget::Float64    # log likelihood
  H::Float64            # Hamiltonian
end
HMCSample(pars::Vector{Float64}) = HMCSample(pars, Float64[], Float64[], NaN, NaN)

calc!(s::HMCSample, ll::Function) = ((s.logTarget, s.grad) = ll(s.pars))
update!(s::HMCSample) = (s.H = -s.logTarget+0.5*dot(s.m, s.m))

function leapfrog(s::HMCSample, ve, ll::Function)
  n = deepcopy(s)  # make a full copy
  n.m += 0.5*n.grad*ve
  n.pars += ve * n.m
  calc!(n, ll)
  n.m += 0.5*n.grad*ve
  update!(n)

  n
end

####  HMC task
function SamplerTask(model::MCMCModel, sampler::HMC, runner::MCMCRunner)
  local state0
  local nLeaps, leapStep

  @assert hasgradient(model) "HMC sampler requires model with gradient function"

  # hook inside Task to allow remote resetting
  task_local_storage(:reset,
             (resetPars::Vector{Float64}) -> (state0 = HMCSample(copy(resetPars));
                                              calc!(state0, model.evalallg)) )

  # initialization
  state0 = HMCSample(copy(model.init))
  calc!(state0, model.evalallg)
  @assert isfinite(state0.logTarget) "Initial values out of model support, try other values"

  if isa(sampler.tuner, EmpMCTuner); tune = EmpiricalHMCTune(sampler.nLeaps, sampler.leapStep, 0, 0); end

  #  main loop
  i = 1
  while true
    local j, state

    if isa(sampler.tuner, EmpMCTuner)
      tune.proposed += 1
      nLeaps, leapStep = tune.nLeaps, tune.leapStep
    else
      nLeaps, leapStep = sampler.nLeaps, sampler.leapStep
    end

    state0.m = randn(model.size)
    update!(state0)
    state = deepcopy(state0)

    for j = 1:nLeaps
      state = leapfrog(state, leapStep, model.evalallg)
    end

    # accept if new is good enough
    if rand() < exp(state0.H-state.H)
      ms = MCMCSample(state.pars, state.logTarget, state.grad, state0.pars, state0.logTarget, state0.grad,
        {"accept" => true})
      produce(ms)
      state0 = deepcopy(state)
      if isa(sampler.tuner, EmpMCTuner); tune.accepted += 1; end
    else
      ms = MCMCSample(state0.pars, state0.logTarget, state0.grad, state0.pars, state0.logTarget, state0.grad,
        {"accept" => false})
      produce(ms)
    end

    if isa(sampler.tuner, EmpMCTuner) && i<= runner.burnin && mod(i, sampler.tuner.adaptStep) == 0
      adapt!(tune, sampler.tuner)
      reset!(tune)
      if sampler.tuner.verbose
        println("Burn-in teration $i of $(runner.burnin): ", round(100*tune.rate, 2), " % acceptance rate")
      end
    end

    i += 1
  end
end


text:
function write_output(	file_path::String,
					model::String,
					number_of_components::Int64,
					number_of_abc_samples::Int64,
					m::Array{Float64, 2},
					c::Array{Float64, 2},
					az::Array{Float64, 2},
					dist::Array{Float64, 1},
					w::Array{Float64, 1},
					weighting_scheme::String,
					gamma::Float64,
					t_exec::Float64,
					number_of_iterations::Int64,
					number_of_simulations::Int64)

	file_stream::IOStream = open(file_path, "w")

	@printf(file_stream, "%s", "<output>\n")

	write_key(file_stream, "model", model)
	write_key(file_stream, "number_of_components", number_of_components)
	write_key(file_stream, "number_of_abc_samples", number_of_abc_samples)
	write_key(file_stream, "m", m[:])
	write_key(file_stream, "c", c[:])
	write_key(file_stream, "az", az[:])
	write_key(file_stream, "dist", dist)
	write_key(file_stream, "w", w)
	write_key(file_stream, "weighting_scheme", weighting_scheme)
	write_key(file_stream, "gamma", gamma)
	write_key(file_stream, "t_exec", t_exec)
	write_key(file_stream, "number_of_iterations", number_of_iterations)
	write_key(file_stream, "number_of_simulations", number_of_simulations)

	@printf(file_stream, "%s", "</output>")

	close(file_stream)

	nothing
end


text:
using GR

r = 1
θ, φ = meshgrid(range(0, stop=2π, length=200), range(0, stop=2π, length=100))

x = r * sin.(θ) .* sin.(φ)
y = r * sin.(θ) .* cos.(φ)
z = r * cos.(θ)

surface(x, y, z)


text:
using JLD2, FileIO, DataFrames, Dates

println("Iniciando redução dos dados para matrizes de 3 colunas.")

root_dir = dirname(pwd())
unzipped_files_path = *(root_dir, "/Datasets/Dataset_Bovespa/Unzipped/")
matrix_jld2_file_name = (unzipped_files_path * "COTAHIST" * ".jld2")

# caminho para matriz reduzida
reduced_matrix_file_name = (unzipped_files_path * "COTAHIST_REDUCED" * ".jld2")

file_begin_counter = 1995
file_end_counter = 2019
file_counter = file_begin_counter

time_inicio = Dates.now()

jldopen(reduced_matrix_file_name, "w") do jldfile
    cotahist_reduced_prefix = "COTAHIST_REDUCED"
    cotahist_reduced = JLD2.Group(jldfile, cotahist_reduced_prefix)
    while (file_counter <= file_end_counter)
        local file_counter_string = string(file_counter)

        # criacao do dataframe
        local df = DataFrame(data = [], codneg = [], preult = [])

        # matriz com dados completos
        local matrix = load(matrix_jld2_file_name, "COTAHIST/" * file_counter_string)

        for acao in eachrow(matrix)
            #if (acao[4] == 10 || acao[4] == 20) #mercado lote padrao ou fracionario
            if (acao[4] == 10)
                push!(df, [acao[1], rstrip(acao[3]), acao[13]]) # acao[1] = data; acao[3] = nome acao; acao[13] = preult
            end
        end

        local matrix_reduced = convert(Matrix{Any}, df)
        cotahist_reduced[file_counter_string] = matrix_reduced
        global file_counter += 1
    end
end

time_fim = Dates.now()
println("TI: " * Dates.format(time_inicio, "S:s"))
println("TF: " * Dates.format(time_fim, "S:s"))
println(time_fim - time_inicio)
println("Finalizada execução de redução dos dados para matrizes de 3 colunas.")