text:
function mcmc_run!(x_current,P_in,p,q_vec,q_mat,u_mat,ind_row::T,ind_chain::T) where {T<:Integer}
    accept_val = 0.0
    P_current = P_in
    x_new = x_current.+(q_mat[ind_row,:,ind_chain].*q_vec)
    p_new = p(x_new)
    A = min(1,p_new/P_current) #calcualte acceptance ratio
    u = u_mat[ind_row,ind_chain]
        if u<=A
            x_current[1:end] = x_new
            P_current = p_new
            accept_val = 1.0
        end

    return P_current,accept_val
end

function n_mcmc_runs!(x_current,P_in,p,q_vec,q_mat,u_mat,ind_row_start::T,ind_chain::T,n_runs::T) where {T<:Integer}
    n_accepts = 0.0
    P_current = P_in
    for i = 1:n_runs
        ind_row=ind_row_start+i
        P_current,accept_val = mcmc_run!(x_current,P_current,p,q_vec,q_mat,u_mat,ind_row,ind_chain)
        n_accepts += accept_val
    end
    return P_current,n_accepts
end

function metropolis_mcmc(x_current_vec,P_current_vec,p,q_vec::Array{Float64,1},file_name::String
    ;n_chains::Integer=7,n_samples::Integer=100,burn_in::Integer=100,sample_freq::Integer=10)

    n_draws = n_samples*sample_freq+burn_in #number of random number per chain
    n_features = length(q_vec)

    u_mat = rand(n_draws,n_chains)
    q_mat = rand(Normal(),(n_draws,n_features,n_chains))

    samples_container = Array{Array{typeof(x_current_vec[1]),1},1}(undef,n_chains)
    P_samples_container = Array{Array{Float64,1},1}(undef,n_chains)

    Threads.@threads for j = 1:n_chains
        core_id = j

        x_current,P_current = copy(x_current_vec[j]), copy(P_current_vec[j])

        #Run burn-in
        println("Core id $(core_id): Start Burn-in")
        P_current,na = n_mcmc_runs!(x_current,P_current,p,q_vec,q_mat,u_mat,0,j,burn_in)
        println("Core id $(core_id): Done Burn-in")

        #Create container for samples
        samples = Array{typeof(x_current),1}(undef,n_samples)
        #Propability value of samples
        P_samples = zeros(n_samples)
        accept_rate = 0.0

        println("Core id $(core_id): Start Sampling")
        #Run samples
        for i = 1:n_samples
            println("Core id $(core_id): ---Sample: $(i)/$(n_samples)")
            ind_row_start = burn_in+(i-1)*sample_freq
            P_current,na = n_mcmc_runs!(x_current,P_current,p,q_vec,q_mat,u_mat,ind_row_start,j,sample_freq)
            accept_rate += na
            #Save sample
            samples[i] = copy(x_current)
            #Save propability values
            P_samples[i] = P_current
        end
        println("Core id $(core_id): Done Sampling")
        accept_rate /= n_samples*sample_freq
        println("Core id $(core_id): accept_rate = $(accept_rate)")

        samples_container[j] = samples
        P_samples_container[j] = P_samples
    end
    save("./output/"*file_name*".jld","samples_container",samples_container,"P_samples_container",P_samples_container)
    return nothing
end


text: # some more code

ext: jl

max_stars_repo_name: tomtrogdon/demo

max_issues_repo_name: tomtrogdon/demo

max_forks_repo_name: tomtrogdon/demo


text:
using ZulipReminderBot

include("configuration.jl")

ZulipReminderBot.up(conn)

# There should be `down` function and some way to control db version...


text:
#This file executes the TD3 algorithm: https://arxiv.org/pdf/1802.09477.pdf
#This is more advanced than the simple DDPG algorithm!!

include("src/init.jl")# Include all relevant functions

include("src/td3/episodes.jl")
include("src/td3/gradients.jl")
include("src/td3/updates.jl")

res  = zeros(params.n_iter)
prec = zeros(params.n_iter)

#=NIPS-challenge
include("src/NIPS2018.jl")
#Global Variable of Environment!
@everywhere @pyimport osim.env as osim
@everywhere env = osim.ProstheticsEnv(visualize=false)
env = osim.ProstheticsEnv(visualize=true)
@everywhere env[:change_model](model="2D", prosthetic=false, difficulty=0)
sleep(1)
# =#

#1. Randomly Initialize
#-------------------------------------------------------------
#1.1 First Critic
Qf_1         = init_θ_single_norm(params.structure_net_qf)
Qf_tar_1     = copy(Qf_1)
opt_para_qf1 = optimizers(Qf_1, Adam; lr = 0.5*params.lr_qf)

#1.1 Second Critic
Qf_2         = init_θ_single_norm(params.structure_net_qf2)
Qf_tar_2     = copy(Qf_2)
opt_para_qf2 = optimizers(Qf_2, Adam; lr = 2.0*params.lr_qf)

#1.3 Policy network
Pl           = init_θ_single_norm(params.structure_net_pl)
Pl_tar       = copy(Pl)
opt_para_pl  = optimizers(Pl, Adam; lr = params.lr_pl)

#1.4 Replay Buffer and Latest Trajectories
buffer           = Traj_Data()
buffer.traj      = Array{Trajectory,1}(undef,0)
buffer.timesteps = 0
buffer.n_traj    = 0

#1.5 Buffer for the most recent trajectories
all_traj  = init_all_traj(params)

#using FileIO
#init_nets = load("Results/Videos/10_20_init_nets/mountaincar.jld2")
#Pl = init_nets["θa"]

#main function running DDPG
function TD3()
    weights_old = zeros(length(net_to_vec(Pl)),2)
    a_t_old = all_traj.traj[1].a_t

    #1.4 Init running mean std
    rms_all = init_rms(params.s_init)

    #Iterate
    for i_iter = 1:params.n_iter
        #Init Progress bar
        progress = Progress(params.n_iter,1)

        if floor(i_iter/100) == ceil(i_iter/100)
            #make_video(all_traj,params, i_iter)
            save_data(res, prec, Pl, Qf_1, all_traj, i_iter)
        end

        #2.Collect Data using current stochastic Policy
        #-------------------------------------------------------------
        res[i_iter] = parallel_episodes_td3(params, Pl, Qf_1, Qf_2, all_traj,rms_all)
        #res[i_iter], dump = parallel_episodes_skeleton(params,Pl,Pl_tar,all_traj,rms_all)

        #3.Pump data in Replay Buffer:
        #-------------------------------------------------------------
        #3.1 Generate minibatches from current data_batch
        mini_batch = build_minibatches(all_traj,params)
        #3.2 Pump the data into the buffer
        pump_buffer!(buffer,mini_batch,params)

        #4.Train all the stuff now:
        #-------------------------------------------------------------
        do_td3_update_parallel!(Qf_1, Qf_tar_1, Qf_2, Qf_tar_2, Pl, Pl_tar, buffer, params, opt_para_pl, opt_para_qf1, opt_para_qf2)

        #5.Visualize Progress
        ProgressMeter.update!(progress,i_iter) #update progress bar

        #6. Plot data of interest
        if floor(i_iter/5) == ceil(i_iter/5)
            print("  new visualization: ")
            a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
            weights_old, i_iter,Pl, Qf_1, res, prec/params.n_traj, Pl, Pl)
            #run_episode_skeleton_demo!(params,Pl,Pl_tar,all_traj.traj[1],rms_all,env)
        end
        print("  Result:  ", res[i_iter])
    end
end

TD3()

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# Exports

export sample_requests_uniform
export sample_requests_file
export SolveSettings
export sat_sim
export mass_sim
export sat_resource_sim
export mass_resource_sim

function sample_requests_uniform(num_samples::Int=0;lon_min::Real=-180, lon_max::Real=180,
                                lat_min::Real=-90, lat_max::Real=90,
                                look_angle_min::Real=10, look_angle_max::Real=55,
                                collect_duration_min::Real=30, collect_duration_max::Real=30,
                                reward_min::Real=1.0, reward_max::Real=1.0)

    if num_samples <= 0
        throw(ErrorException("Number of samples must be positive."))
    end

    if reward_min <= 0
        throw(ErrorException("Minimum reward must be positive."))
    end

    if reward_max < reward_min
        throw(ErrorException("Maximum reward must be greater or equal to minimum reward."))
    end

    if collect_duration_min <= 0
        throw(ErrorException("Minimum collect duration must be positive."))
    end

    if collect_duration_max < collect_duration_min
        throw(ErrorException("Maximum collect duration must be greater or equal to minimum collect duration."))
    end

    # Array of requests to sample
    requests = Request[]

    # Sample locations
    lons = rand(Uniform(lon_min, lon_max), num_samples)
    lats = rand(Uniform(lat_min, lat_max), num_samples)

    if collect_duration_min == collect_duration_max
        cdur = [collect_duration_min for idx in 1:num_samples]
    else
        cdur = rand(Uniform(collect_duration_min, collect_duration_max), num_samples)
    end

    if reward_min == reward_max
        rewards = [reward_min for idx in 1:num_samples]
    else
        rewards = rand(Uniform(reward_min, reward_max), num_samples)
    end

    for i in 1:num_samples
        push!(requests, Request(lons[i], lats[i], look_angle_min=look_angle_min,
                look_angle_max=look_angle_max, collect_duration=cdur[i],
                reward=rewards[i], use_degrees=true, id=i))
    end

    return requests
end

function sample_requests_file(num_samples::Int=0, filepath::String="";
            look_angle_min::Real=10, look_angle_max::Real=55,
            collect_duration_min::Real=30, collect_duration_max::Real=30,
            reward_min::Real=1.0, reward_max::Real=1.0)

    if reward_min <= 0
        throw(ErrorException("Minimum reward must be positive."))
    end

    if reward_max < reward_min
        throw(ErrorException("Maximum reward must be greater or equal to minimum reward."))
    end

    if collect_duration_min <= 0
        throw(ErrorException("Minimum collect duration must be positive."))
    end

    if collect_duration_max < collect_duration_min
        throw(ErrorException("Maximum collect duration must be greater or equal to minimum collect duration."))
    end

    json_requests = JSON.parsefile(filepath)
    num_images = length(json_requests)

    if num_samples > num_images
        throw(ErrorException("Number of requested samples greater than data points in file"))
    end

    @debug "Found $num_images images in file"
    indexes = collect(Int, 1:num_images)

    # Presample collect durations and rewards
    if collect_duration_min == collect_duration_max
        cdur = [collect_duration_min for idx in 1:num_samples]
    else
        cdur = rand(Uniform(collect_duration_min, collect_duration_max), num_samples)
    end

    if reward_min == reward_max
        rewards = [reward_min for idx in 1:num_samples]
    else
        rewards = rand(Uniform(reward_min, reward_max), num_samples)
    end

    # Array of requests to sample
    requests = Request[]

    observed_idx = []
    while length(requests) < num_samples
        idx = rand(indexes, 1)[1]
        if !(idx in observed_idx)
            # Create request
            i = length(requests) + 1
            lon = json_requests[idx]["lon"]
            lat = json_requests[idx]["lat"]
            req = Request(lon, lat, look_angle_min=look_angle_min,
                look_angle_max=look_angle_max, collect_duration=cdur[i],
                reward=rewards[i], use_degrees=true, id=i)

            push!(requests, req)

            # Delete index to prevent resampling
            deleteat!(indexes, findfirst(x -> x == idx, indexes))
            # println(length(indexes))
        end
    end

    return requests
end


text:
__precompile__()

module JondeaUtils

export real_cols_to_complex,
       flatten,
       rel_abs_error,
       max_rel_abs_error,
       offsetmap,
       diffhankelh1,
       diffbesselj,
       fresnels,
       fresnelc,
       reds,
       greens,
       blues,
       angle,
       sqrt,
       round,
       ceil,
       floor

import OffsetArrays: OffsetArray
import SpecialFunctions: hankelh1, besselj, erf
import Colors: RGB
import StaticArrays: SVector

import Base.angle
import Base.sqrt
import Base.round
import Base.ceil
import Base.floor

include("special_functions.jl")
include("array.jl")
include("plot.jl")
include("maths.jl")

end # module


text:
mutable struct Throttle{F}
    f::F
    n::Int
    state::Int
end

function Throttle(f::Function, n::Int)
    return Throttle(f, n, 0)
end

function (t::Throttle)()
    iszero(rem(t.state, t.n)) && t.f()
    t.state += 1
end

mutable struct ModelSaver{M}
    model::M
    n::Int
    path::String
    state::Int
end

function ModelSaver(model, n::Int=500; path=pwd())
    return ModelSaver(model, n, path, 0)
end

function (m::ModelSaver)()
    if iszero(rem(m.state, m.n))
        path = joinpath(m.path, "model_$(m.state).jld2")
        FileIO.save(path, Dict("model" => m.model))
    end
    m.state += 1
