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

mutable struct VirtualMachineScaleSetUpdateVMProfile <: SwaggerModel
    osProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetUpdateOSProfile } # spec name: osProfile
    storageProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetUpdateStorageProfile } # spec name: storageProfile
    networkProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetUpdateNetworkProfile } # spec name: networkProfile
    securityProfile::Any # spec type: Union{ Nothing, SecurityProfile } # spec name: securityProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, DiagnosticsProfile } # spec name: diagnosticsProfile
    extensionProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetExtensionProfile } # spec name: extensionProfile
    licenseType::Any # spec type: Union{ Nothing, String } # spec name: licenseType
    billingProfile::Any # spec type: Union{ Nothing, BillingProfile } # spec name: billingProfile
    scheduledEventsProfile::Any # spec type: Union{ Nothing, ScheduledEventsProfile } # spec name: scheduledEventsProfile

    function VirtualMachineScaleSetUpdateVMProfile(;osProfile=nothing, storageProfile=nothing, networkProfile=nothing, securityProfile=nothing, diagnosticsProfile=nothing, extensionProfile=nothing, licenseType=nothing, billingProfile=nothing, scheduledEventsProfile=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("osProfile"), osProfile)
        setfield!(o, Symbol("osProfile"), osProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("storageProfile"), storageProfile)
        setfield!(o, Symbol("storageProfile"), storageProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("networkProfile"), networkProfile)
        setfield!(o, Symbol("networkProfile"), networkProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("securityProfile"), securityProfile)
        setfield!(o, Symbol("securityProfile"), securityProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("extensionProfile"), extensionProfile)
        setfield!(o, Symbol("extensionProfile"), extensionProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("billingProfile"), billingProfile)
        setfield!(o, Symbol("billingProfile"), billingProfile)
        validate_property(VirtualMachineScaleSetUpdateVMProfile, Symbol("scheduledEventsProfile"), scheduledEventsProfile)
        setfield!(o, Symbol("scheduledEventsProfile"), scheduledEventsProfile)
        o
    end
end # type VirtualMachineScaleSetUpdateVMProfile

const _property_map_VirtualMachineScaleSetUpdateVMProfile = Dict{Symbol,Symbol}(Symbol("osProfile")=>Symbol("osProfile"), Symbol("storageProfile")=>Symbol("storageProfile"), Symbol("networkProfile")=>Symbol("networkProfile"), Symbol("securityProfile")=>Symbol("securityProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"), Symbol("extensionProfile")=>Symbol("extensionProfile"), Symbol("licenseType")=>Symbol("licenseType"), Symbol("billingProfile")=>Symbol("billingProfile"), Symbol("scheduledEventsProfile")=>Symbol("scheduledEventsProfile"))
const _property_types_VirtualMachineScaleSetUpdateVMProfile = Dict{Symbol,String}(Symbol("osProfile")=>"VirtualMachineScaleSetUpdateOSProfile", Symbol("storageProfile")=>"VirtualMachineScaleSetUpdateStorageProfile", Symbol("networkProfile")=>"VirtualMachineScaleSetUpdateNetworkProfile", Symbol("securityProfile")=>"SecurityProfile", Symbol("diagnosticsProfile")=>"DiagnosticsProfile", Symbol("extensionProfile")=>"VirtualMachineScaleSetExtensionProfile", Symbol("licenseType")=>"String", Symbol("billingProfile")=>"BillingProfile", Symbol("scheduledEventsProfile")=>"ScheduledEventsProfile")
Base.propertynames(::Type{ VirtualMachineScaleSetUpdateVMProfile }) = collect(keys(_property_map_VirtualMachineScaleSetUpdateVMProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetUpdateVMProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetUpdateVMProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetUpdateVMProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetUpdateVMProfile[property_name]

function check_required(o::VirtualMachineScaleSetUpdateVMProfile)
    true
end

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


text:
# ----------------- #
# Utility functions #
# ----------------- #

"Tells maximum number of arguments for a generic or anonymous function"
function max_args(f::Function)
    if isgeneric(f)
        return methods(f).max_args
    else
        # anonymous function
        # NOTE: This might be quite fragile, but works on 0.3.6 and 0.4-dev
        return length(Base.uncompressed_ast(f.code).args[1])
    end
end

"""
Inplace version of `row_kron`. Sets the `i`th row of `out` equal to
`kron(A[i, :], B[i, :])`
"""
function row_kron!(A::AbstractMatrix, B::AbstractMatrix, out::AbstractMatrix)
    # get input dimensions
    nobsa, na = size(A)
    nobsb, nb = size(B)

    nobsa != nobsb && error("A and B must have same number of rows")

    # fill in each element. To do this we make sure we access each array
    # consistent with its column major memory layout.
    @inbounds for ia=1:na, ib=1:nb, t=1:nobsa
        out[t, nb*(ia-1) + ib] = A[t, ia] * B[t, ib]
    end
    out
end

"""
Compute the row-wise Kronecker product between two matrices `A` and `B`.
The matrices are assumed to be of size `(r, ca)` and `(r, cb)`,
respectively. The output matrix will be `(r, ca*cb)`

The `i`th row of the output matrix will be equal to
`kron(A[i,:], B[i,:])`
"""
:row_kron

function row_kron(A::AbstractMatrix{S}, B::AbstractMatrix{T}) where {S,T}
    nobsa, na = size(A)
    nobsb, nb = size(B)
    out = Array(promote_type(S, T), nobsa, na*nb)
    row_kron!(A, B, out)
    out
end

# specialized version for sparse matrices
function row_kron(A::SparseMatrixCSC{S}, B::SparseMatrixCSC{T}) where {S,T}
    nobsa, na = size(A)
    nobsb, nb = size(B)
    out = spzeros(promote_type(S, T), nobsa, na*nb)
    row_kron!(A, B, out)
    out
end

function gettril(x::Array{Float64, 2})
  n, m = size(x)
  a = Array{eltype(x)}(undef, convert(Int, n.*(n+1)/2))
  k::Int = 1
  @inbounds for i = 1:n
    for j = 1:i
      a[k] = x[i, j]
      k += 1
    end
  end
  return a
end

ext: jl


text:
OP, python, libpython = ARGS

# Special exit codes for this script.
# (see: https://www.tldp.org/LDP/abs/html/exitcodes.html)
code_no_precompile_needed = 113

if VERSION < v"0.7.0"
    error("Unsupported Julia version $VERSION")
end

using Pkg
using InteractiveUtils

@info "Julia version info"
versioninfo(verbose=true)

@info "Julia executable: $(Base.julia_cmd().exec[1])"

pkgid = Base.PkgId(Base.UUID(0x438e738f_606a_5dbb_bf0a_cddfbfd45ab0), "PyCall")
pycall_is_installed = Base.locate_package(pkgid) !== nothing

@info "Trying to import PyCall..."

module DummyPyCall
python = nothing
libpython = nothing
end

try
    # `import PyCall` cannot be caught?
    global PyCall = Base.require(Main, :PyCall)
catch err
    @error "`import PyCall` failed" exception=(err, catch_backtrace())
    global PyCall = DummyPyCall
end

ENV["PYTHON"] = python

# TODO: warn if some relevant JULIA_* environment variables are set
# TODO: use PackageSpec to specify PyCall's UUID

print_logfile = true

function build_pycall()
    if VERSION >= v"1.1.0-rc1"
        @info """Run `Pkg.build("PyCall"; verbose=true)`"""
        Pkg.build("PyCall"; verbose=true)
        global print_logfile = false
    else
        @info """Run `Pkg.build("PyCall")`"""
        Pkg.build("PyCall")
    end
end

if OP == "build"
    build_pycall()
elseif PyCall.python == python || PyCall.libpython == libpython
    @info """
    PyCall is already installed and compatible with Python executable.

    PyCall:
        python: $(PyCall.python)
        libpython: $(PyCall.libpython)
    Python:
        python: $python
        libpython: $libpython
    """
    exit(code_no_precompile_needed)
else
    if PyCall.python !== nothing
        if isempty(libpython)
            @warn """
            PyCall is already installed.  However, you may have trouble using
            this Python executable because it is statically linked to libpython.

            For more information, see:
                https://pyjulia.readthedocs.io/en/latest/troubleshooting.html

            Python executable:
                $python
            Julia executable:
                $(Base.julia_cmd().exec[1])
            """
            exit(code_no_precompile_needed)
        else
            @info """
            PyCall is already installed but not compatible with this Python
            executable.  Re-building PyCall...
            """
            build_pycall()
        end
    elseif pycall_is_installed
        @info """
        PyCall is already installed but importing it failed.
        Re-building PyCall may fix the issue...
        """
        build_pycall()
    else
        @info "Installing PyCall..."
        Pkg.add("PyCall")
    end
end

modpath = Base.locate_package(pkgid)
pkgdir = joinpath(dirname(modpath), "..")

if print_logfile
    logfile = joinpath(pkgdir, "deps", "build.log")
    if isfile(logfile)
        @info "Build log in $logfile"
        print(stderr, read(logfile, String))
    end
end

max_stars_repo_name: grahamgill/pyjulia

max_issues_repo_name: grahamgill/pyjulia

max_forks_repo_name: grahamgill/pyjulia


text:
"""
    get_values(gradient)

Get q-values and b-values
"""
function get_values(gradient)
    nsequence = length(gradient.sequences)
    if gradient.values_type == "q"
        qvalues = repeat(gradient.values, 1, nsequence)
        bvalues =
            gradient.values .^ 2 .* bvalue_no_q.(gradient.sequences)'
    else
        bvalues = repeat(gradient.values, 1, nsequence)
        qvalues = .√(gradient.values ./ bvalue_no_q.(gradient.sequences)')
    end

    qvalues, bvalues
end


text:
import gmsh

gmsh.initialize()

gmsh.model.add("t1")

lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)

gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)

p4 = gmsh.model.geo.addPoint(0, .3, 0, lc)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, p4, 3)
gmsh.model.geo.addLine(4, 1, p4)

gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)

