text:
export FWI, FWIExample, compute_observation, plot, compute_misfit

@with_kw mutable struct FWI
    nz::Int64 = 134
    nx::Int64 = 384
    dz::Float64 = 24.
    dx::Float64 = 24.
    dt::Float64 = 0.0025
    nSteps::Int64 = 2000
    f0::Float64 = 4.5
    nPml::Int64 = 32
    nPad::Int64 = 32 - mod((nz+2*nPml), 32)
    nz_pad::Int64 = nz + 2*nPml + nPad
    nx_pad::Int64 = nx + 2*nPml
    para_fname::String = "para_file.json"
    survey_fname::String = "survey_file.json"
    data_dir_name::String = "Data"
    WORKSPACE::String = mktempdir()
    mask::Union{Missing, Array{Float64, 2}, PyObject} = missing
    mask_neg::Union{Missing, Array{Float64, 2}, PyObject} = missing
    ind_src_x::Union{Missing, Array{Int64, 1}} = missing
    ind_src_z::Union{Missing, Array{Int64, 1}} = missing
    ind_rec_x::Union{Missing, Array{Int64, 1}} = missing
    ind_rec_z::Union{Missing, Array{Int64, 1}} = missing
end

"""
    FWI(nz::Int64, nx::Int64, dz::Float64, dx::Float64, nSteps::Int64, dt::Float64;
            ind_src_z::Array{Int64, 1}, ind_src_x::Array{Int64, 1}, ind_rec_z::Array{Int64, 1}, ind_rec_x::Array{Int64, 1},
            kwargs...)

Creates a `FWI` structure that holds geometry data and settings.
Here `nz = n + 1`, `nx = m + 1`, `dz` and `dx` are mesh sizes.
`nSteps` is the number of total iterations. `dt` is the time step.
`ind_src_z` and `ind_src_x` are arrays of source locations.
`ind_rec_z` and `ind_rec_x` are arrays of receiver locations.
"""
function FWI(nz::Int64, nx::Int64, dz::Float64, dx::Float64, nSteps::Int64, dt::Float64;
        ind_src_z::Array{Int64, 1}, ind_src_x::Array{Int64, 1}, ind_rec_z::Array{Int64, 1}, ind_rec_x::Array{Int64, 1},
        kwargs...)
    fwi = FWI(nx = nx, nz = nz, dx = dx, dz = dz, nSteps = nSteps, dt = dt;
            ind_src_x = ind_src_x,
            ind_src_z = ind_src_z,
            ind_rec_x = ind_rec_x,
            ind_rec_z = ind_rec_z,
            kwargs...)
    Mask = zeros(fwi.nz_pad, fwi.nx_pad)
    Mask[fwi.nPml+1:fwi.nPml+nz, fwi.nPml+1:fwi.nPml+nx] .= 1.0
    Mask[fwi.nPml+1:fwi.nPml+10,:] .= 0.0
    fwi.mask = constant(Mask)
    fwi.mask_neg = 1 - constant(Mask)
    @assert length(ind_rec_x)==length(ind_rec_z)
    @assert length(ind_src_x)==length(ind_src_z)
    paraGen(fwi.nz_pad, fwi.nx_pad, dz, dx, nSteps, dt, fwi.f0, fwi.nPml, fwi.nPad,
            joinpath(fwi.WORKSPACE,fwi.para_fname),
            joinpath(fwi.WORKSPACE,fwi.survey_fname),
            joinpath(fwi.WORKSPACE,fwi.data_dir_name))
    surveyGen(ind_src_z, ind_src_x, ind_rec_z, ind_rec_x, joinpath(fwi.WORKSPACE,fwi.survey_fname))
    return fwi
end

function PyPlot.:plot(fwi::FWI)
    close("all")
    x1, y1, x2, y2 = 0.0, 0.0, fwi.nx_pad*fwi.dx, fwi.nz_pad*fwi.dz
    plot(LinRange(x1, x2, 100), y1*ones(100), "k")
    plot(LinRange(x1, x2, 100), y2*ones(100), "k")
    plot(x1*ones(100), LinRange(y1, y2, 100), "k")
    plot(x2*ones(100), LinRange(y1, y2, 100), "k")

    x1, y1, x2, y2 = fwi.nPml * fwi.dx, fwi.nPml * fwi.dz, (fwi.nx_pad - fwi.nPml) * fwi.dx, (fwi.nz_pad - fwi.nPml) * fwi.dz
    plot(LinRange(x1, x2, 100), y1*ones(100), "g", label="PML Boundary")
    plot(LinRange(x1, x2, 100), y2*ones(100), "g")
    plot(x1*ones(100), LinRange(y1, y2, 100), "g")
    plot(x2*ones(100), LinRange(y1, y2, 100), "g")

    plot( (fwi.nPml .+ fwi.ind_rec_x .- 1) * fwi.dx, (fwi.nPml + fwi.nPad .+ fwi.ind_rec_z .- 1) * fwi.dz, "r^", label="Receiver", markersize=1)
    plot( (fwi.nPml .+ fwi.ind_src_x .- 1) * fwi.dx, (fwi.nPml + fwi.nPad .+ fwi.ind_src_z .- 1) * fwi.dz, "bv", label="Source", markersize=1)
    gca().invert_yaxis()
    xlabel("Distance")
    ylabel("Depth")
    legend()
    axis("equal")
end

function FWIExample()
    ind_src_x = collect(4:8:384)
    ind_src_z = 2ones(Int64, size(ind_src_x))
    ind_rec_x = collect(3:381)
    ind_rec_z = 2ones(Int64, size(ind_rec_x))
    FWI(134,384, 24., 24., 2000, 0.0025;
        ind_src_x = ind_src_x,
        ind_src_z = ind_src_z,
        ind_rec_x = ind_rec_x,
        ind_rec_z = ind_rec_z)
end