end

mutable struct ExploitabilityCallback{SOL<:DeepCFRSolver, ESOL<:ExploitabilitySolver}
    sol::SOL
    e_sol::ESOL
    n::Int
    state::Int
    hist::CFR.ExploitabilityHistory
end

function ExploitabilityCallback(sol::DeepCFRSolver, n::Int=1; p::Int=1)
    e_sol = ExploitabilitySolver(sol, p)
    return ExploitabilityCallback(sol, e_sol, n, 0, CFR.ExploitabilityHistory())
end

function (cb::ExploitabilityCallback)()
    if iszero(rem(cb.state, cb.n))
        sol = cb.sol
        initialize!(sol.Π)
        train_policy!(sol, sol.strategy_batches, deepcopy(sol.strategy_opt))
        e = CFR.exploitability(cb.e_sol, sol)
        push!(cb.hist, cb.state, e)
    end
    cb.state += 1
end

@recipe f(cb::ExploitabilityCallback) = cb.hist

struct FittingCallback{SOL<:AbstractDeepCFRSolver,N}
    sol::SOL
    verbose::Bool
    io::IO
    hist::NTuple{N,Vector{Float64}}
end

function FittingCallback(sol; verbose=true, io=stderr)
    return FittingCallback(sol, verbose, io, Tuple(Float64[] for _ in eachindex(sol.V)))
end

function (cb::FittingCallback)()
    sol = cb.sol
    io = cb.io
    cb.verbose && println()
    for p in eachindex(sol.Mv)
        d = optimality_distance(sol.V[p], sol.Mv[p])
        push!(cb.hist[p], d)
        if cb.verbose
            println(io, p, ": ", round(d, sigdigits=3))
        end
    end
end

@recipe f(cb::FittingCallback) = hcat(cb.hist...)


text:
# User-defined joint
mutable struct GenericJoint{T,N} <: AbstractJoint{T,N}
    spring::T
    damper::T

    F::SVector{3,T}
    τ::SVector{3,T}

    childid::Int64

    g::Function

    function GenericJoint{N}(body1::AbstractBody{T}, body2::AbstractBody, gfunc; spring = zero(T), damper = zero(T)) where {T,N}
        F = zeros(T,3)
        τ = zeros(T,3)
        childid = body2.id

        new{T,N}(spring, damper, F, τ, childid, gfunc), body1.id, body2.id
    end

end

### General functions
@inline constraintmat(::GenericJoint) = I

### Constraints and derivatives
## Discrete-time position wrappers (for dynamics)
g(joint::GenericJoint, statea::State, stateb::State, Δt) = joint.g(joint, posargsnext(statea, Δt)..., posargsnext(stateb, Δt)...)
g(joint::GenericJoint, stateb::State, Δt) = joint.g(joint, posargsnext(stateb, Δt)...)
g(joint::GenericJoint, statea::State, stateb::State) = joint.g(joint, posargsc(statea)..., posargsc(stateb)...)
g(joint::GenericJoint, stateb::State) = joint.g(joint, posargsc(stateb)...)

## Derivatives NOT accounting for quaternion specialness
@inline function ∂g∂posa(joint::GenericJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    z = [xa;params(qa);xb;params(qb)]

    G = ForwardDiff.jacobian(x -> joint.g(joint,x[1:3],UnitQuaternion(x[4:7]...,false),x[8:10],UnitQuaternion(x[11:14]...,false)), z)
    X = G[:,1:3]
    Q = G[:,4:7]

    return X, Q
end
@inline function ∂g∂posb(joint::GenericJoint, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    z = [xa;params(qa);xb;params(qb)]

    G = ForwardDiff.jacobian(x -> joint.g(joint,x[1:3],UnitQuaternion(x[4:7]...,false),x[8:10],UnitQuaternion(x[11:14]...,false)), z)
    X = G[:,8:10]
    Q = G[:,11:14]

    return X, Q
end
@inline function ∂g∂posb(joint::GenericJoint, xb::AbstractVector, qb::UnitQuaternion)
    z = [xb;params(qb)]

    G = ForwardDiff.jacobian(x -> joint.g(joint,x[1:3],UnitQuaternion(x[4:7]...,false)), z)
    X = G[:,1:3]
    Q = G[:,4:7]

    return X, Q
end


text:
#### CSV Reader ####
# These are temporary methods to interpret CSV files
# We will probably revert to using the DelimitedFiles module and use actual CSV
# But for now, we are using semicolon delimitetd DSV

function parseline(filepath::AbstractString, selected_line::Int64)
    # Function to read a line from a CSV
    # Returns a dictionary with Float64 values where possible
    if !isfile(filepath)
        error("File does not exist at " * filepath)
    end
    open(filepath) do file
        linecount = 1
        for line in eachline(file)
            if linecount == selected_line
                return tryparse_fallback(split(line, ';'))
            end
            linecount += 1
        end
        error("File only contains " * string(linecount-1) * " lines")
    end
end

function findmatches(filepath::AbstractString, match_string::String, selected_col::Int64 = 1)
    # Function to search for a string in the specified column in a CSV
    # Returns an array of line numbers for where string is found
    if !isfile(filepath)
        error("File does not exist")
    end
    found_lines = Array{Int64, 1}()
    open(filepath) do file
        linecount = 1
        for line in eachline(file)
            if match_string == split(line, ';'; limit = selected_col+1)[selected_col]
                push!(found_lines, linecount)
            end
            linecount += 1
        end
        return found_lines
    end
end

function findmatches(filepath::AbstractString, match_string::String, selected_col::String; header_row::Int64 = 1)
    if !isfile(filepath)
        error("File does not exist")
    end
    selected_col_numbers = findall(isequal(selected_col), parseline(filepath, header_row))
    if length(selected_col_numbers) == 0
        error("Header " * selected_col * " does not exist.")
    elseif length(selected_col_numbers) > 1
        error("Header " * selected_col * " is not unique")
    end
    return findmatches(filepath, match_string, selected_col_numbers[1])
end

function findmatches_pair(filepath::AbstractString, match_string1::String, match_string2::String,selected_col1::Int64 = 1, selected_col2::Int64 = 2; ordered = false)
    # Function to search for pairs of strings in two specified columns
    # Returns an array of line numbers for where string is found
    if !isfile(filepath)
        error("File does not exist")
    end
    found_lines = Array{Int64, 1}()
    match_set = Set([match_string1, match_string2])
    open(filepath) do file
        linecount = 1
        for line in eachline(file)
            row_entries = split(line, ';')
            if ordered
                if match_string1 == row_entries[selected_col1] && match_string2 == row_entries[selected_col2]
                    push!(found_lines, linecount)
                end
            else
                if match_set == Set([row_entries[selected_col1], row_entries[selected_col2]])
                    push!(found_lines, linecount)
                end
            end
            linecount += 1
        end
        return found_lines
    end
end

function findmatches_pair(filepath::AbstractString, match_string1::String, match_string2::String, selected_col1::String, selected_col2::String; header_row::Int64 = 1, ordered = false)
    if !isfile(filepath)
        error("File does not exist")
    end
    selected_col_numbers1 = findall(isequal(selected_col1), parseline(filepath, header_row))
    selected_col_numbers2 = findall(isequal(selected_col2), parseline(filepath, header_row))
    if length(selected_col_numbers1) == 0
        error("Header " * selected_col1 * " does not exist.")
    elseif length(selected_col_numbers1) > 1
        error("Header " * selected_col1 * " is not unique")
    end
    if length(selected_col_numbers2) == 0
        error("Header " * selected_col2 * " does not exist.")
    elseif length(selected_col_numbers2) > 1
        error("Header " * selected_col2 * " is not unique")
    end
    return findmatches_pair(filepath, match_string1, match_string2, selected_col_numbers1[1], selected_col_numbers2[1]; ordered = ordered)
end

function tryparse_fallback(list)
    # Wrapper for tryparse to retain string if not parsable to Float, and return missing if empty string
    parsed = [tryparse(Float64, i) for i in list]
    return [list[i] == "" ? missing : parsed[i] == nothing ? list[i] : parsed[i] for i in 1:length(list)]
end


text:
## header to provide surface and context
using Cairo
c = CairoARGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgba(cr,0.0,0.0,0.0,0.0); # transparent black
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);
save(cr);

s1 = CairoRecordingSurface()
c1 = CairoContext(s1)

## original example, following here
move_to(c1, 128.0, 25.6);
line_to(c1, 230.4, 230.4);
rel_line_to(c1, -102.4, 0.0);
curve_to(c1, 51.2, 230.4, 51.2, 128.0, 128.0, 128.0);
close_path(c1);

move_to(c1, 64.0, 25.6);
rel_line_to(c1, 51.2, 51.2);
rel_line_to(c1, -51.2, 51.2);
rel_line_to(c1, -51.2, -51.2);
close_path(c1);

set_line_width(c1, 10.0);
set_source_rgb(c1, 0, 0, 1);
fill_preserve(c1);
set_source_rgb(c1, 0, 0, 0);
stroke(c1);

## now using the recorded painting
set_source(cr, s1, 0.0, 0.0)
paint(cr)

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_record0.png");


text:
# TODO: Hold on to constant parameters such as padding and stride for Convolutions
"""
$(SIGNATURES)

Extract the kernel parameters from `node`.
"""
params(backend::nGraph.Backend{nGraph.CPU}, node::nGraph.NodeLike) = CPUKernelParams(node)
params(backend::nGraph.Backend{nGraph.GPU}, node::nGraph.NodeLike) = GPUKernelParams(node)
params(backend, node::XNode) = params(backend, unx(node))

# Cache object for recording seen kernels.
"""
    CPUKernelParams

Parameter set for nGraph CPU Kernels.

$(METHODLIST)

$(FIELDS)
"""
struct CPUKernelParams{IS, OS, IT, OT, NIF}
    "The description of the op."
    description::String

    # IO Sizes
    "Sizes of the input tensors."
    input_sizes::IS
    "Sizes of the output tensors."
    output_sizes::OS
    "Element types of the input tensors."
    input_types::IT
    "Element types of the output tensors."
    output_types::OT

    # MKLDNN Formats
    "Flag indicating if node is a `mkldnn` op."
    ismkl::Bool
    "Enumj values of the tensor mkldnn input formats."
    input_formats::NTuple{NIF, Int64}
end

# Get the MKL format string for an op
mkldnn_string(x) = last(split(nGraph.Lib.get_mkldnn_string(x), ":"))

function CPUKernelParams(node::nGraph.NodeLike)
    description = nGraph.description(node)

    # Input processing
    num_inputs = nGraph.get_input_size(node)

    input_sizes = ntuple(x -> nGraph.get_input_shape(node, x), num_inputs)
    input_types = ntuple(x -> nGraph.get_input_element_type(node, x), num_inputs)

    # Only
    ismkl = nGraph.is_mkldnn(node)
    input_formats = ntuple(
        # At the moment, still forwarding to nGraph.Lib.
        x -> nGraph.Lib.get_input_format_int(node.ptr, convert(UInt, x-1)),
        num_inputs
    )

    # output processing
    num_outputs = nGraph.get_output_size(node)
    output_sizes = ntuple(x -> nGraph.get_output_shape(node, x), num_outputs)
    output_types = ntuple(x -> nGraph.get_output_element_type(node, x), num_outputs)

    return CPUKernelParams(
        description,
        input_sizes,
        output_sizes,
        input_types,
        output_types,
        ismkl,
        input_formats,
    )
end

## GPU
"""
    GPUKernelParams

Parameters for a GPU Kernel.

$(METHODLIST)

$(FIELDS)
"""
struct GPUKernelParams{IS, OS, IT, OT}
    "The description of the op."
    description::String

    # IO Sizes
    "Sizes of the input tensors."
    input_sizes::IS
    "Sizes of the output tensors."
    output_sizes::OS
    "Element types of the input tensors."
    input_types::IT
    "Element types of the output tensors."
    output_types::OT
end

function GPUKernelParams(node::nGraph.NodeLike)
    description = nGraph.description(node)

    # Input processing
    num_inputs = nGraph.get_input_size(node)
    input_sizes = ntuple(x -> nGraph.get_input_shape(node, x), num_inputs)
    input_types = ntuple(x -> nGraph.get_input_element_type(node, x), num_inputs)

    # output processing
    num_outputs = nGraph.get_output_size(node)
    output_sizes = ntuple(x -> nGraph.get_output_shape(node, x), num_outputs)
    output_types = ntuple(x -> nGraph.get_output_element_type(node, x), num_outputs)

    return GPUKernelParams(
        description,
        input_sizes,
        output_sizes,
        input_types,
        output_types,
    )
end

#####
##### Caches
#####

abstract type AbstractKernelCache end

Base.getindex(cache::AbstractKernelCache, args...) = getindex(cache.cache, args...)
Base.setindex!(cache::AbstractKernelCache, args...) = setindex!(cache.cache, args...)
Base.haskey(cache::AbstractKernelCache, args...) = haskey(cache.cache, args...)
Base.delete!(cache::AbstractKernelCache, args...) = delete!(cache.cache, args...)

# NOTE: Keeping caches as separate types rather than parametric to be backwards compatible
# with older serialized objects.
struct CPUKernelCache <: AbstractKernelCache
    file::String
    cache::Dict{Tuple{CPUKernelParams, IOConfig}, Float64}
end

struct GPUKernelCache <: AbstractKernelCache
    file::String
    cache::Dict{Tuple{GPUKernelParams, IOConfig}, Union{Float64, Vector{CUDNNAlgorithm}}}
end

# Method to determine if we can select the algorithm for this kernel.
can_select_algo(c::GPUKernelCache, p::GPUKernelParams) = _cs(c.cache[p])

CPUKernelCache(file; kw...) = _make_cache(CPUKernelCache, file; kw...)::CPUKernelCache
GPUKernelCache(file; kw...) = _make_cache(GPUKernelCache, file; kw...)::GPUKernelCache

function _make_cache(::Type{T}, file; force_new = false) where {T}
    # If the cache path already exists, just return the existing object.
    # The type assertion for the function will make sure we don't return something weird.
    if (force_new == false) && ispath(file)
        # Deserialize the file.
        cache_db = deserialize(file)

        # Check if we have a key for the current environment context
        # TODO: Don't really need this for the GPU ...
        ctx = _env_context()
        haskey(cache_db, ctx) && return cache_db[ctx]::T
    end

    # Otherwise, create the object.
    return T(
        file,
        last(fieldtypes(T))()
    )
end

function save(cache::AbstractKernelCache)
    # Make the directory for this cache if needed.
    dir = dirname(cache.file)
    ispath(dir) || mkdir(dir)

    # Deserialize the existing cache.
    ctx = _env_context()
    if ispath(cache.file)
        cache_db = deserialize(cache.file)::Dict
        cache_db[ctx] = cache
    else
        cache_db = Dict(ctx => cache)
    end

    # Serialize to a temporary file first
    # Then swap with the original cache.
    #
    # This ensures that if we're interrupted during serialization (which does happen)
    # we don't corrupt the original cache.
    path = tempname()
    open(path; write = true) do io
        serialize(io, cache_db)
    end

    mv(path, cache.file; force = true)
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


max_stars_repo_name: rlerichev/ProjectGitExercise

max_issues_repo_name: rlerichev/ProjectGitExercise

max_forks_repo_name: rlerichev/ProjectGitExercise


text:
# when ############################

when(ta::TimeArray, period::Function, t::Integer) =
    ta[findall(period.(timestamp(ta)) .== t)]
when(ta::TimeArray, period::Function, t::String) =
    ta[findall(period.(timestamp(ta)) .== t)]

# from, to ######################

from(ta::TimeArray{T, N, D}, d::D) where {T, N, D} =
    length(ta) == 0 ? ta :
        d < timestamp(ta)[1] ? ta :
        d > timestamp(ta)[end] ? ta[1:0] :
        ta[searchsortedfirst(timestamp(ta), d):end]

to(ta::TimeArray{T, N, D}, d::D) where {T, N, D} =
    length(ta) == 0 ? ta :
        d < timestamp(ta)[1] ? ta[1:0] :
        d > timestamp(ta)[end] ? ta :
        ta[1:searchsortedlast(timestamp(ta), d)]

###### findall ##################

Base.findall(ta::TimeArray{Bool,1}) = findall(values(ta))

###### findwhen #################

findwhen(ta::TimeArray{Bool,1}) = timestamp(ta)[findall(values(ta))]

###### head, tail ###########

@generated function head(ta::TimeArray{T,N}, n::Int=6) where {T,N}
    new_values = (N == 1) ? :(values(ta)[1:n]) : :(values(ta)[1:n, :])

    quote
        new_timestamp = timestamp(ta)[1:n]
        TimeArray(new_timestamp, $new_values, colnames(ta), meta(ta))
    end
end

 @generated function tail(ta::TimeArray{T,N}, n::Int=6) where {T,N}
    new_values = (N == 1) ? :(values(ta)[start:end]) : :(values(ta)[start:end, :])

    quote
        start = length(ta) - n + 1
        new_timestamp = timestamp(ta)[start:end]
        TimeArray(new_timestamp, $new_values, colnames(ta), meta(ta))
    end
end

###### first, last ###########

Base.first(ta::TimeArray) = head(ta, 1)

Base.last(ta::TimeArray) = tail(ta, 1)


text:
function set_volume(nodes, vecAx, vecAy, vecAz, cellxmax, cellymax, cellzmax)
    volume = zeros(cellxmax, cellymax, cellzmax)
    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                vec_rx = nodes[i+1,j+1,k+1,1] - nodes[i,j,k,1]
                vec_ry = nodes[i+1,j+1,k+1,2] - nodes[i,j,k,2]
                vec_rz = nodes[i+1,j+1,k+1,3] - nodes[i,j,k,3]

                Ax = vecAx[i,j,k,1] + vecAy[i,j,k,1] + vecAz[i,j,k,1]
                Ay = vecAx[i,j,k,2] + vecAy[i,j,k,2] + vecAz[i,j,k,2]
                Az = vecAx[i,j,k,3] + vecAy[i,j,k,3] + vecAz[i,j,k,3]

                volume[i,j,k] = abs(Ax*vec_rx + Ay*vec_ry + Az*vec_rz) /3
            end
        end
    end
    return volume
end

function set_cellcenter(nodes, cellxmax, cellymax, cellzmax)
    cellcenter = zeros(cellxmax, cellymax, cellzmax, 3)
    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:3
                    xl = 0.125*(nodes[i+1,j+1,k+1,l] + nodes[i+1,j,k+1,l] + nodes[i,j+1,k+1,l] + nodes[i,j,k+1,l] +
                                nodes[i+1,j+1,k,l] + nodes[i+1,j,k,l] + nodes[i,j+1,k,l] + nodes[i,j,k,l])

                    cellcenter[i,j,k,l] = xl
                end
            end
        end
    end
    return cellcenter
end

function set_mu(Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio, Rd)
    mu = zeros(cellxmax, cellymax, cellzmax)

    # サザーランドの式
    # https://www.jstage.jst.go.jp/article/jsam1937/37/4/37_4_694/_pdf/-char/ja
    mu0 = 1.82e-5 # 基準粘度[Pa s]
    T0 = 293.15   # 基準温度[K]
    C = 117       # サザーランド定数[K]

    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                # T = p/(rho Rd )
                T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)

                if T <0
                    println("T is munus ( 'A' ) !!")
                    println("i : "*string(i)*"  j : "*string(j)*"  k : "*string(k))
                    println(Qbase[i,j,k,:])
                end

                mu[i,j,k] = mu0 * (T/T0)^1.5 * (T0+C)/(T+C)
            end
        end
    end

    return mu
