text:
using LatexLH
using Test

testDir = joinpath(@__DIR__, "test_files");

function beamer_test()
	@testset "Beamer" begin
		figPath = joinpath(testDir, "afqtMeanByQual.pdf");
		lineV = figure_slide("Title", figPath);
		@test isa(lineV, Vector{String})

		fPath = joinpath(testDir, "beamer_test.tex");
		isfile(fPath)  &&  rm(fPath);
		open(fPath, "w") do io
			write_beamer_header(io);
			for figName ∈ ("afqtMeanByQual.pdf", "fracGradByQual.pdf",
				"parameter_table_test.tex")

				figPath = joinpath(testDir, figName);
				write_figure_slide(io, figName, figPath);
			end
			write_beamer_footer(io);
		end
		@test isfile(fPath)

		success = typeset_file(fPath);
		@test success
	end
end

function param_tb_test()
	@testset "ParameterTable" begin
		pt = ParameterTable();
		add_row!(pt, "p1", "param 1", "1.23");
		# Not robust. Skipping 2 header rows
		@test nrows(pt.tb) == 3
		@test pt.tb.bodyV[3] == "p1 & param 1 & 1.23"
		add_row!(pt, "p2", "param 2", "2.34");
		write_table(pt, joinpath(testDir, "parameter_table_test.tex"))
	end
end

@testset "LatexLH" begin
	@testset "CellColor" begin
		x = LatexLH.CellColor("blue", 40);
		s = color_string(x);
		# println(s)
		@test s == "\\cellcolor{blue!40}"

		x = LatexLH.CellColor("green", 0)
		s = color_string(x);
		# println(s)
		@test s == ""
	end

	@testset "Cell" begin
		c = LatexLH.Cell("abc", 1, 'c', LatexLH.CellColor("blue", 40))
		@test c.text == "abc"
		@test cell_string(c) == "\\cellcolor{blue!40}abc"

		c = LatexLH.Cell("abc");
		@test c.text == "abc"

		c = LatexLH.Cell("abc", 'c');
		@test c.text == "abc"
		@test c.align == 'c'
		@test cell_string(c) == "abc"

		c = LatexLH.Cell("abc", 2, 'c', LatexLH.CellColor("blue", 40))
		@test cell_string(c) == "\\multicolumn{2}{c}{\\cellcolor{blue!40}abc}"
	end

	@testset "Table" begin
		nCols = 5;
		tb = LatexLH.Table(nCols, "lSSSS");
		LatexLH.add_row!(tb, " & \\multicolumn{3}{c}{Heading 1} & ");
		@test length(tb.bodyV) == 1

		cellV = [LatexLH.Cell("cell1"),  LatexLH.Cell("cell2", 4, 'c')]
		rowStr = LatexLH.make_row(tb, cellV);
		@test rowStr == "cell1 & \\multicolumn{4}{c}{cell2}"
		LatexLH.add_row!(tb, rowStr);
		@test length(tb.bodyV) == 2

		LatexLH.add_row!(tb, "\\cline{2-5}")
		LatexLH.add_row!(tb, "x1 & 1.2 & 2.3 & 3.4 & 4.5")

		filePath = joinpath(testDir, "TableTest.tex")
		if isfile(filePath)
			rm(filePath)
		end
		write_table(tb, filePath)
		@test isfile(filePath)
	end

	param_tb_test();
    beamer_test();
    include("symbol_table_test.jl")
end

# -----------


text:
using FinEtools
using add2FinEtools

radius, nperradius = 5.0, 8
fens, fes = Q4spheren(radius, nperradius)
File =  "Sphere_mesh.vtk"
vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4);
@async run(`"paraview.exe" $File`)

ext: jl

max_stars_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_issues_repo_name: PetrKryslUCSD/FinEtoolsUseCase

max_forks_repo_name: PetrKryslUCSD/FinEtoolsUseCase


text:
srand(1)

let
    @testset "pinvfactors $((m, n, k))" for (m, n, k) in ((100, 50, 40), (50, 100, 40))
        X = randsparse(m, n, k, 1000/m/n)
        XM = Matrix(X)
        XMI = pinv(XM)
        pin = pinvfact(X)
        @test rank(pin) == k
        @test pin \ Matrix(1.0I, m, m) ≈ XMI
        @test pin' \ Matrix(1.0I, n, n) ≈ XMI'
        pin = pinvfact(XM)
        @test rank(pin) == k
        @test pin \ Matrix(1.0I, m, m) ≈ XMI
        @test pin' \ Matrix(1.0I, n, n) ≈ XMI'
    end
end


text:
import SparseArrays
import NonlinearEquations
@NonlinearEquations.equations function gw1d(h, Q)
	NonlinearEquations.setnumequations(3)
	NonlinearEquations.addterm(1, -2 * h[1] + h[2] + Q[1])
	NonlinearEquations.addterm(2, h[1] - 2 * h[2] + h[3] + Q[1])
	NonlinearEquations.addterm(3, -2 * h[3] + h[2] + Q[3])
end
function solveforh(Q)
	res(h) = gw1d_residuals(h, Q)
	jac(h) = gw1d_h(h, Q)
	h0 = zeros(N)
	h = jac(h0) \ -res(h0)
	return h
end

N = 3
Q = 8 * ones(N) / N ^ 2
h = solveforh(Q)


text:
abstract AbstractFeedForwardUnit

type RectifiedLinearUnit <: AbstractFeedForwardUnit
    W:: Vector{Float64}
    bias:: Float64
    activationfunction:: Function

    function RectifiedLinearUnit(W:: Vector{Float64}, bias:: Float64)
        function activationfunction(x)
            max(0,x)
        end
        new(W, bias, activationfunction)
    end
end

type LogisticUnit <: AbstractFeedForwardUnit
    W:: Vector{Float64}
    bias:: Float64
    activationfunction:: Function

    function LogisticUnit(W:: Vector{Float64}, bias:: Float64)
        new(W, bias, σ)
    end
end

function forward(u:: AbstractFeedForwardUnit, X::Array{Float64,2})
    u.activationfunction(X*u.W + u.bias)
end

