text:
using ParticleDetection
using Test
using Serialization

@testset "Detect Particles" begin
    ima=deserialize(joinpath(@__DIR__, "particles.bin"))
    @time P=detect_particles(ima,6,10,0.2f0)
    @test find_best_radius(ima,P,10) == 12
    @test P==[(1081, 32), (1112, 92), (1145, 120), (1101, 149), (818, 308), (2205, 322),
             (786, 412), (839, 444), (864, 449), (774, 464), (830, 487), (774, 503),
              (871, 510), (801, 523), (816, 537), (835, 541), (805, 549), (863, 555),
              (848, 558), (895, 570), (916, 574), (850, 576), (863, 609), (997, 674),
              (973, 682), (921, 685), (986, 686), (931, 698), (996, 710), (1403, 714),
              (1417, 723), (1160, 740), (1283, 762), (1288, 778), (909, 784), (932, 802),
              (916, 805), (895, 809), (944, 847), (963, 854), (344, 866), (345, 873),
              (917, 873), (930, 951), (701, 971), (876, 1032), (1524, 1099), (1765, 1226),
              (1778, 1238), (1827, 1278), (944, 1280), (933, 1292), (961, 1294), (1285, 1318),
              (1688, 1329), (1285, 1333), (1748, 1333), (1219, 1341), (1285, 1347), (108, 1423),
              (53, 1464), (1276, 1552), (623, 1578), (1194, 1580), (1177, 1583), (1183, 1609),
              (1235, 1610), (1264, 1650), (2193, 1655), (1924, 1674), (1312, 1684), (1914, 1686),
              (2287, 1690), (1966, 1710), (2225, 1710), (1147, 1715), (2253, 1718), (1113, 1719),
              (1327, 1722), (1396, 1738), (1773, 1749), (2154, 1762), (1201, 1770), (2054, 1775),
              (1892, 1779), (1957, 1782), (2023, 1783), (1126, 1794), (1112, 1802), (2086, 1807),
              (1090, 1833), (1573, 1846), (1759, 1849), (2086, 1850), (1127, 1851), (1457, 1859),
              (1167, 1864), (1510, 1885), (887, 1892), (557, 1914), (1214, 1918), (2066, 1939),
              (1848, 1943), (1914, 1953), (1899, 1955), (905, 1957), (1851, 1957), (1866, 1957),
              (2269, 1958), (998, 1965), (1005, 2000), (1373, 2017), (1029, 2043), (2225, 2059),
              (887, 2075), (1632, 2146), (1249, 2147), (1734, 2152), (879, 2153), (836, 2169),
              (892, 2173), (95, 2302), (162, 2327), (431, 2346), (45, 2364), (161, 2364), (268, 2364)]
end
#
# ima=deserialize(joinpath(@__DIR__, "particles.bin"))
# @time P=detect_particles(ima,6,10,0.2f0)
# b_rad=find_best_radius(ima,P,10)
# @time P2=detect_particles(ima,6,b_rad,0.2f0)
#
#
# using PyPlot
# pygui(true)
#
# imshow(ima)
# scatter(getindex.(P2,2).-1,getindex.(P2,1).-1,color="red")
# scatter(getindex.(P,2).-1,getindex.(P,1).-1,color="blue",alpha=0.5)
#
#


text:
module GLUE
using Pkg.Artifacts

function download()
    artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
    glue_hash = artifact_hash("glue", artifact_toml)
    if glue_hash === nothing || !artifact_exists(glue_hash)
        glue_hash = create_artifact() do artifact_dir
            download_glue_data = tempdir()*"/download_glue_data.py"
            Base.download("https://raw.githubusercontent.com/nyu-mll/jiant/master/scripts/download_glue_data.py", download_glue_data)
            run(`python $download_glue_data --data_dir $artifact_dir`)
            rm(download_glue_data)
        end
        bind_artifact!(artifact_toml, "glue", glue_hash)
    end
    return artifact_path(glue_hash)
end

end

ext: jl


text:
using DimensionalityReduction
using Gadfly

