text:
function crandn(dims...)
    sqrt(0.5) * ( randn(dims) + 1im*randn(dims) )
end
"""
    scm_channel(angles,weights,nAntennas::Integer;nCoherence=1,AS=2.0)

Generate `nCoherence` channel realization according to angular power density (spectrum) `f`
where `f` is given as a superposition of weighted Laplace kernels located at `angles` and
with standard deviation `AS`
The spectrum `f` defines a covariance matrix of the channel `h`.
The first row of this covariance matrix is returned as second output.
"""
function scm_channel(angles,weights,nAntennas::Integer;nCoherence=1,AS=2.0)
    f(x)  = spectrum(x,angles,weights;AS=AS)
    (h,t) = chan_from_spectrum(f,nAntennas,nCoherence=nCoherence)
end

"""
    chan_from_spectrum(f,nAntennas::Integer;nCoherence=1)

Generate `nCoherence` channel realization according to angular power density (spectrum) `f`
(which is a cont. function).
The spectrum `f` defines a covariance matrix of the channel `h`.
The first row of this covariance matrix is returned as second output.
"""
function chan_from_spectrum(f,nAntennas::Integer;nCoherence=1)
    OF = 20 # oversampling factor (ideally, would use continuous freq. spectrum...)
    nFreqSamples = OF*nAntennas
    fs = Vector{Float64}([f(i/nFreqSamples*2*pi) for i=0:nFreqSamples-1])

    # avoid instabilities due to almost infinite energy at some frequencies
    # (this should only happen at "endfire" of a uniform linear array where --
    # because of the arcsin-transform -- the angular psd grows to infinity
    almostInfThreshold = max(1,nFreqSamples) # use nFreqSamples as threshold value...
    almostInfFreqs = abs.(fs) .> almostInfThreshold
    fs[almostInfFreqs] = almostInfThreshold * angle.(fs[almostInfFreqs])

    if sum(fs) > 0
        fs = fs./sum(fs)*nFreqSamples # normalize energy
    end

    x = crandn(nFreqSamples,nCoherence);
    h = ifft(sqrt.(fs).*x,1).*sqrt(nFreqSamples)
    h = h[1:nAntennas,:]

    # t is the first row of the covariance matrix of h (which is Toeplitz and Hermitian)
    t = fft(fs)./nFreqSamples
    t = t[1:nAntennas]
    (h,t)
end

"""
    spectrum(x_rad,angles,weights,AS=2.0)

Generate spectrum at values `asin(x_rad/pi) [deg]` corresponding
to superposition of weighted Laplace kernels with means `angles`
and standard deviations `AS`.
This functions takes care of the folding that happens when transforming
from [-180,180] deg. to [rad]-space...
"""
function spectrum(x_rad,angles,weights;AS=2.0)
    x_rad = mod(x_rad + pi,2*pi)-pi
    x_deg = asin(x_rad/pi)./pi*180

    v =     laplace(    x_deg,angles,weights,AS=AS)
    v = v + laplace(180-x_deg,angles,weights,AS=AS)
    v = v .* (180/pi*2*pi./sqrt(pi^2 - x_rad.^2))
end

"""
    laplace(x_deg::Float64,angles,weights;AS=2.0)

Value at `x_deg` (in deg) of sum of weighted Laplace kernels
(kernels are wrapped to [-180,180] degrees).
"""
function laplace(x_deg::Float64,angles,weights;AS=2.0)
    v = 0.0;
    for i in 1:length(weights)
        xshifted = x_deg - angles[i]
        xshifted = mod(xshifted + 180,360)-180;
        v += weights[i]/(2*AS)*exp(-abs(xshifted)/AS)
    end
    v
end

"""
    laplace(x_deg::Vector{Float64},angles,weights;AS=2.0)

Value at `x_deg` (in deg) of sum of weighted Laplace kernels
(kernels are wrapped to [-180,180] degrees).
"""
function laplace(x_deg::Vector{Float64},angles,weights;AS=2.0)
    v = zeros(x_deg)
    for j in 1:length(x_deg)
        vj = 0.0;
        for i in 1:length(weights)
            xshifted = x_deg[j] - angles[i]
            xshifted = mod(xshifted + 180,360)-180;
            vj += weights[i]/(2*AS)*exp(-abs.(xshifted)/AS)
        end
        v[j] = vj;
    end
    v
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
using StanSample, Test

#set_cmdstan_home!(homedir() * "/Projects/StanSupport/cmdstan_stanc3")

ProjDir = @__DIR__
cd(ProjDir) # do

  bernoulli_model = "
    functions{
      #include model_specific_funcs.stan
      #include shared_funcs.stan // a comment
      //#include shared_funcs.stan // a comment
    }
    data {
      int<lower=1> N;
      int<lower=0,upper=1> y[N];
    }
    parameters {
      real<lower=0,upper=1> theta;
    }
    model {
      model_specific_function();
      theta ~ beta(my_function(),1);
      y ~ bernoulli(theta);
    }
  ";

  sm = SampleModel("bernoulli", bernoulli_model)
  observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
  rc = stan_sample(sm; data=observeddata)

  if success(rc)
    samples = read_samples(sm, :array)
    @test sum(samples)/length(samples) ≈ 0.33 rtol=0.05
  end

#end

ext: jl

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
mutable struct Model
    env::Env # Cplex environment
    lp::Ptr{Cvoid} # Cplex problem (lp)
    has_int::Bool # problem has integer variables?
    has_qc::Bool # problem has quadratic constraints?
    has_sos::Bool # problem has Special Ordered Sets?
    callback::Any
    terminator::Vector{Cint}
end
function Model(env::Env, lp::Ptr{Cvoid})
    notify_new_model(env)
    model = Model(env, lp, false, false, false, nothing, Cint[0])
    function model_finalizer(model)
        if model.lp != C_NULL
            free_problem(model)
        else
            # User must have called `free_problem` directly.
        end
        notify_freed_model(env)
    end
    finalizer(model_finalizer, model)
    set_terminate(model)
    return model
end

function Model(env::Env, name::String="CPLEX.jl")
    @assert is_valid(env)
    stat = Vector{Cint}(undef, 1)
    tmp = @cpx_ccall(createprob, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cchar}), env.ptr, stat, name)
    if tmp == C_NULL
        throw(CplexError(env, stat))
    end
    return Model(env, tmp)
end

function read_model(model::Model, filename::String)
    stat = @cpx_ccall(readcopyprob, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cchar}), model.env.ptr, model.lp, filename, C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    prob_type = get_prob_type(model)
    if prob_type in [:MILP,:MIQP, :MIQCP]
        model.has_int = true
    end
    if prob_type in [:QP, :MIQP, :QCP, :MIQCP]
        model.has_qc = true
    end
end

function write_model(model::Model, filename::String)
    if endswith(filename,".mps")
        filetype = "MPS"
    elseif endswith(filename,".lp")
        filetype = "LP"
    else
        error("Unrecognized file extension: $filename (Only .mps and .lp are supported)")
    end
    stat = @cpx_ccall(writeprob, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cchar}, Ptr{Cchar}), model.env.ptr, model.lp, filename, filetype)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

## TODO: deep copy model, reset model

function c_api_getobjsen(model::Model)
    sense_int = @cpx_ccall(getobjsen, Cint, (
                           Ptr{Cvoid},
                           Ptr{Cvoid},
                           ),
                           model.env.ptr, model.lp)

    return sense_int
end
function get_sense(model::Model)
    sense_int = c_api_getobjsen(model)
    if sense_int == 1
        return :Min
    elseif sense_int == -1
        return :Max
    else
        error("CPLEX: problem object or environment does not exist")
    end
end

function set_sense!(model::Model, sense)
    if sense == :Min
        @cpx_ccall(chgobjsen, Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Cint), model.env.ptr, model.lp, 1)
    elseif sense == :Max
        @cpx_ccall(chgobjsen, Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Cint), model.env.ptr, model.lp, -1)
    else
        error("Unrecognized objective sense $sense")
    end
end

function c_api_chgobjsen(model::Model, sense_int::Cint)
    @cpx_ccall(chgobjsen, Nothing, (Ptr{Cvoid}, Ptr{Cvoid}, Cint),
               model.env.ptr, model.lp, sense_int)
end

function c_api_getobj(model::Model, sized_obj::FVec,
                      col_start::Cint, col_end::Cint)

    nvars = num_var(model)
    stat = @cpx_ccall(getobj, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, sized_obj,
                      col_start - Cint(1), col_end - Cint(1))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function get_obj(model::Model)
    nvars = num_var(model)
    obj = Vector{Cdouble}(undef, nvars)
    stat = @cpx_ccall(getobj, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, obj, 0, nvars-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return obj
end

const type_map = Dict(
     0 => :LP,
     1 => :MILP,
     3 => :FIXEDMILP, # actually fixed milp
     5 => :QP,
     7 => :MIQP,
     8 => :FIXEDMIQP,
    10 => :QCP,
    11 => :MIQCP
)

const rev_prob_type_map = Dict(
    :LP    => 0,
    :MILP  => 1,
    :FIXEDMILP  => 3,
    :QP    => 5,
    :MIQP  => 7,
    :FIXEDMIQP  => 8,
    :QCP   => 10,
    :MIQCP => 11
)

function get_prob_type(model::Model)
  ret = @cpx_ccall(getprobtype, Cint, (
                   Ptr{Cvoid},
                   Ptr{Cvoid}),
                   model.env.ptr, model.lp)
  ret == -1 && error("No problem of environment")
  return type_map[Int(ret)]
end

function set_prob_type!(model::Model, tyint::Int)
    stat = @cpx_ccall(chgprobtype, Cint, (
                     Ptr{Cvoid},
                     Ptr{Cvoid},
                     Cint),
                     model.env.ptr, model.lp, tyint)
     if stat != 0
         throw(CplexError(model.env, stat))
     end
end
set_prob_type!(model::Model, ty::Symbol) = set_prob_type!(model, rev_prob_type_map[ty])

function set_obj!(model::Model, c::Vector)
    nvars = num_var(model)
    stat = @cpx_ccall(chgobj, Cint, (
                        Ptr{Cvoid},
                        Ptr{Cvoid},
                        Cint,
                        Ptr{Cint},
                        Ptr{Cdouble}
                        ),
                        model.env.ptr, model.lp, nvars, Cint[0:nvars-1;], float(c))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function c_api_chgobj(model::Model, indices::IVec, values::FVec)
    nvars = length(indices)
    stat = @cpx_ccall(chgobj, Cint, (
                        Ptr{Cvoid},
                        Ptr{Cvoid},
                        Cint,
                        Ptr{Cint},
                        Ptr{Cdouble}
                        ),
                        model.env.ptr, model.lp, nvars,
                        indices .- Cint(1), values)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function c_api_chgobjoffset(model::Model, offset::Float64)
    stat = @cpx_ccall(
        chgobjoffset,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cdouble),
        model.env.ptr, model.lp, offset
    )
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return
end

function c_api_getobjoffset(model::Model)
    objoffset_p = Ref{Float64}()
    stat = @cpx_ccall(
        getobjoffset,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cdouble}),
        model.env.ptr, model.lp, objoffset_p
    )
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return objoffset_p[]
end

set_warm_start!(model::Model, x::Vector{Float64}, effortlevel::Integer = CPX_MIPSTART_AUTO) = set_warm_start!(model, Cint[1:length(x);], x, effortlevel)

function set_warm_start!(model::Model, indx::IVec, val::FVec, effortlevel::Integer)
    stat = @cpx_ccall(addmipstarts, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Cint,
                      Cint,
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Ptr{Cint},
                      Ptr{Ptr{Cchar}}
                      ),
                      model.env.ptr, model.lp, 1, length(indx), Cint[0], indx .- Cint(1), val, Cint[effortlevel], C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function free_problem(model::Model)
    stat = @cpx_ccall(freeprob, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    model.lp = C_NULL
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function set_terminate(model::Model)
    stat = @cpx_ccall(setterminate, Cint, (Ptr{Cvoid},Ptr{Cint}), model.env.ptr, model.terminator)
    if stat != 0
        throw(CplexError(env, stat))
    end
end

terminate(model::Model) = (model.terminator[1] = 1)

mutable struct ConflictRefinerData
    stat::Int
    nrows::Int # number of rows participating in the conflict
    rowind::Vector{Cint} # index of the rows that participate
    rowstat::Vector{Cint} # state of the rows that participate
    ncols::Int # number of columns participating in the conflict
    colind::Vector{Cint} # index of the columns that participate
    colstat::Vector{Cint} # state of the columns that participates
end

function c_api_getconflict(model::Model)
    # This function always calls refineconflict first, which starts the conflict refiner.
    # In other words, any call to this function is expensive.

    # First, compute the conflict.
    confnumrows_p = Ref{Cint}()
    confnumcols_p = Ref{Cint}()
    stat = @cpx_ccall(
        refineconflict,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}),
        model.env.ptr, model.lp, confnumrows_p, confnumcols_p)
    if stat != 0
        throw(CplexError(model.env, stat))
    end

    # Then, retrieve it.
    confstat_p = Ref{Cint}()
    rowind = Vector{Cint}(undef, confnumrows_p[])
    rowbdstat = Vector{Cint}(undef, confnumrows_p[])
    confnumrows_p = Ref{Cint}()
    colind = Vector{Cint}(undef, confnumcols_p[])
    colbdstat = Vector{Cint}(undef, confnumcols_p[])
    confnumcols_p  = Ref{Cint}()
    stat = @cpx_ccall(
        getconflict,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
        model.env.ptr, model.lp, confstat_p, rowind, rowbdstat, confnumrows_p, colind, colbdstat, confnumcols_p)
    if stat != 0
        throw(CplexError(model.env, stat))
    end

    return ConflictRefinerData(confstat_p[], confnumrows_p[], rowind, rowbdstat, confnumcols_p[], colind, colbdstat)
