text:
using Revise
using AutomotiveDrivingModels
using AutomotivePOMDPs
using AutoViz
using AutonomousMerging
using Distributions
using Random

action_type(d::DriverModel{A}) where A = A

function AutomotiveDrivingModels.propagate(veh::Entity{VehicleState,D,I}, action::LatLonAccel, roadway::Roadway, ΔT::Float64,
                                           nobackup::Bool,
                                           wraparound::Bool) where {D,I}

    a_lat = action.a_lat
    a_lon = action.a_lon

     v = veh.state.v
     ϕ = veh.state.posF.ϕ
    ds = v*cos(ϕ)
     t = veh.state.posF.t
    dt = v*sin(ϕ)

    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²
    Δt = dt*ΔT + 0.5*a_lat*ΔT²

    ds₂ = ds + a_lon*ΔT
    nobackup ? ds₂ = max(0, ds₂) : nothing
    dt₂ = dt + a_lat*ΔT
    speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
    v₂ = sqrt(dt₂*dt₂ + ds₂*ds₂) # v is the magnitude of the velocity vector
    ϕ₂ = atan(dt₂, ds₂)

    if wraparound && veh.posF.s + Δs > roadway[veh.posF.roadind.tag].curve[end].s
        lane = roadway[veh.posF.roadind.tag]
        ds_to_end =  lane.curve[end].s - veh.posF.s
        ds_beginning = Δs - ds_to_end
        curveind = lane.curve[1]
        roadind = RoadIndex{I,T}(curveind, lane.exits[1].target.tag)
        roadind = move_along(roadind, roadway, ds_beginning)
    end
    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    footpoint = roadway[roadind]
    posG = VecE2{Float64}(footpoint.pos.x,footpoint.pos.y) + polar(t + Δt, footpoint.pos.θ + π/2)

    posG = VecSE2{Float64}(posG.x, posG.y, footpoint.pos.θ + ϕ₂)

    veh = VehicleState(posG, roadway, v₂)
    return veh
end

function wrap_around(roadway::Roadway, veh::Vehicle)
    lane = get_lane(roadway, veh)
    s_end = get_end(lane)
    s = veh.state.posF.s
    if s >= s_end
        posF = Frenet(lane, 0.0)
        vehstate = VehicleState(posF, roadway, veh.state.v)
        return Vehicle(vehstate,veh.def, veh.id)
    end
    return veh
end

function AutomotiveDrivingModels.simulate!(scene::Scene, roadway::Roadway, models::Dict{Int64, D}, dt::Float64) where {A, D<:DriverModel}
    s = deepcopy(scene)
    acts = Vector{action_type(models[1])}(undef, length(s))
    get_actions!(acts, s, roadway, models)
    tick!(s, roadway, acts, dt, true)
    for (i, veh) in enumerate(s)
        s[i] = wrap_around(roadway, s[i])
    end
    return s
end

function AutomotiveDrivingModels.simulate!(s0::Scene, roadway::Roadway, models::Dict{Int64, D}, nsteps::Int64, dt::Float64) where D <: DriverModel
    scenes = Vector{Scene}(undef, nsteps+1)
    scenes[1] = s0
    for t=1:nsteps
        scenes[t+1] = simulate!(scenes[t], roadway, models, dt)
    end
    return scenes
end

function initial_scene(n_vehicles, roadway, model;
                       timestep::Float64 = 0.1,
                       n_lanes::Int64 = 3,
                       burn_in::Int64 = 20,
                       init_velocity_dist = Normal(10.0, 3.0),
                       rng::AbstractRNG = Random.GLOBAL_RNG)
    # populate roadway
    scene = Scene()
    models = Dict{Int64, typeof(model)}()
    for i=1:n_vehicles
        v = rand(rng, init_velocity_dist)
        lane_tag = LaneTag(1, rand(rng, 1:n_lanes))
        lane = roadway[lane_tag]
        s = rand(rng, 0.0:lane.curve[end].s)
        posF = Frenet(lane, s)
        vehstate = VehicleState(posF, roadway, v)
        veh = Vehicle(vehstate, VehicleDef(), i)
        models[i] = deepcopy(model)
        set_desired_speed!(models[i], v)
        push!(scene, veh)
    end

    scenes = simulate!(scene, roadway, models, burn_in, timestep)
    return scenes[end], models, scenes
end

##### SIMULATION PARAMETERS #########

roadway = gen_straight_roadway(3, 200.0)

time_step = 0.1

n_vehicles = 20

idm = IntelligentDriverModel(v_des=10.0,
                               d_cmf = 2.0,
                                d_max=2.0,
                                T = 1.5,
                                s_min = 2.0,
                                a_max = 2.0)
model = Tim2DDriver(time_step,
                    mlon = idm,
                    mlat = ProportionalLaneTracker(),
                    mlane = MOBIL(time_step, mlon=idm, politeness=0.35, safe_decel = 2.0))

s0, models, _ = initial_scene(n_vehicles, roadway, model, burn_in=100)

n_steps = 300

@time scenes = simulate!(s0, roadway, models, n_steps, 0.1)

using Blink
using Interact

w = Window()
ui = @manipulate for s=1:n_steps
    render([roadway, scenes[s]])
end
body!(w, ui)


text:
using SeisML
using Test

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

ext: jl


text:
module Evidence

using Nullables
using DataFrames
using CSV

function getEvidence(evidence, idMapping, runDir)
    evidenceMapFile = "$runDir/evidenceIDMapping.tsv"
    sampleNodeValue = Dict()
    if haskey(evidence, "dna")
        evidenceDNA = evidence["dna"]
        if haskey(evidenceDNA, "genomic")
            genomicFile = evidenceDNA["genomic"]
            genomicEvidence = getGenomic(genomicFile, idMapping)
        end
        sampleNodeValue = genomicEvidence["sampleNodeValue"]
        outputToEvidenceMap("genomic", evidenceMapFile, genomicEvidence["geneNodesMap"])
    end

    return sampleNodeValue
end

function outputToEvidenceMap(evidenceType, file, geneNodesMap)
    outfile = open(file, "w")

    for gene in keys(geneNodesMap)
        for node in geneNodesMap[gene]
            id = node[:Database_Identifier]
            nodeName = node[:Node_Name]
            write(outfile, "$evidenceType\t$gene\t$id\t$nodeName\n")
        end
    end

    close(outfile)
end

function getGenomic(file, idMap)
    headers = split(strip(readline(file)),'\t')
    ncols = length(headers)
    datatypes = vcat(String, [Float64 for i=2:ncols])

    df = CSV.read(file, delim="\t", types=datatypes)
    samples = names(df)
    popfirst!(samples)

    geneNodeMap = Dict()
    sampleNodeValue = Dict()
    geneNodesMap = Dict()
    for row in eachrow(df)
        gene = row[:gene]
        if haskey(idMap, gene) == true
            nodes = idMap[gene]
            geneNodesMap[gene] = nodes;
            geneValue = Dict()
            for sample in samples
                value = row[sample]
                if isa(value, Number)
                    if value != -999
                        if haskey(sampleNodeValue, sample) == false
                            sampleNodeValue[sample] = Dict()
                        end
                        for node in nodes
                            sampleNodeValue[sample][node[:Database_Identifier]] = (value < 0.01) ? 0.01 : value
                        end
                    end
                else
                   println("At least one of the values is not a number. Found non number for sample $sample row $gene Use -999 for unknown")
                   exit(1)
                end
            end
        end
    end

    return Dict("sampleNodeValue" => sampleNodeValue,
                "geneNodesMap" => geneNodesMap)
end

function getDnaEvidence(evidence, idMapping)
    divideby = copynumber ? 2 : 1
    data = readlines(fname)

    header = popfirst!(data)
    headerparts = split(chomp(header), "\t")

    samplenodestate = Dict()

    i = 0
    for column in headerparts
        i = i + 1
        if i > 1
            samplenodestate[column] = Dict()
        end
    end

    gene = ""
    for line in data
        lineparts = split(chomp(line), "\t")
        i = 0
        for column in headerparts
            i = i + 1
            if i == 1
                gene = lineparts[1]
            else
                if haskey(genetonodes, gene)
                    for node in Array(genetonodes[gene])
                        samplenodestate[column][node] = oneNormal ?
                                                            parse(Float64, lineparts[i]) :
                                                            copynumber ?
                                                               parse(Float64,lineparts[i]) / 2 :
                                                               parse(Float64,lineparts[i]) - 1
                    end
                end
            end
        end
    end

    return Dict("columns" => headerparts, "samplenodestate" => samplenodestate)
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
function inverse_rotation_matrix(q0::Float64, q1::Float64, q2::Float64, q3::Float64)

	a11::Float64 = q0 * q0 + q1 * q1 - q2 * q2 - q3 * q3
	a12::Float64 = -2.0 * (q0 * q3 - q1 * q2)
	a13::Float64 = 2.0 * (q0 * q2 + q1 * q3)
	a21::Float64 = 2.0 * (q0 * q3 + q1 * q2)
	a22::Float64 = q0 * q0 - q1 * q1 + q2 * q2 - q3 * q3
	a23::Float64 = -2.0 * (q0 * q1 - q2 * q3)
	a31::Float64 = -2.0 * (q0 * q2 - q1 * q3)
	a32::Float64 = 2.0 * (q0 * q1 + q2 * q3)
	a33::Float64 = q0 * q0 - q1 * q1 - q2 * q2 + q3 * q3

	return (a11, a12, a13, a21, a22, a23, a31, a32, a33)

end


text:
#=
In this file we define how agents are plotted and how the plots are updated while stepping
=#
# TODO: I should check whether it is worth to type-parameterize this.
struct ABMStepper # {X, C, M, S, O, AC, AS, AM, HA}
    ac # ::C
    am # ::M
    as # ::S
    offset # ::O
    scheduler # ::X
    pos # ::Observable
    colors # ::AC
    sizes # ::AS
    markers # ::AM
    heatarray # ::HA
    heatobs # ::HO