gmsh.model.geo.addPlaneSurface([1], 1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(0, [1, 2], 1)
gmsh.model.addPhysicalGroup(1, [1, 2], 2)
gmsh.model.addPhysicalGroup(2, [1], 6)

gmsh.model.setPhysicalName(2, 6, "My surface")

gmsh.model.mesh.generate(2)

gmsh.write("t1.msh")

gmsh.finalize()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
"""
    download(
        npd :: NASAPrecipitationDataset,
        geo :: GeoRegion = GeoRegion("GLB")
    ) -> nothing

Downloads a NASA Precipitation dataset specified by `npd` for a GeoRegion specified by `geo`

Arguments
=========
- `npd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range
- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat
"""
function download(
	npd :: IMERGHalfHourly{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the IMERG GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("IMERG")
	else
		isinGeoRegion(geo,GeoRegion("IMERG"))
	end

	fnc  = imergrawfiles()
	lon,lat = gpmlonlat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,48)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Day(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(dt) ..."

		ymdfnc = Dates.format(dt,dateformat"yyyymmdd")
		npddir = joinpath(npd.hroot,"$(year(dt))",@sprintf("%03d",dayofyear(dt)))
		for it = 1 : 48

			@debug "$(modulelog()) - Loading data into temporary array for timestep $(fnc[it])"

			npdfnc = "$(npd.fpref).$ymdfnc-$(fnc[it]).$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitationCal"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitationCal"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitationCal"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"
			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,it] = varii / 3600
				else; var[ilon,ilat,it] = NaN32
				end
			end
		end

		save(var,dt,npd,geo,ginfo)

	end

end

function download(
	npd :: IMERGDaily{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the IMERG GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("IMERG")
	else
		isinGeoRegion(geo,GeoRegion("IMERG"))
	end

	lon,lat = gpmlonlat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,31)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Month(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(year(dt))-$(month(dt)) ..."

		npddir = joinpath(npd.hroot,"$(yrmo2dir(dt))")
		ndy = daysinmonth(dt)

		for dy in 1 : ndy
			dtii   = Date(year(dt),month(dt),dy)
			ymdfnc = "$(ymd2str(dtii))"
			npdfnc = "$(npd.fpref).$ymdfnc-$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitationCal"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitationCal"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitationCal"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"

			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,dy] = varii / 86400
				else; var[ilon,ilat,dy] = NaN32
				end
			end

		end

		save(view(var,:,:,1:ndy),dt,npd,geo,ginfo)

	end

end

function download(
	npd :: IMERGMonthly{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the IMERG GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("IMERG")
	else
		isinGeoRegion(geo,GeoRegion("IMERG"))
	end

	lon,lat = gpmlonlat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,12)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Year(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(year(dt)) ..."

		npddir = joinpath(npd.hroot,"$(year(dt))")

		for mo in 1 : 12
			dtii   = Date(year(dt),mo,1)
			ymdfnc = "$(ymd2str(dtii))-S000000-E235959"
			npdfnc = "$(npd.fpref).$ymdfnc.$(@sprintf("%02d",mo)).$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitation"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitation"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitation"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"

			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,mo] = varii / 3600
				else; var[ilon,ilat,mo] = NaN32
				end
			end

		end

		save(var,dt,npd,geo,ginfo)

	end

end

function download(
	npd :: TRMM3Hourly{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the TRMM GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("TRMM")
	else
		isinGeoRegion(geo,GeoRegion("TRMM"))
	end

	lon,lat = trmmlonloat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,8)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Day(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(dt) ..."

		ymdfnc = Dates.format(dt,dateformat"yyyymmdd")

		for it = 1 : 8

			@debug "$(modulelog()) - Loading data into temporary array for timestep $(fnc[it])"

			if isone(it); di = dt-Day(1)
				  npddir = joinpath(npd.hroot,"$(year(di))",@sprintf("%03d",dayofyear(di)))
			else; npddir = joinpath(npd.hroot,"$(year(dt))",@sprintf("%03d",dayofyear(dt)))
			end

			npdfnc = "$(npd.fpref).$ymdfnc.$((it-1)*3).$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitation"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitation"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitation"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"
			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,it] = varii / 3600
				else; var[ilon,ilat,it] = NaN32
				end
			end
		end

		save(var,dt,npd,geo,ginfo)

	end

end

function download(
	npd :: TRMMDaily{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the TRMM GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("TRMM")
	else
		isinGeoRegion(geo,GeoRegion("TRMM"))
	end

	lon,lat = trmmlonlat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,31)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Month(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(year(dt))-$(month(dt)) ..."

		npddir = joinpath(npd.hroot,"$(yrmo2dir(dt))")
		ndy = daysinmonth(dt)

		for dy in 1 : ndy
			dtii   = Date(year(dt),month(dt),dy)
			ymdfnc = "$(ymd2str(dtii))"
			npdfnc = "$(npd.fpref).$ymdfnc.$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitation"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitation"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitation"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"

			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,dy] = varii / 86400
				else; var[ilon,ilat,dy] = NaN32
				end
			end

		end

		save(view(var,:,:,1:ndy),dt,npd,geo,ginfo)

	end

end

function download(
	npd :: TRMMMonthly{ST,DT},
	geo :: GeoRegion = GeoRegion("GLB")
) where {ST<:AbstractString, DT<:TimeType}

	@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from $(npd.dtbeg) to $(npd.dtend)"

	if geo.regID == "GLB"
		@info "$(modulelog()) - Global dataset request has been detected, switching to the TRMM GeoRegion"
		addNPDGeoRegions(); geo = GeoRegion("TRMM")
	else
		isinGeoRegion(geo,GeoRegion("TRMM"))
	end

	lon,lat = trmmlonlat(); nlon = length(lon); nlat = length(lat)
	ginfo = RegionGrid(geo,lon,lat)

	@info "$(modulelog()) - Preallocating temporary arrays for extraction of $(npd.lname) data for the $(geo.name) GeoRegion from the original gridded dataset"
	glon = ginfo.glon; nglon = length(glon); iglon = ginfo.ilon
	glat = ginfo.glat; nglat = length(glat); iglat = ginfo.ilat
	tmp0 = zeros(Float32,nglat,nglon)
	var  = zeros(Float32,nglon,nglat,12)

	if typeof(geo) <: PolyRegion
		  msk = ginfo.mask
	else; msk = ones(nglon,nglat)
	end

	if iglon[1] > iglon[end]
		shift = true
		iglon1 = iglon[1] : nlon; niglon1 = length(iglon1)
		iglon2 = 1 : iglon[end];  niglon2 = length(iglon2)
		tmp1 = @view tmp0[:,1:niglon1]
		tmp2 = @view tmp0[:,(niglon1+1):niglon2]
	else
		shift = false
		iglon = iglon[1] : iglon[end]
	end

	if iglat[1] > iglat[end]
		iglat = iglat[1] : -1 : iglat[end]
	else
		iglat = iglat[1] : iglat[end]
	end

	for dt in npd.dtbeg : Year(1) : npd.dtend

		@info "$(modulelog()) - Downloading $(npd.lname) data for the $(geo.name) GeoRegion from the NASA Earthdata servers using OPeNDAP protocols for $(year(dt)) ..."

		npddir = joinpath(npd.hroot,"$(year(dt))")

		for mo in 1 : 12
			dtii   = Date(year(dt),mo,1)
			ymdfnc = "$(ymd2str(dtii))"
			npdfnc = "$(npd.fpref).$ymdfnc.$(npd.fsuff)"
			ds = NCDataset(joinpath(npddir,npdfnc))
			if !shift
				NCDatasets.load!(ds["precipitation"].var,tmp0,iglat,iglon,1)
			else
				NCDatasets.load!(ds["precipitation"].var,tmp1,iglat,iglon1,1)
				NCDatasets.load!(ds["precipitation"].var,tmp2,iglat,iglon2,1)
			end
			close(ds)

			@debug "$(modulelog()) - Extraction of data from temporary array for the $(geo.name) GeoRegion"

			for ilat = 1 : nglat, ilon = 1 : nglon
				varii = tmp0[ilat,ilon]
				mskii = msk[ilon,ilat]
				if (varii != -9999.9f0) && !isnan(mskii)
					  var[ilon,ilat,mo] = varii / 3600
				else; var[ilon,ilat,mo] = NaN32
				end
			end

		end

		save(var,dt,npd,geo,ginfo)

	end

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
module TestAqua

using Aqua
using Kaleido
using Test

Aqua.test_all(
    Kaleido;
    ambiguities = false,
    project_extras = true,
    stale_deps = true,
    deps_compat = true,
    project_toml_formatting = true,
)

@testset "Method ambiguity" begin
    # Not including `Base` due to
    # https://github.com/JuliaCollections/DataStructures.jl/pull/511
    Aqua.test_ambiguities(Kaleido)
end

@testset "Compare test/Project.toml and test/environments/main/Project.toml" begin
    @test Text(read(joinpath(@__DIR__, "Project.toml"), String)) ==
          Text(read(joinpath(@__DIR__, "environments", "main", "Project.toml"), String))
end

end  # module


text:
function _add_system_balance_slacks!(
    optimization_container::OptimizationContainer,
    slack_name::String,
    expression::Symbol,
    single_first_axes::Bool = false,
)
    time_steps = model_time_steps(optimization_container)
    expression_array = get_expression(optimization_container, expression)
    var_name_up = make_variable_name(slack_name, SLACK_UP)
    var_name_dn = make_variable_name(slack_name, SLACK_DN)
    single_first_axes && (first_index = [axes(expression_array)[1][1]])
    !single_first_axes && (first_index = axes(expression_array)[1])
    variable_up =
        add_var_container!(optimization_container, var_name_up, first_index, time_steps)
    variable_dn =
        add_var_container!(optimization_container, var_name_dn, first_index, time_steps)
    for ix in first_index, jx in time_steps
        variable_up[ix, jx] = JuMP.@variable(
            optimization_container.JuMPmodel,
            base_name = "$(var_name_up)_{$(ix), $(jx)}",
            lower_bound = 0.0
        )
        variable_dn[ix, jx] = JuMP.@variable(
            optimization_container.JuMPmodel,
            base_name = "$(var_name_dn)_{$(ix), $(jx)}",
            lower_bound = 0.0
        )
        add_to_expression!(expression_array, ix, jx, variable_up[ix, jx], 1.0)
        add_to_expression!(expression_array, ix, jx, variable_dn[ix, jx], -1.0)
        JuMP.add_to_expression!(
            optimization_container.cost_function,
            (variable_dn[ix, jx] + variable_up[ix, jx]) * BALANCE_SLACK_COST,
        )
    end
    return
end

function add_slacks!(
    optimization_container::OptimizationContainer,
    ::Type{CopperPlatePowerModel},
)
    _add_system_balance_slacks!(
        optimization_container,
        ACTIVE_POWER,
        :nodal_balance_active,
        true,
    )
    return
end

function add_slacks!(
    optimization_container::OptimizationContainer,
    ::Type{T},
) where {T <: PM.AbstractActivePowerModel}
    _add_system_balance_slacks!(optimization_container, ACTIVE_POWER, :nodal_balance_active)
    return
end

function add_slacks!(
    optimization_container::OptimizationContainer,
    ::Type{T},
) where {T <: PM.AbstractPowerModel}
    _add_system_balance_slacks!(optimization_container, ACTIVE_POWER, :nodal_balance_active)
    _add_system_balance_slacks!(
        optimization_container,
        REACTIVE_POWER,
        :nodal_balance_reactive,
    )
    return
end


text:
# load code
include("channel.jl")
include("vertex.jl")

abstract type Action end

# get interpolated / extrapolated self energy for general action
function get_Σ(
    w :: Float64,
    m :: Mesh,
    a :: Action
    ) :: Float64

    # init value
    val = 0.0

    # check if in bounds, otherwise extrapolate as 1 / w
    if abs(w) <= m.σ[end]
        p   = get_param(abs(w), m.σ)
        val = sign(w) * (p.lower_weight * a.Σ[p.lower_index] + p.upper_weight * a.Σ[p.upper_index])
    else
        val = m.σ[end] * a.Σ[end] / w
    end

    return val
end

# get interpolated vertex component for general action
function get_Γ_comp(
    comp        :: Int64,
    site        :: Int64,
    bs          :: Buffer,
    bt          :: Buffer,
    bu          :: Buffer,
    r           :: Reduced_lattice,
    a           :: Action,
    apply_flags :: Function
    ;
    ch_s        :: Bool = true,
    ch_t        :: Bool = true,
    ch_u        :: Bool = true
    )           :: Float64

    # init with bare value
    val = a.Γ[comp].bare[site]

    # add s channel
    if ch_s
        # check for site exchange
        site_s = site

        if bs.exchange_flag
            site_s = r.exchange[site_s]
        end

        # apply other flags
        comp_s        = comp
        sgn_s, comp_s = apply_flags(bs, comp_s)

        # check for mapping to u channel and interpolate
        if bs.map_flag
            val += sgn_s * get_vertex(site_s, bs, a.Γ[comp_s], 3)
        else
            val += sgn_s * get_vertex(site_s, bs, a.Γ[comp_s], 1)
        end
    end

    # add t channel
    if ch_t
        # check for site exchange
        site_t = site

        if bt.exchange_flag
            site_t = r.exchange[site_t]
        end

        # apply other flags
        comp_t        = comp
        sgn_t, comp_t = apply_flags(bt, comp_t)

        # interpolate
        val += sgn_t * get_vertex(site_t, bt, a.Γ[comp_t], 2)
    end

    # add u channel
    if ch_u
        # check for site exchange
        site_u = site

        if bu.exchange_flag
            site_u = r.exchange[site_u]
        end

        # apply other flags
        comp_u        = comp
        sgn_u, comp_u = apply_flags(bu, comp_u)

        # check for mapping to s channel and interpolate
        if bu.map_flag
            val += sgn_u * get_vertex(site_u, bu, a.Γ[comp_u], 1)
        else
            val += sgn_u * get_vertex(site_u, bu, a.Γ[comp_u], 3)
        end
    end

    return val
end

# get interpolated vertex component on all lattice sites
function get_Γ_comp_avx!(
    comp        :: Int64,
    r           :: Reduced_lattice,
    bs          :: Buffer,
    bt          :: Buffer,
    bu          :: Buffer,
    a           :: Action,
    apply_flags :: Function,
    temp        :: SubArray{Float64, 1, Array{Float64, 3}}
    ;
    ch_s        :: Bool = true,
    ch_t        :: Bool = true,
    ch_u        :: Bool = true
    )           :: Nothing

    # init with bare value
    @turbo temp .= a.Γ[comp].bare

    # add s channel
    if ch_s
        # apply flags
        comp_s        = comp
        sgn_s, comp_s = apply_flags(bs, comp_s)

        # check for mapping to u channel and interpolate
        if bs.map_flag
            get_vertex_avx!(r, bs, a.Γ[comp_s], 3, temp, bs.exchange_flag, sgn_s)
        else
            get_vertex_avx!(r, bs, a.Γ[comp_s], 1, temp, bs.exchange_flag, sgn_s)
        end
    end

    # add t channel
    if ch_t
        # apply flags
        comp_t        = comp
        sgn_t, comp_t = apply_flags(bt, comp_t)

        # interpolate
        get_vertex_avx!(r, bt, a.Γ[comp_t], 2, temp, bt.exchange_flag, sgn_t)
    end

    # add u channel
    if ch_u
        # apply flags
        comp_u        = comp
        sgn_u, comp_u = apply_flags(bu, comp_u)

        # check for mapping to s channel and interpolate
        if bu.map_flag
            get_vertex_avx!(r, bu, a.Γ[comp_u], 1, temp, bu.exchange_flag, sgn_u)
        else
            get_vertex_avx!(r, bu, a.Γ[comp_u], 3, temp, bu.exchange_flag, sgn_u)
        end
    end

    return nothing
end

# load saving and reading for channels and vertices
include("disk.jl")

# load specialized code for different symmetries
include("action_lib/action_su2.jl")  ; include("checkpoint_lib/checkpoint_su2.jl")
include("action_lib/action_u1_dm.jl"); include("checkpoint_lib/checkpoint_u1_dm.jl")

# interface function to replace action with another action (except for bare)
function replace_with!(
    a1 :: Action,
    a2 :: Action
    )  :: Nothing

    # replace self energy
    @turbo a1.Σ .= a2.Σ

    # replace vertices
    for i in eachindex(a1.Γ)
        replace_with!(a1.Γ[i], a2.Γ[i])
    end

    return nothing
end

# interface function to replace action with another action only on the vertex level (except for bare)
function replace_with_Γ!(
    a1 :: Action,
    a2 :: Action
    )  :: Nothing

    # replace vertices
    for i in eachindex(a1.Γ)
        replace_with!(a1.Γ[i], a2.Γ[i])
    end

    return nothing
end

# interface function to multiply action with factor (except for bare)
function mult_with!(
    a   :: Action,
    fac :: Float64
    )   :: Nothing

    # multiply self energy
    @turbo a.Σ .*= fac

    # multiply vertices
    for i in eachindex(a.Γ)
        mult_with!(a.Γ[i], fac)
    end

    return nothing
end

# interface function to multiply action with factor only on the vertex level (except for bare)
function mult_with_Γ!(
    a   :: Action,
    fac :: Float64
    )   :: Nothing

    # multiply vertices
    for i in eachindex(a.Γ)
        mult_with!(a.Γ[i], fac)
    end

    return nothing
end

# interface function to reset an action to zero (except for bare)
function reset!(
    a :: Action
    ) :: Nothing

    mult_with!(a, 0.0)

    return nothing
end

# interface function to reset an action to zero only on the vertex level (except for bare)
function reset_Γ!(
    a :: Action
    ) :: Nothing

    mult_with_Γ!(a, 0.0)

    return nothing
end

# interface function to multiply action with some factor and add to other action (except for bare)
function mult_with_add_to!(
    a2  :: Action,
    fac :: Float64,
    a1  :: Action
    )   :: Nothing

    # multiply add for the self energy
    @turbo a1.Σ .+= fac .* a2.Σ

    # multiply add for the vertices
    for i in eachindex(a1.Γ)
        mult_with_add_to!(a2.Γ[i], fac, a1.Γ[i])
    end

    return nothing
end

# interface function to multiply action with some factor and add to other action only on the vertex level (except for bare)
function mult_with_add_to_Γ!(
    a2  :: Action,
    fac :: Float64,
    a1  :: Action
    )   :: Nothing

    # multiply add for the vertices
    for i in eachindex(a1.Γ)
        mult_with_add_to!(a2.Γ[i], fac, a1.Γ[i])
    end

    return nothing
end

# interface function to add two actions (except for bare)
function add_to!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to!(a2, 1.0, a1)

    return nothing
end

# interface function to add two actions only on the vertex level (except for bare)
function add_to_Γ!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to_Γ!(a2, 1.0, a1)

    return nothing
end

# interface function to subtract two actions (except for bare)
function subtract_from!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to!(a2, -1.0, a1)

    return nothing
end

# interface function to subtract two actions only on the vertex level (except for bare)
function subtract_from_Γ!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to_Γ!(a2, -1.0, a1)

    return nothing
end

"""
    get_abs_max(
        a :: Action
        ) :: Float64

Returns maximum absolute value of an action across all vertex components.
"""
function get_abs_max(
    a :: Action
    ) :: Float64

    abs_max_Γ = zeros(Float64, length(a.Γ))

    for i in eachindex(a.Γ)
        abs_max_Γ[i] = get_abs_max(a.Γ[i])
    end

    abs_max = maximum(abs_max_Γ)

    return abs_max
end

# scan cut through channel, where x is assumed to be generated by the get_mesh function with linear fraction p0
# in order to optimize the frequency mesh we use the following heuristics:
#   a) the linear extent δ should be large enough to capture sharp peaks at finite frequencies, that is, δ >= 'point where channel has decayed to 10 percent of peak value'
#   b) the linear extent δ should be large enough to capture sign changes at finite frequencies, that is, δ >= zero crossing
#   c) the linear extent δ should be small enough to capture sharp peaks at the origin, that is, Δ should fulfill p1 < Δ < p2,
#      where Δ is the relative deviation between the values at the origin and the first finite frequency
#   d) the linear extent δ should be bounded to avoid overambitious shrinking or broadening of the mesh, that is, num_lin * p3 <= δ <= p4,
#      where num_lin is the number of linear frequencies
function scan(
    Λ  :: Float64,
    x  :: Vector{Float64},
    y  :: Vector{Float64},
    p0 :: Float64,
    p1 :: Float64,
    p2 :: Float64,
    p3 :: Float64,
    p4 :: Float64
    )  :: Float64

    # determine current mesh layout
    num_lin = ceil(Int64, p0 * (length(x) - 1))

    # a) check for peaks at finite frequencies and set δa as 'point where channel has decayed to 10 percent of peak value'
    δa  = Inf
    idx = argmax(abs.(y))

    if idx > 1
        for i in idx + 1 : length(x)
            if abs(y[i] / y[idx]) <= 0.1
                δa = x[i]
                break
            end
        end
    end

    # b) check for sign change at finite frequencies and set δb accordingly
    δb = Inf

    for i in 2 : length(y)
        if sign(y[i]) != sign(y[1])
            δb = x[i]
            break
        end
    end

    # c) check if the value at the origin is finite. If so, set δc such that p1 < Δ < p2
    δc = num_lin * x[2]

    if abs(y[1]) / abs(y[idx]) > 1e-3
        # compute Δ
        Δ = abs(y[2] - y[1]) / max(abs(y[2]), abs(y[1]))

        # adjust δc if Δ is out of required bounds
        while (p1 < Δ < p2) == false
            # if Δ is too large decrease δc by two percent
            if Δ >= p2
                δc *= 0.98
            # if Δ is too small increase δc by two percent
            elseif Δ <= p1
                δc *= 1.02
            end

            # ensure δc is within required bounds
            if δc < num_lin * p3
                break
            elseif δc > p4
                break
            end

            # generate new reference data
            xp = get_mesh(δc, x[end], length(x) - 1, p0)
            yp = similar(y)

            for i in eachindex(yp)
                p     = get_param(xp[i], x)
                yp[i] = p.lower_weight * y[p.lower_index] + p.upper_weight * y[p.upper_index]
            end

            # recompute Δ
            Δ = abs(yp[2] - yp[1]) / max(abs(yp[2]), abs(yp[1]))
        end
    end

    # d) merge results from heuristics a), b) and c) and perform sanity checks
    δ = min(max(num_lin * p3, min(δa, δb, δc)), p4)

    return δ
end

# auxiliary function to scan a single channel
function scan_channel(
    Λ   :: Float64,
    p_Ω :: NTuple{5, Float64},
    p_ν :: NTuple{5, Float64},
    Ω   :: Vector{Float64},
    ν   :: Vector{Float64},
    ch  :: Channel
    )   :: NTuple{2, Float64}

    # deref data
    q3 = ch.q3

    # determine position of the maximum
    idxs = argmax(abs.(q3))

    # get cuts through the maximum
    q3_Ω   = q3[idxs[1], :, idxs[3], idxs[4]]
    q3_ν_1 = Float64[q3[idxs[1], idxs[2],      x,       x] - q3[idxs[1], idxs[2],     end,     end] for x in eachindex(ν)]
    q3_ν_2 = Float64[q3[idxs[1], idxs[2],      x, idxs[4]] - q3[idxs[1], idxs[2],     end, idxs[4]] for x in eachindex(ν)]
    q3_ν_3 = Float64[q3[idxs[1], idxs[2], idxs[3],      x] - q3[idxs[1], idxs[2], idxs[3],     end] for x in eachindex(ν)]

    # scan bosonic cut
    Ω_lin = scan(Λ, Ω, q3_Ω, p_Ω[1], p_Ω[2], p_Ω[3], p_Ω[4] * Λ, p_Ω[5] * Λ)

    # scan fermionic cuts
    ν_lin_1 = scan(Λ, ν, q3_ν_1, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin_2 = scan(Λ, ν, q3_ν_2, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin_3 = scan(Λ, ν, q3_ν_3, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin   = min(ν_lin_1, ν_lin_2, ν_lin_3)

    return Ω_lin, ν_lin
end

# resample an action to new meshes via scanning and trilinear interpolation
function resample_from_to(
    Λ      :: Float64,
    p_σ    :: NTuple{2, Float64},
    p_Ωs   :: NTuple{5, Float64},
    p_νs   :: NTuple{5, Float64},
    p_Ωt   :: NTuple{5, Float64},
    p_νt   :: NTuple{5, Float64},
    p_χ    :: NTuple{5, Float64},
    lins   :: NTuple{5, Float64},
    bounds :: NTuple{5, Float64},
    m_old  :: Mesh,
    a_old  :: Action,
    a_new  :: Action,
    χ      :: Vector{Matrix{Float64}}
    )      :: Mesh

    # scale linear bounds
    σ_lin  = lins[2] * Λ
    Ωs_lin = lins[3] * Λ
    νs_lin = lins[4] * Λ
    Ωt_lin = lins[3] * Λ
    νt_lin = lins[4] * Λ
    χ_lin  = lins[5] * Λ

    # adjust meshes via scanning once required scale is reached
    if Λ < lins[1]
        # scan self energy
        σ_lin = p_σ[2] * m_old.σ[argmax(abs.(a_old.Σ))]

        # scan the channels
        Ωs_lins, νs_lins = zeros(Float64, length(a_old.Γ)), zeros(length(a_old.Γ))
        Ωt_lins, νt_lins = zeros(Float64, length(a_old.Γ)), zeros(length(a_old.Γ))

        for i in eachindex(a_old.Γ)
            Ωs_lins[i], νs_lins[i] = scan_channel(Λ, p_Ωs, p_νs, m_old.Ωs, m_old.νs, a_old.Γ[i].ch_s)
            Ωt_lins[i], νt_lins[i] = scan_channel(Λ, p_Ωt, p_νt, m_old.Ωt, m_old.νt, a_old.Γ[i].ch_t)
        end

        Ωs_lin, νs_lin = minimum(Ωs_lins), minimum(νs_lins)
        Ωt_lin, νt_lin = minimum(Ωt_lins), minimum(νt_lins)

        # scan the correlations
        χ_lins = zeros(Float64, length(χ))

        for i in eachindex(χ_lins)
            χ_lins[i] = scan(Λ, m_old.χ, χ[i][argmax(abs.(χ[i]))[1], :], p_χ[1], p_χ[2], p_χ[3], p_χ[4] * Λ, p_χ[5] * Λ)
        end

        χ_lin = minimum(χ_lins)
    end

    # build new frequency meshes according to scanning results
    σ     = get_mesh( σ_lin, bounds[2] * max(Λ, bounds[1]), m_old.num_σ - 1, p_σ[1])
    Ωs    = get_mesh(Ωs_lin, bounds[3] * max(Λ, bounds[1]), m_old.num_Ω - 1, p_Ωs[1])
    νs    = get_mesh(νs_lin, bounds[4] * max(Λ, bounds[1]), m_old.num_ν - 1, p_νs[1])
    Ωt    = get_mesh(Ωt_lin, bounds[3] * max(Λ, bounds[1]), m_old.num_Ω - 1, p_Ωt[1])
    νt    = get_mesh(νt_lin, bounds[4] * max(Λ, bounds[1]), m_old.num_ν - 1, p_νt[1])
    χ     = get_mesh( χ_lin, bounds[5] * max(Λ, bounds[1]), m_old.num_χ - 1, p_χ[1])
    m_new = Mesh(m_old.num_σ, m_old.num_Ω, m_old.num_ν, m_old.num_χ, σ, Ωs, νs, Ωt, νt, χ)

    # resample self energy
    for w in eachindex(m_new.σ)
        a_new.Σ[w] = get_Σ(m_new.σ[w], m_old, a_old)
    end

    # resample vertices
    for i in eachindex(a_new.Γ)
        resample_from_to!(m_old, a_old.Γ[i], m_new, a_new.Γ[i])
    end

    return m_new
end

# generate action dummy
function get_action_empty(
    symmetry :: String,
    r        :: Reduced_lattice,
    m        :: Mesh
    ;
    S        :: Float64 = 0.5
    )        :: Action

    if symmetry == "su2"
        return get_action_su2_empty(S, r, m)
    elseif symmetry == "u1-dm"
        return get_action_u1_dm_empty(r, m)
    end
end

"""
    read_checkpoint(
        file     :: HDF5.File,
        Λ        :: Float64
        )        :: Tuple{Float64, Float64, Mesh, Action}

Read checkpoint of FRG calculation from HDF5 file.
Returns cutoff Λ, ODE stepwidth dΛ, frequency meshes (wrapped in Mesh struct) and vertices (wrapped in Action struct).
"""
function read_checkpoint(
    file     :: HDF5.File,
    Λ        :: Float64
    )        :: Tuple{Float64, Float64, Mesh, Action}

    # read symmetry group from file
    symmetry = read(file, "symmetry")

    if symmetry == "su2"
        return read_checkpoint_su2(file, Λ)
    elseif symmetry == "u1-dm"
        return read_checkpoint_u1_dm(file, Λ)
    end
end

# load tests and timers
include("test.jl")
include("timers.jl")


text:
"""
    SIS

A type representing a Susceptible-Infected-Susceptible model with infection rate `β` and recovery rate `γ`.
"""
struct SIS <: AbstractCompartimentalModel
    β::Real
    γ::Real
end

SI(sis::SIS) = SI(sis.β)

num_states(::SIS) = 2

# Metapopulation

function rate(i, state, mp::Metapopulation{SIS})
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    h = mp.h
    β*state[i,1]*state[i,2] + γ*state[i,2] + sum(D .* state[i,:]) * (outdegree(h,i) > 0)
end

function update_state_and_rates!(state, a, k, mp::Metapopulation{SIS})
    N = nv(mp.h)
    D = mp.D
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    h = mp.h
    p = vcat(β*state[k,1]*state[k,2], γ*state[k,2], D.*state[k,:]*(outdegree(h,k)>0))
    j = sample(ProbabilityWeights(p, a[k]))
    if j == 1
        state[k,1] -= 1
        state[k,2] += 1
        a[k] = rate(k, state, mp)
    elseif j == 2
        state[k,1] += 1
        state[k,2] -= 1
        a[k] = rate(k, state, mp)
    else
        μ = k
        ν = rand(outneighbors(h, μ))
        state[μ,j-2] -= 1 # NB. j-2 is the state of the individual that migrates
        state[ν,j-2] += 1 # NB. because the first 2 elements of p are the reactions
        a[μ] = rate(μ, state, mp)
        a[ν] = rate(ν, state, mp)
    end
end

function meanfield_fun(mp::Metapopulation{SIS})
    L = normalized_laplacian(mp.h)
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    f! = function(dx, x, p, t)
        tmp = (β.*x[:,1] .- γ).*x[:,2]
        dx[:,1] .= .-tmp .- D[1].*(L'*x[:,1])
        dx[:,2] .= tmp .- D[2].*(L'x[:,2])
    end
    return f!
end

# Contact Process

function rate(k, state, cp::ContactProcess{SIS})
    if state[k] == 1
        l = length(filter(i->state[i]==2, inneighbors(cp.g, k)))
        return cp.dynamics.β * l
    else
        return cp.dynamics.γ
    end
end

function update_state_and_rates!(state, a, k, cp::ContactProcess{SIS})
    g = cp.g
    β = cp.dynamics.β
    if state[k] == 1
        state[k] = 2
        for i in Iterators.filter(j->state[j]==1, outneighbors(g, k))
            a[i] += β
        end
    else
        state[k] = 1
        for i in Iterators.filter(j->state[j]==1, outneighbors(g, k))
            a[i] -= β
        end
    end
    a[k] = rate(k, state, cp)
end

function update_output!(output, state, n, k, cp::ContactProcess{SIS})
    if state[k] == 2 # newly infected node
        output[n] = output[n-1] + [-1, 1]
    else
        output[n] = output[n-1] + [1, -1]
    end
end

function init_state_mf(cp::ContactProcess{SIS}, x0)
    return float(x0 .- 1)
end

function meanfield_fun(cp::ContactProcess{SIS})
    A = adjacency_matrix(cp.g)
    β = cp.dynamics.β
    γ = cp.dynamics.γ
    f! = function(dx, x, p, t)
        dx .= β.*(1.0 .- x).*(A*x) .- γ.*x
    end
    return f!
end

function finalize_meanfield(cp::ContactProcess{SIS}, sol)
    N = nv(cp.g)
    n = length(sol.u)
    output = Array{Float64,2}(undef, n, N)
    for i in 1:n
        output[i,:] .= sol.u[i]
    end
    return sol.t, output
end

# Metaplex

function rate(k, state, mpx::Metaplex{SIS})
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    od = outdegree(mpx.h, x_μ[k]) > 0
    if x_i[k] == 1
        l = length(filter(i->x_i[i]==2 && x_μ[i]==x_μ[k], inneighbors(mpx.g, k)))
        return β*l + D[1]*od
    else
        return γ + D[2]*od
    end
end

function update_state_and_rates!(state, a, k, mpx::Metaplex{SIS})
    N = nv(mpx.g)
    M = nv(mpx.h)
    g = mpx.g
    h = mpx.h
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    popcounts = state[:popcounts]
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    μ = x_μ[k]
    if x_i[k] == 1
        if rand()*a[k] < D[1] # susceptible node migrates
            ν = rand(outneighbors(h, μ))
            x_μ[k] = ν
            popcounts[μ,1] -= 1
            popcounts[ν,1] += 1
            a[k] = rate(k, state, mpx)
        else # susceptible node becomes infected
            x_i[k] = 2
            popcounts[μ,1] -= 1
            popcounts[μ,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g,k))
                a[i] += β
            end
        end
    else
        if rand()*a[k] < D[2] # infected node migrates
            ν = rand(outneighbors(h,μ))
            x_μ[k] = ν
            popcounts[μ,2] -= 1
            popcounts[ν,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1, outneighbors(g,k))
                if x_μ[i] == μ
                    a[i] -= β
                elseif x_μ[i] == ν
                    a[i] += β
                end
            end
        else # infected node recovers
            x_i[k] = 1
            popcounts[μ,2] -= 1
            popcounts[μ,1] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g,k))
                a[i] -= β
            end
        end
    end
end

function meanfield_fun(mpx::Metaplex{SIS})
    N = nv(mpx.g)
    M = nv(mpx.h)
    A = adjacency_matrix(mpx.g)
    L = normalized_laplacian(mpx.h)
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    f! = function(dx, x, p, t)
        for i in 1:N
            dx[1,i,:] .= .-D[1]*L'*x[1,i,:]
            dx[2,i,:] .= .-D[2]*L'*x[2,i,:]
        end
        for μ in 1:M
            infection = β.*x[1,:,μ].*(A*x[2,:,μ]) .- γ.*x[2,:,μ]
            dx[1,:,μ] .-= infection
            dx[2,:,μ] .+= infection
        end
    end
    return f!
end

# Heterogeneous Metaplex

function rate(k, state, mpx::HeterogeneousMetaplex{SIS})
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    μ = x_μ[k]
    od = outdegree(mpx.h, x_μ[k]) > 0
    if x_i[k] == 1
        l = length(filter(i->x_i[i]==2 && x_μ[i]==x_μ[k], inneighbors(mpx.g[μ], k)))
        return β*l + D[1]*od
    else
        return γ + D[2]*od
    end
end

function update_state_and_rates!(state, a, k, mpx::HeterogeneousMetaplex{SIS})
    g = mpx.g
    h = mpx.h
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    popcounts = state[:popcounts]
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    μ = x_μ[k]
    if x_i[k] == 1
        if rand()*a[k] < D[1] # susceptible node migrates
            ν = rand(outneighbors(h, μ))
            x_μ[k] = ν
            popcounts[μ,1] -= 1
            popcounts[ν,1] += 1
            a[k] = rate(k, state, mpx)
        else # susceptible node becomes infected
            x_i[k] = 2
            popcounts[μ,1] -= 1
            popcounts[μ,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g[μ],k))
                a[i] += β
            end
        end
    else
        if rand()*a[k] < D[2] # infected node migrates
            ν = rand(outneighbors(h,μ))
            x_μ[k] = ν
            popcounts[μ,2] -= 1
            popcounts[ν,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1, outneighbors(g[μ],k))
                if x_μ[i] == μ
                    a[i] -= β
                elseif x_μ[i] == ν
                    a[i] += β
                end
            end
        else # infected node recovers
            x_i[k] = 1
            popcounts[μ,2] -= 1
            popcounts[μ,1] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g[μ],k))
                a[i] -= β
            end
        end
    end
end

function meanfield_fun(mpx::HeterogeneousMetaplex{SIS})
    N = nv(mpx.g[1])
    M = nv(mpx.h)
    A = adjacency_matrix.(mpx.g)
    L = normalized_laplacian(mpx.h)
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    f! = function(dx, x, p, t)
        for i in 1:N
            dx[1,i,:] .= .-D[1]*L'*x[1,i,:]
            dx[2,i,:] .= .-D[2]*L'*x[2,i,:]
        end
        for μ in 1:M
            infection = β.*x[1,:,μ].*(A[μ]*x[2,:,μ]) .- γ.*x[2,:,μ]
            dx[1,:,μ] .-= infection
            dx[2,:,μ] .+= infection
        end
    end
    return f!
end

# HeterogeneousMetapopulation

function rate(i, state, mp::HeterogeneousMetapopulation{SIS})
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    h = mp.h
    ks = mp.ks
    N = mp.N
    βᵢ = β * ks[i] / N
    return βᵢ*state[i,1]*state[i,2] + γ*state[i,2] + sum(D .* state[i,:]) * (outdegree(h,i) > 0)
end

function update_state_and_rates!(state, a, k, mp::HeterogeneousMetapopulation{SIS})
    M = nv(mp.h)
    N = mp.N
    ks = mp.ks
    D = mp.D
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    h = mp.h
    βₖ = β * ks[k] / N
    p = vcat(βₖ*state[k,1]*state[k,2], γ*state[k,2], D.*state[k,:]*(outdegree(h,k)>0))
    j = sample(ProbabilityWeights(p, a[k]))
    if j == 1
        state[k,1] -= 1
        state[k,2] += 1
        a[k] = rate(k, state, mp)
    elseif j == 2
        state[k,1] += 1
        state[k,2] -= 1
        a[k] = rate(k, state, mp)
    else
        μ = k
        ν = rand(outneighbors(h, μ))
        state[μ,j-2] -= 1 # NB. j-2 is the state of the individual that migrates
        state[ν,j-2] += 1 # NB. because the first 2 elements of p are the reactions
        a[μ] = rate(μ, state, mp)
        a[ν] = rate(ν, state, mp)
    end
end

function meanfield_fun(mp::HeterogeneousMetapopulation{SIS})
    #L = normalized_laplacian(mp.h)
    L = laplacian_matrix(mp.h)
    ks = mp.ks
    N = mp.N
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    βs = β * ks / N
    f! = function(dx, x, p, t)
        tmp = (βs .* x[:,1] .- γ).*x[:,2]
        dx[:,1] .= .-tmp .- D[1].*(L'*x[:,1])
        dx[:,2] .= tmp .- D[2].*(L'x[:,2])
    end
    return f!
end


text:
import ..UncertainDatasets:
	AbstractUncertainValueDataset,
	UncertainValueDataset,
	UVAL_COLLECTION_TYPES

import ..SamplingConstraints:
	NoConstraint,
	TruncateLowerQuantile,
	TruncateUpperQuantile,
	TruncateQuantiles

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint) -> Vector{T} where T

Resample `x` (a collection of uncertain values) once, applying the provided sampling `constraint`.

Returns a `length(x)`-element vector. The `i`-th element of this vector is generated by
truncating the `i`-th uncertain value by the sampling `constraint`, then drawing a single random
number from the truncated value.

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

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate each distribution at +- 0.5 standard deviations, then resample.
resample(uvals, TruncateStd(0.5))
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint)
	[resample(uv.values[i], constraint) for i in 1:length(uv)]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}) -> Vector{T} where T

Resample `x` (a collection of uncertain values) once, applying the provided sampling `constraint`s.
The number of constraints must match the number of elements in `x`.

Returns a `length(x)`-element vector. The `i`-th element of this vector is generated by
truncating the `i`-th uncertain value by the `i`-th sampling `constraint`, then drawing a single random
number from the truncated value.

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

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate each distribution at +- 0.5 standard deviations, then resample.
resample(uvals, TruncateStd(0.5))
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint})
    [resample(uv.values[i], constraint[i]) for i in 1:length(uv)]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) `n` times, applying the provided sampling `constraint`.