end

function c_api_chgname(model::Model, key::Cchar, ij::Cint, name::String)
    stat = @cpx_ccall(
        chgname,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, Ptr{Cchar}),
        model.env.ptr, model.lp, key, ij, name
    )
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdstatistic

"""
    'feffery'_antdstatistic(;kwargs...)

An AntdStatistic component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `precision` (Real; optional)
- `prefix` (optional): . prefix has the following type: lists containing elements 'mode', 'content'.
Those elements have the following types:
  - `mode` (a value equal to: 'text', 'icon'; optional)
  - `content` (String; optional)
- `showGroupSeparator` (Bool; optional)
- `style` (Dict; optional)
- `suffix` (optional): . suffix has the following type: lists containing elements 'mode', 'content'.
Those elements have the following types:
  - `mode` (a value equal to: 'text', 'icon'; optional)
  - `content` (String; optional)
- `title` (String; optional)
- `titleTooltip` (String; optional)
- `value` (Real; optional)
- `valueStyle` (Dict; optional)
"""
function 'feffery'_antdstatistic(; kwargs...)
        available_props = Symbol[:id, :className, :loading_state, :precision, :prefix, :showGroupSeparator, :style, :suffix, :title, :titleTooltip, :value, :valueStyle]
        wild_props = Symbol[]
        return Component("'feffery'_antdstatistic", "AntdStatistic", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
module SokobanUndirectedModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase

#####
##### game logic
#####

const NUM_OBJECTS = 4
const AGENT = 1
const WALL = 2
const BOX = 3
const TARGET = 4
const NUM_ACTIONS = 4

struct LevelDataset
    characters::Array{Char, 3}
end

function LevelDataset(file::AbstractString)
    lines = readlines(file)

    height = 0
    width = 0
    num_levels = 0

    for i in 1 : length(lines)
        if lines[i] == ""
            height = i - 2
            width = length(lines[i - 1])
            num_levels = length(lines) ÷ (height + 2)
            break
        end
    end

    characters = Array{Char}(undef, height, width, num_levels)

    for k in 1 : num_levels
        for j in 1 : width
            for i in 1 : height
                characters[i, j, k] = lines[(k - 1) * (height + 2) + 1 + i][j]
            end
        end
    end

    return LevelDataset(characters)
end

"""
The dataset (000.txt) used is [boxoban-levels](https://github.com/deepmind/boxoban-levels/blob/master/medium/train/000.txt) (LICENSE file included along with 000.txt)

The following is the citation of the dataset:

@misc{boxobanlevels,
author = {Arthur Guez, Mehdi Mirza, Karol Gregor, Rishabh Kabra, Sebastien Racaniere, Theophane Weber, David Raposo, Adam Santoro, Laurent Orseau, Tom Eccles, Greg Wayne, David Silver, Timothy Lillicrap, Victor Valdes},
title = {An investigation of Model-free planning: boxoban levels},
howpublished= {https://github.com/deepmind/boxoban-levels/},
year = "2018",
}
"""
mutable struct SokobanUndirected{R, RNG} <: GW.AbstractGridWorld
    tile_map::BitArray{3}
    agent_position::CartesianIndex{2}
    reward::R
    rng::RNG
    done::Bool
    box_reward::R
    box_positions::Vector{CartesianIndex{2}}
    level_number::Int
    dataset::LevelDataset
end

function SokobanUndirected(; R = Float32, file = joinpath(dirname(functionloc(GW.eval)[1]), "envs/sokoban/boxoban-levels/medium/train/000.txt"), rng = Random.GLOBAL_RNG)
    dataset = LevelDataset(file)
    height, width, _ = size(dataset.characters)

    tile_map = falses(NUM_OBJECTS, height, width)
    agent_position = CartesianIndex(1, 1)

    reward = zero(R)
    done = false
    box_reward = one(R)
    box_positions = CartesianIndex{2}[]
    level_number = 1

    env = SokobanUndirected(tile_map, agent_position, reward, rng, done, box_reward, box_positions, level_number, dataset)

    GW.reset!(env)

    return env
end

function GW.reset!(env::SokobanUndirected)
    rng = env.rng
    num_levels = size(env.dataset.characters, 3)

    new_level_number = rand(rng, 1 : num_levels)

    set_level!(env, new_level_number)

    env.reward = zero(env.reward)
    env.done = false

    return nothing
end

function GW.act!(env::SokobanUndirected, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    tile_map = env.tile_map
    box_positions = env.box_positions
    box_reward = env.box_reward

    r1 = sum(pos -> tile_map[TARGET, pos], box_positions) * box_reward
    agent_position = env.agent_position

    if action == 1
        dest = GW.move_up(agent_position)
        beyond_dest = GW.move_up(dest)
    elseif action == 2
        dest = GW.move_down(agent_position)
        beyond_dest = GW.move_down(dest)
    elseif action == 3
        dest = GW.move_left(agent_position)
        beyond_dest = GW.move_left(dest)
    else
        dest = GW.move_right(agent_position)
        beyond_dest = GW.move_right(dest)
    end

    if !tile_map[WALL, dest]
        if !tile_map[BOX, dest]
            tile_map[AGENT, agent_position] = false
            env.agent_position = dest
            tile_map[AGENT, dest] = true
        else
            if !tile_map[BOX, beyond_dest] && !tile_map[WALL, beyond_dest]
                tile_map[BOX, dest] = false
                tile_map[BOX, beyond_dest] = true

                box_idx = findfirst(==(dest), box_positions)
                box_positions[box_idx] = beyond_dest
                tile_map[AGENT, agent_position] = false
                env.agent_position = dest
                tile_map[AGENT, dest] = true
            end
        end
    end

    env.done = all(pos -> tile_map[TARGET, pos], box_positions)
    r2 = sum(pos -> tile_map[TARGET, pos], box_positions) * box_reward
    env.reward = r2 - r1

    return nothing
end

function set_level!(env::SokobanUndirected, level_number::Int)
    tile_map = env.tile_map
    characters = env.dataset.characters
    box_positions = env.box_positions

    fill!(tile_map, false)
    empty!(env.box_positions)

    env.level_number = level_number

    height, width, num_levels = size(characters)

    for j in 1 : width
        for i in 1 : height
            char = characters[i, j, level_number]
            if char === '@'
                tile_map[AGENT, i, j] = true
                env.agent_position = CartesianIndex(i, j)
            elseif char === '#'
                tile_map[WALL, i, j] = true
            elseif char === '$'
                tile_map[BOX, i, j] = true
                push!(box_positions, CartesianIndex(i, j))
            elseif char === '.'
                tile_map[TARGET, i, j] = true
            end
        end
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::SokobanUndirected) = size(env.tile_map, 2)
GW.get_width(env::SokobanUndirected) = size(env.tile_map, 3)

GW.get_action_names(env::SokobanUndirected) = (:MOVE_UP, :MOVE_DOWN, :MOVE_LEFT, :MOVE_RIGHT)
GW.get_object_names(env::SokobanUndirected) = (:AGENT, :WALL, :BOX, :TARGET)

function GW.get_pretty_tile_map(env::SokobanUndirected, position::CartesianIndex{2})
    characters = ('☻', '█', '▒', '✖', '⋅')

    object = findfirst(@view env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function GW.get_pretty_sub_tile_map(env::SokobanUndirected, window_size, position::CartesianIndex{2})
    tile_map = env.tile_map
    agent_position = env.agent_position

    characters = ('☻', '█', '▒', '✖', '⋅')

    sub_tile_map = GW.get_sub_tile_map(tile_map, agent_position, window_size)

    object = findfirst(@view sub_tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

function Base.show(io::IO, ::MIME"text/plain", env::SokobanUndirected)
    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nsub_tile_map:\n"
    str = str * GW.get_pretty_sub_tile_map(env, GW.get_window_size(env))
    str = str * "\nreward: $(env.reward)"
    str = str * "\ndone: $(env.done)"
    str = str * "\nlevel_number: $(env.level_number)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::SokobanUndirected) = ('w', 's', 'a', 'd')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: SokobanUndirected} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SokobanUndirected} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: SokobanUndirected} = env.env.tile_map

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: SokobanUndirected} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: SokobanUndirected} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: SokobanUndirected} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: SokobanUndirected} = env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: SokobanUndirected} = env.env.done

end # module


text:
using .DataFrames

function levene(df::DataFrame, observationscolumn::Symbol, factorcolumns::Vector{Symbol})
    observations = df[!, observationscolumn]
    length(observations) > 0 || return
    eltype(observations) <: Number || error("Obervations must be numeric")
    levene(df[!, observationscolumn], [df[!, x] for x ∈ factorcolumns])
end


text:
# # Huber regression

# This example can be found here: <https://web.stanford.edu/~boyd/papers/pdf/cvx_applications.pdf>.
# Here we set `big_example = false` to only generate a small example which takes less time to run.
big_example = false
if big_example
    n = 300
    number_tests = 50
else
    n = 50
    number_tests = 10
end

# Generate data for Huber regression.
using Random
Random.seed!(1);
number_samples = round(Int, 1.5 * n);
beta_true = 5 * randn(n);
X = randn(n, number_samples);
Y = zeros(number_samples);
v = randn(number_samples);

#-

