text:
abstract type AbstractPathSampler end

function _samplepaths!(_npaths, npaths, pmf, semirandom::Bool, canmodifypmf::Bool)
    if semirandom
        sampled = 0
        for i in 1:length(pmf)
            p = pmf[i]
            n = floor(Int, p * npaths + 1e-6)
            if n > 0
                if sampled == 0 && !canmodifypmf
                    pmf = copy(pmf)
                end
                pmf[i] -= n / npaths
                _npaths[i] += n
                sampled += n
            end
        end
        if sampled < npaths
            if sampled > 0
                # sum(pmf) should be equal to (npaths - sampled) / npaths
                # since we have removed sampled / npaths
                # pmf is not modified as it could be the vector of probas
                pmf ./= (npaths - sampled) / npaths
                npaths -= sampled
            end
            _samplepaths!(_npaths, npaths, pmf, false, true)
        else
            _npaths
        end
    else
        cmf = cumsum(pmf)
        @assert isapprox(cmf[end], 1)
        cmf[end] = 1
        samples = rand(Float64, npaths)
        sort!(samples)
        i = 1
        for j in samples
            while j >= cmf[i]
                i += 1
            end
            _npaths[i] += 1
        end
        _npaths
    end
end
function _samplepaths(npaths, pmf, semirandom::Bool, canmodifypmf::Bool)
    _samplepaths!(zeros(Int, length(pmf)), npaths, pmf, semirandom, canmodifypmf)
end
infpaths(g, node) = fill(-1, length(SOI.get(g, SOI.OutTransitions(), node)))

function samplepaths(pathsampler::AbstractPathSampler, sp::SOI.AbstractStochasticProgram, node, npaths::Vector{Int}, t, num_stages)
    npathss = Vector{Int}[similar(npaths) for i in 1:length(SOI.get(sp, SOI.OutTransitions(), node))]
    for i in eachindex(npaths)
        _npaths = samplepaths(pathsampler, sp, node, npaths[i], t, num_stages)
        for c in eachindex(npathss)
            npathss[c][i] = _npaths[c]
        end
    end
    npathss
end

struct ProbaPathSampler <: AbstractPathSampler
    semirandom::Bool
end
function samplepaths(pathsampler::ProbaPathSampler, g::SOI.AbstractStochasticProgram, node, npaths::Int, t, num_stages)
    if npaths == -1
        infpaths(g, node)
    else
        pmf = SOI.get.(g, SOI.Probability(), SOI.get(g, SOI.OutTransitions(), node))
        _samplepaths(npaths, pmf, pathsampler.semirandom, false)
    end
end

struct NumPathsPathSampler <: AbstractPathSampler
    semirandom::Bool
end
function samplepaths(pathsampler::NumPathsPathSampler, g::SOI.AbstractStochasticProgram, node, npaths::Int, t, num_stages)
    if npaths == -1
        infpaths(g, node)
    else
        # TODO check if needs +-1 in argument of NumberOfPathsFrom
        den = SOI.get(g, SOI.NumberOfPathsFrom(num_stages - (t-1)), node)
        pmf = map(tr->SOI.get(g, SOI.NumberOfPathsFrom(num_stages - t), SOI.get(g, SOI.Target(), tr)) / den, SOI.get(g, SOI.OutTransitions(), node))
        _samplepaths(npaths, pmf, pathsampler.semirandom, true)
    end
end


