text:
function train_layer_1_sigmoid_fast2(neurons,vars,obs,n_per_part,inds_all,x,y)

    # a_all=zeros(neurons,vars+1)
    a_all=Vector{Vector{Float64}}()
    ooops=0
    i_keep=Vector{Int64}()
    errs=Array{Float64}(undef,0)
    obs1=13
    for i=1:neurons # num_ranges^vars
        # if n_per_part[i+1]-n_per_part[i]!=0

            try

                itmp=randperm(obs)[1:(obs1)]
                atmp=[ones(obs1) x[itmp,:]]\isigm1.(y[itmp])
                # atmp=[[ones(obs1) x[itmp,:]];[ones(obs1) x[itmp,:]]]\[isigm1.(y[itmp]);zeros(obs1)]
                # obs1=copy(obs1)
                err1=mean(abs.(sigm1.([ones(obs1) x[itmp,:]]*atmp)-y[itmp]))

                # mat1=sigm1.([ones(obs1) x[itmp,:]]*atmp)
                # der1=mat1.*(1.0 .- mat1)
                # err1=mean(abs.(mat1-y[itmp]))+mean(abs.(der1))

                # println(err1)
                # atmp=[ones(n_per_part[i+1]-n_per_part[i]) x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:]]\isigm1.(y[inds_all[n_per_part[i]+1:n_per_part[i+1]]])
                # obs1=length(inds_all[n_per_part[i]+1:n_per_part[i+1]])

                # if sum(isnan.(atmp))==0
                #     # a_all[i,:].=atmp
                #     err=[ones(n_per_part[i+1]-n_per_part[i]) x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:]]*atmp.-isigm1.(y[inds_all[n_per_part[i]+1:n_per_part[i+1]]])
                #     # println(i," ",mean(abs.(err)))
                # else
                #     ooops+=1
                #     # println("local inverse problem, neuron=",i," total ooops=",ooops)
                # end
                # # println(a_all[i,:])

                # println(i," ",obs1)

                # if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 push!(a_all,atmp)
                # else push!(a_all,zeros(vars+1)); ooops+=1 end

                if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 && sum(abs.(atmp[2:end]))>0 # && err1<0.02
                    push!(a_all,atmp)
                    push!(errs,err1)
                    push!(i_keep,i)
                else
                    ooops+=1
                end

            catch ex
                ooops+=1
                println("problem neuron=",i," total ooops=",ooops," ex=",ex)
                # push!(a_all,zeros(length(inds_all[n_per_part[i]+1:n_per_part[i+1]])))
                push!(a_all,zeros(vars+1))
                # continue
            end
        # end
        # println("calculating weights for neuron ",i," of ",neurons)
        if i in 1:1000:neurons println(Dates.format(now(), "HH:MM:SS")," Calculating weights w for ",i," neuron, of ",neurons," Total.") end

    end
    # println("ooops=",ooops)
    # layer1=zeros(obs,neurons)
    # for i=1:neurons
    #     if i in 1:100:neurons println(Dates.format(now(), "HH:MM:SS")," Applying sigmoid on ",i," neuron, of ",neurons," Total.") end
    #     # layer1[:,i].=sigm1.([ones(obs) x]*a_all[i,:])
    #     # println(size(a_all[i]))
    #     layer1[:,i].=sigm1.([ones(obs) x]*a_all[i])
    #     # x[i_so,:]*a_all[i,:]
    # end
    # n_per_part=[n_per_part[1];n_per_part[i_keep.+1]]
    # a_all=a_all[i_keep]
    println(minimum(errs)," ",mean(errs)," ",maximum((errs)))

    println(Dates.format(now(), "HH:MM:SS")," Constructing Layer 1")
    ii1=sortperm(errs)[1:20000]
    a_all=a_all[ii1]
    layer1=sigm1.([ones(obs) x]*reduce(hcat, a_all))
    # layer1_deriv=layer1.*(1.0 .- layer1)

    # layer1=zeros(obs,length(a_all))
    # for i=1:obs
    #     for j=1:length(a_all)
    #         layer1[i,j]=(x[i,:]*a_all[j][2:end]')[1]+a_all[j][1]
    #     end
    #     println(i)
    # end

    println(Dates.format(now(), "HH:MM:SS")," Selecting Best Features")
    cors=zeros(size(layer1,2))
    for i=1:length(cors)
        # cors[i]=cor(layer1[:,i],y)
        cors[i]=mean(abs.(layer1[:,i]-y))# +mean(abs.(layer1_deriv[:,i]))
    end
    ico=sortperm(abs.(cors))[1:20000]
    # ico=cors.<0.25
    layer1=layer1[:,ico]
    a_all=a_all[ico]
    println(mean(abs.(cors[ico])))

    println(Dates.format(now(), "HH:MM:SS")," Solving system (",obs," by ",size(layer1,2)+1,") for weights v in the output Layer")
    # println("calculating linear terms...")
    # a_layer1=[sigm1.(layer1) ones(obs)]\y

    # a_layer1=[layer1 ones(obs)]\y
    # println(layer1'*layer1)
    # a_layer1 = lsmr(([layer1 ones(obs)]), (y),verbose=true)

    a_layer1=inv([layer1 ones(obs)]'*[layer1 ones(obs)])*[layer1 ones(obs)]'*y

    # println(Dates.format(now(), "HH:MM:SS")," All weights computed.")
    # predl1=[sigm1.(layer1) ones(obs)]*a_layer1

    return a_all,a_layer1,layer1

end


text:
struct ClusterCovariance{T} <: CovarianceEstimator
    clusters::T
end

cluster(x::Symbol) = ClusterCovariance((x,))
cluster(args...) = ClusterCovariance(args)

Vcov.completecases(df::AbstractDataFrame, v::ClusterCovariance) = DataFrames.completecases(df, collect(v.clusters))

function materialize(df::AbstractDataFrame, v::ClusterCovariance)
    ClusterCovariance(NamedTuple{v.clusters}(ntuple(i -> group(df[!, v.clusters[i]]), length(v.clusters))))
end

function df_FStat(::RegressionModel, v::ClusterCovariance, ::Bool)
    minimum((length(c.pool) for c in values(v.clusters))) - 1
end

function S_hat(x::RegressionModel, v::ClusterCovariance)
    # Cameron, Gelbach, & Miller (2011): section 2.3
    dim = size(modelmatrix(x), 2) * size(residuals(x), 2)
    S = zeros(dim, dim)
    G = typemax(Int)
    for c in combinations(keys(v.clusters))
        # no need for group in case of one fixed effect, since was already done in VcovMethod
        f = (length(c) == 1) ? v.clusters[c[1]] : group((v.clusters[var] for var in c)...)
        # capture length of smallest dimension of multiway clustering in G
        G = min(G, length(f.pool))
        S += (-1)^(length(c) - 1) * helper_cluster(modelmatrix(x), residuals(x), f)
    end
    # scale total vcov estimate by ((N-1)/(N-K)) * (G/(G-1))
    # another option would be to adjust each matrix given by helper_cluster by number of its categories
    # both methods are presented in Cameron, Gelbach and Miller (2011), section 2.3
    # I choose the first option following reghdfe
    rmul!(S, (size(modelmatrix(x), 1) - 1) / dof_residual(x) * G / (G - 1))
end

# res is a Vector in OLS, Matrix in IV
function helper_cluster(X::Matrix, res::Union{Vector, Matrix}, f::CategoricalVector)
    X2 = zeros(eltype(X), length(f.pool), size(X, 2) * size(res, 2))
    index = 0
    for k in 1:size(res, 2)
        for j in 1:size(X, 2)
            index += 1
            @inbounds @simd for i in 1:size(X, 1)
                X2[f.refs[i], index] += X[i, j] * res[i, k]
            end
        end
    end
    return Symmetric(X2' * X2)
end

function StatsBase.vcov(x::RegressionModel, v::ClusterCovariance)
    xtx = inv(crossmodelmatrix(x))
    pinvertible(Symmetric(xtx * S_hat(x, v) * xtx))
end

ext: jl


text:
import HoTRG_lemon.LatticeModule:
	getTrotterparameter,
	getTrotteriteration,
	getTrotterlayers,
	getBeta,
	setTrotterparameter!,
	setTrotterIteration!

abstract QuantumSimulator{T} <: Simulator{T}
abstract Quantum2dSimulator{T} <: QuantumSimulator{T}

#---
# functions for QuantumSimulator
function initializeCount!(simulator::QuantumSimulator)
	simulator.countiteration = 0
	simulator.trottercount = 0
end

function getTrotterCount(simulator::QuantumSimulator)
	return simulator.trottercount
end

function getSpaceCount(simulator::QuantumSimulator)
	return getCount(simulator) - getTrotterCount(simulator)
end

function countUp!(simulator::QuantumSimulator, whichone::AbstractString)
	countUp!(simulator)
	if whichone == "trotter"
		simulator.trottercount = simulator.trottercount + 1
	end
end

function countDown!(simulator::QuantumSimulator, whichone::AbstractString)
	countDown!(simulator)
	if whichone == "trotter"
		simulator.trottercount = simulator.trottercount - 1
	end
end

function getBeta(simulator::QuantumSimulator)
	return iteration2layer(getTrotteriteration(simulator)) * getTrotterparameter(simulator)
end

function getFirstTerm{T}(simulator::QuantumSimulator{T})
	return one(T)/getTrotterparameter(simulator)
end

#---
# inherit from QuantumSpinModel
single_args_non_preserving = [
							:getTrotterparameter,
							:getTrotteriteration,
							:getTrotterlayers
							]

for f in single_args_non_preserving
	eval(quote

	  function $(f)(simulator::QuantumSimulator)
	  return $(f)(simulator.lattice)
		end

	end)
end

set_functions = [
				:setTrotterparameter!,
				:setTrotterIteration!
				]

for f in set_functions
	eval(quote

		  function $(f){T}(simulator::QuantumSimulator{T}, parameter)
		  $(f)(simulator.lattice, parameter)
		end

	end)
end


text:
####
# Support for Fractional derivatives
#
# defined as
####

export LeftIntegral,LeftDerivative, RightDerivative, RightIntegral

@calculus_operator LeftIntegral
@calculus_operator LeftDerivative

@calculus_operator RightIntegral
@calculus_operator RightDerivative

# DLMF18.17.9 with μ=0.5 and α=β=0

function LeftIntegral(S::Jacobi,k)
    if S.b==0
        ConcreteLeftIntegral(S,k)
    else
        J=Jacobi(0.,S.a,domain(S))
        LeftIntegralWrapper(LeftIntegral(J,k)*Conversion(S,J),k)
    end
end

function LeftIntegral(S::Ultraspherical,k)
    J = Jacobi(S)
    LeftIntegralWrapper(LeftIntegral(J,k)*Conversion(S,J),0.5)
end

LeftIntegral(S::Chebyshev,k) = LeftIntegralWrapper(
    LeftIntegral(Ultraspherical(1//2,domain(S)),k)*Conversion(S,Ultraspherical(1//2,domain(S))),
    0.5)

function rangespace(Q::ConcreteLeftIntegral{Jacobi{DD,RR}}) where {DD<:IntervalOrSegment,RR}
    μ=Q.order
    S=domainspace(Q)

    JacobiWeight(S.b+μ,zero(S.b),Jacobi(S.b+μ,S.a-μ,domain(S)))
end

function RightIntegral(S::Jacobi,k)
    if S.a==0
        ConcreteRightIntegral(S,k)
    else
        J=Jacobi(S.b,zero(S.b),domain(S))
        RightIntegralWrapper(RightIntegral(J,k)*Conversion(S,J),k)
    end
end

function rangespace(Q::ConcreteRightIntegral{Jacobi{DD,RR}}) where {DD<:IntervalOrSegment,RR}
    μ=Q.order
    S=domainspace(Q)
    @assert S.a==0

    JacobiWeight(zero(S.a),S.a+μ,Jacobi(S.b-μ,S.a+μ,domain(S)))
end

for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
    @eval begin
        bandwidths(Q::$TYP{Jacobi{DD,RR},<:Any,T}) where {DD<:IntervalOrSegment,RR,T} = (0,0)
        getindex(Q::$TYP{Jacobi{DD,RR},<:Any,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,T} =
            convert(T, jacobi_frac_getindex(domain(Q),zero(T),Q.order,k,j))::T
    end
end

jacobi_frac_getindex(d::IntervalOrSegment,α,μ,k::Integer,j::Integer) =
    jacobi_frac_getindex((arclength(d)/2)^μ,α,μ,k,j)
jacobi_frac_getindex(c::Number,α,μ,k::Integer,j::Integer) =
    k==j ? c*exp(loggamma(α+k)-loggamma(α+μ+k)) : zero(promote_type(typeof(c),typeof(α),typeof(μ)))

# jacobi_frac_addentries!(d::IntervalOrSegment,α,μ,A,kr::UnitRange)=
#     jacobi_frac_addentries!((arclength(d)/2)^μ,α,μ,A,kr)
# function jacobi_frac_addentries!(c::Number,α,μ,A,kr::UnitRange)
#     γ=c*gamma(α+1)/gamma(α+1+μ)
#     for k=1:first(kr)-1
#         γ*=(α+k)/(α+μ+k)
#     end
#     for k=kr
#        A[k,k]+=γ
#        γ*=(α+k)/(α+μ+k)
#         #should be gamma(S.α+k)/gamma(S.α+μ+k)=
#         #(S.α+k-1)/(S.α+μ+k-1)*gamma(S.α+k-1)/gamma(S.α+μ+k-1)
#     end
#     A
# end

function LeftIntegral(S::JacobiWeight{Jacobi{DD,RR}},k) where {DD,RR}
    J=S.space
    @assert S.α ≈ 0
    @assert S.β ≈ J.b
    ConcreteLeftIntegral(S,k)
end

function RightIntegral(S::JacobiWeight{Jacobi{DD,RR}},k) where {DD,RR}
    J=S.space
    @assert S.α ≈ J.a
    @assert S.β ≈ 0
    ConcreteRightIntegral(S,k)
end

function LeftIntegral(S::JacobiWeight{Chebyshev{DD,RR}},k) where {DD,RR}
    # convert to Jacobi
    Q=LeftIntegral(JacobiWeight(S.β,S.α,Jacobi(S.space)),k)
    LeftIntegralWrapper(Q*Conversion(S,domainspace(Q)),k)
end

function RightIntegral(S::JacobiWeight{Chebyshev{DD,RR}},k) where {DD,RR}
    # convert to Jacobi
    Q=RightIntegral(JacobiWeight(S.β,S.α,Jacobi(S.space)),k)
    RightIntegralWrapper(Q*Conversion(S,domainspace(Q)),k)
end

for (TYP,WRAP) in ((:LeftIntegral,:LeftIntegralWrapper),
                    (:RightIntegral,:RightIntegralWrapper))
    @eval function $TYP(S::JacobiWeight{PS},k) where PS<:PolynomialSpace
        JS = JacobiWeight(S.β,S.α,Jacobi(S.space))
        $WRAP($TYP(JS,k)*Conversion(S,JS),k)
    end
end

#DLMF18.17.9
function rangespace(Q::ConcreteLeftIntegral{JacobiWeight{Jacobi{DD,RR},DD,RR,TT}}) where {DD<:IntervalOrSegment,RR,TT}
    μ=Q.order
    S=domainspace(Q)
    J=S.space
    if isapprox(S.β,-μ)
        Jacobi(J.b+μ,J.a-μ,domain(J))
    else
        JacobiWeight(S.β+μ,zero(S.β),Jacobi(J.b+μ,J.a-μ,domain(J)))
    end
end

function rangespace(Q::ConcreteRightIntegral{JacobiWeight{Jacobi{DD,RR},DD,RR,TT}}) where {DD<:IntervalOrSegment,RR,TT}
    μ=Q.order
    S=domainspace(Q)
    J=S.space
    @assert S.α==J.a
    @assert S.β==0
    if isapprox(S.β,-μ)
        Jacobi(J.b-μ,J.a+μ,domain(J))
    else
        JacobiWeight(zero(S.α),S.α+μ,Jacobi(J.b-μ,J.a+μ,domain(J)))
    end
end

for TYP in (:ConcreteLeftIntegral,:ConcreteRightIntegral)
    @eval bandwidths(Q::$TYP{JacobiWeight{Jacobi{DD,RR},DD,RR,TT}}) where {DD<:IntervalOrSegment,RR,TT} = (0,0)
end

getindex(Q::ConcreteLeftIntegral{JacobiWeight{Jacobi{DD,RR},DD,RR,TT},<:Any,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,TT,T} =
    convert(T,jacobi_frac_getindex(domain(Q),domainspace(Q).β,Q.order,k,j))::T

getindex(Q::ConcreteRightIntegral{JacobiWeight{Jacobi{DD,RR},DD,RR,TT},<:Any,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,TT,T} =
    convert(T,jacobi_frac_getindex(domain(Q),domainspace(Q).α,Q.order,k,j))::T

function choosedomainspace(Q::LeftIntegral{UnsetSpace,T},sp::JacobiWeight) where T
    @assert sp.β>0 && isapproxinteger(sp.α)

    if isapprox(Q.order,sp.β)
        Legendre(domain(sp))
    else
        JacobiWeight(sp.β-Q.order,0.,Jacobi(sp.β-Q.order,0.,domain(sp)))
    end
end

choosedomainspace(Q::LeftIntegral{UnsetSpace,T},sp::PolynomialSpace) where {T}=
                JacobiWeight(-Q.order,zero(T),Jacobi(-Q.order,-Q.order,domain(sp)))

function choosedomainspace(Q::RightIntegral{UnsetSpace,T},sp::JacobiWeight) where T
    @assert sp.α>0  && isapproxinteger(sp.β)

    if isapprox(Q.order,sp.α)
        Legendre(domain(sp))
    else
        JacobiWeight(zero(T),sp.α-Q.order,Jacobi(0.,sp.α-Q.order))
    end
end
choosedomainspace(Q::RightIntegral{UnsetSpace,T},sp::PolynomialSpace) where {T}=
    JacobiWeight(zero(T),-Q.order,Jacobi(-Q.order,-Q.order,domain(sp)))

# Define Left/RightDerivative
for (DTYP,QTYP,DWRAP,QWRAP) in ((:LeftDerivative,:LeftIntegral,:LeftDerivativeWrapper,:LeftIntegralWrapper),
                            (:RightDerivative,:RightIntegral,:RightDerivativeWrapper,:RightIntegralWrapper))
    @eval begin
        function $DTYP(S::Space,k::Real)
            i=ceil(Int,k)
            r=i-k
            $DWRAP(i<0 ? $QTYP(S,-k) : Derivative(i)*$QTYP(S,r),k)
        end
        $QTYP(S::SumSpace,k) = $QWRAP(InterlaceOperator(Diagonal([map(s->$QTYP(s,k),S.spaces)...]),SumSpace),k)
    end
end


text:
# SPDX-License-Identifier: X11
# 2020-08-29 (Live)
# Don't be late (100pt)

D, S, T = parse.(Int, split(readline(), " "))
println(S * T < D ? "No" : "Yes")

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# Now, we need an abstract type representing a numerical discretization scheme
# for the ODE
"""
Represents a map that given (t_n,u_n) returns (t_n+1,u_n+1) and cache for the
corresponding `ODEOperator` and `NonlinearOperator`
"""
abstract type ODESolver <: GridapType end

function solve_step!(
  uF::AbstractVector,solver::ODESolver,op::ODEOperator,u0::AbstractVector,t0::Real,cache) # -> (uF,tF,cache)
  @abstractmethod
end

# Default API

function solve_step!(
  uF::AbstractVector,solver::ODESolver,op::ODEOperator,u0::AbstractVector,t0::Real) # -> (uF,tF,cache)
  solve_step!(uF,solver,op,u0,t0,nothing)
end

function solve(
  solver::ODESolver,op::ODEOperator,u0::AbstractVector,t0::Real,tf::Real)
  GenericODESolution(solver,op,u0,t0,tf)
end

# testers

function test_ode_solver(solver::ODESolver,op::ODEOperator,u0,t0,tf)
  solution = solve(solver,op,u0,t0,tf)
  test_ode_solution(solution)
end

# Specialization

include("ForwardEuler.jl")

include("ThetaMethod.jl")

include("AffineThetaMethod.jl")

include("RungeKutta.jl")


text:
export run_acdcpf

""
function run_acdcpf(file::String, model_constructor, solver; kwargs...)
    data = PowerModels.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcpf(data::Dict{String,Any}, model_constructor, solver; kwargs...)
end

""
function run_acdcpf(data::Dict{String,Any}, model_constructor, solver; kwargs...)
    pm = PowerModels.build_model(data, model_constructor, post_acdcpf; kwargs...)
    #display(pm)
    return PowerModels.optimize_model!(pm, solver; solution_builder = get_solution_acdc)
end

""
function post_acdcpf(pm::AbstractPowerModel)
    add_ref_dcgrid!(pm)
    PowerModels.variable_voltage(pm, bounded = false)
    PowerModels.variable_generation(pm, bounded = false)
    PowerModels.variable_branch_flow(pm, bounded = false)

    # dirty, should be improved in the future TODO
    if typeof(pm) <: PowerModels.SOCBFPowerModel
        PowerModels.variable_branch_current(pm, bounded = false)
    end

    variable_active_dcbranch_flow(pm, bounded = false)
    variable_dcbranch_current(pm, bounded = false)
    variable_dc_converter(pm, bounded = false)
    variable_dcgrid_voltage_magnitude(pm, bounded = false)

    PowerModels.constraint_model_voltage(pm)
    constraint_voltage_dc(pm)

    for (i,bus) in ref(pm, :ref_buses)
        @assert bus["bus_type"] == 3
        PowerModels.constraint_theta_ref(pm, i)
        PowerModels.constraint_voltage_magnitude_setpoint(pm, i)
    end

    for (i, bus) in ref(pm, :bus)# PowerModels.ids(pm, :bus)
        constraint_kcl_shunt(pm, i)
        # PV Bus Constraints
        if length(ref(pm, :bus_gens, i)) > 0 && !(i in ids(pm,:ref_buses))
            # this assumes inactive generators are filtered out of bus_gens
            @assert bus["bus_type"] == 2
            PowerModels.constraint_voltage_magnitude_setpoint(pm, i)
            for j in ref(pm, :bus_gens, i)
                PowerModels.constraint_active_gen_setpoint(pm, j)
            end
        end
    end

    for i in PowerModels.ids(pm, :branch)
        # dirty, should be improved in the future TODO
        if typeof(pm) <: PowerModels.SOCBFPowerModel
            PowerModels.constraint_flow_losses(pm, i)
            PowerModels.constraint_voltage_magnitude_difference(pm, i)
            PowerModels.constraint_branch_current(pm, i)
        else
            PowerModels.constraint_ohms_yt_from(pm, i)
            PowerModels.constraint_ohms_yt_to(pm, i)
        end
    end
    for i in PowerModels.ids(pm, :busdc)
        constraint_kcl_shunt_dcgrid(pm, i)
    end
    for i in PowerModels.ids(pm, :branchdc)
        constraint_ohms_dc_branch(pm, i)
    end
    for (c, conv) in PowerModels.ref(pm, :convdc)
        #constraint_converter_filter_transformer_reactor(pm, c)
        constraint_conv_transformer(pm, c)
        constraint_conv_reactor(pm, c)
        constraint_conv_filter(pm, c)
        if conv["type_dc"] == 2
            constraint_dc_voltage_magnitude_setpoint(pm, c)
            constraint_reactive_conv_setpoint(pm, c)
        else
            if conv["type_ac"] == 2
                constraint_active_conv_setpoint(pm, c)
            else
                constraint_active_conv_setpoint(pm, c)
                constraint_reactive_conv_setpoint(pm, c)
            end
        end
        constraint_converter_losses(pm, c)
        constraint_converter_current(pm, c)
    end
end


text:
using Profile, ChromeProfileFormat

function foo()
    A = randn(1000, 1000)

    return inv(A)
end

foo()
Profile.clear()
@profile foo()

ChromeProfileFormat.save_cpuprofile("test.cpuprofile", from_c=true)


text:
function calculate_prob(n::Int, k::Int)
    p = [1.0]

    for x in 2:k
        new_p = [(1 + (1 - x)/n) * p[1]]
        for y in 1:Int(floor((x - 1)/2))
            push!(new_p, (x + 1 - (2 * y))/n * p[y] + (1 + (1 - x + y)/n) * p[y + 1])
        end

        if x % 2 == 0
            push!(new_p, 1/n * p[end])
        end

        p = new_p
    end

    println(1 - sum(p))
end

function main()
    chips = 1000000
    defects = 20000
    calculate_prob(chips, defects)
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
module DAMMmodel

using DocStringExtensions
using Statistics
using LsqFit
using DataFrames
using SparseArrays
using UnicodeFun
using GLMakie
using Unitful: R, L, mol, K, kJ, °C, m, g, cm, hr, mg, s, μmol
using UnitfulMoles: molC
using Unitful, UnitfulMoles
@compound CO₂

include(joinpath("constructors", "constants.jl"))
include(joinpath("constructors", "constructors.jl"))
include(joinpath("functions", "maths", "DAMM.jl"))
include(joinpath("functions", "maths", "DAMMfit.jl"))
include(joinpath("functions", "maths", "DAMMmat.jl"))
include(joinpath("functions", "viz", "DAMMplot.jl"))
include(joinpath("functions", "viz", "DAMMviz.jl"))
include(joinpath("functions", "maths", "qbins.jl"))
export DAMM, DAMMfit, DAMMmat, DAMMplot, DAMMviz, qbins, sDAMMmat, sDAMMmatq

function __init__()
    Unitful.register(DAMMmodel)
end

end


text:
using FunctionalCollections: PersistentVector, push, assoc

struct UnfoldType end

"""
    gen_fn = Unfold(kernel::GenerativeFunction)

Return a new generative function that applies the kernel in sequence, passing the return value of one application as an input to the next.

The kernel accepts the following arguments:

- The first argument is the `Int` index indicating the position in the sequence (starting from 1).

- The second argument is the *state*.

- The kernel may have additional arguments after the state.

The return type of the kernel must be the same type as the state.

The returned generative function accepts the following arguments:

- The number of times (N) to apply the kernel.

- The initial state.

- The rest of the arguments (not including the state) that will be passed to each kernel application.

The return type of the returned generative function is `FunctionalCollections.PersistentVector{T}` where `T` is the return type of the kernel.

If `kernel` has optional trailing arguments, the corresponding arguments can be omitted from calls to `Unfold(kernel)`.
"""
struct Unfold{T,U} <: GenerativeFunction{PersistentVector{T},VectorTrace{UnfoldType,T,U}}
    kernel::GenerativeFunction{T,U}
end

export Unfold

function has_argument_grads(gen_fn::Unfold)
    (false, has_argument_grads(gen_fn.kernel)[2:end]...)
end

# TODO
accepts_output_grad(gen_fn::Unfold) = false

function (gen_fn::Unfold)(args...)
    (_, _, retval) = propose(gen_fn, args)
    retval
end

function unpack_args(args::Tuple)
    len = args[1]
    init_state = args[2]
    params = args[3:end]
    (len, init_state, params)
end

function check_length(len::Int)
    if len < 0
        error("unfold got length of $len < 0")
    end
end

include("simulate.jl")
include("generate.jl")
include("propose.jl")
include("assess.jl")
include("generic_update.jl")
include("update.jl")
include("regenerate.jl")
include("backprop.jl")


text:
# Invertible network layer from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export CouplingLayerIRIM, CouplingLayerIRIM3D

"""
    IL = CouplingLayerIRIM(C::Conv1x1, RB::ResidualBlock)

or

    IL = CouplingLayerIRIM(n_in, n_hidden; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, logdet=false, ndims=2) (2D)

    IL = CouplingLayerIRIM(n_in, n_hidden; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, logdet=false, ndims=3) (3D)

    IL = CouplingLayerIRIM3D(n_in, n_hidden; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, logdet=false) (3D)

 Create an i-RIM invertible coupling layer based on 1x1 convolutions and a residual block.

 *Input*:

 - `C::Conv1x1`: 1x1 convolution layer

 - `RB::ResidualBlock`: residual block layer consisting of 3 convolutional layers with ReLU activations.

 or

 - `nx`, `ny`, `nz`: spatial dimensions of input

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 *Output*:

 - `IL`: Invertible i-RIM coupling layer.

 *Usage:*

 - Forward mode: `Y = IL.forward(X)`

 - Inverse mode: `X = IL.inverse(Y)`

 - Backward mode: `ΔX, X = IL.backward(ΔY, Y)`

 *Trainable parameters:*

 - None in `IL` itself

 - Trainable parameters in residual block `IL.RB` and 1x1 convolution layer `IL.C`

 See also: [`Conv1x1`](@ref), [`ResidualBlock!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct CouplingLayerIRIM <: NeuralNetLayer
    C::Conv1x1
    RB::Union{ResidualBlock, FluxBlock}
end

@Flux.functor CouplingLayerIRIM

# 2D Constructor from input dimensions
function CouplingLayerIRIM(n_in::Int64, n_hidden::Int64;
                           k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, ndims=2)

    # 1x1 Convolution and residual block for invertible layer
    C = Conv1x1(n_in)
    RB = ResidualBlock(n_in÷2, n_hidden; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, ndims=ndims)

    return CouplingLayerIRIM(C, RB)
end

CouplingLayerIRIM3D(args...;kw...) = CouplingLayerIRIM(args...; kw..., ndims=3)

# 2D Forward pass: Input X, Output Y
function forward(X::AbstractArray{Float32, N}, L::CouplingLayerIRIM) where N

    # Get dimensions
    k = Int(L.C.k/2)

    X_ = L.C.forward(X)
    X1_, X2_ = tensor_split(X_)

    Y1_ = X1_
    Y2_ = X2_ + L.RB.forward(Y1_)

    Y_ = tensor_cat(Y1_, Y2_)
    Y = L.C.inverse(Y_)

    return Y
end

# 2D Inverse pass: Input Y, Output X
function inverse(Y::AbstractArray{Float32, N}, L::CouplingLayerIRIM; save=false) where N

    # Get dimensions
    k = Int(L.C.k/2)

    Y_ = L.C.forward(Y)
    Y1_, Y2_ = tensor_split(Y_)

    X1_ = Y1_
    X2_ = Y2_ - L.RB.forward(Y1_)

    X_ = tensor_cat(X1_, X2_)
    X = L.C.inverse(X_)

    if save == false
        return X
    else
        return X, X_, Y1_
    end
end

# 2D Backward pass: Input (ΔY, Y), Output (ΔX, X)
function backward(ΔY::AbstractArray{Float32, N}, Y::AbstractArray{Float32, N}, L::CouplingLayerIRIM; set_grad::Bool=true) where N

    # Recompute forward state
    k = Int(L.C.k/2)
    X, X_, Y1_ = inverse(Y, L; save=true)

    # Backpropagate residual
    if set_grad
        ΔY_ = L.C.forward((ΔY, Y))[1]
    else
        ΔY_, Δθ_C1 = L.C.forward((ΔY, Y); set_grad=set_grad)[1:2]
    end
    ΔYl_, ΔYr_ = tensor_split(ΔY_)
    if set_grad
        ΔY1_ = L.RB.backward(ΔYr_, Y1_) + ΔYl_
    else
        ΔY1_, Δθ_RB = L.RB.backward(ΔYr_, Y1_; set_grad=set_grad)
        ΔY1_ = ΔY1_ + ΔYl_
    end

    ΔX_ = tensor_cat(ΔY1_, ΔYr_)
    if set_grad
        ΔX = L.C.inverse((ΔX_, X_))[1]
    else
        ΔX, Δθ_C2 = L.C.inverse((ΔX_, X_); set_grad=set_grad)[1:2]
    end

    set_grad ? (return ΔX, X) : (return ΔX, cat(Δθ_C1+Δθ_C2, Δθ_RB; dims=1), X)
end

## Jacobian utilities

# 2D
function jacobian(ΔX::AbstractArray{Float32, N}, Δθ::Array{Parameter, 1}, X::AbstractArray{Float32, N}, L::CouplingLayerIRIM) where N

    # Get dimensions
    k = Int(L.C.k/2)

    ΔX_, X_ = L.C.jacobian(ΔX, Δθ[1:3], X)
    X1_, X2_ = tensor_split(X_)
    ΔX1_, ΔX2_ = tensor_split(ΔX_)

    ΔY1_, Y1__ = L.RB.jacobian(ΔX1_, Δθ[4:end], X1_)
    Y2_ = X2_ + Y1__
    ΔY2_ = ΔX2_ + ΔY1_

    Y_ = tensor_cat(X1_, Y2_)
    ΔY_ = tensor_cat(ΔX1_, ΔY2_)
    ΔY, Y = L.C.jacobianInverse(ΔY_, Δθ[1:3], Y_)

    return ΔY, Y

end

# 2D/3D
function adjointJacobian(ΔY::AbstractArray{Float32, N}, Y::AbstractArray{Float32, N}, L::CouplingLayerIRIM) where N
    return backward(ΔY, Y, L; set_grad=false)
end

## Other utils

# Clear gradients
function clear_grad!(L::CouplingLayerIRIM)
    clear_grad!(L.C)
    clear_grad!(L.RB)
end

# Get parameters
function get_params(L::CouplingLayerIRIM)
    p1 = get_params(L.C)
    p2 = get_params(L.RB)
    return cat(p1, p2; dims=1)
end


text:
using Test

function testapproxequal(a::Vector{Float64}, b::Vector{Float64}, tol::Float64,
  verbose::Bool)
  v::Float64 = maximum(abs.(a-b))
  verbose && (println(a) ; println(b))
  verbose && println("$v < $tol ?")
  @test v < tol
end

function testapproxequal(a::Float64, b::Float64, tol::Float64, verbose::Bool)
  v::Float64 = abs(a-b)
  verbose && println("$v < $tol ?")
  @test v < tol
end

println("Running tests: ", Threads.nthreads(), " threads")

include("finiteFK_test.jl")

@testset "LGModel tests" begin
  @time include("lgModel_test.jl")
end

include("mvlgModel_test.jl")

@testset "SMC Sampler tests" begin
  @time include("smcSampler_test.jl")
end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ bb7a1f18-84c1-4c0a-807d-dc5575a7dd0c
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ c879bad0-7944-11eb-3cc5-6d068e87585a
md"# Arrays"

# ╔═╡ 0a8c911e-9dd2-4f7e-868b-1cd845135c53
md"## Definition"

# ╔═╡ 34558e07-9821-49a5-957b-7ef713fa9a6b
md"An array is an ordered collection of elements. The [Wikipedia webpage](https://en.wikipedia.org/wiki/Array_data_type) tells us that it relates to the mathematical concepts of vectors (1d array), matrices (2d arrays) and tensors (multi-dimensional arrays)."

# ╔═╡ cda72a9e-6ec1-4746-a6f8-e0926b5a50d9
md"## Creating an array"

# ╔═╡ daf15091-45a3-4419-beb3-643d324a1f70
md"### Empty array"

# ╔═╡ 92370118-f98d-453b-99d2-0bf053204278
md"We can create an empty array using the following commands:"

# ╔═╡ 45fcb20b-e604-458c-80d2-d96324d8a63b
a = []

# ╔═╡ 8933fc3c-ccfe-4fa8-98af-270fc8e15280
md"`a` is an empty array, and it is of type `Any` which means it can store any type of values. This is the most general declaration of an array here, because we don't assume the type of the elements that will be stored in the array. It is very practical, but adds some limitations in terms of performance."

# ╔═╡ d3e2cdd5-2dd8-45e2-8f3d-f15d39c18d1b
md"The prefered way when we now the type of the elements we will store is to declare it as follows:"

# ╔═╡ 49090226-6636-4464-b496-1e0b83e5655f
a2 = Array{Float64,1}()

# ╔═╡ 792b67f0-f1f1-4b4a-a720-aecffab167c6
md"Now if we ask Julia what is the type of `a2`:"

# ╔═╡ 7a150c77-11b6-4525-836f-1e371ba91e70
typeof(a2)

# ╔═╡ 94ad54e0-b6a0-4872-94e4-532cfd237acf
md"It returns: $(typeof(a2)). We can also use directly `eltype()` to get the type of the elements in an `Array`:"

# ╔═╡ 9b364eae-7060-4034-9924-1490fd88d275
eltype(a2)

# ╔═╡ 70cb6051-c849-4d14-84b6-3d7b48b42ebc
md"There is also a short hand to `Array{Float64,1}()`:"

# ╔═╡ 6743a46c-4f24-4eca-b4a3-18de4dc9618a
a_float = Float64[]

# ╔═╡ 34f5b3dc-0331-4e32-aaf2-7807fb1855e9
md"And it works with other types, *e.g.*:"

# ╔═╡ 72fd91e2-e83d-4173-8c29-c3f7e3ca7fb2
a3 = Int64[]

# ╔═╡ 189853b7-212b-41ec-95ba-be4993b842ac
eltype(a3)

# ╔═╡ f881789b-9cae-4de3-86e3-624f732aed82
md"### With values"

# ╔═╡ 7e554061-1ed9-4cdc-aa02-f6d7565106b7
md"We can directly instantiate an `Array` by providing values:"

# ╔═╡ 7f9e7e08-6021-4406-a09c-419f4c64b30d
[1, 2]

# ╔═╡ d13f18d8-fa64-4d80-97c7-dd392358f579
md"Or using the function directly:"

# ╔═╡ 78959d1d-b522-4567-8353-73fac5e9d120
Int64[1.0,2.0]

# ╔═╡ c45e1610-c944-40e5-b7c6-c48c8bd1c626
md"Julia tries to promote the values by default to specialise the `Array`, meaning making an `Array` of one specific type:"

# ╔═╡ 0fde7f3d-cde2-4381-89b7-4234d76f3855
[1.0, 1]

# ╔═╡ 9b5093ba-e389-4a7b-a152-40311dcea1e6
md"Here Julia promoted the `Float64` and `Int64` given as inputs into `Float64`. But if it can't, it will make an `Array` with element type `Any`:"

# ╔═╡ 48f5b4e6-83ad-4886-ac57-7aba4e42c7ca
[1.0, 1, "1"]

# ╔═╡ c95603b0-3ddd-4f9a-b59d-12af2b2f492b
md"There is a short-hand to make `Arrays` filled with zeros or ones of a certain type:"

# ╔═╡ 7e2fc580-8ea9-4409-b3e2-2e94a550a81a
zeros(Int64,5)

# ╔═╡ 088555de-40a0-4e84-8c64-8b4c8121dcca
zeros(Float64,5)

# ╔═╡ 2b9c0567-3cf7-47e8-a595-ab4f2f0fcb18
md"Or using any sequence:"

# ╔═╡ 57d80a18-57cb-424c-afca-dba5924bb8f2
collect(1:10)

# ╔═╡ 5e4a882d-608d-47d4-b594-1b33857646ac
md"### Multiple dimensions"

# ╔═╡ 2337fdaf-40d9-4d21-bb19-ac98defef41a
md"""We only made one-dimensional arrays -also called Vectors- until here, but we can make arrays of any dimensions. Let's make an Array of increased dimension:

- one dimension (*i.e.* Vector)
"""

# ╔═╡ 5009bbdb-79fe-49f0-b642-0a866e202136
zeros(Int64, 5)

# ╔═╡ ad0e06ce-7670-4205-94c3-81aa576fa30f
md"- two dimensions (*i.e.* Matrix)"

# ╔═╡ f3aaca99-af7b-468d-b2eb-571a2c664903
zeros(Int64,5,2)

# ╔═╡ fecb5933-8410-4622-98ea-862c8fb8ab59
md"""Note that a matrix can be declared using `[` too using a space to make a new column, and a new line (or `;`) for a new row:

- using a new line:
"""

# ╔═╡ d9f43a57-e061-4f4a-9bd3-c4dcdac4fb25
[1 2
 3 4
 5 6]

# ╔═╡ d3cc12a3-8251-444c-93e9-c5c09527d034
md"- using `;`:"

# ╔═╡ a152378c-db20-45f2-af28-6546b577d9fe
[1 2; 3 4; 5 6]

# ╔═╡ 0abc7154-9a2c-4547-9ebd-02f80d2b8ecf
md"- three dimensions (*i.e.* Tensor)"

# ╔═╡ 6ebef395-d692-448f-ae12-491c385f8dd2
zeros(Int64,5,2,2)

# ╔═╡ a04c4b0c-365c-4032-8884-bd2b2df040bf
md"etc..."

# ╔═╡ 5778335e-0b31-423c-9ca1-8a8b9add3680
md"## Working with arrays"

# ╔═╡ 50fedd78-95c0-4c24-86e0-d6fe7ab99a65
md"Let's make a vector and a matrix first:"

# ╔═╡ 1baf26a7-8af0-4473-86d9-497ab2546051
vec = [1,3,5,7]

# ╔═╡ 2b4dbfbb-173d-43fb-8e3f-0827536814e3
mat = [1 2
	   3 4
	   5 6
	   7 9]

# ╔═╡ 8b7fa713-9b89-46dc-b4fa-da2e348fd18e
md"### Indexing"

# ╔═╡ 8f74faf5-829e-4c08-8463-5695e574af82
md"""Indexing is made using `[`. For example we would type `vec[index]` to get the value in `vec` stored at `index`:"""

# ╔═╡ 75a44b16-5279-4d08-ba6a-dce7ea21db15
vec[2]

# ╔═╡ 8fef9d3e-b549-4a80-8314-d0a19cf6f1c1
md"Similarly, `mat[index_1, index_2]` will return the value in `mat` stored at `index_1` in the first dimension (the rows) and `index_2` in the second dimension (the columns):"

# ╔═╡ ce4bafc7-4ce2-4349-a0ec-d350ef4853f2
mat[3,2]

# ╔═╡ ae116904-1ced-43bc-a007-997123a32ffd
md"""There are some usefull keywords to index. For example use `:` to get all values in a dimension.

- All rows and first column only:
"""

# ╔═╡ 5d3771f0-0b75-4ea2-8fd5-1047e1a35779
mat[:,1]

# ╔═╡ d642b188-9e7c-493f-a9e4-711cf03227e8
md"- All columns and first row only:"

# ╔═╡ a95aaf6c-ff91-4c3c-a83a-d71afdb94878
mat[1,:]

# ╔═╡ 91e738e6-9238-4a48-ae58-0bea688e6c90
md"or to get a range of values:"

# ╔═╡ 5e2387f1-5dd7-4b49-91e4-96d057ed31fa
mat[1:3,:]

# ╔═╡ 60efc230-9ada-4de9-896a-c5d37d7fa9e0
md"or `end` to get the last value:"

# ╔═╡ 4d38a50e-57ff-461d-8af9-9f1f188521e4
mat[1,end]

# ╔═╡ b4a2e4b8-f16f-4fbd-baba-687f7468ca7f
md"We can also index using booleans:"

# ╔═╡ a939313b-1c86-492d-9993-545066bd8cd4
values_we_want = [false,true,false,false]

# ╔═╡ fca126bf-b5d6-47b5-bbe9-994f7acad0fe
vec[values_we_want]

# ╔═╡ 862f073e-8006-4671-a8e0-a6cf2024ae68
md"This is a very convenient feature because it allows us to index with more complex approaches, *e.g.*:"

# ╔═╡ f630778e-1887-4fe3-9f35-c4e286e45cba
vec[vec .> 3]

# ╔═╡ c58f05ad-dd0b-482f-9891-514b960b8bb4
md"We will see those kind of operations in a different notebook."

# ╔═╡ 67208966-c604-4889-b06e-f385647990d1
md"### Add or remove elements"

# ╔═╡ c764ac39-8d2a-4fa1-b3f3-96a24ba5527d
md"#### At the end of the array:"

# ╔═╡ 9fc4cb37-69fd-446c-95fe-30abb7559601
md"We can add elements to the end of an array using `push!`:"

# ╔═╡ 47e8bdbc-ad36-447c-8f79-f78626814be1
begin
	vectmp1 = copy(vec)
	push!(vectmp1,10)
	vectmp1
end

# ╔═╡ b0d77b3e-9674-449a-be5f-a5790a4d4089
md"""

!!! note

	We would simply use `push!(vec,10)` outside of a Pluto notebook (*e.g.* in a Julia script). But Pluto doesn't like when we mutate a variable (meaning changing its value) at several different places in the notebook. So we apply a trick here, were we put the code in a `begin ... end` statement to make it a single statement, and apply `push!` on a copy of the vector (here `vectmp1`). If it is too complicated, no worries, just consider this big blob of code just as if it was only written `push!(vec,10)`.
"""

# ╔═╡ d94f73bb-0e3f-42ea-8e8f-7e3e982836fa
md"And remove the last element using `pop!`:"

# ╔═╡ 23a8f2fa-9153-4eec-a876-db97a999b943
begin
	vectmp2 = copy(vec)
	pop!(vectmp2)
	vectmp2
end

# ╔═╡ f530fe3b-1b0d-472a-b330-d12a46e68c22
md"""
!!! note
	Note the `!` at the end of the name of the functions. Functions ending with a `!` in Julia mutate their first argument, meaning they modify it "in-place". Those functions generally have their non-mutating counterpart without the `!` that computes the result and return it, but don't modify the object given as the first argument. But we'll get into that in a further notebook.
"""

# ╔═╡ 8fbfb20a-5f5e-422e-9b59-90513050b3b8
md"#### At the beginning of the array:"

# ╔═╡ 60608c63-2b59-474e-92a2-3800e6254792
md"We can add or remove elements at the beginning of the array using `pushfirst!` and `popfirst!`:"

# ╔═╡ f153b213-1d4a-4481-ba83-29b1f8131963
begin
	vectmp3 = copy(vec)
	pushfirst!(vectmp3,0)
end

# ╔═╡ a99f7d96-654a-46b0-91a2-9120ab30af44
begin
	vectmp4 = copy(vec)
	popfirst!(vectmp4)
	vectmp4
end

# ╔═╡ 1cb1ec65-48d9-47f1-8a3d-f3767fa03ad2
md"""
!!! note

We are still applying the `begin .. end` trick here because we're inside of a Pluto notebook (and it's the case of many following functions that only exist in a mutating form). Just use `pushfirst!(vec,0)` and `popfirst!(vec)` outside of Pluto.
"""

# ╔═╡ 33c34e8f-c5b6-43a8-a329-12fc55949a2e
md"#### Anywhere in the array:"

# ╔═╡ fad8c97c-7cd6-4c30-bc44-4242653c0e96
md"We can add or remove elements anywhere in the array using `insert!` and `deleteat!`:"

# ╔═╡ 56d63b69-e6bf-4140-8494-9ad36af2382c
begin
	vectmp5 = copy(vec)
	insert!(vectmp5,2,15)
end

# ╔═╡ 4d7f948d-666d-4038-9a68-fe2a8ab09b82
begin
	vectmp6 = copy(vec)
	deleteat!(vectmp6,2)
end

# ╔═╡ 7d06e1cb-ebe9-42cb-a717-f5ee01f3492a
md"### Concatenate"

# ╔═╡ 2cb2f5ad-10e4-4369-ae9e-3ab09665de68
md"We can concatenate several one-dimensional arrays together using `append!` (at the end) and `prepend!` (at the beginning):"

# ╔═╡ 3337c025-6bf9-4fae-b8f9-7fb0bb004451
vec2 = [30,31,32]

# ╔═╡ 4fbc1851-7d17-4907-95c8-cb7862c4e73a
begin
	vectmp7 = copy(vec)
	append!(vectmp7,vec2)
end

# ╔═╡ aab84642-b9fa-48f3-bb97-76d0d8b01d79
begin
	vectmp8 = copy(vec)
	prepend!(vectmp8,vec2)
end

# ╔═╡ 0e38f381-b0b4-4137-86c1-118077f3bf89
md"And you can still use `deleteat!`, but now using a range:"

# ╔═╡ dcb40e07-31cf-449c-b37a-f6ebb89fb65a
begin
	vectmp9 = copy(vec)
	deleteat!(vectmp9,1:3)
end

# ╔═╡ 8d8e2770-ffb7-4c4b-a624-786edac752e4
md"For more complex arrays we can use `hcat` and `vcat` (and it works with vectors too):"

# ╔═╡ 6a8e56f1-f797-45d3-896e-64d5c39f6023
hcat(vec,vec)

# ╔═╡ 7dcbf14a-04a4-4693-9918-b4164853cfaf
hcat(mat,[10,11,12,13])

# ╔═╡ a3883207-0020-4959-9b70-049998963bf0
vcat(mat,[10 11])

# ╔═╡ 4e1cacbc-2824-4823-a9a2-f1736129bf2f
md"We can also simply use `[` to make a new array:"

# ╔═╡ eec05d6a-ae8c-49c0-88e3-2fa19a458dfc
[mat; [10 11]]

# ╔═╡ c886a9f1-ece9-4e30-a42e-819de1e60794
[mat [10,11,12,13]]

# ╔═╡ 8a7c9b3f-e4e8-49d3-ae2f-da6d116963d7
md"### Reshape"

# ╔═╡ 62635c8c-8f21-42a9-9522-3440872aed4d
md"We can change the dimensions of an Array, or reshaping it using `reshape`. For example we can reshape vec to a 2x2 matrix:"

# ╔═╡ ff69b3ec-2c54-48ce-98d7-14fbf197982c
reshape(vec,2,2)

# ╔═╡ 17d44621-6264-466d-9459-a5cda803a7c3
md"""
Or reshape `mat` into a vector:
"""

# ╔═╡ 2c596432-756b-4670-b0b4-79d45664685f
reshape(mat,8)

# ╔═╡ aa7cc737-788c-499c-acc0-a0be06d03e10
md"""
If we are not sure about the size of our matrix, we can also use `:`:
"""

# ╔═╡ aa55dbcf-1491-4183-967a-959f72870fd0
reshape(mat,:)

# ╔═╡ 6aafc7b2-f28b-42a7-b181-d7385ba1359b
md"## Further information"

# ╔═╡ 959ea351-1b13-49ef-9c96-0fe1d4c7b886
md"""
There are some basic functions such as `length` or `size` to get the length of an array, and the size for each dimension:
"""

# ╔═╡ 414af929-a7d3-45b6-ac05-3fe05870bae9
length(mat)

# ╔═╡ e48e5fb2-5cbd-4413-b384-289599a588c1
size(mat)

# ╔═╡ f01961ff-e8ae-4219-9bc4-669c18f12f9c
md"""
We will see in another notebook that arrays are very powerfull, and that we can make more complex operations such as finding a value in an array using a function, sorting, or doing mathematical operations.

If you want more information, you can visit the [official documentation here](https://docs.julialang.org/en/v1/manual/arrays/)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.18"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0ec322186e078db08ea3e7da5b8b2885c099b393"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.0"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "57312c7ecad39566319ccf5aa717a20788eb8c1f"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.18"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─bb7a1f18-84c1-4c0a-807d-dc5575a7dd0c
# ╟─c879bad0-7944-11eb-3cc5-6d068e87585a
# ╟─0a8c911e-9dd2-4f7e-868b-1cd845135c53
# ╟─34558e07-9821-49a5-957b-7ef713fa9a6b
# ╟─cda72a9e-6ec1-4746-a6f8-e0926b5a50d9
# ╟─daf15091-45a3-4419-beb3-643d324a1f70
# ╟─92370118-f98d-453b-99d2-0bf053204278
# ╠═45fcb20b-e604-458c-80d2-d96324d8a63b
# ╟─8933fc3c-ccfe-4fa8-98af-270fc8e15280
# ╟─d3e2cdd5-2dd8-45e2-8f3d-f15d39c18d1b
# ╠═49090226-6636-4464-b496-1e0b83e5655f
# ╟─792b67f0-f1f1-4b4a-a720-aecffab167c6
# ╠═7a150c77-11b6-4525-836f-1e371ba91e70
# ╟─94ad54e0-b6a0-4872-94e4-532cfd237acf
# ╠═9b364eae-7060-4034-9924-1490fd88d275
# ╟─70cb6051-c849-4d14-84b6-3d7b48b42ebc
# ╠═6743a46c-4f24-4eca-b4a3-18de4dc9618a
# ╟─34f5b3dc-0331-4e32-aaf2-7807fb1855e9
# ╠═72fd91e2-e83d-4173-8c29-c3f7e3ca7fb2
# ╠═189853b7-212b-41ec-95ba-be4993b842ac
# ╟─f881789b-9cae-4de3-86e3-624f732aed82
# ╟─7e554061-1ed9-4cdc-aa02-f6d7565106b7
# ╠═7f9e7e08-6021-4406-a09c-419f4c64b30d
# ╟─d13f18d8-fa64-4d80-97c7-dd392358f579
# ╠═78959d1d-b522-4567-8353-73fac5e9d120
# ╟─c45e1610-c944-40e5-b7c6-c48c8bd1c626
# ╠═0fde7f3d-cde2-4381-89b7-4234d76f3855
# ╟─9b5093ba-e389-4a7b-a152-40311dcea1e6
# ╠═48f5b4e6-83ad-4886-ac57-7aba4e42c7ca
# ╟─c95603b0-3ddd-4f9a-b59d-12af2b2f492b
# ╠═7e2fc580-8ea9-4409-b3e2-2e94a550a81a
# ╠═088555de-40a0-4e84-8c64-8b4c8121dcca
# ╟─2b9c0567-3cf7-47e8-a595-ab4f2f0fcb18
# ╠═57d80a18-57cb-424c-afca-dba5924bb8f2
# ╟─5e4a882d-608d-47d4-b594-1b33857646ac
# ╟─2337fdaf-40d9-4d21-bb19-ac98defef41a
# ╠═5009bbdb-79fe-49f0-b642-0a866e202136
# ╟─ad0e06ce-7670-4205-94c3-81aa576fa30f
# ╠═f3aaca99-af7b-468d-b2eb-571a2c664903
# ╟─fecb5933-8410-4622-98ea-862c8fb8ab59
# ╠═d9f43a57-e061-4f4a-9bd3-c4dcdac4fb25
# ╟─d3cc12a3-8251-444c-93e9-c5c09527d034
# ╠═a152378c-db20-45f2-af28-6546b577d9fe
# ╟─0abc7154-9a2c-4547-9ebd-02f80d2b8ecf
# ╠═6ebef395-d692-448f-ae12-491c385f8dd2
# ╟─a04c4b0c-365c-4032-8884-bd2b2df040bf
# ╟─5778335e-0b31-423c-9ca1-8a8b9add3680
# ╟─50fedd78-95c0-4c24-86e0-d6fe7ab99a65
# ╠═1baf26a7-8af0-4473-86d9-497ab2546051
# ╠═2b4dbfbb-173d-43fb-8e3f-0827536814e3
# ╟─8b7fa713-9b89-46dc-b4fa-da2e348fd18e
# ╟─8f74faf5-829e-4c08-8463-5695e574af82
# ╠═75a44b16-5279-4d08-ba6a-dce7ea21db15
# ╟─8fef9d3e-b549-4a80-8314-d0a19cf6f1c1
# ╠═ce4bafc7-4ce2-4349-a0ec-d350ef4853f2
# ╟─ae116904-1ced-43bc-a007-997123a32ffd
# ╠═5d3771f0-0b75-4ea2-8fd5-1047e1a35779
# ╟─d642b188-9e7c-493f-a9e4-711cf03227e8
# ╠═a95aaf6c-ff91-4c3c-a83a-d71afdb94878
# ╟─91e738e6-9238-4a48-ae58-0bea688e6c90
# ╠═5e2387f1-5dd7-4b49-91e4-96d057ed31fa
# ╟─60efc230-9ada-4de9-896a-c5d37d7fa9e0
# ╠═4d38a50e-57ff-461d-8af9-9f1f188521e4
# ╟─b4a2e4b8-f16f-4fbd-baba-687f7468ca7f
# ╠═a939313b-1c86-492d-9993-545066bd8cd4
# ╠═fca126bf-b5d6-47b5-bbe9-994f7acad0fe
# ╟─862f073e-8006-4671-a8e0-a6cf2024ae68
# ╠═f630778e-1887-4fe3-9f35-c4e286e45cba
# ╟─c58f05ad-dd0b-482f-9891-514b960b8bb4
# ╟─67208966-c604-4889-b06e-f385647990d1
# ╟─c764ac39-8d2a-4fa1-b3f3-96a24ba5527d
# ╟─9fc4cb37-69fd-446c-95fe-30abb7559601
# ╠═47e8bdbc-ad36-447c-8f79-f78626814be1
# ╟─b0d77b3e-9674-449a-be5f-a5790a4d4089
# ╟─d94f73bb-0e3f-42ea-8e8f-7e3e982836fa
# ╠═23a8f2fa-9153-4eec-a876-db97a999b943
# ╟─f530fe3b-1b0d-472a-b330-d12a46e68c22
# ╟─8fbfb20a-5f5e-422e-9b59-90513050b3b8
# ╟─60608c63-2b59-474e-92a2-3800e6254792
# ╠═f153b213-1d4a-4481-ba83-29b1f8131963
# ╠═a99f7d96-654a-46b0-91a2-9120ab30af44
# ╟─1cb1ec65-48d9-47f1-8a3d-f3767fa03ad2
# ╟─33c34e8f-c5b6-43a8-a329-12fc55949a2e
# ╟─fad8c97c-7cd6-4c30-bc44-4242653c0e96
# ╠═56d63b69-e6bf-4140-8494-9ad36af2382c
# ╠═4d7f948d-666d-4038-9a68-fe2a8ab09b82
# ╟─7d06e1cb-ebe9-42cb-a717-f5ee01f3492a
# ╟─2cb2f5ad-10e4-4369-ae9e-3ab09665de68
# ╠═3337c025-6bf9-4fae-b8f9-7fb0bb004451
# ╠═4fbc1851-7d17-4907-95c8-cb7862c4e73a
# ╠═aab84642-b9fa-48f3-bb97-76d0d8b01d79
# ╟─0e38f381-b0b4-4137-86c1-118077f3bf89
# ╠═dcb40e07-31cf-449c-b37a-f6ebb89fb65a
# ╟─8d8e2770-ffb7-4c4b-a624-786edac752e4
# ╠═6a8e56f1-f797-45d3-896e-64d5c39f6023
# ╠═7dcbf14a-04a4-4693-9918-b4164853cfaf
# ╠═a3883207-0020-4959-9b70-049998963bf0
# ╟─4e1cacbc-2824-4823-a9a2-f1736129bf2f
# ╠═eec05d6a-ae8c-49c0-88e3-2fa19a458dfc
# ╠═c886a9f1-ece9-4e30-a42e-819de1e60794
# ╟─8a7c9b3f-e4e8-49d3-ae2f-da6d116963d7
# ╟─62635c8c-8f21-42a9-9522-3440872aed4d
# ╠═ff69b3ec-2c54-48ce-98d7-14fbf197982c
# ╟─17d44621-6264-466d-9459-a5cda803a7c3
# ╠═2c596432-756b-4670-b0b4-79d45664685f
# ╟─aa7cc737-788c-499c-acc0-a0be06d03e10
# ╠═aa55dbcf-1491-4183-967a-959f72870fd0
# ╟─6aafc7b2-f28b-42a7-b181-d7385ba1359b
# ╟─959ea351-1b13-49ef-9c96-0fe1d4c7b886
# ╠═414af929-a7d3-45b6-ac05-3fe05870bae9
# ╠═e48e5fb2-5cbd-4413-b384-289599a588c1
# ╟─f01961ff-e8ae-4219-9bc4-669c18f12f9c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
import Base.Broadcast: broadcasted

# ZonedDateTime arithmetic
Base.:(+)(x::ZonedDateTime) = x
Base.:(-)(x::ZonedDateTime, y::ZonedDateTime) = x.utc_datetime - y.utc_datetime

function Base.:(+)(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(DateTime(zdt, Local) + p, timezone(zdt))
end
function Base.:(+)(zdt::ZonedDateTime, p::TimePeriod)
    return ZonedDateTime(DateTime(zdt, UTC) + p, timezone(zdt); from_utc=true)
end
function Base.:(-)(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(DateTime(zdt, Local) - p, timezone(zdt))
end
function Base.:(-)(zdt::ZonedDateTime, p::TimePeriod)
    return ZonedDateTime(DateTime(zdt, UTC) - p, timezone(zdt); from_utc=true)
end

function broadcasted(::typeof(+), r::StepRange{ZonedDateTime{FixedTimeZone}}, p::DatePeriod)
    start, step, stop = first(r), Base.step(r), last(r)
    return StepRange(start + p, step, stop + p)
end

function broadcasted(::typeof(+), r::StepRange{ZonedDateTime{VariableTimeZone}}, p::DatePeriod)
    start, step, stop = first(r), Base.step(r), last(r)

    # Since the local time + period can result in an invalid local datetime when working with
    # VariableTimeZones we will use `first_valid` and `last_valid` which avoids issues with
    # non-existent and ambiguous dates.

    start = first_valid(DateTime(start, Local) + p, timezone(start), step)
    stop = last_valid(DateTime(stop, Local) + p, timezone(stop), step)
    return StepRange(start, step, stop)
end

broadcasted(::typeof(-), r::StepRange{ZonedDateTime}, p::DatePeriod) = broadcast(+, r, -p)


text:
println("\n\n\nStarting runtests.jl $(join(ARGS, " ")) ...")
pushfirst!(LOAD_PATH, joinpath(dirname(@__FILE__), "../src"))

using Test, NamedDimensions, FunctionalData, StatsBase

@testset "basic" begin
    data = [1 2 3;4 5 6]
    n = named(data, :a, :b)
    r = named(n, :a => 2)
    @test r.data == [4,5,6]
    @test r.names == [:b]
    @test array(r) == [4,5,6]

    @test ndims(n) == 2
    @test size(n) == (2,3)
    @test size(n,1) == 2
    @test size(n,:a) == 2
    @test size(n,2) == 3
    @test size(n,:b) == 3
    @test length(n) == 6
    @test len(n) == 3

    r = @p named n :a
    @test r.names == [:b, :a]
    @test size(r.data) == (3,2)

    n = named(data, :c)
    @test n.names == [:dimA, :c]

    b = ones(2,3,4)
    n = named(b, :a, :b, :c)
    r = @p named n :a :b
    @test r.names == [:c, :a, :b]
    @test size(r.data) == (4,2,3)

    r = @p named n :a :b | named :c :a
    @test r.names == [:b, :c, :a]
    @test size(r.data) == (3,4,2)

    r = @p named data :a :b | named (:b=>2:3)
    @test r.names == [:a, :b]
    @test r.data == [2 3; 5 6]
    r = @p named data :a :b | named (:a=>1:2)
    @test r.names == [:b, :a]
    @test r.data == data'

    n = named(data, :a, :b)
    @test n[:a] == at(n,:a)
    @test (@p n :b (2:3)) == named(n, :b => 2:3)
    @test (@p n :b (2:3) :a 2) == named([5,6], :b)

    @testset "linear indexing" begin
        @test n[1] == 1
        @test n[2:3] == [4,2]
    end

    buf = IOBuffer()
    showinfo(buf,n)
    @test String(take!(buf)) == "----  2 a x 3 b\n    type: Array{Int64,2}   size: (2, 3)\n    min:  1   max: 6\n    mean: 3.5   median: 3.5\n"
    showinfo(buf, n, "n in runtests")
    @test String(take!(buf)) == "n in runtests  --  2 a x 3 b\n    type: Array{Int64,2}   size: (2, 3)\n    min:  1   max: 6\n    mean: 3.5   median: 3.5\n"

    @test minimum(n) == 1
    # @test maximum(n) == 6
    # @test mean(n) == 3.5
    # @test minimum(n,1) == named([1,2,3],:b)
    # @test minimum(n,:a) == named([1,2,3],:b)
    # @test minimum(n,2) == named([1,4],:a)
    # @test minimum(n,:b) == named([1,4],:a)
    # @test maximum(n,2) == named([3,6],:a)
    # @test std(n,2) == named(vec(std(n.data,2)), :a)
    # @test std(n,:a) == named(vec(std(n.data,1)), :b)
    # @test n .+ n == named([2 4 6;8 10 12], :a, :b)
    # @test n .+ [1 2 3] == named([2 4 6; 5 7 9], :a, :b)
    # @test array(named(ones(2,3,4)) .+ ones(2)) == 2*ones(2,3,4)
    # @test array(named(ones(2,3,4)) .+ ones(2,1)) == 2*ones(2,3,4)
    # @test array(named(ones(2,3,4)) .+ ones(2,1,1)) == 2*ones(2,3,4)
    # @test n .+ [1 2]' == named([2 3 4; 6 7 8], :a, :b)
    # @test [1 2 3] .+ n == named([2 4 6; 5 7 9], :a, :b)
    # @test [1 2]' .+ n == named([2 3 4; 6 7 8], :a, :b)
    # @test n .* data == named(data.^2, :a, :b)
    # @test n * data' == named([14 32; 32 77], :a, :dimB)
    # @test (@p minimum n :b | minimum :a) == 1

    # @test NamedDimensions.dimnames(named(rand(2,3),:a,:b), 1) == [:dimA]
    # @test NamedDimensions.dimnames(named(rand(2,3),:dimA,:b), 1) == [:dimB]
    # @test NamedDimensions.dimnames(named(rand(2,3),:dimA,:dimB), 1) == [:dimC]
    # @test NamedDimensions.dimnames(named(rand(2,3),:dimC,:dimB), 1) == [:dimA]
    # @test NamedDimensions.dimnames(named(rand(2,3),:dimC,:dimB), 3) == [:dimA, :dimD, :dimE]
end

@testset "FunctionalData" begin
    data = [1 2 3;4 5 6]
    n = named(data, :a, :b)

    @test concat(n,n) == named(concat(data,data), :a, :b)
    @test flatten([n,n]) == named(concat(data,data), :a, :b)
    @test concat(n,n) == named(concat(data,data), :a, :b)
    @test flatten([n,n]) == named(concat(data,data), :a, :b)
    @test stack([n,n], :c) == named(stack(Any[data,data]), :a, :b, :c)
    @test unstack(stack([n,n], :c)) == [named(data, :a, :b), named(data, :a, :b)]

    @test fst(n) == named([1,4], :a)
    @test last(n) == named([3,6], :a)
    @test last(n,:b) == named([3,6], :a)
    @test last(n,:a) == named([4,5,6], :b)
    @test take(n,2) == named([1 2; 4 5], :a, :b)
    @test takelast(n,2) == named([2 3; 5 6], :a, :b)
    @test part(n,3:3) == named(col([3;6]), :a, :b)
    @test drop(n,2) == part(n, 3:3)

    @test named(n,:a=>:end) == named([4,5,6],:b)
    @test named(n,:b=>:end) == named([3,6], :a)

    # @test (@p map n id) == n
    # @test (@p map n array) == n.data

    # @test (@p map2 n n plus) == named(data+data, :a, :b)
    # @test (@p map2 n data plus) == named(data+data, :a, :b)
    # @test (@p map2 data n plus) == named(data+data, :a, :b)
end

println("done!")


text:
module Recurrences

export @lrs
export Seq, ClosedForm

using MacroTools: postwalk, @capture, unblock
using Combinatorics
import AbstractAlgebra: Ring, RingElement, Field, FieldElement, Generic.Frac, RingElem
import AbstractAlgebra: parent_type, elem_type, base_ring, parent, needs_parentheses, displayed_with_minus_in_front, canonical_unit
import AbstractAlgebra: zero!, mul!, add!, addeq!
using Nemo
using SymEngine
using LinearAlgebra
import Base: convert

const RExpr = Union{Expr,Symbol,Number}

include("helper.jl")
include("utils.jl")
include("hyperterm.jl")
include("sequence.jl")
include("zuercher.jl")
include("rectypes.jl")
include("recsystem.jl")
include("petkovsek.jl")
include("macros.jl")

end # module


text:
# a basic test to  debug the spherical vector implementation
# -----------------------------------------------------------------

using ACE, StaticArrays, LinearAlgebra
mat(X::AbstractVector{<:StaticVector}) = hcat(X...)
rotc = ACE.Rotations3D.Rot3DCoeffs()

##
L = 1
Lmax = 3
U = fill( SA[0.0im, 0.0im, 0.0im], ((Lmax+1)^2, (Lmax+1)^2) )
for l = 1:2:Lmax, m = -l:l, μ = -l:l
   u = ACE.Rotations3D.local_cou_coe(rotc, SA[l], SA[m], SA[μ], L, 1)
   row = ACE.SphericalHarmonics.index_y(l, m)
   col = ACE.SphericalHarmonics.index_y(l, μ)
   U[row, col] = u
end

##

# check for  F(Qr) = D F(r)  ⇔   D' F(Qr) = F(r)

sh = ACE.SphericalHarmonics.SHBasis(Lmax)
r = @SVector randn(3); r /= norm(r)
Y = ACE.evaluate(sh, r)
UxY = U * Y

Q, D = ACE.Wigner.rand_QD(1)
Y_Q = ACE.evaluate(sh, Q * r)
DxUxY_D = Ref(D') .* (U * Y_Q)

DxUxY_D ≈ UxY

##

using ACE
using StaticArrays, Random, Printf, Test, LinearAlgebra, ACE.Testing
using ACE: evaluate, evaluate_d, SymmetricBasis, NaiveTotalDegree, PIBasis
using ACE.Random: rand_rot, rand_refl
using ACE.Wigner

# construct the 1p-basis
D = ACE.NaiveTotalDegree()
maxdeg = 6
ord = 1
B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, D = D, pin = 0, pout = 0)
# generate a configuration
nX = 1
Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

##

L = 1
ord = 1
φ = ACE.SphericalVector(L; T = ComplexF64)
pibasis = PIBasis(B1p, ord, 3; property = φ, isreal=false)

spec = ACE.get_spec(pibasis)
# specodd = filter(b -> isodd(b[1].l) && (b[1].n == 1), spec)
specf = [ (n=1, l=l, m=m) for l = 0:Lmax for m = -l:l ]
If = [ findall(isequal([b]), spec)[1] for b in specf ]

U = fill( SA[0.0im, 0.0im, 0.0im], length(specf), length(pibasis))
for (i, b) in enumerate(specf)
   for μ = -b.l:b.l
      bμ = (n = b.n, l = b.l, m = μ)
      iμ = findall(isequal([bμ]), spec)[1]
      u = ACE.Rotations3D.local_cou_coe(rotc, SA[b.l], SA[b.m], SA[μ], L, 1)
      U[i, iμ] = u
   end
end

##

Xs1  = Xs[1:1]
cfg1 = ACEConfig(Xs1)
AA = evaluate(pibasis, cfg1)
Y = evaluate(sh, Xs[1].rr)
 AA[If] / AA[If[1]] ≈ Y/Y[1]
##

Xs1  = Xs[1:1]
cfg1 = ACEConfig(Xs1)
AA = evaluate(pibasis, cfg1)
UxAA = U * AA

Q, D = ACE.Wigner.rand_QD(1)
Qcfg = ACEConfig( Ref(Q) .* Xs1 )
DtxUxAA_Q = Ref(D') .* (U * evaluate(pibasis, Qcfg))
DtxUxAA_Q ≈ UxAA

##

mat(DtxUxAA_Q - UxAA)[:,1]

##


text:
function nummax(A)
    max = A[firstindex(A)]
    count = 1
    for i in (firstindex(A):lastindex(A))
        if A[i] > max
            max = A[i]
            count = 0
        end
        if A[i] == max
            count += 1
        end
    end
    return count
end

print(nummax([3,4,1,2,7,0,1,9,2,3,3,4,5,6]))


text:
"""
Delete named dataset.

The driver will attempt to delete the named dataset in a driver specific
fashion. Full featured drivers will delete all associated files, database
objects, or whatever is appropriate. The default behaviour when no driver
specific behaviour is provided is to attempt to delete the passed name as a
single file.

It is unwise to have open dataset handles on this dataset when it is deleted.
"""
delete(dataset::Dataset, filename::AbstractString) =
    GDAL.deletedataset(dataset, filename)

"""
Rename a dataset.

This may including moving the dataset to a new directory or filesystem.

It is unwise to have open dataset handles on this dataset when it is being
renamed.
"""
rename(dataset::Dataset,
       newname::AbstractString,
       oldname::AbstractString) =
    GDAL.renamedataset(dataset, newname, oldname)

"Copy the files associated with a dataset."
copyfiles(dataset::Dataset,
          newname::AbstractString,
          oldname::AbstractString) =
    GDAL.copydatasetfiles(dataset, newname, oldname)

"""
Copy all dataset raster data.

This function copies the complete raster contents of one dataset to another
similarly configured dataset. The source and destination dataset must have the
same number of bands, and the same width and height. The bands do not have to
have the same data type.

This function is primarily intended to support implementation of driver
specific `CreateCopy()` functions. It implements efficient copying, in
particular "chunking" the copy in substantial blocks and, if appropriate,
performing the transfer in a pixel interleaved fashion.

Currently the only papszOptions value supported are : "INTERLEAVE=PIXEL" to
force pixel interleaved operation and "COMPRESSED=YES" to force alignment on
target dataset block sizes to achieve best compression. More options may be
supported in the future.
"""
copywholeraster(source::Dataset, dest::Dataset) =
    GDAL.datasetcopywholeraster(source, dest, Ptr{Ptr{UInt8}}(pointer(C_NULL)),
                                ProgressFunc(C_NULL), C_NULL)

"""
Create a copy of a dataset.

This method will attempt to create a copy of a raster dataset with the
indicated filename, and in this drivers format. Band number, size, type,
projection, geotransform and so forth are all to be copied from the
provided template dataset.

Note: many sequential write once formats (such as JPEG and PNG) don't implement
the `Create()` method but do implement this `CreateCopy()` method. If the
driver doesn't implement `CreateCopy()`, but does implement `Create()` then the
default `CreateCopy()` mechanism built on calling `Create()` will be used.

It is intended that `CreateCopy()` will often be used with a source dataset
which is a virtual dataset allowing configuration of band types, and other
information without actually duplicating raster data (see the VRT driver). This
is what is done by the gdal_translate utility for example.

This function will validate the creation option list passed to the driver
with the `GDALValidateCreationOptions()` method. This check can be disabled by
defining the configuration option `GDAL_VALIDATE_CREATION_OPTIONS=NO`.

After you have finished working with the returned dataset, it is required to
close it with `GDALClose()`. This does not only close the file handle, but also
ensures that all the data and metadata has been written to the dataset
(`GDALFlushCache()` is not sufficient for that purpose).

In some situations, the new dataset can be created in another process through
the GDAL API Proxy mechanism.

### Parameters
* `filename`      the name for the new dataset. UTF-8 encoded.
* `dataset`       the dataset being duplicated.
* `strict`        `TRUE` if the copy must be strictly equivelent, or more
normally `FALSE` if the copy may adapt as needed for the output format.
* `options`       additional format dependent options controlling creation
of the output file. `The APPEND_SUBDATASET=YES` option can be specified to
avoid prior destruction of existing dataset.

### Returns
a pointer to the newly created dataset (may be read-only access).
"""
createcopy(filename::AbstractString, dataset::Dataset, driver::Driver;
           strict::Bool = false) =
    GDAL.createcopy(driver, filename, dataset, strict, C_NULL,
                    ProgressFunc(C_NULL), C_NULL)

function createcopy(f::Function, args...; kwargs...)
    ds = createcopy(args...; kwargs...)
    try
        f(ds)
    finally
        GDAL.close(ds)
    end
end

create(filename::AbstractString,
       driver::Driver;
       width::Integer = 0,
       height::Integer = 0,
       nbands::Integer = 0,
       dtype::DataType = Any,
       options::Vector{ASCIIString} = Vector{ASCIIString}()) =
    GDAL.create(driver, filename, width, height, nbands,
                _gdaltype[dtype], Ptr{Ptr{UInt8}}(C_NULL))

create(filename::AbstractString,
       drivername::AbstractString;
       width::Integer = 0,
       height::Integer = 0,
       nbands::Integer = 0,
       dtype::DataType = Any) =
    GDAL.create(getdriver(drivername), filename, width, height,
                nbands, _gdaltype[dtype], Ptr{Ptr{UInt8}}(C_NULL))

function create(f::Function, args...; kwargs...)
    ds = create(args...; kwargs...)
    try
        f(ds)
    finally
        GDAL.close(ds)
    end
end

"""
Open a raster file as a GDALDataset.

This function will try to open the passed file, or virtual dataset name by
invoking the Open method of each registered `GDALDriver` in turn. The first
successful open will result in a returned dataset. If all drivers fail then
`NULL` is returned and an error is issued.

Several recommendations:

If you open a dataset object with `GA_Update` access, it is not recommended to
open a new dataset on the same underlying file. The returned dataset should
only be accessed by one thread at a time. To use it from different threads,
you must add all necessary code (mutexes, etc.) to avoid concurrent use of the
object. (Some drivers, such as GeoTIFF, maintain internal state variables that
are updated each time a new block is read, preventing concurrent use.)

Starting with GDAL 1.6.0, if `GDALOpenShared()` is called on the same
`pszFilename` from two different threads, a different `GDALDataset` object will
be returned as it is not safe to use the same dataset from different threads,
unless the user does explicitly use mutexes in its code.

For drivers supporting the VSI virtual file API, it is possible to open a file
in a `.zip` archive (see `VSIInstallZipFileHandler()`), a `.tar/.tar.gz/.tgz`
archive (see `VSIInstallTarFileHandler()`), or a HTTP / FTP server
(see `VSIInstallCurlFileHandler()`)

In some situations (dealing with unverified data), the datasets can be opened
in another process through the GDAL API Proxy mechanism.

See also: `GDALOpenShared()`, `GDALOpenEx()`

### Parameters
* `filename`  the name of the file to access. In the case of exotic drivers
this may not refer to a physical file, but instead contain information for the
driver on how to access a dataset. It should be in UTF-8 encoding.
* `access`      the desired access, either `GA_Update` or `GA_ReadOnly`. Many
drivers support only read only access.
"""
function read(filename::AbstractString; shared::Bool = false)
    if shared
        return GDAL.openshared(filename, GDAL.GA_ReadOnly)
    else
        return GDAL.openex(filename, GDAL.GA_ReadOnly, C_NULL, C_NULL, C_NULL)
    end
end

function read(f::Function, args...; kwargs...)
    ds = read(args...; kwargs...)
    try
        f(ds)
    finally
        GDAL.close(ds)
    end
end

function update(filename::AbstractString; shared::Bool = false)
    if shared
        return GDAL.openshared(filename, GDAL.GA_Update)
    else
        return GDAL.openex(filename, GDAL.GA_Update, C_NULL, C_NULL, C_NULL)
    end
end

function update(f::Function, args...; kwargs...)
    ds = update(args...; kwargs...)
    try
        f(ds)
    finally
        GDAL.close(ds)
    end
end

function write(filename::AbstractString,
               dataset::Dataset;
               strict::Bool = false)
    (dataset == C_NULL) && error("Can't write closed dataset")
    GDAL.close(createcopy(filename, dataset, getdriver(dataset), strict=strict))
end

function write(dataset::Dataset,
               filename::AbstractString,
               driver::Driver;
               strict::Bool = false)
    (dataset == C_NULL) && error("Can't write closed dataset")
    GDAL.close(createcopy(filename, dataset, driver, strict=strict))
end

function write(dataset::Dataset,
               filename::AbstractString,
               drivername::AbstractString;
               strict::Bool = false)
    (dataset == C_NULL) && error("Can't write closed dataset")
    GDAL.close(createcopy(filename, dataset, getdriver(name), strict=strict))
end

function write(f::Function, args...; kwargs...)
    ds = write(args...; kwargs...)
    try
        f(ds)
    finally
        GDAL.close(ds)
    end
end

"Fetch raster width in pixels."
width(dataset::Dataset) = GDAL.getrasterxsize(dataset)

"Fetch raster height in pixels."
height(dataset::Dataset) = GDAL.getrasterysize(dataset)

"Fetch the number of raster bands on this dataset."
nraster(dataset::Dataset) = GDAL.getrastercount(dataset)

"Fetch the number of feature layers on this dataset."
nlayer(dataset::Dataset) = GDAL.datasetgetlayercount(dataset)

"Fetch the driver that the dataset was created with"
getdriver(dataset::Dataset) = GDAL.getdatasetdriver(dataset)

"""
Add a band to a dataset.

This method will add a new band to the dataset if the underlying format
supports this action. Most formats do not.

Note that the new `GDALRasterBand` is not returned. It may be fetched after
successful completion of the method through `GetRasterBand(GetRasterCount())`
as the newest band will always be the last band.

### Parameters
* `eType`           the data type of the pixels in the new band.
* `options`    a list of NAME=VALUE option strings. The supported options
are format specific. `NULL` may be passed by default.
"""
addband(dataset::Dataset, eType::GDAL.GDALDataType) =
    GDAL.addband(dataset, eType, C_NULL)

addband(dataset::Dataset, eType::GDAL.GDALDataType) =
    GDAL.addband(dataset, eType, Ptr{Ptr{UInt8}}(C_NULL))

"""
Fetch files forming dataset.

Returns a list of files believed to be part of this dataset. If it returns an
empty list of files it means there is believed to be no local file system files
associated with the dataset (for instance a virtual dataset). The returned file
list is owned by the caller and should be deallocated with `CSLDestroy()`.

The returned filenames will normally be relative or absolute paths depending on
the path used to originally open the dataset. The strings will be UTF-8 encoded
"""
filelist(dataset::Dataset) =
    loadstringlist(GDAL.C.GDALGetFileList(Ptr{Void}(dataset)))

"Fetch a layer by index (between 0 and GetLayerCount()-1)"
fetchlayer(dataset::Dataset, i::Integer) = GDAL.datasetgetlayer(dataset, i)

"Fetch a feature layer for a dataset from its name"
fetchlayer(dataset::Dataset, name::AbstractString) =
    GDAL.datasetgetlayerbyname(dataset, name)

"""
Delete the indicated layer (at index i) from the datasource.

### Returns
OGRERR_NONE on success, or OGRERR_UNSUPPORTED_OPERATION if deleting layers
is not supported for this datasource.
"""
function deletelayer(dataset::Dataset, i::Integer)
    result = GDAL.datasetdeletelayer(dataset, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete layer")
    end
end

"""
This function attempts to create a new layer on the dataset with the indicated
name, coordinate system, geometry type.

### Parameters
* **dataset**: the dataset
* **name**: the name for the new layer. This should ideally not match any
    existing layer on the datasource.
* **spatialref**: the coordinate system to use for the new layer, or NULL if
    no coordinate system is available.
* **geomtype**: the geometry type for the layer. Use wkbUnknown if there are no
    constraints on the types geometry to be written.
* **papszOptions**: a StringList of name=value (driver-specific) options.
"""
createlayer(dataset::Dataset, name::AbstractString,
            geomtype::GDAL.OGRwkbGeometryType = GDAL.wkbUnknown) =
    GDAL.datasetcreatelayer(dataset, name, SpatialRef(C_NULL), geomtype, C_NULL)

createlayer(dataset::Dataset, name::AbstractString, spatialref::SpatialRef,
            geomtype::GDAL.OGRwkbGeometryType = GDAL.wkbUnknown) =
    GDAL.datasetcreatelayer(dataset, name, spatialref, geomtype, C_NULL)

"""
Duplicate an existing layer.

### Parameters
* **dataset**: the dataset handle.
* **layer**: source layer.
* **name**: the name of the layer to create.
* **papszOptions**: a StringList of name=value (driver-specific) options.
"""
copylayer(dataset::Dataset,layer::FeatureLayer, name::AbstractString) =
    GDAL.datasetcopylayer(dataset, layer, name, C_NULL)

"""
Execute an SQL statement against the data store.

### Parameters
* **dataset**: the dataset handle.
* **query**: the SQL statement to execute.
* **spatialfilter**: geometry which represents a spatial filter. Can be NULL.
* **dialect**: allows control of the statement dialect. If set to NULL, the
    OGR SQL engine will be used, except for RDBMS drivers that will use their
    dedicated SQL engine, unless OGRSQL is explicitly passed as the dialect.
    Starting with OGR 1.10, the SQLITE dialect can also be used.

### Returns
an OGRLayer containing the results of the query.
Deallocate with ReleaseResultSet().
"""
executesql(dataset::Dataset, query::AbstractString) =
    GDAL.datasetexecutesql(dataset, query, Geometry(C_NULL), "")

executesql(dataset::Dataset, query::AbstractString, dialect::AbstractString) =
    GDAL.datasetexecutesql(dataset, query, Geometry(C_NULL), dialect)

executesql(dataset::Dataset, query::AbstractString, spatialfilter::Geometry,
           dialect::AbstractString) =
    GDAL.datasetexecutesql(dataset, query, spatialfilter, dialect)

function executesql(f::Function, dataset::Dataset, args...)
    result = executesql(dataset, args...)
    try
        f(result)
    finally
        releaseresultset(dataset, result)
    end
end

"""
Release results of ExecuteSQL().

### Parameters
* **dataset**: the dataset handle.
* **layer**: the result of a previous ExecuteSQL() call.
"""
releaseresultset(dataset::Dataset, layer::FeatureLayer) =
    GDAL.datasetreleaseresultset(dataset, layer)

"""
Fetch the affine transformation coefficients.

Fetches the coefficients for transforming between pixel/line (P,L) raster
space, and projection coordinates (Xp,Yp) space.

```julia
   Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
   Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];
```

In a north up image, `padfTransform[1]` is the pixel width, and
`padfTransform[5]` is the pixel height. The upper left corner of the upper left
pixel is at position `(padfTransform[0],padfTransform[3])`.

The default transform is `(0,1,0,0,0,1)` and should be returned even when a
`CE_Failure` error is returned, such as for formats that don't support
transformation to projection coordinates.

### Parameters
* `buffer`   a six double buffer into which the transformation will be placed.

### Returns
`CE_None` on success, or `CE_Failure` if no transform can be fetched.
"""
function getgeotransform!(dataset::Dataset, transform::Vector{Cdouble})
    @assert length(transform) == 6
    result = GDAL.getgeotransform(dataset, pointer(transform))
    (result == GDAL.CE_Failure) && error("Failed to get geotransform")
    transform
end

getgeotransform(dataset::Dataset) = getgeotransform!(dataset, Array(Cdouble, 6))

"Set the affine transformation coefficients."
function setgeotransform!(dataset::Dataset, transform::Vector{Cdouble})
    @assert length(transform) == 6
    result = GDAL.setgeotransform(dataset, pointer(transform))
    (result == GDAL.CE_Failure) && error("Failed to transform raster dataset")
end

"Get number of GCPs for this dataset. Zero if there are none."
ngcp(dataset::Dataset) = GDAL.getgcpcount(dataset)

"""
Fetch the projection definition string for this dataset in OpenGIS WKT format.

It should be suitable for use with the OGRSpatialReference class. When a
projection definition is not available an empty (but not `NULL`) string is
returned.

See also: http://www.gdal.org/ogr/osr_tutorial.html
"""
getproj(dataset::Dataset) = GDAL.getprojectionref(dataset)

"Set the projection reference string for this dataset."
function setproj(dataset::Dataset, projstring::AbstractString)
    result = GDAL.setprojection(dataset, projstring)
    (result == GDAL.CE_Failure) && error("Could not set projection")
end

# "Fetch a format specific internally meaningful handle."
# _getinternalhandle(dataset::GDALDatasetH, request::Ptr{UInt8}) =
#     GDALGetInternalHandle(dataset, request)::Ptr{Void}

# """
# Add one to dataset reference count. The reference is one after instantiation.

# ### Returns
# the post-increment reference count.
# """
# _referencedataset(dataset::GDALDatasetH) = GDALReferenceDataset(dataset)::Cint

# """
# Subtract one from dataset reference count. The reference is one after
# instantiation. Generally when the reference count has dropped to zero the
# dataset may be safely deleted (closed).

# ### Returns
# the post-decrement reference count.
# """
# _dereferencedataset(dataset::GDALDatasetH) =
#     GDALDereferenceDataset(dataset)::Cint

"""
Build raster overview(s).

If the operation is unsupported for the indicated dataset, then CE_Failure is
returned, and CPLGetLastErrorNo() will return CPLE_NotSupported.

### Parameters
* `pszResampling`   one of "NEAREST", "GAUSS", "CUBIC", "AVERAGE", "MODE",
"AVERAGE_MAGPHASE" or "NONE" controlling the downsampling method applied.
* `nOverviews`      number of overviews to build.
* `panOverviewList` the list of overview decimation factors to build.
* `nListBands`      number of bands to build overviews for in `panBandList`.
Build for all bands if this is 0.
* `panBandList`     list of band numbers.
* `pfnProgress`     a function to call to report progress, or `NULL`.
* `pProgressData`   application data to pass to the progress function.

### Returns
`CE_None` on success or `CE_Failure` if the operation doesn't work.

For example, to build overview level 2, 4 and 8 on all bands the following
call could be made:

```C
   int       anOverviewList[3] = { 2, 4, 8 };
   poDataset->BuildOverviews( "NEAREST", 3, anOverviewList, 0, NULL,
                              GDALDummyProgress, NULL );
```
"""
function buildoverviews(dataset::Dataset, overviewlist::Vector{Cint};
                        resampling::AbstractString = "NEAREST")
    result = GDAL.buildoverviews(dataset, resampling, length(overviewlist),
                                 overviewlist, 0, C_NULL,
                                 ProgressFunc(C_NULL), C_NULL)
    (result == GDAL.CE_Failure) && error("Failed to build overviews")
end

function buildoverviews(dataset::Dataset,
                        overviewlist::Vector{Cint};
                        bandlist::Vector{Cint} = Cint[],
                        resampling::AbstractString = "NEAREST")
    result = GDAL.buildoverviews(dataset, resampling, length(overviewlist),
                                 overviewlist, length(bandlist), bandlist,
                                 ProgressFunc(C_NULL), C_NULL)
    (result == GDAL.CE_Failure) && error("Failed to build overviews")
end

# """
# Fetch all open GDAL dataset handles.

# NOTE: This method is not thread safe. The returned list may change at any time
# and it should not be freed.

# ### Parameters
# * `pnCount`     integer into which to place the count of dataset pointers being
# returned.

# ### Returns
# a pointer to an array of dataset handles.
# """
# getopendatasets(hDS::Ptr{Ptr{GDALDatasetH}}, pnCount::Ptr{Cint}) =
#     GDALGetOpenDatasets(hdS, pnCount)

"""
Adds a mask band to the dataset.

The default implementation of the `CreateMaskBand()` method is implemented
based on similar rules to the .ovr handling implemented using the
`GDALDefaultOverviews` object. A TIFF file with the extension .msk will be
created with the same basename as the original file, and it will have one band.
The mask images will be deflate compressed tiled images with the same block
size as the original image if possible.

If you got a mask band with a previous call to `GetMaskBand()`, it might be
invalidated by CreateMaskBand(). So you have to call `GetMaskBand()` again.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Parameters
* `nFlags`  ignored. `GMF_PER_DATASET` will be assumed.

### Returns
`CE_None` on success or `CE_Failure` on an error.
"""
function createdatasetmaskband(dataset::Dataset)
    result = GDAL.createdatasetmaskband(dataset)
    (result == GDAL.CE_Failure) && error("Failed to create dataset mask band")
end

"""
Get output projection for GCPs.

The projection string follows the normal rules from `GetProjectionRef()`.

### Returns
internal projection string or `""` if there are no GCPs. It should not be
altered, freed or expected to last for long.
"""
getgcpproj(dataset::Dataset) = GDAL.getgcpprojection(dataset)

# """
# Fetch GCPs.

# ### Returns
# pointer to internal GCP structure list. It should not be modified, and may
# change on the next GDAL call.
# """
# _getgcps(dataset::GDALDatasetH) = GDALGetGCPs(dataset)::Ptr{GDAL_GCP}

# """
# Assign GCPs.

# This method assigns the passed set of GCPs to this dataset, as well as setting
# their coordinate system. Internally copies are made of the coordinate system
# and list of points, so the caller remains responsible for deallocating these
# arguments if appropriate.

# Most formats do not support setting of GCPs, even formats that can handle GCPs.
# These formats will return `CE_Failure`.

# ### Parameters
# * `nGCPCount`           number of GCPs being assigned.
# * `pasGCPList`          array of `nGCPCount` GCP structures to be assigned
# * `pszGCPProjection`    the new OGC WKT coordinate system to assign for the GCP
# output coordinates. Should be `""` if no output coordinate system is known.

# ### Returns
# `CE_None` on success, `CE_Failure` on failure (including if action is not
# supported for this format).
# """
# _setgcps(dataset::GDALDatasetH,
#          nGCPCount::Integer,
#          pasGCPList::Ptr{GDAL_GCP},
#          pszGCPProjection::Ptr{UInt8}) =
#     GDALSetGCPs(dataset, nGCPCount, pasGCPList, pszGCPProjection)::CPLErr

# """
#     GDALDatasetAdviseRead(GDALDatasetH hDS,
#                           int nXOff,
#                           int nYOff,
#                           int nXSize,
#                           int nYSize,
#                           int nBufXSize,
#                           int nBufYSize,
#                           GDALDataType eDT,
#                           int nBandCount,
#                           int * panBandMap,
#                           char ** papszOptions) -> CPLErr
# Advise driver of upcoming read requests.
# """
# function datasetadviseread{T <: GDALDatasetH}(hDS::Ptr{T},nDSXOff::Integer,nDSYOff::Integer,nDSXSize::Integer,nDSYSize::Integer,nBXSize::Integer,nBYSize::Integer,eBDataType::GDALDataType,nBandCount::Integer,panBandCount,papszOptions)
#     ccall((:GDALDatasetAdviseRead,libgdal),CPLErr,(Ptr{GDALDatasetH},Cint,Cint,Cint,Cint,Cint,Cint,GDALDataType,Cint,Ptr{Cint},Ptr{Cstring}),hDS,nDSXOff,nDSYOff,nDSXSize,nDSYSize,nBXSize,nBYSize,eBDataType,nBandCount,panBandCount,papszOptions)
# end

# """
#     GDALDatasetGetStyleTable(GDALDatasetH hDS) -> OGRStyleTableH
# Returns dataset style table.
# ### Parameters
# * **hDS**: the dataset handle
# ### Returns
# handle to a style table which should not be modified or freed by the caller.
# """
# function datasetgetstyletable{T <: GDALDatasetH}(arg1::Ptr{T})
#     checknull(ccall((:GDALDatasetGetStyleTable,libgdal),Ptr{OGRStyleTableH},(Ptr{GDALDatasetH},),arg1))
# end

# """
#     GDALDatasetSetStyleTableDirectly(GDALDatasetH hDS,
#                                      OGRStyleTableH hStyleTable) -> void
# Set dataset style table.
# ### Parameters
# * **hDS**: the dataset handle
# * **hStyleTable**: style table handle to set
# """
# function datasetsetstyletabledirectly{T <: GDALDatasetH}(arg1::Ptr{T},arg2::Ptr{OGRStyleTableH})
#     ccall((:GDALDatasetSetStyleTableDirectly,libgdal),Void,(Ptr{GDALDatasetH},Ptr{OGRStyleTableH}),arg1,arg2)
# end

# """
#     GDALDatasetSetStyleTable(GDALDatasetH hDS,
#                              OGRStyleTableH hStyleTable) -> void
# Set dataset style table.
# ### Parameters
# * **hDS**: the dataset handle
# * **hStyleTable**: style table handle to set
# """
# function datasetsetstyletable{T <: GDALDatasetH}(arg1::Ptr{T},arg2::Ptr{OGRStyleTableH})
#     ccall((:GDALDatasetSetStyleTable,libgdal),Void,(Ptr{GDALDatasetH},Ptr{OGRStyleTableH}),arg1,arg2)
# end

# """
#     GDALDatasetStartTransaction(GDALDatasetH hDS,
#                                 int bForce) -> OGRErr
# For datasources which support transactions, StartTransaction creates a transaction.
# ### Parameters
# * **hDS**: the dataset handle.
# * **bForce**: can be set to TRUE if an emulation, possibly slow, of a transaction mechanism is acceptable.
# ### Returns
# OGRERR_NONE on success.
# """
# function datasetstarttransaction{T <: GDALDatasetH}(hDS::Ptr{T},bForce::Integer)
#     ccall((:GDALDatasetStartTransaction,libgdal),OGRErr,(Ptr{GDALDatasetH},Cint),hDS,bForce)
# end

# """
#     GDALDatasetCommitTransaction(GDALDatasetH hDS) -> OGRErr
# For datasources which support transactions, CommitTransaction commits a transaction.
# ### Returns
# OGRERR_NONE on success.
# """
# function datasetcommittransaction{T <: GDALDatasetH}(hDS::Ptr{T})
#     ccall((:GDALDatasetCommitTransaction,libgdal),OGRErr,(Ptr{GDALDatasetH},),hDS)
# end

# """
#     GDALDatasetRollbackTransaction(GDALDatasetH hDS) -> OGRErr
# For datasources which support transactions, RollbackTransaction will roll back a datasource to its state before the start of the current transaction.
# ### Returns
# OGRERR_NONE on success.
# """
# function datasetrollbacktransaction{T <: GDALDatasetH}(hDS::Ptr{T})
#     ccall((:GDALDatasetRollbackTransaction,libgdal),OGRErr,(Ptr{GDALDatasetH},),hDS)
# end

# """
#     GDALDatasetTestCapability(GDALDatasetH hDS,
#                               const char * pszCap) -> int
# Test if capability is available.
# ### Parameters
# * **hDS**: the dataset handle.
# * **pszCap**: the capability to test.
# ### Returns
# TRUE if capability available otherwise FALSE.
# """
# function datasettestcapability{T <: GDALDatasetH}(arg1::Ptr{T},arg2)
#     ccall((:GDALDatasetTestCapability,libgdal),Cint,(Ptr{GDALDatasetH},Cstring),arg1,arg2)
# end


text:
using StaticCompiler
using StaticTools

function rand_matrix(argc::Int, argv::Ptr{Ptr{UInt8}})
    argc == 3 || return printf(stderrp(), c"Incorrect number of command-line arguments\n")
    rows = parse(Int64, argv, 2)            # First command-line argument
    cols = parse(Int64, argv, 3)            # Second command-line argument

    M = MallocArray{Float64}(undef, rows, cols)
    rng = static_rng()
    @inbounds for i=1:rows
        for j=1:cols
            M[i,j] = rand(rng)
        end
    end
    printf(M)
    free(M)
end

# Attempt to compile
path = compile_executable(rand_matrix, (Int64, Ptr{Ptr{UInt8}}), "./")


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

"""
SoftGlobalScope is a package that simplifies the [variable scoping rules](https://docs.julialang.org/en/stable/manual/variables-and-scoping/)
for code in *global* scope.   It is intended for interactive shells (the REPL, [IJulia](https://github.com/JuliaLang/IJulia.jl),
etcetera) to make it easier to work interactively with Julia, especially for beginners.

In particular, SoftGlobalScope provides a function `softscope` that can transform Julia code from using the default
"hard" (local) scoping rules to simpler "soft" scoping rules in global scope, and a function `softscope_include_string`
that can evaluate a whole string (similar to `include_string`) using these rules.

For example, if `s` is a global variable in the current module (e.g. `Main`), the following code is an error in
Julia 1.0:
```
for i = 1:10
    s += i     # declares a new local variable `s`
end
```
Instead, you can transform the expression using `softscope` to automatically insert the necessary `global` keyword:
```jl
julia> softscope(Main, :(for i = 1:10
           s += i
       end))
:(for i = 1:10
      #= REPL[3]:2 =#
      global s += i
  end)
```
You can then execute the statement with `eval`. Alternatively, you can decorate
the expression with the `@softscope` macro:
```jl
julia> s = 0;

julia> @softscope for i = 1:10
           s += i
       end

julia> s
55
```
This macro should only be used in the global scope (e.g., via the REPL); using
this macro within a function is likely to lead to unintended consequences.

You can execute an entire sequence of statements using "soft" global scoping
rules via `softscope_include_string`:
```jl
julia> softscope_include_string(Main, \"\"\"
       s = 0
       for i = 1:10
           s += i
       end
       s
       \"\"\")
55
```
(This function works like `include_string`, returning the value of the last evaluated expression.)

On Julia 0.6, `softscope` is the identity and `softscope_include_string` is equivalent to
`include_string`, since the `global` keyword is not needed there.
"""
module SoftGlobalScope
export softscope, softscope_include_string, @softscope

if VERSION < v"0.7.0-DEV.2308" # before julia#19324 we don't need to change the ast
    softscope(m::Module, ast) = ast
    softscope_include_string(m::Module, code::AbstractString, filename::AbstractString="string") =
        @static isdefined(Base, Symbol("@__MODULE__")) ? include_string(m, code, filename) : Core.eval(m, :(include_string($code, $filename)))
else
    using Base.Meta: isexpr

    const assignments = Set((:(=), :(+=), :(-=), :(*=), :(/=), :(//=), :(\=), :(^=), :(÷=), :(%=), :(<<=), :(>>=), :(>>>=), :(|=), :(&=), :(⊻=), :($=)))
    const calls = Set((:call, :comparison, :(&&), :(||), :ref, :tuple))

    # extract the local variable names (e.g. `[:x]`) from assignments (e.g. `x=1`) etc.
    function localvars(ex::Expr)
        if isexpr(ex, :(=)) || isexpr(ex, :(::))
            return localvars(ex.args[1])
        elseif isexpr(ex, :tuple) || isexpr(ex, :block)
            return localvars(ex.args)
        else
            return Any[]
        end
    end
    localvars(ex::Symbol) = [ex]
    localvars(ex) = Any[]
    localvars(a::Vector) = vcat(localvars.(a)...)

    # Deal with assignments where the LHS is always local but the RHS might require global statements
    # For example, with previously defined a, let a = (a = 1) ; end -> let a = (global a = 1) ; end
    function localassignment(ex, globals, locals, insertglobal)
        if isexpr(ex, :block)
            args = []
            for arg in ex.args
                if isexpr(arg, :(=))
                    push!(args, Expr(arg.head, arg.args[1], _softscope(arg.args[2], globals, locals, insertglobal)))
                    union!(locals, localvars(arg.args[1]))
                elseif arg isa Symbol
                    push!(args, arg)
                else
                    error("Unknown syntax - please file an issue in the SoftGlobalScope.jl repository")
                end
            end
            return Expr(ex.head, args...)
        elseif isexpr(ex, :(=))
            return Expr(ex.head, ex.args[1], _softscope(ex.args[2], globals, locals, insertglobal))
        elseif ex isa Symbol
            return ex
        else
            error("Unknown syntax - please file an issue in the SoftGlobalScope.jl repository")
        end
    end

    """
        _softscope(ex, globals, locals, insertglobal::Bool=false)

    Transform expression `ex` to "soft" scoping rules, where `globals` is a collection
    (e.g. `Set`) of global-variable symbols to implicitly qualify with `global`, and
    `insertglobal` is whether to insert the `global` keyword at the top level of
    `ex`.  (Usually, you pass `insertglobal=false` to start with and then it is
    recursively set to `true` for local scopes introduced by `for` etcetera.)
    NOTE: `_softscope`` may mutate the `globals` argument (if there are `local` declarations.)
    """
    function _softscope(ex::Expr, globals, locals, insertglobal::Bool=false, noassignment::Bool=false)
        if isexpr(ex, :for)
            return Expr(ex.head, localassignment(ex.args[1], copy(globals), copy(locals), insertglobal),
                _softscope(ex.args[2], copy(globals), copy(locals), true))
        elseif isexpr(ex, :while)
            return Expr(ex.head, _softscope(ex.args[1], copy(globals), copy(locals), insertglobal),
                _softscope(ex.args[2], copy(globals), copy(locals), true))
        elseif isexpr(ex, :try)
            try_clause = _softscope(ex.args[1], copy(globals), copy(locals), true)
            catch_clause = _softscope(ex.args[3], copy(globals), ex.args[2] isa Symbol ? union!(locals, ex.args[2:2]) : copy(locals), true)
            if length(ex.args) == 3
                return Expr(:try, try_clause, ex.args[2], catch_clause)
            else
                finally_clause = _softscope(ex.args[4], copy(globals), copy(locals), true)
                return Expr(:try, try_clause, ex.args[2], catch_clause, finally_clause)
            end
        elseif isexpr(ex, :let)
            letlocals = union(locals, localvars(ex.args[1]))
            return Expr(ex.head, localassignment(ex.args[1], copy(globals), copy(locals), true),
                _softscope(ex.args[2], copy(globals), letlocals, true))
        elseif isexpr(ex, :block) || isexpr(ex, :if) || isexpr(ex, :elseif) || isexpr(ex, :toplevel)
            return Expr(ex.head, _softscope.(ex.args, Ref(globals), Ref(locals), insertglobal)...)
        elseif isexpr(ex, :global)
            union!(globals, localvars(ex.args))
            return ex
        elseif isexpr(ex, :local)
            union!(locals, localvars(ex.args)) # affects globals in surrounding scope!
            return ex
        elseif ex.head in calls
            return Expr(ex.head, _softscope.(ex.args, Ref(globals), Ref(locals), insertglobal, ex.head === :tuple)...)
        elseif isexpr(ex, :kw) || (noassignment && ex.head in assignments)
            return Expr(ex.head, ex.args[1], _softscope(ex.args[2], globals, locals, insertglobal))
        elseif insertglobal && ex.head in assignments
            if isexpr(ex.args[1], :call)
                return ex
            elseif ex.args[1] in globals && !(ex.args[1] in locals) # Simple assignment to global
                return Expr(:global, Expr(ex.head, ex.args[1], _softscope(ex.args[2], globals, locals, insertglobal)))
            end
            softex = Expr(ex.head, _softscope.(ex.args, Ref(globals), Ref(locals), insertglobal)...)
            if isexpr(ex.args[1], :tuple) # Assignment to a tuple
                vars = [var for var in localvars(ex.args[1].args) if (var in globals) && !(var in locals)]
                return isempty(vars) ? softex : Expr(:block, Expr(:global, vars...), softex)
            else
                return softex
            end
        elseif !insertglobal && isexpr(ex, :(=)) # only assignments in the global scope need to be considered
            union!(globals, localvars(ex))
            return ex
        else
            return ex
        end
    end
    _softscope(ex, globals, locals, insertglobal::Bool=false, noassignment::Bool=false) = ex

    softscope(m::Module, ast) = _softscope(ast, Set(@static VERSION < v"0.7.0-DEV.3526" ? names(m, true) : names(m, all=true)), Set{Symbol}())

    # we want to add line numbers to most expressions, but we can
    # only do this by wrapping them in :block (:toplevel doesn't work),
    # and some expressions are toplevel-only.  We also need to shift
    # any existing line numbers by line-1.
    const toplevel_only = (:module, :primitive, :abstract, :struct)
    _add_linenum(ex, line, filesym) = Expr(:block, LineNumberNode(line, filesym), ex)
    add_linenum(ex, line, filesym) = _add_linenum(ex, line, filesym)
    shift_linenum(ex::LineNumberNode, line, filesym) = ex.file === :none ? LineNumberNode(ex.line+line-1, filesym) : ex
    shift_linenum(ex::Expr, line, filesym) = Expr(ex.head, shift_linenum.(ex.args, line, filesym)...)
    shift_linenum(ex, line, filesym) = ex
    add_linenum(ex::LineNumberNode, line, filesym) = shift_linenum(ex, line, filesym)
    function add_linenum(ex::Expr, line, filesym)
        if ex.head == :toplevel
            return Expr(:toplevel, add_linenum.(ex.args, line, filesym)...)
        end
        exshift = shift_linenum(ex, line, filesym)
        if exshift.head in toplevel_only
            return exshift
        else
            return _add_linenum(exshift, line, filesym)
        end
    end

    function softscope_include_string(m::Module, code::AbstractString, filename::AbstractString="string")
        # read through the code line by line, keeping count to preserve line-number information
        pos = 1 # current position in the code
        lastpos = lastindex(code)
        line = 1
        retval = nothing
        filesym = Symbol(filename) # LineNumberNode needs a Symbol
        while pos ≤ lastpos
            startpos = pos
            e, pos = Meta.parse(code, startpos, greedy=true, raise=false)
            if isexpr(e, :incomplete) || isexpr(e, :error)
                throw(LoadError(filename, line,
                                ErrorException("syntax: " * e.args[1])))
            end
            if e !== nothing # ignore blank and comment lines (#12)
                e = softscope(m, add_linenum(e, line, filesym))
                retval = Core.eval(m, e)
            end
            line += count(==('\n'), SubString(code, startpos, prevind(code, pos)))
        end
        return retval
    end
end

"""
    softscope(m::Module, ast)

Transform the abstract syntax tree `ast` (a quoted Julia expression) to use "soft"
scoping rules for the global variables defined in `m`, returning the new expression.
"""
softscope

if VERSION < v"0.7.0-DEV.481" # the version that __module__ was introduced
    macro softscope(ast)
        esc(softscope(current_module(), ast))
    end
else
    macro softscope(ast)
        esc(softscope(__module__, ast))
    end
end

"""
    @softscope(expr)

Apply "soft" scoping rules to the argument of the macro. For example
```jl
julia> s = 0;

julia> @softscope for i = 1:10
           s += i
       end

julia> s
55
```
"""
:(@softscope)

"""
    softscope_include_string(m::Module, code::AbstractString, filename::AbstractString="string")

Like [`include_string`](@ref), but evaluates `code` using "soft"
scoping rules for the global variables defined in `m`.
"""
softscope_include_string

end # module


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

if !isdefined(@__MODULE__, :DefsWorkspaceId)
    const DefsWorkspaceId = String
else
    @warn("Skipping redefinition of DefsWorkspaceId to String")
end


text:
using FooBar
GoodCounter()
for j = 1:7
	jj = j * 10
	@show(jj)
	uu = 6
	archie = 75
	@show(archie)
end

max_stars_repo_name: ArchieCall/FooBar

max_issues_repo_name: ArchieCall/FooBar

max_forks_repo_name: ArchieCall/FooBar