text:
# Tutorial 2.5. Beyond square lattices: graphene
# ==============================================
#
# Physics background
# ------------------
#  Transport through a graphene quantum dot with a pn-junction
#
# Kwant features highlighted
# --------------------------
#  - Application of all the aspects of tutorials 1-3 to a more complicated
#    lattice, namely graphene

# For computing eigenvalues
sla = pyimport("scipy.sparse.linalg")

# For plotting

# Define the graphene lattice
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
                                 [(0, 0), (0, 1 / sqrt(3))])
a, b = kwant.lattice.Monatomic.(graphene.sublattices)

function make_system(;r=10, w=2.0, pot=0.1)

    #### Define the scattering region. ####
    # circular scattering region
    function circle(pos)
        x, y = pos
        return hypot(x,y) < r
    end

    syst = kwant.Builder()

    # w: width and pot: potential maximum of the p-n junction
    py"""
    def potential(site):
        (x, y) = site.pos
        d = y * $cos_30 + x * $sin_30
        return $pot * $tanh(d / $w)
    """

    syst[graphene.shape(circle, (0, 0))] = py"potential"

    # specify the hoppings of the graphene lattice in the
    # format expected by builder.HoppingKind
    hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
    syst[[kwant.builder.HoppingKind(hopping) for hopping in hoppings]] = -1

    # Modify the scattering region
    deleteat!(syst,a(0,0))
    syst[a(-2, 1), b(2, 2)] = -1

    #### Define the leads. ####
    # left lead
    sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))

    function lead0_shape(pos)
        x, y = pos
        return (-0.4 * r < y < 0.4 * r)
    end

    lead0 = kwant.Builder(sym0)
    lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
    lead0[[kwant.builder.HoppingKind(hopping) for hopping in hoppings]] = -1

    # The second lead, going to the top right
    sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))

    function lead1_shape(pos)
        v = pos[2] * sin_30 - pos[1] * cos_30
        return (-0.4 * r < v < 0.4 * r)
    end

    lead1 = kwant.Builder(sym1)
    lead1[graphene.shape(lead1_shape, (0, 0))] = pot
    lead1[[kwant.builder.HoppingKind(hopping) for hopping in hoppings]] = -1

    return syst, [lead0, lead1]
end

function compute_evs(syst)
    # Compute some eigenvalues of the closed system
    sparse_mat = syst.hamiltonian_submatrix(sparse=true)

    evs = sla.eigs(sparse_mat, 2)[1]
    print(real(evs))
end

function plot_conductance(syst, energies)
    # Compute transmission as a function of energy
    data = []
    for energy in energies
        smatrix = kwant.smatrix(syst, energy)
        append!(data,smatrix.transmission(0, 1))
    end

    pyplot.figure()
    pyplot.plot(energies, data)
    pyplot.xlabel("energy [t]")
    pyplot.ylabel("conductance [e^2/h]")
    pyplot.gcf()
end

function plot_bandstructure(flead, momenta)
    bands = kwant.physics.Bands(flead)
    energies = [bands(k) for k in momenta]

    pyplot.figure()
    pyplot.plot(momenta, energies)
    pyplot.xlabel("momentum [(lattice constant)^-1]")
    pyplot.ylabel("energy [t]")
    pyplot.gcf()
end


text:
import Random: CloseOpen01_64, CloseOpen12_64, FloatInterval, MT_CACHE_F

mutable struct MemorizingRNG{S} <: AbstractRNG
    memory::Vector{Float64}
    start::Int
    finish::Int
    idx::Int
    source::S
end

MemorizingRNG(source) = MemorizingRNG(Float64[], 1, 0, 0, source)

# Low level API (copied from MersenneTwister)
@inline mr_avail(r::MemorizingRNG) = r.finish - r.idx
@inline mr_empty(r::MemorizingRNG) = r.idx == r.finish
@inline mr_pop!(r::MemorizingRNG) = @inbounds return r.memory[r.idx+=1]

@inline reserve_1(r::MemorizingRNG) = mr_empty(r) && gen_rand!(r, 1)
@inline reserve(r::MemorizingRNG, n::Integer) = mr_avail(r) < n && gen_rand!(r, n)

# precondition: !mr_empty(r)
@inline rand_inbounds(r::MemorizingRNG, ::Type{CloseOpen12_64}) = mr_pop!(r)
@inline rand_inbounds(r::MemorizingRNG, ::Type{CloseOpen01_64}) = rand_inbounds(r, CloseOpen12_64) - 1.0
@inline rand_inbounds(r::MemorizingRNG) = rand_inbounds(r, CloseOpen01_64)

Random.rng_native_52(rng::MemorizingRNG) = Random.rng_native_52(rng.source)

function rand(r::MemorizingRNG, ::Union{I, Type{I}}) where I <: FloatInterval
    reserve_1(r)
    rand_inbounds(r, I)
end

# rand(r::MemorizingRNG, ::Type{Float64}) = rand(r, CloseOpen01)

function gen_rand!(r::MemorizingRNG{MersenneTwister}, n::Integer)
    len = length(r.memory)
    if len < r.idx + n
        resize!(r.memory, len+MT_CACHE_F)
        Random.gen_rand(r.source) # could be faster to use dsfmt_fill_array_close1_open2
        r.memory[len+1:end] = r.source.vals
        Random.mt_setempty!(r.source)
    end
    r.finish = length(r.memory)
    return r
end

ext: jl


text:
# create dataframe of each group by hour
# calculate probability of each distance

module GrpHrHeatmap
    using DataFrames, CSV, FreqTables, Plots
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_org = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        hours = [Int64(floor(e/10000)) for e in df_org.time]

        # add hour array to df
        df_h = DataFrame(hour=hours)
        df_new = hcat(df_org, df_h)

        # pivot
        pivot = freqtable(df_new, :lidar, :hour)
        println(pivot)

        #heatmap
        x = names(pivot, 2)
        y = names(pivot, 1)
        freqs = [value for(name, value) in enumerate(pivot)]
        probs = freqs ./ length(df_org.lidar)
        heatmap(x, y, probs)

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/group_hour_heatmap/group_hour_heatmap.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
    abstract type Formatted

Wrap a resource (e.g. filename or IO stream), adding format/coding information.

`Formatted` objects can be created with `infer` or `specify` and can be used
with `open`, `close`, `read` and `write`.
"""
abstract type Formatted end

function Base.show(io::IO, ::MIME"text/plain", f::Formatted)
	print(io, "$(typeof(f)): $(f.resource)")
	if f.format == nothing
		if length(f.format_guesses) == 1
			mime = f.format_guesses[1]
			print(io, "\n inferred format: $(mime)")
		elseif length(f.format_guesses) > 1
			print(io, "\n inferred format possibilities:")
			for mime in f.format_guesses
				print(io, "\n  $(mime)")
			end
		else
			print(io, "\n unknown format")
		end
		if length(f.coding_guesses) == 1
			mime = f.coding_guesses[1]
			print(io, "\n inferred coding: $(mime)")
		elseif length(f.coding_guesses) > 1
			print(io, "\n inferred coding possibilities:")
			for mime in f.coding_guesses
				print(io, "\n  $(mime)")
			end
		end
	else
		print(io, "\n specified format: $(f.format)")
		if f.coding != nothing
			print(io, "\n specified coding: $(f.coding)")
		end
	end
end

"""
    resource(f::Formatted) -> resource

Return the resource (e.g. filename or IO stream) wrapped by `f`.
"""
resource(f::Formatted) = f.resource

"""
	struct FormattedFilename <: Formatted

	FormattedFilename(filename::AbstractString) -> formatted

Wrap a `filename`, adding format/coding information.

See also: `Formatted`.
"""
struct FormattedFilename <: Formatted
	resource::String
	format::Union{MIME,Nothing}
	coding::Union{MIME,Nothing}
	format_guesses::Vector{MIME}
	coding_guesses::Vector{MIME}
end

"""
	struct FormattedIO <: Formatted

	FormattedIO(io::IO) -> formatted

Wrap IO stream `io`, adding format/coding information

See also: `Formatted`.
"""
struct FormattedIO <: Formatted
	resource::IO
	format::Union{MIME,Nothing}
	coding::Union{MIME,Nothing}
	format_guesses::Vector{MIME}
	coding_guesses::Vector{MIME}
end

"""
	infer(filename|io|formatted) -> ::Formatted

Guess the format/coding of a resource (filename, IO stream, `Formatted` object).

# Arguments

	filename::AbstractString
	io::IO
	formatted::Formatted

# Examples

Guess format from filename extension:

```julia
julia> f = infer("insulin.pdb")
FormattedFilename: insulin.pdb
 inferred format: structure/x-pdb
```

Guess format from an IO stream:

```julia
julia> io = open("insulin.pdb")
IOStream(<file insulin.pdb>)

julia> f = infer(io)
FormattedIO: IOStream(<file insulin.pdb>)
 inferred format: structure/x-pdb
```

Guess format from a `Formatted` object:

```julia
julia> f1 = infer("insulin.dat")
FormattedFilename: insulin.dat
 unknown format

julia> f2 = open(f1)
FormattedIO: IOStream(<file insulin.dat>)
 unknown format

julia> f3 = infer(f2)
FormattedIO: IOStream(<file insulin.dat>)
 inferred format: structure/x-pdb
```
"""
infer(filename::AbstractString) =
		FormattedFilename(filename, nothing, nothing, infercontent(filename)...)

infer(io::IO) = FormattedIO(io, nothing, nothing, infercontent(io)...)

infer(f::T) where {T<:Formatted} =
        T(f.resource, nothing, nothing, infercontent(f.resource)...)

"""
	specify(filename|io|formatted, format, [coding]) -> ::Formatted

Specify the format/coding of a resource (filename, IO stream, Formatted object).

# Arguments

	filename::AbstractString
	io::IO
	formatted::Formatted
	format::AbstractString
	coding::AbstractString

# Examples

Specify a format for a filename:

```julia
julia> f = specify("insulin.dat", "structure/x-pdb")
FormattedFilename: insulin.dat
 specified format: structure/x-pdb
```

Specify a format and coding for an IO stream:

```julia
julia> io = open("myoglobin.gro.gz")
IOStream(<file myoglobin.gro.gz>)

julia> f = specify(io, "structure/x-gro", "application/gzip")
FormattedIO: IOStream(<file myoglobin.gro.gz>)
 specified format: structure/x-gro
 specified coding: application/gzip
```

Specify a format for an existing Formatted object:

```julia
julia> f1 = infer("lysozyme.dat")
FormattedFilename: lysozyme.dat
 unknown format

julia> f2 = specify(f1, "structure/x-gro")
FormattedFilename: lysozyme.dat
 specified format: structure/x-gro
```
"""
specify(filename::AbstractString, format::AbstractString,
        coding::AbstractString) = FormattedFilename(filename,
        MIME{Symbol(format)}(), MIME{Symbol(coding)}(), [], [])

specify(filename::AbstractString, format::AbstractString) =
		FormattedFilename(filename, MIME{Symbol(format)}(), nothing, [], [])

specify(io::IO, format::AbstractString, coding::AbstractString) =
		FormattedIO(io, MIME{Symbol(format)}(), MIME{Symbol(coding)}(), [], [])

specify(io::IO, format::AbstractString) =
		FormattedIO(io, MIME{Symbol(format)}(), nothing, [], [])

specify(f::T, format::AbstractString, coding::AbstractString) where
        {T<:Formatted} =
        T(f.resource, MIME{Symbol(format)}(), MIME{Symbol(coding)}(), [], [])

specify(f::T, format::AbstractString) where {T<:Formatted} =
		T(f.resource, MIME{Symbol(format)}(), nothing, [], [])

"""
    formatted(filename|io|formatted) -> ::Formatted

Wrap a resource, preserving existing format/coding information, and inferring
format/coding when there is no such information.

# Arguments

	filename::AbstractString
	io::IO
	formatted::Formatted
"""
formatted(filename::AbstractString) = infer(filename)

formatted(io::IO) = infer(io)

formatted(f::Formatted) = f

"""
    isspecified(f::Formatted) -> ::Bool

Check if the format of `f` was explicitly specified.
"""
isspecified(f::Formatted) = (f.format != nothing)

"""
    isinferred(f::Formatted) -> ::Bool

Check if the format of `f` was inferred from a filename extension or stream
signature.
"""
isinferred(f::Formatted) =
        (f.format == nothing) && (length(f.format_guesses) > 0)

"""
    isunknown(f::Formatted) -> ::Bool

Check if the format of `f` was neither specified nor successfully inferred.
"""
isunknown(f::Formatted) =
        (f.format == nothing) && (length(f.format_guesses) == 0)

"""
    isambiguous(f::Formatted) -> ::Bool

