text:
using LinearAlgebra
using Printf
using Random

"""
    ABCDs, umax, S = scaleABCD(ABCD, nlev=2, f=0, xlim=1, ymax=nlev+5, umax, N=100000, N0=10)

Scale the loop filter of a general delta-sigma modulator for dynamic range.

Inputs:
- `ABCD`:   The state-space description of the loop filter.
- `nlev`:   The number of levels in the quantizer.
- `xlim`:   A vector or scalar specifying the limit for each state variable.
- `ymax`:   The stability threshold. Inputs that yield quantizer inputs above `ymax`
            are considered to be beyond the stable range of the modulator.
- `umax`:   The maximum allowable input amplitude. `umax` is calculated if it
            is not supplied.
Outputs:
- `ABCDs`:  The state-space description of the scaled loop filter.
- `S`:      The diagonal scaling matrix, `ABCDs` = [`S`A`S`⁻¹ `S`B; C`S`⁻¹ D]; xs = `S`x;
"""
function scaleABCD(ABCD, nlev=2, f=0, xlim=1, ymax=NaN, umax=NaN, N_sim=100000, N0=10)
    if isnan(ymax)
        ymax = nlev + 5
    end
    order = size(ABCD, 1) - 1
    if size(xlim) == ()
        xlim = xlim * ones(order)   # Convert scalar xlim to a vector
    end

    # Make this function repeatable
    rng = MersenneTwister(0)

    # Envelope for smooth start-up
    raised_cosine = 0.5 * (1 .- cos.(pi/N0 * (0:N0-1)))
    if isnan(umax)
        # Simulate the modulator with DC or sine wave inputs to detect its
        # stable input range.
        # First get a rough estimate of umax.
        ulist = (0.1:0.1:1.0) * (nlev-1)
        umax = nlev - 1
        N = 1000
        u0 = [exp.(2im*pi*f*(-N0:-1)) .* raised_cosine; exp.(2im*pi*f*(0:N-1))]
        u0 = u0 .+ 0.01 * randn(rng, N+N0, 2) * [1, 1im]
        u0 = real(u0)

        for u = ulist
            _, _, _, y = simulateDSM(u*u0, ABCD, nlev)
            if maximum(abs.(y)) > ymax
                umax = u    # umax is the smallest input found which causes 'instability'
                break
            end
        end

        if umax == ulist[1]
            error(@sprintf("Modulator is unstable even with an input amplitude of %.1f.", umax))
        end
    end

    # More detailed simulation
    N = N_sim
    u0 = [exp.(2im*pi*f*(-N0:-1)) .* raised_cosine; exp.(2im*pi*f*(0:N-1))]
    u0 = u0 .+ 0.01 * randn(rng, N+N0, 2) * [1, 1im]
    u0 = real(u0)
    maxima = zeros(order) .- 1
    ulist = range(0.7*umax, umax, length=10)
    for u = ulist
        _, _, xmax, y = simulateDSM(u*u0, ABCD, nlev)
        if maximum(abs.(y)) > ymax
            break
        end
        # We need to do this at least once.
        umax = u    # umax is the largest input which keeps |y| < ymax
        maxima, = findmax([maxima xmax], dims=2)
    end

    # Scale the modulator so that all states are at most xlim.
    scale = maxima[:] ./ xlim
    S = diagm(1 ./ scale)
    Sinv = diagm(scale)
    A, B, C, D = partitionABCD(ABCD)
    ABCDs = [S*A*Sinv S*B; C*Sinv D]

    return ABCDs, umax, S
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
export FDSNevq, FDSNevt

"""
    (H,R) = FDSNevq(ot)

Multi-server query for the events with the closest origin time to `ot`.
Returns an Array{SeisHdr,1} in H with event headers and an Array{SeisSrc,1}
in R in H with corresponding source process info.

Keywords: evw, rad, reg, mag, nev, src, to

See also: SeisIO.KW
"""
function FDSNevq(ot::String;
  evw::Union{Array{Real,1},Array{Float64,1},Array{Int64,1}} = KW.evw,
  rad::Array{Float64,1} = KW.rad,
  reg::Array{Float64,1} = KW.reg,
  mag::Array{Float64,1} = KW.mag,
  nev::Int64 = KW.nev,
  src::String = KW.src,
  to::Int = KW.to,
  v::Int64 = KW.v)

  if isempty(reg) && !isempty(rad)
    if length(rad) == 4
      append!(rad, [-30.0, 700.0])
    end
    search_coords = string( "&latitude=", rad[1], "&longitude=", rad[2],
                            "&minradius=", rad[3], "&maxradius=", rad[4],
                            "&mindepth=", rad[5], "&maxdepth=", rad[6] )
  else
    if isempty(reg)
      reg = Float64[-90.0, 90.0, -180.0, 180.0, -30.0, 700.0]
    elseif length(reg) == 4
      append!(reg, [-30.0, 700.0])
    end
    search_coords = string("&minlat=", reg[1], "&maxlat=", reg[2],
                           "&minlon=", reg[3], "&maxlon=", reg[4],
                           "&mindepth=", reg[5], "&maxdepth=", reg[6])
  end

  # Determine time window
  ot2::Float64 = try
    d2u(DateTime(ot))
  catch
    if length(ot) <= 14
      ot0 = string(ot[1:4],"-",ot[5:6],"-",ot[7:8],"T",ot[9:10],":",ot[11:12])
      if length(ot) > 12
        ot1 = string(ot0, ":", ot[13:14])
      else
        ot1 = string(ot0, ":00")
      end
    end
    d2u(DateTime(ot1))
  end
  d0 = string(u2d(ot2 - abs(evw[1])))
  d1 = string(u2d(ot2 + abs(evw[2])))
  oti = round(Int64, ot2*sμ)

  # multi-server query (most FDSN servers do NOT have an event service)
  if lowercase(src) == "all"
    sources = String["EMSC", "INGV", "IRIS", "LMU", "NCEDC", "NIEP", "ORFEUS", "SCEDC", "USGS", "USP"]
  else
    sources = split(src,",")
  end
  sources = [strip(i) for i in sources]
  EvCat = Array{SeisHdr,1}(undef, 0)
  EvSrc = Array{SeisSrc,1}(undef, 0)
  origin_times = Array{Int64,1}(undef, 0)
  for k in sources
    v > 1 && println(stdout, "Querying ", k)
    url = string(fdsn_uhead(String(k)), "event/1/query?",
                                        "starttime=", d0, "&endtime=", d1,
                                        search_coords,
                                        "&minmag=", mag[1], "&maxmag=", mag[2],
                                        "&format=xml")
    v > 0 && println(stdout, "URL = ", url)
    req_info_str = "\nFDSN event query:"

    (R, parsable) = get_http_req(url, req_info_str, to)
    if parsable
      str_req = String(R)
      v > 1 && println(stdout, "REQUEST BODY:\n", str_req)
      xdoc = parse_string(str_req)
      event_xml!(EvCat, EvSrc, xdoc)
      v > 1 && println(stdout, "CATALOG:\n", EvCat)
    end
  end

  if nev > 0
    # Sort based on earliest origin time
    # sort!(EvCat, by = H -> abs(round(Int64, d2u(getfield(H, :ot))*sμ) - oti))
    for H in EvCat
      # push!(origin_times, round(Int64, d2u(getfield(H, :ot))*sμ))
      push!(origin_times, abs(round(Int64, d2u(getfield(H, :ot))*sμ) - oti))
    end
    # k = sortperm(abs.(origin_times.-oti))
    k = sortperm(origin_times)

    n0 = min(length(EvCat), nev)
    n0 < nev && @warn(string("Catalog only contains ", n0, " events (original request was ", nev,")"))
    return EvCat[k[1:n0]], EvSrc[k[1:n0]]
    # return EvCat[1:n0], EvSrc[1:n0]
  else
    return EvCat, EvSrc
  end
end

"""
    FDSNevt(ot::String, chans::String)

Get trace data for the event closest to origin time `ot` on channels `chans`.
Returns a SeisEvent structure.

Standard keywords: fmt, mag, nd, opts, pha, rad, reg, src, to, v, w

Other keywords:
* len::Real (120.0): desired record length in minutes
* model::String ("iasp91"): Earth velocity model for phase calculations

See also: distaz!, FDSNevq, FDSNsta, SeisKW
"""
function FDSNevt(ot::String, chans::Union{String,Array{String,1},Array{String,2}};
  len::Real             = 120.0,
  evw::Union{Array{Real,1},Array{Float64,1},Array{Int64,1}} = KW.evw,
  fmt::String           = KW.fmt,
  mag::Array{Float64,1} = KW.mag,
  model::String         = "iasp91",
  nd::Real              = KW.nd,
  opts::String          = KW.opts,
  pha::String           = KW.pha,
  rad::Array{Float64,1} = KW.rad,
  reg::Array{Float64,1} = KW.reg,
  src::String           = KW.src,
  to::Int64             = KW.to,
  v::Int64              = KW.v,
  w::Bool               = KW.w
  )

  C = fdsn_chp(chans, v=v)

  # Create header
  v > 0 && println(stdout, now(), ": event query begins.")
  (H,R) = FDSNevq(ot, nev=1,
                  rad=rad,
                  reg=reg,
                  mag=mag,
                  src=src,
                  to=to,
                  evw=evw,
                  v=v
                  )
  H = H[1]
  R = R[1]

  # Create channel data
  v > 0 && println(stdout, now(), ": data query begins.")
  s = H.ot                                      # Start time for FDSNsta is event origin time
  t = u2d(d2u(s) + 60*len)                      # End time is len minutes later
  (d0, d1) = parsetimewin(s,t)
  S = SeisData()
  FDSNget!(S, C, fmt=fmt, nd=nd, rad=rad, reg=reg, s=d0, si=true, src=src, t=d1, to=to, v=v, w=w)

  v > 1 && println(stdout, now(), ": channels initialized.")
  v > 2 && println(stdout, S)

  # Initialize SeisEvent structure
  Ev = SeisEvent(hdr = H, source = R, data = S)
  v > 1 && println(stdout, now(), ": SeisEvent created.")
  v > 1 && println(stdout, S)

  # Update Ev with distance, azimuth
  distaz!(Ev)
  v > 1 && println(stdout, now(), ": Δ,Θ updated.")

  # Add phase arrivals to :data
  v > 0 && println(stdout, now(), ": phase query begins.")
  get_pha!(Ev, pha=pha, model=model, to=to, v=v)
  return Ev
end

ext: jl


text:
export AAFT

"""
    AAFT()

An amplitude-adjusted-fourier-transform (AAFT) surrogate[^Theiler1991].

AAFT surrogates have the same linear correlation, or periodogram, and also
preserves the amplitude distribution of the original data.

AAFT surrogates can be used to test the null hypothesis that the data come from a
monotonic nonlinear transformation of a linear Gaussian process
(also called integrated white noise)[^Theiler1991].

[^Theiler1991]: J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
"""
struct AAFT <: Surrogate end

function surrogenerator(x, method::AAFT, rng = Random.default_rng())
    n = length(x)
    m = mean(x)
    forward = plan_rfft(x)
    inverse = plan_irfft(forward * x, n)
    𝓕 = forward * (x .- m)

    init = (
        x_sorted = sort(x),
        ix = zeros(Int, n),
        inverse = inverse,
        m = m,
        𝓕 = 𝓕,
        r = abs.(𝓕),
        ϕ = angle.(𝓕),
        shuffled𝓕 = similar(𝓕),
        coeffs = zero(𝓕),
        n = n,
    )

    return SurrogateGenerator(method, x, similar(x), init, rng)
end

function (sg::SurrogateGenerator{<:AAFT})()
    s, rng = sg.s, sg.rng

    init_fields = (:x_sorted, :ix, :inverse, :m, :r, :ϕ, :shuffled𝓕, :coeffs, :n)
        x_sorted, ix,  inverse, m, r, ϕ, shuffled𝓕, coeffs, n =
        getfield.(Ref(sg.init), init_fields)

    coeffs .= rand(rng, Uniform(0, 2π), length(shuffled𝓕))
    shuffled𝓕 .= r .* exp.(coeffs .* 1im)
    s .= (inverse * shuffled𝓕) .+ m

    # Rescale back to original values to obtain AAFT surrogate.
    sortperm!(ix, s)
    s[ix] .= x_sorted

    return s