end

function set_lambda(Qbase, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)

    lambda = zeros(cellxmax, cellymax, cellzmax)

    # サザーランドの式
    lam0 = 22.3*10^(-3)  # 基準熱伝導率　[W/mK]
    T0 = 273.15       # 基準温度[K]
    C = 125       # サザーランド定数[K]

    for i in 1:cellxmax
        for j in 1:cellymax
            for k in 1:cellzmax
                T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)
                lambda[i,j,k] = lam0*((T0+C)/(T+C))*(T/T0)^1.5
            end
        end
    end

    return lambda
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


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

name = "Leptonica"
version = v"1.79.0"

# Collection of sources required to build Leptonica
sources = [
    ArchiveSource("http://www.leptonica.org/source/leptonica-$(version).tar.gz",
                  "045966c9c5d60ebded314a9931007a56d9d2f7a6ac39cb5cc077c816f62300d8"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/leptonica-*/
export CPPFLAGS="-I$prefix/include"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
install_license leptonica-license.txt
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("liblept", :liblept),
    ExecutableProduct("convertfilestopdf", :convertfilestopdf),
    ExecutableProduct("convertfilestops", :convertfilestops),
    ExecutableProduct("convertformat", :convertformat),
    ExecutableProduct("convertsegfilestopdf", :convertsegfilestopdf),
    ExecutableProduct("convertsegfilestops", :convertsegfilestops),
    ExecutableProduct("converttopdf", :converttopdf),
    ExecutableProduct("converttops", :converttops),
    ExecutableProduct("fileinfo", :fileinfo),
    ExecutableProduct("imagetops", :imagetops),
    ExecutableProduct("xtractprotos", :xtractprotos),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Giflib_jll", uuid="59f7168a-df46-5410-90c8-f2779963d0ec"))
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8"))
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828"))
    Dependency(PackageSpec(name="libwebp_jll", uuid="c5f90fcd-3b7e-5836-afba-fc50a0988cb2"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="OpenJpeg_jll", uuid="643b3616-a352-519d-856d-80112ee9badc"))
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
############### Define the solver: DeepPolyBounds ###############################

struct DeepPolyBounds <: Solver end

###### Forward propagation

function forward_linear(solver::DeepPolyBounds, L::LayerNegPosIdx, input::SymbolicIntervalBounds)
    output_Low, output_Up = interval_map(L.W_neg, L.W_pos, input.sym.Low, input.sym.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    return SymbolicIntervalBounds(sym, input.lbs, input.ubs)
end

function relaxed_relu_gradient_lower(l::Real, u::Real)
    ((u <= 0) || ((l < 0) && (u <= -l))) && return 0.
    return 1.
end

""" Old implementation, ~4x slower
function relaxed_relu_gradient_lower(l::Real, u::Real)
    u <= 0 && return 0.
    l >= 0 && return 1.
    u <= -l && return 0
    return 1.
end
"""

function forward_act(solver::DeepPolyBounds, L::LayerNegPosIdx{ReLU}, input::SymbolicIntervalBounds)
    n = n_nodes(L)
    output_Low, output_Up = copy(input.sym.Low), copy(input.sym.Up)
    los, his = bounds(input.sym, input.lbs[L.index], input.ubs[L.index])

    slopes = relaxed_relu_gradient.(los, his)
    output_Up .*= slopes
    output_Up[:, end] .+= slopes .* max.(-los, 0)

    slopes .= relaxed_relu_gradient_lower.(los, his)
    output_Low .*= slopes

    sym =  SymbolicInterval(output_Low, output_Up, domain(input))
    output = SymbolicIntervalBounds(sym, input.lbs, input.ubs)
    output.lbs[L.index] = los
    output.ubs[L.index] = his
    return output
end

function forward_act(solver::DeepPolyBounds, L::LayerNegPosIdx{Id}, input::SymbolicIntervalBounds)
    return input
end


text:
# model
g = 9.81

L1 = 0.5
L2 = 0.5
L3 = 1.5
L4 = L3

l1 = L1/2
l2 = L2/2
l3 = L3/2
l4 = L4/2

m1 = 1.0
m2 = 1
m3 = 1
m4 = 1

m = [m1, m2, m3, m4]

J1 = 1
J2 = 1
J3 = 1
J4 = 1

M1 = Diagonal([m1*ones(2);J1])
M2 = Diagonal([m2*ones(2);J2])
M3 = Diagonal([m3*ones(2);J3])
M4 = Diagonal([m4*ones(2);J4])

M = [M1, M2, M3, M4]

leg = DoggoLeg(L1,L3, [m1,m2,m3,m4], [J1,J2,J3,J4])


text:
module Kernels
# super type
abstract type Kernel end

# diferent kernels to use
type RBFKernel <: Kernel
	params::Array{Real,1}
	RBFKernel(s2, l) = new([s2,l])
	RBFKernel(params) = new(params)
end

type BrownianKernel <: Kernel
	params::Array{Float64,1}
	BrownianKernel(s2::Float64) = new([s2])
end

type OrnsteinUhlenbeckKernel <: Kernel
	params::Array{Float64,1}
	OrnsteinUhlenbeckKernel(s2::Float64, l::Float64) = new([s2, l])
end

type WhiteNoiseKernel <: Kernel
	params::Array{Float64,1}
	WhiteNoiseKernel(s2) = new([s2])
end

type PeriodicRBFKernel <: Kernel
	params::Vector
	PeriodicRBFKernel(s2::Float64, l::Float64, p::Float64) = new([s2,l,p])
end

type CauchyKernel <: Kernel
	params::Vector
	CauchyKernel(s2::Float64, scale::Float64) = new([s2, scale])
end

type RationalQuadraticKernel <: Kernel
	params::Vector
	RationalQuadraticKernel(s2::Float64, l::Float64, alpha::Float64) = new([s2, l, alpha])
end

type ConstantKernel <: Kernel
	params::Vector
	ConstantKernel(s2::Float64) = new([s2])
end

export Kernel,
	RBFKernel,
	BrownianKernel,
	OrnsteinUhlenbeckKernel,
	WhiteNoiseKernel,
	PeriodicRBFKernel,
	CauchyKernel,
	computeCovarianceMatrix,
	computeCrossCovarianceMatrix

# covariance function for various kernels
# rbf
computeCovariance(k::RBFKernel, x::Array{Float64, 1}, x_prime::Array{Float64, 1}) = k.params[1] * exp(-.5*norm(x-x_prime)^2/(k.params[2]^2))

# brownian motion
computeCovariance(k::BrownianKernel, x::Array{Float64, 1}, x_prime::Array{Float64, 1}) = float((k.params[1] * min(x, x_prime))[1])

# ornstein-uhlenbeck
computeCovariance(k::OrnsteinUhlenbeckKernel, x::Array{Float64, 1}, x_prime::Array{Float64, 1}) = k.params[1] * exp(-norm(x-x_prime)/k.params[2])

# white noise
computeCovariance(k::WhiteNoiseKernel, x::Array{Float64, 1}, x_prime::Array{Float64, 1}) = x==x_prime ? k.params[1] : 0

#periodic
computeCovariance(k::PeriodicRBFKernel, x::Vector, x_prime::Vector) = k.params[1] * exp(-2*sin(norm(x-x_prime)/k.params[3])^2/(k.params[2]^2) )

#cauchy
computeCovariance(k::CauchyKernel, x::Vector, x_prime::Vector) = k.params[1]/(1 + (norm(x-x_prime)/k.params[2])^2)

#rational quadratic
computeCovariance(k::RationalQuadraticKernel, x::Vector, x_prime::Vector) = k.params[1]*(1 + (norm(x-x_prime)^2 / (2*k.params[3]*k.params[2]^2)))^(-k.params[3])

#constant
computeCovariance(k::ConstantKernel, x::Vector, x_prime::Vector) = k.params[1]

# to_string(k::Kernel) =

function setParameters{T<:Real}(k::Kernel, newParams::Array{T, 1})
	k.params[:] = newParams[:]
end

#compute K_xx
function computeCovarianceMatrix{T}(k::Kernel, points::Array{T})
	n = size(points)[1]
	cov = zeros(n,n)
	for i in 1:n, j in i:n
		cov[i,j] = computeCovariance(k, vec(points[i,:]), vec(points[j,:]))
		cov[j,i] = cov[i,j]
	end
	cov
end

#compute K_xy
function computeCrossCovarianceMatrix{T}(k::Kernel, x::Array{T}, x_prime::Array{T})
	n = size(x)[1]
	m = size(x_prime)[1]
	cov = zeros(n,m)
	for i in 1:n, j in 1:m
		cov[i,j] = computeCovariance(k, vec(x[i,:]), vec(x_prime[j,:]))
	end
	cov
end

end


text:
struct VectorBackedUTF8String <: AbstractString
    buffer::Vector{UInt8}
end

Base.:(==)(x::VectorBackedUTF8String, y::VectorBackedUTF8String) = x.buffer == y.buffer

function Base.show(io::IO, x::VectorBackedUTF8String)
    print(io, '"')
    print(io, string(x))
    print(io, '"')
    return
end

Base.pointer(s::VectorBackedUTF8String) = pointer(s.buffer)

Base.pointer(s::VectorBackedUTF8String, i::Integer) = pointer(s.buffer) + i - 1

Base.pointer(s::SubString{VectorBackedUTF8String}, i::Integer) = pointer(s.string) + s.offset + i - 1

@inline Base.ncodeunits(s::VectorBackedUTF8String) = length(s.buffer)

Base.codeunit(s::VectorBackedUTF8String) = UInt8

@inline function Base.codeunit(s::VectorBackedUTF8String, i::Integer)
    @boundscheck checkbounds(s.buffer, i)
    return @inbounds s.buffer[i]
end

Base.thisind(s::VectorBackedUTF8String, i::Int) = Base._thisind_str(s, i)

Base.nextind(s::VectorBackedUTF8String, i::Int) = Base._nextind_str(s, i)

Base.isvalid(s::VectorBackedUTF8String, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i

Base.@propagate_inbounds function Base.iterate(s::VectorBackedUTF8String, i::Int=firstindex(s))
    i > ncodeunits(s) && return nothing
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return our_next_continued(s, i, u)
end

function our_next_continued(s::VectorBackedUTF8String, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

# The following functions all had implementations in WeakRefStrings (which was taken
# as a starting point for the code in this file), but aren't needed for the
# TextParse.jl use. For now we leave stubs that throw errors around. If this type
# turns out to be useful beyond TextParse.jl, these should be implemented properly.

Base.:(==)(x::String, y::VectorBackedUTF8String) = error("Not yet implemented.")

Base.:(==)(y::VectorBackedUTF8String, x::String) = x == y

Base.hash(s::VectorBackedUTF8String, h::UInt) = error("Not yet implemented.")

Base.print(io::IO, s::VectorBackedUTF8String) = error("Not yet implemented.")

Base.textwidth(s::VectorBackedUTF8String) = error("Not yet implemented.")

Base.string(x::VectorBackedUTF8String) = unsafe_string(pointer(x.buffer), length(x.buffer))

Base.convert(::Type{VectorBackedUTF8String}, x::String) = error("Not yet implemented.")

Base.convert(::Type{String}, x::VectorBackedUTF8String) = error("Not yet implemented.")

Base.String(x::VectorBackedUTF8String) = error("Not yet implemented.")

Base.Symbol(x::VectorBackedUTF8String) = error("Not yet implemented.")


text:
module MultiFieldCellFieldsTests

using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.ReferenceFEs
using Gridap.Fields
using Gridap.Integration
using Gridap.CellData
using Gridap.MultiField
using Gridap.TensorValues
using Test

using Gridap.Arrays: BlockArrayCooMap

domain = (0,1,0,1)
cells = (2,2)
model = CartesianDiscreteModel(domain,cells)

trian = Triangulation(model)

u1(x) = sin(x[1])
cf1 = CellField(u1,trian)

u2(x) = cos(x[2])
cf2 = CellField(u2,trian)

cf = MultiFieldCellField([cf1,cf2])

@test cf1 === cf[1]
@test cf2 === cf[2]

_cf1, _cf2 = cf

@test cf1 === _cf1
@test cf2 === _cf2

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

degree = order

trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)
x = get_cell_points(quad)

trian_Γ = SkeletonTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,degree)
x_Γ = get_cell_points(quad_Γ)

V = TestFESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},order);conformity=:H1)
Q = TestFESpace(model,ReferenceFE(lagrangian,Float64,order-1),conformity=:L2)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

