text:
"""
Implementation of Proximal Bundle Method.

The implementation is based on
  Krzysztof C. Kiwiel, "Proximity control in bundle methods for convex nondifferentiable minimization"
  Mathematical Programming 46(1-3), 1990
"""

using LinearAlgebra
using SparseArrays
using Printf

const disable_purge_cuts = true

mutable struct ProximalMethod <: AbstractMethod
    model::BundleModel

    y::Array{Float64,1}  # current iterate of dimension n
    fy::Array{Float64,1} # objective values at y for N functions
    g::Dict{Int,SparseVector{Float64}} # subgradients of dimension n for N functions

    cuts::Dict{JuMP.ConstraintRef,Dict{String,Any}}

    iter::Int # iteration counter

    u::Float64
    x0::Array{Float64,1}	# current best solution (at iteration k)
    fx0::Array{Float64,1}	# current best objective values
    d::Array{Float64,1}
    v::Array{Float64,1}
    sum_of_v::Float64
    ϵ_v::Float64
    i::Int
    α::Array{Float64,1}

    statistics::Dict{Any,Any} # arbitrary collection of statistics
    eval_time::Float64        # function evaluation time
    start_time::Float64       # start time

    params::Parameters

    function ProximalMethod(n::Int, N::Int, func;
        init::Array{Float64,1} = zeros(n),
        params::Parameters = Parameters())

        pm = new()
        pm.model = BundleModel(n, N, params.ncuts_per_iter, func)

        @assert length(init) == n

        pm.y = copy(init)
        pm.fy = zeros(N)
        pm.g = Dict()

        pm.cuts = Dict()

        pm.iter = 0

        pm.u = params.u
        pm.x0 = copy(init)
        pm.fx0 = zeros(N)
        pm.d = zeros(n)
        pm.v = zeros(N)
        pm.sum_of_v = 0.0
        pm.ϵ_v = Inf
        pm.i = 0
        pm.α = zeros(N)

        pm.statistics = Dict(
            "total_eval_time" => 0.0)
        pm.start_time = time()

        pm.params = params

        return pm
    end
end

ProximalMethod(n::Int, N::Int, func, init::Array{Float64,1}) = ProximalMethod(n, N, func, init = init)

# This returns BundleModel object.
get_model(method::ProximalMethod)::BundleModel = method.model

# This returns solution.
get_solution(method::ProximalMethod) = method.x0

# This returns objective value.
get_objective_value(method::ProximalMethod) = sum(method.fx0)

# This sets the termination tolerance.
set_bundle_tolerance!(method::ProximalMethod, tol::Float64) = set_parameter(method.params, "ϵ_s", tol)

# This creates an objective function to the bundle model.
function add_objective_function!(method::ProximalMethod)
    bundle = get_model(method)
    d = bundle.model[:x]
    v = bundle.model[:θ]
    @objective(bundle.model, Min,
          sum(v[j] for j = 1:bundle.ncuts_per_iter)
        + 0.5 * method.u * sum(d[i]^2 for i = 1:bundle.n))
end

# This may collect solutions from the bundle model.
function collect_model_solution!(method::ProximalMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    if JuMP.termination_status(model) in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.LOCALLY_SOLVED]
        d = model[:x]
        v = model[:θ]
        for i = 1:bundle.n
            method.d[i] = JuMP.value(d[i])
            method.y[i] = method.x0[i] + method.d[i]
        end
        for j = 1:bundle.ncuts_per_iter
            method.v[j] = JuMP.value(v[j])
        end
        method.sum_of_v = sum(method.v)
        for (ref, cut) in method.cuts
            @assert JuMP.is_valid(get_jump_model(method), ref)
            cut["dual"] = JuMP.dual(ref)
            if cut["dual"] > -1e-6 #&& cut["α"] > 0.0
                cut["age"] += 1
            end
        end
    else
        # JuMP.print(model)
        @error "Unexpected model solution status ($(JuMP.termination_status(model)))"
    end
end

# Should the method terminate?
function termination_test(method::ProximalMethod)
    model = get_jump_model(method)
    if JuMP.termination_status(model) ∉ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        return true
    end
    if method.sum_of_v >= -method.params.ϵ_s * (1 + abs(sum(method.fx0)))
        println("TERMINATION: Optimal: v = ", method.sum_of_v)
        return true
    end
    if method.iter >= method.params.maxiter
        println("TERMINATION: Maximum number of iterations reached.")
        return true
    end
    if sum(method.fx0) <= method.params.obj_limit
        println("TERMINATION: Dual objective limit reached.")
        return true
    end
    return false
end

"""
This calls function `BundleModel.evaluate_f(y)` to get function value `fy` and gradient `g`.
This method assumes user-defined function of the form
    `evaluate_f(y::Vector{Float64})::Tuple{Float64,Array{Float64,2}}`
    returning function evaluation value and gradient as first and second outputs, resp.
"""
function evaluate_functions!(method::ProximalMethod)
    stime = time()
    method.fy, method.g = method.model.evaluate_f(method.y)
    method.statistics["total_eval_time"] += time() - stime

    if method.iter == 0
        method.x0 = copy(method.y)
        method.fx0 = copy(method.fy)
    end
end

# This updates the bundle pool by removing and/or adding bundle objects.
function update_bundles!(method::ProximalMethod)
    purge_bundles!(method)

    bundle = get_model(method)
    model = get_model(bundle)

    sumfy = sum(method.fy)
    sumfx0 = sum(method.fx0)
    if sumfy - sumfx0 <= method.params.m_L * method.sum_of_v
        # update bundles first
        for (ref, cut) in method.cuts
            i::Int = cut["cut_index"]
            g::SparseVector{Float64} = cut["g"]
            offset = g' * method.d
            offset += sum(method.fx0[j] for j in bundle.cut_indices[i])
            offset -= sum(method.fy[j] for j in bundle.cut_indices[i])
            JuMP.add_to_function_constant(ref, offset)
            # @show ref
        end

        d = model[:x]
        for i = 1:bundle.n
            if JuMP.has_upper_bound(d[i])
                JuMP.set_upper_bound(d[i], JuMP.upper_bound(d[i]) + method.x0[i] - method.y[i])
            end
            if JuMP.has_lower_bound(d[i])
                JuMP.set_lower_bound(d[i], JuMP.lower_bound(d[i]) + method.x0[i] - method.y[i])
            end
        end

        # serious step
        method.x0 = copy(method.y)
        method.fx0 = copy(method.fy)
        fill!(method.α, 0.0)
    else
        for j in eachindex(method.α)
            method.α[j] = method.fx0[j] - method.fy[j] + method.g[j]' * method.d
        end
    end

    add_bundles!(method)

    u = copy(method.u)
    if sumfy - sumfx0 <= method.params.m_L * method.sum_of_v
        if method.i > 0 && sumfy - sumfx0 <= method.params.m_R * method.sum_of_v
            u = 2 * method.u * (1 - (sumfy - sumfx0) / method.sum_of_v)
        elseif method.i > 3
            u = method.u / 2
        end
        # @show u, method.u/10, method.u_min
        newu = max(u, method.u / 10, method.params.u_min)
        method.ϵ_v = max(method.ϵ_v, -2 * method.sum_of_v)
        method.i = max(method.i + 1, 1)
        if newu != method.u
            method.u = newu
            method.i = 1
        end
        update_objective!(method)
    else
        p = norm(-method.u .* method.d, 1)
        α_tilde = -p^2 / method.u - method.sum_of_v

        method.ϵ_v = min(method.ϵ_v, p + α_tilde)
        if sum(method.α) > max(method.ϵ_v, -10 * method.sum_of_v) && method.i < -3
            u = 2 * method.u * (1 - (sumfy - sumfx0) / method.sum_of_v)
        end
        # @show u, 10 * method.u
        newu = min(u, 10 * method.u)
        method.i = min(method.i - 1, -1)
        if newu != method.u
            method.u = newu
            method.i = -1
            update_objective!(method)
        end
    end
    # update_objective!(method)
end

function purge_bundles!(method::ProximalMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    ncuts_removed = 0
    ncols = JuMP.num_variables(model)
    ncuts = length(method.cuts)
    for (ref, cut) in method.cuts
        if ncuts - ncuts_removed <= ncols
            break
        end
        if cut["age"] >= method.params.max_age
            JuMP.delete(model, ref)
            delete!(method.cuts, ref)
            ncuts_removed += 1
        end
    end
    if ncuts_removed > 0
        @printf("Removed %d inactive cuts.\n", ncuts_removed)
    end
end

function add_bundles!(method::ProximalMethod)
    # add bundles as constraints to the model
    bundle = get_model(method)
    model = get_model(bundle)
    for i = 1:bundle.ncuts_per_iter
        cut_indices = bundle.cut_indices[i]
        g = sum(method.g[j] for j in cut_indices)
        α = sum(method.α[j] for j in cut_indices)

        cut_violation = -α + g' * method.d - method.v[i]
        if method.iter == 0 || cut_violation > method.params.ϵ_float
            d = model[:x]
            v = model[:θ][i]
            ref = @constraint(model, sum(g[k] * d[k] for k = 1:bundle.n) - v <= α)
            method.cuts[ref] = Dict(
                "age" => 0,
                "dual" => 0.0,
                "cut_index" => i,
                "g" => g,
                "α" => α
            )
            # @show ref
        end
    end
end

# This updates the iteration information.
function update_iteration!(method::ProximalMethod)
    method.iter += 1
end

# This displays iteration information.
function display_info!(method::ProximalMethod)
    model = get_jump_model(method)
    nrows = 0
    for tp in [MOI.LessThan{Float64}, MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}]
        nrows += num_constraints(model, AffExpr, tp)
    end
    fx0 = sum(method.fx0)
    @printf("Iter %4d: ncols %5d, nrows %5d, fx0 %+e, fy %+e, m %+e, v %e, u %e, i %+d, master time %6.1fs, eval time %6.1fs, time %6.1fs\n",
        method.iter, num_variables(model), nrows,
        fx0,
        sum(method.fy),
        method.sum_of_v + fx0,
        method.sum_of_v,
        method.u, method.i,
        sum(method.model.time), method.statistics["total_eval_time"], time() - method.start_time)
end

"""
The following functions are specific for this method only.
"""