@doc raw"""
    compute_observation(sess::PyObject, fwi::FWI,
        cp::Union{Array{Float64}, PyObject},
        cs::Union{Array{Float64}, PyObject},
        ρ::Union{Array{Float64}, PyObject},
        stf_array::Union{Array{Float64}, PyObject},
        shot_ids::Union{Missing, Array{<:Integer}} = missing;
        gpu_id::Int64 = 0)

Computes the observations using given parameters. Note that `shot_ids` are 1-based. If it is missing, all shots are used.
"""
function compute_observation(sess::PyObject, fwi::FWI,
    cp::Union{Array{Float64}, PyObject},
    cs::Union{Array{Float64}, PyObject},
    ρ::Union{Array{Float64}, PyObject},
    stf_array::Union{Array{Float64}, PyObject},
    shot_ids::Union{Missing, Array{<:Integer}} = missing;
    gpu_id::Int64 = 0)
    cp_pad, cs_pad, ρ_pad = try_pad(fwi, cp, cs, ρ)
    stf_array = constant(stf_array)
    if length(size(stf_array))==1
        stf_array = repeat(stf_array', length(fwi.ind_src_z), 1)
    end
    λ_pad, μ_pad = velocity_to_moduli(cp_pad, cs_pad, ρ_pad)
    if ismissing(shot_ids)
        shot_ids = collect(1:length(fwi.ind_src_x))
    end
    shot_ids = shot_ids .- 1
    shot_ids_ = constant(shot_ids, dtype=Int32)
    data = fwi_obs_op(λ_pad, μ_pad, ρ_pad, stf_array, gpu_id, shot_ids_, joinpath(fwi.WORKSPACE, fwi.para_fname) )
    run(sess, data)
    data = zeros(length(shot_ids), fwi.nSteps, length(fwi.ind_rec_z))
    for i = 1:length(shot_ids)
        A = read("$(fwi.WORKSPACE)/Data/Shot$(shot_ids[i]).bin")
        data[i,:,:] = reshape(reinterpret(Float32,A),(fwi.nSteps ,length(fwi.ind_rec_z)))
    end
    data
end

@doc raw"""
    compute_misfit(fwi::FWI,
        cp::Union{Array{Float64}, PyObject},
        cs::Union{Array{Float64}, PyObject},
        ρ::Union{Array{Float64}, PyObject},
        stf_array::Union{Array{Float64}, PyObject},
        shot_ids::Union{Array{Int64}, PyObject} = missing;
        gpu_id::Int64 = 0, is_masked::Bool = false,
        cp_ref::Union{Array{Float64}, PyObject, Missing} = missing,
        cs_ref::Union{Array{Float64}, PyObject, Missing} = missing,
        ρ_ref::Union{Array{Float64}, PyObject, Missing} = missing)

Computes the misfit function for the simulation parameters $c_p$, $c_s$, $\rho$, and source time functions `stf_array`

- If `is_masked` is false, `compute_misfit` will add the mask `fwi.mask` to all variables.
- `gpu_id` is an integer in {0,1,2,...,#gpus-1}
- `shot_ids` is 1-based. If it is missing, all shots are used.
"""
function compute_misfit(fwi::FWI,
    cp::Union{Array{Float64}, PyObject},
    cs::Union{Array{Float64}, PyObject},
    ρ::Union{Array{Float64}, PyObject},
    stf_array::Union{Array{Float64}, PyObject},
    shot_ids::Union{Array{Int64}, PyObject} = missing;
    gpu_id::Int64 = 0, is_masked::Bool = false,
    cp_ref::Union{Array{Float64}, PyObject, Missing} = missing,
    cs_ref::Union{Array{Float64}, PyObject, Missing} = missing,
    ρ_ref::Union{Array{Float64}, PyObject, Missing} = missing)

    cp_pad, cs_pad, ρ_pad = try_pad(fwi, cp, cs, ρ)
    if !ismissing(cp_ref)
        cp_ref, cs_ref, ρ_ref = try_pad(fwi, cp_ref, cs_ref, ρ_ref)
    end

    cp_masked, cs_masked,ρ_masked = cp_pad, cs_pad, ρ_pad
    if !is_masked
        cp_masked = cp_pad .* fwi.mask + cp_ref .* fwi.mask_neg
        cs_masked = cs_pad .* fwi.mask + cs_ref .* fwi.mask_neg
        ρ_masked = ρ_pad .* fwi.mask + ρ_ref .* fwi.mask_neg
    end
    λ_masked, μ_masked = velocity_to_moduli(cp_masked, cs_masked,ρ_masked)

    stf_array = constant(stf_array)
    if length(size(stf_array))==1
        stf_array = repeat(stf_array', length(fwi.ind_src_z), 1)
    end
    if ismissing(shot_ids)
        shot_ids = collect(1:length(fwi.ind_src_x))
    end
    shot_ids = constant(shot_ids, dtype=Int32) - 1
    misfit = fwi_op(λ_masked, μ_masked, ρ_masked, stf_array, gpu_id, shot_ids, joinpath(fwi.WORKSPACE, fwi.para_fname))
end

#------------------------------------------------------------------------------------
function padding(fwi::FWI, cp::Union{PyObject, Array{Float64,2}})
    cp = constant(cp)
    nz, nx, nPml, nPad = fwi.nz, fwi.nx, fwi.nPml, fwi.nPad
    nz_orig, nx_orig = size(cp)
    tran_cp = tf.reshape(cp, (1, nz_orig, nx_orig, 1))
    if nz_orig!=nz || nx_orig!=nx
        @info "resizee image to required size"
        tran_cp = squeeze(tf.image.resize_bilinear(tran_cp, (nz, nx)))
    end
	cp_pad = tf.pad(cp, [nPml (nPml+nPad); nPml nPml], "SYMMETRIC")
	cp_pad = cast(cp_pad, Float64)
	return cp_pad
end

function padding(fwi::FWI, cp::Union{PyObject, Array{Float64,2}}, cq...)
    o = Array{PyObject}(undef, 1 + length(cq))
    o[1] = padding(fwi, cp)
    for i = 1:length(cq)
        o[i+1] = padding(fwi, cq[i])
    end
    o
end

function try_pad(fwi::FWI, cp::Union{PyObject, Array{Float64,2}})
    cp = convert_to_tensor(cp, dtype=Float64)
    if size(cp)!=(fwi.nz_pad, fwi.nx_pad)
        return padding(fwi, cp)
    else
        return cp
    end
end

function try_pad(fwi::FWI, cp::Union{PyObject, Array{Float64,2}},cq...)
    o = Array{PyObject}(undef, 1 + length(cq))
    o[1] = try_pad(fwi, cp)
    for i = 1:length(cq)
        o[i+1] = try_pad(fwi, cq[i])
    end
    o
end


text:
# This file is part of InvertedFiles.jl

using Intersections
import SimilaritySearch: search
export search, prepare_posting_lists_for_querying

"""
	prepare_posting_lists_for_querying(idx::AbstractInvertedFile, q, pools=getpools(idx), tol=1e-6)

Fetches and prepares the involved posting lists to solve `q`
"""
function prepare_posting_lists_for_querying(idx::AbstractInvertedFile, q, pools=getpools(idx), tol=1e-6)
	Q = getcachepostinglists(pools)

	@inbounds for (tokenID, val) in sparseiterator(q)
		val < tol && continue
		L = idx.lists[tokenID]
		if length(L) > 0
            push_posting_list!(Q, idx, tokenID, val)
		end
	end

	Q
end

"""
	search(idx::AbstractInvertedFile, q, res::KnnResult; pools=nothing)

Searches `q` in `idx` using the cosine dissimilarity, it computes the full operation on `idx`. `res` specify the query
"""
function search(idx::AbstractInvertedFile, q, res::KnnResult; pools=getpools(idx), tol=1e-6, t=1)
	Q = prepare_posting_lists_for_querying(idx, q, pools, tol)
    P = getcachepositions(length(Q), pools)

    cost = search(idx, Q, P, t) do objID, d
		push!(res, objID, d)
	end

    (res=res, cost=cost)
end

ext: jl


text:
mutable struct History
    logicsets::Array
    current::Integer
end

function Base.push!(x::History, logicset::LogicalCombo)
  x.logicsets = [x.logicsets[1:x.current]..., logicset]
  x.current += 1
  x
end

function Base.pop!(x::History)
  pop!(x.logicsets)
  x.current -= 1
  x
end

update!(x::History, logicset::LogicalCombo) = x.logicsets[x.current] = logicset

activelogicset(x::History) = x.logicsets[x.current]
priorlogicsets(x::History) = x.logicsets[1:x.current]
futurelogicsets(x::History) = x.logicsets[(x.current+1):end]

back!(x::History) = x.current = (x.current > 0 ? x.current - 1 : 0)
next!(x::History) = x.current = (x.current < length(x.logicsets) ? x.current + 1 : x.current)

activecommand(x::History)  = activelogicset(x).commands[end]
priorcommands(x::History)  = priorlogicsets(x) .|> z -> z.commands[end]

function futurecommands(x::History)
  y = futurelogicsets(x)
  (length(y) == 0) && return String[]
  y .|> z -> z.commands[end]
end

History() = History([logicalparse("#Session Started", verbose = false)], 1)

# active = History()
# priorcommands(active)
# futurecommands(active)
#
# push!(active, logicalparse("x in 1:4", activelogicset(active) ) )
# push!(active, logicalparse("y in 1:4", activelogicset(active) ) )
# push!(active, logicalparse("z in 1:4", activelogicset(active) ) )
#
# activelogicset(active)
# priorlogicsets(active)
#
# activecommand(active)
# priorcommands(active)
#
# back!(active)
# back!(active)
# futurecommands(active)
# priorcommands(active)
#
# next!(active)
# push!(active, logicalparse("q in 1:4", activelogicset(active) ) )
# priorcommands(active)


text:
rel_error(a, b) = (a - b) / b
abs_rel_error(a, b) = abs((a - b) / b)

function trimmed_mean(x)
    l, u = quantile(x, [0.2, 0.8])
    mean(x[(x .> l) .& (x .< u)])
end

rng = MersenneTwister(9237429)
parameters = Dict(
    :max_food => 400,
    :food_alpha => 2,
    :food_beta => 10,
    :growth_rate => 1.5,
    :death_rate => 1.05
)

######### testing beverton_holt(x, K, nu, mu)

# If available food is constant, and the population begins under the
# maximum, then it will approach the maximum in the limit but not exceed it.
step_cnt = 200
constant_food = beverton_holt_constant_food(2, step_cnt, parameters)
@test all(constant_food .< 400)
@test abs_rel_error(constant_food[step_cnt], parameters[:max_food]) < 1e-10

# Assert the value is decreasing.
greater_cnt = 0
for diff_idx in 1:(step_cnt - 1)
    if constant_food[diff_idx + 1] < diff_idx
        greater_cnt += 1
    end
end
@test greater_cnt == 0

# If the population begins high, and food is constant, it will fall.
high_pop = 2 * parameters[:max_food]
high_food = beverton_holt_constant_food(high_pop, step_cnt, parameters)
@test all(high_food .> 400)
# Falling slowly because rate is 1.05, so this value is 0.015.
@test abs_rel_error(high_food[step_cnt], parameters[:max_food]) < 0.015

# If the nu and mu are in the right order in the function, then the rise to
# max population will be faster than the fall from high population.
@test abs(constant_food[step_cnt] - 400) < abs(high_food[step_cnt] - 400)

########## testing forecast_population(young, step_cnt, params, rng)

# Increasing the food parameter increases the average value.
lower_food_parameters = copy(parameters)
lower_food_parameters[:max_food] = parameters[:max_food] / 2
high_food = forecast_population(50, 10000, parameters, rng)
low_food = forecast_population(50, 10000, lower_food_parameters, rng)
@test mean(high_food) > mean(low_food)

# Increasing the death rate should lower the average.
# Decreasing the death rate should raise it.
death_rate_trials = collect(0.5:0.1:1.5)
for death_rate in death_rate_trials
    higher_death_parameters = copy(parameters)
    higher_death_parameters[:death_rate] = 1 + (parameters[:death_rate] - 1) * death_rate
    # Copy the rng to get the same draws of random numbers.
    # It's a "variance reduction technique".
    rng2 = copy(rng)
    low_death = forecast_population(50, 10000, parameters, rng)
    high_death = forecast_population(50, 10000, higher_death_parameters, rng2)
    if death_rate < 0.99
        @test trimmed_mean(high_death) > trimmed_mean(low_death)
    elseif death_rate > 1.01
        @test trimmed_mean(high_death) < trimmed_mean(low_death)
    else
        @test abs(trimmed_mean(high_death) - trimmed_mean(low_death)) < 1
    end
end

# The means are near the same value where I last saw them.
long_run_cnt = 100
long_run_max = 10000
long_run_means = zeros(long_run_cnt)
long_run_trimmed = zeros(long_run_cnt)
for long_run_idx in 1:long_run_cnt
    long_run = forecast_population(2, long_run_max, parameters, rng)
    long_run_later = long_run[1000:long_run_max]
    long_run_means[long_run_idx] = mean(long_run)
    long_run_trimmed[long_run_idx] = trimmed_mean(long_run)
end

for lr_arr in [long_run_means, long_run_trimmed]
    @test abs(mean(lr_arr) - 72.5) < 0.5
    @test std(lr_arr) < 0.6
end

# Build an empirical distribution function.
single_sample = forecast_population(50, long_run_max, parameters, rng)
default_cdf = StatsBase.ecdf(single_sample[1000:long_run_max])

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
// declaration and initialization in same statement

int main() {
 int x = 7;
 printInt(x);
 return 0 ;

}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
"""
    tophat(img; dims=coords_spatial(img), r=1)
    tophat(img, se)

Performs morphological top-hat transform for given image, i.e., `img - opening(img, se)`.

$(_docstring_se)

This _white_ top-hat transform can be used to extract small white elements and details from
an image. To extract black details, the _black_ top-hat transform, also known as bottom-hat
transform, [`bothat`](@ref) can be used.

# Examples

```jldoctest; setup = :(using ImageMorphology)
julia> img = falses(5, 5); img[1, 1] = true; img[3:5, 3:5] .= true; img
5×5 BitMatrix:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  1  1
 0  0  1  1  1
 0  0  1  1  1

julia> Int.(tophat(img))
5×5 Matrix{$Int}:
 1  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

julia> Int.(tophat(img, strel_diamond(img))) # use diamond shape SE
5×5 Matrix{$Int}:
 1  0  0  0  0
 0  0  0  0  0
 0  0  1  0  0
 0  0  0  0  0
 0  0  0  0  0
```
"""
tophat(img; kwargs...) = tophat!(similar(img, maybe_floattype(eltype(img))), img, similar(img); kwargs...)
tophat(img, se) = tophat!(similar(img, maybe_floattype(eltype(img))), img, se, similar(img))

"""
    tophat!(out, img, buffer; [dims], [r])
    tophat!(out, img, se, buffer)

The in-place version of [`tophat`](@ref) with input image `img` and output image `out`. The
intermediate erosion result is stored in `buffer`.
"""
function tophat!(out, img, buffer; dims=coords_spatial(img), r=nothing)
    return tophat!(out, img, strel_box(img, dims; r), buffer)
end
function tophat!(out, img, se, buffer)
    opening!(out, img, se, buffer)
    @. out = img - out
    return out
end
tophat!(out::AbstractArray{<:Color3}, img, se, buffer) = throw(ArgumentError("color image is not supported"))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ExpressRouteCircuitRoutesTableApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the currently advertised routes table associated with the express route circuit in a resource group.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: devicePath::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitsRoutesTableListResult
"""
function _swaggerinternal_expressRouteCircuitsListRoutesTable(_api::ExpressRouteCircuitRoutesTableApi, resourceGroupName::String, circuitName::String, peeringName::String, devicePath::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", ExpressRouteCircuitsRoutesTableListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}/routeTables/{devicePath}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "devicePath", devicePath)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function expressRouteCircuitsListRoutesTable(_api::ExpressRouteCircuitRoutesTableApi, resourceGroupName::String, circuitName::String, peeringName::String, devicePath::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitsListRoutesTable(_api, resourceGroupName, circuitName, peeringName, devicePath, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function expressRouteCircuitsListRoutesTable(_api::ExpressRouteCircuitRoutesTableApi, response_stream::Channel, resourceGroupName::String, circuitName::String, peeringName::String, devicePath::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_expressRouteCircuitsListRoutesTable(_api, resourceGroupName, circuitName, peeringName, devicePath, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export expressRouteCircuitsListRoutesTable


text:
# Fix ambiguities on julia 0.4

*(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) * b)

*(a::fmpz, b::ResElem{fmpz}) = b*a

+(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) + b)

+(a::fmpz, b::ResElem{fmpz}) = b + a

-(a::ResElem{fmpz}, b::fmpz) = parent(a)(data(a) - b)

-(a::fmpz, b::ResElem{fmpz}) = parent(b)(a - data(b))

function ==(a::ResElem{fmpz}, b::fmpz)
   z = base_ring(a)(b)
   return data(a) == mod(z, modulus(a))
end

==(a::fmpz, b::ResElem{fmpz}) = b == a

#

*(::fmpz, ::PolyElem{fmpz}) = nothing

*(::PolyElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::PolyElem{fmpz}) = nothing

+(::PolyElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::PolyElem{fmpz}) = nothing

-(::PolyElem{fmpz}, ::fmpz) = nothing

==(::fmpz, ::PolyElem{fmpz}) = nothing

==(::PolyElem{fmpz}, ::fmpz) = nothing

divexact(::PolyElem{fmpz}, ::fmpz) = nothing

evaluate(::PolyElem{fmpz}, ::fmpz) = nothing

#

*(::fmpz, ::SeriesElem{fmpz}) = nothing

*(::SeriesElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::SeriesElem{fmpz}) = nothing

+(::SeriesElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::SeriesElem{fmpz}) = nothing

-(::SeriesElem{fmpz}, ::fmpz) = nothing

==(::fmpz, ::SeriesElem{fmpz}) = nothing

==(::SeriesElem{fmpz}, ::fmpz) = nothing

*(::fmpz, ::RelSeriesElem{fmpz}) = nothing

*(::RelSeriesElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::RelSeriesElem{fmpz}) = nothing

+(::RelSeriesElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::RelSeriesElem{fmpz}) = nothing

-(::RelSeriesElem{fmpz}, ::fmpz) = nothing

==(::fmpz, ::RelSeriesElem{fmpz}) = nothing

==(::RelSeriesElem{fmpz}, ::fmpz) = nothing

*(::fmpz, ::AbsSeriesElem{fmpz}) = nothing

*(::AbsSeriesElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::AbsSeriesElem{fmpz}) = nothing

+(::AbsSeriesElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::AbsSeriesElem{fmpz}) = nothing

-(::AbsSeriesElem{fmpz}, ::fmpz) = nothing

==(::fmpz, ::AbsSeriesElem{fmpz}) = nothing

==(::AbsSeriesElem{fmpz}, ::fmpz) = nothing

*(::fmpz, ::MatElem{fmpz}) = nothing

*(::MatElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::MatElem{fmpz}) = nothing

+(::MatElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::MatElem{fmpz}) = nothing

-(::MatElem{fmpz}, ::fmpz) = nothing

==(::MatElem{fmpz}, ::fmpz) = nothing

divexact(::MatElem{fmpz}, ::fmpz) = nothing

#

setindex_t!(a::nmod_mat, b::GenRes{fmpz}, i::Int, j::Int) = setindex_!(a, data(b), i, j)

*(::FracElem{fmpz}, ::fmpz) = nothing

*(::fmpz, ::FracElem{fmpz}) = nothing

+(::FracElem{fmpz}, ::fmpz) = nothing

+(::fmpz, ::FracElem{fmpz}) = nothing

-(::FracElem{fmpz}, ::fmpz) = nothing

-(::fmpz, ::FracElem{fmpz}) = nothing

==(::FracElem{fmpz}, ::fmpz) = nothing

==(::fmpz, ::FracElem{fmpz}) = nothing

divexact(::FracElem{fmpz}, ::fmpz) = nothing

divexact(::fmpz, ::FracElem{fmpz}) = nothing


text:
using DTU_BrayamValqui_SP2021
using Test
using Gurobi
using GLPK
using CPLEX

println("*****************OPF Solution*****************")
#OPF------------------------------------------------------#
@time (syscost_det,p_G_value,p_D_value,Dual_constraint6,f_value)=solver_OPF(Gurobi.Optimizer,"results/OPF_Solution.txt")

@time (syscost_det,Dual_constraint2,Dual_constraint1,λ_value,f_value)=solver_OPF_dual(Gurobi.Optimizer,"results/OPF_Dual_Solution.txt")

#cd("C:\\Users\\braya\\.julia\\dev\\DTU_BrayamValqui_SP2021\\test")


text:
using testCase
using Test

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


text:
using Dagitty
using Test

tests = [
    "DAG",
    "implied_conditional_independencies",
    "utils",
    "is_d_separated",
    "adj_sets",
]

for t ∈ tests
    @testset "$t" begin
        include("test_$(t).jl")
    end
end


text:
export ArduinoWithExternalTempUnit, getStatus, resetDAQ
# TODO comment relevant Arduino code once added to project
struct ArduinoWithExternalTempUnit <: SurveillanceUnit
  sd::Vector{SerialDevice}
  CommandStart::String
  CommandEnd::String
  delim::String
  numSensors::Number
  maxTemps::Array
  selectSensors::Array
  nameSensors::Array
end

function Base.close(su::ArduinoWithExternalTempUnit)
  for s in su.sd
    close(s.sp)
  end
end

#init from SurveillanceUnit.jl
function ArduinoWithExternalTempUnit(params::Dict)
  # Here we could put more parameters into the TOML file
  su = ArduinoWithExternalTempUnit(params["connection"],
                                   params["numSensors"],
                                   params["maxTemps"],
                                   params["selectSensors"],
                                   params["nameSensors"])
  setMaximumTemps(su, params["maxTemps"])
  return su
end

function ArduinoWithExternalTempUnit(portAdress::Vector{T}, numSensors::Number, maxTemps::Array, selectSensors::Array, nameSensors::Array) where T <: AbstractString
    # general parameters
    pause_ms::Int=30
    timeout_ms = [500, 1000]
    delim::String="#"
    delim_read::String="#"
    delim_write::String="#"
    baudrate = [9600, 115200]
    CommandStart="!"
    CommandEnd="*"
    ndatabits::Integer=8
    parity::SPParity=SP_PARITY_NONE
    nstopbits::Integer=1

    ### build SU arduino ###
    spSU = SerialPort(portAdress[1])
    open(spSU)
	set_speed(spSU, baudrate[1])
	set_frame(spSU,ndatabits=ndatabits,parity=parity,nstopbits=nstopbits)
	#set_flow_control(spSU,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
    sleep(2)
    flush(spSU)
    write(spSU, "!ConnectionEstablished*#")
    response=readuntil(spSU, Vector{Char}(delim_read), timeout_ms[1]);
    @info response
    if(!(response == "ArduinoSurveillanceV1" || response == "ArduinoSurveillanceV2") )
        close(spSU)
        @error "Connected to wrong Device: SU" response portAdress[1]
    else
        @info "Connection to ArduinoSU established."
    end

    ### build Temp arduino ###
    spTU = SerialPort(portAdress[2])
    open(spTU)
	set_speed(spTU, baudrate[2])
	set_frame(spTU,ndatabits=ndatabits,parity=parity,nstopbits=nstopbits)
	#set_flow_control(spTU,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
    sleep(2)
    flush(spTU)
    write(spTU, "!ConnectionEstablished*#")
    response=readuntil(spTU, Vector{Char}(delim_read), timeout_ms[2]);
    @info response
    if(!(response == "ArduinoTemperatureUnitV2") )
        close(spTU)
        @error "Connected to wrong Device: TU" response portAdress[2]
    else
        @info "Connection to ArduinoTU established."
    end

    sds = [SerialDevice(spSU, pause_ms, timeout_ms[1], delim_read, delim_write),
           SerialDevice(spTU, pause_ms, timeout_ms[2], delim_read, delim_write)]

    sp = ArduinoWithExternalTempUnit(sds, CommandStart, CommandEnd, delim, numSensors, maxTemps, selectSensors, nameSensors)

    return sp;
end

### general command structure
function ArduinoCommand(Arduino::ArduinoWithExternalTempUnit, cmd::String, id::Int)
    cmd=Arduino.CommandStart*cmd*Arduino.CommandEnd*Arduino.delim;
    return query(Arduino.sd[id],cmd);
end

###############################################
#----------- Temperature Unit ---------- "id=2"

function getTemperatures(Arduino::ArduinoWithExternalTempUnit; names::Bool=false)
    TempDelim = ","

    Temps = ArduinoCommand(Arduino, "GET:ALLTEMPS", 2)
    Temps = Temps[7:end]  #filter out "TEMPS:" at beginning of answer

    #@info Temps now()
    temp =  tryparse.(Float64,split(Temps,TempDelim))

    if length(temp) == Arduino.numSensors
        tempFloat = zeros(length(temp))
        for i=1:min(length(tempFloat),length(temp))
            if temp[i] != nothing
                tempFloat[i] = temp[i]
            end
        end

        if names
            return [tempFloat[Arduino.selectSensors] Arduino.nameSensors[Arduino.selectSensors]]
        else
            return tempFloat[Arduino.selectSensors]
        end
    else
        return zeros(length(Arduino.selectSensors))
    end
end

export showCommands
function showCommands(Arduino::ArduinoWithExternalTempUnit)
    print(ArduinoCommand(Arduino, "GET:COMMANDS", 2))
end

export setMaximumTemps
function setMaximumTemps(Arduino::ArduinoWithExternalTempUnit, maxTemps::Array)
    if length(maxTemps) == Arduino.numSensors
        maxTempString=""
        for i in 1:Arduino.numSensors
            if i != Arduino.numSensors
                maxTempString *= string(maxTemps[i])*","
            else
                maxTempString *= string(maxTemps[i])
            end
        end
        ack = ArduinoCommand(Arduino, "SET:MAXTEMPS:<"*maxTempString*">", 2)
        @info "acknowledge of MaxTemps from TempUnit."
    else
        @warn "Please parse a maximum temperature for each sensor" Arduino.numSensors
    end
end

export getMaximumTemps
function getMaximumTemps(Arduino::ArduinoWithExternalTempUnit)
    println(ArduinoCommand(Arduino, "GET:MAXTEMPS", 2))
end

export CheckACQ
function CheckACQ(Arduino::ArduinoWithExternalTempUnit,ACQ)
    if ACQ=="ACQ"
        @info "Command Received"
        return ACQ;
    else
        @warn "Error, Unknown response" ACQ
    end
end

function GetStatus(Arduino::ArduinoWithExternalTempUnit)
    status = ArduinoCommand(Arduino,"GET:STATS", 1)
    return status
end

function resetDAQ(Arduino::ArduinoWithExternalTempUnit)
    ACQ = ArduinoCommand(Arduino,"RESET:RP", 1)
    CheckACQ(Arduino,ACQ)
end

hasResetDAQ(su::ArduinoWithExternalTempUnit) = true


text:
# ------------------------------------------------------------------------------------------
# # Julia es rápido
#
# Frecuentemente, se usan benchmarks para comparar lenguajes. Éstos benchmarks pueden llevar
# a largas discusiones; primero, para saber lo que se está midiendo y segundo para explicar
# sus diferencias. Estas preguntas sencillas a veces son mucho más complicadas de lo que uno
# se imagina.
#
# El propósito de este notebook es para que tú puedas hacer un benchmark simple. Uno puede
# leer el notebook y ver que sucedió en la Macbook Pro 4-core Intel Core i7 del autor, o
# correrlo uno mismo.
#
# (Este material empezó como parte de una clase que dió Steven Johnson en MIT:
# https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-
# registers.ipynb.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Esquema de este notebook
# - Definir la función suma
# - Implementación y benchmark de la función suma en ...
#     - C
#     - python (interno)
#     - python (numpy)
#     - python (hecho a mano)
#     - Julia (interno)
#     - Julia (hecho a mano)
# - Resúmenes de resultados
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # `sum`: una función fácil de entender
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Considera la función  **sum** `sum(a)`, la cual calcula
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# con $n$ longitud de `a`.
# ------------------------------------------------------------------------------------------

a = rand(10^7) # Vector 1D uniforme en [0,1)

sum(a)

# ------------------------------------------------------------------------------------------
# El resultado esperado es 0.5 * 10^7, pues el promedio de cada entrada es .5
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Benchmarks de algunas manera en algunos lenguajes
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Julia cuenta con el paquete `BenchmarkTools.jl` para hacer benchmarks fácil y rápidos
# ------------------------------------------------------------------------------------------

#Pkg.add("BenchmarkTools")

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# #  1. El lenguaje C
#
# C es considerado el estándar dorado: difícil para el usuario, fácil para la máquina. Estar
# dentro de un factor de 2 de C puede ser muy satisfactorio. Sin emargo, aún dentro de C,
# existen muchos tipos de optimizaciones posibles que un usuario de C promedio puedo o no
# aprovechar.
#
# Su autor no habla C, entonces no va a leer la celda siguiente, pero será feliz en saber
# que puedes poner código de C en una sesión de Julia, compilarlo, y correrlo. Nota que
# `"""` denota una cadena de varias líneas.
# ------------------------------------------------------------------------------------------

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # Haz un directorio temporario

# compila a una biblioteca compartida pipeando C_code a gcc
# (funciona sólo con gcc instalado):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# define una funcion de Julia que llama a la función de C:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a) # teclea \approx y luego <TAB> para obtener el símbolo ≈

≈  # alias para la función `isapprox`

?isapprox

# ------------------------------------------------------------------------------------------
# ¡Ahora podemos correr el benchmark directo desde Julia!
# ------------------------------------------------------------------------------------------

c_bench = @benchmark c_sum($a)

println("C: Tiempo más rápido fue $(minimum(c_bench.times) / 1e6) msec")

d = Dict()  # un diccionario a.k.a un arreglo asociativo
d["C"] = minimum(c_bench.times) / 1e6  # en milisegundos
d

using Plots
gr()

t = c_bench.times / 1e6 # tiempos en milisegundos
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

# ------------------------------------------------------------------------------------------
# # 2. Python y `sum` interno
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# El paquete `PyCall` provee un interfaz de Julia a Python
# ------------------------------------------------------------------------------------------

#Pkg.add("PyCall")

using PyCall

# Llama una función de bajo nivel de PyCall para obtener una lista de Python
# porque por default PyCall convertirá a un arreglo de NumPy en vez (benchmarkeamos NumPy más abajo)

apy_list = PyCall.array2py(a, 1, 1)

# obtener el "sum" interno de Python
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($apy_list)

d["Python interno"] = minimum(py_list_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 3. Python: `numpy`
#
# ## Aprovechar arquitectura "SIMD" pero sólo cuando funciona
#
# `numpy` es una biblioteca de C optimizada, llamable desde Python, que se puede instalar en
# Julia haciendo:
# ------------------------------------------------------------------------------------------

using Conda
#Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]
apy_numpy = PyObject(a) # convierte a un arreglo de NumPy por default

py_numpy_bench = @benchmark $numpy_sum($apy_numpy)

numpy_sum(apy_list) # python thing

numpy_sum(apy_list) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 4. Python, hecho a mano
# ------------------------------------------------------------------------------------------

py"""
def py_sum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($apy_list)

sum_py(apy_list)

sum_py(apy_list) ≈ sum(a)

d["Python hecho a mano"] = minimum(py_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 5. Julia (interno)
#
#    ## Escrito directo en Julia, ¡no en C!
# ------------------------------------------------------------------------------------------

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia interno"] = minimum(j_bench.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # 6. Julia (hecho a mano)
# ------------------------------------------------------------------------------------------

function mysum(A)
    s = 0.0  # s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hecho a mano"] = minimum(j_bench_hand.times) / 1e6
d

# ------------------------------------------------------------------------------------------
# # Summary
# ------------------------------------------------------------------------------------------

for (key, value) in sort(collect(d))
    println(rpad(key, 20, "."), lpad(round(value, 1), 8, "."))
end

for (key, value) in sort(collect(d), by=x->x[2])
    println(rpad(key, 20, "."), lpad(round(value, 2), 10, "."))
end

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia-ES/12. Julia es rapido.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using BoardGameGeek

id = 75674 # https://boardgamegeek.com/boardgame/75674

game = gameinfo(id)
@test game.name == "Julian: Triumph Before the Storm"
@test game.usersrated > 65

reviews = gamereviews(id)
@test length(reviews) == game.usersrated
@test sum(r -> r.rating, reviews) / length(reviews) ≈ game.average

reviews2 = gamereviews(id; waittime=2.5, pagesize=50)
@test reviews2 == reviews

games = gameinfo([id])
@test first(games).name == game.name


text:
function days(times :: Array{Int64,1}, m :: MITgcmDatas
              ;min = 15, max = 20 , dt = m.dt )
  return (min*86400 .< (times * dt).< max*86400)
end

export days

function filestimes(m:: MITgcmDatas, varfile :: ASCIIString)
  rr = Regex("^$(varfile)\\.0.*\\.data\$")
  filenames = filter( rr, m.files)
  nsteps  = map(filenames) do x
    rr1 = Regex("^$(varfile)\\.0")
    x1 = replace(x, rr1, s"")
    p = match(r"[0-9]+",x1)
    parse(Int64,p.match)
  end
  return (filenames, nsteps )
end

function exactTimeStep(times :: Array{Int64,1}, m :: MITgcmDatas;
                       minstep :: Int64 = 56880, maxstep :: Int64 = 57000)
  return  minstep .< times .<= maxstep
end
export exactTimeStep


text:
import Base.Test: @test
using RNG.MersenneTwisters

stdout_ = STDOUT
pwd_ = pwd()
cd(joinpath(Pkg.dir("RNG"), "test/MersenneTwisters"))
mkpath("./actual")

for mt_name in (:MT19937, )

    outfile = open(string(
        "./actual/check-$(lowercase("$mt_name")).out"
    ), "w")
    redirect_stdout(outfile)

    @eval $mt_name()
    x = @eval $mt_name(123)

    for i in 1:100
        @printf "%.9f\n" rand(x)
    end

    close(outfile)
end
redirect_stdout(stdout_)

@test success(`diff -ru expected actual`)
cd(pwd_)


text:
# Composable Widget API

abstract WidgetMod <: Widget

signal(w::WidgetMod) = signal(w.widget)

immutable Styled{W <: Widget} <: WidgetMod
    widget::W
    style::Dict
end

style(widget, style) = StyledWidget(widget, style)

immutable Labeled{W <: Widget} <: WidgetMod
    widget::W
    label::String
end

label(widget, label) = LabeledWidget(widget, label)

immutable WithClass{class, W <: Widget} <: WidgetMod
    widget::W
end

addclass(widget, class) = WithClass{symbol(class)}(widget)

abstract Container <: Widget

immutable WidgetStack{direction} <: Container
    widgets::Vector{Widget}
end

hstack(widgets::Widget...) = WidgetStack{:x}([widget...])
vstack(widgets::Widget...) = WidgetStack{:y}([widget...])
zstack(widgets::Widget...) = WidgetStack{:z}([widget...])


text:
function digamma(x::Float64) # Use Int instead?
  ccall(("digamma", :libRmath), Float64, (Float64, ), x)
end

function trigamma(x::Float64) # Use Int instead?
  ccall(("trigamma", :libRmath), Float64, (Float64, ), x)
end

function fit(::Type{Bernoulli}, x::Array)
  for i in 1:length(x)
    if !insupport(Bernoulli(), x[i])
      error("Bernoulli observations must be 0/1 values")
    end
  end
  return Bernoulli(mean(x))
end

function fit(::Type{Beta}, x::Array)
  for i in 1:length(x)
    if !insupport(Beta(), x[i])
      error("Bernoulli observations must be in [0,1]")
    end
  end
  x_bar = mean(x)
  v_bar = var(x)
  a = x_bar * (((x_bar * (1. - x_bar)) / v_bar) - 1.)
  b = (1. - x_bar) * (((x_bar * (1. - x_bar)) / v_bar) - 1.)
  return Beta(a, b)
end

function fit(::Type{Binomial}, x::Real, n::Real)
  if x > n || x < 0
    error("For binomial observations, x must lie in [0, n]")
  end
  return Binomial(int(n), x / n)
end

# Categorical

# Cauchy

# Chisq

# Dirichlet

function fit(::Type{DiscreteUniform}, x::Array)
  DiscreteUniform(min(x), max(x))
end

function fit(::Type{Exponential}, x::Array)
  for i in 1:length(x)
    if !insupport(Exponential(), x[i])
      error("Exponential observations must be non-negative values")
    end
  end
  return Exponential(mean(x))
end

# FDist

function fit(::Type{Gamma}, x::Array)
  a_old = 0.5 / (log(mean(x)) - mean(log(x)))
  a_new = 0.5 / (log(mean(x)) - mean(log(x)))

  z = 1.0 / a_old + (mean(log(x)) - log(mean(x)) + log(a_old) - digamma(a_old)) / (a_old^2 * (1.0 / a_old - trigamma(a_old)))
  a_new, a_old = (1.0 / z, a_new)

  iterations = 1

  while abs(a_old - a_new) > 10e-16 && iterations <= 10_000
    z = 1.0 / a_old + (mean(log(x)) - log(mean(x)) + log(a_old) - digamma(a_old)) / (a_old^2 * (1.0 / a_old - trigamma(a_old)))
    a_new, a_old = (1.0 / z, a_new)
  end

  if iterations >= 10_000
    error("Parameter estimation failed to converge in 10,000 iterations")
  end

  Gamma(a_new, mean(x) / a_new)
end

function fit(::Type{Geometric}, x::Array)
  Geometric(1.0 / mean(convert(Array{Int}, x)))
end

# HyperGeometric

# Logistic

function fit(::Type{Laplace}, x::Array)
  a = median(x)
  deviations = 0.0
  for i in 1:length(x)
    deviations += abs(x[i] - a)
  end
  b = deviations / length(x)
  Laplace(a, b)
end

# logNormal

# function fit(Multinomial)
# end

# Assumes observations are rows
# Positive-definite errors
# function fit{T <: Real}(x::Matrix{T}, d::MultivariateNormal)
#   MultivariateNormal(reshape(mean(x, 1), size(x, 2)), cov(x))
# end

# fit([1.0 -0.1; -2.0 0.1], MultivariateNormal())

# NegativeBinomial

# NoncentralBeta

# NoncentralChisq

# NoncentralF

# NoncentralT

# Normal
# Not strict MLE
function fit(::Type{Normal}, x::Array)
  Normal(mean(x), std(x))
end

# Poisson
function fit(::Type{Poisson}, x::Array)
  for i in 1:length(x)
    if !insupport(Poisson(), x[i])
      error("Poisson observations must be non-negative integers")
    end
  end
  Poisson(mean(x))
end

# TDist

# Uniform
function fit(::Type{Uniform}, x::Array)
  Uniform(min(x), max(x))
end

# Weibull


text:
module Space

using ..Util
using ..IDS
using Spec
import Flux     # FIXME: Can we excise this from this submodule?
import Tracker
import ForwardDiff
import Random
import Distributions
using Random: GLOBAL_RNG, AbstractRNG
using DocStringExtensions

export Ω
export ΩBase, memrand, linearize, unlinearize
export ΩProj, parentω, memrand, proj
export TaggedΩ, tag, Tags, hastags, transfertags
export SimpleΩ, LinearΩ, update, nelem
export defΩ, defΩProj

include("omega.jl")         # Sample Space
include("proj.jl")          # Sample Space Projection
include("tagged.jl")        # Space space Tagged with metadata
include("common.jl")
include("seed.jl")          # Setting the seed

# Omega Types
include("simple.jl")        # Simple Ω
include("linear.jl")        # Linear Ω
include("rng.jl")           # Linear Ω
include("defaults.jl")

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    WeightedCrossValidation(weighter, partitioner;
                            lambda=1.0, loss=Dict())

Weighted cross-validation in which samples are split into folds with
`partitioner` method and are weighted with `weighter` method. Weights
are raised to `lambda` power in `[0,1]`. Optionally, specify `loss`
function from `LossFunctions.jl` for some of the variables.

    WeightedCrossValidation(weighter, k; shuffle=true,
                            lambda=1.0, loss=Dict())

Alternatively, specify the desired number of folds `k` and `shuffle`
options for a `UniformPartitioner`.

## References

* Sugiyama et al. 2006. [Importance-weighted cross-validation for
  covariate shift](https://link.springer.com/chapter/10.1007/11861898_36)
* Sugiyama et al. 2007. [Covariate shift adaptation by importance weighted
  cross validation](http://www.jmlr.org/papers/volume8/sugiyama07a/sugiyama07a.pdf)
"""
struct WeightedCrossValidation{W<:AbstractWeighter,
                               P<:AbstractPartitioner,
                               T<:Real} <: AbstractErrorEstimator
  weighter::W
  partitioner::P
  lambda::T
  loss::Dict{Symbol,SupervisedLoss}

  function WeightedCrossValidation{W,P,T}(weighter, partitioner, lambda, loss) where {W,P,T}
    @assert 0 ≤ lambda ≤ 1 "lambda must lie in [0,1]"
    new(weighter, partitioner, lambda, loss)
  end
end

WeightedCrossValidation(weighter::W, partitioner::P;
                        lambda=1.0, loss=Dict()) where {W<:AbstractWeighter,
                                                        P<:AbstractPartitioner} =
  WeightedCrossValidation{W,P,typeof(lambda)}(weighter, partitioner, lambda, loss)

WeightedCrossValidation(weighter::W, k::Int; shuffle=true,
                        lambda=1.0, loss=Dict()) where {W<:AbstractWeighter} =
  WeightedCrossValidation(weighter, UniformPartitioner(k, shuffle), lambda=lambda, loss=loss)

function error(solver::AbstractLearningSolver,
               problem::LearningProblem,
               eestimator::WeightedCrossValidation)
  # retrieve problem info
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  partitioner = eestimator.partitioner
  weighter = eestimator.weighter
  lambda = eestimator.lambda
  loss  = eestimator.loss
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[var][1])
    end
  end

  # weight all samples
  weights = weight(sdata, weighter) .^ lambda

  # folds for cross-validation
  folds  = subsets(partition(sdata, partitioner))
  nfolds = length(folds)

  # error for a fold k
  function ε(k)
    # source and target indices
    sinds = [ind for i in [1:k-1; k+1:nfolds] for ind in folds[i]]
    tinds = folds[k]

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solution   = solve(subproblem, solver)

    # loss for each variable
    losses = map(ovars) do var
      y = hold[var]
      ŷ = solution[var]
      w = view(weights, tinds)
      ℒ = value(loss[var], y, ŷ, AggMode.WeightedSum(w)) / length(y)
      var => ℒ
    end

    Dict(losses)
  end

  # compute error for each fold in parallel
  εs = foldxt(vcat, Map(ε), 1:nfolds)

  # combine error from different folds
  Dict(var => mean(get.(εs, var, 0)) for var in ovars)
end


text:
#!/usr/bin/env julia

using ...ReconstructionEvaluations
using ..HDF5

const RE = ReconstructionEvaluations

function spvec_to_arrays(X::SparseVector)
  X = dropzeros(X)
  find(X), nonzeros(X)
end

function fill_array_indices( is1, vs1, is2, vs2 )

  is = sort(union(is1,is2))

  l1, l2 = length(vs1), length(vs2)
  v1, v2 = zeros(eltype(vs1),length(is)), zeros(eltype(vs2),length(is))

  p1,p2 = 1,1# current "pointer" to vs1&2
  for (i,v) in enumerate(is)
    if p1 <= l1 && v == is1[p1]
      v1[i] = vs1[p1]
      p1 += 1
    end

    if p2 <= l2 && v == is2[p2]
      v2[i] = vs2[p2]
      p2 += 1
    end
  end

  is, v1, v2
end

function append_score( score, score_record )
  scores = HDF5.h5read(score_record, "/main")

  push!(scores, score)

  RE.write_h5(scores, score_record)
end


text:
knn_search(dist_mat, k, metric::Symbol) = knn_search(dist_mat, k, Val(metric))

"""
    knn_search(dist_mat, k, :precomputed) -> knns, dists
Find the `k` nearest neighbors of each point in a precomputed distance
matrix.
"""
knn_search(dist_mat, k, ::Val{:precomputed}) = _knn_from_dists(dist_mat, k)

"""
    knn_search(X, k, metric) -> knns, dists
Find the `k` nearest neighbors of each point in `X` by `metric`.
"""
function knn_search(X,
                    k,
                    metric::SemiMetric)
    if size(X, 2) < 4096
        return knn_search(X, k, metric, Val(:pairwise))
    else
        return knn_search(X, k, metric, Val(:approximate))
    end
end

# compute all pairwise distances
# return the nearest k to each point v, other than v itself
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:pairwise}) where {S <: Real}
    num_points = size(X, 2)
    dist_mat = Array{S}(undef, num_points, num_points)
    pairwise!(dist_mat, metric, X, dims=2)
    # all_dists is symmetric distance matrix
    return _knn_from_dists(dist_mat, k)
end

# find the approximate k nearest neighbors using NNDescent
function knn_search(X::AbstractMatrix{S},
                    k,
                    metric,
                    ::Val{:approximate}) where {S <: Real}
    knngraph = nndescent(X, k, metric)
    return knn_matrices(knngraph)
end

function _knn_from_dists(dist_mat::AbstractMatrix{S}, k) where {S <: Real}
    knns_ = [partialsortperm(view(dist_mat, :, i), 2:k+1) for i in 1:size(dist_mat, 1)]
    dists_ = [dist_mat[:, i][knns_[i]] for i in eachindex(knns_)]
    knns = hcat(knns_...)::Matrix{Int}
    dists = hcat(dists_...)::Matrix{S}
    return knns, dists
end


text: repeatstr(n, s) = s^n

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module MultivariateWrapper

using JavaCall

const AbstractMultivariateRealDistribution = @jimport org.hipparchus.distribution.multivariate.AbstractMultivariateRealDistribution
const Class = @jimport java.lang.Class
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const MixtureMultivariateNormalDistribution = @jimport org.hipparchus.distribution.multivariate.MixtureMultivariateNormalDistribution
const MixtureMultivariateRealDistribution = @jimport org.hipparchus.distribution.multivariate.MixtureMultivariateRealDistribution
const MultivariateNormalDistribution = @jimport org.hipparchus.distribution.multivariate.MultivariateNormalDistribution
const MultivariateRealDistribution = @jimport org.hipparchus.distribution.MultivariateRealDistribution
const Object = @jimport java.lang.Object
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix

export AbstractMultivariateRealDistribution
export MixtureMultivariateNormalDistribution
export MixtureMultivariateRealDistribution
export MultivariateNormalDistribution
export density
export equals
export get_class
export get_components
export get_covariances
export get_dimension
export get_means
export get_standard_deviations
export hash_code
export notify
export notify_all
export reseed_random_generator
export sample
export to_string
export wait

include("abstract_multivariate_real_distribution.jl")
include("mixture_multivariate_normal_distribution.jl")
include("mixture_multivariate_real_distribution.jl")
include("multivariate_normal_distribution.jl")

end


text:
"""
nonfixedIncome
(Temporary income)
any income from household  - housemembers- that is Temporary
"""

mutable struct houseHoldIncome :> nonfixedIncome

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
module IJuliaConverter
using JSON

include("jupyter2code.jl")

export jupyter2code

end # module

max_stars_repo_name: matthew-lake/IJuliaConverter

max_issues_repo_name: matthew-lake/IJuliaConverter

max_forks_repo_name: matthew-lake/IJuliaConverter


text:
struct PreplottedTree
    xnode::Vector{Float64}
    ynode::Vector{Float64}
    xbranch::Vector{Float64}
    ybranch::Vector{Float64}
    labelnode::Vector{String}
    #=
    TODO: Add fields for branch ID data, colour, width
    =#
end

"""
Get X coords of all the nodes, and Y coords of the outer nodes.
"""
function _tree_coords_pass1(p, A, x, y, usebrl)
    id = getid(p)

    if usebrl
        brlen = brlength(p)
        brlen = isnothing(brlen) ? 1.0 : brlen
        x += brlen
    else
        x += 1.0
    end
    A.xnode[id] = x

    if isouter(p)
        A.ynode[id] = y
        y -= 1-0
    end
    for c in children(p)
        y = _tree_coords_pass1(c, A, x, y, usebrl)
    end

    return y
end

function _elbow_brach_coords(A, idp, idc)
    xp, yp = A.xnode[idp], A.ynode[idp]
    xc, yc = A.xnode[idc], A.ynode[idc]
    push!(A.xbranch, xp, xp, xc, NaN)
    push!(A.ybranch, yp, yc, yc, NaN)

    return nothing
end

function _straight_brach_coords(A, idp, idc)
    xp = A.xnode[idp]
    xc, yc = A.xnode[idc], A.ynode[idc]
    push!(A.xbranch, xp, xc, NaN)
    push!(A.ybranch, yc, yc, NaN)

    return nothing
end

function _tree_coords_pass2(p, A, y)
    ychildrensum = 0.0
    for c in children(p)
        y = isinner(c) ? _tree_coords_pass2(c, A, y) : A.ynode[getid(c)]
        ychildrensum += y
    end

    id = getid(p)
    nc = outdegree(p)
    y = ychildrensum / nc
    A.ynode[id] = y

    if nc > 1
        for (i, c) in enumerate(children(p))
            idc = getid(c)
            if ((i == 1) || (i == nc))    # Draw elbow
                _elbow_brach_coords(A, id, idc)
            else                          # Draw straight branch
                _straight_brach_coords(A, id, idc)
            end
        end
    else
        idc = getid(first(children(p)))
        _straight_brach_coords(A, id, idc)
    end

    return y
end

function _preplot_tree(tree, usebrl)
    # ni = number of inner nodes, no = number of outer nodes
    ni, no = mapreduce(x -> isinner(x) ? (1, 0) : (0, 1), .+, getnodes(tree))
    n = ni + no

    xnode = zeros(Float64, n)
    ynode = zeros(Float64, n)
    xbranch = Float64[]
    ybranch = Float64[]
    labelnode = getlabel.(getnodes(tree))

    # Preallocate guessing that all branches are elbows: 3 coords and NaN
    sizehint!(xbranch, 4 * n)
    sizehint!(ybranch, 4 * n)

    A = PreplottedTree(xnode, ynode, xbranch, ybranch, labelnode)

    y = Float64(no)    # Start with a height (`y`) of `no`

    for c in children(getroot(tree))
        y =_tree_coords_pass1(c, A, 0.0, y, usebrl)
    end
    A.xnode[getid(getroot(tree))] = 0.0

    _tree_coords_pass2(getroot(tree), A, NaN)

    return A
end

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
module BipartiteGraphs

export BipartiteEdge, BipartiteGraph

export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
       𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST

using DocStringExtensions
using UnPack
using SparseArrays
using Graphs
using Setfield

###
### Edges & Vertex
###
@enum VertType SRC DST ALL

struct BipartiteEdge{I<:Integer} <: Graphs.AbstractEdge{I}
    src::I
    dst::I
    function BipartiteEdge(src::I, dst::V) where {I,V}
        T = promote_type(I, V)
        new{T}(T(src), T(dst))
    end
end

Graphs.src(edge::BipartiteEdge) = edge.src
Graphs.dst(edge::BipartiteEdge) = edge.dst

function Base.show(io::IO, edge::BipartiteEdge)
    @unpack src, dst = edge
    print(io, "[src: ", src, "] => [dst: ", dst, "]")
end

Base.:(==)(a::BipartiteEdge, b::BipartiteEdge) = src(a) == src(b) && dst(a) == dst(b)

###
### Graph
###
"""
$(TYPEDEF)

A bipartite graph representation between two, possibly distinct, sets of vertices
(source and dependencies). Maps source vertices, labelled `1:N₁`, to vertices
on which they depend (labelled `1:N₂`).

# Fields
$(FIELDS)

# Example
```julia
using ModelingToolkit

ne = 4
srcverts = 1:4
depverts = 1:2

# six source vertices
fadjlist = [[1],[1],[2],[2],[1],[1,2]]

# two vertices they depend on
badjlist = [[1,2,5,6],[3,4,6]]

bg = BipartiteGraph(7, fadjlist, badjlist)
```
"""
mutable struct BipartiteGraph{I<:Integer,F<:Vector{Vector{I}},B<:Union{Vector{Vector{I}},I},M} <: Graphs.AbstractGraph{I}
    ne::Int
    fadjlist::F # `fadjlist[src] => dsts`
    badjlist::B # `badjlist[dst] => srcs` or `ndsts`
    metadata::M
end
BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) = BipartiteGraph(ne, fadj, badj, metadata)

"""
```julia
Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
```

Test whether two [`BipartiteGraph`](@ref)s are equal.
"""
function Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
    iseq = (bg1.ne == bg2.ne)
    iseq &= (bg1.fadjlist == bg2.fadjlist)
    iseq &= (bg1.badjlist == bg2.badjlist)
    iseq
end

"""
$(SIGNATURES)

Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations.
"""
function BipartiteGraph(nsrcs::T, ndsts::T, backedge::Val{B}=Val(true); metadata=nothing) where {T,B}
    fadjlist = map(_->T[], 1:nsrcs)
    badjlist = B ? map(_->T[], 1:ndsts) : ndsts
    BipartiteGraph(0, fadjlist, badjlist, metadata)
end

Base.eltype(::Type{<:BipartiteGraph{I}}) where I = I
function Base.empty!(g::BipartiteGraph)
    foreach(empty!, g.fadjlist)
    g.badjlist isa AbstractVector && foreach(empty!, g.badjlist)
    g.ne = 0
    if g.metadata !== nothing
        foreach(empty!, g.metadata)
    end
    g
end
Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")

@noinline throw_no_back_edges() = throw(ArgumentError("The graph has no back edges."))

if isdefined(Graphs, :has_contiguous_vertices)
    Graphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false
end
Graphs.is_directed(::Type{<:BipartiteGraph}) = false
Graphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1)
𝑑vertices(g::BipartiteGraph) = g.badjlist isa AbstractVector ? axes(g.badjlist, 1) : Base.OneTo(g.badjlist)
has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g)
has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g)
𝑠neighbors(g::BipartiteGraph, i::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i]
function 𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M
    g.badjlist isa AbstractVector || throw_no_back_edges()
    M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
end
Graphs.ne(g::BipartiteGraph) = g.ne
Graphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
Graphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}

nsrcs(g::BipartiteGraph) = length(𝑠vertices(g))
ndsts(g::BipartiteGraph) = length(𝑑vertices(g))

function Graphs.has_edge(g::BipartiteGraph, edge::BipartiteEdge)
    @unpack src, dst = edge
    (src in 𝑠vertices(g) && dst in 𝑑vertices(g)) || return false  # edge out of bounds
    insorted(𝑠neighbors(src), dst)
end

###
### Populate
###
struct NoMetadata
end
const NO_METADATA = NoMetadata()

Graphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md=NO_METADATA) = add_edge!(g, BipartiteEdge(i, j), md)
function Graphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_METADATA)
    @unpack fadjlist, badjlist = g
    s, d = src(edge), dst(edge)
    (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
    @inbounds list = fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)
    if md !== NO_METADATA
        insert!(g.metadata[s], index, md)
    end

    g.ne += 1
    if badjlist isa AbstractVector
        @inbounds list = badjlist[d]
        index = searchsortedfirst(list, s)
        insert!(list, index, s)
    end
    return true  # edge successfully added
end

function Graphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where T
    if type === DST
        if g.badjlist isa AbstractVector
            push!(g.badjlist, T[])
        else
            g.badjlist += 1
        end
    elseif type === SRC
        push!(g.fadjlist, T[])
    else
        error("type ($type) must be either `DST` or `SRC`")
    end
    return true  # vertex successfully added
end

###
### Edges iteration
###
Graphs.edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(ALL))
𝑠edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC))
𝑑edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(DST))

struct BipartiteEdgeIter{T,G} <: Graphs.AbstractEdgeIter
    g::G
    type::Val{T}
end

Base.length(it::BipartiteEdgeIter) = ne(it.g)
Base.length(it::BipartiteEdgeIter{ALL}) = 2ne(it.g)

Base.eltype(it::BipartiteEdgeIter) = edgetype(it.g)

function Base.iterate(it::BipartiteEdgeIter{SRC,<:BipartiteGraph{T}}, state=(1, 1, SRC)) where T
    @unpack g = it
    neqs = nsrcs(g)
    neqs == 0 && return nothing
    eq, jvar = state

    while eq <= neqs
        eq′ = eq
        vars = 𝑠neighbors(g, eq′)
        if jvar > length(vars)
            eq += 1
            jvar = 1
            continue
        end
        edge = BipartiteEdge(eq′, vars[jvar])
        state = (eq, jvar + 1, SRC)
        return edge, state
    end
    return nothing
end

function Base.iterate(it::BipartiteEdgeIter{DST,<:BipartiteGraph{T}}, state=(1, 1, DST)) where T
    @unpack g = it
    nvars = ndsts(g)
    nvars == 0 && return nothing
    ieq, jvar = state

    while jvar <= nvars
        eqs = 𝑑neighbors(g, jvar)
        if ieq > length(eqs)
            ieq = 1
            jvar += 1
            continue
        end
        edge = BipartiteEdge(eqs[ieq], jvar)
        state = (ieq + 1, jvar, DST)
        return edge, state
    end
    return nothing
end

function Base.iterate(it::BipartiteEdgeIter{ALL,<:BipartiteGraph}, state=nothing)
    if state === nothing
        ss = iterate((@set it.type = Val(SRC)))
    elseif state[3] === SRC
        ss = iterate((@set it.type = Val(SRC)), state)
    elseif state[3] == DST
        ss = iterate((@set it.type = Val(DST)), state)
    end
    if ss === nothing && state[3] == SRC
        return iterate((@set it.type = Val(DST)))
    else
        return ss
    end
end

###
### Utils
###
function Graphs.incidence_matrix(g::BipartiteGraph, val=true)
    I = Int[]
    J = Int[]
    for i in 𝑠vertices(g), n in 𝑠neighbors(g, i)
        push!(I, i)
        push!(J, n)
    end
    S = sparse(I, J, val, nsrcs(g), ndsts(g))
end

end # module


text:
export TabularDeterministicPolicy

"""
    TabularDeterministicPolicy(table::Vector{Int}, nactions::Int) -> TabularDeterministicPolicy

Always select the specific action according to the state defined in `table`.
`nactions` is used to define the number of possible actions.
"""
Base.@kwdef struct TabularDeterministicPolicy <: AbstractPolicy
    table::Vector{Int}
    nactions::Int
end

(p::TabularDeterministicPolicy)(s::Int) = p.table[s]
(p::TabularDeterministicPolicy)(obs::Observation) = p(get_state(obs))
get_prob(p::TabularDeterministicPolicy, s::Int, a::Int) = p.table[s] == a ? 1.0 : 0.0

function update!(p::TabularDeterministicPolicy, correction)
    s, a = correction
    p.table[s] = a
end


text:
# Construct contexts for up to 7 inputs if MyInt = UInt128
# Example of context for 2 inputs
# Context[1] = 0x0c = 1100
# Context[2] = 0x0a = 1010
# We consider this as a 2 by 4 bit matrix where Context[i] is the ith row.
# Each column of this matrx corresponds to 2 inputs to the circuit.
# All columns of the matrix correspond to all possible inputs to the circuit.
# Column 1 corresponds to input1 = 1, input2 = 1
# Column 2 corresponds to input1 = 1, input2 = 0
# Column 3 corresponds to input1 = 0, input2 = 1
# Column 4 corresponds to input1 = 0, input2 = 0
# Thus, the 4-bit output string describes the output of the circuit for all inputs.
# This specifies a truth table where columns of the Context matrix correspoing to
#   rows of the truth table.

# When there are 3 inputs, Context has 3 components, one for each input.
# Thus, the matrix has 3 rows and 8 = 2^3 columns and the function computed by the
#   circuit is specified by the 8 bit output string.

# The logic function multi-output circuit is specified by sequence of output bit
#   strings for each output, or by the concatenation of these bit strings into
#   a single bit string.

export construct_ones, construct_contexts, construct_context, get_bits

# Ones[i] is a "mask" with 1 bits in positions 1 to 2^i where position 1 is least significant
function construct_ones( numinputs::Integer )
  Ones = zeros(MyInt,numinputs)
  Ones[1] = 3
  for i = 2:numinputs
    Ones[i] = Ones[i-1] | (Ones[i-1] << 2^(i-1))
  end
  Ones
end

function construct_contexts( numinputs::Integer )
  Ones = construct_ones( numinputs )
  Contexts = [ zeros( MyInt, i ) for i = 1:numinputs ]
  Contexts[1][1] = 0x2
  for i = 2:numinputs
    Contexts[i][1] = Ones[i-1] << 2^(i-1)
    for j = 2:i
      Contexts[i][j] = Contexts[i-1][j-1] << 2^(i-1) | Contexts[i-1][j-1]
      #println("(i,j):",(i,j),Contexts[i][j])
    end
  end
  Contexts
end

function construct_context( numinputs::Integer )
  construct_contexts( numinputs )[numinputs]
end

# v is a vector of outputs from a subset of a circuit
# The length of v is the number of gates in the circuit
# Each element of v is interpreted as a bitstring of length numinputs
# The result is a vector of bitstrings of length 2^numinputs
#   where result[i][j] is
# Example:   let v = [0xe, 0x5, 0xa] = [1110, 0101, 1010]
# get_bits(v,2) = [0x2, 0x5, 0x6, 0x5] = [0010, 0101, 0110, 0101]
# get_bits can be interpreted as a transpose operation on the bit matrix:
#    1110
#    0101
#    1010
# The columns of this matrix are the reversed result
function get_bits( v::Vector{MyInt}, numinputs::Integer )
  result = zeros(MyInt,2^numinputs)
  reverse_v = reverse(v)
  mask = 1
  mask_shift = 0
  for i = 1:2^numinputs
    shift = 0
    for j = 1:length(v)
      result[i] |= (reverse_v[j] & mask) >> (mask_shift - shift)
      shift += 1
    end
    mask <<= 1
    mask_shift += 1
  end
  result
end


text:
"""
# QuakeML

Read and write files in the QuakeML format which describes seismic events.

## User functions

- [`QuakeML.read`](@ref): Read a QuakeML file
- [`QuakeML.readstring`](@ref): Read a QuakeML document from a string
- [`preferred_focal_mechanism`](@ref): Get the preferred focal mechanism for an event
- [`preferred_magnitude`](@ref): Get the preferred magnitude for an event
- [`preferred_origin`](@ref): Get the preferred origin for an event
- [`quakeml`](@ref): Create an XML document from a set of events which can
  be written with `print(io, quakeml(qml))`
"""
module QuakeML

using Dates: DateTime, Time
import UUIDs

import EzXML

export
    EventParameters,
    has_focal_mechanism,
    has_magnitude,
    has_origin,
    preferred_focal_mechanism,
    preferred_focal_mechanisms,
    preferred_magnitude,
    preferred_magnitudes,
    preferred_origin,
    preferred_origins,
    quakeml

include("compat.jl")
include("util.jl")
include("types.jl")
include("io.jl")
include("constructors.jl")
include("accessors.jl")
include("types_module.jl")

end # module

max_stars_repo_name: anowacki/QuakeML


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 21c74f60-fe61-4b5f-8fe5-1840b0a4bcb0
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 705f3bc2-ea0a-11eb-2d81-21de5db26c02
begin
	using CSV
	using DataFrames
	using WordCloud
end

# ╔═╡ 30b58656-7d24-4ff3-908f-c372de38fc9b
md"
## Preprocess
"

# ╔═╡ c9d86708-8ce9-48dd-8d18-ba04b4c4a22d
sms2token(sms::String) = [m.match for m in eachmatch(r"[a-z']+", lowercase(sms))]

# ╔═╡ a7b1da07-fa68-44c9-b01d-b6b0196f02b1
begin
	df_raw = CSV.read("../data/SMSSpamCollection", DataFrame, header=[:Label, :SMS])

	df = copy(df_raw)
	df[!, :SMS] .= sms2token.(df.SMS)

	spam_df = df[df.Label.=="spam", :]
	ham_df = df[df.Label.=="ham", :]

	# df_raw
	df
end

# ╔═╡ 6d356343-6f13-4792-b084-b55aa0af7510
md"
## Text mining
"

# ╔═╡ b4a30d03-5f3c-4fb3-9351-9c8f689032a2
begin
	features = unique(vcat(df.SMS...))

	function calc_tf(words::Vector)
		tf = Dict([(w, Float64(count(x->x==w, words))) for w in features])

		n = sum(values(tf))
		for (k, v) in tf
			tf[k] = (v == 0) ? 1e-1/n : v/n
		end

		return tf
	end
end

# ╔═╡ f10b2755-70d1-4cfe-a759-17f283bb7cb3
begin
	spam_p = nrow(spam_df) / nrow(df)
	spam_tf = calc_tf(vcat(spam_df.SMS...))

	ham_p = nrow(ham_df) / nrow(df)
	ham_tf = calc_tf(vcat(ham_df.SMS...))
end;

# ╔═╡ ac757e08-0b19-4aa8-a791-8a677280b12f
function plot_word_cloud(tf::Dict, top_n::Integer)
	tf = sort(collect(tf), by=x->x[2])

	return wordcloud(
		map(x->x[1], tf)[1:top_n],
		map(x->x[2], tf)[1:top_n],
		font="LiSong Pro",
		mask=shape(box, 400*2, 300*2, 10*2),
		colors=:Dark2_3,
		angles=(0, 90),
		density=0.4
	)
end

# ╔═╡ 7f30886f-32e3-498c-8166-ea5f18b5e725
md"
### Non-spam word cloud
"

# ╔═╡ d9127fa7-e83d-4c21-a2d0-9fa497e1af44
plot_word_cloud(spam_tf, 400)

# ╔═╡ 1fca2c23-5ff9-42cc-a5fa-a872d4cf0cab
md"
### Spam word cloud
"

# ╔═╡ b4d6a0ed-8e41-44d1-b736-e3acc1cd511e
plot_word_cloud(ham_tf, 400)

# ╔═╡ e7997329-45fa-42cb-8d6b-399089099110
md"
## Model
"

# ╔═╡ dcb7da4c-5092-4667-845a-af4d3cd7c9d3
md"
$P(Cs|\vec{x}) = \frac{P(\vec{x}|Cs) P(Cs)}{P(\vec{x})}$

$P(Cs|\vec{x}) = \frac{P(\vec{x}|Cs) P(Cs)}{P(\vec{x}|Cs)P(Cs) + P(\vec{x}|Ch)P(Ch)}$

$P(\vec{x}|C) = P(w1|C)P(w2|C)P(w3|C)...$
"

# ╔═╡ a3f53b05-f2dd-46af-b573-3c5de2c3f709
begin
	function p_x_c(tokens::Vector, tf::Dict)
		p = 1
		for t in tokens
			p *= tf[t]
		end

		return p
	end

	function infer(s::String)
		tokens = filter(x->x in features, sms2token(s))
		(isempty(tokens)) && (return 0.5)

		return p_x_c(tokens, spam_tf)*spam_p / (
			p_x_c(tokens, spam_tf)*spam_p + p_x_c(tokens, ham_tf)*ham_p
		)
	end
end

# ╔═╡ f95f4e29-44ec-48ae-8ef3-878c3f88ba4c
md"
## Validation

Is the following message a spam?
"

# ╔═╡ 60b39235-a767-452e-998a-6f6efa0a7c1e
infer(
	"Congratulations! You have won a free launch today~~~"
)

# ╔═╡ 393e099f-135c-4861-8581-fb5a48066df4
infer(
	"I am the most handsome guy around the world!!"
)

# ╔═╡ Cell order:
# ╟─21c74f60-fe61-4b5f-8fe5-1840b0a4bcb0
# ╠═705f3bc2-ea0a-11eb-2d81-21de5db26c02
# ╟─30b58656-7d24-4ff3-908f-c372de38fc9b
# ╠═c9d86708-8ce9-48dd-8d18-ba04b4c4a22d
# ╠═a7b1da07-fa68-44c9-b01d-b6b0196f02b1
# ╟─6d356343-6f13-4792-b084-b55aa0af7510
# ╠═b4a30d03-5f3c-4fb3-9351-9c8f689032a2
# ╠═f10b2755-70d1-4cfe-a759-17f283bb7cb3
# ╠═ac757e08-0b19-4aa8-a791-8a677280b12f
# ╟─7f30886f-32e3-498c-8166-ea5f18b5e725
# ╠═d9127fa7-e83d-4c21-a2d0-9fa497e1af44
# ╟─1fca2c23-5ff9-42cc-a5fa-a872d4cf0cab
# ╠═b4d6a0ed-8e41-44d1-b736-e3acc1cd511e
# ╟─e7997329-45fa-42cb-8d6b-399089099110
# ╟─dcb7da4c-5092-4667-845a-af4d3cd7c9d3
# ╠═a3f53b05-f2dd-46af-b573-3c5de2c3f709
# ╟─f95f4e29-44ec-48ae-8ef3-878c3f88ba4c
# ╠═60b39235-a767-452e-998a-6f6efa0a7c1e
# ╠═393e099f-135c-4861-8581-fb5a48066df4


text:
### ------------ compute contact forces/torques --------------------------------
# functions setComputationFlag(), selectContactPairsWithEvent!(),
# selectContactPairsNoEvent!() and getDistances!() are called to compute distances
# between objs with mpr algorithm, and use them as zero-crossing functions
# (see file ...\ContactDetectionMPR\handler.jl)
# after computing all distances between colliding shapes response calculations are done
function computeContactForcesAndTorques(sim::ModiaLang.SimulationModel, scene, world, time, file)
  ch = scene.options.contactDetection
  # Compute signed distances of all contact shapes during zero-crossing computation
    #try
        setComputationFlag(ch)
        if ModiaLang.isEvent(sim)             # with Event
            # handle event
            TimerOutputs.@timeit sim.timer "Modia3D_1 selectContactPairsWithEvent!"  selectContactPairsWithEvent!(sim, scene, ch, world)
        elseif ModiaLang.isZeroCrossing(sim)  # no Event
            # compute zero crossings
            TimerOutputs.@timeit sim.timer "Modia3D_1 selectContactPairsNoEvent!" selectContactPairsNoEvent!(sim, scene, ch, world)
        else
            # compute contact distances
            TimerOutputs.@timeit sim.timer "Modia3D_1 getDistances!"  getDistances!(scene, ch, world)
        end
        TimerOutputs.@timeit sim.timer "Modia3D_1 dealWithContacts!"  dealWithContacts!(sim, scene, ch, world, time, file)
    #catch e
        #println("Bug in Modia3D/src/Composition/dynamicCollision.jl: error curing collision detection at time = ", time)
    #end
end

# each shape pair in contact is stored in ch.contactDict
# if there is an event simulation status is updated and for response calculation
# a contact material is defined
# contact and/or support points are visualized
# only if enableContactDetection is true responseCalculation for force and torque computation is called
# (see file ...\Composition\responseCalculation\elasticCollisionResponse.jl)
# further, at an event simulation status is updated, contact material is replaced
# and the actual contactDict is stored
function dealWithContacts!(sim, scene, ch, world, time, file)
  simh = sim.eventHandler
  for (pairID, pair) in ch.contactDict
    obj1 = pair.obj1
    obj2 = pair.obj2
    rContact      = (pair.contactPoint1 + pair.contactPoint2)/2.0
    contactNormal = pair.contactNormal
    if ModiaLang.isEvent(sim)
      # println("$(sim.time): ", obj1.path, " ", obj2.path)
      getMaterialContactStart(scene, ch, simh, pair, pairID, obj1, obj2, rContact, contactNormal)
      visualizeContactAndSupportPoints(ch, world)
    end
    if scene.options.enableContactDetection
      (f1,f2,t1,t2) = responseCalculation(pair.contactPairMaterial, obj1, obj2,
                                          rContact, contactNormal,
                                          pair.distanceWithHysteresis, time, file, sim)

      # Transform forces/torques in local part frames
      obj1.f += obj1.R_abs*f1
      obj1.t += obj1.R_abs*t1
      obj2.f += obj2.R_abs*f2
      obj2.t += obj2.R_abs*t2
  end; end

  if ModiaLang.isEvent(sim)
    deleteMaterialLastContactDictContactEnd(scene, ch, simh)
end; end

# mainly its for assigning contact materials.
# therefore, lastContactDict stores information of shape pairs in contact at last step
# if the actual shape pair was in contact at the last step as well, the already assigned contact material is taken
# otherwise, contact material between these two shapes is choosen and simulation status is set
function getMaterialContactStart(scene, ch, simh, pair, pairID, obj1, obj2, rContact, contactNormal)
  if haskey(ch.lastContactDict, pairID)
    # use material (reference) from previous event
    if scene.options.enableContactDetection
      pair.contactPairMaterial = ch.lastContactDict[pairID].contactPairMaterial    # improve later (should avoid to inquire pairID twice)
    end
  else
    # determine contact pair material
    if scene.options.enableContactDetection
      pair.contactPairMaterial = contactStart(obj1, obj2, rContact, contactNormal,
                                            scene.options.elasticContactReductionFactor)
    end
    simh.restart = max(simh.restart, ModiaLang.Restart)
    simh.newEventIteration = false
    if scene.options.enableContactDetection
      logEvents(simh, pair, obj1, obj2, rContact, contactNormal)
    end
end; end

# important feature is stored at an event
function logEvents(simh, pair, obj1, obj2, rContact, contactNormal)
  if simh.logEvents
    name1 = Modia3D.fullName(obj1)
    name2 = Modia3D.fullName(obj2)
    n     = contactNormal
    println("        distance(", name1, ",", name2, ") = ", pair.distanceWithHysteresis, " became < 0")

    @inbounds println("            contact normal = [", round(n[1], sigdigits=3), ", ", round(n[2], sigdigits=3), ", ", round(n[3], sigdigits=3), "], contact position = [", round(rContact[1], sigdigits=3), ", ", round(rContact[2], sigdigits=3), ", ", round(rContact[3], sigdigits=3),"], c_res = ", round(pair.contactPairMaterial.c_res, sigdigits=3) , " d_res = ", round(pair.contactPairMaterial.d_res, sigdigits=3), "\n")
end; end

# each contact shape pair which was in lastContactDict but isn't in actual contactDict
# (= shape pair was in contact at last step but isn't at the actual step)
# contact material is "deleted", it is set to nothing
# and simulation status is set.
# lastContactDict is deleted and filled with feature of the actual contactDict
function deleteMaterialLastContactDictContactEnd(scene, ch, simh)
  for (pairID, pair) in ch.lastContactDict
    if !haskey(ch.contactDict, pairID)
      if scene.options.enableContactDetection
        contactEnd(pair.contactPairMaterial, pair.obj1, pair.obj2)
      end
      simh.restart = max(simh.restart, ModiaLang.Restart)
      simh.newEventIteration = false
      if !simh.logEvents    # ModiaLang.isLogEvents(simh.logger)
        break
      end
      name1 = Modia3D.fullName(pair.obj1)
      name2 = Modia3D.fullName(pair.obj2)
      println("        distance(", name1, ",", name2, ")  became > 0")
  end; end

  # delete lastContactDict and save contactDict in lastContactDict
  empty!(ch.lastContactDict)
  for (pairID, pair) in ch.contactDict
    push!(ch.lastContactDict, pairID => pair)
end; end

### ------- visualizing contact and support points -----------------------------
# this function visualizes contact and support points if visualizeContactPoints and/or
# visualizeSupportPoints is true
# if a contact point and/or support point is nothing the point is transparent and therefore invisible
# otherwise its visualized as a sphere
function visualizeContactAndSupportPoints(ch, world)
    contactDictCollect = collect(ch.contactDict)
    lengthDict = length(ch.contactDict)

    if ch.visualizeContactPoints
        maxLength = length(world.contactVisuObj1)
        @inbounds for i=1:maxLength
            if i <= lengthDict
                transparency = 0.0
                point1 = contactDictCollect[i][2].contactPoint1
                point2 = contactDictCollect[i][2].contactPoint2
            else
                transparency = 1.0
                point1 = Modia3D.ZeroVector3D
                point2 = Modia3D.ZeroVector3D
            end
            setVisualizationContactProperties!(world.contactVisuObj1[i], transparency, point1)
            setVisualizationContactProperties!(world.contactVisuObj2[i], transparency, point2)
        end
        if lengthDict > maxLength
            nVisualContSupPoints = lengthDict
            printWarnContSupPoints(nVisualContSupPoints)
        end
    end

    if ch.visualizeSupportPoints
        maxLength = length(world.supportVisuObj1A)
        @inbounds for i=1:maxLength
            if i <= lengthDict && contactDictCollect[i][2].supportPointsDefined
                transparency = 0.0
                support1A = contactDictCollect[i][2].support1A
                support2A = contactDictCollect[i][2].support2A
                support3A = contactDictCollect[i][2].support3A
                support1B = contactDictCollect[i][2].support1B
                support2B = contactDictCollect[i][2].support2B
                support3B = contactDictCollect[i][2].support3B
            else
                transparency = 1.0
                support1A = Modia3D.ZeroVector3D
                support2A = Modia3D.ZeroVector3D
                support3A = Modia3D.ZeroVector3D
                support1B = Modia3D.ZeroVector3D
                support2B = Modia3D.ZeroVector3D
                support3B = Modia3D.ZeroVector3D
            end
            setVisualizationContactProperties!(world.supportVisuObj1A[i], transparency, support1A)
            setVisualizationContactProperties!(world.supportVisuObj2A[i], transparency, support2A)
            setVisualizationContactProperties!(world.supportVisuObj3A[i], transparency, support3A)
            setVisualizationContactProperties!(world.supportVisuObj1B[i], transparency, support1B)
            setVisualizationContactProperties!(world.supportVisuObj2B[i], transparency, support2B)
            setVisualizationContactProperties!(world.supportVisuObj3B[i], transparency, support3B)
        end
        if lengthDict > maxLength
            nVisualContSupPoints = lengthDict
            printWarnContSupPoints(nVisualContSupPoints)
        end
    end
end

printWarnContSupPoints(nVisualContSupPoints) = @warn("If all contact points and/or support points should be visualized please set nVisualContSupPoints = $nVisualContSupPoints in SceneOptions.")

# the sphere is visible and its absolute position is updated (this point was computed by the mpr algorithm)
function setVisualizationContactProperties!(obj::Composition.Object3D, transparency::Float64, point::SVector{3,Float64})
    obj.r_abs = point
    obj.feature.visualMaterial.transparency = transparency
end


text:
struct Fibonacci{T>:Real} end
Fibonacci(d::DataType) = d <: Real ? Fibonacci{d}() : error("Not real type")

Base.iterate(x::Fibonacci{T} where {T<:Real} = (zero(T), (one(T), one(T)))
Base.iterate(x::Fibonacci{T}, state::Tuple{T, T}) where {T<:Real} =
(state[1], (state[2], state[1] + state[2]))

#REPL

julia> for e in Fibonacci(Int64)
e > 100 && break
print(e, "")
end
0 1 1 2 3 5 8 13 21 34 55 89

max_stars_repo_name: jocrv/julia-lang

max_issues_repo_name: jocrv/julia-lang

max_forks_repo_name: jocrv/julia-lang


text:
using MPIMeasurements
using PyPlot

params = Dict{String,Any}()
params["studyName"]="TemporalValidation"
params["studyDescription"]="HarmonicValidation over time"
params["scannerOperator"]="FloThieben"
params["dfStrength"]=[20e-3]

scanner = MPIScanner("MPS.toml")
daq = getDAQ(scanner)

delay = 10.0
numRepetitions = 100

# This version does not store the data
#u = measurement(daq, params, controlPhase=true)

# This version does store the data in a custom location
#filename = measurementRepeatability(daq,"/home/labuser/TemporalValidation.mdf", numRepetitions,
#                                     delay, params, controlPhase=true)
filename = "/home/labuser/TemporalValidation.mdf"

# This version does store the data in the MDFStore
#filename = measurementRepeatability(daq, MDFStore, params,  controlPhase=true)
f = MPIFile(filename)

u = getMeasurements(f, false, #frames=measFGFrameIdx(f),
            fourierTransform=true, bgCorrection=true,
             tfCorrection=false)

#showDAQData(daq,u)


text:
#=
Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.

37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49

It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?
=#
using Primes

function calc()
  start, step = 1,2
  iterations, currentsquaresize = 3,3
  cprimes, call = 0,1
  while true
    start += step
    for i in start:step:(start+(iterations*step))
      if isprime(i)
        cprimes = cprimes + 1
      end
      call = call + 1
    end
    ratio = cprimes / call
    if ratio < 0.1
      return currentsquaresize
    end
    start += (iterations*step)
    currentsquaresize += 2
    step += 2
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf

using Atomistic
using InteratomicPotentials
using Molly
using Plots

N = 864
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"
coupling_factor = 10 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_system = generate_atoms_in_cubic_cell(N, :Ar, box_size, reference_temp)

eq_steps = 2000
eq_thermostat = Molly.AndersenThermostat(reference_temp, Δt * coupling_factor)
eq_simulator = MollySimulator(Δt, eq_steps, coupling=eq_thermostat)
potential = InteratomicPotentials.LennardJones(1.657e-21u"J", 0.34u"nm", 0.765u"nm", [:Ar])

eq_result = @time simulate(initial_system, eq_simulator, potential)

temp = plot_temperature(eq_result, 10)
energy = plot_energy(eq_result, 10)

prod_steps = 5000
prod_simulator = MollySimulator(Δt, prod_steps, t₀=get_time(eq_result))

prod_result = @time simulate(get_system(eq_result), prod_simulator, potential)

display(plot_temperature!(temp, prod_result, 10))
display(plot_energy!(energy, prod_result, 10))

rdf = plot_rdf(prod_result, potential.σ, Int(0.95 * prod_steps))
display(rdf)
savefig(rdf, "artifacts/argon_lj_molly_ip_rdf.svg")

;


text:
#=
We want to support these two syntaxes

a < Number
a < Number("Number")
a < Number("Number", "N")
a < Number("Number", "N", "Parameter N")
a < Number("Number", "N", "Parameter N", 0.0)

a < Number(0.0)
a < Number("Number", 0.0)
a < Number("Number", "N", 0.0)
a < Number("Number", "N", "Parameter N", 0.0)

a < String()
a < String("String")
a < String("String")
a < Number(0.0, "a")

=#

#in_gh(sym) = Expr(:., :GH, QuoteNode(sym))
in_gh(sym) = Symbol("GH$(sym)")
kgh_io_function_names = Symbol[]
is_kgh_io_function_name(sym) =
  sym in kgh_io_function_names

macro ghdef(name, init)
  let str = string(name),
      ghname = esc(in_gh(name))
    quote
      push!(kgh_io_function_names, $(QuoteNode(name)))
      $(ghname)(description::Base.String, short_description=description[1:1], message=description*" parameter", value=$init) =
        [$(str), description, short_description, message, value]
      $(ghname)(value::Base.Any, description=$(str), short_description=description[1:1], message=description*" parameter") =
        [$(str), description, short_description, message, value]
    end
  end
end
@ghdef(String, "")
@ghdef(Path, "")
@ghdef(Boolean, false)
@ghdef(Number, 0.0)
@ghdef(Integer, 0)
@ghdef(Point, u0())
@ghdef(Vector, vx(1))
@ghdef(Any, nothing)
@ghdef(Eval, nothing)
@ghdef(JL, nothing)
@ghdef(Strings, [])
@ghdef(Paths, [])
@ghdef(Booleans, [])
@ghdef(Numbers, [])
@ghdef(Integers, [])
@ghdef(Points, [])
@ghdef(Vectors, [])
@ghdef(Many, [])
@ghdef(Evals, [])
@ghdef(JLs, [])
@ghdef(Stringss, [])
@ghdef(Pathss, [])
@ghdef(Booleanss, [])
@ghdef(Numberss, [])
@ghdef(Integerss, [])
@ghdef(Pointss, [])
@ghdef(Vectorss, [])
@ghdef(Manies, [])
@ghdef(Evalss, [])
@ghdef(JLss, [])

export define_kgh_function

kgh_forms(text, idx=1) =
  let (expr, idx) = Meta.parse(text, idx, greedy=true, depwarn=false)
    isnothing(expr) ?
      [] :
      [expr, kgh_forms(text, idx)...]
  end

is_kgh_io_function_call(e) =
  e isa Expr && e.head === :call && is_kgh_io_function_name(e.args[1])

match_expr(e1, e2) =
  e2 === :_ ||
  e1 == e2 ||
  e1 isa Expr && e2 isa Expr &&
  match_expr(e1.head, e2.head) &&
  length(e1.args) == length(e2.args) &&
  all([match_expr(e1, e2) for (e1, e2) in zip(e1.args, e2.args)])

is_kgh_input(form) =
  match_expr(form, :(_ < _)) &&
  is_kgh_io_function_call(form.args[3])
is_kgh_output(form) =
  (match_expr(form, :(_ > _)) &&
   is_kgh_io_function_call(form.args[3])) ||
  is_kgh_io_function_call(form)

kgh_io_param(form) =
  form.args[2] == :_ ? :__result : form.args[2]
kgh_io_call(form) =
  let param = form.args[2],
      form = form.args[3],
      func = Expr(:., :Khepri, QuoteNode(in_gh(form.args[1])))
    match_expr(form, :(_())) ||
    (match_expr(form, :(_(_))) && !(form.args[2] isa String)) ?
        :($(func)($(form.args[2:end]...), $(string(param)))) :
        :($(func)($(form.args[2:end]...)))
  end

kgh_inputs(forms) = filter(is_kgh_input, forms)
kgh_outputs(forms) = filter(is_kgh_output, forms)

create_kgh_function(name::String, body::String) =
  let forms = kgh_forms(body),
      inputs = filter(is_kgh_input, forms),
      outputs = filter(is_kgh_output, forms),
      forms = filter(f -> !(f in inputs || f in outputs), forms),
      inp_params = map(kgh_io_param, inputs),
      out_params = map(kgh_io_param, outputs),
      inp_forms = map(kgh_io_call, inputs),
      out_forms = map(kgh_io_call, outputs),
      inps = Symbol("__inps_$name"),
      inp_inits = [:($(inp) = $(inps)[$(i)]) for (i, inp) in enumerate(inp_params)],
      outs = Symbol("__outs_$name"),
      out_inits = [:($(outs)[$(i)] = $(out)) for (i, out) in enumerate(out_params)],
      inp_docs = Symbol("__doc_inps_$name"),
      out_docs = Symbol("__doc_outs_$name"),
      shapes = Symbol("__shapes_$name")
    if isempty(inp_forms) && isempty(out_forms)
      quote
        $(inp_docs) = Any[$(inp_forms...)]
        $(out_docs) = Any[$(out_forms...)]
        $(inps) = Array{Any,1}(undef, $(length(inp_inits)))
        $(outs) = Array{Any,1}(undef, $(length(out_inits)))
        $(shapes) = Shape[]
        function $(Symbol("__func_$name"))()
          nothing
        end
      end
    else
      quote
        $(inp_docs) = Any[$(inp_forms...)]
        $(out_docs) = Any[$(out_forms...)]
        $(inps) = Array{Any,1}(undef, $(length(inp_inits)))
        $(outs) = Array{Any,1}(undef, $(length(out_inits)))
        $(shapes) = Shape[]
        function $(Symbol("__func_$name"))()
          $(inp_inits...)
          __result = begin
                $(forms...)
              end
          $(out_inits...)
          nothing
        end
      end
    end
  end

define_kgh_function(name::String, body::String) =
  Base.eval(Main, create_kgh_function(name, body))
#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")

f(a + b)

c = a + 1

c > Any()
d > String("Bar")
"""

create_kgh_function("foo", fn)
=#

#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")
e > String("Bar")
_ > Number()

e = f(a + b)
g(b*c)
"""

create_kgh_function("foo", fn)

fn2 = """
e = f(a + b)
g(b*c)
"""

create_kgh_function("foo", fn2)
using InteractiveUtils
InteractiveUtils.@code_lowered Khepri.foo()
=#
#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")

f(a + b)

c = a + 1

c > Any()
d > String("Bar")
"""

eval(create_kgh_function("foo", fn))
=#

#=
create_kgh_function("foo", raw"""
path < String()
autos < Numbers()
costs < Numbers()

open(path, "w") do f
  for (a, c) in zip(autos, costs)
    println(f, "$a $c")
  end
end
""")

dump(:("foo $bar"))

=#

#=
str1 = """# v < Type("Input", "I", "Parameter I", default)
a < Number(2)
b < Number(3)
_ > Number()

sqrt(a^2 + b^2)"""

create_kgh_function("zzz", str1)
=#


text:
module JsonCsvConverter

using JSON
using DataFrames
using DataStructures

type Converter
    convert::Function
    _csv_to_json::Function
    _json_to_csv::Function
    _get_file_type::Function
    _extract_row_values::Function
    _resolve_options::Function

    function Converter()
        this = new()

        """
        Convert supplied JSON/CSV file to a CSV/JSON file. Returns TRUE if successful; FALSE otherwise.
        """
        this.convert = function (input_file::String, output_file::String; parse_nested_json=false)
            input_filetype = this._get_file_type(input_file)
            output_filetype = this._get_file_type(output_file)

            if input_filetype == "csv" && output_filetype == "json"
                return this._csv_to_json(input_file, output_file, parse_nested_json)
            elseif input_filetype == "json" && output_filetype == "csv"
                return this._json_to_csv(input_file, output_file)
            else
                println("File extensions of the input and output files must be either 'csv' or 'json', and must be different.")
                return false
            end
        end

        this._json_to_csv = function (input_file::String, output_file::String)
            if !isfile(input_file)
                println("Cannot find file $input_file")
                return false
            end

            data = JSON.parsefile(input_file; dicttype=DataStructures.OrderedDict)

            if length(data) == 0
                println("No data was found in the JSON file")
                return false
            end

            df_data = Dict()
            for row in data
                headers = keys(row)
                for header in headers
                    if get(df_data, header, false) === false
                        df_data[header] = []
                    end

                    entity = row[header]
                    if !isa(entity, Number) && !isa(entity, String)
                        entity = JSON.json(entity)
                    end

                    push!(df_data[header], entity)
                end
            end

            writetable(output_file, DataFrame(df_data))

            true
        end

        this._csv_to_json = function (input_file::String, output_file::String, parse_nested_json::Bool)
            if !isfile(input_file)
                println("Cannot find file $input_file")
                return false
            end

            df = readtable(input_file)
            headers = names(df)

            data = []
            for row in eachrow(df)
                push!(data, this._extract_row_values(headers, row, parse_nested_json))
            end

            open(output_file, "w") do f
                write(f, JSON.json(data))
            end

            true
        end

        this._extract_row_values = function (headers::Array{Symbol,1}, row::DataFrames.DataFrameRow{DataFrames.DataFrame}, parse_nested_json::Bool)
            formatted_row = OrderedDict()
            for header in headers
                entity = row[Symbol(header)]

                if entity === NA
                    entity = ""
                elseif parse_nested_json
                    try
                        parsed = JSON.parse(entity)
                        entity = parsed
                    catch e
                    end
                end

                formatted_row[header] = entity
            end

            formatted_row
        end

        this._get_file_type = function (filename::String)
            fragments = split(filename, ".")
            fragments[length(fragments)]
        end

        this._resolve_options = function (options::Dict{String,Bool})
            defaults = Dict(
                "parse_nested_json" => false,
            )

            for (key, value) in defaults
                if get(options, key, nothing) === nothing
                    options[key] = value
                end
            end

            options
        end

        this
    end

end

end


text:
immutable EnumerableGroupBySimple{T,TKey,TS,SO,ES} <: Enumerable{T}
    source::SO
    elementSelector::ES
end

immutable Grouping{TKey,T} <: AbstractArray{T,1}
    key::TKey
    elements::Array{T,1}
end

import Base.size
size{TKey,T}(A::Grouping{TKey,T}) = size(A.elements)
import Base.getindex
getindex{TKey,T}(A::Grouping{TKey,T},i) = A.elements[i]
import Base.length
length{TKey,T}(A::Grouping{TKey,T}) = length(A.elements)

function group_by{TS}(source::Enumerable{TS}, f_elementSelector::Function, elementSelector::Expr)
    TKey = Base.return_types(f_elementSelector, (TS,))[1]

    SO = typeof(source)

    T = Grouping{TKey,TS}

    return EnumerableGroupBySimple{T,TKey,TS,SO,FunctionWrapper{TKey,Tuple{TS}}}(source,f_elementSelector)
end

# TODO This should be rewritten as a lazy iterator
function start{T,TKey,TS,SO,ES}(iter::EnumerableGroupBySimple{T,TKey,TS,SO,ES})
    result = OrderedDict{TKey,T}()
    for i in iter.source
        key = iter.elementSelector(i)
        if !haskey(result, key)
            result[key] = Grouping(key,Array(TS,0))
        end
        push!(result[key].elements,i)
    end
    return collect(values(result)),1
end

function next{T,TKey,TS,SO,ES}(iter::EnumerableGroupBySimple{T,TKey,TS,SO,ES}, state)
    results = state[1]
    curr_index = state[2]
    return results[curr_index], (results, curr_index+1)
end

function done{T,TKey,TS,SO,ES}(iter::EnumerableGroupBySimple{T,TKey,TS,SO,ES}, state)
    results = state[1]
    curr_index = state[2]
    return curr_index > length(results)
end

immutable EnumerableGroupBy{T,TKey,TR,SO,ES,RS} <: Enumerable{T}
    source::SO
    elementSelector::ES
    resultSelector::RS
end

function group_by{TS}(source::Enumerable{TS}, f_elementSelector::Function, elementSelector::Expr, f_resultSelector::Function, resultSelector::Expr)
    TKey = Base.return_types(f_elementSelector, (TS,))[1]

    SO = typeof(source)

    TR = Base.return_types(f_resultSelector, (TS,))[1]

    T = Grouping{TKey,TR}

    return EnumerableGroupBy{T,TKey,TR,SO,FunctionWrapper{TKey,Tuple{TS}},FunctionWrapper{TR,Tuple{TS}}}(source,f_elementSelector,f_resultSelector)
end

# TODO This should be rewritten as a lazy iterator
function start{T,TKey,TR,SO,ES}(iter::EnumerableGroupBy{T,TKey,TR,SO,ES})
    result = OrderedDict{TKey,T}()
    for i in iter.source
        key = iter.elementSelector(i)
        if !haskey(result, key)
            result[key] = Grouping(key,Array(TR,0))
        end
        push!(result[key].elements,iter.resultSelector(i))
    end
    return collect(values(result)),1
end

function next{T,TKey,TR,SO,ES}(iter::EnumerableGroupBy{T,TKey,TR,SO,ES}, state)
    results = state[1]
    curr_index = state[2]
    return results[curr_index], (results, curr_index+1)
end

function done{T,TKey,TR,SO,ES}(iter::EnumerableGroupBy{T,TKey,TR,SO,ES}, state)
    results = state[1]
    curr_index = state[2]
    return curr_index > length(results)
end


text:
using FractionalCalculus
using Plots

x=collect(0.1:0.1:5)
y=collect(0.1:0.1:5)
z=zeros(50, 50)

for j in range(1, 50, step=1)
    z[j, :] = fracdiff(sin, 0.02*j, y, 0.01, RLDiff_Approx())
end

plot(x, y, z, st=:surface)

savefig("./3dexample.png")