text:
################################################# FILE DESCRIPTION ##########################################################

# This file contains the combinatorial API for Graph.

################################################# IMPORT/EXPORT #############################################################

export nv, ne, vertices, edges, hasvertex, hasedge, fadj, fadj!, outdegree, indegree

################################################# COMBINATORIAL BASICS ######################################################

""" Return nv(g) x ne(g) """
Base.size(g::Graph) = (g.nv,g.ne)

""" The list of the vertices in the graph """
vertices(g::Graph) = 1 : nv(g)

""" Returns an edge iterator containing all edges in the graph """
edges(g::Graph) = EdgeIter(indxs(g))

""" Check if the vertex(s) exists """
hasvertex(g::Graph, v::VertexID) = 1 <= v <= nv(g)
hasvertex(g::Graph, vs::VertexList) = 1 .<= vs .<= nv(g)

""" Check if the edge(s) exists """
hasedge(g::Graph, e::EdgeID) = indxs(g)[e] > 0
hasedge(g::Graph, es::EdgeList) = indxs(g)[es] .> 0

################################################# ADJACENCY #################################################################

""" Vertex v's out-neighbors in the graph """
fadj(g::Graph, v::VertexID) = fadj(indxs(g), v)

"""
Retrieve a list of vertices connect to vertex v.

This method copies the adjacencies onto the input array,
and is comparitively faster, and causes no mallocs.
"""
fadj!(g::Graph, v::VertexID, adj::Vector{Int}) = fadj!(indxs(g), v, adj)

""" Vertex v's outdegree in the graph """
outdegree(g::Graph, v::VertexID) = outdegree(indxs(g), v)

""" Outdegree of vertex v for v in vs """
outdegree(g::Graph, vs::VertexList) = outdegree(indxs(g), vs)

""" Outdegrees of all the vertices in the graph """
outdegree(g::Graph) = outdegree(indxs(g))

""" Vertex v's indegree in the graph """
indegree(g::Graph, v::VertexID) = indegree(indxs(g), v)

""" Indegree of vertex v for v in vs """
indegree(g::Graph, vs::VertexList) = indegree(indxs(g), vs)

""" Indegrees of all the vertices in the graph """
indegree(g::Graph) = indegree(indxs(g))

""" Shorcut to vertex v's out-neighbors in the graph """
function Base.getindex(g::Graph, v)
   encode(g, fadj(g, decode(g, v)))
end


text:
include("../examples/nonconvex_smooth.jl")
using PyPlot
using Printf
function plot_loss_landscape(s=4.0)
	npts = 10000
	w = rand(npts)
	lw = rand(npts)
	fig, ax = subplots()
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)
	for (i, wi) = enumerate(w)
       lw[i] = loss(wi,s)
	end
	ax.plot(w, lw,".")
end
function plot_bifurcation()
	npts = 50
	wstar = rand(npts)
	nSteps = 50000
	η = LinRange(0.0001, 0.1, 20)
	fig, ax = subplots()
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	for (i, ηi) = enumerate(η)
		@show i, ηi
		wstar .= rand(npts)
		for j = 1:npts
			wstar[j] = run(wstar[j], ηi, nSteps)
		end
		ax.plot(ηi*ones(npts), wstar, "bo", ms=3.0)
	end
	ax.grid(true)
	ax.set_xlabel(L"$ \eta $", fontsize=30)
	ax.set_ylabel(L"$ w^* $", fontsize=30)
end

ext: jl

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
#= inlet bc selector
3: from data
=#
const inlet_BC_switch = 3
const inlet_type = "Q"

# numerical domain
const Ccfl   = 0.9         # Courant number
const cycles = 15

# blood properties
const rho = 1060.    # density [kg/m3]
const mu  = 4.e-3   # dynamic viscosity [Pa⋅s]
const gamma_profile = 9

# vessel properties
#const Pext = 0.			     # external pressure [Pa]
const initial_pressure = 0.


text:
## ExoplanetsSysSim/src/target.jl
## (c) 2015 Eric B. Ford

#using Distributions

struct KeplerTarget
  #sys::PlanetarySystem   # Make array for planetary systems aroud multiple stars in one target?
  sys::Vector{PlanetarySystemAbstract}
  cdpp::Array{Float64,2} # fractional, not ppm; 2D to allow for multiple time scales, months, quarters or seasons/spacecraft rotation angles
                           # QUERY: Do we want a separate CDPP for SC?  Or will CDPP's in different months be for LC/SC depending on this variable?
                            # QUERY: Should this be moved to KeplerTargetObs?
                            # QUERY: Should we not add this to target and use the star id to lookup CDPP from the stellar table?
  contam::Float64     # QUERY: Do we want/need this, since we're able to generate multiple stars in a single target?
  data_span::Float64
  duty_cycle::Float64
  window_function_id::Int64 # Points to the id of the window function for this target
  #channel::Int64         # E.g., if we cared which Kepler channel the target fell on
  #has_sc::Vector{Bool}   # TODO OPT: Make Immutable Vector or BitArray for speed?  QUERY: Should this go in KeplerTargetObs?
  #                       # QUERY: Do we want a separate CDPP for SC?  Or will CDPP's in different months be for LC/SC depending on this variable?
  #ra::Float64           # E.g., if we cared about position on sky     QUERY:  Should we replace with galactic longitude and latitute?
  #dec::Floa64           #
end
num_planets(t::KeplerTarget) = sum( num_planets, t.sys)
flux(t::KeplerTarget) = sum(flux,t.sys)+t.contam

function draw_asymmetric_normal(mu::Real, sig_plus::Real, sig_minus::Real)
  @assert sig_minus >= zero(sig_minus)
  stdn = randn()
  mu + ( (stdn>=zero(stdn)) ? sig_plus*stdn : sig_minus*stdn )
  #mu + ( (stdn>=zero(stdn)) ? sig_plus*stdn : abs(sig_minus)*stdn )
end

function make_cdpp_array(star_id::Integer)
  star_table(id::Integer,sym::Symbol) = StellarTable.star_table(id,sym)::Float64
  cdpp_arr = (1.0e-6*sqrt(1.0/24.0/LC_duration)) .* Float64[star_table(star_id, :rrmscdpp01p5)*sqrt(1.5), star_table(star_id, :rrmscdpp02p0)*sqrt(2.), star_table(star_id,:rrmscdpp02p5)*sqrt(2.5), star_table(star_id,:rrmscdpp03p0)*sqrt(3.), star_table(star_id,:rrmscdpp03p5)*sqrt(3.5), star_table(star_id,:rrmscdpp04p5)*sqrt(4.5), star_table(star_id,:rrmscdpp05p0)*sqrt(5.), star_table(star_id,:rrmscdpp06p0)*sqrt(6.), star_table(star_id,:rrmscdpp07p5)*sqrt(7.5), star_table(star_id,:rrmscdpp09p0)*sqrt(9.), star_table(star_id,:rrmscdpp10p5)*sqrt(10.5), star_table(star_id,:rrmscdpp12p0)*sqrt(12.), star_table(star_id,:rrmscdpp12p5)*sqrt(12.5), star_table(star_id,:rrmscdpp15p0)*sqrt(15.)]
end