function forward(u:: AbstractFeedForwardUnit, x::Vector{Float64})
    u.activationfunction(u.W'*x + u.bias)
end

max_stars_repo_name: helgith/ExampleNeuralNetworks

max_issues_repo_name: helgith/ExampleNeuralNetworks

max_forks_repo_name: helgith/ExampleNeuralNetworks


text:
function rotate(n::Int, target::Char)
    ascii = Int(target)

    if 97 <= ascii <= 122
        # Lower case letters
        code = ascii - 97
        code_1 = (code + n) % 26

        return Char(code_1 + 97)

    elseif 65 <= ascii <= 90
        # Upper case letters
        code = ascii - 65
        code_1 = (code + n) % 26

        return Char(code_1 + 65)

    else
        return target
    end

end

function rotate(n::Int, target::String)
    return String([rotate(n, c) for c in target])
end

for i in 0:26
    macro_name = Symbol("R$(i)_str")
    @eval begin
        macro $macro_name(p)
            rotate($i, p)
        end
    end
end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
module hmm

type HMM{T}
   n_hid
   n_obs
   a :: Array{T,2}
   b :: Array{T,2}
   pi :: Array{T,1}
end

function HMM(n_hid, n_obs)
   HMM(n_hid, n_obs, normalize_stoch_map(rand((n_hid, n_hid))),
       normalize_stoch_map(rand((n_obs, n_hid))),
       normalize_stoch_map(rand(n_hid)))
end

function get_forward(m, obs)

   function rescale!(x)
      fac = 1. / sum(x)
      x[:] *= fac
      fac
   end

   steps = size(obs)
   scale = Array(eltype(m.a), steps)
   alpha = Array(eltype(m.a), (steps, m.n_hid))

   alpha[1, :] = m.pi .* m.b[obs[1], :]
   scale[1] = rescale!(alpha[1, :])

   for t in 2:steps
      alpha[t, :] = (m.a * alpha[t-1, :]) .* m.b[obs[t], :]
      scale[t] = rescale!(alpha[t, :])
   end

   alpha, scale
end

function get_backward(m, obs, scale)
   steps = size(obs)
   beta = Array(eltype(m.a), (steps, m.n_hid))
   beta[steps, :] = scale[steps]

   for t in steps-1:-1:1
      beta[t, :] = scale[t] * (m.b[obs[t+1], :] .* beta[t+1, :]) * m.a
   end
   beta
end

function iterate(m, obs)
   alpha, scale = get_forward(m, obs)
   beta = get_backward(m, obs, scale)

   gamma = alpha .* beta
   gamma /= broadcast(/, gamma, sum(gamma, 2))

   mn =
end

function normalize_stoch_map(a)
   normalization = sum(a, 1)
   broadcast(/, a, normalization)
end

end

max_stars_repo_name: dseuss/stuff

max_issues_repo_name: dseuss/stuff

max_forks_repo_name: dseuss/stuff


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Mesa"
version = v"20.1.5"

# Collection of sources required to complete build
sources = [
    ArchiveSource("ftp://ftp.freedesktop.org/pub/mesa/mesa-$version.tar.xz", "fac1861e6e0bf1aec893f8d86dbfb9d8a0f426ff06b05256df10e3ad7e02c69b"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir build
cd build
apk add py3-mako
meson -D b_ndebug=true -D buildtype=release -D strip=true -D llvm=false ../mesa* --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
mv $prefix/bin/opengl32.dll $prefix/bin/opengl32sw.dll
install_license ../mesa*/docs/license.html
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:x86_64),
    Windows(:i686)
]

# The products that we will ensure are always built
products = [
    LibraryProduct("opengl32sw", :opengl32sw; dont_dlopen=true)
]

# Dependencies that must be installed before this package can be built
dependencies = [
  Dependency("Zlib_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
struct SemiconductorInterval{T}
    semiconductor::SemiconductorData
    S::T
end

function belongs(x,S::SemiconductorInterval )
    DomainSets.SVector(x...) in S.S
end

#∈(x,S::SemiconductorInterval) = belongs(x,S)

const SemiconductorDataInput = AbstractVector{SemiconductorInterval}

struct MOSFETInputDeck{T}
    Source::SemiconductorInterval
    Gate::MOSStructure
    Drain::SemiconductorInterval
    region
    DrainContact::T
    GateContact::T
    SourceContact::T
    Lch::Number
end
getSemi(x,S) = S[belongs.(Ref(x),S)][1].semiconductor

function MOSFETInputDeck(N_d,Nb,T,l,h,xd,ld,tox,type;lscal=1.0,ϕₘ=4.15,gate_ov=nothing)
    if type == :PMOS
        SD = SemiconductorData(T,BoltzmanDist(),PSilicon(N_d))
        SG = SemiconductorData(T,BoltzmanDist(),NSilicon(Nb))
    else
        SD = SemiconductorData(T,BoltzmanDist(),NSilicon(N_d))
        SG = SemiconductorData(T,BoltzmanDist(),PSilicon(Nb))
    end
    L = 2*ld+l
    if isnothing(lscal)
        lscal = L
    end
    if isnothing(gate_ov)
        gate_ov = 0.1*ld/lscal
    else
        gate_ov = gate_ov/lscal
    end
    SourceInt = (0..ld/lscal) × ((h-xd)/lscal..h/lscal)
    DrainInt = ((l+ld)/lscal..(L)/lscal) × ((h-xd)/lscal..h/lscal)
    GateInt = (0..L/lscal) × (0..h/lscal)
    GateMOS = MOSStructure(Metal(ϕₘ),SiO2(),SG,tox,h)
    Source = SemiconductorInterval(SD,SourceInt)
    Drain = SemiconductorInterval(SD,DrainInt)
    SourceContact = (0.0..0.8*ld/lscal)
    GateContact = (ld/lscal-gate_ov..(ld+l)/lscal+gate_ov)
    DrainContact = ((ld+l)/lscal+0.2*ld/lscal..L/lscal)
    MOSFETInputDeck(Source,GateMOS,Drain,GateInt,DrainContact,GateContact,SourceContact,l)
end

function C(x,y,fetIn::MOSFETInputDeck)
    Cv = [C(fetIn.Source.semiconductor),C(fetIn.Gate.semiconductor),C(fetIn.Drain.semiconductor)]
    return dot(Cv,belongVector(x,y,fetIn))
end

L(fetIN::MOSFETInputDeck) = fetIN.region.b[1]
H(fetIN::MOSFETInputDeck) = fetIN.region.b[2]
L_ch(fetIN::MOSFETInputDeck) = fetIN.Lch

function belongVector(x,y,fetIn::MOSFETInputDeck)
    if DomainSets.SVector((x,y)) ∈ fetIn.region
        if  belongs((x,y),fetIn.Source)
            return [1,0,0]
        elseif belongs((x,y),fetIn.Drain)
            return [0,0,1]
        else
            return [0,1,0]
        end
    else
        println("($(x),$(y)) does not belong to the Grid")
        return [0.0,0.0,0.0]
    end
end

function Grid(fetIn::MOSFETInputDeck,nx,ny,lscal)
    l = L(fetIn)/lscal
    h = H(fetIn)/lscal
    x = LinRange(0.0,l,nx)
    y = LinRange(0.0,h,ny)
    dx = diff(x)[1]
    dy = diff(y)[1]
    X = VoronoiFVM.Grid(x,y)
    VoronoiFVM.bfacemask!(X, [(fetIn.SourceContact.left)/lscal,h-dy/2.0],[fetIn.SourceContact.right/lscal,h+dy/2.0],6)
    VoronoiFVM.bfacemask!(X, [fetIn.GateContact.left/lscal,h-dy/2.0],[fetIn.GateContact.right/lscal,h+dy/2.0],7)
    VoronoiFVM.bfacemask!(X, [fetIn.DrainContact.left/lscal,h-dy/2.0],[fetIn.DrainContact.right/lscal,h+dy/2.0],8)

    return X

end

function MOSFETSimulation(fetIn::MOSFETInputDeck,Vg,Vd,nx,ny;verbose=false,ξ₀=1e-10)
    Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(fetIn.Gate.semiconductor,L(fetIn))
    Rf(n,p) = 0.001
    Cf(x,y) = C(x,y,fetIn)
    Jn = q*1400*Nv(fetIn.Source.semiconductor)*Vscal/Jscal
    Jp = q*734*Nc(fetIn.Source.semiconductor)*Vscal/Jscal
    G = Grid(fetIn,nx,ny,lscal)

    #u,nn,pp = DDVoronoiMOSFET(G,Ec(fetIn.Source.semiconductor),Ev(fetIn.Source.semiconductor),Nc(fetIn.Source.semiconductor),Nv(fetIn.Source.semiconductor),Ec(fetIn.Source.semiconductor)-Ec(fetIn.Gate.semiconductor),Vg,0.0,lscal,Rscal,Cscal,Vscal,λ,Jn,Jp,Rf,Cf,fetIn;verbose=verbose,ξ₀=ξ₀)
    u,nn,pp = DDVoronoiMOSFET(G,Ec(fetIn.Source.semiconductor),Ev(fetIn.Source.semiconductor),Nc(fetIn.Source.semiconductor),Nv(fetIn.Source.semiconductor),Ec(fetIn.Source.semiconductor)-Ec(fetIn.Gate.semiconductor),Vg,Vd,lscal,Rscal,Cscal,Vscal,λ,Jn,Jp,Rf,Cf,fetIn;verbose=verbose,ξ₀=ξ₀)
    return G[Coordinates].*lscal*1e4,Ec(fetIn.Source.semiconductor).-u,nn,pp
end

function MOSFETSimulation(N_d,Nb,T,l,h,Vg,Vd,Nx,Ny)
    SD = SemiconductorData(T,BoltzmanDist(),NSilicon(N_d))
    SG = SemiconductorData(T,BoltzmanDist(),PSilicon(Nb))
    #= SourceInt = (0..l/3.0) × (h/2.0..h)
    Source = SemiconductorInverval(SD,SourceInt)
    DrainInt = (2.0*l/3.0..l) × (h/2.0..h)
    Drain = SemiconductorInverval(SD,DrainInt)
    GateInt = setdiff((0..l)×(0..h),SourceInt)
    GateInt = setdiff(GateInt,DrainInt)
    Gate = SemiconductorInverval(SG,GateInt) =#

    #InputDeck = [Drain,Gate,Source]

    Vscal,lscal,Cscal,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(SG,l)
    Rf(n,p) = 0.0
    function Cf(x,y)
        if x < l/3.0 && y > 0.5h
            return C(SD)
        elseif x > 2.0l/3.0 && y > 0.5h
            return C(SD)
        else
            return C(SG)
        end
    end
    Plots.contour(LinRange(0,l,101),LinRange(0,h,101),Cf,fill=true)

    Jn = q*1400*Nv(SD)*Vscal/Jscal
    Jp = q*734*Nc(SD)*Vscal/Jscal

    x = LinRange(0,l/lscal,Nx)
    dx = diff(x)[1]
    y = LinRange(0,h/lscal,Ny)
    dy = diff(y)[1]
    X = VoronoiFVM.Grid(x,y)
    0.2l/lscal-dx
    1.1*h/lscal
    VoronoiFVM.bfacemask!(X, [0.0,h/lscal-dy/2.0],[0.2l/lscal,h/lscal+dy/2.0],6)
    VoronoiFVM.bfacemask!(X, [0.25*l/lscal,h/lscal-dy/2.0],[0.75l/lscal,h/lscal+dy/2.0],7)
    VoronoiFVM.bfacemask!(X, [0.8*l/lscal,h/lscal-dy/2.0],[1.1l/lscal,h/lscal+dy/2.0],8)

    u,nn,pp = DDVoronoiMOSFET(X,Ec(SD),Ev(SD),Nc(SD),Nv(SD),Ec(SD)-Ec(SG),Vg,0.0,lscal,Rscal,Cscal,Vscal,Ld0,l,λ,21,Jn,Jp,21,Rf,Cf)

    u,nn,pp = DDVoronoiMOSFET(X,Ec(SD),Ev(SD),Nc(SD),Nv(SD),Ec(SD)-Ec(SG),Vg,Vd,lscal,Rscal,Cscal,Vscal,Ld0,l,λ,21,Jn,Jp,21,Rf,Cf;u0=u/Vscal,n0=nn,p0=pp)
    Cords = X[Coordinates]
    y1 = unique(Cords[2,:])[2]
    mask = Cords[2,:] .== y1
    u0 = u[mask]
    u0 .-= u0[1]
    return (Cords,Ec(SD).-u0,nn*Cscal,pp*Cscal)
end


text:
"""
For running the fluorescence model. Uses the 'molecularODE' fuction with the Sundials ODE solver
to run the model.
Arguments:  spike_times = the frames when a spike occurs
            cell_radius = the modelled radius of the cell (in metres)
            baseline = the baseline amount of free calcium (M)
            calcium_rate = The rate at which calcium is kicked out of the cell (s^-1)
            indicator = the concentration of indicator in the cell (M)
            endogeneous = the concentration of endogeneous mobile buffer in the cell (M)
            immobile = the concentration of endogeneous immobile buffer in the cell (M)
            b_i = the backward rate of the indicator binding (s^-1)
            f_i = the forward rate of the indicator binding (M^-1 s^-1)
            b_e = the backward rate of the endogeneous binding (s^-1)
            f_e = the forward rate of the endogeneous binding (M^-1 s^-1)
            b_im = the backward rate of the immobile binding (s^-1)
            f_im = the forward rate of the immobile binding (M^-1 s^-1)
            excitation = the proportion of BCa that gets excited at each time step (s^-1)
            release = the proportion of BCa* that releases a photon at each time step (s^-1)
            peak = the amount of free calcium in the cell after a spike (M)
            frequency = sampling frequency (s^-1)
            capture_rate = the probability of capturing a given photon
Returns:    fluorescence = Array{Float64, 1}, the modelled fluorescence
            spiking_sim = T x 5 Array{Float, 2}, the full model
            spiking_sim_time = Array{Float, 1} length T, the time points
"""
function calciumFluorescenceModel(spike_train::Array{Int64,1};  cell_radius::Float64=1e-5,
                                  baseline::Float64=0.045e-6,   calcium_rate::Float64=0.001,
                                  indicator::Float64=1e-4,      endogeneous::Float64=1e-4,
                                  immobile::Float64=78.7e-6,    b_i::Float64=160.0,         f_i::Float64=7.76699e8,
                                  b_e::Float64=10000.0,         f_e::Float64=100e6,         b_im::Float64=524.0,
                                  f_im::Float64=2.47e8,         excitation::Float64=0.15,   release::Float64=0.11,
                                  peak::Float64=2.9e-7,         frequency::Float64=100.0,
                                  capture_rate::Float64=0.62)
  cell_volume, m_baseline, m_indicator, m_endogeneous, m_immobile, m_peak, m_f_i, m_f_e, m_f_im = convertUnits(cell_radius::Float64, baseline::Float64, indicator::Float64, endogeneous::Float64, immobile::Float64, peak::Float64, f_i::Float64, f_e::Float64, f_im::Float64)
  model_time, spike_indices, model_indices = getModelTime(spike_train, frequency)
  # baseline dynamics for 25000 seconds to reach baseline conditions
  # The ODEs are solved for every 100 seconds
  baseline_time = collect(0:100.0:25000)
  baseline_initial_conditions = [m_baseline, 0.0, 0.0, 0.0, 0.0]
  baseline_sim = Sundials.cvode((t, x, xdot) -> molecularODE(t, x, xdot, m_baseline, calcium_rate, m_indicator, m_endogeneous, m_immobile, b_i, m_f_i, b_e, m_f_e, b_im, m_f_im, excitation, release), baseline_initial_conditions, baseline_time)
  # pre-spiking dynamics, ODEs are solved for every millisecond
  if spike_indices[1] > 1 # if there is no spike in the first time bin
    pre_spike_time = model_time[1:spike_indices[1]-1]
    pre_spike_initial_conditions = baseline_sim[end,:]
    pre_spike_sim = Sundials.cvode((t, x, xdot) -> molecularODE(t, x, xdot, m_baseline, calcium_rate, m_indicator, m_endogeneous, m_immobile, b_i, m_f_i, b_e, m_f_e, b_im, m_f_im, excitation, release), pre_spike_initial_conditions, pre_spike_time)
    # spiking dynamics, discontinuous jump at each spike
    spiking_sim = pre_spike_sim
    spiking_sim_time = pre_spike_time
    for sp in 1:length(spike_indices)-1
      spike_time = model_time[spike_indices[sp]:spike_indices[sp+1]-1]
      spiking_initial_conditions = [m_peak; spiking_sim[end, 2:end]]
      spike_sim = Sundials.cvode((t, x, xdot) -> molecularODE(t, x, xdot, m_baseline, calcium_rate, m_indicator, m_endogeneous, m_immobile, b_i, m_f_i, b_e, m_f_e, b_im, m_f_im, excitation, release), spiking_initial_conditions, spike_time)
      spiking_sim = [spiking_sim; spike_sim]
      spiking_sim_time = [spiking_sim_time; spike_time]
    end
  else # if there is a spike in the first time bin
    spiking_sim = zeros(Float64,(0,5))
    spiking_sim_time = zeros(Float64, 0)
    spiking_initial_conditions = [m_peak; baseline_sim[end, 2:end]]
    for sp in 1:length(spike_indices)-1
      spike_time = model_time[spike_indices[sp]:spike_indices[sp+1]-1]
      spike_sim = Sundials.cvode((t, x, xdot) -> molecularODE(t, x, xdot, m_baseline, calcium_rate, m_indicator, m_endogeneous, m_immobile, b_i, m_f_i, b_e, m_f_e, b_im, m_f_im, excitation, release), spiking_initial_conditions, spike_time)
      spiking_sim = [spiking_sim; spike_sim]
      spiking_sim_time = [spiking_sim_time; spike_time]
      spiking_initial_conditions = [m_peak; spiking_sim[end, 2:end]]
    end
  end
  # dynamics for the last spike until the end of the trail
  last_spike_time = model_time[spike_indices[end]:end]
  last_spike_initial_conditions = [m_peak; spiking_sim[end, 2:end]]
  last_spike_sim = Sundials.cvode((t, x, xdot) -> molecularODE(t, x, xdot, m_baseline, calcium_rate, m_indicator, m_endogeneous, m_immobile, b_i, m_f_i, b_e, m_f_e, b_im, m_f_im, excitation, release), last_spike_initial_conditions, last_spike_time)
  spiking_sim = [spiking_sim; last_spike_sim]
  spiking_sim_time = [spiking_sim_time; last_spike_time]
  # photon capturing
  fluorescence = [getPhotonsCollected(ebc, release, capture_rate) for ebc in spiking_sim[model_indices, end]]
  fluorescence = shiftFluorescence(fluorescence)
  spiking_sim_molars = moleculesToMolars(spiking_sim, cell_volume)
  return fluorescence, spiking_sim_molars[model_indices,:], spiking_sim_time[model_indices]
end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d

include("../utils/ortho_gram_schmidt.jl")
include("orthonormalize.jl")
include("test_ortho.jl")

function test_main()
  NN = [25; 25; 25]
  AA = [0.0; 0.0; 0.0]
  BB = [6.0; 6.0; 6.0]

  Npoints = prod(NN)

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB )
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  srand(1234)
  v = rand(Npoints,4)
  orthonormalize!(LF,v)
  test_ortho( ΔV, 4, v )
end

test_main()


text:
function ScanningStrategy2MapInfo(SS::ScanningStrategy, division::Int)
    """
    Nonte:
        This function assumed 2 detectors which is seeing same polarization angle.
    """
    resol = Resolution(SS.nside)
    npix = nside2npix(SS.nside)

    month = Int(SS.duration / division)
    ω_hwp = rpm2angfreq(SS.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    hit_matrix = zeros(npix, 3, 3)
    detmap = zeros(npix)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod, psi_tod, time_array = get_pointing_pixels(SS, BEGIN, END)
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            #psi_tod_jth_det = psi_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array[k]
                ipix = pix_tod_jth_det[k]
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 2 # Assumed 2 detectors
                Cross[ipix,1,1] += 2sin(hwp_ang - psi) # Assumed 2 detectors
                Cross[ipix,1,2] += 2cos(hwp_ang - psi)
                Cross[ipix,2,1] += 2sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += 2cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += 2sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += 2cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += 2sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += 2cos(hwp_ang - 4psi)
                hit_matrix[ipix,:,:] .+= 2HitMatrix(hwp_ang - 2psi)
            end
        end
        BEGIN = END
        next!(p)
    end
    @views @inbounds @threads for i in eachindex(hit_map)
        Inv_hitmat = inv.(hit_matrix[i,:,:])
        hit_matrix[i,:,:] .= Inv_hitmat
        detmap[i] = det(Inv_hitmat)./hit_map[i]
    end
    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    outmap = Dict{String,AbstractArray{Float64}}(
        "hitmap" => hit_map,
        "xl1" => link1,
        "xl2" => link2,
        "xl3" => link3,
        "xl4" => link4,
        "detM" => detmap,
        "dI" => hit_matrix[:,1,1],
        "dQ" => hit_matrix[:,2,2],
        "dU" => hit_matrix[:,3,3]
        )
    return outmap
end

function get_psiDataBase(SS::ScanningStrategy,; division::Int, idx, map_div)
    #=
    divisionはtodの分割計算
    idxはマップ(npix)をRingOrderでmap_div*nside個の領域に分割したときに北極から順に貼られるインデックス
    map_divはマップをmap_div*nside個に分割するためのパラメータ(4くらいがいい)
    =#

    resol = Resolution(SS.nside)
    npix = nside2npix(SS.nside)
    map_division = map_div*SS.nside
    println("To get a psi database in full-sky, specify idx in the range of [1, $map_division] and run the job.")

    split = npix/map_division
    month = Int(SS.duration / division)
    ω_hwp = rpm2angfreq(SS.hwp_rpm)

    psi_db = [Float64[] for i in 1:split]
    under = Int(split*(idx-1))
    upper = Int(split*idx)
    println("ipix_range = [$(under+1),$upper]")

    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod, psi_tod, time_array = get_pointing_pixels(SS, BEGIN, END)
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array[k]
                ipix = pix_tod_jth_det[k]
                psi = 2ω_hwp*t - psi_tod_jth_det[k]

                if under < ipix <= upper
                    push!(psi_db[ipix-under], psi)
                end
            end
        end
        BEGIN = END
        next!(p)
    end
    return (map_div, idx, psi_db)
end

function TwoTelescopes_ScanningStrategy2map(SS1::ScanningStrategy, SS2::ScanningStrategy, division::Int)
    resol = Resolution(SS1.nside)
    npix = nside2npix(SS1.nside)

    month = Int(SS1.duration / division)
    ω_hwp = rpm2angfreq(SS1.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    hit_matrix = zeros(npix, 3, 3)
    detmap = zeros(npix)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod1, psi_tod1, time_array1 = get_pointing_pixels(SS1, BEGIN, END)
        pix_tod2, psi_tod2, time_array2 = get_pointing_pixels(SS2, BEGIN, END)
        pix_tod = @views [pix_tod1 pix_tod2]
        psi_tod = @views [psi_tod1 psi_tod2]
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            #psi_tod_jth_det = psi_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array1[k]
                ipix = pix_tod_jth_det[k]
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 1
                Cross[ipix,1,1] += sin(hwp_ang - psi)
                Cross[ipix,1,2] += cos(hwp_ang - psi)
                Cross[ipix,2,1] += sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += cos(hwp_ang - 4psi)
                hit_matrix[ipix,:,:] .+= HitMatrix(hwp_ang - 2psi)
            end
        end
        BEGIN = END
        next!(p)
    end

    @views @inbounds @threads for i = 1:npix
        detmap[i] = det((hit_matrix[i,:,:])./hit_map[i])
        hit_matrix[i,:,:] .= inv.(hit_matrix[i,:,:])
    end
    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    out_map = @views [hit_map, link1, link2, link3, link4, detmap]
    return out_map
end

function TwoTelescopes_ScanningStrategy2MapInfo(SS1::ScanningStrategy, SS2::ScanningStrategy, division::Int)
    resol = Resolution(SS1.nside)
    npix = nside2npix(SS1.nside)

    month = Int(SS1.duration / division)
    ω_hwp = rpm2angfreq(SS1.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    hit_matrix = zeros(npix, 3, 3)
    detmap = zeros(npix)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod1, psi_tod1, time_array1 = get_pointing_pixels(SS1, BEGIN, END)
        pix_tod2, psi_tod2, time_array2 = get_pointing_pixels(SS2, BEGIN, END)
        pix_tod = @views [pix_tod1 pix_tod2]
        psi_tod = @views [psi_tod1 psi_tod2]
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            #psi_tod_jth_det = psi_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array1[k]
                ipix = pix_tod_jth_det[k]
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 2
                Cross[ipix,1,1] += 2sin(hwp_ang - psi)
                Cross[ipix,1,2] += 2cos(hwp_ang - psi)
                Cross[ipix,2,1] += 2sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += 2cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += 2sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += 2cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += 2sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += 2cos(hwp_ang - 4psi)
                hit_matrix[ipix,:,:] .+= 2HitMatrix(hwp_ang - 2psi)
            end
        end
        BEGIN = END
        next!(p)
    end

    @views @inbounds @threads for i = 1:npix
        detmap[i] = det((hit_matrix[i,:,:])./hit_map[i])
        hit_matrix[i,:,:] .= inv.(hit_matrix[i,:,:])
    end

    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    outmap = Dict{String,AbstractArray{Float64}}(
        "hitmap" => hit_map,
        "xl1" => link1,
        "xl2" => link2,
        "xl3" => link3,
        "xl4" => link4,
        "detM" => detmap,
        "dI" => hit_matrix[:,1,1],
        "dQ" => hit_matrix[:,2,2],
        "dU" => hit_matrix[:,3,3]
        )
    return outmap
end

function ThreeTelescopes_ScanningStrategy2map(SS1::ScanningStrategy, SS2::ScanningStrategy, SS3::ScanningStrategy, division::Int)
    resol = Resolution(SS1.nside)
    npix = nside2npix(SS1.nside)

    month = Int(SS1.duration / division)
    ω_hwp = rpm2angfreq(SS1.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    hit_matrix = zeros(npix, 3, 3)
    detmap = zeros(npix)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        pix_tod1, psi_tod1, time_array1 = get_pointing_pixels(SS1, BEGIN, END)
        pix_tod2, psi_tod2, time_array2 = get_pointing_pixels(SS2, BEGIN, END)
        pix_tod3, psi_tod3, time_array3 = get_pointing_pixels(SS3, BEGIN, END)
        pix_tod = @views [pix_tod1 pix_tod2 pix_tod3]
        psi_tod = @views [psi_tod1 psi_tod2 psi_tod3]
        @views @inbounds for j = eachindex(psi_tod[1,:])
            pix_tod_jth_det = pix_tod[:,j]
            #psi_tod_jth_det = psi_tod[:,j]
            psi_tod_jth_det = ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(psi_tod[:,1])
                t = time_array1[k]
                ipix = pix_tod_jth_det[k]
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 1
                Cross[ipix,1,1] += sin(hwp_ang - psi)
                Cross[ipix,1,2] += cos(hwp_ang - psi)
                Cross[ipix,2,1] += sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += cos(hwp_ang - 4psi)
                hit_matrix[ipix,:,:] .+= HitMatrix(hwp_ang - 2psi)
            end
        end
        BEGIN = END
        next!(p)
    end
    @views @inbounds @threads for i = 1:npix
        detmap[i] = det((hit_matrix[i,:,:])./hit_map[i])
    end
    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    out_map = @views [hit_map, link1, link2, link3, link4, detmap]
    return out_map
end

function Mapmaking(SS::ScanningStrategy, division::Int)
    resol = Resolution(SS.nside)
    npix = nside2npix(SS.nside)

    month = Int(SS.duration / division)
    ω_hwp = rpm2angfreq(SS.hwp_rpm)

    hit_map = zeros(npix)
    Cross = zeros(npix, 4, 2)
    BEGIN = 0
    p = Progress(division)
    @views @inbounds for i = 1:division
        END = i * month
        theta_tod, phi_tod, psi_tod, time_array = get_pointings_tuple(SS, BEGIN, END)
        @inbounds for j = eachindex(psi_tod[1,:])
            theta_tod_jth_det = @views theta_tod[:,j]
            phi_tod_jth_det = @views phi_tod[:,j]
            psi_tod_jth_det = @views ifelse(ω_hwp == 0.0, -psi_tod[:,j], psi_tod[:,j])
            @views @inbounds @simd for k = eachindex(time_array)
                t = @views time_array[k]
                ipix = ang2pixRing(resol, theta_tod_jth_det[k], phi_tod_jth_det[k])
                psi = psi_tod_jth_det[k]
                hwp_ang = 4ω_hwp*t

                hit_map[ipix] += 1
                Cross[ipix,1,1] += sin(hwp_ang - psi)
                Cross[ipix,1,2] += cos(hwp_ang - psi)
                Cross[ipix,2,1] += sin(hwp_ang - 2psi)
                Cross[ipix,2,2] += cos(hwp_ang - 2psi)
                Cross[ipix,3,1] += sin(hwp_ang - 3psi)
                Cross[ipix,3,2] += cos(hwp_ang - 3psi)
                Cross[ipix,4,1] += sin(hwp_ang - 4psi)
                Cross[ipix,4,2] += cos(hwp_ang - 4psi)
            end
        end
        BEGIN = END
        next!(p)
    end
    link1 = @views @. (Cross[:,1,1]/hit_map)^2 + (Cross[:,1,2]/hit_map)^2
    link2 = @views @. (Cross[:,2,1]/hit_map)^2 + (Cross[:,2,2]/hit_map)^2
    link3 = @views @. (Cross[:,3,1]/hit_map)^2 + (Cross[:,3,2]/hit_map)^2
    link4 = @views @. (Cross[:,4,1]/hit_map)^2 + (Cross[:,4,2]/hit_map)^2
    out_map = @views [hit_map, link1, link2, link3, link4]
    return out_map
end


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

using Base.Test

script = """
# Issue #11948
f(x) = x+1
workspace()
@assert @__MODULE__() === Main
@assert isdefined(Main, :f)
@assert !@isdefined LastMain
@eval Core.Main begin
    @assert @__MODULE__() === Main
    @assert !isdefined(Main, :f)
    LastMain.f(2)

    # PR #12990
    io = IOBuffer()
    show(io, Pair)
    @assert String(take!(io)) == "Pair"
    @assert !Base.inbase(LastMain)
end
"""
exename = Base.julia_cmd()
@test success(pipeline(`$exename --startup-file=no -e $script`, stdout=STDOUT, stderr=STDERR))

# issue #17764
script2 = """
mutable struct Foo end
workspace()
@eval Core.Main begin
    mutable struct Foo end
    @assert Tuple{Type{LastMain.Foo}} !== Tuple{Type{Main.Foo}}
end
"""
@test success(pipeline(`$exename --startup-file=no -e $script2`, stdout=STDOUT, stderr=STDERR))

# Issue #22101
mktempdir() do dir
    withenv("JULIA_DEBUG_LOADING" => nothing) do
        # We need to ensure that the module does a nontrivial amount of work during precompilation
        write(joinpath(dir, "Test22101.jl"), """
            __precompile__()
            module Test22101
                export f22101
                f22101() = collect(1:10)
                f22101()
            end
        """)
        write(joinpath(dir, "testdriver.jl"), """
            using Base.Test
            insert!(LOAD_PATH, 1, $(repr(dir)))
            insert!(Base.LOAD_CACHE_PATH, 1, $(repr(dir)))
            @test !isdefined(Main, :f22101)
            @test !isdefined(Main, :LastMain)
            begin
                using Test22101
                @test isdefined(Main, :f22101)
                @test f22101()::Vector{Int} == collect(1:10)
                @eval workspace() using Test22101
                @test f22101()::Vector{Int} == collect(1:10)
                @test isdefined(Main, :f22101)
            end
            @test isdefined(Main, :f22101)
            @test !isdefined(Main, :LastMain)
            @test isdefined(Core.Main, :f22101)
            nothing
        """)
        # Ensure that STDIO doesn't get swallowed (helps with debugging)
        cmd = `$(Base.julia_cmd()) --startup-file=no --precompiled=yes --compilecache=yes $(joinpath(dir, "testdriver.jl"))`
        @test success(pipeline(cmd, stdout=STDOUT, stderr=STDERR))
    end
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
f(a,b) = 333.75*b^6+a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b)

a=77617.0
b=33096.0

println("Float64: ", f(a, b))
println("BigFloat: ", f(BigFloat(a), BigFloat(b)))


text:
module SimGBS

export restrictionEnzyme,
    digestGenome,           # main function for in slico digestion using restriction enzyme
    definePopulation,       # main function for defining population structure
    GBS                    # main function for GBS simulation

using Compat,
    DelimitedFiles,         # for reading and writing text-delimited files
    GZip,                   # for reading and writing gzipped files
    Random,                 # for random permutation
    Distributions,          # for sampling distribution
    StatsBase,              # for weighted vectors
    Plots,                  # for plotting
    LsqFit  # fitting non-linear models

include("genome.jl")          # step 1. in slico digestion using restriction enzyme
include("pop.jl")             # step 2. set up population structure
include("gbs.jl")             # step 3. simultae GBS data

end


text:
using ECHOsched
using Test

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

pth="C:/Users/arch/Dropbox/5-Data/DATA/EnergyMeter/UKDS-7591/DwellingDatasets"
query="GasDemand_kWhPer30min"
N=5

@time InitProject(;N=1:5,pth=pth,query=query)
@code_warntype InitProject(;N=1,pth=pth,query=query)

ScheduleExtraction(;verboseOutput=true)


text:
# # Section Allocation
# Suppose you have $n$ students in a class who need to be assigned to $m$
# discussion sections. Each student needs to be assigned to exactly one section.
# Each discussion section should have between 6 and 10 students. Suppose an
# $n \times m$ preference matrix $P$ is given, where $P_{ij}$ gives student
# $i$'s ranking for section $j$ (1 would mean it is the student's top choice,
# 10,000 or a large number would mean the student can not attend that section).
#
# The goal will be to get an allocation matrix $X$, where $X_{ij} = 1$ if
# student $i$ is assigned to section $j$ and $0$ otherwise.

using Convex, GLPK
aux(str) = joinpath(@__DIR__, "aux_files", str) # path to auxiliary files

# Load our preference matrix, `P`
include(aux("data.jl"))

X = Variable(size(P), :Bin)

# We want every student to be assigned to exactly one section. So, every row
# must have exactly one non-zero entry. In other words, the sum of all the
# columns for every row is 1. We also want each section to have between 6 and 10
# students, so the sum of all the rows for every column should be between these.
constraints = [sum(X, dims=2) == 1, sum(X, dims=1) <= 10, sum(X, dims=1) >= 6]

# Our objective is simple `sum(X .* P)`, which can be more efficiently
# represented as `vec(X)' * vec(P)`. Since each entry of `X` is either 0 or 1,
# this is basically summing up the rankings of students that were assigned to them.
# If all students got their first choice, this value will be the number of
# students since the ranking of the first choice is 1.
p = minimize(vec(X)' * vec(P), constraints)

solve!(p, GLPK.Optimizer)
p.optval


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

mutable struct BackendAddressPoolPropertiesFormat2 <: SwaggerModel
    backendIPConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration} } # spec name: backendIPConfigurations
    loadBalancerBackendAddresses::Any # spec type: Union{ Nothing, Vector{LoadBalancerBackendAddress2} } # spec name: loadBalancerBackendAddresses
    loadBalancingRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: loadBalancingRules
    outboundRule::Any # spec type: Union{ Nothing, SubResource } # spec name: outboundRule
    outboundRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: outboundRules
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function BackendAddressPoolPropertiesFormat2(;backendIPConfigurations=nothing, loadBalancerBackendAddresses=nothing, loadBalancingRules=nothing, outboundRule=nothing, outboundRules=nothing, provisioningState=nothing)
        o = new()
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("backendIPConfigurations"), backendIPConfigurations)
        setfield!(o, Symbol("backendIPConfigurations"), backendIPConfigurations)
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("loadBalancerBackendAddresses"), loadBalancerBackendAddresses)
        setfield!(o, Symbol("loadBalancerBackendAddresses"), loadBalancerBackendAddresses)
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("loadBalancingRules"), loadBalancingRules)
        setfield!(o, Symbol("loadBalancingRules"), loadBalancingRules)
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("outboundRule"), outboundRule)
        setfield!(o, Symbol("outboundRule"), outboundRule)
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("outboundRules"), outboundRules)
        setfield!(o, Symbol("outboundRules"), outboundRules)
        validate_property(BackendAddressPoolPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type BackendAddressPoolPropertiesFormat2

const _property_map_BackendAddressPoolPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("backendIPConfigurations")=>Symbol("backendIPConfigurations"), Symbol("loadBalancerBackendAddresses")=>Symbol("loadBalancerBackendAddresses"), Symbol("loadBalancingRules")=>Symbol("loadBalancingRules"), Symbol("outboundRule")=>Symbol("outboundRule"), Symbol("outboundRules")=>Symbol("outboundRules"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_BackendAddressPoolPropertiesFormat2 = Dict{Symbol,String}(Symbol("backendIPConfigurations")=>"Vector{NetworkInterfaceIPConfiguration}", Symbol("loadBalancerBackendAddresses")=>"Vector{LoadBalancerBackendAddress2}", Symbol("loadBalancingRules")=>"Vector{SubResource}", Symbol("outboundRule")=>"SubResource", Symbol("outboundRules")=>"Vector{SubResource}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ BackendAddressPoolPropertiesFormat2 }) = collect(keys(_property_map_BackendAddressPoolPropertiesFormat2))
Swagger.property_type(::Type{ BackendAddressPoolPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BackendAddressPoolPropertiesFormat2[name]))}
Swagger.field_name(::Type{ BackendAddressPoolPropertiesFormat2 }, property_name::Symbol) =  _property_map_BackendAddressPoolPropertiesFormat2[property_name]

function check_required(o::BackendAddressPoolPropertiesFormat2)
    true
end

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


text:
using LocalCoverage, Test

lockfile = "/tmp/testingLocalCoverage" # prevent infinite recursion when testing

covdir = normpath(joinpath(@__DIR__, "..", "coverage"))

if !isfile(lockfile)
    @test isdir(LocalCoverage.pkgdir("LocalCoverage"))
    tracefile = joinpath(covdir, "lcov.info")
    @test !isfile(tracefile)
    touch(lockfile)
    cov = generate_coverage("LocalCoverage"; genhtml=true, show_summary=false)
    tab = LocalCoverage.coverage_summary_table(cov)
    @test tab[1,1] == "src/LocalCoverage.jl"
    LocalCoverage.coverage_summary(cov)
    rm(lockfile)
    @test isfile(tracefile)
    clean_coverage("LocalCoverage", rm_directory=true)
    @test !isdir(covdir)
end


text:
@testset "Test Conversion" begin
    @info "Test: Ellipse Parameter Conversion"

    @testset "From Geometric to Algebraic and back" begin
        𝛏₀ = [10, 5, 25, 25, π/4]
        geo_to_alg = GeometricToAlgebraic()
        alg_to_geo = AlgebraicToGeometric()
        𝛉 = geo_to_alg(𝛏₀)
        𝛏₁ = alg_to_geo(𝛉)
        @test all(𝛏₀ .≈ 𝛏₁)
    end

    @testset "From Algebraic to Geometric and back" begin
        𝛉₀ = [0.025, -0.03, 0.025, -0.5, -0.5, 11.5]
        geo_to_alg = GeometricToAlgebraic()
        alg_to_geo = AlgebraicToGeometric()
        𝛏 = alg_to_geo(𝛉₀)
        𝛉₁ = geo_to_alg(𝛏)
        @test all(𝛉₀ .≈ 𝛉₁)
    end

    @testset "From  Algebraic to Latent Algebraic and back" begin
        𝛏₀ = [10, 5, 25, 25, π/4]
        geo_to_alg = GeometricToAlgebraic()
        𝛉₀ = geo_to_alg(𝛏₀)
        𝛉₀ = 𝛉₀ / norm(𝛉₀)

        𝛋 = GuaranteedEllipseFitting.to_latent_parameters(GuaranteedEllipseFitting.FirstLatentEllipseParametrisation(), 𝛉₀)
        𝛉₁ = GuaranteedEllipseFitting.from_latent_parameters(GuaranteedEllipseFitting.FirstLatentEllipseParametrisation(), 𝛋)
        @test all(𝛉₀ .≈ 𝛉₁)

        𝛋 = GuaranteedEllipseFitting.to_latent_parameters(GuaranteedEllipseFitting.SecondLatentEllipseParametrisation(), 𝛉₀)
        𝛉₁ = GuaranteedEllipseFitting.from_latent_parameters(GuaranteedEllipseFitting.SecondLatentEllipseParametrisation(), 𝛋)
        @test all(𝛉₀ .≈ 𝛉₁)
    end

end


text:
@generated function backpropagationoptimize{F}(mlmodel::MLAlgorithm{F}, input::AbstractArray, answers::AbstractArray, cost::Function)

  #check to make sure the mlmodel has backpropagation
  if !hasbackpropagation(mlmodel)
    return :(throw(ArgumentError("mlmodel $(typeof(mlmodel)) doesn't support backpropagation")))
  end

  #for analysis purposes, allow the backpropagation of a single layer.
  if mlmodel <: Layer
    quote
      #TODO:  Add support for this.
      throw(ArgumentError("single layer backprop has not been implemented yet."))
    end
  else
    quote
      warn("memory optimized version of backpropagation has not been implemented yet")

      #generate a storage unit for the mlmodel
      storage = BackpropStorage(mlmodel)
      #generate a results unit for the mlmodel
      res = Vector{F}(outputs(mlmodel))

      evaluate!(res, mlmodel, input, storage)

      res_deltas = (d(cost))(answers, res)

      backpropagate!(mlmodel, input, res, res_deltas, storage)
    end
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
@testset "basics" begin

	bc = rand(2,3)
	bcde = rand(2,3,4,5);

	@reduce C[c] := sum(b) bc[b,c]
	@test C == vec(sum(bc, dims=1))

	@reduce B[b] := prod(c:3) bc[b,c]
	@test B == vec(prod(bc, dims=2))

	using Statistics

	@reduce W[b⊗d] := Statistics.mean(c,e) bcde[b,c,d,e]
	V = vec(mean(bcde, dims=(2,4)))
	@test V == W

	@reduce Z[e,b,β] := Statistics.std(c,α) bcde[b,c,α⊗β,e]  (α in 1:2)
	@test size(Z) == (5,2,2)

	@reduce A[_,e,_,b] := sum(c) bcde[b,c,3,e]
	@test size(A) == (1,5,1,2)
	B = sum(bcde[:,:,3,:], dims=2)
	@test A[1,3,1,2] == B[2,1,3]
	C = similar(A);
	@reduce C[_,e,_,b] = sum(c) bcde[b,c,3,e]
	@test A ≈ C

	@reduce D[b][c] := sum(d,e) bcde[b,c,d,e] + 1
	@reduce E[b][c] := sum(e,d) bcde[b,c,d,e] + 1
	@reduce F[c][b] := sum(d:4,e:5) bcde[b,c,d,e] + 1 # doesn't warn
    @reduce F′[c][b] := sum(d,e) bcde[b,c,d,e] + 1  d in 1:4, e in 1:5

	@test D[1][2] ≈ sum(bcde, dims=(3,4))[1,2, 1,1] + 4*5
	@test D[1][2] ≈ E[1][2] ≈ F[2][1]
    @test F′ ≈ F

end
@testset "scalar" begin

	bc = rand(2,3)

	@reduce S[] := sum(i,j) bc[j,i]
	T = @reduce sum(i,j) bc[j,i]
	U = @reduce sum(j,i) bc[j,i]
	@test S[] ≈ sum(bc) ≈ T ≈ U

end
@testset "inference" begin

    # inference for a⊗b⊗c had an (Any[]...) dots problem at first
    B = randn(8,24);
    @reduce A[b,c, y,z] := sum(a:2, x:2) B[a⊗b⊗c, x⊗y⊗z]  b in 1:2, y in 1:3
    @test size(A) == (2,2, 3,4)

    C = similar(A)
    @reduce C[b,c, y,z] = sum(a, x) B[a⊗b⊗c, x⊗y⊗z]
    @test C ≈ A

    # with a:2 given, it doesn't try product inference, just leaves a :, which is OK.
    @reduce C[b,c, y,z] = sum(a:2, x) B[a⊗b⊗c, x⊗y⊗z]
    @test C ≈ A

end
@testset "recursion" begin

    # from readme
    A = rand(4);
    B = randn(4,4);
    R = @reduce sum(i) A[i] * log( @reduce _[i] := sum(j) A[j] * exp(B[i,j]) )
    @reduce inner[i] := sum(j) A[j] * exp(B[i,j])
    S = @reduce sum(i) A[i] * log(inner[i])
    @test S == R

    R2 = @reduce sum(i) A[i] * log( @reduce inner[i] = sum(j) A[j] * exp(B[i,j]) )
    @test S == R2

    # scalar result inside
    tot = @reduce sum(i) A[i]
    @cast N[i] := A[i] / tot
    @cast M[i] := A[i] / @reduce sum(i) A[i]
    @test N == M

    # this needs size of the result of inner macro, failed at first:
    A = rand(2,3); B = rand(3,4); C = rand(4,5);
    @reduce W[i,_,l] := sum(j) A[i,j] * (@matmul _[j,l] := sum(k) B[j,k] * C[k,l])
    @test W ≈ reshape(A * B * C, 2,1,5)

end


text:
# Basic library functionality

#
# API call wrapper
#

macro apicall(fun, rettyp, argtypes, args...)
    isa(fun, QuoteNode) || error("first argument to @apicall should be a symbol")

    target = if startswith(String(fun.value), "LLVMExtra")
        fun
    else
        :($fun, libllvm)
    end

    return quote
        ccall($target, $(esc(rettyp)), $(esc(argtypes)), $(map(esc, args)...))
    end
end


text:
using Simulacrum
using HyperGraphs, Symbolics

@variables t X(t) Y(t) Z(t)
che1 = ChemicalHyperEdge([X], [Y], 2)		# chemical hyperedge unimolecular
che2 = ChemicalHyperEdge([X, Y], [Z], 2)	# chemical hyperedge multimolecular
@test isunimolecular(che1) && ismultimolecular(che2)
@test !isunimolecular(che2) && !ismultimolecular(che1)

# reversible reactions
@variables t A(t) B(t) C(t) k₁ k₂
fwd_bwd = reversible_reaction([A, B], [C], [k₁, k₂])
fwd = ChemicalHyperEdge([A, B], [C], k₁)
bwd = ChemicalHyperEdge([C], [A, B], k₂)
@test [fwd, bwd] == fwd_bwd
@test reversible_reaction(ChemicalHyperEdge([A, B], [C]), [k₁, k₂])[2] == bwd
@test (@reversible ChemicalHyperEdge([A, B], [C]) [k₁, k₂]) == [fwd, bwd]
@test reversible_reaction([A, B], [C], [k₁, k₂])[2] == bwd


text:
import StaticArrays; const SA = StaticArrays
setindex(a::SA.StaticArray, args...) = Base.setindex(a, args...)

set(arr, ::typeof(SA.normalize), val) = SA.norm(arr) * val
set(arr, ::typeof(SA.norm), val)      = val/SA.norm(arr) * arr # shoud we check val is positive?


text:
using ConstrainedDynamics
using ConstrainedDynamics: vrotate, Vmat
using Rotations
using StaticArrays
using LinearAlgebra

Δt = 0.01
length1 = 1.0
width, depth = 1.0, 1.0
box1 = Box(width, depth, length1, length1, color = RGBA(1., 1., 0.))

for i=1:10
    # Links
    origin = Origin{Float64}()
    link1 = Body(box1)
    link1.m = 1.0
    link1.J = diagm(ones(3))

    axis = rand(3)
    axis = axis/norm(axis)
    qoff = one(UnitQuaternion) # rand(UnitQuaternion)

    # Constraints
    joint1 = EqualityConstraint(Revolute(origin, link1, axis; qoffset = qoff))

    links = [link1]
    constraints = [joint1]
    shapes = [box1]

    mech = Mechanism(origin, links, constraints, g = 0., shapes = shapes)

    xθ = rand(1)
    vω = rand(1)
    Fτ = rand(1)

    setPosition!(mech,joint1,xθ;iter=false)
    setVelocity!(mech,joint1,vω)
    function control!(mechanism, k)
        if k==1
            setForce!(mech,joint1,Fτ)
        end
        return
    end
    # setForce!(mech,joint1,Fτ)

    storage = simulate!(mech, 10., control!, record = true)

    angend = mod((xθ + (vω + Fτ*Δt)*10.0)[1],2pi)
    qend = qoff*UnitQuaternion(cos(angend/2), (axis*sin(angend/2))..., false)
    minresult = mod(minimalCoordinates(mech, joint1)[1],2pi)

    @test isapprox(norm(minresult - angend), 0.0; atol = 1e-3)
    @test isapprox(norm(link1.state.xc), 0.0; atol = 1e-3)
    @test isapprox(norm(Vmat(link1.state.qc/qend)), 0.0; atol = 1e-3)
end


text:
# ---
# title: JuliaRL\_MADDPG\_SpeakerListener
# cover: assets/JuliaRL_MADDPG_SpeakerListenerEnv.png
# description: MADDPG applied to SpeakerListenerEnv
# date: 2021-08-28
# author: "[Peter Chen](https://github.com/peterchen96)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Statistics
using Flux

mutable struct MeanRewardHook <: AbstractHook
    episode::Int
    eval_rate::Int
    eval_episode::Int
    episodes::Vector
    mean_rewards::Vector
end

function (hook::MeanRewardHook)(::PostEpisodeStage, policy, env)
    if hook.episode % hook.eval_rate == 0
        ## evaluate policy's performance
        rew = 0
        for _ in 1:hook.eval_episode
            reset!(env)
            while !is_terminated(env)
                env |> policy |> env
                rew += reward(env)
            end
        end

        push!(hook.episodes, hook.episode)
        push!(hook.mean_rewards, rew / hook.eval_episode)
    end
    hook.episode += 1
end

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:MADDPG},
    ::Val{:SpeakerListener},
    ::Nothing;
    seed=123,
)
    rng = StableRNG(seed)
    env = SpeakerListenerEnv(max_steps = 25)

    init = glorot_uniform(rng)
    critic_dim = sum(length(state(env, p)) + length(action_space(env, p)) for p in (:Speaker, :Listener))

    create_actor(player) = Chain(
        Dense(length(state(env, player)), 64, relu; init = init),
        Dense(64, 64, relu; init = init),
        Dense(64, length(action_space(env, player)); init = init)
        )
    create_critic(critic_dim) = Chain(
        Dense(critic_dim, 64, relu; init = init),
        Dense(64, 64, relu; init = init),
        Dense(64, 1; init = init),
        )
    create_policy(player) = DDPGPolicy(
            behavior_actor = NeuralNetworkApproximator(
                model = create_actor(player),
                optimizer = Flux.Optimise.Optimiser(ClipNorm(0.5), ADAM(1e-2)),
            ),
            behavior_critic = NeuralNetworkApproximator(
                model = create_critic(critic_dim),
                optimizer = Flux.Optimise.Optimiser(ClipNorm(0.5), ADAM(1e-2)),
            ),
            target_actor = NeuralNetworkApproximator(
                model = create_actor(player),
            ),
            target_critic = NeuralNetworkApproximator(
                model = create_critic(critic_dim),
            ),
            γ = 0.95f0,
            ρ = 0.99f0,
            na = length(action_space(env, player)),
            start_steps = 0,
            start_policy = nothing,
            update_after = 512 * env.max_steps, # batch_size * env.max_steps
            act_limit = 1.0,
            act_noise = 0.,
        )
    create_trajectory(player) = CircularArraySARTTrajectory(
            capacity = 1_000_000, # replay buffer capacity
            state = Vector{Float64} => (length(state(env, player)), ),
            action = Vector{Float64} => (length(action_space(env, player)), ),
        )

    agents = MADDPGManager(
        Dict(
            player => Agent(
                policy = NamedPolicy(player, create_policy(player)),
                trajectory = create_trajectory(player),
            ) for player in (:Speaker, :Listener)
        ),
        SARTS, # trace's type
        512, # batch_size
        100, # update_freq
        0, # initial update_step
        rng
    )

    stop_condition = StopAfterEpisode(8_000, is_show_progress=!haskey(ENV, "CI"))
    hook = MeanRewardHook(0, 800, 100, [], [])
    Experiment(agents, env, stop_condition, hook, "# play SpeakerListener with MADDPG")
