text:
module TestingZippedArrays

using Test
using ZippedArrays
using ZippedArrays: throw_indices_mismatch, all_match

function generate_array(::Type{T}, dims::Dims{N}) where {T,N}
    A = Array{T,N}(undef, dims)
    for i in 1:length(A)
        A[i] = convert(T, i)
    end
    return A
end

@testset "zipped arrays" begin
    dims = (2,3,4)
    A = generate_array(Int32, dims)
    B = generate_array(Float64, dims)
    C = generate_array(Char, dims)
    D = generate_array(Int16, (1, dims...))
    E = generate_array(Float32, (dims .+ 1))
    V = view(E, map(d -> 2:d+1, dims)...)

    for bool in (true, false)
        @test_throws DimensionMismatch throw_indices_mismatch(bool, A, B, E)
    end
    @test_throws ErrorException ZippedArray()
    @test_throws MethodError ZippedArray(A,D)
    @test_throws DimensionMismatch ZippedArray(A,E)
    @test all_match(nothing, cos)

    # Zip 1 array.
    Z = ZippedArray(D)
    @test IndexStyle(Z) === IndexStyle(D)
    @test eltype(Z) === Tuple{eltype(D)}
    @test ndims(Z) == ndims(D)
    @test size(Z) == size(D)
    @test size(Z,1) == size(D,1)
    @test axes(Z) == axes(D)
    @test axes(Z,2) == axes(D,2)
    @test setindex!(Z, (11,), 1,2,3,1) === Z
    for i in (12, CartesianIndex(1,1,2,3))
        x = Z[i] # save values
        @test Z[i] === (D[i],)
        Z[i] = (D[i]+1,)
        @test D[i] == x[1]+1
    end
    @test map(x -> x[1], Z) == D
    @test_throws BoundsError Z[-1]

    # Zip 2 regular arrays.
    Z = ZippedArray(A,B)
    @test IndexStyle(Z) === IndexLinear()
    @test eltype(Z) === Tuple{eltype(A),eltype(B)}
    @test ndims(Z) == ndims(A)
    @test size(Z) == size(A)
    @test size(Z,1) == size(A,1)
    @test axes(Z) == axes(A)
    @test axes(Z,2) == axes(A,2)
    @test setindex!(Z, (7, 2.3), 2,3,1) === Z
    for i in (12, CartesianIndex(1,2,3))
        x = Z[i] # save values
        @test Z[i] === (A[i], B[i])
        Z[i] = (A[i]+1, B[i]+2)
        @test A[i] == x[1]+1
        @test B[i] == x[2]+2
    end
    @test map(x -> x[1], Z) == A
    @test map(x -> x[2], Z) == B
    @test_throws BoundsError Z[-1]

    # Zip 2 regular arrays and a view.
    Z = ZippedArray(A,V,C)
    @test IndexStyle(Z) === IndexCartesian()
    @test eltype(Z) === Tuple{eltype(A),eltype(V),eltype(C)}
    @test ndims(Z) == ndims(A)
    @test size(Z) == size(A)
    @test size(Z,1) == size(A,1)
    @test axes(Z) == axes(A)
    @test axes(Z,2) == axes(A,2)
    @test setindex!(Z, (-1, 0.3, 'x'), 1,2,3) === Z
    for i in (12, CartesianIndex(1,2,3))
        x = Z[i] # save values
        @test Z[i] === (A[i], V[i], C[i])
        Z[i] = (A[i]+1, V[i]+2, C[i]+3)
        @test A[i] == x[1]+1
        @test V[i] == x[2]+2
        @test C[i] == x[3]+3
    end
    @test map(x -> x[1], Z) == A
    @test map(x -> x[2], Z) == V
    @test map(x -> x[3], Z) == C
    @test_throws BoundsError Z[1,2,-1]

    # This works for vectors.
    n = 23
    a = rand(Bool, n)
    b = rand(Float32, n)
    c = rand(Int16, n)
    @test [(a[i],b[i],c[i]) for i in 1:n] == ZippedArray(a,b,c)

    sort!(ZippedArray(a,b);
          lt = (x,y) -> ifelse(x[1] == y[1], x[2] < y[2], x[1] < y[1]))
    flag = true
    for i in 2:n
        flag &= ((a[i-1] < a[i])|((a[i-1] == a[i])&(b[i-1] <= b[i])))
    end
    @test flag
end

end


text:
using Brillouin, Test
using Brillouin.KPaths: cartesianize, latticize, Bravais

