text:
using CCDReduction: axes_min_length,
                    fits_indices,
                    convert_value

@testset "bias subtraction" begin
    # testing non-mutating version
    @test @inferred(subtract_bias(ones(500, 500), ones(500, 500))) == zeros(500, 500)

    # testing mutating version
    frame = ones(500, 500)
    bias = ones(500, 500)
    @inferred subtract_bias!(frame, bias)
    @test frame == zeros(500, 500)

    # testing error
    @test_throws DimensionMismatch subtract_bias(ones(500, 1), ones(500, 500))
end

@testset "overscan subtraction" begin
    # testing non-mutating version
    @test @inferred(subtract_overscan(ones(500, 600), (1:2, :))) == zeros(500, 600)
    @test @inferred(subtract_overscan(ones(500, 600), "[:, 1:2]")) == zeros(500, 600)

    # testing mutating version
    frame = ones(500, 600)
    frame1 = ones(500, 600)
    @inferred subtract_overscan!(frame, (:, 540:600))
    @inferred subtract_overscan!(frame1, "[540:600, :]")
    @test frame == zeros(500, 600)
    @test frame == frame1

    # testing error
    @test_throws BoundsError subtract_overscan(ones(500, 600), (500:600, :))
end

@testset "flat correction" begin
    #testing non-mutating version
    frame = ones(5, 5)
    flat = fill(2.0, (5, 5))
    @test @inferred(flat_correct(frame, flat)) == ones(5, 5)
    @test @inferred(flat_correct(frame, flat, norm_value = 4.0)) == fill(2.0, (5, 5))

    #testing mutating version
    frame = ones(5, 5)
    flat = rand(5, 5)
    reduced_flat = flat ./ mean(flat)
    reduced_frame = frame ./ reduced_flat

    @inferred flat_correct!(frame, flat)
    @test reduced_frame == frame

    #testing error
    @test_throws DimensionMismatch flat_correct(ones(5, 5), ones(5, 6))
    @test_throws ErrorException flat_correct(ones(5, 5), ones(5, 5), norm_value = -2)
end

@testset "trim" begin
    @test trim(reshape(1:25, 5, 5), (:, 4:5)) == trimview(reshape(1:25, 5, 5), (:, 4:5))
    @test trim(reshape(1:25, 5, 5), (4:5, :)) == trimview(reshape(1:25, 5, 5), (4:5, :))
    @test trim(reshape(1:25, 5, 5), (:, 4:5)) == [1:5 6:10 11:15]
    @test trim(reshape(1:25, 5, 5), (1:2, :)) == [3:5 8:10 13:15 18:20 23:25]
    @test trim(reshape(1:25, 5, 5), "[4:5, :]") == trimview(reshape(1:25, 5, 5), (:, 4:5))
    @test trim(reshape(1:25, 5, 5), "[:, 4:5]") == trimview(reshape(1:25, 5, 5), (4:5, :))
    @test trim(reshape(1:25, 5, 5), "[1:5, 4:5]") == trimview(reshape(1:25, 5, 5), (4:5, :))
    @test trim(reshape(1:25, 5, 5), "[4:5, 1:5]") == trimview(reshape(1:25, 5, 5), (:, 4:5))

    # testing output types
    @test trimview(ones(5, 5), (:, 3:5)) isa SubArray
    @test trim(ones(5, 5), (:, 3:5)) isa Array

    # testing errors
    @test_throws ErrorException trim(ones(5, 5), (4:5, 1:4))
    @test_throws ErrorException trim(ones(5, 5), (:, :))
    @test_throws ErrorException trim(ones(5, 5), (4:6, :))
end

@testset "cropping" begin
    # testing on odd dimension array
    x = reshape(1:25, (5, 5))

    @test crop(x, (3, 3)) == [7:9 12:14 17:19]
    @test_logs (:warn, "dimension 1 changed from 4 to 5") @test crop(x, (4, 3)) == [6:10 11:15 16:20]
    @test crop(x, (4, 3), force_equal = false) == [6:9 11:14 16:19]
    @test crop(x, (1, 1)) == fill(13, (1, 1))
    @test_logs (:warn, "dimension 2 changed from 4 to 5") @test crop(x, (3, 4)) == [2:4 7:9 12:14 17:19 22:24]
    @test crop(x, (3, 4), force_equal = false) == [2:4 7:9 12:14 17:19]

    # testing on even dimension array
    x = reshape(1:16, (4, 4))

    @test cropview(x, (2, 2))  == [6:7 10:11]
    @test_logs (:warn, "dimension 1 changed from 3 to 4") @test cropview(x, (3, 2)) == [5:8 9:12]
    @test cropview(x, (3, 2), force_equal = false) == [5:7 9:11]
    @test cropview(x, (1, 1), force_equal = false) == fill(6, (1, 1))
    @test_logs (:warn, "dimension 2 changed from 3 to 4") @test cropview(x, (2, 3)) == [2:3 6:7 10:11 14:15]
    @test cropview(x, (2, 3), force_equal = false) == [2:3 6:7 10:11]

    # testing with colon
    x = reshape(1:25, (5, 5))

    @test cropview(x, (:, :)) == x
    @test cropview(x, (:, 3)) == [6:10 11:15 16:20]
    @test cropview(x, (3, :)) == [2:4 7:9 12:14 17:19 22:24]
    @test_logs (:warn, "dimension 1 changed from 2 to 3") @test cropview(x, (2, :)) == [2:4 7:9 12:14 17:19 22:24]
    @test cropview(x, (2, :), force_equal = false) == [2:3 7:8 12:13 17:18 22:23]
    @test_logs (:warn, "dimension 2 changed from 2 to 3") @test cropview(x, (:, 2)) == [6:10 11:15 16:20]
    @test cropview(x, (:, 2), force_equal = false) == [6:10 11:15]

    # testing output types
    @test cropview(ones(5, 5), (3, 3)) isa SubArray
    @test crop(ones(5, 5), (3, 3)) isa Array

    # testing error
    @test_throws BoundsError cropview(ones(5, 5), (7, 3))
    @test_throws ErrorException cropview(ones(5, 5), (3, -5))
    @test_throws ErrorException cropview(ones(5, 5), (3, 4, 5))
end

@testset "combining" begin
    x = [reshape(1:4, (2,2)) for i = 1:10]

    @test combine(x) == [1:2 3:4]
    @test combine(x, method = sum) == [1:2 3:4] .* 10

    # testing error
    @test_throws DimensionMismatch combine(rand(5, 5), rand(6, 6))
    @test_throws MethodError combine()
end

@testset "dark subtraction" begin
    # testing non-mutating version
    @test @inferred(subtract_dark(ones(5, 5), ones(5, 5))) == zeros(5, 5)
    @test @inferred(subtract_dark(ones(5, 5), ones(5, 5), data_exposure = 1, dark_exposure = 4)) == fill(0.75, 5, 5)
    @test @inferred(subtract_dark(ones(Float32, 5, 5), ones(Float32, 5, 5), data_exposure = 13, dark_exposure = 17)) ≈ fill(Float32(4/17), 5, 5)

    # testing mutating version
    frame = ones(5, 5)
    dark_frame = ones(5, 5)
    @inferred(subtract_dark!(frame, dark_frame))
    @test frame == zeros(5, 5)

    frame = fill(5.0, 5, 5)
    dark_frame = ones(5, 5)
    @inferred(subtract_dark!(frame, dark_frame, data_exposure = 2, dark_exposure = 1))
    @test frame  == fill(3.0, 5, 5)

    frame = fill(Float32(3), 5, 5)
    dark_frame = fill(Float32(2), 5, 5)
    @inferred(subtract_dark!(frame, dark_frame, data_exposure = 11, dark_exposure = 19))
    @test frame ≈ fill(Float32(35/19), 5, 5)

    # testing error
    @test_throws DimensionMismatch subtract_dark!(ones(5,5), ones(6,6))
end

@testset "helper" begin
    # testing axes_min_length
    @test axes_min_length((:, :)) == 1
    @test axes_min_length((:, 5:10)) == 2
    @test axes_min_length((1:2, :)) == 1
    @test axes_min_length((1:2, 5:12)) == 1
    @test axes_min_length((10:100, 1:5)) == 2

    # testing fits_indices
    @test fits_indices("[1024:2048, 200:300]") == [200:300, 1024:2048]
    @test fits_indices("[:, 200:300]") == [200:300, :]
    @test fits_indices("[1024:2048, :]") == [:, 1024:2048]
    @test fits_indices("[200:300, 1024:2048]") == [1024:2048, 200:300]

    # testing convert_value
    @test convert_value(Int16, 5.4) == 5
    @test convert_value(Int32, 5.4) == 5
    @test convert_value(Int64, -5.4) == -5
    @test convert_value(Float64, -5.4) ≈ -5.4
    @test convert_value(Float32, -5.4) ≈ -5.4
    @test convert_value(Float32, -5.4) ≈ -5.4
end


text:
function read_legacy(io::IO, ver::Float32)
  Z = getfield(BUF, :buf)
  L = getfield(BUF, :int64_buf)

  # read begins ------------------------------------------------------
  N     = fastread(io, Int64)
  checkbuf_strict!(L, 2*N)
  fast_readbytes!(io, Z, 3*N)
  c1    = copy(Z[1:N])
  c2    = copy(Z[N+1:2*N])
  y     = code2typ.(getindex(Z, 2*N+1:3*N))
  cmp   = fastread(io) % Bool
  read!(io, L)
  nx    = getindex(L, N+1:2*N)
  cmp && checkbuf_8!(Z, maximum(nx))

  ver < 0.5f0 && error("No legacy support for SeisIO file format version < 0.5")

  # Get file creation time
  fname = split(io.name)[2][1:end-1]
  isfile(fname) || error(string("Can't stat file ", fname))
  st = stat(fname)
  t0 = st.ctime
  if ver == 0.5f0
    S = SeisData(N)
    setfield!(S, :id, read_string_vec(io, Z))
    setfield!(S, :name, read_string_vec(io, Z))
    setfield!(S, :loc, InstrumentPosition[read(io, code2loctyp(getindex(c1, i))) for i = 1:N])
    read!(io, S.fs)
    read!(io, S.gain)

    # Here is the hard part. For older files, we read a degenerate InstResp
    R = InstrumentResponse[]
    for i = 1:N
      T = code2resptyp(getindex(c2, i))
      (T == GenResp) && (push!(R, read(io, GenResp)); continue)
      et = Complex{T == PZResp ? Float32 : Float64}

      # skip :c, it was never used
      fastskip(io, T == PZResp ? 4 : 8)

      # read poles
      np = fastread(io, Int64)
      p = zeros(et, np)
      read!(io, p)

      # read zeros
      nz = fastread(io, Int64)
      z = zeros(et, nz)
      read!(io, z)

      # push to R
      push!(R, T(p = p, z = z))
    end
    S.resp .= R
    setfield!(S, :units, read_string_vec(io, Z))
    setfield!(S, :src, read_string_vec(io, Z))
    setfield!(S, :misc, [read_misc(io, Z) for i = 1:N])
    setfield!(S, :notes, [read_string_vec(io, Z) for i = 1:N])
    setfield!(S, :t, [read!(io, Array{Int64, 2}(undef, getindex(L, i), 2)) for i = 1:N])
    setfield!(S, :x,
    FloatArray[cmp ?
    (fast_readbytes!(io, Z, getindex(nx, i)); Blosc.decompress(getindex(y,i), Z)) :
    read!(io, Array{getindex(y,i), 1}(undef, getindex(nx, i))) for i = 1:N])
  elseif ver < 0.53
    #=  Read process for 0.51 and 0.52 is identical; 0.52 added two Types <:
        InstrumentResponse that v0.51 doesn't know about, so this is a safe
        procedure =#
    S = SeisData(N)
    setfield!(S, :id, read_string_vec(io, Z))
    setfield!(S, :name, read_string_vec(io, Z))
    setfield!(S, :loc, InstrumentPosition[read(io, code2loctyp(getindex(c1, i))) for i = 1:N])
    read!(io, S.fs)
    read!(io, S.gain)

    # Here is the hard part. For older files, we read a degenerate InstResp
    R = InstrumentResponse[]
    for i = 1:N
      T = code2resptyp(getindex(c2, i))
      if T in (PZResp, PZResp64, GenResp)
        push!(R, read(io, T))
      elseif T == MultiStageResp
        K = fastread(io, Int64)
        codes = fastread(io, K)
        M = MultiStageResp(K)
        A = Array{RespStage,1}(undef, 0)
        for j = 1:K
          c = codes[j]
          if c == 0x03
            units_out = String(fastread(io, fastread(io, Int64)))
            units_in = String(fastread(io, fastread(io, Int64)))
            M.i[j] = units_in
            M.o[j] = units_out
            if j == 2
              M.i[1] = units_out
            end
            CR = CoeffResp(
                       read!(io, Array{Float64,1}(undef, fastread(io, Int64))),
                       read!(io, Array{Float64,1}(undef, fastread(io, Int64)))
                       )
            push!(A, CR)
          else
            push!(A, read(io, code2resptyp(c)))
          end
        end
        M.stage .= A
        read!(io, M.fs)
        read!(io, M.gain)
        read!(io, M.fg)
        read!(io, M.delay)
        read!(io, M.corr)
        read!(io, M.fac)
        read!(io, M.os)
        push!(R, M)
      end
    end
    S.resp .= R
    setfield!(S, :units, read_string_vec(io, Z))
    for i = 1:N
      if typeof(S.resp[i]) == MultiStageResp
        K = length(S.resp[i].fs)
        if K > 0
          S.resp[i].o[1] = S.units[i]
        end
      end
    end
    setfield!(S, :src, read_string_vec(io, Z))
    setfield!(S, :misc, [read_misc(io, Z) for i = 1:N])
    setfield!(S, :notes, [read_string_vec(io, Z) for i = 1:N])
    setfield!(S, :t, [read!(io, Array{Int64, 2}(undef, getindex(L, i), 2)) for i = 1:N])
    setfield!(S, :x,
    FloatArray[cmp ?
    (fast_readbytes!(io, Z, getindex(nx, i)); Blosc.decompress(getindex(y,i), Z)) :
    read!(io, Array{getindex(y,i), 1}(undef, getindex(nx, i))) for i = 1:N])
  end
  return S