end

#+ tangle=false
using Plots
ex = E`JuliaRL_MADDPG_SpeakerListener`
run(ex)
plot(ex.hook.episodes, ex.hook.mean_rewards, xlabel="episode", ylabel="mean episode reward")

savefig("assets/JuliaRL_MADDPG_SpeakerListenerEnv.png") #hide

# ![](assets/JuliaRL_MADDPG_SpeakerListenerEnv.png)


text:
"""
    NexText(nexSerial, name; pid=pid, cid=cid)

A Nextion Text UI control.

A `NexText` is a Nextion UI control to display text on a page.

On Nextion side, it have the following properties:
- bco
- pco
- font
- xcen
- ycen
- pw
- txt
- isbr
"""
struct NexText <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    stringvalued::IStringValued
    fontstyleable::IFontStyleable
    colourable::IColourable
    touchable::ITouchable

    function NexText(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), IStringValued(nid), IFontStyleable(nid), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexText, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexText, property::Symbol)
    _getcommonproperty(obj, property)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

# Created 3/31/2015
# Based on visualization code originally written in Matlab provided by
# Mykel Kochenderfer, mykel@stanford.edu

module Label270

export get_textual_label, get_code

#code is an int with 4 digits corresponding to cc,vc,ua,da
function get_textual_label(code::Int,prev_code::Int,crossing::Bool,full::Bool=true)

  if code == 4010
    s = full ? "Climb" : "CL"
  elseif code == 5001
    s = full ? "Descend" : "DS"
  elseif code == 4110
    s = full ? "Climb, Crossing" : "CLX"
  elseif code == 5101
    s = full ? "Descend, Crossing" : "DSX"
  elseif code == 5002
    pc = prev_code
    if (pc==0 || pc==1000 || #if its a new encounter
        pc==6002 || pc==6003 || pc==6004 || pc==6005 || #if following a preventative limit climb
        pc==6020 || pc==6030 || pc==6040 || pc==6050) #if following a preventative limit descent
      s = full ? "Level-off (DNC)" : "LV_DNC"
    else
      s = full ? "Level-off (Weaken)" : "LV_WKN"
    end
  elseif code == 4020
    pc = prev_code
    if (pc==0 || pc==1000 || #if its a new encounter
        pc==6002 || pc==6003 || pc==6004 || pc==6005 || #if following a preventative limit climb
        pc==6020 || pc==6030 || pc==6040 || pc==6050) #if following a preventative limit descent
      s = full ? "Level-off (DND)" : "LV_DND"
    else
      s = full ? "Level-off (Weaken)" : "LV_WKN"
    end
  elseif code == 4210
    s = full ? "Reverse Climb" : "RV_CL"
  elseif code == 5201
    s = full ? "Reverse Descend" : "RV_DS"
  elseif code == 4310
    s =  full ? "Increase Climb" : "CL+"
  elseif code == 5301
    s = full ? "Increase Descend" : "DS+"
  elseif code == 4410
    if crossing
      s = full ? "Maintain Vertical Speed, Crossing" : "MTVSX"
    else
      s = full ? "Maintain Vertical Speed" : "MTVS"
    end
  elseif code == 5401
    if crossing
      s = full ? "Maintain Vertical Speed, Crossing" : "MTVSX"
    else
      s = full ? "Maintain Vertical Speed" : "MTVS"
    end
  elseif code == 6002
    s = full ? "Monitor Vertical Speed (DNC)" : "MOVS_DNC"
  elseif code == 6003
    s = full ? "Monitor Vertical Speed (500)" : "MOVS500"
  elseif code == 6004
    s = full ? "Monitor Vertical Speed (1000)" : "MOVS1000"
  elseif code == 6005
    s = full ? "Monitor Vertical Speed (2000)" : "MOVS2000"
  elseif code == 6020
    s = full ? "Monitor Vertical Speed (DND)" : "MOVS_DND"
  elseif code == 6030
    s = full ? "Monitor Vertical Speed (-500)" : "MOVS-500"
  elseif code == 6040
    s = full ? "Monitor Vertical Speed (-1000)" : "MOVS-1000"
  elseif code == 6050
    s = full ? "Monitor Vertical Speed (-2000)" : "MOVS-2000"
  elseif code == 6422
    s = full ? "Maintain Vertical Speed, MTE ()" : "MTVS_MTE"
  elseif code == 5022
    s = full ? "Level-off, MTE" : "LV_MTE"
  elseif code == 4022
    s = full ? "Multi-Aircraft Encounter (Descending)" : "MTE_DS"
  elseif code == 6022
    s = full ? "Multi-aircraft Encounter" : "MTE"
  elseif code == 1000 || code == 0000
    s = full ? "Clear of Conflict" : "COC"
  else
    s = full ? "N/A ($code)" : "?$code"
  end

  return s
