text:
@recipe f(x::AbstractVector, gp::AbstractGP) = gp(x)
@recipe f(gp::FiniteGP) = (gp.x, gp)
@recipe function f(x::AbstractVector, gp::FiniteGP)
    length(x) == length(gp.x) ||
        throw(DimensionMismatch("length of `x` and `gp.x` has to be equal"))
    scale::Float64 = pop!(plotattributes, :ribbon_scale, 1.0)
    scale >= 0.0 || error("`ribbon_scale` keyword argument must be non-negative")

    # compute marginals
    μ, σ2 = mean_and_var(gp)

    ribbon := scale .* sqrt.(σ2)
    fillalpha --> 0.3
    linewidth --> 2
    return x, μ
end

"""
    plot(x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)
    plot!([plot, ]x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)

Plot the predictive mean for the projection `f` of a Gaussian process and a ribbon of
`ribbon_scale` standard deviations around it versus `x`.

!!! note
    Make sure to load [Plots.jl](https://github.com/JuliaPlots/Plots.jl) before you use
    this function.

# Examples

Plot the mean and a ribbon of 3 standard deviations:

```julia
using Plots

gp = GP(SqExponentialKernel())
plot(gp(rand(5)); ribbon_scale=3)
```
"""
RecipesBase.plot(::AbstractVector, ::FiniteGP; kwargs...)
@doc (@doc RecipesBase.plot(::AbstractVector, ::FiniteGP)) RecipesBase.plot!(
    ::AbstractVector, ::FiniteGP; kwargs...
)
@doc (@doc RecipesBase.plot(::AbstractVector, ::FiniteGP)) RecipesBase.plot!(
    ::RecipesBase.AbstractPlot, ::AbstractVector, ::FiniteGP; kwargs...
)

"""
    plot(f::FiniteGP; kwargs...)
    plot!([plot, ]f::FiniteGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection `f` of a Gaussian
process versus `f.x`.
"""
RecipesBase.plot(::FiniteGP; kwargs...)
@doc (@doc RecipesBase.plot(::FiniteGP)) RecipesBase.plot!(::FiniteGP; kwargs...)
@doc (@doc RecipesBase.plot(::FiniteGP)) RecipesBase.plot!(
    ::RecipesBase.AbstractPlot, ::FiniteGP; kwargs...
)

"""
    plot(x::AbstractVector, gp::AbstractGP; kwargs...)
    plot!([plot, ]x::AbstractVector, gp::AbstractGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection `gp(x)` of the Gaussian
process `gp`.
"""
RecipesBase.plot(::AbstractVector, ::AbstractGP; kwargs...)
@doc (@doc RecipesBase.plot(::AbstractVector, ::AbstractGP)) RecipesBase.plot!(
    ::AbstractVector, ::AbstractGP; kwargs...
)
@doc (@doc RecipesBase.plot(::AbstractVector, ::AbstractGP)) RecipesBase.plot!(
    ::RecipesBase.AbstractPlot, ::AbstractVector, ::AbstractGP; kwargs...
)

"""
    sampleplot([x::AbstractVector=f.x, ]f::FiniteGP; samples=1, kwargs...)

Plot samples from the projection `f` of a Gaussian process versus `x`.

!!! note
    Make sure to load [Plots.jl](https://github.com/JuliaPlots/Plots.jl) before you use
    this function.

When plotting multiple samples, these are treated as a _single_ series (i.e.,
only a single entry will be added to the legend when providing a `label`).

# Example

```julia
using Plots

gp = GP(SqExponentialKernel())
sampleplot(gp(rand(5)); samples=10, linealpha=1.0)
```
The given example plots 10 samples from the projection of the GP `gp`.
The `linealpha` is modified from default of 0.35 to 1.

---
    sampleplot(x::AbstractVector, gp::AbstractGP; samples=1, kwargs...)

Plot samples from the finite projection `gp(x, 1e-9)` versus `x`.
"""
@userplot struct SamplePlot{X<:AbstractVector,F<:FiniteGP}
    x::X
    f::F
end

# default constructor (recipe forwards arguments as tuple)
SamplePlot((x, f)::Tuple{<:AbstractVector,<:FiniteGP}) = SamplePlot(x, f)

# `FiniteGP`s without explicit `x`
SamplePlot((f,)::Tuple{<:FiniteGP}) = SamplePlot((f.x, f))

# `AbstractGP`s with explicit `x`
# zero mean observation noise with variance 1e-9 to avoid numerical issues in
# Cholesky decomposition
SamplePlot((x, gp)::Tuple{<:AbstractVector,<:AbstractGP}) = SamplePlot((gp(x, 1e-9),))

@recipe function f(sp::SamplePlot)
    nsamples::Int = pop!(plotattributes, :samples, 1)
    samples = rand(sp.f, nsamples)

    flat_x = repeat(vcat(sp.x, NaN), nsamples)
    flat_f = vec(vcat(samples, fill(NaN, 1, nsamples)))

    # Set default attributes
    linealpha --> 0.35
    label --> ""

    return flat_x, flat_f
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using Mimi
include("region_parameters.jl")
include("two-region-model.jl")

run1 = run_my_model()

#Check results
run1[:emissions, :E_Global]

ext: jl


text:
function ecephyscache()
    ecephys_project_cache.EcephysProjectCache.from_warehouse(manifest=ecephysmanifest)
end

# These will take a while to download
function getsessiontable()
    @info "Please wait, this can take a few seconds"
    CSV.read(IOBuffer(ecephyscache().get_session_table().to_csv()), DataFrame);
end
export getsessiontable

function getprobes()
    CSV.read(IOBuffer(ecephyscache().get_probes().to_csv()), DataFrame);
end
export getprobes

function getchannels()
    CSV.read(IOBuffer(ecephyscache().get_channels().to_csv()), DataFrame);
end
export getchannels

function getunits(; filter_by_validity=true, amplitude_cutoff_maximum = 0.1, presence_ratio_minimum = 0.9, isi_violations_maximum = 0.5)
    str = ecephyscache().getunits(filter_by_validity=filter_by_validity,
                            amplitude_cutoff_maximum=amplitude_cutoff_maximum,
                            presence_ratio_minimum=presence_ratio_minimum,
                            isi_violations_maximum=isi_violations_maximum).to_csv()
    CSV.read(IOBuffer(str), DataFrame);
end
export getunits

function getunitanalysismetricsbysessiontype(session_type; filter_by_validity=true, amplitude_cutoff_maximum = 0.1, presence_ratio_minimum = 0.9, isi_violations_maximum = 0.5) # Yeah thats python
    str = ecephyscache().get_unit_analysis_metrics_by_session_type(session_type,
                            filter_by_validity=filter_by_validity,
                            amplitude_cutoff_maximum=amplitude_cutoff_maximum,
                            presence_ratio_minimum=presence_ratio_minimum,
                            isi_violations_maximum=isi_violations_maximum).to_csv()
    CSV.read(IOBuffer(str), DataFrame);
end
export getunitanalysismetricsbysessiontype

function getallunitmetrics() # This one is really slow
    metrics1 = get_unit_analysis_metrics_by_session_type("brain_observatory_1.1",
                            amplitude_cutoff_maximum = Inf,
                            presence_ratio_minimum = -Inf,
                            isi_violations_maximum = Inf)

    metrics2 = get_unit_analysis_metrics_by_session_type("functional_connectivity",
                            amplitude_cutoff_maximum = Inf,
                            presence_ratio_minimum = -Inf,
                            isi_violations_maximum = Inf)

    vcat(analysis_metrics1, analysis_metrics2)
end
export getallunitmetrics

function getsessiondata(session_id::Int; filter_by_validity=true, amplitude_cutoff_maximum = 0.1, presence_ratio_minimum = 0.9, isi_violations_maximum = 0.5)
    ecephyscache().get_session_data(session_id; filter_by_validity=filter_by_validity,
                            amplitude_cutoff_maximum=amplitude_cutoff_maximum,
                            presence_ratio_minimum=presence_ratio_minimum,
                            isi_violations_maximum=isi_violations_maximum)
end
export getsessiondata

abstract type AbstractSession end

struct Session <: AbstractSession
    pyObject
end
export Session
Session(session_id::Int; kwargs...) = Session(getsessiondata(session_id; kwargs...))
getid(S::AbstractSession) = S.pyObject.ecephys_session_id
getprobes(S::AbstractSession) = CSV.read(IOBuffer(S.pyObject.probes.to_csv()), DataFrame)
getprobeids(S::AbstractSession) = getprobes(S)[!, :id]
getchannels(S::AbstractSession) = CSV.read(IOBuffer(S.pyObject.channels.to_csv()), DataFrame)
function getchannels(S::AbstractSession, probeid)
    c = getchannels(S)
    c = subset(c, :probe_id=>ByRow(==(probeid)))
end
function getprobecoordinates(S::AbstractSession)
    c = subset(getchannels(S),              :anterior_posterior_ccf_coordinate => ByRow(!ismissing),
                                            :dorsal_ventral_ccf_coordinate => ByRow(!ismissing),
                                            :left_right_ccf_coordinate => ByRow(!ismissing))
    x = c[!, :anterior_posterior_ccf_coordinate]
    y = c[!, :dorsal_ventral_ccf_coordinate]
    z = c[!, :left_right_ccf_coordinate]
    return (x, y, z)
end
function getprobecoordinates(S::AbstractSession, probeid)
    c = subset(getchannels(S, probeid),              :anterior_posterior_ccf_coordinate => ByRow(!ismissing),
                                            :dorsal_ventral_ccf_coordinate => ByRow(!ismissing),
                                            :left_right_ccf_coordinate => ByRow(!ismissing))
    x = c[!, :anterior_posterior_ccf_coordinate]
    y = c[!, :dorsal_ventral_ccf_coordinate]
    z = c[!, :left_right_ccf_coordinate]
    return (x, y, z)
end
function getstructureacronyms(channelids::Vector{Int})
    channels = getchannels()
    acronyms = Vector{Any}(undef, size(channelids))
    [acronyms[i] = channels[channels.id.==channelids[i], :ecephys_structure_acronym][1] for i ∈ 1:length(channelids)]
    return acronyms
end

function getstructureids(channelids::Vector{Int})
    channels = getchannels()
    acronyms = Vector{Any}(undef, size(channelids))
    [acronyms[i] = channels[channels.id.==channelids[i], :ecephys_structure_id][1] for i ∈ 1:length(channelids)]
    return acronyms
end

function getstimuli(S::Session)
    str =  S.pyObject.stimulus_presentations.to_csv()
    CSV.read(IOBuffer(str), DataFrame);
end

function getunitmetrics(session::AbstractSession)
    str = session.pyObject.units.to_csv()
    CSV.read(IOBuffer(str), DataFrame);
end

function getstimulusname(session::AbstractSession, time::Number; stimulus_table=getstimuli(session))
    idx = findlast(stimulus_table.start_time .< time)
    if isnothing(idx)
        "blank"
    else
        stimulus_table.stimulus_name[idx]
    end
end
getstimulusname(session::AbstractSession, times; stimulus_table=getstimuli(session), kwargs...) = getstimulusname.([session], times; stimulus_table, kwargs...)

function getstimuli(S::Session, stimulusname::String)
    stimulus_table = getstimuli(S)
    df = subset(stimulus_table, :stimulus_name=>ByRow(==(stimulusname)))
end

function getstimuli(session::Session, times::Union{Tuple, UnitRange, LinRange, Vector})
    stimuli = getstimuli(session)
    idxs = [findfirst(time .< stimuli.stop_time) for time ∈ times] # Find first frame that ends after each time point
    return stimuli[idxs, :]
end

function getepochs(S::Session)
    p = S.pyObject.get_stimulus_epochs() # Why is this so slow
    CSV.read(IOBuffer(p.to_csv()), DataFrame);
end

function getepochs(S::Session, stimulusname)
    epoch_table = getepochs(S)
    df = subset(epoch_table, :stimulus_name=>ByRow(==(stimulusname)))