end

Base.show(io::IO, ::ABMStepper) =
println(io, "Helper structure for stepping and updating the plot of an agent based model. ",
"It is outputted by `abm_plot` and can be used in `Agents.step!`, see `abm_plot`.")

"Initialize the abmstepper and the plotted observables. return the stepper"
function abm_init_stepper_and_plot!(ax, fig, model;
        ac = JULIADYNAMICS_COLORS[1],
        as = 10,
        am = :circle,
        scheduler = model.scheduler,
        offset = nothing,
        aspect = DataAspect(),
        scatterkwargs = NamedTuple(),
        heatarray = nothing,
        heatkwargs = NamedTuple(),
        add_colorbar = true,
        static_preplot! = default_static_preplot,
    )

    heatkwargs = merge((colormap=JULIADYNAMICS_CMAP,), heatkwargs)
    o, e = modellims(model)
    is3d = length(o) == 3
    @assert length(o) == 2 || is3d "Only 2D and 3D spaces can be plotted."
    # TODO: once graph plotting is possible, this will be adjusted
    @assert typeof(model.space) <: Union{Agents.ContinuousSpace, Agents.DiscreteSpace}
    xlims!(ax, o[1], e[1])
    ylims!(ax, o[2], e[2])
    is3d && zlims!(ax, o[3], e[3])
    is3d || (ax.aspect = aspect)

    if !isnothing(heatarray)
        # TODO: This is also possible for continuous spaces, we have to
        # get the matrix size, and then make a range for each dimension
        # and do heatmap!(ax, x, y, heatobs)
        #
        # TODO: use surface!(heatobs) here?
        matrix = Agents.get_data(model, heatarray, identity)
        if !(matrix isa AbstractMatrix) || size(matrix) ≠ size(model.space)
            error("The heat array property must yield a matrix of same size as the grid!")
        end
        heatobs = Observable(matrix)
        hmap = heatmap!(ax, heatobs; heatkwargs...)
    else
        heatobs = nothing
    end
    if add_colorbar && !isnothing(heatobs)
        Colorbar(fig[1, 1][1, 2], hmap, width = 20)
        # rowsize!(fig[1,1].fig.layout, 1, ax.scene.px_area[].widths[2]) # Colorbar height = axis height
    end

    static_preplot!(ax, model)

    if is3d && am == :circle
        am = Sphere(Point3f0(0), 1)
    end
    ids = scheduler(model)
    colors  = ac isa Function ? Observable(to_color.([ac(model[i]) for i ∈ ids])) : to_color(ac)
    sizes   = as isa Function ? Observable([as(model[i]) for i ∈ ids]) : as
    markers = am isa Function ? Observable([am(model[i]) for i ∈ ids]) : am
    postype = is3d ? Point3f0 : Point2f0
    if isnothing(offset)
        pos = Observable(postype[model[i].pos for i ∈ ids])
    else
        pos = Observable(postype[model[i].pos .+ offset(model[i]) for i ∈ ids])
    end

    # Here we make the decision of whether the user has provided markers, and thus use
    # `scatter`, or polygons, and thus use `poly`:
    if is3d
        meshscatter!(
            ax, pos;
            color = colors, marker = markers, markersize = sizes,
            scatterkwargs...
        )
    else
        if user_used_polygons(am, markers)
            # For polygons we always need vector, even if all agents are same polygon
            if markers isa Observable
                markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
            else
                markers = Observable([translate(am, p) for p in pos])
            end
            poly!(ax, markers; color = colors, scatterkwargs...)
        else
            scatter!(
                ax, pos;
                color = colors, marker = markers, markersize = sizes,
                scatterkwargs...
            )
        end
    end
    return ABMStepper(
        ac, am, as, offset, scheduler,
        pos, colors, sizes, markers,
        heatarray, heatobs
    )
end

default_static_preplot(ax, model) = nothing

function modellims(model)
    if model.space isa Agents.ContinuousSpace
        e = model.space.extent
    elseif model.space isa Agents.DiscreteSpace
        e = size(model.space.s) .+ 1
    end
    return zero.(e), e
end

function user_used_polygons(am, markers)
    if (am isa Polygon)
        return true
    elseif (am isa Function) && (markers[][1] isa Polygon)
        return true
    else
        return false
    end
end

#=
    Agents.step!(abmstepper, model, agent_step!, model_step!, n::Int)
Step the given `model` for `n` steps while also updating the plot that corresponds to it,
which is produced with the function [`abm_plot`](@ref).

You can still call this function with `n=0` to update the plot for a new `model`,
without doing any stepping.
=#
function Agents.step!(abmstepper::ABMStepper, model, agent_step!, model_step!, n)
    @assert (n isa Int) "Only stepping with integer `n` is possible with `abmstepper`."
    ac, am, as = abmstepper.ac, abmstepper.am, abmstepper.as
    offset = abmstepper.offset
    pos, colors = abmstepper.pos, abmstepper.colors
    sizes, markers =  abmstepper.sizes, abmstepper.markers

    Agents.step!(model, agent_step!, model_step!, n)

    if Agents.nagents(model) == 0
        @warn "The model has no agents"
    end
    ids = abmstepper.scheduler(model)
    if isnothing(offset)
        pos[] = [model[i].pos for i in ids]
    else
        pos[] = [model[i].pos .+ offset(model[i]) for i in ids]
    end
    if ac isa Function; colors[] = to_color.([ac(model[i]) for i in ids]); end
    if as isa Function; sizes[] = [as(model[i]) for i in ids]; end
    if am isa Function; markers[] = [am(model[i]) for i in ids]; end
    # If we use Polygons as markers, do a final update:
    if user_used_polygons(am, markers)
        # translate all polygons according to pos
        markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
    end
    # Finally update the heat array, if any
    if !isnothing(abmstepper.heatarray)
        newmatrix = Agents.get_data(model, abmstepper.heatarray, identity)
        abmstepper.heatobs[] = newmatrix
    end
    return nothing
end


text:
function conv_bn_block(input, output, kernel_size)
    return Chain(
        Conv((kernel_size,), input=>output),
        BatchNorm(output),
        x -> relu.(x))
end

function fc_bn_block(input, output)
    return Chain(
        Dense(input,output),
        BatchNorm(output),
        x -> relu.(x))
end

function CreateSingleKNNGraph(X::AbstractArray, K::Int)
    F, N = size(X)
    kdtree = KDTree(X)
    cat([X[:,(knn(kdtree, X[:,i], K+1, true)[1][2:K+1])] for i in 1:N]..., dims=3)
end

@nograd CreateSingleKNNGraph

struct EdgeConv
    layers::AbstractArray
    K::Int
    mlp
    maxpool_K
end

function EdgeConv(layers::AbstractArray, K::Int)
    m = []

    for i in 1:(length(layers)-1)
        if i == 1
            push!(m, conv_bn_block(2*layers[i], layers[i+1], 1))
        else
            push!(m, conv_bn_block(layers[i], layers[i+1], 1))
        end
    end

    EdgeConv(layers, K, Chain(m...), MaxPool((K,)))
end

function (m::EdgeConv)(X)
    F, N, B = size(X)
    # X: [F, N, B]

    KNNGraph = cat([CreateSingleKNNGraph(X[:,:,i], m.K) for i in 1:B]..., dims=4)
    # KNNGraph: [F, K, N, B]

    X = reshape(X, F, 1, N, B)
    # X: [F, 1, N, B]

    X = cat([X for i in 1:m.K]...,dims=2)
    # X: [F, K, N, B]

    X = cat(X, KNNGraph-X, dims=1)
    # X: [2*F, K, N, B]

    X = PermutedDimsArray(X, (2,3,1,4))
    # X: [K, N, 2*F, B]

    X = reshape(X, N*m.K, 2*F,B)
    # X: [K*N, 2*F, B]

    X = m.mlp(X)
    # X: [K*N, an, B]

    an = size(X)[2]
    X = reshape(X, m.K, an*N, B)
    # X: [K, N*an, B]

    X = m.maxpool_K(X)
    # X: [1, N*an, B]

    X = reshape(X, N, an, B)
    # X: [N, an, B]

    X = PermutedDimsArray(X, (2,1,3))
    # X: [an, N, B]

    return X
end

@functor EdgeConv

struct DGCNN_cls
    EdgeConv1::EdgeConv
    EdgeConv2::EdgeConv
    conv_3
    maxpool_3
    fc_4
    drop_4
    fc_5
    drop_5
    fc_6
end

function DGCNN_cls(num_classes::Int, K::Int, N::Int)
    DGCNN_cls(
        EdgeConv([3,32,64,64], K),
        EdgeConv([64,128,256], K),
        conv_bn_block(256,1024,1),
        MaxPool((N,)),
        fc_bn_block(1024,512),
        Dropout(0.5),
        fc_bn_block(512,256),
        Dropout(0.5),
        Dense(256,num_classes)
        )
end

function (m::DGCNN_cls)(X)

    # X: [N, 3, B]

    X = PermutedDimsArray(X, (2,1,3))
    C,N,B = size(X)
    # X: [3, N, B]

    X = m.EdgeConv1(X)
    # X: [64, N, B]

    X = m.EdgeConv2(X)
    # X: [128, N, B]

    X = PermutedDimsArray(X, (2,1,3))
    # X: [N, 128, B]

    X = m.conv_3(X)
    # X: [N, 1024, B]

    X = m.maxpool_3(X)
    # X: [1, 1024, B]

    X = reshape(X, :, B)
    # X: [1024, B]

    X = m.drop_4(m.fc_4(X))
    # X: [512, B]

    X = m.drop_5(m.fc_5(X))
    # X: [256, B]

    X = m.fc_6(X)
    # X: [num_classes, B]

    return softmax(X, dims=1)

end

@functor DGCNN_cls

max_stars_repo_name: nirmal-suthar/DGCNN

max_issues_repo_name: nirmal-suthar/DGCNN