text:
include("../examples/baker.jl")
using LinearAlgebra
using PyPlot
using PyCall
@pyimport matplotlib.collections as coll
@pyimport matplotlib.colors as cs
@pyimport matplotlib.cm as cm
function plot_w(m)
	s = [0.3, 0.1, 0.1, 1.0]
	n = 100
	x = LinRange(0.,2*pi,n)
	u = [[x[i] + 2*pi*rand(),x[j] + 2*pi*rand()]
		 for i=1:n,j = 1:n][:]
	q = [rand(2) for i=1:n, j=1:n][:]
	w = [zeros(2) for i=1:n, j=1:n][:]
	D2 = [zeros(2,2) for i=1:n, j=1:n][:]
	n = length(u)
	z = zeros(n)
	nw = zeros(n)
	segments = zeros(n, 2, 2)
	eps = 1.e-1
	filepath = string("/home/nishac/Research/PhDThesis/papers/",
					  "Decomposition-of-Linear-Response/figs/")
	for i = 1:20
		@show i
		u .= next.(u, Ref(s))
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
	end

	for i = 1:m
		@show i
		D2 .= pushforward_second_order.(u, q, Ref(s))
		w .= pushforward.(u, w, Ref(s))
		w .+= tensordot(D2, q)
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
		w ./= (z.*z)
		w .= w .- dot.(w, q).*q
		u .= next.(u, Ref(s))
	end
	nw = norm.(w)
	segments .= create_line_colls(u, q, eps)
	lc = coll.LineCollection(segments,
						cmap=plt.get_cmap("cool"),
						norm=cs.Normalize(minimum(nw),
						maximum(nw)))
	lc.set_array(nw)
	lc.set_linewidth(2)

	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_ylim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.set_ylabel(L"$x_2$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	ax.add_collection(lc)
	ax.axis("scaled")

	cbar = fig.colorbar(cm.ScalarMappable(
						norm=cs.Normalize(minimum(nw),
						maximum(nw)),
					   cmap=plt.get_cmap("cool")), ax=ax,
						orientation="horizontal",shrink=0.4,
						pad=0.1)

	cbar.ax.tick_params(labelsize=30)
	cbar.ax.xaxis.get_offset_text().set_fontsize(30)
	plt.tight_layout()
	#savefig(string(filepath,"w_n_",m,".png"))
	return segments, nw
end
function plot_y(m)
	s = [0.3, 0., 0.3, 0.0]
	n = 180
	x = LinRange(0.,2*pi,n)
	u = [[x[i] + 2*pi*rand(),x[j] + 2*pi*rand()]
		 for i=1:n,j = 1:n][:]
	q = [rand(2) for i=1:n, j=1:n][:]
	w = [zeros(2) for i=1:n, j=1:n][:]
	y = [zeros(2) for i=1:n, j=1:n][:]
	v = [zeros(2) for i=1:n, j=1:n][:]
	X = [zeros(2) for i=1:n, j=1:n][:]

	D2 = [zeros(2,2) for i=1:n, j=1:n][:]
	DX = [zeros(2,2) for i=1:n, j=1:n][:]
	n = length(u)
	z = zeros(n)
	a = zeros(n)
	delta = zeros(n)
	ny = zeros(n)
	segments = zeros(n, 2, 2)
	eps = 1.e-1
	filepath = string("/home/nishac/Research/PhDThesis/papers/",
					  "Decomposition-of-Linear-Response/figs/")
	for i = 1:20
		@show i
		u .= next.(u, Ref(s))
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
	end

	for i = 1:m
		@show i
		D2 .= pushforward_second_order.(u, q, Ref(s))
		y .= pushforward.(u, y, Ref(s))
		y .+= tensordot(D2, v)
		X .= pert.(u, 1) .+ pert.(u, 3)
		DX .= dpert.(u,1) .+ dpert.(u,3)
		v .= pushforward.(u, v, Ref(s))
		v .+= X
		w .= pushforward.(u, w, Ref(s))
		w .+= tensordot(D2, q)

		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z

		a .= dot.(v, q)
		v .-= a.*q

		w ./= (z.*z)
		w .= w .- dot.(w, q).*q
		y ./= z
		y .+= -a.*w .+ tensordot(DX, q)
		delta .= dot.(y, q) .+ dot.(v, w)
		y .-= delta.*q
		u .= next.(u, Ref(s))
	end
	qperp = [[-qi[2], qi[1]] for qi in q]
	ny = dot.(y, qperp)
	segments .= create_line_colls(u, q, eps)
	lc = coll.LineCollection(segments,
						cmap=plt.get_cmap("coolwarm"),
						norm=cs.Normalize(minimum(ny),
						maximum(ny)))
	lc.set_array(ny)
	lc.set_linewidth(2)

	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_ylim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.set_ylabel(L"$x_2$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	ax.add_collection(lc)
	ax.axis("scaled")

	cbar = fig.colorbar(cm.ScalarMappable(
						norm=cs.Normalize(minimum(ny),
						maximum(ny)),
					   cmap=plt.get_cmap("coolwarm")), ax=ax,
						orientation="horizontal",shrink=0.4,
						pad=0.1)

	cbar.ax.tick_params(labelsize=30)
	cbar.ax.xaxis.get_offset_text().set_fontsize(30)
	plt.tight_layout()
	#savefig(string(filepath,"w_n_",m,".png"))
	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_ylim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.set_ylabel(L"$x_2$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)
	pts = hcat(u...)
	vecs = hcat(y...)
	x_pts = [pts[1,:] - eps*vecs[1,:] pts[1,:] + eps*vecs[1,:]]'
	y_pts = [pts[2,:] - eps*vecs[2,:] pts[2,:] + eps*vecs[2,:]]'

	ax.plot(x_pts, y_pts, "darkcyan")
	ax.axis("scaled")

	return segments, ny
end

function plot_v(m)
	s = [0.1, 0.1, 0., 0.1]
	n = 180
	x = LinRange(0.,2*pi,n)
	u = [[x[i] + 2*pi*rand(),x[j] + 2*pi*rand()]
		 for i=1:n,j = 1:n][:]
	q = [rand(2) for i=1:n, j=1:n][:]
	v = [zeros(2) for i=1:n, j=1:n][:]
	X = [zeros(2) for i=1:n, j=1:n][:]
	n = length(u)
	z = zeros(n)
	a = zeros(n)
	segments = zeros(n, 2, 2)
	eps = 1.e-1
	filepath = string("/home/nishac/Research/PhDThesis/papers/",
					  "Decomposition-of-Linear-Response/figs/")
	for i = 1:20
		@show i
		u .= next.(u, Ref(s))
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
	end
	for i = 1:m
		@show i
		X .= pert.(u, 1) .+ pert.(u, 2) .+ pert.(u, 4)
		v .= pushforward.(u, v, Ref(s))
		v .+= X
		q .= pushforward.(u, q, Ref(s))
		z .= norm.(q)
		q .= q./z
		a .= dot.(v, q)
		v .-= a.*q
		u .= next.(u, Ref(s))
	end
	a .= norm.(v)
	segments .= create_line_colls(u, q, eps)
	lc = coll.LineCollection(segments,
						cmap=plt.get_cmap("coolwarm"),
						norm=cs.Normalize(minimum(a),
						maximum(a)))
	lc.set_array(a)
	lc.set_linewidth(2)

	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_ylim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.set_ylabel(L"$x_2$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)

	ax.add_collection(lc)
	ax.axis("scaled")

	cbar = fig.colorbar(cm.ScalarMappable(
						norm=cs.Normalize(minimum(a),
						maximum(a)),
					   cmap=plt.get_cmap("coolwarm")), ax=ax,
						orientation="horizontal",shrink=0.4,
						pad=0.1)

	cbar.ax.tick_params(labelsize=30)
	cbar.ax.xaxis.get_offset_text().set_fontsize(30)
	plt.tight_layout()
	#savefig(string(filepath,"w_n_",m,".png"))
	fig, ax = subplots(1,1)
	ax.set_xlim([0,2*pi])
	ax.set_ylim([0,2*pi])
	ax.set_xlabel(L"$x_1$", fontsize=30)
	ax.set_ylabel(L"$x_2$", fontsize=30)
	ax.xaxis.set_tick_params(labelsize=30)
	ax.yaxis.set_tick_params(labelsize=30)
	pts = hcat(u...)
	vecs = hcat(v...)
	x_pts = [pts[1,:] - eps*vecs[1,:] pts[1,:] + eps*vecs[1,:]]'
	y_pts = [pts[2,:] - eps*vecs[2,:] pts[2,:] + eps*vecs[2,:]]'

	ax.plot(x_pts, y_pts, "royalblue")
	ax.axis("scaled")

	return segments, a
end
function tensordot(A, b)
	return [Ai*b[i] for (i, Ai) in enumerate(A)]
end
function pushforward_second_order(u, v1, s)
	d2u = d2step(u, s)
	d2u_v1 = reshape([dot(d2u[:,i], v1) for i=1:4], 2, 2)
	return d2u_v1
end
function create_line_colls(u, q, eps)
	n = length(u)
	lc = zeros(n, 2, 2)
	for i = 1:n
		lc[i,1,:] = u[i] .- eps*q[i]
		lc[i,2,:] = u[i] .+ eps*q[i]
	end
	return lc
end
function pushforward(u, q, s)
	du = dstep(u, s)
	q = du*q
	return q
end

ext: jl

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using GR

   w = 400
   h = 200

#setviewport(0, 1, 0, 1)
#setwindow(1, w, 1, h)
#=
   mwidth, mheight, width, height = GR.inqdspsize()
   if w > h
        ratio = float(h) / w
        msize = mwidth * w / width
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
    else
        ratio = float(w) / h
        msize = mheight * h / height
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
    end
=#
println(msize," ", msize * ratio, ", ratio: ", ratio)
setviewport(0, msize, 0, msize * ratio)
setwindow(1, w, 1, h)

setmarkertype(GR.MARKERTYPE_DOT)
setmarkersize(1)

setcolorrep(1, 0.0, 0.0, 1.0)
setlinecolorind(1)
setfillcolorind(1)

i = 0
while i < 300000

  #clearws()
  #fillrect(1, 1.1, 0, ratio)
  #polymarker([1, 400], [1, 86])
  polymarker([100], [86])
  updatews()
  i += 1
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using Test
using SMTPClient

@testset "Errors" begin
    @testset "Error message for Humans(TM)" begin
        let errmsg = "Couldn't resolve host name"
            server = "smtp://nonexists"
            body = IOBuffer("test")

            try
                send(server, ["nobody@earth"], "nobody@earth", body)
                @assert false, "send should fail"
            catch e
                @test occursin(string(errmsg), string(e))
            end
        end
    end

    @testset "Non-blocking send" begin
        let errmsg = "Couldn't resolve host name"
            opt = SendOptions(blocking = false)
            server = "smtp://nonexists"
            body = IOBuffer("test")

            future = send(server, ["nobody@earth"], "nobody@earth", body, opt)
            e = fetch(future)
            @test occursin(string(errmsg), string(e))
        end
    end
end


text:
function GenVertexArrays(count::Integer)
	a = Array(VertexArray, count)
	ccall( (:glGenVertexArrays, lib), Void, (GLsizei, Ptr{Buffer}), count, a)
	GetError()
	return a
end

function GenVertexArray()
	return GenVertexArrays(1)[1]
end

function BindVertexArray(va::VertexArray)
	ccall( (:glBindVertexArray, lib), Void, (VertexArray,), va)
	GetError() # TODO: Benchmark overhead
end
BindVertexArray(va::Integer) = BindVertexArray(VertexArray(va))


text:
@testset "shapes" begin
    inputsize = 20
    poollength = 10
    datalength = 100
    x = rand(Float32, inputsize, poollength, 1, datalength)
    m = DARNN(inputsize, 10, 10, poollength, 1)
    @test size(m(x)) == (1, datalength)
    if Flux.CUDA.functional()
        @test size(gpu(m)(gpu(x))) == (1, datalength)
    end
end

@testset "misc" begin
    @test repr(DARNN(10, 10, 10, 1, 1)) == "DARNN(10, 10, 10, 1, 1)"
end

# @testset "constructors" begin
#     @test size(Dense(10, 100).weight) == (100, 10)
#     @test size(Dense(10, 100).bias) == (100,)
#     @test Dense(rand(100, 10), rand(100)).σ == identity
#     @test Dense(rand(100, 10)).σ == identity
# end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using ArraysOfArrays, Distributions, StatsBase, PDMats, IntervalSets, ValueShapes

@testset "distribution_density" begin
    mvt = @inferred MvTDist(1.5, PDMat([2.0 1.0; 1.0 3.0]))
    mvdd = @inferred BAT.DistributionDensity(mvt)

    @testset "properties" begin
        @test typeof(mvdd) <: AbstractDensity
        @test parent(mvdd) == mvt
        @test totalndof(mvdd) == 2
    end

    @testset "BAT.eval_logval_unchecked" begin
        @test (@inferred BAT.eval_logval_unchecked(mvdd, [0.0, 0.0])) ≈ -2.64259602
    end

    @testset "BAT.var_bounds" begin
        let
            dist = @inferred NamedTupleDist(a = 5, b = Normal(), c = -4..5, d = MvNormal([1.2 0.5; 0.5 2.1]), e = [Normal(1.1, 0.2)] )
            density = @inferred BAT.DistributionDensity(dist)
            s = BAT.bat_sampler(unshaped(density))
            @test all([rand(s) in BAT.var_bounds(density) for i in 1:10^4])
        end
    end

    @testset "statistics" begin
        mvn = @inferred(product_distribution([Normal(-1.0), Normal(0.0), Normal(1.0)]))
        dist_density = @inferred(BAT.DistributionDensity(mvn))
        @test @inferred(bat_findmode(dist_density)).result == mode(mvn)
    end
end


text:
abstract type Movielens <: AbstractDataset end

function load_dataset(d::Movielens)
    # Movielens dataset has different kinds of meta data,
    # depending on the size of rating.
    # So we leave load function as abstract.
    error("load method is not implemented.")
end

function load_inter(d::Movielens)
    error("load_inter method is not implemented.")
end

function download(d::Movielens; kwargs...)
    download_zip(d.url, d.name, d.dataset_dir; kwargs...)
    return
end

"""
Movielens100k
"""
@with_kw_noshow struct Movielens100k <: Movielens
    name::AbstractString = "movielens100k"
    url::AbstractString = "http://files.grouplens.org/datasets/movielens/ml-100k.zip"
    dataset_dir::AbstractString = joinpath(@__DIR__, "..", "..", "dataset", "movielens100k")
end

function load_dataset(d::Movielens100k)
    return load_inter(d), load_user(d), load_item(d)
end

function load_inter(d::Movielens100k)
    csv = CSV.File(
        joinpath(d.dataset_dir, "u.data"),
        delim = "\t",
        header = [:userid, :movieid, :rating, :timestamp],
    )
    return csv
end

function load_user(d::Movielens100k)
    csv = CSV.File(
        joinpath(d.dataset_dir, "u.user"),
        delim = "|",
        header = [:userid, :age, :gender, :occupation, :zipcode],
    )
    return csv
end

function load_item(d::Movielens100k)
    csv = CSV.File(
        joinpath(d.dataset_dir, "u.item"),
        delim = "|",
        header = [
            :movieid,
            :movie_title,
            :release_date,
            :video_release_date,
            :IMDbURL,
            :unknown,
            :Action,
            :Adventure,
            :Animation,
            :Childrens,
            :Comedy,
            :Crime,
            :Documentary,
            :Drama,
            :Fantasy,
            :FilmNoir,
            :Horror,
            :Musical,
            :Mystery,
            :Romance,
            :SciFi,
            :Thriller,
            :War,
            :Western,
        ],
    )
    return csv
end

"""
Movielens1M
"""

@with_kw_noshow struct Movielens1M <: Movielens
    name::AbstractString = "movielens1m"
    url::AbstractString = "http://files.grouplens.org/datasets/movielens/ml-1m.zip"
    dataset_dir::AbstractString = joinpath(@__DIR__, "..", "..", "dataset", "movielens1m")
end

function load_dataset(d::Movielens1M)
    return load_inter(d), load_user(d), load_item(d)
end

function load_inter(d::Movielens1M)
    csv = CSV.File(
        joinpath(d.dataset_dir, "ratings.dat"),
        delim = "::",
        header = [:userid, :movieid, :rating, :timestamp],
    )
    return csv
end

function load_user(d::Movielens1M)
    csv = CSV.File(
        joinpath(d.dataset_dir, "users.dat"),
        delim = "::",
        header = [:userid, :gender, :age, :occupation, :zipcode],
    )
    return csv
end

function load_item(d::Movielens1M)
    csv = CSV.File(
        joinpath(d.dataset_dir, "movies.dat"),
        delim = "::",
        header = [:movieid, :title, :genres],
    )
    return csv
end


text:
#=
Define data structure to hold results of each iteration of ABC Population Monte Carlo
algorithm.

# Arguments
* `particles::Array` 1 or 2D array containing accepted particles
* `distances::Array` 1 or 2D array containing distances associated with accepted
 particles
* `weights::Array{Float64, 1}` array of weights associated with each accepted particle
* `acceptbw`
* `nsampled`
=#

#ABC Posterior can be over either a single parameter or multivarate
"""
ABC PMC type describing number of parameters in distribution
###### Subtypes: `Univariate`, `Multivariate`
"""
abstract type ParticleDimension end

"""
ABC PMC type for single-parameter distributions
###### Supertype: `ParticleDimension`
"""
type Univariate <: ParticleDimension end

"""
ABC PMC type for multi-parameter distributions

###### Supertype: `ParticleDimension`
"""
type Multivariate <: ParticleDimension end

#Support posteriors based on both singe and multiple distance measures
"""
ABC PMC type describing number of distance metrics used
###### Subtypes: `SingleMeasure`, `MultiMeasure`
"""
abstract type DistanceMeasure end

"""
ABC PMC type describing single distance measure algorithms
###### Supertype: `DistanceMeasure`
"""
type SingleMeasure <: DistanceMeasure end

"""
ABC PMC type describing single distance measure algorithms
###### Supertype: `DistanceMeasure`
"""
type MultiMeasure <: DistanceMeasure end

#Base type for posterior object
"""
Define data structure to hold results of each iteration of ABC Population Monte Carlo
algorithm.

# Fieldnames
* `particles`accepted particles
* `distances::Array{G <: Real}` distances associated with accepted particles
* `weights::StatsBase.AnalyticWeights{Float64, Array{Float64, 1}` importance weights each accepted particle
* `acceptbw::Array{G <: Real, 1} or G <: Real` acceptance acceptbw
* `nsampled::Array{Int64, 1}` number of samples tested before acceptance for each particle
"""
abstract type AbcPmcStep{P <: ParticleDimension, D <: DistanceMeasure} end

#Define alias for easier function definition
UnivariateAbcPmc{D <: DistanceMeasure} = AbcPmcStep{Univariate, D}
MultivariateAbcPmc{D <: DistanceMeasure} = AbcPmcStep{Multivariate, D}

SingleMeasureAbcPmc{P <: ParticleDimension} = AbcPmcStep{P, SingleMeasure}
MultiMeasureAbcPmc{P <: ParticleDimension} = AbcPmcStep{P, MultiMeasure}

#Single parameter and single distance measure
type USAbcPmc{T <: Number, G <: Real} <: AbcPmcStep{Univariate, SingleMeasure}
    particles::Array{T, 1}
    distances::Array{G, 1}
    weights::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}
    acceptbw::G
    nsampled::Array{Int64, 1}
end

#Single parameter and multiple distance measures
type UMAbcPmc{T <: Number, G <: Real} <: AbcPmcStep{Univariate, MultiMeasure}
    particles::Array{T, 1}
    distances::Array{G, 2}
    weights::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}
    acceptbw::Array{G, 1}
    nsampled::Array{Int64, 1}
end

#Multiple parameters and single distance measure
type MSAbcPmc{T <: Number, G <: Real} <: AbcPmcStep{Multivariate, SingleMeasure}
    particles::Array{T, 2}
    distances::Array{G, 1}
    weights::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}
    acceptbw::G
    nsampled::Array{Int64, 1}
end

#Multiple parameters and multiple distance measures
type MMAbcPmc{T <: Number, G <: Real} <: AbcPmcStep{Multivariate, MultiMeasure}
    particles::Array{T, 2}
    distances::Array{G, 2}
    weights::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}
    acceptbw::Array{G, 1}
    nsampled::Array{Int64, 1}