end


text:
#EasyPlotMPL base types & core functions
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

const scalemap = Dict{Symbol, String}(
	:lin => "linear",
	:log => "log",
)

const linestylemap = Dict{Symbol, String}(
	:none    => "none",
	:solid   => "-",
	:dash    => "--",
	:dot     => ":",
	:dashdot => "-.",
)

const markermap = Dict{Symbol, String}(
	:none    => "",
	:square    => "s",
	:diamond   => "D", #Diagsquare; diamond is d
	:uarrow    => "^", :darrow => "v",
	:larrow    => "<", :rarrow => ">",
	:cross     => "+", :+ => "+",
	:diagcross => "x", :x => "x",
	:circle    => "o", :o => "o",
	:star      => "*", :* => "*",
)

#==Base types
===============================================================================#
mutable struct WfrmBuilder <: EasyPlot.AbstractWfrmBuilder
	ref::PyCall.PyObject #"Axes" reference
	theme::EasyPlot.Theme
	fold::Optional{EasyPlot.FoldedAxis}
end

#Keeps key features of Matplotlib/PyPlot state needed for plotting
mutable struct MPLState
	interactive::Bool
	backend::Symbol
	guimode::Bool
end

mutable struct WfrmAttributes
	label
	color #linecolor
	linewidth
	linestyle
	marker
	markersize
	markerfacecolor
	markeredgecolor
	markeredgewidth
	fillstyle
end
WfrmAttributes(;label=nothing,
	color=nothing, linewidth=nothing, linestyle=nothing,
	marker=nothing, markersize=nothing, markerfacecolor=nothing,
	markeredgecolor=nothing) =
	WfrmAttributes(label, color, linewidth, linestyle,
		marker, markersize, markerfacecolor, markeredgecolor, linewidth,
		nothing == markerfacecolor ? "none" : "full"
	)

#==GUI state control
===============================================================================#
function _getstate()
	return MPLState(
		PyPlot.matplotlib.is_interactive(),
		PyPlot.pygui(),
		!PyPlot.isjulia_display[1], #TODO: Hack!... not part of PyPlot interface.
	)
end

function _applystate(s::MPLState)
	PyPlot.pygui(s.backend)
	PyPlot.pygui(s.guimode)
#@show pygui(), PyPlot.isjulia_display[1]
	#Must be applied last (in case backend does not support interactive):
	PyPlot.matplotlib.interactive(s.interactive)
end

#==Helper functions
===============================================================================#
const HEX_CODES = UInt8[
	'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F'
]
function int2mplcolorstr(v::UInt)
	result = Array{UInt8}(undef, 7) #6HEX+Hash symbol
	result[1] = '#'
	for i in length(result):-1:2
		result[i] = HEX_CODES[(v & 0xF)+1]
		v >>= 4
	end
	return String(result)
end

function mapcolor(v::Colorant)
	v = convert(RGB24, v)
	return int2mplcolorstr(UInt(v.color))
end
mapfacecolor(v) = mapcolor(v) #In case we want to diverge

#Linewidth:
maplinewidth(w) = w
maplinewidth(::Nothing) = maplinewidth(1) #default

#Marker size:
mapmarkersize(sz) = 5*sz
mapmarkersize(::Nothing) = mapmarkersize(1)

function maplinestyle(v::Symbol)
	result = get(linestylemap, v, missing)
	if ismissing(result)
		@info("Line style not supported")
		result = maplinestyle(nothing)
	end
	return result
end
maplinestyle(::Nothing) = "-" #default

function mapmarkershape(v::Symbol)
	result = get(markermap, v, missing)
	if ismissing(result)
		@info("Marker shape not supported")
		result = "o" #Use some supported marker
	end
	return result
end
mapmarkershape(::Nothing) = "" #default (no marker)

function WfrmAttributes(id::String, attr::EasyPlot.WfrmAttributes)
	markerfacecolor = attr.glyphfillcolor==colorant"transparent" ?
		nothing : mapfacecolor(attr.glyphfillcolor)

	return WfrmAttributes(label=id,
		color=mapcolor(attr.linecolor),
		linewidth=maplinewidth(attr.linewidth),
		linestyle=maplinestyle(attr.linestyle),
		marker=mapmarkershape(attr.glyphshape),
		markersize=mapmarkersize(attr.glyphsize),
		markerfacecolor=markerfacecolor,
		markeredgecolor = mapcolor(attr.glyphlinecolor),
	)
end

#==AbstractWfrmBuilder implementation
===============================================================================#
EasyPlot.needsfold(b::WfrmBuilder) = b.fold

#Add DataF1 results:
function _addwfrm(ax, d::DataF1, a::WfrmAttributes)
	kwargs = Any[]
	for attrib in fieldnames(typeof(a))
		v = getfield(a,attrib)

		if v != nothing
			push!(kwargs, tuple(attrib, v))
		end
	end

	wfrm = ax.plot(d.x, d.y; kwargs...)
end

#Called by EasyPlot, for each individual DataF1 ∈ DataMD.
function EasyPlot.addwfrm(b::WfrmBuilder, d::DataF1, id::String,
	la::EasyPlot.LineAttributes, ga::EasyPlot.GlyphAttributes, strip::Int)
	attr = EasyPlot.WfrmAttributes(b.theme, la, ga) #Apply theme to attributes
	mplattr = WfrmAttributes(id, attr) #Attributes understood by MPL
	_addwfrm(b.ref, d, mplattr)
end

#Last line


text:
#=
- Open this script from Atom
- Install package uber-Juno
- restart Atom and open again
- cmd+shift+enter will start
=#
using Plots
pyplot()
plot(rand(4,4))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
mutable struct Reaction

    ec_number::String
    kegg_reaction_number::String
    kegg_enzyme_name::String
    kegg_reaction_markup::String

    function Reaction()
		this = new()
	end
end

mutable struct Organism

    organism_id::String
    organism_code::String
    species_description::String
    species_taxonomy::String

    function Organism()
		this = new()
	end
end

mutable struct Sequence

    gene_location::String
    type::Symbol
    body::Array{SubString{String},1}
    header::String

    function Sequence()
		this = new()
	end
end


text:
module PlotPlants

import FileIO: load, save
import ImageIO: load, save
import ImageMagick: load, save

using DataFrames: DataFrame, sort!
using Distributed: @everywhere, addprocs, pmap, workers
using DocStringExtensions: TYPEDFIELDS
using GLM: @formula, adjr2, coef, coeftable, confint, lm, predict
using KernelDensity: kde, pdf
using PkgUtility: read_nc
using ProgressMeter: @showprogress
using PyCall: PyObject, pyimport
using PyPlot: Figure, clf, figure, rc, subplot
using Random: randstring

# export public types
export FormatNC, FormatTIFF, LinearRegressionResult

# export canvas related functions
export create_canvas, save_canvas!, save_gif!, set_titles!, set_xlabels!,
       set_xlims!, set_xticklabels!, set_xticks!, set_xylabels!, set_xylims!,
       set_xyticklabels!, set_xyticks!, set_ylabels!, set_ylims!,
       set_yticklabels!, set_yticks!

# export functions related to plotting
export plot_density!, plot_ellipse!, plot_hexbin!, plot_line_regress!,
       plot_stoma!, preview_data, preview_dataset!

# regression related functions
export line_regress

# statistics related function
export calculate_density, test_slope

# string related functions
export latex_symbol, latex_unit, parse_symbol, subscript, superscript

# text font and render functions
export use_sans, use_sans_tex, use_sans_text, use_serif, use_serif_tex,
       use_serif_text

# include the files
include("types/file_format.jl")
include("types/regression.jl" )

include("canvas/canvas.jl"  )
include("canvas/gif.jl"     )
include("canvas/titles.jl"  )
include("canvas/xylabels.jl")
include("canvas/xylims.jl"  )
include("canvas/xyticks.jl" )

include("plots/density.jl"     )
include("plots/ellipse.jl"     )
include("plots/hexbin.jl"      )
include("plots/line_regress.jl")
include("plots/preview.jl"     )
include("plots/stoma.jl"       )

include("regression/linear.jl")

include("statistics/density.jl")
include("statistics/slope.jl"  )

include("string/latex_symbol.jl")
include("string/latex_unit.jl"  )
include("string/parse_symbol.jl")
include("string/sub_super.jl"   )

include("text/font.jl"  )
include("text/render.jl")

end # module


text:
# complex dense, projected away from origin
using Pseudospectra, Test

@testset "Inverse projection" begin
    opts = Dict{Symbol,Any}(:ax => [-1.2,1.2,-1.1,1.1], :npts => 60,
                            :levels => -5:0.5:-1, :proj_lev => -5.0)
    A = Pseudospectra.landau_fox_li(225,16)

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
end


text:
module ANOVA

include("anova_table.jl")
export ANOVAtest,
  effects,
  anova

end # module


text:
include("rijke.jl")
using PyPlot
N = 18
u0Init = zeros(N)
u0Init[1] = 1.
s = [7.0, 0.2]
nRunUp = 40000
#u0 = Rijke(u0Init, s, nRunUp)
t = LinRange(-1.01,-0.99,500)
heat_sqrt = qfun_sqrt.(t)
heat_poly = qfun_poly.(t)
fig, ax = subplots(1,1)
ax.plot(t, heat_sqrt, label=L"$\dot{q}$, using sqrt", ".")
ax.plot(t, heat_poly, label=L"$\dot{q}$, using poly", ".")
ax.xaxis.set_tick_params(labelsize=18)
ax.yaxis.set_tick_params(labelsize=18)
ax.set_xlabel("delayed velocity",fontsize=18)
fig.legend(fontsize=18)

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
#
# bradi.jl --
#
# Find a global minimum of an univariate function with the "Bracket then Dig"
# method described in:
#
#     Ferréol Soulez, Éric Thiébaut, Michel Tallon, Isabelle Tallon-Bosc
#     and Paulo Garcia, "Optimal a posteriori fringe tracking in optical
#     interferometry", Proc. SPIE 9146, Optical and Infrared Interferometry
#     IV, 91462Y (July 24, 2014); doi:10.1117/12.2056590
#
#-------------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2014-2020, Éric Thiébaut.
#

module BraDi

using ..Brent: fmin_atol, fmin_rtol, fminbrkt

"""
    BraDi.minimize([T=Float64,] f, x) -> (xbest, fbest)
    BraDi.maximize([T=Float64,] f, x) -> (xbest, fbest)

finds the global minimum (resp. maximum) of an univariate function `f`.  The
argument `x` is an abstract vector of coordinates in monotonic order; `x[1]`
and `x[end]` are the endpoints of the global search interval and it is assumed
that no more than a single local minimum lies in any subinterval
`[x[i],x[i+2]]`.  The result is a tuple of 2-values: `xbest` the position of
the global minimum and `fbest` the function value at this position.

Optional argument `T` is the floating-point type used for the computations.

If specified, keywords `atol` and `rtol` set the absolute and relative
tolerances for the precision (see `Brent.fmin`).

For instance:

    BraDi.minimize(f, range(a,stop=b,length=n))

performs a global search in the closed interval `[a,b]` which is sampled by `n`
equalspaced points.

These functions implement the BraDi ("Bracket" then "Dig") algorithm described
in:

> Ferréol Soulez, Éric Thiébaut, Michel Tallon, Isabelle Tallon-Bosc and Paulo
> Garcia, "Optimal a posteriori fringe tracking in optical interferometry",
> Proc. SPIE 9146, Optical and Infrared Interferometry IV, 91462Y (July 24,
> 2014); doi:10.1117/12.2056590

See also `Brent.fmin` and `Step.minimize`.

"""
minimize(f, x::AbstractVector{<:Real}; kdws...) =
    minimize(Float64, f, x; kdws...)

maximize(f, x::AbstractVector{<:Real}; kdws...) =
    maximize(Float64, f, x; kdws...)