@testset "KPath & KPathInterpolant" begin
    # --- `cumdist` ---
    kvs_2d = [[0,0], [0,1], [1,1], [1,0], [0,0], [1,1], [-1,-2], [2,-1]]
    @test cumdists(kvs_2d) ≈ [0,1,2,3,4,4+√2,4+√2+√13, 4+√2+√13+√10]
    kvs_3d = [[2,1.5,3], [4,5,6], [-1.1, 2, -2.3]]
    @test cumdists(kvs_3d) ≈ [0, √25.25, √25.25+√103.9]

    # --- `KPath` ---
    sgnum = 227 # face-centered cubic (cF)
    Rsᶜᶸᵇⁱᶜ = [[1,0,0], [0,1,0], [0,0,1]] # conventional direct basis
    kp = irrfbz_path(sgnum, Rsᶜᶸᵇⁱᶜ)

    # test that `Rs` has no impact in sgnum 227 (cubic face-centered, FCC)
    @test kp == irrfbz_path(sgnum, [rand(3) for _ in 1:3])

    # `getindex`
    @test collect(kp) == [:Γ => [0.0, 0.0, 0.0]
                          :X => [0.5, 0.0, 0.5]
                          :U => [0.625, 0.25, 0.625]
                          :K => [0.375, 0.375, 0.75]
                          :Γ => [0.0, 0.0, 0.0]
                          :L => [0.5, 0.5, 0.5]
                          :W => [0.5, 0.25, 0.75]
                          :X => [0.5, 0.0, 0.5]]

    # `paths`
    @test paths(kp) == [[:Γ, :X, :U], [:K, :Γ, :L, :W, :X]]

    # `show(::IO, ::MIME"text/plain", ::KPath)`
    kp_show_reference = """
    KPath{3} (6 points, 2 paths, 8 points in paths):
     points: :U => [0.625, 0.25, 0.625]
             :W => [0.5, 0.25, 0.75]
             :K => [0.375, 0.375, 0.75]
             :Γ => [0.0, 0.0, 0.0]
             :L => [0.5, 0.5, 0.5]
             :X => [0.5, 0.0, 0.5]
      paths: [:Γ, :X, :U]
             [:K, :Γ, :L, :W, :X]
      basis: [-6.283185, 6.283185, 6.283185]
             [6.283185, -6.283185, 6.283185]
             [6.283185, 6.283185, -6.283185]"""
    test_show(kp, kp_show_reference)

    # `cartesianize` & `latticize`
    pGs = 2π.*[[-1,1,1], [1,-1,1], [1,1,-1]] # primitive reciprocal basis
    kp′ = Brillouin.cartesianize(kp)
    @test all(zip(points(kp), points(kp′))) do ((klab, kv), (klab′, kv′))
        klab == klab′ && kv'pGs ≈ kv′
    end
    kp′′ = Brillouin.latticize(kp′)
    @test keys(points(kp′′)) == keys(points(kp))          # point labels agree
    @test paths(kp′′) == paths(kp)                        # segments agree
    @test all(values(points(kp′′)) .≈ values(points(kp))) # point coordinates approx. agree
    kp′′′ = Brillouin.cartesianize!(Brillouin.latticize(kp′, pGs[[2,3,1]]))
    @test keys(points(kp′′′)) == keys(points(kp′))
    @test paths(kp′′′) == paths(kp′)
    @test all(values(points(kp′′′)) .≈ values(points(kp′)))

    # `extended_bravais` and `irrfbz_path`
    sgnums = (1:2, 1:17, 1:230)
    for D in (1, 2, 3)
        Dᵛ = Val(D)
        for sgnum in sgnums[D]
            Rs = begin
                # the 3D triclinic case (space group 1 & 2) needs more care, since Brillouin
                # only allows all-obtuse or all-acute bases; just hardcode two examples
                if sgnum == 1 && D == 3    # check all-obtuse triclinic `Rs`
                    Bravais.DirectBasis([1.0, 0.0, 0.0], [1.18, 1.54, 0.0], [0.95, 0.44, 1.32])
                elseif sgnum == 2 && D == 3 # check all-acute triclinic `Rs`
                    Bravais.DirectBasis([1.0, 0.0, 0.0], [-0.53, 0.49, 0.0], [0.13, -0.30, 0.74])
                else                        # go with `directbasis` for everything else
                    Bravais.directbasis(sgnum, Dᵛ)
                end
            end
            bt = Bravais.bravaistype(sgnum, D, normalize=false)
            # extended Bravais types
            ebt = Brillouin.KPaths.extended_bravais(sgnum, bt, Rs, Dᵛ)
            @test contains(string(ebt), bt)

            # hard to test output of `irrfbz_path` systematically; just test that it returns
            # a `KPath`, and matches `unshuffle_hpkot_setting! ∘ get_points` and `get_paths`
            kp = irrfbz_path(sgnum, Rs)
            @test kp isa KPath{D}
            @test points(kp) == Brillouin.KPaths.unshuffle_hpkot_setting!(
                                    Brillouin.KPaths.get_points(ebt, Rs, Dᵛ), bt, D)
            @test paths(kp)  == Brillouin.KPaths.get_paths(ebt, Dᵛ)
        end
    end
    @test_throws DomainError Brillouin.KPaths.extended_bravais(110, "Q", nothing, Val(3))  # "undefined bravais type"
    @test_throws DomainError Brillouin.KPaths.extended_bravais(194, "cP", nothing, Val(3)) # `_throw_conflicting_sgnum_and_bravais`
    @test_throws DomainError Brillouin.KPaths.extended_bravais(194, "cF", nothing, Val(3)) # `_throw_conflicting_sgnum_and_bravais`
    @test_throws DomainError Brillouin.KPaths.extended_bravais(38, "tI", nothing, Val(3))  # `_throw_basis_required`
    Rs′ = Bravais.DirectBasis([1, 0, 0], [0.3, 0.8, 0], [-1.6, 0.8, 0.9]) # neither all-obtuse nor all-acute
    @test_throws DomainError Brillouin.KPaths.extended_bravais(1, "aP", Rs′, Val(3))       # `_throw_basis_required`

    # --- `KPathInterpolant` ---
    # `interpolate`
    N = 100
    kpi = interpolate(kp, N)
    @test abs(length(kpi) - N) ≤ 1 # allow deviation by 1 point in total (too tight?)

    # all points in `kp` must be explicitly contained in `kpi`
    @test all(∈(kpi), values(points(kp)))
    @test kpi == interpolate(kp, length=N)

    # iterate vs. indexing
    @test collect(kpi) == [kpi[i] for i in eachindex(kpi)]

    # `splice`
    kps = splice(kp, N)
    highsym_points = sum(path->length(path), paths(kp))
    segments       = sum(path->length(path)-1, paths(kp))
    @test length(kps) == highsym_points + segments*N

    # `cartesianize`
    # test that `cartesianize` commutes with `interpolate` when applied to a `kp`/`kpi`
    cntr = Bravais.centering(5, 2)         # centering type 'c' in plane group 5
    Rs   = Bravais.directbasis(5, Val(2))
    Gs   = Bravais.reciprocalbasis(Rs)
    pGs  = Bravais.primitivize(Gs, cntr)

    kp   = irrfbz_path(5, Rs, Val(2))
    kpi  = interpolate(Brillouin.cartesianize(kp), 100)
    kpi′ = Brillouin.cartesianize!(interpolate(kp, 100))

    @test typeof(kpi) === typeof(kpi′) === KPathInterpolant{2}
    @test kpi′ ≈ kpi # test that `interpolate` commutes with `latticize`/`cartesianize`
    @test Brillouin.cartesianize!(latticize!(kpi)) ≈ kpi′

    @test Brillouin.cartesianize!(Brillouin.latticize(kpi, pGs[[2,1]])) ≈ kpi′

    # interpolate via `density` kwarg
    kp = irrfbz_path(227, Rsᶜᶸᵇⁱᶜ)
    for ρ in [1, 10, 20, 21, 100, 3.5]
        kpi_ρ = interpolate(kp, density=ρ)

        # test a lower bound on how many points must be in kpi
        segments = sum(p -> length(p) - 1, kp.paths)
        tot_dist = sum(last.(cumdists.(Brillouin.cartesianize(kpi_ρ).kpaths)))
        min_points = tot_dist*ρ - segments
        @test length(kpi_ρ) > min_points

        # TODO: test that the density of k-points is _at least_ equal to ρ at every segment
    end

    # interpolate with wrong kwargs
    @test_throws ArgumentError interpolate(kp)
    @test_throws ArgumentError interpolate(kp; length=10, density=2)