Check if several possible formats were inferred for `f`.
"""
isambiguous(f::Formatted) =
        (f.format == nothing) && (length(f.format_guesses) > 1)

function resolveformat(f::Formatted)
	format = f.format
	if format == nothing
		if length(f.format_guesses) == 0
			error("$(f.resource): could not determine format")
		else
			format = f.format_guesses[1]
			if length(f.format_guesses) > 1
				@warn("$(f.resource): ambiguous format, assuming $(format)")
			end
		end
	end
	format
end

"""
	getformat(f::Formatted) -> ::String

Resolve the format associated with the resource wrapped by `f`.

If the format was not specified and could not be inferred, an error is thrown.
If the format was not specified and several possible formats were inferred, the
first guess is returned, and a warning is recorded.
"""
getformat(f::Formatted) = string(resolveformat(f))

function resolvecoding(f::Formatted)
	coding = f.coding
	if f.format == nothing
		if length(f.coding_guesses) != 0
			coding = f.coding_guesses[1]
			if length(f.coding_guesses) > 1
				@warn("$(f.resource): ambiguous coding, assuming $(coding)")
			end
		end
	end
    coding
end

"""
	getcoding(f::Formatted) -> ::String

Resolve the coding associated with the resource wrapped by `f`.

If the coding was not specified and could not be inferred, no coding is assumed
and `nothing` is returned. If the coding was not specified and several possible
codings were inferred, the first guess is returned, and a warning is recorded.
"""
function getcoding(f::Formatted)
    coding = resolvecoding(f)
    if coding == nothing
        nothing
    else
        string(coding)
    end
end


text:
using Test
using BenchmarkTools

using iLQGames:
    iLQGames,
    iLQSolver,
    GeneralGame,
    Unicycle4D,
    NPlayerUnicycleCost,
    AffineStrategy,
    SystemTrajectory,
    dynamics,
    n_states,
    n_controls,
    horizon,
    player_costs,
    quadraticize,
    _quadraticize_ad,
    generate_nplayer_navigation_game,
    solve!,
    ProximityCost

using iLQGames.TestUtils
using StaticArrays
using LinearAlgebra

# generate a game
T_horizon = 10.
ΔT = 0.1

"--------------------------------- Unicycle4D ---------------------------------"

x01 = @SVector [-3., 0., 0., 0.]
x02 = @SVector [0.,  3., -pi/2, 0.]
x03 = @SVector [-3.,  3., -pi/4, 0.]
x0 = vcat(x01, x02, x03)
# goal states (goal position of other player with opposite orientation)
xg1 = @SVector [3., 0., 0., 0.]
xg2 = @SVector [0., -3., -pi/2, 0.]
xg3 = @SVector [3., -3., -pi/4, 0.]
g = generate_nplayer_navigation_game(Unicycle4D, NPlayerUnicycleCost, T_horizon,
                                     ΔT, xg1, xg2, xg3;
                                     proximitycost=ProximityCost([2.0, 2.0, 2.0],
                                                                 [0., 50.0, 50.0]))
dyn = dynamics(g)
nx = n_states(dyn)
nu = n_controls(dyn)
pcs = player_costs(g)
h = horizon(g)
zero_op = zero(SystemTrajectory{h, ΔT, nx, nu})

quad_sanity_check(g)

# solve the lq game
solver = iLQSolver(g; state_regularization=5.0, control_regularization=5.0)
# - setup initial_strategy
steer_init(k::Int) = cos(k/h*pi) * deg2rad(0)
acc_init(k::Int) = -cos(k/h*pi)*0.1
γ_init = SizedVector{h}([AffineStrategy((@SMatrix zeros(nu, nx)),
                                 (@SVector [steer_init(k), 0.7*acc_init(k),
                                            steer_init(k), acc_init(k),
                                            steer_init(k), acc_init(k)])) for k in 1:h])
# generate initial operating point from simulating initial strategy
# solve the game
@benchmark(solve!(copy(zero_op), copy(γ_init), $g, $solver, $x0)) |> display


text:
#
# script to run for the first time
#
# see https://github.com/JuliaLang/IJulia.jl
#

#
# packages to build
#
Pkg.init()

max_stars_repo_name: sdpython/pymyinstall

max_issues_repo_name: sdpython/pymyinstall

max_forks_repo_name: sdpython/pymyinstall


text:
# Letter                           Value
# A, E, I, O, U, L, N, R, S, T       1
# D, G                               2
# B, C, M, P                         3
# F, H, V, W, Y                      4
# K                                  5
# J, X                               8
# Q, Z                               10

function score(ch::Char)
    ch = uppercase(ch)
    if ! (Int('A') <= Int(ch) <= Int('Z'))
        return 0
    end

    for (chars, value) in [
        "DG"    => 2,
        "BCMP"  => 3,
        "FHVWY" => 4,
        "K"     => 5,
        "JX"    => 8,
        "QZ"    => 10,
    ]
        if ch in chars
            return value
        end
    end

    return 1
end

function score(str::AbstractString)
    sum(map(score, collect(str)))
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Somoclu
using Base.Test

import Somoclu: distance

# This is an Euclidean distance estimate computed from Julia,
function distance(p1::Ptr{Cfloat}, p2::Ptr{Cfloat}, d::Cuint)
    s = Cfloat(0.0)
    for i = 1:d
        v = unsafe_load(p1, i) - unsafe_load(p2, i)
        s += v*v
    end
    (isinf(s) || isnan(s)) && throw(ErrorException("Invalid distance computed."))
    return sqrt(s)::Cfloat
end

function deterministic_codebook(useCustomDistance=false)
    ncolumns, nrows = 2, 2
    initialcodebook = Array{Float32}(zeros(2, ncolumns*nrows))
    som = Som(ncolumns, nrows, initialcodebook=initialcodebook,
              useCustomDistance=useCustomDistance)
    println("useCustomDistance: $(som.useCustomDistance)")
    data = Array{Float32}([0.1 0.2; 0.3 0.4]);
    train!(som, data);
    correct_codebook = Array{Float32, 2}([0.15  0.126894  0.173106  0.15;
                                          0.35  0.326894  0.373106  0.35]);
    return sum(som.codebook - correct_codebook) < 10e-6
end

@test deterministic_codebook()
@test deterministic_codebook(true)


text:
__precompile__()

"""
The Search optimization algorithm.
"""
module MendelSearch
#
# Required external modules.
#
using LinearAlgebra
using Printf

export mendel_search
#
# SearchSetup sets up the structures needed for the MendelSearch routines,
# including the Parameter data structure whose entries control optimization.
#
include("SearchSetup.jl")

"""
Computes the values of a function f over a user defined grid of
points or the constrained optimum of f by the method of recursive
quadratic programming. The hessian d2f of f is approximated by Davidon's
secant update unless an exact value is supplied by the user. Users must
supply the objective function fun. Parameter initial values and
constraints are passed through the parameter data structure.
"""
function mendel_search(fun::Function, parameter::Parameter)
  #
  # Set variables controlling convergence and numerical differentiation.
  #
  conv_crit = 1e-6
  conv_tests = 4
  max_iters = 1000
  max_steps = 3
  max_step_length = Inf
  dp = 1e-7
  small = 1e-5
  #
  # Write a header.
  #
  cases = parameter.cases
  io = parameter.output_unit
  travel = parameter.travel
  goal = parameter.goal
  println(io, " \n                Search, Julia Version\n")
  println(io, "            (c) Copyright Kenneth Lange, 2019\n")
  if parameter.title != ""
    println(io, "Title: " * parameter.title)
  end
  println(io, "Grid or search option: " * travel)
  println(io, "Minimize or maximize: " * goal)
  #
  # Abbreviate the names of passed variables.
  #
  constraints = parameter.constraints
  pars = parameter.parameters
  points = parameter.points
  travel = parameter.travel
  pname = parameter.name
  par = parameter.par
  pmin = parameter.min
  pmax = parameter.max
  constraint = parameter.constraint
  level = parameter.constraint_level
  function_value = parameter.function_value
  #
  # Create variables required by searching.
  #
  grad_provided = false
  hess_provided = false
  f = 0.0
  df = zeros(pars)
  grad = zeros(pars)
  d2f = zeros(pars, pars)
  hess = zeros(pars, pars)
  df_old = zeros(pars)
  par_old = zeros(pars)
  #
  # Check for bound and constraint errors.
  #
  error = check_constraints(constraint, level, pname, par, pmin, pmax,
    io, travel)
  if error; return (par, 0.0); end
  #
  # Compute function values over a user defined grid of points.
  # Keep track of the optimal point and optimal function value.
  #
  fmin = Inf
  iter_min = 0
  best_point = zeros(pars)
  if travel == "grid"
    for iteration = 1:points
      parameter.par = vec(parameter.grid[iteration, :])
      (f, df, d2f) = fun(parameter.par)
      iteration_output!(io, f, iteration, parameter.par, pname)
      function_value[iteration] = f
      if goal == "maximize"
        f = - f
      end
      if f < fmin
        fmin = f
        iter_min = iteration
        copyto!(best_point, par)
      end
    end
    #
    # Return the optimal point and optimal function value.
    #
    if goal == "minimize"
      println(io, " \nThe minimum function value of ",
        round(fmin; digits = 5)," occurs at iteration ",iter_min, ".")
      return (best_point, fmin)
    else
      println(io, " \nThe maximum function value of ",
        round(-fmin; digits = 5)," occurs at iteration ",iter_min, ".")
      return (best_point, - fmin)
    end
    #
    # Otherwise minimize the objective function. df and d2f
    # are the first two differentials. ker is the kernel of
    # the constraint matrix.
    #
  else
    df = zeros(pars)
    d2f = Matrix{Float64}(I, pars, pars)
    (matrix_u, d, matrix_v) = svd(constraint, full = true)
    ker = matrix_v[:, constraints + 1:pars]
    conv_count = 0
    forward = true
    #
    # Fetch the first function value and the two differentials.
    # Output the first iteration.
    #
    iteration = 1
    (f, grad_provided, grad, hess) = dfun(fun, parameter, dp, forward)
    iteration_output!(io, f, iteration, par, pname)
    function_value[1] = f
    #
    # Copy the user supplied first and second differentials.
    #
    hess_provided =  typeof(hess) != Nothing
    if goal == "minimize"
      copyto!(df, grad)
      if hess_provided; copyto!(d2f, hess); end
    else
      f = - f
      copyto!(df, - grad)
      if hess_provided; copyto!(d2f, - hess); end
    end
    #
    # Preserve the current point.
    #
    copyto!(best_point, par)
    fmin = f
    iter_min = 1
    #
    # Enter the main iteration loop.
    #
    for iteration = 2:max_iters
    #
    # Solve the quadratic programming problem. If cycles = 0,
    # then the quadratic programming algorithm has failed. Reset
    # the Hessian approximation to the identity and try again.
    #
      for j = 1:2
        tableau = create_tableau(d2f, df, constraint, level, par)
        (cycles, delta) =
        quadratic_program(tableau, par, pmin, pmax, pars, constraints)
        #
        # Check whether quadratic programming succeeded.
        #
        if cycles > 0
          break
        else
          c = Inf
          for i = 1:pars
            if d2f[i, i] > 0.0; c = min(c, d2f[i, i]); end
          end
          #
          # In an emergency revert to steepest descent.
          #
          if c >= Inf
            d2f = Matrix{Float64}(I, pars, pars)
          #
          # Otherwise, bump up the diagonal entries of the hessian.
          #
          else
            for i = 1:pars
              if d2f[i, i] <= c; d2f[i, i] = c * (1.0 + rand()); end
            end
          end
        end
      end
      #
      # Compute a new function value. If it represents a sufficient
      # decrease, then accept the new point. Otherwise backtrack by
      # fitting a quadratic through the two current points with the
      # given slope at the first point. Quit, if necessary, after the
      # maximum number of steps.
      #
      l2_norm = norm(delta)
      if l2_norm > max_step_length
        delta = max_step_length * delta / l2_norm
      end
      d = min(dot(df, delta), 0.0)
      copyto!(df_old, df)
      copyto!(par_old, par)
      f_old = f
      t = 1.0
      steps = -1
      for step = 0:max_steps
        steps = steps + 1
        #
        # Fetch another function value, and if provided, the user
        # supplied second differential.
        #
        par = par_old + t * delta
        parameter.par = par
        (f, grad_provided, grad, hess) = dfun(fun, parameter, dp, forward)
        if goal == "minimize"
          copyto!(df, grad)
          if hess_provided; copyto!(d2f, hess); end
        else
          f = - f
          copyto!(df, - grad)
          if hess_provided; copyto!(d2f, - hess); end
        end
        #
        # Test for a sufficient decrease in the objective.
        #
        if f <= f_old + t * d / 10.0 || step == max_steps
          break
        else
          t = max( -d * t * t / (2.0 * (f - f_old - t * d)), t / 10.0)
        end
      end
      #
      # Output the current iteration, and update the best point.
      #
      if goal == "minimize"
        iteration_output!(io, f, iteration, par, pname, steps)
      else
        iteration_output!(io, - f, iteration, par, pname, steps)
      end
      if f < fmin
        fmin = f
        iter_min = iteration
        copyto!(best_point, par)
      end
      #
      # Check the convergence criterion. If it has been satisfied a
      # sufficient number of times, then exit the main loop. Otherwise,
      # output the current iteration and continue the search.
      #
      if abs(f_old - f) > conv_crit
        forward = true
        conv_count = 0
      else
        forward = false
        conv_count = conv_count + 1
      end
      if conv_count >= conv_tests
        break
      end
      #
      # Update the approximate Hessian by Davidon's secant method.
      #
      if !hess_provided
        delta = par - par_old
        u = df - df_old - d2f * delta
        c = dot(u, delta)
        if c^2 > (1e-8) * dot(u, u) * dot(delta, delta)
          c = -1.0 / c
          #
          # Ensure that the updated Hessian is positive definite on
          # the kernel of the constraint matrix.
          #
          matrix = inv(ker' * d2f * ker)
          v = ker' * u
          d = dot(v, matrix * v)
          c = min(c, 0.8 / d)
          d2f = d2f - c * u * u'
        end
      end
    end
  end
  #
  # Output the optimal function value.
  #
  if goal == "minimize"
    println(io, " \nThe minimum function value of ",
      round(fmin; digits = 5)," occurs at iteration ", iter_min, ".")
    function_value[2] = fmin
  else
    println(io, " \nThe maximum function value of ",
      round(- fmin; digits = 5)," occurs at iteration ", iter_min, ".")
    function_value[2] = - fmin
  end
  #
  # If the asymptotic covariance matrix is desired, then
  # record which parameters occur on a boundary
  #
  boundary = falses(pars)
  if parameter.standard_errors
    for i = 1:pars
      boundary[i] = (par[i] <= pmin[i] + small || par[i] >= pmax[i] - small)
    end
    boundaries = sum(boundary)
    #
    # If the exact first differential is available, approximate the
    # columns of the second differential by forward differences of
    # the first differential.
    #
    if !hess_provided
      if grad_provided
        for i = 1:pars
          if boundary[i]; continue; end
          d = dp * max(abs(par[i]), 1.0)
          if par[i] + d >= pmax[i]; d = -d; end
          parameter.par[i] = par[i] + d
          (g, grad_provided, grad, hess) = dfun(fun, parameter, dp, forward)
          if goal == "maximize"; g = - g; end
          parameter.par[i] = parameter.par[i] - d
          d2f[1:i, i] = (grad[1:i] - df[1:i]) / d
          d2f[i, 1:i] = d2f[1:i, i]
        end
      #
      # Otherwise, use the second difference formula of Dennis and Schnabel
      # to approximate the entries of the second differential. See page
      # on page 321 of their book.
      #
      else
        #
        # If necessary, recompute the first differential using central differences.
        #
        if forward
          forward = false
          parameter.par = par
          (f, grad_provided, df, hess) = dfun(fun, parameter, dp, forward)
        end
        #
        # Compute and store function values and parameter increments.
        #
        u = zeros(pars)
        delta = zeros(pars)
        dp23 = cbrt(dp)^2
        for i = 1:pars
          if boundary[i]; continue; end
          d = dp23 * max(abs(par[i]), 1.0)
          if parameter.par[i] + d >= pmax[i]; d = -d; end
          parameter.par[i] = parameter.par[i] + d
          (g, grad, hess) = fun(parameter.par)
          if goal == "maximize"; g = - g; end
          parameter.par[i] = parameter.par[i] - d
          delta[i] = d
          u[i] = g
        end
        #
        # Apply these to compute the entries of the second differential.
        #
        for i = 1:pars
          for j = 1:i
            if !boundary[i] && !boundary[j]
              parameter.par[i] = parameter.par[i] + delta[i]
              parameter.par[j] = parameter.par[j] + delta[j]
              (g, grad, hess) = fun(parameter.par)
              if goal == "maximize"; g = - g; end
              parameter.par[i] = parameter.par[i] - delta[i]
              parameter.par[j] = parameter.par[j] - delta[j]
              d2f[i, j] = (g - u[i] - u[j] + f) / (delta[i] * delta[j])
              d2f[j, i] = d2f[i, j]
            end
          end
        end
      end
    end
    #
    # Adjust the Hessian for a least squares problem by dividing
    # by the residual mean square.
    #
    reduced = pars - constraints - boundaries
    if cases > 0
      n = cases - reduced
      if n > 0 && f > 0.0
        sigma_sq = 2.0 * f / n
        d2f = d2f / sigma_sq
      else
        fill!(d2f, 0.0)
      end
    end
    #
    # After these preliminaries, compute the asymptotic covariances.
    #
    if reduced > 0
      asy_cov = asymptotic_covariances(io, constraint, boundary, d2f, pname)
    else
      write(io, " \n The asymptotic covariance matrix is undefined.")
      println(io, "")
    end
  end
  #
  # Return the optimal point and optimal function value.
  #
  if goal == "minimize"
    return (best_point, fmin)
  else
    return (best_point, - fmin)
  end
end # function mendel_search

"""
Calculates the tableau used in minimizing the quadratic
0.5 x' Q x + r' x, subject to Ax = b and parameter lower
and upper bounds.
"""
function create_tableau(matrix_q::Matrix{Float64}, r::Vector{Float64},
  matrix_a::Matrix{Float64}, b::Vector{Float64}, x::Vector{Float64})

  m = size(matrix_a, 1)
  #
  # Create the tableau in the absence of constraints.
  #
  if m == 0
    tableau = [matrix_q (-r); -r' 0]
  else
  #
  # In the presence of constraints compute a constant mu via
  # the Gerschgorin circle theorem so that Q + mu * A' * A
  # is positive definite.
  #
    (matrix_u, d, matrix_v) = svd(matrix_a, full = true)
    matrix_p = matrix_v * matrix_q * matrix_v'
    mu = 0.0
    for i = 1:m
      mu = max((norm(matrix_p[:, i], 1) - 2.0 * matrix_p[i, i]) / d[i]^2, mu)
    end
    mu = 2.0 * mu
    #
    # Now create the tableau.
    #
    tableau = [matrix_q + mu * matrix_a' * matrix_a matrix_a' (-r);
               matrix_a zeros(m, m) b - matrix_a * x;
               -r' (b - matrix_a * x)' 0]
  end
  return tableau
end # function create_tableau

"""
Solves the p-dimensional quadratic programming problem

 min [df * delta + 0.5 * delta' * d^2 f * delta]
 subject to: constraint * delta = 0 and pmin <= par + delta <= pmax.

See: Jennrich JI, Sampson PF (1978) "Some problems faced in making
a variance component algorithm into a general mixed model program."
Proceedings of the Eleventh Annual Symposium on the Interface.
Gallant AR, Gerig TM, editors. Institute of Statistics,
North Carolina State University.
"""
function quadratic_program(tableau::Matrix{Float64}, par::Vector{Float64},
  pmin::Vector{Float64}, pmax::Vector{Float64}, p::Int, c::Int)

  delta = zeros(Float64, size(par))
  #
  # See function create_tableau for the construction of the tableau.
  # For checking tolerance, set diag to the diagonal elements of tableau.
  # Begin by sweeping on those diagonal elements of tableau corresponding
  # to the parameters. Then sweep on the diagonal elements corresponding
  # to the constraints. If any parameter fails the tolerance test, then
  # return and reset the approximate Hessian.
  #
  small = 1e-5
  tol = 1e-8
  d = diag(tableau)
  for i = 1:p
    if d[i] <= 0.0 || tableau[i, i] < d[i] * tol
      return (0, delta)
    else
      sweep!(tableau, i, false)
    end
  end
  swept = trues(p)
  for i = p + 1:p + c
    if tableau[i, i] >= 0.0
      return (0, delta)
    else
      sweep!(tableau, i, false)
    end
  end
  #
  # Take a step in the direction tableau(i, end) for the parameters i
  # that are currently swept. If a boundary is encountered, determine
  # the maximal fractional step possible.
  #
  cycle_main_loop = false
  for iteration = 1:1000
    a = 1.0
    for i = 1:p
      if swept[i]
        ui = tableau[i, end]
        if ui > 0.0
          ai = pmax[i] - par[i] - delta[i]
        else
          ai = pmin[i] - par[i] - delta[i]
        end
        if abs(ui) > 1e-10
          a = min(a, ai / ui)
        end
      end
    end
    #
    # Take the fractional step for the currently swept parameters, and
    # reset the transformed partial derivatives for these parameters.
    #
    for i = 1:p
      if swept[i]
        ui = tableau[i, end]
        delta[i] = delta[i] + a * ui
        tableau[i, end] = (1.0 - a) * ui
        tableau[end, i] = tableau[i, end]
      end
    end
    #
    # Find a swept parameter that is critical, and inverse sweep it.
    # Go back and try to take another step or fractional step.
    #
    cycle_main_loop = false
    for i = 1:p
      critical = pmin[i] >= par[i] + delta[i] - small
      critical = critical || pmax[i]<= par[i] + delta[i] + small
      if swept[i] && abs(tableau[i, i])>1e-10 && critical
        sweep!(tableau, i, true)
        swept[i] = false
        cycle_main_loop = true
        break
      end
    end
    if cycle_main_loop; continue; end
    #
    # Find an unswept parameter that violates the KKT condition
    # and sweep it. Go back and try to take a step or fractional step.
    # If no such parameter exists, then the problem is solved.
    #
    for i = 1:p
      ui = tableau[i, end]
      violation = ui > 0.0 && pmin[i] >= par[i] + delta[i] - small
      violation = violation || (ui<0.0 && pmax[i]<= par[i] + delta[i] + small)
      if !swept[i] && violation
        sweep!(tableau, i, false)
        swept[i] = true
        cycle_main_loop = true
        break
      end
    end
    if cycle_main_loop; continue; end
    return (iteration, delta)
  end
  return (0, delta)
end # function quadratic_program

"""
Sweeps or inverse sweeps the symmetric tableau A on its kth diagonal entry.
"""
function sweep!(matrix_a::Matrix{Float64}, k::Int, inverse::Bool = false)

  p = 1.0 / matrix_a[k, k]
  v = matrix_a[:, k]
  matrix_a[:, k] .= 0.0
  matrix_a[k, :] .= 0.0
  if inverse
    v[k] = 1.0
  else
    v[k] = -1.0
  end
  for i = 1:size(matrix_a, 1)
    pv = p * v[i]
    matrix_a[:, i] = matrix_a[:, i] - pv * v
  end
end # function sweep!

"""
Controls the computation of f and its first two differentials.
When the first partials are computed numerically, forward or central
differences are used depending on the logical variable forward. The
numerical differentiation interval is adjusted to take into account
the magnitude of a parameter and whether it lies on its upper bound.
"""
function dfun(fun::Function, parameter::Parameter, dp::Float64, forward::Bool)

  par = parameter.par
  pmax = parameter.max
  (f, df, d2f) = fun(par)
  grad_provided = typeof(df) != Nothing
  if !grad_provided
    df = zeros(length(par))
    for i = 1:length(par)
      d = dp * max(abs(par[i]), 1.0)
      psave = par[i]
      if forward
        if psave + d > pmax[i]; d = -d; end
        par[i] = par[i] + d
        (fplus, grad, hess) = fun(par)
        df[i] = (fplus - f) / d
      else
        par[i] = par[i] + d
        (fplus, grad, hess) = fun(par)
        par[i] = psave - d
        (fminus, grad, hess) = fun(par)
        df[i] = (fplus - fminus) / (d + d)
      end
      par[i] = psave
    end
  end
  return (f, grad_provided, df, d2f)
end # function dfun

"""
Checks for constraint violations. The returned value conveys the
nature of the violation.
"""
function check_constraints(constraint::Matrix{Float64},
  level::Vector{Float64}, pname::Vector{AbstractString}, par::Vector{Float64},
  pmin::Vector{Float64}, pmax::Vector{Float64}, io::IO, travel::AbstractString)

  tol = 1e-4
  error = false
  (constraints, pars) = size(constraint)
  #
  # Check for bound violations.
  #
  for i = 1:pars
    if par[i] < pmin[i]
      println(io,
        " \nError: Parameter ", i, " is less than its minimum.")
      error = true
    elseif par[i] > pmax[i]
      println(io,
        " \nError: Parameter ", i, " is greater than its maximum.")
      error = true
    end
    if pmin[i] > pmax[i] - tol
      println(io, " \nError: The bounds on parameter ", i,
                  " are too close or inconsistent.")
      error = true
    end
  end
  #
  # Check for constraint violations.
  #
  if constraints > 0
    y = constraint * par
    for i = 1:constraints
      if abs(y[i] - level[i]) > tol
        println(io,
          " \nError: Linear equality constraint ", i, " is not satisfied.")
        error = true
      end
      if count(!iszero, constraint[i, :]) != 1
        continue
      else
        for j = 1:pars
          if constraint[i, j] != 0.0 &&
          (abs(par[j] - pmin[j]) < tol || abs(par[j] - pmax[j]) < tol)
            println(io,
            " \nError: parameter ", j, " is constrained to one of its bounds.")
            error = true
          end
        end
      end
    end
    if rank(constraint) < size(constraint, 1)
      println(io, " \nError: Some equality constraints are redundant.")
      error = true
    end
  end
  #
  # Echo the bounds and the constraints.
  #
  if travel == "search"
    println(io, " \nParameter minima and maxima:")
    name = join(pname, "    ")
    println(io, " \n    ", name)
    if pmin[1] == -Inf
      @printf(io, " \n    -Inf     ")
    else
      @printf(io, " \n%12.4e", pmin[1])
    end
    for i = 2:pars
      if pmin[i] == -Inf
        @printf(io, "   -Inf     ")
      else
        @printf(io, "%12.4e", pmin[i])
      end
    end
    if pmax[1] == Inf
      @printf(io, " \n     Inf     ")
    else
      @printf(io, " \n%12.4e", pmax[1])
    end
    for i = 2:pars
      if pmax[i] == Inf
        @printf(io, "    Inf     ")
      else
        @printf(io, "%12.4e", pmax[i])
      end
    end
    println(io, " ")
    if constraints > 0
      println(io, " \nParameter constraints:")
      println(io, " \n    ", name, "    level ")
      println(io, "")
      for i = 1:constraints
        for j = 1:pars
          @printf(io, "%12.4e", constraint[i, j])
        end
        @printf(io, "%12.4e\n", level[i])
      end
    end
  end
  return error
end # function check_constraints

"""
Outputs the current iteration.
"""
function iteration_output!(io::IO, f::Float64, iteration::Int,
  par::Vector{Float64}, pname::Vector{AbstractString}, steps::Int = 0)

  if iteration == 1
    name = join(pname, "    ")
    println(io, " \niter  steps   function    ", name)
  end
  @printf(io, " \n %2i   %2i   %11.4e ", iteration, steps, f)
  for i = 1:length(par)
     @printf(io, "%12.4e", par[i])
  end
  println(io, " ")
end # function iteration_output!

"""
Computes the asymptotic covariance matrix of the
parameter estimates.
"""
function asymptotic_covariances(io::IO, constraint::Matrix{Float64},
  boundary::BitArray{1}, d2f::Matrix{Float64}, pname::Vector{AbstractString})
  #
  # Reparameterize to find the asymptotic covariance matrix.
  # Add an extra constraint for each parameter occurring on a boundary.
  #
  (constraints, pars) = size(constraint)
  asy_cov = zeros(pars, pars)
  asy_cov[1:constraints, :] = constraint
  j = constraints + 1
  for i = 1:pars
    if boundary[i]
      asy_cov[j, i] = 0.0
      j = j + 1
    end
  end
  (matrix_u, d, matrix_v) = svd(asy_cov, full = true)
  ker = matrix_v[:, j:pars]
  matrix = ker' * d2f * ker
  asy_cov = ker * inv(matrix) * ker'
  #
  # Standardize the asymptotic covariance matrix.
  #
  for i = 1:pars
    if asy_cov[i, i] <= 1e-10
      asy_cov[i, i] = 0.0
    else
      asy_cov[i, i] = sqrt(asy_cov[i, i])
    end
    for j = 1:i - 1
      if asy_cov[i, i] > 1e-5 &&  asy_cov[j, j] > 1e-5
        asy_cov[i, j] = asy_cov[i, j] / (asy_cov[i, i] * asy_cov[j, j])
      else
        asy_cov[i, j] = 0.0
        asy_cov[j, i] = 0.0
      end
    end
  end
  #
  # Output the results.
  #
  println(io, " \nThe asymptotic standard errors of the parameters:")
  name = join(pname, "    ")
  println(io, " \n    ", name)
  println(io, " ")
  for i = 1:pars
    @printf(io, "%12.4e", asy_cov[i, i])
  end
  println(io, "")
  println(io, " \nThe asymptotic correlation matrix of the parameters:")
  println(io, " \n    ", name)
  println(io, "")
  @printf(io, "%10.4f\n", 1.0)
  for i = 2:pars
    println(io, "")
    @printf(io, "%10.4f", asy_cov[i, 1])
    for j = 2:i - 1
      @printf(io, "%12.4f", asy_cov[i, j])
    end
    if asy_cov[i, i] > 1e-5
      @printf(io, "%12.4f", 1.0)
    else
      @printf(io, "%12.4f", 0.0)
    end
    println(io, "")
  end
end # function asymptotic_covariances
#
# Method to obtain path to this package's data files
# so they can be used in the documentation and testing routines.
# For example, datadir("Control file.txt") will return
# "/path/to/package/data/Control file.txt"
#
datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)

end # module MendelSearch

max_stars_repo_name: OpenMendel/Search

max_issues_repo_name: OpenMendel/Search

max_forks_repo_name: OpenMendel/Search


text:
centered(C::Type{HyperCube{N,T}}) where {N,T} = C(Vec{N,T}(-0.5), T(1))
function centered(::Type{T}) where T <: HyperCube
    centered(HyperCube{ndims_or(T, 3), eltype_or(T, Float32)})
end

centered(R::Type{HyperRectangle{N,T}}) where {N, T} = R(Vec{N,T}(-0.5), Vec{N,T}(1))
centered(::Type{T}) where {T <: HyperRectangle} = centered(HyperRectangle{ndims_or(T, 3), eltype_or(T, Float32)})

centered(S::Type{HyperSphere{N, T}}) where {N, T} = S(Vec{N,T}(0), T(0.5))
centered(::Type{T}) where {T <: HyperSphere} = centered(HyperSphere{ndims_or(T, 3), eltype_or(T, Float32)})


text:
module runTests
sep_rank_proj_path = dirname(dirname(@__FILE__))

testDir = dirname(@__FILE__)
include(testDir*"\\testUtils_states.jl")
include(testDir*"\\testExamples.jl")
include(testDir*"\\testMoments.jl")
# include(testDir*"\\testSep_Constriants.jl")
include(testDir*"\\testSep_Model.jl")
include(testDir*"\\testSep_Compute.jl")

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
using Dates
using JSON3
using Test
using Logging

using Stonks
using Stonks: JSONContent, APIResponseError, ContentParserError
using Stonks.Models
using Stonks.Parsers

include("test_utils.jl")

@testset "Yahoofinance parsers" begin
  price_content = get_test_data("data/yahoo_prices.json")
  overview_content = get_test_data("data/yahoo_overview.json")
  exchange_content = get_test_data("data/yahoo_exchange.json")
  is_content = get_test_data("data/yahoo_income_statement.json")
  bs_content = get_test_data("data/yahoo_balance_sheet.json")
  cf_content = get_test_data("data/yahoo_cashflow.json")
  eps_content = get_test_data("data/yahoo_earnings.json")
  price_parser = Parsers.YahooPriceParser
  exchange_parser = Parsers.YahooExchangeRateParser
  is_parser = Parsers.YahooIncomeStatementParser
  bs_parser = Parsers.YahooBalanceSheetParser
  cf_parser = Parsers.YahooCashflowParser
  eps_parser = Parsers.YahooEarningsParser

  @testset "Unpack quote summary response" begin
    content = """{
      "quoteSummary": {
        "result": null,
        "error": {
          "code": "Not Found",
          "description": "No fundamentals data found for any of the summaryTypes=incomeStatementHistory"
        }
      }
    }"""
    res = Parsers.unpack_quote_summary_response(JSON3.read(content))
    @test isa(res, APIResponseError)
    content = """{"foo": "bar"}"""
    res = Parsers.unpack_quote_summary_response(JSON3.read(content))
    @test isa(res, ContentParserError)
    content = """{
      "quoteSummary": {
        "result": [{"foo": "bar"}],
        "error": null
      }
    }"""
    res = Parsers.unpack_quote_summary_response(JSON3.read(content))
    @test isa(res, JSONContent)
  end

  @testset "Succesful price response" begin
    data = parse_content(price_parser, price_content)
    @test isa(data, Vector{AssetPrice})
    # 5 datapoint for 4 tickers
    @test length(data) == 5 * 4
    dates = [x.date for x in data]
    @test minimum(dates) == Date("2022-02-09")
    @test maximum(dates) == Date("2022-02-15")
  end

  @testset "Successful price response with from date limit" begin
    from = Date("2022-02-11")
    data = parse_content(price_parser, price_content; from=from)
    @test isa(data, Vector{AssetPrice})
    # 3 datapoint for 4 tickers (filtered using from)
    @test length(data) == 3 * 4
    dates = [x.date for x in data]
    @test minimum(dates) == from
    @test maximum(dates) == Date("2022-02-15")
  end

  @testset "Succesful price response with from and to date limits" begin
    from, to = Date("2022-02-11"), Date("2022-02-14")
    data = parse_content(price_parser, price_content; from=from, to=to)
    @test isa(data, Vector{AssetPrice})
    # 3 datapoint for 4 tickers (filtered using from)
    @test length(data) == 2 * 4
    dates = [x.date for x in data]
    @test minimum(dates) == from
    @test maximum(dates) == to
  end

  @testset "Failed yahoo API response" begin
    content = """{
      "spark": {
        "result": null,
        "error": {
          "code": "Not Found",
          "description": "No data found for spark symbols"
        }
      }
    }"""
    data = parse_content(price_parser, content)
    @test isa(data, APIResponseError)
  end

  @testset "Succesful yahoofinance response on AssetInfo" begin
    parser = Parsers.YahooInfoParser
    data = parse_content(parser, overview_content)
    @test isa(data, Vector{AssetInfo})
    @test length(data) == 1
  end

  @testset "Succesful exchange rate response" begin
    data = parse_content(exchange_parser, exchange_content)
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2022-01-24")
    @test maximum(dates) == Date("2022-02-22")
    @test unique(map(x -> (x.base, x.target), data)) == [("USD", "CAD"), ("EUR", "USD")]
    @test isa(data, Vector{ExchangeRate})
    @test length(data) == 44
  end

  @testset "Succesful exchange rate response with from date filter" begin
    from = Date("2022-02-01")
    data = parse_content(exchange_parser, exchange_content; from=from)
    @test minimum(map(x -> x.date, data)) == from
  end

  @testset "Succesful exchange rate response with from and to date filters" begin
    from, to = Date("2022-02-01"), Date("2022-02-10")
    data = parse_content(exchange_parser, exchange_content; from=from, to=to)
    dates = map(x -> x.date, data)
    @test minimum(dates) == from
    @test maximum(dates) == to
  end

  @testset "Succesful income statement response" begin
    data = parse_content(is_parser, is_content; symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2018-12-31")
    @test maximum(dates) == Date("2021-12-31")
    @test isa(data, Vector{IncomeStatement})
  end

  @testset "Succesful income statement response with from date filter" begin
    data = parse_content(is_parser, is_content; from=Date("2020-01-01"), symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) >= Date("2020-01-01")
  end

  @testset "Failed income statement response" begin
    content = """{
      "quoteSummary": {
        "result": null,
        "error": {
          "code": "Not Found",
          "description": "No fundamentals data found for any of the summaryTypes=incomeStatementHistory"
        }
      }
    }"""
    res = parse_content(is_parser, content)
    @test isa(res, APIResponseError)
  end

  @testset "Succesful balance sheet response" begin
    data = parse_content(bs_parser, bs_content; symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2018-12-31")
    @test maximum(dates) == Date("2021-12-31")
    @test isa(data, Vector{BalanceSheet})
  end

  @testset "Succesful balance sheet response with from date filter" begin
    data = parse_content(bs_parser, bs_content; from=Date("2020-01-01"), symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) >= Date("2020-01-01")
  end

  @testset "Succesful earnings response" begin
    data = parse_content(eps_parser, eps_content; symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2021-03-31")
    @test maximum(dates) == Date("2021-12-31")
    @test isa(data, Vector{Earnings})
  end

  @testset "Succesful earnings response with from date filter" begin
    data = parse_content(eps_parser, eps_content; from=Date("2021-06-01"), symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) >= Date("2020-06-30")
  end

  @testset "Succesful earnings response with yearly frequency (no data)" begin
    data = parse_content(eps_parser, eps_content; symbol="IBM", frequency="yearly")
    @test isa(data, Vector{Earnings})
    @test length(data) == 0
  end

  @testset "Succesful cashflow statement response" begin
    data = parse_content(cf_parser, cf_content; symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2018-12-31")
    @test maximum(dates) == Date("2021-12-31")
    @test isa(data, Vector{Models.CashflowStatement})
  end

  @testset "Succesful cashflow statement response with from date filter" begin
    data = parse_content(cf_parser, cf_content; from=Date("2021-06-01"), symbol="IBM")
    dates = map(x -> x.date, data)
    @test minimum(dates) >= Date("2020-06-30")
  end

  @testset "Succesful cashflow statement response with yearly frequency" begin
    data = parse_content(cf_parser, cf_content; symbol="IBM", frequency="yearly")
    @test isa(data, Vector{Models.CashflowStatement})
    @test length(data) == 4
    @test unique(map(x -> x.frequency, data)) == ["yearly"]
  end

end


text:
"""

Read .csv output files created by Stan's cmdstan executable.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `m`                              : SampleModel
* `output_format`                      : Requested output format
```

### Optional arguments
```julia
* `include_internals`                  : Include internal parameters
* `chains=1:cpp_chains*julia_chains    : Which chains to include in output
* `start=1`                            : First sample to include in output
```
Not exported
"""
function read_csv_files(m::SampleModel, output_format::Symbol;
  include_internals=false,
  chains=1:m.num_chains,
  start=1,
  kwargs...)

  local a3d, monitors, index, idx, indvec, ftype, noofsamples

  # File path components of sample files (missing the "_$(i).csv" part)
  output_base = m.output_base
  name_base ="_chain"

  # How many samples?
  if m.save_warmup
    n_samples = floor(Int,
      (m.num_samples+m.num_warmups)/m.thin)
  else
    n_samples = floor(Int, m.num_samples/m.thin)
  end

  init_a3d = true
  current_chain = 0
  #println("Reading $(m.num_chains) chains.")

  # Read .csv files and return a3d[n_samples, parameters, n_chains]
  for i in 1:m.num_julia_chains   # Number of exec processes
    for k in 1:m.num_cpp_chains   # Number of cpp chains handled in cmdstan

      if (m.use_cpp_chains && m.check_num_chains) ||
        !m.use_cpp_chains || m.num_cpp_chains == 1
        csvfile = output_base*name_base*"_$(i).csv"
      else
        if i == 1
          csvfile = output_base*name_base*"_$(i)_$(k).csv"
        else
          csvfile = output_base*name_base*"_$(i)_$(k + i - 1).csv"
        end
      end
      #println("Reading "*csvfile)

      if isfile(csvfile)
        #println(csvfile*" found!")
        current_chain += 1
        instream = open(csvfile)

        # Skip initial set of commented lines, e.g. containing cmdstan version info, etc.
        skipchars(isspace, instream, linecomment='#')

        # First non-comment line contains names of variables
        line = Unicode.normalize(readline(instream), newline2lf=true)
        idx = split(strip(line), ",")
        index = [idx[k] for k in 1:length(idx)]
        indvec = 1:length(index)
        n_parameters = length(indvec)

        # Allocate a3d as we now know number of parameters
        if init_a3d
          init_a3d = false
          a3d = fill(0.0, n_samples, n_parameters, m.num_chains)
        end

        skipchars(isspace, instream, linecomment='#')
        for j in 1:n_samples
          skipchars(isspace, instream, linecomment='#')
          line = Unicode.normalize(readline(instream), newline2lf=true)
          if eof(instream) && length(line) < 2
            close(instream)
            break
          else
            flds = parse.(Float64, split(strip(line), ","))
            flds = reshape(flds[indvec], 1, length(indvec))
            a3d[j,:,current_chain] = flds
          end
        end   # read in samples
        #println("Filling $(current_chain) of $(size(a3d))")
      end   # read in next file if it exists
    end   # read in all cpp_chains
  end   # read in file all chains

  # Filtering of draws, parameters and chains before further processing

  cnames = convert.(String, idx[indvec])
  if include_internals
    snames = [Symbol(cnames[i]) for i in 1:length(cnames)]
    indices = 1:length(cnames)
  else
    pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
    snames = filter(p -> !(p in  pi), cnames)
    indices = Vector{Int}(indexin(snames, cnames))
  end

  #println(size(a3d))
  res = convert_a3d(a3d[start:end, indices, chains],
    snames, Val(output_format); kwargs...)

  (res, snames)

end   # end of read_samples

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
module LightNLP

include("intdict.jl")
include("depparser/DepParser.jl")
include("ner/NER.jl")
include("tokenization/Tokenization.jl")

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

mainres = ([4, 5, 3],
           [1, 5, 3])
bitres = ([true, true, false],
          [false, true, false])

tskprod(x) = @task for i in x; produce(i); end

for (dest, src, bigsrc, emptysrc, res) in [
    ([1, 2, 3], () -> [4, 5], () -> [1, 2, 3, 4, 5], () -> Int[], mainres),
    ([1, 2, 3], () -> 4:5, () -> 1:5, () -> 1:0, mainres),
    ([1, 2, 3], () -> tskprod(4:5), () -> tskprod(1:5), () -> tskprod(1:0), mainres),
    (falses(3), () -> trues(2), () -> trues(5), () -> trues(0), bitres)]

    @test copy!(copy(dest), src()) == res[1]
    @test copy!(copy(dest), 1, src()) == res[1]
    @test copy!(copy(dest), 2, src(), 2) == res[2]
    @test copy!(copy(dest), 2, src(), 2, 1) == res[2]

    @test copy!(copy(dest), 99, src(), 99, 0) == dest

    @test copy!(copy(dest), 1, emptysrc()) == dest
    x = emptysrc()
    exc = isa(x, AbstractArray) ? BoundsError : ArgumentError
    @test_throws exc copy!(dest, 1, emptysrc(), 1)

    for idx in (0, 4)
        @test_throws BoundsError copy!(dest, idx, src())
        @test_throws BoundsError copy!(dest, idx, src(), 1)
        @test_throws BoundsError copy!(dest, idx, src(), 1, 1)
        x = src()
        exc = isa(x, AbstractArray) ? BoundsError : ArgumentError
        @test_throws exc copy!(dest, 1, x, idx)
        x = src()
        exc = isa(x, AbstractArray) ? BoundsError : ArgumentError
        @test_throws exc copy!(dest, 1, x, idx, 1)
    end

    @test_throws ArgumentError copy!(dest, 1, src(), 1, -1)

    @test_throws BoundsError copy!(dest, bigsrc())

    @test_throws BoundsError copy!(dest, 3, src())
    @test_throws BoundsError copy!(dest, 3, src(), 1)
    @test_throws BoundsError copy!(dest, 3, src(), 1, 2)

    @test_throws BoundsError copy!(dest, 1, src(), 2, 2)
end

# test behavior of shallow and deep copying
let a = Any[[1]], q = QuoteNode([1])
    ca = copy(a); dca = @inferred(deepcopy(a))
    @test ca !== a
    @test ca[1] === a[1]
    @test dca !== a
    @test dca[1] !== a[1]
    @test deepcopy(q).value !== q.value
end

# issue #13124
let a = rand(3, 5)
    b = (a,a)
    c = deepcopy(b)
    @test c[1] === c[2]
end

# issue #14027
@test isnull(deepcopy(Nullable{Array}()))

# issue #15250
let a1 = Core.svec(1, 2, 3, []), a2 = Core.svec(1, 2, 3)
    a3 = Core.svec(a1,a1)
    b1 = deepcopy(a1)
    @test a1 == b1
    @test a1 !== b1
    @test a1[4] !== b1[4]
    b2 = deepcopy(a2)
    @test a2 === b2
    b3 = deepcopy(a3)
    @test a3 == b3
    @test a3 !== b3
    @test a3[1] === a3[2]
end

# issue #16667
let x = BigInt[1:1000;], y = deepcopy(x), v
    # Finalize the original values to make sure the deep copy is indeed
    # independent
    for v in x
        finalize(v)
    end
    # Allocate some memory to make it more likely to trigger an error
    # if `deepcopy` went wrong
    x = BigInt[1:1000;]
    @test y == x
end
let x = BigFloat[1:1000;], y, z, v
    y, z = setprecision(2) do
        deepcopy(x), BigFloat[1:1000;]
    end
    for v in x
        finalize(v)
    end
    x = BigFloat[1:1000;]
    # Make sure the difference in precision doesn't affect deep copy
    @test y == x
    # Check that the setprecision indeed does something
    @test z != x
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module SigmoidalProgramming

using JuMP

include("problems.jl")
include("sp.jl")
include("functions.jl")

end # module


text:
@testsection "Coordinate systems" begin
    s∂ = simplify ∘ ∂
    (x, y, z), (𝐞₁, 𝐞₂, 𝐞₃), ℬ = init_cartesian()
    (θ, ϕ, r), (𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ), ℬˢ = init_spherical()

    @testsection "Usual coordinate systems" begin
        @test components(𝐞ʳ ⊗ 𝐞ᵠ, ℬˢ) == components(𝐞₃ ⊗ 𝐞₂, ℬ) == components_canon(𝐞₃ ⊗ 𝐞₂)
    end

    @testsection "Partial derivatives" begin
        @test s∂(𝐞ʳ, θ) == 𝐞ᶿ
        @test s∂(𝐞ʳ, ϕ) == sin(θ) * 𝐞ᵠ
        @test s∂(𝐞ᵠ ⊗ 𝐞ᶿ, ϕ) == s∂(𝐞ᵠ, ϕ) ⊗ 𝐞ᶿ + 𝐞ᵠ ⊗ s∂(𝐞ᶿ, ϕ)
        @test s∂(𝐞ʳ ⊗ˢ 𝐞ᵠ, ϕ) == s∂(𝐞ʳ, ϕ) ⊗ˢ 𝐞ᵠ + 𝐞ʳ ⊗ˢ s∂(𝐞ᵠ, ϕ)
    end

    @testsection "Coordinate systems" begin
        # Cartesian
        Cartesian = coorsys_cartesian()
        𝐗 = getcoords(Cartesian)
        𝐄 = unitvec(Cartesian)
        ℬ = get_normalized_basis(Cartesian)
        𝛔 = Tens(SymmetricTensor{2,3}((i, j) -> SymFunction("σ$i$j", real = true)(𝐗...)))
        @test DIV(𝛔, Cartesian) ==
              sum([sum([∂(𝛔[i, j], 𝐗[j]) for j ∈ 1:3]) * 𝐄[i] for i ∈ 1:3])

        # Polar
        Polar = coorsys_polar()
        r, θ = getcoords(Polar)
        𝐞ʳ, 𝐞ᶿ = unitvec(Polar)
        ℬᵖ = get_normalized_basis(Polar)
        f = SymFunction("f", real = true)(r, θ)
        @test simplify(LAPLACE(f, Polar)) ==
              simplify(∂(r * ∂(f, r), r) / r + ∂(f, θ, θ) / r^2)

        # Cylindrical
        Cylindrical = coorsys_cylindrical()
        rθz = getcoords(Cylindrical)
        𝐞ʳ, 𝐞ᶿ, 𝐞ᶻ = unitvec(Cylindrical)
        ℬᶜ = get_normalized_basis(Cylindrical)
        r, θ, z = rθz
        𝐯 = Tens(Vec{3}(i -> SymFunction("v$(rθz[i])", real = true)(rθz...)), ℬᶜ)
        vʳ, vᶿ, vᶻ = getarray(𝐯)
        @test simplify(DIV(𝐯, Cylindrical)) ==
              simplify(∂(vʳ, r) + vʳ / r + ∂(vᶿ, θ) / r + ∂(vᶻ, z))

        # Spherical
        Spherical = coorsys_spherical()
        θ, ϕ, r = getcoords(Spherical)
        𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical)
        ℬˢ = get_normalized_basis(Spherical)
        for σⁱʲ ∈ ("σʳʳ", "σᶿᶿ", "σᵠᵠ")
            @eval $(Symbol(σⁱʲ)) = SymFunction($σⁱʲ, real = true)($r)
        end
        𝛔 = σʳʳ * 𝐞ʳ ⊗ 𝐞ʳ + σᶿᶿ * 𝐞ᶿ ⊗ 𝐞ᶿ + σᵠᵠ * 𝐞ᵠ ⊗ 𝐞ᵠ
        div𝛔 = simplify(DIV(𝛔, Spherical))
        @test simplify(div𝛔 ⋅ 𝐞ʳ) == simplify(∂(σʳʳ, r) + (2σʳʳ - σᶿᶿ - σᵠᵠ) / r)

        # Concentric sphere - hydrostatic part
        θ, ϕ, r = getcoords(Spherical)
        𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical)
        ℬˢ = get_normalized_basis(Spherical)
        𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Sym))
        𝟏 = tensId2(Val(3), Val(Sym))
        k, μ = symbols("k μ", positive = true)
        λ = k - 2μ / 3
        ℂ = 3k * 𝕁 + 2μ * 𝕂
        u = SymFunction("u", real = true)(r)
        𝐮 = u * 𝐞ʳ
        𝛆 = simplify(SYMGRAD(𝐮, Spherical))
        𝛔 = simplify(ℂ ⊡ 𝛆)
        # 𝛔 = simplify(λ * tr(𝛆) * 𝟏 + 2μ * 𝛆)
        @test dsolve(factor(simplify(DIV(𝛔, Spherical) ⋅ 𝐞ʳ)), u) ==
              Eq(u, symbols("C1") / r^2 + symbols("C2") * r)

        # Spheroidal
        Spheroidal = coorsys_spheroidal()
        OM = getOM(Spheroidal)
        @test simplify(LAPLACE(OM[1]^2, Spheroidal)) == 2

    end

end


text:
module MultiPoly

using DataStructures
using Combinatorics
using LinearAlgebra

export
    MPoly, terms, vars, nvars, generators, generator, monomials, deg,
    PolyUnion, newexps, oldexps,
    harmonicpolynomialdim, laplaceharmonicpol,
    evaluate, evaluate_basis, diff, integrate

import Base: zero, one,
    show, showcompact, print, length, endof, getindex, setindex!, copy, promote_rule, convert, eltype,
    *, /, //, -, +, ==, ^, divrem, conj, rem, real, imag, diff

include("mpoly.jl")
include("mpolyarithmetic.jl")
include("mpolyprinting.jl")
include("polyunion.jl")
include("mpolyevaluation.jl")
include("mpolycalculus.jl")

end # module


text:
push!(LOAD_PATH, "../src/")

using Documenter, Normalize

# Check docstring examples
DocMeta.setdocmeta!(Normalize, :DocTestSetup, :(using Normalize, DataFrames); recursive=true)

makedocs(modules=[Normalize], sitename="Normalize.jl")

deploydocs(repo="https://github.com/muhadamanhuri/Normalize.jl.git")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

####################
# LU Factorization #
####################
struct LU{T,S<:AbstractMatrix} <: Factorization{T}
    factors::S
    ipiv::Vector{BlasInt}
    info::BlasInt
    LU{T,S}(factors::AbstractMatrix{T}, ipiv::Vector{BlasInt}, info::BlasInt) where {T,S} = new(factors, ipiv, info)
end
LU(factors::AbstractMatrix{T}, ipiv::Vector{BlasInt}, info::BlasInt) where {T} = LU{T,typeof(factors)}(factors, ipiv, info)

# StridedMatrix
function lufact!(A::StridedMatrix{T}, pivot::Union{Val{false}, Val{true}} = Val(true)) where T<:BlasFloat
    if pivot === Val(false)
        return generic_lufact!(A, pivot)
    end
    lpt = LAPACK.getrf!(A)
    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end
function lufact!(A::HermOrSym, pivot::Union{Val{false}, Val{true}} = Val(true))
    copytri!(A.data, A.uplo, isa(A, Hermitian))
    lufact!(A.data, pivot)
end

"""
    lufact!(A, pivot=Val(true)) -> LU

`lufact!` is the same as [`lufact`](@ref), but saves space by overwriting the
input `A`, instead of creating a copy. An [`InexactError`](@ref)
exception is thrown if the factorization produces a number not representable by the
element type of `A`, e.g. for integer types.
"""
lufact!(A::StridedMatrix, pivot::Union{Val{false}, Val{true}} = Val(true)) = generic_lufact!(A, pivot)
function generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot} = Val(true)) where {T,Pivot}
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    ipiv = Vector{BlasInt}(minmn)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = abs(zero(T))
                for i = k:m
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
end

# floating point types doesn't have to be promoted for LU, but should default to pivoting
lufact(A::Union{AbstractMatrix{T}, AbstractMatrix{Complex{T}}},
    pivot::Union{Val{false}, Val{true}} = Val(true)) where {T<:AbstractFloat} =
        lufact!(copy(A), pivot)

# for all other types we must promote to a type which is stable under division
"""
    lufact(A [,pivot=Val(true)]) -> F::LU

Compute the LU factorization of `A`.

In most cases, if `A` is a subtype `S` of `AbstractMatrix{T}` with an element
type `T` supporting `+`, `-`, `*` and `/`, the return type is `LU{T,S{T}}`. If
pivoting is chosen (default) the element type should also support `abs` and
`<`.

The individual components of the factorization `F` can be accessed by indexing:

| Component | Description                         |
|:----------|:------------------------------------|
| `F[:L]`   | `L` (lower triangular) part of `LU` |
| `F[:U]`   | `U` (upper triangular) part of `LU` |
| `F[:p]`   | (right) permutation `Vector`        |
| `F[:P]`   | (right) permutation `Matrix`        |

The relationship between `F` and `A` is

`F[:L]*F[:U] == A[F[:p], :]`

`F` further supports the following functions:

| Supported function               | `LU` | `LU{T,Tridiagonal{T}}` |
|:---------------------------------|:-----|:-----------------------|
| [`/`](@ref)                      | ✓    |                        |
| [`\\`](@ref)                     | ✓    | ✓                      |
| [`inv`](@ref)                    | ✓    | ✓                      |
| [`det`](@ref)                    | ✓    | ✓                      |
| [`logdet`](@ref)                 | ✓    | ✓                      |
| [`logabsdet`](@ref)              | ✓    | ✓                      |
| [`size`](@ref)                   | ✓    | ✓                      |

# Examples
```jldoctest
julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> F = lufact(A)
Base.LinAlg.LU{Float64,Array{Float64,2}} with factors L and U:
[1.0 0.0; 1.5 1.0]
[4.0 3.0; 0.0 -1.5]
successful: true

julia> F[:L] * F[:U] == A[F[:p], :]
true
```
"""
function lufact(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}}) where T
    S = typeof(zero(T)/one(T))
    AA = similar(A, S, size(A))
    copy!(AA, A)
    lufact!(AA, pivot)
end
# We can't assume an ordered field so we first try without pivoting
function lufact(A::AbstractMatrix{T}) where T
    S = typeof(zero(T)/one(T))
    AA = similar(A, S, size(A))
    copy!(AA, A)
    F = lufact!(AA, Val(false))
    if issuccess(F)
        return F
    else
        AA = similar(A, S, size(A))
        copy!(AA, A)
        return lufact!(AA, Val(true))
    end
end

lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt))
lufact(F::LU) = F

lu(x::Number) = (one(x), x, 1)

"""
    lu(A, pivot=Val(true)) -> L, U, p

Compute the LU factorization of `A`, such that `A[p,:] = L*U`.
By default, pivoting is used. This can be overridden by passing
`Val(false)` for the second argument.

See also [`lufact`](@ref).

# Examples
```jldoctest
julia> A = [4. 3.; 6. 3.]
2×2 Array{Float64,2}:
 4.0  3.0
 6.0  3.0

julia> L, U, p = lu(A)
([1.0 0.0; 0.666667 1.0], [6.0 3.0; 0.0 1.0], [2, 1])

julia> A[p, :] == L * U
true
```
"""
function lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true))
    F = lufact(A, pivot)
    F[:L], F[:U], F[:p]
end

function convert(::Type{LU{T}}, F::LU) where T
    M = convert(AbstractMatrix{T}, F.factors)
    LU{T,typeof(M)}(M, F.ipiv, F.info)
end
convert(::Type{LU{T,S}}, F::LU) where {T,S} = LU{T,S}(convert(S, F.factors), F.ipiv, F.info)
convert(::Type{Factorization{T}}, F::LU{T}) where {T} = F
convert(::Type{Factorization{T}}, F::LU) where {T} = convert(LU{T}, F)

copy(A::LU{T,S}) where {T,S} = LU{T,S}(copy(A.factors), copy(A.ipiv), A.info)

size(A::LU) = size(A.factors)
size(A::LU,n) = size(A.factors,n)

function ipiv2perm(v::AbstractVector{T}, maxi::Integer) where T
    p = T[1:maxi;]
    @inbounds for i in 1:length(v)
        p[i], p[v[i]] = p[v[i]], p[i]
    end
    return p
end

function getindex(F::LU{T,<:StridedMatrix}, d::Symbol) where T
    m, n = size(F)
    if d == :L
        L = tril!(F.factors[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        return L
    elseif d == :U
        return triu!(F.factors[1:min(m,n), 1:n])
    elseif d == :p
        return ipiv2perm(F.ipiv, m)
    elseif d == :P
        return eye(T, m)[:,invperm(F[:p])]
    else
        throw(KeyError(d))
    end
end

issuccess(F::LU) = F.info == 0

function show(io::IO, F::LU)
    println(io, "$(typeof(F)) with factors L and U:")
    show(io, F[:L])
    println(io)
    show(io, F[:U])
    print(io, "\nsuccessful: $(issuccess(F))")
end

A_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
    @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info
A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) =
    A_ldiv_B!(UpperTriangular(A.factors),
    A_ldiv_B!(UnitLowerTriangular(A.factors), b[ipiv2perm(A.ipiv, length(b))]))
A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) =
    A_ldiv_B!(UpperTriangular(A.factors),
    A_ldiv_B!(UnitLowerTriangular(A.factors), B[ipiv2perm(A.ipiv, size(B, 1)),:]))

At_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
    @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, B) A.info
At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) =
    At_ldiv_B!(UnitLowerTriangular(A.factors),
    At_ldiv_B!(UpperTriangular(A.factors), b))[invperm(ipiv2perm(A.ipiv, length(b)))]
At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) =
    At_ldiv_B!(UnitLowerTriangular(A.factors),
    At_ldiv_B!(UpperTriangular(A.factors), B))[invperm(ipiv2perm(A.ipiv, size(B,1))),:]

Ac_ldiv_B!(F::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:Real} =
    At_ldiv_B!(F, B)
Ac_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
    @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, B) A.info
Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) =
    Ac_ldiv_B!(UnitLowerTriangular(A.factors),
    Ac_ldiv_B!(UpperTriangular(A.factors), b))[invperm(ipiv2perm(A.ipiv, length(b)))]
Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) =
    Ac_ldiv_B!(UnitLowerTriangular(A.factors),
    Ac_ldiv_B!(UpperTriangular(A.factors), B))[invperm(ipiv2perm(A.ipiv, size(B,1))),:]

At_ldiv_Bt(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
    @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info
At_ldiv_Bt(A::LU, B::StridedVecOrMat) = At_ldiv_B(A, transpose(B))

Ac_ldiv_Bc(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
    @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info
Ac_ldiv_Bc(A::LU, B::StridedVecOrMat) = Ac_ldiv_B(A, ctranspose(B))

function det(F::LU{T}) where T
    n = checksquare(F)
    issuccess(F) || return zero(T)
    P = one(T)
    c = 0
    @inbounds for i = 1:n
        P *= F.factors[i,i]
        if F.ipiv[i] != i
            c += 1
        end
    end
    s = (isodd(c) ? -one(T) : one(T))
    return P * s
end

function logabsdet(F::LU{T}) where T  # return log(abs(det)) and sign(det)
    n = checksquare(F)
    issuccess(F) || return log(zero(real(T))), log(one(T))
    c = 0
    P = one(T)
    abs_det = zero(real(T))
    @inbounds for i = 1:n
        dg_ii = F.factors[i,i]
        P *= sign(dg_ii)
        if F.ipiv[i] != i
            c += 1
        end
        abs_det += log(abs(dg_ii))
    end
    s = ifelse(isodd(c), -one(real(T)), one(real(T))) * P
    abs_det, s
end

inv!(A::LU{<:BlasFloat,<:StridedMatrix}) =
    @assertnonsingular LAPACK.getri!(A.factors, A.ipiv) A.info
inv!(A::LU{T,<:StridedMatrix}) where {T} =
    @assertnonsingular A_ldiv_B!(A.factors, copy(A), eye(T, size(A, 1))) A.info
inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A))

function _cond1Inf(A::LU{<:BlasFloat,<:StridedMatrix}, p::Number, normA::Real)
    if p != 1 && p != Inf
        throw(ArgumentError("p must be either 1 or Inf"))
    end
    return inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, normA))
end

# Tridiagonal

# See dgttrf.f
function lufact!(A::Tridiagonal{T}, pivot::Union{Val{false}, Val{true}} = Val(true)) where T
    n = size(A, 1)
    info = 0
    ipiv = Vector{BlasInt}(n)
    dl = A.dl
    d = A.d
    du = A.du
    du2 = A.du2

    @inbounds begin
        for i = 1:n
            ipiv[i] = i
        end
        for i = 1:n-2
            # pivot or not?
            if pivot === Val(false) || abs(d[i]) >= abs(dl[i])
                # No interchange
                if d[i] != 0
                    fact = dl[i]/d[i]
                    dl[i] = fact
                    d[i+1] -= fact*du[i]
                    du2[i] = 0
                end
            else
                # Interchange
                fact = d[i]/dl[i]
                d[i] = dl[i]
                dl[i] = fact
                tmp = du[i]
                du[i] = d[i+1]
                d[i+1] = tmp - fact*d[i+1]
                du2[i] = du[i+1]
                du[i+1] = -fact*du[i+1]
                ipiv[i] = i+1
            end
        end
        if n > 1
            i = n-1
            if pivot === Val(false) || abs(d[i]) >= abs(dl[i])
                if d[i] != 0
                    fact = dl[i]/d[i]
                    dl[i] = fact
                    d[i+1] -= fact*du[i]
                end
            else
                fact = d[i]/dl[i]
                d[i] = dl[i]
                dl[i] = fact
                tmp = du[i]
                du[i] = d[i+1]
                d[i+1] = tmp - fact*d[i+1]
                ipiv[i] = i+1
            end
        end
        # check for a zero on the diagonal of U
        for i = 1:n
            if d[i] == 0
                info = i
                break
            end
        end
    end
    LU{T,Tridiagonal{T}}(A, ipiv, convert(BlasInt, info))
end

factorize(A::Tridiagonal) = lufact(A)

function getindex(F::Base.LinAlg.LU{T,Tridiagonal{T}}, d::Symbol) where T
    m, n = size(F)
    if d == :L
        L = Array(Bidiagonal(ones(T, n), F.factors.dl, d))
        for i = 2:n
            tmp = L[F.ipiv[i], 1:i - 1]
            L[F.ipiv[i], 1:i - 1] = L[i, 1:i - 1]
            L[i, 1:i - 1] = tmp
        end
        return L
    elseif d == :U
        U = Array(Bidiagonal(F.factors.d, F.factors.du, d))
        for i = 1:n - 2
            U[i,i + 2] = F.factors.du2[i]
        end
        return U
    elseif d == :p
        return ipiv2perm(F.ipiv, m)
    elseif d == :P
        return eye(T, m)[:,invperm(F[:p])]
    end
    throw(KeyError(d))
end

# See dgtts2.f
function A_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where T
    n = size(A,1)
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            for i = 1:n-1
                ip = ipiv[i]
                tmp = B[i+1-ip+i,j] - dl[i]*B[ip,j]
                B[i,j] = B[ip,j]
                B[i+1,j] = tmp
            end
            B[n,j] /= d[n]
            if n > 1
                B[n-1,j] = (B[n-1,j] - du[n-1]*B[n,j])/d[n-1]
            end
            for i = n-2:-1:1
                B[i,j] = (B[i,j] - du[i]*B[i+1,j] - du2[i]*B[i+2,j])/d[i]
            end
        end
    end
    return B
end

function At_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where T
    n = size(A,1)
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            B[1,j] /= d[1]
            if n > 1
                B[2,j] = (B[2,j] - du[1]*B[1,j])/d[2]
            end
            for i = 3:n
                B[i,j] = (B[i,j] - du[i-1]*B[i-1,j] - du2[i-2]*B[i-2,j])/d[i]
            end
            for i = n-1:-1:1
                if ipiv[i] == i
                    B[i,j] = B[i,j] - dl[i]*B[i+1,j]
                else
                    tmp = B[i+1,j]
                    B[i+1,j] = B[i,j] - dl[i]*tmp
                    B[i,j] = tmp
                end
            end
        end
    end
    return B
end

# Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B)
function Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where T
    n = size(A,1)
    if n != size(B,1)
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    end
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            B[1,j] /= conj(d[1])
            if n > 1
                B[2,j] = (B[2,j] - conj(du[1])*B[1,j])/conj(d[2])
            end
            for i = 3:n
                B[i,j] = (B[i,j] - conj(du[i-1])*B[i-1,j] - conj(du2[i-2])*B[i-2,j])/conj(d[i])
            end
            for i = n-1:-1:1
                if ipiv[i] == i
                    B[i,j] = B[i,j] - conj(dl[i])*B[i+1,j]
                else
                    tmp = B[i+1,j]
                    B[i+1,j] = B[i,j] - conj(dl[i])*tmp
                    B[i,j] = tmp
                end
            end
        end
    end
    return B
end

/(B::AbstractMatrix,A::LU) = At_ldiv_Bt(A,B).'

# Conversions
convert(::Type{AbstractMatrix}, F::LU) = (F[:L] * F[:U])[invperm(F[:p]),:]
convert(::Type{AbstractArray}, F::LU) = convert(AbstractMatrix, F)
convert(::Type{Matrix}, F::LU) = convert(Array, convert(AbstractArray, F))
convert(::Type{Array}, F::LU) = convert(Matrix, F)
full(F::LU) = convert(AbstractArray, F)

function convert(::Type{Tridiagonal}, F::Base.LinAlg.LU{T,Tridiagonal{T}}) where T
    n = size(F, 1)

    dl  = copy(F.factors.dl)
    d   = copy(F.factors.d)
    du  = copy(F.factors.du)
    du2 = copy(F.factors.du2)

    for i = n - 1:-1:1
        li         = dl[i]
        dl[i]      = li*d[i]
        d[i + 1]  += li*du[i]
        if i < n - 1
            du[i + 1] += li*du2[i]
        end

        if F.ipiv[i] != i
            tmp   = dl[i]
            dl[i] = d[i]
            d[i]  = tmp

            tmp      = d[i + 1]
            d[i + 1] = du[i]
            du[i]    = tmp

            if i < n - 1
                tmp       = du[i + 1]
                du[i + 1] = du2[i]
                du2[i]    = tmp
            end
        end
    end
    return Tridiagonal(dl, d, du)
end
convert(::Type{AbstractMatrix}, F::Base.LinAlg.LU{T,Tridiagonal{T}}) where {T} =
    convert(Tridiagonal, F)
convert(::Type{AbstractArray}, F::Base.LinAlg.LU{T,Tridiagonal{T}}) where {T} =
    convert(AbstractMatrix, F)
convert(::Type{Matrix}, F::Base.LinAlg.LU{T,Tridiagonal{T}}) where {T} =
    convert(Array, convert(AbstractArray, F))
convert(::Type{Array}, F::Base.LinAlg.LU{T,Tridiagonal{T}}) where {T} =
    convert(Matrix, F)
full(F::Base.LinAlg.LU{T,Tridiagonal{T}}) where {T} = convert(AbstractArray, F)

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
# ---
# title: 1349. Maximum Students Taking Exam
# id: problem1349
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/maximum-students-taking-exam/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix `seats`  that represent seats distributions in a
# classroom. If a seat is broken, it is denoted by `'#'` character otherwise it
# is denoted by a `'.'` character.
#
# Students can see the answers of those sitting next to the left, right, upper
# left and upper right, but he cannot see the answers of the student sitting
# directly in front or behind him. Return the **maximum** number of students
# that can take the exam together without any cheating being possible..
#
# Students must be placed in seats in good condition.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/29/image.png)
#
#
#
#     Input: seats = [["#",".","#","#",".","#"],
#                     [".","#","#","#","#","."],
#                     ["#",".","#","#",".","#"]]
#     Output: 4
#     Explanation: Teacher can place 4 students in available seats so they don't cheat on the exam.
#
#
# **Example 2:**
#
#
#
#     Input: seats = [[".","#"],
#                     ["#","#"],
#                     ["#","."],
#                     ["#","#"],
#                     [".","#"]]
#     Output: 3
#     Explanation: Place all students in available seats.
#
#
#
# **Example 3:**
#
#
#
#     Input: seats = [["#","."," **.** ",".","#"],
#                     [" **.** ","#"," **.** ","#"," **.** "],
#                     [" **.** ",".","#","."," **.** "],
#                     [" **.** ","#"," **.** ","#"," **.** "],
#                     ["#","."," **.** ",".","#"]]
#     Output: 10
#     Explanation: Place students in available seats in column 1, 3 and 5.
#
#
#
#
# **Constraints:**
#
#   * `seats` contains only characters `'.' and``'#'.`
#   * `m == seats.length`
#   * `n == seats[i].length`
#   * `1 <= m <= 8`
#   * `1 <= n <= 8`
#
#
## @lc code=start
using LeetCode

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


text:
md"""
# Cutting
"""

"""
    cut(panel, axis, at; kerf, cost)::(panel1, panel2)

Cut panel at the specified distance along axis.
The first returned panel is the piece cut to that distance.
The second is what remains (accounting for kerf).

An empty Tuple is returned if the cut can't be made.
"""
function cut(panel::CuttablePanel,
             axis::Axis,
             at::LengthType;
             kerf = (1/8)u"inch",
             cost = money(0.0))
    if distance(panel, axis) < at
        return (())
    end
    # We might short circuit the test of the function if
    #   distance(panel, axis) == at
    # but then we might return a BoughtPanel rather than a Panel.
    # It does affect the cost though:
    cost = panel.cost + (distance(panel, axis) > at ? cost : 0)
    mult(v1, v2) = map(*, v1, v2)
    at_vector = at * unit_vector(axis)
    kerf_vector = kerf * unit_vector(axis)
    cross_vector = mult(panel_lw(panel), unit_vector(other(axis)))
    lw1 = at_vector + cross_vector
    prev = at_vector + kerf_vector
    xy2 = panel_xy(panel) + prev
    lw2 = panel_lw(panel) - prev
    if lw2[1] < 0u"inch" || lw2[2]< 0u"inch"
        lw2 = zerovector
    end
    area1 = area(lw1)
    area2 = area(lw2)
    total_area = area1 + area2
    panel1 = Panel(; length=lw1[1], width=lw1[2],
                   cut_from=panel, cut_at=at, cut_axis=axis,
                   x=panel.x, y=panel.y,
	           cost=cost * area1 / total_area)
    # There might be no off-cut if panel was bigger than at by kerf,
    # for example:
    if area2 > zeroarea
        # Maybe cost should be apportioned by area rather than by
        # count?  That would slightly favor shorter cuts.
        (panel1,
         Panel(; length=lw2[1], width=lw2[2],
               cut_from=panel, cut_at=at, cut_axis=axis,
               x=xy2[1], y=xy2[2],
	       cost=cost * area2 / total_area))
    else
        (panel1,)
    end
end

export cut

let
    KERF = (1/8)u"inch"
    panel1 = BoughtPanel(AvailablePanel("30 by 60", 30u"inch", 60u"inch", money(20.00)))

    # Cut panel1 at 25" down LengthAxis:
    panel2, panel3 = cut(panel1, LengthAxis(), 25u"inch"; kerf=KERF)
    @assert panel2.length == 25u"inch"  """got $(panel2.length), expected $(25u"inch")"""
    @assert panel2.width == 30u"inch"
    @assert panel3.length == panel1.length - panel2.length - KERF
    @assert panel3.width == 30u"inch"   """got $(panel3.width), expected $(30u"inch")"""
    @assert panel2.x == panel1.x   "panel2.x: got $(panel2.x), expected $(panel1.x)"
    @assert panel2.y == panel1.y
    @assert panel3.x == panel1.x + panel2.length + KERF
    @assert panel3.y == panel1.y

    # Cut panel2 at 10" down WidthAxis:
    panel4, panel5 = cut(panel2, WidthAxis(), 10u"inch", kerf=KERF)
    @assert panel4.width == 10u"inch"
    @assert panel4.length == panel2.length
    @assert panel4.x == panel2.x
    @assert panel4.y == panel2.y
    @assert panel5.x == panel2.x
    @assert panel5.y == panel4.y + panel4.width + KERF
    @assert panel5.length == panel2.length
    @assert panel5.width == panel2.width - panel4.width - KERF
end

let
    KERF = (1/8)u"inch"
    panel1 = BoughtPanel(AvailablePanel("4 x 8 x 1/2", 4u"ft", 8u"ft", money(95.00)))
    at = 22u"inch"
    cut1, cut2 = cut(panel1, LengthAxis(), at, kerf=KERF)
    @assert cut1.length == at  "got $(cut1.length), expected $(at)"
    @assert cut1.width == panel1.width  "got $(cut1.width), expected $(panel1.width)"
    @assert cut2.width == panel1.width  "got $(cut2.width), expected $(panel1.width)"
    @assert cut2.length ≈ panel1.length - at - KERF (
        "got $(cut2.length), expected $(panel1.length - at - KERF)")
end


text:
using GalerkinSparseGrids
using Test
using HCubature
using ODE
using SparseArrays
using LinearAlgebra

tests = [
    "elementary.jl",
    "basis_constructor.jl",
    "hier_DG.jl",
    "vhier_DG.jl",
    "transformations.jl",
    "differentiation.jl",
    "solvers.jl"
]

@testset "GalerkinSparseGrids" begin
    for filename in tests
        name = first(splitext(filename))
        include(filename)
    end
end

# Later, we will also test the examples:

# include(joinpath(dirname(@__FILE__), "..", "examples", "interpolation.jl"))
# include(joinpath(dirname(@__FILE__), "..", "examples", "differentiation.jl"))
# include(joinpath(dirname(@__FILE__), "..", "examples", "traveling_wave.jl"))
# include(joinpath(dirname(@__FILE__), "..", "examples", "vlasov_evolve.jl"))


text:
"""
    Problem(network, input, output)

Problem definition for neural verification.
- `network` is of type `Network`
- `input` belongs to `AbstractPolytope` in `LazySets.jl`
- `output` belongs to `AbstractPolytope` in `LazySets.jl`

The verification problem consists of: for all  points in the input set,
the corresponding output of the network must belong to the output set.
"""
struct Problem{P<:AbstractPolytope, Q<:AbstractPolytope}
    network::Network
    input::P
    output::Q
end

"""
    Result
Supertype of all result types.

See also: [`BasicResult`](@ref), [`CounterExampleResult`](@ref), [`AdversarialResult`](@ref), [`ReachabilityResult`](@ref)
"""
abstract type Result end

status(result::Result) = result.status

function validate_status(st::Symbol)
    @assert st ∈ (:SAT, :UNSAT, :Unknown) "unexpected status code: `:$st`.\nOnly (:SAT, :UNSAT, :Unknown) are accepted"
    return st
end

"""
    BasicResult(status::Symbol)

Result type that captures whether the input-output constraint is satisfied.
Possible status values:\n
    :SAT (io constraint is satisfied always)\n
    :UNSAT (io constraint is violated)\n
    :Unknown (could not be determined)
"""
struct BasicResult <: Result
    status::Symbol
end

"""
    CounterExampleResult(status, counter_example)

Like `BasicResult`, but also returns a `counter_example` if one is found (if status = :UNSAT).
The `counter_example` is a point in the input set that, after the NN, lies outside the output set.
"""
struct CounterExampleResult <: Result
    status::Symbol
    counter_example::Vector{Float64}
    CounterExampleResult(s, ce) = new(validate_status(s), ce)
end

"""
    AdversarialResult(status, max_disturbance)

Like `BasicResult`, but also returns the maximum allowable disturbance in the input (if status = :UNSAT).
"""
struct AdversarialResult <: Result
	status::Symbol
	max_disturbance::Float64
    AdversarialResult(s, md) = new(validate_status(s), md)
end

"""
    ReachabilityResult(status, reachable)

Like `BasicResult`, but also returns the output reachable set given the input constraint (if status = :UNSAT).
"""
struct ReachabilityResult <: Result
	status::Symbol
	reachable::Vector{<:AbstractPolytope}
    ReachabilityResult(s, r) = new(validate_status(s), r)
end

# Additional constructors:
CounterExampleResult(s) = CounterExampleResult(s, Float64[])
AdversarialResult(s)    = AdversarialResult(s, -1.0)
ReachabilityResult(s)   = ReachabilityResult(s, AbstractPolytope[])


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Represents downward API info for projecting into a projected volume. Note that this is identical to a downwardAPI volume source without the default mode.

    IoK8sApiCoreV1DownwardAPIProjection(;
        items=nothing,
    )

    - items::Vector{IoK8sApiCoreV1DownwardAPIVolumeFile} : Items is a list of DownwardAPIVolume file
"""
mutable struct IoK8sApiCoreV1DownwardAPIProjection <: SwaggerModel
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1DownwardAPIVolumeFile} } # spec name: items

    function IoK8sApiCoreV1DownwardAPIProjection(;items=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1DownwardAPIProjection, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        o
    end
end # type IoK8sApiCoreV1DownwardAPIProjection

const _property_map_IoK8sApiCoreV1DownwardAPIProjection = Dict{Symbol,Symbol}(Symbol("items")=>Symbol("items"))
const _property_types_IoK8sApiCoreV1DownwardAPIProjection = Dict{Symbol,String}(Symbol("items")=>"Vector{IoK8sApiCoreV1DownwardAPIVolumeFile}")
Base.propertynames(::Type{ IoK8sApiCoreV1DownwardAPIProjection }) = collect(keys(_property_map_IoK8sApiCoreV1DownwardAPIProjection))
Swagger.property_type(::Type{ IoK8sApiCoreV1DownwardAPIProjection }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1DownwardAPIProjection[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1DownwardAPIProjection }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1DownwardAPIProjection[property_name]

function check_required(o::IoK8sApiCoreV1DownwardAPIProjection)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1DownwardAPIProjection }, name::Symbol, val)
end


