text:
"""
Bridges `CP.Inverse` to reification.
"""
struct Inverse2ReificationBridge{T} <: MOIBC.AbstractBridge
    vars_first::Matrix{MOI.VariableIndex}
    vars_first_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    vars_second::Matrix{MOI.VariableIndex}
    vars_second_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_first_reif::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}
    cons_second_reif::Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}
    cons_equivalence::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{Inverse2ReificationBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Inverse,
) where {T}
    return MOIBC.bridge_constraint(
        Inverse2ReificationBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Inverse2ReificationBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Inverse,
) where {T}
    f_scalars = MOIU.scalarize(f)
    f_x = f_scalars[1:s.dimension]
    f_y = f_scalars[(s.dimension + 1):(2 * s.dimension)]

    vars_first = Matrix{MOI.VariableIndex}(undef, s.dimension, s.dimension)
    vars_first_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.dimension, s.dimension)
    vars_second = Matrix{MOI.VariableIndex}(undef, s.dimension, s.dimension)
    vars_second_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, s.dimension, s.dimension)

    for i in 1:s.dimension
        vars_first[i, :], vars_first_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.dimension])
        vars_second[i, :], vars_second_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.dimension])
    end

    #     x_i    =    j     ⟺    y_j    =    i
    # \___ vars1[i, j] ___/    \___ vars2[i, j] ___/

    cons_first_reif = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}(undef, s.dimension, s.dimension)
    cons_second_reif = Matrix{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}}(undef, s.dimension, s.dimension)
    cons_equivalence = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}(undef, s.dimension, s.dimension)

    for i in 1:s.dimension
        for j in 1:s.dimension
            cons_first_reif[i, j] = MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_first[i, j],
                        f_x[i] - T(j)
                    ]
                ),
                CP.Reification(MOI.EqualTo(zero(T)))
            )

            cons_second_reif[i, j] = MOI.add_constraint(
                model,
                MOIU.vectorize(
                    [
                        one(T) * vars_second[i, j],
                        f_y[j] - T(i)
                    ]
                ),
                CP.Reification(MOI.EqualTo(zero(T)))
            )

            cons_equivalence[i, j] = MOI.add_constraint(
                model,
                one(T) * vars_first[i, j] - one(T) * vars_second[i, j],
                MOI.EqualTo(zero(T))
            )
        end
    end

    return Inverse2ReificationBridge(vars_first, vars_first_bin, vars_second, vars_second_bin, cons_first_reif, cons_second_reif, cons_equivalence)
end

function MOI.supports_constraint(
    ::Type{Inverse2ReificationBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Inverse},
) where {T}
    return true
    # Ideally, ensure that the underlying solver supports all the needed
    # reified constraints.
end

function MOIB.added_constrained_variable_types(::Type{Inverse2ReificationBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{Inverse2ReificationBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::Inverse2ReificationBridge{T}, ::MOI.NumberOfVariables) where {T}
    return length(b.vars_first) + length(b.vars_second)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
return length(b.vars_first_bin) + length(b.vars_second_bin)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_first_reif) + length(b.cons_second_reif)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons_equivalence)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return vcat(b.vars_first, b.vars_second)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return vcat(b.vars_first_bin, b.vars_second_bin)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return vcat(b.cons_first_reif, b.cons_second_reif)
end