function minimize(::Type{T}, f, x::AbstractVector{<:Real};
                  atol::Real = fmin_atol(T),
                  rtol::Real = fmin_rtol(T)) where {T<:AbstractFloat}
    _minimize(T, f, x, T(atol), T(rtol))
end

function maximize(::Type{T}, f, x::AbstractVector{<:Real};
                  kdws...) where {T<:AbstractFloat}
    (xbest, fbest) = minimize(T, t -> -f(t), x; kdws...)
    return (xbest, -fbest)
end

function _minimize(::Type{T}, f, x::AbstractVector{<:Real},
                   atol::T, rtol::T) where {T<:AbstractFloat}
    xbest = xa = xb = xc = T(x[1])
    fbest = fa = fb = fc = T(f(xc))
    n = length(x)
    for j in 2 : n + 1
        xa, xb = xb, xc
        fa, fb = fb, fc
        if j ≤ n
            xc = T(x[j])
            fc = T(f(xc))
            if fc < fbest
                xbest = xc
                fbest = fc
            end
        end
        if fa ≥ fb ≤ fc
            # A minimum has been bracketed in [XA,XC].
            xm, fm = fminbrkt(f, xb, fb, xa, fa, xc, fc, atol, rtol)
            if fm < fbest
                xbest = xm
                fbest = fm
            end
        end
    end
    return (xbest, fbest)
end

@doc @doc(minimize) maximize

end # module


text:
## This is the OLD version of Table. We are using the version in lists_new.jl

# We only do Table(expr,[i,imax])
# Our Table is rather slow. Slower than Maxima makelist.
# Table(a(i),[i,10000]) is 2 to 3 x slower than makelist( i, i, 1, 10000)
# If we cheat and only do a single evaluation,
# and setfixed() then the speed is the same.
# But, the 3rd party 'table' for Maxima is 4 times faster than makelist
# in this example.
# Mma 3.0  Timing[mm = Table[a[i],{i,10^5}];]  {0.05 Second, Null}
# 3rd pary maxima: table(a(i),[i,10^5]);  0.14 s
# At commit 3478f83ee403238704ad6af81b42e1c7f4d07cc8
# Table(a(i),[i,10000]) takes bout 0.45--0.49 s, makelist reports 0.503.

set_attribute(:TableOld, :HoldAll)

function apprules(mx::Mxpr{:TableOld})
    expr = mx[1]
    iter = mx[2]
    isym = get_localized_symbol(iter[1])
    imax = meval(iter[2])
#    issym = createssym(isym,Int)  ## Trying out Typed variable
    ex = replsym(deepcopy(expr),iter[1],isym) # takes no time, for simple expression
    args = do_table(imax,isym,ex) # creating a symbol is pretty slow
    removesym(isym)
    mx1 = mxpr(:List,args) # takes no time
    # This mergesyms no longer makes this much faster, because overall Table has slowed by more than 5x
    # syms are being merged alot somewhere else
    mergesyms(mx1,:nothing) # not correct, but stops the merging
    setcanon(mx1)
    setfixed(mx1)
    return mx1
end

# These tests concern not Table per se, but the efficiency
# of evaluation, in particular assigning values to symbols.
# commit 3478f83ee403238704ad6af81b42e1c7f4d07cc8
# Testing Table(a(i),[i,10^5])
# changes                   Time
# normal                    0.49
# disable setsymval         0.35
#
# commit 0d2332b17eb7e5c518a86b1b3044b691fddb5b87
# setsymval is now much faster because we changed data structure in SSJSym
# Testing Table(a(i),[i,10^5])
# normal                    0.36
# *sigh* Cannot reproduce the performance in the previous line. Even if
# I return to that commit. Time is closer to 0.5 than 0.35
# setsymval still seems to be an expensive operation
#
# 4f7d9f6dff207ae2c95736ab058627f503fbad26
# With this commit, Time is 1.83 s. For same code as above.  Table is
# becoming slower. We try using a stripped-down meval1 from
# alteval.jl, time is .43 with no gc.

function do_table(imax::T,isym,ex) where T<:Integer
    args = newargs(imax)
    sisym = getssym(isym)
#    setsymval(sisym,1)
    @inbounds for i in 1:imax
        setsymval(sisym,i)
#        args[i] = meval(ex)
        args[i] = doeval(ex)
#        args[i] = meval1(ex)
        setfixed(args[i])  # no difference ?
    end
    return args
end

# The speed penalty for val::Any instead of val::Int, is factor of 10^4.
# Setting isym.val is by far the slowest step in Table(i,[i,10^6])
# Much slower than Mma and SBCL
# If we us val::Array{Any,1}, the speed penalty is better by a couple
# orders of magnitude

# type NSYM
# #    val::Int
#     val::Array{Any,1}
# end

# function testset()
#     ss = NSYM(Any[0])
#     sum0 = 0
#     for i in 1:10^6
#         ss.val[1] = i
#         sum0 += ss.val[1]
#     end
#     sum0
# end


text:
# ndgrid.jl

using MIRT: ndgrid

using Test: @test, @inferred

@inferred ndgrid(1:3, 2:4)
@test ndgrid(1:3, 2:4) == ([1 1 1; 2 2 2; 3 3 3], [2 3 4; 2 3 4; 2 3 4])
@test size(ndgrid(1:3, 4:5, 6:9)[3]) == (3,2,4)


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

"""
    Variogram

A variogram model (e.g. Gaussian variogram).
"""
abstract type Variogram{T,D} end

"""
    isstationary(γ)

Check if variogram `γ` possesses the 2nd-order stationary property.
"""
isstationary(::Variogram) = false

"""
    range(γ)

Return the range of the variogram `γ` when defined.
"""
Base.range(γ::Variogram) = γ.range

"""
    sill(γ)

Return the sill of the variogram `γ` when defined.
"""
sill(γ::Variogram) = γ.sill

"""
    nugget(γ)

Return the nugget of the variogram `γ` when defined.
"""
nugget(γ::Variogram) = γ.nugget

"""
    param_type(γ)

Return the parameter (e.g. sill, range) type of the variogram.
"""
param_type(::Variogram{T,D}) where {T<:Real,D<:Metric} = T

"""
    result_type(γ, x₁, x₂)

Return result type of γ(x₁, x₂).
"""
result_type(γ::Variogram, x₁::AbstractArray, x₂::AbstractArray) =
  promote_type(param_type(γ), Distances.result_type(γ.distance, x₁, x₂))

#------------------
# IMPLEMENTATIONS
#------------------
"""
    GaussianVariogram(sill=s, range=r, nugget=n, distance=d)

A Gaussian variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct GaussianVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 1e-6 # positive nugget for numerical stability
  distance::D = Euclidean()
end
(γ::GaussianVariogram)(h) = (γ.sill - γ.nugget) * (1 - exp(-3(h/γ.range)^2)) + (h > 0) * γ.nugget
(γ::GaussianVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::GaussianVariogram) = true

"""
    ExponentialVariogram(sill=s, range=r, nugget=n, distance=d)

An exponential variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct ExponentialVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  distance::D = Euclidean()
end
(γ::ExponentialVariogram)(h) = (γ.sill - γ.nugget) * (1 - exp(-3(h/γ.range))) + (h > 0) * γ.nugget
(γ::ExponentialVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::ExponentialVariogram) = true

"""
    MaternVariogram(sill=s, range=r, nugget=n, order=ν, distance=d)

A Matérn variogram with sill `s`, range `r` and nugget `n`. The parameter
ν is the order of the Bessel function. Optionally, use a custom distance `d`.
"""
@with_kw struct MaternVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  order::T  = 1.
  distance::D = Euclidean()