Returns an `n`-element vector of `length(x)`-element vectors. Each of these vectors is an independent
draw from `x`. The `i`-th element of each draw is generated by truncating the `i`-th uncertain value by
the sampling `constraint`, then drawing a single random number from the truncated value.

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

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate the first 50 elements at the 90th percentile range, and the
# last 50 elements at the 40th percentile range.
constraints = [i <= 50 ? TruncateQuantiles(0.05, 0.95) : TruncateQuantiles(0.3, 0.7) for i = 1:100]

# Truncate the distributions, then draw ten independent realisations of the collection subject
# to the provided constraints.
resample(uvals, constraints, 10)
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int)
	[[resample(uv.values[i], constraint) for i in 1:length(uv)] for k = 1:n]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) `n` times, applying the provided sampling `constraint`s.

Returns an `n`-element vector of `length(x)`-element vectors. Each of these vectors is an independent
draw from `x`. The `i`-th element of each draw is generated by truncating the `i`-th uncertain value by
the `i`-th sampling `constraint`, then drawing a single random number from the truncated value.

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

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate the first 50 elements at `± 0.5` standard deviations, and the
# last 50 elements at `± 1.2` standar deviations.
constraints = [i <= 50 ? TruncateStd(0.5) : TruncateStd(1.2) for i = 1:100]