function MOI.get(
    b::Inverse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return copy(b.cons_equivalence)
end


text:
# caller_name.jl

using MIRT: caller_name

using Test: @test

function caller_name_test()

	function f2()
		caller_name()
	end

	line = 2 + @__LINE__
	function f1()
		f2() # this is two lines below @__LINE__ above
	end

	@test isa(f1(), String)
	@test f1()[end-12:end] == ".jl $line f1(): "

	true
end

@test caller_name_test()

ext: jl


text:
using HeavyGas
using Interpolations

const dt = 0.01

m = GasData(6.4, 10, 0.003, 20+273.15, 44.10, HeavyGas.D)

x, vfrac = instantaneous(7., 4.12, 1., dt, 28.:20.:128., m, 2)
println(x,"\n",vfrac)


text:
using SparseArrays
using FluxExtensions: mul, multrans

@testset "testing multiplication of sparse matrix with a mask on rows" begin
	A = [1 2 3; 0 1 0]
	B = sparse([0 1 2; 1 0 1; 1 0 0])
	@test mul(A, B) == A*B
	@test mul(A, B) == mul(A,B, [1,2,3])
	Bz = sparse([0 0 0; 1 0 1; 1 0 0])
	@test mul(A, Bz) == mul(A,B, [2,3])
	Bz = sparse([0 1 2; 0 0 0; 1 0 0])
	@test mul(A, Bz) == mul(A,B, [1,3])
	Bz = sparse([0 1 2; 1 0 1; 0 0 0])
	@test mul(A, Bz) == mul(A,B, [1,2])
end

@testset "testing sparse hcat" begin
	x = SparseArrays.sprand(3, 10, 0.5)
	y = SparseArrays.sprand(3, 10, 0.5)
	@test hcat(Matrix(x), Matrix(y)) ≈ reduce(hcat, [x, y])
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
export get_data!, get_data

@doc """
    S = get_data(method, channels, KWs)
    get_data!(S, method, channels, KWs)

## Arguments
* `method`: String retrieval method
  + `"FDSN"`: FDSNWS dataselect
  + `"IRIS"`: IRISWS timeseries
* `channels`: channel selectors. (Union{String,Array{String,1},Array{String,2}})
* `S`: SeisData structure

## Keywords
### Standard keywords
fmt, nd, opts, rad, reg, si, src, to, v, w, y

### Seismic Processing keywords
Certain keywords apply processing steps to seismic data channels (including
pressure sensors). The relevant keywords and order of operations is as follows:

1. unscale: divide gain from data after download?
2. demean: demean data after download?
3. detrend: detrend data after download?
4. taper: taper data after download?
6. ungap: remove gaps in data after download?
5. rr: remove seismic instrument response after download?

None of these keywords are mutually exclusive. Hey, it's your CPU time.

### Other keywords
* `autoname`: Determine file names from channel ID?
* `msr`: Get multi-stage response? (FDSN only)
* `s`: Start time
* `t`: Termination (end) time
* `xf`: Name of XML file to save station metadata

#### Writing to disk and file names
KW `autoname=true` sets `w=true` automatically. In both cases, each valid
request is written directly to disk before parsing; a file persists even if
parsing fails.

If method is "IRIS", file names are uniquely specified if `w=true`. KW `autoname`
is not needed.

If method is "FDSN", then `autoname=true` attempts to emulate IRIS channel file
naming conventions. A major changes to request syntax is needed for this to
work, however: each request must return *exactly one* channel.

Example:
* `get_data("FDSN", "UW.LON..BHZ", autoname=true)` generates IRIS-style
filenames because the channel name is uniquely specified.
* `get_data("FDSN", "UW.LON..BH?", autoname=true)` still writes to disk, but
can't use IRIS-style file names because the request returns three channels.

## Best Practices
* `get_data("IRIS", ...)` is preferable to `get_data("FDSN", ..., src="IRIS")`
for a single-channel request available from both services.
  + Note: not all stations in IRIS metadata are available from both services.
* Avoid wildcards with `autoname=true`.

## Examples
1. `S = get_data("FDSN", "UW.SEP..EHZ,UW.SHW..EHZ,UW.HSR..EHZ", src="IRIS",
t=(-600))`: using FDSNWS, get the last 10 minutes of data from three short-
period vertical-component channels at Mt. St. Helens, USA.
2. `get_data!(S, "IRIS", "CC.PALM..EHN", t=(-120), f="sacbl", demean=true)`:
using IRISWS, fetch the last two minutes of data from component EHN, station
PALM (Palmer Lift (Mt. Hood), OR, USA,), network CC (USGS Cascade Volcano
Observatory, Vancouver, WA, USA), in bigendian SAC format, and append to
SeisData structure `S`. Remove the mean from the request after download and
before appending.
3. `S = get_data("FDSN", "CC.TIMB..EHZ", t=(-600), autoname=true)`: using FDSNWS,
get the last 10 minutes of data from channel EHZ, station TIMB (Timberline Lodge,
OR, USA). Save directly to disk with an automated file name, in addition to
returning data in a new SeisData structure S.
4. `S = get_data("FDSN", "HV.MOKD..HHZ", s="2012-01-01T00:00:00", t=(-3600))`:
using FDSNWS, create a new SeisData structure `S` with an hour of data ending
at 2012-01-01, 00:00:00 UTC, from HV.MOKD..HHZ (USGS Hawai'i Volcano Observatory).

See also: chanspec, parsetimewin, seis_www, SeisIO.KW
""" get_data
function get_data(method_in::String, C="*"::Union{String,Array{String,1},Array{String,2}};
        autoname::Bool = false                               ,  # Auto-generate file names?
          demean::Bool = false                               ,  # Demean data after download?
         detrend::Bool = false                               ,  # Detrend data after download?
           fmt::String = KW.fmt                              ,  # File format
             msr::Bool = false                               ,  # Get multi-stage response?
              nd::Real = KW.nd                               ,  # Number of days per request (in long requests)
          opts::String = KW.opts                             ,  # Options string
 rad::Array{Float64,1} = KW.rad                              ,  # Query radius
 reg::Array{Float64,1} = KW.reg                              ,  # Query region
           prune::Bool = KW.prune                            ,  # Prune empty channels after query?
              rr::Bool = false                               ,  # Remove instrument response after download?
           s::TimeSpec = 0                                   ,  # Start
              si::Bool = KW.si                               ,  # Fill station info?
           src::String = KW.src                              ,  # Data source
           taper::Bool = false                               ,  # Taper data after download?
           t::TimeSpec = (-600)                              ,  # End or Length (s)
             to::Int64 = KW.to                               ,  # Timeout (s)
           ungap::Bool = false                               ,  # Remove data gaps?
         unscale::Bool = false                               ,  # Unscale (divide by :gain)?
              v::Int64 = KW.v                                ,  # Verbosity
               w::Bool = KW.w                                ,  # Write to disc?
            xf::String = "FDSNsta.xml"                       ,  # XML save file
               y::Bool = KW.y                                   # Sync
     )

  # Parse time window
  α, ω = parsetimewin(s, t)

  # KWs that ovewrite other KWs (rare; keep this behavior to an absolute minimum!)
  (autoname == true) && (w = true)

  # Generate SeisData
  S = SeisData()

  # Condense requests as much as possible
  if method_in == "FDSN"
    if isa(C, String)
      C = parse_chstr(C, fdsn = true)
    elseif isa(C, Array{String,1})
      C = parse_charr(C, fdsn = true)
    end
    R = minreq(C)
    parse_err = FDSNget!(S, C,
                          autoname=autoname,
                          fmt=fmt,
                          msr=msr,
                          nd=nd,
                          opts=opts,
                          rad=rad,
                          reg=reg,
                          s=s,
                          si=si,
                          src=src,
                          t=t,
                          to=to,
                          v=v,
                          w=w,
                          xf=xf,
                          y=y)
  elseif method_in == "IRIS"
    if isa(C, String)
      R = String[strip(String(j)) for j in split(C, ',')]
    elseif isa(C, Array{String,2})
      NC = size(C,1)
      R = Array{String,1}(undef, NC)
      for j = 1:NC
        R[j] = join(C[j,:],'.')
      end
    else
      R = deepcopy(C)
    end
    parse_err = IRISget!(S, R, α, ω, fmt = fmt, opts = opts, to = to, v = v, w = w)
  end

  # DND DND DND
  # Wrapper to a generic "get" function -- leave as example code
  # getfield(SeisIO, Symbol(string(method_in, "get!")))(S, C, α, ω,
  #          f = f, opts = opts, si = si, src = src, to = to, v = v, w = w)
  # DND DND DND

  # ========================================================================
  # Viable operations for any data type
  if prune == true
    if parse_err == false
      v > 0 && @info(tnote("Removing empty channels."))
      prune!(S)
    else
      v > 0 && @info(tnote("Can't prune empty channels; web request wasn't fully parsed."))
    end
  end

  # Sync
  if y == true
    v > 0 && @info(tnote("Synchronizing data."))
    sync!(S)
  end

  # ========================================================================
  # Operations that only make sense for seismic or seismoacoustic data

  if any([demean, detrend, rr, taper, unscale])
    # Get list of channels with sane instrument codes
    CC = get_seis_channels(S)

    # Unscale
    if unscale == true
      v > 0 && @info(tnote("dividing out gain of seismic channels"))
      unscale!(S, chans=CC)
    end

    # Demean
    if demean == true
      v > 0 && @info(tnote("removing mean of seismic channels"))
      demean!(S, chans=CC)
    end

    # Detrend
    if detrend == true
      v > 0 && @info(tnote("detrending seismic channels"))
      detrend!(S, chans=CC)
    end

    # Taper
    if taper == true
      v > 0 && @info(tnote("tapering seismic channel data"))
      taper!(S, chans=CC)
    end

    # Ungap
    if ungap == true
      v > 0 && @info(tnote("ungapping seismic channel data"))
      ungap!(S, chans=CC)
    end

    # Remove response
    if rr == true
      v > 0 && @info(tnote("removing instrument response from seismic channels"))
      remove_resp!(S, chans=CC)
    end
  end

  return S
end

@doc (@doc get_data)
function get_data!(S::SeisData, method_in::String, C="*"::Union{String,Array{String,1},Array{String,2}};
        autoname::Bool = false                               ,  # Auto-generate file names?
          demean::Bool = false                               ,  # Demean data after download?
         detrend::Bool = false                               ,  # Detrend data after download?
           fmt::String = KW.fmt                              ,  # File format
             msr::Bool = false                               ,  # Get multi-stage response?
              nd::Real = KW.nd                               ,  # Number of days per request (in long requests)
          opts::String = KW.opts                             ,  # Options string
 rad::Array{Float64,1} = KW.rad                              ,  # Query radius
 reg::Array{Float64,1} = KW.reg                              ,  # Query region
           prune::Bool = KW.prune                            ,  # Prune empty channels after query?
              rr::Bool = false                               ,  # Remove instrument response after download?
           s::TimeSpec = 0                                   ,  # Start
              si::Bool = KW.si                               ,  # Fill station info?
           src::String = KW.src                              ,  # Data source
           taper::Bool = false                               ,  # Taper data after download?
           t::TimeSpec = (-600)                              ,  # End or Length (s)
             to::Int64 = KW.to                               ,  # Timeout (s)
           ungap::Bool = true                                ,  # Remove data gaps?
         unscale::Bool = true                                ,  # Unscale (divide by :gain)?
              v::Int64 = KW.v                                ,  # Verbosity
               w::Bool = KW.w                                ,  # Write to disc?
            xf::String = "FDSNsta.xml"                       ,  # XML save file
               y::Bool = KW.y                                   # Sync
     )

  U = get_data(method_in, C,
                autoname=autoname,
                demean=demean,
                detrend=detrend,
                fmt=fmt,
                msr=msr,
                nd=nd,
                opts=opts,
                prune=prune,
                rad=rad,
                reg=reg,
                rr=rr,
                s=s,
                si=si,
                src=src,
                t=t,
                taper=taper,
                to=to,
                ungap=ungap,
                unscale=unscale,
                v=v,
                w=w,
                xf=xf,
                y=y)
  v > 2 && println(stdout, "S = \n", U)
  append!(S, U)
  return nothing
end


text:
struct LinearPieces
    left::Float64
    right::Float64
    slope::Float64
    intercept::Float64
end

#PiecewiseLinearActivation
struct PiecewiseLinearActivation <: ActivationFunction
    pieces::Vector{LinearPieces}
end

function (f::PiecewiseLinearActivation)(x)
    for p in f.pieces
        if p.left <= x < p.right
            return p.slope * x + p.intercept
        end
    end
end

l1 = LinearPieces(-Inf, 0.0, 0.0, 0.0)
l2 = LinearPieces(0.0, Inf, 1.0, 0.0)
ReLUact = PiecewiseLinearActivation([l1, l2])
IDact = PiecewiseLinearActivation([LinearPieces(-Inf, Inf, 1.0, 0.0)])

struct PointSearch end

# This is the main function
function solve(solver::PointSearch, problem::Problem)
    layers = problem.network.layers
    input = tovrep(problem.input)
    z1_list = Vector{Vector{Float64}}(undef, 0)
    (W, b, act) = (layers[1].weights, layers[1].bias, layers[1].activation)
    for z0 in points(input)
        z1_hat = affine_map(W, b, z0)
        z1 = act.(z1_hat)
        if is_interior_point(forward_network(solver, problem.network, z1, 1), problem.output)
            push!(z1_list, z1)
        else
            return BasicResult(:violated)
        end
    end
    for i in 1:length(layers)
        (W, b, act) = (layers[i].weights, layers[i].bias, layers[i].activation)
        v_list = vertices_list(input) #z_i-1
        zi_hat_list = affine_map.(v_list)
        zi_hat_tope = VPolytope(zi_hat_list)
        zi_list = Vector{Vector{Float64}}(undef, 0)
        inflections = get_inflection(act) #Vector
        if !isempty(inflections)
            C, d = tosimplehrep(zi_hat_tope)
            n = length(d)
            Ci = [C; zeros(Float64, (2, size(C, 2)))]
            di = [d; 0.0; 0.0]
            for j in 1:length(b)
                if j > 1
                    Ci[n+1,j-1] = 0.0
                    Ci[n+2,j-1] = 0.0
                end
                Ci[n+1,j] = 1.0
                Ci[n+2,j] = -1.0
                for k in length(inflections)
                    di[n+1] = inflections[k]
                    di[n+2] = -inflections[k]
                    intersect = HPolytope(Ci, di)
                    if !isempty(intersect)
                        v_inter = tovrep(intersect)
                        for zi_hat in points(v_inter)
                            zi = act.(zi_hat)
                            if is_interior_point(forward_network(solver, problem.network, zi, i), problem.output)
                                push!(zi_list, zi)
                            else
                                return BasicResult(:violated)
                            end
                        end
                    end
                end
            end
        end
        if i < length(layers)
            if i == 1
                append!(zi_list, z1_list)
                z1_list = 0 #free memory
            else
                append!(zi_list, act.(zi_hat_list))
            end
            input = remove_redundant_vertices(VPolytope(zi_list))
        end
    end
    return BasicResult(:holds)
end

function get_inflection(act::PiecewiseLinearActivation)
    n = length(act.pieces)
    inflections = Vector{Float64}(undef, n-1)
    for i in 1:n-1
        inflections[i] = act.pieces[i].right
    end
    return inflections
end

function get_inflection(act::ReLU)
    return [0.0]
end

function get_inflection(act::Id)
    return []
end

function forward_network(solver::PointSearch, nnet::Network, input::Vector{Float64}, last::Int)
    if last == length(nnet.layers)
        return input
    end
    reach = input
    start = last + 1
    for i in start:length(nnet.layers)
        reach = forward_layer(solver, nnet.layers[i], reach)
    end
    return reach
end

function forward_layer(solver::PointSearch, layer::Layer, input::Vector{Float64})
    return layer.activation.(affine_map(layer.weights, layer.bias, input))
end

function affine_map(W::Matrix{Float64}, b::Vector{Float64}, input::Vector{Float64})
    n = length(b)
    output = Vector{Float64}(undef, n)
    for j in 1:n
        output[j] = sum(W[j,:] .* input) + b[j]
    end
    return output
end


text:
__precompile__()

module MixedVolume

    include("cayley.jl")
    include("circuittable.jl")
    include("mixedcell.jl")
end # module


text:
abstract type AbstractMole end

struct Mole <: AbstractMole
    nuclei::Nuclei
    basis::Basis
    net_charge::Int
    multiplicity::Int
end

function Mole(
    nuclei::Nuclei,
    basis::Union{Dict,String},
    net_charge::Int,
    multiplicity::Int,
    )
    basis = Basis(nuclei, basis)
    return Mole(nuclei, basis, net_charge, multiplicity)
end

function Mole(
    nuclei::String,
    basis::Union{Basis,Dict,String},
    net_charge::Int,
    multiplicity::Int,
    )
    nuclei = Nuclei(nuclei)
    return Mole(nuclei, basis, net_charge, multiplicity)
end

function Mole(atoms, basis)
    return Mole(atoms, basis, 0, 1)
end

abstract type AbstractEnv end

struct Env <: AbstractEnv
    pointcharges::PointCharges
end

function get_basis_idx(nuclei::Nuclei, basis::Basis)
    idx = Int[]
    for i = 1:length(nuclei)
        number = nuclei.numbers[i]
        shells_idx = findfirst(==(number), basis.numbers)
        shells = basis.shells[shells_idx]
        for (j, shell) in enumerate(shells)
            angular_momentum = shell.angular_momentum
            for coefficients in shell.coefficients
                for _ = 1:2*angular_momentum+1
                    push!(idx, i)
                end
            end
        end
    end
    return idx
end

get_basis_idx(mole::Mole) = get_basis_idx(mole.nuclei, mole.basis)


text:
using GenomicTiles
using BioFeatures
using DataFrames

function get_sub_regions(region_id, start, stop)
    subregions = DataFrame( gene_id = [region_id,region_id], start = [start, start + 50], stop = [start+60, start + 100])
    #return Features( subregions, region_id_field = :gene_id  )
    return subregions
end

## example where we do genes and some feature that happens in genes
region_lengths = int( rand(100)*1000 )
starts = 1:1000:100000
stops = starts + region_lengths
gene_id = zeros(ASCIIString,100)
for i=1:100
  gene_id[i] = "GENE$i"
end
# Read these into a BioFeature -
#  in the Description field you have GENE_ID --> LENGTH
#  in the features
regions   = DataFrame(gene_id = gene_id, start = starts, stop= stops)
features  = Features( regions, region_id_field=:gene_id, regions=true )
#features.description
## Features in genes
for i = 1:100
  subfeatures = get_sub_regions("GENE$i",starts[i],stops[i])
  add_dataframe_to_features!(features,subfeatures, region_id_field=:gene_id)
end

features.features
# Then get tiles in one step.
# plot_tiles(features, clustering=kmeans)
# plot_cluster_bar_char(Feature)


text:
import Base.push!  # for DisjointMinSets

"""
```
label = label_components(tf, [connectivity])
label = label_components(tf, [region])
```

Find the connected components in a binary array `tf`. There are two forms that
`connectivity` can take:

- It can be a boolean array of the same dimensionality as `tf`, of size 1 or 3
along each dimension. Each entry in the array determines whether a given
neighbor is used for connectivity analyses. For example, `connectivity = trues(3,3)`
would use 8-connectivity and test all pixels that touch the current one, even
the corners.

- You can provide a list indicating which dimensions are used to
determine connectivity. For example, `region = [1,3]` would not test
neighbors along dimension 2 for connectivity. This corresponds to just
the nearest neighbors, i.e., 4-connectivity in 2d and 6-connectivity
in 3d.

The default is `region = 1:ndims(A)`.

The output `label` is an integer array, where 0 is used for background
pixels, and each connected region gets a different integer index.
"""
label_components(A, connectivity = 1:ndims(A), bkg = 0) = label_components!(zeros(Int, size(A)), A, connectivity, bkg)

#### 4-connectivity in 2d, 6-connectivity in 3d, etc.
# But in fact you can choose which dimensions are connected
let _label_components_cache = Dict{Tuple{Int, Vector{Int}}, Function}()
global label_components!
function label_components!(Albl::AbstractArray{Int}, A::Array, region::Union{Dims, AbstractVector{Int}}, bkg = 0)
    uregion = unique(region)
    if isempty(uregion)
        # Each pixel is its own component
        k = 0
        for i = 1:length(A)
            if A[i] != bkg
                k += 1
                Albl[i] = k
            end
        end
        return Albl
    end
    # We're going to compile a version specifically for the chosen region. This should make it very fast.
    key = (ndims(A), uregion)
    if !haskey(_label_components_cache, key)
        # Need to generate the function
        N = length(uregion)
        ND = ndims(A)
        iregion = [Symbol("i_", d) for d in uregion]
        f! = eval(quote
            local lc!
            function lc!(Albl::AbstractArray{Int}, sets, A::Array, bkg)
                offsets = strides(A)[$uregion]
                @nexprs $N d->(offsets_d = offsets[d])
                k = 0
                @nloops $ND i A begin
                    k += 1
                    val = A[k]
                    label = typemax(Int)
                    if val != bkg
                        @nexprs $N d->begin
                            if $iregion[d] > 1  # is this neighbor in-bounds?
                                if A[k-offsets_d] == val  # if the two have the same value...
                                    newlabel = Albl[k-offsets_d]
                                    if label != typemax(Int) && label != newlabel
                                        label = union!(sets, label, newlabel)  # ...merge labels...
                                    else
                                        label = newlabel  # ...and assign the smaller to current pixel
                                    end
                                end
                            end
                        end
                        if label == typemax(Int)
                            label = push!(sets)   # there were no neighbors, create a new label
                        end
                        Albl[k] = label
                    end
                end
                Albl
            end
        end)
        _label_components_cache[key] = f!
    else
        f! = _label_components_cache[key]
    end
    sets = DisjointMinSets()
    f!(Albl, sets, A, bkg)
    # Now parse sets to find the labels
    newlabel = minlabel(sets)
    for i = 1:length(A)
        if A[i] != bkg
            Albl[i] = newlabel[find_root!(sets, Albl[i])]
        end
    end
    Albl
end
end # let
label_components!(Albl::AbstractArray{Int}, A::BitArray, region::Union{Dims, AbstractVector{Int}}, bkg = 0) = label_components!(Albl, convert(Array{Bool}, A), region, bkg)

#### Arbitrary connectivity
for N = 1:4
    @eval begin
        function label_components!(Albl::AbstractArray{Int,$N}, A::AbstractArray, connectivity::Array{Bool,$N}, bkg = 0)
            if isempty(connectivity) || !any(connectivity)
                # Each pixel is its own component
                k = 0
                for i = 1:length(A)
                    if A[i] != bkg
                        k += 1
                        Albl[i] = k
                    end
                end
                return Albl
            end
            for d = 1:ndims(connectivity)
                (isodd(size(connectivity, d)) && connectivity == flipdim(connectivity, d)) || error("connectivity must be symmetric")
            end
            size(Albl) == size(A) || throw(DimensionMismatch("size(Albl) must be equal to size(A)"))
            @nexprs $N d->(halfwidth_d = size(connectivity,d)>>1)
            @nexprs $N d->(offset_d = 1+halfwidth_d)
            cstride_1 = 1
            @nexprs $N d->(cstride_{d+1} = cstride_d * size(connectivity,d))
            @nexprs $N d->(k_d = 0)
            sets = DisjointMinSets()
            @nloops $N i A begin
                val = @nref($N, A, i)
                label = typemax(Int)
                if val != bkg
                    @nloops($N, j,
                            d->(max(-halfwidth_d, 1-i_d):min(halfwidth_d, size(A,d)-i_d)),  # loop range
                            d->(k_{d-1} = k_d + j_d*cstride_d; # break if linear index is >= 0
                                d > 1 ? (if k_{d-1} > 0 break end) : (if k_0 >= 0 break end);
                                jj_d = j_d+offset_d; ii_d = i_d+j_d),    # pre-expression
                            begin
                        if @nref($N, connectivity, jj)
                            if @nref($N, A, ii) == val
                                newlabel = @nref($N, Albl, ii)
                                if label != typemax(Int) && label != newlabel
                                    label = union!(sets, label, newlabel)
                                else
                                    label = newlabel
                                end
                            end
                        end
                    end)
                    if label == typemax(Int)
                        label = push!(sets)
                    end
                    @nref($N, Albl, i) = label
                end
            end
            # Now parse sets to find the labels
            newlabel = minlabel(sets)
            for i = 1:length(A)
                if A[i]!=bkg
                    Albl[i] = newlabel[find_root!(sets, Albl[i])]
                end
            end
            Albl
        end
        label_components!(Albl::AbstractArray{Int,$N}, A::AbstractArray, connectivity::BitArray{$N}, bkg = 0) =
            label_components!(Albl, A, convert(Array{Bool}, connectivity), bkg)
    end
end

# Copied directly from DataStructures.jl, but specialized
# to always make the parent be the smallest label
immutable DisjointMinSets
    parents::Vector{Int}

    DisjointMinSets(n::Integer) = new([1:n;])
end
DisjointMinSets() = DisjointMinSets(0)

function find_root!(sets::DisjointMinSets, m::Integer)
    p = sets.parents[m]   # don't use @inbounds here, it might not be safe
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    p
end

# an unsafe variant of the above
function find_root_unsafe!(sets::DisjointMinSets, m::Int)
    @inbounds p = sets.parents[m]
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    p
end

function union!(sets::DisjointMinSets, m::Integer, n::Integer)
    mp = find_root!(sets, m)
    np = find_root!(sets, n)
    if mp < np
        sets.parents[np] = mp
        return mp
    elseif np < mp
        sets.parents[mp] = np
        return np
    end
    mp
end

function push!(sets::DisjointMinSets)
    m = length(sets.parents) + 1
    push!(sets.parents, m)
    m
end

function minlabel(sets::DisjointMinSets)
    out = Array(Int, length(sets.parents))
    k = 0
    for i = 1:length(sets.parents)
        if sets.parents[i] == i
            k += 1
        end
        out[i] = k
    end
    out
end

"`component_boxes(labeled_array)` -> an array of bounding boxes for each label, including the background label 0"
function component_boxes(img::AbstractArray{Int})
    nd = ndims(img)
    n = [Vector{Int64}[ fill(typemax(Int64),nd), fill(typemin(Int64),nd) ]
            for i=0:maximum(img)]
    s = size(img)
    for i=1:length(img)
        vcur = ind2sub(s,i)
        vmin = n[img[i]+1][1]
        vmax = n[img[i]+1][2]
        for d=1:nd
            vmin[d] = min(vmin[d], vcur[d])
            vmax[d] = max(vmax[d], vcur[d])
        end
    end
    map(x->map(y->tuple(y...),x),n)
end

"`component_lengths(labeled_array)` -> an array of areas (2D), volumes (3D), etc. for each label, including the background label 0"
function component_lengths(img::AbstractArray{Int})
    n = zeros(Int64,maximum(img)+1)
    for i=1:length(img)
        n[img[i]+1]+=1
    end
    n
end

"`component_indices(labeled_array)` -> an array of pixels for each label, including the background label 0"
function component_indices(img::AbstractArray{Int})
    n = [Int64[] for i=0:maximum(img)]
    for i=1:length(img)
      push!(n[img[i]+1],i)
    end
    n
end

"`component_subscripts(labeled_array)` -> an array of pixels for each label, including the background label 0"
function component_subscripts(img::AbstractArray{Int})
    n = [Tuple[] for i=0:maximum(img)]
    s = size(img)
    for i=1:length(img)
      push!(n[img[i]+1],ind2sub(s,i))
    end
    n
end

"`component_centroids(labeled_array)` -> an array of centroids for each label, including the background label 0"
function component_centroids(img::AbstractArray{Int})
    nd = ndims(img)
    n = [zeros(nd+1) for i=0:maximum(img)]
    s = size(img)
    for i=1:length(img)
        vcur = ind2sub(s,i)
        vsum = n[img[i]+1]
        for d=1:nd
            vsum[d] += vcur[d]
        end
        vsum[end] += 1
    end
    map(x->tuple((x[1:nd]./x[nd+1])...),n)
end


text:
include("./SparseBitmappedTrieTest.jl")
include("./PersistentVectorTest.jl")
include("./PersistentMapTest.jl")
include("./PersistentSetTest.jl")
include("./PersistentListTest.jl")
include("./PersistentQueueTest.jl")
include("./PersistentDisjointSetTest.jl")
include("./PersistentMacroTest.jl")


text:
function orthonormalize!( LF::LF3dGrid, v )
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
    ortho_gram_schmidt!( v )
    v[:,:] = (1/sqrt(ΔV))*v  # the [:,:] is significant
end

function orthonormalize( LF::LF3dGrid, v )
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
    return ortho_gram_schmidt( v )/sqrt(ΔV)
end


text:
using DelimitedFiles, Unicode

"""

# read_diagnose

Read diagnose output file created by cmdstan.

### Method
```julia
read_diagnose(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_diagnose(model::Stanmodel)

  ## Collect the results of a chain in an array ##

  cnames = String[]
  res_type = "diagnose"
  tdict = Dict()
  local sstr

  for i in 1:model.nchains
    file_path = joinpath(model.tmpdir, "$(model.name)_$(res_type)_$(i).csv")
    if isfile(file_path)

      ## A result type file for chain i is present ##

      if i == 1

        # Extract cmdstan version

        str = read(file_path, String)
        sstr = split(str)
        tdict[:stan_version] = "$(parse(Int, sstr[4])).$(parse(Int, sstr[8])).$(parse(Int, sstr[12]))"
      end

      # Position sstr at element with with `probability=...`

      indx = findall(x -> length(x)>11 && x[1:11] == "probability", sstr)[1]
      sstr_lp = sstr[indx]
      sstr_lp = parse(Float64, split(sstr_lp, '=')[2])

      if :lp in keys(tdict)
        append!(tdict[:lp], sstr_lp)
        append!(tdict[:var_id], parse(Int, sstr[indx+11]))
        append!(tdict[:value], parse(Float64, sstr[indx+12]))
        append!(tdict[:model], parse(Float64, sstr[indx+13]))
        append!(tdict[:finite_dif], parse(Float64, sstr[indx+14]))
        append!(tdict[:error], parse(Float64, sstr[indx+15]))
      else

        # First time around, create value array

        tdict[:lp] = [sstr_lp]
        tdict[:var_id] = [parse(Int, sstr[indx+11])]
        tdict[:value] = [parse(Float64, sstr[indx+12])]
        tdict[:model] = [parse(Float64, sstr[indx+13])]
        tdict[:finite_dif] = [parse(Float64, sstr[indx+14])]
        tdict[:error] = [parse(Float64, sstr[indx+15])]
      end
    end
  end

  (tdict, cnames)

end


text:
using Pkg
Pkg.add("CSV")
Pkg.add("DelimitedFiles")
Pkg.add("BenchmarkTools")
Pkg.add("DoubleFloats")
Pkg.add("PolynomialRoots")
Pkg.add("Distributed")
Pkg.add("ProgressMeter")
Pkg.add("GRUtils")

max_stars_repo_name: georgeglidden/turnip

max_issues_repo_name: georgeglidden/turnip

max_forks_repo_name: georgeglidden/turnip


text:
using robustml
using Test

#ftest(3, 4)
#ftest2(3, 4)

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

max_stars_repo_name: jstriaukas/robustml

max_issues_repo_name: jstriaukas/robustml

max_forks_repo_name: jstriaukas/robustml


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Qiskit

using PyCall

using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Gates
using QuantumCircuits.QCircuits.OtherGates
using QuantumCircuits.QCircuits.Instructions
using QuantumCircuits.QCircuits.Registers

import QuantumCircuits.QCircuits.QBase: add!, tomatrix, measure!
import Base.show

export QiskitSimulator, QiskitCircuit, getQRegister

const qiskit = PyNULL() # pyimport("qiskit")
const plt = PyNULL()
IS_QISKIT = false

function __init__()
    try
        copy!(qiskit, pyimport("qiskit"))
        copy!(plt, pyimport("matplotlib.pyplot"))

        # imports
        pyimport("qiskit.ignis")
        pyimport("qiskit.ignis.verification")

        global IS_QISKIT = true
    catch e
        println(e)
        println("ERROR: QuantumCircuits need qiskit, qiskit.ignis and matplotlib python package.")
    end
end

"The error throw when there is no qiskit instaled."
struct NoQiskitError <:Exception
    msg::String
end

"Quantum device provided by Qiskit."
abstract type QiskitDevice <: QuantumDevice end

"Nothing function"
const nop = () -> nothing

"Quantum circuit provided by Qiskit."
mutable struct QiskitCircuit <: QuantumCircuit
    qc::PyObject
    qubits::Dict{Qubit, PyObject}
    cbits::Dict{Cbit, PyObject}
    qRegs::Dict{String, PyObject}
    cRegs::Dict{String, PyObject}
    draw::Function

    function QiskitCircuit(qRegs::Vector{QuantumRegister},
                           cRegs::Union{Vector{ClassicalRegister}, Nothing} = nothing)
        qregs = [qiskit.QuantumRegister(length(r.bits), r.name) for r in qRegs]
        qubits = Dict{Qubit, PyObject}()
        for (qr, qqr) in zip(qRegs, qregs)
            for (k, v) in zip(qr, qqr)
                push!(qubits, k => v)
            end
        end
        qRegs = Dict{String, PyObject}()
        for r in qregs
            push!(qRegs, r.name => r)
        end

        if cRegs != nothing
            cregs = [qiskit.ClassicalRegister(length(r.bits), r.name) for r in cRegs]

            cbits = Dict{Cbit, PyObject}()
            for (cr, ccr) in zip(cRegs, cregs)
                for (k, v) in zip(cr, ccr)
                    push!(cbits, k => v)
                end
            end

            cRegs = Dict{String, PyObject}()
            for r in cregs
                push!(cRegs, r.name => r)
            end

            # , name="A"
            qc = new(qiskit.QuantumCircuit(qregs..., cregs...), qubits, cbits, qRegs, cRegs, nop)
        else
            cRegs = Dict{String, PyObject}()
            cbits = Dict{Cbit, PyObject}()
            qc = new(qiskit.QuantumCircuit(qregs...), qubits, cbits, qRegs, cRegs, nop)
        end

        qc.draw = () -> draw(qc)

        return qc
    end
end

"Draw the circuit."
function draw(qc::QiskitCircuit)
    qc.qc.draw()
    plt.show()
end

"Function return the register"
getQRegister(qc::QiskitCircuit, name::String) = qc.qRegs[name]

"Add quantum gate to Qiskit circuit"
addQiskitCode(qc::QiskitCircuit, gate::X) = qc.qc.x(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Sx) = qc.qc.sx(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Y) = qc.qc.y(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Z) = qc.qc.z(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::H) = qc.qc.h(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::S) = qc.qc.s(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Sd) = qc.qc.sdg(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::T) = qc.qc.t(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Td) = qc.qc.tdg(qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, cx::CX) = qc.qc.cx(qc.qubits[cx.control], qc.qubits[cx.target])
addQiskitCode(qc::QiskitCircuit, gate::U3) = qc.qc.u3(getvalue(gate.θ), getvalue(gate.ϕ), getvalue(gate.λ), qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Rx) = qc.qc.rx(getvalue(gate.θ), qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Ry) = qc.qc.ry(getvalue(gate.θ), qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Rz) = qc.qc.rz(getvalue(gate.θ), qc.qubits[gate.qubit])
addQiskitCode(qc::QiskitCircuit, gate::Rzx) = qc.qc.rzx(getvalue(gate.θ), qc.qubits[gate.qubit1], qc.qubits[gate.qubit2])
addQiskitCode(qc::QiskitCircuit, gate::Barrier) = qc.qc.barrier()
function addQiskitCode(qc::QiskitCircuit, gate::U)
    qc.qc.rz(gate.β, qc.qubits[gate.qubit])
    qc.qc.ry(gate.γ, qc.qubits[gate.qubit])
    qc.qc.rz(gate.δ, qc.qubits[gate.qubit])
end

"Add gate to circuit"
function add!(qc::QiskitCircuit, gate::QuantumGate)
    addQiskitCode(qc, gate)
    nothing
end

function measure!(qc::QiskitCircuit, qubit::Qubit, cbit::Cbit)
    qc.qc.measure(qc.qubits[qubit], qc.cbits[cbit])
    nothing
end

"Qiskit simulator"
struct QiskitSimulator <: QiskitDevice
end

"Show method"
function Base.show(io::IO, qc::QiskitCircuit)
    s = string(qc.qc.draw(output="text"))
    print(io, s[10:end])
end

"Get the matrix version of the cirquit."
function tomatrix(qc::QiskitCircuit)
    backend = qiskit.BasicAer.get_backend("unitary_simulator")

    return qiskit.execute(qc.qc, backend).result().get_unitary()
end

end  # module Qiskit


text:
module AR1

export AR1Model, residuals, simulate

@doc """
Construct an AR(1) model with mean `mu`, standard deviation (at zero
lag) `sigma`, correlation time `tau`, and measurement uncertainty scale
`nu`.
""" ->
immutable AR1Model
    mu::Float64
    sigma::Float64
    tau::Float64
    nu::Float64
end

""" Returns the uncorrelated residuals and associated uncertainties
for the AR(1) model applied to the given data. """
function residuals(model::AR1Model, ts::Array{Float64, 1}, ys::Array{Float64, 1}, dys::Array{Float64, 1})
    n = size(ts, 1)

    s2 = model.sigma*model.sigma
    nu2 = model.nu*model.nu

    rs = zeros(n)
    ss = zeros(n)

    t = ts[1]
    x = 0.0
    dx2 = s2
    k = 1.0

    for i in 1:n
        # Advance the filter
        dt = ts[i] - t
        t = ts[i]
        lam = exp(-dt/model.tau)
        x = lam*x
        dx2 = lam*lam*(dx2 - s2) + s2

        # Predict the observation
        rs[i] = ys[i] - x - model.mu
        ss2 = dx2 + nu2*dys[i]*dys[i]
        ss[i] = sqrt(ss2)

        # Update gain
        k = dx2 / ss2

        # Update state
        x += rs[i]*k
        dx2 -= ss2*k*k
    end

    rs, ss
end

""" Returns a sample from the given AR(1) model at times `ts` with
observational uncertainties `dxs` """
function simulate(model::AR1Model, ts::Array{Float64, 1}, dxs::Array{Float64, 1})
    xs = zeros(ts)
    xsout = zeros(ts)

    xs[1] = model.sigma*randn()
    xsout[1] = xs[1] + model.nu*dxs[1]*randn() + model.mu

    for i in 2:size(xs,1)
        dt = ts[i] - ts[i-1]
        alpha = exp(-dt/model.tau)

        beta = model.sigma*sqrt(1.0 - alpha*alpha)

        xs[i] = alpha*xs[i-1] + beta*randn()
        xsout[i] = xs[i] + model.nu*dxs[i]*randn() + model.mu
    end

    xsout
end

end

max_stars_repo_name: farr/RVChallenge

max_issues_repo_name: farr/RVChallenge

max_forks_repo_name: farr/RVChallenge


text:
module KiteUtils

#= MIT License

Copyright (c) 2020, 2021 Uwe Fechner

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. =#

# data structures for the flight state and the flight log
# functions for creating a demo flight state, demo flight log, loading and saving flight logs
# function se() for reading the settings
# the parameter P is the number of points of the tether, equal to segments+1
# in addition helper functions for working with rotations

using Rotations, StaticArrays, StructArrays, RecursiveArrayTools, Arrow, YAML, LinearAlgebra, DocStringExtensions, Parameters
export Settings, SysState, SysLog, MyFloat

export demo_state, demo_syslog, demo_log, load_log, save_log, export_log # functions for logging
export demo_state_4p, initial_kite_ref_frame                             # functions for four point kite model
export rot, rot3d, ground_dist, calc_elevation, azimuth_east, acos2      # geometric functions
export set_data_path, load_settings, copy_settings, se                   # functions for reading and copying parameters

"""
    const MyFloat = Float32

Type used for position components and scalar SysState members.
"""
const MyFloat = Float32             # type to use for position components and scalar SysState members
const DATA_PATH = ["data"]          # path for log files and other data

include("settings.jl")
include("transformations.jl")

"""
    SysState{P}

Basic system state. One of these is saved per time step. P is the number
of tether particles.

$(TYPEDFIELDS)
"""
struct SysState{P}
    "time since start of simulation in seconds"
    time::Float64
    "orientation of the kite (quaternion, order w,x,y,z)"
    orient::MVector{4, Float32}
    "elevation angle in radians"
    elevation::MyFloat
    "azimuth angle in radians"
    azimuth::MyFloat
    "tether length [m]"
    l_tether::MyFloat
    "reel out velocity [m/s]"
    v_reelout::MyFloat
    "tether force [N]"
    force::MyFloat
    "depower settings"
    depower::MyFloat
    "norm of apparent wind speed [m/s]"
    v_app::MyFloat
    "vector of particle positions in x"
    X::MVector{P, MyFloat}
    "vector of particle positions in y"
    Y::MVector{P, MyFloat}
    "vector of particle positions in z"
    Z::MVector{P, MyFloat}
end

function Base.getproperty(st::SysState, sym::Symbol)
    if sym == :pos
        X = getfield(st, :X)
        Y = getfield(st, :Y)
        Z = getfield(st, :Z)
        pos = zeros(SVector{length(X), MVector{3, MyFloat}})
        for i in 1:length(X)
            pos[i] .= MVector(X[i], Y[i], Z[i])
        end
        pos
    else
        getfield(st, sym)
    end
end

function Base.getproperty(st::StructVector{SysState}, sym::Symbol)
    if sym == :x
        last.getfield(st, :X)
    elseif sym == :y
        last.getfield(st, :Y)
    elseif sym == :z
        last.getfield(st, :Z) # last.(st.Z)
    else
        getfield(st, sym)
    end
end

function Base.show(io::IO, st::SysState)
    println(io, "time      [s]:       ", st.time)
    println(io, "orient    [w,x,y,z]: ", st.orient)
    println(io, "elevation [rad]:     ", st.elevation)
    println(io, "azimuth   [rad]:     ", st.azimuth)
    println(io, "l_tether  [m]:       ", st.l_tether)
    println(io, "v_reelout [m/s]:     ", st.v_reelout)
    println(io, "force     [N]:       ", st.force)
    println(io, "depower   [-]:       ", st.depower)
    println(io, "v_app     [m/s]:     ", st.v_app)
    println(io, "X         [m]:       ", st.X)
    println(io, "Y         [m]:       ", st.Y)
    println(io, "Z         [m]:       ", st.Z)
end

"""
    SysLog{P}

Flight log, containing the basic data as struct of vectors which can be accessed as if it would
be an array structs.
In addition an extended view on the data that includes derived/ calculated values for plotting.
Finally it contains meta data like the name of the log file.

$(TYPEDFIELDS)
"""
struct SysLog{P}
    "name of the flight log"
    name::String
    "struct of vectors that can also be accessed like a vector of structs"
    syslog::StructArray{SysState{P}}
end

function Base.getproperty(log::SysLog, sym::Symbol)
    if sym == :x
        last.(getproperty(log.syslog, :X))
    elseif sym == :y
        last.(getproperty(log.syslog, :Y))
    elseif sym == :z
        last.(getproperty(log.syslog, :Z))
    else
        getfield(log, sym)
    end
end

# functions
function __init__()
    SETTINGS.segments=0 # force loading of settings.yaml
end

"""
    demo_state(P, height=6.0, time=0.0)

Create a demo state with a given height and time. P is the number of tether particles.

Returns a SysState instance.
"""
function demo_state(P, height=6.0, time=0.0)
    a = 10
    X = range(0, stop=10, length=P)
    Y = zeros(length(X))
    Z = (a .* cosh.(X./a) .- a) * height/ 5.430806
    r_xyz = RotXYZ(pi/2, -pi/2, 0)
    q = QuatRotation(r_xyz)
    orient = MVector{4, Float32}(Rotations.params(q))
    elevation = calc_elevation([X[end], 0.0, Z[end]])
    return SysState{P}(time, orient, elevation,0.,0.,0.,0.,0.,0.,X, Y, Z)
end

"""
    initial_kite_ref_frame(vec_c, v_app)

Calculate the initial orientation of the kite based on the last tether segment and
the apparent wind speed.

Parameters:
- `vec_c`: (pos_n-2) - (pos_n-1) n: number of particles without the three kite particles
                                    that do not belong to the main thether (P1, P2 and P3).
- v_app: vector of the apparent wind speed

Returns:
x, y, z:  the unit vectors of the kite reference frame in the ENU reference frame
"""
function initial_kite_ref_frame(vec_c, v_app)
    z = normalize(vec_c)
    y = normalize(cross(v_app, vec_c))
    x = normalize(cross(y, vec_c))
    return (x, y, z)
end

"""
    get_particles(height_k, height_b, width, m_k, pos_pod= [ 75., 0., 129.90381057], vec_c=[-15., 0., -25.98076211], v_app=[10.4855, 0, -3.08324])

Calculate the initial positions of the particels representing
a 4-point kite, connected to a kite control unit (KCU).

Parameters:
- height_k: height of the kite itself, not above ground [m]
- height_b: height of the bridle [m]
- width: width of the kite [m]
- mk: relative nose distance
- pos_pod: position of the control pod
- vec_c: vector of the last tether segment
"""
function get_particles(height_k, height_b, width, m_k, pos_pod= [ 75., 0., 129.90381057], vec_c=[-15., 0., -25.98076211], v_app=[10.4855, 0, -3.08324])
    # inclination angle of the kite; beta = atan(-pos_kite[2], pos_kite[1]) ???
    beta = pi/2.0
    x, y, z = initial_kite_ref_frame(vec_c, v_app)

    h_kx = height_k * cos(beta); # print 'h_kx: ', h_kx
    h_kz = height_k * sin(beta); # print 'h_kz: ', h_kz
    h_bx = height_b * cos(beta)
    h_bz = height_b * sin(beta)
    pos_kite = pos_pod - (h_kz + h_bz) * z + (h_kx + h_bx) * x   # top,        poing B in diagram
    pos_C = pos_kite + h_kz * z + 0.5 * width * y + h_kx * x     # side point, point C in diagram
    pos_A = pos_kite + h_kz * z + (h_kx + width * m_k) * x       # nose,       point A in diagram
    pos_D = pos_kite + h_kz * z - 0.5 * width * y + h_kx * x     # side point, point D in diagram
    pos0 = pos_kite + (h_kz + h_bz) * z + (h_kx + h_bx) * x      # equal to pos_pod, P_KCU in diagram
    [zeros(3), pos0, pos_A, pos_kite, pos_C, pos_D] # 0, p7, p8, p9, p10, p11
end

"""
    demo_state_4p(P, height=6.0, time=0.0)

Create a demo state, using the 4 point kite model with a given height and time. P is the number of tether particles.

Returns a SysState instance.
"""
function demo_state_4p(P, height=6.0, time=0.0)
    a = 10
    X = collect(range(0, stop=10, length=P))
    Y = zeros(length(X))
    Z = (a .* cosh.(X./a) .- a) * height/ 5.430806
    # append the kite particles to X, Y and z
    pod_pos = [X[end], Y[end], Z[end]]
    particles = get_particles(se().height_k, se().h_bridle, se().width, se().m_k, pod_pos)[3:end]
    for i in 1:4
        particle=particles[i]
        x, y, z = particle[1], particle[2], particle[3]
        push!(X, x)
        push!(Y, y)
        push!(Z, z)
    end
    r_xyz = RotXYZ(pi/2, -pi/2, 0)
    q = QuatRotation(r_xyz)
    orient = MVector{4, Float32}(Rotations.params(q))
    elevation = calc_elevation([X[end], 0.0, Z[end]])
    return SysState{P+4}(time, orient, elevation,0.,0.,0.,0.,0.,0.,X, Y, Z)
end

"""
    demo_syslog(P, name="Test flight"; duration=10)

Create a demo flight log  with given name [String] and duration [s] as StructArray. P is the number of tether
particles.
"""
function demo_syslog(P, name="Test flight"; duration=10)
    max_height = 6.03
    steps   = Int(duration * se().sample_freq) + 1
    time_vec = Vector{Float64}(undef, steps)
    myzeros = zeros(MyFloat, steps)
    elevation = Vector{Float64}(undef, steps)
    orient_vec = Vector{MVector{4, Float32}}(undef, steps)
    X_vec = Vector{MVector{P, MyFloat}}(undef, steps)
    Y_vec = Vector{MVector{P, MyFloat}}(undef, steps)
    Z_vec = Vector{MVector{P, MyFloat}}(undef, steps)
    for i in range(0, length=steps)
        state = demo_state(P, max_height * i/steps, i/se().sample_freq)
        time_vec[i+1] = state.time
        orient_vec[i+1] = state.orient
        elevation[i+1] = asin(state.Z[end]/state.X[end])
        X_vec[i+1] = state.X
        Y_vec[i+1] = state.Y
        Z_vec[i+1] = state.Z
    end
    return StructArray{SysState{P}}((time_vec, orient_vec, elevation, myzeros,myzeros,myzeros,myzeros,myzeros,myzeros, X_vec, Y_vec, Z_vec))
end

"""
    demo_log(P, name="Test_flight"; duration=10)

Create an artifical SysLog struct for demonstration purposes. P is the number of tether
particles.
"""
function demo_log(P, name="Test_flight"; duration=10)
    syslog = demo_syslog(P, name, duration=duration)
    return SysLog{P}(name, syslog)
end

"""
    save_log(flight_log, compress=true)

Save a fligh log of type SysLog as .arrow file. By default lz4 compression is used,
if you use **false** as second parameter no compression is used.
"""
function save_log(flight_log, compress=true)
    filename = joinpath(DATA_PATH[1], flight_log.name) * ".arrow"
    if compress
        Arrow.write(filename, flight_log.syslog, compress=:lz4)
    else
        Arrow.write(filename, flight_log.syslog)
    end
end

"""
    export_log(flight_log)

Save a fligh log of type SysLog as .csv file.
"""
function export_log(flight_log)
    @eval using CSV
    filename = joinpath(DATA_PATH[1], flight_log.name) * ".csv"
    Base.invokelatest(CSV.write, filename, flight_log.syslog)
end

"""
    load_log(P, filename::String)

Read a log file that was saved as .arrow file.  P is the number of tether
particles.
"""
function load_log(P, filename::String)
    if isnothing(findlast(isequal('.'), filename))
        fullname = joinpath(DATA_PATH[1], filename) * ".arrow"
    else
        fullname = joinpath(DATA_PATH[1], filename)
    end
    table = Arrow.Table(fullname)
    syslog = StructArray{SysState{P}}((table.time, table.orient, table.elevation, table.azimuth, table.l_tether, table.v_reelout, table.force, table.depower, table.v_app, table.X, table.Y, table.Z))
    return SysLog{P}(basename(fullname[1:end-6]), syslog)
end

function test(save=false)
    if save
        log_to_save=demo_log(7)
        save_log(log_to_save)
    end
    return(load_log(7, "Test_flight.arrow"))
end

precompile(load_log, (Int64, String,))

end


text:
using HTTP.Messages: Response

function cors(stream::Stream, next)

    headers = [
        "Access-Control-Allow-Origin" => "*",
        "Access-Control-Allow-Headers" => "*",
        "Access-Control-Allow-Methods" => "*",
    ]

    res = stream.message.response
    # @info "cors middelware" target = stream.message.target method = stream.message.method

    if stream.message.method == "OPTIONS"
        for i in headers
            HTTP.setheader(res, i)
        end
    else
        for i in headers
            HTTP.setheader(res, i)
        end
        next(stream)
    end
end

function combine_middleware(middleware::Vector)
    i = length(middleware)

    if i == 0
        return x -> x
    end

    fns = Function[stream->middleware[i](stream, identity)]
    for i in reverse(1:length(middleware)-1)
        fn = stream -> middleware[i](stream, fns[i+1])
        pushfirst!(fns, fn)
    end
    return fns[1]
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using PyPlot
using DelimitedFiles

if !isdir("figures_summary_channel")
  mkdir("figures_summary_channel")
end

m = 90
n = 180
h = 5.0 # meter
dz = 3.0 # meters
dx = 3.0
nz = Int64(round((m * h) / dz)) + 1
nx = Int64(round((n * h) / dx)) + 1
z_src = (collect(5:10:nz-5) .- 1 ) .* dz .+ dz/2.0
x_src = (5-1)ones(Int64, size(z_src)) .* dx .+ dx/2.0
z_rec = (collect(5:1:nz-5) .- 1) .* dz .+ dz/2.0
x_rec = (nx-5-1) .* ones(Int64, size(z_rec)) .*dx .+ dx/2.0

z_inj = (54-1)*h + h/2.0
x_inj = (18-1)*h + h/2.0
z_prod = (54-1)*h + h/2.0
x_prod = (168-1)*h + h/2.0

rc("axes", titlesize=20)
rc("axes", labelsize=18)
rc("xtick", labelsize=18)
rc("ytick", labelsize=18)
rc("legend", fontsize=20)
# true model
figure()
    K = 20.0 .* ones(m,n) # millidarcy
    ix = 1:n
    y1 = 45. .+ 10. .* sin.(ix./120.0 .* 2.0 .* pi)
    y2 = 55. .+ 10. .* sin.(ix./120.0 .* 2.0 .* pi)
    for j = 1:n
        for i = 1:m
            if (i > y1[j] && i < y2[j])
                K[i, j] = 120;
            end
        end
    end
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">")
scatter(x_prod, z_prod, c="r", marker="<")
savefig("figures_summary_channel/K_true.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

# # init model
# figure()
# K = 20.0 .* ones(m,n)
# imshow(K, extent=[0,n*h,m*h,0]);
# xlabel("Distance (m)")
# ylabel("Depth (m)")
# cb = colorbar()
# clim([20, 120])
# cb.set_label("Permeability (md)")
# shot_inds = collect(1:length(z_src))
# scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
# scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
# scatter(x_inj, z_inj, c="r", marker=">")
# scatter(x_prod, z_prod, c="r", marker="<")
# savefig("figures_summary/K_init.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

m = 45
n = 90
h = 10.0 # meter
dz = 3.0 # meters
dx = 3.0
nz = Int64(round((m * h) / dz)) + 1
nx = Int64(round((n * h) / dx)) + 1
z_src = (collect(5:10:nz-5) .- 1 ) .* dz .+ dz/2.0
x_src = (5-1)ones(Int64, size(z_src)) .* dx .+ dx/2.0
z_rec = (collect(5:1:nz-5) .- 1) .* dz .+ dz/2.0
x_rec = (nx-5-1) .* ones(Int64, size(z_rec)) .*dx .+ dx/2.0

z_inj = (27-1)*h + h/2.0
x_inj = (9-1)*h + h/2.0
z_prod = (27-1)*h + h/2.0
x_prod = (84-1)*h + h/2.0
iter = 100
Prj_names = ["CO2_channel_4590_pgs", "CO2_channel_45_90"]
K_name = "/Stage6/K$iter.txt"
figure()
iPrj = 1
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">")
scatter(x_prod, z_prod, c="r", marker="<")
savefig("figures_summary_channel/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);

K_name = "/K$iter.txt"
figure()
iPrj = 2
K = readdlm(Prj_names[iPrj] * K_name)
imshow(K, extent=[0,n*h,m*h,0]);
xlabel("Distance (m)")
ylabel("Depth (m)")
cb = colorbar()
clim([20, 120])
cb.set_label("Permeability (md)")
shot_inds = collect(1:length(z_src))
scatter(x_src[shot_inds], z_src[shot_inds], c="w", marker="*")
scatter(x_rec, z_rec, s=16.0, c="r", marker="v")
scatter(x_inj, z_inj, c="r", marker=">")
scatter(x_prod, z_prod, c="r", marker="<")
savefig("figures_summary_channel/K_$(Prj_names[iPrj]).pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);