# TODO: approximate (generalized) principal components rather than principal components could be more efficient for N > T
type FactorModel
    # models X_t = lambda*factors_t + e_t  # Here lambda is not equal to lambda(L) --> static factors
    x::Array{Float64, 2}  # variables to calculate factors from
    number_of_factors::Int64  # columns of factors we use (which e.g. capture a certain percentage of the variation or from Bai Ng info criteria)
    number_of_factors_criterion::String  # Bai Ng info criterion to use
    factor_type::String  # "principal components" later maybe also "maximum likelihood"
    factors::Array{Float64, 2}  # matrix of static factors
    loadings::Array{Float64, 2}
    factor_residuals::Array{Float64, 2}  # residuals from the factor estimation  x = factors * lambda

    function FactorModel(x::Matrix{Float64}, number_of_factors::Int64=int(ceil(minimum(size(x))/2)), number_of_factors_criterion::String="", factor_type::String="principal components")
        # (static) factor equation:
        factors, loadings, number_of_factors = calculate_factors(x, factor_type, number_of_factors)
        # TODO: check if correct factors are used (p 1136 on Bai Ng 2006)
        factor_residuals = x .- factors[:, 1:number_of_factors] * loadings[:, 1:number_of_factors]'
        #factor_residuals = x - x*loadings[:, 1:number_of_factors]*inv(loadings[:, 1:number_of_factors]'loadings[:, 1:number_of_factors])*loadings[:, 1:number_of_factors]' # rescaled version from french dude
        # TODO: for lags we have to regress residuals on its lags and use info criterio to choose lag length p. this should be offloaded into another method
        # TODO: solve chicken and egg problem between number_of_factors and number_of_factor_lags. Both require residuals. --> maybe select one after the other repeatedly until convergence?
        # TODO: do the above e.g. factors, loadings, number_of_factors = augment_factors_by_lags(factors, number_of_factor_lags)  # TODO: do we have to

        return(new(x, number_of_factors, number_of_factors_criterion, factor_type, factors, loadings, factor_residuals))
    end

    function FactorModel(x::Matrix{Float64}, number_of_factors_criterion::String, factor_type::String="principal components"; max_factors = int(ceil(minimum(size(x))/2)))
        #println("finding optimal number of factors with maximum number of factors=", max_factors, " and information criterion=", number_of_factors_criterion)
        # number of factors to include is not given but a criterion --> we use
        # the criterion to determine the number of factors using the Bai, Ng
        # information criteria.
        models = [apply(FactorModel, (x, number_of_factors, number_of_factors_criterion, factor_type)) for number_of_factors in 1:max_factors]
        # we keep all the models in memory which can be a problem depending on the dimensions of x. TODO: will refactor later when debugging and testing is done
        criteria = [calculate_criterion(model) for model in models]
        return models[indmin(criteria[1:max_factors])]  # keep the model with the best information criterion
    end

end

