text:
# Create a parabolic velocity inlet
abstract VelocityIOlet <: IOLetStreaming
type ParabolicVelocityIOlet{T} <: VelocityIOlet
    # Direction of the inlet
    normal :: Vector{T}
    # Center of the inlet
    center :: Vector{T}
    # Radius of the inlet
    radius :: T
    # Maximum speed
    maxspeed :: T

    function ParabolicVelocityIOlet(n₀, Γ, r₀, ν_max)
        n = T[n₀...]
        @assert ndims(n₀) == 1 && ndims(Γ) == 1 && length(Γ) == length(n₀)
        @assert dot(n, n) > 1e-12 && r₀ > 1e-12 && ν_max > 1e-12
        new(n / norm(n), T[Γ...], r₀, ν_max)
    end
end

streaming(
    streamer::VelocityIOlet, quantities::LocalQuantities, sim::Simulation,
    from::Integer, to::Integer, direction::Integer
) = streaming(streamer, sim, from, direction)
function streaming(streamer::VelocityIOlet, sim::Simulation, from::Integer, direction::Integer)
    const bb_direction = sim.lattice.inversion[direction]
    const cᵢ = sim.lattice.celerities[:, direction]
    const wᵢ = sim.lattice.weights[direction]
    const halfway = gridcoords(sim.indexing, from) + 0.5cᵢ
    const μ = velocity(streamer, halfway, sim.time)
    const correction = 2wᵢ * dot(cᵢ, μ) / speed_of_sound_squared
    sim.next_populations[bb_direction, from] = sim.populations[direction, from] - correction
end

# Time-independent parabolic inlet
velocity(iolet::ParabolicVelocityIOlet, position, time) = velocity(iolet, position)
function velocity(iolet::ParabolicVelocityIOlet, position)
    const centered = position - iolet.center
    const z = dot(centered, iolet.normal)
    const radial = (dot(centered, centered) - z * z) / iolet.radius / iolet.radius
    ((1 - radial) * iolet.maxspeed) * iolet.normal
end

type ConstantVelocityIOlet{T <: AbstractFloat} <: VelocityIOlet
    velocity :: Vector{T}
end
streaming(
    streamer::ConstantVelocityIOlet, quantities::LocalQuantities, sim::Simulation,
    from::Integer, to::Integer, direction::Integer
) = streaming(streamer, sim, from, direction)
function streaming(iolet::ConstantVelocityIOlet, sim::Simulation, from::Integer, direction::Integer)
    const bb_direction = sim.lattice.inversion[direction]
    const cᵢ = sim.lattice.celerities[:, direction]
    const wᵢ = sim.lattice.weights[direction]
    const correction = 2wᵢ * dot(cᵢ, iolet.velocity) / speed_of_sound_squared
    sim.next_populations[bb_direction, from] = sim.populations[direction, from] - correction
end

type NashZeroOrderPressure{T} <: IOLetStreaming
    # Direction of the inlet
    normal :: Vector{T}
    # Density at the inlet
    density :: T
    function NashZeroOrderPressure(n₀::Vector, ρ::Number)
        n = T[n₀...]
        @assert dot(n, n) > 1e-12 && ρ > 1e-12
        new(n / norm(n), ρ)
    end
end

streaming(
    iolet::NashZeroOrderPressure, quantities::LocalQuantities, sim::Simulation,
    from::Integer, to::Integer, direction::Integer
) = streaming(iolet, quantities.velocity, sim, from, direction)
function streaming(
        iolet::NashZeroOrderPressure, velocity::Vector, sim::Simulation,
        from::Integer, direction::Integer)
    const invdir = sim.lattice.inversion[direction]
    const iolet_μ = iolet.density * dot(iolet.normal, velocity) * iolet.normal
    const cᵢ = sim.lattice.celerities[:, invdir:invdir]
    const wᵢ = sim.lattice.weights[invdir:invdir]
    const feq = equilibrium(iolet.density, iolet_μ, cᵢ, wᵢ)[1]
    sim.next_populations[invdir, from] = feq
end

type ConstantPopulationIOlet{I <: Integer} <: IOLetStreaming
    # Normal to iolet
    # Should be a lattice vector
    # Should points from fluid site to iolet site
    normal :: Vector{I}
end

function streaming(
    iolet::ConstantPopulationIOlet, quantities::LocalQuantities, sim::Simulation,
    from::Integer, to::Integer, direction::Integer)
    # This object can stream to a different object than expected
    const invdir = sim.lattice.inversion[direction]
    const mirror = gridcoords(sim.indexing, from) + iolet.normal + sim.lattice.celerities[:, invdir]
    const mirror_site = index(sim.indexing, mirror)
    const link = to == 0 ? NOTHING: sim.playground[mirror_site]
    const this_kernel = sim.kernels[sim.playground[from]]
    const streamer = get(this_kernel.streamers, link, NULLSTREAMER)
    @assert !(typeof(streamer) <: IOLetStreaming)
    streaming(streamer, quantities, sim, from, mirror_site, invdir)
end


text:
using Documenter, GoogleAnalytics

makedocs(;
    modules=[GoogleAnalytics],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/vly/GoogleAnalytics.jl/blob/{commit}{path}#L{line}",
    sitename="GoogleAnalytics.jl",
    authors="Val Lyashov, Envato",
    assets=String[],
)

deploydocs(;
    repo="github.com/vly/GoogleAnalytics.jl",
)

ext: jl


text:
const DECLARATIVE_PACKAGES_DIR = Pkg.dir("DeclarativePackages")
const NIGHTLY = v"1.3.0-"
const JULIA_BIN_MODE = Mode(user=(READ+WRITE+EXEC), group=(READ+EXEC), other=(READ+EXEC))
const DEFAULT_CONFIG = """
---
# This is just default location to store a new playground.
# This is used by create and activate if no --name or --path.
default_playground_path: .playground

# Uncomment below to change the default shell. Otherwise the SHELL
# Environment variable will be used.
# default_shell: /usr/local/bin/fish

# Default julia registry settings.
default_registry: \"https://github.com/JuliaLang/METADATA.jl.git\"
default_branch: \"metadata-v2\"

# Default git settings when using install build
default_git_address: \"https://github.com/JuliaLang/julia.git\"
default_git_revision: master

# Allows you to isolate shell and julia history to each playground.
isolated_shell_history: true
isolated_julia_history: true
"""


text:
"""
    AST.record(::ASTMDP, sym::Symbol, val)

Recard an ASTMetric specified by `sym`.
"""
function record(mdp::ASTMDP, sym::Symbol, val)
    if mdp.params.debug
        push!(getproperty(mdp.metrics, sym), val)
    end
end

function record(mdp::ASTMDP; prob=1, logprob=exp(prob), miss_distance=Inf, reward=-Inf, event=false, terminal=false, rate=-Inf)
    AST.record(mdp, :prob, prob)
    AST.record(mdp, :logprob, logprob)
    AST.record(mdp, :miss_distance, miss_distance)
    AST.record(mdp, :reward, reward)
    AST.record(mdp, :intermediate_reward, reward)
    AST.record(mdp, :rate, rate)
    AST.record(mdp, :event, event)
    AST.record(mdp, :terminal, terminal)
end

function record_returns(mdp::ASTMDP)
    # compute returns up to now.
    rewards = mdp.metrics.intermediate_reward
    G = returns(rewards, γ=discount(mdp))
    AST.record(mdp, :returns, G)
    mdp.metrics.intermediate_reward = [] # reset
end

function returns(R; γ=1)
    T = length(R)
    G = zeros(T)
    for t in reverse(1:T)
        G[t] = t==T ? R[t] : G[t] = R[t] + γ*G[t+1]
    end
    return G
end

function combine_ast_metrics(plannervec::Vector)
    return ASTMetrics(
        miss_distance=vcat(map(planner->planner.mdp.metrics.miss_distance, plannervec)...),
        rate=vcat(map(planner->planner.mdp.metrics.rate, plannervec)...),
        logprob=vcat(map(planner->planner.mdp.metrics.logprob, plannervec)...),
        prob=vcat(map(planner->planner.mdp.metrics.prob, plannervec)...),
        reward=vcat(map(planner->planner.mdp.metrics.reward, plannervec)...),
        intermediate_reward=vcat(map(planner->planner.mdp.metrics.intermediate_reward, plannervec)...),
        returns=vcat(map(planner->planner.mdp.metrics.returns, plannervec)...),
        event=vcat(map(planner->planner.mdp.metrics.event, plannervec)...),
        terminal=vcat(map(planner->planner.mdp.metrics.terminal, plannervec)...))