# Apply the constraints element-wise, then draw ten independent realisations
# of the collection subject to those constraints.
resample(uvals, constraints, 10)
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}, n::Int)
	[[resample(uv.values[i], constraint[i]) for i in 1:length(uv)] for k = 1:n]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) once by drawing a single random number from
each of the uncertain values in `x`.

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

## Example

```julia
# Generate some uncertain values represented by gamma distributions
uvals = [UncertainValue(Gamma(i, rand())) for i = 1:100]

# Resample the collection once
resample(uvals)
```
"""
function resample(uvd::UVAL_COLLECTION_TYPES)
	L = length(uvd)
	[resample(uvd.values[i]) for i in 1:L]
end

"""
	resample(uvd::UVAL_COLLECTION_TYPES, n::Int) -> Vector{Vector{T}}

Draw `n` realisations of an uncertain value dataset according to the distributions
of the uncertain values comprising it.

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

## Example

```julia
# Generate some uncertain values represented by gamma distributions
uvals = [UncertainValue(Gamma(i, rand())) for i = 1:100]

# Resample the collection once
resample(uvals)
```
"""
function resample(uvd::UVAL_COLLECTION_TYPES, n::Int)
	L = length(uvd)
	[[resample(uvd.values[i]) for i in 1:L] for k in 1:n]
end


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
include("faberschauder.jl")
using ZigZagBoomerang, SparseArrays, LinearAlgebra, Revise
# Drift
α = 1.5
β = 0.1
b(x) = -β*x + α*sin(2pi*x)
# First derivative
b′(x) = -β + 2*α*2pi*cos(2pi*x)
# Second derivative
b″(x) = -2*α*(2pi)^2*sin(2pi*x)
# length diffusion bridge
T = 2.0

L = 7
n = (2 << L) + 1
ξ0 = 0randn(n)
u, v = 1.5, -0.5  # initial and fianl point
ξ0[1] = u/sqrt(T)
ξ0[end] = v/sqrt(T)
c = ones(n)
c[end] = c[1] = 0.0
θ0 = rand((-1.0, 1.0), n)
θ0[end] = θ0[1] =  0.0 # fix final point
T′ = 5000.0 # final clock of the pdmp

Γ = sparse(1.0I, n, n)
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ!, 0.0, ξ0, θ0, T, 10.0, Boomerang(Γ, ξ0*0, 0.1; ρ=0.9), 1, L, adapt=false);
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ, 0.0, ξ0, rand((-1.0, 1.0), n), T, 40.0*ones(n), ZigZag(Γ, ξ0*0), 5, L, adapt=false);
trace, (t, ξ, θ), (acc, num), c = @time spdmp(∇ϕmoving, 0.0, ξ0,
    θ0, T′, c, ZigZag(Γ, ξ0*0), SelfMoving(), L, T, adapt=true);
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ, 0.0, ξ0, rand((-1.0,1.0), n), T, 100.0*ones(n), FactBoomerang(Γ, ξ0*0, 0.1), 5, L, adapt=false);
ts, ξs = splitpairs(discretize(trace, T′/n))
S = T*(0:n)/(n+1)

### change to false is you want to Plot the results
if true
    error("Stop before plotting the results...")
end

using CairoMakie
p1 = lines(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
for ξ in ξs[1:5:end]
    lines!(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
end
display(p1)
save("./diffbridges.png", p1)


text:
using Documenter, MicrobiomePlots

makedocs(
    format = :html,
    sitename = "MicrobiomePlots.jl",
    pages = [
        "Home" => "index.md",
        "Recipes" => "recipes.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD"
)

deploydocs(
    repo = "github.com/BioJulia/MicrobiomePlots.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
Base.:(==)(t1::CXType, t2::CXType) = clang_equalTypes(t1, t2) != 0

"""
    getCanonicalType(t::CXType) -> CXType
    getCanonicalType(t::CLType) -> CLType