## Generate data for different values of p.
## Solve the resulting problems.
using Convex, SCS, Distributions
lsq_data = zeros(number_tests);
huber_data = zeros(number_tests);
prescient_data = zeros(number_tests);
p_vals = range(0, stop = 0.15, length = number_tests);
for i in 1:length(p_vals)
    p = p_vals[i]
    ## Generate the sign changes.
    factor = 2 * rand(Binomial(1, 1 - p), number_samples) .- 1
    Y = factor .* X' * beta_true + v

    ## Form and solve a standard regression problem.
    beta = Variable(n)
    fit = norm(beta - beta_true) / norm(beta_true)
    cost = norm(X' * beta - Y)
    prob = minimize(cost)
    solve!(prob, SCS.Optimizer; silent_solver = true)
    lsq_data[i] = evaluate(fit)

    ## Form and solve a prescient regression problem,
    ## i.e., where the sign changes are known.
    cost = norm(factor .* (X' * beta) - Y)
    solve!(minimize(cost), SCS.Optimizer; silent_solver = true)
    prescient_data[i] = evaluate(fit)

    ## Form and solve the Huber regression problem.
    cost = sum(huber(X' * beta - Y, 1))
    solve!(minimize(cost), SCS.Optimizer; silent_solver = true)
    huber_data[i] = evaluate(fit)
end

#-
using Plots

plot(p_vals, huber_data, label = "Huber", xlabel = "p", ylabel = "Fit")
plot!(p_vals, lsq_data, label = "Least squares")
plot!(p_vals, prescient_data, label = "Prescient")
#-

## Plot the relative reconstruction error for Huber and prescient regression,
## zooming in on smaller values of p.
indices = findall(p_vals .<= 0.08);
plot(p_vals[indices], huber_data[indices], label = "Huber")
plot!(p_vals[indices], prescient_data[indices], label = "Prescient")


text:
using Optim
using PyPlot
#import celerite
include("../src/celerite.jl")

srand(42)

# The input coordinates must be sorted
t = sort(cat(1, 3.8 * rand(57), 5.5 + 4.5 * rand(68)))
yerr = 0.08 + (0.22-0.08)*rand(length(t))
y = 0.2*(t-5.0) + sin.(3.0*t + 0.1*(t-5.0).^2) + yerr .* randn(length(t))

true_t = linspace(0, 10, 1000)
true_y = 0.2*(true_t-5) + sin.(3*true_t + 0.1*(true_t-5).^2)

PyPlot.clf()
PyPlot.plot(true_t, true_y, "k", lw=1.5, alpha=0.3)
PyPlot.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0)
PyPlot.xlabel("x")
PyPlot.ylabel("y")
PyPlot.xlim(0, 10)
PyPlot.ylim(-2.5, 2.5);

Q = 1.0 / sqrt(2.0)
w0 = 3.0
S0 = var(y) / (w0 * Q)
kernel = celerite.SHOTerm(log(S0), log(Q), log(w0))

# A periodic component
Q = 1.0
w0 = 3.0
S0 = var(y) / (w0 * Q)
kernel = kernel + celerite.SHOTerm(log(S0), log(Q), log(w0))

celerite.TermSum((celerite.SHOTerm(-0.7432145976901582,-0.34657359027997275,1.0986122886681096),celerite.SHOTerm(-1.089788187970131,0.0,1.0986122886681096)))

gp = celerite.Celerite(kernel)
celerite.compute!(gp, t, yerr)
celerite.log_likelihood(gp, y)

#mu, variance = celerite.predict_full_ldlt(gp, y, true_t, return_var=true)
#mu, variance = celerite.predict_full(gp, y, true_t, return_var=true)
mu, variance = celerite.predict(gp, y, true_t, return_var=true)
#mu  = celerite.predict_full(gp, y, true_t, return_cov=false, return_var=false)
#mu  = celerite.predict!(gp, gp.x, y, true_t)
sigma = sqrt.(variance)

PyPlot.plot(true_t, true_y, "k", lw=1.5, alpha=0.3,label="Simulated data")
PyPlot.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0)
PyPlot.plot(true_t, mu, "g",label="Initial kernel")
PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
PyPlot.xlabel("x")
PyPlot.ylabel("y")
PyPlot.xlim(0, 10)
PyPlot.ylim(-2.5, 2.5);
read(STDIN,Char)

PyPlot.clf()
vector = celerite.get_parameter_vector(gp.kernel)
mask = ones(Bool, length(vector))
mask[2] = false  # Don't fit for the first Q
function nll(params)
    vector[mask] = params
    celerite.set_parameter_vector!(gp.kernel, vector)
    celerite.compute!(gp, t, yerr)
    return -celerite.log_likelihood(gp, y)
end;

result = Optim.optimize(nll, vector[mask], Optim.LBFGS())
result

vector[mask] = Optim.minimizer(result)
vector

celerite.set_parameter_vector!(gp.kernel, vector)

#mu, variance = celerite.predict(gp, y, true_t, return_var=true)
#mu, variance = celerite.predict_full_ldlt(gp, y, true_t, return_var=true)
#mu = celerite.predict(gp, y, true_t, return_cov = false, return_var=false)
#mu = celerite.predict!(gp, gp.x, y, true_t)
mu, variance = celerite.predict(gp, y, true_t, return_var=true)
sigma = sqrt.(variance)

PyPlot.plot(true_t, true_y, "k", lw=1.5, alpha=0.3)
PyPlot.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0)
PyPlot.plot(true_t, mu, "r",label="Optimized kernel")
PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
PyPlot.xlabel("x")
PyPlot.ylabel("y")
PyPlot.xlim(0, 10)
PyPlot.ylim(-2.5, 2.5);
PyPlot.legend(loc="upper left")

read(STDIN,Char)
PyPlot.clf()

omega = exp.(linspace(log(0.1), log(20), 5000))
psd = celerite.get_psd(gp.kernel, omega)

for term in gp.kernel.terms
    PyPlot.plot(omega, celerite.get_psd(term, omega), "--", color="orange")
end
PyPlot.plot(omega, psd, color="orange")

PyPlot.yscale("log")
PyPlot.xscale("log")
PyPlot.xlim(omega[1], omega[end])
PyPlot.xlabel("omega")
PyPlot.ylabel("S(omega)");
PyPlot.plot(2*pi/2.0*[1.,1.],[1e-8,1e2])


text:
using OteraEngine
using Test

@testset "OteraEngine.jl" begin
    #When initializing variables directly
    @test_nowarn Template("test1.html")
    tmp = Template("test1.html")
    result = ""
    open("test2.html", "r") do f
        result = read(f, String)
    end
    init = Dict("usr"=>"Julia")
    @test result == tmp(init)

    #In case of no params and specified config
    tmp = Template("test3.html", config=Dict("code_block_start"=>"{{", "code_block_stop"=>"}}"))
    regex = r"[\s\S]*?([0-9]{2})[\s\S]*?"
    m = match(regex, tmp())
    @test 0 <= m.offset <= 100

    #include config file
    @test_throws ArgumentError Template("test3.html", config_path="test.conf")
    tmp = Template("test3.html", config_path="test.toml")
    regex = r"[\s\S]*?([0-9]{2})[\s\S]*?"
    m = match(regex, tmp())
    @test 0 <= m.offset <= 100
end


text:
using Recommendation
using Test
using SparseArrays

include("test_base_recommender.jl")
include("test_data_accessor.jl")

include("baseline/test_user_mean.jl")
include("baseline/test_item_mean.jl")
include("baseline/test_most_popular.jl")
include("baseline/test_threshold_percentage.jl")
include("baseline/test_co_occurrence.jl")

include("model/test_tf_idf.jl")
include("model/test_user_knn.jl")
include("model/test_item_knn.jl")
include("model/test_svd.jl")
include("model/test_mf.jl")

include("metric/test_accuracy.jl")
include("metric/test_ranking.jl")

include("evaluation/test_evaluate.jl")
include("evaluation/test_cross_validation.jl")


text:
# function Sem(;
#         semobs::O,
#         imply = nothing,
#         loss = nothing,
#         semdiff = nothing,
#         algorithm = nothing
#         ) where {
#             O <: SemObs
#         }
#     if isnothing(imply) imply = ImplyCommon()
# end

### Two versions. The second one could be easier to construct
function (model::Sem)(par)
    model.imply(par, model)
    F = model.loss(par, model)
    return F
end

# function (model::Sem)(E, G, par)
#     model.imply(par)
#     if G != nothing
#         model.loss(par, model, E, G)
#     end
#     if E != nothing
#         F = model.loss(par, model, E, G)
#         return F
#     end
# end

# function (model::Sem)(E, G, par)
#     model.imply(par)
#     model.loss(par, model, E, G)
# end

function computeloss(model, par)
    model.imply(par)
    F = model.loss(par, model)
    return F
end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, Random, LinearAlgebra
using LevyArea, Statistics

seed = 113213898
Random.seed!(seed)

m = 10
W = WienerProcess(0.0,zeros(m),nothing)

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

# LevyArea.jl algs: ITER_INT_ALGS

@testset "diagonal noise tests" begin
  true_diag = 1//2 .* W.dW .* W.dW

  Jdiag = StochasticDiffEq.JDiagonal_iip(W.dW)
  StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Jdiag, nothing, 1)
  Jdiagoop = StochasticDiffEq.JDiagonal_oop()

  @test StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, Jdiagoop, nothing, 1) == Jdiag.J
  @test Jdiag.J == true_diag

  for alg ∈ LevyArea.ITER_INT_ALGS
    @test diag(StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, alg, nothing, 1)) == true_diag
  end
end

# Commutative noise tests
@testset "Commutative noise tests" begin
  true_commute = 1//2 .* W.dW .* W.dW'

  Jcommute = StochasticDiffEq.JCommute_iip(W.dW)
  StochasticDiffEq.get_iterated_I!(dt, W.dW, W.dZ, Jcommute, nothing, 1)
  Jcommuteoop = StochasticDiffEq.JCommute_oop()

  @test StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, Jcommuteoop, nothing, 1) == Jcommute.J
  @test Jcommute.J == true_commute

  for alg ∈ LevyArea.ITER_INT_ALGS
    @test diag(StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, alg, nothing, 1)) == diag(true_commute)
  end
end

# moment conditions
# E1(I_{j1} I_{j2}) = Δ δ_{j1,j2}
# E2(I_{j1, j2} I_{j1, j2}) = 1/2 Δ^2
# E3(I_{j1} I_{j2} I_{j3, j4}) = {Δ^2 if j1=j2=j3=j4, 1/2 Δ^2 if j3!=j4 and j1=j3, j2=j4 or j1=j4, j2=j3, 0 otherwise}
function test_moments(m, alg, Δ, samples, p=nothing)
  # generate new random dW
  W = WienerProcess!(0.0,zeros(m),nothing)
  E1 = false .* vec(W.dW) .* vec(W.dW)'
  E2 = zero(E1)
  E3 = zero(E1)
  tmp = zero(E1)
  for _ in 1:samples
    calculate_step!(W,Δ,nothing,nothing)
    accept_step!(W,Δ,nothing,nothing)
    mul!(tmp,vec(W.dW),vec(W.dW)')

    I = StochasticDiffEq.get_iterated_I(Δ, W.dW, W.dZ, alg, p, 1)
    @. E1 = E1 + tmp
    @. E2 = E2 + I * I
    @. E3 = E3 + tmp * I
  end
  @. E1 = E1/samples
  @. E2 = E2/samples
  @. E3 = E3/samples

  @test maximum(abs.(E1 - Δ*I)) < 2e-1
  @test maximum(abs.(E2- Δ^2*(zero(E2).+1)/2)) < 1e-2
  @test maximum(abs.(E3 - Δ^2*(one(E2).+1)/2)) < 1e-2
end

"""
Exercise 2.3.1 from
Kloeden, P. E., Platen, E., & Schurz, H. Numerical solution of SDE through computer
experiments. Springer Science & Business Media. (2012)
"""

function test_compare_sample_mean_and_var(alg, Δ, m, samples=Int(1e6), p=Int(1e2))
  W = WienerProcess(0.0,zeros(m),nothing)
  calculate_step!(W,dt,nothing,nothing)
  for i in 1:10
    accept_step!(W,dt,nothing,nothing)
  end

  xs = [StochasticDiffEq.get_iterated_I(Δ, W.dW, W.dZ, alg, p, 1)[1,2]]
  coms = [1//2*W.dW[1]*W.dW[2]]

  for i in 1:samples
    calculate_step!(W,Δ,nothing,nothing)
    accept_step!(W,Δ,nothing,nothing)
    com = 1//2*W.dW[1]*W.dW[2]
    x = StochasticDiffEq.get_iterated_I(Δ, W.dW, W.dZ, alg, p, 1)[1,2]

    push!(xs,x)
    push!(coms,com)
  end

  @test mean(xs) ≈ 0 atol = 1e-2
  @test mean(coms) ≈ 0 atol = 1e-2
  @test var(xs) ≈ 1//2*Δ^2 rtol = 1e-2
  @test var(coms) ≈ 1//4*Δ^2 rtol = 1e-2
end

@testset "General noise tests" begin
  true_commute = 1//2 .* W.dW .* W.dW'
  samples = Int(1e4)
  for alg ∈ LevyArea.ITER_INT_ALGS
    # Test the relations given in Wiktorsson paper Eq.(2.1)
    Random.seed!(seed)
    I = StochasticDiffEq.get_iterated_I(dt, W.dW, W.dZ, alg, 1, 1)
    Random.seed!(seed)
    A = LevyArea.levyarea(W.dW/√dt, 1, alg)
    @test dt*A + true_commute == I # because of correction with \mu term
    @test diag(A) == diag(zero(A))
    @test diag(true_commute) == diag(I)
    @test I + I' ≈ 2*true_commute atol=1e-12
    @test A ≈ -A' atol=1e-12

    # test moment conditions
    test_moments(m, alg, dt, samples)

    # test other StatsJ
    @time test_compare_sample_mean_and_var(alg, 1.0, 2)
  end
end

@testset "Simulate non-com. SDEs tests" begin
  u₀ = [0.0,0.0,0.0]
  f(u,p,t) = [0.0,0.0,0.0]
  g(u,p,t) = [1.0 0.0;0.0 1.0;0.0 u[1]]
  f!(du,u,p,t) = du .*= false
  function g!(du,u,p,t)
    du[1,1] = 1.0
    du[2,2] = 1.0
    du[3,2] = u[1]
    return nothing
  end
  dt = 1//2^8
  tspan = (0.0,1.0)
  prob = SDEProblem(f,g,u₀,(0.0,1.0); noise_rate_prototype=zeros(3,2))
  prob! = SDEProblem(f!,g!,u₀,(0.0,1.0); noise_rate_prototype=zeros(3,2))

  sol = solve(prob, RKMilGeneral(;ii_approx=IICommutative()), dt=dt, adaptive=false)
  @test sol.retcode == :Success
  sol = solve(prob, RKMilGeneral(;ii_approx=IILevyArea()), dt=dt, adaptive=false)
  @test sol.retcode == :Success
  sol = solve(prob, RKMilGeneral(ii_approx=Fourier()), dt=dt, adaptive=false)
  @test sol.retcode == :Success

  sol = solve(prob!, RKMilGeneral(;ii_approx=IICommutative()), dt=dt, adaptive=false)
  @test sol.retcode == :Success
  sol = solve(prob!, RKMilGeneral(;ii_approx=IILevyArea()), dt=dt, adaptive=false)
  @test sol.retcode == :Success
  sol = solve(prob!, RKMilGeneral(ii_approx=Fourier()), dt=dt, adaptive=false)
  @test sol.retcode == :Success
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq

max_issues_repo_name: JuliaDiffEq/StochasticDiffEq

max_forks_repo_name: JuliaDiffEq/StochasticDiffEq


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Sudoku

makedocs(
    sitename = "Sudoku",
    format = Documenter.HTML(),
    modules = [Sudoku]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
srand(1)

context("migration.jl") do
  # We use complex landscapes here to ensure there are no accidental
  # duplicates.
  landscapes = [
    NK.NKLandscape(96, 1),
    NK.NKqLandscape(96, 1, 2),
    NK.NKpLandscape(96, 1, 0.90),
  ]

  function testmigrate(func::Function, ls::NK.Landscape)
    p = rand(NK.MetaPopulation, ls, 5, 10)
    np = func(p, 1.0, 5)
    for popind = 1:NK.popct(p)
      origpop = p.populations[popind]
      newpop = np.populations[popind]
      sharect = 0
      for genoind = 1:NK.popsize(p)
        origgeno = origpop.genotypes[genoind]
        newgeno = newpop.genotypes[genoind]
        if origgeno == newgeno
          sharect += 1
        end
      end
      @fact sharect --> less_than(NK.popsize(p)) "expected migrations"
    end
  end

  function testmigprob(func::Function, ls::NK.Landscape)
    numpops = 10000
    p = rand(NK.MetaPopulation, ls, 1, numpops)
    np = func(p, 0.25, 1)
    sharect = 0
    for popind = 1:NK.popct(p)
      origgeno = p.populations[popind].genotypes[1]
      newgeno = np.populations[popind].genotypes[1]
      if origgeno == newgeno
        sharect += 1
      end
    end
    sharepct = sharect / numpops
    @fact sharepct --> roughly(0.75, 0.15)
  end

  function testmigct(func::Function, ls::NK.Landscape)
    p = rand(NK.MetaPopulation, ls, 2, 3)
    np = func(p, 1.0, 1)
    for popind = 1:NK.popct(p)
      origpop = p.populations[popind]
      newpop = np.populations[popind]
      sharect = sum(origpop.genotypes .== newpop.genotypes)
      @fact sharect --> 1
    end
  end

  for ls = landscapes

    context("$(ls)") do
      context("migrate(...)") do
        testmigrate(NK.migrate, ls)

        context("should respect migration probability") do
          testmigprob(NK.migrate, ls)
        end

        context("should respect migration count") do
          testmigct(NK.migrate, ls)
        end
      end

      context("linmigrate(...)") do
        testmigrate(NK.linmigrate, ls)

        context("should respect migration probability") do
          testmigprob(NK.linmigrate, ls)
        end

        context("should respect migration count") do
          testmigct(NK.linmigrate, ls)
        end

        context("should only migrate left and right") do
          p = rand(NK.MetaPopulation, ls, 1, 3)
          np = NK.linmigrate(p, 1.0, 1)
          for popind = 1:NK.popct(p)
            newpop = np.populations[popind]
            leftorig = if popind == 1
              p.populations[NK.popct(p)]
            else
              p.populations[popind - 1]
            end
            rightorig = if popind == NK.popct(p)
              p.populations[1]
            else
              p.populations[popind + 1]
            end
            newgeno = newpop.genotypes[1]
            origgenos = [leftorig.genotypes[1], rightorig.genotypes[1]]
            @fact newgeno --> anyof(origgenos...)
          end
        end
      end
    end
  end
end


text:
"""
    epsp(t::Real, lastspike, q, taum, taus)
    epsp(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, taum::AbstractArray{<:Real}, taus::AbstractArray{<:Real})
    epsp(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, taum::CuVecOrMat{<:Real}, taus::CuVecOrMat{<:Real})

Evaluate an EPSP synapse. Modeled as `(ϵ₀ / τm - τs) * (exp(-Δ / τm) - exp(-Δ / τs)) Θ(Δ)`
  (where `Θ` is the Heaviside function and `Δ = t - lastspike`).
Specifically, this is the EPSP time course for the SRM0 model introduced by Gerstner.
Details: [Spiking Neuron Models: Single Neurons, Populations, Plasticity]
         (https://icwww.epfl.ch/~gerstner/SPNM/node27.html#SECTION02323400000000000000)

Use `CuVector` instead of `Vector` for GPU support.

# Fields
- `t`: current time
- `lastspike`: last pre-synaptic spike time
- `q`: amplitude
- `taum`: rise time constant
- `taus`: fall time constant
"""
function epsp(t::Real, lastspike, q, taum, taus)
    Δ = t - lastspike

    return (Δ >= 0 && Δ < Inf && taus != taum) * q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))
end
function epsp(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}, taum::AbstractArray{<:Real}, taus::AbstractArray{<:Real})
    Δ = t .- lastspike
    I = @. q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))

    return map((ts, tm, δ, i) -> (δ >= 0) && (δ < Inf) && (ts != tm) ? i : zero(i), taus, taum, Δ, I)