end


text:
@testset "Elliptic curves" begin
  include("EllCrv/EllCrv.jl")
  include("EllCrv/QQ.jl")
#  include("EllCrv/Finite.jl")
end


text:
function action_info(p::POMCPPlanner, b; tree_in_info=false)
    local a::actiontype(p.problem)
    info = Dict{Symbol, Any}()
    try
        tree = POMCPTree(p.problem, b, p.solver.tree_queries)
        a = search(p, b, tree, info)
        p._tree = tree
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end
    catch ex
        # Note: this might not be type stable, but it shouldn't matter too much here
        # a = convert(actiontype(p.problem), default_action(p.solver.default_action, p.problem, b, ex))
        a = convert(actiontype(p.problem), default_action(RandomPolicy(p.problem, rng=MersenneTwister(1)), p.problem, b, ex))
        info[:exception] = ex
    end
    return a, info
end

action(p::POMCPPlanner, b) = first(action_info(p, b))

function search(p::POMCPPlanner, b, t::POMCPTree, info::Dict)
    all_terminal = true
    nquery = 0
    start_us = CPUtime_us()
    info[:tree_depth] = 0
    dep_count=0
    for i in 1:p.solver.tree_queries
        nquery += 1
        if CPUtime_us() - start_us >= 1e6*p.solver.max_time
            break
        end
        s = rand(p.rng, b)
        if !POMDPs.isterminal(p.problem, s)
            reward, depth = simulate(p, s, POMCPObsNode(t, 1), p.solver.max_depth)
            info[:tree_depth] = (info[:tree_depth] * dep_count + depth) / (dep_count + 1)
            dep_count+=1
            # info[:tree_depth] = max(info[:tree_depth], depth)
            all_terminal = false
        end
    end
    info[:search_time_us] = CPUtime_us() - start_us
    info[:tree_queries] = nquery

    if all_terminal
        throw(AllSamplesTerminal(b))
    end

    h = 1
    best_node = first(t.children[h])
    best_v = t.v[best_node]
    @assert !isnan(best_v)
    for node in t.children[h][2:end]
        if t.v[node] >= best_v
            best_v = t.v[node]
            best_node = node
        end
    end

    return t.a_labels[best_node]
end

solve(solver::POMCPSolver, pomdp::POMDP) = POMCPPlanner(solver, pomdp)

function simulate(p::POMCPPlanner, s, hnode::POMCPObsNode, steps::Int)
    dep = 1.0
    if steps == 0 || isterminal(p.problem, s)
        return 0.0, dep
    end
    t = hnode.tree
    h = hnode.node

    ltn = log(t.total_n[h])
    best_nodes = empty!(p._best_node_mem)
    best_criterion_val = -Inf
    for node in t.children[h]
        n = t.n[node]
        if n == 0 && ltn <= 0.0
            criterion_value = t.v[node]
        elseif n == 0 && t.v[node] == -Inf
            criterion_value = Inf
        else
            criterion_value = t.v[node] + p.solver.c*sqrt(ltn/n)
        end
        if criterion_value > best_criterion_val
            best_criterion_val = criterion_value
            empty!(best_nodes)
            push!(best_nodes, node)
        elseif criterion_value == best_criterion_val
            push!(best_nodes, node)
        end
    end
    ha = rand(p.rng, best_nodes)
    a = t.a_labels[ha]

    sp, o, r = @gen(:sp, :o, :r)(p.problem, s, a, p.rng)

    hao = get(t.o_lookup, (ha, o), 0)
    if hao == 0
        dep = dep + 1.0
        hao = insert_obs_node!(t, p.problem, ha, sp, o)
        v = estimate_value(p.solved_estimator,
                           p.problem,
                           sp,
                           POMCPObsNode(t, hao),
                           steps-1)
        R = r + discount(p.problem)*v
    else
        reward, depth = simulate(p, sp, POMCPObsNode(t, hao), steps-1)
        dep = depth + 1
        R = r + discount(p.problem) * reward
    end

    t.total_n[h] += 1
    t.n[ha] += 1
    t.v[ha] += (R-t.v[ha])/t.n[ha]

    return R, dep
end


text:
module MendelIterHardThreshold

using MendelBase
using SnpArrays
using IHT
using DataFrames
using Distances
using StatsBase

export IterHardThreshold

const Float = Union{Float64,Float32}

include("IHT_utilities.jl")

# """
# This is the wrapper function for the Iterative Hard Thresholding analysis option in Open Mendel.
# """
# function MendelIHT(file_name::String, k::Int64)
#     const MENDEL_IHT_VERSION :: VersionNumber = v"0.2.0"
#     #
#     # Print the logo. Store the initial directory.
#     #
#     print(" \n \n")
#     println("     Welcome to OpenMendel's")
#     println("      IHT analysis option")
#     println("        version ", MENDEL_IHT_VERSION)
#     print(" \n \n")

#     snpmatrix = SnpArray(file_name)
#     phenotype = readdlm(file_name * ".fam", header = false)[:, 6]
#     # phenotype = randn(959) #testing GAW data since it has no phenotype

#     return L0_reg(snpmatrix, phenotype, k)
# end #function MendelIHT

"""
This is the wrapper function for the Iterative Hard Thresholding analysis option in Open Mendel.
"""
function IterHardThreshold(control_file = ""; args...)
    const MENDEL_IHT_VERSION :: VersionNumber = v"0.2.0"
    #
    # Print the logo. Store the initial directory.
    #
    print(" \n \n")
    println("     Welcome to OpenMendel's")
    println("      IHT analysis option")
    println("        version ", MENDEL_IHT_VERSION)
    print(" \n \n")
    println("Reading the data.\n")
    initial_directory = pwd()
    #
    # The user specifies the analysis to perform via a set of keywords.
    # Start the keywords at their default values.
    #
    keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
    #
    # Define some keywords unique to this analysis option.
    #
    keyword["data_type"] = ""
    keyword["predictors_per_group"] = ""
    keyword["manhattan_plot_file"] = ""
    keyword["max_groups"] = ""
    keyword["group_membership"] = ""
    #
    # Process the run-time user-specified keywords that will control the analysis.
    # This will also initialize the random number generator.
    #
    process_keywords!(keyword, control_file, args)
    #
    # Check that the correct analysis option was specified.
    #
    lc_analysis_option = lowercase(keyword["analysis_option"])
    if (lc_analysis_option != "" && lc_analysis_option != "iht")
        throw(ArgumentError("An incorrect analysis option was specified.\n \n"))
    end
    keyword["analysis_option"] = "Iterative Hard Thresholding"
    @assert (keyword["max_groups"] != "")           "Need number of groups. Choose 1 to run normal IHT"
    @assert (keyword["predictors_per_group"] != "") "Need number of predictors per group"

    #
    # Read the genetic data from the external files named in the keywords.
    #
    # (pedigree, person, nuclear_family, locus, snpdata, locus_frame, phenotype_frame,
    #     pedigree_frame, snp_definition_frame) = read_external_data_files(keyword)
    #
    # Execute the specified analysis.
    #
    println(" \nAnalyzing the data.\n")
##
    file_name = keyword["plink_input_basename"]
    snpmatrix = SnpArray(file_name)
    phenotype = readdlm(file_name * ".fam", header = false)[:, 6]
    # y_copy = copy(phenotype)
    # y_copy .-= mean(y_copy)
    groups = vec(readdlm(keyword["group_membership"], Int64))
    k = keyword["predictors_per_group"]
    J = keyword["max_groups"]
    return L0_reg(snpmatrix, phenotype, J, k, groups)