max_forks_repo_name: nirmal-suthar/DGCNN


text:
## Example script for constant-strength doublet-source panel method
using AeroMDAO

##
alpha_u = [0.2, 0.3, 0.2, 0.15, 0.2]
alpha_l = [-0.2, -0.1, -0.1, -0.001]
dzs     = (0., 0.)
airfoil = kulfan_CST(alpha_u, alpha_l, dzs, (0., 0.), 60);      # Kulfan CST airfoil
# airfoil = naca4((0,0,1,2)) # NACA 4-digit airfoil
uniform = Uniform2D(angle = 5.0)
system  = @time solve_case(
                     airfoil, uniform;
                     num_panels = 80
                    );

##
panels    = system.surface_panels
@time ues = surface_velocities(system);
@time cl  = lift_coefficient(system)
@time cls, cms, cps = surface_coefficients(system)

## Printing
println("Cl: $cl")
println("Σᵢ Clᵢ: $(sum(cls))")
println("Σᵢ Cmᵢ: $(sum(cms))")

## Plotting libraries
using Plots
gr();

## Pressure coefficients

# Split surfaces and values
cp_upper, cp_lower = get_surface_values(panels, cps)

plot(marker = 2, label = :none, yflip = true, xlabel = "(x/c)", ylabel = "Cp")
plot!(cp_upper, label = "Upper")
plot!(cp_lower, label = "Lower")

## Lift distribution
cl_upper, cl_lower = get_surface_values(panels, cls)

plot(marker = 2, label = :none, yflip = true, xlabel = "(x/c)", ylabel = "Cl")
plot!(cl_upper, label = "Upper")
plot!(cl_lower, label = "Lower")

## Angle of attack sweep example
function alpha_sweep(α, airfoil)
    uniform = Uniform2D(1.0, α)
    lift_coefficient(solve_case(airfoil, uniform, num_panels = 80))
end

αs = 0:10
cls = @time alpha_sweep.(αs, Ref(airfoil))

##
plot(αs, cls, ylabel = "Cl", xlabel = "α (deg)", label = :none)


text: truncated(d::LogUniform, lo::T, hi::T) where {T<:Real} = LogUniform(max(d.a, lo), min(d.b, hi))


text:
using Base.Test
using DataFrames

let
	dv1 = DataVector[9,1,8,NA,3,3,7,NA]
	dv2 = 1.0 * dv1
	dv3 = DataArray([1:8])
	pdv1 = PooledDataArray(dv1)
	idv1 = IndexedVector(dv1)

	@assert sortperm(dv1) == sortperm(dv2)
	@assert sortperm(dv1) == sortperm(pdv1)
	@assert sortperm(dv1) == sortperm(idv1)
	@assert isequal(sort(dv1), DataArray(sort(dv1)))
	@assert isequal(sort(dv1), DataArray(sort(pdv1)))
	@assert isequal(sort(dv1), DataArray(sort(idv1)))

	d = @DataFrame(dv1 => dv1, dv2 => dv2, dv3 => dv3, pdv1 => pdv1, idv1 => idv1)

	@assert sortperm(d) == sortperm(dv1)
	@assert sortperm(d[["dv3","dv1"]]) == sortperm(dv3)
	@assert sortby(d, "dv1")["dv3"] == sortperm(dv1)
	@assert sortby(d, "dv2")["dv3"] == sortperm(dv1)
	@assert sortby(d, "pdv1")["dv3"] == sortperm(dv1)
	@assert sortby(d, "idv1")["dv3"] == sortperm(dv1)
	@assert sortby(d, ["dv1","pdv1"])["dv3"] == sortperm(dv1)
	@assert sortby(d, ["dv1","dv3"])["dv3"] == sortperm(dv1)
end


text:
using SortingLab;
using BenchmarkTools, DataFrames, CSV, Plots, StatPlots
using RCall, IterTools, uCSV

function string_sort_and_perm_vs_r(svec)
    sort_id = @benchmark sort($svec)
    sortperm_id = @benchmark sortperm($svec)

    radixsort_id = @benchmark radixsort($svec)
    fsortperm_id = @benchmark fsortperm($svec)

    @rput svec;

    r_timing = R"""
    memory.limit(40000)
    res = NULL
    for(i in 1:$(length(radixsort_id.times))) {
        pt = proc.time()
        system.time(sort(svec, method="radix"))
        pt2 = proc.time()
        gc()
        res = c(res, (pt2 - pt)[3])
    }

    resorder = NULL
    for(i in 1:$(length(fsortperm_id.times))) {
        pt = proc.time()
        system.time(order(svec, method="radix"))
        pt2 = proc.time()
        gc()
        resorder = c(resorder, (pt2 - pt)[3])
    }
    rm(svec);gc()
    list(res, resorder)
    """
    sort_id, sortperm_id, radixsort_id, fsortperm_id, r_timing
end

sortperf2timings(x) = vcat([mean(x[i].times)/1e9 for i=1:4], [mean(x[5][1]), mean(x[5][2])])

function string_sort_perf(N, K, randomseed = 1)
    tic()
    srand(randomseed);
    svec = rand("id".*dec.(1:N÷K, 10), N);
    @time id_timings = string_sort_and_perm_vs_r(svec);

    srand(randomseed);
    fstr_samplespace = [randstring(8) for i = 1:N÷K];
    fsvec = rand(fstr_samplespace, N);
    @time fstr_timings = string_sort_and_perm_vs_r(fsvec)

    srand(randomseed);
    rstr_samplespace = [randstring(rand(1:8)) for i = 1:N÷K];
    rsvec = rand(rstr_samplespace, N);
    @time rstr_timings = string_sort_and_perm_vs_r(rsvec);
    toc()

    timings = reduce(vcat, sortperf2timings.([id_timings, fstr_timings, rstr_timings]))
    labels1 = ["IDs", "fixed len 8", "var len 8"]
    labels2 = ["sort", "sortperm"]
    labels3 = ["Julia Base", "Julia SortingLab.jl", "R"]

    group = repeat(labels3,  inner = 2, outer = 3)
    xlabel = repeat(labels1, inner = 6)
    sort_or_perm = repeat(labels2, outer = 9)

    df  = DataFrame(timings = timings, group = group, xlabel = xlabel, sort_or_perm = sort_or_perm, N = repeat([N], inner = 18), K = repeat([K], inner = 18))
    fileout = "string_sort_perf_$(N÷1_000_000)m.csv"
    CSV.write(fileout, df)
end

function read_string_sort_perf_plot(N)
    df = CSV.read("string_sort_perf_$(N÷1_000_000)m.csv")

    df1 = df[(df[:sort_or_perm] .== "sort") .& (df[:group] .!= "R"),:]
    gb1 = groupedbar(
        df1[:xlabel],
        df1[:timings],
        group = df1[:group],
        ylabel = "seconds",
        title = "Base.sort vs SortingLab.radixsort ($(N÷1_000_000)m)")
    savefig(gb1, "benchmarks/sort_vs_radixsort_$(N÷1_000_000)m.png")

    df2 = df[(df[:sort_or_perm] .== "sort") .& (df[:group] .!= "R"),:]
    gb2 = groupedbar(
        df2[:xlabel],
        df2[:timings],
        group = df2[:group],
        ylabel = "seconds",
        title = "Base.sortperm vs SortingLab.fsortperm ($(N÷1_000_000)m)")
    savefig(gb2, "benchmarks/sortperm_vs_fsorperm_$(N÷1_000_000)m.png")
end


text:
# Ported from https://github.com/openai/gym/blob/996e5115621bf57b34b9b79941e629a36a709ea1/gym/envs/classic_control/mountain_car.py
# which has header
#  https://webdocs.cs.ualberta.ca/~sutton/MountainCar/MountainCar1.cp

module MountainCarEnv

using Reinforce: AbstractEnvironment
using LearnBase: DiscreteSet
using RecipesBase
using Distributions
@static if VERSION >= v"0.7"
  using Random: seed!
else
  seed! = srand
end

import Reinforce: reset!, actions, finished, step!

export
  MountainCar,
  reset!,
  step!,
  actions,
  finished,
  f

const min_position = -1.2
const max_position = 0.6
const max_speed = 0.07
const goal_position = 0.5
const min_start = -0.6
const max_start = 0.4

const car_width = 0.05
const car_height = car_width/2.0
const clearance = 0.2*car_height
const flag_height = 0.05

mutable struct MountainCarState
  position::Float64
  velocity::Float64
end

mutable struct MountainCar <: AbstractEnvironment
  state::MountainCarState
  reward::Float64
  seed::Int
end
MountainCar(seed=-1) = MountainCar(MountainCarState(0.0, 0.0), 0.0, seed)

function reset!(env::MountainCar)
  if env.seed >= 0
    seed!(env.seed)
    env.seed = -1
  end

  env.state.position = rand(Uniform(min_start, max_start))
  env.state.velocity = 0.0

  env
end

actions(env::MountainCar, s) = DiscreteSet(1:3)
finished(env::MountainCar, s′) = env.state.position >= goal_position

function step!(env::MountainCar, s::MountainCarState, a::Int)
  position = env.state.position
  velocity = env.state.velocity
  velocity += (a - 2)*0.001 + cos(3*position)*(-0.0025)
  velocity = clamp(velocity, -max_speed, max_speed)
  position += velocity
  if position <= min_position && velocity < 0
    velocity = 0
  end
  position = clamp(position, min_position, max_position)
  env.state = MountainCarState(position, velocity)
  env.reward = -1

  return env.reward, env.state
end

# ------------------------------------------------------------------------
height(xs) = sin(3 * xs)*0.45 + 0.55
rotate(xs::Array{Float64}, ys::Array{Float64}, Θ::Float64) =
  xs.*cos(Θ) .- ys.*sin(Θ), ys.*cos(Θ) .+ xs.*sin(Θ)

translate(xs::Array{Float64}, ys::Array{Float64}, t::Array{Float64}) =
  xs .+ t[1], ys .+ t[2]