end

ext: jl


text:
module ImageSegmentationEvaluation
export
    evaluate,

    BoundaryGradient,
    BoundaryShift,
    BoundaryMooreTracing
# package code goes here
include("utils.jl")
include("supervised.jl")

export unsupervised_metrics
include("unsupervised.jl")

end # module


text:
module SimpleDist

using ..Structs: Query

export dist

function dist(x::Query, y::Query)
    return (x.data - y.data) ^ 2
end

module MidDist
    export mid_dist

    function mid_dist(x, y)
        return ((x.data - y.data) ^ 2) / 2
    end
end

end


text:
push!(LOAD_PATH, "Modules")
using LinearAlgebra
using Colors
using JuliRay

Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function Plane2Sphere(p::Array{T,1}) where T<:Real
    if length(p)≠2
        error("No Point on ℝ²")
    end
    return [2p[1]/(1+p[1]^2+p[2]^2),2p[2]/(1+p[1]^2+p[2]^2),(-1+p[1]^2+p[2]^2)/(1+p[1]^2+p[2]^2)]
end

function Plane2Sphere(θ::Real,p::Array{T,1}) where T<:Real
    if length(p)≠2
        error("No Point on ℝ²")
    end
    c=[0,0,tan(θ)]
    r=1/cos(θ)
    P=p*cot(θ/2+π/4)
    return Plane2Sphere(P)*r+c
end

function POLYGON(vertices;r1=1.0,r2=0.05,r3=0.025)
    n=length(vertices)
    V=rgbColor(csgUnion((p->Sphere(p,r2)).(vertices)),RGB(0.1,0.1,0.1))
    E=rgbColor(csgUnion([Cylinder(vertices[i],vertices[mod(i,n)+1],r3) for i in 1:n]),RGB(0.1,0.1,0.1))
    F=rgbftColor(Polygon(vertices),RGB(0.5,0.5,1),FT(0.1,0.1))
    return csgUnion(V,E,F)
end

function Mirror(q,p1,p2,p3)
    e₃=JuliRay.NormalVector(p1,p2,p3)
    e₁=JuliRay.OrthogonalVector(e₃)
    e₂=cross(e₃,e₁)
    A=hcat(e₁,e₂,e₃)*transpose(hcat(e₁,e₂,-e₃))
    fixedpoint=(p1+p2+p3)/3
    b=fixedpoint-A*fixedpoint
    return A*q+b
end

function Put(F,i,O)
    p1=O
    n=length(F)
    p2,p3=F[mod(i-1,n)+1],F[mod(i,n)+1]
    return (p->Mirror(p,p1,p2,p3)).(F)
end

function f(x)
    ifelse(x>0,exp(-1/x),0)
end
function g(x)
    f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a,b,x)
    g((x-a)/(b-a))
end

## Cube

h=1/√3 # 単位球面に内接させた際の内接球半径
n=4 # n角形
m=6 # m面体
M=80 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径

for i ∈ 1:M
    t=3(i-1)/M
    θ=(1-Smooth(0.5,1.4,t)+Smooth(1.6,2.5,t))*π/2 # パラメータθは0 → π/2 → 0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]
    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    NET(F)=Scaling(Rotate(csgUnion([POLYGON(F[i],r2=0.02,r3=0.02) for i ∈ 1:m]),[0,0,1],45°),1.5/√3)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],2,O)
    F[4]=Put(F[1],3,O)
    F[5]=Put(F[1],4,O)
    F[6]=Put(F[4],1,O)
    object01=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[2],2,O)
    F[5]=Put(F[2],4,O)
    F[6]=Put(F[3],1,O)
    object02=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[1],2,O)
    F[5]=Put(F[2],4,O)
    F[6]=Put(F[3],1,O)
    object03=NET(F)

    F[2]=Put(F[1],4,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],2,O)
    F[5]=Put(F[2],1,O)
    F[6]=Put(F[4],1,O)
    object04=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],2,O)
    F[5]=Put(F[1],4,O)
    F[6]=Put(F[3],1,O)
    object05=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],1,O)
    F[5]=Put(F[2],2,O)
    F[6]=Put(F[4],4,O)
    object06=NET(F)

    F[2]=Put(F[1],4,O)
    F[3]=Put(F[1],2,O)
    F[4]=Put(F[1],3,O)
    F[5]=Put(F[2],1,O)
    F[6]=Put(F[4],1,O)
    object07=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],1,O)
    F[5]=Put(F[2],4,O)
    F[6]=Put(F[3],2,O)
    object08=NET(F)

    F[2]=Put(F[1],4,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],1,O)
    F[5]=Put(F[2],1,O)
    F[6]=Put(F[3],2,O)
    object09=NET(F)

    F[2]=Put(F[1],4,O)
    F[3]=Put(F[1],3,O)
    F[4]=Put(F[3],1,O)
    F[5]=Put(F[2],1,O)
    F[6]=Put(F[4],2,O)
    object10=NET(F)

    F[2]=Put(F[1],1,O)
    F[3]=Put(F[2],3,O)
    F[4]=Put(F[1],2,O)
    F[5]=Put(F[4],3,O)
    F[6]=Put(F[5],1,O)
    object11=NET(F)

    object=csgUnion(
        ParallelTranslation(object01,[-8,4,0]),
        ParallelTranslation(object02,[-3,4,0]),
        ParallelTranslation(object03,[2,4,0]),
        ParallelTranslation(object04,[7,4,0]),
        ParallelTranslation(object05,[-8,0,0]),
        ParallelTranslation(object06,[-3,0,0]),
        ParallelTranslation(object07,[2,0,0]),
        ParallelTranslation(object08,[7,0,0]),
        ParallelTranslation(object09,[-8,-4,0]),
        ParallelTranslation(object10,[-3,-4,0]),
        ParallelTranslation(object11,[5,-3.5,0])
    )
    lng = -90°-15°*(Smooth(0.1,0.5,t)-Smooth(2.5,2.9,t))
    lat = 90°-50°*(Smooth(0.1,0.5,t)-Smooth(2.5,2.9,t))
    camera = LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=0.065,width=640*2,height=360*2)
    render(object, camera=camera, name="NetOfCubes", index=i)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
module HW4

using CommonRLInterface
using StaticArrays
using Compose

import ColorSchemes

export gw

const RL = CommonRLInterface

mutable struct GridWorld <: AbstractEnv
    size::SVector{2, Int}
    rewards::Dict{SVector{2, Int}, Float64}
    state::SVector{2, Int}
end

function GridWorld()
    rewards = Dict(SA[7,7]=> 10.0,
                   SA[3,1]=> -2.0,
                   SA[4,3]=>-10.0,
                   SA[7,6]=> -5.0)
    return GridWorld(SA[10, 10], rewards, SA[1,1])
end

RL.reset!(env::GridWorld) = (env.state = SA[1,1])
RL.actions(env::GridWorld) = (SA[1,0], SA[-1,0], SA[0,1], SA[0,-1])
RL.observe(env::GridWorld) = env.state
RL.terminated(env::GridWorld) = haskey(env.rewards, env.state)

function RL.act!(env::GridWorld, a)
    if rand() < 0.4 # 40% chance of going in a random direction (=30% chance of going in a wrong direction)
        a = rand(actions(env))
    end

    env.state = clamp.(env.state + a, SA[1,1], env.size)

    return get(env.rewards, env.state, -0.1)
end

@provide RL.observations(env::GridWorld) = [SA[x, y] for x in 1:env.size[1], y in 1:env.size[2]]
@provide RL.clone(env::GridWorld) = GridWorld(env.size, copy(env.rewards), env.state)
@provide RL.state(env::GridWorld) = env.state
@provide RL.setstate!(env::GridWorld, s) = (env.state = s)

function render(env::GridWorld; color::Function=s->get(env.rewards, s, -0.1), policy::Union{Function,Nothing}=nothing)
    nx, ny = env.size
    cells = []
    for s in observations(env)
        r = get(env.rewards, s, 0.0)
        clr = get(ColorSchemes.redgreensplit, (r+10.0)/20.0)
        cell = context((s[1]-1)/nx, (ny-s[2])/ny, 1/nx, 1/ny)
        if policy !== nothing
            a = policy(s)
            txt = compose(context(), text(0.5, 0.5, aarrow[a], hcenter, vcenter), stroke("black"))
            compose!(cell, txt)
        end
        clr = tocolor(color(s))
        compose!(cell, rectangle(), fill(clr), stroke("gray"))
        push!(cells, cell)
    end
    grid = compose(context(), linewidth(0.5mm), cells...)
    outline = compose(context(), linewidth(1mm), rectangle(), stroke("gray"))

    s = env.state
    agent_ctx = context((s[1]-1)/nx, (ny-s[2])/ny, 1/nx, 1/ny)
    agent = compose(agent_ctx, circle(0.5, 0.5, 0.4), fill("orange"))

    sz = min(w,h)
    return compose(context((w-sz)/2, (h-sz)/2, sz, sz), agent, grid, outline)
end

@provide RL.render(env::GridWorld) = render(env)

tocolor(x) = x
function tocolor(r::Float64)
    minr = -10.0
    maxr = 10.0
    frac = (r-minr)/(maxr-minr)
    return get(ColorSchemes.redgreensplit, frac)
end

const aarrow = Dict(SA[0,1]=>'↑', SA[-1,0]=>'←', SA[0,-1]=>'↓', SA[1,0]=>'→')

gw = GridWorld()

end


text:
export GameRecorder, push!, trim!, setfinalstate!

# used to record what happened in one full game
mutable struct GameRecorder
    states::Vector{Bitboard}
    moves::Vector{Dirs}
    counter::Int
    final_state::Bitboard
    GameRecorder(init_len) = begin
        dummy_boards = Vector{Bitboard}(undef, init_len)
        new(dummy_boards, Vector{Dirs}(init_len), 0, dummy_boards[1])
    end
    GameRecorder() = GameRecorder(2888)
end

import Base

function Base.push!(gr::GameRecorder, board, move)
    # increase the size by 20% if already at limit
    if gr.counter == length(gr.states)
        resize!(gr.states, round(Int, gr.counter*1.2))
        resize!(gr.moves, round(Int, gr.counter*1.2))
    end

    gr.counter += 1
    gr.states[gr.counter] = board
    gr.moves[gr.counter] = move

    gr
end

function trim!(gr::GameRecorder)
    resize!(gr.states, gr.counter)
    resize!(gr.moves, gr.counter)
    gr
end

function setfinalstate!(gr::GameRecorder, board)
    gr.final_state = board
end

Base.length(gr::GameRecorder) = gr.counter


text:
using NPZ
using CodecZlib

"""
Represents an iterable dataset of type AtariDataSet with the following fields:

`dataset`: Dict{Symbol, Any}, representation of the dataset as a Dictionary with style as `style`
`epochs`: Vector{Int}, list of epochs to load
`repo`: String, the repository from which the dataset is taken
`length`: Integer, the length of the dataset
`batch_size`: Integer, the size of the batches returned by `iterate`.
`style`: Tuple, the type of the NamedTuple, for now SARTS and SART is supported.
`rng`<: AbstractRNG.
`meta`: Dict, the metadata provided along with the dataset
`is_shuffle`: Bool, determines if the batches returned by `iterate` are shuffled.
"""
struct AtariDataSet{T<:AbstractRNG} <:RLDataSet
    dataset::Dict{Symbol, Any}
    epochs::Vector{Int}
    repo::String
    length::Integer
    batch_size::Integer
    style::Tuple
    rng::T
    meta::Dict
    is_shuffle::Bool
end

const samples_per_epoch = Int(1e6)
const atari_frame_size = 84
const epochs_per_game = 50