text:
ENV["GKSwstype"] = "100"  # set 'GR environment' to 'no output' (for Travis CI)
using Documenter, LazySets
import Polyhedra, Optim, Expokit, TaylorModels, Distributions, MiniQhull, Symbolics

include("init.jl")

makedocs(
    sitename = "LazySets.jl",
    modules = [LazySets, Approximations, Arrays, LazySets.Parallel],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/juliareach.css"]),
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Getting Started" => "man/getting_started.md",
            "A Tour of LazySets" => "man/tour.md",
            "Introduction to Convex Sets" => "man/convex_sets.md",
            "Polyhedral Approximations" => "man/polyhedral_approximations.md",
            "Decomposing an Affine Map" => "man/decompose_example.md",
            "Fast 2D LPs" => "man/fast_2d_LPs.md",
            "Iterative Refinement" => "man/iterative_refinement.md",
            "Interval Hulls" => "man/interval_hulls.md",
            "Convex Hulls" => "man/convex_hulls.md",
            "Operations on Sets" => "man/set_operations.md",
            "A Reachability Algorithm" => "man/reach_zonotopes.md",
            "A Hybrid Reachability Algorithm" => "man/reach_zonotopes_hybrid.md",
            "Concrete Polyhedra" => "man/concrete_polyhedra.md",
            "Parallel Approximations" => "man/parallel_approximations.md",
            "Lazy Intersections" => "man/lazy_intersections.md"
        ],
        "Library" => Any[
            "Interfaces" => "lib/interfaces.md",
            "Sets" => [
                "Ball1" => "lib/sets/Ball1.md",
                "Ball2" => "lib/sets/Ball2.md",
                "BallInf" => "lib/sets/BallInf.md",
                "Ballp" => "lib/sets/Ballp.md",
                "Ellipsoid" => "lib/sets/Ellipsoid.md",
                "EmptySet" => "lib/sets/EmptySet.md",
                "HalfSpace" => "lib/sets/HalfSpace.md",
                "HParallelotope" => "lib/sets/HParallelotope.md",
                "HPolygon" => "lib/sets/HPolygon.md",
                "HPolygonOpt" => "lib/sets/HPolygonOpt.md",
                "HPolyhedron" => "lib/sets/HPolyhedron.md",
                "HPolytope" => "lib/sets/HPolytope.md",
                "Hyperplane" => "lib/sets/Hyperplane.md",
                "Hyperrectangle" => "lib/sets/Hyperrectangle.md",
                "Interval" => "lib/sets/Interval.md",
                "Line2D" => "lib/sets/Line2D.md",
                "Line" => "lib/sets/Line.md",
                "LineSegment" => "lib/sets/LineSegment.md",
                "PolynomialZonotope" => "lib/sets/PolynomialZonotope.md",
                "RotatedHyperrectangle" => "lib/sets/RotatedHyperrectangle.md",
                "Singleton" => "lib/sets/Singleton.md",
                "Star" => "lib/sets/Star.md",
                "Universe" => "lib/sets/Universe.md",
                "VPolygon" => "lib/sets/VPolygon.md",
                "VPolytope" => "lib/sets/VPolytope.md",
                "ZeroSet" => "lib/sets/ZeroSet.md",
                "Zonotope" => "lib/sets/Zonotope.md",
            ],
            "Lazy Operations" => [
                "AffineMap" => "lib/lazy_operations/AffineMap.md",
                "Bloating" => "lib/lazy_operations/Bloating.md",
                "CartesianProduct" => "lib/lazy_operations/CartesianProduct.md",
                "Complement" => "lib/lazy_operations/Complement.md",
                "ConvexHull" => "lib/lazy_operations/ConvexHull.md",
                "ExponentialMap" => "lib/lazy_operations/ExponentialMap.md",
                "Intersection" => "lib/lazy_operations/Intersection.md",
                "LinearMap" => "lib/lazy_operations/LinearMap.md",
                "InverseLinearMap" => "lib/lazy_operations/InverseLinearMap.md",
                "MinkowskiSum" => "lib/lazy_operations/MinkowskiSum.md",
                "Rectification" => "lib/lazy_operations/Rectification.md",
                "ResetMap" => "lib/lazy_operations/ResetMap.md",
                "SymmetricIntervalHull" => "lib/lazy_operations/SymmetricIntervalHull.md",
                "Translation" => "lib/lazy_operations/Translation.md",
                "UnionSet" => "lib/lazy_operations/UnionSet.md",
            ],
            "Concrete Operations" => "lib/binary_functions.md",
            "Conversions between set representations" => "lib/conversion.md",
            "Comparisons" => "lib/comparisons.md",
            "Approximations" => "lib/approximations.md",
            "Utility Functions" => "lib/utils.md",
            "Parallel" => "lib/parallel.md",
        ],
        "About" => "about.md"
    ],
    doctest = false,
    strict = true
)