end
function epsp(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}, taum::CuVecOrMat{<:Real}, taus::CuVecOrMat{<:Real})
    Δ = t .- lastspike

    return @. (Δ >= 0) * (Δ < Inf) * (taus != taum) * q / (1 - taus / taum) * (exp(-Δ / taum) - exp(-Δ / taus))
end


text:
const CPetscOptions = Ptr{Cvoid}

"""
    AbstractOptions{PetscLib <: PetscLibType}

Abstract type of PETSc solver options.
"""
abstract type AbstractOptions{PetscLib <: PetscLibType} end

"""
    GlobalOptions{PetscLib <: PetscLibType}

The PETSc global options database.
"""
struct GlobalOptions{PetscLib} <: AbstractOptions{PetscLib} end
Base.cconvert(::Type{CPetscOptions}, obj::GlobalOptions) = C_NULL

"""
    Options{PetscLib <: PetscLibType}(kw -> arg, ...)
    Options(petsclib, kw -> arg, ...)

Create a PETSc options data structure for the `petsclib`.

For construction a set of keyword argment pairs should be given. If the option
has no value it should be set to `nothing` or `true`. Setting an option to
`false` will cause the option not to be set on the PETSc options table.

# Examples
```julia-repl
julia> using PETSc

julia> petsclib = PETSc.petsclibs[1];

julia> PETSc.initialize(petsclib)

julia> opt = PETSc.Options(
                         petsclib,
                         ksp_monitor = nothing,
                         ksp_view = true,
                         pc_type = "mg",
                         pc_mg_levels = 1,
                         false_opt = false,
                     )
#PETSc Option Table entries:
-ksp_monitor
-ksp_view
-pc_mg_levels 1
-pc_type mg
#End of PETSc Option Table entries

julia> opt["ksp_monitor"]
""

julia> opt["pc_type"]
"mg"

julia> opt["pc_type"] = "ilu"
"ilu"

julia> opt["pc_type"]
"ilu"

julia> opt["false_opt"]
ERROR: KeyError: key "bad_key" not found

julia> opt["bad_key"]
ERROR: KeyError: key "bad_key" not found
```

# External Links
$(_doc_external("Sys/PetscOptionsCreate"))
"""
mutable struct Options{T} <: AbstractOptions{T}
    ptr::CPetscOptions
end

Options(petsclib; kwargs...) = Options(petsclib, kwargs...)
function Options(petsclib, ps::Pair...)
    opts = Options(petsclib)
    for (k, v) in ps
        opts[k] = v
    end
    return opts
end

@for_petsc function Options(::$UnionPetscLib)
    opts = Options{$PetscLib}(C_NULL)
    @assert initialized($PetscLib)
    @chk ccall(
        (:PetscOptionsCreate, $petsc_library),
        PetscErrorCode,
        (Ptr{CPetscOptions},),
        opts,
    )
    finalizer(finalize, opts)
    return opts
end

@for_petsc function finalize(opts::Options{$PetscLib})
    finalized($PetscLib) || @chk ccall(
        (:PetscOptionsDestroy, $petsc_library),
        PetscErrorCode,
        (Ptr{CPetscOptions},),
        opts,
    )
    return nothing
end

@for_petsc function Base.push!(
    ::GlobalOptions{$PetscLib},
    opts::Options{$PetscLib},
)
    @chk ccall(
        (:PetscOptionsPush, $petsc_library),
        PetscErrorCode,
        (CPetscOptions,),
        opts,
    )
    return nothing
end

@for_petsc function Base.pop!(::GlobalOptions{$PetscLib})
    @chk ccall((:PetscOptionsPop, $petsc_library), PetscErrorCode, ())
    return nothing
end

@for_petsc function Base.setindex!(opts::AbstractOptions{$PetscLib}, val, key)
    val === true && (val = nothing)
    val === false && (return val)

    @chk ccall(
        (:PetscOptionsSetValue, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, Cstring, Cstring),
        opts,
        string('-', key),
        isnothing(val) ? C_NULL : string(val),
    )

    return val
end

@for_petsc function Base.getindex(opts::AbstractOptions{$PetscLib}, key)
    val = Vector{UInt8}(undef, 256)
    set_ref = Ref{PetscBool}()
    @chk ccall(
        (:PetscOptionsGetString, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, Cstring, Cstring, Ptr{UInt8}, Csize_t, Ptr{PetscBool}),
        opts,
        C_NULL,
        string('-', key),
        val,
        sizeof(val),
        set_ref,
    )
    val = GC.@preserve val unsafe_string(pointer(val))
    set_ref[] == PETSC_TRUE || throw(KeyError(key))
    return val
end

@for_petsc function view(
    opts::AbstractOptions{$PetscLib},
    viewer::AbstractViewer{$PetscLib} = ViewerStdout($PetscLib, MPI.COMM_SELF),
)
    @chk ccall(
        (:PetscOptionsView, $petsc_library),
        PetscErrorCode,
        (CPetscOptions, CPetscViewer),
        opts,
        viewer,
    )
    return nothing
end

@for_petsc GlobalOptions(::$UnionPetscLib) = GlobalOptions{$PetscLib}()

Base.show(io::IO, opts::AbstractOptions) = _show(io, opts)

"""
    with(f, opts::Options)

Call `f()` with the [`Options`](@ref) `opts` set temporarily (in addition to any
global options).
"""
function with(f, opts::Options{T}) where {T}
    global_opts = GlobalOptions{T}()
    push!(global_opts, opts)
    try
        f()
    finally
        pop!(global_opts)
    end
end

"""
    parse_options(args::Vector{String})

Parse the `args` vector into a `NamedTuple` that can be used as the options for
the PETSc solvers.

```sh
julia --project file.jl -ksp_monitor -pc_type mg -ksp_view
```
"""
function parse_options(args::Vector{String})
    i = 1
    opts = Dict{Symbol, Union{String, Nothing}}()
    while i <= length(args)
        @assert args[i][1] == '-' && length(args[i]) > 1
        if i == length(args) || args[i + 1][1] == '-'
            opts[Symbol(args[i][2:end])] = nothing
            i = i + 1
        else
            opts[Symbol(args[i][2:end])] = args[i + 1]
            i = i + 2
        end
    end
    return NamedTuple(opts)
end


text:
using ADCME
include("nnutil.jl")
nntype="ae_scaled"
X, Y = prepare_strain_stress_data1D(domain)
x = (constant(X[:,1]), constant(X[:,2]), constant(X[:,3]))
y = squeeze(nn(x...))
# error()
loss = sum((y-Y)^2)
variable_scope("nn") do
    global opt = AdamOptimizer().minimize(loss)
end

sess = Session(); init(sess)

# for i = 1:10000
#     l, _ = run(sess, [loss, opt])
#     @show i,l
# end
ADCME.load(sess, "Data/learned_nn.mat")
BFGS!(sess, loss, 5000)

out = run(sess, y)
close("all")
plot(X[:,2], out,"+", label="NN")
plot(X[:,2], Y, ".", label="Exact")
legend()

ADCME.save(sess, "Data/learned_nn.mat")

# for α = [0.01,0.05,0.1]
# x = constant(X.*(1 .+ α*rand(size(X)...)))
# y = squeeze(ae(x, [20,20,20,20,1], "nn"))
# out = run(sess, y)
# plot(X[:,2], out,".", label="NN")
# end
# # plot(X[:,2], Y, ".", label="Exact")
# legend()

@load "Data/domain.jld2" domain
X, Y = prepare_strain_stress_data1D(domain)
x = constant(X)
x = (constant(X[:,1]), constant(X[:,2]), constant(X[:,3]))
y = squeeze(nn(x...))
sess = Session(); init(sess)
close("all")
ADCME.load(sess, "Data/learned_nn.mat")
out = run(sess, y)
plot(X[:,1], out,"+", label="NN")
plot(X[:,1], Y, ".", label="Exact")
legend()


text:
using DataFrames, MCMCChains, CSV

