text:
#TODO: add "expand" method to go from shorter to greater frequencies
#TODO: include functionality for filling the missing values (approx, locf, etc.)

# Date vector
mondays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Monday
tuesdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Tuesday
wednesdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Wednesday
thursdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Thursday
fridays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Friday
saturdays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Saturday
sundays(t::AbstractArray{Date,1}) = dayofweek.(t) .== Sunday
# DateTime vector
mondays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Monday
tuesdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Tuesday
wednesdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Wednesday
thursdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Thursday
fridays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Friday
saturdays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Saturday
sundays(t::AbstractArray{DateTime,1}) = dayofweek.(t) .== Sunday
# Time series input
mondays(x::TS) = x[mondays(x.index)]
tuesdays(x::TS) = x[tuesdays(x.index)]
wednesdays(x::TS) = x[wednesdays(x.index)]
thursdays(x::TS) = x[thursdays(x.index)]
fridays(x::TS) = x[fridays(x.index)]
saturdays(x::TS) = x[saturdays(x.index)]
sundays(x::TS) = x[sundays(x.index)]

# If `cal` is true, only returns observations on the last calendar day of the month.
# If `cal` is false, returns all observations for which the next index is a new month.
bow(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofweek.(t)) : BitVector([0;diff(week.(t)).!=0])
eow(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofweek.(t)) : BitVector([diff(week.(t)).!=0;0])
bom(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofmonth.(t)) : BitVector([0;diff(month.(t)).!=0])
eom(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofmonth.(t)) : BitVector([diff(month.(t)).!=0;0])
boq(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofquarter.(t)) : BitVector([0;diff(quarterofyear.(t)).!=0])
eoq(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofquarter(t)) : BitVector([diff(quarterofyear.(t)).!=0;0])
boy(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofyear.(t)) : BitVector([0;diff(year.(t)).!=0])
eoy(t::AbstractArray{Date,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofyear.(t)) : BitVector([diff(year.(t)).!=0;0])
bow(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofweek.(t)) : BitVector([0;diff(week.(t)).!=0])
eow(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofweek.(t)) : BitVector([diff(week.(t)).!=0;0])
bom(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofmonth.(t)) : BitVector([0;diff(month.(t)).!=0])
eom(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofmonth.(t)) : BitVector([diff(month.(t)).!=0;0])
boq(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofquarter.(t)) : BitVector([0;diff(quarterofyear.(t)).!=0])
eoq(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofquarter(t)) : BitVector([diff(quarterofyear.(t)).!=0;0])
boy(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== firstdayofyear.(t)) : BitVector([0;diff(year.(t)).!=0])
eoy(t::AbstractArray{DateTime,1}; cal::Bool=false) = cal ? BitVector(t .== lastdayofyear.(t)) : BitVector([diff(year.(t)).!=0;0])

bow(x::TS; cal::Bool=false) = x[bow(x.index, cal=cal)]
eow(x::TS; cal::Bool=false) = x[eow(x.index, cal=cal)]
bom(x::TS; cal::Bool=false) = x[bom(x.index, cal=cal)]
eom(x::TS; cal::Bool=false) = x[eom(x.index, cal=cal)]
boq(x::TS; cal::Bool=false) = x[boq(x.index, cal=cal)]
eoq(x::TS; cal::Bool=false) = x[eoq(x.index, cal=cal)]
boy(x::TS; cal::Bool=false) = x[boy(x.index, cal=cal)]
eoy(x::TS; cal::Bool=false) = x[eoy(x.index, cal=cal)]

function collapse(x::TS{V,T}, at::AbstractArray{Bool,1}; fun::Function=last, args...)::TS{V,T} where {V,T}
    @assert size(at,1) == size(x,1) "Arguments `x` and `at` must have same number of rows."
    idx = findall(at)
    if fun == last
        return x[idx]
    end
    out = zeros(V, (length(idx), size(x,2)))
    @inbounds for j in 1:size(x,2)
        out[1,j] = fun(x.values[1:idx[1], j]; args...)
    end
    @inbounds for j=1:size(x,2), i=2:length(idx)
        a = idx[i-1] + 1
        b = idx[i]
        out[i,j] = fun(x.values[a:b, j]; args...)
    end
    return ts(out, x.index[at], x.fields)
end

# Here `at` is a function of the index of `x` which generates a BitArray to be used as above
function collapse(x::TS{V,T}, at::Function; fun::Function=last, args...)::TS{V,T} where {V,T}
    return collapse(x, at(x.index), fun=fun; args...)
end

function apply(x::TS{V}, dim::Int=1; fun=sum) where {V}
    if dim == 1
        n = size(x,1)
        out = zeros(V, n)
        @inbounds for i in 1:n
            out[i] = fun(x.values[i,:])
        end
        return ts(out, x.index, Symbol(uppercasefirst(string(fun))))
    elseif dim == 2
        k = size(x,2)
        out = zeros(V, k)
        @inbounds for j in 1:k
            out[j] = fun(x.values[:,j])
        end
        return out'
    else
        error("Argument `dim` must be either 1 (rows) or 2 (columns).")
    end
end


text:
using Statistics,QuadGK

using Optim,Roots

mutable struct  EuropeanOption
    St::Float64 #current assset price 当期价格
    K::Float64 #strike price 行权价
    T::Float64 #maturity time 到期时间
    t::Float64 #valuation date 有效日期
    r::Float64 #risk_less rate 无风险收益率
    σ::Float64 #volatility 波动率
    type::String #option type(call,put) 看涨、看跌
    # EuropeanOtion(St,K,T,t,r,σ,type)=new(St,K,T,t,r,σ,type)
end

#正态分布，数值累计积分算法
function CDF_dn1(x)
    a1 = 0.31938153
    a2 = -0.356563782
    a3 = 1.781477937
    a4 = -1.821255978
    a5 = 1.330274429
    l = abs(x)
    k = 1 / (1 + 0.2316419 * l)

    CND = 1 - 1 / sqrt(2 * π) * exp(-l^2 / 2) * (a1 * k + a2 * k^2 + a3 * k^3 + a4 * k^4 + a5 * k^5)

    if x < 0
        return 1 - CND
    end

    return CND
end

#BSM公式定价
#dN(x)标准正太函数显示表达（即高斯函数）
#CDF_dn(x)高斯函数积分
function BSM_value(o::EuropeanOption,greeks::Bool)
    dN(x)=exp(-0.5*x^2)/sqrt(2*pi)
    CDF_dn(x)=quadgk((x)->dN(x),-20,x)[1]

    d1=(log(o.St/o.K)+(o.r+o.σ^2/2)*o.T)/(o.σ*√o.T)
    d2=d1-o.σ*√o.T

    option_value=0.0

    if o.type=="put"
        option_value=CDF_dn(-d2)*o.K*exp(-o.r*o.T)-CDF_dn(-d1)*o.St
    else
        option_value=o.St*CDF_dn(d1)-exp(-o.r*o.T)*o.K*CDF_dn(d2)
    end

    if greeks
        vega=o.St*dN(d1)*√o.T
        δ=CDF_dn(d1)
        γ=dN(d1)/(o.St*o.σ*√o.T)
        ρ=o.K*o.T*exp(-o.r*o.T)*CDF_dn(d2)
        θ=-(o.St*dN(d1)*o.σ/(2*√o.T)+o.r*o.K*exp(-o.r*o.T)*CDF_dn(d2))
        return [option_value,vega,δ,γ,ρ,θ]
    else
        return [option_value]
    end

end
#
#

#格式化函数为，find_zero提供标准输入
function formatfun(o::EuropeanOption,x)
    o.σ=x
    z=BSM_value(option)[1]
end

#隐含波动率计算
function imp_vol(o::EuropeanOption,x,c0)
    fi(x)=formatfun(o,x)-c0
    z=find_zero(fi,(0,1),Bisection())
    return z
end

# option=EuropeanOption(100.0,100.0,1,0,0.05,0.2,"call")
# option_values=BSM_value(option,false)
#
# vega=option_values[2]

ext: jl


text:
"""
    line(x, y)

Creates an `EChart` where region between points connected by a line.

## Methods
```julia
line(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}})
line(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2})
line(df::AbstractDataFrame, x::Symbol, y::Symbol)
line(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol)
```

## Arguments
* `mark::Union{String, AbstractVector} = "line"` : how to display plotted points
* `step::Union{String, Void} = nothing` : one of {"start", "end", "middle", nothing}
* `legend::Bool = false` : display legend?
* `scale::Bool = false` : show full Y-axis or truncated
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `area`, such as displaying a legend when two or more series present.

## Examples
```julia
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ll = line(x, y)
```
"""
function line(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}};
			mark::Union{String, AbstractVector} = "line",
			step::Union{String, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

	return xy_plot(x, y; mark = mark, step = step, legend = legend, scale = scale, kwargs...)

end

#array, multiple series
function line(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2};
			mark::Union{String, AbstractVector} = "line",
			step::Union{String, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

	return xy_plot(x, y; mark = mark, step = step, legend = legend, scale = scale, kwargs...)

end

#DataFrame, single series
function line(df::AbstractDataFrame, x::Symbol, y::Symbol;
			mark::Union{String, AbstractVector} = "line",
			step::Union{String, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

	return xy_plot(df, x, y; mark = mark, step = step, legend = legend, scale = scale, kwargs...)

end

#DataFrame with group argument
function line(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol;
			mark::Union{String, AbstractVector} = "line",
			step::Union{String, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

	return xy_plot(df, x, y, group; mark = mark, step = step, legend = legend, scale = scale, kwargs...)

end


text:
struct Fσ end
struct FGσ end

# type for doing a SE expt
mutable struct ParamExpt{T}
#	paTD::Data.P_misfit #  to measure data misfit
#	acqgeom::Acquisition.Geom # mainly positions of receivers are used from here, for the source see P
	σ::Matrix{T} # see Niels et al.
	Q::Matrix{T} # Q*k/η (see Niels et al.)
	P::Array{T,3} # pressure snapshots (B*P) from seismic modeling
	LP::Array{T,3} # ∇⋅([Q*k/η](x,z)∇B*P)
	data::Array{T,2}
	data_obs::Array{T,2}
	data_misfit::Array{T,1}
#	datP::Data.TD
#	datLP::Data.TD
	ψ::Vector{T} # state variable, we are going to record this state data=ACQ*ψ
	ψ0::Vector{T} # background state variable, only used for born modelling
	adjsrc::Vector{T}
	adjψ::Vector{T} # adjoint state variable
	adjψ2::Vector{T} # some storage of size adjψ
	paQ::Param{T} # Param for each snapshot
	paσ::Param{T} # Param for each snapshot
	paσ0::Param{T} # Param for each snapshot, only used for Born modelling
	tgrid::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} # tgrid of snapshots
	mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}} # model grid of snapshots
	ACQ::SparseMatrixCSC{T,Int64}
	g::Matrix{T}
	gtemp::Vector{T}
end

"""
some dummy ACQ, use `GeoPhyInv.Acquisition.ACQmat` later
"""
function testACQ(nz,nx,nr)

	ACQ=spzeros(nr,nz*nx)
	for i in 1:2:nr*nz*nx
		ACQ[i]=randn()
	end
	return ACQ
end

# allocate and forward modelling in true model
function ParamExpt(snaps, tgrid, mgrid,  Qv, k, η, σ, ACQmat=nothing; σobs=nothing, Qobs=nothing)
	nz=length(mgrid[1])
	nx=length(mgrid[2])
	nt=length(tgrid)
	(size(snaps)≠(nz,nx,nt)) && error("tgrid, mgrid and snaps inconsistent")
	(size(Qv)≠(nz,nx)) && error("tgrid, mgrid and Qv inconsistent")
	(size(k)≠(nz,nx)) && error("tgrid, mgrid and k inconsistent")
	(size(η)≠(nz,nx)) && error("tgrid, mgrid and η inconsistent")
	(size(σ)≠(nz,nx)) && error("tgrid, mgrid and σ inconsistent")

	LP=zero(snaps);
	ψ=zeros(nz*nx);
	ψ0=zeros(nz*nx);

	adjsrc=zeros(nz*nx);
	adjψ=zeros(nz*nx+1);
	adjψ2=zeros(nz*nx+1);

	Q= k .* Qv ./ η # combine all these

#	dobs=Data.TD_zeros([:P],tgrid,acqgeom)
	#Random.randn!(dobs) # dummy dobs, update later
#	paTD=Data.P_misfit(Data.TD_zeros([:P],tgrid,acqgeom),dobs);

	paQ=Param(mgrid, Q)
	paσ=Param(mgrid, σ)
	paσ0=Param(mgrid, σ)

	if(ACQmat===nothing)
		ACQ=testACQ(nz,nx,5)
	else
		ACQ=ACQmat
	end
	data=zeros(nt,size(ACQ,1))
	data_obs=randn(nt,size(ACQ,1))
	data_misfit=randn(size(ACQ,1))

	pa=ParamExpt(#paTD,
	      #acqgeom,
	      σ, Q, snaps,  LP, #
	      data,
	      data_obs,
	      data_misfit,
	      ψ,ψ0,
	      adjsrc,
	      adjψ,
	      adjψ2,
	     # ssnaps,
	     # data,
	     #datP, datLP,
	     #ψ,
	     paQ, paσ, paσ0,
	     tgrid, mgrid,
	     ACQ,
	    zero(σ),zeros(length(σ)))

	# record P and LP
#	record!(pa.datP, pa.P, pa.mgrid)
#	record!(pa.datLP, pa.LP, pa.mgrid)

	#
	if(!(σobs===nothing) && !(Qobs===nothing))
		updateLP!(pa, Qobs)
		update_observed_data!(pa, σobs)
	end
	return pa
end

#=
# given snapshots, record a data obj TD
function record!(data::Data.TD, snaps, mgrid)
	ifield=1;
	iss=1;
	nr=data.acqgeom.nr
	dat=data.d[iss,ifield]
	for ir = 1:nr[iss]
		din=view(dat,:,ir)
		irx=argmin(abs.(mgrid[2].-data.acqgeom.rx[iss][ir]))
		irz=argmin(abs.(mgrid[1].-data.acqgeom.rz[iss][ir]))
		for it in 1:length(data.tgrid)
			din[it]=snaps[irz,irx,it]
		end
	end
end
=#

# generate LP from P, need to optimize this routine, has a loop over snapshots
function updateLP!(pa::ParamExpt, Q=pa.Q)
	updateA!(pa.paQ, Q)
	@showprogress 10 "time loop LP\t" for it in 1:length(pa.tgrid)
		snap_in=view(pa.P,:,:,it)
		snap_out=view(pa.LP,:,:,it)
		applyA!(snap_out, snap_in, pa.paQ; A=pa.paQ.A)
	end
	return nothing
end

function mod!(pa::ParamExpt)
	updateLP!(pa)
	mod!(pa, pa.σ, Fσ())
end

# forward modeling to generate ψ
# update psi using LP, has a loop over snapshots
# most expensive routine
function mod!(pa::ParamExpt, σ, mode)
	updateA!(pa.paσ, σ)
	qrA=factorize(pa.paσ.A)
	qrAt=factorize(pa.paσ.At)
	fill!(pa.g, 0.0)
	@showprogress 10 "time loop ψ\t" for it in 1:length(pa.tgrid)
		snap_in=view(pa.LP,:,:,it)
		forwψ!(pa.ψ,snap_in,pa.paσ,qrA)

		# record
		dat_slice=view(pa.data,it,:)
		mul!(dat_slice,pa.ACQ,pa.ψ) # record

		# wanna do adjoint modeling? depends on mode
		dat_slice_obs=view(pa.data_obs,it,:)
		adjψ!(pa.adjψ, pa.adjψ2, pa.ψ,
		      pa.adjsrc, pa.g, pa.gtemp,
		      pa.ACQ, pa.data_misfit, dat_slice_obs, pa.paσ, qrAt, mode)
	end
	return nothing
end

function forwψ!(ψ, src, pa::Param,qrA)
	applyinvA!(ψ, src, pa; A=qrA)
end
function adjψ!(adjψ, adjψ2, ψ, adjsrc, g, gtemp, ACQ, data_misfit, data_obs, pa, qrAt,   ::Fσ)
end
function adjψ!(adjψ, adjψ2, ψ, adjsrc, g, gtemp, ACQ, data_misfit, data_obs, pa, qrAt,  ::FGσ)

	mul!(data_misfit,ACQ,ψ)
	for i in eachindex(data_misfit)
		data_misfit[i] -= data_obs[i]
	end
	rmul!(data_misfit,-2.0)

	adjψ_core!(adjψ, adjψ2, ψ, adjsrc, g, gtemp, ACQ, data_misfit, pa, qrAt)

end
function adjψ_core!(adjψ, adjψ2, ψ, adjsrc, g, gtemp, ACQ, data_misfit, pa, qrAt)

	mul!(adjsrc, transpose(ACQ), data_misfit)

	applyinvAt!(adjψ, adjsrc, pa, At=qrAt)

	abT_T_dX!(gtemp, ψ, adjψ, adjψ2, pa.dAdx)
	#= requires a lot of memory (use only for testing)
	mul!(ψψ, adjψ, transpose(ψ))
	copyto!(ψψvec,ψψ)
	mul!(gtemp, pa.dAdx', ψψvec)
	=#

	for i in eachindex(gtemp)
		g[i]+=gtemp[i]
	end

end

# imaging condition
function abT_T_dX!(g, ψ, adjψ, adjψ2, A)
	n=length(ψ)
	fill!(g, 0.0)
	for i in 1:n
		AA=A[i]
		mul!(adjψ2,AA,ψ)
		g[i] = dot(adjψ,adjψ2)
		# do everything by hand, was slow...
		#for k in 1:n
		#	for j in 1:n+1
		#		g[i] += ψ[k]*adjψ[j]*A[j+(k-1)*(n+1),i]
		#	end
		#end
	end
end

function func(x,pa::ParamExpt)
	mod!(pa, x, Fσ())
	f=sum(abs2,pa.data.-pa.data_obs)
	return f
end

function grad(x,pa::ParamExpt)
	mod!(pa, x, FGσ())
	return paE.g
end

# do forward modeling to generate observed data in paTD
function update_observed_data!(pa::ParamExpt, σobs)
	# save sigma
	σsave=copy(pa.σ)
	mod!(pa, σobs, Fσ())
	copyto!(pa.data_obs, pa.data)
	fill!(pa.data, 0.0) # clear data
	# put sigma back
	copyto!(pa.σ, σsave)
	# dummy modeling, just to remove traces of σobs
	mod!(pa, σobs, Fσ())
end


text:
dog, Dog, DOG = "Benjamin", "Samba", "Bernie"

if dog === Dog
    println("There is only one dog, ", DOG)
else
    println("The three dogs are:  ", dog, ", ", Dog, " and ", DOG)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module PlutoPublish

import Pluto, SHA

export publisher

export_cache = joinpath(homedir(), ".julia", "pluto_publish")
if isdir(export_cache)
    # this is a temporary _cache_ directory, so remove all the files that were there
    for to_delete ∈ filter(x -> endswith(x, ".statefile"), readdir(export_cache))
        rm(joinpath(export_cache, to_delete))
    end
else
    mkdir(export_cache)
end

include("client.jl")
include("upload.jl")

function _save_statefile(notebook::Pluto.Notebook, statefile_data::Vector{UInt8}, publisher_url::String)
    # TODO: implement this in WYSIWYR instead when that finally gets released
    static_notebook_hash = bytes2hex(SHA.sha256(notebook.path))
    save_path = joinpath(export_cache, static_notebook_hash * ".statefile")

    decoded_statefile = Pluto.unpack(statefile_data)
    publish_configs = map(
        filter(collect(decoded_statefile["cell_inputs"])) do input_cell
            contains(input_cell.second["code"], "PUBLISH")
        end
    ) do input_cell
        Meta.parse(input_cell.second["code"])
    end
    if length(publish_configs) > 1
        @warn "More than one publishing configuration was included in $(basename(notebook.path)). The first will be used"
    end
    if length(publish_configs) > 0
        publish_configs = first(publish_configs)

        @assert typeof(publish_configs) == Expr "PUBLISH must be assigned to either `true` or `false`"

        # PLUTO_LIVE_PUBLISH_CONFIG = (; args...)   <-- optional ;
        publishing_enabled = if first(publish_configs.args) == :PUBLISH  # without ;
            publish_configs.args[2]  # unsafe, but its client-side from the client
        elseif publish_configs.head == :toplevel && first(first(publish_configs.args).args) == :PUBLISH  # with ;
            publish_configs.args[1].args[2]
        else
            nothing
        end

        @assert typeof(publishing_enabled) == Bool "PUBLISH must be assigned to either `true` or `false`"

        if publishing_enabled
            # actually publish the notebook
            open(save_path, "w") do io
                write(io, statefile_data)
            end

            _publish_notebook(notebook, static_notebook_hash, publisher_url)
        end
    end

    # if we don't have publishing configs then don't publish it
end

_uri(se::Pluto.ServerStartEvent) = parse(HTTP.URI, se.address)
_base_url(se::Pluto.ServerStartEvent) = begin
    uri = _uri(se)
    "$(uri.scheme)://$(uri.host):$(uri.port)"
end
_secret(se::Pluto.ServerStartEvent) = get(HTTP.queryparams(_uri(se)), "secret", "")
_statefile_url(notebook::Pluto.Notebook, start_event::Pluto.ServerStartEvent) = "$(_base_url(start_event))/statefile?id=$(notebook.notebook_id)&secret=$(_secret(start_event))"

# cache the server start event once it's sent
server_start_event = nothing

# configured usage
publisher(destination::String = "https://plutopublish.com") = function(event::Pluto.PlutoEvent)
    global server_start_event
    if typeof(event) == Pluto.FileEditEvent
        try
            response = HTTP.get(_statefile_url(event.notebook, server_start_event))
            _save_statefile(event.notebook, response.body, destination)
        catch
            @warn "Failed to fetch notebook statefile"
            for (exc, bt) in current_exceptions()
                showerror(stdout, exc, bt)
                println(stdout)
            end
        end
    elseif typeof(event) == Pluto.ServerStartEvent
        server_start_event = event
    end
end

end # module


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@doc raw"""
`struct PSH1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\begin{aligned}
\phi_{nlm}({\bm r}, z_1, z_0) &= P_{nl}(r) Y_l^m(\hat{\br r}) \\
P_{nl}(r) &= C_{nl,n'} J_{n'}(r)
\end{aligned}
```
where ``J_{n}`` denotes a radial basis.
"""
struct PSH1pBasis{T, NZ, TJ <: ScalarBasis{T}} <: OneParticleBasis{T}
   J::TJ
   SH::SHBasis{T}
   zlist::SZList{NZ}
   C::SMatrix{NZ, NZ, Matrix{T}}
   spec::SMatrix{NZ, NZ, Vector{PSH1pBasisFcn}}
   nlmap::SMatrix{NZ, NZ, Matrix{Int}}
   Aindices::Matrix{UnitRange{Int}}
end

# --------------------------------------------------
# lots of basic book-keeping

"largest (P^{z z0}_nl)_{nl} sub-basis"
_maxnl(basis::PSH1pBasis) = maximum(size(C, 1) for C in basis.C)

"index into BY (spherical harmonics)"
_idxY(b1p::PSH1pBasis, nlm, iz, iz0) = index_y(nlm.l, nlm.m)

"index into BP (Pnl radial basis)"
_idxR(b1p::PSH1pBasis, nlm, iz, iz0) = b1p.nlmap[iz, iz0][nlm.n, nlm.l+1]

get_basis_spec(basis::PSH1pBasis, z0::AtomicNumber, args...) =
      get_basis_spec(basis, z2i(basis, z0), args...)

get_basis_spec(basis::PSH1pBasis, iz0::Integer) =
      vcat( basis.spec[:, iz0]... )

function get_basis_spec(basis::PSH1pBasis, iz0::Integer, i::Integer)
   idx = 0
   for iz = 1:numz(basis)
      if idx + length(basis.spec[iz, iz0]) >= i
         return basis.spec[iz, iz0][i - idx]
      end
   end
   error("index i = $i not found; probably this means that i > length(basis[iz0])")
end

cutoff(basis::PSH1pBasis) = cutoff(basis.J)

Base.length(basis::PSH1pBasis, z::AtomicNumber, z0::AtomicNumber) =
      length(basis, z2i(basis, z), z2i(basis, z0))

Base.length(basis::PSH1pBasis, iz::Integer, iz0::Integer) =
      length(basis.spec[iz, iz0])

Base.length(basis::PSH1pBasis, z0::AtomicNumber) =
      length(basis, z2i(basis, z0))

Base.length(basis::PSH1pBasis, iz0::Integer) =
      sum(length(basis, iz, iz0) for iz = 1:numz(basis))

# TODO: delete after testing
# i2z(basis::PSH1pBasis, i::Integer) = i2z(basis.zlist, i)
# z2i(basis::PSH1pBasis, z::AtomicNumber) = z2i(basis.zlist, z)
# numz(basis::PSH1pBasis) = length(basis.zlist)

# ----------------------------------------
# basis construction

function PSH1pBasis(J::ScalarBasis, maxn;
                         species = :X,
                         D::AbstractDegree = SparsePSHDegree())
   max_p = length(J)
   # convert species list to an actual list
   if species isa Symbol; species = (species,); end
   # convert maxn parameter to a dictionary if needed
   if maxn isa Number
      maxn = Dict( [ (s1, s2) => maxn for s1 in species, s2 in species ] )
   end
   # produce the specifications for all s1, s2 pairs
   zlist = ZList(species...; static=true)
   NZ = length(zlist)
   spec = Matrix{Vector{PSH1pBasisFcn}}(undef, NZ, NZ)
   C = Matrix{Matrix}(undef, NZ, NZ)
   nlmap = Matrix{Matrix{Int}}(undef, NZ, NZ)
   for iz = 1:NZ, iz0 = 1:NZ
      z, z0 = i2z(zlist, iz), i2z(zlist, iz0)
      maxn_zz0 = maxn[ (chemical_symbol(z), chemical_symbol(z0)) ]
      @assert max_p >= maxn_zz0
      spec[iz, iz0] = _build_PSH_1p_spec(maxn_zz0, D, zlist.list)
      nls = sort( unique( [(n = b.n, l = b.l)  for b in spec[iz, iz0]] ) )
      nlmap[iz, iz0] = _build_nlmap(nls)
      C[iz, iz0] = _init_radial_params(fltype(J), nlmap[iz, iz0], max_p)
   end
   # now we can produce the
   return PSH1pBasis(J, zlist, spec, C, nlmap)
end

function _smatrix(A::AbstractMatrix)
   nrows, ncols = size(A)
   A = identity.(A)
   return SMatrix{nrows, ncols, eltype(A)}(A)
end

function PSH1pBasis(J::ScalarBasis{T}, zlist::SZList, spec::AbstractMatrix,
                    C::AbstractMatrix, nlmap::AbstractMatrix) where {T}
   maxL = maximum( maximum(b.l for b in spec12)  for spec12 in spec )
   SH = SHBasis(maxL, fltype(J))
   NZ = length(zlist)
   C = _smatrix(C)
   spec = _smatrix(spec)
   nlmap = _smatrix(nlmap)
   Aindices = Matrix{UnitRange{Int}}(undef, NZ, NZ)
   P1 = PSH1pBasis(J, SH, zlist, C, spec, nlmap, Aindices)
   set_Aindices!(P1)
   return P1
end

function _build_nlmap(nls)
   maxn = maximum( nl.n for nl in nls )
   maxl = maximum( nl.l for nl in nls )
   nlmap = zeros(Int, (maxn, maxl+1))
   idx = 0
   for nl in nls
      idx += 1
      nlmap[nl.n, nl.l+1] = idx
   end
   return nlmap
end

function _init_radial_params(T, nlmap, maxp)
   # initialise the whole matrix to zeros
   C = zeros(T, maximum(nlmap), maxp)
   # for each l ...
   for l = 1:size(nlmap, 2)
      # ... find the n values for which (n,l) defines an R_nl fcn
      ns = nlmap[:, l]
      ns = sort(ns[ findall(ns .> 0) ])
      # this slice of the matrix gives the parameters for (R_nl)_{n ∈ ns}
      # hence we initialize it to I.
      C[ns, :] = Matrix{T}(I, length(ns), maxp)
   end
   return C
end

# function BasicPSH1pBasis(J::ScalarBasis{T}, zlist::SZList,
#                          spec::Vector{PSH1pBasisFcn}) where {T}
#    # now get the maximum L-degree to generate the SH basis
#    maxL = maximum(b.l for b in spec)
#    SH = SHBasis(maxL, T)
#    # construct a preliminary Aindices array to get an incorrect basis
#    NZ = length(zlist)
#    P = BasicPSH1pBasis(J, SH, zlist, spec,
#                        Matrix{UnitRange{Int}}(undef, NZ, NZ))
#    # ... now fix the Aindices array
#    set_Aindices!(P)
#    return P
# end

# function BasicPSH1pBasis(J::ScalarBasis;
#                          species = :X,
#                          D::AbstractDegree = SparsePSHDegree())
#    # get a generic basis spec
#    spec = _get_PSH_1p_spec(J::ScalarBasis, D::AbstractDegree, species)
#    # construct the basis
#    zlist = ZList(species; static=true)
#    return BasicPSH1pBasis(J, zlist, spec)
# end

# # ------------------------------------------------------
# #  FIO code

# ==(P1::BasicPSH1pBasis, P2::BasicPSH1pBasis) =  ACE1._allfieldsequal(P1, P2)

# write_dict(basis::BasicPSH1pBasis{T}) where {T} = Dict(
#       "__id__" => "ACE1_BasicPSH1pBasis",
#            "J" => write_dict(basis.J),
#           "SH" => write_dict(basis.SH),
#         "spec" => write_dict.(basis.spec),
#        "zlist" => write_dict(basis.zlist),
#    )

# function read_dict(::Val{:ACE1_BasicPSH1pBasis}, D::Dict)
#    J = read_dict(D["J"])
#    SH = read_dict(D["SH"])
#    zlist = read_dict(D["zlist"])
#    spec = read_dict.(D["spec"])
#    P = BasicPSH1pBasis(J, SH, zlist, spec,
#                  Matrix{UnitRange{Int}}(undef, length(zlist), length(zlist)))
#    set_Aindices!(P)
#    return P
# end

# ------------------------------------------------------
#  Evaluation code

fltype(basis::PSH1pBasis{T}) where T = Complex{T}
rfltype(basis::PSH1pBasis{T}) where T = T

alloc_temp(basis::PSH1pBasis, args...) =
      (
        BJ = alloc_B(basis.J, args...),
        BP = zeros(fltype(basis), _maxnl(basis)),
        BY = alloc_B(basis.SH, args...),
        tmpJ = alloc_temp(basis.J, args...),
        tmpY = alloc_temp(basis.SH, args...),
       )

function add_into_A!(A, tmp, basis::PSH1pBasis,
                     R, iz::Integer, iz0::Integer)
   # evaluate the r-basis
   C = basis.C[iz, iz0]  # n' -> nl transformation coefficients
   maxnJ = size(C, 2)    # number of J basis functions (pre-basis)
   evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(R); maxn=maxnJ)
   maxn = size(C,1)      # number of transformed basis functions
   BP = mul!((@view tmp.BP[1:maxn]), C, (@view tmp.BJ[1:maxnJ]))

   # evaluate the R̂-basis for the current neighbour at R
   evaluate!(tmp.BY, tmp.tmpY, basis.SH, R)

   # add the contributions to the A_zklm
   @inbounds for (i, nlm) in enumerate(basis.spec[iz, iz0])
      A[i] += (tmp.BP[_idxR(basis, nlm, iz, iz0)]
                * tmp.BY[_idxY(basis, nlm, iz, iz0)] )
   end
   return nothing
end

# alloc_temp_d(basis::BasicPSH1pBasis, args...) =
#       (
#         BJ = alloc_B(basis.J, args...),
#         tmpJ = alloc_temp(basis.J, args...),
#         BY = alloc_B(basis.SH, args...),
#         tmpY = alloc_temp(basis.SH, args...),
#         dBJ = alloc_dB(basis.J, args...),
#         tmpdJ = alloc_temp_d(basis.J, args...),
#         dBY = alloc_dB(basis.SH, args...),
#         tmpdY = alloc_temp_d(basis.SH, args...),
#        )

# function add_into_A_dA!(A, dA, tmpd, basis::BasicPSH1pBasis, R, iz::Integer, iz0::Integer)
#    r = norm(R)
#    R̂ = R / r
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate_d!(tmpd.BJ, tmpd.dBJ, tmpd.tmpdJ, basis.J, r)
#    evaluate_d!(tmpd.BY, tmpd.dBY, tmpd.tmpdY, basis.SH, R)
#    # add the contributions to the A_zklm, ∇A
#    @inbounds for (i, nlm) in enumerate(basis.spec)
#       iY = index_y(nlm.l, nlm.m)
#       A[i] += tmpd.BJ[nlm.n] * tmpd.BY[iY]
#       dA[i] = (tmpd.dBJ[nlm.n] * tmpd.BY[iY]) * R̂ + tmpd.BJ[nlm.n] * tmpd.dBY[iY]
#    end
#    return nothing
# end


text:
@testset "ModelType: parameter compatibility" begin
    for edge_obj in subtypes(EdgeWiseObjectiveFunction)
        for agg_obj in vcat([subtypes(s) for s in subtypes(AggregationObjectiveFunction)]...) # Assumption: only abstract types directly derive from AggregationObjectiveFunction.
            for type in subtypes(FormulationType)
                for cg in [true, false]
                    for algo in subtypes(AlgorithmChoice)
                        for unc in subtypes(UncertaintyHandling)
                            for uncparams in subtypes(UncertainParameters)
                                @testset "ModelType($edge_obj, $agg_obj, $type, $cg, $algo, $unc, $uncparams)" begin
                                    # Failures.
                                    if uncparams == UncertainCapacity && edge_obj in [Load, KleinrockLoad, FortzThorupLoad]
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                    elseif type == FlowFormulation && cg
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                    elseif (unc != ObliviousUncertainty && unc != RobustUncertainty) && algo != Automatic
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                    elseif unc == NoUncertaintyHandling && uncparams != NoUncertainty
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                    elseif unc == NoUncertaintyHandling && uncparams != NoUncertainty
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                    elseif uncparams == UncertainCapacity
                                        @test_throws ErrorException ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())

                                    # Successes.
                                    elseif edge_obj == AlphaFairness
                                        for α in [0.0, 0.25, 1.0]
                                            for pc in [true, false]
                                                mt = nothing

                                                if pc && α == 0.0
                                                    # Warn when forcing the use of a power cone when there is not available.
                                                    msg = "The power-cone formulation is forced to be used, " *
                                                          "but it is not available for α = $(α). " *
                                                          "The parameter `force_power_cone` is therefore ignored."
                                                    @test_logs (:warn, msg) mt = ModelType(AlphaFairness(α, pc), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                                elseif pc && ! (α in [0.5, 1.5, 2.0])
                                                    # Warn when forcing the use of a power cone when there is nothing else available.
                                                    msg = "The power-cone formulation is forced to be used, " *
                                                          "but it is the only available one for α = $(α). " *
                                                          "The parameter `force_power_cone` is therefore ignored."
                                                    @test_logs (:warn, msg) mt = ModelType(AlphaFairness(α, pc), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                                else
                                                    # Don't warn in other cases.
                                                    @test_nowarn mt = ModelType(AlphaFairness(α, pc), agg_obj(), type(), cg, algo(), unc(), uncparams())
                                                end

                                                @test mt !== nothing
                                                @test typeof(mt.edge_obj) == edge_obj
                                                @test mt.edge_obj.α == α
                                                @test mt.edge_obj.force_power_cone == pc
                                                @test mt.agg_obj == agg_obj()
                                                @test mt.type == type()
                                                @test mt.cg == cg
                                                @test mt.unc == unc()
                                                @test mt.uncparams == uncparams()

                                                # Automatic detection of the parameter.
                                                if (unc == ObliviousUncertainty || unc == RobustUncertainty) && algo == Automatic
                                                    @test mt.algo != algo()
                                                else
                                                    @test mt.algo == algo()
                                                end
                                            end
                                        end
                                    else
                                        # Generic case.
                                        mt = ModelType(edge_obj(), agg_obj(), type(), cg, algo(), unc(), uncparams())

                                        # The model type object should match the arguments.
                                        @test mt.edge_obj == edge_obj()
                                        @test mt.agg_obj == agg_obj()
                                        @test mt.type == type()
                                        @test mt.cg == cg
                                        @test mt.unc == unc()
                                        @test mt.uncparams == uncparams()

                                        # Automatic detection of the parameter.
                                        if (unc == ObliviousUncertainty || unc == RobustUncertainty) && algo == Automatic
                                            @test mt.algo != algo()
                                        else
                                            @test mt.algo == algo()
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
using PathDistribution

################################################################################
# Datasets from
# Roberts, B., & Kroese, D. P. (2007). Estimating the Number of st Paths in a Graph. J. Graph Algorithms Appl., 11(1), 195-214.
# http://dx.doi.org/10.7155/jgaa.00142

# Case 1
adj_mtx =[  0 1 1 1 0 1 1 1 ;
            1 0 0 0 1 1 1 0 ;
            1 0 0 1 1 1 1 1 ;
            1 0 1 0 1 1 1 1 ;
            0 1 1 1 0 1 0 0 ;
            1 1 1 1 1 0 1 1 ;
            1 1 1 1 0 1 0 1 ;
            1 0 1 1 0 1 1 0     ]
# true_num = 397

# Full Enumeration
path_enums = path_enumeration(1, size(adj_mtx,1), adj_mtx)
x_data, y_data = actual_cumulative_count(path_enums)

# Monte Carlo Sampling
samples = monte_carlo_path_sampling(1, size(adj_mtx,1), adj_mtx)
x_data_est, y_data_est = estimate_cumulative_count(samples)

println("===== Case 1 of Roberts & Kroese (2007) =====")
println("The total number of paths:")
println("- Full enumeration      : $(length(path_enums))")
println("- Monte Carlo estimation: $(y_data_est[end])")

# Other ways of using Monte Carlo Sampling
N1 = 500
N2 = 1000

samples = monte_carlo_path_sampling(1, size(adj_mtx,1), adj_mtx)
x_data_est, y_data_est = estimate_cumulative_count(samples)
println("- Monte Carlo estimation: $(y_data_est[end])")

samples = monte_carlo_path_sampling(1, size(adj_mtx,1), adj_mtx, N1, N2)
x_data_est, y_data_est = estimate_cumulative_count(samples)
println("- Monte Carlo estimation: $(y_data_est[end])")

samples = monte_carlo_path_sampling(1, size(adj_mtx,1), adj_mtx, N1, N2)
x_data_est, y_data_est = estimate_cumulative_count(samples, :uniform)
println("- Monte Carlo estimation: $(y_data_est[end])")

samples = monte_carlo_path_sampling(1, size(adj_mtx,1), adj_mtx, N1, N2)
x_data_est, y_data_est = estimate_cumulative_count(samples, :unique)
println("- Monte Carlo estimation: $(y_data_est[end])")

# When you just need the path number estimator
no_path_est = monte_carlo_path_number(1, size(adj_mtx,1), adj_mtx, N1, N2)
println("- Monte Carlo estimation: $(no_path_est)")

################################################################################
# Another example in README.md

data = [
 1   4  79.0 ;
 1   2  59.0 ;
 2   4  31.0 ;
 2   3  90.0 ;
 2   5   9.0 ;
 2   6  32.0 ;
 3   9  89.0 ;
 3   8  66.0 ;
 3   6  68.0 ;
 3   7  47.0 ;
 4   3  14.0 ;
 4   9  95.0 ;
 4   8  88.0 ;
 5   3  44.0 ;
 5   6  83.0 ;
 6   7  33.0 ;
 6   8  37.0 ;
 7  11  79.0 ;
 7  12  10.0 ;
 8   7  95.0 ;
 8  10   0.0 ;
 8  12  30.0 ;
 9  10   5.0 ;
 9  11  44.0 ;
10  13  79.0 ;
10  14  91.0 ;
11  14  53.0 ;
11  15  80.0 ;
11  13  56.0 ;
12  15  75.0 ;
12  14   1.0 ;
13  14  48.0 ;
14  15  25.0 ;
]

st = round.(Int, data[:,1]) #first column of data
en = round.(Int, data[:,2]) #second column of data
len = data[:,3] #third

# Double them for two-ways.
start_node = [st; en]
end_node = [en; st]
link_length = [len; len]

origin = 1
destination = 15

# Path enumeration test

# Full Enumeration
path_enums = path_enumeration(origin, destination, start_node, end_node, link_length)
x_data, y_data = actual_cumulative_count(path_enums)

# Monte-Carlo estimation
N1 = 5000
N2 = 10000
samples = monte_carlo_path_sampling(origin, destination, start_node, end_node, link_length)
samples = monte_carlo_path_sampling(origin, destination, start_node, end_node, link_length, N1, N2)
x_data_est, y_data_est = estimate_cumulative_count(samples)

println("===== Another Example =====")
println("The total number of paths:")
println("- Full enumeration      : $(length(path_enums))")
println("- Monte Carlo estimation: $(y_data_est[end])")

# END


text:
module Output

using ..UltraDark
import Dates
import PencilFFTs, MPI
import NPZ
using Statistics

export OutputConfig
export SummaryStatistics, SummaryStatisticsMaxRms

"""
    OutputConfig

struct containing information about what to output.

summary_statistics should be another struct whose constructor generates summary
statistics from `t`, `a`, `Δt` and a `grids` object.
If it is to be used with a `PencilGrids` object, each field must be concrete
and have a binary representation that MPI can handle.
"""
struct OutputConfig
    "where to write output"
    directory::String
    "times at which to output"
    output_times::Array{Float64}

    "whether to output boxes"
    box::Bool
    "whether to output slices"
    slice::Bool
    "whether to output ψ"
    psi::Bool
    "whether to output ρ"
    rho::Bool

    "Type of summary statistics to collect"
    summary_statistics::DataType
end

function OutputConfig(
    directory, output_times;
    box=true, slice=false, psi=true, rho=true,
    summary_statistics=SummaryStatistics,
)

    OutputConfig(directory, output_times, box, slice, psi, rho, summary_statistics)
end

"""
    output_grids(grids, output_config, step)

Write output from `grids` as specified in `output_config`
"""
function output_grids(grids, output_config, step)

    if output_config.box
        if output_config.psi
            NPZ.npzwrite(
                joinpath(output_config.directory, "psi_$step.npy"),
                grids.ψx
            )
        end
        if output_config.rho
            NPZ.npzwrite(
                joinpath(output_config.directory, "rho_$step.npy"),
                grids.ρx
            )
        end
    end

    if output_config.slice
        if output_config.psi
            NPZ.npzwrite(
                joinpath(output_config.directory, "psi_slice_$step.npy"),
                grids.ψx[1, :, :]
            )
        end
        if output_config.rho
            NPZ.npzwrite(
                joinpath(output_config.directory, "rho_slice_$step.npy"),
                grids.ρx[1, :, :]
            )
        end
    end
end

function output_grids(grids::PencilGrids, output_config, step)

    #TODO: don't use gather.  This sends all data to one node.  Should instead use multithreaded HDF5 output
    if output_config.box
        if output_config.psi
            output = PencilFFTs.gather(grids.ψx)
            if MPI.Comm_rank(grids.MPI_COMM) == 0
                NPZ.npzwrite(
                    joinpath(output_config.directory, "psi_$step.npy"),
                    output,
                )
            end
        end

        if output_config.rho
            output = PencilFFTs.gather(grids.ρx)
            if MPI.Comm_rank(grids.MPI_COMM) == 0
                NPZ.npzwrite(
                    joinpath(output_config.directory, "rho_$step.npy"),
                    output,
                )
            end
        end
    end

    if output_config.slice
        if output_config.psi
            output = PencilFFTs.gather(grids.ψx[1, :, :])
            if MPI.Comm_rank(grids.MPI_COMM) == 0
                NPZ.npzwrite(
                    joinpath(output_config.directory, "psi_slice_$step.npy"),
                    output[1, :, :],
                )
            end
        end

        if output_config.rho
            output = PencilFFTs.gather(grids.ρx)
            if MPI.Comm_rank(grids.MPI_COMM) == 0
                NPZ.npzwrite(
                    joinpath(output_config.directory, "rho_slice_$step.npy"),
                    output[1, :, :],
                )
            end
        end
    end
end

"""
    output_summary_header(output_config)

Write a header for a summary file

The header contains labels for each column of the summary CSV file.
This function overwrites the current contents of the file.
"""
function output_summary_header(output_config)
    open(joinpath(output_config.directory, "summary.csv"), "w") do file
        write(file, column_titles(output_config.summary_statistics))
    end
end

function generate_summary_row(summary)::String
    line = "$(getfield(summary, 1))"
    for i in 2:nfields(summary)
        line = line * ",$(getfield(summary, i))"
    end
    line = line * "\n"
end

"""
    output_summary_row(grids, output_config, t, a, Δt)

Write a new row to the summary file
"""
function output_summary_row(grids, output_config, t, a, Δt)
    summary = output_config.summary_statistics(Dates.now(), t, a, Δt, grids)
    line = generate_summary_row(summary)
    open(joinpath(output_config.directory, "summary.csv"), "a") do file
        write(file, line)
    end
end

"""
    output_summary_row(grids::PencilGrids, output_config, t, a, Δt)

Write a new row to the summary file

If the grids is a PencilGrids, this uses `MPI.Reduce` to compute partial
summaries in each task and combine them.
"""
function output_summary_row(grids::PencilGrids, output_config, t, a, Δt)
    root = 0
    summary = MPI.Reduce(
                   output_config.summary_statistics(Dates.now(), t, a, Δt, grids),
                   pool_summarystat,
                   root,
                   grids.MPI_COMM
                  )

    line = generate_summary_row(summary)
    if MPI.Comm_rank(grids.MPI_COMM) == 0
        open(joinpath(output_config.directory, "summary.csv"), "a") do file
            write(file, line)
        end
    end
end

"""
    SummaryStatistics

Summary statistics including the time, scale factor and current time step.
"""
struct SummaryStatistics
    "wall time"
    date::Dates.DateTime
    "time"
    t::Float64
    "scale factor"
    a::Float64
    "time step"
    Δt::Float64
end

function SummaryStatistics(wall_time, sim_time, a, Δt, grids)
    SummaryStatistics(wall_time, sim_time, a, Δt)
end

"""
    SummaryStatisticsMeanMaxRms

Summary statistics including the mean density and RMS density contrast.
"""
struct SummaryStatisticsMeanMaxRms
    "wall time"
    date::Dates.DateTime
    "time"
    t::Float64
    "scale factor"
    a::Float64
    "time step"
    Δt::Float64
    "mean of density"
    ρx_mean::Float64
    "max of density"
    ρx_max::Float64
    "RMS of density contrast"
    δx_rms::Float64
    "number of grid points summarized"
    n::Int64
end

function SummaryStatisticsMeanMaxRms(wall_time, sim_time, a, Δt, grids)
    ρx_mean = mean(grids.ρx)
    ρx_max = maximum(grids.ρx)
    δx_rms = mean(((grids.ρx .- ρx_mean).^2))^0.5
    n = prod(size(grids.ρx))

    SummaryStatisticsMeanMaxRms(wall_time, sim_time, a, Δt, ρx_mean, ρx_max, δx_rms, n)
end

"""
    column_titles

Generate column titles for a CSV file from fields names of a struct
"""
function column_titles(stat_struct)
    reduce(replace, [":"=>"", "("=>"", ")"=>"", " "=>""], init="$(fieldnames(stat_struct))") * "\n"
end

"""
    pool_summarystat(S1::SummaryStatistics, S2::SummaryStatistics)

MPI reduction operator for summary statistics.

Check that t, a, Δt are equal and return them.
"""
function pool_summarystat(S1::SummaryStatistics, S2::SummaryStatistics)::SummaryStatistics

    if (S1.t != S2.t) || (S1.a != S2.a) || (S1.Δt != S2.Δt)
        @error "Summaries incompatible across nodes" S1 S2
    end

    S1
end

"""
    pool_summarystat(S1::SummaryStatisticsMeanMaxRms, S2::SummaryStatisticsMeanMaxRms)

MPI reduction operator for summary statistics.
"""
function pool_summarystat(S1::SummaryStatisticsMeanMaxRms, S2::SummaryStatisticsMeanMaxRms)::SummaryStatisticsMeanMaxRms

    if (S1.t != S2.t) || (S1.a != S2.a) || (S1.Δt != S2.Δt)
        @error "Summaries incompatible across nodes" S1 S2
    end

    n = S1.n + S2.n
    ρx_mean = (S1.ρx_mean * S1.n + S2.ρx_mean * S2.n) / n
    ρx_max = maximum(S1.ρx_max, S2.ρx_max)
    δx_rms = ((S1.n * S1.δx_rms^2 + S2.n * S2.δx_rms^2) / n)^0.5

    SummaryStatisticsMeanMaxRms(S1.date, S1.t, S1.a, S1.Δt, ρx_mean, ρx_max, δx_rms, n)
end

end # module


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

using Test

@testset "Rationals" begin
    @test 1//1 == 1
    @test 2//2 == 1
    @test 1//1 == 1//1
    @test 2//2 == 1//1
    @test 2//4 == 3//6
    @test 1//2 + 1//2 == 1
    @test (-1)//3 == -(1//3)
    @test 1//2 + 3//4 == 5//4
    @test 1//3 * 3//4 == 1//4
    @test 1//2 / 3//4 == 2//3
    @test 1//0 == 1//0
    @test 5//0 == 1//0
    @test -1//0 == -1//0
    @test -7//0 == -1//0
    @test  (-1//2) // (-2//5) == 5//4

    @test_throws OverflowError -(0x01//0x0f)
    @test_throws OverflowError -(typemin(Int)//1)
    @test_throws OverflowError (typemax(Int)//3) + 1
    @test_throws OverflowError (typemax(Int)//3) * 2
    @test (typemax(Int)//1) * (1//typemax(Int)) == 1
    @test (typemax(Int)//1) / (typemax(Int)//1) == 1
    @test (1//typemax(Int)) / (1//typemax(Int)) == 1
    @test_throws OverflowError (1//2)^63
    @test inv((1+typemin(Int))//typemax(Int)) == -1
    @test_throws ArgumentError inv(typemin(Int)//typemax(Int))
    @test_throws ArgumentError Rational(0x1, typemin(Int32))

    @test @inferred(rationalize(Int, 3.0, 0.0)) === 3//1
    @test @inferred(rationalize(Int, 3.0, 0)) === 3//1
    @test_throws OverflowError rationalize(UInt, -2.0)
    @test_throws ArgumentError rationalize(Int, big(3.0), -1.)
    # issue 26823
    @test_throws InexactError rationalize(Int, NaN)
    # issue 32569
    @test_throws ArgumentError 1 // typemin(Int)
    @test_throws ArgumentError 0 // 0
    @test -2 // typemin(Int) == -1 // (typemin(Int) >> 1)
    @test 2 // typemin(Int) == 1 // (typemin(Int) >> 1)

    @test_throws InexactError Rational(UInt(1), typemin(Int32))
    @test iszero(Rational{Int}(UInt(0), 1))
    @test Rational{BigInt}(UInt(1), Int(-1)) == -1
    @test_broken Rational{Int64}(UInt(1), typemin(Int32)) == Int64(1) // Int64(typemin(Int32))

    for a = -5:5, b = -5:5
        if a == b == 0; continue; end
        if ispow2(b)
            @test a//b == a/b
            @test convert(Rational,a/b) == a//b
        end
        @test rationalize(a/b) == a//b
        @test a//b == a//b
        if b == 0
            @test_throws DivideError round(Integer,a//b) == round(Integer,a/b)
        else
            @test round(Integer,a//b) == round(Integer,a/b)
        end
        for c = -5:5
            @test (a//b == c) == (a/b == c)
            @test (a//b != c) == (a/b != c)
            @test (a//b <= c) == (a/b <= c)
            @test (a//b <  c) == (a/b <  c)
            @test (a//b >= c) == (a/b >= c)
            @test (a//b >  c) == (a/b >  c)
            for d = -5:5
                if c == d == 0; continue; end
                @test (a//b == c//d) == (a/b == c/d)
                @test (a//b != c//d) == (a/b != c/d)
                @test (a//b <= c//d) == (a/b <= c/d)
                @test (a//b <  c//d) == (a/b <  c/d)
                @test (a//b >= c//d) == (a/b >= c/d)
                @test (a//b >  c//d) == (a/b >  c/d)
            end
        end
    end

    @test 0.5 == 1//2
    @test 0.1 != 1//10
    @test 0.1 == 3602879701896397//36028797018963968
    @test Inf == 1//0 == 2//0 == typemax(Int)//0
    @test -Inf == -1//0 == -2//0 == -typemax(Int)//0
    @test floatmin() != 1//(BigInt(2)^1022+1)
    @test floatmin() == 1//(BigInt(2)^1022)
    @test floatmin() != 1//(BigInt(2)^1022-1)
    @test floatmin()/2 != 1//(BigInt(2)^1023+1)
    @test floatmin()/2 == 1//(BigInt(2)^1023)
    @test floatmin()/2 != 1//(BigInt(2)^1023-1)
    @test nextfloat(0.0) != 1//(BigInt(2)^1074+1)
    @test nextfloat(0.0) == 1//(BigInt(2)^1074)
    @test nextfloat(0.0) != 1//(BigInt(2)^1074-1)

    @test 1/3 < 1//3
    @test !(1//3 < 1/3)
    @test -1/3 < 1//3
    @test -1/3 > -1//3
    @test 1/3 > -1//3
    @test 1/5 > 1//5
    @test 1//3 < Inf
    @test 0//1 < Inf
    @test 1//0 == Inf
    @test -1//0 == -Inf
    @test -1//0 != Inf
    @test 1//0 != -Inf
    @test !(1//0 < Inf)
    @test !(1//3 < NaN)
    @test !(1//3 == NaN)
    @test !(1//3 > NaN)

    # PR 29561
    @test abs(one(Rational{UInt})) === one(Rational{UInt})
    @test abs(one(Rational{Int})) === one(Rational{Int})
    @test abs(-one(Rational{Int})) === one(Rational{Int})

    # inf addition
    @test 1//0 + 1//0 == 1//0
    @test -1//0 - 1//0 == -1//0
    @test_throws DivideError 1//0 - 1//0
    @test_throws DivideError -1//0 + 1//0
    @test Int128(1)//0 + 1//0 isa Rational{Int128}
    @test 1//0 + Int128(1)//0 isa Rational{Int128}
end

@testset "Rational methods" begin
    rand_int = rand(Int8)

    for T in [Int8, Int16, Int32, Int128, BigInt]
        @test numerator(convert(T, rand_int)) == rand_int
        @test denominator(convert(T, rand_int)) == 1

        @test typemin(Rational{T}) == -one(T)//zero(T)
        @test typemax(Rational{T}) == one(T)//zero(T)
        @test widen(Rational{T}) == Rational{widen(T)}
    end

    @test iszero(typemin(Rational{UInt}))

    @test Rational(Float32(rand_int)) == Rational(rand_int)

    @test Rational(Rational(rand_int)) == Rational(rand_int)

    @test begin
        var = -Rational(UInt32(0))
        var == UInt32(0)
    end

    @test Rational(rand_int, 3)/Complex(3, 2) == Complex(Rational(rand_int, 13), -Rational(rand_int*2, 39))

    @test Complex(rand_int, 0) == Rational(rand_int)
    @test Rational(rand_int) == Complex(rand_int, 0)

    @test (Complex(rand_int, 4) == Rational(rand_int)) == false
    @test (Rational(rand_int) == Complex(rand_int, 4)) == false

    @test trunc(Rational(BigInt(rand_int), BigInt(3))) == Rational(trunc(BigInt, Rational(BigInt(rand_int),BigInt(3))))
    @test  ceil(Rational(BigInt(rand_int), BigInt(3))) == Rational( ceil(BigInt, Rational(BigInt(rand_int),BigInt(3))))
    @test round(Rational(BigInt(rand_int), BigInt(3))) == Rational(round(BigInt, Rational(BigInt(rand_int),BigInt(3))))

    for a = -3:3
        @test Rational(Float32(a)) == Rational(a)
        @test Rational(a)//2 == a//2
        @test a//Rational(2) == Rational(a/2)
        @test a.//[-2, -1, 1, 2] == [-a//2, -a//1, a//1, a//2]
        for b=-3:3, c=1:3
            @test b//(a+c*im) == b*a//(a^2+c^2)-(b*c//(a^2+c^2))*im
            for d=-3:3
                @test (a+b*im)//(c+d*im) == (a*c+b*d+(b*c-a*d)*im)//(c^2+d^2)
                @test Complex(Rational(a)+b*im)//Complex(Rational(c)+d*im) == Complex(a+b*im)//Complex(c+d*im)
            end
        end
    end
end

# check type of constructed rationals
int_types = Base.BitInteger64_types
for N = int_types, D = int_types
    T = promote_type(N,D)
    @test typeof(convert(N,2)//convert(D,3)) <: Rational{T}
end

# issue #7564
@test typeof(convert(Rational{Integer},1)) === Rational{Integer}

@testset "issue #15205" begin
    T = Rational
    x = Complex{T}(1//3 + 1//4*im)
    y = Complex{T}(1//2 + 1//5*im)
    xf = Complex{BigFloat}(1//3 + 1//4*im)
    yf = Complex{BigFloat}(1//2 + 1//5*im)
    yi = 4

    @test x^y ≈ xf^yf
    @test x^yi ≈ xf^yi
    @test x^true ≈ xf^true
    @test x^false == xf^false
    @test x^1 ≈ xf^1
    @test xf^Rational(2, 1) ≈ xf*xf
    @test Complex(1., 1.)^Rational(2,1) == Complex(1., 1.)*Complex(1.,1.) == Complex(0., 2.)

    for Tf = (Float16, Float32, Float64), Ti = (Int16, Int32, Int64)
        almost_half  = Rational(div(typemax(Ti),Ti(2))  , typemax(Ti))
        over_half    = Rational(div(typemax(Ti),Ti(2))+one(Ti), typemax(Ti))
        exactly_half = Rational(one(Ti)  , Ti(2))

        @test round( almost_half) == 0//1
        @test round(-almost_half) == 0//1
        @test round(Tf,  almost_half, RoundNearestTiesUp) == 0.0
        @test round(Tf, -almost_half, RoundNearestTiesUp) == 0.0
        @test round(Tf,  almost_half, RoundNearestTiesAway) == 0.0
        @test round(Tf, -almost_half, RoundNearestTiesAway) == 0.0

        @test round( exactly_half) == 0//1 # rounds to closest _even_ integer
        @test round(-exactly_half) == 0//1 # rounds to closest _even_ integer
        @test round(Tf,  exactly_half, RoundNearestTiesUp) == 1.0
        @test round(Tf, -exactly_half, RoundNearestTiesUp) == 0.0
        @test round(Tf,  exactly_half, RoundNearestTiesAway) == 1.0
        @test round(Tf, -exactly_half, RoundNearestTiesAway) == -1.0

        @test round(over_half) == 1//1
        @test round(-over_half) == -1//1
        @test round(Tf,  over_half, RoundNearestTiesUp) == 1.0
        @test round(Tf,  over_half, RoundNearestTiesAway) == 1.0
        @test round(Tf, -over_half, RoundNearestTiesUp) == -1.0
        @test round(Tf, -over_half, RoundNearestTiesAway) == -1.0

        @test round(Tf, 11//2, RoundNearestTiesUp) == 6.0
        @test round(Tf, -11//2, RoundNearestTiesUp) == -5.0
        @test round(Tf, 11//2, RoundNearestTiesAway) == 6.0
        @test round(Tf, -11//2, RoundNearestTiesAway) == -6.0

        @test round(Tf, Ti(-1)//zero(Ti)) == -Inf
        @test round(Tf, one(1)//zero(Ti)) == Inf
        @test round(Tf, Ti(-1)//zero(Ti), RoundNearestTiesUp) == -Inf
        @test round(Tf, one(1)//zero(Ti), RoundNearestTiesUp) == Inf
        @test round(Tf, Ti(-1)//zero(Ti), RoundNearestTiesAway) == -Inf
        @test round(Tf, one(1)//zero(Ti), RoundNearestTiesAway) == Inf

        @test round(Tf, zero(Ti)//one(Ti)) == 0
        @test round(Tf, zero(Ti)//one(Ti), RoundNearestTiesUp) == 0
        @test round(Tf, zero(Ti)//one(Ti), RoundNearestTiesAway) == 0
    end
end
@testset "show and Rationals" begin
    io = IOBuffer()
    rational1 = Rational(1465, 8593)
    rational2 = Rational(-4500, 9000)
    @test sprint(show, rational1) == "1465//8593"
    @test sprint(show, rational2) == "-1//2"
    let
        io1 = IOBuffer()
        write(io1, rational1)
        io1.ptr = 1
        @test read(io1, typeof(rational1)) == rational1

        io2 = IOBuffer()
        write(io2, rational2)
        io2.ptr = 1
        @test read(io2, typeof(rational2)) == rational2
    end
end
@testset "parse" begin
    # Non-negative Int in which parsing is expected to work
    @test parse(Rational{Int}, string(10)) == 10 // 1
    @test parse(Rational{Int}, "100/10" ) == 10 // 1
    @test parse(Rational{Int}, "100 / 10") == 10 // 1
    @test parse(Rational{Int}, "0 / 10") == 0 // 1
    @test parse(Rational{Int}, "100//10" ) == 10 // 1
    @test parse(Rational{Int}, "100 // 10") == 10 // 1
    @test parse(Rational{Int}, "0 // 10") == 0 // 1

    # Variations of the separator that should throw errors
    @test_throws ArgumentError parse(Rational{Int}, "100\\10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 \\ 10")
    @test_throws ArgumentError parse(Rational{Int}, "100\\\\10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 \\\\ 10")
    @test_throws ArgumentError parse(Rational{Int}, "100/ /10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 / / 10")
    @test_throws ArgumentError parse(Rational{Int}, "100// /10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 // / 10")
    @test_throws ArgumentError parse(Rational{Int}, "100///10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 /// 10")
    @test_throws ArgumentError parse(Rational{Int}, "100÷10" )
    @test_throws ArgumentError parse(Rational{Int}, "100 ÷ 10")
    @test_throws ArgumentError parse(Rational{Int}, "100 10" )
    @test_throws ArgumentError parse(Rational{Int}, "100   10")

    # Zero denominator, negative denominator, and double negative
    @test_throws ArgumentError parse(Rational{Int}, "0//0")
    @test parse(Rational{Int}, "1000//-100") == -10 // 1
    @test parse(Rational{Int}, "-1000//-100") == 10 // 1

    # Negative Int tests in which parsing is expected to work
    @test parse(Rational{Int}, string(-10)) == -10 // 1
    @test parse(Rational{Int}, "-100/10" ) == -10 // 1
    @test parse(Rational{Int}, "-100 / 10") == -10 // 1
    @test parse(Rational{Int}, "-100//10" ) == -10 // 1

    # Variations of the separator that should throw errors (negative version)
    @test_throws ArgumentError parse(Rational{Int}, "-100\\10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 \\ 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100\\\\10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 \\\\ 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100/ /10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 / / 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100// /10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 // / 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100///10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 /// 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100÷10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100 ÷ 10")
    @test_throws ArgumentError parse(Rational{Int}, "-100 10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100   10")
    @test_throws ArgumentError parse(Rational{Int}, "-100 -10" )
    @test_throws ArgumentError parse(Rational{Int}, "-100   -10")
    @test_throws ArgumentError parse(Rational{Int}, "100 -10" )
    @test_throws ArgumentError parse(Rational{Int}, "100   -10")
    try  # issue 44570
       parse(Rational{BigInt}, "100 10")
       @test_broken false
    catch
       @test_broken true
    end

    # A few tests for other Integer types
    @test parse(Rational{Bool}, "true") == true // true
    @test parse(Rational{UInt8}, "0xff/0xf") == UInt8(17) // UInt8(1)
    @test parse(Rational{Int8}, "-0x7e/0xf") == Int8(-126) // Int8(15)
    @test parse(Rational{BigInt}, "$(big(typemax(Int))*16)/8") == (big(typemax(Int))*2) // big(1)
    # Mixed notations
    @test parse(Rational{UInt8}, "0x64//28") == UInt8(25) // UInt8(7)
    @test parse(Rational{UInt8}, "100//0x1c") == UInt8(25) // UInt8(7)

    # Out of the bounds tests
    # 0x100 is 256, Int test works for both Int32 and Int64
    # The error must be throw even if the canonicalized fraction fits
    # (i.e., would be less than typemax after divided by 2 in examples below,
    # both over typemax values are even).
    @test_throws OverflowError parse(Rational{UInt8}, "0x100/0x1")
    @test_throws OverflowError parse(Rational{UInt8}, "0x100/0x2")
    @test_throws OverflowError parse(Rational{Int}, "$(big(typemax(Int)) + 1)/1")
    @test_throws OverflowError parse(Rational{Int}, "$(big(typemax(Int)) + 1)/2")
end # parse

@testset "round" begin
    @test round(11//2) == round(11//2, RoundNearest) == 6//1 # rounds to closest _even_ integer
    @test round(-11//2) == round(-11//2, RoundNearest) == -6//1 # rounds to closest _even_ integer
    @test round(13//2) == round(13//2, RoundNearest) == 6//1 # rounds to closest _even_ integer
    @test round(-13//2) == round(-13//2, RoundNearest) == -6//1 # rounds to closest _even_ integer
    @test round(11//3) == round(11//3, RoundNearest) == 4//1 # rounds to closest _even_ integer
    @test round(-11//3) == round(-11//3, RoundNearest) == -4//1 # rounds to closest _even_ integer

    @test round(11//2, RoundNearestTiesAway) == 6//1
    @test round(-11//2, RoundNearestTiesAway) == -6//1
    @test round(13//2, RoundNearestTiesAway) == 7//1
    @test round(-13//2, RoundNearestTiesAway) == -7//1
    @test round(11//3, RoundNearestTiesAway) == 4//1
    @test round(-11//3, RoundNearestTiesAway) == -4//1

    @test round(11//2, RoundNearestTiesUp) == 6//1
    @test round(-11//2, RoundNearestTiesUp) == -5//1
    @test round(13//2, RoundNearestTiesUp) == 7//1
    @test round(-13//2, RoundNearestTiesUp) == -6//1
    @test round(11//3, RoundNearestTiesUp) == 4//1
    @test round(-11//3, RoundNearestTiesUp) == -4//1

    @test trunc(11//2) == round(11//2, RoundToZero) == 5//1
    @test trunc(-11//2) == round(-11//2, RoundToZero) == -5//1
    @test trunc(13//2) == round(13//2, RoundToZero) == 6//1
    @test trunc(-13//2) == round(-13//2, RoundToZero) == -6//1
    @test trunc(11//3) == round(11//3, RoundToZero) == 3//1
    @test trunc(-11//3) == round(-11//3, RoundToZero) == -3//1

    @test ceil(11//2) == round(11//2, RoundUp) == 6//1
    @test ceil(-11//2) == round(-11//2, RoundUp) == -5//1
    @test ceil(13//2) == round(13//2, RoundUp) == 7//1
    @test ceil(-13//2) == round(-13//2, RoundUp) == -6//1
    @test ceil(11//3) == round(11//3, RoundUp) == 4//1
    @test ceil(-11//3) == round(-11//3, RoundUp) == -3//1

    @test floor(11//2) == round(11//2, RoundDown) == 5//1
    @test floor(-11//2) == round(-11//2, RoundDown) == -6//1
    @test floor(13//2) == round(13//2, RoundDown) == 6//1
    @test floor(-13//2) == round(-13//2, RoundDown) == -7//1
    @test floor(11//3) == round(11//3, RoundDown) == 3//1
    @test floor(-11//3) == round(-11//3, RoundDown) == -4//1

    for T in (Float16, Float32, Float64)
        @test round(T, true//false) === convert(T, Inf)
        @test round(T, true//true) === one(T)
        @test round(T, false//true) === zero(T)
        @test trunc(T, true//false) === convert(T, Inf)
        @test trunc(T, true//true) === one(T)
        @test trunc(T, false//true) === zero(T)
        @test floor(T, true//false) === convert(T, Inf)
        @test floor(T, true//true) === one(T)
        @test floor(T, false//true) === zero(T)
        @test ceil(T, true//false) === convert(T, Inf)
        @test ceil(T, true//true) === one(T)
        @test ceil(T, false//true) === zero(T)
    end

    for T in (Int8, Int16, Int32, Int64, Bool)
        @test_throws DivideError round(T, true//false)
        @test round(T, true//true) === one(T)
        @test round(T, false//true) === zero(T)
        @test_throws DivideError trunc(T, true//false)
        @test trunc(T, true//true) === one(T)
        @test trunc(T, false//true) === zero(T)
        @test_throws DivideError floor(T, true//false)
        @test floor(T, true//true) === one(T)
        @test floor(T, false//true) === zero(T)
        @test_throws DivideError ceil(T, true//false)
        @test ceil(T, true//true) === one(T)
        @test ceil(T, false//true) === zero(T)
    end

    # issue 34657
    @test round(1//0) === round(Rational, 1//0) === 1//0
    @test trunc(1//0) === trunc(Rational, 1//0) === 1//0
    @test floor(1//0) === floor(Rational, 1//0) === 1//0
    @test ceil(1//0) === ceil(Rational, 1//0) === 1//0
    @test round(-1//0) === round(Rational, -1//0) === -1//0
    @test trunc(-1//0) === trunc(Rational, -1//0) === -1//0
    @test floor(-1//0) === floor(Rational, -1//0) === -1//0
    @test ceil(-1//0) === ceil(Rational, -1//0) === -1//0
    for r = [RoundNearest, RoundNearestTiesAway, RoundNearestTiesUp,
             RoundToZero, RoundUp, RoundDown]
        @test round(1//0, r) === 1//0
        @test round(-1//0, r) === -1//0
    end

    @test @inferred(round(1//0, digits=1)) === Inf
    @test @inferred(trunc(1//0, digits=2)) === Inf
    @test @inferred(floor(-1//0, sigdigits=1)) === -Inf
    @test @inferred(ceil(-1//0, sigdigits=2)) === -Inf
end

@testset "issue 1552" begin
    @test isa(rationalize(Int8, float(pi)), Rational{Int8})
    @test rationalize(Int8, float(pi)) == 22//7
    @test rationalize(Int64, 0.957762604052997) == 42499549//44373782
    @test rationalize(Int16, 0.929261477046077) == 11639//12525
    @test rationalize(Int16, 0.2264705884044309) == 77//340
    @test rationalize(Int16, 0.39999899264235683) == 2//5
    @test rationalize(Int16, 1.1264233500618559e-5) == 0//1
    @test rationalize(UInt16, 0.6666652791223875) == 2//3
    @test rationalize(Int8, 0.9374813124660655) == 15//16
    @test rationalize(Int8, 0.003803032342443835) == 0//1
end
# issue 3412
@test convert(Rational{Int32},0.5) === Int32(1)//Int32(2)

@testset "issue 6712" begin
    @test convert(Rational{BigInt},Float64(pi)) == Float64(pi)
    @test convert(Rational{BigInt},big(pi)) == big(pi)

    @test convert(Rational,0.0) == 0
    @test convert(Rational,-0.0) == 0
    @test convert(Rational,zero(BigFloat)) == 0
    @test convert(Rational,-zero(BigFloat)) == 0
    @test convert(Rational{BigInt},0.0) == 0
    @test convert(Rational{BigInt},-0.0) == 0
    @test convert(Rational{BigInt},zero(BigFloat)) == 0
    @test convert(Rational{BigInt},-zero(BigFloat)) == 0
    @test convert(Rational{BigInt},5e-324) == 5e-324
    @test convert(Rational{BigInt},floatmin(Float64)) == floatmin(Float64)
    @test convert(Rational{BigInt},floatmax(Float64)) == floatmax(Float64)

    @test isa(convert(Float64, big(1)//2), Float64)
end
@testset "issue 16513" begin
    @test convert(Rational{Int32}, pi) == 1068966896 // 340262731
    @test convert(Rational{Int64}, pi) == 2646693125139304345 // 842468587426513207
    @test convert(Rational{Int128}, pi) == 60728338969805745700507212595448411044 // 19330430665609526556707216376512714945
    @test_throws ArgumentError convert(Rational{BigInt}, pi)
end
@testset "issue 5935" begin
    @test rationalize(Int8,  nextfloat(0.1)) == 1//10
    @test rationalize(Int64, nextfloat(0.1)) == 300239975158034//3002399751580339
    @test rationalize(Int128,nextfloat(0.1)) == 300239975158034//3002399751580339
    @test rationalize(BigInt,nextfloat(0.1)) == 300239975158034//3002399751580339
    @test rationalize(Int8,  nextfloat(0.1),tol=0.5eps(0.1)) == 1//10
    @test rationalize(Int64, nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629
    @test rationalize(Int128,nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629
    @test rationalize(BigInt,nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629
    @test rationalize(Int8,  nextfloat(0.1),tol=1.5eps(0.1)) == 1//10
    @test rationalize(Int64, nextfloat(0.1),tol=1.5eps(0.1)) == 1//10
    @test rationalize(Int128,nextfloat(0.1),tol=1.5eps(0.1)) == 1//10
    @test rationalize(BigInt,nextfloat(0.1),tol=1.5eps(0.1)) == 1//10
    @test rationalize(BigInt,nextfloat(parse(BigFloat,"0.1")),tol=1.5eps(big(0.1))) == 1//10
    @test rationalize(Int64, nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936
    @test rationalize(Int128,nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936
    @test rationalize(BigInt,nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936

    @test rationalize(Int8,  prevfloat(0.1)) == 1//10
    @test rationalize(Int64, prevfloat(0.1)) == 1//10
    @test rationalize(Int128,prevfloat(0.1)) == 1//10
    @test rationalize(BigInt,prevfloat(0.1)) == 1//10
    @test rationalize(BigInt,prevfloat(parse(BigFloat,"0.1"))) == 1//10
    @test rationalize(Int64, prevfloat(0.1),tol=0) == 7205759403792793//72057594037927936
    @test rationalize(Int128,prevfloat(0.1),tol=0) == 7205759403792793//72057594037927936
    @test rationalize(BigInt,prevfloat(0.1),tol=0) == 7205759403792793//72057594037927936

    @test rationalize(BigInt,nextfloat(parse(BigFloat,"0.1")),tol=0) == 46316835694926478169428394003475163141307993866256225615783033603165251855975//463168356949264781694283940034751631413079938662562256157830336031652518559744

    @test rationalize(Int8, 200f0) == 1//0
    @test rationalize(Int8, -200f0) == -1//0

    @test [rationalize(1pi,tol=0.1^n) for n=1:10] == [
                 16//5
                 22//7
                201//64
                333//106
                355//113
                355//113
              75948//24175
             100798//32085
             103993//33102
             312689//99532 ]
end

@testset "issue #12536" begin
    @test Rational{Int16}(1,2) === Rational(Int16(1),Int16(2))
    @test Rational{Int16}(500000,1000000) === Rational(Int16(1),Int16(2))
end
# issue 16311
rationalize(nextfloat(0.0)) == 0//1

@testset "rational-exponent promotion rules (issue #3155)" begin
    @test 2.0f0^(1//3) == 2.0f0^(1.0f0/3)
    @test 2^(1//3) == 2^(1/3)
end

@testset "overflow in rational comparison" begin
    @test 3//2 < typemax(Int)
    @test 3//2 <= typemax(Int)
end

# issue #15920
@test Rational(0, 1) / Complex(3, 2) == 0

# issue #16282
@test_throws MethodError 3 // 4.5im

# issue #31396
@test round(1//2, RoundNearestTiesUp) === 1//1

@testset "Unary plus on Rational (issue #30749)" begin
   @test +Rational(true) == 1//1
   @test +Rational(false) == 0//1
   @test -Rational(true) == -1//1
   @test -Rational(false) == 0//1
end

# issue #27039
@testset "gcd, lcm, gcdx for Rational" begin
    # TODO: Test gcd, lcm, gcdx for Rational{BigInt}.
    for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128)
        a = T(6) // T(35)
        b = T(10) // T(21)
        @test gcd(a, b) === T(2)//T(105)
        @test gcd(b, a) === T(2)//T(105)
        @test lcm(a, b) === T(30)//T(7)
        if T <: Signed
            @test gcd(-a) === a
            @test lcm(-b) === b
            @test gcdx(a, b) === (T(2)//T(105), T(-11), T(4))
            @test gcd(-a, b) === T(2)//T(105)
            @test gcd(a, -b) === T(2)//T(105)
            @test gcd(-a, -b) === T(2)//T(105)
            @test lcm(-a, b) === T(30)//T(7)
            @test lcm(a, -b) === T(30)//T(7)
            @test lcm(-a, -b) === T(30)//T(7)
            @test gcdx(-a, b) === (T(2)//T(105), T(11), T(4))
            @test gcdx(a, -b) === (T(2)//T(105), T(-11), T(-4))
            @test gcdx(-a, -b) === (T(2)//T(105), T(11), T(-4))
        end

        @test gcd(a, T(0)//T(1)) === a
        @test lcm(a, T(0)//T(1)) === T(0)//T(1)
        @test gcdx(a, T(0)//T(1)) === (a, T(1), T(0))

        @test gcdx(T(1)//T(0), T(1)//T(2)) === (T(1)//T(0), T(1), T(0))
        @test gcdx(T(1)//T(2), T(1)//T(0)) === (T(1)//T(0), T(0), T(1))
        @test gcdx(T(1)//T(0), T(1)//T(1)) === (T(1)//T(0), T(1), T(0))
        @test gcdx(T(1)//T(1), T(1)//T(0)) === (T(1)//T(0), T(0), T(1))
        @test gcdx(T(1)//T(0), T(1)//T(0)) === (T(1)//T(0), T(1), T(1))
        @test gcdx(T(1)//T(0), T(0)//T(1)) === (T(1)//T(0), T(1), T(0))
        @test gcdx(T(0)//T(1), T(0)//T(1)) === (T(0)//T(1), T(1), T(0))

        if T <: Signed
            @test gcdx(T(-1)//T(0), T(1)//T(2)) === (T(1)//T(0), T(1), T(0))
            @test gcdx(T(1)//T(2), T(-1)//T(0)) === (T(1)//T(0), T(0), T(1))
            @test gcdx(T(-1)//T(0), T(1)//T(1)) === (T(1)//T(0), T(1), T(0))
            @test gcdx(T(1)//T(1), T(-1)//T(0)) === (T(1)//T(0), T(0), T(1))
            @test gcdx(T(-1)//T(0), T(1)//T(0)) === (T(1)//T(0), T(1), T(1))
            @test gcdx(T(1)//T(0), T(-1)//T(0)) === (T(1)//T(0), T(1), T(1))
            @test gcdx(T(-1)//T(0), T(-1)//T(0)) === (T(1)//T(0), T(1), T(1))
            @test gcdx(T(-1)//T(0), T(0)//T(1)) === (T(1)//T(0), T(1), T(0))
            @test gcdx(T(0)//T(1), T(-1)//T(0)) === (T(1)//T(0), T(0), T(1))
        end

        @test gcdx(T(1)//T(3), T(2)) === (T(1)//T(3), T(1), T(0))
        @test lcm(T(1)//T(3), T(1)) === T(1)//T(1)
        @test lcm(T(3)//T(1), T(1)//T(0)) === T(3)//T(1)
        @test lcm(T(0)//T(1), T(1)//T(0)) === T(0)//T(1)

        @test lcm(T(1)//T(0), T(1)//T(2)) === T(1)//T(2)
        @test lcm(T(1)//T(2), T(1)//T(0)) === T(1)//T(2)
        @test lcm(T(1)//T(0), T(1)//T(1)) === T(1)//T(1)
        @test lcm(T(1)//T(1), T(1)//T(0)) === T(1)//T(1)
        @test lcm(T(1)//T(0), T(1)//T(0)) === T(1)//T(0)
        @test lcm(T(1)//T(0), T(0)//T(1)) === T(0)//T(1)
        @test lcm(T(0)//T(1), T(0)//T(1)) === T(0)//T(1)

        if T <: Signed
            @test lcm(T(-1)//T(0), T(1)//T(2)) === T(1)//T(2)
            @test lcm(T(1)//T(2), T(-1)//T(0)) === T(1)//T(2)
            @test lcm(T(-1)//T(0), T(1)//T(1)) === T(1)//T(1)
            @test lcm(T(1)//T(1), T(-1)//T(0)) === T(1)//T(1)
            @test lcm(T(-1)//T(0), T(1)//T(0)) === T(1)//T(0)
            @test lcm(T(1)//T(0), T(-1)//T(0)) === T(1)//T(0)
            @test lcm(T(-1)//T(0), T(-1)//T(0)) === T(1)//T(0)
            @test lcm(T(-1)//T(0), T(0)//T(1)) === T(0)//T(1)
            @test lcm(T(0)//T(1), T(-1)//T(0)) === T(0)//T(1)
        end

        @test gcd([T(5), T(2), T(1)//T(2)]) === T(1)//T(2)
        @test gcd(T(5), T(2), T(1)//T(2)) === T(1)//T(2)

        @test lcm([T(5), T(2), T(1)//T(2)]) === T(10)//T(1)
        @test lcm(T(5), T(2), T(1)//T(2)) === T(10)//T(1)
    end
end

@testset "Binary operations with Integer" begin
    @test 1//2 - 1 == -1//2
    @test -1//2 + 1 == 1//2
    @test 1 - 1//2 == 1//2
    @test 1 + 1//2 == 3//2
    for q in (19//3, -4//5), i in (6, -7)
        @test rem(q, i) == q - i*div(q, i)
        @test mod(q, i) == q - i*fld(q, i)
    end
    @test 1//2 * 3 == 3//2
    @test -3 * (1//2) == -3//2
    @test (6//5) // -3 == -2//5
    @test -4 // (-6//5) == 10//3

    @test_throws OverflowError UInt(1)//2 - 1
    @test_throws OverflowError 1 - UInt(5)//2
    @test_throws OverflowError 1//typemax(Int64) + 1
    @test_throws OverflowError Int8(1) + Int8(5)//(Int8(127)-Int8(1))
    @test_throws InexactError UInt(1)//2 * -1
    @test_throws OverflowError typemax(Int64)//1 * 2
    @test_throws OverflowError -1//1 * typemin(Int64)

    @test Int8(1) + Int8(4)//(Int8(127)-Int8(1)) == Int8(65) // Int8(63)
    @test -Int32(1) // typemax(Int32) - Int32(1) == typemin(Int32) // typemax(Int32)
    @test 1 // (typemax(Int128) + BigInt(1)) - 2 == (1 + BigInt(2)*typemin(Int128)) // (BigInt(1) + typemax(Int128))
end

@testset "Promotions on binary operations with Rationals (#36277)" begin
    inttypes = (Base.BitInteger_types..., BigInt)
    for T in inttypes, S in inttypes
        U = Rational{promote_type(T, S)}
        @test typeof(one(Rational{T}) + one(S)) == typeof(one(S) + one(Rational{T})) == typeof(one(Rational{T}) + one(Rational{S})) == U
        @test typeof(one(Rational{T}) - one(S)) == typeof(one(S) - one(Rational{T})) == typeof(one(Rational{T}) - one(Rational{S})) == U
        @test typeof(one(Rational{T}) * one(S)) == typeof(one(S) * one(Rational{T})) == typeof(one(Rational{T}) * one(Rational{S})) == U
        @test typeof(one(Rational{T}) // one(S)) == typeof(one(S) // one(Rational{T})) == typeof(one(Rational{T}) // one(Rational{S})) == U
    end
    @test (-40//3) // 0x5 == 0x5 // (-15//8) == -8//3
    @test (-4//7) // (0x1//0x3) == (0x4//0x7) // (-1//3) == -12//7
    @test -3//2 + 0x1//0x1 == -3//2 + 0x1 == 0x1//0x1 + (-3//2) == 0x1 + (-3//2) == -1//2
    @test 0x3//0x5 - 2//3 == 3//5 - 0x2//0x3 == -1//15
    @test rem(-12//5, 0x2//0x1) == rem(-12//5, 0x2) == -2//5
    @test mod(0x3//0x1, -4//7) == mod(0x3, -4//7) == -3//7
    @test -1//5 * 0x3//0x2 == 0x3//0x2 * -1//5 == -3//10
    @test -2//3 * 0x1 == 0x1 * -2//3 == -2//3
end

@testset "ispow2 and iseven/isodd" begin
    @test ispow2(4//1)
    @test ispow2(1//8)
    @test !ispow2(3//8)
    @test !ispow2(0//1)
    @test iseven(4//1) && !isodd(4//1)
    @test !iseven(3//1) && isodd(3//1)
    @test !iseven(3//8) && !isodd(3//8)
end

@testset "checked_den with different integer types" begin
    @test Base.checked_den(Int8(4), Int32(8)) == Base.checked_den(Int32(4), Int32(8))
end

@testset "Rational{T} with non-concrete T (issue #41222)" begin
    @test @inferred(Rational{Integer}(2,3)) isa Rational{Integer}
end

@testset "issue #41489" begin
    @test Core.Compiler.return_type(+, NTuple{2, Rational}) == Rational
    @test Core.Compiler.return_type(-, NTuple{2, Rational}) == Rational

    A=Rational[1 1 1; 2 2 2; 3 3 3]
    @test @inferred(A*A) isa Matrix{Rational}
end

@testset "issue #42560" begin
    @test rationalize(0.5 + 0.5im) == 1//2 + 1//2*im
    @test rationalize(float(pi)im) == 0//1 + 165707065//52746197*im
    @test rationalize(Int8, float(pi)im) == 0//1 + 22//7*im
    @test rationalize(1.192 + 2.233im) == 149//125 + 2233//1000*im
    @test rationalize(Int8, 1.192 + 2.233im) == 118//99 + 67//30*im
end

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia


text:
function compare_assert(a::Int64, b::Int64)
    @assert(a == b)
    @assert(!(0 == 1))
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(true)
    @assert(!(false))
    compare_assert(1, 1)
    @assert(true)
    @assert(true)
    println("OK")
end


text:
include("DM1_1.jl")
include("DM2_1.jl")
include("DM1_2.jl")

#======================================================================#
function ReactiveGRASP(
    matrix, #matrice de taille m*n representant les condion de notre SPP
    cost, #liste des couts de taille n de notre SPP
    n, #nombre de variable de notre SPP
    m, #nombre de condtion de notre SPP
    ite, #nombre de tirage au sort de alpha avant recalcul des probabilités
    coupe, #nombre de coupe du segment [0,1], si on à pas de liste alpha personalisé
    alphaset, #permet de passer une liste de alpha personnelle
    temps #ressource en seconde alloué a notre reactive-GRASP
    )

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe, alphaset)
    evol_p=Float64[]
    t=time()
    z_global = zeros(Int64, length(p))
    ite_global = zeros(Int64, length(p))
    nb_boucle = 0

    while (time()-t <= temps)
        append!(evol_p,p) #MaJ de notre historique de probabilité
        cpt=1
        nb_boucle += 1
        #On s'assure que chaque alpha s'exprime au moins une fois.
        for i in 1:length(p)
            (SOL,z, crts) = GRASP(cost, matrix, n, m, p[i])
            z_cumul[i] += z
            if z < zWorst
                zWorst = z
            elseif z > zBest
                zBest = z
            end
        end #fin for
        while (cpt <= ite)
            #Selection du alpha
            prob = rand(Float64)
            alpha_choisit = choix_alpha(p,prob)
            #On lance GRASP avec cet alpha
            (SOL, z, crts) = GRASP(cost, matrix, n, m, p[alpha_choisit])
            #Amelioration
            (SOL, z) = exchange1_1(SOL,n,m,cost,crts,matrix)
            #Réinitialisation
            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z
            #Mise à jours
            if z > zBest
                zBest = z
            end #fin if
            if z < zWorst
                zWorst = z
            end #fin if
            #Incrément
            cpt = cpt+1
        end #Fin while
        #On recalcul les probabilitées
        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
        #On sovegarde le nombre d'itération et les valeurs de z
        for i in 1:length(p)
            z_global[i] = z_global[i] + z_cumul[i]
            ite_global[i] = ite_global[i] + nb_iteration[i]
        end

        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end #fin while
    #Calcul de aAvg
    moyenne_global=zeros(Float64,length(p))
    for i in 1:length(p)
        moyenne_global[i] = z_global[i]/ite_global[i]
    end
    zAvg = sum(moyenne_global)/length(p)

    println("zBest: ", zBest,"   zAvg:  ", zAvg, "    zWorst: ", zWorst , " nombre de recalcul de p: ", nb_boucle)
#    println(evol_p) #A decommenter si l'on souhaite afficher l'évolution des probabilités
    return(evol_p)
end #fin reactive-GRASP

#=======================================#
function intialiser(matrix, cost, n, m, ite, coupe, alphaset) #initialisation des variables
    p = Float64[]
    if alphaset == 0 #On verifie si on posséde ou non un set
        if coupe >= 1 #On verifie que les coupes sont admissibles
            for i in 1:coupe
                push!(p,i/coupe)
            end #fin for
        end #fin if
    else #on initialise le alphaset
        p = alphaset
    end #fin if

    nb_iteration = ones(Int64, length(p)) #On s'assurera plus tard que tous les alpha s'esprime au moins une fois
    z_cumul = zeros(Int64, length(p)) #Sera une liste de stockage pour le recalcul des probabilités

    (SOL,z) =Glouton(cost, matrix, n, m) #Initialisation d'une solution de base
    zBest = z
    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end #fin initialiser

#=======================================#
function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    #initialisation
    q=zeros(Float64,length(p))
    somme_q=0
#    println("p:  ", p)
#    println("zcummul:  ", z_cumul)
#    println("nb_ite;   ", nb_iteration)
#    println("zBest: ", zBest, "    zWorst: ", zWorst)

    for i in 1:length(p)
        moyenne = (z_cumul[i]/nb_iteration[i])
        q[i] = (moyenne-zWorst)/(zBest-zWorst)
#        println(i, " moyenne: ", moyenne, "   q[i] : " , q[i])
        somme_q += q[i]
    end
#    println("--------------------")
    for i in 1:length(q)
        p[i] = q[i]/somme_q
    end

    for i in 2:length(p)-1 #On additionne les proabilités entre elles
        p[i] += p[i-1]
    end
    p[length(p)] = 1 #On s'assure qu'il n'y ai pas d'erreur d'arrondie pour 1
end #fin recalcul_p!

#===========================================#
function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end #fin if
    end #fin for
end #fin choix_alpha

max_stars_repo_path: Archives/DM1_Metaheuristique-master, clean avec ca DM2_2.jl


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

function setup_optimizer(method::ADMM, y::Array{Float32,3}, psf::Array{Float32,3}, σsq::Float32, ν::Float32, ρ::Float32)
    num_clip = sum(y .< -σsq)
    if num_clip > 0
        ratio = num_clip / length(y) * 100
        println("Clipped voxels: $(ratio)%")
    else
        println("No clipping performed.")
    end

    y[y .< -σsq] .= -σsq
    y_shape = size(y)
    data_shape = size(psf)
    @assert all(y_shape .<= data_shape)

    M, N, L = size(y)
    ypad = zeros(Float32, data_shape)
    ypad[1:M, 1:N, 1:L] .= y

    mask_in = zeros(Bool, data_shape)
    mask_in[1:M,1:N,1:L] .= true
    mask_out = .!mask_in

    κ = ν / ρ

    psf = fftshift(psf)

    ypad = to_gpu_or_not_to_gpu(ypad)
    mask_in = to_gpu_or_not_to_gpu(mask_in)
    mask_out = to_gpu_or_not_to_gpu(mask_out)
    psf = to_gpu_or_not_to_gpu(psf)

    H = rfft(psf)
    dft_mat = rfft3_operator(ypad)
    idft_mat = dft_mat'

    # Image formation operator
    K_1 = LinearOperator(data_shape, data_shape,
        x -> idft_mat * (H .* (dft_mat * x)),
        x -> idft_mat * (conj.(H) .* (dft_mat * x)), false, Float32)

    # Identity operator
    K_2 = IdentityOperator(data_shape)

    K_3, Dxx = SecondOrderDifferentialOperator3D(data_shape, (2,2))

    K_4, Dyy = SecondOrderDifferentialOperator3D(data_shape, (1,1))

    K_5, Dzz = SecondOrderDifferentialOperator3D(data_shape, (3,3))

    # Dxy + Dyx = √2 * Dxy = Dxy_2
    K_6, Dxy_2 = SecondOrderDifferentialOperator3D(data_shape, (1,2))

    K_7, Dxz_2 = SecondOrderDifferentialOperator3D(data_shape, (1,3))

    K_8, Dyz_2 = SecondOrderDifferentialOperator3D(data_shape, (2,3))

    K = StackedLinearOperator([K_1, K_2, K_3, K_4, K_5, K_6, K_7, K_8])

    update_x! = function(x::T, z::Vector{T}, λ::Vector{T}) where T<:AbstractArray{Float32,3}
        Z = [dft_mat * (z_j .- λ_j) for (z_j, λ_j) in zip(z, λ)]
        tmp = (conj.(H) .* Z[1] .+ Z[2] .+
               conj.(Dxx) .* Z[3] .+ conj.(Dyy) .* Z[4].+
               conj.(Dzz) .* Z[5] .+ conj.(Dxy_2) .* Z[6] .+
               conj.(Dxz_2) .* Z[7] .+ conj.(Dyz_2) .* Z[8]) ./
              (abs2.(H) .+ abs2.(Dxx) .+ abs2.(Dyy) .+ abs2.(Dzz) .+
               abs2.(Dxy_2) .+ abs2.(Dxz_2) .+ abs2.(Dyz_2) .+ 1.0f0)
        x .= idft_mat * tmp
    end

    update_z! = function(Kx::Vector{T}, z::Vector{T}, λ::Vector{T}) where T<:AbstractArray{Float32,3}
        v = [Kx[i] .+ λ[i] for i in 1:length(z)]
        tmp = (-1.0f0 .+ ρ .* (v[1] .- σsq)) ./ 2ρ
        z[1] .= mask_out .* v[1] .+ mask_in .* (tmp .+ sqrt.(max.(tmp.^2 .+ (ypad .+ (ρ * σsq) .* v[1]) ./ ρ, 0.0f0)))
        z[2] .= max.(v[2], 0.0f0)
        v_norm = sqrt.(max.(sum(v[i].^2 for i in 3:8), floatmin(Float32)))
        for i in 1:6
            z[i + 2] .= block_soft_threshold.(v[i + 2], v_norm, κ)
        end
    end

    optimizer = initialize_optimizer(method, K, update_x!, update_z!)
end


text:
function tmp_sol2vec(sol,j)
    vect = Vector{Float64}(undef,length(sol.u))
    for i = 1:length(sol.u)
        vect[i] = sol.u[i][j]
    end
    return vect
end

ho_model1 = @min_reaction_network rn begin
    (1.0,1.0), X + 2X ↔ Y + X
end
addodes!(ho_model1)
addjumps!(ho_model1)
ho_model2 = @min_reaction_network rn begin
    (X^3/6,Y * X), X + 2X ⟺ Y + X
end
addodes!(ho_model2)
addjumps!(ho_model2)
ho_model3 = @min_reaction_network rn begin
    (X*(X-1)*(X-2)/6,Y * X), X + 2X ⟺ Y+ X
end
u0 = [1000.,1000.]
tspan = (0.,100.)

prob_det1 = ODEProblem(ho_model1, u0, tspan)
prob_det2 = ODEProblem(ho_model2, u0, tspan)
sol_det1 = solve(prob_det1,Tsit5())
sol_det2 = solve(prob_det2,Tsit5())
# @test tmp_sol2vec(sol_det1,1) == tmp_sol2vec(sol_det2,1)
# @test tmp_sol2vec(sol_det1,2) == tmp_sol2vec(sol_det2,2)
@test isapprox(tmp_sol2vec(sol_det1,1), tmp_sol2vec(sol_det2,1), atol=1e-7)
@test isapprox(tmp_sol2vec(sol_det1,2), tmp_sol2vec(sol_det2,2), atol=1e-7)

disc_prob =  DiscreteProblem([1000,1000],(0.,1.))
jump_prob1 = JumpProblem(disc_prob,Direct(),ho_model1)
jump_prob2 = JumpProblem(disc_prob,Direct(),ho_model2)
sol_jump1 = solve(jump_prob1,FunctionMap(),maxiters = 1e6)
sol_jump2 = solve(jump_prob2,FunctionMap(),maxiters = 1e6)

@test 0.95 < std(tmp_sol2vec(sol_jump1,1)) / std(tmp_sol2vec(sol_jump2,1)) < 1.05
@test 0.95 < std(tmp_sol2vec(sol_jump1,2)) / std(tmp_sol2vec(sol_jump2,2)) < 1.05


text:
using Mea.Green
using Base.Test

 @testset "Testing Green.jl" begin

    # setup
    fin_gf_to = "./data/self_short_moy.dat"
    (zn, gf_c) = Green.readgreen_c(fin_gf_to, zn_col=1)
    gf_t = [(7.29148945e+00 -1.25688297e+00im) (4.31791692e+00 -7.37366863e-01im) (-4.37618486e+00 +7.44382438e-01im) (6.02348575e+00 -9.90968651e-01im) (4.05606418e+00 -4.07146712e-01im);
            (6.37534240e+00 -3.08656815e+00im) (3.92112500e+00 -1.86394356e+00im) (-3.80134014e+00 +1.78496589e+00im) (5.15813767e+00 -2.33078630e+00im) (3.66679027e+00 -9.42835945e-01im)]

    @testset "test_readgreen_c" begin:
        """ """
        zn_test = [5.235989e-002, 1.570800e-001]
        @test isapprox(zn, zn_test, atol=1e-6)

        gf_c_test = Green.t_to_c(zn, gf_t)
        gf_t_test = Green.c_to_t(zn, gf_c_test)

        @test isapprox(gf_c, gf_c_test)
        @test isapprox(gf_t, gf_t_test)

    end

    @testset "test_c" begin

        gf_c_test = zeros(Complex{Float64}, (size(zn)[1], 4, 4))
        gf_c_test[:, 1, 1] = gf_t[:, 1] ; gf_c_test[:, 1, 2] = gf_t[:, 3] ; gf_c_test[:, 1, 3] = gf_t[:, 4] ; gf_c_test[:, 1, 4] = gf_t[:, 3]
        gf_c_test[:, 2, 1] = gf_c_test[:, 1, 2]; gf_c_test[:, 2, 2] = gf_t[:, 2] ; gf_c_test[:, 2, 3] = gf_c_test[:, 1, 2] ; gf_c_test[:, 2, 4] = gf_t[:, 5]
        gf_c_test[:, 3, 1] = gf_c_test[:, 1, 3] ; gf_c_test[:, 3, 2] = gf_c_test[:, 2, 3] ; gf_c_test[:, 3, 3] = gf_c_test[:, 1, 1] ; gf_c_test[:, 3, 4] = gf_c_test[:, 1, 2]
        gf_c_test[:, 4, 1] = gf_c_test[:, 1, 4] ; gf_c_test[:, 4, 2] = gf_c_test[:, 2, 4] ; gf_c_test[:, 4, 3] = gf_c_test[:, 3, 4] ; gf_c_test[:, 4, 4] = gf_c_test[:, 2, 2]

        gf_c_test_02 = [(6.02348575e+00 -9.90968651e-01im),  (5.15813767e+00 -2.33078630e+00im)]#, (4.28316747e+00 -2.88633370e+00im), (6.68961495e-02 +1.31238862e-02im)]
        gf_c_test_00 = [(7.29148945e+00 -1.25688297e+00im), (6.37534240e+00 -3.08656815e+00im)]#,  (5.42233260e+00 -4.06248226e+00im), (5.55747767e+00 -2.25696764e+00im)]

        @test isapprox(gf_c, gf_c_test, atol=1e-6)
        @test isapprox(gf_c[:, 1, 3], gf_c_test_02, atol=1e-3)
        @test isapprox(gf_c[:, 1, 1], gf_c_test_00,  atol=1e-7)

   end

    @testset "test_t_to_c" begin

        gf_to_test = Green.c_to_to(zn, gf_c)
        gf_t_test = Green.c_to_t(zn, gf_c)

        gf_test1_c = Green.to_to_c(zn, gf_to_test)
        gf_test2_c = Green.t_to_c(zn, gf_t_test)

        @test isapprox(gf_c, gf_test1_c, rtol=1e-7, atol=1e-7)
        @test isapprox(gf_c, gf_test2_c, rtol=1e-7, atol=1e-7)

    end

 end  # testset


text:
#Script for calculating blockade shift coefficients
#26/07/2017

using JLD
include("functions.jl")

nnc6 = collect(35:82)
atom = "87Rb"

c6S1_2 = zeros(length(nnc6))
c6D3_2 = zeros(length(nnc6))
c6D5_2 = zeros(length(nnc6))

for kk = 1:length(nnc6)
  junk1, junk2, junk3, c6S1_2[kk] = BlockadeShift(atom,nnc6[kk],0,0.5,0.5)
  junk1, junk2, junk3, c6D3_2[kk] = BlockadeShift(atom,nnc6[kk],2,1.5,1.5)
  junk1, junk2, junk3, c6D5_2[kk] = BlockadeShift(atom,nnc6[kk],2,2.5,2.5)
end

save("Data/C6data.jld","nnc6",nnc6,"c6S1_2",c6S1_2,"c6D3_2",c6D3_2,"c6D5_2",c6D5_2)

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
module RANSACVisualizer

using LinearAlgebra

using AbstractPlotting: Point3f0, legend, Face
using AbstractPlotting
using Makie: scatter, linesegments!, Scene, scatter!, lines, arrows!, wireframe, lines!
using Makie: mesh!, cam3d!, Sphere, Point
using GeometryTypes: Cylinder, Circle, HomogenousMesh
using RANSAC
using RANSAC: midpoint, contournormal, outwardsnormal, segmentnormal
using RANSAC: AbstractTranslationalSurface
using StaticArrays
using CSGBuilding
using ColorSchemes
using Parameters: @unpack

# for cone
using Makie
# using AbstractPlotting
# using FileIO
using GeometryTypes
#using Reactive

#=
# for raytracer
#using ImageView
using Images
using ColorVectorSpace
using Colors
#using Dates
using SignedDistanceFields
=#

export  showgeometry,
        showgeometry!,
        showcandlength,
        showshapes,
        showshapes!,
        getrest,
        showtype,
        showbytype,
        showbytype!,
        plotshape,
        plotshape!,
        plotimplshape,
        plotimplshape!,
        drawcircles!,
        wframe,
        plotspantree!,
        givelargest

export  ConeMesh

#=
export  AbstractCamera,
        IsometricCamera,
        PinholeCamera,
        RayResults,
        raymarch,
        los,
        DefIsoCam,
        DefPinCam

export  render,
        DefaultShaderArray
=#
export  plotcontour,
        compandcenter,
        translheatmap,
        translheatmap!

include("visualizations.jl")
include("cone3.jl")
include("translational.jl")
#include("raytracer/raymarching.jl")
#include("raytracer/camera.jl");
#include("raytracer/render.jl");

end # module


text:
using Revise
using DynamicEnergyBudgets
using CompositeFieldVectors
using AxisArrays
using ForwardDiff
using DiffEqSensitivity
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

du = [0.0 for i in 1:12]
u = [0.0, 1e-2, 0.0, 1e-2, 1e-2, 1e-2, 0.0, 1e-2, 0.0, 1e-2, 1e-2, 10.0]
organism = DynamicEnergyBudgets.UntypedPlant();
p = [flatten(organism)...];
n = [flatten_fieldnames(organism)...];
np = AxisArray(p, Axis{:parameters}(n));
jac = ForwardDiff.jacobian(p1->organism(du, u, p1, 1), np)
@test sum(jac) != 0.0

organism.records[1].J


text:
using SafeTestsets

@safetestset "Atomic Solutions                                                                " begin include("atomic_solutions_tests.jl") end
@safetestset "Stochastic Solutions                                                            " begin include("stochastic_solutions_tests.jl") end
@safetestset "Stochastic Tableaus                                                             " begin include("stochastic_tableaus_tests.jl") end
@safetestset "Stochastic Integrators                                                          " begin include("stochastic_integrators_tests.jl") end


text:
using CsminWel
using Base.Test

f(x) = sum(x.^2)
function grad(x)
    2*x, false
end

function g!(stor, x)
    stor[:] = 2*x
end

res0 = CsminWel.csminwel(f, grad, [0.1,0.1])
res1 = CsminWel.csminwel(f, [0.1,0.1])

res2 = CsminWel.csminwel(f, grad, [0.1,0.1])
res3 = CsminWel.csminwel(f, [0.1,0.1], H = eye(2))

res4 = optimize(f, g!, [.1, .1], Csminwel())
res5 = optimize(f, [.1, .1], Csminwel())   ## Default to finite difference
res6 = optimize(OnceDifferentiable(f, [.1, .1]; autodiff = :forward), [.1, .1], Csminwel())

@test Optim.minimum(res0) ≈ Optim.minimum(res1) atol=1.0e-9
@test Optim.minimum(res1) ≈ Optim.minimum(res2) atol=1.0e-9
@test Optim.minimum(res2) ≈ Optim.minimum(res3) atol=1.0e-9
@test Optim.minimum(res3) ≈ Optim.minimum(res4) atol=1.0e-9
@test Optim.minimum(res4) ≈ Optim.minimum(res5) atol=1.0e-9
@test Optim.minimum(res6) ≈ Optim.minimum(res6) atol=1.0e-9

#Maximizing loglikelihood logistic models
using StatsFuns
srand(1)
x = [ones(200) randn(200,4)]
y = [rand() < 0.5 ? 1. : 0. for j in 1:200]

function loglik(beta)
    xb = x*beta
    sum(-y.*xb + log1pexp.(xb))
end

function dloglik(beta)
    xb = x*beta
    px = logistic.(xb)
    -x'*(y.-px)
end

function g!(stor, beta)
    copy!(stor, dloglik(beta))
end

res1 = optimize(loglik, g!, zeros(5), BFGS())
res2 = optimize(loglik, g!, zeros(5), Csminwel())
res3 = optimize(loglik, zeros(5), Csminwel())
res4 = optimize(Optim.OnceDifferentiable(loglik, zeros(5), autodiff = :forward), zeros(5), Csminwel())
res5 = optimize(Optim.OnceDifferentiable(loglik, zeros(5), autodiff = :finite), zeros(5), Csminwel())

@test Optim.minimum(res1) ≈ Optim.minimum(res2) atol=1.0e-9
@test Optim.minimum(res2) ≈ Optim.minimum(res3) atol=1.0e-9
@test Optim.minimum(res3) ≈ Optim.minimum(res4) atol=1.0e-9
@test Optim.minimum(res4) ≈ Optim.minimum(res5) atol=1.0e-9

res = optimize(loglik, g!, zeros(5), Csminwel(), Optim.Options(extended_trace = true, show_trace = true, store_trace = true))

@test isa(res.trace, Array{Optim.OptimizationState{Float64,CsminWel.Csminwel},1})


text:
include("load.jl")
SyncBarriersTests.runtests()