deploydocs(
    repo = "github.com/JuliaReach/LazySets.jl.git"
)


text:
module SoftConfidenceWeighted

export classify, update, SCWParameter

import DataStructures.DefaultDict

const MIN_CONFIDENCE = 0.0
const MAX_CONFIDENCE = 1.0
const DEFAULT_CONFIDENCE = 0.7
const MIN_AGGRESSIVENESS = 0.0
const DEFAULT_AGGRESSIVENESS = 1.0
const VALID_LABEL = [1, -1]
const ERF_ORDER = 30

type SCWParameter
  aggressiveness::Float64
  phi::Float64
  psi::Float64
  zeta::Float64
  mu::DefaultDict
  sigma::DefaultDict

  function SCWParameter(confidence::Float64 = 0.7, aggressiveness::Float64 = 1.0)
    if confidence < MIN_CONFIDENCE
      confidence = MIN_CONFIDENCE
    elseif confidence > MAX_CONFIDENCE
      confidence = MAX_CONFIDENCE
    end

    if aggressiveness < MIN_AGGRESSIVENESS
      aggressiveness = MIN_AGGRESSIVENESS
    end
    phi = probit(confidence)
    psi = 1.0 + phi * phi / 2.0
    zeta = 1.0 + phi * phi

    mu = DefaultDict(0.0)
    sigma = DefaultDict(1.0)

    new(aggressiveness, phi, psi, zeta, mu, sigma)
  end