"""

# read_summary

Read summary output file created by stansummary.

### Method
```julia
read_summary(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_summary(m::Stanmodel)

  fname = "$(m.tmpdir)/$(m.name)_summary.csv"
  !isfile(fname) && stan_summary(m)

  df = CSV.read(fname, delim=",", comment="#")

  cnames = lowercase.(convert.(String, String.(names(df))))
  cnames[1] = "parameters"
  cnames[4] = "std"
  cnames[8] = "ess"
  names!(df, Symbol.(cnames))
  df.parameters = Symbol.(df.parameters)

  ChainDataFrame("CmdStan Summary", df)

end   # end of read_samples


text:
type ELSQVariables{T <: Float, V <: DenseVector}
    b        :: V
    nrmsq    :: Vector{T}
    df       :: V
    dotprods :: V
    tempp    :: V
    r        :: V
    tempn    :: V
    tempn2   :: V
    xb       :: V
    perm     :: Vector{Int}
    inner    :: Dict{Int, Vector{T}}
    mask_n   :: Vector{Int}
    idx      :: BitArray{1}

    ELSQVariables(b::DenseVector{T}, nrmsq::DenseVector{T}, df::DenseVector{T}, dotprods::DenseVector{T}, tempp::DenseVector{T}, r::DenseVector{T}, tempn::DenseVector{T}, tempn2::DenseVector{T}, xb::DenseVector{T}, perm::Vector{Int}, inner::Dict{Int, Vector{T}}, mask_n::Vector{Int}, idx::BitArray{1}) = new(b, nrmsq, df, dotprods, tempp, r, tempn, tempn2, xb, perm, inner, mask_n, idx)
end

function ELSQVariables{T <: Float}(
    b        :: DenseVector{T},
    nrmsq    :: Vector{T},
    df       :: DenseVector{T},
    dotprods :: DenseVector{T},
    tempp    :: DenseVector{T},
    r        :: DenseVector{T},
    tempn    :: DenseVector{T},
    tempn2   :: DenseVector{T},
    xb       :: DenseVector{T},
    perm     :: Vector{Int},
    inner    :: Dict{Int, Vector{T}},
    mask_n   :: Vector{Int},
    idx      :: BitArray{1}
)
    ELSQVariables{T, typeof(b)}(b, nrmsq, df, dotprods, tempp, r, tempn, tempn2, xb, perm, inner, mask_n, idx)
end

function ELSQVariables{T <: Float}(
    x :: DenseMatrix{T},
    y :: DenseVector{T}
)
    # dimensions of arrays
    n,p = size(x)

    # form arrays
    b        = zeros(T, p)
    nrmsq    = vec(sumabs2(x,1)) :: Vector{T}
    df       = zeros(T, p)
    dotprods = zeros(T, p)
    tempp    = zeros(T, p)
    r        = zeros(T, n)
    tempn    = zeros(T, n)
    tempn2   = zeros(T, n)
    xb       = zeros(T, n)
    perm     = collect(1:p)
    mask_n   = ones(Int, n)
    idx      = falses(p)

    # form dictionary
    inner = Dict{Int, Vector{T}}()

    # return container object
    ELSQVariables{eltype(y), typeof(y)}(b, nrmsq, df, dotprods, tempp, r, tempn, tempn2, xb, perm, inner, mask_n, idx)
end

function ELSQVariables{T <: Float}(
    x :: BEDFile{T},
    y :: SharedVector{T},
)
    # dimensions of arrays
    n,p = size(x)

    # process ids?
    pids = procs(x)

    # form arrays
    b        = SharedArray(T, (p,), pids=pids) :: typeof(y)
    nrmsq    = (length(y) - 1) * ones(T, p)
    df       = SharedArray(T, (p,), pids=pids) :: typeof(y)
    dotprods = SharedArray(T, (p,), pids=pids) :: typeof(y)
    tempp    = SharedArray(T, (p,), pids=pids) :: typeof(y)
    r        = SharedArray(T, (n,), pids=pids) :: typeof(y)
    tempn    = SharedArray(T, (n,), pids=pids) :: typeof(y)
    tempn2   = SharedArray(T, (n,), pids=pids) :: typeof(y)
    xb       = SharedArray(T, (n,), pids=pids) :: typeof(y)
    perm     = collect(1:p)
    mask_n   = ones(Int, n)
    idx      = falses(p)

    # form dictionary
    inner = Dict{Int, Vector{T}}()

    # return container object
    ELSQVariables{eltype(y), typeof(y)}(b, nrmsq, df, dotprods, tempp, r, tempn, tempn2, xb, perm, inner, mask_n, idx)
end

function ELSQVariables{T <: Float}(
    x :: BEDFile{T},
    y :: SharedVector{T},
    z :: DenseVector{Int} # <-- this should be the bitmask
)
    w = ELSQVariables(x,y)
    copy!(w.mask_n, z)
    return w
end

"A subroutine to automatically set the default number of folds, if the user has not provided them."
set_cv_fold_num(minfolds::Int, maxfolds::Int) = max(minfolds, min(Sys.CPU_CORES, maxfolds)) :: Int

immutable ELSQCrossvalidationResults{T <: Float}
    mses :: Vector{T}
    b    :: Vector{T}
    bidx :: Vector{Int}
    k    :: Int
    path :: Vector{Int}
    bids :: Vector{String}
end

# constructor for when bids are not available
# simply makes vector of "V$i" where $i are drawn from bidx
function ELSQCrossvalidationResults{T <: Float}(
    mses :: Vector{T},
    b    :: Vector{T},
    bidx :: Vector{Int},
    k    :: Int,
    path :: Vector{Int},
)
    bids = convert(Vector{String}, ["V" * "$i" for i in bidx]) :: Vector{String}
    ELSQCrossvalidationResults{eltype(mses)}(mses, b, bidx, k, path, bids)
end

# function to view an ELSQCrossvalidationResults object
function Base.show(io::IO, x::ELSQCrossvalidationResults)
    println(io, "Crossvalidation results:")
    println(io, "Minimum MSE ", minimum(x.mses), " occurs at k = $(x.k).")
    println(io, "Best model β has the following nonzero coefficients:")
    println(io, DataFrame(Predictor=x.bidx, Name=x.bids, Estimated_β=x.b))
    return nothing
end

function Gadfly.plot(x::ELSQCrossvalidationResults)
    df = DataFrame(ModelSize=x.path, MSE=x.mses)
    plot(df, x="ModelSize", y="MSE", xintercept=[x.k], Geom.line, Geom.vline(color=colorant"red"), Guide.xlabel("Model size"), Guide.ylabel("MSE"), Guide.title("MSE versus model size"))
end

function check_finiteness{T <: Float}(x::T)
    isnan(x) && throw(error("Objective function is NaN, aborting..."))
    isinf(x) && throw(error("Objective function is Inf, aborting..."))
end

function print_descent_error{T <: Float}(iter::Int, loss::T, next_loss::T)
    print_with_color(:red, "\nExchange algorithm fails to descend!\n")
    print_with_color(:red, "Iteration: $(iter)\n")
    print_with_color(:red, "Current Objective: $(loss)\n")
    print_with_color(:red, "Next Objective: $(next_loss)\n")
    print_with_color(:red, "Difference in objectives: $(abs(next_loss - loss))\n")
    throw(error("Descent failure!"))
end

function print_maxiter{T <: Float}(max_iter::Int, loss::T)
    print_with_color(:red, "Exchange algorithm has hit maximum iterations $(max_iter)!\n")
    print_with_color(:red, "Return value may be incorrect\n")
    print_with_color(:red, "Current Loss: $(loss)\n")
end

function errorcheck{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    tol      :: T,
    max_iter :: Int,
    window   :: Int,
    p        :: Int = size(x,2)
)
    tol >= eps()          || throw(ArgumentError("Global tolerance must exceed machine precision"))
    max_iter >= 1         || throw(ArgumentError("Maximum number of iterations must exceed 1"))
    return nothing
end

function errorcheck{T <: Float}(
    x        :: BEDFile{T},
    y        :: SharedVector{T},
    tol      :: T,
    max_iter :: Int,
    window   :: Int,
    p        :: Int = size(x,2)
)
    tol >= eps()          || throw(ArgumentError("Global tolerance must exceed machine precision"))
    max_iter >= 1         || throw(ArgumentError("Maximum number of iterations must exceed 1"))
    return nothing
end

function print_cv_results{T <: Float}(errors::DenseVector{T}, path::DenseVector{Int}, k::Int)
    println("\n\nCrossvalidation Results:")
    println("k\tMSE")
    for i = 1:length(errors)
        println(path[i], "\t", errors[i])
    end
    println("\nThe lowest MSE is achieved at k = ", k)
end

#########################################
### subroutines for exchange_leastsq! ###
#########################################

# subroutine compares current predictor i against all predictors k+1, k+2, ..., p
# these predictors are candidates for inclusion in set
# _exlstsq_innerloop! find best new predictor
function _exlstsq_innerloop!{T <: Float}(
    v   :: ELSQVariables{T},
    k   :: Int,
    i   :: Int,
    p   :: Int,
    tol :: T
)
    l     = v.perm[i] :: Int
    betal = v.b[l]    :: T

    # save values to determine best estimate for current predictor
    b   = v.nrmsq[l] :: T
    a   = (v.df[l] + betal*b) :: T
    adb = a / b :: T
    r   = i

    # inner loop compares current predictor j against all remaining predictors j+1,...,p
    for j = (k+1):p
        idx = v.perm[j] :: Int
        c   = (v.df[idx] + betal*v.dotprods[idx]) :: T
        d   = v.nrmsq[idx] :: T

        # if current inactive predictor beats current active predictor,
        # then save info for swapping
        if c*c/d > a*adb + tol
            a   = c
            b   = d
            r   = j
            adb = a / b
        end
    end # end inner loop over remaining predictor set

    return a, b, r, adb
end

# subroutine to swap predictors
# in exchange_leastsq!, will swap best predictor with current predictor
function _swap_predictors!{T <: Float}(
    v   :: ELSQVariables{T},
    i   :: Int,
    r   :: Int,
    m   :: Int,
    adb :: T
)
    # save ith best predictor
    j = v.perm[i]

    # replace index of ith best with new best predictor
    v.perm[i] = v.perm[r]

    # replace new best predictor with ith best
    # at this point, swap of indices is complete
    v.perm[r] = j

    # replace coefficient of mth best predictor with a / b
    v.b[m] = adb

    # if rth and ith predictors coincide,
    # then set ith best predictor coefficient to zero
    if r != i
        v.b[j] = zero(T)
    end
    return nothing
end

"""
    axpymbz!(y,a,x,b,z[, p=length(y)])

