text:
using Test, Random, Zygote, Distributions, SpecialFunctions, FiniteDifferences
import TruncatedNormal as TruncNorm

for a = -10:10
    d = truncated(Normal(), a, Inf)
    @test TruncNorm.tnmean(a) ≈ mean(d)
    @test TruncNorm.tnstd(a)  ≈ std(d)
    @test TruncNorm.tnvar(a)  ≈ var(d)
    @test a ≤ TruncNorm.tnmean(a) < Inf
    @test 0 ≤ TruncNorm.tnvar(a) ≤ 1
end
@test 1e20 ≤ TruncNorm.tnmean(1e20) < Inf
@test_broken 0 ≤ TruncNorm.tnvar(1e80) ≤ 1

for a = -2:2
    @test mean(TruncNorm.randnt'(a) for _ = 1:10^6) ≈ TruncNorm.tnmean'(a) rtol=0.1
    z, da = TruncNorm.∇randnt(Random.GLOBAL_RNG, a)
    @test da ≈ erfc(z/√2) / erfc(a/√2) * exp((z^2 - a^2)/2)
end

@testset "randnt_half gradient" begin
    for μ = -1:1, σ = 1:2
        dμ_ = mean(gradient(TruncNorm.randnt_half, μ, σ)[1] for _ = 1:10^6)
        dσ_ = mean(gradient(TruncNorm.randnt_half, μ, σ)[2] for _ = 1:10^6)
        dμ, dσ = gradient(μ, σ) do μ, σ
            μ + σ * TruncNorm.tnmean(-μ/σ)
        end
        @test dμ ≈ dμ_ rtol=0.1
        @test dσ ≈ dσ_ rtol=0.1
    end
end

@testset "sqrt1half" begin
    @test TruncNorm.sqrt1half(5) ≈ 5.1925824035672520156
    @test TruncNorm.sqrt1half(0) == 1
    @test TruncNorm.sqrt1half(-1) == TruncNorm.sqrt1half(1) ≈ 1.6180339887498948482
    @test isnan(TruncNorm.sqrt1half(NaN))
    @test TruncNorm.sqrt1half(Inf) == TruncNorm.sqrt1half(-Inf) == Inf
    @test TruncNorm.sqrt1half(1e300) ≈ 1e300
end

@testset "randnt" begin
    @test TruncNorm.randnt(0) > 0
    @test TruncNorm.randnt(1e300) == 1e300
    @test TruncNorm.randnt(Inf) == Inf
    @test isnan(TruncNorm.randnt(NaN))
    @test TruncNorm.randnt(floatmax(Float64)) == floatmax(Float64)
end


text:
using Test
using BundleMethod
using JuMP
using Ipopt
using Random

const BM = BundleMethod

@testset "Abstract Method" begin

    mutable struct BaseMethod <: BM.AbstractMethod
        model::BM.BundleModel
        function BaseMethod()
            return new(BM.BundleModel(1, 1, 1, nothing))
        end
    end

    bm = BaseMethod()
    empty_model = BM.get_model(bm)
    @test empty_model.n == 0
    @test empty_model.N == 0

    BM.get_model(method::BaseMethod) = method.model

    BM.set_user_data(BM.get_model(bm), nothing)

    BM.set_optimizer(bm, Ipopt.Optimizer)
    BM.get_solution(bm)
    @test BM.get_objective_value(bm) == Inf
    BM.set_bundle_tolerance!(bm, 0.0)

    BM.build_bundle_model!(bm)
    model = BM.get_jump_model(bm)
    @test JuMP.num_variables(model) == 2
    @test BM.termination_test(bm) == true

    BM.add_constraints!(bm)
    BM.collect_model_solution!(bm)
    BM.evaluate_functions!(bm)
    BM.update_bundles!(bm)
    BM.purge_bundles!(bm)
    BM.add_bundles!(bm)
    BM.update_iteration!(bm)
    BM.display_info!(bm)
end

@testset "Proximal Method" begin

    include("../examples/simple.jl")

    vm = JuMP.Model(Ipopt.Optimizer)
    @variable(vm, -1 <= vm_x[j=1:n] <= 1)
    @objective(vm, Min, sum(b[i] * (vm_x[j] - a[i,j])^2 for i = 1:N, j = 1:n))
    optimize!(vm)
    objval = JuMP.objective_value(vm)
    xval = Dict{Int,Float64}()
    for j in 1:n
        xval[j] = JuMP.value(vm_x[j])
    end

    @show BM.get_objective_value(pm)
    @show BM.get_solution(pm)
    @test isapprox(objval, BM.get_objective_value(pm), rtol=1e-2)
    for j in 1:n
        @test isapprox(xval[j], BM.get_solution(pm)[j], rtol=1e-2)
    end

    params = BM.Parameters()
    BM.set_parameter(params, "maxiter", 3)
    BM.set_parameter(params, "ϵ_s", 1.e-6)

    pm2 = BM.ProximalMethod(n, N, evaluate_f, params = params)

    # Set optimization solver to the internal JuMP.Model
    model2 = BM.get_jump_model(pm2)
    set_optimizer(model2, Ipopt.Optimizer)
    set_optimizer_attribute(model2, "print_level", 0)

    BM.build_bundle_model!(pm2)
    BM.run!(pm2)
end

@testset "Trust Region Method" begin

    include("../examples/tr_simple.jl")
    vm = JuMP.Model(Ipopt.Optimizer)
    @variable(vm, -1 <= vm_x[j=1:n] <= 1)
    @objective(vm, Min, sum(b[i] * (vm_x[j] - a[i,j])^2 for i = 1:N, j = 1:n))
    optimize!(vm)
    objval = JuMP.objective_value(vm)
    xval = Dict{Int,Float64}()
    for j in 1:n
        xval[j] = JuMP.value(vm_x[j])
    end

    @show BM.get_objective_value(pm)
    @show BM.get_solution(pm)
    @test isapprox(objval, BM.get_objective_value(pm), rtol=1e-2)
    for j in 1:n
        @test isapprox(xval[j], BM.get_solution(pm)[j], rtol=1e-2)
    end
end

ext: jl


text:
"""
Prinz bursting neuron script
"""

using NeuronBuilder, ModelingToolkit, OrdinaryDiffEq, Plots

const Area = 0.0628 # Prinz/Liu 0.0628 mm2
const Cm = 10.0 # specific capacitance cₘ is a biological constant (around) 10 nF/mm^2

Prinz_conv = Cm / Area

channels = [Prinz.Na(100 * Prinz_conv), Prinz.CaS(6 * Prinz_conv), Prinz.CaT(2.5 * Prinz_conv), Prinz.H(0.01 * Prinz_conv),
    Prinz.Ka(50 * Prinz_conv), Prinz.KCa(5 * Prinz_conv), Prinz.Kdr(100 * Prinz_conv), Prinz.leak(0.0)]

τCa = 200.0
Ca∞ = 0.05
fxarea = 14.96 * 0.0628

defaults = Dict(Voltage => BasicVoltageDynamics(),
    Calcium => Prinz.CalciumDynamics(τCa, Ca∞, fxarea),
    Reversal{Calcium} => Prinz.CaReversalDynamics())

somatic_parameters = Dict(
    Reversal{Sodium} => 50.0,
    Reversal{Potassium} => -80.0,
    Reversal{Leak} => -50.0,
    Reversal{Proton} => -20.0,
    Voltage => -60.0,
    Calcium => 0.05,
    Reversal{Calcium} => 0.0)

b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, channels, :test_Prinz)

neur = b()

prob = ODEProblem(neur, [], (0.0, 5000.0), [])

sol = solve(prob, AutoTsit5(Rosenbrock23()))
plot(sol)


text:
using Statistics
export prepareSyntheticDipDataFiles
function prepareSyntheticDipDataFiles(m::Array{Float64},Minv::RegularMesh,filename::String,pad::Int64,theta_phi_rad::Array{Float64,2},samplingBinning::Int64, noiseSample::Float64, noiseAngles::Float64,noiseTrans::Float64,method = MATFree)

n = Minv.n;
pFor = getDipParam(Minv,theta_phi_rad + noiseAngles*randn(size(theta_phi_rad)),noiseTrans*randn(size(theta_phi_rad,1),3),samplingBinning,nworkers(),method);
@time Dremote,pFor = getData(m[:],pFor);

Dlocal = arrangeRemoteCallDataIntoLocalData(Dremote);
Data = Dlocal + noiseSample*mean(Dlocal[:])*randn(size(Dlocal));

pFor = getDipParam(Minv,theta_phi_rad + zeros(size(theta_phi_rad)),zeros(size(theta_phi_rad,1),3),samplingBinning,nworkers(),method);
Dremote,pFor = getData(m[:],pFor);

Dlocal = arrangeRemoteCallDataIntoLocalData(Dremote);
DataClean = Dlocal;

file = matopen(string(filename,"_data.mat"), "w");
write(file,"domain",Minv.domain);
write(file,"noiseSample",noiseSample);
write(file,"noiseAngles",noiseAngles);
write(file,"n",Minv.n);
write(file,"samplingBinning",samplingBinning);
write(file,"Data",Data);
write(file,"DataClean",DataClean);
write(file,"theta_phi_rad",theta_phi_rad);
# write(file,"m_true",m);
write(file,"pad",pad);
close(file);
return;
end

export readDipDataFile
function readDipDataFile(filename::String)
file = matopen(string(filename,"_data.mat"), "r");
domain = read(file,"domain");
n = read(file,"n");
samplingBinning = read(file,"samplingBinning");
Data = read(file,"Data");
# m_true = read(file,"m_true")
pad = read(file,"pad");
theta_phi_rad = read(file,"theta_phi_rad");
close(file);
return n,domain,samplingBinning,Data,theta_phi_rad,pad;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
# Example Configuration for nearly neutral infinite sites simuation
# Run with small values of N, ngens for testing
# Command line to run:  julia run.jl examples/is_example0
# Creates output file examples/in_example0.csv

const nn_simtype = 2    # nearly neutral infinite sites
#type_str = "N"
mu_list_flag =true
@everywhere const N_list     = [5]
const average_by_generations = true
const mu_list= [0.3]   # Mutation rate
const L = 1           # number of loci # not used in infinite alleles
const ngens = 3      # Generations after burn-in
const burn_in= 0.0    # generations of burn_in as a multiple of N
dfe=dfe_mixed
dfe_str = "dfe_mixed"

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
## metadata

function init_metadata(sock::TCPSocket)
    header = RequestHeader(METADATA, 0, 0, DEFAULT_ID)
    req = AllTopicsMetadataRequest()
    send_request(sock, header, req)
    header = recv_header(sock)
    resp = readobj(sock, TopicMetadataResponse)
    return process_response(resp)
end

function _metadata(kc::KafkaClient, topics::Vector{S}) where S<:AbstractString
    cor_id, ch = register_request(kc, TopicMetadataResponse)
    header = RequestHeader(METADATA, 0, cor_id, DEFAULT_ID)
    req = TopicMetadataRequest(topics)
    send_request(random_broker(kc), header, req)
    return ch
end

function metadata(kc::KafkaClient, topics::Vector{S}) where {S<:AbstractString}
    ch = _metadata(kc, topics)
    return map(process_response, ch)
end

function _metadata(kc::KafkaClient)
    cor_id, ch = register_request(kc, TopicMetadataResponse)
    header = RequestHeader(METADATA, 0, cor_id, DEFAULT_ID)
    req = AllTopicsMetadataRequest()
    send_request(random_broker(kc), header, req)
    return ch
end

function metadata(kc::KafkaClient)
    ch = _metadata(kc)
    return map(process_response, ch)
end

function parse_metadata(md::PartitionMetadata)
    if md.partition_error_code != 0
        error("TopicMetadataRequest failed with error code " *
              "$(md.partition_error_code)")
    end
    return Dict(md.partition_id =>
                Dict(:leader => md.leader, :replicas => md.replicas,
                     :isr => md.isr))
end

function parse_metadata(md::TopicMetadata)
    if md.topic_error_code != 0
        error("TopicMetadataRequest failed with error code " *
              "$(md.topic_error_code)")
    end
    partition_metadata_dict = merge(map(parse_metadata, md.partition_metadata)...)
    return Dict(md.topic_name => partition_metadata_dict)
end

function process_response(resp::TopicMetadataResponse)
    if !isempty(resp.topic_metadata)
        topic_metadata_dict = merge(map(parse_metadata, resp.topic_metadata)...)
    else
        topic_metadata_dict = Dict{Symbol,Any}()
    end
    return Dict(:brokers => resp.brokers, :topics => topic_metadata_dict)
end

## produce

# using Message version 0
function make_message(key::Vector{UInt8}, value::Vector{UInt8})
    magic_byte = Int8(0)
    # magic_byte = Int8(1)  # current version of Kafka message binary format
    attributes = Int8(0)  # no compression
    buf = IOBuffer()
    writeobj(buf, magic_byte)
    writeobj(buf, attributes)
    writeobj(buf, key)
    writeobj(buf, value)
    crc = reinterpret(Int32, crc32(take!(buf)))
    return Message(Int32(crc), magic_byte, attributes, key, value)
end