end # testset

ext: jl


text:
# read_samples

"""

Read sample output files created by StanSample.jl.

This method is added to StanRun's read_sample function.

### Method
```julia
read_samples(model::SampleModel; start=1)
```

### Required arguments
```julia
* `model`                    : SampleModel
```

### Optional arguments
```julia
* `start=1`                  : First draw saved in MCMCHains.Chains object
```

"""
function read_samples(model::SampleModel;  start=1)

  local a3d, monitors, index, idx, indvec, ftype, noofsamples

  output_base = model.output_base
  name_base ="_chain"
  n_samples = model.method. num_samples
  n_chains = StanBase.get_n_chains(model)

  # Handle save_warmup
  start = model.method.save_warmup ? model.method.num_warmup+1 : start

  # a3d will contain the samples such that a3d[s, i, c] where

  #   s: num_samples
  #   i: variables (from cmdstan .csv file)
  #   c: n_chains

  # Read .csv files created by each chain

  for i in 1:n_chains
    if isfile(output_base*name_base*"_$(i).csv")
      #noofsamples = 0
      instream = open(output_base*name_base*"_$(i).csv")
      #
      # Skip initial set of commented lines, e.g. containing cmdstan version info, etc.
      #
      skipchars(isspace, instream, linecomment='#')
      #
      # First non-comment line contains names of variables
      #
      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]
      indvec = 1:length(index)

      if i == 1
        a3d = fill(0.0, n_samples, length(indvec), n_chains)
      end

      #println(size(a3d))
      skipchars(isspace, instream, linecomment='#')
      for j in 1:n_samples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = parse.(Float64, split(strip(line), ","))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end   # read in samples
    end   # read in next file
  end   # read in file for each chain

  cnames = convert.(String, idx[indvec])
  chns = convert_a3d(a3d, cnames, Val(:mcmcchains); start=start)

end   # end of read_samples

function read_samples(model::StanModel; chain=1)
  read_samples(default_output_base(model)*"_chain_$(chain).csv")
end

