text:
"""
    six_dof_aero_euler_fixed_mass(x, mass, inertia, forces, moments, h=[0.0, 0.0, 0.0])

Six degrees of freedom dynamic system using Euler angles for attitude representation
and assuming fixed mass.

Flat Earth hypothesis is applied and Earth reference frame is considered
inertial.

It is considered that the aircraft xb-zb plane is a plane of symmetry so that
Jxy and Jyz cross-product of inertia are zero and will not be taken into
account.

The effects of the angular momentum produced by spinning rotors is taken into
account with the optional argument `h`.

# Arguments
- `x::12-element Array{Number,1}`: state vector.
    tas: inertial speed along x_wind axis (m/s)
    α: angle of attack (rad)
    β: angle of side-slip (rad)
    p, q, r: inertial rotatinal velocity expressed in body axis. (rad/s)
    ψ, θ, ϕ: attitude given by Euler angles (zyx). Yaw, pitch, roll. (rad)
    xe, ye, ze: position wrt the inertial system origin expressed in Earth Axis. (m)
- `mass::Number`: total mass of the aircraft (kg)
- `inertia::3×3 Array{Number,2}`: inertia tensor (kg·m²)
- `forces::3-element Array{Number,1}`: total forces expressed in body axis. (N)
- `moments::3-element Array{Number,1}`: total moments expressed in body axis.(N·m)
- `h::3-element Array{Number,1}`: Additional angular momentum contributions such
  as those coming from spinning rotors (kg·m²/s).

# Returns
- `x_dot`: state vector derivative according to the equation of motion,
    inertial properties and applied forces and moments.
"""
function six_dof_aero_euler_fixed_mass(x, mass, inertia, forces, moments, h=[0.0, 0.0, 0.0])

    m = mass
    tas, α, β, p, q, r, ψ, θ, ϕ, xe, ye, ze = x

    Fx, Fy, Fz = forces
    L, M, N = moments

    sψ, cψ = sin(ψ), cos(ψ)
    sθ, cθ = sin(θ), cos(θ)
    sϕ, cϕ = sin(ϕ), cos(ϕ)

    u, v, w = wind2body(tas, 0.0, 0.0, α, β)

    # Linear momentum equations
    u_dot = Fx / m + r * v - q * w
    v_dot = Fy / m - r * u + p * w
    w_dot = Fz / m + q * u - p * v

    tas_dot, α_dot, β_dot = tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

    # Angular momentum equations
    Ix = inertia[1, 1]
    Iy = inertia[2, 2]
    Iz = inertia[3, 3]
    Jxz = -inertia[1, 3]

    Jxz2 = Jxz*Jxz
    Γ = (Ix*Iz - Jxz2)
    temp = (Ix + Iz - Iy)

    # Engine angular momentum contribution
    hx, hy, hz = h

    rhy_qhz = (r*hy - q*hz)
    qhx_phy = (q*hx - p*hy)

    pe_dot = Iz * rhy_qhz + Jxz * qhx_phy
    qe_dot = -r*hx + p*hz
    re_dot = Jxz * rhy_qhz + Ix * qhx_phy

    # Angular momentum equations
    p_dot = L*Iz + N*Jxz - q*r*(Iz*Iz - Iz*Iy + Jxz2) + p*q * Jxz * temp + pe_dot
    p_dot /= Γ
    q_dot = M + (Iz - Ix) * p*r - Jxz * (p*p - r*r) + qe_dot
    q_dot /= Iy
    r_dot = L*Jxz + N*Ix + p*q * (Ix*Ix - Ix*Iy + Jxz2) - q*r * Jxz * temp + re_dot
    r_dot /= Γ

    # Angular Kinematic equations
    ψ_dot = (q * sϕ + r * cϕ) / cθ
    θ_dot = q * cϕ - r * sϕ
    # ϕ_dot = p + (q * sϕ + r * cϕ) * tan(θ)
    ϕ_dot = p + ψ_dot * sθ

    # Linear kinematic equations
    xe_dot =  cθ*cψ * u + (sϕ*sθ*cψ - cϕ*sψ) * v + (cϕ*sθ*cψ + sϕ*sψ) * w
    ye_dot =  cθ*sψ * u + (sϕ*sθ*sψ + cϕ*cψ) * v + (cϕ*sθ*sψ - sϕ*cψ) * w
    ze_dot = -sθ    * u +  sϕ*cθ             * v +  cϕ*cθ             * w

    [tas_dot, α_dot, β_dot, p_dot, q_dot, r_dot, ψ_dot, θ_dot, ϕ_dot, xe_dot, ye_dot, ze_dot]
end


text:
"""
    KMATools

Package for parsing various files produced by KMA. Tested on KMA 1.3.22.
"""
module KMATools

using BioSymbols: DNA

imap(f) = x -> Iterators.map(f, x)
ifilter(f) = x -> Iterators.filter(f, x)

const SPA_HEADER = join(
    [
        "#Template",
        "Num",
        "Score",
        "Expected",
        "Template_length",
        "Query_Coverage",
        "Template_Coverage",
        "Depth",
        "tot_query_Coverage",
        "tot_template_Coverage",
        "tot_depth",
        "q_value",
        "p_value",
    ],
    '\t'
)

"""
    parse_spa(io::IO, path::String) -> Vector{NamedTuple}

Parse .spa file, returning a vector of `NamedTuple` with the following names:
`template, num, score, expected, tlen, qcov, tcov, depth, total_qcov, total_tcov,`
`total_depth, qval, pval`.

Coverages and identities are represented as fractions (`Float64`) in [0.0, 1.0]
"""
function parse_spa(io::IO, path::String)
    counted_lines = enumerate(eachline(io))
    header = let
        it = iterate(counted_lines)
        header = if it === nothing
            error("Missing header from file \"$path\"")
        else
            last(first(it))
        end
        if strip(header) != SPA_HEADER
            error("Malformed header in file \"$path\"")
        end
        header
    end
    fields = fill(SubString("", 1:0), 13)
    counted_lines |>
        ifilter(i -> !isempty(strip(last(i)))) |>
        imap() do (line_number, line)
            strip_split!(fields, line, UInt8('\t'))
            return (;
                template    = String(fields[1]),
                num         = parse(UInt, fields[2], base=10),
                score       = parse(UInt, fields[3], base=10),
                expected    = parse(UInt, fields[4], base=10),
                tlen        = parse(UInt, fields[5], base=10),
                qcov        = round(parse(Float64, fields[6]) / 100, digits=6),
                tcov        = round(parse(Float64, fields[7]) / 100, digits=6),
                depth       = parse(Float64, fields[8]),
                total_qcov  = round(parse(Float64, fields[9]) / 100, digits=6),
                total_tcov  = round(parse(Float64, fields[10]) / 100, digits=6),
                total_depth = parse(Float64, fields[11]),
                qval        = parse(Float64, fields[12]),
                pval        = parse(Float64, fields[13]),
            )
        end |>
    collect
end

const RES_HEADER = join(
    [
        "#Template",
        "Score",
        "Expected",
        "Template_length",
        "Template_Identity",
        "Template_Coverage",
        "Query_Identity",
        "Query_Coverage",
        "Depth",
        "q_value",
        "p_value",
    ],
    '\t'
)

"""
    parse_res(io::IO, path::String) -> Vector{NamedTuple}

Parse .res file, returning a vector of `NamedTuple` with the following names:
`template, score, expected, tlen, tid, tcov, qid, qcov, depth qval, pval`.

Coverages and identities are represented as fractions (`Float64`) in [0.0, 1.0]
"""
function parse_res(io::IO, path::String)
    counted_lines = enumerate(eachline(io))
    header = let
        it = iterate(counted_lines)
        header = if it === nothing
            error("Missing header from file \"$path\"")
        else
            last(first(it))
        end
        if strip(header) != RES_HEADER
            error("Malformed header in file \"$path\"")
        end
        header
    end
    fields = fill(SubString("", 1:0), 11)
    counted_lines |>
        ifilter(i -> !isempty(strip(last(i)))) |>
        imap() do (line_number, line)
            strip_split!(fields, line, UInt8('\t'))
            return (;
                template = String(fields[1]),
                score    = parse(UInt, fields[2], base=10),
                expected = parse(UInt, fields[3], base=10),
                tlen     = parse(UInt, fields[4], base=10),
                tid      = round(parse(Float64, fields[5]) / 100, digits=6),
                tcov     = round(parse(Float64, fields[6]) / 100, digits=6),
                qid      = round(parse(Float64, fields[7]) / 100, digits=6),
                qcov     = round(parse(Float64, fields[8]) / 100, digits=6),
                depth    = parse(Float64, fields[9]),
                qval     = parse(Float64, fields[10]),
                pval     = parse(Float64, fields[11]),
            )
        end |>
    collect
end

function strip_split!(v::Vector{SubString{String}}, s::Union{String, SubString{String}}, sep::UInt8)
    n = 0
    start = 1
    @inbounds for i in 1:ncodeunits(s)
        if codeunit(s, i) == sep
            n += 1
            n >= length(v) && throw(BoundsError(v, n+1))
            substr = SubString(s, start, i-1)
            v[n] = strip(substr)
            start = i + 1
        end
    end
    n + 1 != length(v) && error("Incorrect number of fields for strip_split!")
    @inbounds v[n+1] = strip(SubString(s, start, ncodeunits(s)))
    v
end

"""
    parse_map(io::IO, path::String) -> Vector{Tuple{String, Vector{Row}}}

Parse .mat file, returning a vector of sequence matrices. A sequence matrix
is a `Tuple{String, Vector{Row}}` with the sequence name as the first part. The `Row`
is `Tuple{DNA, NTuple{6, UInt32}}`, with the 6 depths being the number of
A, C, G, T, N and gap, respectively.
"""
function parse_mat(io::IO, path::String)
    rowT = Tuple{DNA, NTuple{6, UInt32}}
    result = Vector{Tuple{String, Vector{rowT}}}()
    current = Vector{rowT}()
    fields = Vector{SubString{String}}(undef, 7)
    linedepths = Vector{UInt32}(undef, 6)
    header = nothing
    for line in eachline(io) |> imap(strip) |> ifilter(!isempty)
        if startswith(line, '#')
            if !isempty(current)
                @assert header isa String
                push!(result, (header, copy(current)))
                empty!(current)
            end
            header = String(line[2:end])
            continue
        end
        isnothing(header) && error("Expected header in file \"$path\"")
        strip_split!(fields, line, UInt8('\t'))
        if ncodeunits(first(fields)) != 1
            error("Multi-character reference nucleotide in file \"$path\"")
        end
        refnuc = DNA(first(first(line)))
        depth_tuple = ntuple(i -> parse(UInt32, @inbounds(fields[i+1]), base=10), Val(6))
        @inbounds for i in 1:6
            linedepths[i] = parse(UInt32, fields[i+1], base=10)
        end
        push!(current, (refnuc, depth_tuple))
    end
    isempty(current) || push!(result, (header, current))
    return result
end

export parse_spa, parse_res, parse_mat

end # module

ext: jl


text:
function saveplot(res)
  header = res.header
  df = res.data

  res = []
  time = []

  for col in names(df)
      df[col] = Missings.coalesce.(df[col], 0)
  end

  for row in eachrow(df[3:end])
      push!(res, rootmeansquare(Array(row)))
  end

  for t in df[:2]
    u = Dates.format(Dates.unix2datetime(t), "yyyy-mm-dd HH:MM:SS")
    push!(time, u)
  end

  s = scatter(x=time, y=res, mode="lines")

  layout = Layout(title="Egg Data (Root Mean Square) for " *
                Dates.format(header.start_time, "yyyy-mm-dd HH:MM:SS") * " - " *
                Dates.format(header.end_time, "yyyy-mm-dd HH:MM:SS") *
                " (" * string(header.eggs_reporting) * " Eggs Reporting)",
                width=1200,
                height=700)

  p = plot(s, layout)

  savefig(p, "../out/p.html")
end

function saveplotmultiday(res)
  df = DataFrame()
  for r in res
      df = vcat(df, r.data, cols=:union)
  end

  header = Header()

  header.start_time = res[1].header.start_time
  header.end_time = res[end].header.end_time
  header.eggs_reporting = res[end].header.eggs_reporting

  out = Result(header, df, "")
  saveplot(out)
end


text:
"""
  Provides an include guard mechanism for julia.

  exports:

    * @includeonce
    * @reinclude
"""
module IncludeGuards

export @includeonce, @reinclude

GUARDS = Set{String}()

"""
  Usage: `@includeonce "filename"`

  Use instead of include to ensure that the file is included once and only once
  within this module.
"""
macro includeonce(filename)
	global GUARDS
	locationkey = string( "module:",
			     __module__,
			     "-include:",
			     filename)
	@debug "includeonce::"  __source__ __module__ locationkey
	if locationkey in GUARDS
		@debug "includeonce: skipping include:" filename
		return nothing
	end
	push!(GUARDS, locationkey)
	@debug "includeonce: including:" filename
	return :( Base.include($__module__, $filename) )
end

"""
  Usage: `@reinclude "filename"`

  For use during development from the REPL.  Resets all include guards so that
  modified code gets reloaded.
"""
macro reinclude(filename::String)
	clear()
	return :( Base.include($__module__, $filename) )
end

function clear()
	empty!(GUARDS)
end

end # module

max_stars_repo_name: tgflynn/IncludeGuards