end


text:
function CalculateSampling(in)

	cutoff = 1e-10
	itrace = 1
	wd = zeros(Float32,size(in))
	for itrace = 1 : size(in[:,:],2)
		a = sum(in[:,itrace].*in[:,itrace])
		if (a > cutoff)
			wd[:,itrace] = 1.
		end
	end
	return wd;
end

function CalculateSampling(in,h;cutoff=1e-10)

	itrace = 1
	wd = zeros(Float32,size(in))
	for itrace = 1 : size(in[:,:],2)
		a = sqrt(sum(in[:,itrace].^2))
		if (a > cutoff)
			wd[:,itrace] = 1.
		end
	end
	return wd,h;
end

function CalculateSampling(in::String,wd::String;cutoff=1e-10)
	# calculate sampling operator (1's for live traces, 0's for missing traces)

	@compat parameters = Dict(:cutoff=>cutoff)
	SeisProcess(in,wd,[CalculateSampling],[parameters])

end


text:
using DelimitedFiles
using Statistics

D = Dict()
S = Dict()
for activation in ["tanh", "relu", "leaky_relu", "selu", "elu"]
    D[activation] = zeros(4,3)
    for (k1, nn_depth) in enumerate([3,5,10,20])
        for (k2, nn_width) in enumerate([5,20,40])
            FILEID2 = "nn_width$(nn_width)nn_depth$(nn_depth)sigmaY10000.0dsigma1.0activation$(activation)"
            try
                V = readdlm(FILEID2*".txt")
                D[k1, k2] = mean(V)
                S[k1, k2] = std(V)
                # println(val)
            catch
                println("srun -n 1 -N 1 ./Benchmark_.sh $nn_width $nn_depth 10000.0 1.0 $activation &")
            end
        end
    end
end

using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")

sigmaY = [100.0, 0.001e6, 5000.0, 0.01e6, 50000.0, 0.1e6, 1e6]
dsigma = [0.1, 1.0, 10.0, 100.0]
Mean = zeros(length(sigmaY),length(dsigma))
Std = zeros(length(sigmaY),length(dsigma))
for activation in ["tanh"]#, "relu", "leaky_relu", "selu", "elu"]
    for (k1,σY) in enumerate(sigmaY)
        for (k2,dσ) in enumerate(dsigma)
            FILEID2 = "nn_width20nn_depth3sigmaY$(σY)dsigma$(dσ)activation$(activation).txt"
            # @show isfile(FILEID2)
            try
                V = readdlm(FILEID2)[:]
                Mean[k1,k2] = mean(V)
                Std[k1,k2] = std(V)
            catch
                @show FILEID2
                println("srun -n 1 -N 1 ./Benchmark_.sh 20 3 $(σY) $(dσ) $activation &")
            end
        end
    end
end
close("all")
for i = 1:4
    loglog(sqrt.(sigmaY), Mean[:,i], label="\$d=$(dsigma[i])\$", linewidth=4)
end
xlabel("\$\\sigma_Y\$")
ylabel("Error")
legend()
grid("on")
gca().invert_xaxis()
mpl.save("compare_d.tex")


text:
# TODO(odow):
#
# This script can be used to build the C interface to CPLEX. However, it
# requires you to manually do the following steps first:
#
# 1) Copy cplex.h from CPLEX into this /scripts directory
# 2) Copy cpxconst.h from CPLEX into this /scripts directory
#
# You should also check for an updated version of Clang.jl, and update if one is
# available.

import Clang

const LIBCPX_HEADERS = [
    joinpath(@__DIR__, "cplex.h"),
    joinpath(@__DIR__, "cpxconst.h"),
]

const GEN_DIR = joinpath(dirname(@__DIR__), "src", "gen")

wc = Clang.init(
    headers = LIBCPX_HEADERS,
    output_file = joinpath(GEN_DIR, "libcpx_api.jl"),
    common_file = joinpath(GEN_DIR, "libcpx_common.jl"),
    clang_args = String[
        "-I" * header for header in Clang.find_std_headers()
    ],
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "libcplex",
    clang_diagnostics = true,
)

run(wc)

function manual_corrections_common()
    filename = joinpath(GEN_DIR, "libcpx_common.jl")
    lines = readlines(filename; keep = true)
    for (i, line) in enumerate(lines)
        if occursin("CPXINT_MAX", line)
            lines[i] = replace(line, "= INT_MAX" => "= $(typemax(Cint))")
        elseif occursin("CPXINT_MIN", line)
            lines[i] = replace(line, "= INT_MIN" => "= $(typemin(Cint))")
        elseif occursin("= version", line)
            lines[i] = replace(line, "= version" => "= nothing")
        elseif occursin("= NAN", line)
            lines[i] = replace(line, "= NAN" => "= NaN")
        elseif occursin("# Skipping Typedef: CXType_FunctionProto ", line)
            lines[i] = replace(
                replace(
                    line,
                    "# Skipping Typedef: CXType_FunctionProto " => "const "
                ),
                '\n' => " = Ptr{Cvoid}\n"
            )
        end
    end
    open(filename, "w") do io
        print.(Ref(io), lines)
    end
end
manual_corrections_common()

function manual_corrections_api()
    filename = joinpath(GEN_DIR, "libcpx_api.jl")
    lines = readlines(filename; keep = true)
    for (i, line) in enumerate(lines)
        if occursin("Cstring", line)
            lines[i] = replace(line, "Cstring" => "Ptr{Cchar}")
        end
    end
    open(filename, "w") do io
        print.(Ref(io), lines)
    end
end
manual_corrections_api()

rm(joinpath(GEN_DIR, "LibTemplate.jl"))