"""

# convert_a3d

Convert the output file created by cmdstan to the shape of choice.

### Method
```julia
convert_a3d(a3d_array, cnames, ::Val{Symbol}; start=1)
```
### Required arguments
```julia
* `a3d_array::Array{Float64, 3},`      : Read in from output files created by cmdstan
* `cnames::Vector{AbstractString}`     : Monitored variable names
```

### Optional arguments
```julia
* `::Val{Symbol}`                      : Output format
* `::start=1`                          : First draw for MCMCChains.Chains
```
Method called is based on the output_format defined in the stanmodel, e.g.:

   stanmodel = Stanmodel(`num_samples`=1200, thin=2, name="bernoulli",
     model=bernoullimodel, `output_format`=:mcmcchains);

Current formats supported are:

1. :array (a3d_array format, the default for CmdStan)
2. :namedarray (NamedArrays object)
3. :dataframe (DataFrames object)
4. :mambachains (Mamba.Chains object)
5. :mcmcchains (TuringLang/MCMCChains.Chains object)

Options 3 through 5 are respectively provided by the packages StanDataFrames,
StanMamba, StanMCMCChains and StanMCMCChains.
```

### Return values
```julia
* `res`                       : Draws converted to the specified format.
```
"""
convert_a3d(a3d_array, cnames, ::Val{:array}; start=1) = a3d_array

convert_a3d(a3d_array, cnames, ::Val{:namedarray}; start=1) =
  [NamedArray(a3d_array[:,:,i], (collect(1:size(a3d_array, 1)), Symbol.(cnames)))
    for i in 1:size(a3d_array, 3)]

function convert_a3d(a3d_array, cnames, ::Val{:mcmcchains}; start=1)
  pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
  p = filter(p -> !(p in  pi), cnames)

  MCMCChains.Chains(a3d_array,
    cnames,
    Dict(
      :parameters => p,
      :internals => pi
    );
    start=start
  )
end


text:
# Generates exponential mesh of N elements on [r_mn, r_max]
#
# The domain [r_min, r_max], the mesh will contain both endpoints:
function gen_rmesh_exp(
    r_min::Float64,
    r_max::Float64,
    a::Float64,
    Nr::Int64
)
    N = Nr - 1

    rmesh = zeros(Float64,N+1)
    SMALL = eps()

    if a < 0.0
        error("gen_rmesh_exp: a > 0 required")
    #
    elseif abs(a - 1.0) < SMALL
        α = (r_max - r_min) / N
        for i in 1:N+1
            rmesh[i] = α * (i - 1) + r_min
        end
    #
    else
        if N > 1
            β = log(a)/(N-1)
            α = (r_max - r_min) / (exp(β*N) - 1)
            for i in 1:N+1
                rmesh[i] = α * ( exp(β*(i-1)) - 1 ) + r_min
            end
        elseif N == 1
            rmesh[1] = r_min
            rmesh[2] = r_max
        else
            error("gen_rmesh_exp: N >= 1 required")
        end
    end
    return rmesh
end


text:
using SatelliteSpots
using Test

@testset "SatelliteSpots.jl" begin
    # Write your tests here.
end


text:
using CAGD;
using SparseArrays;
using LinearAlgebraicRepresentation;
Lar = LinearAlgebraicRepresentation;

G = [
	0.5310492999999998 1.0146684 0.3477716 0.8313907882395298 0.6061046999999998 1.0897237999999998 0.42282699999999984 0.9064461979373597 0.5310493 1.0146684000000001 0.6061047 1.0897237772434623 0.3477716 0.8313908 0.422827 0.9064462 0.5310493 0.34777160000000007 0.6061047 0.4228270000000002 1.0146684 0.8313908 1.0897238 0.9064461675456482;
	0.8659989999999999 0.6827212999999999 0.5268921 0.3436144447063971 1.2188832999999994 1.0356056999999999 0.8797763999999998 0.6964987903021808 0.8659989999999999 0.6827213 1.2188833 1.035605657895053 0.5268921 0.3436145 0.8797764000000001 0.6964988 0.8659989999999999 0.5268920999999999 1.2188833 0.8797764 0.6827213 0.3436145 1.0356057 0.6964988122992563;
	0.14191280000000003 0.2169682 0.4947971000000001 0.5698524407571428 0.5200012 0.5950565999999999 0.8728855 0.9479408896095312 0.14191279999999987 0.21696819999999994 0.5200011999999999 0.5950566438156151 0.4947971 0.5698525000000001 0.8728855 0.9479409 0.14191280000000006 0.4947971000000001 0.5200012000000002 0.8728855000000001 0.21696819999999994 0.5698525000000001 0.5950565999999999 0.9479408949632379
];
W = convert(Lar.Points, G')
Delta_0 = SparseArrays.sparse((
	[1, 3, 1, 4, 2, 3, 2, 4, 5, 7, 5, 8, 6, 7, 6, 8, 9, 11, 9, 12, 10, 11, 10, 12, 13, 15, 13, 16, 14, 15, 14, 16, 17, 19, 17, 20, 18, 19, 18, 20, 21, 23, 21, 24, 22, 23, 22, 24],
	[1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24],
	Int8[-1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1]
)...);
Delta_1 = SparseArrays.sparse((
	[1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6],
	[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24],
	Int8[1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1]
)...);
T = [Delta_0, Delta_1];

model = CAGD.Model(G);
CAGD.addModelCells!(model, 1, T[1]);
CAGD.addModelCells!(model, 2, T[2]);

split_model = CAGD.facesplitting(model)
congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)
gift_model = deepcopy(congr_model);
FC, bicon_comps = CAGD.tgw(congr_model, 3)
CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))