function make_produce_request(topic::AbstractString, partition_id::Integer,
                              kvs::Vector{Tuple{Vector{UInt8}, Vector{UInt8}}})
    partition = Int32(partition_id)
    messages = [make_message(k, v) for (k, v) in kvs]
    message_set_elems = [OffsetMessage(0, obj_size(msg), msg)
                         for msg in messages]
    message_set_size = reduce(+, map(obj_size, message_set_elems))
    partition_data = ProduceRequestPartitionData(partition, message_set_size,
                                   MessageSet(message_set_elems))
    topic_data = ProduceRequestTopicData(topic, [partition_data])
    required_acks = 1 # broker will send response after writing to a local log
    timeout = 10_000  # in milliseconds
    req = ProduceRequest(required_acks, timeout, [topic_data])
    return req
end

function _produce(kc::KafkaClient, topic::AbstractString, partition::Integer,
                  kvs::Vector{Tuple{Vector{UInt8}, Vector{UInt8}}})
    pid = Int32(partition)
    cor_id, ch = register_request(kc, ProduceResponse)
    header = RequestHeader(PRODUCE, 0, cor_id, DEFAULT_ID)
    req = make_produce_request(topic, pid, kvs)
    send_request(find_leader(kc, topic, pid), header, req)
    return ch
end

function produce(kc::KafkaClient, topic::AbstractString, partition::Integer,
                 kvs::Vector{Tuple{Vector{UInt8}, Vector{UInt8}}})
    ensure_leader(kc, topic, partition)
    ch = _produce(kc, topic, partition, kvs)
    return map(process_response, ch)
end

function process_response(resp::ProduceResponse)
    partition, error_code, offset = resp.responses[1][2][1]
    if error_code != 0
        error("ProduceRequest failed with error code: $error_code")
    end
    return offset
end

## fetch

function make_fetch_request(topic::AbstractString, prt_id::Integer, offset::Int64,
                            max_wait_time::Int, min_bytes::Int, max_bytes::Int)
    partition = Int32(prt_id)
    partition_fetch = FetchRequestPartitionData(partition, offset, max_bytes)
    topic_fetch = FetchRequestTopicData(topic, [partition_fetch])
    req = FetchRequest(-1, max_wait_time, min_bytes, [topic_fetch])
    return req
end

function _fetch(kc::KafkaClient, topic::AbstractString,
                prt_id::Integer, offset::Int64;
                max_wait_time=100, min_bytes=1024, max_bytes=1024*1024)
    partition = Int32(prt_id)
    cor_id, ch = register_request(kc, FetchResponse)
    header = RequestHeader(FETCH, 0, cor_id, DEFAULT_ID)
    req = make_fetch_request(topic, partition, offset,
                             max_wait_time, min_bytes, max_bytes)
    send_request(find_leader(kc, topic, partition), header, req)
    return ch
end

function fetch(kc::KafkaClient, topic::AbstractString,
               partition::Integer, offset::Int64;
               max_wait_time=100, min_bytes=1024, max_bytes=1024*1024,
               key_type=Vector{UInt8}, value_type=Vector{UInt8})
    ensure_leader(kc, topic, partition)
    ch = _fetch(kc, topic, partition, offset; max_wait_time=max_wait_time,
                min_bytes=min_bytes, max_bytes=max_bytes)
    return map(r -> process_response(r, key_type, value_type), ch)
end

function process_response(resp::FetchResponse)
    pr = resp.topic_results[1].partition_results[1]
    if pr.error_code != 0
        error("FetchRequest failed with error code: $(pr.error_code)")
    end
    elements = pr.message_set.elements
    results = [(elem.offset, elem.message.key, elem.message.value)
               for elem in elements]
    return results
end

function process_response(resp::FetchResponse, ::Type{K}, ::Type{V}) where {K,V}
    return [(offset, convert(K, key), convert(V, value))
            for (offset, key, value) in process_response(resp)]
end

## offset listing

function make_offset_request(replica_id::Int32,
                             topic::AbstractString, partition::Int32,
                             time::Int64, max_number_of_offsets::Int64)
    partition_data = OffsetRequestPartitionData(partition, time,
                                                max_number_of_offsets)
    topic_data = OffsetRequestTopicData(topic, [partition_data])
    req = OffsetRequest(replica_id, [topic_data])
    return req
end

function _list_offsets(kc::KafkaClient, topic::AbstractString, partition::Integer,
                       time::Int64, max_number_of_offsets::Int64)
    pid = Int32(partition)
    cor_id, ch = register_request(kc, OffsetResponse)
    header = RequestHeader(OFFSETS, 0, cor_id, DEFAULT_ID)
    req = make_offset_request(find_leader_id(kc, topic, partition),
                              topic, pid, time, max_number_of_offsets)
    send_request(find_leader(kc, topic, partition), header, req)
    return ch
end

function list_offsets(kc::KafkaClient, topic::AbstractString, partition::Integer,
                      time::Int64, max_number_of_offsets::Int64)
    ensure_leader(kc, topic, partition)
    ch = _list_offsets(kc, topic, partition, time, max_number_of_offsets)
    return map(process_response, ch)
end

function earliest_offset(kc::KafkaClient,
                         topic::AbstractString, partition::Integer)
    ch = list_offsets(kc, topic, partition, -2, 1)
    return map(offsets -> offsets[1], ch)
end

function latest_offset(kc::KafkaClient,
                       topic::AbstractString, partition::Integer)
    ch = list_offsets(kc, topic, partition, -1, 1)
    return map(offsets -> offsets[1], ch)
end

function process_response(resp::OffsetResponse)
    pd = resp.topic_data[1].partition_data[1]
    if pd.error_code != 0
        error("OffsetRequest failed with error code: $(pd.error_code)")
    end
    return pd.offsets
end

## API versions

function api_versions(kc::KafkaClient)
    sock = random_broker(kc)
    cor_id, ch = register_request(kc, ApiVersionsResponse)
    header = RequestHeader(API_VERSIONS, 1, cor_id, DEFAULT_ID)
    send_request(sock, header)
    resp = take!(ch)
    return resp.api_versions
end


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
# %%
module Laser
# modules
using Parameters
using AtomicUnits
export Pulse, Time, getEt, getAt, getEtAt, get_env
# Et => give E(t) and A(t) = -\int dt' E(t')
# At => give A(t) and E(t) = -dA/dt
abstract type FormType end
struct Et <: FormType end
struct At <: FormType end
const dict_formtype = Dict{Symbol,T where T <: FormType}(:Et => Et(), :At => At())
struct Pulse{TE <: FormType,TF <: AbstractFloat,TI <: Integer}
    form::TE
    I::TF
    E::TF
    A::TF
    λ::TF
    ω::TF
    T::TF
    ϕ::TF
    env_type::Symbol
    ncyc::TI
    tL::TF
    t0::TF
    tR::TF
    tprev::Vector{TF}
    Atprev::Vector{TF}
end
function Pulse(fint, wlen, cep, env_type, ncyc_env, form_type::Symbol=:Et)
    @assert haskey(dict_formtype, form_type)
    I = float(fint)
    E = i2e(fint)
    λ = float(wlen)
    ω = wlen2au(wlen)
    A = E / ω
    T = 2π / ω
    ϕ = deg2rad(float(cep))
    tL = 0.0
    t0 = T * ncyc_env / 2
    tR = T * ncyc_env
    tprev = [0.0]
    Atprev = [0.0]
    return Pulse(
        dict_formtype[form_type],
        I,
        E,
        A,
        λ,
        ω,
        T,
        ϕ,
        env_type,
        ncyc_env,
        tL,
        t0,
        tR,
        tprev,
        Atprev,
    )
end
include("getEtAt.jl")
include("time.jl")
include("flux.jl")
end
# %%


text:
using Tables

abstract type AbstractConnector end

Base.length(conn::AbstractConnector) = Inf
hasnext(conn::AbstractConnector) = true
reset!(conn::AbstractConnector) = nothing

mutable struct TablesConnector <: AbstractConnector
    rows
    state::Int
    args::Dict{Symbol, Any}
end

function TablesConnector(data; shuffle::Bool = false)
    if !Tables.istable(data)
        throw(ArgumentError("data must have the Tables.jl interface"))
    end

    if shuffle
        data = data[Random.shuffle(1:size(data,1)), :]
    end

    return TablesConnector(Tables.rows(data), 0, Dict{Symbol, Any}())
end

function TablesConnector(data, orderBy::Symbol; rev::Bool = false)
    if !(orderBy in propertynames(data))
        throw(ArgumentError("data doesn't have the column $orderBy"))
    end

    data = sort(data, orderBy, rev = rev)

    return TablesConnector(data)
end

TablesConnector(filename::String) = TablesConnector(CSV.read(filename; header = false))

Base.length(conn::TablesConnector) = length(conn.rows)
hasnext(conn::TablesConnector) = conn.state < length(conn)

function next(conn::TablesConnector)
    if conn.state >= length(conn)
        return nothing
    end

    conn.state += 1

    return DataFrame([conn.rows[conn.state]])
end

reset!(conn::TablesConnector) = conn.state = 0

mutable struct GeneratorConnector <: AbstractConnector
    generator::Function
    args::Dict{Symbol, Any}
end

function GeneratorConnector(generator::Function;
                            args...)
    return GeneratorConnector(generator, args)
end

function next(conn::GeneratorConnector)
    total = 100
    data = conn.generator(;n_samples = total, conn.args...)

    return DataFrame(data[1 + Int(floor(rand(1,1)[1] .* size(data)[1])), :])
end


text:
using Test
using LevelDB

run(`rm -rf level.db`)
run(`rm -rf level.db.2`)
run(`rm -rf level.db.3`)

@testset "DB basic operations" begin
    @test_throws ErrorException LevelDB.DB("level.db")

    db = LevelDB.DB("level.db", create_if_missing = true)
    close(db)
    @test !isopen(db)

    @test_throws ErrorException LevelDB.DB("level.db", error_if_exists = true)
    db = LevelDB.DB("level.db")

    db[[0x00]] = [0x01]
    @test db[[0x00]] == [0x01]
    delete!(db, [0x00])
    # as in Dict, deleting a non-existing key should not throw an error
    delete!(db, [0x00])
    @test_throws KeyError db[[0x00]]
    close(db)
    @test db.handle        == C_NULL
    @test db.options       == C_NULL
    @test db.write_options == C_NULL
    @test db.read_options  == C_NULL
end

@testset "DB batching and iteration" begin
    db = LevelDB.DB("level.db.2", create_if_missing = true)
    d = Dict([0xa] => [0x1],
             [0xb, 0xb] => [0x2],
             [0xc, 0xc] => [0x3],
             [0xd] => [0x4],
             [0xe] => [0x5])

    db[keys(d)] = values(d)

    @test db[[0xa]] == [0x1]
    @test db[[0xb, 0xb]] == [0x2]
    @test db[[0xc, 0xc]] == [0x3]
    @test db[[0xd]] == [0x4]
    @test db[[0xe]] == [0x5]

    function size_of(db)
        i = 0
        for (k, v) in db
            i += 1
        end
        i
    end
    @test size_of(db) == 5

    d2 = Dict{Vector{UInt8}, Vector{UInt8}}()
    for (k, v) in db
        d2[k] = v
    end
    @test d == d2

    for (k, v) in db
        delete!(db, k)
    end
    @test size_of(db) == 0

    close(db)
    @test db.handle        == C_NULL
    @test db.options       == C_NULL
    @test db.write_options == C_NULL
    @test db.read_options  == C_NULL

    # nothing should happen here
    close(db)
    @test db.handle        == C_NULL
    @test db.options       == C_NULL
    @test db.write_options == C_NULL
    @test db.read_options  == C_NULL
end

@testset "DB Errors" begin
    @test_throws ErrorException LevelDB.DB("level.db.3")
end

@testset "DB key range iterator" begin
    db = LevelDB.DB("level.db", create_if_missing = true)
    d = Dict([0xa] => [0x1],
             [0xb, 0xb] => [0x2],
             [0xc, 0xc] => [0x3],
             [0xd] => [0x4],
             [0xe] => [0x5],
             )

    db[keys(d)] = values(d)
    n = 0
    # Test a range iterator that begins with the second item in DB
    iter = LevelDB.RangeView(db, [0xb, 0xb], [0xd])
    for (k, v) in iter
        if n == 0
            @test v == [0x2]
        elseif n == 1
            @test v ==  [0x3]
        elseif n == 2
            @test v ==  [0x4]
        end
        n +=1
    end
    @test n == 3 # assert the iterator stopped at key 0xd

    # Test a range iterator that begins with a non-existing  item in DB
    n = 0
    iter = LevelDB.RangeView(db, [0xb, 0xa], [0xc, 0xd])
    for (k, v) in iter
        if n == 0
            @test v == [0x2]
        elseif n == 1
            @test v ==  [0x3]
        end
        n +=1
    end
    @test n == 2  # assert the iterator stopped at key 0xc 0xc

    close(db)

end

run(`rm -rf level.db`)
run(`rm -rf level.db.2`)
run(`rm -rf level.db.3`)


text:
## test the main entry point in Celeste: the `infer` function
import JLD

"""
test infer with a single (run, camcol, field).
This is basically just to make sure it runs at all.
"""
function test_infer_single()
    # very small patch of sky that turns out to have 4 sources.
    # We checked that this patch is in the given field.
    box = ParallelRun.BoundingBox(164.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]
    result = ParallelRun.one_node_infer(field_triplets, datadir; box=box)
end

function test_source_division_parallelism()
    box = ParallelRun.BoundingBox(164.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]
    results = ParallelRun.divide_sources_and_infer(box, datadir)
    @test length(results) == 4
end

test_source_division_parallelism()
test_infer_single()