end

function get_code(cc::Int,vc::Int,ua::Int,da::Int)
  dcc = digits(cc)
  dvc = digits(vc)
  dua = digits(ua)
  dda = digits(da)

  @assert length(dcc) == length(dvc) == length(dua) == length(dda) == 1

  s = string(dcc[1],dvc[1],dua[1],dda[1])
  return code = parse(Int64, s)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
"""
This module provides:
- An extension of the `AbstractArray` interface in order to properly deal with mutable caches.
- A mechanism to generate lazy arrays resulting from operations between arrays.
- A collection of concrete implementations of `AbstractArray`.

The exported names in this module are:

$(EXPORTS)
"""
module Arrays

using InplaceArrays.Helpers
using InplaceArrays.Inference

using DocStringExtensions
using Test
using FillArrays
using Base: @propagate_inbounds

export array_cache
export getindex!
export getitems!
export testitem
export uses_hash
export test_array
export testitems
export array_caches

export CachedArray
export CachedMatrix
export CachedVector
export setsize!

export CompressedArray
export LocalToGlobalArray
export Table

export kernel_cache
export kernel_caches
export apply_kernels!
export apply_kernel!
export apply_kernel
export test_kernel
export bcast
export elem
export contract
export kernel_return_type
export kernel_return_types
export Kernel