Return the canonical type for a CXType.

Clang's type system explicitly models typedefs and all the ways a specific type can be
represented. The canonical type is the underlying type with all the "sugar" removed.
For example, if 'T' is a typedef for 'int', the canonical type for 'T' would be 'int'.
Wrapper for libclang's [`clang_getCanonicalType`](@ref).
"""
getCanonicalType(t::CXType) = clang_getCanonicalType(t)
getCanonicalType(t::CLType)::CLType = clang_getCanonicalType(t)

"""
    isConstQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "const" qualifier set, without looking through typedefs
that may have added "const" at a different level.
Wrapper for libclang's [`clang_isConstQualifiedType`](@ref).
"""
isConstQualifiedType(t::Union{CXType,CLType})::Bool = clang_isConstQualifiedType(t)

"""
    isVolatileQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "volatile" qualifier set, without looking through typedefs
that may have added "volatile" at a different level.
Wrapper for libclang's [`clang_isVolatileQualifiedType`](@ref).
"""
isVolatileQualifiedType(t::CXType)::Bool = clang_isVolatileQualifiedType(t)

"""
    isRestrictQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "restrict" qualifier set, without looking through typedefs
that may have added "restrict" at a different level.
Wrapper for libclang's [`clang_isRestrictQualifiedType`](@ref).
"""
isRestrictQualifiedType(t::Union{CXType,CLType})::Bool = clang_isRestrictQualifiedType(t)

"""
    getAddressSpace(t::Union{CXType,CLType})