function generate_kepler_target_from_table(sim_param::SimParam)
  # generate_star = get_function(sim_param,"generate_star")
  generate_planetary_system = get_function(sim_param,"generate_planetary_system")
  ar_table_sigmas = true
  max_draws_star_properties = 20
  min_star_radius = 0.5
  min_star_mass = 0.5
  max_star_radius = 2.0
  max_star_mass = 2.0
  max_star_density = 1000.0
  max_star_id = StellarTable.num_usable_in_star_table()

  star_table(id::Integer,sym::Symbol) = StellarTable.star_table(id,sym)

  @assert(1<=max_star_id)
  star_id = rand(1:max_star_id)
  mass = 0.0
  dens = 0.0
  radius = 0.0
    #if use_star_table_sigmas
    if get(sim_param,"use_star_table_sigmas",false)
        attmpt_num = 0
        while (!(min_star_radius<radius<max_star_radius))# || (!(min_star_mass<mass<max_star_mass)) || (!(0.0<dens<max_star_density))
            if attmpt_num >= max_draws_star_properties
                star_id = rand(1:max_star_id)
                attmpt_num = 0
            end
            radius = draw_asymmetric_normal( star_table(star_id,:radius), star_table(star_id,:radius_err1), abs(star_table(star_id,:radius_err2)) )
            #mass = draw_asymmetric_normal( star_table(star_id,:mass), star_table(star_id,:mass_err1), abs(star_table(star_id,:mass_err2)) )
            #dens = draw_asymmetric_normal( star_table(star_id,:dens), star_table(star_id,:dens_err1), abs(star_table(star_id,:dens_err2)) )
            attmpt_num += 1
        end

        # ZAMS mass-radius relation taken from 15.1.1 of Allen's Astrophysical Quantities (2002)
        if radius > 1.227
            mass = 10^((log10(radius)-0.011)/0.64)
        else
            mass = 10^((log10(radius)+0.02)/0.917)
        end
        dens   = (mass*sun_mass_in_kg_IAU2010*1000.)/(4//3*pi*(radius*sun_radius_in_m_IAU2015*100.)^3)  # Self-consistent density (gm/cm^3)
  else
    radius = star_table(star_id,:radius)
    mass   = star_table(star_id,:mass)
    #dens   = star_table(star_id,:dens)
    dens   = (mass*sun_mass_in_kg_IAU2010*1000.)/(4//3*pi*(radius*sun_radius_in_m_IAU2015*100.)^3)  # Self-consistent density (gm/cm^3)
  end
  star = SingleStar(radius,mass,1.0, star_id)     # TODO SCI: Allow for blends, binaries, etc.
  cdpp_arr = make_cdpp_array(star_id)
  contam = 0.0 # rand(LogNormal(1.0e-3,1.0))      # TODO SCI: Come up with better description of Kepler targets, maybe draw from real contaminations
  data_span = star_table(star_id, :dataspan)
  duty_cycle = star_table(star_id, :dutycycle)
  if StellarTable.star_table_has_key(:wf_id)
     wf_id = star_table(star_id,:wf_id)
  else
     wf_id = WindowFunction.get_window_function_id(star_table(star_id,:kepid))
  end
  # ch = rand(DiscreteUniform(1,84))              # Removed channel in favor of window function id
  ps = generate_planetary_system(star, sim_param)
  return KeplerTarget([ps],repeat(cdpp_arr, outer=[1,1]),contam,data_span,duty_cycle,wf_id)
end

function generate_kepler_target_simple(sim_param::SimParam)
   generate_star = get_function(sim_param,"generate_star")
   generate_planetary_system = get_function(sim_param,"generate_planetary_system")
   star::StarAbstract = generate_star(sim_param)

   mean_log_cdpp = 4.9759601617565465   # mean frmo star table
   stddev_log_cdpp = 0.6704860437536709 # std dev from star table
  rrmscdpp_5hr = exp(mean_log_cdpp+stddev_log_cdpp*randn())
  cdpp_5hr = 1.0e-6 * rrmscdpp_5hr * sqrt(5.0/24.0 / LC_duration )
  contam = 0.0 # rand(LogNormal(1.0e-3,1.0))   # TODO SCI: Come up with better description of Kepler targets, maybe draw from real contaminations
  wf_id = 0
  # ch = rand(DiscreteUniform(1,84))              # Removed channel in favor of window function id
  ps = generate_planetary_system(star, sim_param)
  return KeplerTarget([ps],fill(cdpp_5hr,num_cdpp_timescales,num_quarters),contam,mission_data_span,mission_duty_cycle,0)
end

function test_target(sim_param::SimParam)
  generate_kepler_target_simple(sim_param)
  StellarTable.setup_star_table(sim_param)
  generate_kepler_target_from_table(sim_param)
end


text:
"""
```
estimate(m, data; verbose=:low, proposal_covariance=Matrix())
```

Estimate the DSGE parameter posterior distribution.

### Arguments:
- `m::AbstractModel`: model object

### Optional Arguments:
- `data`: well-formed data as `Matrix` or `DataFrame`. If this is not provided, the `load_data` routine will be executed.

### Keyword Arguments:
- `verbose::Symbol`: The desired frequency of function progress messages printed to standard out.
   - `:none`: No status updates will be reported.
   - `:low`: Status updates will be provided in csminwel and at each block in
     Metropolis-Hastings.
   - `:high`: Status updates provided at each iteration in Metropolis-Hastings.
- `proposal_covariance::Matrix`: Used to test the metropolis_hastings algorithm with a precomputed
  covariance matrix for the proposal distribution. When the Hessian is singular,
  eigenvectors corresponding to zero eigenvectors are not well defined, so eigenvalue
  decomposition can cause problems. Passing a precomputed matrix allows us to ensure that
  the rest of the routine has not broken.
- `mle`: Set to true if parameters should be estimated by maximum likelihood directly.
    If this is set to true, this function will return after estimating parameters.
- `run_MH`: Set to false to disable sampling from the posterior.
"""
function estimate(m::AbstractModel, df::DataFrame;
                  verbose::Symbol=:low,
                  proposal_covariance::Matrix=Matrix(0,0),
                  mle::Bool = false,
                  run_MH::Bool = true)
    data = df_to_matrix(m, df)
    estimate(m, data; verbose=verbose, proposal_covariance=proposal_covariance,
             mle = mle, run_MH = run_MH)
end
function estimate(m::AbstractModel;
                  verbose::Symbol=:low,
                  proposal_covariance::Matrix=Matrix(0,0),
                  mle::Bool = false,
                  run_MH::Bool = true)
    # Load data
    df = load_data(m; verbose=verbose)
    estimate(m, df; verbose=verbose, proposal_covariance=proposal_covariance,
             mle = mle, run_MH = run_MH)
end
function estimate(m::AbstractModel, data::Matrix{Float64};
                  verbose::Symbol=:low,
                  proposal_covariance::Matrix=Matrix(0,0),
                  mle::Bool = false,
                  run_MH::Bool = true)

    ########################################################################################
    ### Step 1: Find posterior/likelihood mode (if reoptimizing, run optimization routine)
    ########################################################################################

    # Specify starting mode

    vint = get_setting(m, :data_vintage)
    if reoptimize(m)
        println("Reoptimizing...")

        # Inputs to optimization algorithm
        n_iterations       = get_setting(m, :optimization_iterations)
        ftol               = get_setting(m, :optimization_ftol)
        xtol               = get_setting(m, :optimization_xtol)
        gtol               = get_setting(m, :optimization_gtol)
        step_size          = get_setting(m, :optimization_step_size)
        converged          = false

        # If the algorithm stops only because we have exceeded the maximum number of
        # iterations, continue improving guess of modal parameters
        total_iterations = 0
        optimization_time = 0
        max_attempts = get_setting(m, :optimization_attempts)
        attempts = 1

        while !converged
            tic()
            out, H = optimize!(m, data;
                               method = get_setting(m, :optimization_method),
                               ftol=ftol, grtol = gtol, xtol = xtol,
                               iterations=n_iterations, show_trace=true, step_size=step_size,
                               verbose=verbose,
                               mle = mle)

            attempts += 1
            total_iterations += out.iterations
            converged = out.converged || attempts > max_attempts

            if VERBOSITY[verbose] >= VERBOSITY[:low]
                @printf "Total iterations completed: %d\n" total_iterations
                @printf "Optimization time elapsed: %5.2f\n" optimization_time += toq()
            end

            # Write params to file after every `n_iterations` iterations
            params = map(θ->θ.value, m.parameters)
            h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
                file["params"] = params
            end
        end

        # write parameters to file one last time so we have the final mode
        h5open(rawpath(m, "estimate", "paramsmode.h5"),"w") do file
            file["params"] = params
        end
    end

    params = map(θ->θ.value, m.parameters)

    # Sampling does not make sense if mle=true
    if mle || !run_MH
        return nothing
    end

    ########################################################################################
    ### Step 2: Compute proposal distribution
    ###
    ### In Metropolis-Hastings, we draw sample parameter vectors from
    ### the proposal distribution, which is a degenerate multivariate
    ### normal centered at the mode. Its variance is the inverse of
    ### the hessian. We find the inverse via eigenvalue decomposition.
    ########################################################################################

    # Calculate the Hessian at the posterior mode
    hessian = if calculate_hessian(m)
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("Recalculating Hessian...")
        end

        hessian, _ = hessian!(m, params, data; verbose=verbose)

        h5open(rawpath(m, "estimate","hessian.h5"),"w") do file
            file["hessian"] = hessian
        end

        hessian

    # Read in a pre-calculated Hessian
    else
        fn = hessian_path(m)
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("Using pre-calculated Hessian from $fn")
        end

        hessian = h5open(fn,"r") do file
            read(file, "hessian")
        end

        hessian
    end

    # Compute inverse hessian and create proposal distribution, or
    # just create it with the given cov matrix if we have it
    propdist = if isempty(proposal_covariance)
        # Make sure the mode and hessian have the same number of parameters
        n = length(params)
        @assert (n, n) == size(hessian)

        # Compute the inverse of the Hessian via eigenvalue decomposition
        S_diag, U = eig(hessian)
        big_eig_vals = find(x -> x > 1e-6, S_diag)
        rank = length(big_eig_vals)

        S_inv = zeros(n, n)
        for i = (n-rank+1):n
            S_inv[i, i] = 1/S_diag[i]
        end

        hessian_inv = U*sqrt.(S_inv) #this is the inverse of the hessian
        DSGE.DegenerateMvNormal(params, hessian_inv)
    else
        DSGE.DegenerateMvNormal(params, proposal_covariance)
    end

    if DSGE.rank(propdist) != n_parameters_free(m)
        println("problem –    shutting down dimensions")
    end

    ########################################################################################
    ### Step 3: Sample from posterior using Metropolis-Hastings algorithm
    ########################################################################################

    # Set the jump size for sampling
    cc0 = get_setting(m, :mh_cc0)
    cc = get_setting(m, :mh_cc)

    metropolis_hastings(propdist, m, data, cc0, cc; verbose=verbose)

    ########################################################################################
    ### Step 4: Calculate and save parameter covariance matrix
    ########################################################################################

    compute_parameter_covariance(m)

    return nothing
end

"""
```
metropolis_hastings{T<:AbstractFloat}(propdist::Distribution, m::AbstractModel,
    data::Matrix{T}, cc0::T, cc::T; verbose::Symbol = :low)
```

Implements the Metropolis-Hastings MCMC algorithm for sampling from the posterior
distribution of the parameters.

### Arguments

- `propdist`: The proposal distribution that Metropolis-Hastings begins sampling from.
- `m`: The model object
- `data`: Data matrix for observables
- `cc0`: Jump size for initializing Metropolis-Hastings.
- `cc`: Jump size for the rest of Metropolis-Hastings.

### Optional Arguments

- `verbose`: The desired frequency of function progress messages printed to
  standard out. One of:

```
   - `:none`: No status updates will be reported.
   - `:low`: Status updates provided at each block.
   - `:high`: Status updates provided at each draw.
```
"""
function metropolis_hastings{T<:AbstractFloat}(propdist::Distribution,
                                               m::AbstractModel,
                                               data::Matrix{T},
                                               cc0::T,
                                               cc::T;
                                               verbose::Symbol=:low)

    # If testing, set the random seeds at fixed numbers
    if m.testing
        srand(m.rng, 654)
    end

    # Set number of draws, how many we will save, and how many we will burn
    # (initialized here for scoping; will re-initialize in the while loop)
    n_blocks = 0
    n_sim    = 0
    mhthin   = 0
    n_burn   = 0
    n_params = n_parameters(m)

    # Initialize algorithm by drawing para_old from a normal distribution centered on the
    # posterior mode until the parameters are within bounds or the posterior value is
    # sufficiently large.
    para_old = rand(propdist, m; cc=cc0)
    post_old = -Inf

    n_blocks = n_mh_blocks(m)
    n_sim    = n_mh_simulations(m)
    n_burn   = n_mh_burn(m)
    mhthin   = mh_thin(m)

    initialized = false
    while !initialized
        post_old = posterior!(m, para_old, data; mh = true)
        if post_old > -Inf
            propdist.μ = para_old
            initialized = true
        else
            para_old = rand(propdist, m; cc=cc0)
        end
    end

    # Report number of blocks that will be used
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Blocks: $n_blocks")
        println("Draws per block: $n_sim")
    end

    # For n_sim*mhthin iterations within each block, generate a new parameter draw.
    # Decide to accept or reject, and save every (mhthin)th draw that is accepted.
    all_rejections = 0

    # Initialize matrices for parameter draws and transition matrices
    mhparams = zeros(n_sim, n_parameters(m))

    # Open HDF5 file for saving parameter draws
    simfile = h5open(rawpath(m,"estimate","mhsave.h5"),"w")
    n_saved_obs = n_sim * (n_blocks - n_burn)
    parasim = d_create(simfile, "mhparams", datatype(Float32),
                       dataspace(n_saved_obs,n_params),
                       "chunk", (n_sim,n_params))

    # Keep track of how long metropolis_hastings has been sampling
    total_sampling_time = 0.

    for block = 1:n_blocks

        tic()
        block_rejections = 0

        for j = 1:(n_sim*mhthin)

            # Draw para_new from the proposal distribution
            para_new = rand(propdist, m; cc = cc)

            # Solve the model (checking that parameters are within bounds and
            # gensys returns a meaningful system) and evaluate the posterior
            post_new = posterior!(m, para_new, data; mh = true)

            if VERBOSITY[verbose] >= VERBOSITY[:high]
                println("Block $block, Iteration $j: posterior = $post_new")
            end

            # Choose to accept or reject the new parameter by calculating the
            # ratio (r) of the new posterior value relative to the old one We
            # compare min(1, r) to a number drawn randomly from a uniform (0, 1)
            # distribution. This allows us to always accept the new draw if its
            # posterior value is greater than the previous draw's, but it gives
            # some probability to accepting a draw with a smaller posterior
            # value, so that we may explore tails and other local modes.
            r = exp(post_new - post_old)
            x = rand(m.rng)

            if x < min(1.0, r)
                # Accept proposed jump
                para_old = para_new
                post_old = post_new
                propdist.μ = para_new

                if VERBOSITY[verbose] >= VERBOSITY[:high]
                    println("Block $block, Iteration $j: accept proposed jump")
                end
            else
                # Reject proposed jump
                block_rejections += 1

                if VERBOSITY[verbose] >= VERBOSITY[:high]
                    println("Block $block, Iteration $j: reject proposed jump")
                end
            end

            # Save every (mhthin)th draw
            if j % mhthin == 0
                draw_index = convert(Int, j/mhthin)
                mhparams[draw_index, :]  = para_old'
            end
        end # of block

        all_rejections += block_rejections
        block_rejection_rate = block_rejections/(n_sim*mhthin)

        ## Once every iblock times, write parameters to a file

        # Calculate starting and ending indices for this block (corresponds to a new chunk in memory)
        block_start = n_sim*(block-n_burn-1)+1
        block_end   = block_start+n_sim-1

        # Write parameters to file if we're past n_burn blocks
        if block > n_burn
            parasim[block_start:block_end, :] = map(Float32, mhparams)
        end

        # Calculate time to complete this block, average block time, and
        # expected time to completion
        if VERBOSITY[verbose] >= VERBOSITY[:low]
            block_time = toq()
            total_sampling_time += block_time
            total_sampling_time_minutes = total_sampling_time/60
            expected_time_remaining_sec     = (total_sampling_time/block)*(n_blocks - block)
            expected_time_remaining_minutes = expected_time_remaining_sec/60

            println("Completed $block of $n_blocks blocks.")
            println("Total time to compute $block blocks: $total_sampling_time_minutes minutes")
            println("Expected time remaining for Metropolis-Hastings: $expected_time_remaining_minutes minutes")
            println("Block $block rejection rate: $block_rejection_rate \n")
        end

    end # of loop over blocks

    close(simfile)

    rejection_rate = all_rejections / (n_blocks*n_sim*mhthin)
    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("Overall rejection rate: $rejection_rate")
    end
end

"""
```
compute_parameter_covariance(m::AbstractModel)
```

Calculates the parameter covariance matrix from saved parameter draws, and writes it to the
parameter_covariance.h5 file in the `workpath(m, "estimate")` directory.

### Arguments
* `m::AbstractModel`: the model object
"""
function compute_parameter_covariance(m::AbstractModel)

    # Read in saved parameter draws
    param_draws_path = rawpath(m,"estimate","mhsave.h5")
    if !isfile(param_draws_path)
        @printf STDERR "Saved parameter draws not found.\n"
        return
    end
    param_draws = h5open(param_draws_path, "r") do f
        read(f, "mhparams")
    end

    # Calculate covariance matrix
    param_covariance = cov(param_draws)

    # Write to file
    h5open(workpath(m, "estimate","parameter_covariance.h5"),"w") do f
        f["mhcov"] = param_covariance
    end

    return param_covariance
end

"""
```
get_estimation_output_files(m)
```

Returns a `Dict{Symbol, String}` with all files created by `estimate(m)`.
"""
function get_estimation_output_files(m::AbstractModel)
    output_files = Dict{Symbol, String}()

    for file in [:paramsmode, :hessian, :mhsave]
        output_files[file] = rawpath(m, "estimate", "$file.h5")
    end

    for file in [:paramsmean, :parameter_covariance]
        output_files[file] = workpath(m, "estimate", "$file.h5")
    end

    for file in [:priors, :prior_posterior_means, :moments]
        output_files[file] = tablespath(m, "estimate", "$file.tex")
    end

    return output_files
end


text:
# Functions for reading FITS files from the SloanDigitalSkySurvey.
import FITSIO
import WCS

# types of things to mask in read_mask().
const DEFAULT_MASK_PLANES = ["S_MASK_INTERP",  # bad pixel (was interpolated)
                             "S_MASK_SATUR",  # saturated
                             "S_MASK_CR",  # cosmic ray
                             "S_MASK_GHOST"]  # electronics artifacts

const BAND_CHAR_TO_NUM = Dict('u'=>1, 'g'=>2, 'r'=>3, 'i'=>4, 'z'=>5)

"""
interp_sky(data, xcoords, ycoords)

Interpolate the 2-d array `data` at the grid of array coordinates spanned
by the vectors `xcoords` and `ycoords` using bilinear interpolation.
The output array will have size `(length(xcoords), length(ycoords))`.
For example, if `x[1] = 3.3` and `y[2] = 4.7`, the element `out[1, 2]`
will be a result of linear interpolation between the values
`data[3:4, 4:5]`.

For coordinates that are out-of-bounds (e.g., `xcoords[i] < 1.0` or
`xcoords[i] > size(data,1)` where the interpolation would index data values
outside the array, the nearest values in the data array are used. (This is
constant extrapolation.)
"""
function interp_sky{T, S}(data::Array{T, 2}, xcoords::Vector{S},
                          ycoords::Vector{S})
    nx, ny = size(data)
    out = Array(T, length(xcoords), length(ycoords))
    for j=1:length(ycoords)
        y0 = floor(Int, ycoords[j])
        y1 = y0 + 1
        yw0 = ycoords[j] - y0
        yw1 = one(S) - yw0

        # modify out-of-bounds indicies to 1 or ny
        y0 = min(max(y0, 1), ny)
        y1 = min(max(y1, 1), ny)

        for i=1:length(xcoords)
            x0 = floor(Int, xcoords[i])
            x1 = x0 + 1
            xw0 = xcoords[i] - x0
            xw1 = one(S) - xw0

            # modify out-of-bounds indicies to 1 or nx
            x0 = min(max(x0, 1), nx)
            x1 = min(max(x1, 1), nx)
            @inbounds out[i, j] = (xw0 * yw0 * data[x0, y0] +
                                   xw1 * yw0 * data[x1, y0] +
                                   xw0 * yw1 * data[x0, y1] +
                                   xw1 * yw1 * data[x1, y1])
        end
    end
    return out
end

"""
read_sky(hdu)

Construct an image of the sky from the sky HDU of a SDSS \"frame\" file.
This is typically the 3rd HDU of the file.
"""
function read_sky(hdu::FITSIO.TableHDU)
    sky_small = squeeze(read(hdu, "ALLSKY"), 3)::Array{Float32, 2}
    sky_x = vec(read(hdu, "XINTERP"))::Vector{Float32}
    sky_y = vec(read(hdu, "YINTERP"))::Vector{Float32}

    # convert sky interpolation coordinates from 0-indexed to 1-indexed
    for i=1:length(sky_x)
        sky_x[i] += 1.0f0
    end
    for i=1:length(sky_y)
        sky_y[i] += 1.0f0
    end

    # interpolate to full sky image
    return interp_sky(sky_small, sky_x, sky_y)
end

"""
read_frame(fname)

Read an SDSS \"frame\" FITS file and return a 4-tuple:

- `image`: sky-subtracted and calibrated 2-d image.
- `calibration`: 1-d \"calibration\".
- `sky`: 2-d sky that was subtracted.
- `wcs`: WCSTransform constructed from image header.
"""
function read_frame(fname)
    f = FITSIO.FITS(fname)
    hdr = FITSIO.read_header(f[1], ASCIIString)::ASCIIString
    image = read(f[1])::Array{Float32, 2}  # sky-subtracted & calibrated data
    calibration = read(f[2])::Vector{Float32}
    sky = read_sky(f[3])
    close(f)

    wcs = WCS.from_header(hdr)[1]

    return image, calibration, sky, wcs
end

"""
decalibrate!(image, sky, calibration, gain)

Convert `image` to raw counts. `image` is modified in-place.

Note that result may not be close to integers due to the analog to digital
conversion process in the instrument.
"""
function decalibrate!{T<:Number}(image::Array{T, 2}, sky::Array{T, 2},
                                 calibration::Vector{T}, gain::Number)
    @assert size(image) == size(sky)
    @assert size(image, 1) == length(calibration)
    for j=1:size(image, 2), i=1:size(image, 1)
        @inbounds image[i, j] = gain * (image[i, j] / calibration[i] +
                                        sky[i, j])
    end
end

"""
read_field_gains(fname, fieldnum)

Return the image gains for field number `fieldnum` in an SDSS
\"photoField\" file `fname`.
"""
function read_field_gains(fname, fieldnum::Integer)

    f = FITSIO.FITS(fname)
    fieldnums = read(f[2], "FIELD")::Vector{Int32}
    gains = read(f[2], "GAIN")::Array{Float32, 2}
    close(f)

    # Find first occurance of `fieldnum` and return the corresponding gain.
    for i=1:length(fieldnums)
        if fieldnums[i] == fieldnum
            bandgains = gains[:, i]
            return Dict(zip("ugriz", bandgains))
        end
    end

    error("field number $fieldnum not found in file: $fname")
end

"""
read_mask(fname[, mask_planes])

Read a \"fpM\"-format SDSS file and return masked image ranges,
based on `mask_planes`. Returns two `Vector{UnitRange{Int}}`,
giving the range of x and y indicies to be masked.
"""
function read_mask(fname, mask_planes=DEFAULT_MASK_PLANES)
    f = FITSIO.FITS(fname)

    # The last (12th) HDU contains a key describing what each of the
    # other HDUs are. Use this to find the indicies of all the relevant
    # HDUs (those with attributeName matching a value in `mask_planes`).
    value = read(f[12], "Value")::Vector{Int32}
    def = read(f[12], "defName")::Vector{ASCIIString}
    attribute = read(f[12], "attributeName")::Vector{ASCIIString}

    # initialize return values
    xranges = UnitRange{Int}[]
    yranges = UnitRange{Int}[]

    # Loop over keys and check if each is a mask plane we're interested in
    # (those with defName == "S_MASKTYPE" and "attributeName" in mask_planes).
    # If so, read from the corresponding HDU and construct ranges to mask.
    for i=1:length(value)
        if (def[i] == "S_MASKTYPE" && attribute[i] in mask_planes)

            # `value` starts from 0, but first table hdu is hdu number 2
            hdunum = value[i] + 2

            cmin = read(f[hdunum], "cmin")::Vector{Int32}
            cmax = read(f[hdunum], "cmax")::Vector{Int32}
            rmin = read(f[hdunum], "rmin")::Vector{Int32}
            rmax = read(f[hdunum], "rmax")::Vector{Int32}

            # "c" ("column") refers to mask's NAXIS1 (x axis);
            # "r" ("row") refers to mask's NAXIS2 (y axis).
            # cmin/cmax and rmin/rmax are 0-based and inclusive, so we add 1
            # to both.
            for j=1:length(cmin)
                push!(xranges, (cmin[j] + 1):(cmax[j] + 1))
                push!(yranges, (rmin[j] + 1):(rmax[j] + 1))
            end
        end
    end

    return xranges, yranges
end

# -----------------------------------------------------------------------------
# PSF-related functions

"""
SDSSPSF

SDSS representation of a spatially variable PSF. The PSF is represented as
a weighted combination of eigenimages (stored in `rrows`), where the weights
vary smoothly across the image as a polynomial of the form

```
weight[k](x, y) = sum_{i,j} cmat[i, j, k] * (rcs * x)^i (rcs * y)^j
```

where `rcs` is a coordinate transformation and `x` and `y` are zero-indexed.
"""
immutable SDSSPSF
    rrows::Array{Float64,2}  # A matrix of flattened eigenimages.
    rnrow::Int  # The number of rows in an eigenimage.
    rncol::Int  # The number of columns in an eigenimage.
    cmat::Array{Float64,3}  # The coefficients of the weight polynomial

    function SDSSPSF(rrows::Array{Float64, 2}, rnrow::Integer, rncol::Integer,
                     cmat::Array{Float64, 3})
        # rrows contains eigen images. Each eigen image is along the first
        # dimension in a flattened form. Check that dimensions match up.
        @assert size(rrows, 1) == rnrow * rncol

        # The second dimension is the number of eigen images, which should
        # match the number of coefficient arrays.
        @assert size(rrows, 2) == size(cmat, 3)

        return new(rrows, Int(rnrow), Int(rncol), cmat)
    end
end

"""
psf(x, y)

Evaluate the PSF at the given image coordinates. The size of the result is
will be `(psf.rnrow, psf.rncol)`, with the PSF (presumably) centered in the
stamp.

This function was originally based on the function sdss_psf_at_points
in astrometry.net:
https://github.com/dstndstn/astrometry.net/blob/master/util/sdss_psf.py
"""
function call(psf::SDSSPSF, x::Real, y::Real)
    const RCS = 0.001  # A coordinate transform to keep polynomial
                       # coefficients to a reasonable size.
    nk = size(psf.rrows, 2)  # number of eigen images.

    # initialize output stamp
    stamp = zeros(psf.rnrow, psf.rncol)

    # Loop over eigen images
    for k=1:nk
        # calculate the weight for the k-th eigen image from psf.cmat.
        # Note that the image coordinates and coefficients are intended
        # to be zero-indexed.
        w = 0.0
        for j=1:size(psf.cmat, 2), i=1:size(psf.cmat, 1)
            w += (psf.cmat[i, j, k] *
                  (RCS * (x - 1.0))^(i-1) * (RCS * (y - 1.0))^(j-1))
        end

        # add the weighted k-th eigen image to the output stamp
        for i=1:length(stamp)
            stamp[i] += w * psf.rrows[i, k]
        end
    end

    return stamp
end

"""
read_psf(fitsfile, band)

Read PSF components for the given band ('u', 'g', 'r', 'i', or 'z')
from the open FITSIO.FITS instance `fitsfile`, which should be a SDSS
\"psField\" file, and return a SDSSPSF instance representing the spatially
variable PSF.
"""
function read_psf(fitsfile::FITSIO.FITS, band::Char)

    # get the HDU, assuming PSF hdus are in order after primary header
    extnum = 1 + BAND_CHAR_TO_NUM[band]
    hdu = fitsfile[extnum]::FITSIO.TableHDU

    nrows = FITSIO.read_key(hdu, "NAXIS2")[1]::Int
    nrow_b = (read(hdu, "nrow_b")::Vector{Int32})[1]
    ncol_b = (read(hdu, "ncol_b")::Vector{Int32})[1]
    rnrow = (read(hdu, "rnrow")::Vector{Int32})[1]
    rncol = (read(hdu, "rncol")::Vector{Int32})[1]
    cmat_raw = read(hdu, "c")::Array{Float32, 3}
    rrows_raw = read(hdu, "rrows")::Array{Array{Float32,1},1}

    # Only the first (nrow_b, ncol_b) submatrix of cmat is used for
    # reasons obscure to the author.
    cmat = Array(Float64, nrow_b, ncol_b, size(cmat_raw, 3))
    for k=1:size(cmat_raw, 3), j=1:nrow_b, i=1:ncol_b
        cmat[i, j, k] = cmat_raw[i, j, k]
    end

    # convert rrows to Array{Float64, 2}, assuming each row is the same length.
    rrows = Array(Float64, length(rrows_raw[1]), length(rrows_raw))
    for i=1:length(rrows_raw)
        rrows[:, i] = rrows_raw[i]
    end

    return SDSSPSF(rrows, rnrow, rncol, cmat)
end

"""
read_photoobj(fname)

Read a source catalog from an SDSS \"photoObj\" FITS file for a given
run/camcol/field combination, returning a Vector of columns and a Vector
of column names.

This is currently pretty specific to Celeste because it returns only columns
that we're interested in, and it does some transformations of some of the
columns.
"""
function read_photoobj(fname, band::Char='r')

    bandnum = BAND_CHAR_TO_NUM[band]

    f = FITSIO.FITS(fname)
    hdu = f[2]::FITSIO.TableHDU

    objid = read(hdu, "objid")::Vector{ASCIIString}
    ra = read(hdu, "ra")::Vector{Float64}
    dec = read(hdu, "dec")::Vector{Float64}
    mode = read(hdu, "mode")::Vector{UInt8}
    thing_id = read(hdu, "thing_id")::Vector{Int32}

    # Get "bright" objects.
    # (In objc_flags, the bit pattern Int32(2) corresponds to bright objects.)
    is_bright = read(hdu, "objc_flags")::Vector{Int32} & Int32(2) .!= 0
    has_child = read(hdu, "nchild")::Vector{Int16} .> 0

    # 6 = star, 3 = galaxy, others can be ignored.
    objc_type = read(hdu, "objc_type")::Vector{Int32}
    is_star = objc_type .== 6
    is_gal = objc_type .== 3
    is_bad_obj = !(is_star | is_gal)

    fracdev = vec((read(hdu, "fracdev")::Matrix{Float32})[bandnum, :])

    # determine mask for rows
    has_dev = fracdev .> 0.
    has_exp = fracdev .< 1.
    is_comp = has_dev & has_exp
    is_bad_fracdev = (fracdev .< 0.) | (fracdev .> 1)
    mask = !(is_bad_fracdev | is_bad_obj | is_bright | has_child)

    # Read the fluxes.
    # Record the cmodelflux if the galaxy is composite, otherwise use
    # the flux for the appropriate type.
    psfflux = read(hdu, "psfflux")::Matrix{Float32}
    cmodelflux = read(hdu, "cmodelflux")::Matrix{Float32}
    devflux = read(hdu, "devflux")::Matrix{Float32}
    expflux = read(hdu, "expflux")::Matrix{Float32}

    # We actually only store the following properties for the given band
    # (see below in Dict construction) NB: the phi quantites in tractor are
    # multiplied by -1.
    phi_dev_deg = read(hdu, "phi_dev_deg")::Matrix{Float32}
    phi_exp_deg = read(hdu, "phi_exp_deg")::Matrix{Float32}

    theta_dev = read(hdu, "theta_dev")::Matrix{Float32}
    theta_exp = read(hdu, "theta_exp")::Matrix{Float32}

    ab_exp = read(hdu, "ab_exp")::Matrix{Float32}
    ab_dev = read(hdu, "ab_dev")::Matrix{Float32}

    close(f)

    # construct result catalog
    catalog = Dict("objid"=>objid[mask],
                   "ra"=>ra[mask],
                   "dec"=>dec[mask],
                   "thing_id"=>thing_id[mask],
                   "mode"=>mode[mask],
                   "is_star"=>is_star[mask],
                   "is_gal"=>is_gal[mask],
                   "frac_dev"=>fracdev[mask],
                   "ab_exp"=>vec(ab_exp[bandnum, mask]),
                   "theta_exp"=>vec(theta_exp[bandnum, mask]),
                   "phi_exp"=>vec(phi_exp_deg[bandnum, mask]),
                   "ab_dev"=>vec(ab_dev[bandnum, mask]),
                   "theta_dev"=>vec(theta_dev[bandnum, mask]),
                   "phi_dev"=>vec(phi_dev_deg[bandnum, mask]))
    for (b, n) in BAND_CHAR_TO_NUM
        catalog["psfflux_$b"] = vec(psfflux[n, mask])
        catalog["compflux_$b"] = vec(cmodelflux[n, mask])
        catalog["expflux_$b"] = vec(expflux[n, mask])
        catalog["devflux_$b"] = vec(devflux[n, mask])
    end

    # left over note about fluxes above: Use the cmodelflux if the
    # galaxy is composite, otherwise use the flux for the appropriate type.

    return catalog
end


text:
module RTE
using CUDA
using Adapt
using ..AtmosphericStates
using DocStringExtensions
import GaussQuadrature
using ..Sources
using ..Fluxes
using ..Optics
using ..BCs

import ..Parameters as RP

export Solver

"""
    struct Solver{
        FT,
        I,
        FTA1D,
        FTA2D,
        AS,
        OP,
        SL,
        SS,
        BCL,
        BCS,
        FXBL,
        FXBS,
        FXL,
        FXS,
    }

The high-level RRTMGP data structure storing the atmospheric state,
optical properties, sources, boundary conditions and fluxes configurations
for a given simulation.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Solver{
    FT <: AbstractFloat,
    I <: Int,
    FTA1D <: AbstractArray{FT, 1},
    FTA2D <: AbstractArray{FT, 2},
    AS <: AbstractAtmosphericState{FT, I, FTA1D},
    OP <: AbstractOpticalProps{FT, FTA2D},
    SL <: Union{AbstractSourceLW{FT, FTA1D, FTA2D}, Nothing},
    SS <: Union{SourceSW2Str{FT, FTA1D, FTA2D}, Nothing},
    BCL <: Union{LwBCs{FT}, Nothing},
    BCS <: Union{SwBCs{FT}, Nothing},
    FXBL <: Union{FluxLW{FT, FTA2D}, Nothing},
    FXBS <: Union{FluxSW{FT, FTA2D}, Nothing},
    FXL <: Union{FluxLW{FT, FTA2D}, Nothing},
    FXS <: Union{FluxSW{FT, FTA2D}, Nothing},
}
    as::AS         # atmospheric state
    op::OP         # optical properties
    src_lw::SL     # source functions
    src_sw::SS     # source functions
    bcs_lw::BCL    # boundary conditions
    bcs_sw::BCS    # boundary conditions
    fluxb_lw::FXBL # temporary storage for bandwise calculations
    fluxb_sw::FXBS # temporary storage for bandwise calculations
    flux_lw::FXL   # fluxes for longwave problem
    flux_sw::FXS   # fluxes for shortwave problem
end

Solver(as, op, src_lw, src_sw, bcs_lw, bcs_sw, fluxb_lw, fluxb_sw, flux_lw, flux_sw) = Solver{
    eltype(as.p_lev),
    typeof(as.ncol),
    typeof(as.t_sfc),
    typeof(as.p_lev),
    typeof(as),
    typeof(op),
    typeof(src_lw),
    typeof(src_sw),
    typeof(bcs_lw),
    typeof(bcs_sw),
    typeof(fluxb_lw),
    typeof(fluxb_sw),
    typeof(flux_lw),
    typeof(flux_sw),
}(
    as,
    op,
    src_lw,
    src_sw,
    bcs_lw,
    bcs_sw,
    fluxb_lw,
    fluxb_sw,
    flux_lw,
    flux_sw,
)

Adapt.@adapt_structure Solver
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
using Statistics, Dates

function test_ach()

	@eval using Plots
	gr()

	nn = 1
	ts = 10
	ACh = fill(0.0, nn)
	synt = 3 # or 8
	Φ = 0.5
	spikesequence = repeat([1, 0, 0, 0, 0], outer=2)
	@show spikesequence
	spt = sim_spikes(nn, ts, spikesequence)
	p = scatter()

	for t = 1:ts

		update_ACh!(ACh, synt, Φ, t, spt[t])
		scatter!(p, [t], [ACh[1]])

	end

	display(p)
end

function sim_spikes(nn, ts, spikesequence)

	lastspike = 1
	out = []

	for i = 1 : ts

		if spikesequence[i] == 1
			push!(out, i)
			lastspike = i
		else
			push!(out, lastspike)
		end

	end
	return out
end

function test_da()

	@eval using Plots
	gr()

	ts = 40
	τ = 20
	da = 0
	BA = repeat([0.5, 0.0, 0.0, 0.0],inner=8, outer=5)
	p = scatter()

	for t = 1:ts
		da = update_da(da, BA[t], τ; δt=1)
		@show da, BA[t]
		scatter!(p, [t], [da], show=true, legend=false)
	end
end

function test_transmission()

	@eval using Plots
	gr()

	# get2 layers, with one synapse between them, and see how activity in the 1st translates to activity in the second.
	nn = [1, 1]
	ts = 60
	ACh = [0.0]
	w = 4.0
	spt1 = sim_spikes(1, 60, repeat([1, 1, 0],inner=20))
	spt2 = 0
	synt = 3 # or 8
	Φ = 0.5
	v = [-60.0]
	rec = [0.0]
	sp = [false]
	C = 100
	a = 0.3
	b = -0.3
	c = -65
	d = 8
	k = 2

	p = scatter()

	for t = 1:ts

		update_ACh!(ACh, synt, Φ, t, spt1[t])
		# calc output.
		input = (-v) .* w .* ACh
		update_activation!(1, v, -60, sp, spt2, t, -40, rec, input, C, a, b, c, d, k;)
		scatter!(p,[t], v, legend=false, color=RGB(0.5, 0.5, 1))

	end

	display(p)
end

function test_weight()

	@eval using Plots
	gr()

	# get2 layers, with one synapse between them, and see how activity in the 1st translates to activity in the second.
	nn = [1, 1]
	ts = 60
	ACh = hcat(0.0)
	w = hcat(4.0)
	γ = hcat(0.0)
	connected = [1]
	spt1 = sim_spikes(1, 60, repeat([1, 1, 0],inner=20))
	spt2 = hcat(0)
	synt = 3 # or 8
	Φ = 0.5
	v = hcat(-60.0)
	rec = hcat(0.0)
	sp = hcat(false)
	C = 100
	a = 0.3
	b = -0.3
	c = -65
	d = 8
	k = 2

	p = scatter()

	for t = 1:ts
		update_ACh!(ACh, synt, Φ, t, spt1[t])
		input = (-v) .* w .* ACh
		update_activation!(1, v, -60, sp, spt2, t, -40, rec, input, C, a, b, c, d, k;)
		update_weights!(w, γ, connected, t, spt1[t], spt2, 1, 40; δt=1, A₋=-1, t₋=15)

		scatter!(p, [t], w, legend=false, color=RGB(0.5, 0.5, 1))
	end
	display(p)
end

function response_before_after_learning2()

	nn = [100, 1000, 5]
	sensory = constructinputsequence((100,), (SparseInput,), stages=[10,100,1], input_bool=Bool[1,1,0], da_bool=Bool[0,1,1])
	numsteps = duration(sensory)
	weights, synapses, trainspikes = train_model(sensory, nn, numsteps, reportvar="spiked")
	_, _, testspikes = test_model(sensory, nn, numsteps, weights, synapses, reportvar="spiked")
	train_plot = plotmeanlayer(trainspikes)
	test_plot = plotmeanlayer(testspikes)
	plot(train_plot,test_plot, layout=(1,2))
end

function many_stimuli()

	nn = [100, 1000, 5]
	numtrain = 4
	numtest = 4

	facenum = tuple(nn[1])
	stimtype = (SparseInput,)
	sstages = [10, 100 , 10]
	input = Bool[1,1,0]
	dastages = Bool[0,1,0]
	reportvar = "spiked"

	p = get_parameters()
	m = MatrixTypes(initialise_matrices(nn, p)...)
	da = 0

	sensory = [constructinputsequence(facenum, stimtype, stages=sstages, input_bool=input, da_bool=dastages) for te in 1:(numtrain + numtest)]

	# Training Period
	for tr in 1:numtrain
		numsteps = duration(sensory[tr])
		reporter = run_all_steps(nn, numsteps, m, p, sensory[tr], da, savevars=nothing, update=true, reportvar=reportvar)
		m = reset(nn, p, m.weights, m.synapses)
	end

	#
	for te in 1:(numtrain + numtest)
		numsteps = duration(sensory[te])
		m = reset(nn, p, m.weights, m.synapses)
		reporter = run_all_steps(nn, numsteps, m, p, sensory[te], da, savevars=nothing, update=false, reportvar=reportvar)
		println(average_layers(reporter,3))
	end
end

function many_observations()
	nn = [100, 100, 1]
	kcn = [100, 1000, 10_000]
	p = []
	for kc in kcn
		nn[2] = kc
		println(nn)
		push!(p,observe_variable("spiked",nn=nn))
	end
	plot(p[1],p[2],p[3])
end

function observe_variable(reportvar::String; nn=[100,1000,1])
	numstim = 1

	facenum = tuple(nn[1])
	stimtype = (SparseInput,)
	sstages = [10, 500 , 10]
	input = Bool[1,1,0]
	dastages = Bool[0,1,0]

	p = get_parameters()
	m = MatrixTypes(initialise_matrices(nn,p)...)
	da = 0

	sensory = [constructinputsequence(facenum, stimtype, stages=sstages, input_bool=input, da_bool=dastages) for te in 1:numstim]
	for i in 1:numstim
		numsteps = duration(sensory[i])
		normalise_layer!(m, p, l=2)
		reporter = run_all_steps(nn, numsteps, m, p, sensory[i], da, savevars=nothing, update=false, reportvar=reportvar)
		p = plot(average_timestep(reporter))
	end
	return p
end

function observe_variable_in_layer(reportvar::String; nn=[100,1000,1], layer=1)

	numstim = 1

	facenum = tuple(nn[1])
	stimtype = (SparseInput,)
	sstages = [10, 500 , 10]
	input = Bool[1,1,0]
	dastages = Bool[0,1,0]

	p = get_parameters()
	m = MatrixTypes(initialise_matrices(nn,p)...)
	da = 0

	sensory = [constructinputsequence(facenum, stimtype, stages=sstages, input_bool=input, da_bool=dastages) for te in 1:numstim]
	for i in 1:numstim
		numsteps = duration(sensory[i])
		normalise_layer!(m, p, l=2)
		reporter = run_all_steps(nn, numsteps, m, p, sensory[i], da, savevars=nothing, update=false, reportvar=reportvar)
		p = plot(average_timestep(reporter,l=layer))
	end
	return p
end

function average_layers(A::Array{<:BrainTypes,1},l)
	holder = Vector{Float64}(undef,length(A))
	for t in 1:length(A)
		holder[t] = mean(A[t].layers[l])
	end
	return mean(holder)
end

function average_timestep(A::Array{<:BrainTypes,1};l=nothing)
	numsteps = length(A)
	numlayers = isnothing(l) ? length(A[1].layers) : l
	out = zeros(numsteps, numlayers)
	for s = 1:numsteps
		for l = 1:numlayers
			out[s,l] = mean(A[s].layers[l])
		end
	end
	return out
end

# this is actually just a less efficient version of initialise_matrices(_,_,weights,synapses)
function reset(nn, p, weights, synapses)

	return MatrixTypes(initialise_matrices(nn, p, weights, synapses)...)
end

function test_many_kenyon_weights()

	totalweights = range(300,500,length=3)
	p = get_parameters()

	numtest=30
	numtrain=30
	nn = [100, 1000, 1]

	for kw in totalweights

		p.weight_target = (kw, kw, kw)
		MBON = test_parameter(p, nn, numtest=numtest, numtrain=numtrain)
		print("\nTrain Average: $(mean(MBON[1:numtest])), Test Average: $(mean(MBON[numtest+1:end]))")

	end
end

function test_parameter(p, nn; numtrain=4, numtest=4, layer_view=3)

	facenum = tuple(nn[1])
	stimtype = (ColorInput,)
	sstages = [20, 50, 0]
	input = Bool[1,1, 0]
	dastages = Bool[0,1, 0]
	reportvar = "spiked"
	da=0

	m = MatrixTypes(initialise_matrices(nn, p)...)
	sensory = [constructinputsequence(facenum, stimtype, stages=sstages, input_bool=input, da_bool=dastages) for te in 1:(numtrain + numtest)]
	# Training Period
	for tr in 1:numtrain

		m = reset(nn, p, m.weights, m.synapses)
		#normalise_layer!(m, p, l=2)
		numsteps = duration(sensory[tr])
		reporter = run_all_steps(nn, numsteps, m, p, sensory[tr], da, savevars=nothing, update=true, normweights=true, reportvar=reportvar)
	end
	#normalise_layer!(m, p, l=2)

	avactivation = zeros(numtrain+numtest)

	for te in 1:(numtrain + numtest)
		numsteps = duration(sensory[te])
		m = reset(nn, p, m.weights, m.synapses)
		reporter = run_all_steps(nn, numsteps, m, p, sensory[te], da, savevars=nothing, update=false, reportvar=reportvar)
		avactivation[te] = average_layers(reporter,layer_view)
	end
	return avactivation
end

function test_learning(p, nn; numtest=16, numtrain=16, numreps=10, progress=false)

	learn_measure = zeros(numreps)

	for i in 1:numreps
		MBON = test_parameter(p, nn, numtest=numtest, numtrain=numtrain, layer_view=3)
		learn_measure[i] = mean(MBON[numtrain+1:end]) - mean(MBON[1:numtrain])
		progress && print("\nRepetition: $i")
	end
	return learn_measure
end

function test_many_parameters(parameter_name::String, parameter_values::Array)

	p = get_parameters()
	fname = "C:\\Users\\MIKKO\\.julia\\dev\\MushroomBody\\src\\output\\Parameter Search\\" * parameter_name * Dates.format(now(), "ddmm_HHMM") * ".txt"

	open(fname, "w") do io

		for (i, parameter) in enumerate(parameter_values)
			setproperty!(p, parameter_name, parameter)

			start = Dates.Time(Dates.now())

			learning_measures = test_learning(p, p.nn)

			finish = Dates.Time(Dates.now())

			print("\nEstimated Time to Completion: $(round((finish-start)*(length(parameter_values)-i),Dates.Second))\n")
			println(io, "$(parameter_name): $(parameter)  Mean Learning: $(mean(learning_measures))  Standard Deviation: $(std(learning_measures))")
		end
		println(io, "\nParameters\n")
		for param in string(p)
			println(io, param)
		end

	end
end

function kcactivation()

	p = get_parameters()
	numtest=4
	nn=[100,1000,1]
	kcnum = [100,1000,4000]
	kc_activation = zeros(length(kcnum))
	for (i, kc) in enumerate(kcnum)
		nn[2] = kc
		kc_activation[i] = mean(test_parameter(p, nn, numtest=numtest, numtrain=0, layer_view=2))
	end
	plot(kcnum, kc_activation)
end


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "helicopter_model.jl")
include(@relpath "helicopter_specifications.jl")

S = helicopter_model()
X0, options = helicopter_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with 1D blocks
partition = [[i] for i in 1:28]

# reachability algorithms
algorithm_reach = BFFPSV18(:δ => 5e-3, :vars => [1], :partition => partition,
                           :assume_homogeneous => true)
algorithm_check = BFFPSV18(:δ => 5e-3, :vars => [1], :partition => partition,
                           :assume_homogeneous => true)

# check property
options[:mode] = "check"
solution = solve(problem, options; op=algorithm_check)
@assert solution.satisfied

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to time and x₁
solution.options[:plot_vars] = [0, 1]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("helicopter")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using Random

Random.seed!(1974)
println("Seed 1974: ",rand(),"\t", rand(), "\t", rand())
Random.seed!(1975)
println("Seed 1975: ",rand(),"\t", rand(), "\t", rand())
Random.seed!(1974)
println("Seed 1974: ",rand(),"\t", rand(), "\t", rand())

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
	gmtspatial(cmd0::String="", arg1=nothing, kwargs...)

Geospatial operations on points, lines and polygons.

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

Parameters
----------

- **A** | **nn** | **nearest_neighbor** :: [Type => Str]     `Arg = [amin_dist][unit]`

    Perform spatial nearest neighbor (NN) analysis: Determine the nearest neighbor of each point
    and report the NN distances and the point IDs involved in each pair.
    ($(GMTdoc)gmtspatial.html#a)
- **C** | **clip** :: [Type => Bool]

    Clips polygons to the map region, including map boundary to the polygon as needed. The result is a closed polygon.
    ($(GMTdoc)gmtspatial.html#c)
- **D** | **duplicates** :: [Type => Str]   `Arg = [+ffile][+aamax][+ddmax][+c|Ccmax][+sfact]`

    Check for duplicates among the input lines or polygons, or, if file is given via +f, check if the
    input features already exist among the features in file.
    ($(GMTdoc)gmtspatial.html#d)
- **E** | **handedness** :: [Type => Str]  `Arg = +|-`

    Reset the handedness of all polygons to match the given + (counter-clockwise) or - (clockwise). Implies Q+
    ($(GMTdoc)gmtspatial.html#e)
- **F** | **force_polygons** :: [Type => Str | []]   `Arg = [l]`

    Force input data to become polygons on output, i.e., close them explicitly if not already closed.
    Optionally, append l to force line geometry.
    ($(GMTdoc)gmtspatial.html#f)
- **I** | **intersections** :: [Type => Str | []]   `Arg = [e|i]`

    Determine the intersection locations between all pairs of polygons.
    ($(GMTdoc)gmtspatial.html#i)
- **N** | **in_polyg** :: [Type => Str]     `Arg = pfile[+a][+pstart][+r][+z]`

    Determine if one (or all, with +a) points of each feature in the input data are inside any of
    the polygons given in the pfile.
    ($(GMTdoc)gmtspatial.html#n)
- **Q** | **area** or **length** :: [Type => Str]      `Arg = [[-|+]unit][+cmin[/max]][+h][+l][+p][+s[a|d]]`

    Measure the area of all polygons or length of line segments.
    ($(GMTdoc)gmtspatial.html#q)
- $(GMT.opt_R)
- **S** | **polyg_process** :: [Type => Int]     `Arg = h|i|j|s|u`

    Spatial processing of polygons.
    ($(GMTdoc)gmtspatial.html#s)
- **T** | **truncate** :: [Type => Str | []]     `Arg = [clippolygon]`

    Truncate polygons against the specified polygon given, possibly resulting in open polygons.
    ($(GMTdoc)gmtspatial.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtspatial(cmd0::String="", arg1=nothing; kwargs...)

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

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:R :V_params :b :d :e :g :h :i :o :yx])
	cmd = parse_these_opts(cmd, d, [[:A :nn :nearest_neighbor], [:C :clip], [:D :duplicates], [:E :handedness],
				[:F :force_polygons], [:I :intersections], [:N :in_polyg], [:Q :area :length],
				[:S :polyg_process], [:T :truncate]])

	common_grd(d, cmd0, cmd, "gmtspatial ", arg1)		# Finish build cmd and run it
end

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


text:
"""
    CircularBuffer{T}(n)

The CircularBuffer type implements a circular buffer of fixed capacity
where new items are pushed to the back of the list, overwriting values
in a circular fashion.

Allocate a buffer of elements of type `T` with maximum capacity `n`.
"""
mutable struct CircularBuffer{T} <: AbstractVector{T}
    capacity::Int
    first::Int
    length::Int
    buffer::Vector{T}

    CircularBuffer{T}(capacity::Int) where {T} = new{T}(capacity, 1, 0, Vector{T}(undef, capacity))
end

CircularBuffer(capacity) = CircularBuffer{Any}(capacity)

"""
    empty!(cb::CircularBuffer)

Reset the buffer.
"""
function Base.empty!(cb::CircularBuffer)
    cb.length = 0
    return cb
end

Base.@propagate_inbounds function _buffer_index_checked(cb::CircularBuffer, i::Int)
    @boundscheck if i < 1 || i > cb.length
        throw(BoundsError(cb, i))
    end
    _buffer_index(cb, i)
end

@inline function _buffer_index(cb::CircularBuffer, i::Int)
    n = cb.capacity
    idx = cb.first + i - 1
    return ifelse(idx > n, idx - n, idx)
end

"""
    cb[i]

Get the i-th element of CircularBuffer.

* `cb[1]` to get the element at the front
* `cb[end]` to get the element at the back
"""
@inline Base.@propagate_inbounds function Base.getindex(cb::CircularBuffer, i::Int)
    cb.buffer[_buffer_index_checked(cb, i)]
end

"""
    cb[i] = data

Store data to the `i`-th element of `CircularBuffer`.
"""
@inline Base.@propagate_inbounds function Base.setindex!(cb::CircularBuffer, data, i::Int)
    cb.buffer[_buffer_index_checked(cb, i)] = data
    return cb
end

"""
    pop!(cb::CircularBuffer)

Remove the element at the back.
"""
@inline function Base.pop!(cb::CircularBuffer)
    @boundscheck (cb.length == 0) && throw(ArgumentError("array must be non-empty"))
    i = _buffer_index(cb, cb.length)
    cb.length -= 1
    return @inbounds cb.buffer[i]
end

"""
    push!(cb::CircularBuffer, data)

Add an element to the back and overwrite front if full.
"""
@inline function Base.push!(cb::CircularBuffer, data)
    # if full, increment and overwrite, otherwise push
    if cb.length == cb.capacity
        cb.first = (cb.first == cb.capacity ? 1 : cb.first + 1)
    else
        cb.length += 1
    end
    @inbounds cb.buffer[_buffer_index(cb, cb.length)] = data
    return cb
end

"""
    popfirst!(cb::CircularBuffer)

Remove the element from the front of the `CircularBuffer`.
"""
function popfirst!(cb::CircularBuffer)
    @boundscheck (cb.length == 0) && throw(ArgumentError("array must be non-empty"))
    i = cb.first
    cb.first = (cb.first + 1 > cb.capacity ? 1 : cb.first + 1)
    cb.length -= 1
    return @inbounds cb.buffer[i]
end

"""
    pushfirst!(cb::CircularBuffer, data)

Insert one or more items at the beginning of CircularBuffer
and overwrite back if full.
"""
function pushfirst!(cb::CircularBuffer, data)
    # if full, decrement and overwrite, otherwise pushfirst
    cb.first = (cb.first == 1 ? cb.capacity : cb.first - 1)
    if length(cb) < cb.capacity
        cb.length += 1
    end
    @inbounds cb.buffer[cb.first] = data
    return cb
end

"""
    append!(cb::CircularBuffer, datavec::AbstractVector)

Push at most last `capacity` items.
"""
function Base.append!(cb::CircularBuffer, datavec::AbstractVector)
    # push at most last `capacity` items
    n = length(datavec)
    for i in max(1, n-capacity(cb)+1):n
        push!(cb, datavec[i])
    end
    return cb
end

"""
    fill!(cb::CircularBuffer, data)

Grows the buffer up-to capacity, and fills it entirely.
It doesn't overwrite existing elements.
"""
function Base.fill!(cb::CircularBuffer, data)
    for i in 1:capacity(cb)-length(cb)
        push!(cb, data)
    end
    return cb
end

"""
    length(cb::CircularBuffer)

Return the number of elements currently in the buffer.
"""
Base.length(cb::CircularBuffer) = cb.length

Base.eltype(::Type{CircularBuffer{T}}) where T = T

"""
    size(cb::CircularBuffer)

Return a tuple with the size of the buffer.
"""
Base.size(cb::CircularBuffer) = (length(cb),)

Base.convert(::Type{Array}, cb::CircularBuffer{T}) where {T} = T[x for x in cb]

"""
    isempty(cb::CircularBuffer)

Test whether the buffer is empty.
"""
Base.isempty(cb::CircularBuffer) = cb.length == 0

""""
    capacity(cb::CircularBuffer)

Return capacity of CircularBuffer.
"""
capacity(cb::CircularBuffer) = cb.capacity

"""
    isfull(cb::CircularBuffer)

Test whether the buffer is full.
"""
isfull(cb::CircularBuffer) = length(cb) == cb.capacity

"""
    first(cb::CircularBuffer)

Get the first element of CircularBuffer.
"""
Base.@propagate_inbounds function Base.first(cb::CircularBuffer)
    @boundscheck (cb.length == 0) && throw(BoundsError(cb, 1))
    return cb.buffer[cb.first]
end
"""
    last(cb::CircularBuffer)

Get the last element of CircularBuffer.
"""
Base.@propagate_inbounds function Base.last(cb::CircularBuffer)
    @boundscheck (cb.length == 0) && throw(BoundsError(cb, 1))
    return cb.buffer[_buffer_index(cb, cb.length)]
end


text:
merge_ext!(S::NodalData, Ω::Int64, rest::Array{Int64, 1}) = nothing

merge!(S::NodalData) = error("NodalData cannot be merged!")


text:
#=

  Euler #38 in Julia.

  """
  Take the number 192 and multiply it by each of 1, 2, and 3:

      192 × 1 = 192
      192 × 2 = 384
      192 × 3 = 576

  By concatenating each product we get the 1 to 9 pandigital,
  192384576. We will call 192384576 the concatenated product of 192
  and (1,2,3)

  The same can be achieved by starting with 9 and multiplying by
  1, 2, 3, 4, and 5, giving the pandigital, 918273645, which is the
  concatenated product of 9 and (1,2,3,4,5).

  What is the largest 1 to 9 pandigital 9-digit number that can be
  formed as the concatenated product of an integer with
  (1,2, ... , n) where n > 1?
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00016527s
function euler38a()
    for n in 9876:-1:9
        println
        s = string(n)
        i = 2
        while length(s) < 9
            s *= string(n*i)
            i += 1
        end
        if length(s) == 9 && is_pandigital(s)
            return parse(Int,s)
        end
    end
    return nothing
end

run_euler(euler38a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@doc raw"""
    onehot(x::Float64, n::Int)
    onehot(x::Vector{Float64}, n::Int)
"""
function onehot(x::Float64, n::Int)
    y = zeros(Bool, n)
    k = trunc(Int, clamp(x, 0.0, 1.0) * n)
    y[k] = true
    return y
end

function onehot(x::Vector{Float64}, n::Int)
    m = length(x)
    y = zeros(Bool, m * n)
    for i = 1:m
        j = (i - 1) * n + 1
        k = trunc(Int, clamp(x[i], 0.0, 1.0) * n) - 1
        y[j + k] = true
    end
    return y
end


text:
function demoHMMinferenceSimple()

    V = 5  # number of visible states
    H = 10 # number of hidden states
    T = 20 # number of timesteps

    h = zeros(Integer,1,T) # holds the state value for the hidden variable at a specific timestep
    v = zeros(Integer,1,T) # holds the state value for the visible variable at a specific timestep

    ph1=condp(rand(H,1)) # probabilities for the states of the hidden variable at timestep 1

    phghm=condp(eye(H,H)) # transition matrix with sum(phghm, 1) = 1 with phghm[i, j] = p(hg=i | hm=j)

    # shuffle the column data in phghm while keeping sum(phghm, 1) = 1
    # below no 2 columns have value 1.0 on the same row
    phghmtmp=condp(eye(H,H))
    r=randperm(H)
    for i=1:H
        phghm[:,i]=phghmtmp[:,r[i]]
    end

    pvgh=condp((rand(V,H))) # emision matrix with sum(pvgh, 1) = 1 with pvgh[i, j] = p(vg = i | h = j)

    h[1]=randgen(ph1) # initialize the hidden variable h(t=1) with a random state based on ph1 distribution
    v[1]=randgen(pvgh[:,h[1]]) # initialize the visible variable v(t=1) with a random state based on pvgh(vg | h)

    for t=2:T
        h[t]=randgen(phghm[:,h[t-1]]) # set the hidden variable state h(t) based on h(t-1) using the transition matrix
        v[t]=randgen(pvgh[:,h[t]]) # set the visible variable state v(t) based on h(t) using the emission matrix
    end

    alpha,loglik=HMMforward(v,phghm,ph1,pvgh) # filtering

    gamma = HMMsmooth(v,phghm,ph1,pvgh,alpha) # smoothing

    maxstate,logprob = HMMviterbi(v,phghm,ph1,pvgh)

    println("Inference log likelihood = $loglik\n")
    println("most likely path (viterbi):\n")
    println(maxstate)
    println("original path (hidden states):")
    println(h)
    println("original path (visible states):")
    println(v)

    PyPlot.figure()
    ax = PyPlot.subplot(2, 1, 1)
    ax[:set_title]("filtering")
    PyPlot.pcolor(alpha)
    ax = PyPlot.subplot(2, 1, 2)
    ax[:set_title]("smoothing")
    PyPlot.pcolor(gamma)

end


text:
# These routines interpolate only within the grid, from cell centers to cell
# edges.

# Interpolations in one direction only, C -> E or E -> C

"""
    grid_interpolate!(out:ScalarGridData,in::ScalarGridData) -> ScalarGridData

Return the 1-d symmetric interpolation of scalar grid data `in` onto scalar grid
data `out`. Either `in` or `out` must be edge component data and the other must
be node data. The direction of interpolation is determined by the relationship
of `in` and `out`. For example, if `in` is dual nodes (cell by cell) and `out` is
primal x-edge components (cell by edge), then the interpolation takes place in the
y direction, since they are different types in this direction. (In the x direction,
they are of the same type (cell), so there is no interpolation in that direction.)
"""
function grid_interpolate! end

# Nodes to edge components

function grid_interpolate!(qu::XEdges{Dual, NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    # E x C <- C x C
    @inbounds for y in 2:NY-1, x in 1:NX-1
        qu[x,y] = (w[x,y] + w[x+1,y])/2
    end
    qu
end

function grid_interpolate!(qv::YEdges{Dual, NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    # C x E <- C x C
    @inbounds for y in 1:NY-1, x in 2:NX-1
        qv[x,y] = (w[x,y] + w[x,y+1])/2
    end
    qv
end

function grid_interpolate!(qu::XEdges{Primal, NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    # C x E <- C x C
    @inbounds for y in 1:NY-1, x in 2:NX-1
        qu[x,y] = (w[x,y] + w[x,y+1])/2
    end
    qu
end

function grid_interpolate!(qv::YEdges{Primal, NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
   # E x C <- C x C
    @inbounds for y in 2:NY-1, x in 1:NX-1
        qv[x,y] = (w[x,y] + w[x+1,y])/2
    end
    qv
end

function grid_interpolate!(qu::XEdges{Primal, NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    # C x E <- E x E
    @inbounds for y in 1:NY-1, x in 2:NX-1
      qu[x,y] = (w[x-1,y] + w[x,y])/2
    end
    qu
end

function grid_interpolate!(qv::YEdges{Primal, NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    # E x C <- E x E
    @inbounds for y in 2:NY-1, x in 1:NX-1
      qv[x,y] = (w[x,y-1] + w[x,y])/2
    end
    qv
end

function grid_interpolate!(qu::XEdges{Dual, NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    # E x C <- E x E
    @inbounds for y in 2:NY-1, x in 1:NX-1
      qu[x,y] = (w[x,y-1] + w[x,y])/2
    end
    qu
end

function grid_interpolate!(qv::YEdges{Dual, NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    # C x E <- E x E
    @inbounds for y in 1:NY-1, x in 2:NX-1
      qv[x,y] = (w[x-1,y] + w[x,y])/2
    end
    qv
end

# Edge components to nodes

function grid_interpolate!(w::Nodes{Dual, NX, NY}, qu::XEdges{Primal,NX, NY}) where {NX, NY}
    # C x C <- C x E
    @inbounds for y in 2:NY-1, x in 1:NX
        w[x,y] = (qu[x,y-1] + qu[x,y])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Dual, NX, NY}, qv::YEdges{Primal,NX, NY}) where {NX, NY}
    # C x C <- E x C
    @inbounds for y in 1:NY, x in 2:NX-1
        w[x,y] = (qv[x-1,y] + qv[x,y])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Dual, NX, NY}, qu::XEdges{Dual,NX, NY}) where {NX, NY}
  # C x C <- E x C
  @inbounds for y in 1:NY, x in 2:NX-1
      w[x,y] = (qu[x-1,y] + qu[x,y])/2
  end
    w
end

function grid_interpolate!(w::Nodes{Dual, NX, NY}, qv::YEdges{Dual,NX, NY}) where {NX, NY}
    # C x C <- C x E
    @inbounds for y in 2:NY-1, x in 1:NX
      w[x,y] = (qv[x,y-1] + qv[x,y])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Primal, NX, NY}, qu::XEdges{Dual,NX, NY}) where {NX, NY}
    # E x E <- E x C
    @inbounds for y in 1:NY-1, x in 1:NX-1
      w[x,y] = (qu[x,y] + qu[x,y+1])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Primal, NX, NY}, qv::YEdges{Dual,NX, NY}) where {NX, NY}
    # E x E <- C x E
    @inbounds for y in 1:NY-1, x in 1:NX-1
      w[x,y] = (qv[x,y] + qv[x+1,y])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Primal, NX, NY}, qu::XEdges{Primal,NX, NY}) where {NX, NY}
    # E x E <- C x E
    @inbounds for y in 1:NY-1, x in 1:NX-1
      w[x,y] = (qu[x,y] + qu[x+1,y])/2
    end
    w
end

function grid_interpolate!(w::Nodes{Primal, NX, NY}, qv::YEdges{Primal,NX, NY}) where {NX, NY}
    # E x E <- E x C
    @inbounds for y in 1:NY-1, x in 1:NX-1
      w[x,y] = (qv[x,y] + qv[x,y+1])/2
    end
    w
end


text:
mutable struct GLBuffer{T} <: GPUArray{T, 1}
    id          ::GLuint
    size        ::NTuple{1, Int}
    buffertype  ::GLenum
    usage       ::GLenum
    context     ::GLContext

    function GLBuffer{T}(ptr::Ptr{T}, buff_length::Int, buffertype::GLenum, usage::GLenum) where T
        id = glGenBuffers()
        glBindBuffer(buffertype, id)
        # size of 0 can segfault it seems
        buff_length = buff_length == 0 ? 1 : buff_length
        glBufferData(buffertype, buff_length * sizeof(T), ptr, usage)
        glBindBuffer(buffertype, 0)

        obj = new(id, (buff_length,), buffertype, usage, current_context())
        finalizer(free, obj)
        obj
    end
end

bind(buffer::GLBuffer) = glBindBuffer(buffer.buffertype, buffer.id)
#used to reset buffer target
bind(buffer::GLBuffer, other_target) = glBindBuffer(buffer.buffertype, other_target)

function similar(x::GLBuffer{T}, buff_length::Int) where T
    GLBuffer{T}(Ptr{T}(C_NULL), buff_length, x.buffertype, x.usage)
end

cardinality(::GLBuffer{T}) where {T} = cardinality(T)

#Function to deal with any Immutable type with Real as Subtype
function GLBuffer(
        buffer::Union{Base.ReinterpretArray{T, 1}, DenseVector{T}};
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T <: GLArrayEltypes
    GLBuffer{T}(pointer(buffer), length(buffer), buffertype, usage)
end

function GLBuffer(
        buffer::DenseVector{T};
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T <: GLArrayEltypes
    GLBuffer{T}(pointer(buffer), length(buffer), buffertype, usage)
end

function GLBuffer(
        buffer::AbstractVector{T};
        kw_args...
    ) where T <: GLArrayEltypes
    GLBuffer(collect(buffer); kw_args...)
end

function GLBuffer{T}(
        buffer::AbstractVector;
        kw_args...
    ) where T <: GLArrayEltypes
    GLBuffer(convert(Vector{T}, buffer); kw_args...)
end

function GLBuffer(
        ::Type{T}, len::Int;
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T <: GLArrayEltypes
    GLBuffer{T}(Ptr{T}(C_NULL), len, buffertype, usage)
end

function indexbuffer(
        buffer::VectorTypes{T};
        usage::GLenum = GL_STATIC_DRAW
    ) where T <: GLArrayEltypes
    GLBuffer(buffer, buffertype = GL_ELEMENT_ARRAY_BUFFER, usage=usage)
end
# GPUArray interface
function gpu_data(b::GLBuffer{T}) where T
    data = Vector{T}(undef, length(b))
    bind(b)
    glGetBufferSubData(b.buffertype, 0, sizeof(data), data)
    bind(b, 0)
    data
end

# Resize buffer
function gpu_resize!(buffer::GLBuffer{T}, newdims::NTuple{1, Int}) where T
    #TODO make this safe!
    newlength = newdims[1]
    oldlen    = length(buffer)
    if oldlen > 0
        old_data = gpu_data(buffer)
    end
    bind(buffer)
    glBufferData(buffer.buffertype, newlength*sizeof(T), C_NULL, buffer.usage)
    bind(buffer, 0)
    buffer.size = newdims
    if oldlen>0
        max_len = min(length(old_data), newlength) #might also shrink
        buffer[1:max_len] = old_data[1:max_len]
    end
    #probably faster, but changes the buffer ID
    # newbuff     = similar(buffer, newdims...)
    # unsafe_copy!(buffer, 1, newbuff, 1, length(buffer))
    # buffer.id   = newbuff.id
    # buffer.size = newbuff.size
    nothing
end

function gpu_setindex!(b::GLBuffer{T}, value::Vector{T}, offset::Integer) where T
    multiplicator = sizeof(T)
    bind(b)
    glBufferSubData(b.buffertype, multiplicator*(offset-1), sizeof(value), value)
    bind(b, 0)
end

function gpu_setindex!(b::GLBuffer{T}, value::Vector{T}, offset::UnitRange{Int}) where T
    multiplicator = sizeof(T)
    bind(b)
    glBufferSubData(b.buffertype, multiplicator*(first(offset)-1), sizeof(value), value)
    bind(b, 0)
    return nothing
end

# copy between two buffers
# could be a setindex! operation, with subarrays for buffers
function unsafe_copy!(a::GLBuffer{T}, readoffset::Int, b::GLBuffer{T}, writeoffset::Int, len::Int) where T
    multiplicator = sizeof(T)
    glBindBuffer(GL_COPY_READ_BUFFER, a.id)
    glBindBuffer(GL_COPY_WRITE_BUFFER, b.id)
    glCopyBufferSubData(
        GL_COPY_READ_BUFFER, GL_COPY_WRITE_BUFFER,
        multiplicator*(readoffset-1),
        multiplicator*(writeoffset-1),
        multiplicator*len
    )
    glBindBuffer(GL_COPY_READ_BUFFER, 0)
    glBindBuffer(GL_COPY_WRITE_BUFFER, 0)
    return nothing
end

function Base.iterate(buffer::GLBuffer{T}) where T
    length(buffer) < 1 && return nothing
    glBindBuffer(buffer.buffertype, buffer.id)
    ptr = Ptr{T}(glMapBuffer(buffer.buffertype, GL_READ_WRITE))
    return (unsafe_load(ptr, i), (ptr, 1))
end

function Base.iterate(buffer::GLBuffer{T}, state::Tuple{Ptr{T}, Int}) where T
    ptr, i = state
    if i > length(buffer)
        glUnmapBuffer(buffer.buffertype)
        return nothing
    end
    val = unsafe_load(ptr, i)
    return (val, (ptr, i + 1))
end

#copy inside one buffer
function unsafe_copy!(buffer::GLBuffer{T}, readoffset::Int, writeoffset::Int, len::Int) where T
    len <= 0 && return nothing
    bind(buffer)
    ptr = Ptr{T}(glMapBuffer(buffer.buffertype, GL_READ_WRITE))
    for i=1:len+1
        unsafe_store!(ptr, unsafe_load(ptr, i+readoffset-1), i+writeoffset-1)
    end
    glUnmapBuffer(buffer.buffertype)
    bind(buffer,0)
    return nothing
end

function unsafe_copy!(a::Vector{T}, readoffset::Int, b::GLBuffer{T}, writeoffset::Int, len::Int) where T
    bind(b)
    ptr = Ptr{T}(glMapBuffer(b.buffertype, GL_WRITE_ONLY))
    for i=1:len
        unsafe_store!(ptr, a[i+readoffset-1], i+writeoffset-1)
    end
    glUnmapBuffer(b.buffertype)
    bind(b,0)
end

function unsafe_copy!(a::GLBuffer{T}, readoffset::Int, b::Vector{T}, writeoffset::Int, len::Int) where T
    bind(a)
    ptr = Ptr{T}(glMapBuffer(a.buffertype, GL_READ_ONLY))
    for i in 1:len
        b[i+writeoffset-1] = unsafe_load(ptr, i+readoffset-2) #-2 => -1 to zero offset, -1 gl indexing starts at 0
    end
    glUnmapBuffer(a.buffertype)
    bind(a,0)
end

function gpu_getindex(b::GLBuffer{T}, range::UnitRange) where T
    multiplicator = sizeof(T)
    offset = first(range)-1
    value = Vector{T}(undef, length(range))
    bind(b)
    glGetBufferSubData(b.buffertype, multiplicator*offset, sizeof(value), value)
    bind(b, 0)
    return value
end


text:
#=
all_tests.jl:
- Julia version:
- Author: bramb
- Date: 2019-01-24
=#
using Test

@testset "all tests" begin
    include("test_hypergraph_visualization.jl")
    include("test_rules_visualization.jl")
    include("test_hypergraph_external_nodes_count.jl")
end


text:
export QBasedPolicy

include("learners.jl")
include("explorers/explorers.jl")

import Functors

Base.@kwdef mutable struct QBasedPolicy{L,E} <: AbstractPolicy
    learner::L
    explorer::E
end

Functors.functor(x::QBasedPolicy) =
    (learner = x.learner,), y -> QBasedPolicy(y.learner, x.explorer)

(p::QBasedPolicy)(env) = p.explorer(p.learner(env), legal_action_space_mask(env))

RLBase.prob(p::QBasedPolicy, env::AbstractEnv) =
    prob(p.explorer, p.learner(env), legal_action_space_mask(env))

RLBase.optimise!(p::QBasedPolicy, x::NamedTuple) = optimise!(p.learner, x)


text:
# using Pkg
# pkg"add Oceananigans GLMakie"

ENV["GKSwstype"] = "100"

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using Printf
using Statistics
using Plots

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.Grids: xnode, ynode, znode

const Lx = 1000kilometers # zonal domain length [m]
const Ly = 2000kilometers # meridional domain length [m]

# number of grid points
Nx = 200
Ny = 400
Nz = 35

# stretched grid
k_center = collect(1:Nz)
Δz_center = @. 10 * 1.104^(Nz - k_center)
const Lz = sum(Δz_center)
z_faces = vcat([-Lz], -Lz .+ cumsum(Δz_center))
z_faces[Nz+1] = 0

arch = CPU()
FT = Float64

grid = VerticallyStretchedRectilinearGrid(architecture = arch,
                                          topology = (Periodic, Bounded, Bounded),
                                          size = (1, Ny, Nz),
                                          halo = (3, 3, 3),
                                          x = (0, Lx),
                                          y = (0, Ly),
                                          z_faces = z_faces)

grid = VerticallyStretchedRectilinearGrid(architecture = CPU(),
                                          topology = (Periodic, Bounded, Bounded),
                                          size = (grid.Nx, grid.Ny, grid.Nz),
                                          halo = (3, 3, 3),
                                          x = (0, grid.Lx),
                                          y = (0, grid.Ly),
                                          z_faces = z_faces)

xu, yu, zu = nodes((Face, Center, Center), grid)
xc, yc, zc = nodes((Center, Center, Center), grid)

u_timeseries = FieldTimeSeries("zonally_averaged_channel.jld2", "u", grid=grid)
b_timeseries = FieldTimeSeries("zonally_averaged_channel.jld2", "b", grid=grid)

@show b_timeseries

anim = @animate for i in length(b_timeseries.times)-15:length(b_timeseries.times)
    b = b_timeseries[i]
    u = u_timeseries[i]

    b_yz = interior(b)[1, :, :]
    u_yz = interior(u)[1, :, :]

    @show umax = max(1e-9, maximum(abs, u_yz))
    @show umax = maximum(abs, u_timeseries[:, :, :, :])
    @show bmax = max(1e-9, maximum(abs, b_yz))

    ulims = (-umax, umax) .* 0.8
    blims = (-bmax, bmax) .* 0.8

    ulevels = vcat([-umax], range(ulims[1], ulims[2], length=31), [umax])
    blevels = vcat([-bmax], range(blims[1], blims[2], length=31), [bmax])

    ylims = (0, grid.Ly) .* 1e-3
    zlims = (-grid.Lz, 0)

    u_yz_plot = contourf(yu * 1e-3, zu, u_yz',
                         xlabel = "y (km)",
                         ylabel = "z (m)",
                         aspectratio = :equal,
                         linewidth = 0,
                         levels = ulevels,
                         clims = ulims,
                         xlims = ylims,
                         ylims = zlims,
                         color = :balance)

    contour!(u_yz_plot,
             yc * 1e-3, zc, b_yz',
             linewidth = 1,
             color = :black,
             levels = blevels,
             legend = :none)
end

mp4(anim, "zonally_averaged_channel_small.mp4", fps = 8) # hide


text:
import Base.valtype
import Distributions.logpdf

"""
    ODESolverType

Types inheriting from abstract type `ODESolverType` declare the numerical
schemes used for finding solutions to backward ODEs
"""
abstract type ODESolverType end

"""
    ODEChangePt

Types inheriting from abstract type `ODEChangePt` decide upon which ODE solvers
are to be used and when to ultimately compute the triplet H,Hν,c.
"""
abstract type ODEChangePt end

"""
    NoChangePt <: ODEChangePt

Struct indicating that only solvers for H,Hν,c are to be employed. The field `λ`
indicates the amount of space (in the units of the number of elements of vector)
that needs to be nevertheless reserved for L,M⁺,μ (even though the latter are
not used).
"""
struct NoChangePt <: ODEChangePt
    λ::Int64
    NoChangePt(λ=0) = new(λ)
end

"""
    SimpleChangePt <: ODEChangePt

Struct indicating that both types of solvers are to be used; for H,Hν,c as well
as for L,M⁺,μ. The ODE solvers for L,M⁺,μ are used first on the terminal part
of the interval, where `λ` gives the lenght of this terminal interval. The ODE
solvers for H,Hν,c are used from then on.
"""
struct SimpleChangePt <: ODEChangePt
    λ::Int64
    SimpleChangePt(λ=0) = new(λ)
end

"""
    getChangePt(changePt::ODEChangePt)

Return the length of terminal interval over which ODE solvers for L,M⁺,μ are
used
"""
getChangePt(changePt::ODEChangePt) = changePt.λ

#TODO implement Jeffrey's priors
"""
    ImproperPrior

Flat prior
"""
struct ImproperPrior end
logpdf(::ImproperPrior, θ) = 0.0

valtype(::Val{T}) where T = T

"""
    idx(::Val{T}) where {T}

Return a tuple containing indices of parameters selected by Val{T}

# Examples
```julia-repl
julia> idx(Val((true, false, false, true, true)))
(1, 4, 5)
```
"""
function idx(::Val{T}) where T
    tuple((i for i in 1:length(T) if T[i])...)
end

"""
    moveToProperPlace(ϑ, θ, ::Val{T}) where {T}

Update parameter vector `θ` at indices specified by `Val{T}` with the
values collected in `ϑ`.

# Examples
```julia-repl
julia> moveToProperPlace([10,20,30], [1,2,3,4,5],
                         Val((true, false, true, false, true)))
5-element Array{Float64,1}:
 10.0
  2.0
 20.0
  4.0
 30.0
```
"""
function moveToProperPlace(ϑ, θ, ::Val{T}) where {T}
    v = valtype(Val{T}())
    m = length(v)
    θᵒ = zero(θ)
    idxNew = [i for i in 1:m if v[i]]
    idxOld = [i for i in 1:m if !v[i]]
    θᵒ[idxNew] .= ϑ
    θᵒ[idxOld] .= θ[idxOld]
    θᵒ
end


text:
function get_panda(;
        timestep=0.01,
        gravity=-9.81,
        friction_coefficient=0.8,
        spring=0.0,
        damper=0.01, # this value comes from the official URDF https://github.com/frankaemika/franka_ros/blob/develop/franka_gazebo/test/launch/panda-gazebo.urdf
        contact=false,
        limits=true,
        model_type=:end_effector,
        # joint_limits=[[-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25],
        #               [ 0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25]],
        # joint_limits=[[-2.9671, -1.8326, -2.9671, -0.4000, -2.9671, -3.8223, -2.9671],
        #               [ 2.9671,  1.8326,  2.9671,  3.1416,  2.9671,  0.0873,  2.9671]],
        joint_limits=[[-2.8973, -1.7628, -2.8973, -0.0698, -2.8973, -3.7525, -2.8973, -0.00],
                      [ 2.8973,  1.7628,  2.8973,  3.0718,  2.8973,  0.0175,  2.8973,  0.04]],
        T=Float64)

    path = joinpath(@__DIR__, "../deps/panda_$(String(model_type)).urdf")
    mech = Mechanism(path, false, T,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    # Adding springs and dampers
    for joint in mech.joints
        joint.damper = true
        joint.spring = true
        joint.translational.spring=spring
        joint.translational.damper=damper
        joint.rotational.spring=spring
        joint.rotational.damper=damper
    end

    # joint limits
    joints = deepcopy(mech.joints)
    if limits
        for (i,joint) in enumerate(joints)
            id = joint.id
            if input_dimension(joint.translational) == 0 && input_dimension(joint.rotational) == 1
                joints[i] = add_limits(mech, joint,
                    rot_limits=[SVector{1}(joint_limits[1][id-1]), SVector{1}(joint_limits[2][id-1])])
            end
            if input_dimension(joint.translational) == 1 && input_dimension(joint.rotational) == 0
                joints[i] = add_limits(mech, joint,
                    tra_limits=[SVector{1}(0.00), SVector{1}(0.04)])
            end
        end
        mech = Mechanism(Origin{T}(), [mech.bodies...], [joints...],
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end

    origin = Origin{T}()
    bodies = mech.bodies
    joints = mech.joints
    contacts = ContactConstraint{T}[]

    if contact
        if model_type == :end_effector
            # TODO place the contact points for each finger of the end-effector
        elseif model_type == :no_end_effector
            # spherical end-effector contact
            location = [-0.01; 0.004; 0.01]
            normal = [0.0; 0.0; 1.0]
            offset = [0.0; 0.0; 0.05]
            contact = contact_constraint(
                get_body(mech, :link7),
                normal,
                friction_coefficient=friction_coefficient,
                contact_point=location,
                offset=offset,
                name=:end_effector)
            push!(contacts, contact)
        end
    end

    mech = Mechanism(origin, bodies, joints, contacts,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    set_minimal_state!(mech, szeros(minimal_dimension(mech)))
    return mech
end

function initialize_panda!(mechanism::Mechanism{T};
    joint_angles=[[0,-0.8,0,1.6,0,-2.4,0]; zeros(input_dimension(mechanism)-7)],
    joint_velocities=zeros(input_dimension(mechanism))) where T

    nu = input_dimension(mechanism)
    zero_velocity!(mechanism)
    y = vcat([[joint_angles[i], joint_velocities[i]] for i=1:nu]...)
    set_minimal_state!(mechanism, y)
    return nothing
end


text:
using Documenter, DiffEqBayesStan

makedocs(
    modules = [DiffEqBayesStan],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "DiffEqBayesStan.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/StanJulia/DiffEqBayesStan.jl.git",
)


text:
#=
The Agfa ColorReference is designed according to the IT8.7 standard.

The spectral reflectance of samples were measured with Minolta spectro-
photometer with option SCE (Specular component excluded). The CM-2002
utilizes d/8 (diffuse illumination/8o viewing) geometry which conforms to
ISO and DIN standards (and also to CIE d/0 recommendations).

Measuring equipment: Minolta spectrofotometer CM-2002.
Wavelength interval: 400 nm - 700 nm.
Wavelength resolution: 10 nm.

Measurer: ela@ee.oulu.fi (University of Oulu, Electrical Eng.Dept)

Further information: ela@ee.oulu.fi (University of Oulu, Electrical Eng.Dept)

References:

   AgfaFotoreference, Trademark of Agfa-Gevaert A.G.,
   Septestraat 27 2640 Mortsel - Belgium,1992.
=#
"""
`agfait_specs()`
"""
function agfait_specs()
    filename = "agfait872.dat"
    fp=@__DIR__
    cd(fp)
    cd("..\\spectra")
    readlines(filename)
end

"""
`agfait(idx::Int)`

Returns the reflectance spectrum of the selected Agfait872 color patch.
Spectral data are from 400 to 700 nm in 10 nm steps.

1st group: 1-264

2nd group: Neutral scale samples (265-286)

3rd group:

* black sample (287)
* white sample (288)
* Minolta white calibration sample (289)

# Examples
```jldoctest
julia> agfait(123)
RSpec(Real[400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, 480.0, 490.0  …  610.0, 620.0, 630.0, 640.0, 650.0, 660.0, 670.0, 680.0, 690.0, 700.0], Real[0.1282, 0.3004, 0.4806, 0.5796, 0.6429, 0.6818000000000001, 0.7078, 0.7159, 0.7071999999999999, 0.6896  …  0.259, 0.2318, 0.21309999999999998, 0.2, 0.1932, 0.19210000000000002, 0.19699999999999998, 0.2095, 0.2311, 0.2594])

julia> agfait(0)
ERROR: DomainError with 0:
Color index out of range. Expecting 0 ≤ index ≤ 289.
```

Measurer: ela@ee.oulu.fi (University of Oulu, Electrical Eng.Dept)
"""
function agfait(idx::Int)
    if 1 ≤ idx ≤ 289
        file = agfait_specs()
        start = 119 * (idx - 1) + 2
        lam = collect(400.0:10.0:700.0)
        spd = map(x -> (parse(Float64, x) / 100.0), file[start:start+30])
        reflectance_spec(lam, spd)
    else
        throw(DomainError(idx, "Color index out of range. Expecting 0 ≤ index ≤ 289."))
    end
end


text:
parseInt(x)=parse(Int, x)
parseMap(x::Array{SubString{String},1})=map(parseInt, x)
pprint(it...;sep=" ",ed="\n")=length(it)>0 ? (@inbounds for i in 1:length(it)-1;print(it[i],sep);end;print(last(it),ed)) : print(ed)

n = readline() |> parseInt
a = readline() |> split |> parseMap
b = zeros(Int,n)
for i in 1:n-1
    b[a[i]]+=1
end
pprint(b...,sep='\n')

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
# Data preprocessing, normalization and all such stuff
# Use `prepare()` to create normalized top N dataset and mappings

using DataFrames

require("data.jl")

function count_by(dataset, col)
    counts = Dict{Any, Int64}()
    for i=1:nrow(dataset)
        val = dataset[i, col]
        if haskey(counts, val)
            counts[val] += 1
        else
            counts[val] = 1
        end
    end
    return counts
end

function most_popular(dataset, col, n=1000)
    # take rows of n most popular (most frequent) col values
    val2count = count_by(dataset, col)
    count2val = [(v, k) for (k, v) in val2count]
    top_n_vals = Set([v for (k, v) in sort(count2val, rev=true)[1:n]])
    rows = Int64[]
    for i=1:nrow(dataset)
        if in(dataset[i, col], top_n_vals)
            push!(rows, i)
        end
    end
    return dataset[rows, :]
end

function filter_low(dataset, threshold)
    return dataset[dataset[:plays] .>= threshold, :]
end

function normalize_dataset(dataset)
    uid2userid = sort(unique(dataset[:userid]))
    userid2uid = Dict([(uid2userid[i], i) for i=1:length(uid2userid)])
    dataset[:uid] = zeros(Int32, nrow(dataset))
    for i=1:nrow(dataset)
        dataset[i, :uid] = userid2uid[dataset[i, :userid]]
    end
    delete!(dataset, :userid)

    aid2artname = sort(unique(dataset[:artname]))
    artname2aid = Dict([(aid2artname[i], i) for i=1:length(aid2artname)])
    dataset[:aid] = zeros(Int32, nrow(dataset))
    for i=1:nrow(dataset)
        dataset[i, :aid] = artname2aid[dataset[i, :artname]]
    end
    delete!(dataset, :artname)
    delete!(dataset, :artmbid)
    return dataset, uid2userid, userid2uid, aid2artname, artname2aid
end

function prepare()
    println("Reading full dataset")
    @time df = read_full_dataset()
    println("Collecting most popular artists")
    @time df = most_popular(df, :artname, TOPN)
    println("Filtering out records with too little plays")
    df = filter_low(df, LOW_THRESHOLD)
    println("Normalizing dataset")
    @time norm_dat = normalize_dataset(df)
    println("Saving")
    @time begin
        writetable(SHORT_DATASET_FILE, norm_dat[1], separator='\t')
        save(ARTIST_FILE, "artists", norm_dat[4])
        save(USER_FILE, "users", norm_dat[2])
    end
end

max_stars_repo_name: dfdx/lastfm-rbm

max_issues_repo_name: dfdx/lastfm-rbm

max_forks_repo_name: dfdx/lastfm-rbm


text:
using DelimitedFiles

function read_vector(doc, n)
    vector = Array{Int, 1}()
    data = readdlm(doc, skipstart=n)
    for i in data[1,:]
        append!(vector, i)
    end
    return(n+1, vector)
end

function read_matrix(P, Q, doc, n)
    matrix = zeros(Int, P, Q)
    data = readdlm(doc, skipstart=n, use_mmap=true)
    for i = 1:P
        for j = 1:Q
            matrix[i,j] = data[i,j]
        end
    end
    return(n+P, matrix)
end

function read_data(doc)
    open(doc, "r") do f
        n = 7
        C = parse(Int, readline(f)) #number of containers
        P = parse(Int, readline(f)) #number of positions
        CP = parse(Int, readline(f)) #number of compatible combinations
        Q = parse(Int, readline(f)) #number of quay cranes
        J = parse(Int, readline(f)) #number of bays
        tt = parse(Int, readline(f)) #quay crane travel time per bay
        d = parse(Int, readline(f)) #clearance, in number of bays

        n, ci = read_vector(doc, n) #container class vector
        n, pj = read_vector(doc, n) #position class vector
        n, bj = read_vector(doc, n) #position class vector

        #n, qj = read_matrix(P,Q,doc,n) #quay cranes position matrix
        #n, sj = read_matrix(P,S,doc,n) #transfer vehicle position matrix
        n, prejj = read_matrix(P, P, doc, n) #position precedence matrix
        n, cpij = read_matrix(CP, 3, doc ,n) #combination c, p, time(cp)
        #return(C,P,CP,S,Q,J,tt,d,mt,beta,EFT,ci,pj,qj,sj,prejj,cpij)
        return(C,P,CP,Q,J,tt,d,ci,pj,bj,prejj,cpij)
    end
end
#doc = "../data/5C_3Type_2QC.txt"
#doc = "../data/QCSP_example.txt"
#doc = "../../data/10C_4Type_2QC.txt"
#doc = "./data/Benchmark/60C_10Type_LessDense_2QC.txt"
#doc = "./data/Benchmark/240C_20Type_LessDense_2QC.txt"
doc = "../data/Benchmark/240C_20Type_LessDense_2QC.txt"
#C,P,CP,S,Q,J,tt,d,mt,beta,EFT,ci,pj,qj,sj,prejj,cpij = read_data(doc)
C,P,CP,Q,J,tt,d,ci,pj,bj,prejj,cpij = read_data(doc)

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
module Actions
    const shift = 1
    const left = 2
    const right = 3
end

type State
    tokens::Vector{Token}
    step::Int
    top::Int
    left::State
    right::Int
    lch::State # left-most child
    rch::State
    score::Float64
    prev::State

    State() = new(Token[], 0)
    State(tokens) = new(tokens, 1, 1, State(), 2, State(), State(), 0.0)
    State(tokens, step, top, left, right, lch, rch) = new(tokens, step, top, left, right, lch, rch, 0.0, left)
end

#State(tokens::Vector{Token}) = State(tokens, 1, 1, State(), 2, State(), State())

function State(s::State, act::Int)
    top, left, right, lch, rch = begin
        if act == Actions.shift
            s.right, s, s.right+1, State(), State()
        elseif act == Actions.left
            s.top, s.left.left, s.right, s.left, s.rch
        elseif act == Actions.right
            s.left.top, s.left.left, s.right, s.left.lch, s
        else
            throw("Invalid action: $(act)")
        end
    end
    State(s.tokens, s.step+1, top, left, right, lch, rch)
end

Base.done(s::State) = s.step == 2 * (length(s.tokens)-1) + 1

Base.isnull(s::State) = s.step == 0

function nextpred(s::State)::Vector{State}
    done(s) && return State[]
    states = State[]
    isnull(s.left) || push!(states, State(s,Actions.left), State(s,Actions.right))
    s.right <= length(s.tokens) && push!(states, State(s,Actions.shift))
    states
end

function nextgold(s::State)::Vector{State}
    done(s) && return State[]
    if isnull(s.left)
        act = Actions.shift
    else
        tokens = s.tokens
        s0, s1 = tokens[s.top], tokens[s.left.top]
        if s1.head == s.top
            act = Actions.left
        elseif s0.head == s.left.top
            reducable = all(i -> tokens[i].head != s.top, s.right:length(tokens))
            act = reducable ? Actions.right : Actions.shift
        elseif s.right <= length(tokens)
            act = Actions.shift
        else
            throw("Invalid data.")
        end
    end
    [State(s,act)]
end

function toheads(s::State)
    heads = fill(0, length(s.tokens))
    while s != nothing
        length(s.lcs) > 0 && (heads[s.lcs[1]] = s.top)
        length(s.rcs) > 0 && (heads[s.rcs[1]] = s.top)
        s = s.prev
    end
    heads
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
abstract type EGOAquisition end

update_parameters!(::EGOAquisition, args...) = error("NOT IMPLEMENTED")
objective(::EGOAquisition, m::Kriging, x) = error("NOT IMPLEMENTED")

import StatsFuns

const φ = StatsFuns.normpdf
const Φ = StatsFuns.normcdf

include("single_obj_aquisitions.jl")
include("constrained_aquisition.jl")


text:
@testset "SAFT methods single phase" begin
    system = PCSAFT(["ethanol"])
    @test Clapeyron.volume(system, 1E5, 298) ≈ 5.9069454773905654e-5 rtol = 0.01 #returns incorrect value
    @test Clapeyron.sat_pure(system, 298)[1] ≈ 7904.2246894075415 rtol = 0.01
    @test Clapeyron.enthalpy_vap(system, 298) ≈ 41720.97771100548 rtol = 0.01
    @test Clapeyron.crit_pure(system)[1] ≈ 1 rtol = 1E10 #T_scale not defined
    @test Clapeyron.pressure(system, 1, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.entropy(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.chemical_potential(system, 1E5, 298)[1] ≈ 1 rtol = 1E10
    @test Clapeyron.internal_energy(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.enthalpy(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.gibbs_free_energy(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.helmholtz_free_energy(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.isochoric_heat_capacity(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.isobaric_heat_capacity(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.isothermal_compressibility(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.isentropic_compressibility(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.speed_of_sound(system, 1E5, 298) ≈ 1 rtol = 1E10 #requires that the model has Mr
    @test Clapeyron.isobaric_expansivity(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.joule_thomson_coefficient(system, 1E5, 298) ≈ 1 rtol = 1E10
    @test Clapeyron.second_virial_coefficient(system, 298) ≈ 1 rtol = 1E10
end


text:
module ComoniconCast

using Pkg
using Markdown
using ExproniconLite
using ComoniconOptions
using ComoniconTypes

export JLArgument, JLOption, JLFlag, JLMD, JLMDFlag, JLMDOption,
    @cast, @main, cast, cast_args, cast_flags, cast_options, default_name,
    get_version, split_leaf_command, split_docstring, read_arguments, read_description,
    read_options, read_flags, split_hint, split_option

include("types.jl")
include("cast.jl")
include("markdown.jl")
include("utils.jl")

end