export apply
export apply_all

import Base: size
import Base: getindex, setindex!
import Base: similar
import Base: IndexStyle

include("Interface.jl")

include("CachedArrays.jl")

include("Kernels.jl")

include("Apply.jl")

include("CompressedArrays.jl")

include("Tables.jl")

include("LocalToGlobalArrays.jl")

end # module


text:
module ILUZero
using LinearAlgebra, SparseArrays

import LinearAlgebra.ldiv!, LinearAlgebra.\, SparseArrays.nnz

export ILU0Precon, \, forward_substitution!, backward_substitution!, nnz, ldiv!, ilu0, ilu0!

# ILU0 type definition
struct ILU0Precon{T <: Real,N <: Integer} <: Factorization{T}
    m::N
    n::N
    l_colptr::Vector{N}
    l_rowval::Vector{N}
    l_nzval::Vector{T}
    u_colptr::Vector{N}
    u_rowval::Vector{N}
    u_nzval::Vector{T}
    l_map::Vector{N}
    u_map::Vector{N}
    wrk::Vector{T}
end

# Allocates ILU0Precon type
function ILU0Precon(A::SparseMatrixCSC{T,N}) where {T <: Real,N <: Integer}
    m, n = size(A)

    # Determine number of elements in lower/upper
    lnz = 0
    unz = 0
    @inbounds for i = 1:n
        for j = A.colptr[i]:A.colptr[i + 1] - 1
            if A.rowval[j] > i
                lnz += 1
            else
                unz += 1
            end
        end
    end

    # Preallocate variables
    l_colptr = zeros(N, n + 1)
    u_colptr = zeros(N, n + 1)
    l_nzval = zeros(T, lnz)
    u_nzval = zeros(T, unz)
    l_rowval = zeros(Int64, lnz)
    u_rowval = zeros(Int64, unz)
    l_map = Vector{N}(undef, lnz)
    u_map = Vector{N}(undef, unz)
    wrk = zeros(T, n)
    l_colptr[1] = 1
    u_colptr[1] = 1

    # Map elements of A to lower and upper triangles, fill out colptr, and fill out rowval
    lit = 1
    uit = 1
    @inbounds for i = 1:n
        l_colptr[i + 1] = l_colptr[i]
        u_colptr[i + 1] = u_colptr[i]
        for j = A.colptr[i]:A.colptr[i + 1] - 1
            if A.rowval[j] > i
                l_colptr[i + 1] += 1
                l_rowval[lit] = A.rowval[j]
                l_map[lit] = j
                lit += 1
            else
                u_colptr[i + 1] += 1
                u_rowval[uit] = A.rowval[j]
                u_map[uit] = j
                uit += 1
            end
        end
    end

    return ILU0Precon(m, n, l_colptr, l_rowval, l_nzval, u_colptr, u_rowval, u_nzval, l_map, u_map, wrk)