end

#Constructor selecting type based on inputs
AbcPmc{T <: Number, G <: Real}(p::Array{T, 1}, d::Array{G, 1}, w::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}, t::G, ts::Array{Int64, 1}) = USAbcPmc(p, d, w, t, ts)
AbcPmc{T <: Number, G <: Real}(p::Array{T, 1}, d::Array{G, 2}, w::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}, t::Array{G, 1}, ts::Array{Int64, 1}) = UMAbcPmc(p, d, w, t, ts)
AbcPmc{T <: Number, G <: Real}(p::Array{T, 2}, d::Array{G, 1}, w::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}, t::G, ts::Array{Int64, 1}) = MSAbcPmc(p, d, w, t, ts)
AbcPmc{T <: Number, G <: Real}(p::Array{T, 2}, d::Array{G, 2}, w::StatsBase.AnalyticWeights{Float64, Float64, Array{Float64, 1}}, t::Array{G, 1}, ts::Array{Int64, 1}) = MMAbcPmc(p, d, w, t, ts)

function copy{A <: AbcPmcStep}(abcpmc::A)
    return AbcPmc(copy(abcpmc.particles), copy(abcpmc.distances), StatsBase.AnalyticWeights(abcpmc.weights.values), copy(abcpmc.acceptbw), copy(abcpmc.nsampled))
end

function StatsBase.sample{A <: UnivariateAbcPmc}(abcpmc::A)
    return StatsBase.sample(abcpmc.particles, abcpmc.weights)
end

function StatsBase.sample{A <: MultivariateAbcPmc}(abcpmc::A)
    row = StatsBase.sample(1:size(abcpmc.particles, 1), abcpmc.weights)
    return abcpmc.particles[row, :]
end


text:
# <md>
# chimera states in a small ring of electro-optic oscillators

# In[]
# loading necessary packages
#using DifferentialEquations
using Plots; pyplot()
using LightGraphs
#using GraphPlot
using LinearAlgebra
using Statistics
using BenchmarkTools
using ProgressMeter
using LaTeXStrings
using DelimitedFiles

# In[]
n = 6   # size of the ring
k = 2    # number of neighbors coupled (in each direction)
σ = -.7    # coupling strength
β = 2*π/3 - 4*σ   # parameter of the electro-optic oscillators
ξ = π/6     # offset parameter
r = 0     # rewiring probability
ring = watts_strogatz(n, 2*k, r)
adj = Matrix(adjacency_matrix(ring))

# In[]
# mapping function for the electro-optic oscillators
function I(x)
    return (1-cos(x))/2
end

# In[]
# simulating the full system to study chimera states
T = 1000000 # total simulation time
x = zeros(n,T) # oscillator states
x[:,1] = 2*π*[.1,.5,.5,.7,.5,.5] + 1e-3*rand(n,1) # initial condition close to a chimera state
p = Progress(T, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow) # progress bar
@time for t in 2:T
    for i in 1:n
        x[i,t] = β*I(x[i,t-1]) + σ*adj[i,:]⋅I.(x[:,t-1]) + ξ
    end
    x[:,t] = mod.(x[:,t],2*π)
    next!(p)
end