text:
#=
    This file includes all the funtions needed to update the operation
    and investment informations
 =#

function update_operation_informations!(h::Int64, y::Int64, s::Int64, mg::Microgrid, ω::Scenarios)
    # Demands
    for (k, a) in enumerate(mg.demands)
        a.timestamp[h,y,s] = ω.demands[k].t[h,y,s]
        a.carrier.power[h,y,s] = ω.demands[k].power[h,y,s]
    end
    # Generations
    for (k, a) in enumerate(mg.generations)
        a.timestamp[h,y,s] = ω.generations[k].t[h,y,s]
        a.carrier.power[h,y,s] = a.powerMax[y,s] * ω.generations[k].power[h,y,s]
    end
    # Grids - We assume the price of electricity is known over the year
    for (k, a) in enumerate(mg.grids)
        if h == 1
            a.cost_in[1:end,y,s] = ω.grids[k].cost_in[1:end,y,s]
            a.cost_out[1:end,y,s] = ω.grids[k].cost_out[1:end,y,s]
        end
    end
end

function update_investment_informations!(y::Int64, s::Int64, mg::Microgrid, ω::Scenarios)
    # Generations
    for (k, a) in enumerate(mg.generations)
        a.cost[y,s] = ω.generations[k].cost[y,s]
    end
    # Storages
    for (k, a) in enumerate(mg.storages)
        a.cost[y,s] = ω.storages[k].cost[y,s]
    end
    # Converters
    for (k, a) in enumerate(mg.converters)
        a.cost[y,s] = ω.converters[k].cost[y,s]
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2016/17")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Test

rate = (u,p,t) -> (0.1/1000.0)*u[1]*u[2]
affect! = function (integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)

rate = (u,p,t) -> 0.01u[2]
affect! = function (integrator)
  integrator.u[2] -= 1
  integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,affect!)

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
integrator = init(jump_prob,FunctionMap())

condition(u,t,integrator) = t == 100
function purge_affect!(integrator)
  integrator.u[2] ÷= 10
  reset_aggregated_jumps!(integrator)
end
cb = DiscreteCallback(condition,purge_affect!,save_positions=(false,false))
sol = solve(jump_prob,FunctionMap(),callback=cb,tstops=[100])
sol = solve(jump_prob,SSAStepper(),callback=cb,tstops=[100])


text:
"""
A function to uniformly select a random agent from a collection
"""
function selectRandomAgent(agents)
    agents[rand(1:end)]
end

"""
A function that unifromly select an agent from a subgroup of agents that have
a minimum amount of tokens. It could be useful to select an agent that is capable
of making a deposit
"""
function selectRandomAgentWithMinBalance(agents, min)
    filteredAgents = filter(a -> a.balance > min, agents)
    selectRandomAgent(filteredAgents)
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
module VideoPlayer

using VideoIO

abstract type AbstractVideoBackend end

include(joinpath(@__DIR__, "..", "deps", "includes.jl"))

play(args...; kwargs...) = play(videobackend(), args...; kwargs...)

export videobackend, play

end # module


text:
using Test
include("../src/fdfd_local_field.jl")

@testset "fdfd_local_field" begin
    @testset "maxwell_fdfd_solver" begin include("test_fdfd_maxwell_solver.jl") end
    @testset "pixel averaging" begin include("test_geometry_code.jl") end
    @testset "finite difference gradient" begin include("test_get_gradient.jl") end
    @testset "parallelization" begin include("test_parallel_vs_serial.jl") end
end


text:
# tests for the orientation code

# packages
using PyPlot
using Debug

# local
using Kernel

include("../src/randvonMisesFisher.jl")

#n = 200
n = 1
v1 = randvonMisesFisher(3,n,.1,[1 1 1])
v2 = randvonMisesFisher(3,n,.1,[1 1 1])
v3 = randvonMisesFisher(3,n,.1,[1 1 1])

v = [v1;v2;v3]

# make sure they are all unit vectors
for c=1:n
	for r=1:3
		nrm = norm( vec( v[(r-1)*3+1:(r-1)*3+3,c] ) )
		v[(r-1)*3+1:(r-1)*3+3,c] /= nrm
	end
end

v = v'

println("size(v) = ", size(v))

#lm = LogMap(v',[1 1 1])

#em = ExpMap(lm,[0.8 0.8 0.8])

#K = SphericalCosKernel(v,v)
Kp = SphericalKernel(v,v)

#println("fro: ", vecnorm(K-Kp))
#println("l2: ", norm(K-Kp))

Kp2 = SphericalKernelLowMem(v,v)

println("fro: ", vecnorm(Kp-Kp2))
println("l2: ", norm(Kp-Kp2))

Cov = SphericalCovariance(v)

println("Spherical Feature and Inverse: ")

println("v: ", v)
Phi = SphericalFeature(v)
println("Phi: ", Phi)

X = SphericalFeatureInverse(Phi)
println("v': ", X)

println("fro: ", vecnorm(v-X))
println("l2: ", norm(v-X))


text:
# Autogenerated wrapper script for libcint_jll for x86_64-apple-darwin
export libcint