text:
using Test
@testset "yeah" begin
    @test true
    sleep(0.5)
    @test true
    sleep(0.5)
end


text:
#Array{Float64,1}
a=[168.3,179.2,165.8,170.5,188.2,174.6,162.8,175.5,178.1,177.1]

#get arithmetic mean
@show mean(a)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function event_is_enabled(obj::EnablingPredicate, arg0::SpacecraftState, arg1::EventDetector, arg2::jdouble)
    return jcall(obj, "eventIsEnabled", jboolean, (SpacecraftState, EventDetector, jdouble), arg0, arg1, arg2)
end


text:
using DSGE, Plots, FileIO
using Plots.PlotMeasures

# This script loads in matrices from the DSGE in Del Negro, Schorfheide, Smets, and Wouters (2007)
# and computes the DSGE impulse responses to exogenous structural shocks and the
# corresponding DSGE-VECM impulse responses with 4 lags using either λ = ∞ or
# the rotation matrix from the DSGE to identify the
# mapping from the innovations in the VECM to the structural shocks.
# The script creates a dictionary called `plots_dict`, whose keys are the names of
# different shocks and values are the plots of the impulse responses of different
# observables to the shock.
#
# Generally, the DSGE-VECM IRFs with λ = 0.33 will look different from the DSGE IRFs, but
# they will yield an overall similar response. See the response of output growth to
# the z shock. In contrast, the λ = ∞ DSGE-VECM IRFs are very close to the DSGE IRFs.
# See the IRF of output growth to the μ or λ_f shocks for examples of DSGE-VECM IRFs
# with λ = ∞ that are nearly identical to the underlying DSGE IRFs. For IRFs that are
# not exactly the same, see the IRF of output growth to the z shock. On impact, the
# IRF is different, but after 5-10 quarters, the impulse responses are almost identical.
fp = dirname(@__FILE__)

matdata     = load(joinpath(fp, "../test/reference/dsgevecm_lambda_irfs.jld2"))
horizon     = 20
nshocks     = size(matdata["RRR"], 2)
obs_dict    = Dict{Symbol, Int}(:output_growth => 1, :consumption_growth => 2,
                                :investment_growth => 3, :hours => 4,
                                :wage_growth => 5, :inflation => 6, :nominalrate => 7)
# The shocks are technology, labor preference, investment price, preference,
# government spending, price markup, and monetary policy
shocks_dict = Dict{Symbol, Int}(:z => 1, :φ => 2, :μ => 3, :b => 4, :g => 5,
                                :λ_f => 6, :r => 7)

# DSGEVECM IRFs in deviations with λ = 0.33
ŷ_dsgevecm_λ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                    vec(matdata["DD"]), matdata["MM"],
                                    matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                    Int(matdata["coint"]), matdata["cct_sim"], matdata["sig_sim"],
                                    matdata["cointvec"], horizon, deviations = true)

# DSGEVECM IRFs with λ = ∞, p = 12
β_λ_∞, Σ_λ_∞ = DSGE.vecm_approx_state_space(matdata["TTT"], matdata["RRR"], matdata["QQ"],
                                       vec(matdata["DD"]), matdata["ZZ"],
                                       zeros(Int(matdata["nvar"]) + Int(matdata["coint"]),
                                             Int(matdata["nvar"]) + Int(matdata["coint"])),
                                       matdata["MM"], Int(matdata["nvar"]), 4, Int(matdata["coint"]),
                                       use_intercept = true)
Σ_λ_∞ = (Σ_λ_∞ + Σ_λ_∞') ./ 2
ŷ_dsgevecm_∞ = DSGE.impulse_responses(matdata["TTT"], matdata["RRR"], matdata["ZZ"],
                                      vec(matdata["DD"]), matdata["MM"],
                                      matdata["QQ"], Int(matdata["k"]), Int(matdata["nvar"]),
                                      Int(matdata["coint"]), β_λ_∞, Σ_λ_∞,
                                      matdata["cointvec"], horizon, deviations = true)

# DSGE IRFs in deviations
dsge_states = zeros(size(matdata["TTT"], 1), horizon, nshocks)
ŷ_dsge = zeros(Int(matdata["nvar"]), horizon, nshocks)
for i in 1:nshocks
    shock = zeros(nshocks)
    shock[i] = -sqrt(matdata["QQ"][i, i])
    dsge_states[:, 1, i] = matdata["RRR"] * shock
    for t in 2:horizon
        dsge_states[:, t, i] = matdata["TTT"] * dsge_states[:, t - 1, i]
    end
    ŷ_dsge[:, :, i] = matdata["ZZ"][1:Int(matdata["nvar"]), :] * dsge_states[:, :, i]
end

# Plot!
plots_dict = Dict{Symbol, Dict{Symbol, Plots.Plot}}()
for (shock, j) in  shocks_dict
    plots_dict[shock] = Dict{Symbol, Plots.Plot}()
    for (obs, i) in obs_dict
        plots_dict[shock][obs] = plot(1:horizon, ŷ_dsgevecm_∞[i, :, j],
                                      label = "DSGE-VECM \\lambda = Inf", color = :black,
                                      linewidth = 3)
        plot!(1:horizon, ŷ_dsge[i, :, j], label = "DSGE",
              color = :red, linewidth = 3, linestyle = :dash)
        plot!(1:horizon, ŷ_dsgevecm_λ[i, :, j], label = "DSGE-VECM \\lambda = 0.33", color = :black,
              linestyle = :dash,
              legend = :bottomright, left_margin = 20px)
    end
end


text:
using Documenter, NetCDF

makedocs(
    modules = [NetCDF],
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliageo.github.io/NetCDF.jl",
        assets = String[],
    ),
    sitename = "NetCDF.jl",
    authors = "Fabian Gans and contributors",
    pages = [
        "Home" => "index.md",
        "Manual" => ["quickstart.md", "highlevel.md", "intermediate.md", "strings.md"],
    ],
)