# In[]
# plot showing the states of the oscillators in a chimera state (coherent cluster: green; incoherent cluster: orange and magenta)
t = 500 # show the final t iterations
scatter(T-t:T,x[1,T-t:T],marker=(:circle,3,.6,:magenta,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
scatter!(T-t:T,x[4,T-t:T],marker=(:circle,3,.6,:orange,stroke(0,.2,:black)))
for i = [2,3,5,6]
    scatter!(T-t:T,x[i,T-t:T],marker=(:circle,3,.6,stroke(1,.2,:white)))
end
title!("Time Series for Chimeras")
xaxis!(L"t",font(15))
yaxis!(L"x_i",font(15),(-.1,2*π))
yticks!([0,π,2*π],[L"0",L"\pi",L"2\pi"])

# In[]
# plot showing the input from the incoherent cluster to the coherent cluster
@show mean(σ*I.(x[1,:])+σ*I.(x[4,:])) # mean of the input
@show std(σ*I.(x[1,:])+σ*I.(x[4,:])) # s.t.d. of the input
t = 500 # show the final t iterations
scatter(T-t:T,σ*I.(x[1,T-t:T])+σ*I.(x[4,T-t:T]),marker=(:circle,3,.6,:red,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Time Series of the Driving")
xaxis!(L"t",font(15))
yaxis!("Input",font(15))

# In[]
# simulating the two-oscillator reduced system
T = 1000000
w = zeros(2,T) # oscillator states
w[:,1] = 2*π*rand(2,1) # random initial condition
p = Progress(T, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
@time for t in 2:T
    # quotient of the 4-node cluster
    w[1,t] = β*I(w[1,t-1]) + 2*σ*(I(w[2,t-1])) + 2*σ*(I(w[1,t-1]))  + ξ
    # quotient of the 2-node cluster
    w[2,t] = β*I(w[2,t-1]) + 4*σ*(I(w[1,t-1])) + ξ
    w[:,t] = mod.(w[:,t],2*π)
    next!(p)
end

# In[]
# plot showing the states of the oscillators in the two-oscillator reduced system (which is unstable in the full system)
t = 500 # show the final t iterations
scatter(T-t:T,w[1,T-t:T],marker=(:circle,3,.6,:orange,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
scatter!(T-t:T,w[2,T-t:T],marker=(:circle,3,.6,:green,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Reduced System")
xaxis!(L"t",font(15))
yaxis!(L"x_i",font(15),(-.1,2*π))
yticks!([0,π,2*π],[L"0",L"\pi",L"2\pi"])

# In[]
# plot showing the input from one cluster to the other
@show mean(2*σ*I.(w[2,:])) # mean of the input
@show std(2*σ*I.(w[2,:])) # s.t.d. of the input
t = 500 # show the final t iterations
scatter(T-t:T,2*σ*I.(w[2,T-t:T]),marker=(:circle,3,.6,:red,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Time Series of the Driving")
xaxis!(L"t",font(15))
yaxis!("Input",font(15))


text:
import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        wait(cond)
        msg = codeunits(unsafe_string(r_ecb[].message))
        println("message ", msg)
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)


text:
@testset "Experiments" begin
	using DrWatson
	@quickactivate
	using Test
	using GenerativeAD

	# first create dummy data
	tstpath = abspath(joinpath(dirname(pathof(GenerativeAD)), "../test/misc"))
	run(`mkdir $tstpath`)
	run(`touch $tstpath/model_a=4_b=0.001.bson`)
	run(`touch $tstpath/a=4_b=0.001.bson`)
	run(`touch $tstpath/a=4_b=0.001_c=bar.bson`)
	run(`touch $tstpath/a=4_b=0.001_d=foo.bson`)

	# now test check_params
	@test !GenerativeAD.check_params(tstpath, (a=4,))
	@test !GenerativeAD.check_params(tstpath, (a=4, b=0.001))
	@test GenerativeAD.check_params(tstpath, (a=4, b=0.002))
	@test !GenerativeAD.check_params(tstpath, (a=4, b=0.001, c="bar"))
	@test GenerativeAD.check_params(tstpath, (a=4, b=0.001, c="foo"))
	@test !GenerativeAD.check_params(tstpath, (a=4, b=0.001, d="foo"))
	@test GenerativeAD.check_params(tstpath, (a=4, b=0.001, g="a"))

	# cleanup
	run(`rm -r $tstpath`)
end


text:
using PolynomialTestSystems

# see https://www.juliahomotopycontinuation.org/examples/tritangents/
f = equations(tritangents())
res = solve(f)
println("# nonsingular: ", nnonsingular(res)) # should be 720

using Bertini
@time bertini(f; MPTYPE=0)


text:
function NewMarkCellsBySteps(r)
    side = Ost
    counter_start = 0
    B = Array{Int, 1}(undef, 1000)
    B, counter_start = movetothecorner(r, West, B)

    counter = movementsandcountmoves(r, Ost)
    movements(r, West)

    markrowbyfixmoves(r, side, counter)
    putmarker!(r)
    movements(r, inverse(side))
    counter -= 1
    while isborder(r, Nord) == false
        move!(r, Nord)
        t_counter = markrowbyfixmoves(r, side, counter)
        if t_counter >= 0
            putmarker!(r)
        end
        movementsbyfixmoves(r, inverse(side), counter - t_counter)
        counter -= 1
    end

    movements(r, Sud)
    movements(r, West)
    returnfromthecorner(r, inverse(West), B, counter_start)
end

function markrowbyfixmoves(r, side, counter)
    while counter > 0
        putmarker!(r)
        counter = trymove!(r, side, counter)
    end
    return counter
end

function movementsbyfixmoves(r, side, counter)
    while counter > 0
        counter = trymove!(r, side, counter)
    end
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function movements(r, side)
    while isborder(r, side) == false
        move!(r, side)
    end
end

function movementsandcountmoves(r, side)
    moves = 0
    while isborder(r, side) == false
        move!(r, side)
        moves += 1
    end
    return moves
end

function trymove!(r::Robot, side::HorizonSide, counter::Int)
    x = 0
    if isborder(r,side)
        while isborder(r, side)
            move!(r,HorizonSide((Int(side) + 1) % 4))
            x = x + 1
        end
        move!(r,side)
        counter -= 1
        while isborder(r,HorizonSide((Int(side) + 3) % 4))
            move!(r,side)
            counter -= 1
        end
        while x != 0
            move!(r,HorizonSide((Int(side) + 3) % 4))
            x = x - 1
        end
        return counter
    else
        move!(r,side)
        return counter - 1
    end
end

function movetothecorner(r, side, B)
    counter = 1
    while ((!isborder(r, side)) || (!isborder(r, HorizonSide((Int(side) + 1) % 4))))
        if !isborder(r, side)
            move!(r, side)
            B[counter] = 1
            counter += 1
        end
        if !isborder(r, HorizonSide((Int(side) + 1) % 4))
            move!(r, HorizonSide((Int(side) + 1) % 4))
            B[counter] = 2
            counter += 1
        end
    end
    return B, counter
end

function returnfromthecorner(r, side, B, counter)
    for i in 1:(counter - 1)
        if B[counter - i] == 1
            move!(r, side)
        end
        if B[counter - i] == 2
            move!(r, HorizonSide((Int(side) + 1) % 4))
        end
    end
end


text:
import DelimitedFiles, DSP

@testset "Rotation correlation" begin
    @testset "Correlation" begin
        let Ts = (Float32, Float64), N = 1000
            for T in Ts
                a, b = rand(T, N), rand(T, N)
                xc = similar(a, 2N-1)
                work_arrays = SeisSplit._correlation_arrays(T, N)
                @test DSP.xcorr(a, b) ≈ SeisSplit._xcorr_fd!(xc, a, b, work_arrays...)
            end
        end
    end

    @testset "Xcorr" begin
        let t = read_test_data(), nphi = 180
            e, n = t
            win1, win2 = n.meta.SAC_user0, n.meta.SAC_user2
            dts, phis, xcorr_map, phi_best, dt_best =
                SeisSplit.rotation_correlation(n, e, win1, win2, nphi=nphi÷2)
            s = splitting(n, e, win1, win2, nphi=nphi)
            @test s.xcorr_phi == phis
            @test s.xcorr_dt == dts
            @test s.xcorr_map == xcorr_map
            @test s.xcorr_phi_best == phi_best
            @test s.xcorr_dt_best == dt_best
            @test phi_best ≈ 42 atol=1
            @test dt_best ≈ 1.35 atol=step(s.xcorr_phi)
        end
    end

    @testset "Xcorr map" begin
        # Compare to SHEBA-computed answer, though SHEBA appears
        # to include an extra point before 0 s dt, and takes the
        # absolute correlation coefficient.
        let t = read_test_data(), nphi = 180
            e, n = t
            dts, phis, xc, phi_best, dt_best =
                rotation_correlation(n, e, n.meta.SAC_user0, n.meta.SAC_user2, nphi=nphi)
            sheba_map = DelimitedFiles.readdlm(joinpath(@__DIR__, "data", "sheba.xc"))
            sheba_map = sheba_map'[2:end,1:180]
            xc = xc[1:2:80,:] .|> abs
            @test sheba_map ≈ xc rtol=0.01
        end
    end
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
function getrice2010parameters(filename)
    p = Dict{Symbol,Any}()

    T = 60
    # p[:timesteps] = 1:T # Time periods (5 years per period); not currently used in constructrice
    regions = ["US", "EU", "Japan", "Russia", "Eurasia", "China", "India", "MidEast", "Africa", "LatAm", "OHI", "OthAsia"]

    # Open RICE_2010 Excel File to Read Parameters
    f = readxlsx(filename)

    # Time Step
    # p[:tstep] = 10 # Years per Period; not currently used in constructrice
    # p[:dt] =  10 # Time step parameter for model equations; not currently used in constructrice

    # If optimal control
    ifopt = true # Indicator where optimized is 1 and base is 0

    # Preferences
    p[:elasmu] =  getparam_single(f, "B18:B18", regions) # Elasticity of MU of consumption
    # p[:prstp] =  getparam_single(f, "B15:B15", regions) # Rate of Social Time Preference; not currently used in constructrice

    # Population and technology

    # Capital elasticity in production function
    p[:gama] = 0.300
    p[:dk]  = getparam_single(f, "B8:B8", regions) # Depreciation rate on capital (per year)
    p[:k0] = getparam_single(f, "B11:B11", regions) #Initial capital
    # p[:miu0] = getparam_single(f, "B103:B103", regions) # Initial emissions control rate for base case 2010; not currently used in constructrice
    # p[:miubase] = getparam_timeseries(f, "B103:BI103", regions, T) # Optimized emission control rate results from RICE2010 (base case); duplicate line to p[:MIU]

    # Carbon cycle

    # Initial Conditions
    p[:mat0] =  787.0 # Initial Concentration in atmosphere 2010 (GtC)
    p[:mat1] =  829.0
    p[:mu0] = 1600. # Initial Concentration in upper strata 2010 (GtC)
    p[:ml0] =  10010. # Initial Concentration in lower strata 2010 (GtC)

    # Carbon cycle transition matrix

    # Flow paramaters
    p[:b12] = 12.0/100 # Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b23] = 0.5/100 # Carbon cycle transition matrix shallow to deep ocean

    # Parameters for long-run consistency of carbon cycle
    p[:b11] = 88.0/100 # Carbon cycle transition matrix atmosphere to atmosphere
    p[:b21] = 4.704/100 # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
    p[:b22] = 94.796/100 # Carbon cycle transition matrix shallow ocean to shallow oceans
    p[:b32] = 0.075/100 # Carbon cycle transition matrix deep ocean to shallow ocean
    p[:b33] = 99.925/100 # Carbon cycle transition matrix deep ocean to deep oceans

    # Climate model parameters
    p[:t2xco2] = 3.2 # Equilibrium temp impact (oC per doubling CO2)
    fex0 = 0.83 # 2010 forcings of non-CO2 GHG (Wm-2)
    fex1 = 0.3 # 2100 forcings of non-CO2 GHG (Wm-2)
    p[:tocean0] = .0068 #  Initial lower stratum temp change (C from 1900)
    p[:tatm0] = 0.83 # Initial atmospheric temp change 2005 (C from 1900)
    p[:tatm1] = 0.98 # Initial atmospheric temp change 2015 (C from 1900)

    # Transient TSC Correction ("Speed of Adjustment Parameter")
    p[:c1] = 0.208 # Climate equation coefficient for upper level
    p[:c3] = 0.310 # Transfer coefficient upper to lower stratum
    p[:c4] = 0.05 # Transfer coefficient for lower level
    p[:fco22x] = 3.8 # Forcings of equilibrium CO2 doubling (Wm-2)

    # Climate damage parameters
    p[:a1] = getparam_single(f, "B24:B24", regions) # Damage intercept
    p[:a2] = getparam_single(f, "B25:B25", regions) # Damage quadratic term
    p[:a3] = getparam_single(f, "B26:B26", regions) # Damage exponent

    # Welfare Weights
    alpha0 = transpose(f["Data"]["B359:BI370"]) # Read in alpha
    p[:alpha] = convert(Array{Float64}, alpha0) # Convert to type used by Mimi

    # Abatement cost
    # Exponent of control cost function
    p[:expcost2] = getparam_single(f, "B38:B38", regions)

    # Availability of fossil fuels
    # Maximum cumulative extraction fossil fuels (GtC)
    # p[:fosslim] = 6000; not currently used in constructrice.

    # Scaling parameters
    # Multiplicative scaling coefficient
    p[:scale1] = getparam_single(f, "B52:B52", regions)
    p[:scale1] = getparam_single(f, "B52:B52", regions)

    # Additive scaling coefficient (combines two additive scaling coefficients from RICE for calculating utility with welfare weights)
    scale2 = Array{Float64}(undef, length(regions))
    for (i,r) in enumerate(regions)
        data = f[r]["B53:C53"]
        scale2[i] = data[1] - data[2]
    end
    p[:scale2] = scale2

    # p[:savebase] = getparam_timeseries(f, "B97:BI97", regions, T) # Optimized savings rate in base case for RICE2010; not currently used in constructrice
    # p[:optlrsav] = getparam_single(f, "BI97:BI97", regions) # Optimized savings rate in base case for RICE2010 for last period (fraction of gross output); not currently used in constructrice
    p[:l] = getparam_timeseries(f, "B56:BI56", regions, T) # Level of population and labor
    p[:al] = getparam_timeseries(f, "B20:BI20", regions, T) # Level of total factor productivity
    p[:sigma] = getparam_timeseries(f, "B40:BI40", regions, T) # CO2-equivalent-emissions output ratio
    p[:pbacktime] = getparam_timeseries(f, "B36:BI36", regions, T) # Backstop price
    p[:cost1] = getparam_timeseries(f, "B31:BI31", regions, T) # Adjusted cost for backstop
    regtree = getparam_timeseries(f, "B43:BI43", regions, T) # Regional Emissions from Land Use Change
    p[:rr] = getparam_timeseries(f, "B17:BI17", regions, T) # Social Time Preference Factor

    # Global Emissions from Land Use Change (Sum of regional emissions for land use change in RICE model)
    etree = Array{Float64}(undef, T)
    for i = 1:T
        etree[i] = sum(regtree[i,:])
    end
    p[:etree] = etree

    # Exogenous forcing for other greenhouse gases
    forcoth =  Array{Float64}(undef, 60)
    data = f["Global"]["B21:BI21"]
    for i=1:T
        forcoth[i] = data[i]
    end
    p[:forcoth] = forcoth

    # Fraction of emissions in control regime
    p[:partfract] = ones(60, length(regions))

    # Savings Rate (base case RICE2010)
    p[:S] = getparam_timeseries(f, "B97:BI97", regions, T)

    # MIU (base case RICE2010)
    p[:MIU] = getparam_timeseries(f, "B103:BI103", regions, T)

    # SEA LEVEL RISE PARAMETERS
    p[:slrmultiplier] = getparam_single(f, "B49:B49", regions) # Multiplier for SLR
    p[:slrelasticity] = getparam_single(f, "C49:C49", regions) # SLR elasticity of substitution
    p[:slrdamlinear] = getparam_single(f, "B48:B48", regions) # SLR damage parameter (linear)
    p[:slrdamquadratic] = getparam_single(f, "C48:C48", regions) # SLR damage parameter (quadratic)

    # Thermal Expansion
    p[:therm0] = f["SLR"]["B9:B9"][1] # Thermal Expansion initial conditions (SLR per decade)
    p[:thermadj] = f["SLR"]["B8:B8"][1] # Thermal Expansion adjustment rate/calibration (per decade)
    p[:thermeq] = f["SLR"]["B7:B7"][1] # Thermal Expansion equilibrium (m/degree C)

    # Glaciers and Small Ice Caps (GSIC)
    p[:gsictotal] = f["SLR"]["B12:B12"][1] # GSIC total ice (SLR equivalent in meters)
    p[:gsicmelt] = f["SLR"]["B13:B13"][1] # GSIC melt rate (meters/year/degree C)
    p[:gsicexp] = f["SLR"]["B11:B11"][1] # GSIC exponent (assumed)
    # p[:gsiceq] = f["SLR"]["B14:B14"][1] # GSIC equilibrium temperature (degrees C relative to global T of -1 degree C from 2000)

    # Greenland Ice Sheet (GIS)
    p[:gis0] = f["SLR"]["B18:B18"][1] # GIS initial ice volume (meters)
    p[:gismelt0] = f["SLR"]["B17:B17"][1] # GIS initial melt rate (mm per year)
    p[:gismeltabove] = f["SLR"]["B20:B20"][1] # GIS melt rate above threshold (mm/year/degree C)
    p[:gismineq] = f["SLR"]["B19:B19"][1] # GIS minimum equilibrium temperature (degrees C)
    p[:gisexp] = f["SLR"]["B21:B21"][1] # GIS exponent on remaining

    # Antarctic Ice Sheet (AIS)
    p[:aismelt0] = f["SLR"]["B24:B24"][1] # AIS initial melt rate (mm/year)
    p[:aismeltlow] = f["SLR"]["B26:B26"][1] # AIS melt rate lower (mm/year/degrees C) [T < 3 degrees C over Antarctic)
    p[:aismeltup] = f["SLR"]["B27:B27"][1] # AIS melt rate upper (mm/year/degrees C) [T = 8 degrees C over Antarctic]
    p[:aisratio] = f["SLR"]["B29:B29"][1] # AIS ratio t-ant/t-glob
    p[:aisinflection] = f["SLR"]["B28:B28"][1] # AIS inflection point (degrees C)
    p[:aisintercept] = f["SLR"]["B25:B25"][1] # AIS intercept (mm/year)
    p[:aiswais] = f["SLR"]["B31:B31"][1] # AIS total remaining ice volume (m) for WAIS
    p[:aisother] = f["SLR"]["B32:B32"][1] # AIS total remaining ice volume (m) for other AIS

    return p
end


text:
using Knet: Knet, AutoGrad

module AutoML

using ArgParse

export
    Knet,
    AutoGrad,
    Model,
    Config,
    savemodel,
    loadmodel,
    preprocess,
    preprocess2,
    preparedata,
    build,
    hyperoptimization,
    train,
    partialtrain,
    predictdata,
    crossvalidate,
    getbatches,
    csv2data,
    readtext,
    splitdata,
    house_rentals,
    splice_junction,
    cifar100,
    imdb_movie_review,
    quora_questions,
    default_of_credit,
    prediction_of_return,
    train_house_rentals,
    train_gene_sequences,
    train_cifar_100,
    train_imdb,
    train_quora,
    train_default_of_credit,
    train_prediction_of_return

const DIR = @__DIR__
const SAVEDIR = abspath(joinpath(DIR, "..", "models"))
const DATADIR = abspath(joinpath(DIR, "..", "data"))

include("config.jl")
include("layers.jl")
include("data.jl")
include("utils.jl")
include("model.jl")
include("build.jl")
include("clean.jl")
include("preprocess.jl")
include("hyperopt.jl")
include("train.jl")
include("examples.jl")
include("resnet.jl")
include("bert.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "Automated Machine Learning Tool"

    @add_arg_table s begin
        ("--config"; default=nothing; help="Yaml file for config")
        ("--train"; default=nothing; help="Training dataset file.")
        ("--model_name"; default="model"; help="Model name")
        ("--model_path"; default=nothing; help="Pretrained model path")
        ("--predict"; default=nothing; help="Test dataset file.")
        ("--epochs"; default="1"; help="Number of epochs")
        ("--optimize"; default="true"; help="Hyperparameter optimization")
    end

    isa(args, AbstractString) && (args=split(args))
    println(s.description)
    o = parse_args(args, s; as_symbols=true)

    istrain = o[:train] != nothing
    ispredict = o[:predict] != nothing

    if istrain
        o[:config] == nothing && throw("Config file is missing")

        yamlpath = isabspath(o[:config]) ? o[:config] :
            joinpath(dirname(DIR), o[:config])
        trainpath = isabspath(o[:train]) ? o[:train] :
            joinpath(dirname(DIR), o[:train])
        epochs = parse(Int, o[:epochs])
        optimize = parse(Bool, o[:optimize])
        modelname = o[:model_name]

        inputfeatures, outputfeatures = readyaml(yamlpath)
        model = Model(inputfeatures, outputfeatures; name=modelname)
        model, dtrn = train(model, trainpath;
                            epochs=epochs,
                            showprogress=false,
                            issave=true,
                            optimize=optimize)

        println("Error: ", model.model(dtrn.x, dtrn.y))
        println("Accuracy: ", accuracy(model.model, dtrn))

        return
    elseif ispredict
        o[:model_path] == nothing && throw("Model path is missing")

        modelpath = isabspath(o[:model_path]) ? o[:model_path] :
            joinpath(dirname(DIR), o[:model_path])
        testpath = isabspath(o[:predict]) ? o[:predict] :
            joinpath(dirname(DIR), o[:predict])

        model = loadmodel(modelpath)
        dtst = getbatches(model, testpath; showtime=false)
        println("Error: ", model.model(dtst.x, dtst.y))
        println("Accuracy: ", accuracy(model.model, dtst))
    else
        throw("Use --train or --predict parameters to use package")
    end
end

basename(PROGRAM_FILE) == "AutoML.jl" && main(ARGS)

end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
using DEMC
using Test

@testset "DEMC.jl" begin
    # Write your own tests here.
    println("test with norm pdf in serial computation")
    include("example_normpdf.jl")
end


text:
using PyPlot
using NetCDF

print("Examing mapplot_something in a global map")
plot_problem = AModel_Import(filename = "data/example.nc")
sst = ncread(plot_problem.data_based_on, "sst")
mapplot_something(plot_problem, sst)
savefig("An-Example-of-global-plot.png")

print("Examing mapplot_something in a ENSO area")
ENSO_area= Give_An_ENSO_area(plot_problem)
mapplot_something(ENSO_area, sst)
savefig("An-Example-of-ENSO-plot.png")


text:
function iscyclic(a::Vector)
    for i=1:length(a)-1
      if a == circshift(a, i)
        return true
      end
    end
    # else
    return false
end
export iscyclic

function isunique(a)::Bool
    return length(a) == length(unique(a))
end
export isunique

function mapArray(arr::AbstractArray, low, high)
"""
  map array values to region [0,1] where:
  * arr<=low -> 0
  * arr>=high -> 1
  * arr∈[low,high] -> values are linearily mapped from within [low,high] to [0,1]

  * quick test:
  *   mapArray(rand(3,3,3),0.25,0.75)
"""
  _arr, _low, _high = promote(arr, low, high)
  _arr = clamp.(_arr, _low, _high)
  return (_arr-_low)/(_high-_low)
end
export mapArray


text:
type RenderWindow
    ptr::Ptr{Void}
    _view::View

    function RenderWindow(ptr::Ptr{Void})
        w = new(ptr)
        w
        # finalizer(w, destroy)
        # w
    end
end


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

name = "OptimPack"
version = v"3.1.0"

# Collection of sources required to build OptimPAck
sources = [
    "https://github.com/emmt/OptimPack/releases/download/v$(version)/optimpack-$(version).tar.gz" =>
    "fa1efdbebf6efa42a1d6b4f6223c7ca0b871fad5c0b53f7dc9be296e4c766190",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/optimpack-*
./configure --prefix=$prefix --host=$target --enable-shared=yes --enable-static=no
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    # This is needed in order to build the shared library on Windows
    FLAGS+=(LDFLAGS="-no-undefined")
fi
make -j${nproc} "${FLAGS[@]}"
make install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libopk",    :libopk),
    LibraryProduct("libcobyla", :libcobyla),
    LibraryProduct("libbobyqa", :libbobyqa),
    LibraryProduct("libnewuoa", :libnewuoa)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using Zygote

LEARNING_RATE = 0.5

miditofreq(n) = 440 * 2^((n-69)/12)

overtones(f, n = 4) = map((x) -> (f * x, 1 / x), 1:n)

function loudness(amp)
    p_ref = 20 * 10^-6 # FIXME
    p_e = amp / sqrt(2)
    spl = 20 * log10(p_e / p_ref)
    (1/16) * 2^(spl/10)
end

# https://sethares.engr.wisc.edu/comprog.html
function dissonance(freq1, freq2, amp1, amp2)
    b_1 = 3.6
    b_2 = 5.75
    d_max = 0.24
    s_1 = 0.0207
    s_2 = 18.96

    l12 = min(loudness(amp1), loudness(amp2))

    f_1, f_2 = freq1 < freq2 ? [freq1, freq2] : [freq2, freq1]
    s = d_max / (s_1 * f_1 + s_2)
    x = s * (f_2 - f_1)
    l12 * (exp(-b_1 * x) - exp(-b_2 * x))
end

function total_dissonance(freq1, freq2)
    overtones1 = overtones(freq1)
    overtones2 = overtones(freq2)
    s = 0
    for (f1, a1) in overtones1
        for (f2, a2) in overtones2
            s += dissonance(f1, f2, a1, a1)
        end
    end
    s
end

function adjust(pitches)
    notes, bends = map(collect, zip(pitches...))
    freqs = miditofreq.(notes)
    function loss(xs)
        adjusted_freqs = freqs .* centtoratio.(bends + xs)
        sum(
            total_dissonance(f1, f2)
            for f1 in adjusted_freqs
            for f2 in adjusted_freqs
            if f1 != f2
        )
    end
    guess = -LEARNING_RATE * loss'(zero(bends))

    bends + guess
end

max_stars_repo_name: miksu/FluidIntonation

max_issues_repo_name: miksu/FluidIntonation

max_forks_repo_name: miksu/FluidIntonation


text:
const _index = [1,2,3,4,5,6,0,0,7,8,9,10,11,12,13,14,0,0,0,0,15,16,17,18]

const _amino = AminoAcid[0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46;
 0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46  0x46;
 0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x2a  0x4c;
 0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x2a  0x53  0x53;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53;
 0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59;
 0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59  0x59;
 0x2a  0x2a  0x2a  0x2a  0x2a  0x51  0x2a  0x2a  0x2a  0x2a  0x2a  0x59  0x2a  0x2a  0x2a  0x2a  0x2a  0x2a;
 0x2a  0x2a  0x2a  0x2a  0x2a  0x51  0x2a  0x2a  0x2a  0x2a  0x2a  0x2a  0x51  0x4c  0x2a  0x4c  0x2a  0x2a;
 0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43;
 0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43  0x43;
 0x2a  0x57  0x57  0x57  0x57  0x2a  0x57  0x43  0x2a  0x2a  0x57  0x57  0x2a  0x2a  0x57  0x2a  0x2a  0x57;
 0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57  0x57;
 0x4c  0x4c  0x54  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c;
 0x4c  0x4c  0x54  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c;
 0x4c  0x4c  0x54  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c;
 0x4c  0x4c  0x54  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x53  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c  0x4c;
 0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50;
 0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50;
 0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50;
 0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50  0x50;
 0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48;
 0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48  0x48;
 0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51;
 0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51  0x51;
 0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52;
 0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52;
 0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52;
 0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52  0x52;
 0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49;
 0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49  0x49;
 0x49  0x4d  0x4d  0x49  0x4d  0x49  0x49  0x49  0x49  0x49  0x4d  0x49  0x49  0x49  0x4d  0x49  0x49  0x49;
 0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d  0x4d;
 0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54;
 0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54;
 0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54;
 0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54  0x54;
 0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e;
 0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e  0x4e;
 0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4e  0x4b  0x4b  0x4b  0x4b  0x4e  0x4b  0x4b  0x4e  0x4b  0x4b  0x4b;
 0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b  0x4b;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53;
 0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53  0x53;
 0x52  0x2a  0x52  0x52  0x53  0x52  0x53  0x52  0x52  0x52  0x47  0x53  0x52  0x52  0x53  0x52  0x52  0x53;
 0x52  0x2a  0x52  0x52  0x53  0x52  0x53  0x52  0x52  0x52  0x47  0x53  0x52  0x52  0x53  0x52  0x52  0x4b;
 0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56;
 0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56;
 0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56;
 0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56  0x56;
 0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41;
 0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41;
 0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41;
 0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41  0x41;
 0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44;
 0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44  0x44;
 0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45;
 0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45  0x45;
 0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47;
 0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47;
 0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47;
 0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47  0x47]

const _start = convert(BitArray,[false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false   true  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
  true  false  false   true   true  false  false  false   true  false   true  false  false  false  false  false  false   true;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
  true  false  false   true  false  false  false  false   true   true  false  false  false  false  false  false  false   true;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false   true  false   true   true  false  false  false   true  false  false  false  false  false  false  false   true  false;
 false   true  false   true   true  false  false  false   true  false  false  false  false  false  false  false  false  false;
 false   true   true   true   true  false  false  false   true  false   true  false  false  false  false  false  false  false;
  true   true   true   true   true   true   true   true   true   true   true   true   true   true   true   true   true   true;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false   true  false   true   true  false   true  false   true  false   true  false  false  false   true  false   true   true;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false;
 false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false  false])

const _codons = Dict{Vector{Nucleotide},Int}( nucleotide(['G','C','U'])=>53,nucleotide(['G','C','t'])=>53,nucleotide(['u','U','G'])=>4,nucleotide(['C','C','t'])=>21,nucleotide(['U','G','c'])=>14,
nucleotide(['u','U','U'])=>1,nucleotide(['C','u','U'])=>17,nucleotide(['C','A','T'])=>25,nucleotide(['c','T','T'])=>17,nucleotide(['c','T','G'])=>20,
nucleotide(['A','T','A'])=>35,nucleotide(['A','C','g'])=>40,nucleotide(['T','C','g'])=>8,nucleotide(['G','A','t'])=>57,nucleotide(['U','C','a'])=>7,
nucleotide(['U','C','G'])=>8,nucleotide(['C','G','c'])=>30,nucleotide(['G','A','U'])=>57,nucleotide(['t','C','A'])=>7,nucleotide(['t','G','C'])=>14,
nucleotide(['U','G','G'])=>16,nucleotide(['C','g','G'])=>32,nucleotide(['G','t','A'])=>51,nucleotide(['T','G','C'])=>14,nucleotide(['U','G','C'])=>14,
nucleotide(['g','C','T'])=>53,nucleotide(['T','G','a'])=>15,nucleotide(['C','C','a'])=>23,nucleotide(['G','C','c'])=>54,nucleotide(['C','A','a'])=>27,
nucleotide(['C','t','T'])=>17,nucleotide(['A','A','C'])=>42,nucleotide(['A','u','A'])=>35,nucleotide(['c','A','C'])=>26,nucleotide(['a','T','G'])=>36,
nucleotide(['A','A','a'])=>43,nucleotide(['u','U','C'])=>2,nucleotide(['G','A','G'])=>60,nucleotide(['G','a','G'])=>60,nucleotide(['g','A','C'])=>58,
nucleotide(['c','C','C'])=>22,nucleotide(['C','U','u'])=>17,nucleotide(['T','g','T'])=>13,nucleotide(['c','C','U'])=>21,nucleotide(['C','G','G'])=>32,
nucleotide(['a','T','C'])=>34,nucleotide(['A','G','t'])=>45,nucleotide(['A','C','a'])=>39,nucleotide(['c','A','U'])=>25,nucleotide(['a','C','G'])=>40,
nucleotide(['G','T','A'])=>51,nucleotide(['G','C','a'])=>55,nucleotide(['A','C','u'])=>37,nucleotide(['c','T','C'])=>18,nucleotide(['c','G','U'])=>29,
nucleotide(['a','T','A'])=>35,nucleotide(['T','t','T'])=>1,nucleotide(['A','U','U'])=>33,nucleotide(['g','T','G'])=>52,nucleotide(['C','T','G'])=>20,
nucleotide(['C','a','U'])=>25,nucleotide(['a','C','C'])=>38,nucleotide(['G','g','T'])=>61,nucleotide(['u','A','C'])=>10,nucleotide(['a','C','A'])=>39,
nucleotide(['U','c','U'])=>5,nucleotide(['G','U','u'])=>49,nucleotide(['A','A','A'])=>43,nucleotide(['G','T','a'])=>51,nucleotide(['G','T','G'])=>52,
nucleotide(['a','U','U'])=>33,nucleotide(['C','C','C'])=>22,nucleotide(['U','C','g'])=>8,nucleotide(['C','G','t'])=>29,nucleotide(['A','t','G'])=>36,
nucleotide(['G','G','t'])=>61,nucleotide(['A','C','G'])=>40,nucleotide(['a','A','A'])=>43,nucleotide(['A','T','T'])=>33,nucleotide(['G','G','u'])=>61,
nucleotide(['C','U','A'])=>19,nucleotide(['g','C','A'])=>55,nucleotide(['G','G','c'])=>62,nucleotide(['U','U','C'])=>2,nucleotide(['A','a','G'])=>44,
nucleotide(['G','a','C'])=>58,nucleotide(['U','C','A'])=>7,nucleotide(['A','u','G'])=>36,nucleotide(['A','U','u'])=>33,nucleotide(['a','C','U'])=>37,
nucleotide(['C','T','g'])=>20,nucleotide(['T','T','g'])=>4,nucleotide(['T','C','G'])=>8,nucleotide(['T','t','A'])=>3,nucleotide(['U','a','U'])=>9,
nucleotide(['C','C','U'])=>21,nucleotide(['A','a','C'])=>42,nucleotide(['g','C','U'])=>53,nucleotide(['c','G','T'])=>29,nucleotide(['C','g','U'])=>29,
nucleotide(['C','t','A'])=>19,nucleotide(['A','C','c'])=>38,nucleotide(['t','A','G'])=>12,nucleotide(['A','A','T'])=>41,nucleotide(['A','g','A'])=>47,
nucleotide(['C','g','C'])=>30,nucleotide(['G','A','a'])=>59,nucleotide(['G','G','g'])=>64,nucleotide(['T','T','G'])=>4,nucleotide(['G','c','A'])=>55,
nucleotide(['A','T','t'])=>33,nucleotide(['C','C','u'])=>21,nucleotide(['A','T','c'])=>34,nucleotide(['A','U','c'])=>34,nucleotide(['G','C','C'])=>54,
nucleotide(['U','G','a'])=>15,nucleotide(['t','C','G'])=>8,nucleotide(['t','A','C'])=>10,nucleotide(['G','U','g'])=>52,nucleotide(['C','c','T'])=>21,
nucleotide(['C','C','g'])=>24,nucleotide(['U','G','u'])=>13,nucleotide(['A','T','g'])=>36,nucleotide(['G','U','c'])=>50,nucleotide(['A','G','T'])=>45,
nucleotide(['U','G','A'])=>15,nucleotide(['G','C','u'])=>53,nucleotide(['A','t','C'])=>34,nucleotide(['c','G','A'])=>31,nucleotide(['g','G','U'])=>61,
nucleotide(['G','u','C'])=>50,nucleotide(['T','t','G'])=>4,nucleotide(['C','G','a'])=>31,nucleotide(['A','c','G'])=>40,nucleotide(['C','A','C'])=>26,
nucleotide(['G','U','U'])=>49,nucleotide(['C','C','T'])=>21,nucleotide(['C','a','C'])=>26,nucleotide(['C','a','A'])=>27,nucleotide(['A','A','u'])=>41,
nucleotide(['G','U','G'])=>52,nucleotide(['a','U','C'])=>34,nucleotide(['T','a','A'])=>11,nucleotide(['a','A','C'])=>42,nucleotide(['G','U','A'])=>51,
nucleotide(['g','G','A'])=>63,nucleotide(['G','A','g'])=>60,nucleotide(['C','u','G'])=>20,nucleotide(['A','g','T'])=>45,nucleotide(['U','c','G'])=>8,
nucleotide(['G','T','C'])=>50,nucleotide(['C','g','A'])=>31,nucleotide(['a','G','C'])=>46,nucleotide(['U','U','G'])=>4,nucleotide(['G','A','A'])=>59,
nucleotide(['U','U','U'])=>1,nucleotide(['U','C','C'])=>6,nucleotide(['T','a','T'])=>9,nucleotide(['u','A','A'])=>11,nucleotide(['G','g','G'])=>64,
nucleotide(['g','A','A'])=>59,nucleotide(['c','G','G'])=>32,nucleotide(['G','c','G'])=>56,nucleotide(['U','c','A'])=>7,nucleotide(['A','G','G'])=>48,
nucleotide(['G','g','A'])=>63,nucleotide(['G','C','A'])=>55,nucleotide(['a','A','U'])=>41,nucleotide(['u','G','A'])=>15,nucleotide(['c','U','A'])=>19,
nucleotide(['g','G','G'])=>64,nucleotide(['A','a','U'])=>41,nucleotide(['T','C','A'])=>7,nucleotide(['G','C','G'])=>56,nucleotide(['U','A','u'])=>9,
nucleotide(['u','G','U'])=>13,nucleotide(['c','A','G'])=>28,nucleotide(['C','A','u'])=>25,nucleotide(['c','U','C'])=>18,nucleotide(['G','u','G'])=>52,
nucleotide(['a','U','A'])=>35,nucleotide(['a','G','U'])=>45,nucleotide(['G','G','C'])=>62,nucleotide(['t','T','A'])=>3,nucleotide(['T','A','G'])=>12,
nucleotide(['A','U','C'])=>34,nucleotide(['C','a','T'])=>25,nucleotide(['A','U','a'])=>35,nucleotide(['a','A','G'])=>44,nucleotide(['G','A','c'])=>58,
nucleotide(['A','A','g'])=>44,nucleotide(['G','G','G'])=>64,nucleotide(['C','g','T'])=>29,nucleotide(['U','C','U'])=>5,nucleotide(['T','T','t'])=>1,
nucleotide(['t','G','A'])=>15,nucleotide(['T','a','C'])=>10,nucleotide(['T','A','g'])=>12,nucleotide(['U','A','C'])=>10,nucleotide(['c','C','T'])=>21,
nucleotide(['T','T','T'])=>1,nucleotide(['C','G','T'])=>29,nucleotide(['G','U','a'])=>51,nucleotide(['T','c','G'])=>8,nucleotide(['c','G','C'])=>30,
nucleotide(['U','u','C'])=>2,nucleotide(['U','A','G'])=>12,nucleotide(['A','g','G'])=>48,nucleotide(['T','T','C'])=>2,nucleotide(['g','G','T'])=>61,
nucleotide(['U','a','C'])=>10,nucleotide(['t','G','T'])=>13,nucleotide(['c','C','G'])=>24,nucleotide(['A','G','c'])=>46,nucleotide(['A','G','g'])=>48,
nucleotide(['U','A','c'])=>10,nucleotide(['u','C','C'])=>6,nucleotide(['A','A','t'])=>41,nucleotide(['T','G','c'])=>14,nucleotide(['U','u','U'])=>1,
nucleotide(['c','U','G'])=>20,nucleotide(['G','a','A'])=>59,nucleotide(['g','T','T'])=>49,nucleotide(['G','u','U'])=>49,nucleotide(['C','G','C'])=>30,
nucleotide(['C','U','G'])=>20,nucleotide(['g','G','C'])=>62,nucleotide(['C','T','T'])=>17,nucleotide(['A','A','c'])=>42,nucleotide(['G','u','A'])=>51,
nucleotide(['U','U','u'])=>1,nucleotide(['g','U','G'])=>52,nucleotide(['a','U','G'])=>36,nucleotide(['t','C','C'])=>6,nucleotide(['u','C','U'])=>5,
nucleotide(['C','A','U'])=>25,nucleotide(['g','A','U'])=>57,nucleotide(['A','g','U'])=>45,nucleotide(['T','C','C'])=>6,nucleotide(['T','g','C'])=>14,
nucleotide(['C','c','G'])=>24,nucleotide(['a','C','T'])=>37,nucleotide(['g','C','C'])=>54,nucleotide(['C','G','A'])=>31,nucleotide(['u','A','G'])=>12,
nucleotide(['U','U','A'])=>3,nucleotide(['A','G','A'])=>47,nucleotide(['t','G','G'])=>16,nucleotide(['U','G','g'])=>16,nucleotide(['G','a','U'])=>57,
nucleotide(['c','U','U'])=>17,nucleotide(['T','G','g'])=>16,nucleotide(['C','T','c'])=>18,nucleotide(['G','c','T'])=>53,nucleotide(['C','u','C'])=>18,
nucleotide(['G','T','c'])=>50,nucleotide(['T','A','t'])=>9,nucleotide(['C','c','U'])=>21,nucleotide(['C','U','C'])=>18,nucleotide(['U','G','U'])=>13,
nucleotide(['G','A','C'])=>58,nucleotide(['u','C','A'])=>7,nucleotide(['G','c','C'])=>54,nucleotide(['T','A','T'])=>9,nucleotide(['u','A','U'])=>9,
nucleotide(['U','a','G'])=>12,nucleotide(['T','g','A'])=>15,nucleotide(['c','T','A'])=>19,nucleotide(['A','u','U'])=>33,nucleotide(['c','A','T'])=>25,
nucleotide(['T','a','G'])=>12,nucleotide(['a','G','G'])=>48,nucleotide(['a','G','T'])=>45,nucleotide(['T','t','C'])=>2,nucleotide(['G','T','g'])=>52,
nucleotide(['t','A','A'])=>11,nucleotide(['G','t','C'])=>50,nucleotide(['T','g','G'])=>16,nucleotide(['T','T','c'])=>2,nucleotide(['A','c','T'])=>37,
nucleotide(['A','t','T'])=>33,nucleotide(['g','T','C'])=>50,nucleotide(['t','T','C'])=>2,nucleotide(['a','T','T'])=>33,nucleotide(['T','T','A'])=>3,
nucleotide(['T','c','C'])=>6,nucleotide(['T','T','a'])=>3,nucleotide(['T','G','G'])=>16,nucleotide(['T','G','T'])=>13,nucleotide(['A','G','U'])=>45,
nucleotide(['g','C','G'])=>56,nucleotide(['c','C','A'])=>23,nucleotide(['u','G','C'])=>14,nucleotide(['C','C','G'])=>24,nucleotide(['T','A','a'])=>11,
nucleotide(['g','A','T'])=>57,nucleotide(['T','G','A'])=>15,nucleotide(['g','A','G'])=>60,nucleotide(['A','c','C'])=>38,nucleotide(['C','u','A'])=>19,
nucleotide(['t','A','T'])=>9,nucleotide(['A','U','A'])=>35,nucleotide(['A','T','G'])=>36,nucleotide(['T','c','A'])=>7,nucleotide(['G','A','u'])=>57,
nucleotide(['A','U','G'])=>36,nucleotide(['G','c','U'])=>53,nucleotide(['G','G','T'])=>61,nucleotide(['A','G','a'])=>47,nucleotide(['T','C','t'])=>5,
nucleotide(['U','a','A'])=>11,nucleotide(['G','C','g'])=>56,nucleotide(['C','G','u'])=>29,nucleotide(['A','g','C'])=>46,nucleotide(['g','U','C'])=>50,
nucleotide(['U','A','U'])=>9,nucleotide(['c','A','A'])=>27,nucleotide(['A','A','G'])=>44,nucleotide(['C','G','g'])=>32,nucleotide(['C','a','G'])=>28,
nucleotide(['g','T','A'])=>51,nucleotide(['u','G','G'])=>16,nucleotide(['C','A','t'])=>25,nucleotide(['U','U','a'])=>3,nucleotide(['T','c','T'])=>5,
nucleotide(['U','g','U'])=>13,nucleotide(['t','T','G'])=>4,nucleotide(['C','T','t'])=>17,nucleotide(['A','A','U'])=>41,nucleotide(['a','A','T'])=>41,
nucleotide(['A','u','C'])=>34,nucleotide(['T','C','c'])=>6,nucleotide(['g','U','A'])=>51,nucleotide(['C','U','U'])=>17,nucleotide(['A','U','g'])=>36,
nucleotide(['C','U','a'])=>19,nucleotide(['A','a','T'])=>41,nucleotide(['U','c','C'])=>6,nucleotide(['C','A','g'])=>28,nucleotide(['T','C','a'])=>7,
nucleotide(['A','T','a'])=>35,nucleotide(['A','G','C'])=>46,nucleotide(['G','A','T'])=>57,nucleotide(['U','g','G'])=>16,nucleotide(['A','G','u'])=>45,
nucleotide(['U','C','u'])=>5,nucleotide(['A','C','t'])=>37,nucleotide(['U','g','C'])=>14,nucleotide(['G','T','t'])=>49,nucleotide(['C','A','c'])=>26,
nucleotide(['C','U','c'])=>18,nucleotide(['U','U','g'])=>4,nucleotide(['C','t','G'])=>20,nucleotide(['A','a','A'])=>43,nucleotide(['A','T','C'])=>34,
nucleotide(['C','A','A'])=>27,nucleotide(['U','A','a'])=>11,nucleotide(['t','C','T'])=>5,nucleotide(['a','G','A'])=>47,nucleotide(['C','c','A'])=>23,
nucleotide(['G','t','T'])=>49,nucleotide(['U','u','G'])=>4,nucleotide(['A','c','A'])=>39,nucleotide(['U','C','c'])=>6,nucleotide(['C','C','A'])=>23,
nucleotide(['u','U','A'])=>3,nucleotide(['G','U','C'])=>50,nucleotide(['T','A','A'])=>11,nucleotide(['U','U','c'])=>2,nucleotide(['C','T','a'])=>19,
nucleotide(['G','G','U'])=>61,nucleotide(['u','C','G'])=>8,nucleotide(['U','A','g'])=>12,nucleotide(['A','C','T'])=>37,nucleotide(['C','A','G'])=>28,
nucleotide(['C','G','U'])=>29,nucleotide(['G','g','C'])=>62,nucleotide(['U','g','A'])=>15,nucleotide(['G','a','T'])=>57,nucleotide(['G','g','U'])=>61,
nucleotide(['G','G','a'])=>63,nucleotide(['A','c','U'])=>37,nucleotide(['C','c','C'])=>22,nucleotide(['A','C','U'])=>37,nucleotide(['A','C','C'])=>38,
nucleotide(['U','u','A'])=>3,nucleotide(['G','G','A'])=>63,nucleotide(['T','C','T'])=>5,nucleotide(['U','A','A'])=>11,nucleotide(['C','C','c'])=>22,
nucleotide(['G','C','T'])=>53,nucleotide(['C','T','C'])=>18,nucleotide(['T','A','c'])=>10,nucleotide(['T','A','C'])=>10,nucleotide(['g','U','U'])=>49,
nucleotide(['T','G','t'])=>13,nucleotide(['G','t','G'])=>52,nucleotide(['C','T','A'])=>19,nucleotide(['A','t','A'])=>35,nucleotide(['C','U','g'])=>20,
nucleotide(['C','t','C'])=>18,nucleotide(['A','C','A'])=>39,nucleotide(['t','T','T'])=>1,nucleotide(['G','T','T'])=>49 )

immutable CodonTables
  ids::Vector{Int}
  amino::Matrix{AminoAcid}
  start::BitArray{2}
  codons::Dict{Vector{Nucleotide},Int}

  CodonTables(t::Vector{Int},a::Matrix{AminoAcid},s::BitArray{2},c::Dict{Vector{Nucleotide},Int}) = new(t,a,s,c)
end

const CODON_TABLES = CodonTables(_index,_amino,_start, _codons );

function show(io::IO,ct::CodonTables)
  println("Table IDs:")
  println([1:length(ct.ids)][ct.ids .!= 0])
  println(CODON_TABLES.amino)
  println("Codons IDs:")
  print(collect(CODON_TABLES.codons))
end

codon2aa(codon,tableid::Int) = CODON_TABLES.amino[ CODON_TABLES.codons[codon] , CODON_TABLES.ids[tableid] ]

function isstop(seq,tableid::Int)
  idx = CODON_TABLES.ids[tableid]
  len = length(seq)
  if len==3
    return( CODON_TABLES.amino[ CODON_TABLES.codons[ seq ] , idx ] == '*' )
  else
    ncod = Int(floor(len/3))
    out = falses(ncod)
    for i in 1:ncod
      out[i] = CODON_TABLES.amino[ CODON_TABLES.codons[ seq[(3i - 2):(3i)] ] , idx ] == '*'
    end
  return(out)
  end
end

isstop(tableid::Int) = CODON_TABLES.amino[ : , CODON_TABLES.ids[tableid] ] .== '*'

function isstart(seq,tableid::Int)
  idx = CODON_TABLES.ids[tableid]
  len = length(seq)
  if len==3
    return( CODON_TABLES.start[ CODON_TABLES.codons[ seq ] , idx ] )
  else
    ncod = Int(floor(len/3))
    out = falses(ncod)
    for i in 1:ncod
      out[i] = CODON_TABLES.start[ CODON_TABLES.codons[ seq[(3i - 2):(3i)] ] , idx ]
    end
  return(out)
  end
end

isstart(tableid::Int) = CODON_TABLES.start[ : , CODON_TABLES.ids[tableid] ]

function translateCDS(seq,tableid::Int)
  len = length(seq)
  idx = CODON_TABLES.ids[tableid]
  if len%3 == 0
    if CODON_TABLES.start[ CODON_TABLES.codons[  seq[1:3] ] , idx ]
      if CODON_TABLES.amino[ CODON_TABLES.codons[ seq[len-2:len] ] , idx ] == '*'
	out = Array( AminoAcid, div(len,3)-1 )
	ind = 1
	out[ ind ] = 'M'
	for i in 4:3:len-5
	  ind += 1
	  out[ ind ] = CODON_TABLES.amino[ CODON_TABLES.codons[ seq[i:i+2] ] , idx ]
	  if out[ ind ] == '*'
	    throw("$(seq[i:i+2]) at $i is a stop codon on Table $tableid")
	  end
	end
	return(out)
      else
	throw("$(seq[len-2:len]) is not a valid stop codon on Table $tableid")
      end
    else
      throw("$(seq[1:3]) is not a valid start codon on Table $tableid")
    end
  else
    throw("CDS length must be a multiple of three")
  end
end

function translatetostop(seq,tableid::Int)
  idx = CODON_TABLES.ids[tableid]
  out = AminoAcid[]
  for i in 1:3:length(seq)-2
    aa = CODON_TABLES.amino[ CODON_TABLES.codons[ seq[i:i+2] ] , idx ]
    if aa != '*'
      push!( out , aa )
    else
      return(out)
    end
  end
  return(out)
end

function translate(seq,tableid::Int)
  idx = CODON_TABLES.ids[tableid]
  out = AminoAcid[]
  for i in 1:3:length(seq)-2
    push!( out , CODON_TABLES.amino[ CODON_TABLES.codons[ seq[i:i+2] ] , idx ] )
  end
  return(out)
end


text:
export toOdeFile, fromOdeFile, parseOutputFile

newline = os[OS_NAME].newline

#-------------------------------------------------------------------------------
#
#                       Parsing XPPjl-model to ODE-file
#
#-------------------------------------------------------------------------------

"""
Function for writing a Model-instance to an .ode file

    toOdeFile(ModelInstance)

...will write the model instance to an ode file with the name ModelInstance.name
"""
function toOdeFile(M::Model; path = "")
    file = "#" * M.name * newline * "#generated using XPPjl" * newline * newline * "#ODEs:$newline"
    for r in M.odes
        file *= r[1] * "\'=" * r[2] * newline
    end
    #Algebraic equation automatically also auxvars
    file *= newline * newline * "#Algebraic and auxilliary equations:" * newline
    for a in M.alg
        file *=  a[1] * "=" * a[2] * newline
    end
    for x in M.aux
        file *= "aux " * x[1] * "=" * x[2] * newline
    end
    #Parameters
    file *= newline * newline * "#Parameters:\n"
    for p in sort(collect(M.pars))
        file *= "p " * p[1] * "=" * string(p[2]) *newline
    end
    #Initials
    file *= newline * newline * "#Initials:\n"
    for i in M.init
        file *= "init " * i[1] * "=" * string(i[2]) *newline
    end
    #Settings
    file *= newline * newline * "#Settings:\n"
    for s in M.spec
        file *= "@ " * s[1] * "=" * string(s[2]) *newline
    end
    file *= "done" * newline
    f = open(joinpath(path, M.name) * ".ode", "w")
    write(f, file)
    close(f)
end

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

"""
High-level routine for generating a model instance from an existing .ode-file

    ModelName = fromOdeFile(\"odeFilenameAsString.ode\")

Syntax requirements:

    - odes are specified with ' instead of dVar/dt

    - parameters are specified with the 'p '-prefix

    - settings are specified with the '@ '-prefix

    - initial conditions are specified with the 'init '-prefix

    - every line can only contain a single specification, i.e. no concatentation of multiple parameter/initial/setting specifiactions in a single line

"""
function fromOdeFile(filename::AbstractString)
     f = open(filename)
     name = split(filename, ".ode")[1] * "_xppjl"
     M = parseOdeFile(f, name)
     M.originalState = Dict(:pars => deepcopy(M.pars), :init => deepcopy(M.init), :spec => deepcopy(M.spec))
     close(f)
     return(M)
end

#-------------------------------------------------------------------------------
#
#                       Parsing ODE-file to XPPjl
#
#-------------------------------------------------------------------------------

#                       Parsing subfunctions
#-------------------------------------------------------------------------------

type ParsedLine
    b::Bool
    name::AbstractString
    value::Any
end
function ignore(l)
    b = length(l) < 3 || l[1] == '#' || contains(l, "done")
    name = ""
    value = ""
    return(ParsedLine(b, name, value))
end

function parameter(l)
    b = l[1:2] == "p "
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = parts[1][3:end]
        value = float(parts[2])
    end
    return(ParsedLine(b, name, value))
end

function initial(l)
    b = l[1:5] == "init "
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = parts[1][6:end]
        value = float(parts[2])
    end
    return(ParsedLine(b, name, value))
end

function method(l)
    b = l[1:2] == "@ "
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = parts[1][3:end]
        value = parts[2]
    end
    return(ParsedLine(b, name, value))
end

function variable(l)
    b = contains(split(l, "=")[1], "\'")
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = parts[1][1:end-1]
        value = parts[2]
    end
    return(ParsedLine(b, name, value))
end

function auxilliary(l)
    b = l[1:4] == "aux "
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = split(parts[1], "aux ")[2]
        value = parts[2]
    end
    return(ParsedLine(b, name, value))
end

function algebraic(l)
    b = length(split(l, "=")) >= 2
    name = ""
    value = ""
    if b
        parts = split(l, "=")
        name = parts[1]
        value = split(l, name*"=")[2]
    end
    return(ParsedLine(b, name, value))
end

#                       High-level routine
#-------------------------------------------------------------------------------

"""
Function for parsing .ode files that obey the following rules:

    - odes are specified with \' instead of dVar/dt

    - parameters are specified with the 'p '-prefix

    - settings are specified with the '@ '-prefix

    - initial conditions are specified with the 'init '-prefix

    - every line can only contain a single specification, i.e. no concatentation of multiple parameter/initial/setting specifiactions in a single line

"""
function parseOdeFile(f::IOStream, modelname::AbstractString)
    odes = Dict()
    init = Dict()
    pars = Dict()
    aux = Dict()
    alg =Dict()
    spec = Dict()
    vars = Any[]
    for l in eachline(f)
        l = split(string(l), newline)[1];
        l = split(l, "\r")[1]
        if ignore(l).b
            #comment or empty line: do nothing
            #both auxilliary and algebraic equation treated as the same
        elseif parameter(l).b
            pars[parameter(l).name] = parameter(l).value
        elseif initial(l).b
            init[initial(l).name] = initial(l).value
        elseif method(l).b
            spec[method(l).name] = method(l).value
        elseif variable(l).b
            push!(vars, variable(l).name)
            odes[variable(l).name] = variable(l).value
        elseif auxilliary(l).b
            push!(vars, auxilliary(l).name)
            aux[auxilliary(l).name] = auxilliary(l).value
        elseif algebraic(l).b
            alg[algebraic(l).name] = algebraic(l).value
        end
    end
    M = Model(odes, init, pars; name = modelname, aux = aux, alg = alg, spec = spec, vars = vars)
    return(M)
end

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

#-------------------------------------------------------------------------------
#
#                       Parsing Output-file to XPPjl
#
#-------------------------------------------------------------------------------

"""
Parse output file and store it as new SimulationData-instance in the Model.sims-dict

    SimulationDataInstance = parseOutputFile(file-IOstream, ModelInstance)
"""
function parseOutputFile(f::IOStream, M::Model, name = false)
    if name == false
        #Get the new key for the dict
        k = length(M.sims) + 1
    else
        #Overwrite the last simulation
        k = name
    end
    #Instantiate new SimulationData-structure
    M.sims[k] = SimulationData(M, name)
    #loop over lines in file
    for l in eachline(f)
        #Remove newline and space at the end of each line
        l = split(string(l), " \n")[1];
        #Split line to get data points
        pts = split(l, " ");
        #Loop over var-list & append to SimulationData.D-instance based on index
        for v in M.vars
            i = findfirst(M.vars, v)
            push!(M.sims[k].D[v], float(pts[i]))
        end
    end
    return(M)
end

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

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
module TestLDict

using FactCheck
using DataCubes

facts("LDict tests") do
  context("constructor tests") do
    d1 = LDict(Dict(:a=>1, :c=>3, :b=>2), [:c, :b, :a])
    d2 = LDict([:c, :b, :a], [3, 2, 1])
    d3 = LDict(:c=>3, :b=>2, :a=>1)
    d4 = LDict((:c,3), (:b,2), (:a,1))
    @fact d1 --> d2
    @fact d2 --> d3
    @fact d3 --> d4
    @fact_throws LDict([:a,:b,:c], [1,2])
    @fact DataCubes.create_ldict_nocheck() --> isempty
    @fact DataCubes.create_ldict_nocheck((:a,3),(:b,2)) --> LDict(:a=>3,:b=>2)
    @fact DataCubes.create_ldict_nocheck(Dict(:a=>3,:b=>2)) --> Dict(:a=>3,:b=>2)
    @fact LDict(Dict(:a=>3,:b=>2)) --> Dict(:a=>3,:b=>2)
    @fact eltype(LDict(:a=>3,:b=>2)) --> Pair{Symbol,Int}
    @fact keys(LDict(:a=>3,:b=>2)) --> [:a,:b]
    @fact values(LDict(:a=>3,:b=>2)) --> [3,2]
    @fact_throws LDict(:a=>3,:b=>'a')[:x]
    d5 = LDict(:a=>3, :b=>'a')
    @fact get!(d5, :a, 1.0) --> 3
    @fact get!(d5, :c, 1.0) --> 1.0
  end
  context("method tests") do
    @fact merge(LDict(:a=>1, :b=>2, :c=>3), LDict(:a=>10, :d=>15)) -->
      LDict(:a=>10, :b=>2, :c=>3, :d=>15)
    @fact merge(LDict(:a=>1, :b=>2, :c=>3), LDict(:a=>10, :d=>"15")) -->
      LDict(:a=>10, :b=>2, :c=>3, :d=>"15")
    @fact delete(LDict(:a=>[10,11,12], 10=>3, :c=>5.0, :d=>1), :a, 10) -->
      LDict(:c=>5.0, :d=>1)
    @fact pick(LDict(:a=>[10,11,12], 10=>3, :c=>5.0, :d=>1), [:d, :c]) -->
      LDict(:d=>1, :c=>5.0)
  end
end

end


text:
using Documenter, SomervilleCouncilParser

DocMeta.setdocmeta!(SomervilleCouncilParser, :DocTestSetup,
                    :(using SomervilleCouncilParser); recursive=true)

makedocs(; format=Documenter.HTML(; prettyurls=true), modules=[SomervilleCouncilParser],
         sitename="SomervilleCouncilParser.jl",
         pages=["Home" => "index.md", "Getting Started" => "new_users.md",
                "Examples" => "examples.md", "API reference" => "api.md"], strict=true)

deploydocs(; repo="github.com/hannahilea/SomervilleCouncilParser.jl.git", push_preview=true,
           devbranch="main", devurl="dev")


text:
using Documenter, POMDPStressTesting, PyPlot, Seaborn

# TODO: Literate.jl to import notebook

makedocs(
    modules = [POMDPStressTesting, AST, BlackBox, GrayBox],
    sitename = "POMDPStressTesting.jl",
    authors = "Robert Moss",
    clean = false,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Installation" => "install.md",
            "Guide" => "guide.md",
            "Solvers" => "solvers.md",
        ],
        "Example" => "example.md",
        "Library/Interface" => "interface.md",
        "Contributing" => "contrib.md"
    ]
)

deploydocs(
    repo = "github.com/sisl/POMDPStressTesting.jl.git",
)


text:
"A type for the effective wave ansatz."
mutable struct EffectiveWave{T<:AbstractFloat}
    hankel_order::Int
    amplitudes::Array{Complex{T}} # the effective ampliudes
    k_eff::Complex{T} # the effective wavenumber
    θ_eff::Complex{T} # the effective transmission angle
end

zero(W::Type{EffectiveWave{T}}) where {T<:AbstractFloat} = EffectiveWave(0,[zero(Complex{T})],zero(Complex{T}),zero(Complex{T}))

effective_waves(ω::T, medium::Medium{T}, specie::Specie{T}; kws...) where T<:AbstractFloat =  effective_waves(ω, medium, [specie]; kws...)

"Calculates the effective wavenumbers and return Vector{EffectiveWave}."
function effective_waves(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; tol::T = 1e-6,
    extinction_rescale::Bool = false, kws...) where T<:AbstractFloat
    # as there will be likely more than 1 k_eff we set extinction to false.

    k_effs = wavenumbers(ω, medium, species; tol = tol, kws... )
    wave_effs = [
        EffectiveWave(ω, k_eff, medium, species; tol = tol, extinction_rescale=extinction_rescale, kws...)
    for k_eff in k_effs]

    return wave_effs
end

function EffectiveWave(ω::T, k_eff::Complex{T}, medium::Medium{T}, species::Vector{Specie{T}};
        θin::T = 0.0, tol::T = 1e-7,
        hankel_order::Int = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier::T = 1.005,
        method::Symbol = :none,
        extinction_rescale::Bool = true,
        kws...
    ) where T<:AbstractFloat

    k = ω/medium.c
    θ_eff = transmission_angle(k, k_eff, θin; tol = tol)
    if method == :WienerHopf
        if hankel_order > 0
            error("the Wiener Hopf method has only been implemented for monopole scatterers, i.e. hankel order = 0. ")
        end
        amps = wienerhopf_wavevectors(ω, k_eff, medium, species;
            hankel_order=hankel_order, tol = tol,
            radius_multiplier=radius_multiplier,
            θin = θin, kws...
        )
    else
        amps = effective_wavevectors(ω, k_eff, medium, species;
            hankel_order=hankel_order, tol = tol,
            radius_multiplier=radius_multiplier
        )
    end
    wave_eff = EffectiveWave(hankel_order, amps, k_eff, θ_eff)
    if extinction_rescale && method != :WienerHopf
        amps = amps.*scale_amplitudes_effective(ω, wave_eff, medium, species; tol = tol, θin=θin)
    end

    return EffectiveWave(hankel_order, amps, k_eff, θ_eff)
end