"""
    dataset(dataset::String, epochs::Vector{Int}; repo::String, style::Tuple, rng<:AbstractRNG, is_shuffle::Bool, max_iters::Int64, batch_size::Int64)

Creates a dataset of enclosed in a AtariDataSet type and other related metadata for the `dataset` that is passed.
The `AtariDataSet` type is an iterable that fetches batches when used in a for loop for convenience during offline training.

`dataset`: String, name of the datset.
`index`: Int, analogous to v
`epochs`: Vector{Int}, list of epochs to load
`repo`: Name of the repository of the dataset
`style`: the style of the iterator and the Dict inside AtariDataSet that is returned.
`rng`: StableRNG
`max_iters`: maximum number of iterations for the iterator.
`is_shuffle`: whether the dataset is shuffled or not. `true` by default.
`batch_size`: batch_size that is yielded by the iterator. Defaults to 256.

The returned type is an infinite iterator which can be called using `iterate` and will return batches as specified in the dataset.
"""
function dataset(
    game::String,
    index::Int,
    epochs::Vector{Int};
    style=SARTS,
    repo = "atari-replay-datasets",
    rng = StableRNG(123),
    is_shuffle = true,
    batch_size=256
)

    try
        @datadep_str "$repo-$game-$index"
    catch
        throw("The provided dataset is not available")
    end

    path = @datadep_str "$repo-$game-$index"

    @assert length(readdir(path)) == 1
    folder_name = readdir(path)[1]

    folder_path = "$path/$folder_name"
    files = readdir(folder_path)
    file_prefixes = collect(Set(map(x->join(split(x,"_")[1:2], "_"), files)))
    fields = map(collect(file_prefixes)) do x
        if split(x, "_")[1] == "\$store\$"
            x = split(x, "_")[2]
        else
            x = x
        end
    end

    s_epochs = Set(epochs)

    dataset = Dict()

    for (prefix, field) in zip(file_prefixes, fields)
        for epoch in s_epochs
            @assert epoch <= epochs_per_game
            data = open("$folder_path/$(prefix)_ckpt.$epoch.gz") do file
                stream = GzipDecompressorStream(file)
                NPZ.npzreadarray(stream)
            end

            if field == "observation"
                data = permutedims(data, [3, 2, 1]) # not sure about the orientation of the frame
            end

            if haskey(dataset, field)
                if field == "observation"
                    dataset[field] = cat(dataset[field], data, dims=3)
                else
                    dataset[field] = cat(dataset[field], data, dims=1)
                end
            else
                dataset[field] = data
            end
        end
    end

    num_epochs = length(s_epochs)

    atari_verify(dataset, num_epochs)

    N_samples = size(dataset["observation"])[3]

    final_dataset = Dict{Symbol, Any}()
    meta = Dict{String, Any}()

    for (key, d_key) in zip(["observation", "action", "reward", "terminal"], Symbol.(["state", "action", "reward", "terminal"]))
            final_dataset[d_key] = dataset[key]
    end

    for key in keys(dataset)
        if !(key in ["observation", "action", "reward", "terminal"])
            meta[key] = dataset[key]
        end
    end

    return AtariDataSet(final_dataset, epochs, repo, N_samples, batch_size, style, rng, meta, is_shuffle)

end

function iterate(ds::AtariDataSet, state = 0)
    rng = ds.rng
    batch_size = ds.batch_size
    length = ds.length
    is_shuffle = ds.is_shuffle
    style = ds.style

    if is_shuffle
        inds = rand(rng, 1:length-1, batch_size)
    else
        if (state+1) * batch_size <= length
            inds = state*batch_size+1:(state+1)*batch_size
        else
            return nothing
        end
        state += 1
    end

    batch = (state = view(ds.dataset[:state], :, :, inds),
    action = view(ds.dataset[:action], inds),
    reward = view(ds.dataset[:reward], inds),
    terminal = view(ds.dataset[:terminal], inds))

    if style == SARTS
        batch = merge(batch, (next_state = view(ds.dataset[:state], :, :, (1).+(inds)),))
    end

    return batch, state
end

take(ds::AtariDataSet, n::Integer) = take(ds.dataset, n)
length(ds::AtariDataSet) = ds.length
IteratorEltype(::Type{AtariDataSet}) = EltypeUnknown() # see if eltype can be known (not sure about carla and adroit)

function atari_verify(dataset::Dict, num_epochs::Int)
    @assert size(dataset["observation"]) == (atari_frame_size, atari_frame_size, num_epochs*samples_per_epoch)
    @assert size(dataset["action"]) == (num_epochs * samples_per_epoch,)
    @assert size(dataset["reward"]) == (num_epochs * samples_per_epoch,)
    @assert size(dataset["terminal"]) == (num_epochs * samples_per_epoch,)
end


text:
α = 0.025

@testset "nested archimedean copulas helpers" begin
  Random.seed!(43)
  u = nestedcopulag("clayton", [[1,2],[3,4]], [2., 3.], 1.1, [0.1 0.2 0.3 0.4 0.5; 0.2 0.3 0.4 0.5 0.6])
  @test u ≈ [0.153282 0.182421 0.636606 0.679396; 0.381051 0.430175 0.254842 0.279192] atol=1.0e-5
  Random.seed!(43)
  n = nestedstep("clayton", [0.2 0.8; 0.1 0.7], [0.2, 0.4], 2., 1.5)
  @test n ≈ [0.0504023 0.545041; 0.0736747 0.58235] atol=1.0e-5
  @test_throws AssertionError nestedstep("clayto", [0.2 0.8; 0.1 0.7], [0.2, 0.4], 2., 1.5)
end

@testset "nested archimedean copulas exceptions" begin
  nestedarchcopulagen(100000, [2, 2], [2., 2.], 0.5, "frank")
  @test_throws DomainError testnestedθϕ([2, 2], [2.1, 2.2], 0.5, "gumbel")
  @test_throws DomainError testnestedθϕ([2, 2], [2.1, 2.2], 3.5, "gumbel")
  @test_throws DomainError testnestedθϕ([2, 2], [0.8, 1.1], 0.5, "amh")
  @test_throws AssertionError testnestedθϕ([2, 2], [0.8], 0.5, "amh")
  @test_throws AssertionError nestedarchcopulagen(100000, [2, 2], [2., 2.], 0.5, "fran")
  @test_throws DomainError nestedarchcopulagen(500000, [2.2, 3.6, 1.1], "gumbel")
  @test_throws DomainError nestedarchcopulagen(500000, [4.2, 3.6, 0.1], "gumbel")
end

@testset "nested gumbel copula" begin
  @testset "single nested" begin
    Random.seed!(44)
    x = nestedarchcopulagen(500000, [2,2], [4.2, 6.1], 2.1, "gumbel", 1)
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
    M = [1. 0.7619 0.52380 0.52380 0.52380; 0.7619 1. 0.52380 0.52380 0.52380; 0.52380 0.52380 1. 0.83606 0.52380]
    @test corkendall(x)[1:3,:] ≈ M atol=1.0e-2
    @test tail(x[:,1], x[:,2], "r", 0.01) ≈ 2-2^(1/4.2) atol=1.0e-1
    @test tail(x[:,2], x[:,3], "r", 0.01) ≈ 2-2^(1/2.1) atol=1.0e-1
    @test tail(x[:,1], x[:,5], "r", 0.01) ≈ 2-2^(1/2.1) atol=1.0e-1
    @test tail(x[:,3], x[:,4], "r", 0.01) ≈ 2-2^(1/6.1) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0
    @test tail(x[:,1], x[:,3], "l", 0.00001) ≈ 0
  end
  @testset "double nested" begin
    Random.seed!(43)
    x = nestedarchcopulagen(200000, [[2,2], [2,2]], [[4.1, 3.8],[5.1, 6.1]], [1.9, 2.4], 1.2, "gumbel")
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
    @test corkendall(x)[1,:] ≈ [1., 0.7560, 0.47368, 0.47368, 1/6, 1/6, 1/6, 1/6] atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r", 0.01) ≈ 2-2^(1/4.1) atol=1.0e-1
    @test tail(x[:,1], x[:,3], "r", 0.01) ≈ 2-2^(1/1.9) atol=1.0e-1
    @test tail(x[:,1], x[:,5], "r", 0.01) ≈ 2-2^(1/1.2) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0
    @test tail(x[:,1], x[:,3], "l", 0.00001) ≈ 0
  end
  @testset "hierarchical" begin
    Random.seed!(42)
    x = nestedarchcopulagen(500000, [4.2, 3.6, 1.1], "gumbel")
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
    @test corkendall(x)[1:2,:] ≈ [1. 0.7619 0.72222 0.0909; 0.7619 1. 0.72222 0.0909] atol=1.0e-2
    @test tail(x[:,2], x[:,3], "r", 0.01) ≈ 2-2^(1/3.6) atol=1.0e-1
    @test tail(x[:,3], x[:,4], "r", 0.01) ≈ 2-2^(1/1.1) atol=1.0e-2
    @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0
  end
end

@testset "nested Ali-Mikhail-Haq copula" begin
  Random.seed!(43)
  x = nestedarchcopulagen(200000, [3, 2], [0.8, 0.7], 0.5, "amh", 2)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  c = corkendall(x)
  @test c[1:3,1] ≈ [1., 0.23373, 0.23373] atol=1.0e-2
  @test c[3:6,4] ≈ [0.1288, 1., 0.19505, 0.1288] atol=1.0e-2
  @test tail(x[:,4], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,4], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "l", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
end

@testset "nested Frank copula" begin
  Random.seed!(44)
  x = nestedarchcopulagen(250000, [3, 2],  [8., 10.], 2., "frank", 2)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  @test corkendall(x)[1:4,1] ≈ [1.0, 0.60262, 0.60262, 0.2139] atol=1.0e-2
  @test corkendall(x)[3:5,4] ≈ [0.2139, 1.0, 0.6658] atol=1.0e-2
  @test corkendall(x)[6:7,6] ≈ [1.0, 0.2139] atol=1.0e-2
  @test tail(x[:,4], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,4], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "l", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "r", 0.0001) ≈ 0
  @test tail(x[:,6], x[:,7], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
end

@testset "nested Clayton copula" begin
  Random.seed!(42)
  x = nestedarchcopulagen(100000, [2, 3],  [3., 4.], 1.5, "clayton", 2)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,7], Uniform(0,1))) > α
  cc = corkendall(x)
  @test cc[:,1] ≈ [1.0, 3/5, 1.5/3.5, 1.5/3.5, 1.5/3.5, 1.5/3.5, 1.5/3.5] atol=1.0e-2
  @test cc[:,5] ≈ [1.5/3.5, 1.5/3.5, 2/3, 2/3, 1.0, 1.5/3.5, 1.5/3.5] atol=1.0e-2
  @test cc[6,7] ≈ 1.5/3.5 atol=1.0e-2
  @test tail(x[:,4], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
  @test tail(x[:,1], x[:,5], "l", 0.01) ≈ 2^(-1/(1.5)) atol=1.0e-1
  @test tail(x[:,1], x[:,2], "l", 0.01) ≈ 2^(-1/3) atol=1.0e-1
  @test tail(x[:,4], x[:,5], "l", 0.01) ≈ 2^(-1/4) atol=1.0e-1
  @test tail(x[:,6], x[:,7], "l", 0.01) ≈ 2^(-1/(1.5)) atol=1.0e-1
end


text:
using Diffusion
import MPI
using CUDA

MPI.Init()
nx = ny                        = Csize_t(8)
nt                             = Csize_t(100)
lx = ly = lam = cp_min = amp_T = Cdouble(42.0)
Cp_d                           = pointer(CUDA.ones(Cdouble, nx, ny))
comm                           = MPI.COMM_WORLD.val

Diffusion.julia_diffusion(nx, ny, nt, lx, ly, lam, cp_min, amp_T, Cp_d, comm)

MPI.Finalize()

max_stars_repo_name: omlins/libdiffusion

max_issues_repo_name: omlins/libdiffusion

max_forks_repo_name: omlins/libdiffusion


text:
############################################################
# judiVector # #################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiVector, judiVectorException, subsample, judiVector_to_SeisBlock, src_to_SeisBlock
export time_resample, time_resample!, judiTimeInterpolation
export write_shot_record, get_data, convert_to_array, rebuild_jv

############################################################

# structure for seismic data as an abstract vector
mutable struct judiVector{vDT<:Number, AT} <: joAbstractLinearOperator{vDT,vDT}
    name::String
    m::Integer
    n::Integer
    nsrc::Integer
    geometry::Geometry
    data::Array{AT, 1}
end

mutable struct judiVectorException <: Exception
    msg :: String
end

############################################################

## outer constructors

"""
    judiVector
        name::String
        m::Integer
        n::Integer
        nsrc::Integer
        geometry::Geometry
        data

Abstract vector for seismic data. This vector-like structure contains the geometry and data for either\\
receiver data (shot records) or source data (wavelets).

Constructors
============

Construct vector from `Geometry` structure and cell array of shot records or wavelets. The `data` keyword\\
can also be a single (non-cell) array, in which case the data is the same for all source positions:

    judiVector(geometry, data)

Construct vector for observed data from `SegyIO.SeisBlock`. `segy_depth_key` is the `SegyIO` keyword \\
that contains the receiver depth coordinate:

    judiVector(SegyIO.SeisBlock; segy_depth_key="RecGroupElevation")

Construct vector for observed data from out-of-core data container of type `SegyIO.SeisCon`:

    judiVector(SegyIO.SeisCon; segy_depth_key="RecGroupElevation")

Examples
========

(1) Construct data vector from `Geometry` structure and a cell array of shot records:

    dobs = judiVector(rec_geometry, shot_records)

(2) Construct data vector for a seismic wavelet (can be either cell arrays of individual\\
wavelets or a single wavelet as an array):

    q = judiVector(src_geometry, wavelet)

(3) Construct data vector from `SegyIO.SeisBlock` object:

    using SegyIO
    seis_block = segy_read("test_file.segy")
    dobs = judiVector(seis_block; segy_depth_key="RecGroupElevation")

(4) Construct out-of-core data vector from `SegyIO.SeisCon` object (for large SEG-Y files):

    using SegyIO
    seis_container = segy_scan("/path/to/data/directory","filenames",["GroupX","GroupY","RecGroupElevation","SourceDepth","dt"])
    dobs = judiVector(seis_container; segy_depth_key="RecGroupElevation")

"""
function judiVector(geometry::Geometry, data::Array{T, N}) where {T, N}
    T == Float32 || (data = tof32(data))
    N < 3 || throw(judiVectorException("Only 1D (trace) and 2D (record) input data supported"))
    # length of vector
    n = 1
    nsrc = get_nsrc(geometry)
    m = n_samples(geometry, nsrc)

    dataCell = Array{Array{T, N}, 1}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = deepcopy(data)
    end
    return judiVector{T, Array{T, N}}("Seismic data vector", m, n, nsrc, geometry, dataCell)
end

# constructor if data is passed as a cell array
function judiVector(geometry::Geometry, data::Array{Array{T, N}, 1}) where {T, N}
    T == Float32 || (data = tof32.(data))

    # length of vector
    nsrc = get_nsrc(geometry)
    m = n_samples(geometry, nsrc)

    n = 1
    return judiVector{T, Array{T, N}}("Seismic data vector",m,n,nsrc,geometry,data)
end

############################################################
# constructors from SEGY files or out-of-core containers

# contructor for in-core data container
function judiVector(data::SegyIO.SeisBlock; segy_depth_key="RecGroupElevation")
    # length of data vector
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))

    numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
    numSamples = get_header(data,"ns")[end]
    m = numTraces*numSamples
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{Array{Float32, 2}, 1}(undef, nsrc)
    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        dataCell[j] = convert(Array{Float32, 2}, data.data[:,traces])
    end

    return judiVector{Float32, Array{Float32, 2}}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for in-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisBlock)

    # length of data vector
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))

    numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
    numSamples = get_header(data,"ns")[end]
    m = numTraces*numSamples
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{Array{Float32, 2}, 1}(undef, nsrc)
    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        dataCell[j] = convert(Array{Float32, 2}, data.data[:,traces])
    end

    return judiVector{Float32, Array{Float32, 2}}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from single container