@recipe function f(env::MountainCar)
  legend := false
  xlims := (min_position, max_position)
  ylims := (0, 1.1)
  grid := false
  ticks := nothing

  # Mountain
  @series begin
    xs = range(min_position, stop = max_position, length = 100)
    ys = height.(xs)
    seriestype := :path
    linecolor --> :blue
    xs, ys
  end

  # Car
  @series begin
    fillcolor --> :black
    seriestype := :shape

    θ = cos(3 * env.state.position)
    xs = [-car_width/2, -car_width/2, car_width/2, car_width/2]
    ys = [0, car_height, car_height, 0] .+ clearance
    xs, ys = rotate(xs, ys, θ)
    translate(xs, ys, [env.state.position, height(env.state.position)])
  end

  # Flag
  @series begin
    linecolor --> :red
    seriestype := :path

    [goal_position, goal_position], [height(goal_position), height(goal_position) + flag_height]
  end
end
end


text:
module TestLearningCurves

# using Revise
using Test
using MLJ
import Random.seed!
seed!(1234)

@load KNNRegressor

include("foobarmodel.jl")

x1 = rand(100);
x2 = rand(100);
x3 = rand(100);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

@testset "learning curves" begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=50, rng=1)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda),
                     lower=0.0001, upper=0.1, scale=:log10)
    curve = @test_logs((:info, r"Using measure=rms"),
                       (:info, r"Training"),
                       (:info, r"Training of best"),
                       MLJ.learning_curve!(mach; range=r_lambda))
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curve2 = MLJ.learning_curve!(mach; range=r_n)
    curve3 = learning_curve(ensemble, X, y; range=r_n)
    @test curve2.measurements ≈ curve3.measurements
end

end # module
true


text:
struct AType
  other_obj::SType
end


text:
using wgregseq, FASTX, DataFrames, CSV, BioSequences, CairoMakie,  StatsBase

# Set path
dir = @__DIR__
home_dir = joinpath(split(dir, "/")[1:end-2])

# Import promoter list and infer types that can be infered automatically

# Promoters from Ecocyc
promoter_list_ecocyc = CSV.read(
    "/$home_dir/data/promoter_list_ecocyc.csv",
    DataFrames.DataFrame,
    types=Dict(
        "promoter"=>String,
        "tss"=>Float64,
        "direction"=>String
    )
)
promoter_list_ecocyc
# Promoters from RegulonDB
promoter_list_regulonDB = CSV.read(
    "/$home_dir/data/promoter_list_regulon_DB.csv",
    DataFrames.DataFrame,
    types=Dict(
        "promoter"=>String,
        "tss"=>Float64,
        "direction"=>String
    )
)

# Replace columns by nicer types
promoter_list_ecocyc.genes = parse.(Vector{String}, promoter_list_ecocyc.genes)
promoter_list_ecocyc.evidence = parse.(Vector{String}, promoter_list_ecocyc.evidence)
promoter_list_ecocyc.gene_position = parse.(Vector{Float64}, promoter_list_ecocyc.gene_position)
promoter_list_regulonDB.genes = parse.(Vector{String}, promoter_list_regulonDB.genes)
promoter_list_regulonDB.gene_position = parse.(Vector{Float64}, promoter_list_regulonDB.gene_position)

# Join the datasets
df = promoter_list_ecocyc[(.~ isnan.(promoter_list_ecocyc.tss)) .& (map(x-> x != ["None"], promoter_list_ecocyc.genes)), :]
promoter_list_regulonDB.tss = coalesce.(promoter_list_regulonDB.tss, NaN)
#append!(df, promoter_list_regulonDB)

##

# Drop missing entries
dropmissing!(df)
df.direction = convert(Vector{String}, df.direction)

# Sort gene entries to remove duplicates
for i in 1:nrow(df)
    genes, positions = df[i, ["genes", "gene_position"]]
    index = sortperm(genes)
    df[i, ["genes", "gene_position"]] = genes[index], positions[index]
end

# Alphabetically sort dataframe and remove duplicates
sort!(df, "genes")
df = unique(df)

# Last 5 entries are artifacts
df = df[5:end, :]

# Store list
CSV.write("/$home_dir/data/promoter_list_processed.csv", df)

##
df_genes = CSV.read("/$home_dir/data/all_genes_table.csv", DataFrames.DataFrame)[1:end-1, :]
# Import promoters with genes annotated
df = df[map(x -> length(x)>0, df.genes), :]

# Find genes which have a promoter annotated and which don't
gdf = groupby(df, "promoter")
gene_promoter = []
for group in gdf
    push!(gene_promoter, group.genes...)