##
    # execution_error = iht_gwas(person, snpdata, pedigree_frame, keyword)
    # if execution_error
    #   println(" \n \nERROR: Mendel terminated prematurely!\n")
    # else
    #   println(" \n \nMendel's analysis is finished.\n")
    # end
    # #
    # # Finish up by closing, and thus flushing, any output files.
    # # Return to the initial directory.
    # #
    # close(keyword["output_unit"])
    # cd(initial_directory)
    # return nothing
end #function MendelIHT

"""
Calculates the IHT step β+ = P_k(β - μ ∇f(β)).
Returns step size (μ), and number of times line search was done (μ_step).
This function updates: b, xb, xk, gk, xgk, idx
"""
function iht!(
    v        :: IHTVariable{T},
    x        :: SnpLike{2},
    y        :: Vector{T},
    J        :: Int,
    k        :: Int,
    mean_vec :: Vector{T},
    std_vec  :: Vector{T},
    iter     :: Int = 1,
    nstep    :: Int = 50,
) where {T <: Float}

    # compute indices of nonzeroes in beta
    v.idx .= v.b .!= 0
    if sum(v.idx) == 0
        _init_iht_indices(v, J, k)
    end

    # store relevant columns of x. Need to do this on 1st iteration.
    # afterwards, only do if support changes
    if !isequal(v.idx, v.idx0) || iter < 2
        copy!(v.xk, view(x, :, v.idx))
    end

    # store relevant components of gradient
    v.gk .= v.df[v.idx]

    # now compute subset of x*g
    SnpArrays.A_mul_B!(v.xgk, v.xk, v.gk, mean_vec[v.idx], std_vec[v.idx]) # v.df = X'(y - Xβ)

    # warn if xgk only contains zeros
    all(v.xgk .== zero(T)) && warn("Entire active set has values equal to 0")

    # compute step size
#    mu = _iht_stepsize(v, k) :: T # does not yield conformable arrays?!
    μ = (sum(abs2, v.gk) / sum(abs2, v.xgk)) :: T

    # check for finite stepsize
    isfinite(μ) || throw(error("Step size is not finite, is active set all zero?"))

    # compute gradient step
    _iht_gradstep(v, μ, J, k)

    # update xb
    v.xk .= view(x, :, v.idx)
    SnpArrays.A_mul_B!(v.xb, v.xk, v.b[v.idx], mean_vec[v.idx], std_vec[v.idx])

    # calculate omega
    ω_top, ω_bot = _iht_omega(v)

    # backtrack until mu < omega and until support stabilizes
    μ_step = 0
    while _iht_backtrack(v, ω_top, ω_bot, μ, μ_step, nstep)

        # stephalving
        μ /= 2

        # recompute gradient step
        copy!(v.b,v.b0)
        # v.itc = v.itc0
        _iht_gradstep(v, μ, J, k)

        # recompute xb
        v.xk .= view(x, :, v.idx)
        SnpArrays.A_mul_B!(v.xb, v.xk, v.b[v.idx], mean_vec[v.idx], std_vec[v.idx])

        # calculate omega
        ω_top, ω_bot = _iht_omega(v)

        # increment the counter
        μ_step += 1
    end

    return μ::T, μ_step::Int
end

"""
This function performs IHT on GWAS data.
"""
function L0_reg(
    x        :: SnpLike{2},
    y        :: Vector{T},
    J        :: Int,
    k        :: Int,
    group    :: Vector{Int};
    v        :: IHTVariable = IHTVariables(x, y, J, k),
    # v        :: IHTVariables = IHTVariables(x, y, J, k),
    mask_n   :: Vector{Int} = ones(Int, size(y)),
    tol      :: T = 1e-4,
    max_iter :: Int = 100,
    max_step :: Int = 50,
) where {T <: Float}

    # start timer
    tic()

    # first handle errors
    @assert J >= 0        "Value of J (max number of groups) must be nonnegative!\n"
    @assert k >= 0        "Value of k (max predictors per group) must be nonnegative!\n"
    @assert max_iter >= 0 "Value of max_iter must be nonnegative!\n"
    @assert max_step >= 0 "Value of max_step must be nonnegative!\n"
    @assert tol > eps(T)  "Value of global tol must exceed machine precision!\n"

    # initialize return values
    mm_iter   = 0                 # number of iterations of L0_reg
    mm_time   = 0.0               # compute time *within* L0_reg
    next_loss = oftype(tol,Inf)   # loss function value

    # initialize floats
    current_obj = oftype(tol,Inf) # tracks previous objective function value
    the_norm    = 0.0             # norm(b - b0)
    scaled_norm = 0.0             # the_norm / (norm(b0) + 1)
    μ           = 0.0             # Landweber step size, 0 < tau < 2/rho_max^2

    # initialize integers
    mu_step = 0                   # counts number of backtracking steps for mu

    # initialize booleans
    converged = false             # scaled_norm < tol?

    # compute some summary statistics for our snpmatrix
    maf, minor_allele, missings_per_snp, missings_per_person = summarize(x)
    people, snps = size(x)

    #precompute mean and standard deviations for each snp. Note that (1) the mean is
    #given by 2 * maf, and (2) based on which allele is the minor allele, might need to do
    #2.0 - the maf for the mean vector.
    mean_vec = zeros(snps)
    for i in 1:snps
        minor_allele[i] ? mean_vec[i] = 2.0 - 2.0maf[i] : mean_vec[i] = 2.0maf[i]
    end
    std_vec = std_reciprocal(x, mean_vec)

    #add intercept (at the end)
    x        = [x SnpArray(ones(people))] #NOTE this creates A LOT of extra memory!!!!
    mean_vec = [mean_vec; zero(T)]
    std_vec  = [std_vec; one(T)]

    #
    # Begin IHT calculations
    #
    fill!(v.xb, 0.0)       #initialize β = 0 vector, so Xβ = 0
    copy!(v.r, y)          #redisual = y-Xβ-intercept = y  CONSIDER BLASCOPY!
    v.r[mask_n .== 0] .= 0 #bit masking? idk why we need this yet
    v.group .= group       #assign the groups in the beginning

    # Calculate the gradient v.df = -X'(y - Xβ) = X'(-1*(Y-Xb)). All future gradient
    # calculations are done in iht!. Note the negative sign will be cancelled afterwards
    # when we do b+ = P_k( b - μ∇f(b)) = P_k( b + μ(-∇f(b))) = P_k( b + μ*v.df)
    SnpArrays.At_mul_B!(v.df, x, v.r, mean_vec, std_vec)

    for mm_iter = 1:max_iter
        # save values from previous iterate
        copy!(v.b0, v.b)   # b0 = b    CONSIDER BLASCOPY!
        copy!(v.xb0, v.xb) # Xb0 = Xb  CONSIDER BLASCOPY!
        # v.itc0 = v.itc     # update intercept as well
        loss = next_loss

        #calculate the step size μ. TODO: check how adding intercept affects this
        (μ, μ_step) = iht!(v, x, y, J, k, mean_vec, std_vec, max_step, mm_iter)

        # iht! gives us an updated x*b. Use it to recompute residuals and gradient
        # v.r .= y .- v.xb - v.itc
        v.r .= y .- v.xb
        v.r[mask_n .== 0] .= 0 #bit masking, idk why we need this yet

        # v.df = X'(y - Xβ - intercept)
        SnpArrays.At_mul_B!(v.df, x, v.r, mean_vec, std_vec, similar(v.df))

        # update loss, objective, gradient, and check objective is not NaN or Inf
        next_loss = sum(abs2, v.r) / 2
        !isnan(next_loss) || throw(error("Objective function is NaN, aborting..."))
        !isinf(next_loss) || throw(error("Objective function is Inf, aborting..."))

        # track convergence
        the_norm    = chebyshev(v.b, v.b0) #max(abs(x - y))
        scaled_norm = the_norm / (norm(v.b0, Inf) + 1)
        converged   = scaled_norm < tol

        if converged
            mm_time = toq()   # stop time
            return gIHTResults(mm_time, next_loss, mm_iter, copy(v.b), J, k, group)
        end

        if mm_iter == max_iter
            mm_time = toq() # stop time
            throw(error("Did not converge!!!!! The run time for IHT was " *
                string(mm_time) * "seconds"))
        end
    end
end #function L0_reg