Returns the address space of the given type.
Wrapper for libclang's [`clang_getAddressSpace`](@ref).
"""
getAddressSpace(t::Union{CXType,CLType}) = clang_getAddressSpace(t)

"""
    getTypedefName(t::Union{CXType,CLType}) -> String
Return the typedef name of the given type.
Wrapper for libclang's [`clang_getTypedefName`](@ref).
"""
function getTypedefName(t::Union{CXType,CLType})
    cxstr = clang_getTypedefName(t)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getPointeeType(t::CXType) -> CXType
    getPointeeType(t::CLType) -> CLType
Return the type of the pointee for pointer types.
Wrapper for libclang's [`clang_getPointeeType`](@ref).
"""
getPointeeType(t::CXType) = clang_getPointeeType(t)
getPointeeType(t::CLType)::CLType = clang_getPointeeType(t)

"""
    getTypeDeclaration(t::CXType) -> CXCursor
    getTypeDeclaration(t::CLType) -> CLCursor
Return the cursor for the declaration of the given type. To get the type of the cursor,
see [`type`](@ref). Wrapper for libclang's [`clang_getTypeDeclaration`](@ref).
"""
getTypeDeclaration(t::CXType) = clang_getTypeDeclaration(t)
getTypeDeclaration(t::CLType)::CLCursor = clang_getTypeDeclaration(t)

"""
    getFunctionTypeCallingConv(t::Union{CXType,CLFunctionNoProto,CLFunctionProto}) -> CXCallingConv
Return the calling convention associated with a function type.
Wrapper for libclang's [`clang_getFunctionTypeCallingConv`](@ref).
"""
getFunctionTypeCallingConv(t::Union{CXType,CLFunctionNoProto,CLFunctionProto})::CXCallingConv = clang_getFunctionTypeCallingConv(t)

"""
    getCursorResultType(t::CXType) -> CXType
    getCursorResultType(t::Union{CLFunctionNoProto,CLFunctionProto}) -> CLType
Return the return type associated with a function type.
Wrapper for libclang's [`clang_getResultType`](@ref).
"""
getCursorResultType(t::CXType) = clang_getResultType(t)
getCursorResultType(t::Union{CLFunctionNoProto,CLFunctionProto})::CLType = clang_getResultType(t)

## TODO:
# clang_getExceptionSpecificationType

"""
    getNumArguments(t::Union{CXType,CLFunctionNoProto,CLFunctionProto}) -> Int
Return the number of non-variadic parameters associated with a function type.
Wrapper for libclang's [`clang_getNumArgTypes`](@ref).
"""
getNumArguments(t::Union{CXType,CLFunctionNoProto,CLFunctionProto,CLUnexposed})::Int = clang_getNumArgTypes(t)

"""
    getArgType(t::CXType, i::Unsigned) -> CXType
    getArgType(t::Union{CLFunctionNoProto,CLFunctionProto}, i::Integer) -> CLType
Return the type of a parameter of a function type.
Wrapper for libclang's [`clang_getArgType`](@ref).
"""
getArgType(t::CXType, i::Unsigned) = clang_getArgType(t, Unsigned(i))
getArgType(t::Union{CLFunctionNoProto,CLFunctionProto,CLUnexposed}, i::Integer)::CLType = getArgType(t.type, Unsigned(i))

"""
    isVariadic(t::Union{CXType,CLType}) -> Bool
Return true if the CXType is a variadic function type.
Wrapper for libclang's [`clang_isFunctionTypeVariadic`](@ref).
"""
isVariadic(t::Union{CXType,CLType})::Bool = clang_isFunctionTypeVariadic(t)

"""
    isPODType(t::Union{CXType,CLType}) -> Bool
Return true if the CXType is a plain old data type.
Wrapper for libclang's [`clang_isPODType`](@ref).
"""
isPODType(t::CXType)::Bool = clang_isPODType(t)

"""
    getElementType(t::CXType) -> CXType
    getElementType(t::Union{CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray,CLComplex}) -> CLType
Return the element type of an array, complex, or vector type.
Wrapper for libclang's [`clang_getElementType`](@ref).
"""
getElementType(t::CXType) = clang_getElementType(t)
getElementType(t::Union{CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray,CLComplex})::CLType = clang_getElementType(t)
# old API: clang_getArrayElementType

"""
    getNumElements(t::Union{CXType,CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray}) -> Int
Return the number of elements of an array or vector type.
Wrapper for libclang's [`clang_getNumElements`](@ref).
"""
getNumElements(t::Union{CXType,CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray})::Int = clang_getNumElements(t)
# old API: clang_getArraySize

"""
    getNamedType(t::CXType) -> CXType
    getNamedType(t::CLElaborated) -> CLType
Return the type named by the qualified-id.
Wrapper for libclang's [`clang_Type_getNamedType`](@ref).
"""
getNamedType(t::CXType) = clang_Type_getNamedType(t)
getNamedType(t::CLElaborated)::CLType = clang_Type_getNamedType(t)

"""
    getSizeOf(t::CXType) -> Int
    getSizeOf(t::CLType) -> Int
Return the size of a type in bytes as per C++[expr.sizeof] standard,

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getSizeOf(t::CXType)::Int = clang_Type_getSizeOf(t)
getSizeOf(t::CLType) = getSizeOf(t.type)

"""
    getAlignOf(t::CXType) -> Int
    getAlignOf(t::CLType) -> Int
Return the alignment of a type in bytes as per C++[expr.alignof] standard.

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getAlignOf(t::CXType)::Int = clang_Type_getAlignOf(t)
getAlignOf(t::CLType) = getAlignOf(t.type)

"""
    getOffsetOf(t::CXType, s) -> Int
    getOffsetOf(t::CLType, s::AbstractString) -> Int
Return the offset of a field named S in a record of type T in bits as it
would be returned by __offsetof__ as per C++11[18.2p4].

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getOffsetOf(t::CXType, s)::Int = clang_Type_getOffsetOf(t, s)
getOffsetOf(t::CLType, s::AbstractString) = getOffsetOf(t.type, s)

"""
    isInvalid(t::CXType) -> Bool
    isInvalid(t::CLType) -> Bool
Return true if the type is a valid type.
"""
isInvalid(t::CXType) = kind(t) != CXType_Invalid
isInvalid(t::CLType) = isInvalid(t.type)

## TODO:
# clang_Type_isTransparentTagTypedef (ObjectiveC)
# clang_Type_getClassType (C++)
# clang_Type_getNumTemplateArguments (C++)
# clang_Type_getTemplateArgumentAsType (C++)
# clang_Type_getCXXRefQualifier (C++)
#

# helper
"""
    kind(t::CXType) -> CXTypeKind
    kind(t::CLType) -> CXTypeKind
Return the kind of the given type.
"""
kind(t::CXType) = t.kind
kind(t::CLType) = kind(t.type)

"""
    spelling(t::Union{CXType,CLType}) -> String