end
(γ::MaternVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget
  ν = γ.order

  # shift lag by machine precision to
  # avoid explosion at the origin
  h2 = h + eps(eltype(h))
  h3 = sqrt(2.0ν)h2/r

  (s - n) * (1 - 2.0^(1 - ν)/gamma(ν) * h3^ν * besselk(ν, h3)) + (h > 0) * n
end
(γ::MaternVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::MaternVariogram) = true

"""
    SphericalVariogram(sill=s, range=r, nugget=n, distance=d)

A spherical variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct SphericalVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  distance::D = Euclidean()
end
(γ::SphericalVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * (1.5(h/r) - 0.5(h/r)^3) + (h ≥ r) * (s - n) + (h > 0) * n
end
(γ::SphericalVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::SphericalVariogram) = true

"""
    CubicVariogram(sill=s, range=r, nugget=n, distance=d)

A cubic variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct CubicVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  distance::D = Euclidean()
end
(γ::CubicVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * (7*(h/r)^2 - (35/4)*(h/r)^3 + (7/2)*(h/r)^5 - (3/4)*(h/r)^7) +
  (h ≥ r) * (s - n) + (h > 0) * n
end
(γ::CubicVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::CubicVariogram) = true

"""
    PentasphericalVariogram(sill=s, range=r, nugget=n, distance=d)

A pentaspherical variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct PentasphericalVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  distance::D = Euclidean()
end
(γ::PentasphericalVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * ((15/8)*(h/r) - (5/4)*(h/r)^3 + (3/8)*(h/r)^5) +
  (h ≥ r) * (s - n) + (h > 0) * n
end
(γ::PentasphericalVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::PentasphericalVariogram) = true

"""
    PowerVariogram(scaling=s, exponent=a, nugget=n, distance=d)

A power variogram with scaling `s`, exponent `a` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct PowerVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  scaling::T  = 1.
  nugget::T   = 0.
  exponent::T = 1.
  distance::D = Euclidean()
end
(γ::PowerVariogram)(h) = @. γ.scaling*h^γ.exponent + (h > 0) * γ.nugget
(γ::PowerVariogram)(x, y) = γ(evaluate(γ.distance, x, y))

"""
    SineHoleVariogram(sill=s, range=r, nugget=n, distance=d)

A sine hole variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct SineHoleVariogram{T<:Real,D<:Metric} <: Variogram{T,D}
  range::T  = 1.
  sill::T   = 1.
  nugget::T = 0.
  distance::D = Euclidean()
end
(γ::SineHoleVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  # shift lag by machine precision to
  # avoid explosion at the origin
  h = h + eps(eltype(h))

  (s - n) * (1 - sin(π*h/r)/(π*h/r)) + (h > 0) * n
end
(γ::SineHoleVariogram)(x, y) = γ(evaluate(γ.distance, x, y))
isstationary(::SineHoleVariogram) = true

#------------------------
# COMPOSITIVE VARIOGRAMS
#------------------------
"""
    CompositeVariogram(γ₁, γ₂, ..., γₙ)

A composite (additive) model of variograms γ(h) = γ₁(h) + γ₂(h) + ⋯ + γₙ(h).
"""
struct CompositeVariogram <: Variogram{Real,Metric}
  γs::Vector{Variogram}
  CompositeVariogram(g, gs...) = new([g, gs...])
end
(c::CompositeVariogram)(h) = sum(γ(h) for γ in c.γs)
(c::CompositeVariogram)(x, y) = sum(γ(x,y) for γ in c.γs)
isstationary(c::CompositeVariogram) = all(isstationary(γ) for γ in c.γs)
sill(c::CompositeVariogram) = sum(sill(γ) for γ in c.γs)
nugget(c::CompositeVariogram) = sum(nugget(γ) for γ in c.γs)
param_type(c::CompositeVariogram) = promote_type([param_type(γ) for γ in c.γs]...)

"""
    γ₁ + γ₂

Return compositive (additive) model of variograms `γ₁` and `γ₂`.
"""
+(γ₁::Variogram, γ₂::Variogram) = CompositeVariogram(γ₁, γ₂)
+(c::CompositeVariogram, γ::Variogram) = CompositeVariogram(c.γs..., γ)
+(γ::Variogram, c::CompositeVariogram) = CompositeVariogram(γ, c.γs...)
+(c₁::CompositeVariogram, c₂::CompositeVariogram) = CompositeVariogram(c₁.γs..., c₂.γs...)


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/SAG.jl")
include("../opt/SAGA.jl")
include("../opt/SAGA_PROX.jl")
include("../opt/SGD.jl")
include("../data/data.jl")

using ProximalOperators
using Gadfly
using DataFrames
using ArgParse

#for lambda in [1e-2,1e-3,1e-4,1e-5]
#    for step in [0.1, 0.05, 0.01, 0.005, 0.001]

function main(args)

    lambda = 0
    step = 0
    iter = 300
    mode = "logreg"

    s = ArgParseSettings(description = "Example 1 for argparse.jl: minimal usage.")

    @add_arg_table s begin
        "--step","-s"
        "--lam", "-l"
        "--iter"
        "--mode"
    end
    parsed_args = parse_args(s)
    println("Parsed args:")
    for (key,val) in parsed_args
        if key == "lam"
            lambda = val
        elseif key == "step"
            step = val
        elseif key == "iter"
            iter = val
        elseif key == "mode"
            mode = val
        end
        println("  $key  =>  $(repr(val))")
    end
    println(lambda)
    println(step)
    lambda = parse(Float64,lambda)
    step = parse(Float64,step)
    iter = parse(Int32, iter)
    X, y = getFullData()
    (X_test,y_test) = getTestData()

    # init initial point x_0
    x0 = zeros(size(X,2),1)

    plotname = string("comparison-lam=", lambda,  "-step=",step, "-iter=",iter, "-mode=", mode ,".svg" )
    println(plotname)

    if mode == "logreg"
    f = f_logreg(X,y,x0,lambda)
    g = sub_g_logreg()
    elseif mode == "svm"
        f = f_svm(X,y,x0,lambda)
        g = sub_g_svm()
    else
        error("wrong mode!")
    end

    n = size(X,1)

    prox_operator = NormL2(lambda)

    @time (theta_sag, strings, iter,val_array_sag) = SAG(X,y,x0,f,g,1e-8,lambda,iter,step,n)

    @time (theta_saga, strings, iter,val_array_saga) = SAGA(X,y,x0,f,g,1e-8,lambda,iter,step,n)

    @time (theta_saga_prox, strings, iter,val_array_saga_prox) = SAGA(X,y,x0,f,g,1e-8,lambda,iter,step,n,prox_operator)

    @time (theta_sgd,state,val_array_sgd,stop_array_sgd) = sgd(X,y,x0,f,g,1e-8,lambda,iter,step,n)

    # EVALUATION #

    # evaluate on training set #
    @printf("Accuracy of sag on trainingdata = %f.\n",evaluate(X,y,theta_sag,mode,false))
    @printf("Accuracy of saga on trainingdata = %f.\n",evaluate(X,y,theta_saga,mode,false))
    @printf("Accuracy of saga_prox on trainingdata = %f.\n",evaluate(X,y,theta_saga_prox,mode,false))
    @printf("Accuracy of sgd on trainingdata = %f.\n",evaluate(X,y,theta_sgd,mode,false))

    println()

    # evalute test dat #
    @printf("Accuracy of sag on testdata = %f.\n",evaluate(X_test,y_test,theta_sag,mode,false))
    @printf("Accuracy of saga on testdata = %f.\n",evaluate(X_test,y_test,theta_saga,mode,false))
    @printf("Accuracy of saga_prox on testdata = %f.\n",evaluate(X_test,y_test,theta_saga_prox,mode,false))
    @printf("Accuracy of sgd on testdata = %f.\n",evaluate(X_test,y_test,theta_sgd,mode,false))

    # plot values

    df = DataFrame(iter = [identity(x) for x in 1:iter],
                    obj = val_array_saga,
                    obj2 = val_array_saga_prox,
                    obj3 = val_array_sgd,
                    obj4 = val_array_sag)

    p = plot(layer(df,x=:iter,y=:obj4,Geom.line,Theme(default_color=colorant"green")),
            layer(df,x=:iter,y=:obj,Geom.line,Theme(default_color=colorant"red")),
            layer(df,x=:iter,y=:obj2,Geom.line,Theme(default_color=colorant"blue")),
            layer(df,x=:iter,y=:obj3,Geom.line,Theme(default_color=colorant"orange")))
    img = SVG(plotname, 6inch, 6inch)
    draw(img, p)

end

main(ARGS)

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
using MaxwellBase
using SimpleConstants
using Test
using Statistics: mean
using AbbreviatedTypes

Base.isapprox(a::Tuple, b::Tuple; kws...) = all(p -> isapprox(p...; kws...), zip(a,b))

# @testset "MaxwellBase" begin

include("enumtype.jl")
include("util.jl")
include("field.jl")
include("material.jl")
include("object.jl")
include("smoothing.jl")
include("param.jl")

# end  # @testset "MaxwellBase"


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of line breaks.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Line breaks inside cells" begin
    data = ["This line contains\nthe velocity [m/s]" 10.0;
            "This line contains\nthe acceleration [m/s^2]" 1.0;
            "This line contains\nthe time from the\nbeginning of the simulation" 10;]

    header = ["Information", "Value"]

    expected = """
┌─────────────────────────────┬───────┐
│                 Information │ Value │
├─────────────────────────────┼───────┤
│          This line contains │  10.0 │
│          the velocity [m/s] │       │
│          This line contains │   1.0 │
│    the acceleration [m/s^2] │       │
│          This line contains │    10 │
│           the time from the │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(String, data; header = header, linebreaks = true)
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│         Information         │ Value │
├─────────────────────────────┼───────┤
│     This line contains      │ 10.0  │
│     the velocity [m/s]      │       │
│     This line contains      │  1.0  │
│  the acceleration [m/s^2]   │       │
│     This line contains      │  10   │
│      the time from the      │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(
        String,
        data;
        header = header,
        alignment = :c,
        linebreaks = true
    )
    @test result == expected

    expected = """
┌─────────────────────────────┬───────┐
│ Information                 │ Value │
├─────────────────────────────┼───────┤
│ This line contains          │ 10.0  │
│ the velocity [m/s]          │       │
│ This line contains          │ 1.0   │
│ the acceleration [m/s^2]    │       │
│ This line contains          │ 10    │
│ the time from the           │       │
│ beginning of the simulation │       │
└─────────────────────────────┴───────┘
"""

    result = pretty_table(
        String,
        data;
        header = header,
        alignment = :l,
        linebreaks = true
    )
    @test result == expected

    expected = """
┌────────────────────────────────────────────────────────────────────┬───────┐
│                                                        Information │ Value │
├────────────────────────────────────────────────────────────────────┼───────┤
│                             This line contains\\nthe velocity [m/s] │  10.0 │
│                       This line contains\\nthe acceleration [m/s^2] │   1.0 │
│ This line contains\\nthe time from the\\nbeginning of the simulation │    10 │
└────────────────────────────────────────────────────────────────────┴───────┘
"""

    result = pretty_table(String, data; header = header)
    @test result == expected

    # Show only the first line
    # --------------------------------------------------------------------------

    expected = """
┌────────────────────┬───────┐
│        Information │ Value │
├────────────────────┼───────┤
│ This line contains │  10.0 │
│ This line contains │   1.0 │
│ This line contains │    10 │
└────────────────────┴───────┘
"""

    result = pretty_table(
        String,
        data;
        header = header,
        cell_first_line_only = true
    )
    @test result == expected
end

@testset "Auto wrapping" begin
    table = [1 """Ouviram do Ipiranga as margens plácidas
                  De um povo heróico o brado retumbante,
                  E o sol da Liberdade, em raios fúlgidos,
                  Brilhou no céu da Pátria nesse instante.""";
             2 """Se o penhor dessa igualdade
                  Conseguimos conquistar com braço forte,
                  Em teu seio, ó Liberdade,
                  Desafia o nosso peito a própria morte!""";
             3 """Ó Pátria amada, Idolatrada, Salve! Salve!
                  Brasil, um sonho intenso, um raio vívido
                  De amor e de esperança à terra desce,
                  Se em teu formoso céu, risonho e límpido,"""]

    header = ["Verse number", "Verse"]

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │                          Verse │
├──────────────┼────────────────────────────────┤
│            1 │         Ouviram do Ipiranga as │
│              │               margens plácidas │
│              │     De um povo heróico o brado │
│              │                    retumbante, │
│              │       E o sol da Liberdade, em │
│              │                raios fúlgidos, │
│              │       Brilhou no céu da Pátria │
│              │                nesse instante. │
├──────────────┼────────────────────────────────┤
│            2 │    Se o penhor dessa igualdade │
│              │     Conseguimos conquistar com │
│              │                   braço forte, │
│              │      Em teu seio, ó Liberdade, │
│              │        Desafia o nosso peito a │
│              │                 própria morte! │
├──────────────┼────────────────────────────────┤
│            3 │    Ó Pátria amada, Idolatrada, │
│              │                  Salve! Salve! │
│              │   Brasil, um sonho intenso, um │
│              │                    raio vívido │
│              │       De amor e de esperança à │
│              │                   terra desce, │
│              │         Se em teu formoso céu, │
│              │             risonho e límpido, │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(
        String,
        table;
        header        = header,
        autowrap      = true,
        linebreaks    = true,
        body_hlines   = [1, 2],
        columns_width = [-1, 30]
    )

    @test result == expected

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │             Verse              │
├──────────────┼────────────────────────────────┤
│      1       │     Ouviram do Ipiranga as     │
│              │        margens plácidas        │
│              │   De um povo heróico o brado   │
│              │          retumbante,           │
│              │    E o sol da Liberdade, em    │
│              │        raios fúlgidos,         │
│              │    Brilhou no céu da Pátria    │
│              │        nesse instante.         │
├──────────────┼────────────────────────────────┤
│      2       │  Se o penhor dessa igualdade   │
│              │   Conseguimos conquistar com   │
│              │          braço forte,          │
│              │   Em teu seio, ó Liberdade,    │
│              │    Desafia o nosso peito a     │
│              │         própria morte!         │
├──────────────┼────────────────────────────────┤
│      3       │  Ó Pátria amada, Idolatrada,   │
│              │         Salve! Salve!          │
│              │  Brasil, um sonho intenso, um  │
│              │          raio vívido           │
│              │    De amor e de esperança à    │
│              │          terra desce,          │
│              │     Se em teu formoso céu,     │
│              │       risonho e límpido,       │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(
        String,
        table;
        header        = header,
        alignment     = :c,
        autowrap      = true,
        linebreaks    = true,
        body_hlines   = [1, 2],
        columns_width = [-1, 30]
    )

    @test result == expected

    expected = """
┌──────────────┬────────────────────────────────┐
│ Verse number │ Verse                          │
├──────────────┼────────────────────────────────┤
│ 1            │ Ouviram do Ipiranga as         │
│              │ margens plácidas               │
│              │ De um povo heróico o brado     │
│              │ retumbante,                    │
│              │ E o sol da Liberdade, em       │
│              │ raios fúlgidos,                │
│              │ Brilhou no céu da Pátria       │
│              │ nesse instante.                │
├──────────────┼────────────────────────────────┤
│ 2            │ Se o penhor dessa igualdade    │
│              │ Conseguimos conquistar com     │
│              │ braço forte,                   │
│              │ Em teu seio, ó Liberdade,      │
│              │ Desafia o nosso peito a        │
│              │ própria morte!                 │
├──────────────┼────────────────────────────────┤
│ 3            │ Ó Pátria amada, Idolatrada,    │
│              │ Salve! Salve!                  │
│              │ Brasil, um sonho intenso, um   │
│              │ raio vívido                    │
│              │ De amor e de esperança à       │
│              │ terra desce,                   │
│              │ Se em teu formoso céu,         │
│              │ risonho e límpido,             │
└──────────────┴────────────────────────────────┘
"""

    result = pretty_table(
        String,
        table;
        header        = header,
        alignment     = :l,
        autowrap      = true,
        linebreaks    = true,
        body_hlines   = [1, 2],
        columns_width = [-1, 30]
    )

    @test result == expected

    # Test with additional rows
    # --------------------------------------------------------------------------

    expected = """
┌─────┬──────────┬──────────────┬────────────────────────────────┐
│ Row │    Verso │ Verse number │ Verse                          │
├─────┼──────────┼──────────────┼────────────────────────────────┤
│   1 │ Primeiro │ 1            │ Ouviram do Ipiranga as         │
│     │          │              │ margens plácidas               │
│     │          │              │ De um povo heróico o brado     │
│     │          │              │ retumbante,                    │
│     │          │              │ E o sol da Liberdade, em       │
│     │          │              │ raios fúlgidos,                │
│     │          │              │ Brilhou no céu da Pátria       │
│     │          │              │ nesse instante.                │
├─────┼──────────┼──────────────┼────────────────────────────────┤
│   2 │  Segundo │ 2            │ Se o penhor dessa igualdade    │
│     │          │              │ Conseguimos conquistar com     │
│     │          │              │ braço forte,                   │
│     │          │              │ Em teu seio, ó Liberdade,      │
│     │          │              │ Desafia o nosso peito a        │
│     │          │              │ própria morte!                 │
├─────┼──────────┼──────────────┼────────────────────────────────┤
│   3 │ Terceiro │ 3            │ Ó Pátria amada, Idolatrada,    │
│     │          │              │ Salve! Salve!                  │
│     │          │              │ Brasil, um sonho intenso, um   │
│     │          │              │ raio vívido                    │
│     │          │              │ De amor e de esperança à       │
│     │          │              │ terra desce,                   │
│     │          │              │ Se em teu formoso céu,         │
│     │          │              │ risonho e límpido,             │
└─────┴──────────┴──────────────┴────────────────────────────────┘
"""

    result = pretty_table(
        String, table;
        header                = header,
        alignment             = :l,
        autowrap              = true,
        linebreaks            = true,
        body_hlines           = [1, 2],
        columns_width         = [-1, 30],
        show_row_number       = true,
        row_names             = ["Primeiro", "Segundo", "Terceiro"],
        row_name_column_title = "Verso"
    )

    @test result == expected
end


text:
function module_from_string(mod)
    ms = split(mod, '.')

    out = Main

    loaded_module = findfirst(==(first(ms)), string.(Base.loaded_modules_array()))

    if loaded_module !== nothing
        out = Base.loaded_modules_array()[loaded_module]
        popfirst!(ms)
    end

    for m in Symbol.(ms)
        if isdefined(out, m)
            resolved = getfield(out, m)

            if resolved isa Module
                out = resolved
            else
                return out
            end
        end
    end

    return out
end

is_module_loaded(mod) = mod == "Main" || module_from_string(mod) !== Main

function get_modules(toplevel=nothing, mods=Set(Module[]))
    top_mods = toplevel === nothing ? Base.loaded_modules_array() : [toplevel]

    for mod in top_mods
        push!(mods, mod)

        for name in names(mod, all=true)
            if !Base.isdeprecated(mod, name) && isdefined(mod, name)
                thismod = getfield(mod, name)
                if thismod isa Module && thismod !== mod && !(thismod in mods)
                    push!(mods, thismod)
                    get_modules(thismod, mods)
                end
            end
        end
    end
    mods
end

max_stars_repo_name: fonsp/julia-vscode

max_issues_repo_name: fonsp/julia-vscode

max_forks_repo_name: fonsp/julia-vscode


text:
__precompile__(true)
module SatelliteTasking

# Julia Packages
using Printf
using JSON
using UUIDs
using Distributed
using LinearAlgebra
using Statistics
using Distributions
using Discretizers
using Parameters

# POMDP Solvers
using Random
using POMDPs
using MCTS

# Package imports
using PGFPlots
using JuMP
using Gurobi

# Other dependencies
using SatelliteDynamics

# Includes
include("data_structures.jl")
include("models.jl")
include("utils.jl")
include("access.jl")
include("graph.jl")
include("milp.jl")
include("mdpcore.jl")
include("mdpfs.jl")
include("mdpmcts.jl")
include("baseline.jl")
include("analysis.jl")
include("visualization.jl")
include("simulation.jl")

# Export satellite planning submodule
# include(joinpath(".", "satellite_planning", "satellite_planning.jl"))

# Export analysis submodule
# include(joinpath(".", "analysis", "analysis.jl"))
# @reexport using SatelliteTasking.Analysis

end # End module


text:
module CMPUtils
    using Taro
    Taro.init()
println("Howdy!!!, Sit back and relax! This might take some time...")
    println("Loading the mighty plotting functions")
    using GRUtils
    println("Let me get some tools to handle images. ")
    using Images
    using ImageIO
    using ImageMagick
    using LinearAlgebra
    using Pipe: @pipe
    using ProgressMeter: @showprogress
    using CSV
    using DataFrames: DataFrame
    using StatsBase
    println("Did you hear something? Tuning up my audio powers!!!!")
    using WAV
    using MFCC: mfcc, deltas
    using Glob
    println("Oh, I almost forgot about the PDFs!")
    using TextAnalysis
    using Languages
    using StatsBase
    include("plotting.jl")
    include("recode.jl")
    include("classdemo.jl")

    export recodeimage, recodeaudio, recodetext, img_similarity, recodeaudio2
    precompile(recodeimage, (String,Int),)
    precompile(recodeaudio, (String,Int),)
    precompile(recodetext, (String,Int),)
    println("CMPUtils reporting for duty!")
end


text:
raw_data = readlines("../inputs/day03.txt")
num_bits = length(raw_data[1])

function line2ints(line)
    return map(x -> parse(Int, x), collect(line))
end

function get_power_consumption(data)

    function get_gamma(average_vals)
        output_string = ""
        for val in average_vals
            if val < 0.5
                output_string = output_string * "0"
            else
                output_string = output_string * "1"
            end
        end
        return parse(Int, output_string, base = 2)
    end

    function get_epsilon(average_vals)
        output_string = ""
        for val in average_vals
            if val < 0.5
                output_string = output_string * "1"
            else
                output_string = output_string * "0"
            end
        end
        return parse(Int, output_string, base = 2)
    end

    column_sums = sum(data)
    average_value = column_sums / length(data)

    gamma = get_gamma(average_value)
    epsilon = get_epsilon(average_value)

    return gamma * epsilon

end

int_data = map(x -> line2ints(x), raw_data)

power_consumption = get_power_consumption(int_data)
println("Day 03 Part 1 Solution: ", power_consumption)

function get_most_common_in_col(data, index)
    average_val = sum(data) / length(data)
    if average_val[index] >= 0.5
        return 1
    else
        return 0
    end
end

function get_least_common_in_col(data, index)
    average_val = sum(data) / length(data)
    if average_val[index] >= 0.5
        return 0
    else
        return 1
    end
end

function filter_by_most_common_index(data, index)
    filter_val = get_most_common_in_col(data, index)
    return data[map(x -> x[index] == filter_val, data)]
end

function filter_by_least_common_index(data, index)
    filter_val = get_least_common_in_col(data, index)
    return data[map(x -> x[index] == filter_val, data)]
end

function get_oxygen_rating(data)
    output_data = data
    for ii in range(1, num_bits)
        if length(output_data) == 1
            break
        else
            output_data = filter_by_most_common_index(output_data, ii)
        end
    end

    output_string = ""
    for bit in output_data[1]
        output_string = output_string * string(bit)
    end
    return parse(Int, output_string, base = 2)
end

function get_co2_rating(data)
    output_data = data
    for ii in range(1, num_bits)
        _backup_data = output_data
        if length(output_data) == 1
            break
        else
            output_data = filter_by_least_common_index(output_data, ii)
        end

        if length(output_data) == 0
            output_data = _backup_data
        end
    end

    output_string = ""
    for bit in output_data[1]
        output_string = output_string * string(bit)
    end
    return parse(Int, output_string, base = 2)
end

oxygen_rating = get_oxygen_rating(int_data)
co2_rating = get_co2_rating(int_data)
life_support = oxygen_rating * co2_rating

println("Day 03 Part 2 Solution: ", life_support)


text:
function maxlength(svec, lo = 1, hi = length(svec))
    if lo > hi
        return -1
    end
    ml = length(svec[lo])
    for i = lo+1:hi
        ml = max(ml, length(svec[i]))
    end
    ml
end

const CHAR0 = UInt8(0)
function code_unit0(str, pos)
    @inbounds return sizeof(str) < pos ? CHAR0 : codeunit(str, pos)
end

# inspiration http://www.cs.princeton.edu/courses/archive/spring03/cs226/lectures/radix
# radix 8 bytes at a time because 8 bytes fit into a word
function three_way_radix_qsort!(svec, lo = 1, hi = length(svec), skipbytes = 0, lb = UInt[])
    if hi <= lo
        # println("terminate")
        return svec
    elseif hi - lo == 1
        if svec[lo] > svec[hi]
            svec[lo], svec[hi] = svec[hi], svec[lo]
        end
        return svec
    elseif hi - lo < 16
        svec = sort!(svec, lo, hi, InsertionSort, Base.Forward)
        return svec
    end

    #  pick a pivot
    # choose median

    # cmppos the position of letter to compare

    # i is the lower cursor, j is the upper cursor
    # p, q is used to keep track of equal elements in 4 way partition
    i, p = lo-1, lo
    j, q = hi+1, hi

    if length(lb) == 0
        # lb[lo:hi] = FastGroupBy.load_bits.(UInt, @view(svec[lo:hi]), skipbytes)
        lb = FastGroupBy.load_bits.(UInt, svec, skipbytes)
    end

    pivotl = sort!(lb[rand(lo:hi,3)])[2]

    # after this loop the data should have been partition 4 ways
    while i < j
        while true
            i += 1
            lb[i] < pivotl || break
        end
        while true
            j -= 1
            pivotl < lb[j] || break

            if j == lo
                break;
            end
        end
        if i > j
            break;
        end
        # swap them now
        svec[i], svec[j] = svec[j], svec[i]
        lb[i], lb[j] = lb[j], lb[i]
        if lb[i] == pivotl
            svec[p], svec[i] = svec[i], svec[p]
            lb[p], lb[i] = lb[i], lb[p]
            p += 1
        end
        if lb[j] == pivotl
            svec[q], svec[j] = svec[j], svec[q]
            lb[q], lb[j] = lb[j], lb[q]
            q -= 1
        end
    end

    # if all elements are equal on that position then sort again
    if p >= q
        # println("p >= q")
        if FastGroupBy.maxlength(svec, lo, hi) > skipbytes + 8
            lb[lo:hi] = FastGroupBy.load_bits.(UInt, @view(svec[lo:hi]), skipbytes+8)
            three_way_radix_qsort!(svec, lo, hi, skipbytes + 8, lb)
            return svec
        end
    end

    # @show svec
    # println(string(p:q),pivotl,string(j:i),"lo:Hi",string(lo,hi))

    if lb[i] < pivotl
        i += 1
    end

    for k = lo:p-1
        svec[k], svec[j] = svec[j], svec[k]
        lb[k], lb[j] = lb[j], lb[k]
        j -= 1
    end

    for k = hi:-1:q+1
        svec[k], svec[i] = svec[i], svec[k]
        lb[k], lb[i] = lb[i], lb[k]
        i += 1
    end

    # recursive sort

    # println("lessthan: ", lo,j,cmppos)
    three_way_radix_qsort!(svec, lo, j, skipbytes, lb)

    if i == hi && lb[i] == pivotl
        i += 1
    end
    if FastGroupBy.maxlength(svec, j+1, i-1) > skipbytes + 8
        # println("mid: ", j+1,i-1,cmppos+1)
        lb[j+1:i-1] = FastGroupBy.load_bits.(UInt, @view(svec[j+1:i-1]), skipbytes+8)
        three_way_radix_qsort!(svec, j+1, i-1, skipbytes + 8, lb) # need to compute lb again
    end

    # println("gt: ", i,hi,cmppos+1)
    three_way_radix_qsort!(svec, i, hi, skipbytes, lb)

    return svec
end


text:
# Training a deeputoencoder on MNIST
# The method uses Nesterov's accelerated gradient, with minibatches
# (c) David Barber, University College London 2015

PlotResults=false

useproc("GPU") # GPU about 4 times faster than CPU
#useproc("CPU") # GPU about 4 times faster than CPU

using MAT

Ntrain=60000
BatchSize=500
TrainingIts=20000 # number of Nesterov updates (takes about 6 minutes on a Titan GPU)
include("loadmnist.jl")
images,label=loadmnist()
r=randperm(size(images,2))
data=images[:,r]
# bound away from 0 and 1 to avoid log(0) problems:
tol=0.000001
data[data.>(1-tol)]=1-tol
data[data.<tol]=tol

H=[784 1000 500 250 30 250 500 1000 784] # number of units in each layer

L=length(H) # number of hidden layers
# node indices:
w=Array(ADnode,L) # weights
bias=Array(ADnode,L) # biases
h=Array(ADnode,L) # hidden layer index (input layer is h[1])

StartCode()
ytrain=h[1]=ADnode()
for layer=2:L-1
    w[layer]=ADvariable()
    bias[layer]=ADvariable()
#    h[layer]=kinklinAXplusBias(w[layer],h[layer-1],bias[layer])
#    h[layer]=absAXplusBias(w[layer],h[layer-1],bias[layer])
#    h[layer]=rectlinAXplusBias(w[layer],h[layer-1],bias[layer])
    h[layer]=alex(w[layer]*h[layer-1]) # a log-linear transfer
#    h[layer]=rectlin(w[layer]*h[layer-1])
end
w[L]=ADvariable()
bias[L]=ADvariable()
h[L]=AXplusBias(w[L],h[L-1],bias[L])
ypred=sigmoid(h[L]) # just use for testing, not training
meanSqloss=meanSquareLoss(ypred,ytrain)
loss=BinaryKullbackLeiblerLossXsigmoidY(ytrain,h[L])
net=EndCode()

#instantiate root node values:
net.value[h[1]]=data[:,1:BatchSize]

for i=2:L
    net.value[w[i]]=.5*sign(randn(H[i],H[i-1]))/sqrt(H[i-1])
    net.value[bias[i]]=.005*sign(randn(H[i],1))+0.1
end
#initvalue=deepcopy(net.value[w[L]])

net=compile(net) # compile the DAG and preallocate memory

@gpu CUDArt.init([0]) # let the user do device management
@gpu net=convert(net,PROC,Float32) # use 32 bits -- much faster than 64 bit

#gradcheck(net)

println("Training: using $(net.gpu==true? "GPU" : "CPU") $(net.eltype)")
tstart=time()
error=Array(net.eltype,0)
SquareError=Array(net.eltype,0)
ParsToUpdate=Parameters(net)
velo=NesterovInit(net)
#avgrad=GradientDescentMomentumInit(net); Momentum=0.5
minibatchstart=1 # starting datapoint for the minibatch
#ForwardPassList!(net,ExcludeNodes=[ypred])
    for iter=1:TrainingIts
        LearningRate=(0.25/(1+iter/10000))
        minibatchstart,minibatch=GetBatch(minibatchstart,BatchSize,Ntrain)
        net.value[ytrain]=cArray(net.gpu,net.eltype,data[:,minibatch]) # select batch of data
        ADforward!(net)
        ADbackward!(net)
        push!(error,extract(net.value[net.FunctionNode])[1])
        push!(SquareError,784*extract(net.value[meanSqloss])[1])
        printover("iteration $iter: training loss = $(error[iter]) : meanSqLoss = $(784*to_host(net.value[meanSqloss])[1])")
        for par in ParsToUpdate
            #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
            #GradientDescentMomentumUpdate!(net.value[par],net.gradient[par],avgrad[par],Momentum,LearningRate)
            NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,iter)
        end
    end
tend=time();
println("\nTraining took $(tend-tstart) seconds\n")

if PlotResults
   net=convert(net,"CPU") # to make analysis easier
   figure(1)
   plot(error); title("training loss")
   figure(2)
    for i=1:10 # plot the reconstructions for a few datapoints
            p=imshow([reshape(net.value[h[1]][:,i],28,28)'  reshape(net.value[ypred][:,i],28,28)'],interpolation="none",cmap=ColorMap("gray"))
            display(p)
            println("press key to continue")
            readline(STDIN)
    end
end

end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
print_infix(io::IO, op, xs...) = print_infix(io, op, xs)

function print_infix(io::IO, op, xs)
    compact = get(io, :compact, false)
    if length(xs) == 1
        # infix operator must have at least 2 args
        return print_func(io, op, xs[1])
    end

    for (k, x) in enumerate(xs)
        print(io, x)
        if k != lastindex(xs)
            if compact
                print(io, op)
            else
                print(io, " ", op, " ")
            end
        end
    end
end

print_func(io::IO, f, xs...) = print_func(io, f, xs)
print_func(io::IO, ::typeof(-), x) = print(io, "-", x)
function print_func(io::IO, f, xs)
    compact = get(io, :compact, false)
    print(io, f, "(")
    for (k, each) in enumerate(xs)
        print(io, each)
        if k != lastindex(xs)
            if compact
                print(io, ",")
            else
                print(io, ", ")
            end
        end
    end
    print(io, ")")
end

function print_term(io::IO, t::Term)
    if isinfix(t.head)
        print_infix(io, t.head, t.args)
    else
        print_func(io, t.head, t.args)
    end
end

const UNIPOW = Dict()

for (k, each) in enumerate("¹²³⁴⁵⁶⁷⁸⁹")
    UNIPOW[k] = each
end

function print_term(io::IO, t::Term{typeof(^)})
    compact = get(io, :compact, false)
    x, y = t.args
    if y in 1:9
        print(io, x, UNIPOW[Int(y)])
    else
        print(io, x, "^", y)
    end
end

Base.show(io::IO, x::Numeric) = print(io, x.value)
Base.show(io::IO, x::Variable) = print(io, x.name)
Base.show(io::IO, ::Im) = print(io, "im")
Base.show(io::IO, ::Constant{sym}) where sym = print(io, sym)
Base.show(io::IO, t::Term) = print_term(io, t)


text:
function _elementary_path_lazy_callback(
    graph,
    source,
    destination,
    m,
    x,
    cb_data,
)
    # Callback to generate subtour-elimination constraints.

    # Based on the hypothesis that the strengthening constraints have been
    # added beforehand! I.e. at most one edge incoming and at most one outgoing
    # (except for source/destination).
    x = Dict(e => callback_value(cb_data, x[e]) for e in edges(graph))

    # This callback is sometimes called for noninteger nodes! Just let the
    # solver go on for these nodes (finding a loop does not really make sense).
    # Use this loop to count nonzero values in the solution, this will be
    # useful later to get out of a loop sooner.
    n_x_nonzero = 0
    for e in edges(graph)
        # Test for a quite forgiving tolerance.
        ε = 0.01
        if x[e] >= ε && x[e] <= 1.0 - ε
            # The value for x[e] is clearly not integer.
            return
        end

        # Count nonzero values.
        if x[e] >= 0.5
            n_x_nonzero += 1
        end
    end

    # Helper methods.
    inedges(g, v) = (edgetype(g)(x, v) for x in inneighbors(g, v))
    outedges(g, v) = (edgetype(g)(v, x) for x in outneighbors(g, v))

    # Find a subtour. Start from the source, use edges until the destination
    # is found.
    # This is ensured to be correct, because only one outgoing edge from the
    # current node, as per the constraints (and this value is integer).
    current_node = source # TODO: factor this loop out? Also useful for sorting edges.
    while current_node != destination
        for e in outedges(graph, current_node)
            if x[e] >= 0.5
                x[e] = 0.0
                current_node = dst(e)
                n_x_nonzero -= 1
                break
            end
        end

        # If all nonzero values have been found, stop iterating.
        if n_x_nonzero == 0
            break
        end
    end

    # Find the edges that remain after the path is removed.
    lhs_edges = Edge[]
    for e in edges(graph)
        # Don't consider edges adjacent to the source or the destination
        if src(e) == source || dst(e) == destination
            continue
        end

        if x[e] >= 0.5
            push!(lhs_edges, e)
        end
    end

    # Nothing left? No subtour!
    if length(lhs_edges) == 0
        return
    end

    # Add lazy constraints (subtour elimination). Separate all subtours to
    # generate constraints that are as tight as possible.
    while length(lhs_edges) > 0
        # Find one subtour among these edges. Same technique as above: start
        # from a node of the tour, use edges until you find the starting point.
        first_edge = lhs_edges[1]
        deleteat!(lhs_edges, 1)

        # No loop of just one edge.
        if length(lhs_edges) == 0
            break
        end

        src_node = src(first_edge) # Memorise the source for this iteration.
        cur_node = dst(lhs_edges[1]) # First node to explore.
        con_edges = Edge[first_edge] # List of edges in the loop.

        while cur_node != src_node # While not back to the source...
            for (i, e) in enumerate(lhs_edges) # Go through all edges...
                if src(e) == cur_node # If this one starts from the current node:
                    cur_node = dst(e)
                    push!(con_edges, e)
                    deleteat!(lhs_edges, i)
                    break
                end
            end
        end

        # Build the corresponding constraint.
        con = @build_constraint(
            sum(x[e] for e in con_edges) <= length(con_edges) - 1
        )
        MOI.submit(m, MOI.LazyConstraint(cb_data), con)
    end
end

function _sort_path(
    path::Vector{Tuple{Int, Int}},
    source::Int,
    destination::Int,
    n_vertices::Int,
)
    sorted = Tuple{Int, Int}[]
    current_node = source
    edges = copy(path)

    i = 0
    while length(edges) > 0
        # Find the corresponding edge (or the first one that matches).
        for i in 1:length(edges)
            e = edges[i]
            if e[1] == current_node
                push!(sorted, e)
                current_node = e[2]
                deleteat!(edges, i)
                break
            end
        end

        if current_node == destination
            break
        end

        # Safety: ensure this loop does not make too many iterations.
        i += 1
        if i > n_vertices
            error(
                "Assertion failed: infinite loop when sorting the edges of the path $path",
            )
        end
    end

    if length(edges) > 0
        error(
            "Edges remaining after sorting the edges: the solution is likely to contain at least a subtour, $edges",
        )
    end

    return sorted
end

function formulation(
    i::ElementaryPathInstance{Int, Maximise},
    f::DefaultLinearFormulation
)
    n = nv(i.graph)

    # Helper methods.
    inedges(g, v) = (edgetype(g)(x, v) for x in inneighbors(g, v))
    outedges(g, v) = (edgetype(g)(v, x) for x in outneighbors(g, v))

    # List all the vertices that are not the source nor the destination.
    other_nodes = collect(filter(vertices(i.graph)) do v
        return v != i.src && v != i.dst
    end)

    # Build the optimisation model behind solve_linear.
    model = Model(f.solver)
    set_silent(model)

    x = @variable(model, [e in edges(i.graph)], binary = true)
    # Normally, with this formulation, explicitly having binary variables
    # is not necessary, but this is only valid if minimising a linear function.
    # In all other cases, there is no guarantee to have an integer-feasible
    # solution. Moreover, adding new constraints may make this definition
    # required (especially budget constraints).

    # Set the names, due to the use of anonymous variables with JuMP.
    for e in edges(i.graph)
        set_name(x[e], "x_$(src(e))_$(dst(e))")
    end

    # For a vertex in: ∑ input flow - ∑ output flow.
    function edge_incidence(v)
        if length(inedges(i.graph, v)) > 0
            ins = sum(x[e] for e in inedges(i.graph, v))
        else
            ins = 0.0
        end

        if length(outedges(i.graph, v)) > 0
            outs = sum(x[e] for e in outedges(i.graph, v))
        else
            outs = 0.0
        end

        return ins - outs
    end

    # Ensure that the unit flow goes from the source to the destination.
    @constraint(model, edge_incidence(i.src) == -1)
    @constraint(model, edge_incidence(i.dst) == 1)

    # Flow conservation at each node (except source and destination).
    @constraint(model, [v in other_nodes], edge_incidence(v) == 0)

    # Eliminate a large number of simple subtours.
    @constraint(model, sum(x[e] for e in inedges(i.graph, i.src)) == 0)
    @constraint(model, sum(x[e] for e in outedges(i.graph, i.src)) == 1)
    @constraint(model, sum(x[e] for e in inedges(i.graph, i.dst)) == 1)
    @constraint(model, sum(x[e] for e in outedges(i.graph, i.dst)) == 0)

    for v in other_nodes
        @constraint(model, sum(x[e] for e in inedges(i.graph, v)) <= 1)
        @constraint(model, sum(x[e] for e in outedges(i.graph, v)) <= 1)
    end

    # Objective function.
    @objective(model, Max, sum(i.rewards[e] * x[e] for e in keys(i.rewards)))

    # Set the lazy-constraint callback to ensure the solution is always an elementary path.
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> _elementary_path_lazy_callback(
            i.graph,
            i.src,
            i.dst,
            model,
            x,
            cb_data,
        ),
    )

    return model, x