#"""
#Calculates the IHT step β+ = P_k(β - μ ∇f(β)).
#Returns step size (μ), and number of times line search was done (μ_step).
#
#This function updates: b, xb, xk, gk, xgk, idx
#"""
#function iht!(
#    v         :: IHTVariable,
#    snpmatrix :: Matrix{Float64},
#    y         :: Vector{Float64},
#    k         :: Int;
#    iter      :: Int = 1,
#    nstep     :: Int = 50,
#)
#    # compute indices of nonzeroes in beta and store them in v.idx (also sets size of v.gk)
#    _iht_indices(v, k)
#
#    # fill v.xk, which stores columns of snpmatrix corresponding to non-0's of b
#    v.xk[:, :] .= snpmatrix[:, v.idx]
#
#    # fill v.gk, which store only k largest components of gradient (v.df)
#    # fill_perm!(v.gk, v.df, v.idx)  # gk = g[v.idx]
#    v.gk .= v.df[v.idx]
#
#    # now compute X_k β_k and store result in v.xgk
#    A_mul_B!(v.xgk, v.xk, v.gk)
#
#    # warn if xgk only contains zeros
#    all(v.xgk .≈ 0.0) && warn("Entire active set has values equal to 0")
#
#    #compute step size and notify if step size too small
#    #μ = norm(v.gk, 2)^2 / norm(v.xgk, 2)^2
#    μ = sum(abs2, v.gk) / sum(abs2, v.xgk)
#    isfinite(μ) || throw(error("Step size is not finite, is active set all zero?"))
#    μ <= eps(typeof(μ)) && warn("Step size $(μ) is below machine precision, algorithm may not converge correctly")
#
#    #Take the gradient step and compute ω. Note in order to compute ω, need β^{m+1} and xβ^{m+1} (eq5)
#    _iht_gradstep(v, μ, k)
#    ω = compute_ω!(v, snpmatrix) #is snpmatrix required? Or can I just use v.x
#
#    #compute ω and check if μ < ω. If not, do line search by halving μ and checking again.
#    μ_step = 0
#    for i = 1:nstep
#        #exit loop if μ < ω where c = 0.01 for now
#        if _iht_backtrack(v, ω, μ); break; end
#
#        #if μ >= ω, step half and warn if μ falls below machine epsilon
#        μ /= 2
#        μ <= eps(typeof(μ)) && warn("Step size equals zero, algorithm may not converge correctly")
#
#        # recompute gradient step
#        copy!(v.b, v.b0)
#        _iht_gradstep(v, μ, k)
#
#        # re-compute ω based on xβ^{m+1}
#        A_mul_B!(v.xb, snpmatrix, v.b)
#        ω = sqeuclidean(v.b, v.b0) / sqeuclidean(v.xb, v.xb0)
#
#        μ_step += 1
#    end
#
#    return (μ, μ_step)
#end
#

end #module


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1], URLs[2]])

dec = 64
modulus = 4800
base_frequency = 125000000
samples_per_period = div(modulus, dec)
periods_per_frame = 3

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)

modeDAC(rp, "STANDARD")
frequencyDAC(rp, 1, 1, base_frequency / modulus)

println(" frequency = $(frequencyDAC(rp,1,1))")
amplitudeDAC(rp, 1, 1, 4000)
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

triggerMode(rp, "EXTERNAL")
ramWriterMode(rp, "TRIGGERED")
masterTrigger(rp, false)

startADC(rp)
masterTrigger(rp, true)

sleep(1.0)

uFirstPeriod = readData(rp, 0, 1)
uCurrentPeriod = readData(rp, currentFrame(rp), 1)
#RedPitayaDAQServer.disconnect(rp)

figure(1)
clf()
subplot(1,2,1)
plot(vec(uFirstPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,1,:,:]))
subplot(1,2,2)
plot(vec(uFirstPeriod[:,3,:,:]))
plot(vec(uCurrentPeriod[:,3,:,:]))
legend(("first period", "current period"))

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
#=
Converted files 04.04.2021:
- attribute.feature
- entity.feature

=#

# from attribute.feature -> not implemented in python
@given("put entity type: person") do context
    @fail "Implement me"
end

# implemented in thing_steps.py but not in java
@when("delete attribute: \$x") do context
    @fail "Implement me"
end


text:
const MAX_RAND_PAD_LENGTH_MOD8 = 3

randpad() = rand(UInt8, rand(8*(1:MAX_RAND_PAD_LENGTH_MOD8)))

function rand_primitive_buffer(;lpad=randpad(), rpad=randpad())
    dtype = rand(PRIMITIVE_ELTYPES)
    len = rand(1:MAX_VECTOR_LENGTH)
    v = rand(dtype, len)
    b = convert(Vector{UInt8}, reinterpret(UInt8, v))
    b = vcat(lpad, b, rpad)
    len, dtype, length(lpad), b, v
end

function rand_nullableprimitive_buffer(;lpad=randpad(), rpad=randpad())
    dtype = rand(PRIMITIVE_ELTYPES)
    len = rand(1:MAX_VECTOR_LENGTH)
    pres = rand(Bool, len)
    mask = Arrow.bitpackpadded(pres)
    vraw = rand(dtype, len)
    b = vcat(lpad, mask, convert(Vector{UInt8}, reinterpret(UInt8, vraw)), rpad)
    v = Union{dtype,Missing}[pres[i] ? vraw[i] : missing for i ∈ 1:len]
    len, dtype, length(mask), length(lpad), b, v
end


text:
using Documenter, Econometrics

makedocs(modules = [Econometrics], sitename = "Econometrics.jl")
deploydocs(repo = "github.com/Nosferican/Econometrics.jl.git")


text:
"""
$(TYPEDEF)

Defines interpolating DenseHamiltonian object

# Fields

$(FIELDS)
"""
struct InterpDenseHamiltonian{T} <: AbstractDenseHamiltonian{T}
    "Interpolating object"
    interp_obj::Any
    "Size"
    size::Any
end

"""
$(TYPEDEF)

Defines interpolating SparseHamiltonian object

# Fields

$(FIELDS)
"""
struct InterpSparseHamiltonian{T} <: AbstractSparseHamiltonian{T}
    "Interpolating object"
    interp_obj::Any
    "Size"
    size::Any
end

function InterpDenseHamiltonian(
    s,
    hmat;
    method = "bspline",
    order = 1,
    unit = :h,
)
    if ndims(hmat) == 3
        hsize = size(hmat)[1:2]
        htype = eltype(hmat)
    elseif ndims(hmat) == 1
        hsize = size(hmat[1])
        htype = eltype(sum(hmat))
    else
        throw(ArgumentError("Invalid input data dimension."))
    end

    interp_obj = construct_interpolations(
        s,
        unit_scale(unit) * hmat,
        method = method,
        order = order,
    )
    InterpDenseHamiltonian{htype}(interp_obj, hsize)
end

function (H::InterpDenseHamiltonian)(s)
    H.interp_obj(1:size(H, 1), 1:size(H, 1), s)
end

# The argument `p` is not essential for `InterpDenseHamiltonian`
# It exists to keep the `update_cache!` interface consistent across
# all `AbstractHamiltonian` types
function update_cache!(cache, H::InterpDenseHamiltonian, tf, s::Real)
    for i = 1:size(H, 1)
        for j = 1:size(H, 1)
            @inbounds cache[i, j] = -1.0im * H.interp_obj(i, j, s)
        end
    end
end

function update_vectorized_cache!(cache, H::InterpDenseHamiltonian, tf, s)
    hmat = H(s)
    iden = Matrix{eltype(H)}(I, size(H))
    cache .= 1.0im * (transpose(hmat) ⊗ iden - iden ⊗ hmat)
end

function get_cache(H::InterpDenseHamiltonian{T}, vectorize) where {T}
    if vectorize == true
        hsize = size(H, 1) * size(H, 1)
        Matrix{T}(undef, hsize, hsize)
    else
        Matrix{T}(undef, size(H))
    end
end

get_cache(H::InterpDenseHamiltonian{T}) where {T} = Matrix{T}(undef, size(H))