deploydocs(repo = "github.com/JuliaGeo/NetCDF.jl")


text:
using AngleBetweenVectors
using Test

point1 = (1.0, 0.0)
point2 = (0.0, 1.0)
point3 = (-1.0, 0.0)

@test angle(point1, point2) == angle(point2, point1)
@test angle(point1, point3) == angle(point3, point1)
@test angle(point3, point2) == angle(point2, point3)

@test angle(point1, point2) == pi/2.0
@test angle(point1, point3) == pi/1.0

point1 = (1.0, -1.0, 0.0,  0.0)
point2 = (0.0,  1.0, 0.0, -1.0)

# note: 2pi/3 !== Float64(2*BigFloat(pi)/3)
@test angle(point1, point2) == Float64(2*BigFloat(pi)/3)
@test angle(point1, point2) == angle(point2, point1)


text:
# # HMM with continuous observations

using Distributions
using HMMBase
using Plots

a = [0.6, 0.4]
A = [0.7 0.3; 0.4 0.6]
B = [MvNormal([0.0,5.0],[1.0,1.0]), MvNormal([5.0,10.0],[1.0,1.0])]
hmm = HMM(a, A, B)
#-

z, y = rand(hmm, 250, seq = true)
z_viterbi = viterbi(hmm, y);
#-

plot(y, linetype=:steppre, label=["Observations (1)", "Observations (2)"], size=(600,200))
#-

plot(z, linetype=:steppre, label="True hidden state", size=(600,200))
#-

plot(z_viterbi, linetype=:steppre, label="Viterbi decoded hidden state", size=(600,200))


text:
using ABCDistances
Libdl.dlopen("/usr/lib/liblapack.so.3", Libdl.RTLD_GLOBAL); ##Needed to avoid PyPlot problems on my work machine
##I sometimes also need run "export LD_LIBRARY_PATH=$HOME/.julia/v0.6/Conda/deps/usr/lib/:$LD_LIBRARY_PATH" - see https://github.com/JuliaPy/Conda.jl/issues/105
using PyPlot
using Distributions
import Distributions.length, Distributions._rand!, Distributions._pdf ##So that these can be extended

P = [ 1 0; ##Prey growth; requires 1 unit of prey
      1 1; ##Predation; requires 1 predator and 1 prey
      0 1] ##Predator death; requires 1 unit of predator
Q = [ 2 0; ##Prey growth; 2 units of prey produced
      0 2; ##Predation; 2 units of predator produced
      0 0] ##Predator death; nothing produced

stoichiometry_LV = Stoichiometry(P, Q)

##An illustration of the complete data simulation function
srand(1)
(LV_times, LV_states) = gillespie_sim(stoichiometry_LV, [50, 100], [0.5,0.0025,0.3], 100.0, 10000)

PyPlot.figure();
PyPlot.subplot(121);
plot(LV_states[1,:], LV_states[2,:])
PyPlot.title("Trajectory")
PyPlot.xlabel("Prey");
PyPlot.ylabel("Predators");
PyPlot.subplot(122);
plot(LV_times, LV_states[1,:]);
plot(LV_times, LV_states[2,:]);
PyPlot.legend(["Prey", "Predators"]);
PyPlot.ylim([0.,510.]); ##Stops legend blocking lines
PyPlot.title("Trace plots");
PyPlot.xlabel("Time");
PyPlot.ylabel("Population");
PyPlot.tight_layout();
PyPlot.savefig("LV_full_data.pdf");