using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("libcint")
JLLWrappers.@declare_library_product(libcint, "@rpath/libcint.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libcint,
        "lib/libcint.5.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# stc_load.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.6.0
# 16.12.2016
# Last Edit: 07.07.2021

"""# stc_load

Load data from Aerodyne STC files generated by TDLWintel

---

### Examples

`time,data = stc_load(\"K:\\\\Data\\\\140113_135354.stc\")`\n
`time,data = stc_load([\"K:\\\\Data\\\\140113_135354.stc\",\"K:\\\\Data\\\\140114_000000.stc\"])`\n
`time,data = stc_load(\"K:\\\\Data\\\\\")` # Load all STC files in the given directory and subdirectories\n
`time,data = stc_load(\"K:\\\\Data\\\\\",verbose=true,cols=[\"Traw\",\"X1\"])` # Load only Traw and X1 from all files in the given directory and display information as it is processed\n
`time,data = stc_load(\"K:\\\\Data\\\\\",DateTime(2014,6,27,11,58,0))` # Load files starting at the given timestamp\n
`time,data = stc_load(\"K:\\\\Data\\\\\",DateTime(2014,6,27,11,58,0),DateTime(2014,7,4,12,0,0))` # Load files between the given timestamps\n

---

`time, D = stc_load(F)`\n
* **time**::Array{DateTime,1} = Parsed time values from the STC file
* **D**::DataFrame = Data from the STC file
* **F**::Array{String,1} = Array of STC files to load

`time, D = stc_load(F;verbose=false,cols=[])`\n
* **time**::Array{DateTime,1} =
* **F**::String = File name and location
* **verbose**::Bool (optional) = Display what is happening, FALSE is default
* **cols**::Array{String,1} or Array{Symbol,1} (optional) = List of columns to return, [] is default

`time, D = stc_load(F,mindate;verbose=false,cols=[])`\n
`time, D = stc_load(F,mindate,maxdate;verbose=false,cols=[])`\n
* **mindate**::DateTime = Load files including and after this date and time
* **maxdate**::DateTime = Load files up to this date and time
"""
function stc_load(Dr::String,mindate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)
	## Load STC files in the given directory (Dr) including and beyond the given date ##
	return stc_load(Dr,mindate,DateTime(9999,1,1,0,0,0);verbose=verbose,cols=cols,recur=recur)
end # End of stc_load(Dr::String,mindate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)

function stc_load(Dr::String,mindate::DateTime,maxdate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)
	## Load STC files in the given directory (Dr) between the given dates ##

	# Check for Directory
	if isdir(Dr) == false
		error("First input should be a directory")
	end

	# List Files
	verbose ? println("Listing Files") : nothing
	(Fstc,folders) = dirlist(Dr,regex=r"\d{6}_\d{6}\.stc$",recur=recur) # List STC files

	# Parse Times
	(Tstc,Fstc) = aerodyne_parsetime(Fstc)

	# Sort Info
	f1 = findlast(Tstc .<= mindate)
	if isnothing(f1)
		f1 = 1
	end
	f2 = findlast(Tstc .< maxdate)
	f = collect(f1:f2)
	Tstc = Tstc[f]
	Fstc = Fstc[f]

	# Remove Files Out of Range
	begin
		# STC Files
		f = findall(mindate .<= Tstc .< maxdate)
		Tstc = Tstc[f]
		Fstc = Fstc[f]
	end

	# Load Data
	verbose ? println("Loading:") : nothing
	(Tstc,Dstc) = stc_load(Fstc,verbose=verbose,cols=cols)

	# Remove Time Values Out of Range
	begin
		# Ensure all rows are defined in SPEFile
		if in(true,[isequal(:SPEFile,i) for i in names(Dstc)]) # See if :SPEFile exists first
			f = ismissing.(Dstc.SPEFile)
			Dstc.SPEFile[f] .= ""
		end

		f = findall(mindate .<= Tstc .< maxdate)
		Tstc = Tstc[f]
		Dstc = Dstc[f,:]
	end

	return Tstc,Dstc
end # End of stc_load(Dr::String,mindate::DateTime,maxdate::DateTime;verbose::Bool=false,cols::Array{String,1}=String[],recur::Int=9999)

function stc_load(F::Array{T,1};verbose::Bool=false,cols::Array{String,1}=String[]) where T <: String
	## Load a list of file names ##

	# Sort files by dates
	(Ftime,F) = aerodyne_parsetime(F)

	# Load and concatinate data
	if !isempty(F)
		(t,D) = stc_load(F[1],verbose=verbose,cols=cols) # Initial load
		for i=2:1:length(F)
			(tempT,tempD) = stc_load(F[i],verbose=verbose,cols=cols)
			t = [t;tempT]
			D = [D;tempD]
		end
	else
		# No Files Found
		t = Vector{DateTime}(undef,0)
		D = DataFrame()
	end

	return t,D
end # End stc_load(F::Array{T,1};verbose::Bool=false,cols::Array{String,1}=String[]) where T <: String

function stc_load(F::String;verbose::Bool=false,cols::Array{String,1}=String[])
	ext = splitext(F)[2][2:end] # Save the extension

	verbose ? println("  " * F) : nothing

	## Check for a proper file
	if isdir(F) == true
		return stc_load(F,DateTime(0,1,1,0,0,0),DateTime(9999,1,1,0,0,0),verbose=verbose,cols=cols)
	elseif isempty(ext) == true
		error("No file extension")
	elseif ext != "stc"
		error("Extension is not STC. Returning nothing...")
		return
	end

	#################
	##  Load Data  ##
	#################
	## Load Header Information
	fid = open(F,"r")
	h1 = readline(fid)
	h2 = readline(fid)

	h3 = permutedims(readdlm(IOBuffer("\"" * replace(h2[1:end],"," => "\",\"") * "\""),','),[2,1])
	h = Array{String}(undef,length(h3))
	for i=1:1:length(h3)
		h[i] = strip(String(h3[i])) # Remove leading and trailing whitespace
		h[i] = replace(h[i]," " => "_") # Replace remaining whitespace with _
	end
	close(fid)

	## Check for duplicate column names and adjust them
	unames = unique(h) # unique column names
	for i in unames
		check = 0
		for j = 1:length(h)
			if i == h[j]
				check += 1

				if check > 1
					h[j] = h[j] * "-" * string(check)
				end
			end
		end
	end

	## List column types
	coltypes = Any[Float64 for i=1:length(h)]
	for i=1:1:length(h)
		if h[i] == "SPEFile"
			coltypes[i] = Union{Missing,String}
		end

		if h[i] == "StatusW"
			coltypes[i] = Int
		end
	end

	## Todo: Replace column names with reasonable names

	## Load data
	D = DataFrame()
	try
		D = CSV.read(F,DataFrame;delim=",",types=coltypes,header=h,datarow=3,missingstring="-")
	catch
		println("Cannot load " * F)
		error("ERROR loading file")
	end

	#########################
	##  Parse Time Format  ##
	#########################
	time = Array{DateTime}(undef,length(D.time)) # Preallocate time column
	secs = Dates.Second # Initialize so it doesn't have to do it every time in the loop
	millisecs = Dates.Millisecond # Initialize so it doesn't have to do it every time in the loop
	for i=1:1:length(D.time)
		secs = Dates.Second(floor(D.time[i]))
		millisecs = Dates.Millisecond(floor(1000(D.time[i] - floor(D.time[i]))))
		time[i] = DateTime(1904,1,1,0,0,0) + secs + millisecs
	end

	##################################
	##  Keep only selected columns  ##
	##################################
	if !isempty(cols)
		# Check cols' type
		if typeof(cols) != Array{Symbol,1} && typeof(cols) != Symbol
			temp = Array{Symbol}(undef,length(cols))
			for i=1:1:length(cols)
				temp[i] = Symbol(cols[i]) # Convert all the values to symbols
			end
			cols = temp
		end

		# Make sure each entry exists
		fields = Symbol.(names(D))
		cols_bool = fill!(Array{Bool}(undef,length(cols)),false) # Preallocate false
		for i=1:1:length(cols)
			for j=1:1:length(fields)
				if fields[j] == cols[i]
					cols_bool[i] = true
				end
			end
		end
		cols = cols[cols_bool]

		# Remove Unwanted column
		if isempty(cols)
			D = DataFrame()
			time = Array{DateTime}(undef,0)
		else
			D = D[:,cols]
		end
	end

	return time,D
end # End of stc_load(F::String;verbose::Bool=false,cols::Array{String,1}=String[])


text:
test_dir = dirname(@__FILE__)
tests = [
  "boolean",
  "code",
  "configuration",
  "exec",
  "float",
  "integer",
  "name",
  "stack"
]

println("Running tests:")

for t in tests
  f = "test_$t.jl"
  println("  * $f...")
  include(joinpath(test_dir, f))
end


text:
# Lebesgue measure

export Lebesgue

struct Lebesgue{X} <: PrimitiveMeasure end

function Base.show(io::IO, μ::Lebesgue{X}) where X
    io = IOContext(io, :compact => true)
    print(io, "Lebesgue(")
    print(io, X)
    print(io, ")")
end

Lebesgue(X) = Lebesgue{X}()

sampletype(::Lebesgue{ℝ}) = Float64
sampletype(::Lebesgue{ℝ₊}) = Float64
sampletype(::Lebesgue{𝕀}) = Float64

testvalue(::Lebesgue{ℝ}) = 0.0
testvalue(::Lebesgue{𝕀}) = 0.5
testvalue(::Lebesgue{ℝ₊}) = 1.0
testvalue(::Lebesgue{<:Real}) = 0.0

logdensity(::Lebesgue, x) = zero(x)

Base.:∘(::typeof(basemeasure), ::Type{Lebesgue}) = Lebesgue


text:
@inline function _det(r11, r12, r13, r21, r22, r23, r31, r32, r33)
    r11 * r22 * r33 -
    r11 * r32 * r23 -
    r21 * r12 * r33 +
    r21 * r32 * r13 +
    r31 * r12 * r23 -
    r31 * r22 * r13
end

@inline function _mul_trace(
    x11, x12, x13, x21, x22, x23, x31, x32, x33,
    y11, y12, y13, y21, y22, y23, y31, y32, y33
)

    return (x11 * y11 + x12 * y21 + x13 * y31) +  # z11
           (x21 * y12 + x22 * y22 + x23 * y32) +  # z22
           (x31 * y13 + x32 * y23 + x33 * y33)    # z33
end

get_sform(x::NIVolume) = get_sform(x.header)
function get_sform(hdr::NIfTI1Header)
    if hdr.sform_code > 0
        return @inbounds Float32[
            hdr.srow_x[1]  hdr.srow_x[2]  hdr.srow_x[3]  hdr.srow_x[4]
            hdr.srow_y[1]  hdr.srow_y[2]  hdr.srow_y[3]  hdr.srow_y[4]
            hdr.srow_z[1]  hdr.srow_z[2]  hdr.srow_z[3]  hdr.srow_z[4]
            0              0              0              1
        ]
    else
        # TODO there's no sform in this case, should we just throw an error?
        return nothing
    end
end

get_qform(x::NIVolume) = get_qform(x.header)
function get_qform(hdr::NIfTI1Header)
    if hdr.qform_code <= 0
        return @inbounds Float32[
            hdr.pixdim[2]   0              0              0
            0               hdr.pixdim[3]  0              0
            0               0              hdr.pixdim[4]  0
            0               0              0              1
        ]
    else
        dx = hdr.pixdim[2]
        dy = hdr.pixdim[3]
        # aka qfac left handedness
        if hdr.pixdim[1] < 0
            dz = -hdr.pixdim[4]
        else
            dz = hdr.pixdim[4]
        end
        b = hdr.quatern_b
        c = hdr.quatern_c
        d = hdr.quatern_d
        b2 = b*b
        c2 = c*c
        d2 = d*d
        a = 1 - b2 - c2 - d2
        if a < 1.e-7
            a = 1 / sqrt(b2 + c2 + d2)
            b *= a
            c *= a
            d *= a       # normalize (b,c,d) vector
            a = zero(a)  # a = 0 ==> 180 degree rotation
        else
            a = sqrt(a)   # angle = 2*arccos(a)
        end
        return Float32[
            (a*a+b*b-c*c-d*d)*dx   (2*b*c-2*a*d)*dy       (2*b*d+2*a*c)*dz      hdr.qoffset_x
            (2*b*c+2*a*d)*dx       (a*a+c*c-b*b-d*d)*dy   (2*c*d-2*a*b)*dz      hdr.qoffset_y
            (2*b*d-2*a*c)*dx       (2*c*d+2*a*b)*dy       (a*a+d*d-c*c-b*b)*dz  hdr.qoffset_z
            0                      0                      0                     1
        ]
    end
end
"""
    getaffine(x::NIVolume)

Gets a 4x4 affine transformation volume's header's sform
"""
getaffine(x::NIVolume) = getaffine(x.header)

# Convert a NIfTI header to a 4x4 affine transformation matrix
"""
    getaffine(hdr::NIfTI1Header)

Gets a 4x4 affine transformation matrix from a header's sform
"""
function getaffine(hdr::NIfTI1Header)
    if hdr.sform_code > 0
        return get_sform(hdr)
    else
        return get_qform(hdr)
    end
end

"""
    function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}

Set the affine of a `NIfTI1Header` to 4x4 affine matrix `affine`
"""
function setaffine(h::NIfTI1Header, affine::Array{T,2}) where {T}
    size(affine, 1) == size(affine, 2) == 4 ||
        error("affine matrix must be 4x4")
    affine[4, 1] == affine[4, 2] == affine[4, 3] == 0 && affine[4, 4] == 1 ||
        error("last row of affine matrix must be [0 0 0 1]")
    h.qform_code = one(Int16)
    h.sform_code = one(Int16)
    h.pixdim = (zero(Float32), h.pixdim[2:end]...,)
    h.quatern_b = zero(Float32)
    h.quatern_c = zero(Float32)
    h.quatern_d = zero(Float32)
    h.qoffset_x = zero(Float32)
    h.qoffset_y = zero(Float32)
    h.qoffset_z = zero(Float32)
    h.srow_x = convert(Tuple{Vararg{Float32}}, (affine[1, :]...,))
    h.srow_y = convert(Tuple{Vararg{Float32}}, (affine[2, :]...,))
    h.srow_z = convert(Tuple{Vararg{Float32}}, (affine[3, :]...,))
    h
end

"""
    orientation(img)::Tuple{Symbol,Symbol,Symbol}

Returns a tuple providing the orientation of a NIfTI image.
"""
orientation(x) = orientation(x.header)
function orientation(hdr::NIfTI1Header)
    if hdr.sform_code > 0
        return @inbounds _dir2ori(
            hdr.srow_x[1], hdr.srow_x[2], hdr.srow_x[3],
            hdr.srow_y[1], hdr.srow_y[2], hdr.srow_y[3],
            hdr.srow_z[1], hdr.srow_z[2], hdr.srow_z[3]
        )
    elseif hdr.qform_code <= 0
        return @inbounds _dir2ori(
            hdr.pixdim[2], 0, 0,
            0, hdr.pixdim[3], 0,
            0, 0, hdr.pixdim[4]
        )
    else
        dx = hdr.pixdim[2]
        dy = hdr.pixdim[3]
        # aka qfac left handedness
        if hdr.pixdim[1] < 0
            dz = -hdr.pixdim[4]
        else
            dz = hdr.pixdim[4]
        end
        b = hdr.quatern_b
        c = hdr.quatern_c
        d = hdr.quatern_d
        b2 = b*b
        c2 = c*c
        d2 = d*d
        a = 1 - b2 - c2 - d2
        if a < 1.e-7
            a = 1 / sqrt(b2 + c2 + d2)
            b *= a
            c *= a
            d *= a       # normalize (b,c,d) vector
            a = zero(a)  # a = 0 ==> 180 degree rotation
        else
            a = sqrt(a)   # angle = 2*arccos(a)
        end
        return _dir2ori(
            (a*a+b*b-c*c-d*d)*dx,  (2*b*c-2*a*d)*dy,      (2*b*d+2*a*c)*dz,
            (2*b*c+2*a*d)*dx,      (a*a+c*c-b*b-d*d)*dy,  (2*c*d-2*a*b)*dz,
            (2*b*d-2*a*c)*dx,      (2*c*d+2*a*b)*dy,      (a*a+d*d-c*c-b*b)*dz,
        )
    end
end

_encoding_name(x) = _encoding_name(Int(x))
@inline function _encoding_name(x::Int)
    if x === 1
        return :left
    elseif x === -1
        return :right
    elseif x === 2
        return :posterior
    elseif x === -2
        return :anterior
    elseif x === 3
        return :inferior
    elseif x === -3
        return :superior
    else
        error("$x does not map to a dimension name.")
    end
end

function _dir2ori(xi, xj, xk, yi, yj, yk, zi, zj, zk)
    # Normalize column vectors to get unit vectors along each ijk-axis
    # normalize i axis
    val = sqrt(xi*xi + yi*yi + zi*zi)
    if val == 0
        error("Invalid rotation directions.")
    end
    xi /= val
    yi /= val
    zi /= val

    # normalize j axis
    val = sqrt(xj*xj + yj*yj + zj*zj)
    if val == 0
        error("Invalid rotation directions.")
    end
    xj /= val
    yj /= val
    zj /= val

    # orthogonalize j axis to i axis, if needed
    val = xi*xj + yi*yj + zi* zj  # dot product between i and j
    if abs(val) > .0001
        xj -= val*xi
        yj -= val*yi
        zj -= val*zi

        val = sqrt(xj*xj + yj*yj + zj*zj)  # must renormalize
        if val == 0
            error("The first and second dimensions cannot be parallel.")
        end
        xj /= val
        yj /= val
        zj /= val
    end

    # normalize k axis; if it is zero, make it the cross product i x j
    val = sqrt(xk*xk + yk*yk + zk*zk)
    if val == 0
        xk = yi*zj-zi*yj
        yk = zi*xj-zj*xi
        zk = xi*yj-yi*xj
    else
        xk = xk/val
        yk = yk/val
        zk = zk/val
    end

    # orthogonalize k to i
    val = xi*xk + yi*yk + zi*zk  # dot product between i and k
    if abs(val) > 0.0001
        xk -= val*xi
        yk -= val*yi
        zk -= val*zi

        # must renormalize
        val = sqrt(xk*xk + yk*yk + zk*zk)
        if val == 0
            return 0  # I think this is suppose to be an error output
        end
        xk /= val
        yk /= val
        zk /= val
    end

    # orthogonalize k to j */
    val = xj*xk + yj*yk + zj*zk  # dot product between j and k
    if abs(val) > 0.0001
        xk -= val*xj
        yk -= val*yj
        zk -= val*zj

        val = sqrt(xk*xk + yk*yk + zk*zk)
        if val == 0
            return 0  # bad
        end
        xk /= val
        yk /= val
        zk /= val
    end

    # at this point Q is the rotation matrix from the (i,j,k) to (x,y,z) axes
    detQ = _det(xi, xj, xk, yi, yj, yk, zi, zj, zk)
    # if( detQ == 0.0 ) return ; /* shouldn't happen unless user is a DUFIS */

    # Build and test all possible +1/-1 coordinate permutation matrices P;
    # then find the P such that the rotation matrix M=PQ is closest to the
    # identity, in the sense of M having the smallest total rotation angle.

    # Despite the formidable looking 6 nested loops, there are
    # only 3*3*3*2*2*2 = 216 passes, which will run very quickly.
    vbest = -666
    ibest = pbest=qbest=rbest= 1
    jbest = 2
    kbest = 3
    for (i, j, k) in ((1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1))
        for p in (-1, 1)           # p,q,r are -1 or +1
            for q in (-1, 1)       # and go into rows 1,2,3
                for r in (-1, 1)
                    p11, p12, p13 = _nval_other_zero(i, p)
                    p21, p22, p23 = _nval_other_zero(j, q)
                    p31, p32, p33 = _nval_other_zero(k, r)
                    #=
                    P[1,i] = p
                    P[2,j] = q
                    P[3,k] = r
                    detP = det(P)  # sign of permutation
                    =#
                    detP = _det(p11, p12, p13, p21, p22, p23, p31, p32, p33)
                    # doesn't match sign of Q
                    if detP * detQ >= 0.0
                        # angle of M rotation = 2.0 * acos(0.5 * sqrt(1.0 + trace(M)))
                        # we want largest trace(M) == smallest angle == M nearest to I
                        val = _mul_trace(
                            p11, p12, p13, p21, p22, p23, p31, p32, p33,
                            xi, xj, xk, yi, yj, yk, zi, zj, zk
                        )
                        if val > vbest
                            vbest = val
                            ibest = i
                            jbest = j
                            kbest = k
                            pbest = p
                            qbest = q
                            rbest = r
                        end
                    end
                end
            end
        end
    end
    # At this point ibest is 1 or 2 or 3; pbest is -1 or +1; etc.

    # The matrix P that corresponds is the best permutation approximation
    # to Q-inverse; that is, P (approximately) takes (x,y,z) coordinates
    # to the (i,j,k) axes.

    # For example, the first row of P (which contains pbest in column ibest)
    # determines the way the i axis points relative to the anatomical
    # (x,y,z) axes.  If ibest is 2, then the i axis is along the y axis,
    # which is direction P2A (if pbest > 0) or A2P (if pbest < 0).

    # So, using ibest and pbest, we can assign the output code for
    # the i axis.  Mutatis mutandis for the j and k axes, of course.

    return (_encoding_name(ibest*pbest), _encoding_name(jbest*qbest), _encoding_name(kbest*rbest))
end

@inline function _nval_other_zero(n, val)
    if n === 1
        return val, 0, 0
    elseif n === 2
        return 0, val, 0
    else
        return 0, 0, val
    end
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.DictEncoding

Represents the "pool" of possible values for a [`DictEncoded`](@ref)
array type. Whether the order of values is significant can be checked
by looking at the `isOrdered` boolean field.

The `S` type parameter, while not tied directly to any field, is the
signed integer "index type" of the parent DictEncoded. We keep track
of this in the DictEncoding in order to validate the length of the pool
doesn't exceed the index type limit. The general workflow of writing arrow
data means the initial schema will typically be based off the data in the
first record batch, and subsequent record batches need to match the same
schema exactly. For example, if a non-first record batch dict encoded column
were to cause a DictEncoding pool to overflow on unique values, a fatal error
should be thrown.
"""
mutable struct DictEncoding{T, S, A} <: ArrowVector{T}
    id::Int64
    data::A
    isOrdered::Bool
    metadata::Union{Nothing, Dict{String, String}}
end

indextype(::Type{DictEncoding{T, S, A}}) where {T, S, A} = S
indextype(::T) where {T <: DictEncoding} = indextype(T)

Base.size(d::DictEncoding) = size(d.data)

@propagate_inbounds function Base.getindex(d::DictEncoding{T}, i::Integer) where {T}
    @boundscheck checkbounds(d, i)
    return @inbounds ArrowTypes.arrowconvert(T, d.data[i])
end

# convenience wrapper to signal that an input column should be
# dict encoded when written to the arrow format
struct DictEncodeType{T} end
getT(::Type{DictEncodeType{T}}) where {T} = T

"""
    Arrow.DictEncode(::AbstractVector, id::Integer=nothing)

Signals that a column/array should be dictionary encoded when serialized
to the arrow streaming/file format. An optional `id` number may be provided
to signal that multiple columns should use the same pool when being
dictionary encoded.
"""
struct DictEncode{T, A} <: AbstractVector{DictEncodeType{T}}
    id::Int64
    data::A
end

DictEncode(x::A, id=-1) where {A} = DictEncode{eltype(A), A}(id, x)
Base.IndexStyle(::Type{<:DictEncode}) = Base.IndexLinear()
Base.size(x::DictEncode) = (length(x.data),)
Base.iterate(x::DictEncode, st...) = iterate(x.data, st...)
Base.getindex(x::DictEncode, i::Int) = getindex(x.data, i)
ArrowTypes.ArrowKind(::Type{<:DictEncodeType}) = DictEncodedKind()
Base.copy(x::DictEncode) = DictEncode(x.data, x.id)

"""
    Arrow.DictEncoded

A dictionary encoded array type (similar to a `PooledArray`). Behaves just
like a normal array in most respects; internally, possible values are stored
in the `encoding::DictEncoding` field, while the `indices::Vector{<:Integer}`
field holds the "codes" of each element for indexing into the encoding pool.
Any column/array can be dict encoding when serializing to the arrow format
either by passing the `dictencode=true` keyword argument to [`Arrow.write`](@ref)
(which causes _all_ columns to be dict encoded), or wrapping individual columns/
arrays in [`Arrow.DictEncode(x)`](@ref).
"""
struct DictEncoded{T, S, A} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    indices::Vector{S}
    encoding::DictEncoding{T, S, A}
    metadata::Union{Nothing, Dict{String, String}}
end

DictEncoded(b::Vector{UInt8}, v::ValidityBitmap, inds::Vector{S}, encoding::DictEncoding{T, S, A}, meta) where {S, T, A} =
    DictEncoded{T, S, A}(b, v, inds, encoding, meta)

Base.size(d::DictEncoded) = size(d.indices)

isdictencoded(d::DictEncoded) = true
isdictencoded(x) = false
isdictencoded(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = true

function signedtype(n::Integer)
    typs = (Int8, Int16, Int32, Int64)
    typs[something(findfirst(n .≤ typemax.(typs)), 4)]
end

signedtype(::Type{UInt8}) = Int8
signedtype(::Type{UInt16}) = Int16
signedtype(::Type{UInt32}) = Int32
signedtype(::Type{UInt64}) = Int64
signedtype(::Type{T}) where {T <: Signed} = T

indtype(d::DictEncoded{T, S, A}) where {T, S, A} = S
indtype(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = indtype(c.data)

dictencodeid(colidx, nestedlevel, fieldid) = (Int64(nestedlevel) << 48) | (Int64(fieldid) << 32) | Int64(colidx)

getid(d::DictEncoded) = d.encoding.id
getid(c::Compressed{Z, A}) where {Z, A <: DictEncoded} = c.data.encoding.id

function arrowvector(::DictEncodedKind, x::DictEncoded, i, nl, fi, de, ded, meta; dictencode::Bool=false, dictencodenested::Bool=false, kw...)
    id = x.encoding.id
    if !haskey(de, id)
        de[id] = Lockable(x.encoding)
    else
        encodinglockable = de[id]
        @lock encodinglockable begin
            encoding = encodinglockable.x
            # in this case, we just need to check if any values in our local pool need to be delta dicationary serialized
            deltas = setdiff(x.encoding, encoding)
            if !isempty(deltas)
                ET = indextype(encoding)
                if length(deltas) + length(encoding) > typemax(ET)
                    error("fatal error serializing dict encoded column with ref index type of $ET; subsequent record batch unique values resulted in $(length(deltas) + length(encoding)) unique values, which exceeds possible index values in $ET")
                end
                data = arrowvector(deltas, i, nl, fi, de, ded, nothing; dictencode=dictencodenested, dictencodenested=dictencodenested, dictencoding=true, kw...)
                push!(ded, DictEncoding{eltype(data), ET, typeof(data)}(id, data, false, getmetadata(data)))
                if typeof(encoding.data) <: ChainedVector
                    append!(encoding.data, data)
                else
                    data2 = ChainedVector([encoding.data, data])
                    encoding = DictEncoding{eltype(data2), ET, typeof(data2)}(id, data2, false, getmetadata(encoding))
                    de[id].x = encoding
                end
            end
        end
    end
    return x
end

function arrowvector(::DictEncodedKind, x, i, nl, fi, de, ded, meta; dictencode::Bool=false, dictencodenested::Bool=false, kw...)
    @assert x isa DictEncode
    id = x.id == -1 ? dictencodeid(i, nl, fi) : x.id
    x = x.data
    len = length(x)
    validity = ValidityBitmap(x)
    if !haskey(de, id)
        # dict encoding doesn't exist yet, so create for 1st time
        if DataAPI.refarray(x) === x || DataAPI.refpool(x) === nothing
            # need to encode ourselves
            x = PooledArray(x; signed=true, compress=true)
            inds = DataAPI.refarray(x)
            pool = DataAPI.refpool(x)
        else
            pool = DataAPI.refpool(x)
            refa = DataAPI.refarray(x)
            inds = copyto!(similar(Vector{signedtype(length(pool))}, length(refa)), refa)
        end
        # horrible hack? yes. better than taking CategoricalArrays dependency? also yes.
        if typeof(pool).name.name == :CategoricalRefPool
            if eltype(x) >: Missing
                pool = vcat(missing, DataAPI.levels(x))
            else
                pool = DataAPI.levels(x)
                for i = 1:length(inds)
                    @inbounds inds[i] -= 1
                end
            end
        else
            # adjust to "offset" instead of index
            for i = 1:length(inds)
                @inbounds inds[i] -= 1
            end
        end
        data = arrowvector(pool, i, nl, fi, de, ded, nothing; dictencode=dictencodenested, dictencodenested=dictencodenested, dictencoding=true, kw...)
        encoding = DictEncoding{eltype(data), eltype(inds), typeof(data)}(id, data, false, getmetadata(data))
        de[id] = Lockable(encoding)
    else
        # encoding already exists
          # compute inds based on it
          # if value doesn't exist in encoding, push! it
          # also add to deltas updates
        encodinglockable = de[id]
        @lock encodinglockable begin
            encoding = encodinglockable.x
            len = length(x)
            ET = indextype(encoding)
            pool = Dict{Union{eltype(encoding), eltype(x)}, ET}(a => (b - 1) for (b, a) in enumerate(encoding))
            deltas = eltype(x)[]
            inds = Vector{ET}(undef, len)
            categorical = typeof(x).name.name == :CategoricalArray
            for (j, val) in enumerate(x)
                if categorical
                    val = get(val)
                end
                @inbounds inds[j] = get!(pool, val) do
                    push!(deltas, val)
                    return length(pool)
                end
            end
            if !isempty(deltas)
                if length(deltas) + length(encoding) > typemax(ET)
                    error("fatal error serializing dict encoded column with ref index type of $ET; subsequent record batch unique values resulted in $(length(deltas) + length(encoding)) unique values, which exceeds possible index values in $ET")
                end
                data = arrowvector(deltas, i, nl, fi, de, ded, nothing; dictencode=dictencodenested, dictencodenested=dictencodenested, dictencoding=true, kw...)
                push!(ded, DictEncoding{eltype(data), ET, typeof(data)}(id, data, false, getmetadata(data)))
                if typeof(encoding.data) <: ChainedVector
                    append!(encoding.data, data)
                else
                    data2 = ChainedVector([encoding.data, data])
                    encoding = DictEncoding{eltype(data2), ET, typeof(data2)}(id, data2, false, getmetadata(encoding))
                    de[id].x = encoding
                end
            end
        end
    end
    if meta !== nothing && getmetadata(encoding) !== nothing
        merge!(meta, getmetadata(encoding))
    elseif getmetadata(encoding) !== nothing
        meta = getmetadata(encoding)
    end
    return DictEncoded(UInt8[], validity, inds, encoding, meta)
end

@propagate_inbounds function Base.getindex(d::DictEncoded, i::Integer)
    @boundscheck checkbounds(d, i)
    @inbounds valid = d.validity[i]
    !valid && return missing
    @inbounds idx = d.indices[i]
    return @inbounds d.encoding[idx + 1]
end

@propagate_inbounds function Base.setindex!(d::DictEncoded{T}, v, i::Integer) where {T}
    @boundscheck checkbounds(d, i)
    if v === missing
        @inbounds d.validity[i] = false
    else
        ix = findfirst(d.encoding.data, v)
        if ix === nothing
            push!(d.encoding.data, v)
            @inbounds d.indices[i] = length(d.encoding.data) - 1
        else
            @inbounds d.indices[i] = ix - 1
        end
    end
    return v
end

function Base.copy(x::DictEncoded{T, S}) where {T, S}
    pool = copy(x.encoding.data)
    valid = x.validity
    inds = x.indices
    refs = copy(inds)
    @inbounds for i = 1:length(inds)
        refs[i] = refs[i] + one(S)
    end
    return PooledArray(PooledArrays.RefArray(refs), Dict{T, S}(val => i for (i, val) in enumerate(pool)), pool)
end

function compress(Z::Meta.CompressionType, comp, x::A) where {A <: DictEncoded}
    len = length(x)
    nc = nullcount(x)
    validity = compress(Z, comp, x.validity)
    inds = compress(Z, comp, x.indices)
    return Compressed{Z, A}(x, [validity, inds], len, nc, Compressed[])
end

function DataAPI.levels(x::DictEncoded)
    rp = DataAPI.refpool(x)   # may contain missing values
    Missing <: eltype(rp) || return rp
    convert(AbstractArray{nonmissingtype(eltype(rp))}, deleteat!(rp, ismissing.(rp)))
end

function makenodesbuffers!(col::DictEncoded{T, S}, fieldnodes, fieldbuffers, bufferoffset, alignment) where {T, S}
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += blen
    # indices
    blen = sizeof(S) * len
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += padding(blen, alignment)
    return bufferoffset
end

DataAPI.refarray(x::DictEncoded{T, S}) where {T, S} = x.indices .+ one(S)

DataAPI.refpool(x::DictEncoded) = copy(x.encoding.data)

function writebuffer(io, col::DictEncoded, alignment)
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    # write indices
    n = writearray(io, col.indices)
    @debug 1 "writing array: col = $(typeof(col.indices)), n = $n, padded = $(padding(n, alignment))"
    writezeros(io, paddinglength(n, alignment))
    return
end


text:
#include("../ReadWrite/Header.jl")

"""
**SeisBin**

*Bin seismic data onto regular grid using trace headers*

**IN**

* in: filename of input, irregularly sampled data
* out: filename of output, regularly sampled data
* style="sxsygxgy"
* ang=90
* gamma=1
* osx=0
* osy=0
* ogx=0
* ogy=0
* omx=0
* omy=0
* ohx=0
* ohy=0
* oh=0
* oaz=0
* dsx=1
* dsy=1
* dgx=1
* dgy=1
* dmx=1
* dmy=1
* dhx=1
* dhy=1
* dh=1
* daz=1
* min_isx=0
* max_isx=0
* min_isy=0
* max_isy=0
* min_igx=0
* max_igx=0
* min_igy=0
* max_igy=0
* min_imx=0
* max_imx=0
* min_imy=0
* max_imy=0
* min_ihx=0
* max_ihx=0
* min_ihy=0
* max_ihy=0
* min_ih=0
* max_ih=0
* min_iaz=0
* max_iaz=0

**OUT**

*Credits: Aaron Stanton, 2015*

"""

function SeisBin(in,out;style="sxsygxgy",ang=90,gamma=1,osx=0,osy=0,ogx=0,ogy=0,omx=0,omy=0,ohx=0,ohy=0,oh=0,oaz=0,dsx=1,dsy=1,dgx=1,dgy=1,dmx=1,dmy=1,dhx=1,dhy=1,dh=1,daz=1,min_isx=0,max_isx=0,min_isy=0,max_isy=0,min_igx=0,max_igx=0,min_igy=0,max_igy=0,min_imx=0,max_imx=0,min_imy=0,max_imy=0,min_ihx=0,max_ihx=0,min_ihy=0,max_ihy=0,min_ih=0,max_ih=0,min_iaz=0,max_iaz=0,ntrace=10000)

	rad2deg = 180/pi;
	deg2rad = pi/180;
	gammainv = 1/gamma;
	if (ang > 90)
		ang2=-deg2rad*(ang-90)
	else
		ang2=deg2rad*(90-ang)
	end

	nsx = max_isx - min_isx + 1
	nsy = max_isy - min_isy + 1
	ngx = max_igx - min_igx + 1
	ngy = max_igy - min_igy + 1
	nmx = max_imx - min_imx + 1
	nmy = max_imy - min_imy + 1
	nhx = max_ihx - min_ihx + 1
	nhy = max_ihy - min_ihy + 1
	nh  = max_ih  - min_ih  + 1
	naz = max_iaz - min_iaz + 1

	if (style=="sxsygxgy")
		nx1=nsx;nx2=nsy;nx3=ngx;nx4=ngy;
		ox1=osx;ox2=osy;ox3=ogx;ox4=ogy;
		dx1=dsx;dx2=dsy;dx3=dgx;dx4=dgy;
		label2="sx";label3="sy";label4="gx";label5="gy";
		unit2="m";unit3="m";unit4="m";unit5="m";
	elseif (style=="mxmyhxhy")
		nx1=nmx;nx2=nmy;nx3=nhx;nx4=nhy;
		ox1=omx;ox2=omy;ox3=ohx;ox4=ohy;
		dx1=dmx;dx2=dmy;dx3=dhx;dx4=dhy;
		label2="mx";label3="my";label4="hx";label5="hy";
		unit2="m";unit3="m";unit4="m";unit5="m";
	elseif (style=="mxmyhaz")
		nx1=nmx;nx2=nmy;nx3=nh;nx4=naz;
		ox1=omx;ox2=omy;ox3=oh;ox4=oaz;
		dx1=dmx;dx2=dmy;dx3=dh;dx4=daz;
		label2="mx";label3="my";label4="h";label5="az";
		unit2="m";unit3="m";unit4="m";unit5="Degrees";
	elseif (style=="sxsyhxhy")
		nx1=nsx;nx2=nsy;nx3=nhx;nx4=nhy;
		ox1=osx;ox2=osy;ox3=ohx;ox4=ohy;
		dx1=dsx;dx2=dsy;dx3=dhx;dx4=dhy;
		label2="sx";label3="sy";label4="hx";label5="hy";
		unit2="m";unit3="m";unit4="m";unit5="m";
	elseif (style=="gxgyhxhy")
		nx1=ngx;nx2=ngy;nx3=nhx;nx4=nhy;
		ox1=ogx;ox2=ogy;ox3=ohx;ox4=ohy;
		dx1=dgx;dx2=dgy;dx3=dhx;dx4=dhy;
		label2="gx";label3="gy";label4="hx";label5="hy";
		unit2="m";unit3="m";unit4="m";unit5="m";
	elseif (style=="sxsyhaz")
		nx1=nsx;nx2=nsy;nx3=nh;nx4=naz;
		ox1=osx;ox2=osy;ox3=oh;ox4=oaz;
		dx1=dsx;dx2=dsy;dx3=dh;dx4=daz;
		label2="sx";label3="sy";label4="hx";label5="az";
		unit2="m";unit3="m";unit4="m";unit5="Degrees";
	elseif (style=="gxgyhaz")
		nx1=ngx;nx2=ngy;nx3=nh;nx4=naz;
		ox1=ogx;ox2=ogy;ox3=oh;ox4=oaz;
		dx1=dgx;dx2=dgy;dx3=dh;dx4=daz;
		label2="gx";label3="gy";label4="h";label5="az";
		unit2="m";unit3="m";unit4="m";unit5="Degrees";
	else
		error("style not recognized.")
	end
	nx_out = nx1*nx2*nx3*nx4

	stream = open(ParseHeaderName(in))
	@compat nx_in = int(filesize(stream)/(4*length(fieldnames(Header))))
	seek(stream, header_count["n1"])
	nt = read(stream,Int32)
	seek(stream, header_count["d1"])
	dt = read(stream,Float32)
	close(stream)
	d = zeros(Float32,nt,1)
	h = Array(Header,1)
	h[1] = InitSeisHeader()

	extent = Extent(convert(Int32,nt),convert(Int32,nx1),convert(Int32,nx2),convert(Int32,nx3),convert(Int32,nx4),
		   convert(Float32,0),convert(Float32,ox1),convert(Float32,ox2),convert(Float32,ox3),convert(Float32,ox4),
		   convert(Float32,dt),convert(Float32,dx1),convert(Float32,dx2),convert(Float32,dx3),convert(Float32,dx4),
		   "Time",label2,label3,label4,label5,
		   "s",unit2,unit3,unit4,unit5,
		   "")

	j = 1
	for ix1 = 1 : nx1
		for ix2 = 1 : nx2
			for ix3 = 1 : nx3
				for ix4 = 1 : nx4
					h[1].tracenum = convert(typeof(h[1].tracenum),j)
					h[1].o1 = convert(typeof(h[1].o1),0)
					h[1].n1 = convert(typeof(h[1].n1),nt)
					h[1].d1 = convert(typeof(h[1].d1),dt)
					if (style=="sxsygxgy")
						h[1].isx = convert(typeof(h[1].isx),ix1 - 1 + min_isx)
						h[1].isy = convert(typeof(h[1].isy),ix2 - 1 + min_isy)
						h[1].igx = convert(typeof(h[1].igx),ix3 - 1 + min_igx)
						h[1].igy = convert(typeof(h[1].igy),ix4 - 1 + min_igy)
						sx_rot = convert(Float32,(ix1 - 1 + min_isx)*dsx + osx);
						sy_rot = convert(Float32,(ix2 - 1 + min_isy)*dsy + osy);
						gx_rot = convert(Float32,(ix3 - 1 + min_igx)*dgx + ogx);
						gy_rot = convert(Float32,(ix4 - 1 + min_igy)*dgy + ogy);
						h[1].sx =  (sx_rot-osx)*cos(ang2) + (sy_rot-osy)*sin(ang2) + osx;
						h[1].sy = -(sx_rot-osx)*sin(ang2) + (sy_rot-osy)*cos(ang2) + osy;
						h[1].gx =  (gx_rot-ogx)*cos(ang2) + (gy_rot-ogy)*sin(ang2) + ogx;
						h[1].gy = -(gx_rot-ogx)*sin(ang2) + (gy_rot-ogy)*cos(ang2) + ogy;
						h[1].hx = h[1].gx - h[1].sx
						h[1].hy = h[1].gy - h[1].sy
						h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
						h[1].az = rad2deg*atan2((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
						if (h[1].az < 0)
							h[1].az += 360.0
						end
						h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
						h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
						mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
						my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
						hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
						hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
						h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
						h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
						h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
						h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
						h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
						h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))
					elseif (style=="mxmyhxhy")
						h[1].imx = convert(typeof(h[1].imx),ix1 - 1 + min_imx)
						h[1].imy = convert(typeof(h[1].imy),ix2 - 1 + min_imy)
						h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
						h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
						mx_rot = convert(Float32,(ix1 - 1 + min_imx)*dmx + omx);
						my_rot = convert(Float32,(ix2 - 1 + min_imy)*dmy + omy);
						hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
						hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
						h[1].mx =  (mx_rot-omx)*cos(ang2) + (my_rot-omy)*sin(ang2) + omx;
						h[1].my = -(mx_rot-omx)*sin(ang2) + (my_rot-omy)*cos(ang2) + omy;
						h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
						h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
						h[1].sx = h[1].mx - h[1].hx/(1 + gammainv);
						h[1].sy = h[1].my - h[1].hy/(1 + gammainv);
						h[1].gx = h[1].mx + h[1].hx*(1-(1/(1 + gammainv)));
						h[1].gy = h[1].my + h[1].hy*(1-(1/(1 + gammainv)));
						sx_rot = (h[1].sx-osx)*cos(ang2) - (h[1].sy-osy)*sin(ang2) + osx;
						sy_rot = (h[1].sx-osx)*sin(ang2) + (h[1].sy-osy)*cos(ang2) + osy;
						gx_rot = (h[1].gx-ogx)*cos(ang2) - (h[1].gy-ogy)*sin(ang2) + ogx;
						gy_rot = (h[1].gx-ogx)*sin(ang2) + (h[1].gy-ogy)*cos(ang2) + ogy;
						h[1].isx = convert(Int32,round((sx_rot-osx)/dsx))
						h[1].isy = convert(Int32,round((sy_rot-osy)/dsy))
						h[1].igx = convert(Int32,round((gx_rot-ogx)/dgx))
						h[1].igy = convert(Int32,round((gy_rot-ogy)/dgy))
						h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
						h[1].az = rad2deg*atan2((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
						if (h[1].az < 0)
							h[1].az += 360.0
						end
						h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
						h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))
					elseif (style=="mxmyhaz")
						h[1].imx = convert(typeof(h[1].imx),ix1 - 1 + min_imx)
						h[1].imy = convert(typeof(h[1].imy),ix2 - 1 + min_imy)
						h[1].ih = convert(typeof(h[1].ih),  ix3 - 1 + min_ih)
						h[1].iaz = convert(typeof(h[1].iaz),ix4 - 1 + min_iaz)
						mx_rot = convert(Float32,(ix1 - 1 + min_imx)*dmx + omx);
						my_rot = convert(Float32,(ix2 - 1 + min_imy)*dmy + omy);
						h[1].mx =  (mx_rot-omx)*cos(ang2) + (my_rot-omy)*sin(ang2) + omx;
						h[1].my = -(mx_rot-omx)*sin(ang2) + (my_rot-omy)*cos(ang2) + omy;
						h[1].h = convert(Float32,(ix3 - 1 + min_ih)*dh + oh);
						h[1].az = convert(Float32,(ix4 - 1 + min_iaz)*daz + oaz);
						if (h[1].az <= 90)
							h[1].hx = h[1].h*cos(deg2rad*h[1].az);
							h[1].hy = h[1].h*sin(deg2rad*h[1].az);
						elseif (h[1].az > 90 && h[1].az <= 180)
							h[1].hx =-h[1].h*cos(pi-(deg2rad*h[1].az));
							h[1].hy = h[1].h*sin(pi-(deg2rad*h[1].az));
						elseif (h[1].az > 180 && h[1].az <= 270)
							h[1].hx =-h[1].h*cos((deg2rad*h[1].az)-pi);
							h[1].hy =-h[1].h*sin((deg2rad*h[1].az)-pi);
						else
							h[1].hx = h[1].h*cos(2*pi-(deg2rad*h[1].az));
							h[1].hy =-h[1].h*sin(2*pi-(deg2rad*h[1].az));
						end
						h[1].sx = h[1].mx - h[1].hx/(1 + gammainv);
						h[1].sy = h[1].my - h[1].hy/(1 + gammainv);
						h[1].gx = h[1].mx + h[1].hx*(1-(1/(1 + gammainv)));
						h[1].gy = h[1].my + h[1].hy*(1-(1/(1 + gammainv)));
						sx_rot = (h[1].sx-osx)*cos(ang2) - (h[1].sy-osy)*sin(ang2) + osx;
						sy_rot = (h[1].sx-osx)*sin(ang2) + (h[1].sy-osy)*cos(ang2) + osy;
						gx_rot = (h[1].gx-ogx)*cos(ang2) - (h[1].gy-ogy)*sin(ang2) + ogx;
						gy_rot = (h[1].gx-ogx)*sin(ang2) + (h[1].gy-ogy)*cos(ang2) + ogy;
						h[1].isx = convert(Int32,round((sx_rot-osx)/dsx))
						h[1].isy = convert(Int32,round((sy_rot-osy)/dsy))
						h[1].igx = convert(Int32,round((gx_rot-ogx)/dgx))
						h[1].igy = convert(Int32,round((gy_rot-ogy)/dgy))
						hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
						hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
						h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
						h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
					elseif (style=="sxsyhxhy")
						h[1].isx = convert(typeof(h[1].isx),ix1 - 1 + min_isx)
						h[1].isy = convert(typeof(h[1].isy),ix2 - 1 + min_isy)
						h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
						h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
						sx_rot = convert(Float32,(ix1 - 1 + min_isx)*dsx + osx);
						sy_rot = convert(Float32,(ix2 - 1 + min_isy)*dsy + osy);
						hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
						hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
						h[1].sx =  (sx_rot-osx)*cos(ang2) + (sy_rot-osy)*sin(ang2) + osx;
						h[1].sy = -(sx_rot-osx)*sin(ang2) + (sy_rot-osy)*cos(ang2) + osy;
						h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
						h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
						h[1].gx = h[1].sx + h[1].hx;
						h[1].gy = h[1].sy + h[1].hy;
						h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
						h[1].az = rad2deg*atan2((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
						if (h[1].az < 0)
							h[1].az += 360.0
						end
						h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
						h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
						mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
						my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
						hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
						hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
						h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
						h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
						h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
						h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
						h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
						h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))
					elseif (style=="gxgyhxhy")
						h[1].igx = convert(typeof(h[1].igx),ix1 - 1 + min_igx)
						h[1].igy = convert(typeof(h[1].igy),ix2 - 1 + min_igy)
						h[1].ihx = convert(typeof(h[1].ihx),ix3 - 1 + min_ihx)
						h[1].ihy = convert(typeof(h[1].ihy),ix4 - 1 + min_ihy)
						gx_rot = convert(Float32,(ix1 - 1 + min_igx)*dgx + ogx);
						gy_rot = convert(Float32,(ix2 - 1 + min_igy)*dgy + ogy);
						hx_rot = convert(Float32,(ix3 - 1 + min_ihx)*dhx + ohx);
						hy_rot = convert(Float32,(ix4 - 1 + min_ihy)*dhy + ohy);
						h[1].gx =  (gx_rot-ogx)*cos(ang2) + (gy_rot-ogy)*sin(ang2) + ogx;
						h[1].gy = -(gx_rot-ogx)*sin(ang2) + (gy_rot-ogy)*cos(ang2) + ogy;
						h[1].hx =  (hx_rot-ohx)*cos(ang2) + (hy_rot-ohy)*sin(ang2) + ohx;
						h[1].hy = -(hx_rot-ohx)*sin(ang2) + (hy_rot-ohy)*cos(ang2) + ohy;
						h[1].sx = h[1].gx - h[1].hx;
						h[1].sy = h[1].gy - h[1].hy;
						h[1].h = sqrt((h[1].hx^2) + (h[1].hy^2))
						h[1].az = rad2deg*atan2((h[1].gy-h[1].sy),(h[1].gx-h[1].sx))
						if (h[1].az < 0)
							h[1].az += 360.0
						end
						h[1].mx = h[1].sx + h[1].hx/(1 + gammainv);
						h[1].my = h[1].sy + h[1].hy/(1 + gammainv);
						mx_rot = (h[1].mx-omx)*cos(ang2) - (h[1].my-omy)*sin(ang2) + omx;
						my_rot = (h[1].mx-omx)*sin(ang2) + (h[1].my-omy)*cos(ang2) + omy;
						hx_rot = (h[1].hx-ohx)*cos(ang2) - (h[1].hy-ohy)*sin(ang2) + ohx;
						hy_rot = (h[1].hx-ohx)*sin(ang2) + (h[1].hy-ohy)*cos(ang2) + ohy;
						h[1].imx = convert(Int32,round((mx_rot-omx)/dmx))
						h[1].imy = convert(Int32,round((my_rot-omy)/dmy))
						h[1].ihx = convert(Int32,round((hx_rot-ohx)/dhx))
						h[1].ihy = convert(Int32,round((hy_rot-ohy)/dhy))
						h[1].ih = convert(Int32,round((h[1].h-oh)/dh))
						h[1].iaz = convert(Int32,round((h[1].az-oaz)/daz))
					elseif (style=="sxsyhaz")
						error("sxsyhaz not developed yet.")
					elseif (style=="gxgyhaz")
						error("gxgyhaz not developed yet.")
					end
					h[1].selev = convert(typeof(h[1].selev),0)
					h[1].gelev = convert(typeof(h[1].gelev),0)
					h[1].trid = convert(typeof(h[1].trid),0)
					SeisWrite(out,d,h,extent,itrace=j)
					j += 1
				end
			end
		end
	end

	out_d = ParseDataName(out)
	out_h = ParseHeaderName(out)
	stream_d = open(out_d,"a")
	stream_h = open(out_h,"a")

	if (style=="sxsygxgy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].isx - min_isx)*nx2*nx3*nx4 + (h[k].isy - min_isy)*nx3*nx4 + (h[k].igx - min_igx)*nx4 + h[k].igy - min_igy + 1
				if (h[k].isx >= min_isx && h[k].isx <= max_isx && h[k].isy >= min_isy && h[k].isy <= max_isy && h[k].igx >= min_igx && h[k].igx <= max_igx && h[k].igy >= min_igy && h[k].igy <= max_igy)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="mxmyhxhy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].imx - min_imx)*nx2*nx3*nx4 + (h[k].imy - min_imy)*nx3*nx4 + (h[k].ihx - min_ihx)*nx4 + h[k].ihy - min_ihy + 1
				if (h[k].imx >= min_imx && h[k].imx <= max_imx && h[k].imy >= min_imy && h[k].imy <= max_imy && h[k].ihx >= min_ihx && h[k].ihx <= max_ihx && h[k].ihy >= min_ihy && h[k].ihy <= max_ihy)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="mxmyhaz")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].imx - min_imx)*nx2*nx3*nx4 + (h[k].imy - min_imy)*nx3*nx4 + (h[k].ih - min_ih)*nx4 + h[k].iaz - min_iaz + 1
				if (h[k].imx >= min_imx && h[k].imx <= max_imx && h[k].imy >= min_imy && h[k].imy <= max_imy && h[k].ih >= min_ih && h[k].ih <= max_ih && h[k].iaz >= min_iaz && h[k].iaz <= max_iaz)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="sxsyhxhy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].isx - min_isx)*nx2*nx3*nx4 + (h[k].isy - min_isy)*nx3*nx4 + (h[k].ihx - min_ihx)*nx4 + h[k].ihy - min_ihy + 1
				if (h[k].isx >= min_isx && h[k].isx <= max_isx && h[k].isy >= min_isy && h[k].isy <= max_isy && h[k].ihx >= min_ihx && h[k].ihx <= max_ihx && h[k].ihy >= min_ihy && h[k].ihy <= max_ihy)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="gxgyhxhy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].igx - min_igx)*nx2*nx3*nx4 + (h[k].igy - min_igy)*nx3*nx4 + (h[k].ihx - min_ihx)*nx4 + h[k].ihy - min_ihy + 1
				if (h[k].igx >= min_igx && h[k].igx <= max_igx && h[k].igy >= min_igy && h[k].igy <= max_igy && h[k].ihx >= min_ihx && h[k].ihx <= max_ihx && h[k].ihy >= min_ihy && h[k].ihy <= max_ihy)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="sxsyhaz")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].isx - min_isx)*nx2*nx3*nx4 + (h[k].isy - min_isy)*nx3*nx4 + (h[k].ih - min_ih)*nx4 + h[k].iaz - min_iaz + 1
				if (h[k].isx >= min_isx && h[k].isx <= max_isx && h[k].isy >= min_isy && h[k].isy <= max_isy && h[k].ih >= min_ih && h[k].ih <= max_ih && h[k].iaz >= min_iaz && h[k].iaz <= max_iaz)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	elseif (style=="gxgyhaz")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h[k].igx - min_igx)*nx2*nx3*nx4 + (h[k].igy - min_igy)*nx3*nx4 + (h[k].ih - min_ih)*nx4 + h[k].iaz - min_iaz + 1
				if (h[k].igx >= min_igx && h[k].igx <= max_igx && h[k].igy >= min_igy && h[k].igy <= max_igy && h[k].ih >= min_ih && h[k].ih <= max_ih && h[k].iaz >= min_iaz && h[k].iaz <= max_iaz)
					h[k].tracenum = convert(Int32,itrace)
					if (itrace > 0 && itrace <= nx_out)
						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
						PutHeader(stream_h,h[k],itrace)
					end
				end
			end
			j += num_traces_in
		end
	end
	close(stream_d)
	close(stream_h)

end


text:
# Testing part 1

# Initialize Globals
const USE_GPU = false
const BENCHMARK = false
const VISUALIZE = false

# Run the XPU script
include("../scripts-part1/diffusion3D_xpu.jl")
xc, H = diffusion_3D(32)

# Run the Multi XPU script
ParallelStencil.@reset_parallel_stencil()
include("../scripts-part1/diffusion3D_multixpu.jl")
xc2, H2 = diffusion_3D(32)

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2; atol = 1e-5) for (v1,v2) in zip(values(d1), values(d2))])
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d_xpu = Dict(:X=> xc[inds], :H=>H[inds, inds, 15])
d_mxpu = Dict(:X=> xc2[inds], :H=>H2[inds, inds, 15])

@testset "Reference Test - Diffusion 3D" begin
    @test_reference "reftest-files/test_1.bson" d_xpu by=comp
    @test_reference "reftest-files/test_1.bson" d_mxpu by=comp
end


text:
@testset "Permute Tuner" begin
    funcvec=[random_decorrelate,fit_mix]
    instruct=Permute_Instruct(funcvec, [.5,.5],100,100)
    tuner=Permute_Tuner(instruct)
    #need to test update_weights functionality
    #want to supply some fake data to induce a .8 .2 categorical
    tuner.successes[:,1]=falses(TUNING_MEMORY*instruct.func_limit)
    tuner.successes[:,2]=falses(TUNING_MEMORY*instruct.func_limit)
    tuner.successes[1:Int(floor(TUNING_MEMORY*instruct.func_limit*.8)),1].=true
    tuner.successes[1:Int(floor(TUNING_MEMORY*instruct.func_limit*.2)),2].=true
    update_weights!(tuner)
    @test tuner.inst.weights==[.8,.2]    #check clamping
    tuner.successes[:,1]=falses(TUNING_MEMORY*instruct.func_limit)
    @test tuner.inst.min_clmps==[.01,.01]
    @test tuner.inst.max_clmps==[1.,1.]
    update_weights!(tuner)
    @test tuner.inst.weights==[.01,1-.01]

    #test override
    instruct=Permute_Instruct(funcvec, [.5,.5],100,100,override_time=2.,override_weights=[.75,.25])
    tuner=Permute_Tuner(instruct)
    tune_weights!(tuner, [(1,1.,1.)])
    @test tuner.inst.weights==[.5,.5] #no override
    tuner.time_history=fill(3.,CONVERGENCE_MEMORY)
    tune_weights!(tuner, [(1,1.,1.)])
    @test tuner.inst.weights==[.75,.25]

    #more clamping tests
    testvec=[0.02693244970132842, 0.031512823560878485, 0.050111382705776294, 0.6893314065796903, 0.04206537140157707, 0.02013161919125878, 0.026535815205008566, 0.051758654139053166, 0.03497967993105292, 0.013999262917669668, 0.01264153466670629]
    target=[0.02527900179828276, 0.029859375657832827, 0.04845793480273063, 0.6876779586766446, 0.040411923498531406, 0.02, 0.02488236730196291, 0.050105206236007505, 0.03332623202800726, 0.02, 0.02]

    clamp_pvec!(testvec,fill(.02,length(testvec)),fill(1.,length(testvec)))
    @test isprobvec(testvec)
    @test testvec==target

    #high clamp
    testvec=[.9,.08,.02]
    clamp_pvec!(testvec, fill(.05,length(testvec)), [.45,.45,.45])
    @test testvec==[.45, .2825, .2675]
end


text:
lp_solvers = Any[]
ip_solvers = Any[]
ip_dual_solvers = Any[]
semi_solvers = Any[]
sos_solvers = Any[]
lazy_solvers, lazy_soc_solvers, lazylocal_solvers, cut_solvers, cutlocal_solvers, heur_solvers, info_solvers = Any[], Any[], Any[], Any[], Any[], Any[], Any[]
quad_solvers = Any[]
quad_soc_solvers = Any[]
rsoc_solvers = Any[]
nlp_solvers = Any[]
convex_nlp_solvers = Any[]
minlp_solvers = Any[]
sdp_solvers = Any[]
global_solvers = Any[]

#lp_solvers = ["cbc", "conopt", "cplex", "minos", "pathnlp", "mosek", "xpress"]
lp_solvers = ["conopt", "minos", "pathnlp", "mosek", "xpress"]

lp_exclude = Dict()

if Sys.islinux() || Sys.iswindows()
    push!(lp_solvers, "xa")
end

quad_solvers = [("conopt", "local"),
#                ("cplex", "nodual"),
                ("minos", "local"),
                ("pathnlp", "local"),
                ("mosek", ""),
                ("xpress", "")]
quad_exclude = Dict( "conopt" => ["socp1", "qcp1"],
#                    "cplex" => ["socp1"],
                    "minos" => ["socp1"],
                    "pathnlp" => ["socp1"],
                    "xpress" => ["qcp1"])
# mosek fails?
mip_solvers = ["xpress", "scip"]#, "cplex", "cbc"]

# Investigate if we can solve those
#soc_solvers = copy(quad_solvers)
soc_solvers = Any[]

push!(convex_nlp_solvers, "")
push!(convex_nlp_solvers, "conopt")
# buggy on my system
#push!(convex_nlp_solvers, "ipopt")
push!(convex_nlp_solvers, "minos")
push!(convex_nlp_solvers, "snopt")
push!(convex_nlp_solvers, "mosek")

append!(convex_nlp_solvers, nlp_solvers)

push!(minlp_solvers, "")

# Default configuration.
const config_mosek = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                               optimal_status=MOI.OPTIMAL,
                               query=false,
                               infeas_certificates=false, # Do not ask for infeasibility certificates.
                               modify_lhs=false)

const config_mosek_noduals = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                                       optimal_status=MOI.OPTIMAL,
                                       query=false,
                                       duals=false,
                                       infeas_certificates=false,
                                       modify_lhs=false)

# Default configuration.
const config_knitro = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                               optimal_status=MOI.LOCALLY_SOLVED,
                               query=false,
                               infeas_certificates=false, # Do not ask for infeasibility certificates.
                               modify_lhs=false)

const config_knitro_noduals = MOIT.TestConfig(atol=1e-5, rtol=1e-8,
                                       optimal_status=MOI.LOCALLY_SOLVED,
                                       query=false,
                                       duals=false,
                                       infeas_certificates=false,
                                       modify_lhs=false)

config_solver = Dict(
"mosek"   => config_mosek,
)

config_solver_noduals = Dict(
"mosek"   => config_mosek_noduals,
)


text:
using Catch22
using StatsBase

catch2 = FeatureSet([StatsBase.mean, StatsBase.std], [:mean, :standard_deviation], [["distribution","centre"], ["distribution","spread"]], ["Arithmetic mean", "Sample standard deviation"])
export catch2

catch24Names = vcat(:mean, :standard_deviation, Catch22.featurenames...)
catch24 = catch2+catch22
export catch24

catch23 = catch24\FeatureSet(FC_LocalSimple_mean1_tauresrat)
catch23Names = getnames(catch23)
export catch23

catch21 = catch22\FeatureSet(FC_LocalSimple_mean1_tauresrat)
catch21Names = getnames(catch21)
export catch21

IntFeatures = FeatureSet([CO_f1ecac,
                        CO_FirstMin_ac, IN_AutoMutualInfoStats_40_gaussian_fmmi,  SB_BinaryStats_diff_longstretch0,  SB_BinaryStats_mean_longstretch1,
                        PD_PeriodicityWang_th0_01,
                        FC_LocalSimple_mean1_tauresrat,
                        CO_Embed2_Dist_tau_d_expfit_meandiff,
                        SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1])
catch24sansInt = catch24\IntFeatures
catch24sansIntnames = getnames(catch24sansInt)
catch22sansInt = catch22\IntFeatures
catch22sansIntnames = getnames(catch22sansInt)
export catch24sansInt, catch22sansInt

function featureRepresentation(X, featureFunc::Function)
    F = featureFunc(X)
end
export featureRepresentation

function NonstationaryProcesses.forcemat(x::AbstractFeatureArray)
    if x isa AbstractFeatureVector
        x = Catch22.featureMatrix(forcemat(Array(x)), Catch22.featureDims(x))
    end
    return x
end

function intersectFeatures(X::AbstractFeatureArray, Y::AbstractFeatureArray)
    # Will be stable in the order of features in X
    fx, fy = Catch22.featureDims(X), Catch22.featureDims(Y)
    if any(.!in.(fx, (fy,)))
        error("The features of Y are neither the same as nor a superset of the features in X")
    end
    fs = intersect(Catch22.featureDims(X), Catch22.featureDims(Y))
    return (X[fs, :], forcemat(Y)[fs, :])
end
function intersectFeatures(X::AbstractFeatureArray, Y::Array)
    # Label Y with the same features as X
    if size(X, 1) != size(Y, 1)
        @error "X and Y do not have the same number of rows"
    end
    fs = Catch22.featureDims(X)
    Y = FeatureMatrix(Y, fs)
    return (X, Y)
end
function intersectFeatures(X, Y, Z, otherarrays...)
    otherarrays = [Y, Z, otherarrays...]
    for i ∈ 1:lastindex(otherarrays)
        X, otherarrays[i] = intersectFeatures(X, otherarrays[i])
    end
    # Then do it again in case we missed any
    for i ∈ 1:lastindex(otherarrays)
        X, otherarrays[i] = intersectFeatures(X, otherarrays[i])
    end
    return (X, otherarrays...)
end
export intersectFeatures


text:
# All the word, context and language model related files

############ MODEL CREATION ##############

# character based model utilities
cembed(m) = m[1]
wchar(m) = m[2]; bchar(m) = m[3];
wforw(m) = m[4]; bforw(m) = m[5];
wback(m) = m[6]; bback(m) = m[7];
wsoft(m) = m[8]; bsoft(m) = m[9];

# To-load uniquely trained word-based language model
makewmodel1(d)=[ d["cembed"],
                 d["char"][1],
                 d["char"][2],
                 d["forw"][1],
                 d["forw"][2],
                 d["back"][1],
                 d["back"][2],
                 d["soft"][1],
                 d["soft"][2] ]

# In order to load the word-model either as gpu or normal Array
function makewmodel(d)
    d1 = makewmodel1(d)
    if gpu() >= 0
        return map(KnetArray, d1)
    else
        return map(Array, d1)
    end
end

############ PREPROCESSING ##############

function maptoint(sentences, v::Vocab)
    MAXWORD = 32
    wdict = empty!(v.idict) # it is already empty ?
    cdict = v.cdict
    unkcid = cdict[v.unkchar]
    words = Vector{Int}[]
    sents = Vector{Int}[]

    maxwordlen = 0; maxsentlen = 0;
    for w in (v.sosword, v.eosword)
        wid = get!(wdict, w, 1+length(wdict))
        word = Array{Int}(length(w)) #Array(Int, length(w))
        wordi = 0 # to check 2 byte characters
        for c in w
            word[wordi+=1] = get(cdict, c, unkcid)
        end
        (wordi != length(w)) && error("Missing in single word process")
        (wordi > maxwordlen) && (maxwordlen = wordi)
        push!(words, word)
    end

    for s in sentences
        sent = Array{Int}(length(s.word))  #Array(Int, length(s.word))
        senti = 0
        for w in s.word
            ndict = length(wdict)
            wid = get!(wdict, w, 1+ndict)
            sent[senti+=1] = wid
            if wid == 1+ndict
                word = Array{Int}(length(w)) #Array(Int, length(w))
                wordi = 0
                for c in w
                    word[wordi+=1] = get(cdict, c, unkcid)
                end
                (wordi != length(w)) && error("Missing in single word process")
                if wordi > MAXWORD; wordi=MAXWORD; word = word[1:wordi]; end;
                (wordi > maxwordlen) && (maxwordlen = wordi)
                push!(words, word)
            end
        end
        @assert(senti == length(s.word))
        (senti > maxsentlen) && (maxsentlen = senti)
        push!(sents, sent)
    end
    @assert(length(wdict) == length(words))
    return words, sents, maxwordlen, maxsentlen
end

# To make ready for character based lstm, data[i] corresponds to ith time step input, to charlstm
function tokenbatch(words, maxlen, sos, eos, pad=eos)
    B = length(words) # batchsize
    T = maxlen + 2
    data = [ Array{Int}(B) for t in 1:T ]
    mask = [ Array{Float32}(B) for t in 1:T ]
    @inbounds for t in 1:T
        for b in 1:B
            N = length(words[b]) # wordlen
            n = t - T + N + 1 # cursor
            if n < 0
                mask[t][b] = 0
                data[t][b] = pad
            else
                mask[t][b] = 1
                if n == 0
                    data[t][b] = sos
                elseif n <= N
                    data[t][b] = words[b][n]
                elseif n == N+1
                    data[t][b] = eos
                else
                    error()
                end
            end
        end
    end
    return data, mask
end

# Gold-batch for final layer of the bilstm implementation
function goldbatch(sentences, maxlen, wdict, unkwid, pad=unkwid)
    B = length(sentences)
    T = maxlen
    data = [ Array{Int}(B) for t in 1:T ]
    mask = [ Array{Float32}(B) for t in 1:T ]
    for t in 1:T
        for b in 1:B
            N = length(sentences[b])
            n = t - T + N
            if n <= 0
                mask[t][b] = 0
                data[t][b] = pad
            else
                mask[t][b] = 1
                data[t][b] = get(wdict, sentences[b].word[n], unkwid)
            end
        end
    end
    return data, mask
end

############# MODEL FORWARDING FUNCTIONALITY ############

# col-major lstm
function _lstm(weight, bias, hidden, cell, input; mask=nothing)
    gates = weight * vcat(input, hidden) .+ bias
    H = size(hidden, 1)
    forget = sigm.(gates[1:H, :])
    ingate = sigm.(gates[1+H:2H, :])
    outgate = sigm.(gates[1+2H:3H, :])
    change = tanh.(gates[1+3H:4H, :])
    (mask != nothing) && (mask = reshape(mask, 1, length(mask)))

    cell = cell .* forget + ingate .* change
    hidden = outgate .* tanh.(cell)
    if mask != nothing
        hidden = hidden .* mask
        cell = cell .* mask
    end
    return (hidden, cell)
end

# Token-batched words are coming as input to the system, col-major lstm
function charlstm(model, data, mask)
    weight, bias, embeddings = wchar(model), bchar(model), cembed(model)
    T = length(data)
    B = length(data[1])
    H = div(length(bias), 4)

    if isa(weight, KnetArray)
        mask = map(KnetArray, mask)
    end

    czero = fill!(similar(bias, H, B), 0)
    hidden = cell = czero
    for t in 1:T
        (hidden, cell) = _lstm(weight, bias, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
    end
    return hidden
end

# col-major bilstm implementation
function wordlstm(model, data, mask, embeddings)
    weight, bias = wforw(model), bforw(model)
    T = length(data)
    B = length(data[1])
    H = div(length(bias), 4)

    if isa(weight, KnetArray)
        mask = map(KnetArray, mask)
    end

    wzero = fill!(similar(bias, H, B), 0)

    # forward lstm
    hidden = cell = wzero
    fhiddens = Array{Any}(T-2)  # fhiddens = Array(Any, T-2) : deprecated
    for t in 1:T-2
        (hidden, cell) = _lstm(weight, bias, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
        fhiddens[t] = hidden
    end

    # backward lstm
    weight_b, bias_b = wback(model), bback(model)
    hidden = cell = wzero
    bhiddens = Array{Any}(T-2)  # bhiddens = Array(Any, T-2) : deprecated
    for t in T:-1:3
        (hidden, cell) = _lstm(weight_b, bias_b, hidden, cell, embeddings[:, data[t]]; mask=mask[t])
        bhiddens[t-2] = hidden
    end
    return fhiddens, bhiddens
end

# col-major loss function implementation, as explained in paper
function lmloss(model, data, mask, forw, back; result=nothing)
    T = length(data)
    B = length(data[1])
    weight, bias = wsoft(model), bsoft(model)
    idx(t,b,n) = data[t][b] + (b-1)*n

    total = count = 0
    for t in 1:T
        ypred = weight * vcat(forw[t], back[t]) .+ bias
        nrows,ncols = size(ypred)
        index = Int[]
        for b=1:B
            if mask[t][b]==1
                push!(index, idx(t,b,nrows))
            end
        end
        o1 = logp(ypred, 1)
        o2 = o1[index]
        total += sum(o2)
        count += length(o2)
    end

    if result != nothing
        result[1] += AutoGrad.getval(total)
        result[2] += count
    end
    return total
end

################# To COMPUTE RELATED EMBEDDINGS ###################

# To fill the constant word embeddings
function fillwvecs!(sentences, isents, wembed; GPUFEATS=false)
    for (s, isents) in zip(sentences, isents)
        empty!(s.wvec)
        for w in isents
            if GPUFEATS
                push!(s.wvec, wembed[:, w])
            else
                push!(s.wvec, Array(wembed[:, w]))
            end
        end
    end
end

# To fill the context embeddings of words in given sentences one by one
function fillcvecs!(sentences, forw, back; GPUFEATS=false)
    T = length(forw)
    for i in 1:length(sentences)
        s = sentences[i]
        empty!(s.fvec)
        empty!(s.bvec)
        N = length(s)
        for n in 1:N
            t = T-N+n
            if GPUFEATS #GPU
                push!(s.fvec, forw[t][:,i])
                push!(s.bvec, back[t][:,i])
            else #CPU
                push!(s.fvec, Array(forw[t][:,i]))
                push!(s.bvec, Array(back[t][:,i]))
            end
        end
    end
end

# fill both word and context embeddings
function fillvecs!(wmodel, sentences, vocab; batchsize=128)

    words, sents, maxwordlen, maxsentlen = maptoint(sentences, vocab)
    sow = vocab.cdict[vocab.sowchar]
    eow = vocab.cdict[vocab.eowchar]
    paw = vocab.cdict[vocab.unkchar]

    # word-embeddings calcutation
    wembed = Any[]
    #free_KnetArray();
    for i=1:batchsize:length(words)
        j = min(i+batchsize-1,length(words))
        wij = view(words,i:j)
        maxij = maximum(map(length, wij))
        cdata, cmask = tokenbatch(wij, maxij, sow, eow)
        push!(wembed, charlstm(wmodel, cdata, cmask))
    end
    wembed =hcat(wembed...) # Here I applied the changes from hcatn-> hcat in newer version
    fillwvecs!(sentences, sents, wembed)

    sos,eos,unk = vocab.idict[vocab.sosword], vocab.idict[vocab.eosword], vocab.odict[vocab.unkword]
    result = zeros(2)
    #free_KnetArray()
    for i=1:batchsize:length(sents)
        j = min(i+batchsize-1, length(sents))
        isentij = view(sents, i:j)
        maxij = maximum(map(length, isentij))
        wdata, wmask = tokenbatch(isentij, maxij, sos, eos)
        forw, back = wordlstm(wmodel, wdata, wmask, wembed)
        sentij = view(sentences, i:j)
        fillcvecs!(sentij, forw, back)
        odata, omask = goldbatch(sentij, maxij, vocab.odict, unk)
        lmloss(wmodel,odata,omask,forw,back; result=result)
    end
    return exp(-result[1]/result[2])
end

# calculate the coverage
function unkrate(sentences)
    total = known = 0
    for s in sentences
        for w in s.word
            total += 1
            known += haskey(s.vocab.odict, w)
        end
    end
    return (total-known,total)
end