function (h::InterpDenseHamiltonian)(
    du,
    u::Matrix{T},
    p,
    t::Real,
) where {T<:Complex}
    fill!(du, 0.0 + 0.0im)
    H = h(t)
    gemm!('N', 'N', -1.0im, H, u, 1.0 + 0.0im, du)
    gemm!('N', 'N', 1.0im, u, H, 1.0 + 0.0im, du)
end

function InterpSparseHamiltonian(
    s_axis,
    H_list::AbstractArray{SparseMatrixCSC{T,Int},1};
    unit = :h,
) where {T<:Number}
    interp_obj = construct_interpolations(
        collect(s_axis),
        unit_scale(unit) * H_list,
        method = "gridded",
        order = 1,
    )
    InterpSparseHamiltonian{T}(interp_obj, size(H_list[1]))
end

function (H::InterpSparseHamiltonian)(s)
    H.interp_obj(s)
end

function get_cache(H::InterpSparseHamiltonian{T}, vectorize) where {T}
    if vectorize == true
        hsize = size(H, 1) * size(H, 1)
        spzeros(T, hsize, hsize)
    else
        spzeros(T, size(H)...)
    end
end

get_cache(H::InterpSparseHamiltonian{T}) where {T} = spzeros(T, size(H)...)

# The argument `p` is not essential for `InterpSparseHamiltonian`
# It exists to keep the `update_cache!` interface consistent across
# all `AbstractHamiltonian` types
update_cache!(cache, H::InterpSparseHamiltonian, p, s::Real) =
    cache .= -1.0im * H(s)

function update_vectorized_cache!(
    cache,
    H::InterpSparseHamiltonian,
    tf,
    s::Real,
)
    hmat = H(s)
    iden = sparse(I, size(H))
    cache .= 1.0im * (transpose(hmat) ⊗ iden - iden ⊗ hmat)
end

function (h::InterpSparseHamiltonian)(
    du,
    u::Matrix{T},
    tf,
    s::Real,
) where {T<:Number}
    H = h(s)
    du .= -1.0im * (H * u - u * H)
end