function judiVector(data::SegyIO.SeisCon; segy_depth_key="RecGroupElevation")
    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
    end
    m = numTraces*data.ns
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = split(data,j)
    end

    return judiVector{Float32, SegyIO.SeisCon}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for single out-of-core data container and given geometry
function judiVector(geometry::Geometry, data::SegyIO.SeisCon; vDT::DataType=Float32)
    vDT == Float32 || throw(judiVectorException("Domain and range types not supported"))

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
    end
    m = numTraces*data.ns
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = split(data,j)
    end

    return judiVector{Float32, SegyIO.SeisCon}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers
function judiVector(data::Array{SegyIO.SeisCon,1}; segy_depth_key="RecGroupElevation")

    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4))
    end
    m = numTraces*data[1].ns    # SEGY requires same number of samples for every trace
    n = 1

    # extract geometry from data container
    geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = data[j]
    end

    return judiVector{Float32, SegyIO.SeisCon}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers and given geometry
function judiVector(geometry::Geometry, data::Array{SegyIO.SeisCon})
    # length of data vector
    nsrc = length(data)
    numTraces = 0
    for j=1:nsrc
        numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4))
    end
    m = numTraces*data[1].ns
    n = 1

    # fill data vector with pointers to data location
    dataCell = Array{SegyIO.SeisCon}(undef, nsrc)
    for j=1:nsrc
        dataCell[j] = data[j]
    end

    return judiVector{Float32, SegyIO.SeisCon}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

############################################################
## overloaded Base functions

# conj(jo)
conj(a::judiVector{vDT, AT}) where {vDT, AT} =
    judiVector{vDT, AT}("conj("*a.name*")",a.m,a.n,a.nsrc,a.geometry,a.data)

# transpose(jo)
transpose(a::judiVector{vDT, AT}) where {vDT, AT} =
    judiVector{vDT, AT}(""*a.name*".'",a.n,a.m,a.nsrc,a.geometry,a.data)

# adjoint(jo)
adjoint(a::judiVector{vDT, AT}) where {vDT, AT} =
        judiVector{vDT, AT}(""*a.name*".'",a.n,a.m,a.nsrc,a.geometry,a.data)

##########################################################

# Overload needed base function for SegyIO objects
vec(x::SegyIO.SeisCon) = vec(x[1].data)
abs(x::SegyIO.IBMFloat32) = abs(Float32(x))

##########################################################

for opo=[:+, :-, :*, :/]
    @eval begin
        $opo(a::judiVector{avDT, SeisCon}, ::T) where {avDT, T<:Real} = throw(DomainError(a, "Addition for OOC judiVectors not supported."))
        $opo(::T, b::judiVector{bvDT, SeisCon}) where {bvDT, T<:Real} = throw(DomainError(b, "Addition for OOC judiVectors not supported."))
        $opo(::judiVector{avDT, SeisCon}, b::judiVector{bvDT, SeisCon}) where {avDT, bvDT} = throw(DomainError(b, "Addition for OOC judiVectors not supported."))
    end
end

for ipop=[:lmul!, :rmul!, :rdiv!, :ldiv!]
    @eval begin
        $ipop(a::judiVector{avDT, SeisCon}, ::T) where {avDT, T<:Real} = throw(DomainError(a, "Addition for OOC judiVectors not supported."))
        $ipop(::T, b::judiVector{bvDT, SeisCon}) where {bvDT, T<:Real} = throw(DomainError(b, "Addition for OOC judiVectors not supported."))
        $ipop(::judiVector{avDT, SeisCon}, b::judiVector{bvDT, SeisCon}) where {avDT, bvDT} = throw(DomainError(b, "Addition for OOC judiVectors not supported."))
    end
end

# *(joLinearFunction, judiVector)
function *(A::joLinearFunction{ADDT,ARDT},v::judiVector{avDT, AT}) where {ADDT, ARDT, avDT, AT}
    A.n == size(v,1) || throw(judiVectorException("Shape mismatch: A:$((A.m, A.n)), v: $(size(v))"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiVector):",A.name,typeof(A),avDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(joLinearOperator, judiVector)
function *(A::joLinearOperator{ADDT,ARDT},v::judiVector{avDT, AT}) where {ADDT, ARDT, avDT, AT}
    A.n == size(v,1) || throw(judiVectorException("Shape mismatch: A:$((A.m, A.n)), v: $(size(v))"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiVector):",A.name,typeof(A),avDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# vcat
function vcat(ai::Vararg{judiVector{avDT, AT}, N}) where {avDT, AT, N}
    N == 1 && (return ai[1])
    N > 2 && (return vcat(ai[1], vcat(ai[2:end]...)))
    a, b = ai
    typeof(a.geometry) == typeof(b.geometry) || throw(judiVectorException("Geometry type mismatch"))
    m = a.m + b.m
    n = 1
    nsrc = a.nsrc + b.nsrc

    data = vcat(a.data, b.data)

    if AT == SegyIO.SeisCon
        nsamples = vcat(a.geometry.nsamples, b.geometry.nsamples)
    else
        xloc = vcat(a.geometry.xloc, b.geometry.xloc)
        yloc = vcat(a.geometry.yloc, b.geometry.yloc)
        zloc = vcat(a.geometry.zloc, b.geometry.zloc)
    end

    dt = vcat(a.geometry.dt, b.geometry.dt)
    nt = vcat(a.geometry.nt, b.geometry.nt)
    t = vcat(a.geometry.t, b.geometry.t)

    if AT == SegyIO.SeisCon
        geometry = GeometryOOC{Float32}(data,dt,nt,t,nsamples,a.geometry.key,a.geometry.segy_depth_key)
    else
        geometry = GeometryIC{Float32}(xloc, yloc, zloc, dt, nt, t)
    end

    return judiVector{avDT, AT}(a.name, m, n, nsrc, geometry, data)
end

# push!
function push!(a::judiVector{T, mT}, b::judiVector{T, mT}) where {T, mT}
    append!(a.data, b.data)
    a.nsrc += b.nsrc
    a.m += b.m
    push!(a.geometry, b.geometry)
end

# dot product
function dot(a::judiVector{avDT, AT}, b::judiVector{bvDT, AT}) where {avDT, bvDT, AT}
# Dot product for data containers
    size(a) == size(b) || throw(judiVectorException("Size mismatch: a:$(size(a)) b:$(size(b))"))
    compareGeometry(a.geometry, b.geometry) == 1 || throw(judiVectorException("Geometry mismatch"))
    dotprod = 0f0
    for j=1:a.nsrc
        dotprod += a.geometry.dt[j] * dot(vec(a.data[j]),vec(b.data[j]))
    end
    return dotprod
end

# norm
function norm(a::judiVector{avDT, AT}, p::Real=2) where {avDT, AT}
    if p == Inf
        return max([maximum(abs.(a.data[i])) for i=1:a.nsrc]...)
    end
    x = 0.f0
    for j=1:a.nsrc
        x += a.geometry.dt[j] * sum(abs.(vec(a.data[j])).^p)
    end
    return x^(1.f0/p)
end

# abs
function abs(a::judiVector{avDT, AT}) where {avDT, AT}
    b = deepcopy(a)
    for j=1:a.nsrc
        b.data[j] = abs.(a.data[j])
    end
    return b
end

# Subsample data container
"""
    subsample(x,source_numbers)

Subsample seismic data vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `judiVector`, `judiModeling`, \\
`judiProjection`, `judiJacobian`, `Geometry`, `judiRHS`, `judiPDE`, `judiPDEfull`.

Examples
========

(1) Extract 2 shots from `judiVector` vector:

    dsub = subsample(dobs,[1,2])

(2) Extract geometry for shot location 100:

    geometry_sub = subsample(dobs.geometry,100)

(3) Extract Jacobian for shots 10 and 20:

    Jsub = subsample(J,[10,20])

"""
function subsample(a::judiVector{avDT, AT},srcnum) where {avDT, AT}
    geometry = subsample(a.geometry,srcnum)     # Geometry of subsampled data container
    return judiVector(geometry,a.data[srcnum])
end

getindex(x::judiVector,a) = subsample(x,a)

# Create SeisBlock from judiVector container to write to file
function judiVector_to_SeisBlock(d::judiVector{avDT, AT}, q::judiVector{avDT, QT};
                                 source_depth_key="SourceSurfaceElevation",
                                 receiver_depth_key="RecGroupElevation") where {avDT, AT, QT}

    typeof(d.geometry) <: GeometryOOC && (d.geometry = Geometry(d.geometry))
    typeof(q.geometry) <: GeometryOOC && (q.geometry = Geometry(q.geometry))

    blocks = Array{Any}(undef, d.nsrc)
    count = 0
    for j=1:d.nsrc

        # create SeisBlock
        blocks[j] = SeisBlock(d.data[j])
        numTraces = size(d.data[j],2)
        traceNumbers = convert(Array{Integer,1},count+1:count+numTraces)

        # set headers
        set_header!(blocks[j], "GroupX", Int.(round.(d.geometry.xloc[j]*1f3)))
        if length(d.geometry.yloc[j]) == 1
            set_header!(blocks[j], "GroupY", Int(round(d.geometry.yloc[j][1]*1f3)))
        else
            set_header!(blocks[j], "GroupY", convert(Array{Integer,1},round.(d.geometry.yloc[j]*1f3)))
        end
        set_header!(blocks[j], receiver_depth_key, Int.(round.(d.geometry.zloc[j]*1f3)))
        set_header!(blocks[j], "SourceX", Int.(round.(q.geometry.xloc[j][1]*1f3)))
        set_header!(blocks[j], "SourceY", Int.(round.(q.geometry.yloc[j][1]*1f3)))
        set_header!(blocks[j], source_depth_key, Int.(round.(q.geometry.zloc[j][1]*1f3)))

        set_header!(blocks[j], "dt", Int(d.geometry.dt[j]*1f3))
        set_header!(blocks[j], "FieldRecord",j)
        set_header!(blocks[j], "TraceNumWithinLine", traceNumbers)
        set_header!(blocks[j], "TraceNumWithinFile", traceNumbers)
        set_header!(blocks[j], "TraceNumber", traceNumbers)
        set_header!(blocks[j], "ElevationScalar", -1000)
        set_header!(blocks[j], "RecSourceScalar", -1000)
        count += numTraces
    end

    # merge into single block
    fullblock = blocks[1]
    for j=2:d.nsrc
        fullblock = merge(fullblock,blocks[j])
        blocks[j] = []
    end
    return fullblock
end

function src_to_SeisBlock(q::judiVector{avDT, QT};
                                 source_depth_key="SourceSurfaceElevation",
                                 receiver_depth_key="RecGroupElevation") where {avDT, QT}

    return judiVector_to_SeisBlock(q, q;
        source_depth_key=source_depth_key,
        receiver_depth_key=receiver_depth_key)
end

function write_shot_record(srcGeometry, srcData, recGeometry, recData,options)
    q = judiVector(srcGeometry, srcData)
    d = judiVector(recGeometry, recData)
    pos = [srcGeometry.xloc[1][1], srcGeometry.yloc[1][1],  srcGeometry.zloc[1][1]]
    pos = join(["_"*string(trunc(p; digits=2)) for p in pos])
    file = join([string(options.file_name), pos,".segy"])
    block_out = judiVector_to_SeisBlock(d, q)
    segy_write(join([options.file_path,"/",file]), block_out)
    container = scan_file(join([options.file_path,"/",file]),
                          ["GroupX", "GroupY", "dt", "SourceSurfaceElevation", "RecGroupElevation"];
                          chunksize=256)
    return container
end

function time_resample!(x::judiVector, dt_new; order=2)
    x.m = 0
    for j=1:x.nsrc
        dataInterp, geom = time_resample(x.data[j], subsample(x.geometry, j), dt_new)
        x.data[j] = dataInterp
        x.geometry.dt[j] = dt_new
        x.geometry.nt[j] = geom.nt[1]
        x.geometry.t[j] = geom.t[1]
        x.m += prod(size(dataInterp))
    end
    return x
end

function time_resample(x::judiVector, dt_new; order=2)
    xout = deepcopy(x)
    time_resample!(xout, dt_new; order=order)
    return xout
end

function judiTimeInterpolation(geometry::Geometry, dt_coarse, dt_fine)
# Time interpolation as a linear operator (copies input data)

    nsrc = length(geometry.xloc)
    m = 0
    n = 0
    for j=1:nsrc
        nt_coarse = Int(trunc(geometry.t[j]/dt_coarse + 1))
        nt_fine = Int(trunc(geometry.t[j]/dt_fine + 1))
        n += length(geometry.xloc[j])*nt_coarse
        m += length(geometry.xloc[j])*nt_fine
    end
    I = joLinearFunctionFwd_T(m,n,
                              v -> time_resample(v, dt_fine),
                              w -> time_resample(w, dt_coarse),
                              Float32,Float32,name="Time resampling")
    return I
end

####################################################################################################
# Indexing

setindex!(x::judiVector, y, i) = x.data[i][:] = y

firstindex(x::judiVector) = 1

lastindex(x::judiVector) = x.nsrc

axes(x::judiVector) = Base.OneTo(x.nsrc)

ndims(x::judiVector) = length(size(x))

similar(x::judiVector) = 0f0*x

similar(x::judiVector, element_type::DataType, dims::Union{AbstractUnitRange, Integer}...) = 0f0*x

function fill!(x::judiVector{vDT, AT}, val) where {vDT, AT}
    for j=1:length(x.data)
        fill!(x.data[j], val)
    end
end

function sum(x::judiVector)
    s = 0f0
    for j=1:length(x.data)
        s += sum(vec(x.data[j]))
    end
    return s
end

Base.IteratorSize(d::judiVector) = Base.SizeUnknown()
isfinite(v::judiVector) = all(all(isfinite.(v.data[i])) for i=1:v.nsrc)
iterate(S::judiVector, state::Integer=1) = state > S.nsrc ? nothing : (S.data[state], state+1)

####################################################################################################

# Integration/differentiation of shot records

function cumsum(x::judiVector;dims=1)
    y = deepcopy(x)
    cumsum!(y, x; dims=dims)
    return y
end

function cumsum!(y::judiVector, x::judiVector;dims=1)
    dims == 1 || dims == 2 || throw(judiVectorException("Dimension $(dims) is out of range for a 2D array"))
    h = dims == 1 ? x.geometry.dt[1] : 1f0              # receiver dimension is non-dimensional
    for i = 1:x.nsrc
        cumsum!(y.data[i], x.data[i], dims=dims)
    end
    rmul!(y, h)
    return y
end

function diff(x::judiVector;dims=1)
    # note that this is not the same as default diff in julia, the first entry stays in the diff result
    dims == 1 || dims == 2 || throw(judiVectorException("Dimension $(dims) is out of range for a 2D array"))
    y = (dims == 1 ? 1f0 / x.geometry.dt[1] : 1f0) * x        # receiver dimension is non-dimensional
    for i = 1:x.nsrc
        copy!(selectdim(y.data[i], dims, 2:size(y.data[i],dims)), diff(y.data[i],dims=dims))
    end
    return y
end

####################################################################################################

BroadcastStyle(::Type{judiVector}) = Base.Broadcast.DefaultArrayStyle{1}()

ndims(::Type{judiVector{Float32, Array{Float32, 2}}}) = 1

### +/- ####
broadcasted(::typeof(+), x::judiVector, y::judiVector) = x + y
broadcasted(::typeof(-), x::judiVector, y::judiVector) = x - y

broadcasted(::typeof(+), x::judiVector, y::Number) = x + y
broadcasted(::typeof(-), x::judiVector, y::Number) = x - y

broadcasted(::typeof(+), x::Number, y::judiVector) = x + y
broadcasted(::typeof(-), x::Number, y::judiVector) = x - y

### * ####
function broadcasted(::typeof(*), x::judiVector, y::judiVector)
    size(x) == size(y) || throw(judiVectorException("Size mismatch: x:$(size(x)), y:$(size(y))"))
    compareGeometry(x.geometry, y.geometry) == 1 || throw(judiVectorException("Geometry mismatch"))
    typeof(x.data[1]) == SeisCon && throw("Addition for OOC judiVectors not supported.")
    typeof(y.data[1]) == SeisCon && throw("Addition for OOC judiVectors not supported.")
    z = deepcopy(x)
    for j=1:length(x.data)
        z.data[j] = x.data[j] .* y.data[j]
    end
    return z
end

function broadcasted(::typeof(*), x::judiVector, y::Number)
    z = deepcopy(x)
    for j=1:length(x.data)
        z.data[j] .*= y
    end
    return z
end

broadcasted(::typeof(*), x::Number, y::judiVector) = broadcasted(*, y, x)

### / ####
function broadcasted(::typeof(/), x::judiVector, y::judiVector)
    size(x) == size(y) || throw(judiVectorException("Size mismatch: x:$(size(x)), y:$(size(y))"))
    compareGeometry(x.geometry, y.geometry) == 1 || throw(judiVectorException("geometry mismatch"))
    typeof(x.data[1]) == SeisCon && throw("Addition for OOC judiVectors not supported.")
    typeof(y.data[1]) == SeisCon && throw("Addition for OOC judiVectors not supported.")
    z = deepcopy(x)
    for j=1:length(x.data)
        z.data[j] = x.data[j] ./ y.data[j]
    end
    return z
end

broadcasted(::typeof(/), x::judiVector, y::Number) = broadcasted(*, x, 1/y)

# Materialize for broadcasting
function materialize!(x::judiVector, y::judiVector)
    for j=1:length(x.data)
        x.data[j] .= y.data[j]
    end
end

function broadcast!(identity, x::judiVector, y::judiVector)
    copy!(x,y)
end

function broadcasted(identity, x::judiVector)
    return x
end

function copy!(x::judiVector, y::judiVector)
    for j=1:x.nsrc
        x.data[j] .= y.data[j]
    end
    x.geometry = deepcopy(y.geometry)
end

copy(x::judiVector) = 1f0 * x

function get_data(x::judiVector{T, SeisCon}) where T
    shots = Array{Array{Float32, 2}, 1}(undef, x.nsrc)
    rec_geometry = Geometry(x.geometry)
    for j=1:x.nsrc
        shots[j] = convert(Array{Float32, 2}, x.data[j][1].data)
    end
    return judiVector(rec_geometry, shots)
end

get_data(x::judiVector{T, Array{Float32, 2}}) where T = x

function convert_to_array(x::judiVector)
    y = vec(x.data[1])
    if x.nsrc > 1
        for j=2:x.nsrc
            y = [y; vec(x.data[j])]
        end
    end
    return y
end

function isapprox(x::judiVector, y::judiVector; rtol::Real=sqrt(eps()), atol::Real=0)
    compareGeometry(x.geometry, y.geometry) == 1 || throw(judiVectorException("Geometry mismatch"))
    isapprox(x.data, y.data; rtol=rtol, atol=atol)
end

############################################################

function A_mul_B!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.weights)
        x.weights[j] .= z.weights[j]
    end
end

function A_mul_B!(x::judiVector, F::Union{joAbstractLinearOperator, joLinearFunction}, y::Union{Array, PhysicalParameter})
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.data)
        x.data[j] .= z.data[j]
    end