end

function classify(params::SCWParameter, data::Dict)
  margin = 0.0
  for (feature, weight) = data
    if haskey(params.mu, feature)
      margin += params.mu[feature] * weight
    end
  end
  if margin > 0.0
    return 1
  end
  -1
end

function update(params::SCWParameter, data::Dict, label::Int64)
  if !in(label, VALID_LABEL)
    error("invalid label: $label")
  end

  sigma_x = get_sigma_x(params, data)
  margin_mean, variance = get_margin_mean_and_variance(params, label, data, sigma_x)

  if params.phi * sqrt(variance) <= margin_mean
    return params
  end

  alpha, beta = get_alpha_and_beta(params, margin_mean, variance)
  if alpha == 0.0 || beta == 0.0
    return params
  end

  for (feature, weight) = sigma_x
    params.mu[feature] += alpha * label * weight
    params.sigma[feature] -= beta * weight * weight
  end

  params
end

function get_sigma_x(params::SCWParameter, data::Dict)
  sigma_x = DefaultDict(1.0)
  for (feature, weight) = data
    sigma_x[feature] *= params.sigma[feature] * weight
  end
  sigma_x
end

function get_margin_mean_and_variance(params::SCWParameter, label::Int64, data::Dict, sigma_x::DefaultDict)
  margin_mean = 0.0
  variance    = 0.0
  for (feature, weight) = data
    margin_mean += params.mu[feature] * weight
    variance    += sigma_x[feature] * weight
  end
  (margin_mean * label, variance)