text:
"""
	grdcut(cmd0::String="", arg1=[], kwargs...)

Produce a new outgrid file which is a subregion of ingrid. The subregion is specified with
``limits`` (the -R); the specified range must not exceed the range of ingrid (but see ``extend``).

Full option list at [`grdcut`]($(GMTdoc)grdcut.html)

Parameters
----------

- **F** | **clip** | **cutline** :: [Type => Str | GMTdaset | Mx2 array | NamedTuple]	`Arg = array|fname[+c] | (polygon=Array|Str, crop2cutline=Bool, invert=Bool)`

    Specify a closed polygon (either a file or a dataset). All grid nodes outside the
    polygon will be set to NaN (>= GMT6.2).
    ($(GMTdoc)grdcut.html#f)
- **G** | **outgrid** | **outfile** | **save** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdcut(....) form.
    ($(GMTdoc)grdcut.html#g)
- **img** | **usegdal** :: [Type => Any]

    Force the cut operation to be done by GDAL. Works for images where GMT fails or even crash.
- $(GMT.opt_J)
- **N** | **extend** :: [Type => Str or []]

    Allow grid to be extended if new region exceeds existing boundaries. Append nodata value
    to initialize nodes outside current region [Default is NaN].
    ($(GMTdoc)grdcut.html#n)
- $(GMT.opt_R)
- **S** | **circ_subregion** :: [Type => Str]    ``Arg = [n]lon/lat/radius[unit]``

    Specify an origin and radius; append a distance unit and we determine the corresponding
    rectangular region so that all grid nodes on or inside the circle are contained in the subset.
    ($(GMTdoc)grdcut.html#s)
- $(GMT.opt_V)
- **Z** | **z_subregion** :: [Type => Str]       ``Arg = [n|N |r][min/max]``

    Determine a new rectangular region so that all nodes outside this region are also outside
    the given z-range.
    ($(GMTdoc)grdcut.html#z)
- $(GMT.opt_f)
"""
function grdcut(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdcut", cmd0, arg1)

	arg2 = nothing
	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
    cmd, = parse_common_opts(d, "", [:R :V_params :f])
    opt_J, = parse_J(d, "")
    if (!startswith(opt_J, " -JX"))  cmd *= opt_J  end
	cmd = parse_these_opts(cmd, d, [[:G :outgrid :outfile :save], [:N :extend], [:S :circ_subregion], [:Z :z_subregion]])
	cmd, args, n, = add_opt(d, cmd, 'F', [:F :clip :cutline], :polygon, Array{Any,1}([arg1, arg2]),
	                        (crop2cutline="_+c", invert="_+i"))
	if (n > 0)  arg1, arg2 = args[:]  end
	(show_kwargs[1]) && return print_kwarg_opts([:img :usegdal], "Any")		# Just print the options

	if (cmd0 != "" && ((find_in_dict(d, [:img :usegdal])[1]) !== nothing))
		(cmd0[1] == '@') && (cmd0 = gmtwhich(cmd0)[1].text[1])	# A remote file
		ds = Gdal.read(cmd0)
		t = split(scan_opt(cmd, "-R"), '/')
		opts = ["-projwin", t[1], t[4], t[2], t[3]]		# -projwin <ulx> <uly> <lrx> <lry>
		if ((outname = scan_opt(cmd, "-G")) == "")
			gd2gmt(gdaltranslate(ds, opts))
		else
			gdaltranslate(ds, opts; dest=outname)
			return nothing				# Since it wrote a file so nothing to return
		end
	else
		common_grd(d, cmd0, cmd, "grdcut ", arg1, arg2)		# Finish build cmd and run it
	end
end

# ---------------------------------------------------------------------------------------------------
grdcut(arg1, cmd0::String=""; kw...) = grdcut(cmd0, arg1; kw...)


text:
using BasisFunctions
using MoreBasisFunctions
using DomainSets
using FrameFun

f(x) = cos(π*x)

L = LagrangeGaussLegendre(21)
FL = Fun(f, L, support(L))
println(FL(0.0))

B = Bernstein(21)
FB = Fun(f, B, support(B))
println(FB(0.5))


text:
BUCH_DISCARDED = 0
BUCH_NF        = 0
BUCH_SIZE      = 0
BUCH_REDUCED   = 0

function zerobuch()
    global BUCH_DISCARDED
    global BUCH_NF
    global BUCH_SIZE
    global BUCH_REDUCED
    BUCH_DISCARDED = 0
    BUCH_NF    = 0
    BUCH_SIZE      = 0
    BUCH_REDUCED   = 0
end

#-----------------------------------------------------------------------------
# Regular reductions

# try to reduce f with g in a regular way and return result
function reduction_step(f, g)
    @info "reducing $(f) with $(g)" f g

    for t in terms(f)
        # if divides
        success, u = divides(t, leading_monomial(g))
        if success
            u = divexact(t, leading_term(g))
            @info "$(leading_term(g)) | $t"

            newf = f - u*g
            @info "reduction! $(f) - $(u)*$(g) --> $newf"
            return true, newf
        end
    end
    return false, f
end

# // same as above, but in terms of set G //
# try to reduce f with G in a regular way and return result
function reduction_step(f, G::Vector{T}) where {T}
    for g in G
        success, newf = reduction_step(f, g)
        if success
            return true, newf
        end
    end
    return false, f
end

# regular normal form of f w.r.t. G
function normal_form(f, G)
    @info "computing normalform of $f w.r.t $G.."

    success = true
    newf = f
    while success
        success, newf = reduction_step(newf, G)
        @info "reduction $success"
    end
    newf
end

#-----------------------------------------------------------------------------
# S-polynomial

# returns u, v,  such that
# u*a - v*b is spoly(a, b)
function mults(a, b)
    u = lcm(leading_monomial(a), leading_monomial(b))
    u = divexact(u, leading_term(a))

    v = lcm(leading_monomial(a), leading_monomial(b))
    v = divexact(v, leading_term(b))

    u, v
end

# S-polynomial of a and b
function spoly(a, b)
    u, v = mults(a, b)
    u*a - v*b
end

function spoly(ab)
    spoly(ab...)
end

#-----------------------------------------------------------------------------

function buchberger_criterion1(f)
    f1, f2 = f
    if isone(gcd(leading_monomial(f1), leading_monomial(f2)))
        @warn "DISCARDED $((f1, f2)) by 1st Criterion"
        global BUCH_DISCARDED
        BUCH_DISCARDED += 1
        return true
    end
    return false
end

#-----------------------------------------------------------------------------
# Groebner basis things

function groebner_basis(F)
    G = deepcopy(F)

    P = []
    for i in 1:length(G)
        for j in i+1:length(G)

            if buchberger_criterion1((G[i], G[j]))
                continue
            end

            push!(P, (G[i], G[j]))
        end
    end

    @info "generated initial G and P:"
    @info "F = $F"
    @info "G = $G"
    @info "P = $P"

    while !isempty(P)
        sort!(P, by=leading_monomial ∘ spoly)
        f = popfirst!(P)
        @info "selected $f"

        global BUCH_NF
        BUCH_NF += 1
        fNF = normal_form(spoly(f), G)

        if iszero(fNF)
            global BUCH_REDUCED
            BUCH_REDUCED += 1
        else
            for fj in G
                if buchberger_criterion1((fNF, fj)) || iszero(spoly(fNF, fj))
                    continue
                end

                push!(P, (fNF, fj))
                @info "$fNF  $fj"
                @warn "SPOLY $(last(P))"
            end

            # update G
            fNFnormalzed = fNF * inv(leading_coefficient(fNF))
            push!(G, fNFnormalzed)
        end

        @info "updated G and P"
        @info "G = $G"
        @info "P = $P"
    end

    global BUCH_SIZE
    BUCH_SIZE = length(G) - length(F)

    global BUCH_REDUCED
    global BUCH_DISCARDED

    G
end

#-----------------------------------------------------------------------------

R, (x,y,z) = PolynomialRing(QQ, ["x","y", "z"], ordering=:degrevlex)

F = [x + y + z, x*y + y*z + x*z, x*y*z - 1]

# evaluation
G = groebner_basis(F)
println("############")
println(G)


text:
function convert_func(op, T1, W1, T2, W2 = W1)
    typ1 = LLVM_TYPES[T1]
    typ2 = LLVM_TYPES[T2]
    vtyp1 = vtype(W1, typ1)
    vtyp2 = vtype(W2, typ2)
    instrs = """
    %res = $op $vtyp2 %0 to $vtyp1
    ret $vtyp1 %res
    """
    quote
        $(Expr(:meta, :inline))
        Vec(llvmcall($instrs, _Vec{$W1,$T1}, Tuple{_Vec{$W2,$T2}}, data(v)))
    end
end
function identity_func(W, T1, T2)
    vtyp1 = vtype(W, LLVM_TYPES[T1])
    instrs = """
    ret $vtyp1 %0
    """
    quote
        $(Expr(:meta, :inline))
        Vec(llvmcall($instrs, _Vec{$W,$T1}, Tuple{_Vec{$W,$T2}}, data(v)))
    end
end

@generated function Vec{W,F}(v::Vec{W,T}) where {W,F<:Union{Float32,Float64},T<:Integer}
    convert_func(T <: Signed ? "sitofp" : "uitofp", F, W, T)
end

@generated function Vec{W,T}(v::Vec{W,F}) where {W,F<:Union{Float32,Float64},T<:Integer}
    convert_func(T <: Signed ? "fptosi" : "fptoui", T, W, F)
end
@generated function Vec{W,T1}(v::Vec{W,T2}) where {W,T1<:Integer,T2<:Integer}
    sz1 = sizeof(T1)::Int; sz2 = sizeof(T2)::Int
    if sz1 < sz2
        convert_func("trunc", T1, W, T2)
    elseif sz1 == sz2
        identity_func(W, T1, T2)
    else
        convert_func(((T1 <: Signed) && (T2 <: Signed)) ? "sext" : "zext", T1, W, T2)
    end
end

@inline Base.float(v::Vec{W,I}) where {W, I <: Union{UInt64, Int64}} = Vec{W,Float64}(v)
@generated function Base.float(v::Vec{W,I}) where {W, I <: Integer}
    ex = if 8W ≤ REGISTER_SIZE
        :(Vec{$W,Float64}(v))
    else
        :(Vec{$W,Float32}(v))
    end
    Expr(:block, Expr(:meta, :inline), ex)
end
@inline Base.float(v::AbstractSIMD{W,T}) where {W,T <: Union{Float32,Float64}} = v
@inline Base.float(vu::VecUnroll) = VecUnroll(fmap(float, vu.data))

    # for bytes2 in (1,2,4)
    #     bytes2 == bytes1 && break
    #     llvmint2 = "i$(8bytes2)"
    #     jint2 = Symbol(:Int, 8bytes2)
    #     juint2 = Symbol(:UInt, 8bytes2)
    #     expr = :(convert_func("trunc", W, $llvmint2, $jint2, $llvmint, $jint))
    #     @eval @generated Vec{W,$jint2}(v::Vec{W,$jint}) where {W} = $expr
    #     expr = :(convert_func("trunc", W, $llvmint2, $juint2, $llvmint, $jint))
    #     @eval @generated Vec{W,$juint2}(v::Vec{W,$jint}) where {W} = $expr
    #     expr = :(convert_func("trunc", W, $llvmint2, $jint2, $llvmint, $juint))
    #     @eval @generated Vec{W,$jint2}(v::Vec{W,$juint}) where {W} = $expr
    #     expr = :(convert_func("trunc", W, $llvmint2, $juint2, $llvmint, $juint))
    #     @eval @generated Vec{W,$juint2}(v::Vec{W,$juint}) where {W} = $expr

    #     expr = :(convert_func("sext", W, $llvmint, $jint, $llvmint2, $jint2))
    #     @eval @generated Vec{W,$jint}(v::Vec{W,$jint2}) where {W} = $expr
    #     expr = :(convert_func("zext", W, $llvmint, $juint, $llvmint2, $jint2))
    #     @eval @generated Vec{W,$juint}(v::Vec{W,$jint2}) where {W} = $expr
    #     expr = :(convert_func("zext", W, $llvmint, $jint, $llvmint2, $juint2))
    #     @eval @generated Vec{W,$jint}(v::Vec{W,$juint2}) where {W} = $expr
    #     expr = :(convert_func("zext", W, $llvmint, $juint, $llvmint2, $juint2))
    #     @eval @generated Vec{W,$juint}(v::Vec{W,$juint2}) where {W} = $expr
    # end
    # expr = :(identity_func(W, $llvmint, $juint, $jint))
    # @eval @generated Vec{W,$juint}(v::Vec{W,$jint}) where {W} = $expr
    # expr = :(identity_func(W, $llvmint, $jint, $juint))
    # @eval @generated Vec{W,$jint}(v::Vec{W,$juint}) where {W} = $expr
# end

@generated function Vec{W,Float32}(v::Vec{W,Float64}) where {W}
    convert_func("fptrunc", Float32, W, Float64, W)
end
@generated function Vec{W,Float64}(v::Vec{W,Float32}) where {W}
    convert_func("fpext", Float64, W, Float32, W)
end

@inline Base.convert(::Vec{W,T}, s::T) where {W, T <: NativeTypes} = vbroadcast(Val{W}(), s)
@inline Base.convert(::Vec{W,T}, s::T) where {W, T <: Integer} = vbroadcast(Val{W}(), s)
@inline Base.convert(::Vec{W,T}, s::NativeTypes) where {W, T} = vbroadcast(Val{W}(), T(s))
@inline Base.convert(::Vec{W,T1}, s::T2) where {W, T1 <: Integer, T2 <: Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64}} = vbroadcast(Val{W}(), s % T1)
@inline function Base.convert(::Type{T}, v::Vec{W,S}) where {T<:Number,S,W}
    if S <: T
        v
    else
        Vec{W,T}(v)
    end
end
@inline function Base.convert(::Type{U}, v::Vec{W,S}) where {N,W,T,U<:VecUnroll{N,W,T},S}
    VecUnroll{N}(convert(T, v))