end
genes_w_promoter = unique(vcat(gene_promoter...))
genes_wo_promoter = filter(x -> x ∉ genes_w_promoter, df_genes.gene)
println("Percentage of genes without a promoter: $(length(genes_wo_promoter)/(length(genes_wo_promoter)+length(genes_w_promoter)))")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -0.9705716506753843 <= q <= 1.6199082189427347)

                     add_NL_constraint(m, :((0.27812717705086953 + -0.4967460082986306*(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))) + -0.883304587160544*(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))))/(1 + exp(-(0.27812717705086953 + -0.4967460082986306*(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))) + -0.883304587160544*(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))))))) + (0.48450040123340576 + -0.32897580116872405*(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))) + 0.43090905819531766*(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))))/(1 + exp(-(0.48450040123340576 + -0.32897580116872405*(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(-0.6740257057257448 + -0.0007917455797805673*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.5330540695702357*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))))))) + 0.43090905819531766*(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))/(1 + exp(-(0.3709718635886192 + -0.47514710764751555*(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.08584163188401472 + 0.6683752749724832*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + -0.39230803394698555*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))) + 0.6627028336583534*(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))))/(1 + exp(-(-0.7769463282105753 + 0.25325339074048747*(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2]))/(1 + exp(-(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])))) + 0.5840452662855768*(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))/(1 + exp(-(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Parameters
include("common.jl")

# Na-K ATPase (NKA) parameters
@with_kw struct NKAParams
    k_o = 5.4
    na_o = 140.0
    # Max Na-K ATPase current (μA/cm²)
    IMAX = 3.147
    # Na half-saturate constant of Na-K ATPase (mM)
    KM_NA = 10.0
    # K half-saturate constant of Na-K ATPase (mM)
    KM_K = 1.5
    # ATP half-saturate constant of Na-K ATPase (mM)
    KM_ATP = 8E-3
    # ADP half-saturate constant of Na-K ATPase (mM)
    KI_ADP = 0.1
    # Factor of extracellular potassium
    FAC_K = _mm(k_o, KM_K)
    # Factor of extracellular sodium
    FAC_NA = expm1(na_o / 67.3) / 7
end

# Na/K ATPase current (iNaK) (uA/cm^2), scalar version
function inak(na_i, atp_i, adp_i, vm, evfrtm1, KM_ATP, KI_ADP, KM_NA, IMAX, FAC_K, FAC_NA)
    evfrt = evfrtm1 + 1
    f_na = _hills(na_i, KM_NA, 1.5)
    f_nak_inv = 1.0 + 0.1245 * _ra(-0.1 * vm) + 0.0365 * FAC_NA / evfrt
    f_atp = _mm(atp_i * _mm(KI_ADP, adp_i), KM_ATP)
    iNak = IMAX * FAC_K * f_na * f_atp / f_nak_inv
end

function inak(na_i, atp_i, adp_i, vm, evfrtm1, pNKA::NKAParams)
    @unpack KM_ATP, KI_ADP, KM_NA, IMAX, FAC_K, FAC_NA = pNKA
    inak(na_i, atp_i, adp_i, vm, evfrtm1, KM_ATP, KI_ADP, KM_NA, IMAX, FAC_K, FAC_NA)
end

# Plasma membrane calcium ATPase (PMCA) parameters
@with_kw struct PMCAParams
    # Max PMCA current (μA/cm²)
    IMAX = 0.575
    # Ca half-saturation constant (mM)
    KM_CA = 5E-4
    # ATP 1st half-saturation constant (mM)
    KM1_ATP = 0.012
    # ATP 2nd half-saturation constant (mM)
    KM2_ATP = 0.23
    # ADP half-inhibition constant (mM)
    KI_ADP = 1.0
end

# PMCA current (iPCa) (uA/cm^2), scalar version
function ipca(ca_i, atp_i, adp_i, KM1_ATP, KI_ADP, KM2_ATP, KM_CA, IMAX)
    f_atp = _mm(atp_i * _mm(KI_ADP, adp_i), KM1_ATP) + _mm(atp_i, KM2_ATP)
    f_ca = _mm(ca_i, KM_CA)
    iPca = IMAX * f_atp * f_ca
end

function ipca(ca_i, atp_i, adp_i, pPMCA::PMCAParams)
    @unpack KM1_ATP, KI_ADP, KM2_ATP, KM_CA, IMAX = pPMCA
    ipca(ca_i, atp_i, adp_i, KM1_ATP, KI_ADP, KM2_ATP, KM_CA, IMAX)
end

# SER ca ATPase (SERCA) parameters
# Params from Li et al. (2015)
@with_kw struct SERCAParams
    # Max forward rate (mM/ms)
    VF = 2.989E-4
    # Max reverse rate (mM/ms)
    VR = 3.179E-4
    # Michaelis constant for Ca of forward reaction (mM)
    KMF_CA = 1.5E-4
    # Michaelis constant for Ca of reverse reaction (mM)
    KMR_CA = 3.3
    # Cooperativity of forward reaction
    NFB = 0.55
    # Cooperativity of reverse reaction
    NRB = 0.5
    # Michaelis constant for ATP (mM)
    KM_ATP = 0.01
    # 1st michaelis constant for ADP (mM)
    KI1_ADP = 0.14
    # 2nd michaelis constant for ADP (mM)
    KI2_ADP = 5.1
end

# SERCA calcium flux (jUp) (mM/ms), scalar version
function jup(ca_i, ca_nsr, atp_i, adp_i,
             NFB, NRB, KMF_CA, KMR_CA, KM_ATP, KI1_ADP, KI2_ADP, VF, VR)
    fb = (ca_i / KMF_CA) ^ NFB
    rb = (ca_nsr / KMR_CA) ^ NRB
    f_atp = inv(KM_ATP / atp_i * _mm_reci(KI1_ADP, adp_i) + _mm_reci(KI2_ADP, adp_i))
    jUp = f_atp * (VF * fb - VR * rb) / (1.0 + fb + rb)
end

function jup(ca_i, ca_nsr, atp_i, adp_i, pSERCA::SERCAParams)
    @unpack NFB, NRB, KMF_CA, KMR_CA, KM_ATP, KI1_ADP, KI2_ADP, VF, VR = pSERCA
    jup(ca_i, ca_nsr, atp_i, adp_i,
        NFB, NRB, KMF_CA, KMR_CA, KM_ATP, KI1_ADP, KI2_ADP, VF, VR)
end


text:
using RDatasets
using Random
using Statistics
using LinearAlgebra
using Plots

iris = dataset("datasets", "iris")

X = Matrix(iris[:, 1:4])
y = iris.Species

# Exercise

# Exercise

# Exercise

# Preparing data

function prepare_data(X, y; do_normal=true, do_onehot=true, kwargs...)
    X_train, y_train, X_test, y_test = split(X, y; kwargs...)

    if do_normal
        X_train, X_test = normalize(X_train, X_test; kwargs...)
    end

    classes = unique(y)

    if do_onehot
        y_train = onehot(y_train, classes)
        y_test = onehot(y_test, classes)
    end

    return X_train, y_train, X_test, y_test, classes
end

Random.seed!(666)

iris = dataset("datasets", "iris")

X = Matrix(iris[:, 1:4])
y = iris.Species

X_train, y_train, X_test, y_test, classes = prepare_data(X', y; dims=2)

# Testing dims in prepare_data

Random.seed!(666)
aux1 = prepare_data(X, y; dims=1)
Random.seed!(666)
aux2 = prepare_data(X', y; dims=2)

norm(aux1[1] - aux2[1]')

# SimpleNet

struct SimpleNet{T<:Real}
    W1::Matrix{T}
    b1::Vector{T}
    W2::Matrix{T}
    b2::Vector{T}
end

# Exercise

# Initialize SimpleNet

Random.seed!(666)

m = SimpleNet(size(X_train,1), 5, size(y_train,1))

# Exercise

# Evaluate SimpleNet

m(X_train[:,1:2])

# Gradient computation

function grad(m::SimpleNet, x::AbstractVector, y; ϵ=1e-10)
    z1 = m.W1*x .+ m.b1
    a1 = max.(z1, 0)
    z2 = m.W2*a1 .+ m.b2
    a2 = exp.(z2) ./ sum(exp.(z2), dims=1)
    l = -sum(y .* log.(a2 .+ ϵ))

    e_z2 = exp.(z2)
    l_part = (- e_z2 * e_z2' + Diagonal(e_z2 .* sum(e_z2))) / sum(e_z2)^2

    l_a2 = - y ./ (a2 .+ ϵ)
    l_z2 = l_part * l_a2
    l_a1 = m.W2' * l_z2
    l_z1 = l_a1 .* (a1 .> 0)
    l_x = m.W1' * l_z1

    l_W2 = l_z2 * a1'
    l_b2 = l_z2
    l_W1 = l_z1 * x'
    l_b1 = l_z1

    return l, l_W1, l_b1, l_W2, l_b2
end

# Few comments about the gradient

g_all = [grad(m, X_train[:,k], y_train[:,k]) for k in 1:size(X_train,2)]

typeof(g_all)

mean_tuple(d::AbstractArray{<:Tuple}) = Tuple([mean([d[k][i] for k in 1:length(d)]) for i in 1:length(d[1])])

g_mean = mean_tuple(g_all)

typeof(g_mean)

# Exercise

# Exercise

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
module operators_sparse

export SparseOperator, diagonaloperator

import Base: ==, *, /, +, -
import ..operators

using ..bases, ..states, ..operators, ..operators_dense, ..sparsematrix

"""
    SparseOperator(b1[, b2, data])

Sparse array implementation of Operator.

The matrix is stored as the julia built-in type `SparseMatrixCSC`
in the `data` field.
"""
mutable struct SparseOperator <: Operator
    basis_l::Basis
    basis_r::Basis
    data::SparseMatrixCSC{Complex128, Int}
    function SparseOperator(b1::Basis, b2::Basis, data)
        if length(b1) != size(data, 1) || length(b2) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(b1, b2, data)
    end
end

SparseOperator(b::Basis, data::SparseMatrixCSC{Complex128, Int}) = SparseOperator(b, b, data)
SparseOperator(b::Basis, data::Matrix{Complex128}) = SparseOperator(b, sparse(data))
SparseOperator(op::DenseOperator) = SparseOperator(op.basis_l, op.basis_r, sparse(op.data))

SparseOperator(b1::Basis, b2::Basis) = SparseOperator(b1, b2, spzeros(Complex128, length(b1), length(b2)))
SparseOperator(b::Basis) = SparseOperator(b, b)

Base.copy(x::SparseOperator) = SparseOperator(x.basis_l, x.basis_r, copy(x.data))
Base.full(a::SparseOperator) = DenseOperator(a.basis_l, a.basis_r, full(a.data))

"""
    sparse(op::Operator)

Convert an arbitrary operator into a [`SparseOperator`](@ref).
"""
Base.sparse(a::Operator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead."))
Base.sparse(a::SparseOperator) = copy(a)
Base.sparse(a::DenseOperator) = SparseOperator(a.basis_l, a.basis_r, sparse(a.data))

==(x::SparseOperator, y::SparseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && (x.data == y.data)

# Arithmetic operations
+(a::SparseOperator, b::SparseOperator) = (check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data+b.data))
+(a::SparseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data+b.data))
+(a::DenseOperator, b::SparseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data+b.data))

-(a::SparseOperator) = SparseOperator(a.basis_l, a.basis_r, -a.data)
-(a::SparseOperator, b::SparseOperator) = (check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data-b.data))
-(a::SparseOperator, b::DenseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data-b.data))
-(a::DenseOperator, b::SparseOperator) = (check_samebases(a,b); DenseOperator(a.basis_l, a.basis_r, a.data-b.data))

*(a::SparseOperator, b::SparseOperator) = (check_multiplicable(a, b); SparseOperator(a.basis_l, b.basis_r, a.data*b.data))
*(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, complex(b)*a.data)
*(a::Number, b::SparseOperator) = SparseOperator(b.basis_l, b.basis_r, complex(a)*b.data)

/(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, a.data/complex(b))

operators.dagger(x::SparseOperator) = SparseOperator(x.basis_r, x.basis_l, x.data')
operators.ishermitian(A::SparseOperator) = ishermitian(A.data)

operators.tensor(a::SparseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))
operators.tensor(a::DenseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))
operators.tensor(a::SparseOperator, b::DenseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))

operators.trace(op::SparseOperator) = (check_samebases(op); trace(op.data))

function operators.ptrace(op::SparseOperator, indices::Vector{Int})
    operators.check_ptrace_arguments(op, indices)
    rank = length(op.basis_l.shape) - length(indices)
    shape = [op.basis_l.shape; op.basis_r.shape]
    data = sparsematrix.ptrace(op.data, shape, indices)
    b_l = ptrace(op.basis_l, indices)
    b_r = ptrace(op.basis_r, indices)
    SparseOperator(b_l, b_r, data)
end

function operators.expect(op::SparseOperator, state::DenseOperator)
    check_samebases(op.basis_r, state.basis_l)
    check_samebases(op.basis_l, state.basis_r)
    result = Complex128(0.)
    @inbounds for colindex = 1:op.data.n
        for i=op.data.colptr[colindex]:op.data.colptr[colindex+1]-1
            result += op.data.nzval[i]*state.data[colindex, op.data.rowval[i]]
        end
    end
    result
end

function operators.permutesystems(rho::SparseOperator, perm::Vector{Int})
    @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
    @assert isperm(perm)
    shape = [rho.basis_l.shape; rho.basis_r.shape]
    data = sparsematrix.permutedims(rho.data, shape, [perm; perm + length(perm)])
    SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
end

operators.identityoperator(::Type{SparseOperator}, b1::Basis, b2::Basis) = SparseOperator(b1, b2, speye(Complex128, length(b1), length(b2)))
operators.identityoperator(b1::Basis, b2::Basis) = identityoperator(SparseOperator, b1, b2)
operators.identityoperator(b::Basis) = identityoperator(b, b)

"""
    diagonaloperator(b::Basis)

Create a diagonal operator of type [`SparseOperator`](@ref).
"""
function diagonaloperator(b::Basis, diag::Vector{T}) where T <: Number
  @assert 1 <= length(diag) <= prod(b.shape)
  SparseOperator(b, spdiagm(convert(Vector{Complex128}, diag)))
end

# Fast in-place multiplication implementations
operators.gemm!(alpha, M::SparseOperator, b::DenseOperator, beta, result::DenseOperator) = sparsematrix.gemm!(convert(Complex128, alpha), M.data, b.data, convert(Complex128, beta), result.data)
operators.gemm!(alpha, a::DenseOperator, M::SparseOperator, beta, result::DenseOperator) = sparsematrix.gemm!(convert(Complex128, alpha), a.data, M.data, convert(Complex128, beta), result.data)
operators.gemv!(alpha, M::SparseOperator, b::Ket, beta, result::Ket) = sparsematrix.gemv!(convert(Complex128, alpha), M.data, b.data, convert(Complex128, beta), result.data)
operators.gemv!(alpha, b::Bra, M::SparseOperator, beta, result::Bra) = sparsematrix.gemv!(convert(Complex128, alpha), b.data, M.data, convert(Complex128, beta), result.data)

end # module


text:
module FileSessionAdapter

using Genie.Sessions, Genie, Genie.Configuration
import Serialization, Logging

const SESSIONS_PATH = "sessions"

"""
    write(session::Sessions.Session) :: Sessions.Session

Persists the `Session` object to the file system, using the configured sessions folder and returns it.
"""
function write(session::Sessions.Session) :: Sessions.Session
  if ! isdir(joinpath(SESSIONS_PATH))
    @warn "Sessions folder $(abspath(SESSIONS_PATH)) does not exist"
    @info "Creating sessions folder at $(abspath(SESSIONS_PATH))"
    mkpath(SESSIONS_PATH)
  end

  try
    open(joinpath(SESSIONS_PATH, session.id), "w") do (io)
      Serialization.serialize(io, session)
    end
  catch ex
    @error "Error serializing session"

    rethrow(ex)
  end

  session
end

"""
    read(session_id::Union{String,Symbol}) :: Union{Nothing,Sessions.Session}
    read(session::Sessions.Session) :: Union{Nothing,Sessions.Session}

Attempts to read from file the session object serialized as `session_id`.
"""
function read(session_id::Union{String,Symbol}) :: Union{Nothing,Sessions.Session}
  try
    isfile(joinpath(SESSIONS_PATH, session_id)) || return write(Session(session_id))
  catch ex
    @error "Can't check session file"
    @error ex

    write(Session(session_id))
  end

  try
    open(joinpath(SESSIONS_PATH, session_id), "r") do (io)
      Serialization.deserialize(io)
    end
  catch ex
    @error "Can't read session"
    @error ex

    write(Session(session_id))
  end
end
function read(session::Sessions.Session) :: Union{Nothing,Sessions.Session}
  read(session.id)
end

end


text:
#####
##### Activate the right project (usually `docs`)
#####

using Pkg
project_path = get(ENV, "INPUT_PROJECT", "")
if isempty(project_path)
    project_path = mktempdir()
end
Pkg.activate(project_path)

#####
##### `dev` the package code into our project
#####

package_path = get(ENV, "INPUT_PACKAGE_PATH", "")
if isempty(package_path)
    package_path = pwd()
end

Pkg.develop(PackageSpec(path=package_path))
Pkg.instantiate()

#####
##### parse the package name out of its Project.toml
#####

try
    using TOML # stdlib on Julia v1.6+
catch
    @info "Adding TOML..."
    Pkg.add("TOML")
    using TOML
end

name_str = TOML.parsefile(joinpath(package_path, "Project.toml"))["name"]
name_sym = Symbol(name_str)
@info "Package name: $(name_sym)"

#####
##### Load Documenter
#####

try
    using Documenter
catch
    @warn("Documenter.jl was not found in the `$(project_path)` environment. Adding latest release...")
    Pkg.add("Documenter")
    using Documenter
end

#####
##### Run the doctests with `fix=true`
#####

@eval begin
    using $(name_sym)
    DocMeta.setdocmeta!($(name_sym), :DocTestSetup, :(using $name_sym); recursive=true)
    doctest($(name_sym); fix=true)
end

#####
##### Reset the docs Project + Manifest so that they don't show up in the diff
#####

function reset(file)
    try
        path = joinpath(project_path, file)
        read(`git checkout HEAD -- $path`) # `read` instead of `run` to suppress output
        @info "Reset $path"
    catch
        # We may not be in a git repo, and the file may not be committed,
        # so we will ignore errors here.
    end
end

reset("Project.toml")
reset("Manifest.toml")

max_stars_repo_name: julia-actions/julia-fix-doctests

max_issues_repo_name: julia-actions/julia-fix-doctests

max_forks_repo_name: julia-actions/julia-fix-doctests


text:
using DataFrames
using CSV
using Plots
#%%
data_1 = open("ode05.dat")
data_2 = open("ode025.dat")
data_3 = open("ode01.dat")
data_4 = open("ode001.dat")

df1=CSV.read(data_1, header=false, delim=" ")
df2=CSV.read(data_2, header=false, delim=" ")
df3=CSV.read(data_3, header=false, delim=" ")
df4=CSV.read(data_4, header=false, delim=" ")

plot(df1[1],df1[2], label="h=0.5")
plot!(df2[1],df2[2], label="h=0.25")
plot!(df3[1],df3[2], label="h=0.01")
plot!(df4[1],df4[2], label="h=0.001")
# xaxis!(:log10)
# yaxis!(:log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
########################################################
#   Experimento del algoritmo inicialización aprendida aprendida con base de datos de air self noise
########################################################
using Random
using Plots
using TOML
using CSV
using DataFrames
using StatsBase
using HypothesisTests
using TimerOutputs

FICHERO_CONFIGURACION = "Experimentos/.config.toml"
config = TOML.parsefile(FICHERO_CONFIGURACION)["air-self-noise"]
FILE = config["FICHERO_DATOS"]
DIRECTORIO_RESULTADOS = config["DIRECTORIO_RESULTADOS"]
NUMERO_EJECUCIONES = config["NUMERO_EJECUCIONES"]

include("../../OptimizedNeuralNetwork.jl/src/OptimizedNeuralNetwork.jl")

using .OptimizedNeuralNetwork

#------------------------------------------------------
#                Data preprocesing
#------------------------------------------------------
input_dimension = 5
output_dimension = 1
atributes = 1:input_dimension
label = 6
df = DataFrame(CSV.File(FILE, header=false))
# Miramos si hay valores perdidos
display(describe(df)) # No hay
# Transformamos en matrices y vectores
X = Matrix(df[:, atributes])
Y = Vector(df[:,label])
len = length(Y)
index = Integer(2*len/3) #separador de conjunto de entrenamiento y test
α = 0.9 # heurística ortogonalidad comentada en memoria
n = ceil(Integer, index*α) # numbers of nodes
valor_estancamiento = 5
tol = 0.001
η = 0.005 # valor heurístico tras varias pruebas (0.1 era demasiado grande)

println("Se va a entrenar con $n neuronas")
println("El tamaño de test es de $(len-index)")
println("El conjunto de entrenamiento $(index)")

### Valores donde almacenar los datos
# índeces del array donde se almacenan los datos
indice_algoritmo_inicializcion = 1
indice_aleatorio_y_backpropagation = 2

dfTiempo = [
        Array{Float64}(undef, 2 )
        for _ in 1:NUMERO_EJECUCIONES
        ]
dfNombre = ["algoritmo inicialización aprendida", "Aleatorio y Backpropagation"]
dfErrorEntrenamiento = [
        Array{Float64}(undef, 2 )
        for _ in 1:NUMERO_EJECUCIONES
        ]
dfErrorTest = [
        Array{Float64}(undef, 2 )
        for _ in 1:NUMERO_EJECUCIONES
]
for i in 1:NUMERO_EJECUCIONES
        # suffle
        sort = randperm(len)
        Xs = X[sort, :]
        Ys = Y[sort]
        # separate train from test
        X_train = Xs[1:index,:]
        Y_train = Ys[1:index]
        X_test = Xs[index+1:len, :]
        Y_test = Ys[index+1:len]

        # Data normalization
        #dt = fit(ZScoreTransform, X_train, dims=1)
        dt = fit(UnitRangeTransform, X_train, dims=1)
        X_test_normalized = StatsBase.transform(dt, X_test)
        X_train_normalized = StatsBase.transform(dt, X_train)

        dt_y = fit(UnitRangeTransform, Y_train, dims=1)
        Y_test_normalized = StatsBase.transform(dt_y, Y_test)
        Y_train_normalized = StatsBase.transform(dt_y, Y_train)

        #------------------------------------------------------
        #                Get neuronal network
        #------------------------------------------------------
        println("\nEjecución $i")
        # Experimentamos con nuestro algoritmo
        M = 1
        dfTiempo[i][indice_algoritmo_inicializcion]= @elapsed h_initialized = nn_from_data(X_train_normalized, Y_train_normalized, n, M)
        # Función de evaluación por forward propagation
        evaluate_initialized(x) = forward_propagation(h_initialized,
                RampFunction,x
                )

        println("Resultados con h ajustada")
        error_in_train_initialize = error_in_data_set(X_train_normalized, Y_train_normalized, evaluate_initialized)
        println("Ha tardado un tiempo de $(dfTiempo[i][indice_algoritmo_inicializcion])")
        println("El error en el conjunto de entrenamiento es de $error_in_train_initialize")

        # Almacenamos datos en el csv
        dfErrorEntrenamiento[i][indice_algoritmo_inicializcion] = error_in_train_initialize
        dfErrorTest[i][indice_algoritmo_inicializcion] = error_in_data_set(X_test_normalized, Y_test_normalized, evaluate_initialized)

        ####### Test con backpropagation
        # Entrenaremos hasta que el error sea igual o menor
        # Ajuste con comienzo de datos inicial
        println("\n--- Resultados con h aleatoria ---")
        time_backpropagation = @elapsed h_random = RandomWeightsNN(input_dimension, n, output_dimension)
        evaluate_random(x) = forward_propagation(h_random,RampFunction,x)
        error_in_train_backpropagation = error_in_data_set(X_train_normalized, Y_train_normalized, evaluate_random)
        iterations = 0
        last_error = error_in_train_backpropagation
        stoped_iterations = 0
        println("El error inicial en backpropagation es de $error_in_train_backpropagation")
        while( error_in_train_initialize < error_in_train_backpropagation
                &&
                stoped_iterations < valor_estancamiento
        )
                println("El error en la iteración $iterations: $error_in_train_backpropagation")
                time_backpropagation += @elapsed backpropagation!(h_random,
                        X_train_normalized, Y_train_normalized,
                        RampFunction, derivativeRampFunction,
                        n,
                        η)
                iterations += 1
                evaluate_random(x) = forward_propagation(h_random,
                RampFunction,x
                )
                error_in_train_backpropagation = error_in_data_set(X_train_normalized, Y_train_normalized, evaluate_random)
                if(abs(error_in_train_backpropagation - last_error) < tol ||  error_in_train_backpropagation > tol + last_error)
                        stoped_iterations += 1
                else
                        stoped_iterations = 0
                end
                last_error = error_in_train_backpropagation
        end
        evaluate_random(x) = forward_propagation(h_random,RampFunction,x)
        dfTiempo[i][indice_aleatorio_y_backpropagation] = time_backpropagation
        dfErrorEntrenamiento[i][indice_aleatorio_y_backpropagation] = error_in_train_backpropagation
        dfErrorTest[i][indice_aleatorio_y_backpropagation] = error_in_data_set(X_test_normalized, Y_test_normalized, evaluate_random)
        println(regression(X_test_normalized, Y_test_normalized, evaluate_random))
        println("Durante $time_backpropagation")
end

# Mostramos resultados y los guardamos en le fichero de resultados
#Mostramos en pantalla resultados
DF_NOMBRES = DataFrame(
        Método = dfNombre,
)
# Añadimos los tiempos
DF_TIEMPOS = hcat(
        DF_NOMBRES,
        DataFrame(dfTiempo, ["Tiempo $(i)" for i in 1:NUMERO_EJECUCIONES])
)
# Añadimos error en entrenamiento
DF_ERROR_ENTRENAMIENTO = hcat(
        DF_NOMBRES,
        DataFrame(dfErrorEntrenamiento, ["Error entrenamiento $(i)" for i in 1:NUMERO_EJECUCIONES])
)

DF_ERROR_TEST = hcat(
                DF_NOMBRES,
                DataFrame(dfErrorTest, ["Error test $(i)" for i in 1:NUMERO_EJECUCIONES])
)

# Guardamos los datos en el directorio respectivo
CSV.write(DIRECTORIO_RESULTADOS*"tiempos.csv", DF_TIEMPOS)
CSV.write(DIRECTORIO_RESULTADOS*"error_entrenamiento.csv", DF_ERROR_ENTRENAMIENTO)
CSV.write(DIRECTORIO_RESULTADOS*"error_test.csv", DF_ERROR_TEST)

# Test de los signos de Wilcoxon

resultados = "\n
        $(SignedRankTest(map(x-> x[1]-x[2],dfTiempo)))
        "
println(resultados)
write(DIRECTORIO_RESULTADOS*"TEST_WILCOXON", resultados)


text:
#Julia Gadfly Histogram

using Gadfly
plot(x=randn(113), Geom.histogram(bincount=10))


text:
module GalacticMOI

using MathOptInterface, GalacticOptim, GalacticOptim.SciMLBase
const MOI = MathOptInterface

struct MOIOptimizationProblem{T,F<:OptimizationFunction,uType,P} <: MOI.AbstractNLPEvaluator
    f::F
    u0::uType
    p::P
    J::Matrix{T}
    H::Matrix{T}
    cons_H::Vector{Matrix{T}}
end

function MOIOptimizationProblem(prob::OptimizationProblem)
    num_cons = prob.ucons === nothing ? 0 : length(prob.ucons)
    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p, num_cons)
    T = eltype(prob.u0)
    n = length(prob.u0)
    return MOIOptimizationProblem(
        f,
        prob.u0,
        prob.p,
        zeros(T, num_cons, n),
        zeros(T, n, n),
        Matrix{T}[zeros(T, n, n) for i in 1:num_cons],
    )
end

MOI.features_available(::MOIOptimizationProblem) = [:Grad, :Hess, :Jac]

function MOI.initialize(
    moiproblem::MOIOptimizationProblem,
    requested_features::Vector{Symbol},
)
    available_features = MOI.features_available(moiproblem)
    for feat in requested_features
        if !(feat in available_features)
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
    return
end

function MOI.eval_objective(moiproblem::MOIOptimizationProblem, x)
    return moiproblem.f(x, moiproblem.p)
end

function MOI.eval_constraint(moiproblem::MOIOptimizationProblem, g, x)
    g .= moiproblem.f.cons(x)
    return
end

function MOI.eval_objective_gradient(moiproblem::MOIOptimizationProblem, G, x)
    moiproblem.f.grad(G, x)
    return
end

# This structure assumes the calculation of moiproblem.J is dense.
function MOI.jacobian_structure(moiproblem::MOIOptimizationProblem)
    rows, cols = size(moiproblem.J)
    return Tuple{Int,Int}[(i, j) for j in 1:cols for i in 1:rows]
end

function MOI.eval_constraint_jacobian(moiproblem::MOIOptimizationProblem, j, x)
    if isempty(j)
        return
    elseif moiproblem.f.cons_j === nothing
        error(
            "Use OptimizationFunction to pass the derivatives or " *
            "automatically generate them with one of the autodiff backends",
        )
    end
    moiproblem.f.cons_j(moiproblem.J, x)
    for i in eachindex(j)
        j[i] = moiproblem.J[i]
    end
    return
end

# Because the Hessian is symmetrical, we choose to store the upper-triangular
# component. We also assume that it is dense.
function MOI.hessian_lagrangian_structure(moiproblem::MOIOptimizationProblem)
    num_vars = length(moiproblem.u0)
    return Tuple{Int,Int}[(row, col) for col in 1:num_vars for row in 1:col]
end

function MOI.eval_hessian_lagrangian(
    moiproblem::MOIOptimizationProblem{T},
    h,
    x,
    σ,
    μ,
) where {T}
    n = length(moiproblem.u0)
    if iszero(σ)
        fill!(h, zero(T))
    else
        moiproblem.f.hess(moiproblem.H, x)
        k = 0
        for col in 1:n
            for row in 1:col
                k += 1
                h[k] = σ * moiproblem.H[row, col]
            end
        end
    end
    if !isempty(μ) && !all(iszero, μ)
        moiproblem.f.cons_h(moiproblem.cons_H, x)
        for (μi, Hi) in zip(μ, moiproblem.cons_H)
            k = 0
            for col in 1:n
                for row in 1:col
                    k += 1
                    h[k] += μi * Hi[row, col]
                end
            end
        end
    end
    return
end

_create_new_optimizer(opt::MOI.AbstractOptimizer) = opt
_create_new_optimizer(opt::MOI.OptimizerWithAttributes) = MOI.instantiate(opt)

function __map_optimizer_args(
    prob::OptimizationProblem,
    opt::Union{MOI.AbstractOptimizer,MOI.OptimizerWithAttributes};
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...
)
    optimizer = _create_new_optimizer(opt)
    for (key, value) in kwargs
        MOI.set(optimizer, MOI.RawOptimizerAttribute("$(key)"), value)
    end
    if !isnothing(maxtime)
        MOI.set(optimizer, MOI.TimeLimitSec(), maxtime)
    end
    if !isnothing(reltol)
        @warn "common reltol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments."
    end
    if !isnothing(abstol)
        @warn "common abstol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments."
    end
    if !isnothing(maxiters)
        @warn "common maxiters argument is currently not used by $(optimizer). Set number of interations via optimizer specific keyword aguments."
    end
    return optimizer
end

function SciMLBase.__solve(
    prob::OptimizationProblem,
    opt::Union{MOI.AbstractOptimizer,MOI.OptimizerWithAttributes};
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...
)
    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)
    opt_setup = __map_optimizer_args(
        prob,
        opt;
        abstol=abstol,
        reltol=reltol,
        maxiters=maxiters,
        maxtime=maxtime,
        kwargs...
    )
    num_variables = length(prob.u0)
    θ = MOI.add_variables(opt_setup, num_variables)
    if prob.lb !== nothing
        @assert eachindex(prob.lb) == Base.OneTo(num_variables)
        for i in 1:num_variables
            if prob.lb[i] > -Inf
                MOI.add_constraint(opt_setup, θ[i], MOI.GreaterThan(prob.lb[i]))
            end
        end
    end
    if prob.ub !== nothing
        @assert eachindex(prob.ub) == Base.OneTo(num_variables)
        for i in 1:num_variables
            if prob.ub[i] < Inf
                MOI.add_constraint(opt_setup, θ[i], MOI.LessThan(prob.ub[i]))
            end
        end
    end
    if MOI.supports(opt_setup, MOI.VariablePrimalStart(), MOI.VariableIndex)
        @assert eachindex(prob.u0) == Base.OneTo(num_variables)
        for i in 1:num_variables
            MOI.set(opt_setup, MOI.VariablePrimalStart(), θ[i], prob.u0[i])
        end
    end
    MOI.set(
        opt_setup,
        MOI.ObjectiveSense(),
        prob.sense === GalacticOptim.MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE,
    )
    if prob.lcons === nothing
        @assert prob.ucons === nothing
        con_bounds = MOI.NLPBoundsPair[]
    else
        @assert prob.ucons !== nothing
        con_bounds = MOI.NLPBoundsPair.(prob.lcons, prob.ucons)
    end
    MOI.set(
        opt_setup,
        MOI.NLPBlock(),
        MOI.NLPBlockData(con_bounds, MOIOptimizationProblem(prob), true),
    )
    MOI.optimize!(opt_setup)
    if MOI.get(opt_setup, MOI.ResultCount()) >= 1
        minimizer = MOI.get(opt_setup, MOI.VariablePrimal(), θ)
        minimum = MOI.get(opt_setup, MOI.ObjectiveValue())
        opt_ret = Symbol(string(MOI.get(opt_setup, MOI.TerminationStatus())))
    else
        minimizer = fill(NaN, num_variables)
        minimum = NaN
        opt_ret = :Default
    end
    return SciMLBase.build_solution(
        prob,
        opt,
        minimizer,
        minimum;
        original=opt_setup,
        retcode=opt_ret
    )