end

function solve(
    i::ElementaryPathInstance{Int, Maximise},
    f::DefaultLinearFormulation
)
    m, x = formulation(i, f)
    optimize!(m)
    return ElementaryPathSolution(i, _extract_lp_solution(i, x))
end


text:
using Oceananigans.TimeSteppers: ab2_step_field!

import Oceananigans.TimeSteppers: ab2_step!

combine_events(::Nothing, tracer_events) = MultiEvent(tuple(tracer_events...))
combine_events(free_surface_event, tracer_events) = MultiEvent(tuple(free_surface_event, tracer_events...))

function ab2_step!(model::HydrostaticFreeSurfaceModel, Δt, χ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    step_field_kernel! = ab2_step_field!(device(model.architecture), workgroup, worksize)

    # Launch velocity update kernels

    velocities_events = []

    for name in (:u, :v)
        Gⁿ = model.timestepper.Gⁿ[name]
        G⁻ = model.timestepper.G⁻[name]
        velocity_field = model.velocities[name]

        event = step_field_kernel!(velocity_field, Δt, χ, Gⁿ, G⁻,
                                   dependencies=Event(device(model.architecture)))

        push!(velocities_events, event)
    end

    # Launch tracer update kernels

    tracer_events = []

    for name in tracernames(model.tracers)
        Gⁿ = model.timestepper.Gⁿ[name]
        G⁻ = model.timestepper.G⁻[name]
        tracer_field = model.tracers[name]

        event = step_field_kernel!(tracer_field, Δt, χ, Gⁿ, G⁻,
                                   dependencies=Event(device(model.architecture)))

        push!(tracer_events, event)
    end

    velocities_update = MultiEvent(Tuple(velocities_events))

    # Update the free surface if not using a rigid lid once the velocities have finished updating.
    free_surface_event = ab2_step_free_surface!(model.free_surface, velocities_update, model, χ, Δt)

    tracer_and_free_surface_events = combine_events(free_surface_event, tracer_events)

    wait(device(model.architecture), tracer_and_free_surface_events)

    return nothing
end

#####
##### Free surface time-stepping: explicit, implicit, rigid lid ?
#####

ab2_step_free_surface!(free_surface::ExplicitFreeSurface, velocities_update, model, χ, Δt) =
    explicit_ab2_step_free_surface!(free_surface, velocities_update, model, χ, Δt)

function explicit_ab2_step_free_surface!(free_surface, velocities_update, model, χ, Δt)

    event = launch!(model.architecture, model.grid, :xy,
                    _ab2_step_free_surface!,
                    model.free_surface.η,
                    χ,
                    Δt,
                    model.timestepper.Gⁿ.η,
                    model.timestepper.G⁻.η,
                    dependencies=Event(device(model.architecture)))

    return event
end

@kernel function _ab2_step_free_surface!(η, χ::FT, Δt, Gηⁿ, Gη⁻) where FT
    i, j = @index(Global, NTuple)

    @inbounds begin
        η[i, j, 1] += Δt * ((FT(1.5) + χ) * Gηⁿ[i, j, 1] - (FT(0.5) + χ) * Gη⁻[i, j, 1])
    end
end

function ab2_step_free_surface!(free_surface::ImplicitFreeSurface, velocities_update, model, χ, Δt)

    ##### Implicit solver for η

    ## Need to wait for u* and v* to finish
    wait(device(model.architecture), velocities_update)
    fill_halo_regions!(model.velocities, model.architecture, model.clock, fields(model) )

    ## Leaving this here for now. There may be some scenarios where stepping forward η and then using
    ## the stepped forward value as a guess is helpful.
    ## η_save = deepcopy(free_surface.η)
    ## event = explicit_ab2_step_free_surface!(free_surface, velocities_update, model, χ, Δt)
    ## wait(device(model.architecture), event)

    ## We need vertically integrated U,V
    event = compute_vertically_integrated_volume_flux!(free_surface, model)
    wait(device(model.architecture), event)
    u=free_surface.barotropic_volume_flux.u
    v=free_surface.barotropic_volume_flux.v
    fill_halo_regions!(u.data ,u.boundary_conditions, model.architecture, model.grid, model.clock, fields(model) )
    fill_halo_regions!(v.data ,v.boundary_conditions, model.architecture, model.grid, model.clock, fields(model) )

    ## Compute volume scaled divergence of the barotropic transport and put into solver RHS
    event = compute_volume_scaled_divergence!(free_surface, model)
    wait(device(model.architecture), event)

    ## Include surface pressure term into RHS
    RHS = free_surface.implicit_step_solver.solver.settings.RHS
    RHS .= RHS/(model.free_surface.gravitational_acceleration*Δt)
    η = free_surface.η
    fill_halo_regions!(RHS   , η.boundary_conditions, model.architecture, model.grid)
    fill_halo_regions!(η.data, η.boundary_conditions, model.architecture, model.grid)
    ##  need to subtract Azᵃᵃᵃ(i, j, 1, grid)*η[i,j, 1]/(g*Δt^2)
    event = add_previous_free_surface_contribution(free_surface, model, Δt )
    wait(device(model.architecture), event)
    fill_halo_regions!(RHS   , η.boundary_conditions, model.architecture, model.grid)
    ## RHS .= RHS .+ free_surface.η.data/Δt

    ## Then we can invoke solve_for_pressure! on the right type via calculate_pressure_correction!
    x  = free_surface.implicit_step_solver.solver.settings.x
    x .= η.data
    fill_halo_regions!(x ,η.boundary_conditions, model.architecture, model.grid)
    solve_poisson_equation!(free_surface.implicit_step_solver.solver, RHS, x; Δt=Δt, g=free_surface.gravitational_acceleration)
    fill_halo_regions!(x ,η.boundary_conditions, model.architecture, model.grid)
    free_surface.η.data .= x

    ## The explicit form of this function defaults to returning an event, we do the same for now.
    return event
end


text:
"""
This module provides a set of helper function to safely infer return types of functions.

In Gridap, we rely as less as possible in type inference. But, when needed, we adopt
the following mechanism in order to compute returned types. We do not rely on
the `Base._return_type` function.

This module exports following functions:
$(EXPORTS)

"""
module Inference

using DocStringExtensions
using FillArrays

export testvalue
export testvalues
export testargs
export testargs_broadcast
export return_type
export return_type_broadcast

"""
$(TYPEDSIGNATURES)

Returns the type returned by function `f` when called with arguments
of the types in `Ts`.

The underlying implementation uses the function [`testargs`](@ref) to generate
some test values in order to call the function and determine the returned type.
This mechanism does not use `Base._return_type`. One of the advantages is
that the given function `f` is called, and thus, meaningful error messages
will be displayed if there is any error in `f`.

"""
function return_type(f::Function,Ts...)
  args = testargs(f,Ts...)
  #try
    typeof(f(args...))
  #catch e
  #  if isa(e,DomainError)
  #    s = "Function $(nameof(f)) cannot be evaluated at $args, its not in the domain.\n"
  #    s *= " Define function `testargs(::typeof{$(nameof(f))},Ts...)`\n"
  #    s *= " which sould return an argument tuple in the function domain."
  #    error(s)
  #  else
  #    throw(e)
  #  end
  #end
  # TODO try - catch block makes function very in-efficient when
  # called multiple times
end

"""
    return_type_broadcast(f::Function,Ts::DataType...) -> DataType

Like [`return_type`](@ref), but when function `f` is used in a broadcast operation.
"""
function return_type_broadcast(f::Function,Ts...)
  args = testargs_broadcast(f,Ts...)
  r = broadcast(f,args...)
  typeof(r)
end

function return_type_broadcast(f::Function,Ts::Type{<:Number}...)
  return_type(f,Ts...)
end

"""
$(SIGNATURES)
"""
function testargs_broadcast(f::Function,Ts...)
  testvalues(Ts...)
end

function testargs_broadcast(f::Function,Ts::Type{<:Number}...)
  testargs(f,Ts...)
end

"""
    testargs(f::Function,Ts::DataType...) -> Tuple

Returns a tuple with valid arguments of the types in `Ts` in order to call
function `f`. It defaults to `testvalues(Ts...)`, see the [`testvalues`](@ref)
function.
The user can overload the `testargs`
function for particular functions if the default test arguments are not in the domain
of the function and a `DomainError` is raised.

# Examples

For the following function, the default test argument (which is a zero)
is not in the domain. We can overload the `testargs` function to provide
a valid test argument.

```jldoctests
using InplaceArrays.Inference
import InplaceArrays.Inference: testargs
foo(x) = sqrt(x-1)
testargs(::typeof(foo),T::DataType) = (one(T),)
return_type(foo, Int)
# output
Float64
```

"""
testargs(f::Function,Ts...) = testvalues(Ts...)

"""
    testvalue(::Type{T}) where T

Returns an arbitrary instance of type `T`. It defaults to `zero(T)` for
non-array types and to an empty array for array types.
This function is used to compute the default test arguments in
[`testargs`](@ref).
It can be overloaded for new types `T` if `zero(T)` does not makes sense.
"""
function testvalue end

testvalue(::Type{T}) where T = zero(T)

function testvalue(::Type{T}) where T<:AbstractArray{E,N} where {E,N}
   similar(T,fill(0,N)...)
end

function testvalue(::Type{T}) where T<:Fill{E,N} where {E,N}
  Fill(zero(E),fill(0,N)...)
end

"""
    testvalues(Ts::DataType...) -> Tuple

Returns a tuple with test values for each of the types in `Ts`.
Equivalent to `map(testvalue,Ts)`.
"""
function testvalues(a,b...)
  ta = testvalue(a)
  tb = testvalues(b...)
  (ta,tb...)
end

function testvalues(a)
  ta = testvalue(a)
  (ta,)
end

end # module


text:
# for convenience
const Seq = Tuple{String,String}
const Reference = Vector{Seq}

function savefasta(filename::AbstractString, ref::Vector{Tuple{String,String}})
	io = open(filename,"w")

	for (name,seq) in ref
		println(io, '>', name)
		println(io, seq)
	end

	close(io)
end
savefasta(filename::AbstractString, ref::Tuple{String,String}) = savefasta(filename,[ref])

# # FastaIO version
# function loadfasta(ref::AbstractString)
# 	out = Reference()
# 	for x in FastaReader(ref)
# 		push!(out, (x[1],uppercase(x[2])))
# 	end
# 	out
# end

# Bio.Seq version
function loadfasta(ref::AbstractString)
	# a = collect(open(FASTAReader{DNASequence},ref))
	# Seq[(seqname(s),convert(String,sequence(s))) for s in a]
	open(FASTA.Reader,ref) do io
		Seq[(FASTA.identifier(record),string(FASTA.sequence(record))) for record in io]
	end
end

function seq2codons(reference::AbstractString)
    Symbol[Symbol(reference[i:i+2]) for i in 1:length(reference)-2]
end
seq2codons(reference::Tuple{S,T}) where {S<:AbstractString,T<:AbstractString} = (reference[1], seq2codons(reference[2]))
seq2codons(references::AbstractArray) = map(seq2codons, references)

const codon2aaDict = Dict{Symbol,Symbol}(:AAA=>:K, :AAC=>:N, :AAG=>:K, :AAT=>:N, :ACA=>:T, :ACC=>:T, :ACG=>:T, :ACT=>:T, :AGA=>:R, :AGC=>:S, :AGG=>:R, :AGT=>:S, :ATA=>:I, :ATC=>:I, :ATG=>:M, :ATT=>:I, :CAA=>:Q, :CAC=>:H, :CAG=>:Q, :CAT=>:H, :CCA=>:P, :CCC=>:P, :CCG=>:P, :CCT=>:P, :CGA=>:R, :CGC=>:R, :CGG=>:R, :CGT=>:R, :CTA=>:L, :CTC=>:L, :CTG=>:L, :CTT=>:L, :GAA=>:E, :GAC=>:D, :GAG=>:E, :GAT=>:D, :GCA=>:A, :GCC=>:A, :GCG=>:A, :GCT=>:A, :GGA=>:G, :GGC=>:G, :GGG=>:G, :GGT=>:G, :GTA=>:V, :GTC=>:V, :GTG=>:V, :GTT=>:V, :TAA=>:*, :TAC=>:Y, :TAG=>:*, :TAT=>:Y, :TCA=>:S, :TCC=>:S, :TCG=>:S, :TCT=>:S, :TGA=>:*, :TGC=>:C, :TGG=>:W, :TGT=>:C, :TTA=>:L, :TTC=>:F, :TTG=>:L, :TTT=>:F)

codon2aa(c::Symbol) = get(codon2aaDict,c,:?)
codon2aa(c::AbstractString) = codon2aa(convert(Symbol,c))

codon2aa(A::AbstractArray) = map!(codon2aa,Array{Symbol}(undef,size(A)),A)


text:
include("area.jl")
include("vectorout.jl")
include("volume.jl")

function Assembleout(IK,JK,VK,Node, Ele, Cen, Face_out, EE, NU, FF)
    @inbounds @sync @distributed for i = 1:size(Face_out,1)
        v = Volume(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                   Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                   Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                   Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
        A1 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                  Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                  Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3])
        A2 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                  Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                  Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
        A3 = area(Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                  Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                  Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
        A4 = area(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                  Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                  Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
        n = Vectorout(Node[Face_out[i,1],1],Node[Face_out[i,1],2],Node[Face_out[i,1],3],
                      Node[Face_out[i,2],1],Node[Face_out[i,2],2],Node[Face_out[i,2],3],
                      Node[Face_out[i,3],1],Node[Face_out[i,3],2],Node[Face_out[i,3],3],
                      Cen[Face_out[i,4],1],Cen[Face_out[i,4],2],Cen[Face_out[i,4],3])
        bAx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(5/12)*A2+n[3,1]*(1/12)*A3+n[4,1]*(5/12)*A4)
        bAy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(5/12)*A2+n[3,2]*(1/12)*A3+n[4,2]*(5/12)*A4)
        bAz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(5/12)*A2+n[3,3]*(1/12)*A3+n[4,3]*(5/12)*A4)
        bA = [bAx 0 0; 0 bAy 0; 0 0 bAz; bAy bAx 0; 0 bAz bAy; bAz 0 bAx]
        bBx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(5/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(1/12)*A4)
        bBy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(5/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(1/12)*A4)
        bBz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(5/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(1/12)*A4)
        bB = [bBx 0 0; 0 bBy 0; 0 0 bBz; bBy bBx 0; 0 bBz bBy; bBz 0 bBx]
        bCx = (1/v)*(n[1,1]*(1/3)*A1+n[2,1]*(1/12)*A2+n[3,1]*(5/12)*A3+n[4,1]*(5/12)*A4)
        bCy = (1/v)*(n[1,2]*(1/3)*A1+n[2,2]*(1/12)*A2+n[3,2]*(5/12)*A3+n[4,2]*(5/12)*A4)
        bCz = (1/v)*(n[1,3]*(1/3)*A1+n[2,3]*(1/12)*A2+n[3,3]*(5/12)*A3+n[4,3]*(5/12)*A4)
        bC = [bCx 0 0; 0 bCy 0; 0 0 bCz; bCy bCx 0; 0 bCz bCy; bCz 0 bCx]
        bDx = (1/v)*(n[2,1]*(1/12)*A2+n[3,1]*(1/12)*A3+n[4,1]*(1/12)*A4)
        bDy = (1/v)*(n[2,2]*(1/12)*A2+n[3,2]*(1/12)*A3+n[4,2]*(1/12)*A4)
        bDz = (1/v)*(n[2,3]*(1/12)*A2+n[3,3]*(1/12)*A3+n[4,3]*(1/12)*A4)
        bD = [bDx 0 0; 0 bDy 0; 0 0 bDz; bDy bDx 0; 0 bDz bDy; bDz 0 bDx]
        B = hcat(bA,bB,bC,bD)
        D = (EE/((1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0; NU 1-NU NU 0 0 0; NU NU 1-NU 0 0 0; 0 0 0 (1-2*NU)/2 0 0; 0 0 0 0 (1-2*NU)/2 0; 0 0 0 0 0 (1-2*NU)/2]
        ke = B'*D*B*v
        DOF = [3*Face_out[i, 1]-2,3*Face_out[i, 1]-1,3*Face_out[i, 1],
               3*Face_out[i, 2]-2,3*Face_out[i, 2]-1,3*Face_out[i, 2],
               3*Face_out[i, 3]-2,3*Face_out[i, 3]-1,3*Face_out[i, 3],
               3*Face_out[i, 5]-2,3*Face_out[i, 5]-1,3*Face_out[i, 5]]
        xi = i*144-143
        @inbounds for n1 = 1:12
            for n2 = 1:12
                IK[xi] = DOF[n1]
                JK[xi] = DOF[n2]
                VK[xi] = ke[n1,n2]
                xi += 1
            end
        end
    end
    return IK, JK, VK
end

max_stars_repo_path: src/multi-core juSFEM/assembleout.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/multi-core juSFEM/assembleout.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/multi-core juSFEM/assembleout.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
# String Decorators
# -----------------

remove_newlines(s) = replace(s, r"\r|\n" => "")

macro dna_str(seq, flag)
    if flag == "s"
        return DNASequence(remove_newlines(seq))
    elseif flag == "d"
        return quote
            DNASequence($(remove_newlines(seq)))
        end
    end
    error("Invalid DNA flag: '$(flag)'")
end

macro dna_str(seq)
    return DNASequence(remove_newlines(seq))
end

macro rna_str(seq, flag)
    if flag == "s"
        return RNASequence(remove_newlines(seq))
    elseif flag == "d"
        return quote
            RNASequence($(remove_newlines(seq)))
        end
    end
    error("Invalid RNA flag: '$(flag)'")
end

macro rna_str(seq)
    return RNASequence(remove_newlines(seq))
end

macro aa_str(seq, flag)
    if flag == "s"
        return AminoAcidSequence(remove_newlines(seq))
    elseif flag == "d"
        return quote
            AminoAcidSequence($(remove_newlines(seq)))
        end
    end
    error("Invalid Amino Acid flag: '$(flag)'")
end

macro aa_str(seq)
    return AminoAcidSequence(remove_newlines(seq))
end

macro char_str(seq, flag)
    if flag == "s"
        return CharSequence(remove_newlines(seq))
    elseif flag == "d"
        return quote
            CharSequence($(remove_newlines(seq)))
        end
    end
    error("Invalid Char flag: '$(flag)'")
end

macro char_str(seq)
    return CharSequence(remove_newlines(seq))
end


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

function configure(version_offset, min_julia_version)
    # The version of this JLL is decoupled from the upstream version.
    # Whenever we package a new upstream release, we initially map its
    # version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
    # So for example version 2.6.3 would become 200.600.300.

    name = "NetCDF"
    upstream_version = v"4.8.1"
    version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                            upstream_version.minor * 100 + version_offset.minor,
                            upstream_version.patch * 100 + version_offset.patch)

    # Collection of sources required to build NetCDF
    sources = [
        ArchiveSource("https://github.com/Unidata/netcdf-c/archive/v$(upstream_version).zip",
                      "e4e75523466de4187cff29784ff12755925f17e753bff0c9c46cd670ca63c6b2")
    ]

    # HDF5.h in /workspace/artifacts/805ccba77cd286c1afc127d1e45aae324b507973/include
    # Bash recipe for building across all platforms
    script = raw"""
cd $WORKSPACE/srcdir/netcdf-c-*

export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"
export LDFLAGS_MAKE="${LDFLAGS}"
CONFIGURE_OPTIONS=""

if [[ ${target} == *-mingw* ]]; then
    export LIBS="-lhdf5-0 -lhdf5_hl-0 -lcurl-4 -lz"
    # linking fails with: "libtool:   error: can't build x86_64-w64-mingw32 shared library unless -no-undefined is specified"
    # unless -no-undefined is added to LDFLAGS
    LDFLAGS_MAKE="${LDFLAGS} ${LIBS} -no-undefined"

    # testset fails on mingw (NetCDF 4.8.1)
    # libtool: link: cc -fno-strict-aliasing -o .libs/pathcvt.exe pathcvt.o  -L/workspace/destdir/bin ../liblib/.libs/libnetcdf.dll.a -lhdf5-0 -lhdf5_hl-0 -lcurl-4 -lz -L/workspace/destdir/lib
    # pathcvt.o:pathcvt.c:(.text+0x15c): undefined reference to `NCpathcvt'
    CONFIGURE_OPTIONS="--disable-testsets"

elif [[ "${target}" == *-apple-* ]]; then
    # this file is referenced by hdf.h by not installed
    touch ${includedir}/features.h
fi

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-utilities \
    --enable-shared \
    --disable-static \
    $CONFIGURE_OPTIONS
make LDFLAGS="${LDFLAGS_MAKE}" -j${nproc}
make install
nc-config --all
"""

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    # Set equal to the supported platforms in HDF5
    platforms = [
        Platform("x86_64", "linux"),
        # HDF5_jll on armv7l should use the same glibc as the root filesystem
        # before it can be used
        # https://github.com/JuliaPackaging/Yggdrasil/pull/1090#discussion_r432683488
        # Platform("armv7l", "linux"; libc="glibc"),
        Platform("aarch64", "linux"; libc="glibc"),
        Platform("x86_64", "macos"),
        Platform("x86_64", "windows"),
        Platform("i686", "windows"),
    ]
    if min_julia_version == v"1.6"
        push!(platforms, Platform("aarch64","macos"))
    end
    # The products that we will ensure are always built
    products = [
        LibraryProduct("libnetcdf", :libnetcdf),
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = [
        Dependency(PackageSpec(name="HDF5_jll"), compat="1.12.1"),
        Dependency("Zlib_jll"),
    ]

    jll_stdlibs = Dict(
        v"1.3" => [
            Dependency("LibCURL_jll", v"7.71.1"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll", v"1.9.0"),
            Dependency("MbedTLS_jll", v"2.16.8"),
            Dependency("nghttp2_jll", v"1.40.0"),
        ],
        v"1.6" => [
            Dependency("LibCURL_jll"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll"),
            Dependency("MbedTLS_jll", v"2.24.0"),
            Dependency("nghttp2_jll"),
        ]
    )

    append!(dependencies, jll_stdlibs[min_julia_version])

    # Build the tarballs, and possibly a `build.jl` as well.
    return name, version, sources, script, platforms, products, dependencies
end

max_stars_repo_name: jwahlstrand/Yggdrasil

max_issues_repo_name: jwahlstrand/Yggdrasil

max_forks_repo_name: jwahlstrand/Yggdrasil