text:
# ---
# title: 1332. Remove Palindromic Subsequences
# id: problem1332
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/remove-palindromic-subsequences/description/>
# hidden: true
# ---
#
# Given a string `s` consisting only of letters `'a'` and `'b'`. In a single
# step you can remove one palindromic **subsequence**  from `s`.
#
# Return the minimum number of steps to make the given string empty.
#
# A string is a subsequence of a given string, if it is generated by deleting
# some characters of a given string without changing its order.
#
# A string is called palindrome if is one that reads the same backward as well
# as forward.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "ababa"
#     Output: 1
#     Explanation: String is already palindrome
#
#
# **Example 2:**
#
#
#
#     Input: s = "abb"
#     Output: 2
#     Explanation: " **a** bb" -> " **bb** " -> "".
#     Remove palindromic subsequence "a" then "bb".
#
#
# **Example 3:**
#
#
#
#     Input: s = "baabb"
#     Output: 2
#     Explanation: " **baa** b **b** " -> "b" -> "".
#     Remove palindromic subsequence "baab" then "b".
#
#
# **Example 4:**
#
#
#
#     Input: s = ""
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 1000`
#   * `s` only consists of letters 'a' and 'b'
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Copyright 2013 Peter Chronz
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

macro length(typ)
    # Check whether this type contains any fields whose size cannot be
    # inferred right now. Assuming we can only get the size of Numbers (is this
    # true?) reliably.
    tayp::Type = eval(typ)
    typs = tayp.types
    dirty = false
    for t in typs
        if !(t<:Number)
            dirty = true
            break
        end
    end
    if !dirty
        # Return both a generic version and a specific one that delegates to the
        # generic one.
        quote
            function $(esc(:give_length))(::Type{$(typ)})
                # TODO XXX Is there a 1-to-1 correspondence between the literal
                # type definition of fields, what names returns and what types
                # contains? The ordering in specific is of importance!
                typs = $(tayp).types
                len = 0
                for t in typs
                    len += sizeof(t)
                end
                len
            end
            function $(esc(:give_length))(obj::$(typ))
                $(esc(:give_length))($(typ))
            end
        end
    else
        # Expression to access a certain field of the object. This is spliced in
        # below.
        fields = names(tayp)
        typs = tayp.types
        exs::Array{Expr} = Array(Expr, length(fields))
        for i = 1:Base.length(typs)
            accex = Expr(:quote, :($(fields[i])))
            exs[i] = typs[i] <: Number ? :(len += sizeof(obj.($(accex)))) : :(len += give_length(obj.($(accex))))
        end
        # Return only the specific one.
        quote
            function $(esc(:give_length))(obj::$(typ))
                len = 0
                $(exs...)
                len
            end
        end
    end
end
# Create a function that will convert your given type to bytes.
macro bytes(typ)
    t::DataType = eval(typ)
    fields = names(t)
    typs = t.types
    exs = Array(Expr, length(fields))
    for i = 1:length(fields)
        ex = Expr(:., :obj, Expr(:quote, fields[i]))
        exs[i] = quote
            bs = bytes($ex)
            bs_len = length(bs)
            byts[pos:pos + bs_len - 1] = bs
            pos += bs_len
        end
    end
    quote
        function $(esc(:bytes))(obj::$(typ))
            objsize = give_length(obj)
            byts::Bytes = Array(Uint8, objsize)
            bs::Bytes
            pos::Integer = 1
            $(exs...)
            byts
        end
    end
end
# A @string macro; probably quite similar to the bytes macro
macro string(typ, args...)
    t::DataType = eval(typ)
    fields::Vector{Symbol} = names(t)
    strfields::Vector{Symbol}
    if length(args) > 0
        strfields = eval(args[1])
    else
        strfields = Symbol[]
    end
    typs::Tuple = t.types
    exs::Vector{Expr} = Array(Expr, 0)
    push!(exs, :(strs[1] = string("\n$(repeat("\t", tabs))Type: ", $(typ), "\n")))
    push!(exs, :(pos::Int = 2))
    for i = 2:length(typs)
        objacc = Expr(:., :obj, Expr(:quote, fields[i]))
        if ismatch(r"^pad.*", string(fields[i]))
            continue
        elseif typs[i] <: OfpStruct || eltype(typs[i]) <: OfpStruct
            push!(exs, :(strs[pos] = string($(Expr(:quote, fields[i])), " = ",
                string($(objacc), tabs + 1), "\n")))
            push!(exs, :(pos += 1))
        elseif typs[i] <: Bytes && contains(strfields, fields[i])
            push!(exs, :(strs[pos] = string(repeat("\t", tabs), $(Expr(:quote,
                fields[i])), " = ", ascii($(objacc)), "\n")))
            push!(exs, :(pos += 1))
        else
            push!(exs, :(strs[pos] = string(repeat("\t", tabs), $(Expr(:quote,
                fields[i])), " = ", string($(objacc)), "\n")))
            push!(exs, :(pos += 1))
        end
    end
    quote
        import Base.string
        $(esc(:string))(obj::$(typ)) = $(esc(:string))(obj, 0)
        function $(esc(:string))(obj::$(typ), tabs::Int)
            strs::Vector{String} = Array(String, length($(typ).types))
            $(exs...)
            "$(strs[1:(pos - 1)]...)"
        end
    end
end
# Create a byte-based constructor for the given type.
function bytesconstructor(t::Symbol, typ::Type, fields::Vector{Symbol},
    typs::Tuple, sizes::Dict{Symbol, Int}, varsizes::Dict{Symbol, Expr})
    # Expanding the array dynamically since we need to filter out the padding.
    # The penalty should be quite low, since there are rather few fields in each
    # type, and the construction macro is going to be called just once per type.
    exs::Vector{Expr} = Expr[]
    hasheader::Bool = false
    for i = 1:length(fields)
        # Padding fields are set by the constructors. We do not have any access
        # to those fields from here.
        # The OfpHeader needs to be created by the inner constructor implicitly.
        ex::Expr
        if ismatch(r"^pad.*", string(fields[i]))
            if typs[i] <: Number
                ex = :(pos += $(sizeof(typs[i])))
            elseif typs[i] <: Array
                ex = :(pos += $(sizes[fields[i]]))
            else
                error("Unrecognized type for padding.")
            end
        elseif typs[i] == Uint8
            ex = quote
                val = bytes[pos]
                local $(fields[i]) = val
                field_len += 1
                fieldvals[field_len] = val
                pos += 1
            end
        elseif typs[i] == Uint16
            ex = quote
                val = btoui(bytes[pos], bytes[pos + 1])
                local $(fields[i]) = val
                field_len += 1
                fieldvals[field_len] = val
                pos += 2
            end
        elseif typs[i] == Uint32
            ex = quote
                val = btoui(bytes[pos], bytes[pos + 1], bytes[pos + 2], bytes[pos + 3])
                local $(fields[i]) = val
                field_len += 1
                fieldvals[field_len] = val
                pos += 4
            end
        elseif typs[i] == Uint64
            ex = quote
                val = btoui(bytes[pos], bytes[pos + 1], bytes[pos + 2], bytes[pos + 3], bytes[pos + 4], bytes[pos + 5], bytes[pos + 6], bytes[pos + 7])
                local $(fields[i]) = val
                field_len += 1
                fieldvals[field_len] = val
                pos += 8
            end
        elseif typs[i] == OfpHeader
            hasheader = true
            ex = :((field_len += 1; fieldvals[field_len] = header))
        elseif typs[i] <: OfpStruct
            ex = quote
                o = $(typs[i])(bytes[pos:end])
                isa(o, $(typs[i]))
                field_len += 1
                fieldvals[field_len] = o
                pos += give_length(o)
            end
        elseif typs[i] <: Bytes
            lenex::Expr = if haskey(sizes, fields[i])
                :(pos + $(sizes[fields[i]]) - 1)
            elseif haskey(varsizes, fields[i])
                :(pos + $(varsizes[fields[i]]) - 1)
            elseif i == length(fields)
                :(length(bytes))
            else
                error("Cannot process a byte array of unknown length if it is not the last field!")
            end
            ex = quote
                field_len += 1
                fieldvals[field_len] = bytes[pos:$(lenex)]
                pos = $(lenex) + 1
            end
        elseif typs[i] <: Array
            eltyp = eltype(typs[i])
            if !(eltype(typs[i]) <: OfpStruct)
                error("Handling array fields works only with element types that
                    subtype OfpStruct and not with $(typs[i])!")
            end
            ex = quote
                arr = Array($(eltyp), 0)
                # XXX If I change the following 4 lines to a conditional
                # statement, where "constructor" is declared and defined within
                # the body of the if-statement, then I get a strange looking
                # error which tells me that constructor is not defined. Bug in
                # Julia?
                constructor = $(eltyp)
                if !isleaftype($(eltyp))
                    constructor = $(symbol(string(eltyp, "Factory")))
                end
                while pos <= Base.length(bytes)
                    o = constructor(bytes[pos:end])
                    arr = [arr..., o]
                    len = give_length(o)
                    pos += len
                end
                field_len += 1
                fieldvals[field_len] = arr
            end
        end
        exs = [exs..., ex]
    end
    args::Vector{Expr} = hasheader ? [:(header::OfpHeader), :(bytes::Bytes)] :
        [:(bytes::Bytes)]
    quote
        function $(esc(t))($(args...))
            # Values of the fields to be passed in to the inner constructor
            fieldvals = cell(16)
            field_len = 0
            pos = 1
            $(exs...)
            $(t)(fieldvals[1:field_len]...)
        end
    end
end
# szs[1]-->Dict{Symbol, Int}: length of a field
# szs[2]-->Dict{Symbol, Expr}: expression to compute the length of a field
macro bytesconstructor(t, szs...)
    typ::Type = eval(t)
    fields::Vector{Symbol} = names(typ)
    typs::Tuple = typ.types
    sizes::Dict{Symbol, Int} = length(szs) > 0 ? (length(eval(szs[1])) > 0 ?
        eval(szs[1]) : Dict{Symbol, Int}()) : Dict{Symbol, Int}()
    varsizes::Dict{Symbol, Expr} = length(szs) > 1 ? (length(eval(szs[2])) > 0 ?
        eval(szs[2]) : Dict{Symbol, Expr}()) : Dict{Symbol, Expr}()
    bytesconstructor(t, typ, fields, typs, sizes, varsizes)
end

max_stars_repo_name: pchronz/juliaflow

max_issues_repo_name: pchronz/juliaflow

max_forks_repo_name: pchronz/juliaflow


text:
using EzXML
using LaserFileSplitter

svgfile = "data/exp_raw/DeathMoonDevil.svg"

begin
    svgdoc = readxml(svgfile);
    svgfile
end

begin
    elms = collect(eachelement(root(svgdoc)))
    svgfile
end

exnode = elms[1000]

exnode["class"]

elements_with_class = [node for node in elms if haskey(node, "class")]

classes = Set(node["class"] for node in elements_with_class)
classes_nodedict = Dict{String, Vector{EzXML.Node}}(class=>[] for class in classes)

for elm in elements_with_class
    push!(classes_nodedict[elm["class"]], elm)
end

class_lengths = length.(values(classes_nodedict))
sorted_class_lengths = sort(class_lengths)

exnode["d"]

#named_elms = filter(n->haskey(n, "name"), elms)
paths = [node for node in elms if node.name == "path"]
paths = filter(n->n.name == "path", elms)

#write(datadir("processed/DeathMoonDevil_same.svg"), svgdoc)

EzXML.parentelement(elms[100])

parents = Set(EzXML.parentelement.(elms))

EzXML.nodename(elms[100])

collect(eachattribute(root(svgdoc)))
nodetype.(elms)

k = 1
splitby = 2000

for i in 2001:length(paths)
    unlink!(paths[i])
end
#for i in 1:2000
    #unlink!(paths[i])
#end

#for i

num_per_chunk = 1000

filename = svgfile
outfolder = "data/processed/chunked_deathmoondevil2"
max_num_chunks = 10
save_svg_in_chunks(filename, outfolder, max_num_chunks)

newelms = collect(eachelement(root(svgdoc)))

write(datadir("processed/DeathMoonDevil_fewpaths_0.svg"), svgdoc)

for i in 2001:length(paths)
    link!(root(svgdoc), paths[i])
end

write(datadir("processed/DeathMoonDevil_shouldbesame.svg"), svgdoc)

max_stars_repo_name: zoemcc/LaserFileSplitter

max_issues_repo_name: zoemcc/LaserFileSplitter

max_forks_repo_name: zoemcc/LaserFileSplitter


text:
# This is code for computing a transit model integrated over
# a time step.

# First the version without derivatives:
function integrate_timestep(param::Array{T,1},trans::Transit_Struct{T},time::T,dt::T,tol::T,maxdepth::Int64) where {T <: Real}

  function solver(param::Array{T,1},time::T) where {T <: Real}
  # This predicts the impact parameter as a function of time in
  # terms of the transit_struct.
  # For now, transit is approximated as a straight line, where params are given as:
  # param[1] = t0 (mid-transit time; units of days/JD)
  # param[2] = v/R_* (sky velocity in terms of radius of the star - units are day^{-1}
  # param[3] = r (radius-ratio)
  # param[4] = b (impact parameter at time t0)
  btime = sqrt(param[4]^2+(param[2]*(time-param[1]))^2)
  return
  end

  # Now, adopt Dan Foreman-Mackey's recursive integration approach:
  function integrate(param::Array{T,1}, f1::T, f2::T, t1::T, t2::T, trans::Transit_Struct{T}) where {T <: Real}
  tmid = 0.5*(t1+t2);
  bmid = solver(param,tmid)
  trans.b = bmid
  fmid = transit_poly!(trans)
  fapprox = 0.5*(f1+f2)
  d = abs(fmid-fapprox)
  if d > tol && depth < maxdepth
     a = integrate(param,f1,fmid,t1,tmid,depth+1)
     b = integrate(param,fmid,f2,tmid,t2,depth+1)
     return a+b
  end
  return fapprox*(t2-t1)
  end

  t1 = time - 0.5*dt
  b1 = solver(param,t1)
  trans.b = b1
  f1 = transit_poly!(trans)
  t2 = time + 0.5*dt
  b2 = solver(param,t2)
  trans.b = b2
  f2 = transit_poly!(trans)
  fint = integrate(param,f1,f2,t1,t2,trans)
return fint
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# 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.
using Test
using Tesseract

# -------------------------------------------------------------------------------------------------
# Test tess_init()
@testset "tess_init" begin
    local inst = TessInst("eng", datadir)
    @test tess_init(inst, "eng", datadir) == true

    @suppress begin
        local inst = TessInst("eng", datadir)
        @test tess_init(inst, "eng", ".") == false

        local inst = TessInst("eng", datadir)
        @test tess_init(inst, "zzzz", datadir) == false

        local inst = TessInst("eng", datadir)
        @test tess_init(inst, "zzzz", ".") == false
    end

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_init(inst, "eng", datadir)) == false
end

# -------------------------------------------------------------------------------------------------
# Test tess_initialized_languages()
@testset "tess_initialized_languages" begin
    local inst = TessInst("eng", datadir)

    @test tess_initialized_languages(inst) == "eng"

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_initialized_languages(inst)) == nothing
end

# -------------------------------------------------------------------------------------------------
# Test tess_loaded_languages()
@testset "tess_loaded_languages" begin
    local inst = TessInst("eng", datadir)

    @test tess_loaded_languages(inst) == [ "eng" ]

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_loaded_languages(inst)) == nothing
end

# -------------------------------------------------------------------------------------------------
# Test tess_available_languages()
@testset "tess_available_languages" begin
    local inst = TessInst("eng", datadir)

    @test length(tess_available_languages(inst)) > 0

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_available_languages(inst)) == nothing
end

# -------------------------------------------------------------------------------------------------
# Test tess_params()
@testset "tess_params()" begin
    local inst = TessInst("eng", datadir)

    @test isempty(tess_params(inst)) == false

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_params(inst)) == nothing
end

# -------------------------------------------------------------------------------------------------
# Test tess_params(filename)
@testset "tess_params(filename)" begin
    local inst = TessInst("eng", datadir)
    local file = safe_tmp_file()

    @test tess_params(inst, file) == true
    @test filesize(file) > 0
    rm(file)

    local err  = "Instance has been freed."
    local file = safe_tmp_file()
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_params(inst, file)) == false
    @test filesize(file) == 0
    rm(file)
end

# -------------------------------------------------------------------------------------------------
# Test tess_params(IO)
@testset "tess_default_params(IO)" begin
    local inst = TessInst("eng", datadir)

    local buffer = IOBuffer()
    @test tess_params(inst, buffer) == true
    @test length(take!(buffer)) > 0

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    local buffer = IOBuffer()
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_params(inst, buffer)) == false
    @test length(take!(buffer)) == 0
end

# -------------------------------------------------------------------------------------------------
# Test tess_params_parsed()
@testset "tess_params_parsed()" begin
    local inst = TessInst("eng", datadir)

    @test isempty(tess_params_parsed(inst)) == false

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test (@test_logs (:error, err) isempty(tess_params_parsed(inst))) == true
end

# -------------------------------------------------------------------------------------------------
# Test tess_read_config()
@testset "tess_read_config" begin
    local inst = TessInst("eng", datadir)

    local filename, io = mktemp(;cleanup=false)

    try
        println(io, "tessedit_char_blacklist amz")
        println(io, "dawg_debug_level 3")
        close(io)
        @test_nowarn tess_read_config(inst, filename)
    finally
        rm(filename)
    end

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test_logs (:error, err) tess_read_config(inst, filename)
end

# -------------------------------------------------------------------------------------------------
# Test tess_read_debug_config()
@testset "tess_read_debug_config" begin
    local inst = TessInst("eng", datadir)

    local filename, io = mktemp(;cleanup=false)

    try
        println(io, "tessedit_char_blacklist amz")
        println(io, "dawg_debug_level 3")
        close(io)
        @test_nowarn tess_read_debug_config(inst, filename)
    finally
        rm(filename)
    end

    local err  = "Instance has been freed."
    local inst = TessInst("eng", datadir)
    tess_delete!(inst)
    @test_logs (:error, err) tess_read_debug_config(inst, filename)
end


text:
using LinearAlgebra
using LoopVectorization
using BenchmarkTools
using StaticArrays
using RecursiveFactorization
function cholfact!(A)
    # fills in lower triangular portion
    # total flops:
    #   2*n*(n+1)*(n-1)/6 + n + n*(n+1)/2
    # total assignments
    #   n*(n+1)*(n-1)/6 + (n+1)/2
    n = size(A,1)
    for j = 1:n

        # (n+1)*(n-1)/6
        for k = 1:j-1  # loop over columns
            for i = j:n
                @inbounds A[i,j] -= A[i,k] * A[j,k]
            end
        end
        ajj = sqrt(A[j,j])  # 1 flop

        # (n+1)/2 flops
        for i = j:n
            @inbounds A[i,j] /= ajj
        end
    end
    return A
end

function avxchol!(A)
    n = size(A,1)
    for j = 1:n

        # (n+1)*(n-1)/6
        for k = 1:j-1  # loop over columns
            @turbo for i = j:n
                A[i,j] -= A[i,k] * A[j,k]
            end
        end
        ajj = sqrt(A[j,j])  # 1 flop

        # (n+1)/2 flops
        @turbo for i = j:n
            A[i,j] /= ajj
        end
    end
end

function _chol_expr(j,n)
    quote
        for k = 1:$j-1  # loop over columns
            @turbo for i = $j:$n
                A[i,$j] -= A[i,k] * A[$j,k]
            end
        end
        ajj = sqrt(A[$j,$j])

        @turbo for i = $j:$n
            A[i,$j] /= ajj
        end
    end
end

@generated function avxchol_unrolled!(A, ::Val{n}) where n
    expr = [_chol_expr(j,n) for j = 1:n]
    quote
        $(Expr(:block, expr...))
    end
end

n = 50
A = randn(n,n); A = A'A;
C = zero(A)
println("Raw Loop")
@btime begin
    $C .= $A
    cholfact!($C)
end

C .= A
avxchol!(C)
LowerTriangular(C) ≈ cholesky(A).L
println("avxchol")
@btime begin
    $C .= $A
    cholfact!($C)
end

println("LAPACK")
@btime begin
    $C .= $A
    LAPACK.potrf!('L', $C)
end

println("Static Arrays")
SM = similar_type(SMatrix{n,n,Float64})
@btime cholesky($(SM(A)))
@btime begin
    Cs = $SM($A)
    $C .= cholesky(Cs).L
end

p = zeros(Int,n)
lu = copy(A)
RecursiveFactorization.lu!(lu, p)
println("RF lu")
@btime RecursiveFactorization.lu!($lu, $p)


text:
#########################
# All RefTetrahedron   #
#########################
struct RefTetrahedron <: AbstractRefShape end
@inline get_num_faces(::Type{RefTetrahedron}, ::Type{Val{dim}}) where {dim} = dim + 1
@inline get_num_vertices(::Type{RefTetrahedron}, ::Type{Val{dim}}) where {dim} = dim + 1

function reference_coordinates(::Type{RefTetrahedron}, ::Type{Val{1}})
    return [Vec{1, Float64}((0.0,)),Vec{1, Float64}((1.0,))]
end
function gettopology(::Type{RefTetrahedron}, ::Type{Val{1}})
    return Dict(0=>2,1=>1)
end
function reference_coordinates(::Type{RefTetrahedron}, ::Type{Val{2}})
    [Vec{2, Float64}((0.0, 0.0)),
     Vec{2, Float64}((1.0, 0.0)),
     Vec{2, Float64}((0.0, 1.0))]
end
function reference_edges(::Type{RefTetrahedron},::Type{Val{2}})
    [[Vec{2, Float64}((1.0, 0.0)),Vec{2, Float64}((0.0, 1.0))],
     [Vec{2, Float64}((0.0, 1.0)),Vec{2, Float64}((0.0, 0.0))],
     [Vec{2, Float64}((0.0, 0.0)),Vec{2, Float64}((1.0, 0.0))]]
end
@inline reference_edge_nodes(::Type{RefTetrahedron},::Type{Val{2}}) = ((2,3),(3,1),(1,2))

function gettopology(::Type{RefTetrahedron}, ::Type{Val{2}})
    return Dict(0=>3,1=>3,2=>1)
end

"""
get_nodal_points(shp::RefTetrahedron, dim, order)
get points for a nodal basis of order `order` on a `dim`
    dimensional simplex
"""
function get_nodal_points(::Type{RefTetrahedron}, ::Type{Val{1}}, order)
    points = Vector{Vec{1,Float64}}()
    vertices = reference_coordinates(RefTetrahedron, Val{1})
    topology = Dict{Int, Int}()
    append!(points, vertices)
    push!(topology, 0=>length(points))
    append!(points, _interior_points(vertices, order))
    push!(topology, 1=>length(points)-topology[0])
    points, topology
end

function get_nodal_points(::Type{RefTetrahedron}, ::Type{Val{2}}, order)
    points = Vector{Vec{2,Float64}}()
    vertices = reference_coordinates(RefTetrahedron, Val{2})
    topology = Dict{Int, Int}()
    append!(points, vertices)
    push!(topology, 0=>length(points))
    [append!(points, _interior_points(verts, order)) for verts in reference_edges(RefTetrahedron, Val{2})]
    push!(topology, 1=>length(points)-topology[0])
    append!(points, _interior_points(vertices, order))
    push!(topology, 2=>length(points)-topology[0]-topology[1])
    points, topology
end

function _interior_points(verts, order)
    n = length(verts)
    ls = [(verts[i] - verts[1])/order for i in 2:n]
    m = length(ls)
    grid_indices =  []
    if m == 1
        grid_indices = [[i] for i in 1:order-1]
    elseif m == 2 && order > 2
        grid_indices = [[i,j] for i in 1:order-1 for j in 1:order-i-1]
    end
    pts = Vector{typeof(verts[1])}()
    for indices in grid_indices
        res = verts[1]
        for (i,ii) in enumerate(indices)
            res += (ii) * ls[m - i+1]
        end
        push!(pts,res)
    end
    pts
end

function weighted_normal(J::Tensor{2,2}, face::Int, ::Type{RefTetrahedron}, ::Type{Val{2}})
    @inbounds begin
        face == 1 && return Vec{2}((-(J[2,1] - J[2,2]), J[1,1] - J[1,2]))
        face == 2 && return Vec{2}((-J[2,2], J[1,2]))
        face == 3 && return Vec{2}((J[2,1], -J[1,1]))
    end
    throw(ArgumentError("unknown face number: $face"))
end

""" Compute volume of a simplex spanned by vertices `verts` """
function volume(verts::Vector{Vec{N, T}}) where {N,T}
    # Volume of reference simplex element is 1/n!
    n = length(verts) - 1
    ref_verts = reference_coordinates(RefTetrahedron, Val{n})
    A, b = get_affine_map(ref_verts, verts)
    F = svd(A)
    return *(F.S...)/factorial(n)
end


text:
# using Pkg
# Pkg.add(PackageSpec(url="https://github.com/RoboticExplorationLab/MatrixCalculus.jl.git"))


text:
"""
Flatten.jl converts data between nested and flat structures, using `flatten()`,
`reconstruct()` and `update!()` functions. This facilitates building modular,
composable structs while allowing access to solvers and optimisers that require flat
vectors of parameters -- or any other use case that requires extraction or modification
of a list of values from a nested type. Importantly it's type-stable and _very_ fast.

Flatten is also flexible. The types to use and ignore can be specified, and fields can be
ignored using field level traits like `flattenable` from FieldMetadata.jl. Method overrides
can also be defined for custom types.

Flatten allows 'querying' to extract some types and ignore others, here using `flatten`:

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> struct Bar{X}
           x::X
       end

julia> obj = Foo(1, :two, Foo(Bar(3), 4.0, 5.0f0));

julia> use = Union{Int, Float32}; # Immediately return Int and AbstractFloat fields

julia> ignore = Bar;  # Dont return Bar or iterate over Bar fields

julia> flatten(obj, use, ignore) # Flatten all Int and Float32 except fields of Bar
(1, 5.0f0)

julia> flatten(Foo(:one, :two, Foo(Bar(:three), 4.0, :five)), Symbol, Bar) # Return all symbols, except in Bar
(:one, :two, :five)
```