dv, dq = get_cell_shapefuns(Y)
du, dp = get_cell_shapefuns_trial(X)

n = VectorValue(1,2)

cellmat = integrate( (n⋅dv)*dp + dq*dp, quad)
cellvec = integrate( n⋅dv, quad)
@test isa(cellvec,LazyArray{<:Fill{BlockArrayCooMap{1}}})
@test isa(cellmat,LazyArray{<:Fill{BlockArrayCooMap{2}}})

cellmat1 = integrate( ((n⋅dv) - dq)*((n⋅du) + dp), quad)
cellmat2 = integrate( (n⋅dv)*(n⋅du) + (n⋅dv)*dp - dq*(n⋅du) - dq*dp, quad)
test_array(cellmat1,cellmat2,≈)

cellmat1 = integrate( (n⋅dv)*2, quad)
cellmat2 = integrate( (n⋅dv)*fill(2,num_cells(trian)), quad)
test_array(cellmat1,cellmat2,≈)

α = CellField(2,trian)
op(u,∇u,v,∇v,α) = α*(u⋅v) + ∇u⊙∇v
cellmat1 = integrate( op∘(du,∇(du),dv,∇(dv),α) , quad)
cellmat2 = integrate( α*(du⋅dv) + ∇(du)⊙∇(dv) , quad)
test_array(cellmat1,cellmat2,≈)

@test isa(cellmat2,LazyArray{<:Fill{BlockArrayCooMap{2}}})
#@test_broken isa(cellmat1,LazyArray{<:Fill{BlockArrayCooMap{2}}})

α = CellField(2,trian)
op2(u,∇u,α) = α*(∇u⋅u)
cellmat1 = integrate( dv⋅(op2∘(du,∇(du),α)),quad)
cellmat2 = integrate( dv⋅(α*(∇(du)⋅du)),quad)
test_array(cellmat1,cellmat2,≈)
@test isa(cellmat2,LazyArray{<:Fill{BlockArrayCooMap{2}}})
#@test_broken isa(cellmat1,LazyArray{<:Fill{BlockArrayCooMap{2}}})

conv(u,∇u,α) = α*(u⋅∇u)
dconv(du,∇du,u,∇u,α) = conv(u,∇du,α)+conv(du,∇u,α)

u = zero(U)
cellvec2 = integrate(dv⊙(α*(u⋅∇(u))),quad)
cellvec1 = integrate(dv⊙(conv∘(u,∇(u),α)),quad)
test_array(cellvec1,cellvec2,≈)
@test isa(cellvec2,LazyArray{<:Fill{BlockArrayCooMap{1}}})
@test isa(cellvec1,LazyArray{<:Fill{BlockArrayCooMap{1}}})

cellmat1 = integrate( dv⋅(dconv∘(du,∇(du),u,∇(u),α)) , quad)
cellmat2 = integrate( dv⋅( α*(du⋅∇(u)) + α*(u⋅∇(du))), quad)
test_array(cellmat1,cellmat2,≈)
@test isa(cellmat2,LazyArray{<:Fill{BlockArrayCooMap{2}}})
#@test_broken isa(cellmat1,LazyArray{<:Fill{BlockArrayCooMap{2}}})

cellmat_Γ = integrate(  jump(n⋅dv)*dp.⁺ + mean(dq)*jump(dp), quad_Γ)
cellvec_Γ = integrate(  jump(n⋅dv) + mean(dq), quad_Γ)
L = 1
R = 2
@test isa(cellmat_Γ,LazyArray{<:Fill{BlockArrayCooMap{2}}})
@test isa(cellvec_Γ,LazyArray{<:Fill{BlockArrayCooMap{1}}})

cell = 1
@test isa(cellmat_Γ[cell][Block(L,R)],BlockArrayCoo)
@test isa(cellvec_Γ[cell][Block(L)],BlockArrayCoo)

cellmat1_Γ = integrate(((n⋅dv.⁺)-dq.⁻)*((n⋅du.⁺)+dp.⁻),quad_Γ)
cellmat2_Γ = integrate((n⋅dv.⁺)*(n⋅du.⁺)+(n⋅dv.⁺)*dp.⁻-dq.⁻*(n⋅du.⁺)-dq.⁻*dp.⁻,quad_Γ)
test_array(cellmat1_Γ,cellmat2_Γ,≈)

#a = cellmat_Γ
#using BenchmarkTools
#cache = array_cache(a)
#@btime getindex!($cache,$a,2)

end # module


text:
using Pkg
Pkg.update()
Pkg.add("Plots")
Pkg.add("CSV")
Pkg.add("DataFrames")


text:
#dof_contain.jl
function zmatrix(lmm::LMM, i)
    l   = zeros(Int, lmm.covstr.sn[i])'
    rzm = Matrix{eltype(lmm.covstr.z)}(undef, nobs(lmm), length(lmm.covstr.zrndur[i]) * lmm.covstr.sn[i])
    sn = 1
    for b = 1:length(lmm.covstr.vcovblock)
        zblock    = view(lmm.covstr.z, lmm.covstr.vcovblock[b], lmm.covstr.zrndur[i])
        for s = 1:length(lmm.covstr.sblock[b][i])
            zi    = view(zblock, lmm.covstr.sblock[b][i][s], :)
            l[sn] = 1
            copyto!(view(rzm, lmm.covstr.vcovblock[b][lmm.covstr.sblock[b][i][s]], :), kron(l, zi))
            l[sn] = 0
            sn   += 1
        end
    end
    rzm
end
function fullzmatrix(lmm)
    fzm = Matrix{Int}(undef, nobs(lmm), 0)
    for r = 1:length(lmm.covstr.random)
        rzm = zmatrix(lmm, r)
        fzm = hcat(fzm, rzm)
    end
    fzm
end

function rankxz(lmm::LMM)
    rank(hcat(lmm.data.xv, fullzmatrix(lmm)))
end
function rankxz(lmm::LMM, i)
    rank(hcat(lmm.data.xv, zmatrix(lmm, i)))
end