end

# Updates ILU0Precon type in-place based on matrix A
function ilu0!(LU::ILU0Precon{T,N}, A::SparseMatrixCSC{T,N}) where {T <: Real,N <: Integer}
    m = LU.m
    n = LU.n
    l_colptr = LU.l_colptr
    l_rowval = LU.l_rowval
    l_nzval = LU.l_nzval
    u_colptr = LU.u_colptr
    u_rowval = LU.u_rowval
    u_nzval = LU.u_nzval
    l_map = LU.l_map
    u_map = LU.u_map

    # Redundant data or better speed... speed is chosen, but this might be changed.
    # This shouldn't be inbounded either.
    for i = 1:length(l_map)
        l_nzval[i] = A.nzval[l_map[i]]
    end
    for i = 1:length(u_map)
        u_nzval[i] = A.nzval[u_map[i]]
    end

    @inbounds for i = 1:m - 1
        multiplier = u_nzval[u_colptr[i + 1] - 1]
        for j = l_colptr[i]:l_colptr[i + 1] - 1
            l_nzval[j] /= multiplier
        end
        for j = u_colptr[i + 1]:u_colptr[i + 2] - 2
            multiplier = u_nzval[j]
            qn = j + 1
            rn = l_colptr[i + 1]
            pn = l_colptr[u_rowval[j]]
            while pn < l_colptr[u_rowval[j] + 1] && l_rowval[pn] <= i + 1
                while qn < u_colptr[i + 2] && u_rowval[qn] < l_rowval[pn]
                    qn += 1
                end
                if qn < u_colptr[i + 2] && l_rowval[pn] == u_rowval[qn]
                    u_nzval[qn] -= multiplier * l_nzval[pn]
                end
                pn += 1
            end
            while pn < l_colptr[u_rowval[j] + 1]
                while rn < l_colptr[i + 2] && l_rowval[rn] < l_rowval[pn]
                    rn += 1
                end
                if rn < l_colptr[i + 2] && l_rowval[pn] == l_rowval[rn]
                    l_nzval[rn] -= multiplier * l_nzval[pn]
                end
                pn += 1
            end
        end
    end
    return
end

# Constructs ILU0Precon type based on matrix A
function ilu0(A::SparseMatrixCSC{T,N}) where {T <: Real,N <: Integer}
    LU = ILU0Precon(A)
    ilu0!(LU, A)
    return LU
end

# Solves L\b and stores the solution in y
function forward_substitution!(y::AbstractVector{T}, LU::ILU0Precon{T,N}, b::AbstractVector{T}) where {T <: Real,N <: Integer}
    n = LU.n
    l_colptr = LU.l_colptr
    l_rowval = LU.l_rowval
    l_nzval  = LU.l_nzval

    y .= zero(T)

    @inbounds for i = 1:n
        y[i] += b[i]
        for j = l_colptr[i]:l_colptr[i + 1] - 1
            y[l_rowval[j]] -= l_nzval[j] * y[i]
        end
    end
    return y
end

# Solves U\y and stores the solution in x
function backward_substitution!(x::AbstractVector{T}, LU::ILU0Precon{T,N}, y::AbstractVector{T}) where {T <: Real,N <: Integer}
    n = LU.n
    u_colptr = LU.u_colptr
    u_rowval = LU.u_rowval
    u_nzval = LU.u_nzval

    wrk = LU.wrk
    wrk .= y

    @inbounds for i = n:-1:1
        x[i] = wrk[i] / u_nzval[u_colptr[i + 1] - 1]
        for j = u_colptr[i]:u_colptr[i + 1] - 2
            wrk[u_rowval[j]] -= u_nzval[j] * x[i]
        end
    end
    return x
end

# Solves LU\b overwriting x
function ldiv!(x::AbstractVector{T}, LU::ILU0Precon{T,N}, b::AbstractVector{T}) where {T <: Real,N <: Integer}
    (length(b) == LU.n) || throw(DimensionMismatch())
    n = LU.n
    l_colptr = LU.l_colptr
    l_rowval = LU.l_rowval
    l_nzval = LU.l_nzval
    u_colptr = LU.u_colptr
    u_rowval = LU.u_rowval
    u_nzval = LU.u_nzval
    wrk = LU.wrk

    wrk .= zero(T)

    @inbounds for i = 1:n
        wrk[i] += b[i]
        for j = l_colptr[i]:l_colptr[i + 1] - 1
            wrk[l_rowval[j]] -= l_nzval[j] * wrk[i]
        end
    end
    @inbounds for i = n:-1:1
        x[i] = wrk[i] / u_nzval[u_colptr[i + 1] - 1]
        for j = u_colptr[i]:u_colptr[i + 1] - 2
            wrk[u_rowval[j]] -= u_nzval[j] * x[i]
        end
    end
    return x
end

# Solves LU\b overwriting b
function ldiv!(LU::ILU0Precon{T,N}, b::AbstractVector{T}) where {T <: Real,N <: Integer}
    ldiv!(b, LU, b)
end

# Returns LU\b
function \(LU::ILU0Precon{T,N}, b::Vector{T}) where {T <: Real,N <: Integer}
    length(b) == LU.n || throw(DimensionMismatch())
    x = zeros(T, length(b))
    ldiv!(x, LU, b)
end

# Returns the number of nonzero
function nnz(LU::ILU0Precon{T,N}) where {T <: Real,N <: Integer}
    return length(LU.l_nzval) + length(LU.u_nzval)
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Cylinder(radius, bottom, top)

A solid circular cylinder embedded in R³ with given `radius`,
delimited by `bottom` and `top` planes.

    Cylinder(radius, segment)

Alternatively, construct a right circular cylinder with given `radius`
and `segment` between origin of `bottom` and `top` planes.

    Cylinder(radius)

Finally, construct a right vertical circular cylinder with given `radius`.

See https://en.wikipedia.org/wiki/Cylinder.
"""
struct Cylinder{T} <: Primitive{3,T}
  radius::T
  bot::Plane{T}
  top::Plane{T}
end

function Cylinder(radius::T, segment::Segment{3,T}) where {T}
  a, b = extrema(segment)
  v    = b - a
  bot  = Plane(a, v)
  top  = Plane(b, v)
  Cylinder(radius, bot, top)
end

function Cylinder(radius::T) where {T}
  _0 = (T(0), T(0), T(0))
  _1 = (T(0), T(0), T(1))
  segment = Segment(_0, _1)
  Cylinder(radius, segment)
end

paramdim(::Type{<:Cylinder}) = 3

isconvex(::Type{<:Cylinder}) = true

radius(c::Cylinder) = c.radius

axis(c::Cylinder) = Line(origin(c.bot), origin(c.top))

planes(c::Cylinder) = (c.bot, c.top)

isright(c::Cylinder) = isright(boundary(c))

boundary(c::Cylinder) = CylinderSurface(c.radius, c.bot, c.top)


text:
## names.jl retrieve and set dimension names
## (c) 2013 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

import Base.names

## `names` is somewhat loaded, as some of the semantics were transferred tp `fieldnames`
names(dict::Associative) = collect(keys(dict))
allnames(a::NamedArray) = [names(dict) for dict in a.dicts]
allnames(a::AbstractArray) = [[string(i) for i in 1:size(a, d)] for d in 1:ndims(a)] ## default dimension names
names(a::NamedArray, d::Integer) = names(a.dicts[d])
dimnames(a::NamedArray) = [dn for dn in a.dimnames]
dimnames(a::NamedArray, d::Integer) = a.dimnames[d]
dimnames(a::AbstractArray) = [letter(i) for i in 1:ndims(a)]
dimnames(a::AbstractArray, d::Integer) = letter(i)

## string versions of the above
strnames(dict::Associative) = map(string, names(dict))
strnames(a::NamedArray) = [strnames(d) for d in a.dicts]
strnames(a::NamedArray, d::Integer) = strnames(a.dicts[d])
strdimnames(a::NamedArray) = [string(dn) for dn in a.dimnames]
strdimnames(a::NamedArray, d::Integer) = string(a.dimnames[d])

## seting names, dimnames
function setnames!(a::NamedArray, v::Vector, d::Integer)
    size(a.array,d) == length(v) || error("inconsistent vector length")
    eltype(keys(a.dicts[d])) == eltype(v) || error("inconsistent name type")
    ## a.dicts is a tuple, so we need to replace it as a whole...
    vdicts = OrderedDict[]
    for i = 1:length(a.dicts)
        if i==d
            push!(vdicts, OrderedDict(zip(v, 1:length(v))))
        else
            push!(vdicts, a.dicts[i])
        end
    end
    a.dicts = tuple(vdicts...)
end

function setnames!(a::NamedArray, v, d::Integer, i::Integer)
    @assert 1 <= d <= ndims(a)
    @assert 1 <= i <= size(a, d)
    filter!((k,v) -> v!=i, a.dicts[d]) # remove old name
    a.dicts[d][v] = i
end

function setdimnames!{T,N}(a::NamedArray{T,N}, dn::NTuple{N})
    a.dimnames = dn
end

setdimnames!(a::NamedArray, dn::Vector) = setdimnames!(a, tuple(dn...))

function setdimnames!{T,N}(n::NamedArray{T,N}, v, d::Integer)
    1 <= d <= N || throw(BoundsError(size(n), d))
    vdimnames = Array(Any, N)
    for i=1:N
        if i==d
            vdimnames[i] = v
        else
            vdimnames[i] = n.dimnames[i]
        end
    end
    n.dimnames = tuple(vdimnames...)
end

defaultnames(a::AbstractArray, dim::Integer) = [string(i) for i=1:size(a, dim)]
defaultnames(a::AbstractArray) = [defaultnames(d) for d in size(a)]


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
@doc doc"""
    TruncatedDensity