end
@inline Base.convert(::Type{Vec{W,T}}, v::Vec{W,S}) where {T<:Number,S,W} = Vec{W,T}(v)
@inline Base.unsafe_trunc(::Type{I}, v::Vec{W,T}) where {W,I,T} = convert(Vec{W,I}, v)
@inline Base.:(%)(v::Vec{W,T}, ::Type{I}) where {W,I,T} = convert(Vec{W,I}, v)

@inline Base.convert(::Type{T}, v::Union{Mask,VecUnroll{<:Any, <:Any, Bool, <: Mask}}) where {T <: Union{Base.HWReal,Bool}} = ifelse(v, one(T), zero(T))
@inline Base.convert(::Type{<:AbstractSIMD{W,T}}, v::Union{Mask{W},VecUnroll{<:Any, W, Bool, <: Mask}}) where {W, T <: Union{Base.HWReal,Bool}} = ifelse(v, one(T), zero(T))
@inline Base.convert(::Type{Bit}, v::Union{Mask,VecUnroll{<:Any, <:Any, Bool, <: Mask}}) = v
@inline Base.convert(::Type{<:AbstractSIMD{W,Bit}}, v::Union{Mask{W},VecUnroll{<:Any, W, Bool, <: Mask}}) where {W} = v
# @inline Base.convert(::typeof{T}, v::T)

@generated function Base.reinterpret(::Type{T1}, v::Vec{W2,T2}) where {W2, T1 <: NativeTypes, T2}
    W1 = W2 * sizeof(T2) ÷ sizeof(T1)
    Expr(:block, Expr(:meta, :inline), :(reinterpret(Vec{$W1,$T1}, v)))
end
@generated function Base.reinterpret(::Type{Vec{W1,T1}}, v::Vec{W2,T2}) where {W1, W2, T1, T2}
    @assert sizeof(T1) * W1 == W2 * sizeof(T2)
    convert_func("bitcast", T1, W1, T2, W2)
end


text:
AbstractTrees.children(x::AbstractXMLNode) = collect(child_elements(x));
AbstractTrees.printnode(io::IO, x::AbstractXMLNode) = print(io, name(x));

#Processes XML codebook tree into a usable JSON table
function process_codebook()
    xdoc = parse_file("PSIDCodebook.xml");
    r = root(xdoc);
    t = Tree(r)
    # Travel down the tree until we reach the list of variables
    name(t[1])
    current_node = children(r)[1]
    name(current_node)
    current_node = children(current_node)[1]
    name(current_node)
    current_node = children(current_node)[2]
    name(current_node)
    ce = current_node["VARIABLE"]
    # Ingest the variables into a DataFrame
    outdf = DataFrame()
    outdf.NAME =  [content(ce[i]["NAME"][1]) for i in 1:length(ce)]
    outdf.YEAR =  [content(ce[i]["YEAR"][1]) for i in 1:length(ce)]
    outdf.QTEXT =  [content(ce[i]["QTEXT"][1]) for i in 1:length(ce)]
    outdf.ETEXT =  [content(ce[i]["ETEXT"][1]) for i in 1:length(ce)]
    outdf.TYPE_ID =  [content(ce[i]["TYPE_ID"][1]) for i in 1:length(ce)]
    outdf.LABEL =  [content(ce[i]["LABEL"][1]) for i in 1:length(ce)]

    # codes
    list_codes = [ce[i]["LIST_CODE"][1] for i in 1:length(ce)]

    list_codes[1]

    # Take a codexml vector of some length containing (value, text) pairs
    # Return a dict
    function process_codes(codexml)
        codes = codexml["CODE"]
        vals = [content(c["VALUE"][1]) for c in codes]
        texts = [content(c["TEXT"][1]) for c in codes]
        Dict(v => t for (v, t) in zip(vals, texts))
    end
    process_codes(list_codes[1])

    morecodes = process_codes.(list_codes)
    outdf.codedict = morecodes

    #j = objecttable(outdf)
    j = arraytable(outdf)
    write("output/codebook.json", j)
end


text:
"""
    tear_graph(sys) -> sys

Tear the bipartite graph in a system.
"""
function tear_graph(sys)
    find_solvables!(sys)
    s = structure(sys)
    @unpack graph, solvable_graph, assign, inv_assign, scc = s

    partitions = map(scc) do c
        ieqs = filter(eq->isalgeq(s, eq), c)
        vars = inv_assign[ieqs]

        td = TraverseDAG(graph.fadjlist, length(assign))
        e_solved, v_solved, e_residue, v_tear = tearEquations!(td, solvable_graph.fadjlist, ieqs, vars)
    end

    @set! sys.structure.partitions = partitions
    return sys
end

function tearing_sub(expr, dict, s)
    expr = ModelingToolkit.fixpoint_sub(expr, dict)
    s ? simplify(expr) : expr
end

function tearing_reassemble(sys; simplify=false)
    s = structure(sys)
    @unpack fullvars, partitions, assign, inv_assign, graph, scc = s
    eqs = equations(sys)

    ### extract partition information
    rhss = []
    solvars = []
    ns, nd = nsrcs(graph), ndsts(graph)
    active_eqs  = trues(ns)
    active_vars = trues(nd)
    rvar2req = Vector{Int}(undef, nd)
    for (ith_scc, (e_solved, v_solved, e_residue, v_tear)) in enumerate(partitions)
        for ii in eachindex(e_solved)
            ieq = e_solved[ii]; ns -= 1
            iv = v_solved[ii]; nd -= 1
            rvar2req[iv] = ieq

            active_eqs[ieq] = false
            active_vars[iv] = false

            eq = eqs[ieq]
            var = fullvars[iv]
            rhs = solve_for(eq, var; simplify=simplify, check=false)
            # if we don't simplify the rhs and the `eq` is not solved properly
            (!simplify && var in vars(rhs)) && (rhs = SymbolicUtils.polynormalize(rhs))
            # Since we know `eq` is linear wrt `var`, so the round off must be a
            # linear term. We can correct the round off error by a linear
            # correction.
            rhs -= expand_derivatives(Differential(var)(rhs))*var
            @assert !(var in vars(rhs)) """
            When solving
            $eq
            $var remainded in
            $rhs.
            """
            push!(rhss, rhs)
            push!(solvars, var)
        end
        # DEBUG:
        #@show ith_scc solvars .~ rhss
        #Main._nlsys[] = eqs[e_solved], fullvars[v_solved]
        #ModelingToolkit.topsort_equations(solvars .~ rhss, fullvars)
        #empty!(solvars); empty!(rhss)
    end

    ### update SCC
    eq_reidx = Vector{Int}(undef, nsrcs(graph))
    idx = 0
    for (i, active) in enumerate(active_eqs)
        eq_reidx[i] = active ? (idx += 1) : -1
    end

    rmidxs = Int[]
    newscc = Vector{Int}[]; sizehint!(newscc, length(scc))
    for component′ in newscc
        component = copy(component′)
        for (idx, eq) in enumerate(component)
            if active_eqs[eq]
                component[idx] = eq_reidx[eq]
            else
                push!(rmidxs, idx)
            end
        end
        push!(newscc, component)
        deleteat!(component, rmidxs)
        empty!(rmidxs)
    end

    ### update graph
    var_reidx = Vector{Int}(undef, ndsts(graph))
    idx = 0
    for (i, active) in enumerate(active_vars)
        var_reidx[i] = active ? (idx += 1) : -1
    end

    newgraph = BipartiteGraph(ns, nd)

    function visit!(ii, gidx, basecase=true)
        ieq = basecase ? ii : rvar2req[ii]
        for ivar in 𝑠neighbors(graph, ieq)
            # Note that we need to check `ii` against the rhs states to make
            # sure we don't run in circles.
            (!basecase && ivar === ii) && continue
            if active_vars[ivar]
                add_edge!(newgraph, gidx, var_reidx[ivar])
            else
                # If a state is reduced, then we go to the rhs and collect
                # its states.
                visit!(ivar, gidx, false)
            end
        end
        return nothing
    end

    ### update equations
    newstates = setdiff([fullvars[diffvars_range(s)]; fullvars[algvars_range(s)]], solvars)
    varidxmap = Dict(newstates .=> 1:length(newstates))
    neweqs = Vector{Equation}(undef, ns)
    newalgeqs = falses(ns)

    dict = Dict(value.(solvars) .=> value.(rhss))

    for ieq in Iterators.flatten(scc); active_eqs[ieq] || continue
        eq = eqs[ieq]
        ridx = eq_reidx[ieq]

        visit!(ieq, ridx)

        if isdiffeq(eq)
            neweqs[ridx] = eq.lhs ~ tearing_sub(eq.rhs, dict, simplify)
        else
            newalgeqs[ridx] = true
            if !(eq.lhs isa Number && eq.lhs != 0)
                eq = 0 ~ eq.rhs - eq.lhs
            end
            rhs = tearing_sub(eq.rhs, dict, simplify)
            if rhs isa Symbolic
                neweqs[ridx] = 0 ~ rhs
            else # a number
                if abs(rhs) > 100eps(float(rhs))
                    @warn "The equation $eq is not consistent. It simplifed to 0 == $rhs."
                end
                neweqs[ridx] = 0 ~ fullvars[inv_assign[ieq]]
            end
        end
    end

    ### update partitions
    newpartitions = similar(partitions, 0)
    emptyintvec = Int[]
    for ii in eachindex(partitions)
        _, _, og_e_residue, og_v_tear = partitions[ii]
        isempty(og_v_tear) && continue
        e_residue = similar(og_e_residue)
        v_tear = similar(og_v_tear)
        for ii in eachindex(og_e_residue)
            e_residue[ii] = eq_reidx[og_e_residue[ii]]
            v_tear[ii] = var_reidx[og_v_tear[ii]]
        end
        # `emptyintvec` is aliased to save memory
        # We need them for type stability
        push!(newpartitions, (emptyintvec, emptyintvec, e_residue, v_tear))
    end

    obseqs = solvars .~ rhss

    @set! sys.structure.graph = newgraph
    @set! sys.structure.scc = newscc
    @set! sys.structure.fullvars = fullvars[active_vars]
    @set! sys.structure.partitions = newpartitions
    @set! sys.structure.algeqs = newalgeqs
    @set! sys.eqs = neweqs
    @set! sys.states = newstates
    @set! sys.reduced_states = [get_reduced_states(sys); solvars]
    @set! sys.observed = vcat(observed(sys), obseqs)
    return sys