"""
    dof_contain(lmm, i)

!!! warning
    Experimental!
    Compute rank(XZi) for each random effect that syntactically contain factor assigned for β[i] element (Where Zi - Z matrix for random effect i).
    Minimum returned. If no random effect found  N - rank(XZ) returned.
"""
function dof_contain(lmm, i)
    ind  = lmm.mm.assign[i]
    sym  = get_symb(lmm.mf.f.rhs.terms[ind])
    rr   = Vector{Int}(undef, 0)
    for r = 1:length(lmm.covstr.random)
        if length(intersect(sym, get_symb(lmm.covstr.random[r].model))) > 0
            push!(rr, rankxz(lmm, r))
        end
    end
    if length(rr) > 0
        return minimum(rr)
    else
        return nobs(lmm) - rankxz(lmm)
    end
end

function dof_contain(lmm)
    dof   = zeros(Int, length(lmm.mm.assign))
    rrt   = zeros(Int, length(lmm.covstr.random))
    rz   = 0
    for i = 1:length(lmm.mm.assign)
        ind  = lmm.mm.assign[i]
        sym  = get_symb(lmm.mf.f.rhs.terms[ind])
        rr   = Vector{Int}(undef, 0)
        for r = 1:length(lmm.covstr.random)
            if length(intersect(sym, get_symb(lmm.covstr.random[r].model))) > 0
                if rrt[r] == 0 rrt[r] = rankxz(lmm, r) end
                push!(rr, rrt[r])
            end
        end
        if length(rr) > 0
            dof[i] = minimum(rr)
        else
            if rz == 0 rz = nobs(lmm) - rankxz(lmm) end
            dof[i] = rz
        end
    end
    dof
end
"""
    dof_contain_f(lmm, i)

!!! warning
    Experimental

"""
function dof_contain_f(lmm, i)
    sym  = get_symb(lmm.mf.f.rhs.terms[i])
    rr   = Vector{Int}(undef, 0)
    for r = 1:length(lmm.covstr.random)
        if length(intersect(sym, get_symb(lmm.covstr.random[r].model))) > 0
            push!(rr, rankxz(lmm, r))
        end
    end
    if length(rr) > 0
        return minimum(rr)
    else
        return nobs(lmm) - rankxz(lmm)
    end
end


text:
export clip_by_global_norm!, global_norm

using Zygote

Zygote.@adjoint argmax(xs; dims = :) = argmax(xs; dims = dims), _ -> nothing

global_norm(gs::Zygote.Grads, ps::Zygote.Params) =
    sqrt(sum(mapreduce(x -> x^2, +, gs[p]) for p in ps))

function clip_by_global_norm!(gs::Zygote.Grads, ps::Zygote.Params, clip_norm::Float32)
    gn = global_norm(gs, ps)
    if clip_norm <= gn
        for p in ps
            gs[p] .*= clip_norm / max(clip_norm, gn)
        end
    end
    gn
end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8, 6))
hist(data.Solar_gen, color="blue", edgecolor="white", bins=100,
     density=true, alpha=0.5)