function principal_components(x)  # calculate factors and loadings for a given x matrix
    T, N = size(x)  # T: time dimension, N cross-sectional dimension
    # see Stock, Watson (1998)
    if T >= N
        eigen_values, eigen_vectors = eig(x'x)  # x'x is NxN
        eigen_values, eigen_vectors = reverse(eigen_values), eigen_vectors[:, size(eigen_vectors, 2):-1:1]  # reverse the order from highest to lowest eigenvalue
        loadings = sqrt(N) * eigen_vectors  # we may rescale the eigenvectors say Bai, Ng 2002
        factors = x*loadings/N  # this is from Bai, Ng 2002 except that for me the inverse of the loadings matrix primed is not the same as the loadings matrix
        #factors = factors*(factors'factors/T)^(1/2)  # TODO: not sure it is correct to rescale like this (see Bai, Ng 2002 p.198)
        #factors = 1/N*x*loadings'  # this is from Stock, Watson 2010
        #loadings = (x'x)*loadings/(N*T)  # and this is from C. Hurlin from University of Orléans... rescales factor loadings to have the same variance as x
    end
    if N > T
        eigen_values, eigen_vectors = eig(x*x')  # x*x' is TxT
        eigen_values, eigen_vectors = reverse(eigen_values), eigen_vectors[:, size(eigen_vectors, 2):-1:1]  # reverse the order from highest to lowest eigenvalue
        factors = sqrt(T) * eigen_vectors  # sqrt(T) comes from normalization F'F/T = eye(r) where r is number of factors (see Bai 2003), factors are Txr
        loadings = x'factors/T  # see e.g. Bai 2003 p.6 or Breitung, Eickmeier 2011 p. 80. Dimension of loadings is Nxr where r here is also N
        # betahat = loadings * chol(loadings'loadings/N)  # C. Hurlin from University of Orléans rescales like this
    end
        # resultingly factors*loadings' estimates x
    return factors, loadings, eigen_values, eigen_vectors
end

# calculate the factors and the rotation matrix to transform data into the space spanned by the factors
function calculate_factors(x::Matrix, factor_type::String="principal components", number_of_factors=ceil(minimum(size(x))/2))
    if factor_type == "principal components"
        factors, loadings = principal_components(x)
        max_factor_number = int(ceil(minimum(size(x))/2))
    elseif factor_type == "squared principal components"  # include squares of X
        factors, loadings = principal_components([x x.^2])
        max_factor_number = minimum([size(x, 1), size(x,2)*2])
    elseif factor_type == "quadratic principal components"  # include squares of X and interaction terms - better only use in combination with targeted_predictors
        pca_cols = x  # columns to use for principal components
        for i in 1:size(x, 2)
            for j in 1:size(x, 2)
                pca_cols = hcat(pca_cols, x[:, i].*x[:, j])  # not very efficient!
            end
        end
        max_factor_number = int(ceil(minimum(size(pca_cols))/2))
        factors, loadings = principal_components(pca_cols)
    end
    if number_of_factors > max_factor_number
        number_of_factors = max_factor_number
        warn("can not estimate more than `minimum(size(x))` factors with $factor_type. Number of factors set to $number_of_factors")
    end
    return factors, loadings, number_of_factors
end

function scree_plot(x::Array{Float64, 2}, max_factors::Int64=size(x, 2); file_name::String="")
    eigen_values = principal_components(x)[3][1:max_factors]
    plt = plot(
        layer(x=[1:length(eigen_values)], y=eigen_values, Geom.point, Theme(default_color=color("green"))),
        layer(x=[1:length(eigen_values)], y=eigen_values, Geom.line),
        Guide.XLabel("Number of Factors"), Guide.YLabel("Eigenvalue"),
        Scale.y_continuous(format=:plain),
        Scale.x_discrete
    )
    if length(file_name) > 0
        draw(PNG(file_name, 20cm, 15cm), plt)
    end
    plt
end

function Base.show(io::IO, fm::FactorModel)
    @printf io "Static Factor Model\n"
    @printf io "Dimensions of X: %s\n" size(fm.x)
    @printf io "Number of factors used: %s\n" fm.number_of_factors
    @printf io "Factors calculated by: %s\n" fm.factor_type
    @printf io "Factor model residual variance: %s\n" sum(fm.factor_residuals.^2)/apply(*, size(fm.x))
end

function calculate_criterion(fm::FactorModel)
    if fm.number_of_factors_criterion != ""
        number_of_factors_criterion = fm.number_of_factors_criterion
        number_of_factors_criterion_value::Float64 = eval(symbol("criterion_$number_of_factors_criterion"))(fm)
    end
    return(number_of_factors_criterion_value)
end

function predict(fm::FactorModel, y::Array{Float64, 1}, h::Int64=1, number_of_lags::Int64=5, number_of_factors::Int64=0)
    if number_of_factors == 0  # number of factors can be given to set them to a different value in the forecasting equation than in the factor equation
        #println("using the same number of factors for forecasting as in the factor equation: ", fm.number_of_factors)
        number_of_factors = fm.number_of_factors
    end
    # makes a h step ahead forecast of y (which can also be in the factor model)
    w = lagged_matrix(y, [0, [h:number_of_lags+h-1]])
    # estimate y_{t+h} = alpha'w_t + Gamma'x_t

    y, x = w[1:end-1, 1], hcat(w[:, 2:end], fm.factors[end-size(w,1)+1:end, 1:number_of_factors])  # last observation of y is reserved for prediction
    x = hcat(ones(size(x, 1)), x)  # add a constant term to the regression
    new_x, design_matrix = x[end,:], x[1:end-1, :]  # we reserve the last observation for prediction and dont use it for learning
    coefficients = inv(design_matrix'design_matrix)*design_matrix'y
    prediction = new_x*coefficients
    residuals = y - design_matrix*coefficients
    #hat_matrix = design_matrix*inv(design_matrix'design_matrix)*design_matrix'
    #residual_variance = (residuals.^2)./(1.-diag(hat_matrix))  # HC 2
    #coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'diagm(residual_variance)design_matrix)*inv(design_matrix'design_matrix)

    return(residuals, prediction[1])
end


text:
function solve_tsp(dist_mtx::Matrix{Int})
    if !issymmetric(dist_mtx)
        error("Asymmetric TSP is not supported.")
    end

    n_nodes = size(dist_mtx, 1)

    lower_diag_row = Int[]
    for i in 1:n_nodes
        for j in 1:i
            push!(lower_diag_row, dist_mtx[i, j])
        end
    end
    buf = 10
    n_rows = length(lower_diag_row) / buf |> ceil |> Int
    rows = String[]
    for i in 1:n_rows
        s = buf * (i-1) + 1
        t = min(buf * i, length(lower_diag_row))
        push!(rows, join(lower_diag_row[s:t], " "))
    end

    name = randstring(10)
    filename = name * ".tsp"
    open(filename, "w") do io
        write(io, "NAME: $(name)\n")
        write(io, "TYPE: TSP\n")
        write(io, "COMMENT: $(name)\n")
        write(io, "DIMENSION: $(n_nodes)\n")
        write(io, "EDGE_WEIGHT_TYPE: EXPLICIT\n")
        write(io, "EDGE_WEIGHT_FORMAT: LOWER_DIAG_ROW \n")
        write(io, "EDGE_WEIGHT_SECTION\n")
        for r in rows
            write(io, "$r\n")
        end
        write(io, "EOF\n")
    end

    return __solve_tsp__(filename)
end

function solve_tsp(x::Vector{Float64}, y::Vector{Float64}; dist="EUC_2D")
    n_nodes = length(x)
    @assert length(x) == length(y)

    name = randstring(10)
    filename = name * ".tsp"

    open(filename, "w") do io
        write(io, "NAME: $(name)\n")
        write(io, "TYPE: TSP\n")
        write(io, "COMMENT: $(name)\n")
        write(io, "DIMENSION: $(n_nodes)\n")
        write(io, "EDGE_WEIGHT_TYPE: $(dist)\n")
        write(io, "EDGE_WEIGHT_FORMAT: FUNCTION \n")
        write(io, "NODE_COORD_TYPE: TWOD_COORDS \n")
        write(io, "NODE_COORD_SECTION\n")
        for i in 1:n_nodes
            write(io, "$i $(x[i]) $(y[i])\n")
        end
        write(io, "EOF\n")
    end

    return __solve_tsp__(filename)
end

function solve_tsp(org_tsp_file::String)
    if !isfile(org_tsp_file)
        error("$(org_tsp_file) is not a file.")
    end

    name = randstring(10)
    filename = name * ".tsp"
    cp(org_tsp_file, filename)

    return __solve_tsp__(filename)
end

function __solve_tsp__(tsp_file::String)
    status = run(`$(Concorde.CONCORDE_EXECUTABLE) $(tsp_file)`, wait=false)
    while !success(status)
        #
    end

    name = split(basename(tsp_file), ".")[1]
    sol_filepath =  name * ".sol"
    opt_tour = read_solution(sol_filepath)

    tsp = readTSP(tsp_file)
    M = Int.(tsp.weights)
    opt_len = tour_length(opt_tour, M)

    cleanup(name)

    return opt_tour, opt_len
end


text:
#__precompile__()

module ExperimentalDesign

using Random
using Logging
using Primes
using DataFrames
using Distributions
using StatsModels
using DocStringExtensions
using LinearAlgebra

import Random: rand

# Types

export    AbstractDesign,   AbstractScreeningDesign,    AbstractFactorialDesign,
    AbstractOptimalDesign,  PlackettBurman, FullFactorial,  FractionalFactorial,
    OptimalDesign, RandomDesign

# Methods

export fullfactorial,  explicit_fullfactorial, plackettburman, isplackettburman,
    paley,     optimize_design,     expanded_design,     generate_random_design,
    rand, random_design, d_criterion, next_offset_divisible_prime

include("design.jl")
include("plackett_burman.jl")
include("factorial.jl")
include("random.jl")
include("d_optimal/variance_predictions.jl")

end


text:
abstract type Basis end

struct Oscillator <: Basis
    n::Int64
    ω::Float64
end

Oscillator(n::Int64) = Oscillator(n,1.0)
index(b::Oscillator) = 1:b.n
spectrum(b::Oscillator) = (0:b.n-1)*b.ω
qnumbers(b::Oscillator) = 0:b.n-1

function hermite(x,n)
    @assert n >= 1
    h1 = exp(-x^2 / 2) * (1 / π)^(1 / 4)
    n == 1 && (return h1)
    h2 = sqrt(2) * x * h1
    n == 2 && (return h2)
    h3 = 0.0
    for j ∈ 3:n
        h3 = sqrt(2/(j - 1)) * x * h2 - sqrt((j-2)/(j-1)) * h1
        h1 = h2; h2 = h3
    end
    n > 2 && (return h3)
end
hermite(x,n,ω) = hermite(√ω * x,n) * ω ^ ( 1 / 4 )

function (b::Oscillator)(x::Array{Float64,1})
    @unpack ω,n = b
    T = x * ones(1, n) |> zero
    T[:,1] = @. hermite(x,1,ω)
    n >= 2 && (T[:,2] = @. sqrt(2 * ω) * x * T[:,1])
    for j ∈ index(b)[3:end]
        m = qnumbers(b)[j]
        T[:,j] = @. sqrt(2 / m) * (√ω * x) * T[:,j-1] - sqrt((m-1) / m) * T[:,j-2]
    end
    return T
end

(b::Oscillator)(x::LinRange{Float64}) = b(x |> collect)
(b::Oscillator)(x::Float64,n::Int64) = hermite(x,n,b.ω)

function hermite_polar(x,n,ω)
    p = Progress(n,.1,"Precomputing...",20)   # minimum update interval: 1 second
    T = zeros(size(x)...,n)
    T[:,:,1] = @. hermite(x,1,ω)
    next!(p)
    n >= 2 && (T[:,:,2] = @. sqrt(2 * ω) * x * T[:,:,1])
    next!(p)
    for j in 3:n
        T[:,:,j] = @. sqrt(2 / (j-1)) * (√ω * x) * T[:,:,j-1] - sqrt((j-2) / (j-1)) * T[:,:,j-2]
    next!(p)
    end
    return T
end
hermite_polar(x,n) = hermite_polar(x,n,1.)

# filter cartesian data
function filterH(psi,x,N,ω)
    p = Progress(N,.1,"Filtering...",20)
    H = Oscillator(N,ω)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T
    psiF = zero(psi)
        for m = 1:N
            next!(p)
            for n = 1:(N+1-m)
            @. psiF += H[:,m]*F̂[m,n]*H[:,n]'
            end
        end
    return psiF
end

# covert to polar coords
function slowpolar(psi,x,N,ω)
    H = Oscillator(N,ω)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T

    Nx = length(x)
    θ = LinRange(0,2*pi,2*Nx)
    r = LinRange(0,last(x),Nx/2 |> Int)'

    psiP = zero(θ*r)
    for m = 1:N, n = 1:(N + 1 - m)
        @. psiP += hermite(r*cos(θ),m) * F̂[m,n] * hermite(r*sin(θ),n)
    end
    return psiP
end

function init_polar(x,N,ω)
    Nx = length(x)
    θ = LinRange(0,2*pi,2*Nx)
    r = LinRange(0,last(x),Nx/2 |> Int)'
    hx = hermite_polar((@. r*cos(θ)),N,ω)
    hy = hermite_polar((@. r*sin(θ)),N,ω)
    return hx,hy
end

function polar(psi,x,ω,Hx,Hy)
    N = size(Hx)[3]
    p = Progress(N,1,"Polar transform...",20)
    H = Oscillator(N,ω)(x)
    T = inv(H'*H)
    F̂ = T*H'*psi*H*T

    Nx = length(x)
    Nθ = 2*Nx
    Nr = Nx/2 |> Int
    θ = LinRange(0,2*pi,Nθ+1)[1:Nθ]
    r = LinRange(0,last(x),Nr)'

    psiFP = zeros(eltype(psi),Nθ,Nr)
    for m = 1:N
        next!(p)
        for n = 1:(N + 1 - m)
        hx = @view Hx[:,:,m]
        hy = @view Hy[:,:,n]
        @. psiFP += hx * F̂[m,n] * hy
    end
    end
    return psiFP
end


text:
using Distributions
using JSON

r = 3
p = 0.25
d = NegativeBinomial( r, p )

x = [ -0.5, -0.25, 0, 0.25, 0.5 ]

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("r", r),
	("p", p),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/number.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

max_stars_repo_name: distributions-io/negative-binomial-mgf

max_issues_repo_name: distributions-io/negative-binomial-mgf

max_forks_repo_name: distributions-io/negative-binomial-mgf


text:
export get_response
export plot_response

"Get the response matrix corresponding to `res`.
Any substitution rules not specified in `res` can be supplied in `rules`."
function ResponseMatrix(res::Result; rules=Dict())

    # get the symbolic response matrix
    @variables Δ
    M = get_response_matrix(res.problem.eom.natural_equation, Δ, order=2)
    M = substitute_all(M, merge(res.fixed_parameters, rules))
    symbols = cat(res.problem.variables, collect(keys(res.swept_parameters)), dims=1)
    compiled_M = [build_function(el, cat(symbols, [Δ], dims=1)) for el in M]

    ResponseMatrix(eval.(compiled_M), symbols, res.problem.eom.variables)
end

"""Evaluate the response matrix `resp` for the steady state `s` at (lab-frame) frequency `Ω`."""
function evaluate_response_matrix(resp::ResponseMatrix, s::StateDict, Ω)
    values = cat([s[var] for var in resp.symbols], [Ω], dims=1)
    f = resp.matrix
    return [Base.invokelatest(el, values) for el in f]
end

## THIS NEEDS REVISING
function _evaluate_response_vector(rmat::ResponseMatrix, s::StateDict, Ω)
    m=evaluate_response_matrix(rmat, s, Ω)
    force_pert = cat([[1.0, 1.0*im] for n in 1:size(m)[1]/2]..., dims=1)
    return inv(m) * force_pert
end

"""
$(TYPEDSIGNATURES)

For `rmat` and a solution dictionary `s`,
calculate the total response to a perturbative force at frequency `Ω`.

"""
function get_response(rmat::ResponseMatrix, s::StateDict, Ω)
    resp = 0

    # uv-type
    for pair in _get_uv_pairs(rmat.variables)
        u,v = rmat.variables[pair]
        this_ω = Float64(substitute_all(u.ω, s))
        uv1 = _evaluate_response_vector(rmat, s, Ω-this_ω)[pair]
        uv2 = _evaluate_response_vector(rmat, s, -Ω+this_ω)[pair]
        resp += sqrt(_plusamp(uv1)^2 + _minusamp(uv2)^2)
    end

    # a-type variables
    for a_idx in _get_as(rmat.variables)
        a = rmat.variables[a_idx]
        uv1 = _evaluate_response_vector(rmat, s, Ω)[a]
        uv2 = _evaluate_response_vector(rmat, s, -Ω)[a]
        resp += sqrt(_plusamp(uv1)^2 + _minusamp(uv2)^2)
    end
    resp
end

# formulas to obtain up- and down- converted frequency components when going from the
# rotating frame into the lab frame
_plusamp(uv) = norm(uv)^2 - 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))
_minusamp(uv) = norm(uv)^2 + 2*(imag(uv[1])*real(uv[2]) - real(uv[1])*imag(uv[2]))

"""
$(TYPEDSIGNATURES)

Plot the linear response of a solution `branch` of `res` for the frequencies `Ω_range`.

This method takes n^2 matrix inversions for n elements of `Ω_range`. It therefore takes longer than `plot_jacobian_spectrum` (which is (O(n)))
but is also valid far from resonance.

"""
function plot_response(res::Result, Ω_range; branch::Int, logscale=false)
    _set_plotting_settings()
    length(size(res.solutions)) != 1 && error("1D plots of not-1D datasets are usually a bad idea.")
    stability = classify_branch(res, branch, "stable") # boolean array
    !any(stability) && error("Cannot generate a spectrum - no stable solutions!")

    X = Vector{Float64}(collect(values(res.swept_parameters))[1][stability])

    response = ResponseMatrix(res)

    Y = Array{Float64, 2}(undef,  length(Ω_range), length(X))

    # this could be optimized by not grabbing the entire huge dictionary every time

    bar = Progress(length(Y), 1, "Solving the linear response ODE for each solution and input frequency ... ", 50)
    for j in 1:(size(Y)[2])
        next!(bar)
        s = get_single_solution(res, branch=branch, index=j)
        for i in 1:(size(Y)[1])
            Y[i,j] = get_response(response, s, reverse(Ω_range)[i])
        end
    end

    Y = logscale ? log.(Y) : Y
    plt = PyPlot.imshow(Y, aspect="auto", extent=[X[1], X[end], Ω_range[1], Ω_range[end]])
    xlabel(Latexify.latexify(string(first(keys(res.swept_parameters)))),fontsize=20); ylabel("noise angular frequency",fontsize=20);
    return plt
end


text:
### A Pluto.jl notebook ###
# v0.12.11

using Markdown
using InteractiveUtils

# ╔═╡ 1c7e00f8-330d-11eb-36a3-fb25098f61f0
begin
	import Pkg; Pkg.add("RollingFunctions")
	using LinearAlgebra
	using Plots
	using QuantEcon
	using RollingFunctions
end

# ╔═╡ 9bb7e430-3309-11eb-045b-45984b90b650
md"""
# Arrays, Tuples, Ranges, and Other Fundamental Types
### Solutions to the Exercise
"""

# ╔═╡ 6b36bba0-330a-11eb-2f33-6d4b9f648302

html"""Exercise 1 (click)"""

# ╔═╡ c9846d96-330b-11eb-05ec-7774cfa3e09e
function ComputedUnconditionVariables(A::AbstractArray,Σ::AbstractArray;max_err =1e-7,maxitr=1000)

	S=zeros(size(A))
	i=0
	while true
		i+=1
		newS = A*S*A' + Σ*Σ'
		if norm(newS-S) < max_err && i>maxitr
			break
		end
		S = newS
	end
	S
end

# ╔═╡ 60fe778c-330e-11eb-3efd-37fe60e5becf
let
	A =[0.8 -0.2;-0.1 0.7]
	Σ = [0.5 0.4;0.4 0.6]
	norm(ComputedUnconditionVariables(A,Σ)-solve_discrete_lyapunov(A,Σ*Σ'))
end

# ╔═╡ 3ea21b48-330f-11eb-1fab-c396a34b1089
html"""
Exercise 2 (click)"""

# ╔═╡ 382c3762-3312-11eb-3371-f71f6337cd13
Θ =[0.8,0.9,0.98]

# ╔═╡ 75bf8eb2-330f-11eb-2313-fb73eb133cbb
function CalculateY(T::Int64,θ::Float64;σ=1,γ=1,y=0)
	A=Array{Float64}(undef,0)
	for i in 1:T
		y = γ + θ*y +σ*randn()
		push!(A,y)
	end
	A
end

# ╔═╡ c8e9e352-3310-11eb-35dd-f72acb3b5373
function CalculateY(T::Int64;σ=1,γ=1,y=0)
	ans=[]
	for θ in Θ
		push!(ans,CalculateY(T,θ;σ=1,γ=1,y=0))
	end
	ans
end

# ╔═╡ f53d440c-3311-11eb-217a-2f5e52af604b
let
	res = CalculateY(150)
	p=plot(1:150,res)
	rollingres = rolling.(mean,res,10)
	plot!(p,1:150,cat.(Ref(zeros(9)),rollingres,dims=1),legend=false)
end

# ╔═╡ 75978c48-3317-11eb-09b8-d33d9c09cb4f
let
	answers = Dict()
	for θ in Θ
		answers[θ]=Array{Float64}(undef,0)
	end
	for N in 1:200
		p=CalculateY(150)
		for (i,v) in enumerate(Θ)
			push!(answers[v],p[i][end])
		end
	end
	all=cat(collect(values(answers))...,dims=1)
	t1=mean(all)# answer for last part doing it for values of Θ collectively
	t2=mean(x->(x^2),all)# answer of last part
	@show t1
	@show t2-t1^2
	histogram(answers;alpha=0.5,label = keys(answers))
end

# ╔═╡ 1365fe9c-3321-11eb-1b03-1bcb8ec30dc8
html"""
Exercise 3 (click)"""

# ╔═╡ 2df0c788-3321-11eb-08c3-cf636bb7d128
begin
	parameters=[0.1,0.2,0.5,1.0,1.0]
	Y=[]
	X=[]
	for M in 1:20
		x1=randn(50)
		x2=randn(50)
		w=randn(50)
		push!(X,cat(x1,x1.^2,x2,ones(50),w,dims=2))
		push!(Y,X[end]*parameters)
	end
	X,Y
end

# ╔═╡ d44bb630-3324-11eb-1dd8-47e7ce22a550
function find_parameters(x::AbstractArray,y::AbstractArray)
	(x'*x)\(x'*y)
end

# ╔═╡ 5f790744-3325-11eb-2629-f52ace14e7bb
para=[ find_parameters(X[M],Y[M]) for M in 1:20]

# ╔═╡ 108c296c-3326-11eb-04a8-3bdb958ddeef
let
	para2=Matrix(reshape(cat(para...,dims=1),(5,20))')# for display
	histogram(para2,bins=0:0.1:1.1;alpha=.2,label=['a','b','c','d','σ'])
end

# ╔═╡ 518146e4-330a-11eb-1720-d35a5e832af2
bigbreak=html"""

"""

# ╔═╡ 697a48ea-330f-11eb-1650-1fb81c301c71
bigbreak

# ╔═╡ 2fe0b666-3321-11eb-0598-eb28cd43af60
bigbreak

# ╔═╡ Cell order:
# ╟─9bb7e430-3309-11eb-045b-45984b90b650
# ╟─6b36bba0-330a-11eb-2f33-6d4b9f648302
# ╠═c9846d96-330b-11eb-05ec-7774cfa3e09e
# ╠═60fe778c-330e-11eb-3efd-37fe60e5becf
# ╟─697a48ea-330f-11eb-1650-1fb81c301c71
# ╟─3ea21b48-330f-11eb-1fab-c396a34b1089
# ╠═382c3762-3312-11eb-3371-f71f6337cd13
# ╠═75bf8eb2-330f-11eb-2313-fb73eb133cbb
# ╠═c8e9e352-3310-11eb-35dd-f72acb3b5373
# ╟─f53d440c-3311-11eb-217a-2f5e52af604b
# ╠═75978c48-3317-11eb-09b8-d33d9c09cb4f
# ╟─2fe0b666-3321-11eb-0598-eb28cd43af60
# ╟─1365fe9c-3321-11eb-1b03-1bcb8ec30dc8
# ╠═2df0c788-3321-11eb-08c3-cf636bb7d128
# ╠═d44bb630-3324-11eb-1dd8-47e7ce22a550
# ╠═5f790744-3325-11eb-2629-f52ace14e7bb
# ╠═108c296c-3326-11eb-04a8-3bdb958ddeef
# ╟─1c7e00f8-330d-11eb-36a3-fb25098f61f0
# ╟─518146e4-330a-11eb-1720-d35a5e832af2

max_stars_repo_path: getting started with julia/Arrays, Tuples, Ranges, and Other Fundamental Types.jl

max_stars_repo_name: nehamsoni/Quantitative-Economics

max_issues_repo_name: nehamsoni/Quantitative-Economics

max_forks_repo_name: nehamsoni/Quantitative-Economics