The silly name is based on BLAS `axpy()` (A*X Plus Y), except that this function performs *A**X* *P*lus *Y* *M*inus *B**Z*.
The idea behind `axpymz!()` is to perform the computation in one pass over the arrays. The output is the same as `y = y + a*x - b*z`.
"""
function axpymbz!{T <: Float}(
    y :: DenseVector{T},
    a :: T,
    x :: DenseVector{T},
    b :: T,
    z :: DenseVector{T};
)
    @inbounds @simd for i in eachindex(y)
        y[i] = y[i] + a*x[i] - b*z[i]
    end
    return nothing
end

# subroutine to refit preditors after crossvalidation
function refit_exlstsq{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    k        :: Int;
    models   :: DenseVector{Int} = collect(1:min(20,size(x,2))),
    tol      :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    quiet    :: Bool = true,
)
    # initialize β vector and temporary arrays
    v = ELSQVariables(x, y)

    # first use exchange algorithm to extract model
    exchange_leastsq!(v, x, y, k, max_iter=max_iter, quiet=quiet, tol=tol, window=k)

    # which components of β are nonzero?
    # cannot use binary indices here since we need to return Int indices
    bidx = find(v.b)

    # allocate the submatrix of x corresponding to the inferred model
    # cannot use SubArray since result is not StridedArray?
    # issue is that bidx is Vector{Int} and not a Range object
    # use of SubArray is more memory efficient; a pity that it doesn't work!
#    x_inferred = sub(sdata(x), :, bidx)
    x_inferred = x[:,bidx]

    # now estimate β with the ordinary least squares estimator β = inv(x'x)x'y
    # return it with the vector of MSEs
    xty = BLAS.gemv('T', x_inferred, sdata(y)) :: Vector{T}
    xtx = BLAS.gemm('T', 'N', x_inferred, x_inferred) :: Matrix{T}
    b   = zeros(T, length(bidx))
    try
        b = (xtx \ xty) :: Vector{T}
    catch e
        warn("caught error: ", e, "\nSetting returned values of b to -Inf")
        fill!(b, -Inf)
    end

    return b, bidx
end

# refitting routine for GWAS data with x', mean, prec files
function refit_exlstsq(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String,
    k        :: Int;
    models   :: DenseVector{Int} = collect(1:20),
    pids     :: DenseVector{Int} = procs(),
    tol      :: Float = convert(T, 1e-6),
    max_iter :: Int   = 100,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)

    # initialize all variables
    x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, pids=pids, header=header)
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}
    v = ELSQVariables(x, y, ones(Int, length(y)))

    # first use exchange algorithm to extract model
    exchange_leastsq!(v, x, y, k, max_iter=max_iter, quiet=quiet, tol=tol, window=k)

    # which components of β are nonzero?
    inferred_model = v.b .!= zero(T)
    bidx = find(inferred_model)

    # allocate the submatrix of x corresponding to the inferred model
    x_inferred = zeros(T, x.geno.n, sum(inferred_model))
    decompress_genotypes!(x_inferred, x, inferred_model)

    # now estimate b with the ordinary least squares estimator b = inv(x'x)x'y
    # return it with the vector of MSEs
    xty = BLAS.gemv('T', x_inferred, sdata(y)) :: Vector{T}
    xtx = BLAS.gemm('T', 'N', x_inferred, x_inferred) :: Matrix{T}
    b   = zeros(T, length(bidx))
    try
        b = (xtx \ xty) :: Vector{T}
    catch e
        warn("caught error: ", e, "\nSetting returned values of b to -Inf")
        fill!(b, -Inf)
    end

    # get predictor names before returning
    bids = prednames(x)[bidx]

    return b, bidx
end

# refitting routine for GWAS data with just genotypes, covariates, y
function refit_exlstsq(
    T        :: Type,
    xfile    :: String,
    x2file   :: String,
    yfile    :: String,
    k        :: Int;
    models   :: DenseVector{Int} = collect(1:20),
    pids     :: DenseVector{Int} = procs(),
    tol      :: Float = convert(T, 1e-6),
    max_iter :: Int   = 100,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)

    # initialize all variables
    x = BEDFile(T, xfile, x2file, pids=pids, header=header)
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}
    v = ELSQVariables(x, y, ones(Int, length(y)))

    # first use exchange algorithm to extract model
    exchange_leastsq!(v, x, y, k, max_iter=max_iter, quiet=quiet, tol=tol, window=k)

    # which components of β are nonzero?
    inferred_model = v.b .!= zero(T)
    bidx = find(inferred_model)

    # allocate the submatrix of x corresponding to the inferred model
    x_inferred = zeros(T, x.geno.n, sum(inferred_model))
    decompress_genotypes!(x_inferred, x, inferred_model)

    # now estimate β with the ordinary least squares estimator b = inv(x'x)x'y
    # return it with the vector of MSEs
    xty = BLAS.gemv('T', x_inferred, sdata(y)) :: Vector{T}
    xtx = BLAS.gemm('T', 'N', x_inferred, x_inferred) :: Matrix{T}
    b   = zeros(T, length(bidx))
    try
        b = (xtx \ xty) :: Vector{T}
    catch e
        warn("caught error: ", e, "\nSetting returned values of b to -Inf")
        fill!(b, -Inf)
    end

    # get predictor names before returning
    bids = prednames(x)[bidx]

    return b, bidx, bids
end

function update_current_best_predictor!{T <: Float}(
    v     :: ELSQVariables{T},
    x     :: DenseMatrix{T},
    betal :: T,
    adb   :: T,
    r     :: Int
)
    # first get index of current best predictor
    m = v.perm[r] :: Int

    # tempn2 = x[:,m]
    copy!(v.tempn2, view(x, :, m))

    # v.r = betal*v.tempn + adb*v.tempn2
    axpymbz!(v.r, betal, v.tempn, adb, v.tempn2)

    return m
end

function update_current_best_predictor!{T <: Float}(
    v     :: ELSQVariables{T},
    x     :: BEDFile{T},
    betal :: T,
    adb   :: T,
    r     :: Int
)
    # first get index of current best predictor
    m = v.perm[r] :: Int

    # v.tempn2 = x[:,m]
    decompress_genotypes!(v.tempn2, x, m)

    # v.tempn2[v.mask_n .== 0] = 0.0
    mask!(v.tempn2, v.mask_n, 0, zero(T))

    # v.r = betal*v.tempn + adb*v.tempn2
    axpymbz!(v.r, betal, v.tempn, adb, v.tempn2)

    # v.r[v.mask_n .== 0] = 0.0
    mask!(v.r, v.mask_n, 0, zero(T))

    return m
end

function get_inner_product!{T <: Float}(
    z :: DenseVector{T},
    w :: DenseVector{T},
    v :: ELSQVariables{T},
    x :: DenseMatrix{T},
    i :: Int
)
    if !haskey(v.inner, i)
        v.inner[i] = BLAS.gemv('T', one(T), x, w)
    end
    copy!(z, v.inner[i])
end

function get_inner_product!{T <: Float}(
    z :: DenseVector{T},
    w :: DenseVector{T},
    v :: ELSQVariables{T},
    x :: BEDFile{T},
    i :: Int;
    pids :: DenseVector{Int} = procs(x)
)
    if !haskey(v.inner, i)
        At_mul_B!(z, x, w, v.mask_n, pids=pids)
        v.inner[i] = copy(z)
    end
    copy!(z, v.inner[i])
end

function get_inner_product!{T <: Float}(
    z :: DenseVector{T},
    w :: DenseVector{T},
    v :: ELSQVariables{T},
    x :: BEDFile{T},
    a :: PlinkGPUVariables{T},
    i :: Int;
    pids :: DenseVector{Int} = procs(x)
)
    if !haskey(v.inner, i)
        At_mul_B!(z, x, w, v.mask_n, a)
        v.inner[i] = copy(z)
    end
    copy!(z, v.inner[i])
end


text:
module FourierOperations

export divergence, curl!, sqnorm

using PencilFFTs.PencilArrays
using Reexport

include("Grids.jl")
@reexport using .Grids

const VectorField{T} = NTuple{N, PencilArray{T,3}} where {N}

"""
    divergence(u::AbstractArray{<:Complex}, grid::FourierGrid)

Compute total divergence ``∑|∇⋅u|²`` in Fourier space, in the local process.
"""
function divergence(uF_local::VectorField{T},
                    grid::FourierGrid) where {T <: Complex}
    div2 = real(zero(T))

    uF = map(global_view, uF_local)
    ux = first(uF)

    @inbounds for (l, I) in enumerate(CartesianIndices(ux))
        k = grid[I]  # (kx, ky, kz)
        div = zero(T)
        for n in eachindex(k)
            v = 1im * k[n] * uF[n][l]
            div += v
        end
        div2 += abs2(div)
    end

    div2
end

# Local grid variant (faster -- with linear indexing!)
function divergence(uF::VectorField{T},
                    grid::FourierGridIterator) where {T <: Complex}
    div2 = real(zero(T))
    @inbounds for (i, k) in enumerate(grid)
        div = zero(T)
        for n in eachindex(k)
            v = 1im * k[n] * uF[n][i]
            div += v
        end
        div2 += abs2(div)
    end
    div2
end

"""
    curl!(ω, u, grid::FourierGrid)

Compute ``ω = ∇×u`` in Fourier space.
"""
function curl!(ωF_local::VectorField{T},
               uF_local::VectorField{T},
               grid::FourierGrid) where {T <: Complex}
    u = map(global_view, uF_local)
    ω = map(global_view, ωF_local)

    @inbounds for I in CartesianIndices(u[1])
        k = grid[I]  # (kx, ky, kz)
        l = LinearIndices(u[1])[I]
        v = (u[1][l], u[2][l], u[3][l])
        ω[1][l] = 1im * (k[2] * v[3] - k[3] * v[2])
        ω[2][l] = 1im * (k[3] * v[1] - k[1] * v[3])
        ω[3][l] = 1im * (k[1] * v[2] - k[2] * v[1])
    end

    ωF_local
end

function curl!(ω::VectorField{T},
               u::VectorField{T},
               grid::FourierGridIterator) where {T <: Complex}
    @inbounds for (i, k) in enumerate(grid)
        v = (u[1][i], u[2][i], u[3][i])
        ω[1][i] = 1im * (k[2] * v[3] - k[3] * v[2])
        ω[2][i] = 1im * (k[3] * v[1] - k[1] * v[3])
        ω[3][i] = 1im * (k[1] * v[2] - k[2] * v[1])
    end
    ω
end

"""
    index_r2c(u::PencilArray)

Return index associated to dimension of real-to-complex transform.

This is assumed to be the first *logical* dimension of the array.
Since indices in the array may be permuted, the actual dimension may be other
than the first.
"""
index_r2c(u::PencilArray) = index_r2c(permutation(u))
index_r2c(::Nothing) = 1
index_r2c(::Val{p}) where {p} = findfirst(==(1), p) :: Int

"""
    sqnorm(u::AbstractArray{<:Complex}, grid::FourierGridIterator)