function update_objective!(method::ProximalMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    d = model[:x]
    v = model[:θ]
    @objective(model, Min,
          sum(v[j] for j = 1:bundle.ncuts_per_iter)
        + 0.5 * method.u * sum(d[i]^2 for i = 1:bundle.n))
end


text:
module Feather

if !isdefined(Core, :String)
    typealias String UTF8String
end
if !isdefined(Base, :view)
    view = sub
end
if Base.VERSION < v"0.5.0-dev+4631"
    unsafe_wrap{A<:Array}(::Type{A}, ptr, len) = pointer_to_array(ptr, len)
    unsafe_string(ptr, len) = utf8(ptr, len)
end

using FlatBuffers, DataStreams, DataFrames, NullableArrays, CategoricalArrays, WeakRefStrings

export Data, DataFrame

# sync with feather
const VERSION = 2

# Arrow type definitions
include("Arrow.jl")

# flatbuffer defintions
include("metadata.jl")

# wesm/feather/cpp/src/common.h
const FEATHER_MAGIC_BYTES = "FEA1".data

ceil_byte(size) = (size + 7) & ~7
bytes_for_bits(size) = div(((size + 7) & ~7), 8)
const BITMASK = UInt8[1, 2, 4, 8, 16, 32, 64, 128]
getbit(byte::UInt8, i) = (byte & BITMASK[i]) == 0
const ALIGNMENT = 8
paddedlength(x) = div((x + ALIGNMENT - 1), ALIGNMENT) * ALIGNMENT
getoutputlength(version, x) = version < VERSION ? x : paddedlength(x)

function writepadded(io, x)
    bw = Base.write(io, x)
    diff = paddedlength(bw) - bw
    Base.write(io, zeros(UInt8, diff))
    return bw + diff
end

juliastoragetype(meta::Void, values_type) = Type_2julia[values_type]
function juliastoragetype(meta::Metadata.CategoryMetadata, values_type)
    R = Type_2julia[values_type]
    return CategoricalValue{String,R}
end
juliastoragetype(meta::Metadata.TimestampMetadata, values_type) = Arrow.Timestamp{TimeUnit2julia[meta.unit],meta.timezone == "" ? :UTC : Symbol(meta.timezone)}
juliastoragetype(meta::Metadata.DateMetadata, values_type) = Arrow.Date
juliastoragetype(meta::Metadata.TimeMetadata, values_type) = Arrow.Time{TimeUnit2julia[meta.unit]}

juliatype{T<:Arrow.Timestamp}(::Type{T}) = DateTime
juliatype{T<:Arrow.Date}(::Type{T}) = Date
juliatype{T}(::Type{T}) = T
juliatype(::Type{Arrow.Bool}) = Bool

addlevels!{T}(::Type{T}, levels, orders, i, meta, values_type, data, version) = return
function addlevels!{T <: CategoricalValue}(::Type{T}, catlevels, orders, i, meta, values_type, data, version)
    ptr = pointer(data) + meta.levels.offset
    offsets = unsafe_wrap(Array, convert(Ptr{Int32}, ptr), meta.levels.length + 1)
    ptr += getoutputlength(version, sizeof(offsets))
    catlevels[i] = map(x->unsafe_string(ptr + offsets[x], offsets[x+1] - offsets[x]), 1:meta.levels.length)
    orders[i] = meta.ordered
    return
end

# DataStreams interface types
type Source <: Data.Source
    path::String
    schema::Data.Schema
    ctable::Metadata.CTable
    data::Vector{UInt8}
    feathertypes::Vector{DataType} # separate from the types in schema, since we need to convert between feather storage types & julia types
    levels::Dict{Int,Vector{String}}
    orders::Dict{Int,Bool}
    columns::Vector{Any} # holds references to pre-fetched columns for Data.getfield
end

# reading feather files
function Source(file::AbstractString)
    # validity checks
    isfile(file) || throw(ArgumentError("'$file' is not a valid file"))
    m = Mmap.mmap(file)
    length(m) < 12 && throw(ArgumentError("'$file' is not in the feather format"))
    (m[1:4] == FEATHER_MAGIC_BYTES && m[end-3:end] == FEATHER_MAGIC_BYTES) ||
        throw(ArgumentError("'$file' is not in the feather format"))
    # read file metadata using FlatBuffers
    metalength = Base.read(IOBuffer(m[length(m)-7:length(m)-4]), Int32)
    metapos = length(m) - (metalength + 7)
    rootpos = Base.read(IOBuffer(m[metapos:metapos+4]), Int32)
    ctable = FlatBuffers.read(Metadata.CTable, m, metapos + rootpos - 1)
    ctable.version < VERSION && warn("This Feather file is old and will not be readable beyond the 0.3.0 release")
    header = String[]
    types = DataType[]
    juliatypes = DataType[]
    columns = ctable.columns
    levels = Dict{Int,Vector{String}}()
    orders = Dict{Int,Bool}()
    for (i, col) in enumerate(columns)
        push!(header, col.name)
        push!(types, juliastoragetype(col.metadata, col.values.type_))
        jl = juliatype(types[end])
        addlevels!(jl, levels, orders, i, col.metadata, col.values.type_, m, ctable.version)
        push!(juliatypes, col.values.null_count == 0 ? (jl <: AbstractString ? String : jl) : Nullable{jl})
    end
    # construct Data.Schema and Feather.Source
    return Source(file, Data.Schema(header, juliatypes, ctable.num_rows), ctable, m, types, levels, orders, Array{Any}(length(columns)))
end

# DataStreams interface
function Data.isdone(io::Feather.Source, row, col)
    rows, cols = size(Data.schema(io))
    return col > cols || row > rows
end
Data.streamtype{T<:Feather.Source}(::Type{T}, ::Type{Data.Column}) = true
Data.streamtype{T<:Feather.Source}(::Type{T}, ::Type{Data.Field}) = true

function Data.getfield{T}(source::Source, ::Type{T}, row, col)
    !isdefined(source.columns, col) && (source.columns[col] = Data.getcolumn(source, T, col))
    return source.columns[col][row]::T
end

checknonull(source, col) = source.ctable.columns[col].values.null_count > 0 && throw(NullException)
getbools(s::Source, col) = Bool[getbit(s.data[s.ctable.columns[col].values.offset + bytes_for_bits(x)], mod1(x, 8)) for x = 1:s.ctable.num_rows]

function unwrap{T}(s::Source, ::Type{T}, col, rows, off=0)
    bitmask_bytes = s.ctable.columns[col].values.null_count > 0 ? bytes_for_bits(rows) : 0
    return unsafe_wrap(Array, convert(Ptr{T}, pointer(s.data) + s.ctable.columns[col].values.offset + getoutputlength(s.ctable.version, bitmask_bytes) + off), rows)::Vector{T}
end

transform!{T}(::Type{T}, A, len) = T[x for x in A]
transform!(::Type{Date}, A, len) = map(x->Arrow.unix2date(x), A)
transform!{P,Z}(::Type{DateTime}, A::Vector{Arrow.Timestamp{P,Z}}, len) = map(x->Arrow.unix2datetime(P, x), A)
transform!{T,R}(::Type{CategoricalValue{T,R}}, A, len) = map(x->x + R(1), A)
function transform!(::Type{Bool}, A, len)
    B = falses(len)
    Base.copy_chunks!(B.chunks, 1, map(x->x.value, A), 1, length(A) * 64)
    return convert(Vector{Bool}, B)
end

function Data.getcolumn{T}(source::Source, ::Type{T}, col)
    checknonull(source, col)
    A = unwrap(source, source.feathertypes[col], col, source.ctable.num_rows)
    return transform!(T, A, source.ctable.num_rows)::Vector{T}
end
function Data.getcolumn(source::Source, ::Type{Bool}, col)
    checknonull(source, col)
    A = unwrap(source, source.feathertypes[col], col, max(1,div(bytes_for_bits(source.ctable.num_rows),8)))
    return transform!(Bool, A, source.ctable.num_rows)::Vector{Bool}
end
function Data.getcolumn{T <: AbstractString}(source::Source, ::Type{T}, col)
    checknonull(source, col)
    offsets = unwrap(source, Int32, col, source.ctable.num_rows + 1)
    values = unwrap(source, UInt8, col, offsets[end], getoutputlength(source.ctable.version, sizeof(offsets)))
    return [unsafe_string(pointer(values, offsets[i]+1), Int(offsets[i+1] - offsets[i])) for i = 1:source.ctable.num_rows]
end
function Data.getcolumn{T}(source::Source, ::Type{Nullable{T}}, col)
    A = transform!(T, unwrap(source, source.feathertypes[col], col, source.ctable.num_rows), source.ctable.num_rows)::Vector{T}
    bools = getbools(source, col)
    return NullableArray{T,1}(A, bools, source.data)
end
function Data.getcolumn{T <: AbstractString}(source::Source, ::Type{Nullable{T}}, col)
    offsets = unwrap(source, Int32, col, source.ctable.num_rows + 1)
    values = unwrap(source, UInt8, col, offsets[end], sizeof(offsets))
    A = [WeakRefString(pointer(values, offsets[i]+1), Int(offsets[i+1] - offsets[i])) for i = 1:source.ctable.num_rows]
    bools = getbools(source, col)
    return NullableArray{WeakRefString{UInt8},1}(A, bools, source.data)
end
function Data.getcolumn{T,R}(source::Source, ::Type{CategoricalValue{T,R}}, col)
    checknonull(source, col)
    refs = transform!(CategoricalValue{T,R}, unwrap(source, R, col, source.ctable.num_rows), source.ctable.num_rows)
    pool = CategoricalPool{String, R}(source.levels[col], source.orders[col])
    return CategoricalArray{String,1,R}(refs, pool)
end
function Data.getcolumn{T,R}(source::Source, ::Type{Nullable{CategoricalValue{T,R}}}, col)
    refs = transform!(CategoricalValue{T,R}, unwrap(source, R, col, source.ctable.num_rows), source.ctable.num_rows)
    bools = getbools(source, col)
    refs = R[ifelse(bools[i], R(0), refs[i]) for i = 1:source.ctable.num_rows]
    pool = CategoricalPool{String, R}(source.levels[col], source.orders[col])
    return NullableCategoricalArray{String,1,R}(refs, pool)
end

"""
`Feather.read{T <: Data.Sink}(file, sink_type::Type{T}, sink_args...)` => `T`
`Feather.read(file, sink::Data.Sink)` => `Data.Sink`

`Feather.read` takes a feather-formatted binary `file` argument and "streams" the data to the
provided `sink` argument. A fully constructed `sink` can be provided as the 2nd argument (the 2nd method above),
or a Sink can be constructed "on the fly" by providing the type of Sink and any necessary positional arguments
(the 1st method above). By default, a `DataFrame` is returned.

Examples:

```julia
# default read method, returns a DataFrame
df = Feather.read("cool_feather_file.feather")

# read a feather file directly into a SQLite database table
db = SQLite.DB()
Feather.read("cool_feather_file.feather", SQLite.Sink, db, "cool_feather_table")
```
"""
function read end

function read(file::AbstractString, sink=DataFrame, args...; append::Bool=false)
    source = Source(file)
    sink = Data.stream!(source, sink, append, args...)
    Data.close!(sink)
    return sink
end

function read{T}(file::AbstractString, sink::T; append::Bool=false)
    source = Source(file)
    sink = Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end

read(source::Feather.Source, sink=DataFrame, args...; append::Bool=false) = (sink = Data.stream!(source, sink, append, args...); Data.close!(sink); return sink)
read{T}(source::Feather.Source, sink::T; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)

# writing feather files
feathertype{T}(::Type{T}) = Feather.julia2Type_[T]
feathertype{S,R}(::Type{CategoricalValue{S,R}}) = julia2Type_[R]
feathertype{P}(::Type{Arrow.Time{P}}) = Metadata.INT64
feathertype(::Type{Date}) = Metadata.INT32
feathertype(::Type{DateTime}) = Metadata.INT64
feathertype{T<:AbstractString}(::Type{T}) = Metadata.UTF8

getmetadata{T}(io, ::Type{T}, A) = nothing
getmetadata(io, ::Type{Date}, A) = Metadata.DateMetadata()
getmetadata{T}(io, ::Type{Arrow.Time{T}}, A) = Metadata.TimeMetadata(julia2TimeUnit[T])
getmetadata(io, ::Type{DateTime}, A) = Metadata.TimestampMetadata(julia2TimeUnit[Arrow.Millisecond], "")
function getmetadata{S,R}(io, T::Type{CategoricalValue{S,R}}, A)
    lvls = CategoricalArrays.levels(A)
    len = length(lvls)
    offsets = zeros(Int32, len+1)
    offsets[1] = off = 0
    for (i,v) in enumerate(lvls)
        off += length(v)
        offsets[i + 1] = off
    end
    offset = position(io)
    total_bytes = writepadded(io, view(reinterpret(UInt8, offsets), 1:(sizeof(Int32) * (len + 1))))
    total_bytes += writepadded(io, lvls)
    return Metadata.CategoryMetadata(Metadata.PrimitiveArray(Metadata.UTF8, Metadata.PLAIN, offset, len, 0, total_bytes), ordered(A))
end

values(A::Vector) = A
values(A::NullableVector) = A.values
values{S,R}(A::Union{CategoricalArray{S,1,R},NullableCategoricalArray{S,1,R}}) = map(x-> x - R(1), A.refs)

# Bool
function writecolumn(io, ::Type{Bool}, o, b, A)
    return writepadded(io, view(reinterpret(UInt8, convert(BitVector, values(A)).chunks), 1:bytes_for_bits(length(A))))
end
# Category
function writecolumn{S,R}(io, ::Type{CategoricalValue{S,R}}, o, b, A)
    return writepadded(io, view(reinterpret(UInt8, values(A)), 1:(length(A) * sizeof(R))))
end
# Date
function writecolumn(io, ::Type{Date}, o, b, A)
    return writepadded(io, view(reinterpret(UInt8, map(Arrow.date2unix, values(A))), 1:(length(A) * sizeof(Int32))))
end
# Timestamp
function writecolumn(io, ::Type{DateTime}, o, b, A)
    return writepadded(io, view(reinterpret(UInt8, map(Arrow.datetime2unix, values(A))), 1:(length(A) * sizeof(Int64))))
end
# Date, Timestamp, Time and other primitive T
function writecolumn{T}(io, ::Type{T}, o, b, A)
    return writepadded(io, view(reinterpret(UInt8, values(A)), 1:(length(A) * sizeof(T))))
end
# List types
valuelength{T}(val::T) = length(String(val))
valuelength{T}(val::Nullable{T}) = isnull(val) ? 0 : length(get(val))

writevalue{T}(io, val::T) = Base.write(io, String(val).data)
writevalue{T}(io, val::Nullable{T}) = isnull(val) ? 0 : Base.write(io, String(get(val)).data)

function writecolumn{T<:Union{Vector{UInt8},AbstractString}}(io, ::Type{T}, offsets, writeoffset, arr)
    len = length(arr)
    off = isempty(offsets) ? 0 : offsets[end]
    ind = isempty(offsets) ? 1 : length(offsets)
    # total_bytes = sizeof(Int32) * (isempty(offsets) ? len + 1 : len)
    append!(offsets, zeros(Int32, isempty(offsets) ? len + 1 : len))
    offsets[1] = isempty(offsets) ? 0 : offsets[1]
    for v in arr
        off += valuelength(v)
        offsets[ind + 1] = off
        ind += 1
    end
    writeoffset && (total_bytes = writepadded(io, view(reinterpret(UInt8, offsets), 1:length(offsets) * sizeof(Int32))))
    total_bytes += offsets[len+1]
    for val in arr
        writevalue(io, val)
    end
    diff = paddedlength(offsets[len+1]) - offsets[len+1]
    if diff > 0
        Base.write(io, zeros(UInt8, diff))
        total_bytes += diff
    end
    return total_bytes
end

writenulls(io, A, null_count, len, total_bytes) = return total_bytes
function writenulls(io, A::NullableVector, null_count, len, total_bytes)
    # write out null bitmask
    if null_count > 0
        null_bytes = Feather.bytes_for_bits(len)
        bytes = BitArray(!A.isnull)
        total_bytes = writepadded(io, view(reinterpret(UInt8, bytes.chunks), 1:null_bytes))
    end
    return total_bytes
end
function writenulls{T <: NullableCategoricalArray}(io, A::T, null_count, len, total_bytes)
    # write out null bitmask
    if null_count > 0
        null_bytes = Feather.bytes_for_bits(len)
        bytes = BitArray(!(A.refs .== 0))
        total_bytes = writepadded(io, view(reinterpret(UInt8, bytes.chunks), 1:null_bytes))
    end
    return total_bytes
end

"DataStreams Sink implementation for feather-formatted binary files"
type Sink <: Data.Sink
    schema::Data.Schema
    ctable::Metadata.CTable
    io::IO
    description::String
    metadata::String
    df::DataFrame
end

# construct a Sink from an existing Feather Source
# function Sink(s::Feather.Source; append::Bool=false)
#     data = open(s.path, append ? "a" : "w")
#     seek(data, append ? filesize(s.path) : s.datapos)
#     return Sink(s.schema, s.options, data, s.datapos, !append)
# end

function Sink{T<:Data.StreamType}(file::Union{IO,AbstractString}, schema::Data.Schema=Data.EMPTYSCHEMA, ::Type{T}=Data.Column;
              description::AbstractString=String(""), metadata::AbstractString=String(""), append::Bool=false)
    df = DataFrame(schema, T)
    if isa(file, AbstractString)
        if append && isfile(file)
            df = Feather.read(file)
            schema.rows += size(df, 1)
        end
        io = open(file, "w")
    else
        if append
            df = Feather.read(file)
            schema.rows += size(df, 1)
        end
        io = seekstart(file)
    end
    writepadded(io, FEATHER_MAGIC_BYTES)
    return Sink(schema, Metadata.CTable("", 0, Metadata.Column[], VERSION, ""), io, description, metadata, df)
end

Base.close(s::Sink) = (applicable(close, s.io) && close(s.io); return nothing)
Base.flush(s::Sink) = (applicable(flush, s.io) && flush(s.io); return nothing)

# DataStreams interface
function Sink{T}(source, ::Type{T}, append::Bool, file::AbstractString)
    sink = Sink(file, Data.schema(source), T; append=append)
    return sink
end

function Sink{T}(sink, source, ::Type{T}, append::Bool)
    if !append
        for col in sink.df.columns
            empty!(col)
        end
    end
    sink.schema = Data.schema(source)
    sink.df = DataFrame(sink.schema, T)
    return sink
end

Data.streamtypes{T<:Feather.Sink}(::Type{T}) = [Data.Column, Data.Field]

Data.streamfield!{T}(sink::Feather.Sink, source, ::Type{T}, row, col, cols, sinkrows) = Data.streamfield!(sink.df, source, T, row, col, cols, sinkrows)

Data.streamcolumn!{T}(sink::Feather.Sink, source, ::Type{T}, col, row) = Data.streamcolumn!(sink.df, source, T, col, row)

function Data.close!(sink::Feather.Sink)
    sch = Data.schema(sink)
    header = Data.header(sch)
    types = Data.types(sch)
    data = sink.df
    io = sink.io
    # write out arrays, building each array's metadata as we go
    rows = size(sch, 1)
    columns = Metadata.Column[]
    for (i, name) in enumerate(header)
        @inbounds T = types[i]
        arr = data[i]
        total_bytes = 0
        offset = position(io)
        null_count = Data.nullcount(arr)
        len = length(arr)
        total_bytes = writenulls(io, arr, null_count, len, total_bytes)
        # write out array values
        TT = eltype(arr) <: Nullable ? eltype(eltype(arr)) : eltype(arr)
        total_bytes += writecolumn(io, TT, Int32[], true, arr)
        values = Metadata.PrimitiveArray(feathertype(TT), Metadata.PLAIN, offset, len, null_count, total_bytes)
        push!(columns, Metadata.Column(String(name), values, getmetadata(io, TT, arr), String("")))
    end
    # write out metadata
    ctable = Metadata.CTable(sink.description, rows, columns, VERSION, sink.metadata)
    meta = FlatBuffers.build!(ctable)
    rng = (meta.head + 1):length(meta.bytes)
    writepadded(io, view(meta.bytes, rng))
    # write out metadata size
    Base.write(io, Int32(length(rng)))
    # write out final magic bytes
    Base.write(io, FEATHER_MAGIC_BYTES)
    close(io)
    sink.ctable = ctable
    return sink
end

"""
`Feather.write{T <: Data.Source}(io, source::Type{T}, source_args...)` => `Feather.Sink`
`Feather.write(io, source::Data.Source)` => `Feather.Sink`

Write a `Data.Source` out to disk as a feather-formatted binary file. The two methods allow the passing of a
fully constructed `Data.Source` (2nd method), or the type of Source and any necessary positional arguments (1st method).

Examples:

```julia
df = DataFrame(...)
Feather.write("shiny_new_feather_file.feather", df)

Feather.write("sqlite_query_result.feather", SQLite.Source, db, "select * from cool_table")
```
"""
function write end

function write{T}(io::Union{AbstractString,IO}, ::Type{T}, args...; append::Bool=false,
                    description::AbstractString=String(""), metadata::AbstractString=String(""))
    source = T(args...)
    sink = Sink(io; description=description, metadata=metadata, append=append)
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end
function write(io::Union{AbstractString,IO}, source; append::Bool=false,
                    description::AbstractString=String(""), metadata::AbstractString=String(""))
    sink = Sink(io; description=description, metadata=metadata, append=append)
    Data.stream!(source, sink, append)
    Data.close!(sink)
    return sink
end

write{T}(sink::Sink, ::Type{T}, args...; append::Bool=false) = (sink = Data.stream!(T(args...), sink, append); Data.close!(sink); return sink)
write(sink::Sink, source; append::Bool=false) = (sink = Data.stream!(source, sink, append); Data.close!(sink); return sink)

end # module

ext: jl


text:
using Flatten
using Base.Test

type Foo{T}
    a::T
    b::T
    c::T
end

type Nested{T1, T2}
    f::Foo{T1}
    b::T2
    c::T2
end

foo = Foo(1.0, 2.0, 3.0)
nested = Nested(Foo(1,2,3), 4.0, 5.0)

@test to_vector(from_vector(Foo, to_vector(foo))) == to_vector(foo)
@test to_tuple(from_tuple(Foo, to_tuple(foo))) == to_tuple(foo)
@test to_tuple(Nested(Foo(1,2,3),4,5)) == (1,2,3,4,5)
@test to_tuple((Nested(Foo(1,2,3),4,5), Nested(Foo(6,7,8), 9, 10))) == (1,2,3,4,5,6,7,8,9,10)
@test to_tuple(Nested(Foo(1,2,3), (4,5), (6,7))) == (1,2,3,4,5,6,7)
@test to_vector(Foo(1,2,3)) == Int[1,2,3]
@test typeof(to_vector(Foo(1,2,3))) == Array{Int, 1}

# Test nested types and tuples
@test to_vector((Nested(Foo(1,2,3),4.0,5.0), Nested(Foo(6,7,8), 9, 10))) == Float64[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
@test typeof(to_vector((Nested(Foo(1,2,3),4.0,5.0), Nested(Foo(6,7,8), 9, 10)))) == Array{Float64, 1}
@test to_tuple((Nested(Foo(1,2,3),4.0,5.0), Nested(Foo(6,7,8), 9, 10))) == (1,2,3,4.0,5.0,6,7,8,9,10)
@test to_tuple(from_tuple(Nested, to_tuple(nested))) == to_tuple(nested)
@test to_vector(from_vector(Nested, to_vector(nested))) == to_vector(nested)
@test to_tuple(from_tuple(Tuple{Nested, Nested}, to_tuple((nested, nested)))) == to_tuple((nested, nested))

# Test non-parametric types
type AnyPoint
    x
    y
end
@test to_tuple(AnyPoint(1,2)) == (1,2)
@test to_tuple(from_tuple(AnyPoint, (1,2))) == (1,2)

# Test function wrapping
type Point
	x
	y
end

function distance(p::Point)
	sqrt(p.x^2 + p.y^2)
end

wrapped_distance = wrap(distance, Point)
@test wrapped_distance([1,2]) == [norm([1,2])]

# Test performance
function to_vector_naive(obj)
    v = Array(Float64, length(fieldnames(obj)))
    for (i, field) in enumerate(fieldnames(obj))
        v[i] = getfield(obj, field)
    end
    v
end

function to_tuple_naive(obj)
    v = (map(field -> getfield(obj, field), fieldnames(obj))...)
end

function from_vector_naive(T, data)
    T(data...)
end

@test to_vector_naive(foo) == to_vector(foo)
@test to_tuple_naive(foo) == to_tuple(foo)

function test_to_vector()
    foo = Foo(1.0, 2.0, 3.0)
    for i = 1:1e3
        data = to_vector(foo)
    end
end

function test_to_vector_naive()
    foo = Foo(1.0, 2.0, 3.0)
    for i = 1:1e3
        data = to_vector_naive(foo)
    end
end

function test_to_tuple()
    foo = Foo(1.0, 2.0, 3.0)
    for i = 1:1e3
        data = to_tuple(foo)
    end
end

function test_to_tuple_naive()
    foo = Foo(1.0, 2.0, 3.0)
    for i = 1:1e3
        data = to_tuple_naive(foo)
    end
end

function test_from_vector()
    foo = Foo(1.0, 2.0, 3.0)
    data = to_vector(foo)
    for i = 1:1e3
        foo2 = from_vector(Foo, data)
    end
end

function test_from_vector_naive()
    foo = Foo(1.0, 2.0, 3.0)
    data = to_vector(foo)
    for i = 1:1e3
        foo2 = from_vector_naive(Foo, data)
    end
end

gc_enable(false)
println("to vector")
test_to_vector()
@time test_to_vector()
println("to vector naive")
test_to_vector_naive()
@time test_to_vector_naive()
println("to tuple")
test_to_tuple()
@time test_to_tuple()
println("to tuple naive")
test_to_tuple_naive()
@time test_to_tuple_naive()
println("from vector")
test_from_vector()
@time test_from_vector()
println("from vector naive")
test_from_vector_naive()
@time test_from_vector_naive()
gc_enable(true)


text:
using LinearAlgebra

Cep = Matrix(0.5I, 2, 2);                                                           # quadratic position cost matrix
Cu = Matrix(0.2I, 2, 2);                                                            # quadratic control cost matrix
β_pos = 0.1;                                                                        # relative weight between instant and terminal pos cost
α_col = 100.0;                                                                      # scale parameter for exponential collision cost
β_col = 0.1;                                                                        # relative weight between instant and terminal pos cost
λ_col = 0.2;                                                                        # bandwidth parameter for exponential collision cost
σ_risk = 0.0;                                                                       # risk-sensitiveness parameter


text:
using GitCommand
using Test

function with_temp_dir(f::Function)
    original_directory = pwd()
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)
    result = f(tmp_dir)
    cd(original_directory)
    rm(tmp_dir; force = true, recursive = true)
    return result
end

@testset "GitCommand.jl" begin
    if Sys.iswindows()
        @test GitCommand._separator() == ';'
    else
        @test GitCommand._separator() == ':'
    end
    with_temp_dir() do tmp_dir
        @test !isdir("General")
        @test !isfile(joinpath("General", "Registry.toml"))
        git() do git
            run(`$git clone https://github.com/JuliaRegistries/General`)
        end
        @test isdir("General")
        @test isfile(joinpath("General", "Registry.toml"))
    end
    with_temp_dir() do tmp_dir
        @test !isdir("General")
        @test !isfile(joinpath("General", "Registry.toml"))
        cmd = git`clone https://github.com/JuliaRegistries/General`
        @test cmd isa Cmd
        run(cmd)
        @test isdir("General")
        @test isfile(joinpath("General", "Registry.toml"))
    end
end


text:
poll_events{T <: BaseExperiment{NullWindow}}(callback,exp::T,time::Float64) = nothing

const SDL_KEYDOWN = 0x00000300
const SDL_KEYUP = 0x00000301
const SDL_WINDOWEVENT = 0x00000200
const SDL_QUIT = 0x00000100

const SDL_WINDOWEVENT_FOCUS_GAINED = 0x0000000c
const SDL_WINDOWEVENT_FOCUS_LOST = 0x0000000d

const type_ptr = 0x0000000000000000       # offsetof(SDL_Event,type)
const keysym_ptr = 0x0000000000000010     # offsetof(SDL_KeyboardEvent,keysym)
const sym_ptr = 0x0000000000000004        # offsetof(SDL_Keysym,sym)
const mod_ptr = 0x0000000000000008        # offsetof(SDL_Keysym,mod)
const win_event_ptr = 0x000000000000000c  # offsetof(SDL_WindowEvent,event)
const event_size = 0x0000000000000038     # sizeof(SDL_Event)

function poll_events(callback,exp::ExtendedExperiment,time::Float64)
  poll_events(callback,next(exp),time)
end

"""
     Weber.poll_events(callback,experiment,time)

Call the function `callback`, possibility multiple times, passing it an event
object each time. The time at which the events are polled is passed,
allowing this time to be stored with the event.

!!! warning

    This function should never be called directly by user code. A new
    method of this function can be implemented to extend Weber,
    allowing it to report new kinds events.

"""
function poll_events{T <: BaseExperiment{SDLWindow}}(callback,exp::T,time::Float64)
  event_bytes = Array{Int8}(event_size)
  event = reinterpret(Ptr{Void},pointer(event_bytes))

  while ccall((:SDL_PollEvent,weber_SDL2),Cint,(Ptr{Void},),event) != 0
    etype = at(event,UInt32,type_ptr)
    if etype == SDL_KEYDOWN
      code = at(event,UInt32,keysym_ptr + sym_ptr)
      mod = at(event,UInt16,keysym_ptr + mod_ptr)
      callback(exp,KeyDownEvent(code,mod,time))
    elseif etype == SDL_KEYUP
      code = at(event,UInt32,keysym_ptr + sym_ptr)
      mod = at(event,UInt16,keysym_ptr + mod_ptr)
      callback(exp,KeyUpEvent(code,mod,time))
    elseif etype == SDL_WINDOWEVENT
      wevent = at(event,UInt8,win_event_ptr)
      if wevent == SDL_WINDOWEVENT_FOCUS_GAINED
        callback(exp,WindowFocused(time))
      elseif wevent == SDL_WINDOWEVENT_FOCUS_LOST
        callback(exp,WindowUnfocused(time))
      end
    elseif etype == SDL_QUIT
      callback(exp,QuitEvent())
    end
  end
end


text:
module SampleProjectTest

using BenchmarkTools
using Statistics
using Printf
using Test

export @bench, @onlybench,
       displaylevel, set_displaylevel,
       set_benchmark_duration_seconds, enable_benchmark

displaylevel = 1
set_displaylevel(level) = global displaylevel = level

benchmark_duration_seconds = 1.0
set_benchmark_duration_seconds(duration) = global benchmark_duration_seconds = duration

run_benchmark = false
enable_benchmark(enabled = true) = global run_benchmark = enabled

benchmark_results = Dict()

"""
This macro can be used in front of `@testset` macros. It will only take effect if the
`run_benchmark` global variable is set to `true`. If so, it will do a simple
benchmark of the code within the `@testset`, and store the results within this module,
which can later be reported by calling `save_benchmark`.

Since all benchmark results are aggregated hierarchically, it is advised not to nest
`@bench` expressions, since that will end up double counting elapsed time and memory.
"""
macro bench(ex)
    if run_benchmark
        test = ex.args[end]
        if test.head == :for
            # rewrite "@bench @testset for ... end" as "for ... @bench @testset begin ... end; end"
            rewritten_expr =
                Expr(test.head, test.args[1],
                    Expr(:macrocall, Symbol("@bench"), LineNumberNode(0),
                        Expr(ex.head, ex.args[1:end-1]..., test.args[2])))
            return :($(esc(rewritten_expr)))
        else
            :(benchit(() -> $(esc(ex))))
        end
    else
        :($(esc(ex)))
    end
end

"Like `@bench`, but if `run_benchmark` is false, the `@testset` won't be executed at all."
macro onlybench(ex)
    if run_benchmark
        :(@bench($(esc(ex))))
    end
end

function benchit(f)
    test = Vector(undef, 1)

    warmup(@benchmarkable $test[1] = $f())

    # "Infinite" number of samples, to enforce the configured benchmark duration.
    # There will be a single eval per sample, which at least should be fine for benchmark
    # expressions longer than ~1 μs.
    benchmark_results[test[1]] = run((@benchmarkable $f()), samples = typemax(Int), seconds = benchmark_duration_seconds)
end

end

max_stars_repo_name: maxbennedich/julia-regression-analysis

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
  export seedlink, seedlink!, sl_info, has_sta, has_stream

# ========================================================================
# Utility functions not for export
function timed_wait!(conn::TCPSocket, to::Real, b::Bool)
  n = bytesavailable(conn)
  m = copy(n)
  ts = time()
  te = copy(ts)
  while (te-ts) < to && m == n
    sleep(1.0)
    b = eof(conn)
    te = time()
    n = bytesavailable(conn)
  end
  if m == bytesavailable(conn)
    close(conn)
    error("Connection timed out.")
  end
  return nothing
end

# This was deprecated in Julia 0.6; hard-copied here, still works
function sync_add(r::Task)
    spawns = get(task_local_storage(), :SPAWNS, ())
    if spawns != ()
        push!(spawns[1], r)
        tls_r = get_task_tls(r)
        tls_r[:SUPPRESS_EXCEPTION_PRINTING] = true
    end
    r
end

function get_sl_ver(vline::String)
  # Versioning will break if SeedLink switches to VV.PPP.NNN format
  ver = 0.0
  vinfo = split(vline)
  for i in vinfo
      if startswith(i, 'v')
          try
              ver = Meta.parse(i[2:end])
          catch
              continue
          end
      end
  end
  return ver
end

function check_sta_exists(sta::Array{String,1}, xstr::String)

  xstreams = get_elements_by_tagname(root(parse_string(xstr)), "station")
  xid = join([join([attribute(xstreams[i], "network"),
                    attribute(xstreams[i], "name")],'.') for i=1:length(xstreams)], ' ')
  N = length(sta)
  x = falses(N)
  for i = 1:N
    id = split(sta[i], '.', keepempty=true)
    sid = join(id[1:2],'.')
    if occursin(sid, xid)
      x[i] = true
    end
  end
  return x
end

function check_stream_exists(S::Array{String,1}, xstr::String;
                             gap::Real=KW.SL.gap,
                             to::Int=KW.to
                             )

  a = ["seedname","location","type"]
  N = length(S)
  x = falses(N)

  xstreams = get_elements_by_tagname(root(parse_string(xstr)), "station")
  xid = String[join([attribute(xstreams[i], "network"),
                     attribute(xstreams[i], "name")],'.') for i=1:length(xstreams)]
  for i = 1:N
    # Assumes the combination of network name and station name is unique
    id = split(S[i], '.', keepempty=true)
    sid = join(id[1:2], '.')
    K = findid(sid, xid)
    if K > 0
      t = Inf

      # Syntax requires that contains(string, "") returns true for any string
      p = ["","",""]
      for j = 3:length(id)
        p[j-2] = replace(id[j], "?" => "")
      end
      R = get_elements_by_tagname(xstreams[K], "stream")
      if !isempty(R)
        for j = 1:length(R)
          if prod([occursin(p[i], attribute(R[j], a[i])) for i=1:length(p)]) == true
            te = replace(attribute(R[j], "end_time"), " " => "T")
            t = min(t, time()-d2u(Dates.DateTime(te)))
          end
        end
      end

      # Treat station as "present" if there's a match
      if minimum(t) < gap
        x[i] = true
      end
    end
  end
  return x
end
# ========================================================================

"""
    info_xml = sl_info(level=LEVEL::String; u=URL::String, port=PORT::Integer)

Retrieve XML output of SeedLink command "INFO `LEVEL`" from server `URL:PORT`.
Returns formatted XML. `LEVEL` must be one of "ID", "CAPABILITIES",
"STATIONS", "STREAMS", "GAPS", "CONNECTIONS", "ALL".

"""
function sl_info(level::String;                                 # verbosity
                to::Int64     = KW.to,                          # Timeout (s)
                 u::String    = "rtserve.iris.washington.edu",  # url
                 port::Int64  = KW.SL.port                      # port
                 )
  conn = connect(TCPSocket(), u, port)
  write(conn, string("INFO ", level, "\r"))
  b = false
  timed_wait!(conn, to, b)
  buf = BUF.buf
  if (bytesavailable(conn) == 0) || (isopen(conn) == false)
    @warn(("Connection sent no data after ", to, " seconds; closing."))
    (isopen(conn)) && close(conn)
    x = Array{UInt8, 1}(undef, 0)
  else
    eflg = false
    x = Array{UInt8, 1}(undef, 1048576)
    i = 0x0000000000000001
    c = 0x01
    k = 0x0000
    while true
      N = conn.buffer.size
      checkbuf!(buf, N)

      # buffer to :buf from conn
      unsafe_read(conn.buffer, pointer(buf), N)
      conn.buffer.ptr = 1
      conn.buffer.size = 0

      # copy to x after filtering out 0x00 and first 64 vals of each 520-byte packet
      j = 0x0000000000000000
      while j < N
        j += 0x0000000000000001
        k += 0x0001
        if i + 0x00000000000001c8 > length(x)
          resize!(x, length(x) + 524288)
        end
        if k >= 0x0040
          c = getindex(buf, j)
          if c != 0x00
            setindex!(x, c, i)
            i += 0x0000000000000001
          end
        end
        if k == 0x0208
          if buf[j] == 0x00
            eflg = true
            break
          end
          k = 0x0000
        end
      end
      eflg && break
      timed_wait!(conn, to, b)
    end
    close(conn)
    deleteat!(x, i:length(x))
  end
  resize!(buf, 65535)
  return String(x)
end

"""
    has_sta(sta[, u=url, port=n])

Check that streams exist at `url` for stations `sta`, formatted
NET.STA. Use "?" to match any single character. Returns `true` for
stations that exist. `sta` can also be the name of a valid config
file or a 1d string array.

Returns a BitArray with one value per entry in `sta.`

SeedLink keywords: gap, port
"""
function has_sta(C::String;
  u::String     = "rtserve.iris.washington.edu",
  port::Int64   = KW.SL.port
  )

  sta, pat = parse_sl(parse_chstr(C))

  # This exists to convert SeedLink syntax (SSS NN) to FDSN (NN.SSS)
  for i = 1:length(sta)
    s = split(sta[i], ' ')
    sta[i] = join([s[2], s[1]], '.')
  end
  # This exists to convert SeedLink syntax (SSS NN) to FDSN (NN.SSS)

  return check_sta_exists(sta, sl_info("STATIONS", u=u, port=port))
end

has_sta(sta::Array{String,1};
  u::String     = "rtserve.iris.washington.edu",
  port::Int64   = KW.SL.port
  ) = check_sta_exists(sta, sl_info("STATIONS", u=u, port=port))

has_sta(sta::Array{String,2};
  u::String   = "rtserve.iris.washington.edu",
  port::Int64 = KW.SL.port
  ) = check_sta_exists([join(sta[i,:], '.') for i=1:size(sta,1)],
                       sl_info("STATIONS", u=u, port=port))

"""
    has_stream(cha[, u=url, port=N, gap=G)

Check that streams with recent data exist at url `u` for channel spec
`cha`, formatted NET.STA.LOC.CHA.DFLAG, e.g. "UW.TDH..EHZ.D,
CC.HOOD..BH?.E". Use "?" to match any single character. Returns `true`
for streams with recent data.

`cha` can also be the name of a valid config file.

    has_stream(sta::Array{String,1}, sel::Array{String,1}, u::String, port=N::Int, gap=G::Real)

If two arrays are passed to has_stream, the first should be
formatted as SeedLink STATION patterns (formated "SSSSS NN", e.g.
["TDH UW", "VALT CC"]); the second be an array of SeedLink selector
patterns (formatted LLCCC.D, e.g. ["??EHZ.D", "??BH?.?"]).

SeedLink keywords: gap, port
"""
function has_stream(sta::Array{String,1}, pat::Array{String,1};
  u::String   = "rtserve.iris.washington.edu",
  to::Int64   = KW.to                               ,  # Timeout (s)
  port::Int64 = KW.SL.port,
  gap::Real   = KW.SL.gap,
  d::Char     = ' '
  )

  L = length(sta)
  cha = Array{String,1}(undef, L)
  for i = 1:L
    s = split(sta[i], d)
    c = split(pat[i], '.')
    cha[i] = join([s[1], s[2], c[1][1:2], c[1][3:5], c[2]], '.')
  end
  return check_stream_exists(cha, sl_info("STREAMS", u=u, port=port), gap=gap, to=to)
end

has_stream(sta::String;
u::String   = "rtserve.iris.washington.edu",
to::Int64   = KW.to                               ,  # Timeout (s)
port::Int64 = KW.SL.port,
gap::Real   = KW.SL.gap,
d::Char     = ','
) = check_stream_exists(String.(split(sta, d)), sl_info("STREAMS", u=u, port=port), gap=gap, to=to)

has_stream(sta::Array{String,1};
  u::String   = "rtserve.iris.washington.edu",
  to::Int64   = KW.to                               ,  # Timeout (s)
  port::Int64 = KW.SL.port,
  gap::Real   = KW.SL.gap
  ) = check_stream_exists(sta, sl_info("STREAMS", u=u, port=port), gap=gap, to=to)

has_stream(sta::Array{String,2};
  u::String   = "rtserve.iris.washington.edu",
  to::Int64   = KW.to                               ,  # Timeout (s)
  port::Int64 = KW.SL.port,
  gap::Real   = KW.SL.gap
  ) = check_stream_exists([join(sta[i,:], '.') for i=1:size(sta,1)],
                          sl_info("STREAMS", u=u, port=port), gap=gap, to=to)

@doc """
    seedlink!(S, sta)

Begin acquiring seedlink data to SeisData structure `S`. New channels
are added to `S` automatically based on `sta`. Connections are added
to S.c.

    S = seedlink(sta)

Create a new SeisData structure `S` to acquire seedlink data. Connection will be in S.c[1].

### INPUTS
* `S`: SeisData object
* `sta`: Array{String, 1} formatted NET.STA.LOC.CHA.DFLAG, e.g. ["UW.TDH..EHZ.D", "CC.HOOD..BH?.E"].
Use "?" to match any single character; leave LOC and CHA fields blank to select all. Don't use "*" for wildcards,
it isn't valid seedlink syntax.

When finished, close connection manually with `close(S.c[n])` where n is connection #.

* Standard keywords: fmt, opts, q, si, to, v, w, y
* SL keywords: gap, kai, mode, port, refresh, safety, x_on_err
""" seedlink!
function seedlink!(S::SeisData, sta::Array{String,1}, patts::Array{String,1};
                    gap::Real=KW.SL.gap,
                    kai::Real=KW.SL.kai,
                    mode::String=KW.SL.mode,
                    port::Int64=KW.SL.port,
                    refresh::Real=KW.SL.refresh,
                    s::TimeSpec=0,
                    t::TimeSpec=300,
                    u::String="rtserve.iris.washington.edu",
                    v::Int64=KW.v,
                    w::Bool=KW.w,
                    x_on_err::Bool=KW.SL.x_on_err
                  )

  # ==========================================================================
  # init, warnings, sanity checks
  Ns = size(sta,1)
  setfield!(BUF, :swap, false)

  # Refresh interval
  refresh = maximum([refresh, eps()])
  refresh < 10 && @warn(string("refresh = ", refresh, " < 10 s; Julia may freeze if no packets arrive between consecutive read attempts."))

  # keepalive interval
  if kai < 240
    @warn("KeepAlive interval < 240s violates IRIS netiquette guidelines. YOU are responsible if you get banned for hammering.")
  end

  # Source for logging
  src = join([u,port],':')

  # ==========================================================================
  # connection and server info retrieval
  push!(S.c, connect(TCPSocket(),u,port))
  q = length(S.c)

  # version, server info
  write(S.c[q],"HELLO\r")
  vline = readline(S.c[q])
  sline = readline(S.c[q])
  ver = get_sl_ver(vline)

  (v > 1) && println("Version = ", ver)
  (v > 1) && println("Server = ", strip(sline,['\r','\n']))
  # ==========================================================================

  # ==========================================================================
  # handshaking

  # create mode string and filename for -w
  (d0,d1) = parsetimewin(s,t)
  s = join(split(d0,r"[\-T\:\.]")[1:6],',')
  t = join(split(d1,r"[\-T\:\.]")[1:6],',')
  if mode in ["TIME", "FETCH"]
    if mode == "TIME"
      m_str = string("TIME ", s, " ", t, "\r")
    else
      m_str = string("FETCH\r")
    end
  else
    m_str = string("DATA\r")
  end

  if w
    fname = hashfname([join(sta,','), join(patts,','), s, t, m_str], "mseed")
    (v > 0) && println(string("Raw packets will be written to file ", fname, " in dir ", realpath(pwd())))
    fid = open(fname, "w")
  end

  # pass strings to server; check responses carefully
  for i = 1:Ns
    # pattern selector
    sel_str = string("SELECT ", patts[i], "\r")
    (v > 1) && println("Sending: ", sel_str)
    write(S.c[q], sel_str)
    sel_resp = readline(S.c[q])
    (v > 1) && println("Response: ", sel_resp)
    if occursin("ERROR", sel_resp)
      if x_on_err
        write(S.c[q], "BYE\r")
        close(S.c[q])
        deleteat!(S.c, q)
        error(string("Error in select string ", patts[i], "; connection closed and deleted, exit with error."))
      else
        @warn(string("Error in select string ", patts[i], "; previous selector, ", i==1 ? "*" : patts[i-1], " used)."))
      end
    end

    # station selector
    sta_str = string("STATION ", sta[i], "\r")
    (v > 1) && println("Sending: ", sta_str)
    write(S.c[q], sta_str)
    sta_resp = readline(S.c[q])
    (v > 1) && println("Response: ", sta_resp)
    if occursin("ERROR", sta_resp)
      if x_on_err
        write(S.c[q], "BYE\r")
        close(S.c[q])
        deleteat!(S.c, q)
        error(string("Error in station string ", sta[i], "; connection closed and deleted, exit with error."))
      else
        @warn(string("Error in station string ", sta[i], " (station excluded)."))
      end
    end

    # mode
    (v > 1) && println("Sending: ", m_str)
    write(S.c[q], m_str)
    m_resp = readline(S.c[q])
    (v > 1) && println("Response: ", m_resp)
  end
  write(S.c[q],"END\r")
  # ==========================================================================

  # ==========================================================================
  # data transfer
  k = @task begin
    j = 0
    while true
      if !isopen(S.c[q])
          (v > 0) && @info(string(timestamp(), ": SeedLink connection closed."))
        w && close(fid)
        break
      else

        #= use of rand() makes it almost impossible for multiple SeedLink
        connections to result in one sleeping indefinitely. =#
        τ = ceil(Int, refresh*(0.8 + 0.2*rand()))
        sleep(τ)
        eof(S.c[q])
        N = div(bytesavailable(S.c[q]), 520)
        if N > 0
          #=
            This syntax consumes 2x memory but is "safe" in that a partially
            downloaded packet will be left in S.c[q], rather than causing the
            SEED parser to throw errors and/or hang the REPL. (From experience,
            it is NOT guaranteed that full packets are transmitted)
          =#
          buf = read(S.c[q], 520*N)
          io = IOBuffer(buf)
          if w
            write(fid, copy(io))
          end
          (v > 1) && printstyled(stdout, now(), ": Processing packets ", color=:green)
          while !eof(io)
            pkt_id = String(read(io, 8))
            (v > 1) && @printf(stdout, "%s, ", pkt_id)
            parserec!(S, BUF, io, v, 65535, 65535)
          end
          (v > 1) && @printf(stdout, "\b\b...done current packet dump.\n")
          seed_cleanup!(S, BUF)
        end

        # SeedLink (non-standard) keep-alive gets sent every kai seconds
        j += τ
        if j ≥ kai
          # Secondary "isopen" loop avoids possible error from race condition
          # maybe a Julia bug? First encountered 2017-01-03
          if isopen(S.c[q])
            j -= kai
            write(S.c[q],"INFO ID\r")
          end
        end

      end
    end
  end
  sync_add(k)
  Base.enq_work(k)
  # ========================================================================

  return S
end
function seedlink!(S::SeisData, C::Union{String,Array{String,1},Array{String,2}};
                    gap::Real=KW.SL.gap,
                    kai::Real=KW.SL.kai,
                    mode::String=KW.SL.mode,
                    port::Int64=KW.SL.port,
                    refresh::Real=KW.SL.refresh,
                    s::TimeSpec=0,
                    t::TimeSpec=300,
                    v::Int64=KW.v,
                    u::String="rtserve.iris.washington.edu",
                    w::Bool=KW.w,
                    x_on_err::Bool=KW.SL.x_on_err)

  if isa(C, String)
    sta, pat = parse_sl(parse_chstr(C))
  elseif ndims(C) == 1
    sta, pat = parse_sl(parse_charr(C))
  else
    sta, pat = parse_sl(C)
  end
  seedlink!(S, sta, pat, u=u, port=port, mode=mode, refresh=refresh, kai=kai, s=s, t=t, x_on_err=x_on_err, v=v, w=w)
  return S
end

@doc (@doc seedlink!)
function seedlink(sta::Array{String,1}, pat::Array{String,1};
  gap::Real=KW.SL.gap,
  kai::Real=KW.SL.kai,
  mode::String=KW.SL.mode,
  port::Int64=KW.SL.port,
  refresh::Real=KW.SL.refresh,
  s::TimeSpec=0,
  t::TimeSpec=300,
  v::Int64=KW.v,
  u::String="rtserve.iris.washington.edu",
  w::Bool=KW.w,
  x_on_err::Bool=KW.SL.x_on_err)

  S = SeisData()
  seedlink!(S, sta, pat, u=u, port=port, mode=mode, refresh=refresh, kai=kai, s=s, t=t, x_on_err=x_on_err, v=v, w=w)
  return S
end

function seedlink(C::Union{String,Array{String,1},Array{String,2}};
  gap::Real=KW.SL.gap,
  kai::Real=KW.SL.kai,
  mode::String=KW.SL.mode,
  port::Int64=KW.SL.port,
  refresh::Real=KW.SL.refresh,
  s::TimeSpec=0,
  t::TimeSpec=300,
  v::Int64=KW.v,
  u::String="rtserve.iris.washington.edu",
  w::Bool=KW.w,
  x_on_err::Bool=KW.SL.x_on_err)

  S = SeisData()
  if isa(C, String)
    sta,pat = parse_sl(parse_chstr(C))
  elseif ndims(C) == 1
    sta,pat = parse_sl(parse_charr(C))
  else
    sta, pat = parse_sl(C)
  end
  seedlink!(S, sta, pat, u=u, port=port, mode=mode, refresh=refresh, kai=kai, s=s, t=t, x_on_err=x_on_err, v=v, w=w)
  return S
end


text:
using Random: seed!
seed!(4231)
using Flux
using Zygote
using BSON: @save, @load
ENV["PLOTS_DEFAULT_BACKEND"] = "PyPlot"
using Plots
pyplot()
using StatsPlots

include("NN.jl")
using .NN
include("Solver.jl")
using .Solver

function RSSapproximate(model; ε::Real, n::Real, hoursperlayer::Real = 0, experimentname::String)
    @assert 0 < ε
    @assert n > 0

    # Make sure the original model is not altered and the copy is in CPU
    model = cpu(deepcopy(model))
    filename = "saves/$experimentname-RSS-n$n-h$hoursperlayer.bson"
    # Check for cached model
    if isfile(filename)
        @info "Found a saved approximation in file '$filename'. Loading from it..."
        @load filename parameters
        try
            Flux.loadparams!(model, parameters)
        catch e
            @warn "Current model is not compatible with the weights in file '$filename'"
            @info "$model\nReturning"
            return nothing
        end
    else
        nparams = countparams(model)
        @info "Subset Summing $nparams parameters\n  n = $n\n  ε = $(ε/nparams)"
        for l ∈ model
            if l isa Dense || l isa Conv
                @info "Subset Summing $l layer with $(sum(length, params(l))) parameters"
                for p ∈ params(l)
                    @time p .= randomsubsetsum.(p, ε = ε/nparams, n = n, verbose = false, minsecs = hoursperlayer*60*60/nparams)
                end
            elseif l isa MaxPool || l isa typeof(flatten)
                # Nothing to be done
            else
                @warn "TODO: Support layer of type `$(typeof(l))`"
                @info "Ignoring"
            end
        end
        @save filename parameters=params(cpu(model))
    end
    return model
end

function buildlenet5(; inputsize::Tuple{Int, Int, Int}, nclasses::Int, init::Function = Flux.kaiming_uniform)
    cnn = Chain(
        Conv((5, 5), inputsize[end]=>6, pad = SamePad(), relu, init = init),
        MaxPool((2,2), stride=2),

        Conv((5, 5), 6=>16, pad = SamePad(), relu, init = init),
        MaxPool((2,2), stride=2),
    )

    cnnoutsize = Flux.outputsize(cnn, inputsize, padbatch=true)

    return Chain(
        cnn...,
        flatten,
        Dense(prod(cnnoutsize), 120, relu, init = init),
        Dense(120, 84, relu, init = init),
        Dense(84, nclasses, init = init),
    )
end

countparams(model) = sum(length, params(model))

function run()
    for dataparams ∈ [
        (datasetname = "MNIST", dataset = getprocessedMNIST()),
        (datasetname = "FashionMNIST", dataset = getprocessedFashionMNIST()),
    ]
        experimentname = "lenet5-$(dataparams.datasetname)"
        model = trainmodel(
            name = experimentname,
            model = buildlenet5(inputsize=(28,28,1), nclasses=10),
            dataset = dataparams.dataset,
            loss = (ŷ, y)->Flux.logitcrossentropy(ŷ, y),
            opt = ADAM(0.001),
            epochs = 50,
            batchsize = 64,
        )
        println("l1 = ", map(W->round(sum(abs, W), digits = 2), params(model)))
        println("l∞ = ", map(W->round(maximum(abs, W), digits = 2), params(model)))
        model = cpu(model[1:4])  # Just the convolutional layers
        ns = 5:5:30
        # Large ε is effectively ignored
        approxmodels = [n => RSSapproximate(model, ε = 10*countparams(model), n = n, hoursperlayer = Inf, experimentname = experimentname) for n ∈ ns]

        maxweightdeltas = [
            map(params(model), params(modelRSS)) do p, p′
                maximum(abs, (p .- p′) ./ maximum(abs, p))
            end
            for modelRSS ∈ last.(approxmodels)
        ]
        plot(
            first.(approxmodels),
            maximum.(maxweightdeltas),
            legends = false,
            marker = :circle,
            xlabel = "sample size",
            ylabel = "maximum relative weight error",
            yticks = 10.0 .^ (0:-1:-4),
            yaxis = :log,
            size = (600, 350),
        )
        savefig("plots/maxweightdelta-$(dataparams.datasetname).pdf")

        @info "Computing outputs"
        # Use GPU if available
        model = gpu(model)
        approxmodels = @. first(approxmodels) => gpu(last(approxmodels))
        inputs_train = gpu(dataparams.dataset.train.inputs)
        inputs_test = gpu(dataparams.dataset.test.inputs)
        referenceoutputs_train = model(inputs_train)
        referenceoutputs_test = model(inputs_test)
        outputerrors_train = [cpu(abs.(referenceoutputs_train .- modelRSS(inputs_train)) ./ maximum(abs, referenceoutputs_train)) for modelRSS ∈ last.(approxmodels)]
        outputerrors_test = [cpu(abs.(referenceoutputs_test .- modelRSS(inputs_test)) ./ maximum(abs, referenceoutputs_test)) for modelRSS ∈ last.(approxmodels)]
        plot(
            first.(approxmodels),
            maximum.(outputerrors_train),
            label = "train",
            marker = :circle,
            line = :solid,
        )
        plot!(
            first.(approxmodels),
            maximum.(outputerrors_test),
            label = "test",
            marker = :circle,
            line = :dash,
        )
        plot!(size = (600, 350))
        xlabel!("sample size")
        ylabel!("maximum relative output error")
        yaxis!(:log)
        yticks!(10.0 .^ (0:-1:-5))
        savefig("plots/maxoutputerrors-$(dataparams.datasetname).pdf")
    end
    return nothing
end

run()

max_stars_repo_name: ArthurWalraven/cnnslth

max_issues_repo_name: ArthurWalraven/cnnslth

max_forks_repo_name: ArthurWalraven/cnnslth


text:
"""
    get_title(pdf::String) -> String

Extract the title from the PDF.

# Arguments
- `pdf::String`: $(PDF_ARGUMENT)

# Returns
- `String`: the title

# Throws
- Exceptions from: [`get_author_title`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf = joinpath(
    pathof(PDFHighlights) |> dirname |> dirname,
    "test",
    "pdf",
    "TestPDF.pdf",
)

get_title(path_to_pdf) == "A dummy PDF for tests"

# output

true
```
"""
function get_title(pdf::String)::String
    _, title = get_author_title(pdf)
    return title
end


text:
#Libraries used
using YAML
using Parameters
using Interpolations
using LinearAlgebra
using ProgressBars
using Plots
using NPZ
using DifferentialEquations

"""
    read_inp(debug=0)

Read the input from the input file, creates all the global numbering variables used in later stages

Inputs:
	inp_file- Name of the input YAML file with .yaml extension
    debug   - Integer 0,1 to switch on or off debugging

Outputs:
	sref      - reference surface area
	bref      - reference span
	cref      - reference chord length
	aseq      - sequence of alpha to investigate of form [starting,ending,number of steps] of variable types [float,float,integer]
	vinf      - value of vinf
	ρ         - density
	alpha     - angle of attack
	nwing     - number of wings
	ntot_patch- total number of patches
	npatch_perwing
	          - patches per wing (used for global patches numbering)
	ntot_lat  - total number of lattices (used for global lattice numbering)
	nspan_perpatch
	          - number of spanwise strips per patch (used to loop over spanwise strips)
	nlat_perpatch
			  - number of total lattices per patch
	xyz_qc_patch
			  - quarter chord locations for each patch of form [(x,y,z),(start,end),global patch number]
	chord_patch
		      - chord lengths for the pacthes of form [(start,end),global patch number]
	twist_patch
			  - twists for the pacthces of the form [(start,end),gloabl patch number]
	 α_zl_patch
	 		  - zero lift alpha values of the form [(start,end),gloabl patch number]

"""
function read_inp(inp_file,debug=0)
        inp_data = YAML.load_file(inp_file)
         @info "Reading Input File"
        name = inp_data["name"]
        sref = inp_data["sref"]
        bref = inp_data["bref"]
        cref = inp_data["cref"]
        nwing= inp_data["nwing"]

        #create arrays required to store
        npatch_perwing = zeros(Int64,nwing)
         for i in 1:nwing
                str = "wing"*string(i)
                npatch_perwing[i] = inp_data[str]["npatch"]
        end

        ntot_patch = sum(npatch_perwing)
        #get some arrays to store root and tip info
        xyz_qc_patch = zeros(3,2,ntot_patch)
        chord_patch  = zeros(2,ntot_patch)
        twist_patch  = zeros(2,ntot_patch)
        α_zl_patch   = zeros(2,ntot_patch)

        nlat_perpatch = zeros(Int64,ntot_patch)
		nspan_perpatch = zeros(Int64,ntot_patch)
        ntot_lat   = 0
        counter    = 0
        for i in 1:nwing
               str = "wing"*string(i)
                for j in 1:npatch_perwing[i]
                        counter = counter + 1
						nspan_perpatch[counter]   = inp_data[str]["patch"*string(j)]["nspan"]
                        nlat_perpatch[counter]    = inp_data[str]["patch"*string(j)]["nlat"]*nspan_perpatch[counter]

                        # at root
                        xyz_qc_patch[:,1,counter] = inp_data[str]["patch"*string(j)]["root"]["xyz_qc"]
                        chord_patch[1,counter]    = inp_data[str]["patch"*string(j)]["root"]["chord"]
                        twist_patch[1,counter]    = inp_data[str]["patch"*string(j)]["root"]["twist"]*pi/180
                        α_zl_patch[1,counter]     = inp_data[str]["patch"*string(j)]["root"]["alpha_zl"]*pi/180

                        # at tip
                        xyz_qc_patch[:,2,counter] = inp_data[str]["patch"*string(j)]["tip"]["xyz_qc"]
                        chord_patch[2,counter]    = inp_data[str]["patch"*string(j)]["tip"]["chord"]
                        twist_patch[2,counter]    = inp_data[str]["patch"*string(j)]["tip"]["twist"]*pi/180
                        α_zl_patch[2,counter]     = inp_data[str]["patch"*string(j)]["tip"]["alpha_zl"]*pi/180

                        ntot_lat = ntot_lat + nlat_perpatch[counter]
                end
        end

		aseq = inp_data["aseq"]
		aseq = [aseq[1]*pi/180,aseq[2]*pi/180,aseq[3]]
		vinf = inp_data["vinf"]
		ρ    = inp_data["rho"]
		alpha= inp_data["alpha"]*pi/180

        if (debug==1)
                #checks

                #print some stuff if needed
                display("ntot_lat:"*string(ntot_lat))
                display(xyz_qc_patch)

        end

        return (sref,bref,cref,aseq,vinf,ρ,alpha,nwing,ntot_patch,npatch_perwing,ntot_lat,nspan_perpatch,nlat_perpatch,xyz_qc_patch,chord_patch,twist_patch,α_zl_patch)
end

"""
    geom_calc(nwing,npatch,npatch_perwing,ntot_lat,nlat_perpatch,xyz_qc_patch,chord_patch,twist_patch,α_zl_patch,Λ,debug)

Calculate the geometric propoerties for all the lattice elements

	Inputs:
		nwing      - number of wings
		ntot_patch - total number of patches
		npatch_perwing
				   - number of patches per wing, form : [Integer value]
		ntot_lat   - total numer of lattices
		nspan_perpatch
				   - number of spanwise strips per patch
		nlat_perpatch
				   - number of lattices per patch (uses global patch numbering)
		xyz_qc_patch
				   - quarter chord locations for each patch of form [(x,y,z),(start,end),global patch number]
		 chord_patch
				   - chord lengths for the pacthes of form [(start,end),global patch number]
		 twist_patch
				   - twists for the pacthces of the form [(start,end),gloabl patch number]
		α_zl_patch
	 			   - zero lift alpha values of the form [(start,end),gloabl patch numb
		debug      - debug flag to enable some debugging options

	Output:
		nspn       - number of span wise locations to plt local distributions like Cl
						: If more than one wing, each wings span wise locations are stored in that order
						  Example - Wing + Tail
						  nspn = span locations on wing + span locations on tail
		spn_map    - mapping for all the lattices to span locations
						: form - Array{Integers},size "ntot_lat" pointing to position in "spn_loc" array
		spn_loc    - locations along span (y-coordinates)
						: form - Array{Float,1}, size "nspn"
		sbar       - locations of starting point of horse-shoe lattice for all lattices
						: form - Array{Float,2}, size {3,ntot_lat}
		ebar       - locations of ending point of horse-shoe lattice for all lattices
						: form - Array{Float,2}, size {3,ntot_lat}
		mbar       - locations of middle point of horse-shoe lattice for all lattices
						: form - Array{Float,2}, size {3,ntot_lat}
		nbar       - perpendicular unit vectors of horse-shoe lattice for all lattices
						: form - Array{Float,2}, size {3,ntot_lat}
		cbar       - locations of control points or quarter chord points of horse-shoe lattice for all lattices
						: form - Array{Float,2}, size {3,ntot_lat}
		dbar       - locations of middle points of horse-shoe lattice for all lattices at Treffts Plane
						: form - Array{Float,2}, size {3,ntot_lat}
		tbar       - unit vector for pointing the tail of horse shoe vortices
						: form - Array{Float,2}, size {3,ntot_lat}
		chord      - chord length values at all lattices
						: form - Array{Float,1}, size {ntot_lat}
		twist      - twist values at all lattices
						: form - Array{Float,1}, size {ntot_lat}
		α_zl       - alpha zero lift values at all lattices
						: form - Array{Float,1}, size {ntot_lat}
		θ          - dihedral angle at all lattices
						: form - Array{Float,1}, size {ntot_lat}
		Λ          - sweep angle at all lattices
						: form - Array{Float,1}, size {ntot_lat}
		ds         - length of filament for all lattices
						: form - Array{Float,1}, size {ntot_lat}

"""
function geom_calc(nwing,ntot_patch,npatch_perwing,ntot_lat,nspan_perpatch,nlat_perpatch,xyz_qc_patch,chord_patch,twist_patch,α_zl_patch,debug=0)

        #create arrays for the lattices
        sbar = zeros(3,ntot_lat)        # starting of bv
        ebar = zeros(3,ntot_lat)        # ending of bv
        mbar = zeros(3,ntot_lat)        # mid of bv
        cbar = zeros(3,ntot_lat)        # ctrl poi of bv
        dbar = zeros(3,ntot_lat)        # trefftz plane mid poi
        nbar = zeros(3,ntot_lat)		# normal for patch/wing/lattice
        tbar = zeros(3,ntot_lat)        # tail direction
		univ = zeros(3,ntot_lat)        #used for calc nbar

        chord = zeros(ntot_lat)         # chord length
        twist = zeros(ntot_lat)			# twist
        α_zl  = zeros(ntot_lat)			# α_zero-lift
        θ     = zeros(ntot_lat)			# dihedral angle
		Λ     = zeros(ntot_lat)         # sweep angle
		ds    = zeros(ntot_lat)         # lattice length

		spn_map 	 = zeros(Int64,ntot_lat) #Spanwise mapping
		lat2wing_map = zeros(Int64,ntot_lat)  #Mapping from lattice to wing

		nlatspn_perwing = zeros(Int64,nwing) #Span per wing
		nlat_perwing    = zeros(Int64,nwing) #Number of lattices in a wing

		# Loop over to find total span wise locations (Could be eliminated/improved)
		nspn  = 0
		ipatch=0
		for i in 1:nwing
			nlatw = 0
			for j in 1:npatch_perwing[i]
				ipatch = ipatch + 1
				nlat_perspan = Int64(nlat_perpatch[ipatch]/nspan_perpatch[ipatch])
				for k in 1:nlat_perspan
					nspn = nspn + 1
				end
				nlatw = nlatw + nlat_perpatch[ipatch]
			end
			nlat_perwing[i] = nlatw
		end

		spn_loc = zeros(nspn) #Spanwise locations

		#Initialise counters
        ilat  = 0
        ipatch= 0
		nspn  = 0

        for i in  tqdm(1:nwing,leave=false)
                for j in 1:npatch_perwing[i]
                        ipatch = ipatch + 1 #for global patch number

						# Interpolate chord, twist and α_zl values
                        itp_chord = LinearInterpolation(xyz_qc_patch[2,:,ipatch],chord_patch[:,ipatch])
                        itp_twist = LinearInterpolation(xyz_qc_patch[2,:,ipatch],twist_patch[:,ipatch])
                        itp_α_zl  = LinearInterpolation(xyz_qc_patch[2,:,ipatch],α_zl_patch[:,ipatch])

						# Get the starting and ending locations along the chord for each patch
						nlat_perspan = Int64(nlat_perpatch[ipatch]/nspan_perpatch[ipatch])
						a = copy(xyz_qc_patch[:,1,ipatch])
                        b = copy(xyz_qc_patch[:,2,ipatch])

						#
						# Computing starting root and tip locations for each span
						#
						a_span = zeros(3,nspan_perpatch[ipatch])
						b_span = zeros(3,nspan_perpatch[ipatch])

						#pvec is perpendicular to (b-a)
						pvec   = (b[1:2] - a[1:2])/norm(b[1:2]-a[1:2])
						pvec   = [pvec[2],-pvec[1]]

						#loop over strips along the chord wise
						ns     = nspan_perpatch[ipatch]
						for ispan in 1:ns
							a_span[:,ispan] = a + 0.25*chord_patch[1,ipatch]*(4*ispan-3-ns)/ns * [pvec[1],pvec[2],0.0]
							b_span[:,ispan] = b + 0.25*chord_patch[2,ipatch]*(4*ispan-3-ns)/ns * [pvec[1],pvec[2],0.0]
						end

						#Sanity check to see the pvec vector is in the right direction
						@info pvec
						if (dot(pvec,[1.0,0.0])<=0.0)
							println("Error in geom_calc span wise variables")
							println("Perp Vector Calc wrong")
							return
						end

						#Loop over the chord wise strips and lattices in each strip
						for ispan in 1:nspan_perpatch[ipatch]
						for k in 1:nlat_perspan

                                ilat = ilat + 1

                                a = copy(a_span[:,ispan])
                                b = copy(b_span[:,ispan])

								#Find start,end,mid,normal vec, dihedral, sweep, α_zl, ...
								sbar[:,ilat] =  a + (b - a)/nlat_perspan .*(k-1)
                                ebar[:,ilat] =  a + (b - a)/nlat_perspan .*(k)
                                mbar[:,ilat] = 0.5*(sbar[:,ilat]+ebar[:,ilat])
                                univ[:,ilat] = cross((ebar[:,ilat] - sbar[:,ilat]),[1.0,0.0,0.0])/norm(cross((ebar[:,ilat] - sbar[:,ilat]),[-1.0,0.0,0.0]))

                                θ[ilat]      = acos(nbar[3,ilat])
								Λ[ilat]      = acos(dot(b[1:2]-a[1:2],[0.0,1.0])/norm(b[1:2]-a[1:2]))
                                chord[ilat]  = itp_chord(mbar[2,ilat])
                                twist[ilat]  = itp_twist(mbar[2,ilat])
                                α_zl[ilat]   = itp_α_zl(mbar[2,ilat])
								ds[ilat]	 = norm(ebar[2:3,ilat]-sbar[2:3,ilat])#*cos(Λ[ilat])

                                cbar[:,ilat] = mbar[:,ilat] - [chord[ilat]/(2*nspan_perpatch[ipatch]),0.0,0.0]
                                dbar[:,ilat] = [-10^8,mbar[2,ilat],mbar[3,ilat]]
								tbar[:,ilat] = [1.0,0.0,0.0]   #Trailing edge vector
								rot          = [cos(twist[ilat]) 0.0 sin(twist[ilat]);0.0 1.0 0.0;-sin(twist[ilat]) 0.0 cos(twist[ilat])]
								nbar[:,ilat] = rot*univ[:,ilat]

								#To get the spn_loc locations only when using the first strip
								if (ispan==1)
									nspn = nspn + 1
									spn_loc[nspn] = mbar[2,ilat]
								end

								#Compute ispn_bef to locate where in spn_map this patch is starting
								ispn_bef = 0
								for ip in 1:ipatch-1
									ispn_bef = ispn_bef + Int64(nlat_perpatch[ip]/nspan_perpatch[ip])
								end

								#Span map
								spn_map[ilat] = ispn_bef + k
								#Lat2Wing map
								lat2wing_map[ilat] = i

						end
						end
					#finding numer of lattices in spn wise per each wing
					nlatspn_perwing[i] = nlatspn_perwing[i] + nlat_perspan
				end
        end

		#Debug stuff
        if debug==1
                display("sbar")
                display(sbar)
                display("xyz_qc_patch")
                display(xyz_qc_patch)
        end

        return nspn,spn_map,spn_loc,lat2wing_map,sbar,ebar,mbar,nbar,cbar,dbar,tbar,chord,twist,α_zl,θ,Λ,ds
end
"""
    calc_vind_sinf(rbar,sbar,tbar)

    Function to influence for semi-infinte vortex element given starting,ending and tail direction

    Input :
		rbar - vector for pointing the location of evaluation
				: form - Array{Float,1}, size {3}
		sbar - vector for pointing the start of filament
				: form - Array{Float,1}, size {3}
		tbar - vector for pointing the trailing of vortex
				: form - Array{Float,1}, size {3}

    Output:
		q    - vector of influence at the given location
				: form - Array{Float,1}, size {3}

"""
function calc_q_sinf(rbar,sbar,tbar)

	(x1,y1,z1) = sbar
	(x,y,z)    = rbar

	avec = [x1-x,y1-y,z1-z]

	amag = norm(avec)

	cprd = cross(avec,tbar)
	dprd = dot(avec,tbar)

	num = cprd* (1.0-(dprd/(amag)))
	den = dot(cprd,cprd)

	#for the point coinciding with the element
	if (den<1e-6)
		q = [0.0,0.0,0.0]
	else
		q =  num/den
	end
	return q
end
"""
    calc_vind_finite(rbar,sbar,ebar)

    Function to find the three velocity componenets for downwash given the line segment and the point

	Input :
		rbar - vector for pointing the location of evaluation
				: form - Array{Float,1}, size {3}
		sbar - vector for pointing the start of filament
				: form - Array{Float,1}, size {3}
		ebar - vector for pointing the end of filament
				: form - Array{Float,1}, size {3}

    Output:
		q    - vector of influence at the given location
				: form - Array{Float,1}, size {3}
"""
function calc_q_finite(rbar,sbar,ebar)

	(x1,y1,z1) = sbar
	(x2,y2,z2) = ebar
	(x,y,z)    = rbar

	avec = [x1-x,y1-y,z1-z]
	bvec = [x2-x,y2-y,z2-z]

	amag = norm(avec)
	bmag = norm(bvec)

	cprd = cross(avec,bvec)
	dprd = dot(avec,bvec)

	num = cprd* (amag+bmag) * (1.0-(dprd/(amag*bmag)))
	den = dot(cprd,cprd)

	#for the point coinciding with the element
	if (den<1e-6)
		q = [0.0,0.0,0.0]
	else
		q =  num/den
	end
	return q
end
"""
    calc_AICs()

    Calculate the AIC matrices required
"""
function calc_AICs(ntot_lat,sbar,ebar,cbar,mbar,nbar,tbar,debug=0)

	pi_inv = 1.0/pi

        AIC  = zeros(ntot_lat,ntot_lat)

		#Other matrices for CVLM
		Gxgamma = zeros(ntot_lat,ntot_lat)
		Gygamma = zeros(ntot_lat,ntot_lat)
		Gzgamma = zeros(ntot_lat,ntot_lat)
		Hxgamma = zeros(ntot_lat,ntot_lat)
		Hygamma = zeros(ntot_lat,ntot_lat)
		Hzgamma = zeros(ntot_lat,ntot_lat)

        AICₘ = zeros(ntot_lat,ntot_lat)
        AICₜ = zeros(ntot_lat,ntot_lat)

	xyz_ref = zeros(3,1)
	ubar = ebar - sbar
	r_bk = cbar - (ones(ntot_lat,1)*xyz_ref')'
	r_bj = mbar - (ones(ntot_lat,1)*xyz_ref')'
	for i in 1:ntot_lat
		for j in 1:ntot_lat

			#Calculating AIC
			#bounded vortex contribution
			qbv     = calc_q_finite(cbar[:,i],sbar[:,j],ebar[:,j])
			qbv_comp= dot(qbv,nbar[:,i])
			#left leg contribution
			qlv     = -calc_q_sinf(cbar[:,i],sbar[:,j],tbar[:,j])
			qlv_comp= dot(qlv,nbar[:,i])
			#right leg contribution
			qrv     = calc_q_sinf(cbar[:,i],ebar[:,j],tbar[:,j])
			qrv_comp= dot(qrv,nbar[:,i])

			AIC[i,j] = 0.25*pi_inv*(qbv_comp+qlv_comp+qrv_comp)
			qbv .= 0.0
			qbv_comp = 0.0
			qlv .= 0.0
			qlv_comp = 0.0
			qrv .= 0.0
			qrv_comp = 0.0

			#Calculating AICₘ
			#bounded vortex contribution
			qbv     = calc_q_finite(mbar[:,i],sbar[:,j],ebar[:,j])
			qbv_comp= dot(qbv,-nbar[:,i])
			#left leg contribution
			qlv     = -calc_q_sinf(mbar[:,i],sbar[:,j],tbar[:,j])
			qlv_comp= dot(qlv,-nbar[:,i])
			#right leg contribution
			qrv     = calc_q_sinf(mbar[:,i],ebar[:,j],tbar[:,j])
			qrv_comp= dot(qrv,-nbar[:,i])

			h         = 0.25*pi_inv*(qbv+qlv+qrv)
			AICₘ[i,j] = 0.25*pi_inv*(qbv_comp+qlv_comp+qrv_comp)

			Gxgamma[i,j] = dot(cross(h,ubar[:,j]),[1.0 0.0 0.0])
			Gygamma[i,j] = dot(cross(h,ubar[:,j]),[0.0 1.0 0.0])
			Gzgamma[i,j] = dot(cross(h,ubar[:,j]),[0.0 0.0 1.0])
			Hxgamma[i,j] = dot(cross(r_bj[:,j],cross(h,ubar[:,j])),[1.0 0.0 0.0])
			Hygamma[i,j] = dot(cross(r_bj[:,j],cross(h,ubar[:,j])),[0.0 1.0 0.0])
			Hzgamma[i,j] = dot(cross(r_bj[:,j],cross(h,ubar[:,j])),[0.0 0.0 1.0])

			qbv .= 0.0
			qbv_comp = 0.0
			qlv .= 0.0
			qlv_comp = 0.0
			qrv .= 0.0
			qrv_comp = 0.0

			#Calculating AICₜ
			#coordinate transformation for calculation purposes
			c0 = cbar[:,i]
			s0 = sbar[:,j]
			e0 = ebar[:,j]

			c0[1] = s0[1] = e0[1] = 0.0

			#bounded vortex contribution
			qbv     = 0.0
			#left leg contribution
			qlv     = -calc_q_sinf(c0,s0,tbar[:,j])
			qlv_comp= dot(qlv,nbar[:,i])
			#right leg contribution
			qrv     = calc_q_sinf(c0,e0,tbar[:,j])
			qrv_comp= dot(qrv,nbar[:,i])

			AICₜ[i,j] = 0.25*pi_inv*(qbv_comp+2.0*qlv_comp+2.0*qrv_comp)
		end
	end

	Gxv = zeros(ntot_lat,3)
	Gyv = zeros(ntot_lat,3)
	Gzv = zeros(ntot_lat,3)
	Hxv = zeros(ntot_lat,3)
	Hyv = zeros(ntot_lat,3)
	Hzv = zeros(ntot_lat,3)

	Gxomega = zeros(ntot_lat,3)
	Gyomega = zeros(ntot_lat,3)
	Gzomega = zeros(ntot_lat,3)
	Hxomega = zeros(ntot_lat,3)
	Hyomega = zeros(ntot_lat,3)
	Hzomega = zeros(ntot_lat,3)

	crsbar  = zeros(3,ntot_lat)

	for i in 1:ntot_lat
		crsbar[:,i] = cross(r_bk[:,i],nbar[:,i])

		Gxv[i,:] = cross(ubar[:,i],[1.0;0.0;0.0])
		Gyv[i,:] = cross(ubar[:,i],[0.0;1.0;0.0])
		Gzv[i,:] = cross(ubar[:,i],[0.0;0.0;1.0])

		Gxomega[i,:] = cross(r_bj[:,i],Gxv[i,:])
		Gyomega[i,:] = cross(r_bj[:,i],Gyv[i,:])
		Gzomega[i,:] = cross(r_bj[:,i],Gzv[i,:])

		Hxv[i,:] = cross(-ubar[:,i],cross(r_bj[:,i],[1.0;0.0;0.0]))
		Hyv[i,:] = cross(-ubar[:,i],cross(r_bj[:,i],[0.0;1.0;0.0]))
		Hzv[i,:] = cross(-ubar[:,i],cross(r_bj[:,i],[0.0;0.0;1.0]))

		Hxomega[i,:] = cross(-r_bj[:,i],cross(ubar[:,i],cross(r_bj[:,i],[1.0;0.0;0.0])))
		Hyomega[i,:] = cross(-r_bj[:,i],cross(ubar[:,i],cross(r_bj[:,i],[0.0;1.0;0.0])))
		Hzomega[i,:] = cross(-r_bj[:,i],cross(ubar[:,i],cross(r_bj[:,i],[0.0;0.0;1.0])))
	end

	J = AIC\[nbar' crsbar']

	P_x = J'*(-[Gxv Gxomega]+Gxgamma*J)
	P_y = J'*(-[Gyv Gyomega]+Gygamma*J)
	P_z = J'*(-[Gzv Gzomega]+Gzgamma*J)

	Q_x = J'*(-[Hxv Hxomega]+Hxgamma*J)
	Q_y = J'*(-[Hyv Hyomega]+Hygamma*J)
	Q_z = J'*(-[Hzv Hzomega]+Hzgamma*J)

    return AIC,AICₘ,AICₜ,J,P_x,P_y,P_z,Q_x,Q_y,Q_z,Gxgamma,Gygamma,
	Gzgamma,Hxgamma,Hygamma,Hzgamma,Gxv,Gyv,Gzv,Hxv,Hyv,Hzv,Gxomega,
	Gyomega,Gzomega,Hxomega,Hyomega,Hzomega,crsbar
end
"""
    motion_controller()

	Calculating the deflection on wing2/tail for stabilising

	Output:
		def - deflection in radian
"""
function motion_controller()

	def = 0.01
	return def
end
"""
    calc_fm()

	Calculating the forces and moments from CVLM

	Input:
		Px -
		Py -
		Pz -
		xst - state vector of u,v,w,p,q,r,phi,theta,psi,x,y,z
		rho - density

	Output:
		Fx -
		Fy -
		Fz -

"""
function calc_fm(Px,Py,Pz,xst,rho)

	u,v,w,p,q,r,phi,theta,psi,x,y,z = xst

	X = [u;v;w;p;q;r]

	Fx = rho*X'*Px*X
	Fy = rho*X'*Py*X
	Fz = rho*X'*Pz*X

	return Fx,Fy,Fz
end

"""
    calc_motion()

	Calculating the 6DOF of motion

	Input:

	Output:
"""
function calc_motion(Px,Py,Pz,Qx,Qy,Qz,param_file,lat2wing_map,ntot_lat,sbar,ebar,cbar,mbar,nbar,tbar)

	inp_data = YAML.load_file(param_file)
	@info "Reading Motion Parameters Input File"

	#condition to determine control
	ictrl = inp_data["ictrl"]

	# time span
	t_begin = inp_data["t_begin"]
	t_end   = inp_data["t_end"]
	tsteps  = inp_data["tsteps"]
	tspan = (t_begin,t_end)

	##########################################
	# initial conditions
	##########################################
	u0 = inp_data["u0"]
	v0 = inp_data["v0"]
	w0 = inp_data["w0"]

	p0 = inp_data["p0"]
	q0 = inp_data["q0"]
	r0 = inp_data["r0"]

	phi0 = inp_data["phi0"]
	theta0 = deg2rad(inp_data["theta0"])
	psi0 = inp_data["psi0"]

	x0 = inp_data["x0"]
	y0 = inp_data["y0"]
	z0 = inp_data["z0"]

	###############
	# Let's get it running with CVLM in the loop - then, we will need to come up with a trim routine here for steady level flight. I've done this before, so I can add it.
	##############

	# inital condition vector [tranlational velocities | angular velocities | euler angles | cartesian position ]
	xst0 = [u0;v0;w0;p0;q0;r0;phi0;theta0;psi0;x0;y0;z0]

	#############################################
	# parameter definitions
	#############################################
	m = inp_data["m"]
	g = inp_data["g"]
	rho = inp_data["rho"]
	T = inp_data["T"] # thrust (will need to determine from induced drag at trim?)
	# Inertia constants for calculating GAMMA, C1-C9
	Ix  = inp_data["Ix"]
	Iy  = inp_data["Iy"]
	Iz  = inp_data["Iz"]
	Ixz = inp_data["Ixz"]

	# constants to deal with inertial coupling
	GAMMA = Ix*Iz-Ixz^2
	c1 = ((Iy-Iz)*Iz-Ixz^2)/GAMMA
	c2 = ((Ix-Iy+Iz)*Ixz)/GAMMA
	c3 =  Iz/GAMMA
	c4 = Ixz/GAMMA
	c5 = (Iz-Ix)/Iy
	c6 = Ixz/Iy
	c7 = 1/Iy
	c8 = ((Ix-Iy)*Ix-Ixz^2)/GAMMA
	c9 = Ix/GAMMA

	# rotating mass along longitudinal axis (spinning propellor or turbomachinery) - these will probably always be zero for using
	Ip     = inp_data["Ip"]  # rotating inertia
	OmegaP = inp_data["OmegaP"]  # rpm of rotating machinery

	##################################################
	# Define parameter vector (would be more convienient as a tuple)
	##################################################
	pq = [m,g,rho,c1,c2,c3,c4,c5,c6,c7,c8,c9,T,Ip,OmegaP]

	# Controls On
	if ictrl == 1

		def = motion_controller()
		#Rotate the tail by some specification
		rot          = [cos(def) 0.0 sin(def);0.0 1.0 0.0;-sin(def) 0.0 cos(def)]

		nbar_new = copy(nbar)
		#
		# Can be improved Immensely
		#
		for ilat in 1:ntot_lat

			if lat2wing_map[ilat] == 2
				nbar_new[:,ilat] = rot*nbar[:,ilat]
			end
		end

		s,s,s,s,Px,Py,Pz,Qx,Qy,Qz, = calc_AICs(ntot_lat,sbar,ebar,cbar,mbar,nbar_new,tbar)

	end

	# ***This is a function inside function***
	## ODE Function
	function calc_dx!(dxst,xst,pq,t)

	    #X = [u,v,w,p,q,r,phi,theta,psi,x,y,z]

	    # parameters:
	    m,g,rho,c1,c2,c3,c4,c5,c6,c7,c8,c9,T,Ip,OmegaP = pq
	    # current state values:
	    u,v,w,p,q,r,phi,theta,psi,x,y,z = xst

	    # Aerodynamic Force/Moment Calculation (without controls {for now})
	    X = [u;v;w;p;q;r]

	    Fx = rho*X'*Px*X
	    Fy = rho*X'*Py*X
	    Fz = rho*X'*Pz*X

	    Mx = rho*X'*Qx*X
	    My = rho*X'*Qy*X
	    Mz = rho*X'*Qz*X

	    # state derivative calculations
	    #dxst = zeros(12)
	    dxst[1] = r*v-q*w+Fx/m-g*sin(theta)+T/m
	    dxst[2] = p*w-r*u+Fy/m+g*cos(theta)*sin(phi)  # Force Equations
	    dxst[3] = q*u-p*v+Fz/m+g*cos(theta)*cos(phi)

	    dxst[4] = (c1*r+c2*p-c4*Ip*OmegaP)*q+c3*Mx+c4*Mz
	    dxst[5] = (c5*p+c7*Ip*OmegaP)*r-c6*(p^2-r^2)+c7*My   # moment equations
	    dxst[6] = (c8*p-c2*r-c9*Ip*OmegaP)*q+c9*Mz+c4*Mx

	    dxst[7] = p+tan(theta)*(q*sin(phi)+r*cos(phi))
	    dxst[8] = q*cos(phi)-r*sin(phi)
	    dxst[9] = (q*sin(phi)+r*cos(phi))/cos(theta)

	    dxst[10:12] = [1 0 0; 0 cos(phi) sin(phi);0 -sin(phi) cos(phi)]*[cos(theta) 0 -sin(theta);0 1 0; sin(theta) 0 cos(theta)]*[cos(phi) sin(phi) 0;-sin(phi) cos(phi) 0; 0 0 1]*X[1:3]

	end

	prob = ODEProblem(calc_dx!,xst0,tspan,pq)

	#nsteps = 111
	#tstops = LinRange(tspan[1],tspan[2],nsteps)
	#sol = solve(prob,tstops=tstops,adaptive=false)
	#sol  = solve(prob,Rodas4())
	#sol  = solve(prob,Rosenbrock23())
	sol = solve(prob,Rodas4(),dtmin=1e-6,saveat=LinRange(t_begin,t_end,tsteps))
	return sol
end
"""
    main(inp_file)

	Main driver function to calculate all the values required

	Input :
	inp_file - Input File

	Output:
	sol - solution for motion
"""
function main(inp_file,param_file,iseq=0)

    sref,bref,cref,aseq,vinf,ρ,alpha,nwing,
	ntot_patch,npatch_perwing,
	ntot_lat,nspan_perpatch,nlat_perpatch,xyz_qc_patch,
	chord_patch,twist_patch,α_zl_patch = read_inp(inp_file)

	nspn,spn_map,spn_loc,lat2wing_map,sbar,ebar,mbar,nbar,cbar,dbar,tbar,chord,twist,α_zl,θ,Λ,ds = geom_calc(nwing,ntot_patch,npatch_perwing,ntot_lat,nspan_perpatch,nlat_perpatch,xyz_qc_patch,chord_patch,twist_patch,α_zl_patch)

	AIC,AICₘ,AICₜ,J,Px,Py,Pz,Qx,Qy,Qz,Gxgamma,Gygamma,
	Gzgamma,Hxgamma,Hygamma,Hzgamma,Gxv,Gyv,Gzv,Hxv,Hyv,Hzv,Gxomega,
	Gyomega,Gzomega,Hxomega,Hyomega,Hzomega,crsbar = calc_AICs(ntot_lat,sbar,ebar,cbar,mbar,nbar,tbar)

	sol = calc_motion(Px,Py,Pz,Qx,Qy,Qz,param_file,lat2wing_map,ntot_lat,sbar,ebar,cbar,mbar,nbar,tbar)

	return sol,Px,Py,Pz,Qx,Qy,Qz,lat2wing_map
end

max_stars_repo_name: rohitvuppala/CVLM

max_issues_repo_name: rohitvuppala/CVLM

max_forks_repo_name: rohitvuppala/CVLM


text:
"""
Thread-safe logging.
"""
module Log

import Base.Threads.threadid

# logging levels in increasing verbosity
@enum LogLevel ERROR WARN INFO DEBUG

const LEVEL = Ref{LogLevel}(INFO)
const VERBOSE = Ref{Bool}(true)  # print stack traces from errors

# Rank for multinode functionality (can be set on startup)
const rank = Ref{Int}(1)
grank() = rank[]

# thread-safe print function
@inline function puts(s...)
    data = string(s..., '\n')
    ccall(:write, Cint, (Cint, Cstring, Csize_t), 1, data, sizeof(data))
end
@inline rtputs(s...) = puts("[$(grank())]<$(threadid())> ", s...)

# logging functions
@inline error(msg...) = LEVEL[] >= ERROR && rtputs("ERROR: ", msg...)
@inline warn(msg...) = LEVEL[] >= WARN && rtputs("WARN: ", msg...)
@inline info(msg...) = LEVEL[] >= INFO && rtputs("INFO: ", msg...)
@inline debug(msg...) = LEVEL[] >= DEBUG && rtputs("DEBUG: ", msg...)

"""
    exception(ex::Exception, msg...)

If the logging level is ERROR or greater, log the message `msg` and
show the exception `ex`. If `Log.VERBOSE[]` is true, also log the
stack trace. Should only be called from a `catch` block, e.g.,

```
try
    ...
catch ex
    Log.exception(ex, "Something", " happened.")
end
```
"""
function exception(exception::Exception, msg...)
    LEVEL[] >= ERROR || return
    if length(msg) > 0
        error(msg...)
    end
    if VERBOSE[]
        stack_trace = catch_stacktrace()
    end
    buf = IOBuffer()
    Base.showerror(buf, exception)
    error(String(take!(buf)))
    if VERBOSE[]
        error("Stack trace:")
        if length(stack_trace) > 100
            stack_trace = vcat(
                [string(line) for line in stack_trace[1:50]],
                @sprintf("...(removed %d frames)...", length(stack_trace) - 100),
                [string(line) for line in stack_trace[(length(stack_trace) - 50):length(stack_trace)]],
            )
        end
        for stack_line in stack_trace
            error(@sprintf("  %s", stack_line))
        end
    end
end

end


text:
# ---
# title: 565. Array Nesting
# id: problem565
# author: Indigo
# date: 2021-06-26
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/array-nesting/description/>
# hidden: true
# ---
#
# A zero-indexed array A of length N contains all integers from 0 to N-1. Find
# and return the longest length of set S, where S[i] = {A[i], A[A[i]],
# A[A[A[i]]], ... } subjected to the rule below.
#
# Suppose the first element in S starts with the selection of element A[i] of
# index = i, the next element in S should be A[A[i]], and then A[A[A[i]]]… By
# that analogy, we stop adding right before a duplicate element occurs in S.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [5,4,0,3,1,6,2]
#     Output: 4
#     Explanation:
#     A[0] = 5, A[1] = 4, A[2] = 0, A[3] = 3, A[4] = 1, A[5] = 6, A[6] = 2.
#
#     One of the longest S[K]:
#     S[0] = {A[0], A[5], A[6], A[2]} = {5, 6, 2, 0}
#
#
#
#
# **Note:**
#
#   1. N is an integer within the range [1, 20,000].
#   2. The elements of A are all distinct.
#   3. Each element of A is an integer within the range [0, N-1].
#
#
## @lc code=start
using LeetCode

function array_nesting(nums::Vector{Int})
    visited = fill(false, length(nums))
    res = 0
    for (idx, num) in enumerate(nums)
        num += 1
        visited[idx] && continue
        tmp = 1
        visited[idx] = true
        while !visited[num]
            visited[num] = true
            num = nums[num] + 1
            tmp += 1
        end
        res = max(res, tmp)
    end
    return res
end
## @lc code=end


text:
module UniformScalingMatrixFunctions
using LinearAlgebra

# Matrix/operator functions f(A) are defined as either the Taylor
# series of f evaluated using the matrix A, or equivalently, as f(A) =
# Q*f(λ)*Qᴴ, where Q is the matrix of eigenvectors of A and λ the
# associated eigenvalues. UniformScaling represents linear operators
# of the form α*I, where α is a scalar, which simultaneously is the
# only eigenvalue. A matrix/operator function f(α*I) is thus simply
# f(α)*I, which is also valid from the Taylor series perspective,
# since I^k==I.
for f in [:exp,
          :sin, :cos, :tan,
          :asin, :acos, :atan,
          :sinh, :cosh, :tanh,
          :asinh, :acosh, :atanh,
          :log, :sqrt]
    @eval begin
        Base.$f(J::UniformScaling) = UniformScaling($(f)(J.λ))
    end
end

export I

end # module


text:
using SemialgebraicGames
using Test

tests = ["system_univariate.jl", "system_multivariate.jl"]

for t in tests
    include(t)
end