##Simulate observed data as in Owen et al 2014 (their dataset D2), but without time 0
obs_times = [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0]
state0 = [50, 100]
theta0 = [1.0, 0.005, 0.6]
srand(1)
(LV_success, LV_obs) = gillespie_partial_sim(stoichiometry_LV, state0, theta0, obs_times, 100000)
σ0 = exp(2.3) ##Scale of observation noise
nobs = 2*length(obs_times)
x0 = vec(LV_obs') + σ0*randn(nobs); ##Noisy observations
PyPlot.figure();
plot(obs_times, x0[1:16], "-o");
plot(obs_times, x0[17:32], "-o");
PyPlot.legend(["Prey", "Predators"]);
PyPlot.xlabel("Time");
PyPlot.ylabel("Noisy observation");
PyPlot.tight_layout();
PyPlot.savefig("LV_obs.pdf");

####################################
##ABC ON A SINGLE DATASET
##nb initial state and σ known
##(parameters are the log rates)
####################################
##Define prior
marg_rate_prior = Uniform(-6.0,2.0)

type LVprior <: ContinuousMultivariateDistribution
end

function length(d::LVprior)
    3
end

function _rand!{T<:Real}(d::LVprior, x::AbstractVector{T})
    x = rand(marg_rate_prior, 3)
end

function _pdf{T<:Real}(d::LVprior, x::AbstractVector{T})
    if (all(-6.0 .<= x .<= 2.0))
        return 8.0^-3
    else
        return 0.0
    end
end

##Define model
function sample_sumstats(pars::Array{Float64,1})
    (success, x) = gillespie_partial_sim(stoichiometry_LV, state0, exp.(pars), obs_times, 100000)
    if (success)
        stats = vec(x') + σ0*randn(nobs)
    else
        stats = zeros(x)
    end
    (success, stats)
end

abcinput = ABCInput();
abcinput.prior = LVprior();
abcinput.sample_sumstats = sample_sumstats;
abcinput.abcdist = WeightedEuclidean(x0)
abcinput.nsumstats = nobs;

srand(1);
pmcoutput_alg3 = abcPMC3(abcinput, 200, 1/2, 50000, store_init=true);
pmcoutput_alg3V = abcPMC3V(abcinput, 200, 1/2, 50000);
pmcoutput_alg5 = abcPMC5(abcinput, 200, 1/2, 50000, store_init=true);
pmcoutput_alg4 = abcPMC4(abcinput, 200, 1/2, 50000, store_init=true);
abcinput.abcdist = MahalanobisEmp(x0)
pmcoutput_alg5_Mahalanobis = abcPMC5(abcinput, 200, 1/2, 50000, store_init=true);

##Plot weights
w1 = pmcoutput_alg3.abcdists[1].w;
w2 = pmcoutput_alg5.abcdists[pmcoutput_alg5.niterations].w;
w3 = pmcoutput_alg4.abcdists[pmcoutput_alg4.niterations].w;
PyPlot.figure(figsize=(12,6));
PyPlot.subplot(121);
plot(obs_times, w1[1:16]/sum(w1), "b-o");
plot(obs_times, w3[1:16]/sum(w3), "g-^");
plot(obs_times, w2[1:16]/sum(w2), "y-*");
PyPlot.ylim([0.,0.1])
PyPlot.xlabel("Time");
PyPlot.ylabel("Relative weight");
PyPlot.title("Prey");
PyPlot.legend(["Algorithm 3","Algorithm 4\n(last iteration)","Algorithm 5\n(last iteration)"]);
PyPlot.subplot(122);
plot(obs_times, w1[17:32]/sum(w1), "b-o");
plot(obs_times, w3[17:32]/sum(w3), "g-^");
plot(obs_times, w2[17:32]/sum(w2), "y-*");
PyPlot.ylim([0.,0.1])
PyPlot.ylabel("Relative weight");
PyPlot.xlabel("Time");
PyPlot.title("Predators");
PyPlot.legend(["Algorithm 3","Algorithm 4\n(last iteration)","Algorithm 5\n(last iteration)"]);
PyPlot.tight_layout();
PyPlot.savefig("LV_weights.pdf");

##Plot MSEs. First only those algorithms used in paper.
outputs = (pmcoutput_alg3, pmcoutput_alg4, pmcoutput_alg5);
pnames = ("Prey growth", "Predation", "Predator death");
leg_code = ("b-o", "g-^", "y-*");
PyPlot.figure(figsize=(12,6))
for i in 1:length(outputs)
    s = outputs[i]
    m = parameter_means(s)
    v = parameter_vars(s)
    c = s.cusims ./ 1000;
    for j in 1:3
        PyPlot.subplot(1, 3, j)
        PyPlot.plot(c, log10.(v[j,:] .+ (m[j,:]-log(theta0[j])).^2), leg_code[i])
    end
end
for i in 1:3
    PyPlot.subplot(1, 3, i)
    PyPlot.title(pnames[i])
    PyPlot.xlabel("Number of simulations (000s)")
    PyPlot.ylabel("log₁₀(MSE)")
    PyPlot.legend(["Algorithm 3", "Algorithm 4", "Algorithm 5"])
end
PyPlot.tight_layout();
PyPlot.savefig("LV_mse.pdf");

##Add MSEs for other algorithms
outputs = (pmcoutput_alg3V, pmcoutput_alg5_Mahalanobis);
leg_code = ("r-x", "k-|");
for i in 1:length(outputs)
    s = outputs[i]
    m = parameter_means(s)
    v = parameter_vars(s)
    c = s.cusims ./ 1000;
    for j in 1:3
        PyPlot.subplot(1, 3, j)
        PyPlot.plot(c, log10.(v[j,:] .+ (m[j,:]-log(theta0[j])).^2), leg_code[i])
    end
end

##Simulations used for distance initialisation

##Uncomment these rows to compare all 3 outputs
##ss_toplot = Array[pmcoutput_alg3.init_sims[1],
##                  pmcoutput_alg4.init_sims[pmcoutput_alg4.niterations],
##                  pmcoutput_alg5.init_sims[pmcoutput_alg5.niterations]];
##ss_names = ["Algorithm 3\n(first iteration)", "Algorithm 4\n(last iteration)", "Algorithm 5\n(last iteration)"];
ss_toplot = Array[pmcoutput_alg3.init_sims[1],
                  pmcoutput_alg5.init_sims[pmcoutput_alg5.niterations]];
ss_names = ["Algorithm 3\n(first iteration)", "Algorithm 5\n(last iteration)"];

PyPlot.figure(figsize=(12,12));
plotcounter = 1;
for ss in ss_toplot
    PyPlot.subplot(length(ss_toplot), 2, plotcounter)
    for i in 1:20
        plot(obs_times, ss[1:16,i])
    end
    plot(obs_times, x0[1:16], "ko", markersize=5)
    PyPlot.ylim([0,1000])
    plotcounter += 1
    PyPlot.subplot(length(ss_toplot), 2, plotcounter)
    for i in 1:20
        plot(obs_times, ss[17:32,i])
    end
    plot(obs_times, x0[17:32], "ko", markersize=5)
    PyPlot.ylim([0,1000])
    plotcounter += 1
end
PyPlot.subplot(length(ss_toplot), 2, 1)
PyPlot.title("Prey");
PyPlot.subplot(length(ss_toplot), 2, 2)
PyPlot.title("Predators");
PyPlot.subplot(length(ss_toplot), 2, 2*length(ss_toplot)-1)
PyPlot.xlabel("Time");
PyPlot.subplot(length(ss_toplot), 2, 2*length(ss_toplot))
PyPlot.xlabel("Time");
for i in 1:length(ss_toplot)
    PyPlot.subplot(length(ss_toplot), 2, 2i-1)
    PyPlot.ylabel("Population")
    line_invis = plot((0,0),(0,0),"w-")
    PyPlot.legend([line_invis], [ss_names[i]], handlelength=0)
    PyPlot.subplot(length(ss_toplot), 2, 2i)
    PyPlot.legend([line_invis], [ss_names[i]], handlelength=0)
end
PyPlot.tight_layout();
PyPlot.savefig("LV_paths.pdf");

##Marginal posterior plots
PyPlot.figure();
samplesABC = (pmcoutput_alg3.parameters[:,:,pmcoutput_alg3.niterations], pmcoutput_alg4.parameters[:,:,pmcoutput_alg4.niterations], pmcoutput_alg5.parameters[:,:,pmcoutput_alg5.niterations]);
weightsABC = (pmcoutput_alg3.weights[:,pmcoutput_alg3.niterations], pmcoutput_alg4.weights[:,pmcoutput_alg4.niterations], pmcoutput_alg5.weights[:,pmcoutput_alg5.niterations]);
for i in 1:3 ##Loop over algorithms
    ww = weightsABC[i]
    for j in 1:3 ##Loop over parameters
        pp = samplesABC[i][j,:]
        PyPlot.subplot(130+j)
        PyPlot.plt[:hist](pp, weights=ww, normed=true, alpha=0.5)
    end
end

##Posterior summaries
outputs = (pmcoutput_alg3, pmcoutput_alg4, pmcoutput_alg5);
mapreduce(p -> parameter_means(p)[:,p.niterations], hcat, outputs)' ##Posterior means
sqrt.(mapreduce(p -> parameter_vars(p)[:,p.niterations], hcat, outputs))' ##Posterior sds

##Investigate best choice of alpha (as requested by reviewers)
function get_mse(pars::Array{Float64, 2}, w::Array{Float64, 1})
    mse = 0.0
    for i in 1:size(pars)[2]
        mse += sum((pars[:, i] .- log(theta0)).^2) * w[i]
    end
    mse
end

alphas = 0.05:0.1:0.95;
MSEs_alg5 = zeros(alphas);
MSEs_alg4 = zeros(alphas);
MSEs_alg3 = zeros(alphas);
abcinput.abcdist = WeightedEuclidean(x0);
srand(1);
for i in 1:length(alphas)
    x = abcPMC5(abcinput, 200, alphas[i], 50000)
    nits = x.niterations
    MSEs_alg5[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
    x = abcPMC3(abcinput, 200, alphas[i], 50000)
    nits = x.niterations
    MSEs_alg3[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
    x = abcPMC4(abcinput, 200, alphas[i], 50000)
    nits = x.niterations
    MSEs_alg4[i] = (nits == 0) ? Inf : get_mse(x.parameters[:,:,nits], x.weights[:,nits])
end

PyPlot.figure();
plot(alphas, log10.(MSEs_alg5), "r-x");
plot(alphas, log10.(MSEs_alg3), "b-o");
plot(alphas, log10.(MSEs_alg4), "g-*");
##Best alpha value around 0.5


text:
tests = ["test-bootsample-non-parametric",
         "test-bootsample-parametric",
         "test-utils",
         "test-dependencies"]

for t in tests
    test_file = "$(t).jl"
    include(test_file)
end