end

function A_mul_B!(x::Union{Array, PhysicalParameter}, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector)
    F.m == size(y, 1) ? x[:] .= (adjoint(F)*y)[:] : x[:] .= (F*y)[:]
end

function A_mul_B!(x::judiVector, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.data)
        x.data[j] .= z.data[j]
    end
end

function A_mul_B!(x::judiVector, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector)
    F.m == size(y, 1) ? z = adjoint(F)*y : z = F*y
    for j=1:length(x.data)
        x.data[j] .= z.data[j]
    end
end

mul!(x::judiWeights, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector) = A_mul_B!(x, F, y)
mul!(x::judiVector, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiWeights) = A_mul_B!(x, F, y)
mul!(x::judiVector, F::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector) = A_mul_B!(x, F, y)
mul!(x::Union{Array, PhysicalParameter}, J::Union{joAbstractLinearOperator, joLinearFunction}, y::judiVector) = A_mul_B!(x, J, y)
mul!(x::judiVector, J::Union{joAbstractLinearOperator, joLinearFunction}, y::Union{Array, PhysicalParameter}) = A_mul_B!(x, J, y)

# Rebuild for backward compatinility
convgeom(x) = GeometryIC{Float32}([getfield(x.geometry, s) for s=fieldnames(GeometryIC)]...)
convdata(x) = convert(Array{Array{Float32, 2}, 1}, x.data)

"""
    reuild_jv(v)
rebuild a judiVector from previous version type or JLD2 reconstructed type
"""
rebuild_jv(v::judiVector{T, AT}) where {T, AT} = v
rebuild_jv(v) = judiVector(convgeom(v), convdata(v))


text:
using HTTP
url = "https://julialang.org"
r = HTTP.get(url)
# HTTP.Messages.Response:
# """
# HTTP/1.1 200 OK
# Connection: keep-alive
# Content-Length: 36344
# Server: GitHub.com
# Content-Type: text/html; charset=utf-8
# x-origin-cache: HIT
# Last-Modified: Thu, 31 Mar 2022 15:26:14 GMT
# Access-Control-Allow-Origin: *
# ETag: "6245c816-8df8"
# expires: Mon, 04 Apr 2022 08:53:51 GMT
# ...
println(r.status) # 200
println(String(r.body))
# <!doctype html> <html lang=en > <meta charset=utf-8 >
# <meta name=viewport  content="width=device-width, initial-scale=1, shrink-to-fit=no">
# <meta http-equiv=x-ua-compatible  content="ie=edge">
# <meta name=author  content="Jeff Bezanson, Stefan Karpinski, Viral Shah, Alan Edelman, et al.">
# <meta name=description  content="The official website for the Julia Language. Julia is a language that is fast,
# dynamic, easy to use, and open source.
println(r.headers)
# Pair{SubString{String}, SubString{String}}["Connection" => "keep-alive", "Content-Length" => "36344", "Server" => "GitHub.com", "Content-Type" => "text/html; charset=utf-8", "x-origin-cache" => "HIT", "Last-Modified" => "Thu, 31 Mar 2022 15:26:14 GMT", "Access-Control-Allow-Origin" => "*", "ETag" => "\"6245c816-8df8\"", "expires" => "Mon, 04 Apr 2022 08:53:51 GMT", "Cache-Control" => "max-age=600", "x-proxy-cache" => "MISS", "X-GitHub-Request-Id" => "F5A0:3CB5:DC5B19:E5B057:624AAFC5", "Via" => "1.1 varnish, 1.1 varnish", "X-Fastly-Request-ID" => "7ce247d3564f7d97a13fb235789763f70b46fc0a", "Fastly-Debug-States" => "DELIVER", "Accept-Ranges" => "bytes", "Date" => "Mon, 04 Apr 2022 08:43:51 GMT", "Age" => "0", "X-Served-By" => "cache-bru1480020-BRU, cache-bru1480075-BRU", "X-Cache" => "MISS, MISS", "X-Cache-Hits" => "0, 0", "X-Timer" => "S1649061831.488585,VS0,VE101", "Vary" => "Accept-Encoding", "Strict-Transport-Security" => "max-age=5356800"]

HTTP.post("http://httpbin.org/post", [], "post body data")
# HTTP.Messages.Response:
# """
# HTTP/1.1 200 OK
# Date: Mon, 04 Apr 2022 08:56:03 GMT
# Content-Type: application/json
# Content-Length: 362
# Connection: keep-alive
# Server: gunicorn/19.9.0
# Access-Control-Allow-Origin: *
# Access-Control-Allow-Credentials: true

# {
#   "args": {},
#   "data": "post body data",
#   "files": {},
#   "form": {},
#   "headers": {
#     "Accept": "*/*",
#     "Content-Length": "14",
#     "Host": "httpbin.org",
#     "User-Agent": "HTTP.jl/1.7.0",
#     "X-Amzn-Trace-Id": "Root=1-624ab2a3-0d8397831936742012f8dd3f"
#   },
#   "json": null,
#   "origin": "81.83.67.231",
#   "url": "http://httpbin.org/post"
# }
# """


text:
export execute, initialise_sampler, timer_output

using FileIO
using TimerOutputs
using CUDA

using QXContexts.Param
using QXContexts.ComputeGraphs
using QXContexts.Sampling

# const timer_output = TimerOutput()
# if haskey(ENV, "QXRUN_TIMER")
#     timeit_debug_enabled() = return true # Function fails to be defined by enable_debug_timings
#     TimerOutputs.enable_debug_timings(Execution)
# end

"""
    write_results(results, output_file)

Save results from calculations for the given
"""
function write_results(results, output_file)
    @assert splitext(output_file)[end] == ".jld2" "Output file should have jld2 suffix"
    if results !== nothing
        save(output_file, "results", results)
    end
end

"""
    initialise_sampler(dsl_file::String,
                       input_file::String,
                       param_file::String;
                       use_mpi::Bool=false,
                       sub_comm_size::Int=1,
                       use_gpu::Bool=false,
                       elt::Type=ComplexF32)

Initialise the sampler
"""
function initialise_sampler(dsl_file::String,
                            input_file::String,
                            param_file::String;
                            use_mpi::Bool=false,
                            sub_comm_size::Int=1,
                            use_gpu::Bool=false,
                            elt::Type=ComplexF32)
    # read dsl file of commands and data file with initial tensors
    @timeit "Parse input files" cg, _ = parse_dsl_files(dsl_file, input_file)

    # Create a context to execute the commands in
    T = if use_gpu
        @assert CUDA.functional() "CUDA installation is not functional, ensure you have a GPU and appropriate drivers"
        CuArray{elt}
    else
        Array{elt}
    end
    @timeit "Create Context" begin
        ctx = QXContext{T}(cg)
        if use_mpi
            ctx = QXMPIContext(ctx, sub_comm_size=sub_comm_size)
        end
    end

    # read sampler parameters from parameter file
    sampler_args = parse_parameters(param_file)

    # Create a sampler to produce bitstrings to get amplitudes for and a variable to store
    # the results.
    @timeit "Create sampler" create_sampler(ctx, sampler_args)
end