Pretty-print the underlying type using the rules of the language of the translation unit
from which it came. If the type is invalid, an empty string is returned.
Wrapper for libclang's [`clang_getTypeSpelling`](@ref).
"""
function spelling(t::Union{CXType,CLType})
    cxstr = clang_getTypeSpelling(t)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    spelling(kind::CXTypeKind) -> String
Return the spelling of a given CXTypeKind.
"""
function spelling(kind::CXTypeKind)
    cxstr = clang_getTypeKindSpelling(kind)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

# visitor
function ct_field_visitor(cursor::CXCursor, list)::CXVisitorResult
    push!(list, cursor)
    return CXVisit_Continue
end

"""
    fields(ty::CXType) -> Vector{CXCursor}
    fields(ty::CLType) -> Vector{CLCursor}
Return a child cursor vector of the given cursor.
"""
function fields(ty::CXType)
    list = CXCursor[]
    ct_visitor_cb = @cfunction(
        ct_field_visitor, CXVisitorResult, (CXCursor, Ref{Vector{CXCursor}})
    )
    GC.@preserve list ccall(
        (:clang_Type_visitFields, LibClang.libclang),
        UInt32,
        (CXType, CXFieldVisitor, Any),
        ty,
        ct_visitor_cb,
        list,
    )
    return list
end
fields(ty::CLType)::Vector{CLCursor} = fields(ty.type)

is_elaborated(ty::CXType) = kind(ty) == CXType_Elaborated
is_elaborated(ty::CLType) = is_elaborated(ty.type)