end

end


text:
function functionalize_fieldnames(x)
  for field in fieldnames(x)
    @eval $field(y::$x) = y.$field
    @eval export $field
  end
end

methods_to_functionalize = [PCB_Specification,
                            WindingGeometry,
                            Windings,
                            Transformer,
                            TransformerPowerAnalysis,
                            CoreGeometry,
                            FerriteProperties,
                            Conductor,
                            Dielectric,
                            Stackup
                            ]

import Base.length
for m in methods_to_functionalize
  functionalize_fieldnames(m)
end


text:
module utilities

export c_0, ϵ_0, peak_width

c_0 = 299792458
ϵ_0 = 8.854187812813e-12

function peak_width(data::Vector{<:Real}, step::T, inverted::Bool)::T where T<:Real
    data .-= minimum(data)
    if inverted
        data = maximum(data) .- data
    end

    target = maximum(data) / 2

    previous::Real = 0
    left::UInt = 0
    right::UInt = 0
    for (i, current) in zip(1:1:length(data), data)
        if left == 0
            if (previous < target) & (current >= target)
                left = i
            end
        elseif right == 0
            if (previous > target) & (current <= target)
                right = i
            end
        else
            break
        end
        previous = current
    end

    (right - left) * step
end

end # utilities


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    trivialhydrogens,
    allhydrogens,
    makehydrogensimplicit,
    makehydrogensexplicit,
    largestcomponentnodes,
    largestcomponentgraph,
    neutralizeacids!,
    neutralizeoniums!,
    depolarize!,
    triplebondanion!,
    canonicalize!,
    canonicalize