The default type used is `Real`. These rules apply to all Flatten.jl functions.

Flatten.jl also uses [FieldMetadata.jl](https://github.com/rafaqz/FieldMetadata.jl)
to provide a `@flattenable` macro, allowing you to choose struct fields to include
and remove from flattening -- defaulting to `true` for all fields.

```jldoctest
using Flatten
import Flatten: flattenable

@flattenable struct Bar{X,Y}
    x::X | true
    y::Y | false
end

flatten(Bar(1, 2))

# output
(1,)
```
Custom `@metadata` methods from FieldMetadata can be used, if they return a Bool.
You can also use cusom functions that follow the following form, returning a boolean:

```julia
f(::Type, ::Type{Var{:fieldname}}) = false
```

Flatten also provides `metaflatten()` to flatten any FieldMetadata.jl
metadata for the same fields `flatten()` returns. This can be useful for attaching
information like descriptions or prior propability distributions to each field.
Regular field data can also be collected with convenience versions of metaflatten:
`fieldnameflatten`, `parentflatten`, `fieldtypeflatten` and `parenttypeflatten`
functions provide lists of fieldnames and types that may be useful for building
parameter display tables.

This package was started by Robin Deits (@rdeits), and its early development
owes much to discussions and ideas from Jan Weidner (@jw3126) and Robin Deits.
"""
module Flatten

using FieldMetadata, Requires
import FieldMetadata: @flattenable, @reflattenable, flattenable

export @flattenable, @reflattenable, flattenable, flatten, reconstruct, update!, modify,
       metaflatten, fieldnameflatten, parentnameflatten, fieldtypeflatten, parenttypeflatten

# Optionally load Unitful and unlittless falttening
function __init__()
    @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" include("unitful.jl")
end

struct EmptyIgnore end

# Default behaviour when no flattentrait/use/ignore args are given
const USE = Real
const IGNORE = EmptyIgnore
const FLATTENTRAIT = flattenable

# Generalised nested struct walker
nested(T::Type, expr_builder, expr_combiner, action) =
    nested(T, Nothing, expr_builder, expr_combiner, action)
nested(T::Type, P::Type, expr_builder, expr_combiner, action) =
    expr_combiner(T, [Expr(:..., expr_builder(T, fn, action)) for fn in fieldnames(T)])

"""
    constructor_of(::Type)
Add methods to define constructors for types with custom type parameters.
"""
@generated constructor_of(::Type{T}) where T = getfield(T.name.module, Symbol(T.name.name))
constructor_of(::Type{T}) where T<:Tuple = tuple
constructor_of(T::UnionAll) = constructor_of(T.body)

"""
    flatten(obj, [flattentrait::Function], [use::Type], [ignore::Type])

Flattening. Flattens a nested struct or tuple to a flat tuple.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `data`: Replacement data - an `AbstractArray`, `Tuple` or type that defines `getfield`.
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
```julia
f(::Type, ::Type{Val{:fieldname}}) = true
```
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> foo = Foo(1, 2, 3)
Foo{Int64,Int64,Int64}(1, 2, 3)

julia> flatten(foo)
(1, 2, 3)

julia> nested = Foo(Foo(1, 2, 3), 4.0, 5.0)
Foo{Foo{Int64,Int64,Int64},Float64,Float64}(Foo{Int64,Int64,Int64}(1, 2, 3), 4.0, 5.0)

julia> flatten(nested)
(1, 2, 3, 4.0, 5.0)
```

Fields can be excluded from flattening with the `flattenable(struct, field)`
method. These are easily defined using `@flattenable` from
[FieldMetadata.jl](https://github.com/rafaqz/FieldMetadata.jl), or defining your own
custom function with FieldMetadata, or manually with the form:

```julia
julia> import Flatten: flattenable

julia> @flattenable struct Partial{A,B,C}
           a::A | true
           b::B | true
           c::C | false
       end

julia> nestedpartial = Partial(Partial(1.0, 2.0, 3.0), 4, 5)
Partial{Partial{Float64,Float64,Float64},Int64,Int64}(Partial{Float64,Float64,Float64}(1.0, 2.0, 3.0), 4, 5)

julia> nestedpartial = Partial(Partial(1.0, 2.0, 3.0), 4, 5)
Partial{Partial{Float64,Float64,Float64},Int64,Int64}(Partial{Float64,Float64,Float64}(1.0, 2.0, 3.0), 4, 5)

julia> flatten(nestedpartial)
(1.0, 2.0, 4)
```

"""
function flatten end

flatten_builder(T, fname, action) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        $action(getfield(obj, $(QuoteNode(fname))), flattentrait, use, ignore)
    else
        ()
    end
end

flatten_combiner(T, expressions) = Expr(:tuple, expressions...)

flatten_inner(T, action) =
    nested(T, flatten_builder, flatten_combiner, action)

flatten(obj) = flatten(obj, flattenable)
flatten(obj, args...) = flatten(obj, flattenable, args...)
flatten(obj, ft::Function) = flatten(obj, ft, USE)
flatten(obj, ft::Function, use) = flatten(obj, ft, use, IGNORE)
flatten(x::I, ft::Function, use::Type{U}, ignore::Type{I}) where {U,I} = ()
flatten(x::U, ft::Function, use::Type{U}, ignore::Type{I}) where {U,I} = (_flatten(x),)
@generated flatten(obj, flattentrait::Function, use, ignore) = flatten_inner(obj, flatten)

_flatten(x) = x

"""
    reconstruct(obj, data, [flattentrait::Function], [use::Type], [ignore::Type])

Reconstruct an object from Tuple or Vector data and an existing object.
Data should be at least as long as the queried fields in the obj.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments
- `obj`: The target type to be reconstructed
- `data`: Replacement data - an `AbstractArray`, `Tuple` or type that defines `getfield`.
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
```julia
f(::Type, ::Type{Val{:fieldname}}) = true
```
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

# Examples

```julia
julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> reconstruct(Foo(1, 2, 3), (1, :two, 3.0))
Foo{Int64,Symbol,Float64}(1, :two, 3.0)
```
"""
function reconstruct end

reconstruct_builder(T, fname, action) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        val, n = $action(x, data, flattentrait, use, ignore, n)
        val
    else
        (getfield(obj, $(QuoteNode(fname))),)
    end
end

reconstruct_combiner(T, expressions) = :($(Expr(:tuple, expressions...)), n)

reconstruct_inner(::Type{T}, action) where T =
    nested(T, reconstruct_builder, reconstruct_combiner, action)

# Run from first data index and extract the final return value from the nested tuple
reconstruct(obj, data) = reconstruct(obj, data, flattenable)
reconstruct(obj, data, args...) = reconstruct(obj, data, flattenable, args...)
reconstruct(obj, data, ft::Function) = reconstruct(obj, data, ft, USE)
reconstruct(obj, data, ft::Function, use) = reconstruct(obj, data, ft, use, IGNORE)
# Need to extract the final return value from the nested tuple
reconstruct(obj, data, ft::Function, use, ignore) =
    reconstruct(obj, data, ft, use, ignore, firstindex(data))[1][1]
# Return value unmodified
reconstruct(x::I, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (x, ), n
# Return value from data. Increment position counter -  the returned n + 1 becomes n
reconstruct(x::U, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} =
    (_reconstruct(x, data[n]),), n + 1
@generated reconstruct(obj, data, flattentrait::Function, use, ignore, n) =
    quote
        args, n = $(reconstruct_inner(obj, reconstruct))
        (constructor_of(typeof(obj))(args...),), n
    end

_reconstruct(x, data) = data

apply(func, data::Tuple{T, Vararg}) where T = (func(data[1]), apply(func, Base.tail(data))...)
apply(func, data::Tuple{}) = ()

"""
    modify(func, obj, args...)

Modify field in a type with a function

"""
modify(func, obj, args...) = reconstruct(obj, apply(func, (flatten(obj, args...))), args...)

"""
    update!(obj, data, [flattentrait::Function], [use::Type], [ignore::Type])
Update a mutable object with a `Tuple` or `Vector` of data.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `data`: Replacement data - an `AbstractArray`, `Tuple` or type that defines `getfield`.
- `flattentrait`: A function returning a Bool, such as a FielMetadat method. With the form:
```julia
f(::Type, ::Type{Val{:fieldname}}) = true
```
- `use`: Types to return in the tuple.
- `ignore`: Types to ignore completly. These are not reurned or recursed over.

# Examples

```jldoctest
julia> using Flatten

julia> mutable struct MutableFoo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> mufoo = MutableFoo(1, 2, 3)
MutableFoo{Int64,Int64,Int64}(1, 2, 3)

julia> update!(mufoo, (2, 4, 6))
MutableFoo{Int64,Int64,Int64}(2, 4, 6)

julia> mufoo = MutableFoo(1, 2, :three)
MutableFoo{Int64,Int64,Symbol}(1, 2, :three)

julia> update!(mufoo, (:four,), Symbol)
MutableFoo{Int64,Int64,Symbol}(1, 2, :four)
```
"""
function update! end

update_builder(T, fname, action) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        val, n = $action(x, data, flattentrait, use, ignore, n)
        setfield!(obj, $(QuoteNode(fname)), val[1])
    end
    ()
end

# Use the reconstruct builder for tuples, they're immutable
update_builder(T::Type{<:Tuple}, fname, action) = reconstruct_builder(T, fname, action)

update_combiner(T, expressions) = :($(Expr(:tuple, expressions...)); ((obj,), n))
update_combiner(T::Type{<:Tuple}, expressions) = reconstruct_combiner(T, expressions)

update_inner(::Type{T}, action) where T =
    nested(T, update_builder, update_combiner, action)

update!(obj, data) = update!(obj, data, flattenable)
update!(obj, data, args...) = update!(obj, data, flattenable, args...)
update!(obj, data, ft::Function) = update!(obj, data, ft, USE)
update!(obj, data, ft::Function, use) = update!(obj, data, ft, use, IGNORE)
update!(obj, data, ft::Function, use, ignore) = begin
    update!(obj, data, ft, use, ignore, firstindex(data))[1][1]
    obj
end
update!(x::I, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (x,), n
update!(x::U, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (data[n],), n + 1
@generated update!(obj, data, flattentrait::Function, use, ignore, n) =
    quote
        x, n = $(update_inner(obj, update!))
        (obj,), n
    end

"""
    metaflatten(obj, func, [flattentrait::Function], [use::Type], [ignore::Type])

Metadata flattening. Flattens data attached to a field using a passed in function
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `func`: A function with the form:
```julia
f(::Type, ::Type{Val{:fieldname}}) = metadata
```
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
```julia
f(::Type, ::Type{Val{:fieldname}}) = true
```
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

We can flatten the @foobar metadata defined earlier:

```jldoctest
julia> using Flatten, FieldMetadata

julia> import Flatten: flattenable

julia> @metadata foobar :foo;

julia> @foobar struct Foo{A,B,C}
           a::A | :bar
           b::B | :foobar
           c::C | :foo
       end;

julia> @foobar struct Bar{X,Y}
           x::X | :foobar
           y::Y | :bar
       end;

julia> metaflatten(Foo(1, 2, Bar(3, 4)), foobar)
(:bar, :foobar, :foobar, :bar)
```
"""
function metaflatten end

metaflatten_builder(T, fname, action) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        $action(x, func, flattentrait, use, ignore, $T, Val{$(QuoteNode(fname))})
    else
        ()
    end
end

metaflatten_inner(T::Type, action) =
    nested(T, metaflatten_builder, flatten_combiner, action)

metaflatten(obj, func::Function) = metaflatten(obj, func, flattenable)
metaflatten(obj, func::Function, args...) = metaflatten(obj, func, flattenable, args...)
metaflatten(obj, func::Function, ft::Function) = metaflatten(obj, func, ft, USE)
metaflatten(obj, func::Function, ft::Function, use) = metaflatten(obj, func, ft, use, IGNORE)
metaflatten(obj, func::Function, ft::Function, use, ignore) =
    metaflatten(obj, func::Function, ft::Function, use, ignore, Nothing, Val{:none})

metaflatten(x::I, func::Function, ft::Function, use::Type{U}, ignore::Type{I}, P, fname) where {U,I} = ()
metaflatten(x::U, func::Function, ft::Function, use::Type{U}, ignore::Type{I}, P, fname) where {U,I} =
    (func(P, fname),)
# Better field names for tuples. TODO what about mixed type tuples?
# Also this could be confusing, consider removing.
metaflatten(xs::NTuple{N,Number}, func::Function, ft::Function, use, ignore, P, fname) where {N} =
    map(x -> func(P, fname), xs)
@generated metaflatten(obj, func::Function, flattentrait::Function, use, ignore, P, fname) =
    metaflatten_inner(obj, metaflatten)

# Helper functions to get field data with metaflatten

fieldname_meta(T, ::Type{Val{N}}) where N = N

"""jldoctest
    fieldnameflatten(obj, args...)

Flatten the field names of an object. Args are passed to metaflatten.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> fieldnameflatten(Foo(1, 2, 3))
(:a, :b, :c)
```
"""
fieldnameflatten(obj, args...) = metaflatten(obj, fieldname_meta, args...)

fieldtype_meta(T, ::Type{Val{N}}) where N = fieldtype(T, N)

"""
    fieldtypeflatten(obj, args...)

Flatten the field types of an object. Args are passed to metaflatten.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> fieldtypeflatten(Foo(1.0, :two, "Three"), Union{Real,String})
(Float64, String)
```
"""
fieldtypeflatten(obj, args...) = metaflatten(obj, fieldtype_meta, args...)

fieldparent_meta(T, ::Type{Val{N}}) where N = T.name.name

"""
    parentnameflatten(obj, args...)

Flatten the name of the parent type of an object. Args are passed to metaflatten.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> struct Bar{X,Y}
           x::X
           y::Y
       end

julia> parentnameflatten(Foo(1, 2, Bar(3, 4)))
(:Foo, :Foo, :Bar, :Bar)
```
"""
parentnameflatten(obj, args...) = metaflatten(obj, fieldparent_meta, args...)

fieldparenttype_meta(T, ::Type{Val{N}}) where N = T

"""
    parenttypeflatten(obj, args...)

Flatten the parent type of an object. Args are passed to metaflatten.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> struct Bar{X,Y}
           x::X
           y::Y
       end

julia> parenttypeflatten(Foo(1, 2, Bar(3, 4)))
(Foo{Int64,Int64,Bar{Int64,Int64}}, Foo{Int64,Int64,Bar{Int64,Int64}}, Bar{Int64,Int64}, Bar{Int64,Int64})
```
"""
parenttypeflatten(obj, args...) = metaflatten(obj, fieldparenttype_meta, args...)

end # module


text:
struct Metadata
    cell_dofs::Matrix{Int}
    dof_cells::Vector{Tuple{Int,Int}}
    dof_cells_offset::Vector{Int}
    node_cells::Vector{Tuple{Int,Int}}
    node_cells_offset::Vector{Int}
    node_dofs::Matrix{Int}
end

function Metadata(dh::DofHandler{dim}) where dim
    cell_dofs = get_cell_dofs_matrix(dh)
    dof_cells, dof_cells_offset = get_dof_cells_matrix(dh, cell_dofs)
    node_cells, node_cells_offset = get_node_cells(dh)
    node_dofs = get_node_dofs(dh)

    meta = Metadata(cell_dofs, dof_cells, dof_cells_offset, node_cells, node_cells_offset, node_dofs)
end

function get_cell_dofs_matrix(dh)
    cell_dofs = zeros(Int, ndofs_per_cell(dh), getncells(dh.grid))
    for i in 1:size(cell_dofs, 2)
        r = dh.cell_dofs_offset[i]:dh.cell_dofs_offset[i+1]-1
        for j in 1:length(r)
            cell_dofs[j,i] = dh.cell_dofs[r[j]]
        end
    end
    cell_dofs
end

function get_dof_cells_matrix(dh, cell_dofs)
    dof_cells_vecofvecs = [Vector{Tuple{Int,Int}}() for i in 1:ndofs(dh)]
    l = 0
    for cellidx in 1:size(cell_dofs, 2)
        for localidx in 1:size(cell_dofs, 1)
            dofidx = cell_dofs[localidx, cellidx]
            push!(dof_cells_vecofvecs[dofidx], (cellidx, localidx))
            l += 1
        end
    end

    dof_cells = Tuple{Int,Int}[]
    sizehint!(dof_cells, l)
    dof_cells_offset = Int[1]
    sizehint!(dof_cells_offset, ndofs(dh)+1)

    for (dofidx, indsincells) in enumerate(dof_cells_vecofvecs)
        append!(dof_cells, indsincells)
        push!(dof_cells_offset, length(dof_cells)+1)
    end

    dof_cells, dof_cells_offset
end

function get_node_cells(dh)
    node_cells_vecofvecs = [Vector{Tuple{Int,Int}}() for i in 1:ndofs(dh)]
    l = 0
    for (cellidx, cell) in enumerate(CellIterator(dh))
        for (localidx, nodeidx) in enumerate(cell.nodes)
            push!(node_cells_vecofvecs[nodeidx], (cellidx, localidx))
            l += 1
        end
    end

    node_cells = Tuple{Int,Int}[]
    node_cells_offset = Int[1]
    sizehint!(node_cells, l)
    sizehint!(node_cells_offset, getnnodes(dh.grid) + 1)

    for (nodeidx, indsincells) in enumerate(node_cells_vecofvecs)
        append!(node_cells, indsincells)
        push!(node_cells_offset, length(node_cells)+1)
    end

    return node_cells, node_cells_offset
end

node_field_offset(dh, f) = sum(view(dh.field_dims, 1:f-1))

function get_node_dofs(dh::DofHandler)
    ndofspernode = sum(dh.field_dims)
    nfields = length(dh.field_dims)
    nnodes = getnnodes(dh.grid)
    _celldofs = fill(0, ndofs_per_cell(dh))
    node_dofs = zeros(Int, ndofspernode, nnodes)
    visited = BitVector(nnodes)
    visited .= false
    for field in 1:nfields
        interpol_points = getnbasefunctions(dh.field_interpolations[field])
        field_dim = dh.field_dims[field]
        node_offset = node_field_offset(dh, field)
        offset = JuAFEM.field_offset(dh, dh.field_names[field])
        for (cellidx, cell) in enumerate(dh.grid.cells)
            celldofs!(_celldofs, dh, cellidx) # update the dofs for this cell
            for idx in 1:min(interpol_points, length(cell.nodes))
                node = cell.nodes[idx]
                if !visited[node]
                    noderange = (offset + (idx-1)*field_dim + 1):(offset + idx*field_dim) # the dofs in this node
                    for i in 1:field_dim
                        node_dofs[node_offset+i,node] = _celldofs[noderange[i]]
                    end
                    visited[node] = true
                end
            end
        end
    end

    return node_dofs
end


text:
# some basic functionality testing
x = 1234.56789

@test fmt(x) == "1234.567890"
@test fmt(x;prec=2) == "1234.57"
@test fmt(x,10,3) == "  1234.568"
@test fmt(x,10,3,:left) == "1234.568  "
@test fmt(x,10,3,:ljust) == "1234.568  "
@test fmt(x,10,3,:right) == "  1234.568"
@test fmt(x,10,3,:lrjust) == "  1234.568"
@test fmt(x,10,3,:zpad) == "001234.568"
@test fmt(x,10,3,:zeropad) == "001234.568"
@test fmt(x,:commas) == "1,234.567890"
@test fmt(x,10,3,:left,:commas) == "1,234.568 "
@test fmt(x,:ipre) == "1234.567890"
@test fmt(x,12) == " 1234.567890"

i = 1234567

@test fmt(i) == "1234567"
@test fmt(i,:commas) == "1,234,567"

@test_throws ErrorException fmt_default(Real)
@test_throws ErrorException fmt_default(Complex)

fmt_default!(Int, :commas, width = 12)
@test fmt(i) == "   1,234,567"
@test fmt(x) == "1234.567890"  # default hasn't changed

fmt_default!(:commas)
@test fmt(i) == "   1,234,567"
@test fmt(x) == "1,234.567890"  # width hasn't changed, but added commas

fmt_default!(Int) # resets Integer defaults
@test fmt(i) == "1234567"
@test fmt(i,:commas) == "1,234,567"

reset!(Int)
fmt_default!(UInt16, 'd', :commas)
@test fmt(0xffff) == "65,535"
fmt_default!(UInt32, UInt16, width=20)
@test fmt(0xfffff) == "           1,048,575"


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# 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

# only returns real roots
function quadraticroots(a::T, b::T, c::T) where {T<:Real}
    temp = b^2 - 4 * a * c

    if temp < zero(T)
        return nothing # no real roots so no ray cylinder intersection
    end

    radical = sqrt(temp)
    if b >= zero(T)
        radicalterm = -b - radical
        x1 = radicalterm / (2a)
        x2 = (2c) / radicalterm
    else
        radicalterm = -b + radical
        x1 = (2c) / radicalterm
        x2 = radicalterm / (2a)
    end

    return SVector{2,T}(x1, x2)
end

replprint(a) = show(IOContext(stdout), "text/plain", a)

@inline samepoint(pt1, pt2) = isapprox(pt1, pt2, atol = 1e-10)

password() = String(collect(rand(('1':'9'..., 'A':'Z'..., 'a':'z'...)) for i in 1:14))

# 1 if val>= 0 and -1 otherwise. Unlike Math.Sign which returns 0 if val==0

@inline sign(val::T) where {T<:Real} = val >= zero(T) ? 1 : -1

@inline function NaNsafeatan(x::T, y::T)::T where {T<:Real}
    if y == zero(T)
        if x == zero(T)
            return zero(T)
        elseif x > zero(T)
            return T(π / 2)
        else
            return -T(π / 2)
        end
    else
        return atan(x, y)
    end
end

@inline function NaNsafeasin(x::T)::T where {T<:Real}
    if x == zero(T)
        return zero(T)
    elseif x == one(T)
        return T(π / 2)
    elseif x == -one(T)
        return -T(π / 2)
    else
        return asin(x)
    end
end

@inline function NaNsafeacos(x::T)::T where {T<:Real}
    if x == zero(T)
        return T(π / 2)
    elseif x == one(T)
        return zero(T)
    elseif x == -one(T)
        return T(π)
    else
        return acos(x)
    end
end


text:
using CairoMakie
using LinearAlgebra

# Vector of NamedTuples, ensemble mean at each iteration
ensemble_means(eki) = getproperty.(eki.iteration_summaries, :ensemble_mean)

# N_param x N_iter matrix, ensemble covariance at each iteration
ensemble_std(eki) = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...))

parameter_names(eki) = eki.inverse_problem.free_parameters.names

function plot_parameter_convergence!(eki, directory; true_parameters=nothing, n_columns=3)

    means = ensemble_means(eki)
    θθ_std_arr = ensemble_std(eki)

    pnames = parameter_names(eki)
    N_param = length(pnames)
    N_iter = length(eki.iteration_summaries) - 1 # exclude 0th element

    n_rows = Int(ceil(N_param / n_columns))
    ax_coords = [(i, j) for i = 1:n_rows, j = 1:n_columns]

    fig = Figure(resolution = (500n_columns, 200n_rows))
    for (i, pname) in enumerate(pnames)
        coords = ax_coords[i]
        ax = Axis(fig[coords...],
            xlabel = "Iteration",
            xticks = 0:N_iter,
            ylabel = string(pname))
        ax.ylabelsize = 20

        mean_values = [getproperty.(means, pname)...]
        lines!(ax, 0:N_iter, mean_values)
        band!(ax, 0:N_iter, mean_values .+ θθ_std_arr[i, :], mean_values .- θθ_std_arr[i, :])
        isnothing(true_parameters) || hlines!(ax, [true_parameters[pname]], color = :red)
    end

    save(joinpath(directory, "parameter_convergence.png"), fig)
end

function plot_pairwise_ensembles!(eki, directory, true_parameters=nothing)

    path = joinpath(directory, "pairwise_ensembles")
    isdir(path) || mkdir(path)

    pnames = parameter_names(eki)
    N_param = length(pnames)
    N_iter = length(eki.iteration_summaries) - 1 # exclude 0th element
    for (i1, pname1) in enumerate(pnames), (i2, pname2) in enumerate(pnames)
        if i1 < i2

            f = Figure()
            axtop = Axis(f[1, 1])
            axmain = Axis(f[2, 1], xlabel = string(pname1), ylabel = string(pname2))
            axright = Axis(f[2, 2])
            scatters = []
            for iteration in [0, 1, N_iter]
                ensemble = eki.iteration_summaries[iteration].parameters
                ensemble = [[particle[pname1], particle[pname2]] for particle in ensemble]
                ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2
                push!(scatters, scatter!(axmain, ensemble))
                density!(axtop, ensemble[:, 1])
                density!(axright, ensemble[:, 2], direction = :y)
            end
            isnothing(true_parameters) || begin
                vlines!(axmain,  [true_parameters[pname1]], color = :red)
                vlines!(axtop,   [true_parameters[pname1]], color = :red)
                hlines!(axmain,  [true_parameters[pname2]], color = :red)
                hlines!(axright, [true_parameters[pname2]], color = :red)
            end
            colsize!(f.layout, 1, Fixed(300))
            colsize!(f.layout, 2, Fixed(200))
            rowsize!(f.layout, 1, Fixed(200))
            rowsize!(f.layout, 2, Fixed(300))
            Legend(f[1, 2], scatters,
                ["Initial ensemble", "Iteration 1", "Iteration $N_iter"],
                position = :lb)
            hidedecorations!(axtop, grid = false)
            hidedecorations!(axright, grid = false)
            xlims!(axright, 0, 10)
            ylims!(axtop, 0, 10)
            save(joinpath(path, "pairwise_ensembles_$(pname1)_$(pname2).png"), f)
        end
    end
end

function plot_error_convergence!(f, eki, directory; true_parameters=nothing)

    means = ensemble_means(eki)
    N_iter = length(eki.iteration_summaries) - 1 # exclude 0th element
    y = eki.mapped_observations

    output_distances = [mapslices(norm, (forward_map(eki.inverse_problem, [means...])[:, 1:(N_iter+1)] .- y), dims = 1)...]
    scatterlines(f[1, 1], 0:N_iter, output_distances, color = :blue, linewidth = 2,
        axis = (title = "Output distance",
            xlabel = "Iteration",
            ylabel = "|G(θ̅ₙ) - y|",
            yscale = log10))

    isnothing(true_parameters) || begin
        weight_distances = [norm(collect(means[iter]) .- collect(true_parameters)) for iter = 0:N_iter]
        scatterlines(f[1, 2], 0:N_iter, weight_distances, color = :red, linewidth = 2,
            axis = (title = "Parameter distance",
                xlabel = "Iteration",
                ylabel = "|θ̅ₙ - θ⋆|",
                yscale = log10))
    end

    nothing
end

function plot_error_convergence!(eki, directory; true_parameters=nothing)

    f = Figure()

    plot_error_convergence!(f, eki, directory; true_parameters)

    save(joinpath(directory, "error_convergence_summary.png"), f);
end


text:
using OpenNI2
using Cxx
using Base.Test

const ni2 = OpenNI2

@testset "Version" begin
    ver = ni2.getVersionNumber()
    @test ver >= v"2.2.0-"
end

@testset "Basics" begin
    ni2.initialize()
    device = ni2.DevicePtr()
    @test ni2.open(device) == ni2.STATUS_OK
    @test ni2.isValid(device)
    @test ni2.hasSensor(device, ni2.SENSOR_DEPTH)

    depth = ni2.VideoStreamPtr()
    @test ni2.create(depth, device, ni2.SENSOR_DEPTH) == ni2.STATUS_OK
    @test ni2.isValid(device)

    rc = ni2.setImageRegistrationMode(device, ni2.IMAGE_REGISTRATION_OFF)
    @test rc == ni2.STATUS_OK

    mode = ni2.getVideoMode(depth)
    w, h = ni2.getResolutionX(mode), ni2.getResolutionY(mode)

    @test ni2.start(depth) == ni2.STATUS_OK
    frame = ni2.VideoFrameRef()
    @test !ni2.isValid(frame)
    @test ni2.waitForAnyStream([depth]) == 0
    ni2.readFrame(depth, frame)
    @test ni2.isValid(frame)
    arr = convert(Array{Cushort,2}, frame)
    @test size(arr) == (w, h)
    ni2.stop(depth)
    ni2.destroy(depth)

    ni2.close(device)
    ni2.shutdown()

    @test !ni2.isValid(device)
end

@testset "DeviceInfo" begin
    ni2.initialize()
    device = ni2.DevicePtr()
    @test ni2.open(device) == ni2.STATUS_OK
    @test ni2.isValid(device)

    di = ni2.getDeviceInfo(device)
    @test !isempty(ni2.getUri(di))
    @test !isempty(ni2.getVendor(di))
    @test !isempty(ni2.getName(di))
    @test ni2.getUsbVendorId(di) != 0x00
    @test ni2.getUsbProductId(di) != 0x00

    ni2.close(device)
    ni2.shutdown()
end

@testset "setVideoMode" begin
    ni2.initialize()
    device = ni2.DevicePtr()
    @test ni2.open(device) == ni2.STATUS_OK
    depth = ni2.VideoStreamPtr()
    @test ni2.create(depth, device, ni2.SENSOR_DEPTH) == ni2.STATUS_OK

    si = ni2.getSensorInfo(device, ni2.SENSOR_DEPTH)
    modes = ni2.getSupportedVideoModes(si)
    @test length(modes) > 0

    # target
    w = 640
    h = 480
    pxfmt = ni2.PIXEL_FORMAT_DEPTH_1_MM

    found = false
    for mode in modes
        if Int(ni2.getResolutionX(mode)) == w &&
                Int(ni2.getResolutionY(mode)) == h &&
                ni2.getPixelFormat(mode) == pxfmt
            found = true
            ni2.setVideoMode(depth, mode)
            break
        end
    end
    @test found

    @test ni2.start(depth) == ni2.STATUS_OK
    frame = ni2.VideoFrameRef()
    @test ni2.waitForAnyStream([depth]) == 0
    ni2.readFrame(depth, frame)
    arr = convert(Array{Cushort,2}, frame)
    @test size(arr) == (w, h)
    ni2.stop(depth)
    ni2.destroy(depth)

    ni2.close(device)
    ni2.shutdown()
end


text:
module CameraManager

export Camera

using CoreExtended
using MatrixMath
using JLGEngine

using ..TransformManager

const Transform = TransformManager.Typ

""" TODO """
type Camera <: JLGEngine.IComponent
	id::Symbol
	transform::Transform

	ortho::Bool
	fov::Float32
	ratio::Float32
	depth::Vec2f

  viewport::Vec4f

  viewMat::Mat4x4f
  projectionMat::Mat4x4f
	modelMat::Mat4x4f

	translateMat::Mat4x4f
	rotationMat::Mat4x4f
	scalingMat::Mat4x4f

	mvpMat::Mat4x4f

	Camera(id=:_) = new(id,nothing,false,0,0,Vec2f(0,1),zeros(Vec4f),IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f,IMat4x4f)
end

""" TODO """
function init(this::Camera)
	this.transform=TransformManager.create(this.id)
end

# preset for manager: includes code here
presetManager(Camera)

""" TODO """
function zoom(this::Camera, zf::Float32)
	this.transform.scaling = ones(Vec3f) * (1+zoom*0.01f0)
end

""" TODO """
function updateTranslation(this::Camera)
	this.translateMat = MatrixMath.translationmatrix(this.transform.position)
end

""" TODO """
function updateRotation(this::Camera)
	this.rotationMat = MatrixMath.computeRotation(this.transform.rotation)
end

""" TODO """
function updateScaling(this::Camera)
	this.scalingMat = MatrixMath.scalingmatrix(this.transform.scaling)
end

""" TODO """
function updateView(this::Camera)
	updateTranslation(this)
	updateRotation(this)
	updateScaling(this)
	this.viewMat = this.rotationMat * this.translateMat * this.scalingMat
end

""" TODO """
function setRatio(this::Camera, ratio::Float32)
	this.ratio=ratio
end

""" TODO """
function setPerspective(this::Camera, fov::Float32,ratio::Float32,near::Float32,far::Float32, ortho=false) #(60.0f0, Float32(windowsize[1]/windowsize[2]), 0.1f0, 10.0f0)
	this.fov=fov
	this.ratio=ratio
	this.depth=Vec2f(near,far)
	this.ortho=ortho
end

""" TODO """
function updateProjection(this::Camera)
	if !this.ortho this.projectionMat = MatrixMath.perspectiveprojection(this.fov, this.ratio, this.depth[1], this.depth[2])
	else this.projectionMat = MatrixMath.orthographicprojection(this.fov, this.ratio, this.depth[1], this.depth[2])
	end
end

""" TODO """
function updateMVP(this::Camera)
	this.mvpMat = this.projectionMat * this.viewMat * this.modelMat
end

""" TODO """
function updateModel(this::Camera, transform::Transform)
	if transform == nothing return end
	this.modelMat = MatrixMath.translationmatrix(transform.position) * MatrixMath.computeRotation(transform.rotation) * MatrixMath.scalingmatrix(transform.scaling)
end

""" TODO """
function update(this::Camera, transform::Transform)
	updateProjection(this)
	updateView(this)
	updateModel(this, transform)
	updateMVP(this)
end

end #CameraManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using VLConstraintBasedModelGenerationUtilities
using BioSymbols

# setup path to sequence file -
path_to_gene_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/G-MZ373340.fasta"

# load -
gene_table = build_gene_table(path_to_gene_sequence_file) |> check

# get the sequence -
gene_seq = gene_table[!,:gene_sequence][1]

# setup local complementOp function -
function complementOp(nucleotide::BioSymbols.DNA)::BioSymbols.RNA

    if nucleotide == DNA_T
        return RNA_U
    elseif nucleotide == DNA_A
        return RNA_A
    elseif nucleotide == DNA_C
        return RNA_C
    elseif nucleotide == DNA_G
        return RNA_G
    else
        return RNA_N
    end
end

# table -
transcription_table = transcribe(gene_seq; complementOperation=complementOp) |> check


text:
# Test naive implementation of entropy calculations towards efficient calculation of entropy on a manifold

export
  logmap_SO2,
  difftheta,
  addtheta,
  rbfAccAt!,
  rbf!,
  rbf,
  evaluateManifoldNaive1D!,
  manifoldLooCrossValidation,
  kde!_CircularNaiveCV,
  getCircMu,
  getCircLambda

global const reci_s2pi=1.0/sqrt(2.0*pi) # 1.0/2.5066282746310002

# On-manifold circular product callbacks

# manifold distance, add, and subtract
function logmap_SO2(Rl::Matrix{Float64})
  ct = abs(Rl[1,1]) > 1.0 ? 1.0 : Rl[1,1]  # reinserting the sign below
  -sign(Rl[2,1])*acos(ct)
end
difftheta(wth1, wth2)::Float64 = logmap_SO2(TUs.R(wth1)'*TUs.R(wth2))
addtheta(wth1, wth2)::Float64 = TUs.wrapRad( wth1+wth2 )

# manifold get Gaussian products mean
getCircMu = (m::Vector{Float64}, s::Vector{Float64}, dummy::Float64) -> TUs.wrapRad(get2DMu(m, s, diffop=difftheta, initrange=(-pi+0.0,pi+0.0)))
# getCircMu = (m::Vector{Float64}, s::Vector{Float64}, dummy::Float64) -> TUs.wrapRad(get2DMuMin(m, s, diffop=difftheta, initrange=(-pi+0.0,pi+0.0)))

getCircLambda(x) = getEuclidLambda(x)

"""
    $SIGNATURES

Probability density function `p(x)`, as estimated by kernels
```math
hatp_{-j}(x) = 1/(N-1) Σ_{i != j}^N frac{1}{sqrt{2pi}σ } exp{ -frac{(x-μ)^2}{2 σ^2} }
```
"""
function normDistAccAt!(ret::AV,
                        idx::Int,
                        x::Float64,
                        sigma::Float64,
                        w::Float64=1.0  )::Nothing where {AV <: AbstractVector}
    global reci_s2pi
    @fastmath ret[idx] += w*reci_s2pi/sigma * exp( - (x^2)/(2.0*(sigma^2)) )
    return nothing
end

function rbfAccAt!(ret::AV,
                   idx::Int,
                   x::Float64,
                   μ::Float64=0.0,
                   σ::Float64=1.0,
                   w::Float64=1.0,
                   diffop::Function=-)::Nothing  where {AV <: AbstractVector}
    #
    normDistAccAt!(ret, idx, diffop(x, μ), σ, w)
    nothing
end
function rbf!(ret::AV,
              x::Float64,
              μ::Float64=0.0,
              σ::Float64=1.0,
              diffop::Function=-)::Nothing  where {AV <: AbstractVector}
    #
    ret[1] = 0.0
    normDistAccAt!(ret, 1, diffop(x,μ), σ)
    nothing
end

function rbf(x::Float64, μ::Float64=0.0, σ::Float64=1.0)
    ret = Vector{Float64}(undef, 1) # initialized in rbf!(..)
    rbf!(ret, x, μ, σ)
    return ret[1]
end

"""
    $SIGNATURES

Evalute the KDE naively as equally weighted Gaussian kernels with common bandwidth.
This function does, however, allow on-manifold evaluations.
"""
function evaluateManifoldNaive1D!(ret::Vector{Float64},
                                  idx::Int,
                                  pts::Array{Float64,1},
                                  bw::Float64,
                                  x::Array{Float64,1},
                                  loo::Int=-1,
                                  diffop=-  )::Nothing
    #
    dontskip = loo == -1
    N = length(pts)
    reci_N = dontskip ? 1.0/N : 1.0/(N-1)
    for j in 1:N
        if dontskip || loo != j
            manifolddist = diffop(pts[loo], pts[j])
            normDistAccAt!(ret, idx, manifolddist, bw, reci_N)
        end
    end

    return nothing
end
function evaluateManifoldNaive1D!(ret::Vector{Float64},
                                  idx::Int,
                                  bd::BallTreeDensity,
                                  x::Array{Float64,1},
                                  loo::Int=-1,
                                  diffop=-  )::Nothing
    #
    evaluateManifoldNaive1D!(ret, idx, getPoints(bd)[:], getBW(bd)[1,1], x, loo, diffop )
end

"""
    $SIGNATURES

Calculate negative entropy with leave one out (j'th element) cross validation.

Background
==========

From: Silverman, B.: Density Estimation for Statistics and Data Analysis, 1986, p.52

Probability density function `p(x)`, as estimated by kernels
```math
hatp_{-j}(x) = 1/(N-1) Σ_{i != j}^N frac{1}{sqrt{2pi}σ } exp{ -frac{(x-μ)^2}{2 σ^2} }
```
and has Cross Validation number as the average log evaluations of leave one out `hatp_{-j}(x)`:
```math
CV(p) = 1/N Σ_i^N log hat{p}_{-j}(x_i)
```

This quantity `CV` is related to an entropy `H(p)` estimate via:
```math
H(p) = -CV(p)
```
"""
function manifoldLooCrossValidation(pts::Array,
                            bw::Float64;
                            own::Bool=true,
                            diffop::Function=-  )
    #
    N = maximum(size(pts))
    h = [bw;]
    loo = zeros(N)
    @inbounds for i in 1:N
        if !own
            # validation version
            pts99 = pts[[1:(i-1);(i+1):end]]
            p99 = kde!(pts99, h)
            loo[i] = log(p99([pts[i];])[1])
        else
            # own naive entropy calculation
            loo[i] = 0.0
            evaluateManifoldNaive1D!(loo, i, pts, bw, pts, i, diffop)
            loo[i] = log(loo[i])
        end
    end
    return  sum(loo)/N
end

function kde!_CircularNaiveCV(points::A) where {A <: AbstractArray{Float64,1}}

  # initial setup parameters
  dims = 1 # size(points,1)
  bwds = zeros(dims)
  # initial testing values
  lower = 0.001
  upper = 2pi

  # excessive for loop for leave one out likelihiood cross validation (Silverman 1986, p.52)
  for i in 1:dims
    minEntropyLOOCV = (bw) -> -manifoldLooCrossValidation(points, bw, own=true, diffop=difftheta)
    res = optimize(minEntropyLOOCV, lower, upper, GoldenSection(), x_tol=0.001)
    bwds[i] = res.minimizer
  end

  # cosntruct the kde with CV optimized bandwidth
  p = kde!( points, bwds, (addtheta,), (difftheta,) )

  return p
end

#


text:
# Like @test_throws, but for macros.
# Julia-0.7 and later always wrap the exception in a LoadError.

@static if VERSION >= v"0.7.0-alpha.0"
   macro test_macro_throws(typ, expr)
      quote
         @test_throws $typ begin
            try
               $expr
            catch e
               rethrow(e.error)
            end
         end
      end
   end
else
   macro test_macro_throws(typ, expr)
      quote
         @test_throws $typ $expr
      end
   end
end


text:
function emit_node(program::Program, symbol_map, call_node::CallNode)
    call_expr = Expr(:call, call_node.f)

    for arg in call_node.args
        if arg isa TrackedValue
            push!(call_expr.args, symbol_map[arg.id])
        else
            push!(call_expr.args, arg)
        end
    end

    return call_expr
end

function emit_node(program::Program, symbol_map, zero_node::ZeroNode)
    return :(zero($(symbol_map[zero_node.node_id])))
end

function emit_node(program::Program, symbol_map, tuple_node::TupleNode)
    tuple_expr = Expr(:tuple)

    for arg in tuple_node.args
        if arg isa TrackedValue
            push!(tuple_expr.args, symbol_map[arg.id])
        else
            push!(tuple_expr.args, arg)
        end
    end

    return tuple_expr
end

function emit_subprogram(program::Program, subprograms, symbol_map, i::Int)
    sp = subprograms[i]

    block = Expr(:block)

    for j in length(program.args[i])+1:length(sp)
        node = program.nodes[sp[j]]
        node_expr = emit_node(program, symbol_map, node)
        push!(block.args, Expr(:(=), symbol_map[sp[j]], node_expr))
    end

    if length(subprograms) > i
        push!(block.args, emit_subprogram(program, subprograms, symbol_map, i+1))
    else
        push!(block.args, Expr(:return, symbol_map[program.return_node]))
    end

    subprogram_args = program.args[i]

    args_expr = if length(subprogram_args) > 1
        args_tuple_expr = Expr(:tuple)

        for node_id in subprogram_args
            push!(args_tuple_expr.args, symbol_map[node_id])
        end

        args_tuple_expr
    elseif length(subprogram_args) == 1
        symbol_map[subprogram_args[1]]
    else
        error("we do not support subprograms with zero arguments")
    end

    return Expr(:->, args_expr, block)
end

function emit_code(program::Program)
    subprograms = create_subprograms(program)
    executable_order = reduce(vcat, subprograms)
    symbol_map = Dict{Int,Symbol}()

    for i in 1:length(executable_order)
        symbol_map[executable_order[i]] = Symbol("s", i)
    end

    return emit_subprogram(program, subprograms, symbol_map, 1)
end


text:
export move_maker

"""
`move_maker(T::Tableau)` generates the list of all possible next positions.
"""
function move_maker(T::Tableau)::Vector{Tableau}
    out = Set{Tableau}()   # place to put the answer

    # PART 1: Move last card in cascade piles

    for k = 1:8
        p = T.Casc.piles[k]
        if length(p) == 0
            continue
        end
        C = p[end]  # last card in this pile

        # 1a Move to another cascade pile

        for j = 1:8
            if j != k
                if add_check(T.Casc.piles[j], C)  # see if we can add C to pile j
                    TT = deepcopy(T)
                    rm_card(TT.Casc, k)
                    add_card(TT.Casc, j, C)
                    push!(out, TT)
                end

            end
        end

        # 1b Move to foundation

        if add_check(T.Found, C)
            TT = deepcopy(T)
            rm_card(TT.Casc, k)
            add_card(TT.Found, C)
            push!(out, TT)
        end

        # 1c Move to a freecell

        if add_check(T.Free, C)
            TT = deepcopy(T)
            rm_card(TT.Casc, k)
            add_card(TT.Free, C)
            push!(out, TT)
        end

    end

    # Part 2: Move a card from a Free Cell

    for C in T.Free

        # 2a move to Tableau

        for k = 1:8
            if add_check(T.Casc.piles[k], C)
                TT = deepcopy(T)
                rm_card(TT.Free, C)
                add_card(TT.Casc, k, C)
                push!(out, TT)
            end
        end

        # 2b move to Foundation
        if add_check(T.Found, C)
            TT = deepcopy(T)
            rm_card(TT.Free, C)
            add_card(TT.Found, C)
            push!(out, TT)
        end

    end
    return collect(out) # report as a list for ImplicitGraph
end


text:
using Ipopt

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

function eval_f(x::Vector{Float64})
  return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
  # Bad: g    = zeros(2)  # Allocates new array
  # OK:  g[:] = zeros(2)  # Modifies 'in place'
  g[1] = x[1]   * x[2]   * x[3]   * x[4]
  g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
  # Bad: grad_f    = zeros(4)  # Allocates new array
  # OK:  grad_f[:] = zeros(4)  # Modifies 'in place'
  grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
  grad_f[2] = x[1] * x[4]
  grad_f[3] = x[1] * x[4] + 1
  grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
  if mode == :Structure
    # Constraint (row) 1
    rows[1] = 1; cols[1] = 1
    rows[2] = 1; cols[2] = 2
    rows[3] = 1; cols[3] = 3
    rows[4] = 1; cols[4] = 4
    # Constraint (row) 2
    rows[5] = 2; cols[5] = 1
    rows[6] = 2; cols[6] = 2
    rows[7] = 2; cols[7] = 3
    rows[8] = 2; cols[8] = 4
  else
    # Constraint (row) 1
    values[1] = x[2]*x[3]*x[4]  # 1,1
    values[2] = x[1]*x[3]*x[4]  # 1,2
    values[3] = x[1]*x[2]*x[4]  # 1,3
    values[4] = x[1]*x[2]*x[3]  # 1,4
    # Constraint (row) 2
    values[5] = 2*x[1]  # 2,1
    values[6] = 2*x[2]  # 2,2
    values[7] = 2*x[3]  # 2,3
    values[8] = 2*x[4]  # 2,4
  end
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
  if mode == :Structure
    # Symmetric matrix, fill the lower left triangle only
    idx = 1
    for row = 1:4
      for col = 1:row
        rows[idx] = row
        cols[idx] = col
        idx += 1
      end
    end
  else
    # Again, only lower left triangle
    # Objective
    values[1] = obj_factor * (2*x[4])  # 1,1
    values[2] = obj_factor * (  x[4])  # 2,1
    values[3] = 0                      # 2,2
    values[4] = obj_factor * (  x[4])  # 3,1
    values[5] = 0                      # 3,2
    values[6] = 0                      # 3,3
    values[7] = obj_factor * (2*x[1] + x[2] + x[3])  # 4,1
    values[8] = obj_factor * (  x[1])  # 4,2
    values[9] = obj_factor * (  x[1])  # 4,3
    values[10] = 0                     # 4,4

    # First constraint
    values[2] += lambda[1] * (x[3] * x[4])  # 2,1
    values[4] += lambda[1] * (x[2] * x[4])  # 3,1
    values[5] += lambda[1] * (x[1] * x[4])  # 3,2
    values[7] += lambda[1] * (x[2] * x[3])  # 4,1
    values[8] += lambda[1] * (x[1] * x[3])  # 4,2
    values[9] += lambda[1] * (x[1] * x[2])  # 4,3

    # Second constraint
    values[1]  += lambda[2] * 2  # 1,1
    values[3]  += lambda[2] * 2  # 2,2
    values[6]  += lambda[2] * 2  # 3,3
    values[10] += lambda[2] * 2  # 4,4
  end
end

function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)
  println("Iteration $iter_count, objective value is $obj_value.")
  return true
end

n = 4
x_L = [1.0, 1.0, 1.0, 1.0]
x_U = [5.0, 5.0, 5.0, 5.0]

m = 2
g_L = [25.0, 40.0]
g_U = [2.0e19, 40.0]

prob = createProblem(n, x_L, x_U, m, g_L, g_U, 8, 10,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)

prob.x = [1.0, 5.0, 5.0, 1.0]
stat = solveProblem(prob)

@test Ipopt.ApplicationReturnStatus[stat] == :Solve_Succeeded
@test_approx_eq_eps prob.x[1] 1.0000000000000000 1e-5
@test_approx_eq_eps prob.x[2] 4.7429996418092970 1e-5
@test_approx_eq_eps prob.x[3] 3.8211499817883077 1e-5
@test_approx_eq_eps prob.x[4] 1.3794082897556983 1e-5
@test_approx_eq_eps prob.obj_val 17.014017145179164 1e-5


text:
using  MarketData, FactCheck

include("timevector.jl")

exitstatus()


text:
using Base.Test

function runtests(name)
    println("     \033[1m*\033[0m \033[31m$(name)\033[0m")
    Core.include("$name.jl")
    nothing
end

function propagate_errors(a,b)
    if isa(a,Exception)
        rethrow(a)
    end
    if isa(b,Exception)
        rethrow(b)
    end
    nothing
end

# looking in . messes things up badly
filter!(x->x!=".", LOAD_PATH)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Documenter, Thebes, Luxor, Rotations

makedocs(
    modules = [Thebes],
    sitename = "Thebes",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true", sidebar_sitename=true),
    pages    = Any[
        "Introduction to Thebes"  => "index.md",
        "The basics"              => "basics.md",
        "Views"                   => "views.md",
        "Polygons and planes"     => "polys.md",
        "Objects"                 => "objects.md",
        "Text"                    => "text.md",
        "Tools"                   => "tools.md",
        "Index"                   => "functionindex.md"
        ]
    )

deploydocs(
    repo = "github.com/cormullion/Thebes.jl.git",
    target = "build"
)


text:
# Adapted from Knet's src/data.jl (author: Deniz Yuret)
using Random: AbstractRNG, shuffle!, GLOBAL_RNG
using Base: @propagate_inbounds

struct MyDataLoader{D,R<:AbstractRNG}
    data::D
    batchsize::Int
    nobs::Int
    partial::Bool
    imax::Int
    indices::Vector{Int}
    s_indices::Vector{Int}
    shuffle::Bool
    rng::R
end

_nobs(data::AbstractArray) = size(data)[end]
_nsim(data::AbstractArray) = size(data)[1]
_getobs(data::Union{Tuple, NamedTuple}, i) = map(Base.Fix2(_getobs, i), data)

function MyDataLoader(data; batchsize=1, shuffle=false, partial=true, rng=GLOBAL_RNG)
    batchsize > 0 || throw(ArgumentError("Need positive batchsize"))

    n = _nobs(data)

    if ndims(data[1]) > 2
        nsim = size(data[1])[1]
    else
        nsim = 1
    end

    if n < batchsize
        @warn "Number of observations less than batchsize, decreasing the batchsize to $n"
        batchsize = n
    end
    imax = partial ? n : n - batchsize + 1

    MyDataLoader(data, batchsize, n, partial, imax, [1:n;], [1:nsim;], shuffle, rng)
end

@propagate_inbounds function Base.iterate(d::MyDataLoader, i=0)     # returns data in d.indices[i+1:i+batchsize]
    i >= d.imax && return nothing
    if d.shuffle && i == 0
        shuffle!(d.rng, d.indices)
        shuffle!(d.rng, d.s_indices)
    end

    nexti = min(i + d.batchsize, d.nobs)
    ids = d.indices[i+1:nexti]
    next_sim = d.s_indices[1]

    batch = d.data[1][next_sim,:, ids]
    time_batch = d.data[2][next_sim,1, ids]

    return ((batch, time_batch), nexti)

end

function Base.length(d::MyDataLoader)
    n = d.nobs / d.batchsize
    d.partial ? ceil(Int,n) : floor(Int,n)
end

function _nobs(data::Union{Tuple, NamedTuple})
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))
    n = _nobs(data[1])
    for i in keys(data)
        ni = _nobs(data[i])
        n == ni || throw(DimensionMismatch("All data inputs should have the same number of observations, i.e. size in the last dimension. " *
            "But data[$(repr(first(keys(data))))] ($(summary(data[1]))) has $n, while data[$(repr(i))] ($(summary(data[i]))) has $ni."))
    end
    return n
end

Base.eltype(::MyDataLoader{D}) where D = D

max_stars_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_issues_repo_name: Romit-Maulik/Tutorials-Demos-Practice

max_forks_repo_name: Romit-Maulik/Tutorials-Demos-Practice