max_issues_repo_name: tgflynn/IncludeGuards

max_forks_repo_name: tgflynn/IncludeGuards


text:
struct Simulation{T1, T2, T3, T4}
    imodel::T1 # input model
    smodel::T2 # synapse model
    omodel::T3 # output model
    filter::T4 # filtering algorithm
end

struct SimulationState{T1, T2}
    hstate::T1 # hidden state
    fstate::T2 # filter state
end

# initialize a state for a fresh simulation based on an input model and a filtering algo
function SimulationState(imodel::InputModel, filter::SF)
    hstate = State(imodel)
    fstate = FilterState(filter)
    return SimulationState(hstate, fstate)
end

function update!(sstate::SimulationState, sim::Simulation, dt, t)
    hstate = sstate.hstate
    fstate = sstate.fstate

    imodel = sim.imodel
    smodel = sim.smodel
    omodel = sim.omodel
    filter = sim.filter

    update!(hstate, imodel, dt, t)    # sample new input
    update!(hstate, smodel, dt, t)    # update weights
    obs = update!(hstate, omodel, dt) # generate observation
    update!(fstate, filter, obs, dt)  # update filter state

    return sstate
end

function run(sim, num_timesteps, timestep, errormeasure, burnin)
    if burnin >= num_timesteps
        throw(DomainError(burnin, "burnin must be smaller than total number of timesteps"))
    end
    sstate = SimulationState(sim.imodel, sim.filter)
    total_error = 0.
    time = 0.
    for k in 1:burnin
        update!(sstate, sim, timestep, time)
        time += timestep
    end
    for k in burnin+1:num_timesteps
        update!(sstate, sim, timestep, time)
        total_error += errormeasure(sstate)
        time += timestep
    end
    return total_error / (num_timesteps - burnin)
end


text:
"Z score a single cell data set

The data will be centered around +10 to make use of the sparse matrix data compression
Not expressed gene/cell combination will get the 0 value and genes lost in the normalization process will get a -1
"
function Zscore(data::SingleCellExpr)
  function scale(v::SparseVector)
    z = copy(v)
    z[z.nzind] = 10 + z[v.nzind] -mean(z[z.nzind]) / std(z[z.nzind])
    return(z)
  end
  function zScore(x::SparseVector)
    ##n = find(x .>= 0) ## there are no 0 vals ;-)
    if (length(find(x .== -1)) == 0 )
      x =   scale(x)
    else
      ok = find(x .>= 0)
      x[ok] = scale(x[ok])
    end
    return x
  end
  ## now I need to dig through the data in row direction (bad)
  nm = length( nonzeros(data.data) )
  pos = 0
  #define max result values
  v = zeros(nm)
  r = zeros(nm)
  c = zeros(nm)
  cells = size(data.data)[2]
  for i =1:size(data.data)[1]
    val = zScore( data.data[i,1:cells] );
    ## collect the Int64 > 0 values
    ok = (pos+1):(pos+length(val.nzind));
    v[ok] = collect(val.nzval);
    c[ok] = val.nzind;
    r[ok] = fill(i, length(val.nzind));
    pos = pos+length(val.nzind);
  end
  data.zscored = sparse( r[1:pos], c[1:pos] , v[1:pos] )

  return data
end


text:
function include_folder(cur_package::Module, cur_folder::AbstractString="."; is_sorted::Bool=false, except_for::AbstractArray=[], verbose::Bool=false, is_test::Bool=false)

  Docs.initmeta(cur_package)

  # allow loading folder by index file

  if isempty(except_for) && endswith(cur_folder,".jl")
    cur_file = cur_folder

    cur_folder = dirname(cur_file)
    except_for = [cur_file]
  end

  # load all files

  all_files = get_all_files(cur_folder, except_for=except_for, is_sorted=true)

  cur_cabinet = make_initial_load(cur_package, all_files, verbose=verbose, is_test=is_test)

  is_test || do_main_load_loop!(cur_package, cur_cabinet, is_sorted=is_sorted, verbose=verbose)

  # raise errors for undefined variables

  isempty(cur_cabinet.file_infos) ||
    start_load_failure!(cur_package, cur_cabinet, verbose=verbose)

  # print verbose messages

  if verbose
    print_load_info(cur_cabinet)
    println("\ndone.\n")
  end

  # return the file cabinet

  cur_cabinet

end


text:
module Shark

include("constants.jl")
include("utils.jl")
include("engine.jl")

end # module

max_stars_repo_name: Joshua-Robison/Shark

max_issues_repo_name: Joshua-Robison/Shark

max_forks_repo_name: Joshua-Robison/Shark


text:
struct Player
    id::String
    ships::Dict{String, Ship}

    function Player(id::String)
        new(id, Dict{String, Ship}())
    end
end

"""
    get_ship(player, ship_id::String)

Returns the ship designated by `ship_id` belonging to this user.
"""
get_ship(player::Player, ship_id::String) = player.ships[ship_id]

"""
    all_ships(player, collect = false)

    Returns either iterator over all ships which belong to the user (collect = false), or a list of all ships which belong to the user (collect = true).
"""
function all_ships(player::Player, coll::Bool=false)
    if coll
        return collect(values(player.ships))
    else
        return values(player.ships)
    end
end

max_stars_repo_name: xinnosuke/Halite-II

max_issues_repo_name: xinnosuke/Halite-II

max_forks_repo_name: xinnosuke/Halite-II


text:
include("straight_roadways.jl")
include("vehicles.jl")


text:
VERSION >= v"0.0.2" && __precompile__(true)

module MarketCycles

export
    SuperSmoother, Decycler, Decycle_OSC, BandPassFilter, DominantCycle, HurstCoefficient, HPLPRoofingFilter,
    ZeroMeanRoofingFilterK0, ZeroMeanRoofingFilterK1, RoofingFilterIndicator,
    ModifiedStochastic, ModifiedRSI, AutoCorrelationIndicator, SingleLagAutoCorrelationIndicator, AutoCorrelationPeriodogram,
    AutoCorrelationReversals, DFTS, AdaptiveRSI, AdaptiveStochastic

    include("ehlers_cycles.jl")

end


text:
using BOMoD
using Random
using StatsBase
using LinearAlgebra
#using Optim
#import Optim: minimizer
include("Binary_Quadratic.jl")
# random
rng = MersenneTwister(1)
# set up BOMoD use
# make binary moduels
mod0 = Mod(Symbol(0))
mod1 = Mod(Symbol(1))
lenght_c = 5
# no constrains model of lenght 10, not all constructs are generated
all_moduels =Group_Mod([mod0,mod1])
σ²_mean = 10^-6
# the design
design = construct_design(all_moduels,lenght_c,order = true)
# the space
a_space = BOMoD.getspace(design)

# settings of the toy datatsettings

lenght_c = 10
λ = 0
Lc = 10
# setting up the toydata
const Q = compute_decay(rng,lenght_c,Lc)

# first sampling step

points = sample(rng,a_space,10,with_index = true )
unpackconstruct(points[:,1][1])
#get y_train values
# zet up the generated data
y_train = map(x -> unpackconstruct(x) |> (x -> quad(x,Q,λ)),points[:,1])
constructs_train = points[:,1]
MyKernel = LevStehnoexp(1)
Model = gp_optimised(constructs_train,y_train,MyKernel)

# Needs improvement to prevent the explicit generation of all constructs.
all_constructs = map(x -> unpackconstruct(x),collect(space))
x_test = [con for (index,con) in enumerate(all_constructs) if index ∉ points[:,2]]
x_test_space = [con for (index,con) in enumerate(space) if index ∉ points[:,2]]
constructs_test = [construct for construct in space if sum([isequal(construct,train_con) for train_con in points[:,1]]) == 0 ]

in points[:,1])), space.space)
new_predictions = save_thompson_sampling(Model.GP_model,x_test_space,10,10^-6)

a_space[138] in points[:,1]

myfilter(x::Ordered_Construct) = sum([isequal(i,x) for i in points[:,1]]) == 0 ? x : nothing


text:
# This runs the forward model to make an ice thickness and then inverts for it.

# TODO:
# - run several tests
# - make into a test suite
# - figure out whether dropping temp or fsl is better

println("Loading BITEModel.jl")
@everywhere begin
    using Plots, StatPlots
    if strip(readstring(`hostname`))=="x230"
        pyplot()
    else
        ENV["GKSwstype"]="100" # for headless/no display
        gr()
    end
    using Optim
    using BITEModel
    const BM=BITEModel
    import JLD2,JLD
    using FileIO
#    import NLopt
    import KissMCMC
    using Base.Iterators

    plotyes = false
    verbose = false
    gid = [BM.SyntheticGlacier(:bench), BM.ITMIXGlacier(:Unteraar)][1]

    # Run it
    ########
    mcmc_niter = 2*10^5
    # parameters to play with
    @show error_fac = [0,1][1]
    seed = 5

    @show meas_dens_s = [:low, :normal, :high]
    @show fsl0s = [0.1, 0.5, 0.9]
    @show temp0s = [-10, -0.5]
    @show sigma_smalls = [5.0, 0.5]
    n_1d_theta = 3
    pm0 = BM.MPara()
    theta0_dict_def = Dict(:btilde=> zeros(BM.F,n_1d_theta), # NOTE in m/y (error)
                           :fsl=> zeros(BM.F,n_1d_theta), # error
                           :temp=> zeros(BM.F,n_1d_theta), # error
                           :dist_exp=> pm0.dist_exp,
                           :iv_h_exp=> pm0.iv_h_exp,
                           :iv_dist_exp1=> pm0.iv_dist_exp1)
    #:iv_dist_exp2=> pm.iv_dist_exp2, this does not do much!?
    vars2pop = [:na] #[:temp, :fsl, :na]

    runs=[:h, :iv, :h_iv]

    function var2filename(run, meas_dens, fsl0, temp0, sigma_small, popit)
        if error_fac==0
            out = "output/fwd-inv/$(BM.getname(gid))_run-$(run)_md-$(meas_dens)_fsl-$(fsl0)_temp-$(temp0)_sigma-$(sigma_small)_sans-$(popit)"
        else
            out = "output/fwd-inv/$(BM.getname(gid))_with-err_run-$(run)_md-$(meas_dens)_fsl-$(fsl0)_temp-$(temp0)_sigma-$(sigma_small)_sans-$(popit)"
        end
    end
end

# Use Core.println
# https://discourse.julialang.org/t/status-of-threads/3729
# Threads.@threads does not work due to plotting and probably printing