"""
    execute(dsl_file::String,
            input_file::Union{String, Nothing}=nothing,
            param_file::Union{String, Nothing}=nothing,
            output_file::String="";
            max_amplitudes::Union{Int, Nothing}=nothing,
            max_slices::Union{Int, Nothing}=nothing,
            timings::Bool=false,
            kwargs...)

Main entry point for running calculations. Loads input data, runs computations and
writes results to output files.
"""
function execute(dsl_file::String,
                 input_file::Union{String, Nothing}=nothing,
                 param_file::Union{String, Nothing}=nothing,
                 output_file::String="";
                 max_amplitudes::Union{Int, Nothing}=nothing,
                 max_slices::Union{Int, Nothing}=nothing,
                 timings::Bool=false,
                 kwargs...)

    if input_file === nothing
        input_file = splitext(dsl_file)[1] * ".jld2"
    end

    if param_file === nothing
        param_file = splitext(dsl_file)[1] * ".yml"
    end
    @timeit "Init sampler" sampler = initialise_sampler(dsl_file, input_file, param_file;
                                 kwargs...)

    @timeit "Simulation" results = sampler(max_amplitudes=max_amplitudes,
                                           max_slices=max_slices)

    if output_file != "" && results !== nothing
        @timeit "Write results" write_results(results, output_file)
        if timings
            print_timer()
        end
    end
    results
end


text:
Base.promote_rule(::Type{PhysicalConstant{name,S,D,SU}},
                  ::Type{Quantity{T,D,TU}}) where {name,S,T,D,SU,TU} =
                      promote_type(Quantity{S,D,SU}, Quantity{T,D,TU})
Base.promote_rule(::Type{Quantity{T,D,TU}},
                  ::Type{PhysicalConstant{name,S,D,SU}}) where {name,S,T,D,SU,TU} =
                      promote_type(Quantity{S,D,SU}, Quantity{T,D,TU})

Base.promote_rule(::Type{PhysicalConstant{name,S,NoDims,U}},
                  ::Type{T}) where {name,S,U,T<:Number} =
                      promote_type(S, T)
Base.promote_rule(::Type{PhysicalConstant{name,S,NoDims,U}},
                  ::Type{Quantity{T,NoDims,TU}}) where {name,S,U,T,TU} =
                      promote_type(S, T)

Base.convert(::Type{T},
             x::PhysicalConstant{name,S,NoDims,U}) where {T<:AbstractFloat,name,S,U} =
                 float(T, x)

Unitful.uconvert(u::Unitful.Units, c::PhysicalConstant{name,T,D,TU}) where {name,T,D,TU} =
    uconvert(u, float(T, c))

max_issues_repo_name: JuliaPhysics/Constants

max_forks_repo_name: JuliaPhysics/Constants


text:
@testset "Test inverse of MapFromFunc" begin
  f = MapFromFunc(x -> x+1, x -> x-1, ZZ, ZZ)
  x = ZZ(1)
  y = f(x)
  @test preimage(f, y) == x

  finv = inv(f)
  y = finv(x)
  @test preimage(finv, y) == x
end


text:
@with_kw mutable struct BatchRegularizer
    buffers # Datasets to use in the regularization
    𝒟s = nothing# Batches of those datasets
    batch_size::Int # Batch size of the data sets
    λ::Float32 = 1f0# regularization coefficients
    loss # Loss function that takes args (π, 𝒟)
end

value_regularization(π, 𝒟) = Flux.mse(value(π, 𝒟[:s]), 𝒟[:value])
action_regularization(π, 𝒟) = Flux.mse(action(π, 𝒟[:s]), 𝒟[:a])
action_value_regularization(π, 𝒟) = begin
    v = value(π, 𝒟[:s], 𝒟[:a])
    v isa Tuple && (v = v[1])
    Flux.mse(v, 𝒟[:value])
end

function (R::BatchRegularizer)(π)
    # sample a random batch for each buffer
    ret0 = false
    ignore() do
        if any([length(b) == 0 for b in R.buffers])
            ret0 = true
        else
            isnothing(R.𝒟s) && (R.𝒟s = [buffer_like(b, capacity=R.batch_size, device=device(π)) for b in R.buffers])
            for (𝒟, buffer) in zip(R.𝒟s, R.buffers)
                rand!(𝒟, buffer)
            end
        end
    end

    # Return the mean
    ret0 ? 0f0 : R.λ*mean([R.loss(π, 𝒟) for 𝒟 in R.𝒟s])
end


text:
@testset "black box variational inference" begin

    Random.seed!(1)

    @gen function model()
        slope = @trace(normal(-1, exp(0.5)), :slope)
        intercept = @trace(normal(1, exp(2.0)), :intercept)
    end

    @gen function approx()
        @param slope_mu::Float64
        @param slope_log_std::Float64
        @param intercept_mu::Float64
        @param intercept_log_std::Float64
        @trace(normal(slope_mu, exp(slope_log_std)), :slope)
        @trace(normal(intercept_mu, exp(intercept_log_std)), :intercept)
    end
    init_param!(approx, :slope_mu, 0.)
    init_param!(approx, :slope_log_std, 0.)
    init_param!(approx, :intercept_mu, 0.)
    init_param!(approx, :intercept_log_std, 0.)

    observations = choicemap()
    update = ParamUpdate(GradientDescent(0.001, 100000), approx)
    black_box_vi!(model, (), observations, approx, (), update;
        iters=500, samples_per_iter=1000, verbose=false)
    slope_mu = get_param(approx, :slope_mu)
    slope_log_std = get_param(approx, :slope_log_std)
    intercept_mu = get_param(approx, :intercept_mu)
    intercept_log_std = get_param(approx, :intercept_log_std)
    @test isapprox(slope_mu, -1., atol=0.001)
    @test isapprox(slope_log_std, 0.5, atol=0.001)
    @test isapprox(intercept_mu, 1., atol=0.001)
    @test isapprox(intercept_log_std, 2.0, atol=0.001)
end

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
module SymFinel

    using Calculus

    include("CoordExpr.jl")
    include("Interpolator.jl")
    include("Regressor.jl")
    include("Domain.jl")
    include("Finel.jl")
    include("Residuals.jl")

end # module

module SymFinelTests

    using Calculus

    include("CoordExpr.jl")
    include("Interpolator.jl")
    include("Regressor.jl")
    include("Domain.jl")
    include("Finel.jl")
    include("Residuals.jl")

    export get_basefun, get_basefun_derivative, get_linpar_contribution, get_interp_linmat,
        _get_gradfunc, regress, finel_get_regr, get_abs_coords!, from_finel_jac!, get_punctual_resfunc,
            vals_to_polys, polys_to_vars

end # module


text:
#=
Using a function rand5() that returns an integer from 1 to 5 (inclusive) with uniform probability, implement a function rand7() that returns an integer from 1 to 7 (inclusive).
=#

function rand5()
    return rand(1:5)
end

function rand7()
    x = 100
    # Create a number between 1 to 21
    while x > 21
        x = rand5() * 5 + rand5() - 5 #uniform probability of having 5, 10, 15, 20 and 25 minus a number between 0 and 4
    end
    return x % 7 + 1
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
@doc doc"""
```Julia
VSIVoltagePT1(;τ_v,τ_P,τ_Q,K_P,K_Q,E_r,P,Q)
```
A node type that applies the frequency and voltage droop control to control the frequency and
voltage dynamics. Implemented according to  Schiffer et. al., "Conditions for stability of
droop-controlled inverter-based microgrids", Automatica, 2014.
`VSIVoltagePT1` models an inverters as AC voltage source which means the amplitude and
frequency can defined by the designer (often called grid-forming inverter mode).
The frequency regulation is assumed to be instantaneous, but the voltage control happens with a delay
`τ_v` that is represented by a first order filter.
In addition simple proportional controllers are implemented for frequency and voltage such that the frequency `ω` and
voltage amplitudes `v` of the inverters are modified depending on the deviations (with respect to a desired value) of the active and reactive powers, respectively.
ift is assumed that active and reactive power are measured via low pass fileters with time constant `τ_P` and `τ_Q`, respectively.
Hence, additionally to `u`, it has the internal dynamic variables
* `ω` representing the frequency of the inverter relative to the grid frequency ``Ω=2π50``Hz, i.e. the real frequency ``ω_r`` of the inverter is given as ``ω_r = Ω + ω``.
* `q_m` is the measured reactive power at the grid connection point.
# Keyword Arguments
- `τ_v`: time constant voltage control delay
- `τ_P`: time constant active power measurement
- `τ_Q`: time constant reactive power measurement
- `K_P`: droop constant frequency droop
- `K_Q`: droop constant voltage droop
- `V_r`: reference/ desired voltage
- `P`: active (real) power infeed
- `Q`: reactive (imag) power infeed
# Mathematical Representation
Using `VSIVoltagePT1` for node ``a`` (according to J. Schiffer et. al., eq. (6)) gives the equations
```math
\dot{\phi}_a=\omega_a\\
 \dot{\omega}_a=\frac{1}{\tau_{P,a}}[-\omega_a-K_{P,a} (\Re\left(u_a \cdot i_a^*\right)-P_{ref,a})]\\
 \tau_v\dot{v}_{a}=-v_a+V_{ref}-K_{Q,a}(q_{m,a}-Q_{ref,a})\\
 \tau_Q \dot{q}_{m,a}=-q_{m,a}+\Im\left(u_a \cdot i_a^*\right)\\
 \dot{u}_a=\dot{v_a}e^{j\phi}+j\omega_a u_a\\
```
In general ``τ_V ≪ τ_P``, assuming ``τ_V = 0`` would then lead to `VSIMinimal`.
"""
@DynamicNode dVOC(τ_v,τ_P,τ_Q,K_P,K_Q,V_r,P,Q) begin
    @assert τ_v > 0 "time constant voltage droop delay should be >0"
    @assert τ_P > 0 "time constant active power measurement should be >0"
    @assert τ_Q > 0 "time constant reactive power measurement should be >0"
    @assert K_Q > 0 "reactive power droop constant should be >0"
    @assert K_P > 0 "active power droop constant reactive power measurement should be >0"
end [[ω, dω],[q_m,dq_m]] begin
    i_c = 1im * i * exp(-1im * θ)
    e_c = 1im * u * exp(-1im * θ)
    p = real(u * conj(i))
    e_d = 0                         # simplification from K.Schmietendorf paper
    e_q = imag(e_c)
    i_d = real(i_c)
    i_q = imag(i_c)

    dθ = ω
    dω = (P - D * ω - p - (X_q_dash - X_d_dash) * i_d * i_q) * Ω_H
    de_q = (1 / T_d_dash) * (E_f - e_q + i_d * (X_d - X_d_dash))
    du = de_q * exp(1im * θ) + u * 1im * ω
end

export dVOC

max_stars_repo_path: example_cases/new node types/dVOC.jl

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_path: example_cases/new node types/dVOC.jl

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_path: example_cases/new node types/dVOC.jl

max_forks_repo_name: luap-pik/PowerDynIC


text:
using Test
using Statistics
using MIPVerify

@testset "import_datasets.jl" begin
    @testset "NamedTrainTestDataset" begin
        @testset "Base.show" begin
            io = IOBuffer()
            mnist = read_datasets("mnist")
            Base.show(io, mnist)
            @test String(take!(io)) == """
                mnist:
                  `train`: {LabelledImageDataset}
                    `images`: 60000 images of size (28, 28, 1), with pixels in [0.0, 1.0].
                    `labels`: 60000 corresponding labels, with 10 unique labels in [0, 9].
                  `test`: {LabelledImageDataset}
                    `images`: 10000 images of size (28, 28, 1), with pixels in [0.0, 1.0].
                    `labels`: 10000 corresponding labels, with 10 unique labels in [0, 9]."""
        end
    end

    @testset "LabelledImageDataset" begin
        @testset "Base.show" begin
            io = IOBuffer()
            mnist = read_datasets("mnist")
            Base.show(io, mnist.train)
            @test String(take!(io)) == """
                {LabelledImageDataset}
                    `images`: 60000 images of size (28, 28, 1), with pixels in [0.0, 1.0].
                    `labels`: 60000 corresponding labels, with 10 unique labels in [0, 9]."""
        end
    end

    @testset "read_datasets" begin
        @testset "MNIST" begin
            d = read_datasets("mnist")
            @test MIPVerify.num_samples(d.test) == 10000
            @test MIPVerify.num_samples(d.train) == 60000
            # Some sanity checks to make sure that the data imported looks like we expect.
            @test d.test.labels[1:10] == [7, 2, 1, 0, 4, 1, 4, 9, 5, 9]
            @test d.train.labels[1:10] == [5, 0, 4, 1, 9, 2, 1, 3, 1, 4]
            @test d.test.images |> mean ≈ 0.13251463f0
            @test d.train.images |> mean ≈ 0.1307003f0
            @test size(d.test.images) == (10000, 28, 28, 1)
            @test size(d.train.images) == (60000, 28, 28, 1)
            @test d.test.images[1, 8, 9, 1] ≈ 0.62352943f0
            @test d.test.images[1, :, :, :] |> mean ≈ 0.09230693f0
        end

        @testset "CIFAR10" begin
            d = read_datasets("cifar10")
            @test MIPVerify.num_samples(d.test) == 10000
            @test MIPVerify.num_samples(d.train) == 50000
            # Some sanity checks to make sure that the data imported looks like we expect.
            @test d.test.labels[1:10] == [3, 8, 8, 0, 6, 6, 1, 6, 3, 1]
            @test d.train.labels[1:10] == [6, 9, 9, 4, 1, 1, 2, 7, 8, 3]
            @test d.test.images |> mean ≈ 0.4765849205984478
            @test d.train.images |> mean ≈ 0.4733630004850899
            @test size(d.test.images) == (10000, 32, 32, 3)
            @test size(d.train.images) == (50000, 32, 32, 3)
            @test d.test.images[1, 8, 9, 1] ≈ 0.6705882352941176
            @test d.test.images[1, :, :, :] |> mean ≈ 0.42504340277777786
        end

        @testset "Unsupported dataset" begin
            @test_throws DomainError read_datasets("the angry hippo of the nile")
        end
    end