Constructor:

    TruncatedDensity(D<:AbstractDensity, B<:VarVolumeBounds)

*BAT-internal, not part of stable public API.*

A density truncated to given bounds, forced to effectively zero outside of
those bounds. In contrast to a truncated distribution, the density is
*not* renormalized.
"""
struct TruncatedDensity{D<:AbstractDensity,B<:AbstractVarBounds,T<:Real} <: AbstractDensity
    density::D
    bounds::B
    logscalecorr::T
end

Base.parent(density::TruncatedDensity) = density.density

var_bounds(density::TruncatedDensity) = density.bounds

ValueShapes.varshape(density::TruncatedDensity) = varshape(parent(density))

function eval_logval(
    density::TruncatedDensity,
    v::Any,
    T::Type{<:Real} = density_logval_type(v);
    use_bounds::Bool = true,
    strict::Bool = false
)
    v_shaped = reshape_variate(varshape(density), v)
    if use_bounds && !variate_is_inbounds(density, v_shaped, strict)
        return log_zero_density(T)
    end

    parent_logval = eval_logval(
        parent(density), v_shaped,
        use_bounds = false, strict = false
    )

    parent_logval + density.logscalecorr
end

function eval_logval_unchecked(density::TruncatedDensity, v::Any)
    eval_logval(
        density, v,
        use_bounds = false, strict = false
    )
end

@doc doc"""
    BAT.truncate_density(density::AbstractDensity, bounds::AbstractVarBounds)::AbstractDensity

*BAT-internal, not part of stable public API.*

Truncate `density` to `bounds`, the resulting density will be effectively
zero outside of those bounds. In contrast `Distributions.truncated`, `truncate_density`
does *not* renormalize the density.

Currently implemented for `BAT.HyperRectBounds` only.
"""
function truncate_density end

function truncate_density(density::AbstractPosteriorDensity, bounds::AbstractVarBounds)
    old_prior = getprior(density)
    old_bounds = var_bounds(old_prior)
    new_bounds = ismissing(old_bounds) ? bounds : var_bounds(old_prior) ∩ bounds
    newprior = truncate_density(getprior(density), new_bounds)
    PosteriorDensity(getlikelihood(density), newprior)
end

function truncate_density(density::ConstDensity, bounds::AbstractVarBounds)
    old_bounds = var_bounds(density)
    new_bounds = ismissing(old_bounds) ? bounds : var_bounds(density) ∩ bounds
    ConstDensity(density.value, new_bounds)
end

function truncate_density(density::DistributionDensity, bounds::HyperRectBounds)
    old_bounds = var_bounds(density)
    new_bounds = ismissing(old_bounds) ? bounds : var_bounds(density) ∩ bounds
    interval_bounds = Interval.(new_bounds.vol.lo, new_bounds.vol.hi)
    dist = density.dist
    r = truncate_dist_hard(dist, interval_bounds)
    TruncatedDensity(DistributionDensity(r.dist), new_bounds, r.logscalecorr)
end

function truncate_density(density::TruncatedDensity{<:DistributionDensity}, bounds::AbstractVarBounds)
    # Tricky: How to decide whether to propagate renormalization of inner densities/distributions ?
    @assert false # Not implemented yet
end

# ToDo: Reject out-of-bounds samples, in case density contains distributions that can't be truncated (e.g. multivariate)
Distributions.sampler(density::TruncatedDensity{<:DistributionDensity}) = Distributions.sampler(parent(density))

function Statistics.cov(density::TruncatedDensity{<:DistributionDensity})
    cov(nestedview(rand(bat_determ_rng(), sampler(density), 10^5)))
end

@doc doc"""
    BAT.truncate_dist_hard(dist::Distribution{Univariate}, bounds::Interval)::Distribution{Univariate}
    BAT.truncate_dist_hard(dist::Distribution{Multivariate}, bounds::AbstractArray{<:Interval})::Distribution{Multivariate}

*BAT-internal, not part of stable public API.*

Generalized variant of `Distributions.truncated` - also handles multivariate
distributions and operates on a best-effort basis: If distributions cannot
be truncated, may return the original distribution.

Returns a `NamedTuple`

```julia
    (dist = trunc_dist, logscalecorr = logscalecorr)
```

with the truncated distribution and the log-PDF amplitude difference to
the original (see [`BAT.trunc_logpdf_ratio`](@ref)).

Mainly used to implement [`BAT.truncate_density`](@ref).
"""
function truncate_dist_hard end

function truncate_dist_hard(dist::Distribution{Univariate}, bounds::AbstractArray{<:Interval})
    @argcheck length(eachindex(bounds)) == 1
    truncate_dist_hard(dist, first(bounds))
end

function truncate_dist_hard(dist::Distribution{Univariate}, bounds::Interval)
    min_lo = minimum(dist)
    max_hi = maximum(dist)
    lo = clamp(minimum(bounds), min_lo, max_hi)
    hi = clamp(max(lo, maximum(bounds)), min_lo, max_hi)

    trunc_dist = truncated(dist, lo, hi)
    logscalecorr = trunc_logpdf_ratio(dist, trunc_dist)
    return (dist = trunc_dist, logscalecorr = logscalecorr)
end

function truncate_dist_hard(dist::Distributions.Truncated, bounds::Interval)
    # Note: Assumes that dist is result of trunctate, not of truncate_dist_hard:
    untrunc_dist = dist.untruncated

    min_lo = minimum(untrunc_dist)
    max_hi = maximum(untrunc_dist)
    lo = clamp(max(minimum(bounds), dist.lower), min_lo, max_hi)
    hi = clamp(max(lo, min(maximum(bounds), dist.upper)), min_lo, max_hi)
    trunc_dist = truncated(untrunc_dist, lo, hi)
    logscalecorr = trunc_logpdf_ratio(dist, trunc_dist)
    return (dist = trunc_dist, logscalecorr = logscalecorr)
end

function truncate_dist_hard(d::Product, bounds::AbstractArray{<:Interval})
    @argcheck length(eachindex(bounds)) == length(d)
    r = truncate_dist_hard.(d.v, bounds)
    trunc_dists = map(x -> x.dist, r)
    logscalecorr = sum(x.logscalecorr for x in r)

    return (dist = Product(trunc_dists), logscalecorr = logscalecorr)
end

function truncate_dist_hard(dist::ConstValueDist, bounds::AbstractVector{<:Interval})
    @argcheck length(eachindex(bounds)) == 0
    (dist = dist, logscalecorr = 0)
end

function truncate_dist_hard(dist::NamedTupleDist{names}, bounds::AbstractArray{<:Interval}) where names
    @argcheck length(eachindex(bounds)) == totalndof(varshape(dist))
    distributions = values(dist)
    accessors = values(varshape(dist))

    r = map((dist, acc) -> truncate_dist_hard(dist, view(bounds, ValueShapes.view_idxs(eachindex(bounds), acc))), distributions, accessors)
    trunc_dist = NamedTupleDist(NamedTuple{names}(map(x -> x.dist, r)))
    logscalecorr = sum(map(x -> x.logscalecorr, r))
    (dist = trunc_dist, logscalecorr = logscalecorr)
end

@doc doc"""
    BAT.trunc_logpdf_ratio(orig_dist::Distribution{TP}, trunc_dist::Distribution{TP})::AbstractFloat

*BAT-internal, not part of stable public API.*

Computes the log-ratio between the amplitude of the PDF of a truncated
distribution and the original (untruncted) distribution, within the support
of the truncated one.

The PDF of both distributions must have the same shape within the support of
`trunc_dist` and may only differ in amplitude.

Mainly used to implement [`BAT.truncate_density`](@ref), in conjunction with
[`BAT.truncate_dist_hard`](@ref). The result contributes to the `logscalecorr`
factor of a [`TruncatedDensity`] that uses truncated distributions internally,
to ensure the density does not get renormalized.
"""
function trunc_logpdf_ratio end

function trunc_logpdf_ratio(orig_dist::Distribution, trunc_dist::Distribution)
    _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist)
end

function _trunc_logpdf_ratio_fallback(orig_dist::Distribution, trunc_dist::Distribution)
    x = mean(trunc_dist)
    +logpdf(orig_dist, x) - logpdf(trunc_dist, x)
end

function trunc_logpdf_ratio(orig_dist::Distributions.Truncated, trunc_dist::Distributions.Truncated)
    T = promote_type(typeof(trunc_dist.logtp), typeof(orig_dist.logtp))
    if orig_dist.untruncated == trunc_dist.untruncated
        convert(T, trunc_dist.logtp - orig_dist.logtp)
    else
        convert(T, _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist))
    end
end

function trunc_logpdf_ratio(orig_dist::Distribution, trunc_dist::Distributions.Truncated)
    T = typeof(trunc_dist.logtp)
    if orig_dist == trunc_dist.untruncated
        trunc_dist.logtp
    else
        convert(T, _trunc_logpdf_ratio_fallback(orig_dist, trunc_dist))
    end
end


text:
# Autogenerated wrapper script for MAGEMin_jll for x86_64-w64-mingw32-libgfortran4
export MAGEMin, libMAGEMin

using MicrosoftMPI_jll
using NLopt_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MAGEMin")
JLLWrappers.@declare_library_product(libMAGEMin, "libMAGEMin.dll")
JLLWrappers.@declare_executable_product(MAGEMin)
function __init__()
    JLLWrappers.@generate_init_header(MicrosoftMPI_jll, NLopt_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libMAGEMin,
        "bin\\libMAGEMin.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        MAGEMin,
        "bin\\MAGEMin.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
pogm_restart.jl
2017-03-31, Donghwan Kim and Jeff Fessler, University of Michigan
=#

export pogm_restart

using LinearAlgebra: norm, opnorm
using Random: seed!
using Test: @test, @test_throws

function gr_restart(Fgrad, ynew_yold, restart_cutoff)
	return sum(Float64, real(-Fgrad .* ynew_yold)) <=
		restart_cutoff * norm(Fgrad) * norm(ynew_yold)
end

"""
`x, out = pogm_restart(x0, Fcost, f_grad, f_L ;
	f_mu=0, mom=:pogm, restart=:gr, restart_cutoff=0.,
	bsig=1, niter=10, g_prox=(z,c)->z, fun=...)`

Iterative proximal algorithms (PGM=ISTA, FPGM=FISTA, POGM) with restart.

# in
- `x0` initial guess
- `Fcost` function for computing the cost function value ``F(x)``
  - (needed only if `restart === :fr`)
- `f_grad` function for computing the gradient of ``f(x)``
- `f_L` Lipschitz constant of the gradient of ``f(x)``

# option
- `f_mu` strong convexity parameter of ``f(x)``; default 0.
  - if `f_mu > 0`, ``(\\alpha, \\beta_k, \\gamma_k)`` is chosen by Table 1 in [KF18]
- `g_prox` function `g_prox(z,c)` for the proximal operator for ``g(x)``
  - `g_prox(z,c)` computes ``argmin_x 1/2 \\|z-x\\|^2 + c \\, g(x)``
- `mom`	momentum option
  - `:pogm` POGM (fastest); default!
  - `:fpgm` (FISTA), ``\\gamma_k = 0``
  - `:pgm` PGM (ISTA), ``\\beta_k = \\gamma_k = 0``
- `restart` restart option
  - `:gr` gradient restart; default!
  - `:fr` function restart
  - `:none` no restart
- `restart_cutoff` for `:gr` restart if cos(angle) < this; default 0.
- `bsig` gradient "gamma" decrease option (value within [0 1]); default 1
  - see ``\\bar{\\sigma}`` in [KF18]