# to use up to 27 workers use product:
@sync @parallel for (run,meas_dens,fsl0) in collect(product(runs,meas_dens_s,fsl0s))
    for temp0 in temp0s, sigma_small in sigma_smalls, popit in vars2pop
        tic()
        fl = var2filename(run, meas_dens, fsl0, temp0, sigma_small, popit)
        Core.println(fl)
        theta0_dict = deepcopy(theta0_dict_def)
        pop!(theta0_dict, popit, nothing)

        # TODO: This is bad as gl,pp,pm,pn,pl get loaded again below!
        # I am not sure how the parameters flow...
        # However, I also need gl in the
        pl = BM.LoadPara(gid)
        dt, pl = BM.make_datatable(gid,pl)

        # Change defaults for :temp and/or :fsl
        dt.fsl.opts[:fsl] = fsl0
        dt.temp.opts[:temp] = temp0
        gl,pp,pm,pn,pl = load_glacier(dt,pl)

        # add the parameter to the options so the
        # forward model is used to produce the thickness:

        # reseed to make the same error for each run:
        srand(seed)

        h_opts = pl.dataset_opts[BM.ThicknessData]
        h_opts[:use_forward_model] = true
        if meas_dens==:low
            h_opts[:profile_xlocs] = [2000.0]
            h_opts[:profile_ylocs] = [0.0]
        elseif meas_dens==:high
            h_opts[:profile_xlocs] = [100:100:5900;]
            h_opts[:profile_ylocs] = [-250:250:250;]
        end
        h_opts[:fwd_para] = (gl, pp, pm, pn)
        h_opts[:err_sigma] = 10.0 * error_fac
        h_opts[:err_cor_length] = 300.0
        h_opts[:profile_margin_dist] = 30.0

        iv_opts = pl.dataset_opts[BM.IVData]
        iv_opts[:use_forward_model] = true
        iv_opts[:fwd_para] = (gl, pp, pm, pn)
        iv_opts[:err_sigma] = 1.0 * error_fac
        iv_opts[:err_cor_length] = 300.0
        if meas_dens==:high
            iv_opts[:iv_grid] = (gl.dem.x, gl.dem.y)
        elseif meas_dens==:low
            iv_opts[:iv_grid] = (gl.dem.x[1:10:end], gl.dem.y[1:10:end])
        end

        # if ITMIX swap out IV and thickness
        if gid isa BM.ITMIXGlacier
            pl.dataset_LOADERS[BM.ThicknessData] = BM.SyntheticBenchLoader
            pl.dataset_LOADERS[BM.IVData] = BM.SyntheticBenchLoader
        end

        # Now update all
        dt, pl = BM.make_datatable(gid,pl)
        # for convenience:
        h_opts_dt = dt.thickness.opts
        iv_opts_dt = dt.iv.opts

        # Again change defaults for :temp and/or :fsl
        dt.fsl.opts[:fsl] = fsl0
        dt.temp.opts[:temp] = temp0
        gl,pp,pm,pn,pl = load_glacier(dt,pl)
        gb,pm = BM.make_bands(gl, pp, pm, pn);

        # Now, finally update pm (only benign things!!!)
        ## set model uncertainty to almost 0
        pm = BM.MPara(pm, sigma_h_model=sigma_small, sigma_iv_model=sigma_small)

        #print("Time to run forward model:"); gc()
        hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb =
            BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

        qtot1d, qd1d = BM.postproc1d(gb, pp, pm)

        fig = BM.plot_run(gb, hs1d, taus1d, ivs1d, hs2d, taus2d,
                          ivs2d, taus1d_l, pp, pm,
                          reuse=false, plot2d=[:h,:iv,:ivlog][1])
        Plots.savefig(fig, fl*"--map.png")

        # Inversion
        ###########
        (theta0, logposterior, logposterior1d, logprior, logprior_debug,
         loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
         pmcmc_defaults, fit_target) =
            BM.init_inverse(gb, pp, pm, pn, n_1d_theta=n_1d_theta, run=run,
                            theta0_dict = theta0_dict     )

        th0_true = deepcopy(theta0.th0)

        if h_opts[:err_sigma]==0 && iv_opts[:err_sigma]==0
            # check that th0_true is indeed minimal
            good = true
            maxll = loglikelihood(th0_true)[1]
            for delta = [1e-8, 1e-4, 1e-3,1e-2,1,10,100,1000]
                for i=1:length(th0_true)
                    th = copy(th0_true)
                    th[i] += delta
                    if logprior(th)>-Inf && loglikelihood(th)[1]>maxll
                        good = false
                        warn("loglikelihood < for $th, delta=$delta, i=$i")
                        Core.println(loglikelihood(th)[1]-maxll)
                    end
                    th[i] -= 2delta
                    if logprior(th)>-Inf && loglikelihood(th)[1]>maxll
                        good = false
                        warn("loglikelihood < for $th, delta=-$delta, i=$i")
                        Core.println(loglikelihood(th)[1]-maxll)
                    end
                end
            end
        end

        # some tests:
        fwdm_fn(theta0.th0);
        # print("Time to run forward model:")
        hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l = fwdm_fn(theta0.th0);
        # print("Time to calculate posterior:")
        logposterior(theta0.th0)[1]

        #         ###########
        #         # using NLopt
        #         #
        #         # Much faster than Optim.jl. LN_SBPLX seems better than LN_NELDERMEAD

        #         ic = th0_true + rand(length(th0_true)) * 0.01

        #         od_ll = (th, grad) -> (logprior(th)[1]>-Inf) ? loglikelihood(th)[1] : -Inf
        #         od_post = (th, grad) -> logposterior(th)[1]

        #         # https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#local-derivative-free-optimization
        #         alg = [:LN_COBYLA, :LN_NELDERMEAD, :LN_BOBYQA, :LN_PRAXIS, :LN_SBPLX][end]
        #         opt = NLopt.Opt(alg, length(ic))

        #         NLopt.max_objective!(opt, od_post)
        #         tol = 1e-6
        #         NLopt.ftol_rel!(opt, tol)
        #         NLopt.ftol_abs!(opt, tol)
        # NLopt.xtol_abs!(opt, tol)
        # NLopt.xtol_rel!(opt, tol)
        # NLopt.maxeval!(opt, 5000)

        # maxf, theta_start, ret = NLopt.optimize(opt, ic)
        # good_returns = [:SUCCESS, :STOPVAL_REACHED, :FTOL_REACHED, :XTOL_REACHED]
        # if !(ret in good_returns)
        #     Core.println("Optimization not successful: $ret\n   Best max found=$maxf vs true $(logposterior(th0_true)[1])")
        #     # else
        #     #     println("Optimization successful.\n     Max found=$maxf vs true $(logposterior(th0_true)[1])")
        # end
        # theta0.th0[:] = theta_start

        theta0.th0[:] = th0_true + rand(length(th0_true)) * 0.01

        #############
        # MCMC
        ############
        pmcmc = BM.MCMCNum(niter = mcmc_niter,
                           nthin = 20,
                           sigma_ppdf = ones(length(theta0.th0))*3e-2,
                           alg=[:emcee, :metro][1]) # metro needs many more iterations:
        # - emcee good with 10^5
        # - metro needs 10^?

        # sample posterior
        thetas, theta_expect, accept_ratio, blobs, mc_hs1d, mc_ivs1d,
        mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d = mcmc(logposterior,
                                                                              theta0, pmcmc;
                                                                              verbose=false,
                                                                              use_progress_meter=false)

        # Figure out errors of h and iv versus the truth
        # 1D
        true_h1d = h_opts_dt[:fwd_true_h1d]
        true_iv1d = iv_opts_dt[:fwd_true_iv1d]
        true_h_th, true_iv_th = BM.h_iv_at_thetas(theta0, true_h1d, true_iv1d)
        h_th, iv_th = BM.h_iv_at_thetas(theta0, mc_hs1d, mc_ivs1d)

        # do some scatter plots of these:
        err_h_th = h_th .- true_h_th
        err_iv_th = iv_th .- true_iv_th

        rel_err_h_th = (h_th .- true_h_th)./true_h_th
        rel_err_iv_th = (iv_th .- true_iv_th)./true_h_th

        # # 2D
        # hfn = h_opts_dt[:fwd_thick_fn]
        # ivfn = iv_opts_dt[:fwd_iv_fn]
        # true_h2d = Gridded(gl.dem.x, gl.dem.y, hfn.(gl.dem.x, gl.dem.y'));
        # true_iv2d = Gridded(gl.dem.x, gl.dem.y, ivfn.(gl.dem.x, gl.dem.y'));

        # # sample the prior to get that distribution too:
        # pmcmc_prior = BM.MCMCNum(niter=5*10^4,
        #                          nthin = 10)
        # thetas_prior, theta_expect_prior, accept_ratio_prior = mcmc(logprior, theta0, pmcmc_prior; verbose=true)

        # pmcmc_prior2 = BM.MCMCNum(niter=10^6,
        #                           nthin = 1, alg=:metro,
        #                           sigma_ppdf=ones(length(theta0.th0))/5  )
        # thetas_prior2, theta_expect_prior, accept_ratio_prior = mcmc(logprior, theta0, pmcmc_prior2; verbose=true)

        varnames = BM.get_varnames(theta0)
        #print_results(thetas, accept_ratio, names=varnames)

        fig = BM.plottheta(thetas, theta0, toplot=:all, size=(1200,1200))
        Plots.savefig(fig, fl*"--thetas.png")

        # plot of h-errs against other variables
        tmp = copy(thetas)
        tmp[theta0.names[:fsl],:] = err_h_th
        va = strip.(varnames)
        va[theta0.names[:fsl]] = ["h_top","h_mid","h_low"]
        fig = StatPlots.cornerplot(tmp'; size=(1200,1200), label=va)
        Plots.savefig(fig, fl*"--thetas-h-temp.png")

        tmp = copy(thetas)
        tmp[theta0.names[:temp],:] = err_h_th
        va = strip.(varnames)
        va[theta0.names[:temp]] = ["h_top","h_mid","h_low"]
        fig = StatPlots.cornerplot(tmp'; size=(1200,1200), label=va)
        Plots.savefig(fig, fl*"--thetas-h-fsl.png")

        tmp = copy(thetas)
        tmp[theta0.names[:fsl],:] = err_iv_th
        va = strip.(varnames)
        va[theta0.names[:fsl]] = ["iv_top","iv_mid","iv_low"]
        fig = StatPlots.cornerplot(tmp'; size=(1200,1200), label=va)
        Plots.savefig(fig, fl*"--thetas-iv-temp.png")

        tmp = copy(thetas)
        tmp[theta0.names[:temp],:] = err_iv_th
        va = strip.(varnames)
        va[theta0.names[:temp]] = ["iv_top","iv_mid","iv_low"]
        fig = StatPlots.cornerplot(tmp'; size=(1200,1200), label=va)
        Plots.savefig(fig, fl*"--thetas-iv-fsl.png")

        # summaries
        # as JLD
        res = KissMCMC.summarize_run(thetas, theta_true=th0_true, names=varnames)
        save(fl*".jld2", "res", res)

        # as text
        dir_,fll = splitdir(fl)
        ar = hcat([[strip(string(v)),e] for (v,e) in zip(res[:var],res[:err])]...)
        l1,l2 = split(reprmime("text/plain", ar),'\n')[2:3]
        write(fl*"--summary", "$fll\n    $l1\n    $l2\n")
        toc()
    end
end


text:
include("../examples/lorenz.jl")
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function sens(s,nSteps)
    d = 3
    n_runup = 2000
    u = rand(d)
    for i = 1:n_runup
    u = next(u, s)
    end
    u_trj = step(u, s, nSteps)
    x, y, z = view(u_trj,1,:), view(u_trj,2,:), view(u_trj,3,:)
    du_trj = dstep(u_trj, s)
    ddu_trj = d2step(u_trj,s)
    f_trj = flow(u_trj,s)

    vs = zeros(d)
    q = rand(d)
    q /= norm(q)

    pp = pert(u_trj, s)
    dpp = dpert(u_trj, s)
    q1 = rand(d)
    vs1 = zeros(d)

    dJds_st = 0.
    dJds_ust = 0.
    dJds_c = 0.
    N = 2000
    nSteps = nSteps + 1
    g = zeros(nSteps)
    J = (z .- 28.0).^2.0
    #dJdu = [0., 0., 1.]
    # Df represents the derivative of f
    # on the
    # unstable manifold.
    a = zeros(nSteps)
    c = zeros(nSteps)
    beta = zeros(nSteps)
    dJf = zeros(nSteps)
    Da = zeros(nSteps)
    Dc = zeros(nSteps)
    Dq = zeros(d)
    Dq1 = zeros(d)
    Dvs = zeros(d)
    Dvs1 = zeros(d)
    d2varphi = zeros(d,d,d)
    d2q = zeros(d,d)
    rho = s[2]

    # nSteps large number.
    for i = 1:nSteps-1
        dui = du_trj[:,:,i] # for large systems, can't store Jacobian.
        ppi = pp[:,i]
        dppi = dpp[:,:,i]
    	beta[i+1] = norm(f_trj[:,i+1])
    	f1 = f_trj[:,i+1]/beta[i+1]

        q1 .= dui*q
        alpha = norm(q1)
        q1 ./= alpha
        alpha2 = alpha*alpha
    	th = dot(q1, f1)
        th2 = 1.0 - th*th
        r1 = dunit_flow([x[i+1],y[i+1],z[i+1]],s)*q1

        vs1 .= dui*vs + ppi
        a[i+1] = dot(vs1, q1 - th*f1)/th2
    	c[i+1] = dot(vs1, f1 - th*q1)/th2

        d2varphi .= permutedims(ddu_trj[:,:,:,i],
				 [1,3,2])
    	for j = 1:d
    		d2q[:,j] .= d2varphi[:,:,j]*q
    	end

        Dq1 .= d2q*q/alpha2 + dui*Dq/alpha2
        dalphadx = dot(Dq1, q1)*alpha
	Dq1 .= Dq1 .- dalphadx*q1/alpha

        Dvs1 .= d2q*vs/alpha + dui*Dvs/alpha + dppi*q1
    	Dth = dot(q1, r1) + dot(f1, Dq1)
    	constant = 2*th*Dth/th2/th2
    	Da[i+1] = constant*dot(vs1, q1 - th*f1)
    	Da[i+1] += 1/th2*(dot(vs1, Dq1 - Dth*f1 - th*r1) +
    				 dot(Dvs1, q1 - th*f1))

    	Dc1 = constant*dot(vs1, f1 - th*q1)
    	Dc1 += 1/th2*(dot(vs1, r1 - Dth*q1 - th*Dq1) +
    				 dot(Dvs1, f1 - th*q1))

        g[i+1] = g[i]/alpha - dalphadx/alpha
        vs1 .= vs1 .- c[i+1]*f1 .- a[i+1]*q1

    	Dvs1 .= Dvs1 .- a[i+1]*Dq1 .- Da[i+1]*q1 .- c[i+1]*r1 .- Dc1*f1

	dJdu = [0.,0., 2*(z[i] - 28.0)]
        dJds_st += dot(dJdu, vs)/nSteps
    	dJf[i+1] = dot(dJdu, f1)

        vs .= vs1
        q .= q1
        Dvs .= Dvs1
	Dq .= Dq1
    end

    for n = 1:N
        J_shift = J[n+1:end]
    	dJf_shift = dJf[n+1:end]
    	betan_shift = beta[n+1:end]

        nJ = length(J_shift)
        g_shift = g[2:nJ+1]
        Da_shift = Da[2:nJ+1]
        a_shift = a[2:nJ+1]
    	c_shift = c[2:nJ+1]
    	beta_shift = beta[2:nJ+1]

        dJds_ust -= dot(J_shift,Da_shift)/nJ
        dJds_ust -= dot(J_shift,g_shift.*a_shift)/nJ
    	dJds_c += dot(c_shift./beta_shift, dJf_shift.*betan_shift)/nJ
    end
    return dJds_st + dJds_ust + dJds_c

end
function get_sens(s)
    nSteps = 200000
    n_exps = size(s)[2]
    dJds = zeros(n_exps)
    var_dJds = zeros(n_exps)
    n_rep = 16
    dJds_proc = SharedArray{Float64}(n_rep)
    for k=1:n_exps
        sk = s[:,k]
        dJds_proc .= 0.
        t = @distributed for j=1:n_rep
            dJds_proc[j] = sens(sk,nSteps)/n_rep
        end
        wait(t)
        dJds[k] = sum(dJds_proc)
        var_dJds[k] = n_rep*sum(dJds_proc.*dJds_proc) - (dJds[k])^2
    end
    save("../data/sens/lorenz/dvarz_ds2_2000.jld", "rho",
         s[2,:], "dJds", dJds, "var_dJds", var_dJds)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using Clang.cindex
using Clang.wrap_c
using Printf

# Set these to correspond to your local filesystem's curl and clang include paths
const CURL_PATH = "/usr/local/opt/curl/include/curl"
const CLANG_INCLUDES = [
    "/usr/local/opt/llvm/include",
]

const SRC_DIR = abspath(@__DIR__, "..", "src")

headers = map(x -> joinpath(CURL_PATH, x), ["curl.h", "easy.h", "multi.h"])
context = wrap_c.init(;
    headers = headers,
    clang_args = String[],
    common_file = joinpath(SRC_DIR, "lC_common_h.jl"),
    clang_includes = CLANG_INCLUDES,
    clang_diagnostics = true,
    header_wrapped = (top_header, cursor_header) -> in(cursor_header, headers),
    header_library = header -> "libcurl",
    header_outputfile = header -> joinpath(SRC_DIR, "lC_") * replace(basename(header), "." => "_") * ".jl",
)
context.options.wrap_structs = true

begin
    context.headers = [joinpath(CURL_PATH, "curl.h")]
    run(context)
end

function write_constants(filename::AbstractString, startswith_identifier::AbstractString, exports_file)
    open(filename, "r") do file
        lines = split(read(file, String), "\n")

        for line in lines
            if startswith(line, startswith_identifier)
                @printf exports_file "export %s\n" split(line, r" |\(")[2]
            end
        end
    end
end

# Generate export statements
open(joinpath(SRC_DIR, "lC_exports_h.jl"), "w+") do exports_file
    println(exports_file, "#   Generating exports")

    write_constants(joinpath(SRC_DIR, "lC_curl_h.jl"), "function", exports_file)
    write_constants(joinpath(SRC_DIR, "lC_common_h.jl"), "const", exports_file)

    # Generate define constants
    open(joinpath(SRC_DIR, "lC_defines_h.jl"), "w+") do defines_file
        println(defines_file, "#   Generating #define constants")

        hashdefs = split(read(`gcc -E -dD -P $(joinpath(CURL_PATH, "curl.h"))`, String), "\n")

        for line in hashdefs
            m = match(r"^\s*#define\s+CURL(\w+)\s+(.+)", line)

            if m !== nothing
                c2 = replace(m.captures[2], "(unsigned long)" => "")
                @printf defines_file "const CURL%-30s = %s\n"  m.captures[1]  c2
                @printf exports_file "export CURL%s\n"  m.captures[1]
            end
        end
    end
end


text:
using Base.Threads
using SequentialLOB
using Plots
pyplot()
using PyPlot
pygui(true)

using AdaptiveABC
using QuantEcon
using Random
using Distributions
using ProgressMeter

# Set Simulator Parameters and Seed
Random.seed!(37162)
ϕ = 0.7
θ = 0.3
σ = 0.01
parameter_true = [ϕ, θ, σ]
n = 1000
p0 = 100
parameter_names = ["phi", "theta", "sigma"]

# Create model, simulate to obtain observation vector and then calculate summary statistics
true_arma = ARMA(ϕ, θ, σ)
obs = simulation(true_arma, ts_length=n)
obs_summary_stats = get_summary_stats(obs)
n_summary_stats = size(obs_summary_stats, 1)

function summary_fn(parameters, n_summary_stats, n_replications)
    success = true
    model = ARMA(parameters...)
    sim_obs = Array{Float64,2}(undef, n, n_replications)
    summary_stats = Array{Float64,2}(undef, n_summary_stats, n_replications)
    for i in 1:n_replications
        sim_obs[:, i] = simulation(model, ts_length=n)
        summary_stats[:, i] = get_summary_stats(obs, sim_obs[:, i])
    end
    return success, summary_stats
end

prior = Prior([
    Uniform(0.05, 0.95),
    Uniform(0.05, 0.95),
    Uniform(0.005, 0.1),
])

weighted_mad = WeightedEuclidean(obs, get_summary_stats, ones(n_summary_stats), "MAD")
sim_sum_stats = zeros(n_summary_stats, 10000, 1)

prog = Progress(10000, 1)

seeds = Int.(rand(MersenneTwister(651411613), UInt32, 10000))
for i in 1:10000
    params = rand(seeds[i], prior)
    success, sum_stats = summary_fn(params, n_summary_stats, 1)
    sim_sum_stats[:,i,:] = sum_stats
    next!(prog)
end

inv_weights = [AdaptiveABC.MAD(sim_sum_stats[i, :, :]) for i in 1:n_summary_stats]  # For each statistic, calculate the MAD across all n_sim, n_rep values
weights = 1.0 ./ inv_weights

weighted_mad.weights = weights

function f(x, weight, reps)
    success, sum_stats = summary_fn(x, n_summary_stats, reps)
    return weight(sum_stats)
end

f_mad(x, grad) = f(x, weighted_mad, 1000)

# phi (x) vs theta (y)

num_points = 20
x_s = repeat(collect(range(0.6, stop=0.8, length=num_points)), inner=num_points)
y_s = repeat(collect(range(0.2, stop=0.4, length=num_points)), outer=num_points)

z_s = zeros(num_points^2)

prog = Progress(num_points^2, 1)

Threads.@threads for i in 1:num_points^2
    z_s[i] = f_mad([x_s[i], y_s[i], 0.01], [0])
    next!(prog)
end

param_true_vals = [ϕ, θ]

z_true = 0
for i in 1:size(x_s, 1)
    if x_s[i] >= param_true_vals[1] && y_s[i] >= param_true_vals[2]
        z_true = z_s[i]
        break
    end
end
z_range = maximum(z_s) - minimum(z_s)
true_plot_z = 0.01 * z_range

plt = Plots.plot(
    x_s, y_s, z_s, st=:surface,
    size=(1500, 800),
        xlabel="Phi (ϕ)" , ylabel="Theta (θ)",
        zlabel="Distance Value", camera=(-39, 15),
        xguidefontsize=14,yguidefontsize=14,zguidefontsize=14,
        xtickfontsize=8,ytickfontsize=8,ztickfontsize=8,
        dpi=200, zguidefontrotation=90, formatter=:plain, legend=:none
)

Plots.plot!(
    plt,
    [param_true_vals[1],param_true_vals[1]],
    [param_true_vals[2],param_true_vals[2]],
    [z_true - true_plot_z,z_true + true_plot_z],
    label="True Parameter Values",
    camera=(-39, 15),
    st=:line, c=:red, linewidth=2, formatter=:plain,);

Plots.savefig(plt, joinpath(@__DIR__, "phi_theta_surface.pdf"))

max_stars_repo_name: GantZA/CALOBMTMD

max_issues_repo_name: GantZA/CALOBMTMD

max_forks_repo_name: GantZA/CALOBMTMD


text:
function detect_cross_pnt(arr, thr, way=:up, gap=1)
    _idx_repo = []
    try
        idx = findall(arr .> thr)
        idx_diff = findall(diff(idx) .> 1)
        _idx_repo = [idx[1]; idx[idx_diff]; idx[idx_diff .+ 1]; idx[end]]
    catch BoundError
        return nothing
    end

    if way == :up
        _check = (x, i) -> x[i-1] < thr < x[i] #< x[i+1]
    elseif way == :down
        _check = (x, i) -> x[i-1] > x[i] > x[i+1]
    else
        return nothing
    end

    sort!(_idx_repo)
    _result = Vector{Int}()
    _previous = -9999

    for idx in _idx_repo
            if _check(arr, idx) && (idx - _previous > gap)
                _result = [_result; idx]
                _previous = idx
            end
    end

    _result
end


text:
# * Auxilliary type definitions for restricted bases

const RestrictionMatrix = BandedMatrix{<:Int, <:FillArrays.Ones}
const RestrictedQuasiArray{T,N,B<:Basis} = SubQuasiArray{T,N,B}
const AdjointRestrictedQuasiArray{T,N,B<:Basis} = QuasiAdjoint{T,<:RestrictedQuasiArray{T,N,B}}

const BasisOrRestricted{B<:Basis} = Union{B,<:RestrictedQuasiArray{<:Any,<:Any,<:B}}
const AdjointBasisOrRestricted{B<:Basis} = Union{<:QuasiAdjoint{<:Any,B},<:AdjointRestrictedQuasiArray{<:Any,<:Any,<:B}}

indices(B::Basis) = axes(B)
indices(B::RestrictedQuasiArray) = B.indices
indices(B, i) = indices(B)[i]

unrestricted_basis(R::AbstractQuasiMatrix) = R
unrestricted_basis(R::RestrictedQuasiArray) = parent(R)

==(A::BasisOrRestricted, B::BasisOrRestricted) =
    unrestricted_basis(A) == unrestricted_basis(B) && indices(A,2) == indices(B,2)

Base.hash(B::RestrictedQuasiArray, h::UInt) = hash(indices(B,2), hash(parent(B), h))

restriction_extents(::Basis) = 0,0
function restriction_extents(B̃::RestrictedQuasiArray)
    B = parent(B̃)
    a,b = B̃.indices[2][[1,end]]
    a-1,size(B,2)-b
end

restriction(B) = Diagonal(Ones{Int}(size(B,2)))
restriction(B̃::RestrictedQuasiArray) = last(LazyArrays.arguments(LazyArrays.ApplyLayout{typeof(*)}(), B̃))

function combined_restriction_selection(A,B)
    parent(A) == parent(B) ||
        throw(ArgumentError("Cannot multiply functions on different grids"))

    la,lb = restriction_extents(A)
    ra,rb = restriction_extents(B)

    lsel = 1+ra:(size(A,2)-max(0,rb-lb))
    rsel = 1+la:(size(B,2)-max(0,lb-rb))
    lsel, rsel
end

function combined_restriction(A,B)
    lsel,rsel = combined_restriction_selection(A,B)
    l,u,r = if !isempty(lsel)
        l = lsel[1]-rsel[1]
        l, -l, 1
    else
        -1, -1, 0
    end
    BandedMatrices._BandedMatrix(Ones{Int}(r,size(B,2)), axes(A,2), l,u)
end

function show(io::IO, B̃::RestrictedQuasiArray{<:Any,2})
    B = parent(B̃)
    a,b = B̃.indices[2][[1,end]]
    N = size(B,2)
    show(io, B)
    write(io, ", restricted to basis functions $(a)..$(b) $(a>1 || b<N ? "⊂" : "⊆") 1..$(N)")
end

IntervalSets.leftendpoint(B::RestrictedQuasiArray) =
    leftendpoint(parent(B))
IntervalSets.rightendpoint(B::RestrictedQuasiArray) =
    rightendpoint(parent(B))

# TODO: This is invalid for non-orthogonal bases such as B-splines,
# since it selects as many quadrature nodes as there are basis
# functions in the restriction. Each B-spline is actually associated
# with a range of quadrature nodes.
locs(B::RestrictedQuasiArray) = locs(parent(B))[indices(B,2)]
real_locs(B::RestrictedQuasiArray) = real_locs(parent(B))[indices(B,2)]

centers(B::RestrictedQuasiArray) = centers(parent(B))[indices(B,2)]

weights(B::RestrictedQuasiArray) = weights(parent(B))[indices(B,2)]
inverse_weights(B::RestrictedQuasiArray) = inverse_weights(parent(B))[indices(B,2)]

distribution(B::RestrictedQuasiArray) = distribution(parent(B))

vandermonde(B::RestrictedQuasiArray) = vandermonde(parent(B))[:,indices(B,2)]


text:
function get_parameter_driver(obj::ForceModel, arg0::JString)
    return jcall(obj, "getParameterDriver", ParameterDriver, (JString,), arg0)
end

function is_supported(obj::ForceModel, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end


text:
import AdversarialPrediction: solve_p_given_abk, marginal_projection, marginal_projection_with_constraint,
    prox_max_sumlargest, prox_max_sumlargest_with_constraint
using LinearAlgebra

tol = 1e-3

@testset "test solve_p_given_rho" begin
    n = 30
    a = rand(n)
    rho = 3 * rand()
    k = rand(1:n)
    b = ones(n) * rho / k

    p = solve_p_given_abk(a, b, k)

    @test all(p .>= 0)
    @test all(p .<= (sum(p) / k + tol))

end

@testset "test projection" begin
    n = 30
    A = (rand(n,n) .- 0.4) ./ 0.01

    P = marginal_projection(A)
    ks = 1:n

    @test all(P .>= 0)
    @test all(sum(P, dims = 2) .<= 1 + tol)
    @test all(P .<= sum(P, dims = 1) ./ ks' .+ tol)
    @test sum(P ./ ks') <= 1 + tol
end

@testset "test projection2" begin

    A = [0.3 0.2; 0.4 0.21]
    X = [0.3 0.205; 0.4 0.205]

    P = marginal_projection(A)

    @test P ≈ X
end


text:
# types.jl

## Exceptions
"""
Exception for ununifiable terms in subst
"""
struct ICMPIn
 left
 right
 subst
 op
end

"""
Exception for ununifiable terms
"""
struct ICMP
 left
 right
 op
end

struct Loop
 left
 right
 op
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
####################################
#
# computes row ind of inverse of (1/2) ⋅ Q (Σy ⊗ Σx + Σx ⊗ Σy) Q'
#
####################################
struct CDInverseReducedSymKroneckerLoss{T<:AbstractFloat, S} <: CoordinateDifferentiableFunction
  Σx::Symmetric{T, S}
  Σy::Symmetric{T, S}
  A::Matrix{T}    # stores Σx⋅Θ⋅Σy
  ind::Int64
  p::Int64
end

function CDInverseReducedSymKroneckerLoss(Σx::Symmetric{T,S}, Σy::Symmetric{T,S}, ind::Int) where {T<:AbstractFloat} where S
  (issymmetric(Σx) && issymmetric(Σy)) || throw(DimensionMismatch())
  (p = size(Σx, 1)) == size(Σy, 1) || throw(DimensionMismatch())
  CDInverseReducedSymKroneckerLoss{T,S}(Σx, Σy, zeros(T, p, p), ind, p)
end

CoordinateDescent.numCoordinates(f::CDInverseReducedSymKroneckerLoss) = div((f.p + 1) * f.p, 2)

function _A_mul_symX_mul_B_rc(
  A::Symmetric{T},
  X::SparseIterate{T},
  B::Symmetric{T},
  r::Int,
  c::Int
  ) where {T<:AbstractFloat}

  _c = sqrt(2.)
  p = size(A, 1)
  v = zero(T)
  for j=1:nnz(X)
    ind = X.nzval2ind[j]
    ri, ci = ind2subLowerTriangular(p, ind)
    if ri == ci
      @inbounds v += A[ri, r] * B[ci, c] * X.nzval[j]
    else
      @inbounds v += (A[ri, r] * B[ci, c] + A[ci, r] * B[ri, c]) * X.nzval[j] / _c
    end
  end
  v
end

function CoordinateDescent.initialize!(f::CDInverseReducedSymKroneckerLoss, x::SparseIterate)
# compute residuals for the loss

  Σx = f.Σx
  Σy = f.Σy
  A = f.A
  p = f.p

  for ac=1:p, ar=1:p
      @inbounds A[ar,ac] = _A_mul_symX_mul_B_rc(Σx, x, Σy, ar, ac)
  end

  nothing
end

function CoordinateDescent.gradient(
  f::CDInverseReducedSymKroneckerLoss{T},
  x::SparseIterate{T},
  j::Int64) where {T <: AbstractFloat}

  A = f.A
  B = f.B
  ri, ci = ind2subLowerTriangular(f.p, j)
  @inbounds v = (A[ri,ci] + A[ci,ri]) / 2.
  v = ri == ci ? v : sqrt(2.) * v
  return j == f.ind ? v - 1. : v
end

function CoordinateDescent.descendCoordinate!(
  f::CDInverseReducedSymKroneckerLoss{T},
  g::ProxL1{T},
  x::SparseIterate{T},
  j::Int64) where {T <: AbstractFloat}

  Σx = f.Σx
  Σy = f.Σy
  A = f.A
  p = size(Σx, 1)

  ri, ci = ind2subLowerTriangular(f.p, j)

  a = zero(T)
  b = zero(T)
  if ri == ci
    @inbounds a = Σx[ri,ri] * Σy[ci,ci]
    @inbounds b = A[ri,ri]
  else
    @inbounds a = (Σx[ri,ri] * Σy[ci,ci] + Σx[ci,ci] * Σy[ri,ri]) / 2. + Σx[ri,ci] * Σy[ri,ci]
    @inbounds b = (A[ri,ci] + A[ci,ri]) / sqrt(2.)
  end
  b = j == f.ind ? b - 1. : b

  @inbounds oldVal = x[j]
  a = one(T) / a
  @inbounds x[j] -= b * a
  newVal = cdprox!(g, x, j, a)
  h = newVal - oldVal

  # update internals
  if ri == ci
    for ac=1:p, ar=1:p
      @inbounds A[ar, ac] += h * Σx[ar, ri] * Σy[ri, ac]
    end
  else
    _c = h / sqrt(2.)
    for ac=1:p, ar=1:p
      @inbounds A[ar, ac] += _c * (Σx[ar, ri] * Σy[ci, ac] + Σx[ar, ci] * Σy[ri, ac])
    end
  end
  h
end

############################################
#
#   compressed version of invHessian
#
############################################

function _mul_symX_S(θ::SparseIterate, S::Symmetric)
    p = size(S, 1)
    rp = length(θ)
    _c = 1. / sqrt(2.)

    out = zeros(p, p)

    i = 1
    @inbounds while i <= θ.nnz
      v = θ.nzval[i]
      ind = θ.nzval2ind[i]
      r, c = ind2subLowerTriangular(p, ind)
      if r == c
          for l=1:p
              @inbounds out[r, l] += S[c, l] * v
          end
      else
          for l=1:p
              t = v * _c
              @inbounds out[r, l] += t * S[c, l]
              @inbounds out[c, l] += t * S[r, l]
          end
      end
      i += 1
    end
    out
end

function _computeReducedVar!(ω, f::CDInverseReducedSymKroneckerLoss, X, Y, Θ)
    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    A = f.A
    ind = f.ind
    _c = 1. / sqrt(2.)

    ΘSx = _mul_symX_S(Θ, f.Σx)
    ΘSy = _mul_symX_S(Θ, f.Σy)

    V = X * ΘSy
    W = Y * ΘSx

    j = 0
    for col=1:p
        for row=col:p
            j += 1

            δ = ind == j ? - 1. : 0.

            fill!(q, 0.)
            fill!(r, 0.)
            t = 0.

            if row == col
                @inbounds t = A[row, row] + δ

                for k=1:nx
                    @inbounds q[k] = X[k, row] * V[k, row] + δ
                end
                for k=1:ny
                    @inbounds r[k] = Y[k, row] * W[k, row] + δ
                end
            else
                @inbounds t = (A[row, col] + A[col, row]) * _c + δ

                for k=1:nx
                    @inbounds q[k] = (X[k, row] * V[k, col] + X[k, col] * V[k, row]) * _c  + δ
                end
                for k=1:ny
                    @inbounds r[k] = (Y[k, row] * W[k, col] + Y[k, col] * W[k, row]) * _c + δ
                end
            end

            σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2
            σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2

            ω[j] = sqrt((σ1/nx + σ2/ny))
        end
    end

end

function invQSymHessian(Sx::Symmetric, Sy::Symmetric, ind, X, Y, λ, options=CDOptions())
    nx, p = size(X)
    ny = size(Y, 1)

    rp = div((p+1)*p, 2)
    f = CDInverseReducedSymKroneckerLoss(Sx, Sy, ind)
    x = SparseIterate(rp)
    x[ind] = 1. / skron(Sx, Sy, ind, ind)
    CoordinateDescent.initialize!(f, x)

    ##################0
    #
    #  first stage
    #
    ##################
    # compute initial variance
    ω = Array{eltype(Sx)}(undef, length(x))
    _computeReducedVar!(ω, f, X, Y, x)
    ω[ind] = 0.

    # compute initial estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, options)

    ##################
    #
    #  second stage
    #
    ##################

    # recompute variance
    opt1 = CDOptions(;
        maxIter=options.maxIter,
        optTol=options.optTol,
        randomize=options.randomize,
        warmStart=true,
        numSteps=options.numSteps)

    _computeReducedVar!(ω, f, X, Y, x)
    ω[ind] = 0.

    # recompute estimate
    g = ProxL1(λ, ω)
    coordinateDescent!(x, f, g, opt1)

    return x
end


text:
include("../src/api.jl")

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

using Plots
using Nemo
using Distributions

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

# basis_small = []
# basis_average = []
# basis_big = []

reductions = []
function uwuwu()

    for (i, (filename, x, y, system)) in enumerate(load_COO_if(from_dim=100, to_dim=130))

        @info "$i-th, loaded a system $filename of size $(length(system))"

        matrices = [
            from_COO(system[1][1], system[1][2], ms[3:end]..., QQ)
            for ms in system
        ]

        sz1 = system[1][1]

        # matrices
        V1 = invariant_subspace(deepcopy(matrices))

        sz2 = length(V1)

        push!(reductions, [filename, sz1, sz2])
        if i > 30
            break
        end

    end
end

uwuwu()

# "BIOMD0000000147.json", 26, 5
# "BIOMD0000000407.json", 50, 4
# "MODEL1108260014.json", 82, 24
# "BIOMD0000000559.json", 101, 14

# plot(1:length, times_W_det, label="W. det", ylabel="time, s",  xlabel="order, n", title="rand vs det, 4 nnz per row")
#plot!(sizes, times_W_rand, label="W. rand")

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


text:
using CUDA
using AgentBasedModels
using Test
using Distributions
using CUDA
using DataFrames
import MacroTools: prettify

if CUDA.has_cuda()
    testplatforms = ["cpu","gpu"]
else
    testplatforms = ["cpu"]
    println("CUDA was not found, only checking cpu.")
end

include("testAgent.jl")
# include("testSubstitution.jl")
# include("testAuxiliar.jl")
# include("testModel.jl")
# include("testRandom.jl")
# include("testCommunity.jl")
# include("testUpdates.jl")
# include("testIntegrator.jl")
# include("testEvent.jl")
# include("testSimulationSpace.jl")
# include("testBound.jl")
# include("testSave.jl")
# include("testMedium.jl")
# include("testConstructors.jl")

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using FDModBiotCyl
using Test

@testset "FDModBiotCyl.jl" begin
    # Write your tests here.
    include("./forward.jl")
end


text:
const FEATUREEXTRACTION_EXTRAKEY = :featureextraction
const REMOVEMISSINGS_DEFAULT = false
const KEEPDATA_DEFAULT = false


text:
using DifferentialEquations
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

include("Dimer_Eq_of_M.jl")
Ham2ndSHO(p, q, t) = 8 * (p[1]^2 + p[2]^2) + (q[1]^2 + q[2]^2);
ω_0 = [1.,1e4]
tspan = (0., 100);
p0 = [1.0,0]
q0 = [0,1.0]
params_1 = [1,2]
params_2 = [1,sqrt(2)]
param = sqrt(2)
# prob = HamiltonianProblem{true}(Ham2ndSHO, p0, q0, tspan);
# sol = solve(prob, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10)

# prob_1 = HamiltonianProblem{true}(Ham2ndSHO, p0, q0, params_1, tspan);
# sol_1 = solve(prob_1, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10)

# prob_2 = HamiltonianProblem{true}(Ham2ndSHO, p0, q0, tspan);
# sol_2 = solve(prob_2, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10)

function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function affect_update_iterator_q1!(integrator)
    integrator.u[5] = integrator.u[5] + 1
end

u0=[0.,0,1,0]
callback_hits_PSS_q1_forward = ContinuousCallback(condition_hits_PSS_q1, nothing,  affect_update_iterator_q1!, rootfind = true)
cb_forward = CallbackSet(callback_hits_PSS_q1_forward)
prob=ODEProblem(Dimer_Eq_of_M,u0,tspan )
sol = solve(prob, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10, callback = cb_forward, save_start = false, save_end = true, save_everystep = false)
plot(sol[2,:],sol[4,:],".")


text:
using Documenter, LyapunovExponents

# Generate src/gallery/examples/*.{md,png}.  In principle, I can
# generate those files in build/gallery/examples/ after `makedocs()`
# below.  However, resolving @ref etc. are done by `makedocs()` so
# it's better to generate Markdown files here.
cd(dirname(@__FILE__)) do
    run(`make gallery_md`)
    include("plot_gallery.jl")
end

makedocs()

if get(ENV, "REQUIRE_VERSIONS", "") == "latest"
    deploydocs(
        deps   = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
        repo   = "github.com/tkf/LyapunovExponents.jl.git",
        julia  = "0.6",
    )
end


text:
@testset "Testing PCI" begin
    # Make sure we find busses for each socket.
    #
    # This isn't a great test but will at least allow for testing of regressions on
    # other systems.
    @test length(CounterTools.findbusses()) == 2
end

@testset "Testing PMU" begin
    # Do some experimenting with a PMU unit.
    busses = CounterTools.findbusses()

    # Get the bus for Socket 0
    bus = busses[1]

    # Open up the PCI path to this counter
    device = CounterTools.SKYLAKE_IMC_REGISTERS[1][1].device
    fn = CounterTools.SKYLAKE_IMC_REGISTERS[1][1].fn

    handle = CounterTools.Handle(bus, device, fn)
    pmu = CounterTools.IMCUncorePMU(handle)

    # First things first - reset the PMU
    CounterTools.reset!(pmu)

    # Read from all registers - entries should all be zero at this point
    @test CounterTools.getunitstatus(pmu) == zero(UInt32)

    for i in 1:4
        @test iszero(CounterTools.getcontrol(pmu, i))
        @test iszero(CounterTools.getcounter(pmu, i))
    end

    # Now, we come up with an event
    # - try to program each of the counters:
    # - make sure our programming was successful
    # - make sure the counters are returning some non-zero value

    # The event we choose simply monitors DRAM read and write traffic.
    dram_reads = CounterTools.UncoreSelectRegister(; event = 0x4, umask = 0x3)
    dram_writes = CounterTools.UncoreSelectRegister(; event = 0x4, umask = 0xC)

    for i in 1:4
        CounterTools.setcontrol!(pmu, i, dram_reads)
        @test CounterTools.getcontrol(pmu, i) == CounterTools.value(dram_reads)

        # TODO: Simple test function
        sleep(1)
        pre = CounterTools.getcounter(pmu, i)
        @test !iszero(pre)
    end

    # Try resetting again
    CounterTools.reset!(pmu)
    for i in 1:4
        @test iszero(CounterTools.getcontrol(pmu, i))
        @test iszero(CounterTools.getcounter(pmu, i))
    end
end

@testset "Testing iMC Monitor" begin
    # Monitor DRAM reads and writes
    dram_reads = CounterTools.UncoreSelectRegister(; event = 0x4, umask = 0x3)
    dram_writes = CounterTools.UncoreSelectRegister(; event = 0x4, umask = 0xC)

    events = (
        dram_reads,
        dram_writes,
        dram_writes,
        dram_reads,
    )

    # Construct a iMC Monitor - attach it to socket 1 (IndexZero(0)) since that is the node
    # where we use `numactl` to bind the progras.
    socket = 1
    monitor = CounterTools.IMCMonitor(events, socket)

    # Test that creating another monitor throws an error
    @test_throws ErrorException CounterTools.IMCMonitor(events, socket)

    # Compile the reading test program.
    ntimes = 10
    array_size = 10^7
    compile_test_program("read", ntimes = ntimes, array_size = array_size)

    pre = read(monitor)
    # NOTE: This is running on numanode zero, so we'll look at the results
    # from Socket 0 for testing purposes
    run_test_program()
    post = read(monitor)

    # Reset the monitor so we don't clobber anything on subsequent runs.
    CounterTools.reset!(monitor)

    # Accumulate and diff
    socket_0 = CounterTools.aggregate(post - pre)
    @show socket_0

    # Compute the expected amount of read traffic.
    # Remeber that the iMC counters count in transactions and that each transaction
    # is 64 bytes
    expected_read_bytes =  array_size * 64 * ntimes
    expected_read_actions = div(expected_read_bytes, 64)

    # We write to the array once to initialize it, so compute the expected number
    # of read actions as well.
    expected_write_bytes = array_size * 64 * 2
    expected_write_actions = div(expected_write_bytes, 64)

    @show socket_0
    @show expected_read_actions
    @show expected_write_actions

    # The number of DRAM reads measured should strictly be greater than the minimum,
    # especially since we exceed the size of the L3 cache
    @test socket_0[1] >= expected_read_actions
    @test socket_0[1] <= 1.2 * expected_read_actions
    @test socket_0[4] >= expected_read_actions
    @test socket_0[4] <= 1.2 * expected_read_actions

    # The number of write actions should be well less than half the number of read actions.
    #
    # The lower and upper bounds are mainly derived from heuristics
    @test socket_0[2] >= 0.8 * expected_write_actions
    @test socket_0[2] <= 1.2 * expected_write_actions
    @test socket_0[3] >= 0.8 * expected_write_actions
    @test socket_0[3] <= 1.2 * expected_write_actions

    ##### Write Test
    compile_test_program("write", ntimes = ntimes, array_size = array_size)
    CounterTools.program!(monitor)

    pre = read(monitor)
    run_test_program()
    post = read(monitor)
    CounterTools.reset!(monitor)

    # Accumulate and diff
    socket_0 = CounterTools.aggregate(post - pre)

    # We write to the array once to initialize it, so compute the expected number
    # of read actions as well.
    expected_write_bytes = array_size * 64 * ntimes
    expected_write_actions = div(expected_write_bytes, 64)

    # In this case, we don't have a really good handle on the number of expected read
    # actions, so lets just set it for 1/5th the number of write actions

    @show socket_0
    @show expected_read_actions
    @show expected_write_actions

    @test socket_0[1] <= expected_write_actions / 5
    @test socket_0[4] <= expected_write_actions / 5

    @test socket_0[2] >= expected_write_actions
    @test socket_0[2] <= 1.2 * expected_write_actions
    @test socket_0[3] >= expected_write_actions
    @test socket_0[3] <= 1.2 * expected_write_actions

    #####
    ##### STREAM
    #####

    # Finally, lets run our modified version of the STREAM benchmark
    compile_test_program("copy", ntimes = ntimes, array_size = array_size)
    CounterTools.program!(monitor)

    pre = read(monitor)
    run_test_program()
    post = read(monitor)
    CounterTools.reset!(monitor)

    socket_0 = CounterTools.aggregate(post - pre)

    # Expect 5 times through the array for both read and write
    expected_read_bytes =  array_size * 64 * ntimes
    expected_read_actions = div(expected_read_bytes, 64)
    expected_write_actions = expected_read_actions

    # Reads
    @test socket_0[1] >= expected_read_actions
    @test socket_0[1] <= 1.2 * expected_read_actions
    @test socket_0[4] >= expected_read_actions
    @test socket_0[4] <= 1.2 * expected_read_actions

    # Writes
    @test socket_0[2] >= expected_write_actions
    @test socket_0[2] <= 1.2 * expected_write_actions
    @test socket_0[3] >= expected_write_actions
    @test socket_0[3] <= 1.2 * expected_write_actions
end


text:
# The license for this code is available at https://github.com/avik-pal/FastStyleTransfer.jl/blob/master/LICENSE.md

#-----------------------------VGG19-----------------------------------

mutable struct vgg19
    slice1
    slice2
    slice3
    slice4
end

function vgg19()
    vgg = VGG19().layers[1]  # get only the first set of layers in the chain
    slice1 = Chain(vgg[1:5]...)
    slice2 = Chain(vgg[6:10]...)
    slice3 = Chain(vgg[11:15]...)
    slice4 = Chain(vgg[16:20]...)
    vgg19(slice1, slice2, slice3, slice4)
end

function (layer::vgg19)(x)
    res1 = layer.slice1(x)
    res2 = layer.slice2(res1)
    res3 = layer.slice3(res2)
    res4 = layer.slice4(res3)
    (res1, res2, res3, res4)
end


text:
using Compose, Gadfly

set_default_plot_size(6inch, 3inch)

srand(123)
plot(x=rand(20), y=rand(20), color=repeat(["A","B"], inner=10),
     Guide.colorkey(title="Species", labels=["Name1","Name2"], pos=[0.85w,-0.3h])
     )


text:
function AccelerationNoise(highway; vmax::Array{Int8, 1} = Int8[3, 5], R::Float32 = 0.2f0, left_boundary::Int64 = 5)

    for v in highway[left_boundary+1:end]
        if v.tipo == 1 || v.tipo == 2

            #Acceleration
            if v.speed != -1
                  Cambiar_Vehiculo!(v, min(v.speed+1, vmax[v.tipo]), v.position, v.tipo, v.change, v.num, v.len)
            end

            #Noise
            if v.speed > 0
                if rand() <= R
                    Change_Vehicle!(v, max(v.speed-1, 0), v.position, v.tipo, v.change, v.num, v.len)
                end
            end
        end
    end
end

function DecelerationMove(C::Highway1D)

    alfa = [0.8f0, 0.75f0]
    left_boundary = Int8(5)

    i = Pos_left(C.highway)
    i_1 = ( i > 1 ? Pos_left(C.highway[1:i-1]) : 0 )

    # This while is for the next steps
    while (i > 0)

        if i_1 != 0
        ############################ change ###################################
          distance = i - i_1 - C.highway[i].len
          distance = (distance < 1 ? 0 : distance)

          #The vehicle i-1 deaccelarate with respect to vehicle i
          estimated_speed = round(Int8, floor( C.highway[i].speed * (1.0 - alfa[C.highway[i].tipo]) + 0.5 ) )
          Change_Vehicle!(C.highway[i_1], min( (estimated_speed+distance), C.highway[i_1].speed),
                              C.highway[i_1].position, C.highway[i_1].tipo, C.highway[i_1].change,
                              C.highway[i_1].num, C.highway[i_1].len)
        end

        #Vehicle i is moved

        next_pos = C.highway[i].speed + C.highway[i].position
        if C.highway[i].position != next_pos  && next_pos > 0 #&& C.carretera[pos_sig].tipo == -1

            if next_pos > C.N
                Empty_Cell!(C.highway[i])
                C.count -= 1
            end

            if next_pos <= C.N && C.highway[next_pos].tipo == -1

                if C.highway[i].position <= left_boundary && next_pos > left_boundary
                    C.count += 1
                end

                Change_Vehicle!(C.highway[next_pos], C.highway[i].speed, next_pos, C.highway[i].tipo, 0
                                    , C.highway[i].num, C.highway[i].len)

                Empty_Cell!(C.highway[i])
            elseif next_pos <= C.N && C.highway[next_pos].tipo == -2 && next_pos > 0
                while C.highway[next_pos].tipo == -2 && next_pos > 1
                   next_pos -= 1
                end
               Change_Vehicle!(C.highway[next_pos], C.highway[i].speed, next_pos, C.highway[i].tipo, 0,
                                C.highway[i].num, C.highway[i].len)
               Empty_Cell!(C.highway[i])
            end
        end
        #Reindexing the vehicle i-1 as the vehicle i

        i = i_1
        i_1 = ( i > 1 ? Pos_left(C.highway[1:i-1]) : 0 )
    end

    alfa = left_boundary = i = i_1 = 0
end

################################## Special Functions for the simulations with curves ############################

######### Sense 1 ###########

function AccelerationNoiseS1(highway, vmax::Array{Int8, 1} = Int8[3, 5], R::Float32 = 0.2f0, left_boundary::Int64 = 5)

    for v in highway
        if v.tipo == 1 || v.tipo == 2

            #Acceleration
            if v.speed != -1
              #Curves zone
                if v.position > left_boundary && ((v.position < 349) || ( v.position > 620 && v.position < 663))
                  vmax[2] = Int8(4)
                end
                Change_Vehicle!(v, min(v.speed+1, vmax[v.tipo]), v.position, v.tipo, v.change, v.num, v.len)

            end

            #Noise
            if v.speed > 0
                if rand() <= R
                    Change_Vehicle!(v, max(v.speed-1, 0), v.position, v.tipo, v.change, v.num, v.len)
                end
            end
            vmax[2] = Int8(5)
        end
    end
  vmax = 0
end

##### Sense 2 #######

function AccelerationNoiseS2(highway, vmax::Array{Int8, 1} = Int8[3, 5], R::Float32 = 0.2f0, left_boundary::Int64 = 5)

    for v in highway[left_boundary+1:end]
        if v.tipo == 1 || v.tipo == 2

            #Acceleration
            #if v.speed != -1
              #Curves zone
                if (v.position > 3518) || ( v.position < 3408 && v.position > 3265)
                    vmax[2] = Int8(4)
                end

                #Hill
                if (v.position > 2500) && v.tipo == 1
                    vmax[1] = Int8(2)
                end
                Change_Vehicle!(v, min(v.speed+1, vmax[v.tipo]), v.position, v.tipo, v.change, v.num, v.len)
            #end

            #Noise
            if v.speed > 0
                if rand() <= R
                    Change_Vehicle!(v, max(v.speed-1, 0), v.position, v.tipo, v.change, v.num, v.len)
                end
            end
        vmax = Int8[3, 5]
        end
    end
  vmax = 0
end


text:
export TransitionProperties

mutable struct TransitionProperties
    pause_for::Float64
    pause_for_cnt::Float64
    can_transition::Bool # true = transition completed.

    function TransitionProperties()
        new(0.0,
            0.0,
            false)
    end
end

function reset_pause(prop::TransitionProperties)
    prop.pause_for_cnt = 0.0
    prop.can_transition = false;
end

function inc_pause_cnt(prop::TransitionProperties, dt::Float64)
    prop.pause_for_cnt += dt;
end

function update(prop::TransitionProperties, dt::Float64)
    inc_pause_cnt(prop, dt)
    if prop.pause_for_cnt >= prop.pause_for
        prop.can_transition = true
    end
end

function ready(prop::TransitionProperties)
    prop.can_transition
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -5.483106975000851 <= q <= 6.193029789058127)

                     add_NL_constraint(m, :(swish(0.7547513043737739 + -0.48113681936364117*swish(-0.10262045764995253 + -0.12743000480166566*swish(0.5778005478095731 + 0.8852363638482865*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*swish(0.07154777431849713 + -0.5862593464659436*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.9820460002458025*swish(0.1263157800600383 + 0.1285425867505281*swish(0.5778005478095731 + 0.8852363638482865*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*swish(0.07154777431849713 + -0.5862593464659436*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) + swish(-0.6861294460251406 + 0.36167652828680286*swish(-0.10262045764995253 + -0.12743000480166566*swish(0.5778005478095731 + 0.8852363638482865*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.1298782680512347*swish(0.07154777431849713 + -0.5862593464659436*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))) + -0.5440699588803759*swish(0.1263157800600383 + 0.1285425867505281*swish(0.5778005478095731 + 0.8852363638482865*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + -0.8793054881197468*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))) + -0.6128023058992138*swish(0.07154777431849713 + -0.5862593464659436*swish(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])) + 0.9939436563912345*swish(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Algorithm 6.1
# Gaussian Elimination

function swap_rows(mat, row1, row2)
	for i in 1:size(mat)[2]
		p1 = mat[row1, i]
		p2 = mat[row2, i]
		tmp = p1
		mat[row1, i] = p2
		mat[row2, i] = tmp
	end
	return mat
end

function solve_gauss(mat)
	n, m = size(mat)
	for k in 1:n
		for i in (k + 1):n
			m = mat[i, k]/mat[k, k]
			for j in k:(n + 1)
				mat[i, j] -= m*mat[k, j]
			end
		end
	end
	println(mat)
	x = zeros(1, n)
	x[n] = mat[n, n + 1]/mat[n ,n]
	for i in (n - 1):-1:1
		sum = 0
		for j in (i + 1):n
			sum = mat[i, j]*x[j]
		end
		x[i] = (mat[i, n + 1] - sum)*mat[i, i]
	end
	return x
end

A = [[1.0 1.0 2]; [2.0 -2.0 3.0];]
println(A)
A = solve_gauss(A)
println(A)


text:
function factoredFastMarching3D(pEik::EikonalParam, mem::EikonalTempMemory)

kappaSquared 	= pEik.kappaSquared;
h 				= pEik.Mesh.h;
src 			= pEik.src;
n 				= pEik.Mesh.n.+1;
HO 				= pEik.HO;

N = prod(n);

if isempty(pEik.T1)
	pEik.T1 = zeros(Float64,tuple(n...));
end
T = pEik.T1;
if isempty(pEik.ordering)
	pEik.ordering = zeros(EikIdxType,N);
end

ordering = pEik.ordering;
if isempty(pEik.OP)
	pEik.OP = zeros(Int8,tuple(n...));
end
OP       = pEik.OP;
OP[:] 	 .= 0;

Done = mem.Done;
Done[:] .= false;

heapArr = mem.V;
heapIdx = mem.J;
frontHeap = initHeap(heapArr,heapIdx);

# time1 = 0.0;
# time2 = 0.0;
# time3 = 0.0;

monoThresh = max(h[1],h[2],h[3])^2;
T[:] .= Inf;

hinv = 1.0./h;
src_cs = loc2cs3D(src,n);

kappaSrc = sqrt(kappaSquared[src_cs]);
T[src_cs] = kappaSrc;

insertToHeap(frontHeap,0.0,src_cs);
k_order = 1;
Neighbors3D = zeros(Int64,6,3);
notBoundary3D = zeros(Bool,6);
curr_loc = zeros(Int64,3);
GradLoc = zeros(Float64,3);
A = zeros(Float64,3);
B = zeros(Float64,3);
OpLoc = zeros(Int8,3);
curr_val = zeros(Float64,1);

while frontHeap.size > 0
	# tic();
	curr = getMin(frontHeap,curr_val);
	cs2loc3D(curr_loc,curr,n);
	@inbounds while getPaddedArr3D(Done,curr_loc[1],curr_loc[2],curr_loc[3]) && frontHeap.size>0
		curr = getMin(frontHeap,curr_val);
		cs2loc3D(curr_loc,curr,n);
	end
	@inbounds if getPaddedArr3D(Done,curr_loc[1],curr_loc[2],curr_loc[3]) && frontHeap.size==0
        break;
    end
	# time1 += toq();
	ordering[k_order] = curr;
    k_order += 1;
	@inbounds setPaddedArr3D(Done,curr_loc[1],curr_loc[2],curr_loc[3],true);
	@inbounds get3DNeighbors(curr_loc,n,Neighbors3D,notBoundary3D);

	for k = 1:6
		@inbounds if notBoundary3D[k]
			@inbounds nei_loc1 = Neighbors3D[k,1];
			@inbounds nei_loc2 = Neighbors3D[k,2];
			@inbounds nei_loc3 = Neighbors3D[k,3];
			if !getPaddedArr3D(Done,nei_loc1,nei_loc2,nei_loc3)
				# tic();
				T0loc = analyticLocal3D(h,src,nei_loc1,nei_loc2,nei_loc3,GradLoc);
				m_loc = kappaSquared[nei_loc1,nei_loc2,nei_loc3];
				@inbounds T_k = getUpdatedFactoredVal_3D(HO,m_loc,T,Done,nei_loc1,
									nei_loc2,nei_loc3,h,hinv,n,src,T0loc,GradLoc,OpLoc,A,B);

				# @inbounds if T_k*T0loc < curr_val[1]
					# @inbounds if abs(T_k*T0loc-curr_val[1]) > monoThresh
						# println("Monotonicity not fulfilled by ",(T_k*T0loc-curr_val),", not corrected.");
					# end
				# end
				# time2 += toq();
				# tic();
				@inbounds if T_k < T[nei_loc1,nei_loc2,nei_loc3]
					@inbounds T[nei_loc1,nei_loc2,nei_loc3] = T_k;
					@inbounds OP[nei_loc1,nei_loc2,nei_loc3] = getOpCode(OpLoc[1],OpLoc[2],OpLoc[3]) ;
					@inbounds insertToHeap(frontHeap,T0loc*T_k,loc2cs3D(nei_loc1,nei_loc2,nei_loc3,n));
				end
				# time3 += toq();
			end
		end
	end
end
# println("time heap: ",time1, " time quad: ",time2, " time heapInsert: ",time3);
return T,ordering;
end;


text:
module Sodium

typealias size_t Csize_t
typealias uint32_t Uint32
typealias uint8_t Uint8

include("ls_exports_h.jl")
include("ls_common_h.jl")
include("ls_sodium_h.jl")
include("ls_defines_h.jl")

end


text:
module IOs

using ..Types
using HDF5
import FileIO

#import FileIO: save

export imread, imsave, readimg, saveimg, readseg, saveseg, readaff, saveaff
export issgmfile, readsgm, savesgm, readec, saveec, readecs, saveecs, save

function imread(fname::AbstractString)
    print("reading file: $(fname) ......")
    if ishdf5(fname)
        ret =  h5read(fname, "/main")
        println("done :)")
        return ret
    else
        # handled by FileIO
        return FileIO.load(fname)
    end
end

function imsave(fname::AbstractString, vol::Array, is_overwrite=true)
    print("saving file: $(fname); ......")
    # remove existing file
    if isfile(fname) && is_overwrite
        rm(fname)
    end

    if contains(fname, ".h5") || contains(fname, ".hdf5")
        h5write(fname, "/main", vol)
    else
        # handled by FileIO
        FileIO.save(fname, vol)
    end
    println("done!")
end

"""
read raw image
"""
function readimg(fimg::AbstractString)
    if ishdf5(fimg)
        f = h5open(fimg)
        if has(f, "img")
          img = read(f["img"])
        elseif has(f, "image")
          img = read(f["image"])
        else
          img = read(f["main"])
        end
        close(f)
    else
        img = reinterpret(UInt8, load(fimg).data)
        if contains(fimg, ".tif")
            # transpose the X and Y
            perm = Vector{Int64}( 1:ndims(img) )
            perm[1:2] = [2,1]
            img = permutedims(img, perm)
        end
    end
    return img
end

"""
save raw image
"""
function save(fimg::AbstractString, img::EMImage, dname::AbstractString="image")
    if isfile(fimg)
        rm(fimg)
    end
    f = h5open(fimg, "w")
    f[dname,"chunk", (8,8,2), "shuffle", (), "deflate", 3] = img
    close(f)
end

function saveimg(fimg::AbstractString, img::EMImage, dname::AbstractString="image")
    save(fimg, img, dname)
end

"""
permute dims of x and y for tif images
"""
function permutetifdims!(arr::Array)
    # transpose the X and Y
    p = Vector{Int64}( 1:ndims(arr) )
    p[1] = 2
    p[2] = 1
    arr = permutedims(arr, p)
    arr
end

"""
read segmentation
for directory, currently only tested with VAST output
"""
function readseg(fseg::AbstractString)
    if isdir(fseg)
        error("reading from directory is not correct now.")
        files_all = readdir(fseg)
        # only collect the tif images
        files = Vector{eltype(files_all)}()
        # get full path
        for file in files_all
            if contains(file,".tif")
                push!(files, joinpath(fseg, file))
            end
        end
        @assert length(files) > 0
        # read one tif and get size
        sz = length(files)
        image2d = load(files[1])
        if contains(image2d.properties["colorspace"], "RGB")
            tmp = reinterpret(UInt8, image2d.data)
            sc,sx,sy = size(tmp)
            seg = zeros(UInt32, (sx,sy,sz))
            for z in 1:sz
                tmp = reinterpret(UInt8, load(files[z]).data)
                tmp = Array{UInt32, 3}(tmp)
                im = tmp[1,:,:] .* 256 .* 256 + tmp[2,:,:] .* 256 + tmp[3,:,:]
                seg[:,:,z] = reshape(im,(sx,sy))
            end
        else
            tmp = reinterpret(UInt32, image2d.data)
            sx,sy = size(tmp)
            seg = zeros(UInt32,(sx,sy,sz))
            for z in 1:sz
                seg[:,:,z] = reinterpret(UInt32, load(files[z]).data)
            end
        end
        seg = permutetifdims!(seg)
    elseif ishdf5(fseg)
        f = h5open(fseg)
        if has(f, "seg")
          seg = read(f["seg"])
        elseif has(f, "segmentation")
          seg = read(f["segmentation"])
        else
          @assert has(f, "main")
          seg = read(f["main"])
        end
        close(f)
    elseif contains(fseg, ".tif")
        error("reading tif is not correct now..")
        image = load(fseg)
        if contains( image.properties["colorspace"], "RGB")
            tmp = reinterpret(UInt8,image.data)
            tmp = Array{UInt32, 4}(tmp)
            sc,sx,sy,sz = size(tmp)
            seg = tmp[1,:,:,:] .* 256 .* 256 + tmp[2,:,:,:] .* 256 + tmp[3,:,:,:]
            seg = reshape(seg, (sx,sy,sz))
        else
            seg = reinterpret(UInt32, image.data)
        end
        seg = permutetifdims!(seg)
    else
        error("unsupported file format!")
    end
    return Segmentation(seg)
end

"""
"""
function save(fseg::AbstractString, seg::Segmentation, dname::AbstractString="segmentation")
    if isfile(fseg)
        rm(fseg)
    end
    f = h5open(fseg, "w")
    f[dname,"chunk", (8,8,2), "shuffle", (), "deflate", 3] = seg
    close(f)
end

function saveseg(fseg::AbstractString, seg::Segmentation, dname::AbstractString="segmentation")
    save(fseg,seg, dname)
end

"""
read affinity map
"""
function readaff(faff::AbstractString)
    f = h5open(faff)
    if has(f, "aff")
      aff = read(f["aff"])
    elseif has(f, "affinityMap")
      aff = read(f["affinityMap"])
    else
      @assert has(f, "main")
      aff = read(f["main"])
    end
    close(f)
    return AffinityMap(aff)
end

"""
save affinity map
"""
function save(faff::AbstractString, aff::AffinityMap, dname::AbstractString="affinityMap")
    if isfile(faff)
        rm(faff)
    end
    f = h5open(faff, "w")
    f[dname,"chunk", (8,8,2,3), "shuffle", (), "deflate", 3] = aff
    close(f)
end

function saveaff(faff::AbstractString, aff::AffinityMap, dname::AbstractString="affinityMap")
    save(faff, aff, dname)
end

"""
whether a file is a sgm file
if the file do not exist, reture false
"""
function issgmfile(fname::AbstractString)
    if !isfile(fname)
        return false
    else
        f = h5open(fname)
        if has(f, "segmentPairs") || has(f, "segmentPair")
            return true
        else
            return false
        end
        close(f)
    end
end

"""
read segmentation with maximum spanning tree
"""
function readsgm(fname::AbstractString)
    f = h5open(fname)
    if has(f, "seg")
      seg = read(f["seg"])
      segmentPairs = read(f["segmentPairs"])
      segmentPairAffinities = read(f["segmentPairAffinities"])
    elseif has(f, "segmentation")
      seg = read(f["segmentation"])
      segmentPairs = read(f["segmentPairs"])
      segmentPairAffinities = read(f["segmentPairAffinities"])
    else
      @assert has(f, "main")
      seg = read(f["main"])
      segmentPairs = read(f["segmentPairs"])
      segmentPairAffinities = read(f["segmentPairAffinities"])
    end
    close(f)
    return SegMST(seg, segmentPairs, segmentPairAffinities)
end

"""
save segmentation with segmentPairsrogram
"""
function save(fsgm::AbstractString, sgm::SegMST)
    f = h5open(fsgm, "w")
    f["segmentation", "chunk", (64,64,8), "shuffle", (), "deflate", 3] = sgm.segmentation
    f["segmentPairs"] = sgm.segmentPairs
    f["segmentPairAffinities"] = sgm.segmentPairAffinities
    close(f)
end
function save(fsgm::AbstractString, seg::Segmentation, segmentPairs::SegmentPairs, segmentPairAffinities::SegmentPairAffinities)
    savesgm( fsgm, SegMST(seg,segmentPairs,segmentPairAffinities) )
end

function savesgm(fsgm::AbstractString, sgm)
    save(fsgm, sgm)
end

end # end of module


text:
using SimplePlotPerf
using Base.Test

# just make sure it runs
prplot(Dict(
    "test" => rand(100)
), bitrand(100))


text:
import LinearAlgebra: mul!, ldiv!
import Base: getindex, iterate

using SparseArrays

struct DiagonalIndices{Tv, Ti <: Integer}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::Vector{Ti}

    function DiagonalIndices{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
        # Check square?
        diag = Vector{Ti}(undef, A.n)

        for col = 1 : A.n
            r1 = A.colptr[col]
            r2 = A.colptr[col + 1] - 1
            r1 = searchsortedfirst(A.rowval, col, r1, r2, Base.Order.Forward)
            if r1 > r2 || A.rowval[r1] != col || iszero(A.nzval[r1])
                throw(LinearAlgebra.SingularException(col))
            end
            diag[col] = r1
        end

        new(A, diag)
    end
end

DiagonalIndices(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = DiagonalIndices{Tv,Ti}(A)

function ldiv!(y::AbstractVector{Tv}, D::DiagonalIndices{Tv,Ti}, x::AbstractVector{Tv}) where {Tv,Ti}
    @inbounds for row = 1 : D.matrix.n
        y[row] = x[row] / D.matrix.nzval[D.diag[row]]
    end
    y
end

@inline getindex(d::DiagonalIndices, i::Int) = d.diag[i]

struct FastLowerTriangular{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::DiagonalIndices{Tv,Ti}
end

struct FastUpperTriangular{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::DiagonalIndices{Tv,Ti}
end

struct StrictlyUpperTriangular{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::DiagonalIndices{Tv,Ti}
end

struct StrictlyLowerTriangular{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::DiagonalIndices{Tv,Ti}
end

struct OffDiagonal{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
    diag::DiagonalIndices{Tv,Ti}
end

"""
Forward substitution for the FastLowerTriangular type
"""
function forward_sub!(F::FastLowerTriangular, x::AbstractVector)
    A = F.matrix

    @inbounds for col = 1 : A.n

        # Solve for diagonal element
        idx = F.diag[col]
        x[col] /= A.nzval[idx]

        # Substitute next values involving x[col]
        for i = idx + 1 : (A.colptr[col + 1] - 1)
            x[A.rowval[i]] -= A.nzval[i] * x[col]
        end
    end

    x
end

"""
Forward substitution
"""
function forward_sub!(α, F::FastLowerTriangular, x::AbstractVector, β, y::AbstractVector)
    A = F.matrix

    @inbounds for col = 1 : A.n

        # Solve for diagonal element
        idx = F.diag[col]
        x[col] = α * x[col] / A.nzval[idx] + β * y[col]

        # Substitute next values involving x[col]
        for i = idx + 1 : (A.colptr[col + 1] - 1)
            x[A.rowval[i]] -= A.nzval[i] * x[col]
        end
    end

    x
end

"""
Backward substitution for the FastUpperTriangular type
"""
function backward_sub!(F::FastUpperTriangular, x::AbstractVector)
    A = F.matrix

    @inbounds for col = A.n : -1 : 1

        # Solve for diagonal element
        idx = F.diag[col]
        x[col] = x[col] / A.nzval[idx]

        # Substitute next values involving x[col]
        for i = A.colptr[col] : idx - 1
            x[A.rowval[i]] -= A.nzval[i] * x[col]
        end
    end

    x
end

function backward_sub!(α, F::FastUpperTriangular, x::AbstractVector, β, y::AbstractVector)
    A = F.matrix

    @inbounds for col = A.n : -1 : 1

        # Solve for diagonal element
        idx = F.diag[col]
        x[col] = α * x[col] / A.nzval[idx] + β * y[col]

        # Substitute next values involving x[col]
        for i = A.colptr[col] : idx - 1
            x[A.rowval[i]] -= A.nzval[i] * x[col]
        end
    end

    x
end

"""
Do mul! with the off-diagonal elements of a matrix.
"""
function mul!(α::T, O::OffDiagonal, x::AbstractVector, β::T, y::AbstractVector) where {T}
    # Specialize for β = 0 and β = 1
    A = O.matrix

    if β != one(T)
        if iszero(β)
            fill!(y, zero(T))
        else
            lmul!(β, y)
        end
    end

    @inbounds for col = 1 : A.n
        αx = α * x[col]
        diag_index = O.diag[col]
        for j = A.colptr[col] : diag_index - 1
            y[A.rowval[j]] += A.nzval[j] * αx
        end
        for j = diag_index + 1 : A.colptr[col + 1] - 1
            y[A.rowval[j]] += A.nzval[j] * αx
        end
    end

    y
end

"""
Computes z := α * U * x + β * y. Because U is StrictlyUpperTriangular
one can set z = x and update x in-place as x := α * U * x + β * y.
"""
function gauss_seidel_multiply!(α, U::StrictlyUpperTriangular, x::AbstractVector, β, y::AbstractVector, z::AbstractVector)
    A = U.matrix

    for col = 1 : A.n
        αx = α * x[col]
        diag_index = U.diag[col]
        @inbounds for j = A.colptr[col] : diag_index - 1
            z[A.rowval[j]] += A.nzval[j] * αx
        end
        z[col] = β * y[col]
    end
    z
end

"""
Computes z := α * L * x + β * y. Because A is StrictlyLowerTriangular
one can set z = x and update x in-place as x := α * L * x + β * y.
"""
function gauss_seidel_multiply!(α, L::StrictlyLowerTriangular, x::AbstractVector, β, y::AbstractVector, z::AbstractVector)
    A = L.matrix

    for col = A.n : -1 : 1
        αx = α * x[col]
        z[col] = β * y[col]
        @inbounds for j = L.diag[col] + 1 : (A.colptr[col + 1] - 1)
            z[A.rowval[j]] += A.nzval[j] * αx
        end
    end
    z
end

##
## Jacobi
##

mutable struct JacobiIterable{T <: Number, solT, vecT, rhsT}
    O::OffDiagonal

    x::solT
    next::vecT
    b::rhsT

    maxiter::Int
end

start(::JacobiIterable) = 1
done(j::JacobiIterable, iteration::Int) = iteration > j.maxiter
function iterate(j::JacobiIterable{T}, iteration::Int=start(j)) where {T}
    if done(j, iteration) return nothing end
    # tmp = D \ (b - (A - D) * x)
    copyto!(j.next, j.b)
    mul!(-one(T), j.O, j.x, one(T), j.next)
    ldiv!(j.x, j.O.diag, j.next)

    nothing, iteration + 1
end

jacobi_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int = 10) =
    JacobiIterable{eltype(x), typeof(x), typeof(x), typeof(b)}(
        OffDiagonal(A, DiagonalIndices(A)), x, similar(x), b, maxiter)

"""
    jacobi!(x, A::SparseMatrixCSC, b; maxiter=10) -> x

Performs exactly `maxiter` Jacobi iterations.

Allocates a temporary vector and precomputes the diagonal indices.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function jacobi!(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int=10)
    iterable = jacobi_iterable(x, A, b, maxiter = maxiter)
    for item = iterable end
    iterable.x
end

##
## Gauss-Seidel
##

mutable struct GaussSeidelIterable{solT, rhsT}
    U::StrictlyUpperTriangular
    L::FastLowerTriangular

    x::solT
    b::rhsT

    maxiter::Int
end

function gauss_seidel_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int = 10)
    D = DiagonalIndices(A)
    GaussSeidelIterable(StrictlyUpperTriangular(A, D), FastLowerTriangular(A, D), x, b, maxiter)
end

start(::GaussSeidelIterable) = 1
done(g::GaussSeidelIterable, iteration::Int) = iteration > g.maxiter
function iterate(g::GaussSeidelIterable, iteration::Int=start(g))
    if done(g, iteration) return nothing end
    # x ← L \ (-U * x + b)
    T = eltype(g.x)
    gauss_seidel_multiply!(-one(T), g.U, g.x, one(T), g.b, g.x)
    forward_sub!(g.L, g.x)

    nothing, iteration + 1
end

"""
    gauss_seidel!(x, A::SparseMatrixCSC, b; maxiter=10) -> x

Performs exactly `maxiter` Gauss-Seidel iterations.

Works fully in-place, but precomputes the diagonal indices.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function gauss_seidel!(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector; maxiter::Int = 10)
    iterable = gauss_seidel_iterable(x, A, b, maxiter = maxiter)
    for item = iterable end
    iterable.x
end

##
## SOR
##

mutable struct SORIterable{T, solT, vecT, rhsT, numT <: Real}
    U::StrictlyUpperTriangular
    L::FastLowerTriangular
    ω::numT

    x::solT
    next::vecT
    b::rhsT

    maxiter::Int
end

start(::SORIterable) = 1
done(s::SORIterable, iteration::Int) = iteration > s.maxiter
function iterate(s::SORIterable{T}, iteration::Int=start(s)) where {T}
    if done(s, iteration) return nothing end

    # next = b - U * x
    gauss_seidel_multiply!(-one(T), s.U, s.x, one(T), s.b, s.next)

    # next = ω * inv(L) * next + (1 - ω) * x
    forward_sub!(s.ω, s.L, s.next, one(T) - s.ω, s.x)

    # Switch current and next iterate
    s.x, s.next = s.next, s.x

    nothing, iteration + 1
end

function sor_iterable(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10)
    D = DiagonalIndices(A)
    T = eltype(x)
    SORIterable{T,typeof(x),typeof(x),typeof(b),eltype(ω)}(
        StrictlyUpperTriangular(A, D), FastLowerTriangular(A, D), ω,
        x, similar(x), b, maxiter
    )
end

"""
    sor!(x, A::SparseMatrixCSC, b, ω::Real; maxiter=10)

Performs exactly `maxiter` SOR iterations with relaxation parameter `ω`.

Allocates a temporary vector and precomputes the diagonal indices.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function sor!(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10)
    iterable = sor_iterable(x, A, b, ω, maxiter = maxiter)
    for item = iterable end
    iterable.x
end

##
## SSOR
##

mutable struct SSORIterable{T, solT, vecT, rhsT, numT <: Real}
    sL::StrictlyLowerTriangular
    sU::StrictlyUpperTriangular
    L::FastLowerTriangular
    U::FastUpperTriangular
    ω::numT
    x::solT
    tmp::vecT
    b::rhsT
    maxiter::Int
end

function ssor_iterable(x::AbstractVector{T}, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10) where {T}
    D = DiagonalIndices(A)
    SSORIterable{T,typeof(x),typeof(x),typeof(b),typeof(ω)}(
        StrictlyLowerTriangular(A, D),
        StrictlyUpperTriangular(A, D),
        FastLowerTriangular(A, D),
        FastUpperTriangular(A, D),
        ω, x, similar(x), b, maxiter
    )
end

start(s::SSORIterable) = 1
done(s::SSORIterable, iteration::Int) = iteration > s.maxiter

function iterate(s::SSORIterable{T}, iteration::Int=start(s)) where {T}
    if done(s, iteration) return nothing end

    # tmp = b - U * x
    gauss_seidel_multiply!(-one(T), s.sU, s.x, one(T), s.b, s.tmp)

    # tmp = ω * inv(L) * tmp + (1 - ω) * x
    forward_sub!(s.ω, s.L, s.tmp, one(T) - s.ω, s.x)

    # x = b - L * tmp
    gauss_seidel_multiply!(-one(T), s.sL, s.tmp, one(T), s.b, s.x)

    # x = ω * inv(U) * x + (1 - ω) * tmp
    backward_sub!(s.ω, s.U, s.x, one(T) - s.ω, s.tmp)

    nothing, iteration + 1
end

"""
    ssor!(x, A::SparseMatrixCSC, b, ω::Real; maxiter=10)

Performs exactly `maxiter` SSOR iterations with relaxation parameter `ω`. Each iteration
is basically a forward *and* backward sweep of SOR.

Allocates a temporary vector and precomputes the diagonal indices.

Throws `LinearAlgebra.SingularException` when the diagonal has a zero. This check
is performed once beforehand.
"""
function ssor!(x::AbstractVector, A::SparseMatrixCSC, b::AbstractVector, ω::Real; maxiter::Int = 10)
    iterable = ssor_iterable(x, A, b, ω, maxiter = maxiter)
    for item = iterable end
    iterable.x
end


text:
#!/usr/bin/env julia

export pe15

pe15(n::Integer=20)=binomial(2n,n)

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe15()|>print
end

max_stars_repo_path: 15. Lattice paths/15.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 15. Lattice paths/15.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 15. Lattice paths/15.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
module SIMDHints

export simd, ivdep

using Base.Broadcast: Broadcasted, BroadcastStyle, throwdm, preprocess

abstract type SIMDHint end

struct SIMD{T} <: SIMDHint
    value::T
end

struct IVDEP{T} <: SIMDHint
    value::T
end

simd(x) = SIMD(x)
ivdep(x) = IVDEP(x)

struct SIMDStyle{f} <: BroadcastStyle end

const HintFuncType = Union{typeof(simd), typeof(ivdep)}

Broadcast.broadcasted(f::HintFuncType, x) =
    Broadcasted{SIMDStyle{f}}(identity, (x,))

@inline function Base.copyto!(dest::AbstractArray, bc::Broadcasted{SIMDStyle{f}}) where f
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    @assert bc.f === identity
    @assert bc.args isa Tuple{Any}
    bc′ = preprocess(dest, bc.args[1])
    if f === ivdep
        @simd ivdep for I in eachindex(bc′)
            @inbounds dest[I] = bc′[I]
        end
    else
        @simd for I in eachindex(bc′)
            @inbounds dest[I] = bc′[I]
        end
    end
    return dest
end

end # module