# TODO: large conjugated system
# TODO: salts and waters should detected by functional group analysis
# TODO: Phosphate, diphosphate, sulfate, nitrate, acetate,
# maleate, fumarate, succinate, citrate, tartrate, oxalate,
# mesylate, tosylate, besylate,
# benzoate, gluconate

"""
    trivialhydrogens(mol::GraphMol) -> Set{Int}

Return a set of trivial hydrogen nodes (light hydrogens which are uncharged,
non-radical, non-stereospecific and attached to organic heavy atoms)
"""
function trivialhydrogens(mol::GraphMol)
    hs = Set{Int}()
    organic_heavy = (
        :B, :C, :N, :O, :F, :Si, :P, :S, :Cl, :As, :Se, :Br, :I)
    for (i, atom) in enumerate(nodeattrs(mol))
        if (atom.symbol != :H || atom.charge != 0 || atom.multiplicity != 1
                || atom.mass !== nothing)
            continue
        elseif atom isa SmilesAtom && atom.stereo !== nothing
            continue
        end
        degree(mol, i) == 1 || continue
        (inc, adj) = iterate(neighbors(mol, i))[1]
        bond = edgeattr(mol, inc)
        nbratom = nodeattr(mol, adj)
        if bond isa SDFileBond && (bond.order != 1 || bond.notation != 0)
            continue
        elseif !(nbratom.symbol in organic_heavy)
            continue
        end
        push!(hs, i)
    end
    return hs