- `niter` number of iterations; default 10
- `fun` function`(iter, xk, yk, is_restart)` user-defined function evaluated each `iter` with secondary `xk`, primary `yk`, and boolean `is_restart` indicating whether this iteration was a restart

# out
- `x`	final iterate
  - for PGM (ISTA): ``x_N = y_N``
  - for FPGM (FISTA): primary iterate ``y_N``
  - for POGM: secondary iterate ``x_N``, see [KF18]
- `out [fun(0, x0, x0, false), fun(1, x1, y1, is_restart), ...]` array of length `[niter+1]`

Optimization Problem: Nonsmooth Composite Convex Minimization
* ``argmin_x F(x),  F(x) := f(x) + g(x))``
  - ``f(x)`` smooth convex function
  - ``g(x)`` convex function, possibly nonsmooth and "proximal-friendly" [CP11]

# Optimization Algorithms:

Accelerated First-order Algorithms when ``g(x) = 0`` [KF18]
iterate as below for given coefficients ``(\\alpha, \\beta_k, \\gamma_k)``
* For k = 0,1,...
  - ``y_{k+1} = x_k - \\alpha  f'(x_k)`` : gradient update
  - ``x_{k+1} = y_{k+1} + \\beta_k  (y_{k+1} - y_k) + \\gamma_k  (y_{k+1} - x_k)`` : momentum update

Proximal versions of the above for ``g(x) \\neq 0`` are in the below references,
and use the proximal operater
``prox_g(z) = argmin_x {1/2\\|z-x\\|^2 + g(x)}``.

- Proximal Gradient method (PGM or ISTA) - ``\\beta_k = \\gamma_k = 0``. [BT09]
- Fast Proximal Gradient Method (FPGM or FISTA) - ``\\gamma_k = 0``. [BT09]
- Proximal Optimized Gradient Method (POGM) - [THG15]
- FPGM(FISTA) with Restart - [OC15]
- POGM with Restart - [KF18]

# references

- [CP11] P. L. Combettes, J. C. Pesquet,
 "Proximal splitting methods in signal processing,"
 Fixed-Point Algorithms for Inverse Problems in Science and Engineering,
 Springer, Optimization and Its Applications, 2011.
- [KF18] D. Kim, J.A. Fessler,
 "Adaptive restart of the optimized gradient method for convex optimization," 2018
 Arxiv:1703.04641,
 [http://doi.org/10.1007/s10957-018-1287-4]
- [BT09] A. Beck, M. Teboulle:
 "A fast iterative shrinkage-thresholding algorithm for linear inverse problems,"
 SIAM J. Imaging Sci., 2009.
- [THG15] A.B. Taylor, J.M. Hendrickx, F. Glineur,
 "Exact worst-case performance of first-order algorithms
 for composite convex optimization," Arxiv:1512.07516, 2015,
 SIAM J. Opt. 2017
 [http://doi.org/10.1137/16m108104x]

Copyright 2017-3-31, Donghwan Kim and Jeff Fessler, University of Michigan
2018-08-13 Julia 0.7.0
2019-02-24 interface redesign
"""
function pogm_restart(x0, Fcost::Function, f_grad::Function, f_L::Real ;
		f_mu::Real = 0.,
		mom::Symbol = :pogm, # :ogm :gm
		restart::Symbol = :gr, # :fr :none
		restart_cutoff::Real = 0.,
		bsig::Real = 1,
		niter::Int = 10,
		g_prox::Function = (z, c::Real) -> z,
		fun::Function = (iter::Int, xk, yk, is_restart::Bool) -> undef,
	)

	!in(mom, (:pgm, :fpgm, :pogm)) && throw(ArgumentError("mom $mom"))
	!in(restart, (:none, :gr, :fr)) && throw(ArgumentError("restart $restart"))
	f_L < 0 && throw(ArgumentError("f_L=$f_L < 0"))
	f_mu < 0 && throw(ArgumentError("f_mu=$f_mu < 0"))
	bsig < 0 && throw(ArgumentError("bsig=$bsig < 0"))
	!((-1 < restart_cutoff) && (restart_cutoff < 1)) &&
		throw(ArgumentError("restart_cutoff=$restart_cutoff"))

	L = f_L
	mu = f_mu
	q = mu/L

	# initialize parameters
	told = 1
	sig = 1
	zetaold = 1 # dummy

	# initialize x
	xold = x0
	yold = x0
	uold = x0
	zold = x0
	Fcostold = Fcost(x0)
	Fgradold = zeros(size(x0)) # dummy

	# save initial
	out = Array{Any}(undef, niter+1)
	out[1] = fun(0, x0, x0, false)

	xnew = []
	ynew = []

# iterations
for iter=1:niter

	# proximal gradient (PGM) update
	if mom === :pgm && mu != 0
		alpha = 2. / (L+mu)
	else
		alpha = 1. / L
	end

	fgrad = f_grad(xold)

	is_restart = false

	if mom === :pgm || mom === :fpgm
		ynew = g_prox(xold - alpha * fgrad, alpha) # standard PG update
		Fgrad = -(1. / alpha) * (ynew - xold) # standard composite gradient mapping
		Fcostnew = Fcost(ynew)

		# restart condition
		if restart != :none
			# function/gradient restart
			if ((restart === :fr && Fcostnew > Fcostold)
			|| (restart === :gr && gr_restart(Fgrad, ynew-yold, restart_cutoff)))
				told = 1
				is_restart = true
			end
			Fcostold = Fcostnew
		end

	elseif mom === :pogm # POGM
		# gradient update for POGM [see KF18]
		unew = xold - alpha * fgrad
		# restart + "gamma" decrease conditions checked later for POGM,
		# unlike PGM, FPGM above

#	else
#		throw("bad mom $mom")
	end

	# momentum coefficient "beta"
	if mom === :fpgm && mu != 0 # known μ > 0
		beta = (1 - sqrt(q)) / (1 + sqrt(q))
	elseif mom === :pogm && mu != 0
		beta = (2 + q - sqrt(q^2+8*q))^2 / 4. / (1-q)
	# for "mu" = 0 or for unknown "mu"
	elseif mom != :pgm
		if mom === :pogm && iter == niter # && restart == 0
			tnew = 0.5 * (1 + sqrt(1 + 8 * told^2))
		else
			tnew = 0.5 * (1 + sqrt(1 + 4 * told^2))
		end

		beta = (told - 1) / tnew
	end

	# momentum update
	if mom === :pgm
		xnew = ynew
	elseif mom === :fpgm
		xnew = ynew + beta * (ynew - yold)
	elseif mom === :pogm # see [KF18]
		# momentum coefficient "gamma"
		if mu != 0
			gamma = (2 + q - sqrt(q^2+8*q)) / 2.
		else
			gamma = sig * told / tnew
		end

		znew = (unew + beta * (unew - uold) + gamma * (unew - xold)
				- beta * alpha / zetaold * (xold - zold))
		zetanew = alpha * (1 + beta + gamma)
		xnew = g_prox(znew, zetanew) # non-standard PG update for POGM

		# non-standard composite gradient mapping for POGM:
		Fgrad = fgrad - 1/zetanew * (xnew - znew)
		ynew = xold - alpha * Fgrad
		Fcostnew = Fcost(xnew)

		# restart + "gamma" decrease conditions for POGM
		if restart != :none
			# function/gradient restart
			if ((restart === :fr && Fcostnew > Fcostold)
			|| (restart === :gr && gr_restart(Fgrad, ynew-yold, restart_cutoff)))
				tnew = 1
				sig = 1
				is_restart = true

			# gradient "gamma" decrease
			elseif sum(Float64, real(Fgrad .* Fgradold)) < 0
				sig = bsig * sig
			end

			Fcostold = Fcostnew
			Fgradold = Fgrad
		end

		uold = unew
		zold = znew
		zetaold = zetanew
	end

	out[iter+1] = fun(iter, xnew, ynew, is_restart) # save

	xold = xnew
	yold = ynew

	if mom != :pgm && mu == 0
		told = tnew
	end
end # for iter

	return ((mom === :pogm) ? xnew : ynew), out
end # pogm_restart()

"""
`pogm_restart_test()`
self test
"""
function pogm_restart_test()
	# Tikhonov regularized LS problem
	seed!(0); M = 30; N = 6; A = randn(M,N); y = randn(M)
	a2 = opnorm(A)^2
	reg = 0.1 * a2
	xh = (A'A + reg*I) \ A'y
	Fcost = (x) -> 1/2 * norm(A * x - y)^2
	cost = (x) -> Fcost(x) + reg/2 * norm(x)^2
	fun = (iter, x, y, is_restart) ->
		(cost(x) - cost(xh), norm(x - xh) / norm(xh), time())
	f_grad = (x) -> A' * (A * x - y)
	g_prox = (z,c) -> z / (1 + reg * c) # proximal operator for 2-norm
	grad = (x) -> A' * (A * x - y) + reg * x
	x0 = A \ y

	x, out = pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=0, mom=:pogm, restart=:gr, restart_cutoff=0.,
		bsig=1, niter=100, g_prox=g_prox, fun=fun)
	@test isapprox(x, xh)
	x, out = pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=reg, mom=:pogm, restart=:gr, restart_cutoff=0., # with mu
		bsig=1, niter=100, g_prox=g_prox, fun=fun)
	@test isapprox(x, xh)
	x, _ = pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=0, mom=:fpgm, niter=100, g_prox=g_prox, fun=fun)
	pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=0, mom=:pgm, niter=2, g_prox=g_prox)
	pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=reg, mom=:fpgm, niter=2, g_prox=g_prox) # with mu
	pogm_restart(x0, Fcost, f_grad, a2;
		f_mu=reg, mom=:pgm, niter=2, g_prox=g_prox) # with mu
	@test_throws ArgumentError pogm_restart(x0, Fcost, f_grad, a2; mom=:bad)
	true
end

"""
`pogm_restart(:test)` self test
"""
function pogm_restart(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	@test pogm_restart_test()
	true
end


text:
using DSGE, JLD2, DifferentialEquations

# Get transition equations
m = KrusellSmithCT()
T, R, C, inverse_basis, basis = solvect(m)

# Make initial guess
n_vars = n_states(m)
n_expectation_errors = n_shocks_expectational(m)
n_shocks = n_shocks_exogenous(m)
x = zeros(Float64, n_vars)

# Set up SDE
f(u,p,t) = T*u
g(u,p,t) = R
dt = 1./90 # daily frequency when one period is one month
tspan = (0.0,4. * 50.) # 50 years of data
u0 = reshape(basis * x, size(inverse_basis, 2), 1)
W = WienerProcess(0.0, 0., 0.)
prob = SDEProblem(f, g, u0, tspan, W)
sol = solve(prob, EM(), dt = dt)

save("fine_data_red_basis.jld2", "states", sol.u)
states = sol.u
data = zeros(1, (length(states) - 1) / 90)

# created 2 years worth of daily data -> transform into quarterly data for discrete observations
for i = 1:length(data)
    full_state = inverse_basis * states[1+ (i-1)*90]
    data[i] = full_state[m.endogenous_states[:output][1]]
end
save("simulated_data.jld2", "data", data)


text:
function test_mathops(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Unary operations for ", T, "                          ") )" begin
        x = randn()

        @test to_julia( -T(x) ) ≈ -x
    end

    @testset "$( string("Binary operations for ", T, "                         ") )" begin
        x = rand()
        y = rand()

        @test to_julia( T(x) + T(y) ) ≈ x + y
        @test to_julia( T(x) * T(y) ) ≈ x * y
        @test to_julia( T(x) - T(y) ) ≈ x - y
        @test to_julia( T(x) / T(y) ) ≈ x / y
        @test to_julia( T(x) ^ T(y) ) ≈ x ^ y
        @test to_julia( T(x) \ T(y) ) ≈ x \ y
    end

    @testset "$( string("Comparisons for ", T, "                               ") )" begin
        x = randn()
        y = x + rand()

        @test Bool( to_julia(T(y) >= T(x)) )
        @test Bool( to_julia(T(y) >  T(x)) )
        @test Bool( to_julia(T(x) <= T(y)) )
        @test Bool( to_julia(T(x) <  T(y)) )
        @test Bool( to_julia(T(x) == T(x)) )
    end

    @testset "$( string("Cross product for ", T, "                             ") )" begin
        n₁ = rand(3)
        n₂ = rand(3)
        c₁ = T(n₁)
        c₂ = T(n₂)

        n₁₂ = n₁×n₂
        c₁₂ = c₁×c₂
        @test casadi.evalf(c₁₂).toarray()[:] == n₁₂

        c₁₂ = cross(c₁, c₂)
        @test casadi.evalf(c₁₂).toarray()[:] == n₁₂
    end
end