Compute squared norm of array in Fourier space, in the local process.
"""
sqnorm(u::PencilArray, grid) = sqnorm((u, ), grid)

function sqnorm(u::VectorField{T}, grid::FourierGridIterator) where {T <: Complex}
    gp = parent(grid) :: FourierGrid
    kx = gp[1]  # global wave numbers along r2c dimension

    # Note: when Nx (size of real input data) is even, the Nyquist frequency is
    # also counted twice.
    Nx = size(gp, 1)
    @assert length(kx) == div(Nx, 2) + 1

    k_zero = kx[1]  # zero mode

    kx_lims = if iseven(Nx)
        (k_zero, kx[end])  # kx = 0 or Nx/2 (Nyquist frequency)
    else
        # We repeat k_zero for type inference reasons.
        (k_zero, k_zero)  # only kx = 0
    end

    s = zero(real(T))

    @inbounds for (i, k) in enumerate(grid)
        # Account for Hermitian symmetry implied by r2c transform along the
        # first logical dimension. Note that `k` is "unpermuted", meaning that
        # k[1] is the first *logical* wave number.
        factor = k[1] in kx_lims ? 1 : 2
        s += factor * sum(v -> abs2(v[i]), u)
    end

    s
end

# Add a variant for real arrays, for completeness.
sqnorm(u::AbstractArray{T} where {T <: Real}) = sum(abs2, u)
sqnorm(u::PencilArray{T} where {T <: Real}) = sqnorm(parent(u))

sqnorm(u::Tuple, args...) = mapreduce(v -> sqnorm(v, args...), +, u)

end


text:
struct Dynamics{T}
    evaluate
    jacobian_state
    jacobian_action
    num_next_state::Int
    num_state::Int
    num_action::Int
    num_parameter::Int
    evaluate_cache::Vector{T}
    jacobian_state_cache::Matrix{T}
    jacobian_action_cache::Matrix{T}
end

Model{T} = Vector{Dynamics{T}} where T

function Dynamics(f::Function, num_state::Int, num_action::Int; num_parameter::Int=0)
    #TODO: option to load/save methods
    @variables x[1:num_state], u[1:num_action], w[1:num_parameter]
    y = f(x, u, w)
    jacobian_state = Symbolics.jacobian(y, x);
    jacobian_action = Symbolics.jacobian(y, u);
    evaluate_func = eval(Symbolics.build_function(y, x, u, w)[2]);
    jacobian_state_func = eval(Symbolics.build_function(jacobian_state, x, u, w)[2]);
    jacobian_action_func = eval(Symbolics.build_function(jacobian_action, x, u, w)[2]);
    num_next_state = length(y)

    return Dynamics(evaluate_func, jacobian_state_func, jacobian_action_func,
                    num_next_state, num_state, num_action, num_parameter,
                    zeros(num_next_state), zeros(num_next_state, num_state), zeros(num_next_state, num_action))
end

function dynamics!(d::Dynamics, state, action, parameter)
    d.evaluate(d.evaluate_cache, state, action, parameter)
    return d.evaluate_cache
end

function jacobian!(jacobian_states, jacobian_actions, dynamics::Vector{Dynamics{T}}, states, actions, parameters) where T
    for (t, d) in enumerate(dynamics)
        d.jacobian_state(d.jacobian_state_cache, states[t], actions[t], parameters[t])
        d.jacobian_action(d.jacobian_action_cache, states[t], actions[t], parameters[t])
        @views jacobian_states[t] .= d.jacobian_state_cache
        @views jacobian_actions[t] .= d.jacobian_action_cache
        fill!(d.jacobian_state_cache, 0.0) # TODO: confirm this is necessary
        fill!(d.jacobian_action_cache, 0.0) # TODO: confirm this is necessary
    end
end

num_trajectory(dynamics::Vector{Dynamics{T}}) where T = sum([d.num_state + d.num_action for d in dynamics]) + dynamics[end].num_next_state

# user-provided dynamics and gradients
function Dynamics(f::Function, fx::Function, fu::Function, num_next_state::Int, num_state::Int, num_action::Int,
    num_parameter::Int=0)
    return Dynamics(f, fx, fu,
                    num_next_state, num_state, num_action, num_parameter,
                    zeros(num_next_state), zeros(num_next_state, num_state), zeros(num_next_state, num_action))
end


text:
# Methods

if VERSION < v"1.2-"
const syms_fieldname = :sparam_syms
else
const syms_fieldname = :slot_syms
end

if VERSION < v"1.2-"
const _uncompress = Base.uncompressed_ast
else
const _uncompress = Base._uncompressed_ast
end

structdata(meth::Method) =
  [meth.module, meth.name, meth.file, meth.line, meth.sig, getfield(meth, syms_fieldname),
   meth.ambig, meth.nargs, meth.isva, meth.nospecialize,
   _uncompress(meth, meth.source)]

initstruct(::Type{Method}) = ccall(:jl_new_method_uninit, Ref{Method}, (Any,), Main)

function newstruct!(meth::Method, mod, name, file, line, sig,
                    syms, ambig, nargs, isva, nospecialize, ast)
  meth.module = mod
  meth.name = name
  meth.file = file
  meth.line = line
  meth.sig = sig
  setfield!(meth, syms_fieldname, syms)
  meth.ambig = ambig
  meth.nospecialize = nospecialize
  meth.nargs = nargs
  meth.isva = isva
  meth.source = ast
  meth.pure = ast.pure
  return meth
end

function structdata(t::TypeName)
  primary = Base.unwrap_unionall(t.wrapper)
  mt = !isdefined(t, :mt) ? nothing :
    [t.mt.name, collect(Base.MethodList(t.mt)), t.mt.max_args,
     isdefined(t.mt, :kwsorter) ? t.mt.kwsorter : nothing]
  [Base.string(VERSION), t.name, t.names, primary.super, primary.parameters,
   primary.types, isdefined(primary, :instance), primary.abstract,
   primary.mutable, primary.ninitialized, mt]
end

# Type Names

baremodule __deserialized_types__ end

function newstruct_raw(cache, ::Type{TypeName}, d)
  name = raise_recursive(d[:data][2], cache)
  name = isdefined(__deserialized_types__, name) ? gensym() : name
  tn = ccall(:jl_new_typename_in, Ref{Core.TypeName}, (Any, Any),
             name, __deserialized_types__)
  cache[d] = tn
  names, super, parameters, types, has_instance,
    abstr, mutabl, ninitialized = map(x -> raise_recursive(x, cache), d[:data][3:end-1])
  tn.names = names
  ndt = ccall(:jl_new_datatype, Any, (Any, Any, Any, Any, Any, Any, Cint, Cint, Cint),
              tn, tn.module, super, parameters, names, types,
              abstr, mutabl, ninitialized)
  ty = tn.wrapper = ndt.name.wrapper
  ccall(:jl_set_const, Cvoid, (Any, Any, Any), tn.module, tn.name, ty)
  if has_instance && !isdefined(ty, :instance)
    # use setfield! directly to avoid `fieldtype` lowering expecting to see a Singleton object already on ty
    Core.setfield!(ty, :instance, ccall(:jl_new_struct, Any, (Any, Any...), ty))
  end
  mt = raise_recursive(d[:data][end], cache)
  if mt != nothing
    mtname, defs, maxa, kwsorter = mt
    tn.mt = ccall(:jl_new_method_table, Any, (Any, Any), name, tn.module)
    tn.mt.name = mtname
    tn.mt.max_args = maxa
    for def in defs
      isdefined(def, :sig) || continue
      ccall(:jl_method_table_insert, Cvoid, (Any, Any, Ptr{Cvoid}), tn.mt, def, C_NULL)
    end
  end
  return tn
end

# Function Types

# Modelled on Base.Serialize
function isanon(t::DataType)
  tn = t.name
  if isdefined(tn, :mt)
    name = tn.mt.name
    mod = tn.module
    return t.super === Function &&
    unsafe_load(Base.unsafe_convert(Ptr{UInt8}, tn.name)) == UInt8('#') &&
    (!isdefined(mod, name) || t != typeof(getfield(mod, name)))
  end
  return false
end

# Called externally from lower(::DataType)
function lower_anon(T::DataType)
  BSONDict(:tag => "jl_anonymous",
           :typename => T.name,
           :params => [T.parameters...])
end

tags[:jl_anonymous] = function (d)
  constructtype(d[:typename].wrapper, d[:params])
end


text:
struct WhiteBlackLists
    whitelist_module::Indices{Module}
    whitelist_function::Indices{Symbol}
    whitelist_globalref::Indices{Core.GlobalRef}
    blacklist_module::Indices{Module}
    blacklist_function::Indices{Symbol}
    blacklist_globalref::Indices{Core.GlobalRef}
end

function WhiteBlackLists()
    WhiteBlackLists(Indices{Module}(), Indices{Symbol}(), Indices{Core.GlobalRef}(),
        Indices{Module}(), Indices{Symbol}(), Indices{Core.GlobalRef}())
end

function iswhite(l::WhiteBlackLists, ex::GlobalRef)
    ex.mod ∈ l.whitelist_module && return(true)
    ex ∈ l.whitelist_globalref && return(true)
    ex.name ∈ l.whitelist_function && return(true)
    return(false)
end

function isblack(l::WhiteBlackLists, ex::GlobalRef)
    ex.mod ∈ l.blacklist_module && return(true)
    ex ∈ l.blacklist_globalref && return(true)
    ex.name ∈ l.blacklist_function && return(true)
    return(false)
end

function Base.isempty(l::WhiteBlackLists)
    isempty(l.whitelist_module) && isempty(l.whitelist_function) && isempty(l.whitelist_function)
end

whitelist!(l::WhiteBlackLists, ex::GlobalRef) = insert!(l.whitelist_globalref, ex)
whitelist!(l::WhiteBlackLists, ex::Symbol) = insert!(l.whitelist_function, ex)
whitelist!(l::WhiteBlackLists, ex::Module) = insert!(l.whitelist_module, ex)
blacklist!(l::WhiteBlackLists, ex::GlobalRef) = insert!(l.blacklist_globalref, ex)
blacklist!(l::WhiteBlackLists, ex::Symbol) = insert!(l.blacklist_function, ex)
blacklist!(l::WhiteBlackLists, ex::Module) = insert!(l.blacklist_module, ex)

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler


text:
using StochasticSearch, FactCheck

facts("[NumberParameter]") do
    context("constructor") do
        p = NumberParameter{Int8}(convert(Int8, 0), convert(Int8, 2),
                                  convert(Int8, 1), "test")
        @fact (typeof(p.min)   == Int8)               --> true
        @fact (typeof(p.max)   == Int8)               --> true
        @fact (typeof(p.value) == Int8)               --> true
        @fact (p.name          == "test")             --> true
        @fact (typeof(p) <: NumberParameter{Int8})    --> true
        p = NumberParameter{Float32}(convert(Float32, 0), convert(Float32, 2),
                                     convert(Float32, 1), "test")
        @fact (typeof(p.min)   == Float32)            --> true
        @fact (typeof(p.max)   == Float32)            --> true
        @fact (typeof(p.value) == Float32)            --> true
        @fact (typeof(p) <: NumberParameter{Float32}) --> true
        p = NumberParameter{Int64}(1, 7, 2, "test")
        @fact (typeof(p) == NumberParameter{Int64})   --> true
        p = NumberParameter{Float64}(1.2, 7.2, 2.34, "test")
        @fact (typeof(p) == NumberParameter{Float64}) --> true
        p = NumberParameter{Int64}(1, "test")
        @fact (p.min == 1)                            --> true
        @fact (p.max == 1)                            --> true
        p = NumberParameter{Int64}(0, 10, "test")
        @fact (0 <= p.value <= 10)                    --> true
    end
    context("neighbor!") do
        p = NumberParameter{Int64}(1, 200, 100, "test")
        n = p.value
        neighbor!(p)
        @fact n                                       --> not(exactly(p.value))
        n = p.value
        neighbor!(p, interval = 8)
        @fact n                                       --> not(exactly(p.value))
        n = p.value
        neighbor!(p, interval = 8, distance = 1)
        @fact n                                       --> not(exactly(p.value))
        p = NumberParameter{Float64}(1.332, 60.2, 44.3, "test")
        n = p.value
        neighbor!(p)
        @fact n                                       --> not(exactly(p.value))
        n = p.value
        neighbor!(p, interval = 3.2231)
        @fact n                                       --> not(exactly(p.value))
        n = p.value
        neighbor!(p, interval = 8.332, distance = 1)
        @fact n                                       --> not(exactly(p.value))
        @fact_throws MethodError neighbor!(p, 8.332, 20.2)
    end
    context("constructor") do
        @fact (IntegerParameter <: NumberParameter)   --> true
        p = IntegerParameter(0, 10, 3, "test")
        @fact (typeof(p) <: NumberParameter{Integer}) --> true
        @fact (typeof(p) <: IntegerParameter)         --> true
        @fact (p.min   == 0 )                         --> true
        @fact (p.max   == 10)                         --> true
        @fact (p.value == 3 )                         --> true
        @fact (p.name  == "test")                     --> true
        @fact_throws ErrorException IntegerParameter(3, 1, 2, "test")
        @fact_throws ErrorException IntegerParameter(1, 3, 0, "test")
        @fact_throws ErrorException IntegerParameter(1, 3, 4, "test")
        @fact_throws MethodError    IntegerParameter("a", 3, 4, "test")
        @fact_throws MethodError    IntegerParameter(1, 3, 2)
        @fact_throws MethodError    IntegerParameter(1.0, 3, 4, "test")
        p = IntegerParameter(convert(Int8, 0), convert(Int8, 2),
                             convert(Int8, 1), "test")
        @fact (typeof(p) == NumberParameter{Integer}) --> true
        @fact (typeof(p) <: IntegerParameter)         --> true
        @fact (typeof(p.min)   == Int8)               --> true
        @fact (typeof(p.max)   == Int8)               --> true
        @fact (typeof(p.value) == Int8)               --> true
        @fact (p.min   == 0)                          --> true
        @fact (p.max   == 2)                          --> true
        @fact (p.value == 1)                          --> true
    end
    context("perturb!") do
        value    = 3
        interval = 10
        p = IntegerParameter(0, 100, value, "rand_test")
        perturb!(p, 10)
        @fact (typeof(p.value) <: Integer)  --> true
        @fact (p.value <= p.max)            --> true
        @fact (p.value >= p.min)            --> true
        perturb!(p, interval)
        value = p.value
        @fact (typeof(p.value) <: Integer)  --> true
        @fact (p.value <= p.max)            --> true
        @fact (p.value <= value + interval) --> true
        @fact (p.value >= p.min)            --> true
        @fact (p.value >= value - interval) --> true
        interval = 103
        perturb!(p, interval)
        @fact (typeof(p.value) <: Integer)  --> true
        @fact (p.value <= p.max)            --> true
        @fact (p.value >= p.min)            --> true
        interval = -1
        @fact_throws ErrorException perturb!(p, interval)
        interval = 0
        @fact_throws ErrorException perturb!(p, interval)
    end
    context("constructor") do
        p = FloatParameter(0.223, 10.122, 3.12, "test")
        @fact (p.min   == 0.223 ) --> true
        @fact (p.max   == 10.122) --> true
        @fact (p.value == 3.12  ) --> true
        @fact (p.name  == "test") --> true
        @fact_throws ErrorException FloatParameter(3.3, 1.223, 2.4, "test")
        @fact_throws ErrorException FloatParameter(1.443, 3.2332, 1.442, "test")
        @fact_throws ErrorException FloatParameter(1.23, 3.2, 3.23, "test")
        @fact_throws MethodError    FloatParameter("a", 3, 4, "test")
        @fact_throws MethodError    FloatParameter(1, 3, 2, "test")
        @fact_throws MethodError    FloatParameter(1, 3, 2)
        p = FloatParameter(convert(Float32, 0), convert(Float32, 2),
                           convert(Float32, 1), "test")
        @fact (typeof(p.min)   == Float32)                  --> true
        @fact (typeof(p.max)   == Float32)                  --> true
        @fact (typeof(p.value) == Float32)                  --> true
        @fact (typeof(p) == NumberParameter{AbstractFloat}) --> true
        @fact (typeof(p) <: FloatParameter)                 --> true
        @fact (p.min   == 0)                                --> true
        @fact (p.max   == 2)                                --> true
        @fact (p.value == 1)                                --> true
    end
    context("perturb!") do
        value    = 3.2
        interval = 50.6
        p = FloatParameter(0., 1000., value, "rand_test")
        perturb!(p, 10.)
        @fact (typeof(p.value) <: AbstractFloat) --> true
        @fact (p.value <= p.max)                 --> true
        @fact (p.value >= p.min)                 --> true
        perturb!(p, interval)
        perturb!(p, interval)
        perturb!(p, interval)
        value = p.value
        @fact (typeof(p.value) <: AbstractFloat) --> true
        @fact (p.value <= value + interval)      --> true
        @fact (p.value <= p.max)                 --> true
        @fact (p.value >= p.min)                 --> true
        @fact (p.value >= value - interval)      --> true
        interval = 103.0
        perturb!(p, interval)
        perturb!(p, interval)
        perturb!(p, interval)
        @fact (typeof(p.value) <: AbstractFloat) --> true
        @fact (p.value <= p.max)                 --> true
        @fact (p.value >= p.min)                 --> true
        interval = -1.2
        @fact_throws ErrorException perturb!(p, interval)
        interval = 0
        @fact_throws ErrorException perturb!(p, interval)
    end
end

facts("[EnumParameter]") do
    p = EnumParameter([IntegerParameter(1, 3, 2, "test")], "test")
    @fact (typeof(p.current) == IntegerParameter)               --> true
    @fact (typeof(p)         <: EnumParameter)                  --> true
    p = EnumParameter([IntegerParameter(1, 3, 2, "test"),
                       FloatParameter(1.1,2.3,1.4, "test")],
                       "test")
    @fact (typeof(p.current) <: NumberParameter)                --> true
    @fact (typeof(p)         <: EnumParameter)                  --> true
    p = EnumParameter([IntegerParameter(1, 3, 2, "test"),
                       FloatParameter(1.1,2.3,1.4, "test"),
                       StringParameter("value", "test")],
                       "test")
    @fact (typeof(p.current) <: Parameter)                      --> true
    @fact (typeof(p)         <: EnumParameter)                  --> true
    p = EnumParameter([EnumParameter([StringParameter("a", "test")], "test"),
                       EnumParameter([StringParameter("b", "test")], "test"),
                       EnumParameter([StringParameter("c", "test")], "test")],
                       "test")
    @fact (typeof(p)       <: EnumParameter)                    --> true
    p = EnumParameter([EnumParameter([StringParameter("a", "test")], "test"),
                       StringParameter("b", "test"),
                       EnumParameter([StringParameter("c", "test")], "test")],
                       "test")
    println(p)
    @fact (typeof(p)       <: EnumParameter)                    --> true
    context("constructors") do
        p = EnumParameter([IntegerParameter(1, 4, 3, "test"),
                           IntegerParameter(1, 6, 3, "test")], "test")
        @fact (typeof(p)         <: EnumParameter)              --> true
        @fact (typeof(p.values)  <: AbstractArray)              --> true
        @fact (typeof(p.values)  == Array{IntegerParameter, 1}) --> true
        @fact (p.values[1].value == 3)                          --> true
        @fact (p.values[2].value == 3)                          --> true
        @fact (p.values[1].name  == "test")                     --> true
        @fact (p.values[2].name  == "test")                     --> true
        @fact (p.name            == "test")                     --> true
        p = EnumParameter([IntegerParameter(1, 4, 3, "test"),
                           IntegerParameter(1, 6, 2, "test")], 1, "test")
        @fact (p.value           == 1)                          --> true
        @fact (p.current.value   == 3)                          --> true
        p = EnumParameter([FloatParameter(1.2, 4.3, 3.2, "test"),
                           IntegerParameter(1, 6, 2, "test")], 1, "test")
        @fact (typeof(p)         <: EnumParameter)              --> true
        p = EnumParameter([FloatParameter(1.2, 4.3, 3.2, "test"),
                           IntegerParameter(1, 6, 2, "test"),
                           StringParameter("value", "test")], "test")
        @fact (typeof(p)         <: EnumParameter)              --> true
        @fact_throws MethodError    EnumParameter([3, 4], "test")
        @fact_throws MethodError    EnumParameter([IntegerParameter(1, 4, 3, "test")])
        @fact_throws ErrorException EnumParameter([IntegerParameter(1, 4, 3, "test"),
                                                   IntegerParameter(1, 6, 3, "test")], 3, "test")
    end
    context("perturb! and perturb_elements!") do
        p = EnumParameter([IntegerParameter(1, 4, 3, "test"),
                           StringParameter("b", "test"),
                           IntegerParameter(1, 6, 3, "test")], 2, "test")
        perturb_elements!(p)
        perturb_elements!(p)
        perturb_elements!(p)
        @fact (p.value == 2)                                --> true
        p = EnumParameter([FloatParameter(1.1, 4.2, 3.221, "test"),
                           FloatParameter(2.3, 4.4, 3.1, "test")], "test")
        perturb!(p)
        perturb!(p)
        perturb!(p)
        @fact (p.values[1].value == 3.221)                  --> true
        p = EnumParameter([StringParameter("valuea", "test"),
                           StringParameter("valueb", "test")], "test")
        perturb!(p)
        perturb!(p)
        perturb!(p)
        @fact (p.values[1].value == "valuea")               --> true
        p = EnumParameter([EnumParameter([StringParameter("a", "test")], "test"),
                           StringParameter("b", "test"),
                           IntegerParameter(1, 4, 3, "test"),
                           EnumParameter([StringParameter("c", "test")], "test")],
                           "test")
        v = p.current
        perturb_elements!(p, 3, 2)
        perturb_elements!(p, 3, 2)
        perturb_elements!(p, 3, 2)
        @fact (p.current == v)                              --> true
        v = p.values[3]
        @fact (v.value <= v.max)                            --> true
        @fact (v.value >= v.min)                            --> true
        perturb_elements!(p, [2, 3, 4])
        v = p.values[3]
        @fact (v.value <= v.max)                            --> true
        @fact (v.value >= v.min)                            --> true
        @fact_throws MethodError    perturb_elements!(p, 2, 2)
        @fact_throws MethodError    perturb_elements!(p, 2)
        @fact_throws ErrorException perturb_elements!(p, [1,3,4,5,6,3])
    end
    context("neighbor!") do
        p = EnumParameter([IntegerParameter(1, 4, 3, "test"),
                           StringParameter("b", "test"),
                           IntegerParameter(1, 6, 3, "test")], 2, "test")
        n = p.current
        neighbor!(p)
        @fact n                                       --> not(exactly(p.value))
        @fact_throws MethodError neighbor!(p, 8, 20)
        @fact_throws MethodError neighbor!(p, 3.2231)
    end
end

facts("[PermutationParameter]") do
    context("constructor") do
        a = [1, 2, 3, 4, 5, 6]
        p = PermutationParameter(a, "test")
        @fact (typeof(p) <: PermutationParameter) --> true
        @fact (p.value == a)                      --> true
        @fact (p.size == length(a))               --> true
    end
    context("neighbor!") do
        a = [1, 2, 3, 4, 5, 6]
        b = [1, 2, 3, 4, 5, 6]
        p = PermutationParameter(a, "test")
        @fact p.value == b                        --> true
        print(p)
        neighbor!(p)
        print(p)
        print(p)
        neighbor!(p, interval = 2)
        print(p)
        @fact p.value != b                        --> true
        @fact p.size == length(b)                 --> true
    end
    context("perturb!") do
        a = [1, 2, 3, 4, 5, 6]
        b = [1, 2, 3, 4, 5, 6]
        p = PermutationParameter(a, "test")
        @fact p.value == b                        --> true
        print(p)
        perturb!(p)
        print(p)
        @fact p.value != b                        --> true
        @fact p.size == length(b)                 --> true
    end
end

facts("[StringParameter]") do
    p = StringParameter("value", "test")
    @fact (typeof(p) == StringParameter) --> true
    @fact (typeof(p) <: Parameter)       --> true
    @fact (p.value == "value")           --> true
    @fact (p.name  == "test")            --> true
    @fact_throws MethodError perturb!(p)
    @fact_throws MethodError StringParameter(2, "test")
    @fact_throws MethodError StringParameter("value")
end

facts("[BoolParameter]") do
    p = BoolParameter(false, "test")
    @fact (typeof(p.value) == Bool)      --> true
    p = BoolParameter(0, "test")
    @fact (p.value == false)             --> true
    p = BoolParameter(1, "test")
    @fact (p.value == true)              --> true
    @fact_throws MethodError BoolParameter("must have start value")
    context("perturb!") do
        p = BoolParameter(false, "test")
        perturb!(p)
        @fact (p.value == true)          --> true
    end
    context("neighbor!") do
        p = BoolParameter(false, "test")
        neighbor!(p)
        @fact (p.value == true)          --> true
        print(p)
    end
end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

using JSON

export get_credentials

function scan_for_credentials(path)

    # Get storage and batch accounts
    batch_accounts = []
    storage_accounts = []
    files = readdir(path)

    for entry in files
        if ~isempty(findall("aad_", entry))
            push!(batch_accounts, entry)
        elseif ~isempty(findall("key_", entry))
            push!(storage_accounts, entry)
        end
    end
    return storage_accounts, batch_accounts
end

function get_storage_credentials(path, credential_file)

    io = open(joinpath(path, credential_file), "r")
    creds = JSON.parse(io)
    storage_account_name = creds[3]["accountName"]
    storage_secret_key = creds[1]["value"]

    return storage_account_name, storage_secret_key
end

function get_batch_credentials(path, credential_file)

    io = open(joinpath(path, credential_file), "r")
    creds = JSON.parse(io)
    batch_account_name = creds[2]["accountName"]
    region = creds[2]["region"]
    tenant = creds[1]["tenant"]
    appId = creds[1]["appId"]
    password = creds[1]["password"]

    return batch_account_name, region, tenant, appId, password
end

# Scan directory for given path for credential files and return list of credentials
function get_credentials(path)

    storage_accounts, batch_accounts = scan_for_credentials(path)
    credentials = []

    if length(storage_accounts) == length(batch_accounts)
        for (storage, batch) in zip(storage_accounts, batch_accounts)
            storage_account, storage_secret = get_storage_credentials(path, storage)
            batch_account, region, aad_tenant, app_id, aad_secret = get_batch_credentials(path, batch)

            credential_dict = Dict{String, Any}(
                "_AD_TENANT" => aad_tenant,
                "_BATCH_ACCOUNT_NAME" => batch_account,
                "_REGION" => region,
                "_AD_BATCH_CLIENT_ID" => app_id,
                "_AD_SECRET_BATCH" => aad_secret,
                "_BATCH_ACCOUNT_URL" => join(["https://", batch_account, ".", region, ".batch.azure.com"]),
                "_BATCH_RESOURCE" => "https://batch.core.windows.net/",
                "_STORAGE_ACCOUNT_NAME" => storage_account,
                "_STORAGE_ACCOUNT_KEY" => storage_secret
            )
            push!(credentials, credential_dict)
        end
    elseif length(batch_accounts) == 1 && length(storage_accounts) >= 1

        batch_account, region, aad_tenant, app_id, aad_secret = get_batch_credentials(path, batch_accounts[1])
        for storage in storage_accounts
            storage_account, storage_secret = get_storage_credentials(path, storage)

            credential_dict = Dict{String, Any}(
                "_AD_TENANT" => aad_tenant,
                "_BATCH_ACCOUNT_NAME" => batch_account,
                "_REGION" => region,
                "_AD_BATCH_CLIENT_ID" => app_id,
                "_AD_SECRET_BATCH" => aad_secret,
                "_BATCH_ACCOUNT_URL" => join(["https://", batch_account, ".", region, ".batch.azure.com"]),
                "_BATCH_RESOURCE" => "https://batch.core.windows.net/",
                "_STORAGE_ACCOUNT_NAME" => storage_account,
                "_STORAGE_ACCOUNT_KEY" => storage_secret
            )
            push!(credentials, credential_dict)
        end
    else
        throw("Need to provide at least one storage account per batch account.")
    end

    return credentials
end


text:
using SolidStateDetectors

#define the charge drift model in your custom json file
cdm = ADLChargeDriftModel()

#define a range of relevant electric field strengths in units of V/m
E = Float32.(0:100:1000000)

#print(cdm.masses)
using StaticArrays
using LinearAlgebra

function SolidStateDetectors.Vl(Emag::T, params::SolidStateDetectors.VelocityParameters{T})::T where {T <: SolidStateDetectors.SSDFloat}
    #define your own parametrization here: e.g. vd = 500 m^2/(Vs) * E
    #beware of units!!
    return 500 * Emag
end

newE = [T(i)*SVector{3,T}(sin(cdm.phi110+π/4),cos(cdm.phi110+π/4),0) for i in E]
ve100a = map(x -> norm(SolidStateDetectors.getVe(x, cdm)), newE)
ve100b = map(x -> SolidStateDetectors.Vl(x,cdm.electrons.axis100), T.(E))

plot(E, [ve100a, ve100b],
    xlabel = "Electric field strength [V/m]", ylabel = "Drift velocity [m/s]")

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
# Construct a sample using Reservoir sampling

type Sample{T} <: StreamStat
    sample::Vector{T}
    n::Int
end

Sample(T::DataType, k::Integer) = Sample(Array(T, k), 0)

function update!(stat::Sample, x::Any)
    k = length(stat.sample)

    stat.n += 1

    if stat.n <= k
        stat.sample[stat.n] = x
    else
       i = rand(1:stat.n)
       if i <= k
           stat.sample[i] = x
       end
    end

    return
end

# This should really return a Nullable
function state{T}(stat::Sample{T})
   if stat.n < length(stat.sample)
       return T[]
   else
       return stat.sample
   end
end

StatsBase.sample{T}(stat::Sample{T}) = state(stat)

nobs(stat::Sample) = stat.n

Base.copy(stat::Sample) = Sample(stat.sample, stat.n)

# TODO: Implement random mixture
# Base.merge(a::Std, b::Std) = Std(merge(a.var, b.var))

function Base.empty!(stat::Sample)
    empty!(stat.sample)
    stat.n = 0
    return
end

function Base.show(io::IO, stat::Sample)
    @printf(io, "Online Random Sample\n")
    @printf(io, " * K: %d\n", length(stat.sample))
    @printf(io, " * N: %d\n", stat.n)
    return
end