end


text:
#### This example shows the growth of the number of quadratic functions in the constrained case
#### The figure in the article was generated with "OPTIMIZE = false"

using EllZeroTrendFiltering
using DSP
using Plots
using Interpolations

### START TEST FUNCTIONS

t = 1:200
srand(1)

circle_segment(N) = sin.(acos.(range(-1, stop=1, length=N)))
linear_trend(N) = (0:N-1)

g_cll = [0.2*linear_trend(120) + circle_segment(120);
     -6 + 8*linear_trend(40);
     -6 + 8*linear_trend(40)]

g_cllc = [5*circle_segment(70);
    -5 + 4*linear_trend(30);
    -5 + 4*linear_trend(30);
     5*circle_segment(70)]

g_rand = () -> randn(200)

g_walk = () -> cumsum(randn(200))

g_square = range(0, stop=4, length=200).^2

g_cube = range(0, stop=4, length=200).^3

g_4 = range(0, stop=4, length=200).^4

g_5 = range(0, stop=4, length=200).^5

g_exp = exp.(range(0, stop=4, length=200))

g_super_exp = exp.(range(0, stop=3, length=200).^2)

g_super_exp2 = [exp.(range(0, stop=2, length=100).^2)   ;  exp.(range(2, stop=0.5, length=100).^2)]

g_super_exp3 = [exp.(range(0, stop=2, length=100).^2)   ; -exp.(range(2, stop=0.5, length=100).^2)]

g_super_exp4 = [exp.(range(2, stop=0.5, length=100).^2) ;  exp.(range(0, stop=2, length=100).^2)]

g_super_exp5 = [exp.(range(2, stop=0, length=100).^2)   ; 2-exp.(range(0, stop=0, length=100).^2)]

function randlin(N)
    t_grid = cumsum( rand(10:25, N) )
    y_grid = randn(N)
    g = interpolate((t_grid,), y_grid, Gridded(Linear()))[t]
    reverse!(g)
end

function randlinfilt1(N)
    t_grid = cumsum( rand(10:25, N) )
    y_grid = randn(N)
    g = interpolate((t_grid,), y_grid, Gridded(Linear()))[t]
    g = reverse(filt(ones(10), 1, reverse(filt(ones(10), 1, g))))
    reverse!(g)
end

function randlinfilt2(N)
    t_grid = cumsum( rand(10:25, N) )
    y_grid = randn(N)
    g = interpolate((t_grid,), y_grid, Gridded(Linear()))[t]
    g =filt(ones(40), 1, g)
end

function randlinfilt3(N)
    t_grid = cumsum( rand(10:25, N) )
    y_grid = randn(N)
    g = interpolate((t_grid,), y_grid, Gridded(Linear()))[t]
    g = filtfilt(ones(10), g)
end

g_data = [g_cll, g_cllc,
        g_exp, g_super_exp, g_super_exp2,
        g_super_exp3, g_super_exp4, g_super_exp5,
        g_square, g_cube, g_4, g_5,
        reverse!(g_cll), reverse!(g_cllc),
        reverse!(g_exp), reverse!(g_super_exp), reverse!(g_super_exp2),
        reverse!(g_super_exp3), reverse!(g_super_exp4), reverse!(g_super_exp5),
        reverse!(g_square), reverse!(g_cube), reverse!(g_4), reverse!(g_5),
        g_rand(), g_rand(), g_rand(), g_rand(), g_rand(), g_rand(), g_rand(), g_rand(), g_rand(), g_rand(),
        g_walk(), g_walk(), g_walk(), g_walk(), g_walk(), g_walk(), g_walk(), g_walk(), g_walk(), g_walk(),
        randlin(5), randlinfilt1(5), randlinfilt2(5), randlinfilt3(5),
        randlin(10), randlinfilt1(10), randlinfilt2(10), randlinfilt3(10),
        randlin(20), randlinfilt1(20), randlinfilt2(20), randlinfilt3(20),
        randlin(40), randlinfilt1(40), randlinfilt2(40), randlinfilt3(40),
        reverse!(randlin(5)), reverse!(randlinfilt1(5)), reverse!(randlinfilt2(5)), reverse!(randlinfilt3(5)),
        reverse!(randlin(10)), reverse!(randlinfilt1(10)), reverse!(randlinfilt2(10)), reverse!(randlinfilt3(10)),
        reverse!(randlin(20)), reverse!(randlinfilt1(20)), reverse!(randlinfilt2(20)), reverse!(randlinfilt3(20)),
        reverse!(randlin(40)), reverse!(randlinfilt1(40)), reverse!(randlinfilt2(40)), reverse!(randlinfilt3(40)),
        cumsum(randlin(5)), cumsum(randlinfilt1(5)), cumsum(randlinfilt2(5)), cumsum(randlinfilt3(5)),
        cumsum(randlin(10)), cumsum(randlinfilt1(10)), cumsum(randlinfilt2(10)), cumsum(randlinfilt3(10)),
        cumsum(randlin(20)), cumsum(randlinfilt1(20)), cumsum(randlinfilt2(20)), cumsum(randlinfilt3(20)),
        cumsum(randlin(40)), cumsum(randlinfilt1(40)), cumsum(randlinfilt2(40)), cumsum(randlinfilt3(40)),
        ]
### END TEST FUNCTIONS

nbr_prob = length(g_data)

V = zeros(length(t), nbr_prob)

g_save = zeros(length(t), nbr_prob)

for prob_i=1:nbr_prob
    g = g_data[prob_i]

    l = compute_discrete_transition_costs(g, t);
    cost_last = QuadraticPolynomial(0.0, 0.0, 0.0)

    Λ =  pwq_dp_constrained(l, cost_last, 15);

    l = [(isassigned(Λ,m,i) ? length(Λ[m,i]) : 0) for m = 1:size(Λ,1), i = 1:size(Λ,2)]

    V[:, prob_i] .= maximum(l, 1)[:]
    g_save[:, prob_i] .= g
end

gr()
plot(size=(800,400));
plot(t, V, lab="", linealpha=0.2, c=:black);
plot!(t, length(t) - t + 1, l=(:red, :dash), lab="")
savefig("complexity.png")
writecsv("complexity.csv", V)


text:
export ParamFcn, AddedParamFcn, reduction, null, Null, constant, Cst, linear, Lin, rbf, rbf′,
        ψ₀, ψj, ψpp1, ψ₀′, ψpp1′

import Base: +, *, show, zero, size, getindex, setindex!

# Type to hold 1D functions
abstract type ParamFcn end

zero(::Type{T}) where {T <: ParamFcn} = null()

## Need tools to identify

# struct AddedParamFcn{T <: Tuple} <: ParamFcn
#     list::T
# end

struct AddedParamFcn{T <: Tuple} <: ParamFcn
    list::T
end
function Base.show(io::IO, Σg::AddedParamFcn)
    println(io, "AddedParamFcn:")
    for g in Σg.list
        println(io, "  $g")
    end
end

# """
#     g₁::ParamFcn + g₂::ParamFcn
#
# Add the parameterized functions so that `(g₁ + g₂)(z) = g₁(z) + g₂(z)`.
# """

+(g::ParamFcn, Σg::AddedParamFcn) = AddedParamFcn((g, Σg.list...))
+(Σg::AddedParamFcn, g::ParamFcn) = AddedParamFcn((Σg.list..., g))
function +(Σg₁::AddedParamFcn, Σg₂::AddedParamFcn)
    AddedParamFcn((Σg₁..., Σg₂...))
end

+(g::ParamFcn...) = AddedParamFcn(g)

(Σg::AddedParamFcn)(z::T) where {T<:Real} = map(g->g(z), Σg.list)
 # mapreduce(g->g(z),+, Σg.list;init=0.0)#
reduction(Σg::AddedParamFcn, z::T) where {T<:Real} = mapreduce(g->g(z),+, Σg.list;init=0.0)

# Add tools to access elements of an AddedParamFcn
# Base.size(A::AddedParamFcn{T}) where {T<:Tuple} = size(A.list)
# Base.getindex(A::AddedParamFcn{T}, i::Vararg{Int,N}) where {T,N} = A.list[i...]
# setindex!(::AddedParamFcn{T}, ::S, ::Int64)

# (*)(α::T)(f::F)(z::U) where {T<:Real, F<:paramfunction, U<:Real} = α*f(z)
struct null <:ParamFcn
end

const Null = null()

(C::null)(z::T) where {T<:Real} = 0.0

function Base.show(io::IO, N::null)
println(io,"Zero")
end

struct constant <:ParamFcn
end

const Cst = constant()

(C::constant)(z::T) where {T<:Real} = 1.0

function Base.show(io::IO, C::constant)
println(io,"Constant")
end

struct linear <:ParamFcn
end

const Lin = linear()

(L::linear)(z::T) where {T<:Real} = z

function Base.show(io::IO, L::linear)
println(io,"Linear")
end

struct rbf <:ParamFcn
        μ::Float64
        σ::Float64
end

function Base.show(io::IO, N::rbf)
println(io,"Gaussian kernel with mean $(N.μ) and std $(N.σ)")
end

(N::rbf)(z::T)  where {T<:Real} = 1/(N.σ*√(2*π))*exp(-(0.5/N.σ^2)*(z-N.μ)^2)

struct rbf′ <:ParamFcn
        μ::Float64
        σ::Float64
end

function Base.show(io::IO, N::rbf′)
println(io,"Derivative of a Gaussian kernel with mean $(N.μ) and std $(N.σ)")
end

(N::rbf′)(z::T)  where {T<:Real} = (N.μ-z)/(N.σ^3*√(2*π))*exp(-(0.5/N.σ^2)*(z-N.μ)^2)

## Notation of Couplings for nonlinear ensemble filtering, Spantini, Baptista, Marzouk
struct ψ₀ <:ParamFcn
        ξ₀::Float64
        σ₀::Float64
        function ψ₀(ξ₀::Float64, σ₀::Float64)
            @assert σ₀ >= 0.0 "std must be >= 0"
            new(ξ₀, σ₀)
        end
        # function ψ₀(ξ₀::T, σ₀::S) where {T<:Real, S<:Real}
        #     @assert σ₀ >= 0.0 "std must be >= 0"
        #     new(ξ₀, σ₀)
        # end
end

function Base.show(io::IO, F::ψ₀)
println(io,"ψ₀ function with mean $(F.ξ₀) and std $(F.σ₀)")
end

function (F::ψ₀)(z::T) where {T<:Real}
    Δ₀ = (1/(F.σ₀*√2))*(z-F.ξ₀)
    return 0.5*((z-F.ξ₀)*(1-erf(Δ₀))-F.σ₀*√(2/π)*exp(-Δ₀^2))
end

## Notation of Couplings for nonlinear ensemble filtering, Spantini, Baptista, Marzouk
struct ψj <:ParamFcn
        ξⱼ::Float64
        σⱼ::Float64
        function ψj(ξⱼ::Float64, σⱼ::Float64)
            @assert σⱼ >= 0.0 "std must be >= 0"
            new(ξⱼ, σⱼ)
        end
end

function Base.show(io::IO, F::ψj)
println(io,"ψj function with mean $(F.ξⱼ) and std $(F.σⱼ)")
end

function (F::ψj)(z::T) where {T<:Real}
    Δⱼ = (1/(F.σⱼ*√2))*(z-F.ξⱼ)
    return 0.5*(1+erf(Δⱼ))
end

## Notation of Couplings for nonlinear ensemble filtering, Spantini, Baptista, Marzouk
#  ψₚ₊₁
struct ψpp1 <:ParamFcn
        ξpp1::Float64
        σpp1::Float64
        function ψpp1(ξpp1::Float64, σpp1::Float64)
            @assert σpp1 >= 0.0 "std must be >= 0"
            new(ξpp1, σpp1)
        end
end

function Base.show(io::IO, F::ψpp1)
println(io,"ψpp1 function with mean $(F.ξpp1) and std $(F.σpp1)")
end

function (F::ψpp1)(z::T) where {T<:Real}
    Δpp1 = (1/(F.σpp1*√2))*(z-F.ξpp1)
    return 0.5*((z-F.ξpp1)*(1+erf(Δpp1)) +F.σpp1*√(2/π)*exp(-Δpp1^2))
end

# " Derivatives ψ′ of ψ for j=0:p+1"
## Notation of Couplings for nonlinear ensemble filtering, Spantini, Baptista, Marzouk
#####  ψ₀′
struct ψ₀′ <:ParamFcn
        ξ₀::Float64
        σ₀::Float64
        function ψ₀′(ξ₀::Float64, σ₀::Float64)
            @assert σ₀ >= 0.0 "std must be >= 0"
            new(ξ₀, σ₀)
        end
end

function Base.show(io::IO, F::ψ₀′)
println(io,"ψ₀′ function with mean $(F.ξ₀) and std $(F.σ₀)")
end

function (F::ψ₀′)(z::T) where {T<:Real}
    Δ₀ = (1/(F.σ₀*√2))*(z-F.ξ₀)
    return 0.5*(1-erf(Δ₀))
end

### ψₚ₊₁′
struct ψpp1′ <:ParamFcn
        ξpp1::Float64
        σpp1::Float64
        function ψpp1′(ξpp1::Float64, σpp1::Float64)
            @assert σpp1 >= 0.0 "std must be >= 0"
            new(ξpp1, σpp1)
        end
end

function Base.show(io::IO, F::ψpp1′)
println(io,"ψpp1′ function with mean $(F.ξpp1) and std $(F.σpp1)")
end