end

"""
    algebraic_equations_scc(sys)

Find strongly connected components of algebraic equations in a system.
"""
function algebraic_equations_scc(sys)
    s = get_structure(sys)
    if !(s isa SystemStructure)
        sys = initialize_system_structure(sys)
        s = structure(sys)
    end

    # skip over differential equations
    algvars = isalgvar.(Ref(s), 1:ndsts(s.graph))
    eqs = equations(sys)
    assign = matching(s, algvars, s.algeqs)

    components = find_scc(s.graph, assign)
    inv_assign = inverse_mapping(assign)

    @set! sys.structure.assign = assign
    @set! sys.structure.inv_assign = inv_assign
    @set! sys.structure.scc = components
    return sys
end

"""
    tearing(sys; simplify=false)

Tear the nonlinear equations in system. When `simplify=true`, we simplify the
new residual residual equations after tearing.
"""
tearing(sys; simplify=false) = tearing_reassemble(tear_graph(algebraic_equations_scc(sys)); simplify=simplify)


text:
"""
play_sound(filename::String, loops::Integer)

Plays a sound effect from the `sounds` subdirctory. It will play the specified number of times. If not specified, it will default to once.
"""
function play_sound(name, loops=0, ticks=-1)
    sound_file = file_path(name, :sounds)
    sample=Mix_LoadWAV_RW(SDL_RWFromFile(sound_file, "rb"), 1);
    if sample == C_NULL
        @warn "Could not load sound file: $sound_file\n$(getSDLError())"
        return
    end
    r = Mix_PlayChannelTimed(Int32(-1), sample, loops, ticks)
    if r == -1
        @warn "Unable to play sound $sound_file\n$(getSDLError())"
    end
end

"""
play_music(name::String, loops::Integer)

Plays music from the `sounds` subdirectory. It will play the file the specified number of times. If not specified, it will default to infinitely.
"""
function play_music(name, loops=-1)
    music_file = file_path(name, :music)
    music = Mix_LoadMUS(music_file)
    Mix_PlayMusic( music, Int32(loops) )
end

const resource_ext = Dict(
        :images=>"[png|jpg|jpeg]",
        :sounds=>"[mp3|ogg|wav]",
        :music=>"[mp3|ogg|wav]",
        :fonts=>"[ttf]")

function image_surface(image::String)
    image_file = file_path(image, :images)
    sf = IMG_Load(image_file)
    if sf == C_NULL
        throw("Error loading $image_file")
    end
    return sf
end

function file_path(name::String, subdir::Symbol)
    path = joinpath(game[].location, String(subdir))
    @assert isdir(path)
    allfiles = readdir(path)
    allexts = resource_ext[subdir]
    validate_name(name)
    for x in allfiles
        if occursin(Regex("$(name)(\\.$(allexts))?", "i"), x)
            return joinpath(path, x)
        end
    end
    # We try to return helpful messages if the file could not be found
    for x in allfiles
        if basename(x) == name
            @warn "Did you mean $x? We can only handle the follwing extensions: $allexts"
        end
        if edit_distance(x, name) / length(name) <= .5
            @warn "Did you mean $x instead of $names. Please check your spelling."
        end
    throw(ArgumentError("No file: $name in $path")); end

end

"""
Simplistic string edit distance method
"""
function edit_distance(x, y)
    #Convert strings to char arrays so that we can index into it
    xx = [i for i in x]
    yy = [i for i in y]

    m=length(xx)
    n=length(yy)

    r = zeros(Int, m+1, n+1)

    # Iterate through substrings
    for i in 1:(m + 1)
        for j in 1:(n + 1)
            if i == 1
                r[i, j] = j
            elseif j == 1
                r[i, j] = i
            elseif xx[i-1] == yy[j-1]
                r[i, j] = r[i-1, j-1]
            else
                r[i, j] = 1 + min(r[i, j-1], r[i-1, j],  r[i-1, j-1])
            end
        end
    end
    return r[m+1, n+1]
end

function validate_name(name::String)
    if occursin(' ', name)
        @warn("Do not use spaces in resource names. It may cause problems when moving accross platforms: $name")
    end
    if lowercase(name) != name
        @warn("Use lowercases names for resource files. It is safer when moving between windows and unix: $name")
    end
end


text:
using RAFF
using DelimitedFiles
using PyPlot

#
# This example shows how solve a problem given by data in a
# file. Also, we show how to use the `model_list` utility structure to
# retrieve pre-defined models in RAFF.
#

# Retrieve the number of parameters, the model function and the model
# function string (not used) for the "Cubic" model.
n, model, modelstr = RAFF.model_list["cubic"]

# Read data from file
open("C1.txt") do fp

    global data = readdlm(fp)

end

# Set starting point
initguess = [0.0, 0.0, 0.0, 0.0]

# Set the number of multistart iterations.
maxms = 10

# Call RAFF. In this case, the model is not a multivariate one. The
# last column of the file indicates which observation is the outlier.
rsol = raff(model, data[:, 1:end - 1], n; MAXMS=maxms, initguess=initguess)

println(rsol)

# Now we plot the solution of the problem

x = data[:, 1]
y = data[:, 2]
c = data[:, 3]

# Plot the obtained solution
modl1 = (x) -> model(x, rsol.solution)

t = minimum(x):0.01:maximum(x)
PyPlot.plot(t, modl1.(t), "r", label="RAFF")

# Plot the 'true' solutios, i. e., the one use to generate the
# problem.
θSol = [2.0, 0, -4.0, -10]
modl2 = (x) -> model(x, θSol)

PyPlot.plot(t, modl2.(t), "b--", label="True Solution")

PyPlot.legend(loc=4)

# Color the obtained outliers
c[rsol.outliers] .= 5.0

# Plot all the points
PyPlot.scatter(x, y, c=c, marker="o", s=50.0, linewidths=0.2,
               cmap=PyPlot.cm."Paired", alpha=0.6)

PyPlot.show()