xlabel("Solar power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)
xticks(fontsize=12)
yticks([0, 1e-3, 2e-3], fontsize=12)
savefig(joinpath(@OUTPUT, "hist_solar.svg")) # hide


text:
module EventStreamGLM

import GLM

import Base: rand
import GLM: BlasReal, DensePred, delbeta!, cholesky!, cholesky, LinPred, invchol, linpred!
import EventStreamMatrix: AbstractEventStreamMatrix, AbstractEventStreamVector

export DensePredConjGrad
export delbeta!
export cholesky!, cholesky
export CGControl
export EventStreamPredConjGrad
export rand
export greet
export EventStreamProcess
export other_intensity
export HomogPoissonProcess
export PoissonProcess
export HawkesProcess

using EventStreamMatrix
using GLM
using IterativeSolvers
using LinearAlgebra
using StatsFuns
using BSplines

greet() = "Damns"

include("ConjGradLinPred.jl")
include("simulation.jl")

end


text:
#=
- Julia's arrays default to column vector behaviour
- Julia's indexes start at 1 and you can access elements using brackets
- Slices are inclusive of endpoints, and `end` means the last component.
- Steps can be informed in the middle
=#
x = [1, 2, 3, 4, 5]
x * 2
ones(5) + x
x[2]
x[1:2]
x[2:end]
x[1:2:end]
length(x)
# - Matrices are first-class citizens.
# - Create a matrix with brackets, spaces to change columns, and `;` to change lines.
# - Only in vectors `,` can be used instead of spaces (compat. Python).
# - You can concatenate matrices, vectors, and numbers with the same syntax.
A = [3 2 4; -1 0 1]
A * ones(3)
B = rand(3, 3)
A * B
B * A'
A[1:2, 1:2]
[
  B A';
  A ones(2,2)
]
# - Linear systems can be solved with `\ `
# - The LinearAlgebra stdlib module contains useful linear algebra tools, functions, and elements.
A = randn(50, 50)
b = A * ones(50)
x = A \ b
using LinearAlgebra
norm(b - A * x)
norm(x .- 1)
# - Everything has a concrete type, every type derives from some abstract type.
# - There is no implicit conversion, but sometimes it looks like there is.
typeof(2)
typeof(2.0)
Float64 |> supertype
Float64 |> supertype |> supertype
typeof(2 + 2.0)
typeof([2, 2.0])
typeof([2, true])
typeof([2, "a"])
# - There are 3 approaches to write functions
function foo(a, b)
  return (a + 1) * (b - 1)
end
foo(2, 3)
bar(x) = x^2 - 2
bar(5)
anon = (a, b, c) -> sqrt(a^2 + b^2 + c^2)
anon(-1, 0, 1)
# - `if CONDITION ... elseif ... else ... end`
# - `while CONDITION ... end`
# - `for VARIABLE in SET|RANGE|ITERATOR ... end`
# - You can specify types
function myfactorial(n::Int)
  if n < 0
    error("Nay nay")
  end
  prod = 1
  for i = 2:n
    prod *= i
  end
  return prod
end
myfactorial(4)
# - You can insert unicode characters with `\something[TAB]`.
function machine_eps(x::AbstractFloat)
  ϵ = x
  while x + ϵ != x
    ϵ /= 2
  end
  return 2ϵ
end
machine_eps(1.0), machine_eps(1e4), machine_eps(Float16(1.0)), machine_eps(1.0f0)
# - you can create a different signature
function machine_eps(x::Integer)
  return zero(x)
end
machine_eps(1), machine_eps(0x05), machine_eps(typemax(1))
# - you can set default positional arguments and keyword arguments
function newton(f, fder, x = 0.0; atol = 1e-6, rtol = 1e-6, max_iter = 1000)
  fx = f(x)
  iter = 0
  ϵ = atol + rtol * abs(fx)
  solved = abs(fx) < ϵ
  tired = iter ≥ max_iter
  while !(solved || tired)
    slope = fder(x)
    if abs(slope) < 1e-12
      error("0 derivative at x = $x")
    end
    x -= fx / slope
    fx = f(x)
    iter += 1
    solved = abs(fx) < ϵ
    tired = iter ≥ max_iter
  end
  return x, fx, iter
end
newton(x -> x - exp(-x), x -> 1 + exp(-x), 1.0)
# - You can use parametric types
# - The type must be inferred from somewhere, so `x` is not optional anymore
function newton_parametric(f, fder, x :: T; atol :: T = √eps(T), rtol :: T = √eps(T), max_iter :: Int = 1000) where T
  fx = f(x)
  iter = 0
  ϵ = atol + rtol * abs(fx)
  slope_ϵ = (√√eps(T))^3 # ϵₘ^(3/4)
  solved = abs(fx) < ϵ
  tired = iter ≥ max_iter
  while !(solved || tired)
    slope = fder(x)
    if abs(slope) < slope_ϵ
      error("0 derivative at x = $x")
    end
    x -= fx / slope
    fx = f(x)
    iter += 1
    solved = abs(fx) < ϵ
    tired = iter ≥ max_iter
  end
  return x, fx, iter
end
for T in [Float16, Float32, Float64, BigFloat]
  x, fx, iter = newton_parametric(x -> x^2 - 2, x -> 2x, one(T))
  println("√2 ≈ $x")
end
# - Formatted printing is possible with the `Printf` stdlib
# - `@printf` is a macro. Interpolation (`"$x"`) is not possible.
using Printf
for T in [Float16, Float32, Float64, BigFloat]
  x, fx, iter = newton_parametric(x -> x^2 - 2, x -> 2x, one(T))
  @printf("√2 ≈ %20.16e, (√2)² - 2 ≈ %20.16e\n", x, x^2 - 2)
end
# - In the REPL, `?` access the help mode.
# - `apropos` can be used to search among the docstrings
#=
---

## Exercises

1. Given `v` and `w` vectors, verify that `[v w]` and `[v, w]` are different. Notice the type of the second object.
1. Search for "dot product" using the `apropos` function.
1. Create a function that orthogonalizes a _list of vectors_ following the given theorem:
**Theorem**: Given $n$ linear independent vector $v_1, v_2, \dots, v_n$, define
$\displaystyle q_i = v_i - \sum_{j = 1}^{i-1} (v_i \cdot q_j) q_j, \ \forall \ i = 1,\dots,n$.
The vectors $q_1, q_2, \dots, q_n$ are orthogonal.
=#

max_stars_repo_name: abelsiqueira/julia-tutorial-at-nlesc

max_issues_repo_name: abelsiqueira/julia-tutorial-at-nlesc

max_forks_repo_name: abelsiqueira/julia-tutorial-at-nlesc


text:
#convertposit.jl - posit conversions for FastSigmoid

function __UT{T <: Posit}(::Type{T})
  (T == Posit) && throw(ArgumentError("__UT cannot operate on general posit"))
  S = sizeof(T) * 8
  Dict(8 => UInt8, 16 => UInt16, 32 => UInt32, 64 => UInt64)[S]
end
function __ST{T <: Posit}(::Type{T})
  (T == Posit) && throw(ArgumentError("__ST cannot operate on general posit"))
  S = sizeof(T) * 8
  Dict(8 => Int8, 16 => Int16, 32 => Int32, 64 => Int64)[S]
end
function __PT{T <: Posit}(::Type{T})
  isbits(T) && return T
  (T == Posit) && throw(ArgumentError("__PT cannot operate on general posit"))
  S = sizeof(T) * 8
  (Posit{S} <: T) && throw(ArgumentError("__PT requires a type with an ES"))
  __es{N,ES}(::Type{Posit{N,ES}}) = ES
  es = __es(T)
  Dict(8 => Posit8{es}, 16 => Posit16{es}, 32 => Posit32{es}, 64 => Posit64{es})[S]
end

#convenience functions to convert to signed and unsigned representations.
@generated function __u{N,ES}(x::Posit{N,ES})
  UT = __UT(Posit{N,ES})
  :(reinterpret($UT, x))
end
@generated function __s{N,ES}(x::Posit{N,ES})
  ST = __ST(Posit{N,ES})
  :(reinterpret($ST, x))
end

doc"""
  `FastSigmoid.__decomp(x)`
  takes a fastsigmoid number and breaks it up into three parts:
  1) sign (bool)
  2) exponent (int, same as backing int type)
  3) fraction (uint)
"""
@generated function __decomp{N,ES}(x::Posit{N,ES})
  BackingFloat = Dict(8 => Float32, 16 => Float32, 32 => Float64)[N]
  BackingSInt  = Dict(8 => Int32, 16 => Int32, 32 => Int64)[N]
  BackingUInt  = Dict(8 => UInt32, 16 => UInt32, 32 => UInt64)[N]
  IntType = Dict(8 => UInt8, 16 => UInt16, 32 => UInt32)[N]
  topbit = one(IntType) << (N - 1)
  invbit = one(IntType) << (N - 2)
  #exp_mask done in Int64
  exp_mask = BackingSInt((1 << ES) - 1)
  #size delta.
  sizedelta = sizeof(BackingUInt) * 8 - N
  quote
    intdecomp = reinterpret($IntType, x)
    sign = false
    #decide the sign by checking if the topbit is zero.
    if (intdecomp & ($topbit) != 0)
      intdecomp = -intdecomp
      sign = true
    end
    #set up some local values.
    local fraction::($BackingUInt)
    local shift::Int64
    local regime::Int64
    #ascertain the shift and regime.
    if (intdecomp & ($invbit) != 0)
      #note that shift and regime are going to have Int64 data type.
      shift = leading_ones(intdecomp | $topbit) + 1
      regime = shift - 3
    else
      shift = leading_zeros(intdecomp) + 1
      regime = -shift + 2
    end
    #exp_mask is going to have the correct size.
    exponent = $exp_mask & (intdecomp >> (N - shift - ES))
    exponent |= ($BackingSInt(regime) << ES)

    #shift over the fraction, boss.
    fraction = ($BackingUInt(intdecomp) << (shift + ES + $sizedelta))

    (sign, exponent, fraction)
  end
end

doc"""
  `FastSigmoid.__p2f(x)`
  converts a posit(x) to the backing float type.
"""
@generated function __p2f{N,ES}(x::Posit{N,ES})

  BackingFloat = Dict(8 => Float32, 16 => Float32, 32 => Float64)[N]
  BackingUInt  = Dict(8 => UInt32, 16 => UInt32, 32 => UInt64)[N]

  #need to implement promotion to higher types if the ES requires it.
  fracshift = Dict(8 => 9, 16 => 9, 32 => 12)[N]
  bias      = Dict(8 => 127, 16 => 127, 32 => 1023)[N]
  expshift  = Dict(8 => 23, 16 => 23, 32 => 52)[N]
  signbit   = Dict(8 => 0x8000_0000, 16 => 0x8000_0000, 32 => 0x8000_0000_0000_0000)[N]

  quote
    #check for the dropout values
    if iszeroinf(x)
      if __u(x) == 0
        return zero($BackingFloat)
      else
        return $BackingFloat(Inf)
      end
    end

    #next pull parts out of the posit value
    (sign, exponent, fraction) = __decomp(x)
    #put them together.
    intresult::($BackingUInt) = fraction >> $fracshift
    intresult |= (exponent + $bias) << $expshift
    #set the sign.
    intresult |= ($signbit) * sign
    reinterpret($BackingFloat, intresult)
  end
end

typealias IEEEFloat Union{Float16, Float32, Float64}

#for conversions to other IEEE float types, don't bother directly converting.
function Base.convert{T <: IEEEFloat, U <: Posit}(::Type{T}, x::U)
  T(__p2f(x))
end

#for conversions back to a posit type, go ahead and do this:
@generated function Base.convert{N, ES, T <: IEEEFloat}(::Type{Posit{N,ES}}, x::T)

  fullsize = sizeof(x) * 8
  fracsize = Dict(Float16 => 10, Float32 => 23, Float64 => 52)[T]
  IntType = Dict(Float16 => UInt16, Float32 => UInt32, Float64 => UInt64)[T]
  SIntType = Dict(Float16 => Int16, Float32 => Int32, Float64 => Int64)[T]
  bias = SIntType((1 << (fullsize - fracsize - 2)) - 1)
  topbit = one(IntType) << (fullsize - 1)
  topmask = ~topbit
  PositType = __PT(Posit{N,ES})
  pinf = one(__UT(Posit{N, ES})) << (N - 1)

  quote
    intform = reinterpret($IntType, x)
    if isnan(x)
      throw(NaNError())
    elseif isinf(x)
      return reinterpret($PositType, $pinf)
    elseif x == zero(T)
      return zero($PositType)
    end

    sign = ((intform & $topbit) != 0)

    #signed exponent value.
    exponent = ($SIntType((intform & $topmask) >> $fracsize)) - $bias

    #right-shifted fraction value.
    fraction = intform << ($fullsize - $fracsize)

    __build(Posit{N,ES}, sign, exponent, fraction)
  end
end

function regime_exponent{T <: Signed,ES}(exponent::T, ::Type{Val{ES}})
  low_mask = T((1 << ES) - 1)
  (exponent >> ES, unsigned(exponent & low_mask))
end

@generated function __build{N, ES}(::Type{Posit{N, ES}}, sign::Bool, exponent::Signed, fraction::Unsigned)

  sizeof(exponent) == sizeof(fraction) || throw(ArgumentError("mismatched integer sizes"))

  #in the __build function, we will promote the fraction to be the size of larger
  #of fraction / UInt(N).  This will be trimmed later.

  frac_width = sizeof(fraction) * 8

  if (frac_width) > N
    #note that in the unquoted parts of generated functions the variable name is
    #assigned to the *type* of the variable.
    UIntType = fraction
    IntType  = exponent

    topbit = one(UIntType) << (frac_width - 1)
    minreal = one(UIntType) << (frac_width - N)
    buffer_width = frac_width
  else
    UIntType = __UT(Posit{N, ES})
    IntType  = __ST(Posit{N, ES})

    topbit = one(UIntType) << (N - 1)
    minreal = one(UIntType)
    buffer_width = N
  end

  #next create code which shifts the fraction correctly.  Fraction
  #should start 3 bits to the right of the desired value, plus making some
  #room for the exponent bits.
  create_data_buffer = if frac_width < N
    :( data_buffer = $UIntType(fraction) << (N - frac_width - 3 - ES) )
  else
    :( data_buffer = fraction >> (3 + ES) )
  end

  #the code for summary mask is different for the two cases as well.
  create_summary_mask = if frac_width > N
    :( value_mask -= one($IntType) )
  else
    :( value_mask = (shift > 0) * value_mask - one($IntType) )
  end

  #finally create "trimmer" code which trims the tail end of the fraction.
  trim_buffer_to_result = if frac_width > N
    ResultType = __UT(Posit{N,ES})
    :( result_value = $ResultType(data_buffer >> $(frac_width - N)) )
  else
    :( result_value = data_buffer )
  end

  #create a series of useful constants.
  toptwo = -(topbit >> 1)
  thirdbit = topbit >> 2
  topmask = ~topbit
  maxreal = topbit - minreal

  #this allows us to directly cast to the non-abstract posit type
  PositType = __PT(Posit{N,ES})

  quote
    #adjust the integers to be left-aligned.
    $create_data_buffer

    (regime, subexponent) = regime_exponent(exponent, Val{ES})

    #generate a "signature" - 001 for less than one, 110 for more than one.
    if exponent < 0
      signature = $thirdbit | ($UIntType(subexponent) << ($buffer_width - ES - 3))
      shift = (-regime - 1)
    else
      signature = $toptwo | ($UIntType(subexponent) << ($buffer_width - ES - 3))
      shift = regime
    end

    data_buffer = signature | data_buffer

    #check the inner bit.
    value_mask = $minreal << shift
    innerbit = ((value_mask) & data_buffer) != 0

    #check the guard bit
    value_mask >>= 1
    guardbit = ((value_mask) & data_buffer) != 0

    #check the summary bit.
    $create_summary_mask
    summarybit = ((value_mask) & data_buffer) != 0

    #do an arithmetic shift on the signed data buffer.
    data_buffer = unsigned(signed(data_buffer) >> shift)

    #mask out the top bit.
    data_buffer &= $topmask
    #perform rounding.
    data_buffer += ((summarybit & guardbit) | (innerbit & guardbit)) * $minreal

    #intercept infinities and convert to maxreal.
    data_buffer = clamp(data_buffer, $minreal, $maxreal)

    $trim_buffer_to_result

    reinterpret($PositType, sign ? -result_value : result_value)
  end
end

#reimplement the reinterpret directive to support non-leaf posit types.
Base.reinterpret{ES}(::Type{Posit{8,ES}}, i::UInt8) = reinterpret(Posit8{ES}, i)
Base.reinterpret{ES}(::Type{Posit{16,ES}}, i::UInt16) = reinterpret(Posit16{ES}, i)
Base.reinterpret{ES}(::Type{Posit{32,ES}}, i::UInt32) = reinterpret(Posit32{ES}, i)
Base.reinterpret{ES}(::Type{Posit{64,ES}}, i::UInt64) = reinterpret(Posit64{ES}, i)


text:
# Generic evaluation arguments

@autohasheq struct Argument{T} <: Evaluable{T}
    expression :: Union{Symbol, Expr}
    isconstant :: Bool
    iselconstant :: Bool
end

Base.show(io::IO, self::Argument) = print(io, "Argument($(self.expression))")
isconstant(self::Argument) = self.isconstant
iselconstant(self::Argument) = self.iselconstant

arguments(::Argument) where T = ()
@destructure prealloc(::Argument, []) = []
@destructure codegen(self::Argument, [], []) = self.expression

@autohasheq struct ArrayArgument{T,N} <: ArrayEvaluable{T,N}
    expression :: Union{Symbol, Expr}
    isconstant :: Bool
    iselconstant :: Bool
    shape :: Shape
end

isconstant(self::ArrayArgument) = self.isconstant
iselconstant(self::ArrayArgument) = self.iselconstant
Base.size(self::ArrayArgument) = self.shape

arguments(::ArrayArgument) where T = ()
@destructure prealloc(::ArrayArgument, []) = []
@destructure codegen(self::ArrayArgument, [], []) = self.expression

# ApplyTransform

@autohasheq struct ApplyTransform <: ArrayEvaluable{Float64,1}
    trans :: Evaluable{TransformChain}
    arg :: Evaluable{Vector{Float64}}
    dims :: Int
end

Base.size(self::ApplyTransform) = (self.dims,)

arguments(self::ApplyTransform) = (self.trans, self.arg)
@destructure prealloc(self::ApplyTransform, []) = [
    :(Vector{Float64}(undef, $(self.dims))),
    :(Vector{Float64}(undef, $(self.dims))),
]
@destructure codegen(self::ApplyTransform, [trans, arg], [ws, out]) = quote
    applytrans($trans, $arg, $ws, $out) :: $(restype(self))
end

# ApplyTransformGrad

@autohasheq struct ApplyTransformGrad <: ArrayEvaluable{Float64,2}
    trans :: Evaluable{TransformChain}
    arg :: Evaluable{Vector{Float64}}
    dims :: Int
end

Base.size(self::ApplyTransformGrad) = (self.dims, self.dims)

arguments(self::ApplyTransformGrad) = (self.trans, self.arg)
@destructure prealloc(self::ApplyTransformGrad, []) = [
    :(Vector{Float64}(undef, $(self.dims))),
    :(Vector{Float64}(undef, $(self.dims))),
    :(Matrix{Float64}(undef, $(self.dims), $(self.dims))),
    :(Matrix{Float64}(undef, $(self.dims), $(self.dims))),
]
@destructure codegen(self::ApplyTransformGrad, [trans, arg], [ptws, ptout, mxws, mxout]) = quote
    applytrans_grad($trans, $arg, $ptws, $ptout, $mxws, $mxout) :: $(restype(self))
end

# Constant

@autohasheq struct Constant{T,N} <: ArrayEvaluable{T,N}
    value :: Array{T,N}
end

# Scalars are wrapped as zero-dimensional arrays
Constant(v::T) where T <: Number = Constant(fill(v, ()))

Base.size(self::Constant) = size(self.value)

arguments(::Constant) = ()
@destructure prealloc(self::Constant, []) = [self.value]
@destructure codegen(::Constant, [], [alloc]) = alloc

# Contract

@autohasheq struct Contract{T,N} <: ArrayEvaluable{T,N}
    terms :: Vector{ArrayEvaluable}
    indices :: Vector{Vector{Int}}
    target :: Vector{Int}
    shape :: Shape

    function Contract(terms::Vector{ArrayEvaluable}, indices::Vector{Vector{Int}}, target::Vector{Int})
        @assert length(terms) == length(indices)
        @assert all(length(ind) == ndims(term) for (ind, term) in zip(indices, terms))

        # Enumerate all index occurences and sizes
        occurences = Dict{Int,Vector{Int}}()
        for (term, inds) in zip(terms, indices)
            for (axid, sz) in zip(inds, size(term))
                push!(get!(occurences, axid, Int[]), sz)
            end
        end

        # Every occurence has consistent size, and appears exactly twice, or once + in target
        for (axid, occ) in occurences
            @assert length(occ) == 2 || length(occ) == 1 && axid in target
            @assert all(occ[1] == sz for sz in occ)
        end

        newshape = Tuple(occurences[axid][1] for axid in target)
        newtype = reduce(promote_type, (eltype(term) for term in terms))
        newdims = length(target)
        new{newtype, newdims}(ArrayEvaluable[Normalize(term) for term in terms], indices, target, newshape)
    end
end

newindex(self::Contract) = maximum(flatten((self.indices..., self.target))) + 1
Base.size(self::Contract) = self.shape

arguments(self::Contract) = Tuple(self.terms)
@destructure prealloc(self::Contract{T}, []) where T = [:((Array{$T})(undef, $(size(self)...)))]

# Here, we essentially roll our own einsum macro. Options are:
# - TensorOperations.jl: optimized for large arrays, which is typically not the case for us.
#   Generated code always allocates.
# - Einsum.jl: Generated code doesn't allocate, but not possible to avoid bounds checking,
#   which we can statically guarantee is not needed.
# May be changed in the future.
@destructure function codegen(self::Contract, terms, [target])
    counters = Dict{Int,Symbol}(i => gensym(string(i)) for i in flatten((self.indices..., self.target)))
    lengths = Dict{Int,Int}(i => size(term, j)
                            for (term, inds) in zip(self.terms, self.indices)
                            for (j, i) in enumerate(inds))
    caxes = Set{Int}(k for k in keys(lengths) if k ∉ self.target)
    taxes = Set{Int}(k for k in keys(lengths) if k in self.target)

    c = i -> (counters[j] for j in i)

    # Code for incrementing an element in the target array
    getcodes = (:($t[$(c(i)...)]) for (t, i) in zip(terms, self.indices))
    code = :(@inbounds $target[$(c(self.target)...)] += *($(getcodes...)))

    # Code for looping over all contracted axes
    for i in caxes
        code = :(for $(counters[i]) in 1:$(lengths[i]); $code end)
    end

    # Initialize accumulator to zero
    code = quote
        $target[$(c(self.target)...)] = $(zero(eltype(self)))
        $code
    end

    # Code for looping over all non-contracted axes
    for i in taxes
        code = :(for $(counters[i]) in 1:$(lengths[i]); $code end)
    end

    :($code; $target)
end

# DropDims

@autohasheq struct DropDims{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable
    dims :: Tuple{Vararg{Int}}

    function DropDims(source::ArrayEvaluable, dims::Tuple{Vararg{Int}})
        @assert all(size(source, d) == 1 for d in dims)
        new{eltype(source), ndims(source)-length(dims)}(source, dims)
    end
end

Base.size(self::DropDims) = Tuple(s for (i, s) in enumerate(size(self.source)) if i ∉ self.dims)

isnormal(::DropDims) = false
arguments(self::DropDims) = (self.source,)
@destructure prealloc(::DropDims, []) = []
@destructure codegen(self::DropDims, [source], []) = :(dropdims($source, dims=$(self.dims)))

# GetIndex
# TODO: More general indexing expressions?

@autohasheq struct GetIndex{T,N} <: ArrayEvaluable{T,N}
    value :: ArrayEvaluable
    indices :: Indices
    map :: IndexMap

    function GetIndex(value::ArrayEvaluable{T,N}, map::IndexMap) where {T,N}
        dimcheck(map, 0, 1) || error("Multidimensional indices not yet supported")
        indices = Index[get(map, i, :) for i in 1:ndims(value)]
        newdims = resultndims(map, ndims(value))
        new{T,newdims}(value, indices, map)
    end
end

Base.size(self::GetIndex) = resultsize(self.map, size(self.value))

isnormal(::GetIndex) = false
arguments(self::GetIndex) = (self.value, values(self.map)...)
@destructure prealloc(::GetIndex, []) = []

@destructure function codegen(self::GetIndex, [value, varindices...], [])
    :(view($value, $(codegen(self.map, varindices, ndims(self.value))...)))
end

# Inflate

@autohasheq struct Inflate{T,N} <: ArrayEvaluable{T,N}
    data :: ArrayEvaluable{T}
    shape :: Shape
    map :: IndexMap

    function Inflate(data::ArrayEvaluable, shape::Shape, map::IndexMap)
        ndims(data) == length(shape) || error("Dimension mismatch")
        resultsize(map, shape) == size(data) || error("Size mismatch")
        dimcheck(map, 1, 1) || error("Only linear indices supported in inflation")
        new{eltype(data), length(shape)}(data, shape, map)
    end
end

Inflate(data::ArrayEvaluable, shape::Shape, indices::Pair...) = Inflate(data, shape, IndexMap(indices...))

Base.size(self::Inflate) = self.shape
separate(self::Inflate) = [
    (Tupl((axis in keys(self.map) ? getindex(self.map[axis], ind) : ind)  for (axis, ind) in enumerate(inds)), data)
    for (inds, data) in separate(self.data)
]

arguments(self::Inflate) = (self.data, values(self.map)...)
@destructure prealloc(self::Inflate{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]

@destructure function codegen(self::Inflate{T}, [data, indices...], [target]) where T
    weaved = codegen(self.map, indices, ndims(self))
    quote
        $target[:] .= $(zero(T))
        $target[$(weaved...)] = $data
        $target
    end
end

# InsertAxis

@autohasheq struct InsertAxis{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable{T}
    axes :: Vector{Int}

    function InsertAxis(source::ArrayEvaluable, axes::Vector{Int})
        isempty(axes) && return source
        all(1 <= i <= ndims(source) + 1 for i in axes) || error("Axes out of range")
        new{eltype(source), ndims(source)+length(axes)}(Normalize(source), sort(axes))
    end
end

function Base.size(self::InsertAxis)
    shape = collect(Int, size(self.source))
    insertmany!(shape, self.axes, 1)
    Tuple(shape)
end

isnormal(::InsertAxis) = false
arguments(self::InsertAxis) = (self.source,)
@destructure prealloc(self::InsertAxis, [dep]) = [:(reshape($dep, $(size(self)...)))]
@destructure codegen(self::InsertAxis, [], [pre]) = pre

# Inv

@autohasheq struct Inv{T,N} <: ArrayEvaluable{T,2}
    source :: ArrayEvaluable{T,2}

    function Inv(source::ArrayEvaluable{T,2}) where T
        size(source, 1) == size(source, 2) || error("Expected square matrix")
        new{T, size(source,1)}(source)
    end
end

Base.size(self::Inv) = size(self.source)

arguments(self::Inv) = (self.source,)
@destructure prealloc(self::Inv{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]

# Would be useful with an inv! function here
@destructure codegen(::Inv, [source], [target]) = :($target .= inv($source); $target)
@destructure codegen(::Inv{T,1}, [source], [target]) where T = :($target[1,1] = one($T) / $source[1,1]; $target)
@destructure codegen(::Inv{T,2}, [source], [target]) where T = quote
    $target[1,1] = $source[2,2]
    $target[2,2] = $source[1,1]
    $target[1,2] = -$source[1,2]
    $target[2,1] = -$source[2,1]
    $target ./= ($source[1,1] * $source[2,2] - $source[1,2] * $source[2,1])
    $target
end
@destructure codegen(::Inv{T,3}, [source], [target]) where T = quote
    $target[1,1] = $source[2,2] * $source[3,3] - $source[2,3] * $source[3,2]
    $target[2,1] = $source[2,3] * $source[3,1] - $source[2,1] * $source[3,3]
    $target[3,1] = $source[2,1] * $source[3,2] - $source[2,2] * $source[3,1]
    $target[1,2] = $source[1,3] * $source[3,2] - $source[1,2] * $source[3,3]
    $target[2,2] = $source[1,1] * $source[3,3] - $source[1,3] * $source[3,1]
    $target[3,2] = $source[1,2] * $source[3,1] - $source[1,1] * $source[3,2]
    $target[1,3] = $source[1,2] * $source[2,3] - $source[1,3] * $source[2,2]
    $target[2,3] = $source[1,3] * $source[2,1] - $source[1,1] * $source[2,3]
    $target[3,3] = $source[1,1] * $source[2,2] - $source[1,2] * $source[2,1]
    $target ./= $source[1,1] * $target[1,1] + $source[1,2] * $target[2,1] + $source[1,3] * $target[3,1]
    $target
end

# Monomials

@autohasheq struct Monomials{T,N} <: ArrayEvaluable{T,N}
    points :: ArrayEvaluable{T}
    degree :: Int
    padding :: Int

    Monomials(points::ArrayEvaluable{T,M}, degree::Int, padding::Int) where {T,M} = new{T,M+1}(points, degree, padding)
    Monomials(points::ArrayEvaluable, degree::Int) = Monomials(points, degree, 0)
end

Base.size(self::Monomials) = (self.degree + self.padding + 1, size(self.points)...)

arguments(self::Monomials) = (self.points,)
@destructure prealloc(self::Monomials{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]

@destructure function codegen(self::Monomials{T}, [points], [target]) where T
    j = gensym("j")

    loopcode = Any[:(@inbounds $target[$(1+self.padding),$j] = 1.0)]
    for i in self.padding .+ (1:self.degree)
        push!(loopcode, :(@inbounds $target[$(i+1),$j] = $target[$i,$j] * $points[$j]))
    end

    quote
        $target .= zero($T)
        for $j = CartesianIndices($(size(self.points)))
            $(loopcode...)
        end
        $target
    end
end

# Neg

@autohasheq struct Neg{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable{T,N}
end

Base.size(self::Neg) = size(self.source)

arguments(self::Neg) = (self.source,)
@destructure prealloc(::Neg, []) = []
@destructure codegen(::Neg, [source], []) = :(-$source)

# Normalize

@autohasheq struct Normalize{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable{T,N}

    function Normalize(source::ArrayEvaluable)
        isnormal(source) && return source
        new{eltype(source), ndims(source)}(source)
    end
end

Base.size(self::Normalize) = size(self.source)

arguments(self::Normalize) = (self.source,)
@destructure prealloc(self::Normalize{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]
@destructure codegen(::Normalize, [source], [target]) = :($target .= $source; $target)

# Product

@autohasheq struct Product{T,N} <: ArrayEvaluable{T,N}
    terms :: Tuple{Vararg{ArrayEvaluable}}

    function Product(terms::Tuple{Vararg{ArrayEvaluable}})
        newshape = broadcast_shape((size(term) for term in terms)...)
        newtype = reduce(promote_type, (eltype(term) for term in terms))
        new{newtype, length(newshape)}(terms)
    end
end

Base.size(self::Product) = broadcast_shape((size(term) for term in self.terms)...)

arguments(self::Product) = self.terms
@destructure prealloc(self::Product{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]

# Note: element-wise product!
@destructure function codegen(self::Product, args, [target])
    code = args[1]
    for arg in args[2:end]
        code = :($code .* $arg)
    end
    :($target .= $code; $target)
end

# Reshape

@autohasheq struct Reshape{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable{T}
    newshape :: Shape

    function Reshape(source::ArrayEvaluable, newshape::Shape)
        @assert prod(newshape) == length(source)
        new{eltype(source), length(newshape)}(Normalize(source), newshape)
    end
end

Base.size(self::Reshape) = self.newshape

isnormal(::Reshape) = false
arguments(self::Reshape) = (self.source,)
@destructure prealloc(self::Reshape, [dep]) = [:(reshape($dep, $(size(self)...)))]
@destructure codegen(self::Reshape, [], [pre]) = pre

# Sum

@autohasheq struct Sum{T,N} <: ArrayEvaluable{T,N}
    source :: ArrayEvaluable{T}
    dims :: Tuple{Vararg{Int}}

    function Sum(source::ArrayEvaluable, dims::Tuple{Vararg{Int}})
        @assert all(1 <= d <= ndims(source) for d in dims)
        new{eltype(source), ndims(source)}(source, dims)
    end
end

Base.size(self::Sum) = Tuple(i in self.dims ? 1 : k for (i, k) in enumerate(size(self.source)))

arguments(self::Sum) = (self.source,)
@destructure prealloc(self::Sum{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]
@destructure codegen(self::Sum, [source], [target]) = :($target .= sum($source, dims=$(self.dims)))

# Add

@autohasheq struct Add{T,N} <: ArrayEvaluable{T,N}
    terms :: Tuple{Vararg{ArrayEvaluable}}

    function Add(terms::Tuple{Vararg{ArrayEvaluable}})
        newshape = broadcast_shape((size(term) for term in terms)...)
        newtype = reduce(promote_type, (eltype(term) for term in terms))
        new{newtype, length(newshape)}(terms)
    end
end

Add(terms...) = Add(terms)

Base.size(self::Add) = broadcast_shape((size(term) for term in self.terms)...)

arguments(self::Add) = self.terms
@destructure prealloc(self::Add{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]

@destructure function codegen(self::Add, args, [target])
    code = args[1]
    for arg in args[2:end]
        code = :($code .+ $arg)
    end
    :($target .= $code; $target)
end

# Tupl(e)

@autohasheq struct Tupl{T} <: Evaluable{T}
    terms :: Tuple{Vararg{Evaluable}}

    function Tupl(terms)
        terms = Tuple(terms)
        # TODO: Is there a better way to do this?
        tupletype = eval(:(Tuple{$((restype(term) for term in terms)...)}))
        new{tupletype}(terms)
    end
end

Base.iterate(self::Tupl) = iterate(self.terms)
Base.iterate(self::Tupl, i::Int) = iterate(self.terms, i)
Base.length(self::Tupl) = length(self.terms)

arguments(self::Tupl) = self.terms
@destructure prealloc(self::Tupl, []) = []
@destructure codegen(self::Tupl, terms, []) = :($(terms...),)

# Zeros

@autohasheq struct Zeros{T,N} <: ArrayEvaluable{T,N}
    shape :: Shape
    Zeros(T, shape::Int...) = new{T,length(shape)}(shape)
end

Base.size(self::Zeros) = self.shape
separate(::Zeros) = []

arguments(::Zeros) = ()
@destructure prealloc(self::Zeros{T}, []) where T = [:(Array{$T}(undef, $(size(self)...)))]
@destructure codegen(self::Zeros{T}, [], [target]) where T = :($target[:] .= zero($T); $target)

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
    using DifferentialEquations
    using OrdinaryDiffEq
    using DiffEqBase
    using Sundials
    using ODEInterfaceDiffEq

        # du -> Derivative
        # u -> State at time t
        # p-> paramter vector
        # t-> time as tuple (init, end)

        function Hes1ExampleODE!(du,u,p,t)
            m, p1, p2 = u;

            k1, P0, v, h, inps = p;

        du[1] = dm = -0.03*m + (1/(1+(p2/P0)^h))*inps;
        du[2] = dp1 = -0.03*p1 + v*m - k1*p1;
        du[3] = dp2 = -0.03*p2 + k1*p1;

        end

        # p -> parameter vector (plus inducer vector at the end)
        # I -> initial Y0 vector

        function Hes1ExampleSteadyState(p,I)

            m, p1, p2 = I;
               return(I)

        end

        # ts -> time vector going from t=0 to t=end every 1
        # p -> parameter vector
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # pre -> Vector with the inducers in the ON. It can be empty

        function Hes1Example_solvecoupledODE(ts, p, sp, inputs, ivss, pre=[])

            maxtime = length(ts);
            Nsp = length(sp);
            Nevents = length(sp)-1;
            Neq = 3;
            if pre != []
                pSS = vcat(p, pre);
            else
                pSS = p;
            end

            final = zeros(maxtime, Neq);

            y_al = Hes1ExampleSteadyState(pSS,ivss) # Calculation of initial guesses for steady state

            y0 = y_al;
            initialV = y0;
            i = 1;

            for q in collect(1:Nevents)
                lts = length(ts[(sp[q]+1):sp[q+1]+1]);  # General way to define the number of elements in each event series
                Tevent = ts[(sp[q]+1):sp[q+1]+1];  # General way to extract the times of each event
                I = inputs[i:(i+(0))];
                pSte = vcat(p, I);

                if q == 1
                    prob = ODEProblem(Hes1ExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                    part1 = DifferentialEquations.solve(prob, Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                else
                    prob = ODEProblem(Hes1ExampleODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                    part1 = DifferentialEquations.solve(prob, Tsit5(),reltol=1.0e-9,abstol=1.0e-9,saveat=1);
                end

                initialV = part1[lts];
                i+=1;

                for d in collect((sp[q]+1):(sp[q]+lts))
                    final[d,:] = part1[d-sp[q]];
                end
            end
            return(final)

        end

        # ts -> time vector going from t=0 to t=end every 1
        # pD -> Matrix containing all the draws for the parameters (can be a single vector if nly 1 draw of parameters is considered)
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # samps -> Vector of Time points that want to be extracted from the simulations (it considers the initial time point as 0)
        # pre -> Vector with the inducers in the ON. It can be empty

        function Hes1Example_SolveAll(ts, pD, sp, inputs, ivss, samps, pre=[])

            if length(size(pD)) == 1
                pD = reshape(pD,size(pD)[1],1);
            end

            if size(pD)[2] != 4
                pD = pD';
            end

            if length(ivss)/3 > 1
                if size(ivss)[2] != 3
                    ivss = ivss';
                end
            end

            AllSolTest = zeros(length(samps), 3, length(pD[:,1]));

            for drawInd in collect(1:length(pD[:,1]))
                p = pD[drawInd,:];
                if length(ivss)/3 > 1
                    ivss2 = ivss[drawInd,:];
                else
                    ivss2 = ivss;
                end
                temp = Hes1Example_solvecoupledODE(ts, p, sp, inputs, ivss2, pre);
                AllSolTest[:,:,drawInd] = temp[convert.(Int,samps.+1),:];
            end

            return(AllSolTest);
        end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct MarconatoMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Xd_pp::Float64
        Xq_pp::Float64
        Td0_p::Float64
        Tq0_p::Float64
        Td0_pp::Float64
        Tq0_pp::Float64
        T_AA::Float64
        ext::Dict{String, Any}
        γd::Float64
        γq::Float64
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of 6-states synchronous machine: Marconato model

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: `(0, nothing)`
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit, validation range: `(0, nothing)`
- `Xq_pp::Float64`: Sub-Transient reactance after EMF in q-axis per unit, validation range: `(0, nothing)`
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: `(0, nothing)`
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: `(0, nothing)`
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: `(0, nothing)`
- `T_AA::Float64`: Time constant of d-axis additional leakage, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `γd::Float64`
- `γq::Float64`
- `states::Vector{Symbol}`: The states are:
	ψq: q-axis stator flux,
	ψd: d-axis stator flux,
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage,
	eq_pp: q-axis subtransient voltage,
	ed_pp: d-axis subtransient voltage
- `n_states::Int`: MarconatoMachine has 6 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct MarconatoMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit"
    Xd_pp::Float64
    "Sub-Transient reactance after EMF in q-axis per unit"
    Xq_pp::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Time constant of d-axis additional leakage"
    T_AA::Float64
    ext::Dict{String, Any}
    γd::Float64
    γq::Float64
    "The states are:
	ψq: q-axis stator flux,
	ψd: d-axis stator flux,
	eq_p: q-axis transient voltage,
	ed_p: d-axis transient voltage,
	eq_pp: q-axis subtransient voltage,
	ed_pp: d-axis subtransient voltage"
    states::Vector{Symbol}
    "MarconatoMachine has 6 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function MarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, ext=Dict{String, Any}(), )
    MarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, ext, ((Td0_pp*Xd_pp)/(Td0_p*Xd_p) )*(Xd-Xd_p), ((Tq0_pp*Xq_pp)/(Tq0_p*Xq_p) )*(Xq-Xq_p), [:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp], 6, InfrastructureSystemsInternal(), )
end

function MarconatoMachine(; R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, ext=Dict{String, Any}(), γd=((Td0_pp*Xd_pp)/(Td0_p*Xd_p) )*(Xd-Xd_p), γq=((Tq0_pp*Xq_pp)/(Tq0_p*Xq_p) )*(Xq-Xq_p), states=[:ψq, :ψd, :eq_p, :ed_p, :eq_pp, :ed_pp], n_states=6, internal=InfrastructureSystemsInternal(), )
    MarconatoMachine(R, Xd, Xq, Xd_p, Xq_p, Xd_pp, Xq_pp, Td0_p, Tq0_p, Td0_pp, Tq0_pp, T_AA, ext, γd, γq, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function MarconatoMachine(::Nothing)
    MarconatoMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Xd_pp=0,
        Xq_pp=0,
        Td0_p=0,
        Tq0_p=0,
        Td0_pp=0,
        Tq0_pp=0,
        T_AA=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`MarconatoMachine`](@ref) `R`."""
get_R(value::MarconatoMachine) = value.R
"""Get [`MarconatoMachine`](@ref) `Xd`."""
get_Xd(value::MarconatoMachine) = value.Xd
"""Get [`MarconatoMachine`](@ref) `Xq`."""
get_Xq(value::MarconatoMachine) = value.Xq
"""Get [`MarconatoMachine`](@ref) `Xd_p`."""
get_Xd_p(value::MarconatoMachine) = value.Xd_p
"""Get [`MarconatoMachine`](@ref) `Xq_p`."""
get_Xq_p(value::MarconatoMachine) = value.Xq_p
"""Get [`MarconatoMachine`](@ref) `Xd_pp`."""
get_Xd_pp(value::MarconatoMachine) = value.Xd_pp
"""Get [`MarconatoMachine`](@ref) `Xq_pp`."""
get_Xq_pp(value::MarconatoMachine) = value.Xq_pp
"""Get [`MarconatoMachine`](@ref) `Td0_p`."""
get_Td0_p(value::MarconatoMachine) = value.Td0_p
"""Get [`MarconatoMachine`](@ref) `Tq0_p`."""
get_Tq0_p(value::MarconatoMachine) = value.Tq0_p
"""Get [`MarconatoMachine`](@ref) `Td0_pp`."""
get_Td0_pp(value::MarconatoMachine) = value.Td0_pp
"""Get [`MarconatoMachine`](@ref) `Tq0_pp`."""
get_Tq0_pp(value::MarconatoMachine) = value.Tq0_pp
"""Get [`MarconatoMachine`](@ref) `T_AA`."""
get_T_AA(value::MarconatoMachine) = value.T_AA
"""Get [`MarconatoMachine`](@ref) `ext`."""
get_ext(value::MarconatoMachine) = value.ext
"""Get [`MarconatoMachine`](@ref) `γd`."""
get_γd(value::MarconatoMachine) = value.γd
"""Get [`MarconatoMachine`](@ref) `γq`."""
get_γq(value::MarconatoMachine) = value.γq
"""Get [`MarconatoMachine`](@ref) `states`."""
get_states(value::MarconatoMachine) = value.states
"""Get [`MarconatoMachine`](@ref) `n_states`."""
get_n_states(value::MarconatoMachine) = value.n_states
"""Get [`MarconatoMachine`](@ref) `internal`."""
get_internal(value::MarconatoMachine) = value.internal

"""Set [`MarconatoMachine`](@ref) `R`."""
set_R!(value::MarconatoMachine, val) = value.R = val
"""Set [`MarconatoMachine`](@ref) `Xd`."""
set_Xd!(value::MarconatoMachine, val) = value.Xd = val
"""Set [`MarconatoMachine`](@ref) `Xq`."""
set_Xq!(value::MarconatoMachine, val) = value.Xq = val
"""Set [`MarconatoMachine`](@ref) `Xd_p`."""
set_Xd_p!(value::MarconatoMachine, val) = value.Xd_p = val
"""Set [`MarconatoMachine`](@ref) `Xq_p`."""
set_Xq_p!(value::MarconatoMachine, val) = value.Xq_p = val
"""Set [`MarconatoMachine`](@ref) `Xd_pp`."""
set_Xd_pp!(value::MarconatoMachine, val) = value.Xd_pp = val
"""Set [`MarconatoMachine`](@ref) `Xq_pp`."""
set_Xq_pp!(value::MarconatoMachine, val) = value.Xq_pp = val
"""Set [`MarconatoMachine`](@ref) `Td0_p`."""
set_Td0_p!(value::MarconatoMachine, val) = value.Td0_p = val
"""Set [`MarconatoMachine`](@ref) `Tq0_p`."""
set_Tq0_p!(value::MarconatoMachine, val) = value.Tq0_p = val
"""Set [`MarconatoMachine`](@ref) `Td0_pp`."""
set_Td0_pp!(value::MarconatoMachine, val) = value.Td0_pp = val
"""Set [`MarconatoMachine`](@ref) `Tq0_pp`."""
set_Tq0_pp!(value::MarconatoMachine, val) = value.Tq0_pp = val
"""Set [`MarconatoMachine`](@ref) `T_AA`."""
set_T_AA!(value::MarconatoMachine, val) = value.T_AA = val
"""Set [`MarconatoMachine`](@ref) `ext`."""
set_ext!(value::MarconatoMachine, val) = value.ext = val
"""Set [`MarconatoMachine`](@ref) `γd`."""
set_γd!(value::MarconatoMachine, val) = value.γd = val
"""Set [`MarconatoMachine`](@ref) `γq`."""
set_γq!(value::MarconatoMachine, val) = value.γq = val


text:
export segy_read

"""
block = segy_read(aws::AWSCore.AWSConfig, bucket::String, path::String)

# Optional Keyword Arguments

- `multipart=false`: read entire object in single REST request (beware of memory overhead for large objects)
- `buffer_size=1024`: buffer size [MB] if multipart=true
- `warn_user=true`: explict warnings

"""
function segy_read(aws::AWSCore.AWSConfig, bucket::String, path::String; multipart::Bool=false, buffer_size::Int=1024, warn_user::Bool=true)

    s3_exists(aws, bucket, path) || @error "AWSS3SegyIO/segy_read: file $path does not exist in $bucket."
    if multipart
        return read_file(aws, bucket, path, warn_user)
    else
        s = IOBuffer(s3_get(aws, bucket, path);)
        return read_file(s, warn_user)
    end

end

"""
block = segy_read(aws::AWSCore.AWSConfig, bucket::String, path::String, keys::Array{String,1})

# Optional Keyword Arguments

- `multipart=false`: read entire object in single REST request (beware of memory overhead for large objects)
- `buffer_size=1024`: buffer size [MB] if multipart=true
- `warn_user=true`: explict warnings

"""
function segy_read(aws::AWSCore.AWSConfig, bucket::String, path::String, keys::Array{String,1}; multipart::Bool=false, buffer_size::Int=1024, warn_user::Bool=true)

    s3_exists(aws, bucket, path) || @error "AWSS3SegyIO/segy_read: file $path does not exist in $bucket."
    if multipart
        return read_file(aws, bucket, path, keys, warn_user)
    else
        s = IOBuffer(s3_get(aws, bucket, path);)
        return read_file(s, keys, warn_user)
    end

end