function (F::ψpp1′)(z::T) where {T<:Real}
    Δpp1 = (1/(F.σpp1*√2))*(z-F.ξpp1)
    return 0.5*(1+erf(Δpp1))
end


text:
using GR
import SpecialFunctions

x = 0:0.1:20
for order in 0:3
  plot(x, SpecialFunctions.besselj.(order, x))
  hold(true)
end
legend("J_0", "J_1", "J_2", "J_3")


text:
"""
Update inversion variable with the initial model.
At the same time, update the lower and upper bounds on the inversion variable.
"""
function initialize!(pa::PFWI)

	randn!(pa.mx.last_x)  # reset last_x

	# use mod_initial as starting model according to parameterization
	copyto!(pa.mx.x, pa.mod_initial, pa.parameterization)

	# bounds
	Seismic_xbound!(pa.mx.lower_x, pa.mx.upper_x, pa)

end

# finalize xfwi or wfwi on pa
function finalize!(pa::PFWI, xminimizer)
	# update mediumm
	x_to_mediumm!(pa, xminimizer)

	# update calculated data using the minimizer
	F!(pa, nothing)

	# put minimizer into mediumi
	copyto!(pa.mediumi, xminimizer, pa.parameterization)

	# update initial model using minimizer, so that running xfwi! next time will start with updated model
	copyto!(pa.mod_initial, pa.mediumi)
end

"""
```julia
update!(pa)
```
Full Waveform Inversion using `Optim` and `Ipopt` packages.
This method updates `pa.mediumm` and `pa.dcal`. More details about the
optional parameters can be found in the documentation of the `Optim`
package.
The gradient is computed using the adjoint state method.
`pa.mediumi` is used as initial model if non-zero.

# Arguments

* `pa::PFWI` : inversion object (updated inside method)

# Keywords Arguments
* `solver` Choose `Symbol`
  * `=:ipopt` for `Ipopt.jl`
  * `=:optim` for `Optim.jl`

# Optional Arguments

* `optim_options` : see `Optim.jl` package for optimization options
* `optim_scheme=LBFGS()` : see `Optim.jl` for other options
* `bounded_flag=true` : constrained inversion using medium parameter bounds (`Ipopt.jl` only has bounded inversion)
* `ipopt_options` : see `Ipopt.jl`
"""
function update!(pa::PFWI{T1,T2,T3};
	   optim_scheme=LBFGS(),
	   optim_options=Optim.Options(show_trace=true,
	   store_trace=true,
	   extended_trace=false,
	   outer_iterations=2,
	   iterations=5,
	   f_tol=1e-5,
	   g_tol=1e-8,
	   x_tol=1e-5, ),
	   bounded_flag=true, solver=nothing,
	   ipopt_options=[["max_iter", 5]]) where {T1,T2,T3<:Union{LS,LS_prior}}

	global fwi_to
	reset_timer!(fwi_to)

	println("updating mediumm and mediumi...")
	println("> xfwi: number of inversion variables:\t", xfwi_ninv(pa))

	initialize!(pa)

	f(x) = ζfunc(x, pa.mx.last_x,  pa)
	g!(storage, x) = ζgrad!(storage, x, pa.mx.last_x, pa)
	if(!bounded_flag)
		(solver==:ipopt) && @info("bounded_flag not functional for Ipopt, switching to Optim")
		"""
		Unbounded LBFGS inversion, only for testing
		"""
		@timeit_debug fwi_to "xfwi!" begin
			res = optimize(f, g!, pa.mx.x, optim_scheme, optim_options)
		end
	else
		"""
		Bounded LBFGS inversion
		"""
		if (solver == :ipopt)
			function eval_f(x)
				return ζfunc(x, pa.mx.last_x,  pa)
			end

			function eval_grad_f(x, grad_f)
				ζgrad!(grad_f, x, pa.mx.last_x, pa)
			end

			function void_g(x, g)

			end

			function void_g_jac(x, mode, rows, cols, values)

			end

			global prob = createProblem(size(pa.mx.x)[1], pa.mx.lower_x, pa.mx.upper_x, 0,
				Array{Float64}(undef,0), Array{Float64}(undef,0), 0, 0,
					eval_f, void_g, eval_grad_f, void_g_jac, nothing)

			addOption(prob, "hessian_approximation", "limited-memory")

			if !(ipopt_options === nothing)
			    for i in 1:size(ipopt_options)[1]
				addOption(prob, ipopt_options[i][1], ipopt_options[i][2])
			    end
			end

			prob.x = pa.mx.x

			@timeit_debug fwi_to "xfwi!" begin
				res = solveProblem(prob)
			end
		else
			@timeit_debug fwi_to "xfwi!" begin
				res = optimize(f, g!, pa.mx.lower_x, pa.mx.upper_x,pa.mx.x, Fminbox(optim_scheme), optim_options)
			end
                end
	end
	println(fwi_to)

	# print some stuff after optimization...
	if (solver == :ipopt)
		pa.verbose && println(ApplicationReturnStatus[res])
	else
		pa.verbose && println(res)
	end

	# finalize optimization, using the optimizer...
	if solver == :ipopt
		finalize!(pa, prob.x)
	else
		finalize!(pa, Optim.minimizer(res))
	end

	return res
end

"""
```julia
update!(pa)
```
Return gradient at the first iteration, i.e., a migration image, when `attrib_inv==Migr()`.
"""
function update!(pa::PFWI{T1,T2,Migr}) where {T1,T2}

	global fwi_to
	reset_timer!(fwi_to)

	initialize!(pa)

	g1=pa.mx.gm[1]
	@timeit_debug fwi_to "xfwi!" begin
		grad!(g1, pa.mx.x, pa.mx.last_x,  pa)
	end
	println(fwi_to)

	# convert gradient vector to model
	visualize_gx!(pa.gmediumm, pa.mediumm, pa.gmediumi, pa.mediumi, g1, pa)

	return pa.gmediumi, g1
end

"""
```julia
update!(pa)
```
Return gradient at the first iteration, i.e., a migration image, without using
adjoint state method, when `attrib_inv==Migr_FD()`.
"""
function update!(pa::PFWI{T1,T2,Migr_FD}) where {T1,T2}
	global fwi_to
	reset_timer!(fwi_to)

	initialize!(pa)
	gx=pa.mx.gm[1]

	println("> number of functional evaluations:\t", xfwi_ninv(pa))

	f(x) = ζfunc(x, pa.mx.last_x,  pa)

	gx = Calculus.gradient(f,pa.mx.x) # allocating new gradient vector here

	println(fwi_to)

	# convert gradient vector to model
	visualize_gx!(pa.gmediumm, pa.mediumm, pa.gmediumi, pa.mediumi, gx, pa)

	return pa.gmediumi, gx
end

"""
Return posterior covariance matrix and estimate as given by Best Linear Unbiased Analysis.
    # Arguments
    - H_obs : observation operator
    - C_obs : observation error covariance matrix
    - C_prior : prior covariance matrix
    - est_0 : prior estimate
    - y_obs : observation vector
    - compute_update : if true returns posterior mean and covariance, if false only covariance
    - diag : if true, assumes that observation covariance matrix is diag and uses it to siply inverse calculation

    # Output
    Posterior covariance and mean (if compute_update = true)
"""
function blue(H_obs, C_obs, C_prior, est_0, y_obs; compute_update = false, diag = false)
    C_prior_inv = inv(C_prior)
    C_prior_inv = 0.5 * (C_prior_inv + C_prior_inv')

    if diag
        C_obs_inv = 1 ./ C_obs
    else
        C_obs_inv = inv(C_obs)
        C_obs_inv = 0.5 * (C_obs_inv + C_obs_inv')
    end

    C_post = inv((H_obs' * C_obs_inv * H_obs + C_prior_inv))

    if compute_update
        est = C_post * (C_prior_inv * est_0 + H_obs' * C_obs_inv * y_obs)
        return C_post, est;
    else
        return C_post;
    end
end

		#=
		Harvest the Optim result to plot functional and gradient --- to be updated later
		if(extended_trace)
			# convert gradient vector to model
			gmediumi = [Medium_zeros(pa.mediumi.mgrid) for itr=1:Optim.iterations(res)]
			gmediumm = [Medium_zeros(pa.mediumm.mgrid) for itr=1:Optim.iterations(res)]
			mediumi = [Medium_zeros(pa.mediumi.mgrid) for itr=1:Optim.iterations(res)]
			mediumm = [Medium_zeros(pa.mediumm.mgrid) for itr=1:Optim.iterations(res)]
			for itr=1:Optim.iterations(res)
				# update mediumm and mediumi
				Seismic_x!(mediumm[itr], mediumi[itr], Optim.x_trace(res)[itr], pa, -1)

				# update gmediumm and gmediumi
				Seismic_gx!(gmediumm[itr],mediumm[itr],gmediumi[itr],mediumi[itr],Optim.trace(res)[itr].metadata["g(x)"],pa,-1)
			end

		else
			f = Optim.minimum(res)
		end
		=#


text:
using DynamicHMCExamples: rel_path
using Literate
using Documenter
using Random
Random.seed!(UInt32[0x57a97f0d, 0x1a38664c, 0x0dddb228, 0x7dbba96f])

####
#### generate using Literate
####

DOCROOT = rel_path("../docs")
DOCSOURCE = joinpath(DOCROOT, "src")
EXAMPLES = ["independent_bernoulli", "linear_regression", "logistic_regression", "multinomial_logistic_regression"]

for example in EXAMPLES
    Literate.markdown(rel_path("example_$(example).jl"), DOCSOURCE)
end

####
#### render & deploy using Documenter
####

makedocs(root = DOCROOT,
         format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
         modules = Module[],
         sitename = "DynamicHMCExamples.jl",
         authors = "Tamás K. Papp",
         pages = vcat(Any["index.md"],
                      Any["example_$(example).md" for example in EXAMPLES]))

deploydocs(root = DOCROOT,
           repo = "github.com/tpapp/DynamicHMCExamples.jl.git")


text:
using ML
using Base.Test

x,y = rand(10),rand(10)

k = GaussianKernel(σ=0.3)
@test k(x,y) ≈ exp(-norm(x-y)^2/(2*k.σ^2))

k = PolynomialKernel(α=10.0,c=1.5,ξ=2.5)
@test k(x,y) ≈ (k.α*dot(x,y) + k.c)^k.ξ

x_train = linspace(0,1,30)[:,:]'
y_train = sin(vec(x_train))
x_test = linspace(0,1,100)[:,:]'
y_test = sin(vec(x_test))

for k in ( GaussianKernel(σ=0.5),
		   PolynomialKernel(α=5.,c=0.1,ξ=4.),
		   ExponentialKernel(σ=10) )
    m = KRR(x_train,y_train,k=k,λ=1e-8)
    fit!(m)
    @test mae(m(x_test),y_test) < 1e-4
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined highlighters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export hl_lt, hl_leq, hl_gt, hl_geq, hl_value

"""
    function hl_lt(n::Number)

Highlight all elements that < `n`.

"""
hl_lt(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(<,data[i,j],n) && data[i,j] < n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :red)
)

"""
    function hl_leq(n::Number)

Highlight all elements that ≤ `n`.

"""
hl_leq(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(≤,data[i,j],n) && data[i,j] ≤ n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :red)
)

"""
    function hl_gt(n::Number)

Highlight all elements that > `n`.

"""
hl_gt(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(>,data[i,j],n) && data[i,j] > n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :blue)
)

"""
    function hl_geq(n::Number)

Highlight all elements that ≥ `n`.

"""
hl_geq(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(≥,data[i,j],n) && data[i,j] ≥ n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :blue)
)

"""
    function hl_value(v::Any)

Highlight all the values that matches `data[i,j] == v`.

"""
hl_value(v) = Highlighter(
    f = (data,i,j)->data[i,j] == v,
    crayon = Crayon(bold       = true,
                    foreground = :yellow)
)


text:
THRESHOLD = 1.2

function fwhm_axis(img, axis, b_mean)
    x, y = size(img)
    if axis == 1
        axis_1 = x
        axis_2 = y
    else
        axis_1 = y
        axis_2 = x
    end
    y_fwhm = []
    y_start = 1
    y_end = axis_2

    for idx = 1:axis_1
        if axis == 1
            temp = img[idx, :]
        else
            temp = img[:, idx]
        end
        if isempty(temp)
            continue
        end
        maxi = maximum(temp)
        half = (maxi + b_mean) / 2
        mid = argmax(temp)
        if half <= b_mean * THRESHOLD # take away lines of background
            continue
        end
        for idy = 1:mid-1
            if temp[idy] < half && temp[idy+1] >= half
                y_start = mid - idy
                break
            end
        end
        for idy = 1:axis_2-mid
            if temp[axis_2-idy] < half && temp[axis_2-idy+1] >= half
                y_end = idy
                break
            end
        end
        temp_fwhm = abs(y_end - y_start)
        append!(y_fwhm, temp_fwhm)
    end
    skew = 0
    if isempty(y_fwhm)
        append!(y_fwhm, 1)
    end
    return (maximum(Float64.(y_fwhm)), skew)
end

function psf_cal(img, b_mean)
    #background = append!(img[1,:],img[end,:])
    #append!(background,img[:,1])
    #append!(background,img[:,end])
    #b_mean = mean(background)

    (y_fwhm, y_skew) = fwhm_axis(img, 1, b_mean)
    (x_fwhm, x_skew) = fwhm_axis(img, 2, b_mean)
    x_fwhm = Float64(x_fwhm)
    y_fwhm = Float64(y_fwhm)
    x_fwhm *= 0.65
    y_fwhm *= 0.65
    return (x_fwhm, y_fwhm, x_skew, y_skew)
end