"""
    has_elaborated_reference(ty::CLType) -> Bool
    has_elaborated_reference(ty::CXType) -> Bool
Return true if the type is an elaborated type or the type indirectly refers to an
elaborated type.
"""
function has_elaborated_reference(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return has_elaborated_reference(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return has_elaborated_reference(elty)
    else
        return is_elaborated(ty)
    end
end
has_elaborated_reference(ty::CLType) = has_elaborated_reference(ty.type)

"""
    get_elaborated_cursor(ty::CLType) -> CLCursor
    get_elaborated_cursor(ty::CXType) -> CXCursor
Return the cursor of the elaborated type that is referenced by the input type.
The input type can be a pointer/array.
This function returns [`clang_getNullCursor`](@ref) if the input type is not refer to an
elaborated type.
"""
function get_elaborated_cursor(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return get_elaborated_cursor(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return get_elaborated_cursor(elty)
    elseif is_elaborated(ty)
        return getTypeDeclaration(ty)
    else
        return getNullCursor()
    end
end
get_elaborated_cursor(ty::CLType)::CLCursor = get_elaborated_cursor(ty.type)

function is_function(ty::CXType)
    canonical_ty = kind(ty) == CXType_Typedef ? getCanonicalType(ty) : ty
    k = kind(canonical_ty)
    return k == CXType_FunctionProto || k == CXType_FunctionNoProto
end
is_function(ty::CLType) = is_function(ty.type)

"""
    has_function_reference(ty::CLType) -> Bool
    has_function_reference(ty::CXType) -> Bool
Return true if the type is a function or the type indirectly refers to a function.
"""
function has_function_reference(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return has_function_reference(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return has_function_reference(elty)
    else
        return is_function(ty)
    end
end
has_function_reference(ty::CLType) = has_function_reference(ty.type)


text:
using Gnuplot
x = 0:0.01:2π
files = Vector{String}()
for i in 1:50
    @gp xlab = "x" ylab = "sin"  x sin.(x .+ i / 10.0) "w l lc 'black' notit"
    save(term="pngcairo size 600,400", output="advancedAnim$(i).png")
    push!(files, "advancedAnim$(i).png")
end
run(`convert -limit memory 1 -limit map 1 $files advancedAnim.gif`)
run(`rm -f $files`)

# https://github.com/gcalderone/Gnuplot.jl/issues/9

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
include("./params.jl")
include("./utils.jl")

function compute_optimal_bw{T<:Real}(C::Function, ws::AbstractArray{T,1})::T

    ################# WRITE SOMETHING #################

    return wopt

end

function get_costfunc(data)
    # Cost function

    return C
end

function sskernel(data; ws=nothing)
    # Example usage:
    # x = 0.5-0.5*log(rand(1,1e3));
    # ws = linspace(0.0001, 1.0, 100);
    # rs = sskernel(x, ws=ws); # or simply sskernel(x)
    #
    # Input argument
    # data: Spike-rate.
    # w:    Bandwidth for radial basis function
    #
    # Output argument
    # rs:   Estimated intensity function

    ################# WRITE SOMETHING #################

    return rs
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
using Millboard
using Base.Test

board = ()
@test """
+==+""" == table(board) |> string

board = ([],)
@test """
+===+===+
|   | 1 |
+===+===+
| 1 |   |
+---+---+""" == table(board) |> string

board = ((1,2,3),)
@test """
+===+=========+
|   |       1 |
+===+=========+
| 1 | (1,2,3) |
+---+---------+""" == table(board) |> string

board = ([1,2,3],)
@test """
+===+===+
|   | 1 |
+===+===+
|   | 1 |
|   | 2 |
| 1 | 3 |
+---+---+""" == table(board) |> string

board = [()]
@test """
+===+====+
|   |  1 |
+===+====+
| 1 | () |
+---+----+""" == table(board) |> string

board = [(1,2,3)]
@test """
+===+=========+
|   |       1 |
+===+=========+
| 1 | (1,2,3) |
+---+---------+""" == table(board) |> string

board = [(1,2,3) (4,5,6)]
@test """
+===+=========+=========+
|   |       1 |       2 |
+===+=========+=========+
| 1 | (1,2,3) | (4,5,6) |
+---+---------+---------+""" == table(board) |> string


text:
function expa(stoichiometric_matrix::Array{Float64,2})

    # initialize -
	(ℳ,ℛ) = size(stoichiometric_matrix)
    T = [Matrix{Float64}(I,ℛ,ℛ) transpose(stoichiometric_matrix)];

    # outer loop -
    for col_index = 1:ℳ

        # adjust the col index (from the org matrix -)
	    local_col_index = col_index + ℛ;

        # get the current RHS -
        RHS = T[:,(ℛ + 1):end] # this gets the cols belonging to the stoichiometric matrix -

        # Pick a col -
	    TMP_COL = RHS[:,col_index]; # this is grabbing a metabolite -

        # how many zeros, + and - in this col -
        idx_zero_array = findall(x->x==0.0,TMP_COL)
        idx_positive_array = findall(x->x>0.0,TMP_COL)
        idx_negative_array = findall(x->x<0.0,TMP_COL)

        # make new T -
        TNEW = T[idx_zero_array,:]

        # zero out non-zero elements -
        TMP_ARR = Array{Array{Float64,1},1}()
        for idx_positive ∈ idx_positive_array
            for idx_negative ∈ idx_negative_array

                # setup α and β -
			    α = abs(T[idx_positive,local_col_index]);
			    β = abs(T[idx_negative,local_col_index]);

                # compute a tmp row -
                TMP_ROW = α*T[idx_negative,:] .+ β*T[idx_positive,:];

                # grab the tmp row -
                push!(TMP_ARR,TMP_ROW)
            end
        end

		# ok, so we have a new set of possible pathways, check for independent pathways -
        P = transpose(hcat(TMP_ARR...))
		if (isempty(P) == false)
			# Update T -
			T = vcat(TNEW,P);
		end
    end # main -

	# ok, check for indepedent pathways -
	rows_to_keep_set = compute_convex_independent_rows(T)

	# grab: which rows should we keep?
    T_IND = T[(rows_to_keep_set |> collect),:]

    # return -
    return T_IND

end # function -

function compute_convex_independent_rows(P)::Set{Int64}

	# ok, so we have a new set of possible pathways, check for independent pathways -
	(NRP,_) = size(P)
	IPA = Matrix{Float64}(I,NRP,NRP)
	for i ∈ 1:NRP

		# get rᵢ -
		rᵢ = P[i,:]

		# compare rᵢ against rⱼ i neq j
		for j ∈ 1:NRP
			if (i != j)

				# get rⱼ -
				rⱼ = P[j,:]

				# check: is convex independent?
				if (is_convex_independent(rᵢ,rⱼ) == true)
					IPA[i,j] = 1.0
				else
					IPA[i,j] = 0.0
				end
			end
		end # inner
	end # outer

	# rows we should keep -
	rows_to_keep_set = Set{Int64}()
	for i ∈ 1:NRP
		sᵢ = sum(IPA[i,:])
		if (sᵢ == NRP)
			push!(rows_to_keep_set,i)
		end
	end

	# return -
	return rows_to_keep_set
end

function is_convex_independent(rᵢ,rⱼ)::Bool
	idx_zero_rᵢ = findall(x->x==0.0,rᵢ)
	idx_zero_rⱼ = findall(x->x==0.0,rⱼ)
	return issubset(idx_zero_rᵢ,idx_zero_rⱼ) ? false : true
end


text:
# gravity
get_gravity(g::T) where T <: Real = SVector{3,T}([0.0; 0.0; g])
get_gravity(g::Vector{T}) where T = SVector{3,T}(g)
get_gravity(g::SVector) = g


text:
doc"""
    FDist(ν1,ν2)

The *F distribution* has probability density function

$f(x; \nu_1, \nu_2) = \frac{1}{x B(\nu_1/2, \nu_2/2)}
\sqrt{\frac{(\nu_1 x)^{\nu_1} \cdot \nu_2^{\nu_2}}{(\nu_1 x + \nu_2)^{\nu_1 + \nu_2}}},
\quad x>0$

It is related to the [`Chisq`](:func:`Chisq`) distribution via the property that if $X_1
\sim \operatorname{Chisq}(\nu_1)$ and $X_2 \sim \operatorname{Chisq}(\nu_2)$, then
$(X_1/\nu_1) / (X_2 / \nu_2) \sim FDist(\nu_1, \nu_2)`.

```julia
FDist(d1, d2)     # F-Distribution with parameters d1 and d2

params(d)         # Get the parameters, i.e. (d1, d2)
```

External links

* [F distribution on Wikipedia](http://en.wikipedia.org/wiki/F-distribution)
    """
immutable FDist{T<:Real} <: ContinuousUnivariateDistribution
    ν1::T
    ν2::T

    function FDist(ν1::T, ν2::T)
        @check_args(FDist, ν1 > zero(ν1) && ν2 > zero(ν2))
        new(ν1, ν2)
    end
end

FDist{T<:Real}(ν1::T, ν2::T) = FDist{T}(ν1, ν2)
FDist(ν1::Integer, ν2::Integer) = FDist(Float64(ν1), Float64(ν2))
FDist(ν1::Real, ν2::Real) = FDist(promote(ν1, ν2)...)

@distr_support FDist 0.0 Inf

#### Conversions
function convert{T <: Real, S <: Real}(::Type{FDist{T}}, ν1::S, ν2::S)
    FDist(T(ν1), T(ν2))
end
function convert{T <: Real, S <: Real}(::Type{FDist{T}}, d::FDist{S})
    FDist(T(d.ν1), T(d.ν2))
end

#### Parameters

params(d::FDist) = (d.ν1, d.ν2)
@inline partype{T<:Real}(d::FDist{T}) = T

#### Statistics

mean{T<:Real}(d::FDist{T}) = (ν2 = d.ν2; ν2 > 2 ? ν2 / (ν2 - 2) : T(NaN))

function mode{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    ν1 > 2 ? ((ν1 - 2)/ν1) * (ν2 / (ν2 + 2)) : zero(T)
end

function var{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    ν2 > 4 ? 2ν2^2 * (ν1 + ν2 - 2) / (ν1 * (ν2 - 2)^2 * (ν2 - 4)) : T(NaN)
end

function skewness{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    if ν2 > 6
        return (2ν1 + ν2 - 2) * sqrt(8(ν2 - 4)) / ((ν2 - 6) * sqrt(ν1 * (ν1 + ν2 - 2)))
    else
        return T(NaN)
    end
end

function kurtosis{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    if ν2 > 8
        a = ν1 * (5ν2 - 22) * (ν1 + ν2 - 2) + (ν2 - 4) * (ν2 - 2)^2
        b = ν1 * (ν2 - 6) * (ν2 - 8) * (ν2 - 2)
        return 12a / b
    else
        return T(NaN)
    end
end

function entropy(d::FDist)
    (ν1, ν2) = params(d)
    hν1 = ν1/2
    hν2 = ν2/2
    hs = (ν1 + ν2)/2
    return log(ν2 / ν1) + lgamma(hν1) + lgamma(hν2) - lgamma(hs) +
        (1 - hν1) * digamma(hν1) + (-1 - hν2) * digamma(hν2) +
        hs * digamma(hs)
end

#### Evaluation & Sampling

@_delegate_statsfuns FDist fdist ν1 ν2

rand(d::FDist) = StatsFuns.RFunctions.fdistrand(d.ν1, d.ν2)


text:
"""
    solve_hadc(problem::HADC, gradient::ScalarGradient)

Compute the ADC using a homogenized ADC model (HADC). This is currently only implemented for
scalar gradients.
"""
function solve(problem::HADC, gradient::ScalarGradient)
    (; model, matrices, odesolver, reltol, abstol) = problem
    (; mesh, D) = model
    (; M_cmpts, S_cmpts, G, volumes) = matrices

    f = gradient.profile
    dir = gradient.dir
    TE = echotime(f)

    # Deduce sizes
    ncompartment = length(D)

    # Initial conditions
    ω₀ = zeros.(size.(mesh.points, 2))

    # Time dependent ODE function
    function Mdω!(dω, ω, p, t)
        # @show t
        (; mS, f, surfint) = p
        mul!(dω, mS, ω)
        dω .+= integral(f, t) .* surfint
    end

    # Create ODE function and Jacobian from matrices
    Jac!(J, _, p, t) = (J .= p.mS)

    # Allocate output array
    adc = zeros(ncompartment)

    # Iterate over compartments and gradient sequences and directions
    for icmpt = 1:ncompartment
        # Free diffusivity in gradient direction
        D₀ = dir' * D[icmpt] * dir

        # Surface integrals in gradient direction
        surfint = G[icmpt] * (D[icmpt] * dir)

        p = (; mS = -S_cmpts[icmpt], f, surfint)

        odefunction = ODEFunction(
            Mdω!,
            jac = Jac!,
            jac_prototype = p.mS,
            mass_matrix = M_cmpts[icmpt],
        )
        odeproblem = ODEProblem(odefunction, ω₀[icmpt], (0, TE), p, progress = false)

        # Solve ODE, keep all time steps (for integral)
        sol = OrdinaryDiffEq.solve(odeproblem, odesolver, reltol = reltol, abstol = abstol)

        # Integral over compartment boundary
        a, = quadgk(t -> integral(f, t) * (surfint' * sol(t)), 0, TE)

        # HADC (free diffusivity minus correction)
        adc[icmpt] = D₀ - a / volumes[icmpt] / int_F²(f)
    end

    adc
end


text:
@doc doc"""
    GroupOperationAction(group::AbstractGroupManifold, AD::ActionDirection = LeftAction())

Action of a group upon itself via left or right translation.
"""
struct GroupOperationAction{G,AD} <: AbstractGroupAction{AD}
    group::G
end

function GroupOperationAction(
    G::AbstractGroupManifold,
    ::TAD = LeftAction(),
) where {TAD<:ActionDirection}
    return GroupOperationAction{typeof(G),TAD}(G)
end

function show(io::IO, A::GroupOperationAction)
    print(io, "GroupOperationAction($(A.group), $(direction(A)))")
end

base_group(A::GroupOperationAction) = A.group

g_manifold(A::GroupOperationAction) = A.group

function switch_direction(A::GroupOperationAction)
    return GroupOperationAction(A.group, switch_direction(direction(A)))
end

apply(A::GroupOperationAction, a, x) = translate(A.group, a, x, direction(A))

apply!(A::GroupOperationAction, y, a, x) = translate!(A.group, y, a, x, direction(A))

function inverse_apply(A::GroupOperationAction, a, x)
    return inverse_translate(A.group, a, x, direction(A))
end

function inverse_apply!(A::GroupOperationAction, y, a, x)
    return inverse_translate!(A.group, y, a, x, direction(A))
end

function apply_diff(A::GroupOperationAction, a, x, v)
    return translate_diff(A.group, a, x, v, direction(A))
end

function apply_diff!(A::GroupOperationAction, vout, a, x, v)
    return translate_diff!(A.group, vout, a, x, v, direction(A))
end

function inverse_apply_diff(A::GroupOperationAction, a, x, v)
    return inverse_translate_diff(A.group, a, x, v, direction(A))
end

function inverse_apply_diff!(A::GroupOperationAction, vout, a, x, v)
    return inverse_translate_diff!(A.group, vout, a, x, v, direction(A))
end

function optimal_alignment(A::GroupOperationAction, x1, x2)
    return inverse_apply(switch_direction(A), x1, x2)
end

function optimal_alignment!(A::GroupOperationAction, y, x1, x2)
    return inverse_apply!(switch_direction(A), y, x1, x2)
end

function center_of_orbit(
    A::GroupOperationAction,
    pts::AbstractVector,
    q,
    mean_method::AbstractEstimationMethod,
)
    m = mean(A.group, pts, mean_method)
    return inverse_apply(switch_direction(A), q, m)
end
function center_of_orbit(A::GroupOperationAction, pts::AbstractVector, q)
    m = mean(A.group, pts)
    return inverse_apply(switch_direction(A), q, m)
end


text:
#
# Advent of Code
# Day 4: Security Through Obscurity
#

println("Advent of Code")
println("Day 4: Security Through Obscurity")
println("")

function pairSorter(pair::Pair{Char,Int})
  fullString = string(pair[2]) * string(pair[1])
  return fullString
end

function parseRooms(inputDataLines)
  total = 0

  for line in inputDataLines
    letterDict = Dict{Char, Int}()

    inSector = false
    sectorID = ""

    inChecksum = false
    checksum = ""

    # evaluate all the letters one by one
    for letter in line
      # check if we have [A-za-z]
      if isalpha(letter)
        # if we've started the checksum, add the letter to the checksum
        if inChecksum
          checksum = checksum * string(letter)
        # if not, and we haven't seen this letter before, start the count
        elseif !haskey(letterDict, letter)
          letterDict[letter] = 1
        # otherwise, increment the count
        else
          letterDict[letter] += 1
        end
      # but if we have a number, we must be in the sectorID
      elseif isnumber(letter)
        # check if we've been in the sectorID yet
        if !inSector
          # we're just starting building the sector string
          sectorID = string(letter)
          inSector = true
        else
          # or, add to the sector string
          sectorID = sectorID * string(letter)
        end
      # If it's not a letter or number, maybe we're starting the checksum
      elseif letter == '['
        if !inChecksum
          inChecksum = true
        end
      end
    end

    # generate the calculated checksum
    generatedChecksum = ""
    letterList = Vector{Pair{Char,Int}}()

    # we need to reverse the numbers for the sort
    lineLength = length(line)

    # get the key, value pairs into a list
    for kv in letterDict
      kv = Pair(kv[1], lineLength - kv[2])
      push!(letterList, kv)
    end

    # sort the list by the values, then the key
    sort!(letterList,by=pairSorter)

    # grab the first five (if they exist)
    if length(letterList) >= 5
      for item in letterList[1:5]
        generatedChecksum = generatedChecksum * string(item[1])
      end

      # compare the generatedChecksum to the checksum
      # and if it matches, then grab the value of the sectorID
      # and then add it to the total

      if generatedChecksum == checksum
        shift = parse(sectorID)

        total += shift

        # part 2 is to decrypt the room name with a shift cypher
        roomName = ""

        for letter in line
          if isnumber(letter)
            break
          end

          if letter == '-'
            newLetter = ' '
          else
            # all letters are lowercase, thank goodness
            newLetter = Char(rem((Int(letter) - 97 + shift), 26) + 97)
          end

          roomName *= string(newLetter)
        end

        println("Room: $roomName ($sectorID)")

      end
    end
  end

  return total

end

# Data file
# Input data has been saved to this file
inputDataFileName = "day04-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

#inputData = "aaaaa-bbb-z-y-x-123[abxyz]\n"
#inputData *= "a-b-c-d-e-f-g-h-987[abcde]\n"
#inputData *= "not-a-real-room-404[oarel]\n"
#inputData *= "totally-real-room-200[decoy]\n"

inputDataLines = split(strip(inputData), '\n')

part1Total = parseRooms(inputDataLines)

println("Part 1 total: $part1Total")


text:
module JuliaOrgs

using GitHub
using TOML
#using ProgressMeter

export get_orglist

ghauth = GitHub.authenticate(ENV["GITHUB_TOKEN"])

function get_orglist()
    orgnames = Dict{String, Any}()

    registrypath = ENV["HOME"] * "/.julia/registries/General"
    packages = TOML.parsefile(registrypath * "/Registry.toml")["packages"]

    for p in packages
        package_path = p[2]["path"]
        repo_url = TOML.parsefile(registrypath * "/" * package_path * "/Package.toml")["repo"]
        name = split(repo_url, "/")[4]
        orgnames[name] = nothing
    end

    function ghowner(x)
        try
            GitHub.owner(x, auth=myauth)
        catch
            x
        end
    end

    return asyncmap(x->ghowner(x), collect(keys(orgnames)))

end

end