end

function get_alpha_and_beta(params::SCWParameter, margin_mean::Float64, variance::Float64)
  alpha_den = variance * params.zeta
  if alpha_den == 0.0
    return (0.0, 0.0)
  end

  term1 = margin_mean * params.phi / 2.0
  alpha = (-1.0 * margin_mean * params.psi + params.phi * sqrt(term1^2 + alpha_den)) / alpha_den
  if alpha <= 0.0
    return (0.0, 0.0)
  end

  if alpha >= params.aggressiveness
    alpha = params.aggressiveness
  end

  beta_num = alpha * params.phi
  term2    = variance * beta_num
  beta_den = term2 + (-1.0 * term2 + sqrt(term2^2 + 4.0 * variance)) / 2.0
  if beta_den == 0.0
    return (0.0, 0.0)
  end

  (alpha, beta_num / beta_den)
end

function probit(p)
  sqrt(2.0) * erf_inv(2.0 * p - 1.0)
end

function erf_inv(z)
  value = 1.0
  term  = 1.0
  c_memo = [1.0]

  for n in 1:ERF_ORDER
    term *= pi * z^2 / 4.0
    c = 0.0
    for m in 0:n-1
      c += c_memo[m + 1] * c_memo[n - m] / (m + 1.0) / (2.0 * m + 1.0)
    end
    push!(c_memo, c)
    value += c * term / (2.0 * n + 1.0)
  end
  sqrt(pi) * z * value / 2.0
end
end