end

"""
    allhydrogens(mol::GraphMol) -> Set{Int}

Return a set of hydrogen nodes.
"""
function allhydrogens(mol::GraphMol)
    hs = Set{Int}()
    for (i, atom) in enumerate(nodeattrs(mol))
        if atom.symbol == :H
            push!(hs, i)
        end
    end
    return hs
end

"""
    makehydrogensimplicit(mol::GraphMol) -> GraphMol

Return molecule whose hydrogen nodes are removed. If option `all` is set to
false, only trivial hydrogens are removed (see [`trivialhydrogens`](@ref)).
"""
function makehydrogensimplicit(mol::GraphMol; all=true)
    hydrogens = all ? allhydrogens : trivialhydrogens
    ns = setdiff(nodeset(mol), hydrogens(mol))
    return graphmol(nodesubgraph(mol, ns))
end

"""
    makehydrogensexplicit(mol::GraphMol) -> GraphMol

Return molecule whose hydrogens are fully attached. If option `all` is set to
false, only trivial hydrogens are removed (see [`trivialhydrogens`](@ref)).
"""
function makehydrogensexplicit(mol::GraphMol)
    newmol = graphmol(mol)
    implicithcount_ = implicithcount(mol)
    for (i, atom) in enumerate(nodeattrs(mol))
        for j in 1:implicithcount_[i]
            n = addnode!(newmol, nodeattrtype(mol)(:H))
            addedge!(newmol, i, n, edgeattrtype(mol)())
        end
    end
    return newmol
end

"""
    largestcomponentnodes(mol::GraphMol) -> Set{Int}

Return a set of nodes in the largest connected component.
"""
largestcomponentnodes(mol::GraphMol
    ) = sortstablemax(connected_components(mol), by=length, init=Set{Int}())

"""
    largestcomponentgraph(mol::GraphMol) -> GraphMol

Return largest connected component of the molecular graph.
"""
largestcomponentgraph(mol::GraphMol
    ) = nodesubgraph(mol, largestcomponentnodes(mol))

"""
    neutralizeacids!(mol::GraphMol)

Neutralize oxo(thio) acids.

Note that this function edits `Atom` object fields directly. The molecular
property vector needs recalculation to apply the changes.
see [`canonicalize!`](@ref).
"""
function neutralizeacids!(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    connectivity_ = connectivity(mol)
    pielectron_ = pielectron(mol)
    for o in findall(
            (atomsymbol_ .== :O) .* (charge_ .== -1) .* (connectivity_ .== 1))
        nbr = iterate(adjacencies(mol, o))[1]
        if pielectron_[nbr] == 1
            cnbrs = adjacencies(mol, nbr)
            pop!(cnbrs, o)
            for cn in cnbrs
                if (atomsymbol_[cn] in (:O, :S) && pielectron_[cn] == 1
                        && connectivity_[cn] == 1)
                    setnodeattr!(mol, o, setcharge(nodeattr(mol, o), 0))
                    break
                end
            end
        end
    end
end

"""
    neutralizeoniums!(mol::GraphMol)

Neutralize 1-3° oniums. Permanently charged quart-oniums are not neutralized.

Note that this function edits `Atom` object fields directly. The molecular
property vector needs recalculation to apply the changes.
see [`canonicalize!`](@ref).
"""
function neutralizeoniums!(mol::GraphMol)
    for o in findall((charge(mol) .== 1) .* (hcount(mol) .> 0))
        setnodeattr!(mol, o, setcharge(nodeattr(mol, o), 0))
    end
end

"""
    depolarize!(mol::GraphMol)

Depolarize oxo groups except in the case that polarization is required for
aromaticity.

Note that this function edits `Atom` object fields directly. The molecular
property vector needs recalculation to apply the changes.
see [`canonicalize!`](@ref).
"""
function depolarize!(mol::GraphMol)
    charge_ = charge(mol)
    isaromatic_ = isaromatic(mol)
    for o in findall((atomsymbol(mol) .== :O) .* (charge_ .== -1))
        @assert degree(mol, o) == 1 "unexpected oxygen degree $(length(nbrs))"
        (inc, adj) = iterate(neighbors(mol, o))[1]
        if charge_[adj] == 1 && !isaromatic_[adj]
            setnodeattr!(mol, o, setcharge(nodeattr(mol, o), 0))
            setnodeattr!(mol, adj, setcharge(nodeattr(mol, adj), 0))
            setedgeattr!(mol, inc, setorder(edgeattr(mol, inc), 2))
        end
    end
end

"""
    triplebondanion!(mol::GraphMol)

Canonicalize anions next to triple bonds (ex. [C-][N+]#N -> C=[N+]=[N-]).

Note that this function edits `Atom` object fields directly. The molecular
property vector needs recalculation to apply the changes.
see [`canonicalize!`](@ref).
"""
function triplebondanion!(mol::GraphMol)
    # TODO: better function name
    charge_ = charge(mol)
    for tb in findall(bondorder(mol) .== 3)
        tbond = edgeattr(mol, tb)
        (u, v) = getedge(mol, tb)
        for (f, s) in ((u, v), (v, u))
            nbrs = copy(neighbors(mol, f))
            pop!(nbrs, findedgekey(mol, f, s))
            length(nbrs) == 1 || continue
            (inc, adj) = iterate(nbrs)[1]
            if charge_[adj] == -1
                setnodeattr!(mol, adj, setcharge(nodeattr(mol, adj), 0))
                setnodeattr!(mol, s, setcharge(nodeattr(mol, s), -1))
                setedgeattr!(mol, inc, setorder(edgeattr(mol, inc), 2))
                setedgeattr!(mol, tb, setorder(edgeattr(mol, tb), 2))
            end
        end
    end
end

"""
    canonicalize!(mol::GraphMol)

Canonicalize molecule notation and apply the changes to the molecular property
vector.

- Neutralize oxo acid, 1-3° ammonium and polarized carbonyls except in the
  case that polarization is required for aromaticity.
- Canonicalize anions next to triple bonds (ex. [C-][N+]#N -> C=[N+]=[N-])
"""
function canonicalize!(mol::GraphMol)
    neutralizeacids!(mol)
    neutralizeoniums!(mol)
    depolarize!(mol)
    triplebondanion!(mol)
    return
end

function canonicalize(mol::GraphMol)
    newmol = graphmol(mol)
    canonicalize!(newmol)
    return newmol
end


text:
using CSV, Distributions, HypothesisTests

data = CSV.read("../data/machine1.csv", header=false)[:,1]
n = length(data)
mu0A, mu0B = 52.2, 53

xPositiveA = sum(data .> mu0A)
testStatisticA = max(xPositiveA, n-xPositiveA)

xPositiveB = sum(data .> mu0B)
testStatisticB = max(xPositiveB, n-xPositiveB)

binom = Binomial(n,0.5)
pValA = 2*ccdf(binom, testStatisticA)
pValB = 2*ccdf(binom, testStatisticB)

println("Binomial mean: ", mean(binom))

println("Case A: mu0: ", mu0A)
println("\tTest statistc: ", testStatisticA)
println("\tP-value: ", pValA)

println("Case B: mu0: ", mu0B)
println("\tTest statistc: ", testStatisticB)
println("\tP-value: ", pValB)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook