text:
# TimeZones.jl supports parsing a direct time zone specification as provided by tzset(3).
# Unfortunately, different platforms only follow parts of the specification which results in
# differences from platform-to-platform.
#
# TimeZones.jl follows the tzset(3) specification and avoids supporting additional features
# which may not be present on all platforms. As the platform specific functionality only
# really deals with corner cases which are not required by any actual time zone definition
# supporting for these features is purely academic.
#
# Here are some examples of the platform specific differences. Under macOS time zone offset
# only supports hours under ±7 days (167 hours) and hours beyond this range result in
# falling back to UTC:
#
# ```bash
# $ TZ="FOO-0" date -r 0
# Thu  1 Jan 1970 00:00:00 UTC
#
# $ TZ="FOO-167" date -r 0
# Wed  7 Jan 1970 23:00:00 FOO
#
# $ TZ="FOO-168" date -r 0
# Thu  1 Jan 1970 00:00:00 UTC
# ```
#
# Under Linux the time zone offset only supports hours up to ±1 day (24 hours) and hours
# beyond that are clamped to exactly 24 hours.
#
# ```bash
# $ TZ="FOO" date --date @0
# Thu Jan  1 00:00:00 FOO 1970
#
# $ TZ="FOO-24" date --date @0
# Fri Jan  2 00:00:00 FOO 1970
#
# $ TZ="FOO-25" date --date @0
# Fri Jan  2 00:00:00 FOO 1970
# ```
#
# Under macOS the transition time only supports hours under ±7 days (167 hours) and hours
# beyond this range result in falling back to UTC:
#
# ```bash
# $ TZ="FOO-1BAR-2,1/0,1/167" date -r 1609455600
# Fri  1 Jan 2021 00:00:00 FOO
#
# $ TZ="FOO-1BAR-2,1/0,1/168" date -r 1609455600
# Thu 31 Dec 2020 23:00:00 UTC
# ```
#
# Under Linux the transition time appears to have no limit:
#
# ```bash
# $ TZ="FOO-1BAR-2,1/0,1/8761" date --date @1609455600
# Fri Jan  1 00:00:00 FOO 2021
# $ TZ="FOO-1BAR-2,1/0,1/8762" date --date @1609455600
# Fri Jan  1 01:00:00 BAR 2021
# ```

# Below is a listing of corner cases with TZ interpretation and behaviours on various
# platforms. We have tests cases for these assist in discovering any behaviour changes that
# may occur over time or when encountering a new platform.
#
# TZ=ABC
# - tzset(3): Invalid, an standard offset is required
# - Linux: Alternate way of writing `TZ=ABC+0`
# - macOS: Invalid, falls back to: `TZ=UTC`
#
# TZ="<>+1"
# - tzset(3): Invalid, an offset name must be at least three characters
# - Linux: Invalid, falls back to: `TZ="<>0"`
# - macOS: Supported
#
# TZ=A1
# - tzset(3): Invalid, an offset name must be at least three characters
# - Linux: Invalid, falls back to: `TZ="<>0"`
# - macOS: Supported
#
# TZ="<αβc>1"
# - tzset(3): Undefined behaviour, unicode support isn't mentioned
# - Linux: Invalid, falls back to: `TZ="<>0"`
# - macOS: Mostly supported, equivalent of `TZ="<____c>1"`
#
# TZ=ABC25
# - tzset(3): Invalid, offset hours must be between 0 and 24.
# - Linux: Clamped to ±24 hours and equivalent to `TZ=ABC24`
# - macOS: Supported
#
# TZ=ABC-168
# - tzset(3): Invalid, offset hours must be between 0 and 24.
# - Linux: Clamped to ±24 hours and equivalent to `TZ=ABC-24`
# - macOS: Invalid, only supports ±167 hours and falls back to: `TZ=UTC`
#
# TZ=FOO0BAR0,M1.5.6,1
# - tzset(3): Supported, DST transition occurs on the last Saturday in January. For
#   January's that only contain 4 Saturdays this is equivalent to M1.4.6.
# - Linux: Supported
# - macOS: Supported
#
# TZ=FOO0BAR0,365,1
# - tzset(3): Undefined behaviour, using a zero-based Julian day of 365 on a non-leap year
#   would be the first day of the following year.
# - Linux: Clamped non-leap year transition to first instant of the following year.
# - macOS: Supported
#
# TZ=FOO-1BAR-2,1/0,1/8762
# - tzset(3): Undefined behaviour, no specification is provided for transition time
# - Linux: Supported
# - macOS: Invalid, falls back to `TZ=UTC`
#
# TZ=FOO0BAR01/+1,2/-1
# - tzset(3): Undefined behaviour, no specification is provided for +/- time
# - Linux: Supported
# - macOS: Invalid, falls back to `TZ=UTC`

function tzdate(dt::DateTime, tz::AbstractString)
    iszero(millisecond(dt)) || throw(ArgumentError("Milliseconds not supported"))
    str = Dates.format(dt, dateformat"yyyy-mm-dd\THH:MM:SS")
    out_format = "%Y-%m-%dT%H:%M:%S%z (%Z)"
    withenv("TZ" => tz) do
        if Sys.isapple()
            readchomp(`date -j -f %Y-%m-%dT%H:%M:%S $str +$out_format`)
        elseif Sys.isunix()
            readchomp(`date --date=$str +$out_format`)
        else
            error("Unsupported platform")
        end
    end
end

@testset "System TZ behaviour" begin
    if Sys.isapple()
        @test tzdate(DateTime(1970), "ABC") == "1970-01-01T00:00:00+0000 (UTC)"
        @test tzdate(DateTime(1970), "<>+1") == "1970-01-01T00:00:00-0100 ()"
        @test tzdate(DateTime(1970), "A1") == "1970-01-01T00:00:00-0100 (A)"
        @test tzdate(DateTime(1970), "<αβc>1") == "1970-01-01T00:00:00-0100 (____c)"

        @test tzdate(DateTime(1970), "ABC24") == "1970-01-01T00:00:00-2400 (ABC)"
        @test tzdate(DateTime(1970), "ABC25") == "1970-01-01T00:00:00-2500 (ABC)"

        @test tzdate(DateTime(1970), "ABC-167") == "1970-01-01T00:00:00+16700 (ABC)"
        @test tzdate(DateTime(1970), "ABC-168") == "1970-01-01T00:00:00+0000 (UTC)"

        tz = "FOO0BAR0,M1.5.6,1"
        @test tzdate(DateTime(2020, 1, 25, 1, 59, 59), tz) == "2020-01-25T01:59:59+0000 (FOO)"
        @test tzdate(DateTime(2020, 1, 25, 2,  0,  0), tz) == "2020-01-25T02:00:00+0000 (BAR)"

        tz = "FOO0BAR0,365,1"
        @test tzdate(DateTime(2020,  1,  2, 1, 59, 59), tz) == "2020-01-02T01:59:59+0000 (BAR)"
        @test tzdate(DateTime(2020,  1,  2, 2,  0,  0), tz) == "2020-01-02T02:00:00+0000 (FOO)"
        @test tzdate(DateTime(2020, 12, 31, 1, 59, 59), tz) == "2020-12-31T01:59:59+0000 (FOO)"
        @test tzdate(DateTime(2020, 12, 31, 2,  0,  0), tz) == "2020-12-31T02:00:00+0000 (BAR)"

        @test tzdate(DateTime(2021, 1, 2, 1, 59, 59), tz) == "2021-01-02T01:59:59+0000 (BAR)"
        @test tzdate(DateTime(2021, 1, 2, 2,  0,  0), tz) == "2021-01-02T02:00:00+0000 (FOO)"
        @test tzdate(DateTime(2022, 1, 1, 1, 59, 59), tz) == "2022-01-01T01:59:59+0000 (FOO)"
        @test tzdate(DateTime(2022, 1, 1, 2,  0,  0), tz) == "2022-01-01T02:00:00+0000 (BAR)"

        @test tzdate(DateTime(2021), "FOO-1BAR-2,1/0,1/8761") == "2021-01-01T00:00:00+0000 (UTC)"
        @test tzdate(DateTime(2021), "FOO-1BAR-2,1/0,1/8762") == "2021-01-01T00:00:00+0000 (UTC)"

        tz = "FOO0BAR0,1/+1,2/-1"
        @test tzdate(DateTime(2020, 1, 2,  0, 59, 59), tz) == "2020-01-02T00:59:59+0000 (UTC)"
        @test tzdate(DateTime(2020, 1, 2,  1,  0,  0), tz) == "2020-01-02T01:00:00+0000 (UTC)"
        @test tzdate(DateTime(2020, 1, 2, 22, 59, 59), tz) == "2020-01-02T22:59:59+0000 (UTC)"
        @test tzdate(DateTime(2020, 1, 2, 23,  0,  0), tz) == "2020-01-02T23:00:00+0000 (UTC)"

    elseif Sys.isunix()
        @test tzdate(DateTime(1970), "ABC") == "1970-01-01T00:00:00+0000 (ABC)"
        @test tzdate(DateTime(1970), "<>+1") == "1970-01-01T00:00:00+0000 ()"
        @test tzdate(DateTime(1970), "A1") == "1970-01-01T00:00:00+0000 ()"
        @test tzdate(DateTime(1970), "<αβc>1") == "1970-01-01T00:00:00+0000 ()"

        @test tzdate(DateTime(1970), "ABC24") == "1970-01-01T00:00:00-2400 (ABC)"
        @test tzdate(DateTime(1970), "ABC25") == "1970-01-01T00:00:00-2400 (ABC)"

        @test tzdate(DateTime(1970), "ABC-167") == "1970-01-01T00:00:00+2400 (ABC)"
        @test tzdate(DateTime(1970), "ABC-168") == "1970-01-01T00:00:00+2400 (ABC)"

        tz = "FOO0BAR0,M1.5.6,1"
        @test tzdate(DateTime(2020, 1, 25), tz) == "2020-01-25T00:00:00+0000 (FOO)"
        @test tzdate(DateTime(2020, 1, 26), tz) == "2020-01-26T00:00:00+0000 (BAR)"

        tz = "FOO0BAR0,365,1"
        @test tzdate(DateTime(2020,  1,  2, 1, 59, 59), tz) == "2020-01-02T01:59:59+0000 (BAR)"
        @test tzdate(DateTime(2020,  1,  2, 2,  0,  0), tz) == "2020-01-02T02:00:00+0000 (FOO)"
        @test tzdate(DateTime(2020, 12, 31, 1, 59, 59), tz) == "2020-12-31T01:59:59+0000 (FOO)"
        @test tzdate(DateTime(2020, 12, 31, 2,  0,  0), tz) == "2020-12-31T02:00:00+0000 (BAR)"

        @test tzdate(DateTime(2021,  1,  2,  1, 59, 59), tz) == "2021-01-02T01:59:59+0000 (BAR)"
        @test tzdate(DateTime(2021,  1,  2,  2,  0,  0), tz) == "2021-01-02T02:00:00+0000 (FOO)"
        @test tzdate(DateTime(2021, 12, 31, 23, 59, 59), tz) == "2021-12-31T23:59:59+0000 (FOO)"
        @test tzdate(DateTime(2022,  1,  1,  0,  0,  0), tz) == "2022-01-01T00:00:00+0000 (BAR)"

        @test tzdate(DateTime(2021), "FOO-1BAR-2,1/0,1/8761") == "2021-01-01T00:00:00+0100 (FOO)"
        @test tzdate(DateTime(2021), "FOO-1BAR-2,1/0,1/8762") == "2021-01-01T00:00:00+0200 (BAR)"

        tz = "FOO0BAR0,1/+1,2/-1"
        @test tzdate(DateTime(2020, 1, 2,  0, 59, 59), tz) == "2020-01-02T00:59:59+0000 (FOO)"
        @test tzdate(DateTime(2020, 1, 2,  1,  0,  0), tz) == "2020-01-02T01:00:00+0000 (BAR)"
        @test tzdate(DateTime(2020, 1, 2, 22, 59, 59), tz) == "2020-01-02T22:59:59+0000 (BAR)"
        @test tzdate(DateTime(2020, 1, 2, 23,  0,  0), tz) == "2020-01-02T23:00:00+0000 (FOO)"
    end
end


text:
module UnderwaterAcoustics

using Requires
using DocStringExtensions

using DSP: amp2db, db2amp, pow2db, db2pow
using Printf

export °

const ° = π/180.0

# basic underwater acoustics
include("uw_basic.jl")

# propagation modeling
include("pm_core.jl")
include("pm_basic.jl")
include("pm_pekeris.jl")
include("pm_rays.jl")
include("pm_bellhop.jl")
include("pm_all.jl")

# plot recipes
function __init__()
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    include("plot.jl")
  end
end

end # module

ext: jl


text:
#---------------------------------------------------------------------#
#This routine advances the solution in time using SSP Osher-Shu RK2 or RK3
#Written by F.X. Giraldo / P.R. Mugg on 9/20/19
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function ssp_rk(q0,c,kstages,Δt,time)

   qp = q0
   q1 = q0
   for ik = 1:kstages
      if (kstages == 2)     #RK2
         if (ik == 1)
            α0 = 1.0; α1 = 0.0; β = 1.0
            t = time
         elseif (ik == 2)
            α0 = 0.5; α1 = 0.5; β = 0.5
            t = time + Δt
         end # ik
      elseif (kstages == 3)  #RK3
         if (ik == 1)
            α0 = 1.0; α1 = 0.0; β = 1.0
            t = time
         elseif (ik == 2)
            α0 = 3.0/4.0; α1 = 1.0/4.0; β = 1.0/4.0
            t = time + Δt
         elseif (ik == 3)
            α0 = 1.0/3.0; α1 = 2.0/3.0; β = 2.0/3.0
            t = time + (Δt/2)
         end # ik
      end #kstages

      #Construct RHS
      (rhs,f,g) = rhs_function(qp,c,t)

      #Update
      qp = α0*q0 + α1*q1 + Δt*β*rhs
      q1 = qp
   end #for ik

   return qp
end #end function


text:
module ROC

	using ROCAnalysis, MLLabelUtils
	using StatsBase: sample

	export AbstractOP, ComplexOP, SimpleOP, findop, changeop!, simpleop,
		AbstractPerfMetric, TPr, FPr, TNr, FNr

	abstract type AbstractOP end

	"""
	Operating point object with additional information.
	"""
	mutable struct ComplexOP <: AbstractOP
		targetclass::Int		# index of the target class
		weights::Matrix{Float64}	# weight matrix (class weights present only on the diagonal)
		rocdata::Matrix{Float64}	# Three column data matrix with the columns: θ i.e. threshold, FPr and FNr
		pos::Int			# position of chosen treshold in r.θ

		function ComplexOP(targetclass::Int, weights::Matrix{Float64}, rocdata::Matrix{Float64}, pos::Int)
			@assert isdiag(weights) "Weight matrix for operating point objects has to be diagonal."
			@assert pos>0 && pos<=size(rocdata,1) "Position of threhsold does not match length of the threshold vector"
			@assert targetclass>0 && targetclass<=size(weights,2) "The target class has to be an integer between 0 and $(size(weights,2))"
			return new(targetclass, weights, rocdata, pos)
		end
	end

	ComplexOP(targetclass::Int, weights::T where T<:AbstractVector{<:Real}, rocdata::Matrix{Float64}, pos::Int) =
		ComplexOP(targetclass, Matrix{Float64}(diagm(weights)), rocdata, pos)

	Base.show(io::IO, op::ComplexOP) = print(io, "Complex operating point object, pos=$(op.pos)/$(size(op.rocdata,1)),"*
					  		" targetclass=$(op.targetclass), weights=$(diag(op.weights))")

	"""
	Operating point object with no additional information.
	"""
	mutable struct SimpleOP <: AbstractOP
		weights::Matrix{Float64}	# weight matrix (class weights present only on the diagonal)
	end

	SimpleOP(w::T where T<:AbstractVector{<:Real}) = SimpleOP(Matrix{Float64}(diagm(w)))
	SimpleOP(cop::ComplexOP) = SimpleOP(cop.weights)

	Base.show(io::IO, op::SimpleOP) = print(io, "Simple Operating point object, weights=$(diag(op.weights))")

	# Define performance metric types
	abstract type AbstractPerfMetric end
	struct TPr <: AbstractPerfMetric end 	# true positive rate (sensitivity)
	struct TNr <: AbstractPerfMetric end 	# true negative rate (specificity)
	struct FPr <: AbstractPerfMetric end 	# false positive rate (fall-out)
	struct FNr <: AbstractPerfMetric end 	# false negative rate

	Base.show(io::IO, apm::TPr) = print(io, "TPr")
	Base.show(io::IO, apm::TNr) = print(io, "TNr")
	Base.show(io::IO, apm::FPr) = print(io, "FPr")
	Base.show(io::IO, apm::FNr) = print(io, "FNr")

	# Define threshold search functions
	_θ_search_(rocdata::AbstractMatrix{Float64}, pm::TPr, threshold::Float64) = begin
		pos = findlast(view(rocdata,:,3).<=1-threshold) # 1-rocdata[:,3].>=threshold
		pos = ifelse(pos==0, 1, pos)
		return rocdata[pos,1], pos
	end

	_θ_search_(rocdata::AbstractMatrix{Float64}, pm::FNr, threshold::Float64) = begin
		pos = findlast(view(rocdata,:,3).<=threshold)
		return rocdata[pos,1], pos
	end

	_θ_search_(rocdata::AbstractMatrix{Float64}, pm::FPr, threshold::Float64) = begin
		pos = findlast(view(rocdata,:,2).>=threshold)
		pos = ifelse(pos==0, 1, pos)
		return rocdata[pos,1], pos
	end

	_θ_search_(rocdata::AbstractMatrix{Float64}, pm::TNr, threshold::Float64) = begin
		pos = findlast(rocdata[:,2].>=1-threshold) # 1-rocdata[:,2].<=threshold
		return rocdata[pos,1], pos
	end

	# Define functions that generate the 'rocdata' i.e. threhsholds, false positive and false negative rates
	_get_rocdata_(x::AbstractMatrix{T}, y::AbstractVector{S}, yu::AbstractVector{S}, targetclass::S, maxops::Int=100) where {T,S} = begin
		idx = findin(yu, [targetclass])[1] # find index corresponding to the targetclass
		n = size(x,2)

		# Define maximum number of operating points
		# and sample thresholds if necessary
		positives = view(x, idx, y.==targetclass)
		upositives::Vector{T} = sort(unique(positives))
		nops::Int = length(upositives)+2
		if nops > maxops
			nops = maxops
			upositives = sort(sample(upositives, maxops-2, replace=false))
		end

		# Pre-allocate output and calculate false positives, false negatives for the data
		rocdata = zeros(Float64, nops, 3)
		rocdata[1,:] = [0.0, 1.0, 0.0]						# minimum threshold case
		rocdata[nops,:] = [1.0, 0.0, 1.0]					# maximum threshold case
		nP = length(positives)
		nN = n-nP

		fp::Float64 = 0.0
		fn::Float64 = 0.0
		@inbounds for (i,θ) in enumerate(upositives)
			fp = 0.0
			fn = 0.0
			# Count false positives and false negatives for the current threshold
			@simd for j in 1:n
				fp = fp + ifelse( (y[j]!=targetclass) && (x[idx,j]*(1.0-θ)/θ >= 1.0-x[idx,j]), 1.0, 0.0)
				fn = fn + ifelse( (y[j]==targetclass) && (x[idx,j]*(1.0-θ)/θ < 1.0-x[idx,j]), 1.0, 0.0)
			end
			rocdata[i+1,1] = upositives[i]
			rocdata[i+1,2] = fp/nP
			rocdata[i+1,3] = fn/nN
		end
		return idx, rocdata
	end

	_get_rocdata_ra_(x::AbstractMatrix, y::AbstractVector{S}, yu::AbstractVector{S}, targetclass::S, maxops::Int=100) where S = begin
		# Preprocess the matrix of probabilities (transform from probability to log-likelihood ration scores)
		idx = findin(yu, [targetclass])[1] 					# find index corresponding to the targetclass
		notidx = setdiff(1:length(yu), idx)     				# find indexes corresponding to the other classes
		positives = 2 * x[idx, y.==targetclass] - 1 				# process probabilities for ROCAnalysis.jl
		negatives = 2* x[idx, y.!=targetclass] - 1

		# Call ROCAnalysis
		r = ROCAnalysis.roc(positives, negatives)

		# If more ops than needed, limit their number
		nops = length(r.θ)+1
		opidx = collect(1:min(nops,maxops))
		if nops > maxops
			opidx[1] = 1; opidx[end] = nops
			opidx[2:end-1] = sort(sample(2:nops-1, maxops-2, replace=false))
		end
		rocdata = [[(1.0+r.θ)/2.0;1.0][opidx] r.pfa[opidx] r.pmiss[opidx]]
		return idx, rocdata
	end

	# Define op generation functions
	findop(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractVector{N}, targetclass::N,
			pm::AbstractPerfMetric, threshold::Float64, method::Symbol=:j4pr, maxops::Int=100) where N =
		findop(x[1], x[2], targetclass, pm, threshold, method, maxops)

	findop(x::T where T<:AbstractMatrix, y::S  where S<:AbstractVector{N}, targetclass::N,
			pm::AbstractPerfMetric, threshold::Float64, method::Symbol=:j4pr, maxops::Int=100) where N = begin

		# Checks
		yu = sort(unique(y))
		@assert size(x,1) == length(yu) "The number of unique labels does not match dimensionality of class-wise probabilities"
		@assert size(x,2) == length(y) "The number of labels does not match the number of observations"
		@assert targetclass in yu "Target class does not seem to be present in list of classes"

		# Apply appropriate ROC analysis method
		if method == :j4pr
			idx, rocdata = _get_rocdata_(x, y, yu, targetclass, maxops)
		elseif method == :ra
			idx, rocdata = _get_rocdata_ra_(x, y, yu, targetclass, maxops)
		else
			warn("[findop] Unknown method for op selection, using :j4pr.")
			idx, rocdata = _get_rocdata_(x, y, yu, targetclass, maxops)
		end

		# Generate op data and create object
		op = ComplexOP(idx, eye(length(yu)), rocdata, 1)

		# Search for best operating point
		changeop!(op, pm, clamp(threshold, 0.0, 1.0))

		return op
	end

	# Define op changing function
	changeop!(op::ComplexOP, pm::AbstractPerfMetric, threshold::Float64) = begin

		# Call threshold search
		θ::Float64, pos::Int = _θ_search_(op.rocdata, pm, threshold)

		# Update position in op
		op.pos = pos

		# Construct the weights
		op.weights[diagind(op.weights)] = 1.0; 					# update diagonal elements
		op.weights[op.targetclass, op.targetclass] = (1-θ+eps())/(θ+eps()) 	# update target class weights
		nothing;
	end

	changeop!(op::ComplexOP, pos::Int64=1) = begin

		# Check position
		@assert pos >= 1 && pos <= size(op.rocdata,1)
			"[changeop!] OP position has to be an integer between 1 and $(size(op.rocdata,1))."

		# Get threshold directly
		θ::Float64 = op.rocdata[pos,1]

		# Update position in op
		op.pos = pos

		# Construct the weights
		op.weights[diagind(op.weights)] = 1.0; 					# update diagonal elements
		op.weights[op.targetclass, op.targetclass] = (1-θ+eps())/(θ+eps()) 	# update target class weights
		nothing;
	end

	# Conversion function from a complex op to a simple op
	simpleop(op::ComplexOP) = SimpleOP(op)
	simpleop(weights::T where T<:AbstractVector{<:Real}) = SimpleOP(diagm(float.(weights)))
	simpleop(weights::T where T<:AbstractMatrix{<:Real}) = begin
		@assert isdiag(weights) && (size(weights,1)==size(weights,2)) "[simpleop] Weight matrix has to be square and diagonal"
		SimpleOP(float.(weights))
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	findop(class, pm, val, method, maxops)

Constructs an untrained cell that when piped data inside, returns the operating point
trying to optimize the performance metric `pm` around the value `val` for class `class`.

# Arguments
  * `class` can be either a specific class i.e. "setosa" or an index in the sorted class list
  * `pm::ROC.AbstractPerfMetric` is the performance metric. Can be `ROC.TPr()`, `ROC.TNr()`, `ROC.FPr()` or `ROC.FNr()`
  * `val::Float64` is the value desired for the metric
  * `method::Symbol` method through which to find the op; available `:j4pr` and `:ra` i.e. `ROCAnalysis.jl` (default `:j4pr`)
  * `maxops::Int` is the maximum number of operating points to consider (default `100`)

# Examples
```
julia> D = DataGenerator.normal2d(2000);        # generate dataset
       (tr,ts) = splitobs(shuffleobs(D),0.1);   # split in training and testing 10%/90%
       w = knn(5,smooth=:dist);                 # define untrained classifier
       wt = w(tr);                              # train classifier
       out = ts|>wt;                            # calculate probabilities on test data
       vop = findop(1,ROC.FPr(),0.025,:j4pr);   # define an operating point search for the first class
       top = vop(out)                           # search for op where for the first class in FPr=0.025
       pt=wt + top;                             # make pipeline
       # pass test data through pipe and calculate the confusion matrix
       +ts |> pt |> x->j4pr.confusionmatrix(targets(indmax,x)-1, Int.(-ts), normalize=true, showmatrix=true)
       # plot ROC
       rocplot(top)

reference labels (columns):
  "0"   "1"
------------------
0.72   0.025555555555555557
0.28   0.9744444444444444
------------------
       ┌────────────────────────────────────────┐
     1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡠⠴⠋⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉⠉│
       │⠀⠀⠀⠀⠀⠀⠀⠀⢀⡴⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⢠⠞⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⢀⡴⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⢀⣠⠞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⡎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⢸⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⡞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⢰⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⡞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   0.5 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
       └────────────────────────────────────────┘
       0                                        1
```
"""
findop(class, pm::ROC.AbstractPerfMetric, val::Float64, method::Symbol=:j4pr, maxops::Int=100) =
	FunctionCell(findop, (class, pm, clamp(val,0.0,1.0), method, maxops), ModelProperties(),
	      		"OP finder, on \"$class\", targeting $pm@$val, method=$method, maxops=$maxops")

############################
# DataCell/Array Interface #
############################

"""
	findop(x, class, pm, val, method, maxops)

Searches for the operating point trying to optimize the performance metric `pm` around the value
`val` for class `class`. The data `x` is assumed to be either a `Tuple` or a `DataCell` that contain
class probabilities and true labels.
"""
# Training
findop(x::T where T<:CellDataL, class, pm::ROC.AbstractPerfMetric, val::Float64, method::Symbol=:j4pr, maxops::Int=100) =
	findop((getx!(x), gety(x)), class, pm, val, method, maxops)
findop(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, class, pm::ROC.AbstractPerfMetric, val::Float64,
       method::Symbol=:j4pr, maxops::Int=100) = findop((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), class, pm, val, method, maxops)
findop(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, class, pm::ROC.AbstractPerfMetric, val::Float64,
       method::Symbol=:j4pr, maxops::Int=100) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[findop] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	origclass = (class isa Int) ? ind2label(class,enc) : class

	# Get op
	op = ROC.findop(x, origclass, pm, clamp(val,0.0,1.0), method, maxops)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), nvars(x[1]))

	FunctionCell(applyop, Model(op, modelprops), "OP, complex, on \"$origclass\"")
end

# Execution (operating point application)
applyop(x::T where T<:CellData, model::Model{<:ROC.AbstractOP}) =
	datacell(applyop(getx!(x), model), gety(x))

applyop(x::T where T<:AbstractVector, model::Model{<:ROC.AbstractOP}) =
	applyop(mat(x, LearnBase.ObsDim.Constant{2}()), model)

applyop(x::T where T<:AbstractMatrix, model::Model{<:ROC.AbstractOP}) =
	model.data.weights * x

##########################################
# Function interface for OP modification #
##########################################

# Change the op of a trained function cell
changeop!(op::T where T<:CellFunT{<:Model{<:ROC.ComplexOP}}, pm::ROC.AbstractPerfMetric, threshold::Float64) = begin
	ROC.changeop!(getx!(op).data, pm, threshold)
end

changeop!(op::T where T<:CellFunT{<:Model{<:ROC.ComplexOP}}, pos::Int) = begin
	ROC.changeop!(getx!(op).data, pos)
end

# Create a simple op
simpleop(weights::T where T<:AbstractArray{<:Real}) =
	FunctionCell(applyop, Model(ROC.simpleop(float.(weights)), ModelProperties(length(weights), length(weights))),
	      "OP, simple")

simpleop(op::T where T<:CellFunT{<:Model{<:ROC.ComplexOP}}) =
	FunctionCell(applyop, Model(ROC.simpleop(getx(op).data), getx(op).properties),
	      "OP, simple")


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

## 12 steps to Navier–Stokes
# ======
# ***

# This Jupyter notebook continues the presentation of the **12 steps to Navier–Stokes**,
# the practical module taught in the interactive CFD class of [Prof. Lorena Barba](http://lorenabarba.com).
# You should have completed [Step 1](./01_Step_1.jl) before continuing, having written your own Julia script
# or notebook and having experimented with varying the parameters of the discretization and observing what happens.

## Step 2: Nonlinear Convection
# -----
# ***

# Now we're going to implement nonlinear convection using the same methods as in step 1.  The 1D convection equation is:

# $$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0$$

# Instead of a constant factor $c$ multiplying the second term, now we have the solution $u$ multiplying it. Thus, the second term of the equation is now *nonlinear*. We're going to use the same discretization as in Step 1 — forward difference in time and backward difference in space. Here is the discretized equation.

# $$\frac{u_i^{n+1}-u_i^n}{\Delta t} + u_i^n \frac{u_i^n-u_{i-1}^n}{\Delta x} = 0$$

# Solving for the only unknown term, $u_i^{n+1}$, yields:

# $$u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n)$$

# As before, the Julia code starts by loading the necessary libraries. Then, we declare some variables
# that determine the discretization in space and time (you should experiment by changing these parameters to see what happens).
# Then, we create the initial condition $u_0$ by initializing the array for the
# solution using $u = 2\ @\ 0.5 \leq x \leq 1$  and $u = 1$ everywhere else in $(0,2)$ (i.e., a hat function).
# **Note**
# We're now going to be exploting Julia's amazing capabilities for Unicode. This lets us write code that is closer to the mathematics we want to model.
# To type `Δx`, write `\Delta<TAB>` and Julia should autocomplete it. You can also autocomplete incomplete characters so there's
# no need to memeorize them all, although they mirror how they're written in LaTeX so that helps a lot.

import Plots    # and our 2D plotting library
Plots.plotly()

nₓ = 41 # n\_<TAB>x
Δx = 2 / (nₓ - 1)
nₜ = 20    #nₜ is the number of timesteps we want to calculate
Δt = .025  #Δt is the amount of time each timestep covers (delta t)

u = ones(nₓ)      #as before, we initialize u with every value equal to 1.
#then set u = 2 between 0.5 and 1 as per our I.C.s
u[Int.((.5/Δx):(1/Δx)] .= 2 # This is a very numpy-ish way of writing this

# An equivalent way would be just
for i in length((.5/Δx):(1/Δx))
      u[i] = 2
end

uₙ = ones(nₓ) #initialize our placeholder array un, to hold the time-stepped solution

# The code snippet below is *unfinished*. We have copied over the line from [Step 1](./01_Step_1.jl)
# that executes the time-stepping update. Can you edit this code to execute the nonlinear convection instead?

for n in nₜ  #iterate through time
    uₙ = copy(u) ##copy the existing values of u into uₙ
    for i in 1:nₓ  ##now we'll iterate through the u array

     ###This is the line from Step 1, copied exactly.  Edit it for our new equation.
     ###then uncomment it and run the cell to evaluate Step 2
           ###u[i] = uₙ[i] - c * dt / dx * (uₙ[i] - uₙ[i-1])
    end
end

plot(range(start = 0, stop = 2, length = nₓ), u) ##Plot the results

# What do you observe about the evolution of the hat function under the nonlinear convection equation? What happens when you change the numerical parameters and run again?

## Learn More

# For a careful walk-through of the discretization of the convection equation with finite differences (and all steps from 1 to 4), watch **Video Lesson 4** by Prof. Barba on YouTube.

# TODO
# from IPython.display import YouTubeVideo
# YouTubeVideo('y2WaK7_iMRI')

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
using GridapGPU
using Documenter

makedocs(;
    modules=[GridapGPU],
    authors="Santiago Badia",
    repo="https://github.com/gridap/GridapGPU.jl/blob/{commit}{path}#L{line}",
    sitename="GridapGPU.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gridap.github.io/GridapGPU.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/gridap/GridapGPU.jl",
)


text:
struct Weight <: EdgeModifier
    w
end

struct Start <: EdgeModifier
    s
end

@testset "Declare Communative" begin
    StenoGraphs.@communative foo(a::Weight, b::Start) = a.w * b.s
    @test foo(Weight(2), Start(3)) == foo(Start(3), Weight(2))

    StenoGraphs.@communative function bar(x::Vector{m} where {m <: EdgeModifier}, y::Int)
        length(x)*y
    end
    @test bar([Weight(1), Start(2)], 2) == 4
end


text:
# ---
# title: Custom Title
# id: custom_id
# ---

# This is the content


text:
# test data generation

using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using IntervalTrees,SparseArrays,NearestNeighbors
using ViewerGL
GL = ViewerGL

function input_collection(data::Array)::Lar.LAR
	assembly = Lar.Struct(data)
	return Lar.struct2lar(assembly)
end

V,FV = Lar.sphere(2)([9,12])
#V,FV = Lar.sphere(2)([12,16])
EV = Lar.simplexFacets(FV)
#EV = Lar.quads2triangles(FV)
mysphere = V,FV,EV

data3d1 = mysphere
data3d2 = Lar.Struct([ Lar.t(0,1,0), mysphere ])
data3d3 = Lar.Struct([ Lar.t(0,0.5,0), Lar.s(0.4,0.4,0.4), mysphere ])
data3d4 = Lar.Struct([ Lar.t(4,0,0), Lar.s(0.8,0.8,0.8), mysphere ])
data3d5 = Lar.Struct([ Lar.t(4,0,0), Lar.s(0.4,0.4,0.4), mysphere ])

model3d = input_collection([ data3d5 ])#, data3d2, data3d3, data3d4, data3d5 ])
V,FV,EV = model3d
VV = [[k] for k in 1:size(V,2)];

GL.VIEW( GL.numbering(0.125)((V,[VV, EV, FV]),GL.COLORS[1],0.25) );

# copEV = Lar.coboundary_0(EV)
# copFE = Lar.coboundary_1(V,FV,EV)
# W = convert(Lar.Points, V')


text:
push!(LOAD_PATH,"../src/")
using Documenter
using ShapeFromShading

makedocs(
    sitename = "Documentation",
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md"
    ],
    modules = [ShapeFromShading]
)
deploydocs(repo = "github.com/betttris13/ShapeFromShading.jl.git")


text:
using Arrow, Chain, CSV, DataFrameMacros, DataFrames

axb_data = CSV.read("../../Patricia_Amillos/full_axb.csv", DataFrame; missingstrings="NA"); # view in table viewer

      # check for columns that are constant - drop them and Column1, which is row numbers
constants = names(axb_data, nm -> isone(length(unique(axb_data[:, nm]))))
axb_reduced = @chain axb_data begin
    select(Not(push!(constants, "Column1")))
    @transform!(
        :accuracy = :response == :answer_key, # evaluate accuracy as Boolean
        :item = string('I', lpad(:item, 2, '0')),
        :trained = :trained == 1,              # convert from number to Boolean
        :x_dental = :x_dental == 1,
    )
end
Arrow.write(
    "./data/axb_data.arrow",
    select(   # re-order the columns so the important ones are first
        axb_reduced,
        :ID,
        :item,
        :accuracy,
        :Task,
        :trained,
        :);
    compress = :zstd
    )


text:
import PyPlot.plt

# #############################################################################
#Spectrum datafile I/O

#Define spectrum rawdata type
immutable JLabSpectrum
    date::String
    real_time::Int64
    live_time::Int64
    bins::Array{Int64,1}
    counts::Array{Int64,1}
    rois::Array{Int64,2}
end

#Read ascii .spe file
function readspe(filename::String)
    #Read file
    file::IOStream = open(filename,"r")
    lines::Array{String,1} = readlines(file)

    #Find fields
    date_mea_idx::Int64 = findfirst(lines,"\$DATE_MEA:\r\n")
    meas_tim_idx::Int64 = findfirst(lines,"\$MEAS_TIM:\r\n")
    data_idx::Int64 = findfirst(lines,"\$DATA:\r\n")
    roi_idx::Int64 = findfirst(lines,"\$ROI:\r\n")

    #Extract data -------------------------------------------------------------
    #Date string
    date_mea::String = chomp(lines[date_mea_idx+1])
    #Measurement times (Live, Real)
    meas_tim::Array{Int64,1} = int64(split(lines[meas_tim_idx+1]," "))
    #Data
    binlim::Array{Int64,1} = int64(split(lines[data_idx+1]," "))
    bins::Array{Int64,1} = [binlim[1]:binlim[2]]
    counts::Array{Int64,1} = int64(lines[(data_idx+2):(data_idx+2+binlim[2]-binlim[1])])
    start = 1
    while counts[start] == 0; start += 1; end
    bins = bins[start:end]
    counts = counts[start:end]
    #ROI
    roinum::Int64 = int64(lines[roi_idx+1])
    rois::Array{Int64,2} = Array(Int64,roinum,2)
    for i = 1:roinum
        rois[i,:] = int64(split(lines[roi_idx+1+i]))
    end

    return JLabSpectrum(date_mea,meas_tim[1],meas_tim[2],bins,counts,rois)
end

#Show bins
function plotspec(spec::JLabSpectrum;logged::Bool=true)
    if logged
        y::Array{Float64,1} = zeros(Float64,length(spec.counts))
        nonzs::Array{Int64,1} = find(spec.counts)
        y[nonzs] = log(spec.counts[nonzs])
        ylabel::String = "log(Counts)"
    else
        y = float64(spec.counts)
        ylabel = "Counts"
    end

    for i=1:length(spec.counts)
        color = [0,0,0.75]
        for j=1:size(spec.rois,1)
            if (spec.rois[j,1] <= spec.bins[i] <= spec.rois[j,2])
                color = [0,0.75,0.75]
                break
            end
        end
        plt.plot([spec.bins[i],spec.bins[i]],[0,y[i]],color=color)
    end
    plt.xlabel("Bin Number")
    plt.ylabel(ylabel)
    plt.axis([spec.bins[1],spec.bins[end],0,1.1*maximum(y)])
    return nothing
end


text:
using CUDAdrv, CUDAnative

function hello_world()
    @cuprintf("Greetings from block %u, thread %u!\n", blockIdx().x, threadIdx().x)
    return
end

dev = CuDevice(0)
ctx = CuContext(dev)

@cuda (2,2) hello_world()

synchronize()

destroy(ctx)


text:
@testset "filter" begin
    i = HashIndices([1,2,3,4,5])

    @test issetequal(filter(iseven, i)::HashIndices, [2, 4])
    @test issetequal(filter(isodd, i)::HashIndices, [1, 3, 5])

    @test issetequal(filterview(iseven, i)::Dictionaries.FilteredIndices, [2, 4])
    @test issetequal(filterview(isodd, i)::Dictionaries.FilteredIndices, [1, 3, 5])

    d = HashDictionary([1,2,3,4,5], [1,3,2,4,5])

    @test issetequal(pairs(filter(iseven, d)::HashDictionary), [3=>2, 4=>4])
    @test issetequal(pairs(filter(isodd, d)::HashDictionary), [1=>1, 2=>3, 5=>5])

    @test issetequal(pairs(filterview(iseven, d)::Dictionaries.FilteredDictionary), [3=>2, 4=>4])
    @test issetequal(pairs(filterview(isodd, d)::Dictionaries.FilteredDictionary), [1=>1, 2=>3, 5=>5])
end


text:
using StaticUnivariatePolynomials
using Test
using ForwardDiff

import StaticUnivariatePolynomials: constant, derivative, integral, coefficient_gradient, exponential_integral

include("monomial_basis.jl")
include("bernstein_basis.jl")


text:
module MacroTools

using DataStructures, Compat
using Compat.Markdown
export @match, @capture, sourcewalk, sourcemap

include("match/match.jl")
include("match/types.jl")
include("match/union.jl")
include("match/macro.jl")

include("utils.jl")
include("structdef.jl")

include("examples/destruct.jl")
include("examples/threading.jl")
include("examples/forward.jl")

include("cstparser/CSTParser.jl")

module hack; using CSTParser; end

using .CSTParser

include("patch/diff.jl")
include("patch/cst.jl")

const animals = Symbol[]
const animals_file = joinpath(@__DIR__, "..", "animals.txt")

# Load and initialize animals symbols.
_animals = split(read(animals_file, String))
resize!(animals, length(_animals))
animals .= Symbol.(lowercase.(_animals))
Compat.Random.shuffle!(animals)

end # module


text:
using HetaSimulator, Plots

p = load_platform("./jak2stat5");

### test run

sim(p.models[:nameless], tspan = (0., 100.), alg = CVODE_BDF()) |> plot

### Load conditions

cond_df = read_conditions("./jak2stat5/Data/conditions.csv")
add_conditions!(p, cond_df)

# chapter 2.3.1, page 21, figure S11 => task1
dts1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_Long_log10-mod.csv");
add_measurements!(p, dts1_csv)

# chapter 2.3.2, page 23, figure S12
dts2_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_Concentrations_log10-mod.csv");
add_measurements!(p, dts2_csv)

# chapter 2.3.3, page 25, figure S13
dts3_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_RNA_log10-mod.csv");
add_measurements!(p, dts3_csv)

# chapter 2.3.4, page 28, figure S14
dts4_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_ActD0_log10-mod.csv");
dts4_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_ActD1_log10-mod.csv");
add_measurements!(p, dts4_0_csv)
add_measurements!(p, dts4_1_csv)

# chapter 2.3.5, page 31, figure S15
dts5_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_Fine_log10-mod.csv");
add_measurements!(p, dts5_csv)

# chapter 2.3.6, page 34, figure S16
dts6_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_CISoe0_log10-mod.csv");
dts6_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_CISoe1_log10-mod.csv");
add_measurements!(p, dts6_0_csv)
add_measurements!(p, dts6_1_csv)

# chapter 2.3.7, page 36, figure S17
dts7_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_CISoe_pEpoR0_log10-mod.csv");
dts7_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_CISoe_pEpoR1_log10-mod.csv");
add_measurements!(p, dts7_0_csv)
add_measurements!(p, dts7_1_csv)

# chapter 2.3.8, page 38, figure S18
dts8_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_SOCS3oe0_log10-mod.csv");
dts8_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_SOCS3oe1_log10-mod.csv");
add_measurements!(p, dts8_0_csv)
add_measurements!(p, dts8_1_csv)

# chapter 2.3.9, page 40, figure S19
dts9_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_SHP1oe0_log10-mod.csv");
dts9_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_SHP1oe1_log10-mod.csv");
add_measurements!(p, dts9_0_csv)
add_measurements!(p, dts9_1_csv)

# chapter 2.3.10, page 42, figure S20
dts10_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR72.5e-9_log10-mod.csv");
dts10_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR72.5e-8_log10-mod.csv");
dts10_2_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR72.5e-7_log10-mod.csv");
dts10_3_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR72.5e-6_log10-mod.csv");
dts10_4_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR72.5e-5_log10-mod.csv");
add_measurements!(p, dts10_0_csv)
add_measurements!(p, dts10_1_csv)
add_measurements!(p, dts10_2_csv)
add_measurements!(p, dts10_3_csv)
add_measurements!(p, dts10_4_csv)

# chapter 2.3.11, page 44, figure S21
dts11_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR302.5e-9_log10-mod.csv"); # 0
dts11_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR302.5e-8_log10-mod.csv"); # 1
dts11_2_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR301.25e-7_log10-mod.csv"); # 2
dts11_3_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR302.5e-7_log10-mod.csv"); # 3
dts11_4_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR301.25e-6_log10-mod.csv"); # 4
dts11_5_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR302.5e-6_log10-mod.csv"); # 5
add_measurements!(p, dts11_0_csv)
add_measurements!(p, dts11_1_csv)
add_measurements!(p, dts11_2_csv)
add_measurements!(p, dts11_3_csv)
add_measurements!(p, dts11_4_csv)
add_measurements!(p, dts11_5_csv)

# chapter 2.3.12, page 46, figure S22
dts12_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR101.25e-8_log10-mod.csv"); # 1.25e-8 # 0
dts12_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR101.75e-8_log10-mod.csv"); # 1.75e-8 # 1
dts12_2_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR101.25e-7_log10-mod.csv"); # 1.25e-7 # 2
dts12_3_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR101.7675e-7_log10-mod.csv"); # 1.7675e-7 # 3 # maybe mistake
dts12_4_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR102.5e-7_log10-mod.csv"); # 2.5e-7 # 4
dts12_5_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR102.5e-6_log10-mod.csv"); # 2.5e-6 # 5
add_measurements!(p, dts12_0_csv)
add_measurements!(p, dts12_1_csv)
add_measurements!(p, dts12_2_csv)
add_measurements!(p, dts12_3_csv)
add_measurements!(p, dts12_4_csv)
add_measurements!(p, dts12_5_csv)

# chapter 2.3.13, page 48, figure S23
dts13_0_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR902.5e-9_log10-mod.csv"); # 0
dts13_1_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR902.5e-8_log10-mod.csv"); # 1
dts13_2_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR901.25e-7_log10-mod.csv"); # 2
dts13_3_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR902.5e-7_log10-mod.csv"); # 3
dts13_4_csv = read_measurements("./jak2stat5/Data/JAK2STAT5_CFUE_DR902.5e-6_log10-mod.csv"); # 4
add_measurements!(p, dts13_0_csv)
add_measurements!(p, dts13_1_csv)
add_measurements!(p, dts13_2_csv)
add_measurements!(p, dts13_3_csv)
add_measurements!(p, dts13_4_csv)

res = sim(p; alg = CVODE_BDF())
plot(res) # very bad picture

max_stars_repo_name: hetalang/qsp-benchmark

max_issues_repo_name: hetalang/qsp-benchmark

max_forks_repo_name: hetalang/qsp-benchmark


text:
import Knet.Ops21: pool
using Statistics: mean

using CUDA.CUDNN:
    cudnnPoolingForward,
    cudnnPoolingForward!,
    cudnnPoolingBackward,
    cudnnGetPoolingNdForwardOutputDim,
    cudnnPoolingDescriptor,
        cudnnPoolingDescriptor_t,
        cudnnCreatePoolingDescriptor,
        cudnnSetPoolingNdDescriptor,
        cudnnDestroyPoolingDescriptor,
    cudnnPoolingMode_t,
        CUDNN_POOLING_MAX,                           # 0,
        CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING, # 1, /* count for average includes padded values */
        CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING, # 2, /* count for average does not include padded values */
        CUDNN_POOLING_MAX_DETERMINISTIC,             # 3
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    pooldims

function pool(x::GPUVal; op=maximum, window=2, padding=0, stride=window, propagateNaN=false, includePadding=false, channelmajor=false, alpha=1)
    mode = (op==maximum ? CUDNN_POOLING_MAX :
            op==mean ? (includePadding ?
                        CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING :
                        CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING) :
            error("Unknown pooling op=$op"))
    nanOpt = (propagateNaN ? CUDNN_PROPAGATE_NAN : CUDNN_NOT_PROPAGATE_NAN)
    format = (channelmajor ? CUDNN_TENSOR_NHWC : CUDNN_TENSOR_NCHW)
    xsize = (channelmajor ? size(x)[2:end-1] : size(x)[1:end-2])
    window = min.(window, xsize .+ 2 .* padding)
    if ndims(x) == 3  # cudnnPoolingForward converts to 4D
        window, padding, stride = (x->x[1]).((window, padding, stride))
    end
    y = cudnnPoolingForward(x; mode, nanOpt, window, padding, stride, format, alpha)
end


text:
include("conjugategradient/root.jl")
include("quasinewton/root.jl")


text:
const KEY = :OOLONG

"""
Similar to `Future`, but we added some customized methods.
"""
struct Promise
    f::Future
    function Promise(args...)
        new(Future(args...))
    end
end

Base.getindex(p::Promise) = getindex(p.f)
Base.wait(p::Promise) = wait(p.f)

"Recursively fetch inner value"
function Base.getindex(p::Promise, ::typeof(!))
    x = p.f[]
    while x isa Promise || x isa Future
        x = x[]
    end
    x
end

function Base.getindex(ps::Vector{Promise})
    res = Vector(undef, length(ps))
    @sync for (i, p) in enumerate(ps)
        Threads.@spawn begin
            res[i] = p[]
        end
    end
    res
end

struct TimeOutError{T} <: Exception
    t::T
end

Base.showerror(io::IO, err::TimeOutError) = print(io, "failed to complete in $(err.t) seconds")

"""
    p::Promise[t::Number]

Try to fetch value during a period of `t`.
A [`TimeOutError`](@ref) is thrown if the underlying data is still not ready after `t`.
"""
function Base.getindex(p::Promise, t::Number, pollint=0.1)
    res = timedwait(t;pollint=pollint) do
        isready(p)
    end
    if res === :ok
        p[]
    else
        throw(TimeOutError(t))
    end
end

Base.put!(p::Promise, x) = put!(p.f, x)
Base.isready(p::Promise) = isready(p.f)

#####

struct PotNotRegisteredError <: Exception
    pid::PotID
end

Base.showerror(io::IO, err::PotNotRegisteredError) = print(io, "can not find any pot associated with the pid: $(err.pid)")

#####

const RESOURCE_REGISTRY = Dict{Symbol, UInt}(
    :cpu => () -> Threads.nthreads(),
    :gpu => () -> length(CUDA.devices())
)

struct ResourceInfo{I<:NamedTuple}
    info::I
end

ResourceInfo() = ResourceInfo(NamedTuple(k=>v() for (k,v) in RESOURCE_REGISTRY))

ResourceInfo(;kw...) = ResourceInfo(kw.data)
Base.keys(r::ResourceInfo) = keys(r.info)
Base.getindex(r::ResourceInfo, x) = getindex(r.info, x)
Base.haskey(r::ResourceInfo, x) = haskey(r.info, x)

function Base.:(<=)(x::ResourceInfo, y::ResourceInfo)
    le = true
    for k in keys(x)
        if haskey(y, k) && x[k] <= y[k]
            continue
        else
            le = false
            break
        end
    end
    le
end

function Base.:(-)(x::ResourceInfo, y::ResourceInfo)
    merge(x, (k => x[k]-v for (k,v) in pairs(y)))
end

struct RequirementNotSatisfiedError <: Exception
    required::ResourceInfo
    remaining::ResourceInfo
end

Base.showerror(io::IO, err::RequirementNotSatisfiedError) = print(io, "required: $(err.required), remaining: $(err.remaining)")

#####

"System level messages are processed immediately"
abstract type AbstractSysMsg end

is_prioritized(msg) = false
is_prioritized(msg::AbstractSysMsg) = true

# !!! force system level messages to be executed immediately
# directly copied from
# https://github.com/JuliaLang/julia/blob/6aaedecc447e3d8226d5027fb13d0c3cbfbfea2a/base/channels.jl#L13-L31
# with minor modification
function Base.put_buffered(c::Channel, v)
    lock(c)
    try
        while length(c.data) == c.sz_max
            Base.check_channel_state(c)
            wait(c.cond_put)
        end
        if is_prioritized(v)
            pushfirst!(c.data, v)  # !!! force sys msg to be handled immediately
        else
            push!(c.data, v)  # !!! force sys msg to be handled immediately
        end
        # notify all, since some of the waiters may be on a "fetch" call.
        notify(c.cond_take, nothing, true, false)
    finally
        unlock(c)
    end
    return v
end

#####

"Similar to `RemoteException`, except that we need the `PotID` info."
struct Failure <: Exception
    pid::PotID
    captured::CapturedException
end

Failure(captured) = Failure(self(), captured)

is_prioritized(::Failure) = true

function Base.showerror(io::IO, f::Failure)
    println(io, "In pot $(f.pid) :")
    showerror(io, re.captured)
end


text:
# Communicator

const NCCL_UNIQUE_ID_BYTES = 128
const ncclUniqueId_t = NTuple{NCCL_UNIQUE_ID_BYTES, Cchar}

"""
    UniqueID()

Generates an Id to be used in `Communicator`. `UniqueID` should be
called once and the Id should be distributed to all ranks in the
communicator before calling `Communicator(nranks, uid, rank)`.
"""
struct UniqueID
    internal::ncclUniqueId_t

    function UniqueID()
        buf = zeros(Cchar, NCCL_UNIQUE_ID_BYTES)
        @apicall(:ncclGetUniqueId, (Ptr{Cchar},), buf)
        new(Tuple(buf))
    end
end

Base.convert(::Type{ncclUniqueId_t}, id::UniqueID) = id.internal

const ncclComm_t = Ptr{Cvoid}

struct Communicator
    handle::ncclComm_t
    Communicator(handle::ncclComm_t) = new(handle)
end

# creates a new communicator (multi thread/process version)
"""
   Communicator(nranks, uid, rank)

Creates a new Communicator (multi thread/process version)
`rank` must be between `0` and `nranks-1` and unique within a communicator
clique. Each rank is associated to a CUDA device which has to be set before
calling `Communicator`. Implicitly synchroniszed with other ranks so it must
be called by different threads/processes or used within `group`.
"""
function Communicator(nranks, comm_id, rank)
    handle_ref = Ref{ncclComm_t}()
    @apicall(:ncclCommInitRank, (Ptr{ncclComm_t}, Cint, ncclUniqueId_t, Cint),
             handle_ref, nranks, comm_id, rank)

    Communicator(handle_ref[])
end

Base.cconvert(::Type{ncclComm_t}, comm::Communicator) = comm.handle

# creates a clique of communicators (single process version)
function Communicator(devices)
    ndev = length(devices)
    comms = Vector{ncclComm_t}(undef, ndev)
    devlist = Cint[dev.handle for dev in devices]
    @apicall(:ncclCommInitAll, (Ptr{ncclComm_t}, Cint, Ptr{Cint}), comms, ndev, devlist)
    Communicator.(comms)
end

function destroy(comm::Communicator)
    @apicall(:ncclCommDestroy, (ncclComm_t,), comm)
end

function abort(comm::Communicator)
    @apicall(:ncclCommAbort, (ncclComm_t,), comm)
end

function getError(comm::Communicator)
    ref = Ref{ncclResult_t}()
    @apicall(:ncclCommGetAsyncError, (ncclComm_t, Ref{ncclResult_t}), comm, ref)
    return NCCLError(ref[])
end

function count(comm::Communicator)
    ref = Ref{Cint}()
    @apicall(:ncclCommCount, (ncclComm_t, Ref{Cint}), comm, ref)
    return ref[]
end

function rank(comm::Communicator)
    ref = Ref{Cint}()
    @apicall(:ncclUserRank, (ncclComm_t, Ref{Cint}), comm, ref)
    return ref[]
end

function device(comm::Communicator)
    ref = Ref{Cint}()
    @apicall(:ncclCommCuDevice, (ncclComm_t, Ref{Cint}), comm, ref)
    return CUDAdrv.CuDevice(ref[])
end


text:
export stretched, shifted, rotated, renormed

"""
    $(TYPEDEF)
Abstract type for image modifiers. These are some model wrappers
that can transform any model using simple Fourier transform properties.
To see the implemented modifier
"""
abstract type AbstractModifier{M<:AbstractModel} <: AbstractModel end

"""
    $(SIGNATURES)
Returns the base model from a modified model.
"""
basemodel(model::AbstractModifier) = model.model
basemodel(model::AbstractModel) = model

unmodified(model::AbstractModel) = basemodel(model)
unmodified(model::AbstractModifier) = unmodified(basemodel(model))

flux(m::AbstractModifier) = flux(m.model)

@inline visanalytic(::Type{<:AbstractModifier{M}}) where {M} = visanalytic(M)
@inline imanalytic(::Type{<:AbstractModifier{M}}) where {M} = imanalytic(M)

radialextent(m::AbstractModifier) = radialextent(basemodel(m))

# @inline function apply_uv_transform(m::AbstractModifier, u::Number, v::Number, scale::Number)
#     ut, vt = transform_uv(m, u, v)
#     scale *= scale_uv(m, u, v)
#     return apply_uv_transform(basemodel(m), ut, vt, scale)
# end

# @inline function apply_uv_transform(::AbstractModel, u::Number, v::Number, scale::Number)
#     return u, v, scale
# end

# function apply_uv_transform(m::AbstractModifier, u::AbstractVector, v::AbstractVector)
#     ut = similar(u)
#     vt = similar(v)
#     T = typeof(scale_uv(m, 0.0, 0.0))
#     scale = ones(T, size(ut))
#     @inbounds for i in eachindex(u,v)
#         up, vp, sp = apply_uv_transform(m, u[i], v[i], scale[i])
#         ut[i] = up
#         vt[i] = vp
#         scale[i] = sp
#     end
#     return ut, vt, scale
# end

#@inline function _visibilities(m::AbstractModifier, u::AbstractArray, v::AbstractArray, args...)
#    ut, vt, scales = apply_uv_transform(m, u, v)
#    scales.*visibilities(unmodified(m), ut, vt, args...)
#end

function modelimage(::NotAnalytic,
    model::AbstractModifier,
    image::ComradeBase.AbstractIntensityMap, alg::FFTAlg, executor=SequentialEx())

    @set model.model = modelimage(model.model, image, alg, executor)
end

function modelimage(::NotAnalytic,
    model::AbstractModifier,
    image::ComradeBase.AbstractIntensityMap, alg::NUFT,
    executor=SequentialEx())
    _modelimage(model, image, alg, executor)
end

@inline function visibility_point(m::AbstractModifier, u, v, args...)
    ut, vt = transform_uv(m, u, v)
    scale = scale_uv(m, u, v)
    scale*visibility(basemodel(m), ut, vt, args...)
end

@inline function ComradeBase.intensity_point(m::AbstractModifier, x, y)
    xt, yt = transform_image(m, x, y)
    scale = scale_image(m, x, y)
    scale*ComradeBase.intensity_point(basemodel(m), xt, yt)
end

"""
    $(TYPEDEF)
Shifts the model by `Δx` units in the x-direction and `Δy` units
in the y-direction.
"""
struct ShiftedModel{T,M<:AbstractModel} <: AbstractModifier{M}
    model::M
    Δx::T
    Δy::T
end

"""
    $(SIGNATURES)
Shifts the model `m` in the image domain by an amount `Δx,Δy`.
"""
shifted(model, Δx, Δy) = ShiftedModel(model, Δx, Δy)
# This is a simple overload to simplify the type system
shifted(model::ShiftedModel, Δx, Δy) = ShiftedModel(basemodel(model), Δx+model.Δx, Δy+model.Δy)
radialextent(model::ShiftedModel, Δx, Δy) = radialextent(model.model) + hypot(abs(Δx), abs(Δy))

@inline transform_image(model::ShiftedModel, x, y) = (x-model.Δx, y-model.Δy)
@inline transform_uv(model::ShiftedModel, u, v) = (u, v)

@inline scale_image(model::ShiftedModel, x, y) = 1.0
@inline scale_uv(model::ShiftedModel, u, v) = exp(2im*π*(u*model.Δx + v*model.Δy))

"""
    $(TYPEDEF)
Renormalizes the flux of the model to the new value `flux`.
We have also overloaded the Base.:* operator as syntactic sugar
although I may get rid of this.
"""
struct RenormalizedModel{M<:AbstractModel,T} <: AbstractModifier{M}
    model::M
    scale::T
    RenormalizedModel(model::M, f::T) where {M,T} = new{M,T}(model, f)
end

"""
    $(SIGNATURES)
Renormalizes the model `m` to have total flux `flux`.
"""
renormed(model::M, f) where {M<:AbstractModel} = RenormalizedModel(model, f)
Base.:*(model::AbstractModel, f::Real) = renormed(model, f)
Base.:*(f::Real, model::AbstractModel) = renormed(model, f)
Base.:/(f::Real, model::AbstractModel) = renormed(model, inv(f))
Base.:/(model::AbstractModel, f::Real) = renormed(model, inv(f))
# Dispatch on RenormalizedModel so that I just make a new RenormalizedModel with a different f
# This will make it easier on the compiler.
Base.:*(model::RenormalizedModel, f::Real) = renormed(model.model, model.scale*f)
# Overload the unary negation operator to be the same model with negative flux
Base.:-(model::AbstractModel) = renormed(model, -1.0)
flux(m::RenormalizedModel) = m.scale*flux(m.model)

@inline transform_image(::RenormalizedModel, x, y) = (x, y)
@inline transform_uv(::RenormalizedModel, u, v) = (u, v)

@inline scale_image(model::RenormalizedModel, x, y) = model.scale
@inline scale_uv(model::RenormalizedModel, u, v) = model.scale

function _visibilities(m::RenormalizedModel, u::AbstractArray, v::AbstractArray)
    m.scale*visibilities(basemodel(m), u, v)
end

"""
    $(TYPEDEF)
Stretched the model in the x and y directions, i.e. the new intensity is
```math
    I_s(x,y) = 1/(αβ) I(x/α, y/β),
```
where were renormalize the intensity to preserve the models flux.
"""
struct StretchedModel{M<:AbstractModel,T} <: AbstractModifier{M}
    model::M
    α::T
    β::T
end

"""
    $(SIGNATURES)
Stretches the model `m` according to the formula
```math
    I_s(x,y) = 1/(αβ) I(x/α, y/β),
```
where were renormalize the intensity to preserve the models flux.
"""
stretched(model, α, β) = StretchedModel(model, α, β)
radialextent(model::StretchedModel) = hypot(model.α, model.β)*radialextent(basemodel(model))

@inline transform_image(model::StretchedModel, x, y) = (x/model.α, y/model.β)
@inline transform_uv(model::StretchedModel, u, v) = (u*model.α, v*model.β)

@inline scale_image(model::StretchedModel, x, y) = inv(model.α*model.β)
@inline scale_uv(::StretchedModel, u, v) = one(eltype(u))

"""
    $(TYPEDEF)
Type for the rotated model. This is more fine grained constrol of
rotated model. In most use cases the end-user should be using
the `rotate` method e.g.

```julia
rotate(model, ξ)
```
"""
struct RotatedModel{M<:AbstractModel,T} <: AbstractModifier{M}
    model::M
    s::T
    c::T
end
function RotatedModel(model::T, ξ::F) where {T, F}
    s,c = sincos(ξ)
    return RotatedModel(model, s, c)
end

"""
    $(SIGNATURES)
Rotates the model by an amount `ξ` in radians.
"""
rotated(model, ξ) = RotatedModel(model, ξ)
posangle(model::RotatedModel) = atan(model.s, model.c)

@inline function transform_image(model::RotatedModel, x, y)
    s,c = model.s, model.c
    return c*x + s*y, -s*x + c*y
end

@inline function transform_uv(model::RotatedModel, u, v)
    s,c = model.s, model.c
    return c*u + s*v, -s*u + c*v
end

@inline scale_image(model::RotatedModel, x, y) = 1.0
@inline scale_uv(model::RotatedModel, u, v) = 1.0


text:
#
# Orbit descriptions.
#

"""
A supertype for all single-point orbit descriptions.
Parameterized by coordinate frame, floating point type,
mass unit, lenth unit, time unit, epoch type, state type,
and parameter type (in order).
"""
abstract type AbstractOrbit{FR, F, MU, LU, TU, AU, E, S<:AbstractState{F, LU, TU, AU}, P<:ParameterVector{F, MU, LU, TU, AU}} end

"""
An orbit, described by a `StateVector`, with parameters described by a `ParameterVector`.
"""
struct Orbit{FR, F, MU, LU, TU, AU, E, S, P} <: AbstractOrbit{FR, F, MU, LU, TU, AU, E, S, P}
    epoch::E
    state::S
    system::P
end

"""
Returns a default frame.
"""
function defaultframe(system::ParameterVector)
    if system isa R2BPParameters
        return CoordinateFrames.BodycentricInertial
    elseif system isa CR3BPParameters
        return CoordinateFrames.BarycentricRotating
    else
        return Inertial
    end
end

"""
Outer constructor for `Orbit`s.
"""
function Orbit(state::AbstractState, system::ParameterVector, epoch=TAIEpoch(now()); frame = defaultframe(system))

    S = eval(Base.typename(typeof(state)).name)
    P = eval(Base.typename(typeof(system)).name)

    if state isa KeplerianState && !(system isa R2BPParameters)
        throw(ArgumentError("Orbits with $P must have CartesianState state vectors."))
    end

    FR = frame
    F  = promote_type(eltype(state), eltype(system))
    MU = massunit(system)

    if system isa CR3BPParameters
        LU = lengthunit(system)
        TU = timeunit(system)
    else
        LU = lengthunit(state)
        TU = timeunit(state)
    end

    AU = angularunit(state)
    E  = typeof(epoch)

    return Orbit{FR, F, MU, LU, TU, AU, E, S{F, LU, TU, AU}, P{F, MU, LU, TU, AU}}(
        epoch,
        convert(S{F, LU, TU, AU}, state),
        convert(P{F, MU, LU, TU, AU}, system)
    )
end

"""
Shows all `Orbit` instances.
"""
function Base.show(io::IO, orbit::Orbit)
    if system(orbit) isa R2BPParameters
        println(io, conic(orbit), " Restricted Two-body Orbit with eltype $(typeof(orbit).parameters[2])")
    elseif system(orbit) isa CR3BPParameters
        println(io, "Circular Restricted Three-body Orbit with eltype $(typeof(orbit).parameters[2])")
    end

    println(io, "")
    Base.show(io, state(orbit); showfloats=false, space="  ")
    println(io, "")
    Base.show(io, system(orbit); showfloats=false, space="  ")

end
Base.show(io::IO, ::MIME"text/plain", orbit::Orbit) = show(io, orbit)

"""
Returns the epoch (timestamp) for the `Orbit`.
"""
epoch(orbit::Orbit) = orbit.epoch

"""
Returns the state vector for the `Orbit`.
"""
state(orbit::Orbit) = orbit.state

"""
Returns the parameter vector for the `Orbit`.
"""
system(orbit::Orbit) = orbit.system

"""
Returns the `lengthunit` for an `Orbit`.
"""
Base.@pure lengthunit(::Orbit{FR, F, MU, LU, TU, AU}) where {FR, F, MU, LU, TU, AU} = LU

"""
Returns the `timeunit` for an `Orbit`.
"""
Base.@pure timeunit(::Orbit{FR, F, MU, LU, TU, AU}) where {FR, F, MU, LU, TU, AU} = TU

"""
Returns the `angularunit` for an `Orbit`.
"""
Base.@pure angularunit(::Orbit{FR, F, MU, LU, TU, AU}) where {FR, F, MU, LU, TU, AU} = AU

"""
Returns the `massunit` for an `Orbit`.
"""
Base.@pure massunit(::Orbit{FR, F, MU, LU, TU, AU}) where {FR, F, MU, LU, TU, AU} = MU

"""
Returns the `velocityunit` for an `Orbit`.
"""
velocityunit(orbit::Orbit) = lengthunit(orbit) / timeunit(orbit)

"""
Returns the `massparamunit` for an `Orbit`.
"""
massparamunit(orbit::Orbit) = lengthunit(orbit)^3 / timeunit(orbit)^2

"""
Returns the `OrbitalFrame` for the `Orbit`.
"""
Base.@pure frame(::Orbit{FR}) where FR = FR

"""
An alias for `Orbit` instances about `R2BP` systems.
"""
const R2BPOrbit = Orbit{FR, F, MU, LU, TU, AU, E, S, P} where {FR, F, MU, LU, TU, AU, E, S<:Union{CartesianStateVector, KeplerianState}, P<:R2BPParameters}

"""
An alias for `Orbit` instances about `R2BP` systems with `KeplerianState` descriptions.
"""
const KeplerianR2BPOrbit = Orbit{FR, F, MU, LU, TU, AU, E, <:KeplerianState, <:R2BPParameters} where {FR, F, MU, LU, TU, AU, E}

"""
An alias for `Orbit` instances about `R2BP` systems with `CartesianState` descriptions.
"""
const CartesianR2BPOrbit = Orbit{FR, F, MU, LU, TU, AU, E, <:CartesianStateVector, <:R2BPParameters} where {FR, F, MU, LU, TU, AU, E}

"""
An alias for `Orbit` instances about `CR3BP` systems.
"""
const CR3BPOrbit = Orbit{FR, F, MU, LU, TU, AU, E, <:CartesianStateVector, <:CR3BPParameters} where {FR, F, MU, LU, TU, AU, E}

"""
An alias for `Orbit` instances about any systems with `CartesianState` descriptions.
"""
const CartesianOrbit = Orbit{FR, F, MU, LU, TU, AU, E, <:CartesianStateVector, <:ParameterVector} where {FR, F, MU, LU, TU, AU, E}

"""
An alias for `Orbit` instances about any systems with `CartesianStateWithSTM` descriptions.
"""
const CartesianOrbitWithSTM = Orbit{FR, F, MU, LU, TU, AU, E, <:CartesianStateWithSTM, <:ParameterVector} where {FR, F, MU, LU, TU, AU, E}


text:
using MAT
using LinearAlgebra
using BlockArrays

file = matopen("./examples/MultiscaleLowRank/MultiscaleLowRankRandom.mat")
X_decom = read(file, "X_decom")
close(file)

x1 = X_decom[:,:,1]
x2 = X_decom[:,:,2]
x3 = X_decom[:,:,3]
x4 = X_decom[:,:,4]
m, n = size(x1)

# first component
n = 64
kr1 = count(x->x!=0, x1)
x1 = vec(x1)

# second component
bn2 = 4; l = floor(Int, n/bn2)
Xb2 = BlockArray(x2, repeat([bn2], l), repeat([bn2], l))
idx2 = CartesianIndices((l, l))
kr2 = map(i->rank(Xb2[Block(i[1], i[2])]), idx2)
kr2 = vec(kr2)
x2 = vec(x2)

# third component
bn3 = 16; l = floor(Int, n/bn3)
Xb3 = BlockArray(x3, repeat([bn3], l), repeat([bn3], l))
idx3 = CartesianIndices((l, l))
kr3 = map(i->rank(Xb3[Block(i[1], i[2])]), idx3)
kr3 = vec(kr3)
x3 = vec(x3)

# fourth component
kr4 = rank(x4)
x4 = vec(x4)

# observation
b = x1 + x2 + x3 + x4

# write to file
f = open("/home/zhenan/Github/AtomicOpt.jl/examples/MultiscaleLowRank/MultiscaleData.txt", "w")
l = length(b)
@printf(f, "x1 = ")
for i in 1:l
    @printf(f, "%.3f ", x1[i])
end
println(f, "\n")
@printf(f, "x2 = ")
for i in 1:l
    @printf(f, "%.3f ", x2[i])
end
println(f, "\n")
@printf(f, "x3 = ")
for i in 1:l
    @printf(f, "%.3f ", x3[i])
end
println(f, "\n")
@printf(f, "x4 = ")
for i in 1:l
    @printf(f, "%.3f ", x4[i])
end
println(f, "\n")
@printf(f, "b = ")
for i in 1:l
    @printf(f, "%.3f ", b[i])
end
println(f, "\n")
@printf(f, "kr1 = %d\n", kr1)
@printf(f, "kr2 = ")
for t in kr2
    @printf(f, "%d ", t)
end
println(f, "\n")
@printf(f, "kr3 = ")
for t in kr3
    @printf(f, "%d ", t)
end
println(f, "\n")
@printf(f, "kr4 = %d\n", kr4)
@printf(f, "m = %d\n", m)
@printf(f, "n = %d\n", n)
close(f)


text:
using ClimateMachine.BalanceLaws

abstract type AbstractFluidModel <: BalanceLaw end

"""
    ModelSetup <: AbstractFluidModel
"""
struct ModelSetup{𝒯, 𝒰, 𝒱, 𝒲, 𝒳} <: AbstractFluidModel
    physics::𝒯
    boundary_conditions::𝒰
    initial_conditions::𝒱
    numerics::𝒲
    parameters::𝒳
end

function ModelSetup(; physics, boundary_conditions, initial_conditions, numerics, parameters)
    return ModelSetup(
        physics,
        unpack_boundary_conditions(boundary_conditions),
        initial_conditions,
        numerics,
        parameters,
    )
end

"""
    DryAtmosModel <: AbstractFluidModel

    temporarily use this struct
"""
Base.@kwdef struct DryAtmosModel{𝒯, 𝒰, 𝒱, 𝒲, 𝒳} <: AbstractFluidModel
    physics::𝒯
    boundary_conditions::𝒰
    initial_conditions::𝒱
    numerics::𝒲
    parameters::𝒳
end

"""
    DryAtmosLinearModel <: AbstractFluidModel

    temporarily use this struct
"""
Base.@kwdef struct DryAtmosLinearModel{𝒯, 𝒰, 𝒱, 𝒲, 𝒳} <: AbstractFluidModel
    physics::𝒯
    boundary_conditions::𝒰
    initial_conditions::𝒱
    numerics::𝒲
    parameters::𝒳
end

function unpack_boundary_conditions(bcs)
    # We need to repackage the boundary conditions to match the
    # boundary conditions interface of the Balance Law and DGModel
    boundaries = (:west, :east, :south, :north, :bottom, :top)
    repackaged_bcs = []

    for boundary in boundaries
        fields = get(bcs, boundary, nothing)
        new_bc = isnothing(fields) ? FluidBC() : FluidBC(fields...)
        push!(repackaged_bcs, new_bc)
    end

    return Tuple(repackaged_bcs)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package RadiationDetectorDSP" begin

# test code goes here

end # testset


text:
# Simply supported square plate with center concentrated load
module simply_supp_square_plate_conc_examples

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute_dsg_model(formul, n = 2, visualize = true)
    E = 30e6;
    nu = 0.3;
    force = 40;
    thickness = 0.1;
    L= 10;
    # analytical solution for the vertical deflection under the load
    analyt_sol=-0.0168;

    tolerance = L/n/1000
    fens, fes = T3block(L/2,L/2,n,n);
    fens.xyz = xyz3(fens)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

        sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 0 L/2 -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[0 L/2 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # simple support
    l1 = selectnode(fens; box = Float64[L/2 L/2 0 L/2 -Inf Inf], inflate = tolerance)
    for i in [3,4,6]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[0 L/2 L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [3,5,6]
        setebc!(dchi, l1, true, i)
    end
    # in-plane, rotations
    l1 = selectnode(fens; box = Float64[0 L/2 0 L/2 -Inf Inf], inflate = tolerance)
    for i in [1, 2, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[0 0 0 0 -Inf Inf], tolerance = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    fi = ForceIntensity(FFlt[0, 0, -force/4, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    @show dchi.values[nl, 3]/analyt_sol*100

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return true
end

function test_convergence()
    formul = FEMMShellT3FFModule
    @info "Simply supported square plated with concentrated force,"
    @info " formulation=$(formul)"

    for n in [2, 4, 8, 16, 32, 64]
        _execute_dsg_model(formul, n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
module Cholesky

export DenseBlockAngular, FactorBlockAngular

"""
    cholesky(A, d)
    Form the matrix S = A*Diag(d)*A' and computes its Cholesky factorization.
    The factorization overwrites F.
"""
function cholesky(A::AbstractMatrix{Tv}, d::AbstractVector{Td}) where{Tv<:Real, Td<:Real}
    warn("cholesky! must be implemented for $(typeof(A))")
    F = cholfact(Symmetric(A*spdiagm(d)*A'))
    return F
end

function cholesky(A::Matrix{Ta}, d::AbstractVector{Td}) where{Ta<:Real, Td<:Real}
    F = cholfact(Symmetric(A*Diagonal(d)*A'))
    return F
end

function cholesky(A::SparseMatrixCSC{Ta, Int64}, d::AbstractVector{Td}) where{Ta<:Real, Td<:Real}
    F = cholfact(Symmetric(A*spdiagm(d)*A'))
    return F
end

"""
    cholesky!(A, d, F)
    Form the matrix S = A*Diag(d)*A' and computes its Cholesky factorization.
    The factorization overwrites F.
"""
function cholesky!(A::AbstractMatrix{Ta}, d::AbstractVector{Td}, F::Factorization{Ta}) where{Ta<:Real, Td<:Real}
    warn("cholesky! must be implemented for $(typeof(A))")
    return F
end

function cholesky!(A::Matrix{Ta}, d::AbstractVector{Td}, F::Base.LinAlg.Cholesky{Ta,Matrix{Ta}}) where{Ta<:Real, Td<:Real}
    F = cholfact!(F, A*Diagonal(d)*A')
    return F
end

function cholesky!(A::SparseMatrixCSC{Ta, Int64}, d::AbstractVector{Td}, F::Base.SparseArrays.CHOLMOD.Factor{Ta}) where{Ta<:Real, Td<:Real}
    # update CHolesky factor
    F = cholfact!(F, Symmetric(A*spdiagm(d)*A'))
    return F
end

"""
    solve!(y, F, b)
    Solves the system F*y = b and overwrites y.
"""
function solve!(y::AbstractVector{T1}, F::Factorization{T2}, b::AbstractVector{T3}) where {T1<:Real, T2<:Real, T3<:Real}
    warn("Implement solve! for concrete type implementations!")
    return y
end

# include("denseCholesky.jl")
# include("sparseCholesky.jl")
include("denseBlockAngular.jl")

end  # module


text:
import PGFPlotsX

function get_underlying(p::PGFPlotsXPlot)
    result = p.underlying_object
    return result
end

function Base.display(p::PGFPlotsXPlot)::Nothing
    if is_runtests() && !open_plots_during_tests()
        @debug(
            string(
                "PREDICTMD_OPEN_PLOTS_DURING_TESTS is false, therefore ",
                "the plot will not be opened.",
                )
            )
    else
        @debug(string("Attempting to display plot..."))
        try
            Base.display(get_underlying(p))
            @debug(string("Displayed plot."))
        catch e
            handle_plotting_error(e)
        end
    end
    return nothing
end

function PGFPlotsX.save(
        filename::String,
        p::PGFPlotsXPlot;
        kwargs...,
        )::Nothing
    underlying_object = get_underlying(p)
    try
        @debug(string("Attempting to save plot..."))
        mkpath(dirname(filename))
        PGFPlotsX.save(
            filename,
            underlying_object;
            kwargs...,
            )
        @debug(string("Saved plot to file: \"", filename, "\"",))
    catch e
        handle_plotting_error(e)
    end
    return nothing
end

function save_plot(
        filename::String,
        p::PGFPlotsXPlot;
        kwargs...,
        )::Nothing
    PGFPlotsX.save(
        filename,
        p;
        kwargs...,
        )
    return nothing
end


text:
########################################################################
# This package is based on the neo4j-contrib/boltkit package on GitHub #
# It is also greatly influence by the JuliaIO/MsgPack.jl Julia package #
########################################################################

__precompile__()
module PackageStream

# Imports/Exports
import Base.showerror
export pack, unpack

# Type of packable values
Atomic = Union{Void, Bool, Integer, AbstractFloat, AbstractString}
"""
    Packable type of data values
    1) Atomic values: `Union{Void, Bool, Integer, AbstractFloat, AbstractString}`
    2) Packable values: `Union{Atomic, Vector, Dict, Tuple}`
"""
Packable = Union{Atomic, Vector, Dict, Tuple}

# Dictionary to handle markers and dispatch associated with packed items
const markers = Dict{Symbol, Tuple{Vector{UInt8}, Function}}(
    :Nothing    => (b"\xC0", s -> nothing),
    :False      => (b"\xC2", s -> false),
    :True       => (b"\xC3", s -> true),
    :Int8       => (b"\xC8", s -> readi(s, Int8)),
    :Int16      => (b"\xC9", s -> readi(s, Int16)),
    :Int32      => (b"\xCA", s -> readi(s, Int32)),
    :Int64      => (b"\xCB", s -> readi(s, Int64)),
    :Float      => (b"\xC1", s -> readn(s, Float64)),
    :String8    => (b"\xD0", s -> unpackstring(s, readn(s, UInt8))),
    :String16   => (b"\xD1", s -> unpackstring(s, readn(s, UInt16))),
    :String32   => (b"\xD2", s -> unpackstring(s, readn(s, UInt32))),
    :List8      => (b"\xD4", s -> unpacklist(s, readn(s, UInt8))),
    :List16     => (b"\xD5", s -> unpacklist(s, readn(s, UInt16))),
    :List32     => (b"\xD6", s -> unpacklist(s, readn(s, UInt32))),
    :Map8       => (b"\xD8", s -> unpackmap(s, readn(s, UInt8))),
    :Map16      => (b"\xD9", s -> unpackmap(s, readn(s, UInt16))),
    :Map32      => (b"\xDA", s -> unpackmap(s, readn(s, UInt32))),
    :Struct8    => (b"\xDC", s -> unpackstruct(s, readn(s, UInt8))),
    :Struct16   => (b"\xDD", s -> unpackstruct(s, readn(s, UInt16))),
)

# Error specific to the package
struct ValueError <: Exception
    content::AbstractString
end
Base.showerror(io::IO, e::ValueError) = print(io, e.content)

# Function to write an item that is big enough to need a head marker
function writehead(
    stream::IOBuffer,
    head::Tuple{Vector{UInt8}, Function},
    value::Packable
    )
    writehead(stream, head[1], value)
end
function writehead(
    stream::IOBuffer,
    head::Vector{UInt8},
    value::Packable
    )
    write(stream, head)
    write(stream, hton(value))
end

# Function pack with dispatch over Void value
function pack(
    stream::IOBuffer,
    value::Void
    )
    write(stream, markers[:Nothing][1])
end

# Function pack with dispatch over Bool value
function pack(
    stream::IOBuffer,
    value::Bool
    )
    value ? write(stream, markers[:True][1]) : write(stream, markers[:False][1])
end

# Function pack with dispatch over Integer value
function pack(
    stream::IOBuffer,
    value::Integer
    )
    if -2^4 ≤ value ≤ typemax(Int8) # Tiny int
        write(stream, Int8(value))
    elseif typemin(Int8) ≤ value ≤ typemax(Int8) # Int8
        writehead(stream, markers[:Int8], Int8(value))
    elseif typemin(Int16) ≤ value ≤ typemax(Int16) # Int16
        writehead(stream, markers[:Int16], Int16(value))
    elseif typemin(Int32) ≤ value ≤ typemax(Int32) # Int32
        writehead(stream, markers[:Int32], Int32(value))
    elseif typemin(Int64) ≤ value ≤ typemax(Int64) # Int64
        writehead(stream, markers[:Int64], Int64(value))
    else
        throw(ValueError("Integer value out of packable range"))
    end
end

# Function pack with dispatch over Float value
function pack(
    stream::IOBuffer,
    value::AbstractFloat
    )
    writehead(stream, markers[:Float], value)
end

# Function pack with dispatch over String value
function pack(
    stream::IOBuffer,
    value::AbstractString
    )
    utf_8 = transcode(String, value)
    l = length(utf_8)
    if l < 0x10
        write(stream, UInt8(0x80 + l))
    elseif l < 0x100
        writehead(stream, markers[:String8], UInt8(l))
    elseif l < 0x10000
        writehead(stream, markers[:String16], UInt16(l))
    elseif l < 0x100000000
        writehead(stream, markers[:String32], UInt32(l))
    else
        throw(ValueError("String too long to pack"))
    end
    write(stream, utf_8)
end

# Function pack with dispatch over list (Vector) value
function pack(
    stream::IOBuffer,
    value::Vector
    )
    l = length(value)
    if l < 0x10
        write(stream, UInt8(0x90 + l))
    elseif l < 0x100
        writehead(stream, markers[:List8], UInt8(l))
    elseif l < 0x10000
        writehead(stream, markers[:List16], UInt16(l))
    elseif l < 0x100000000
        writehead(stream, markers[:List32], UInt32(l))
    else
        throw(ValueError("List too long to pack"))
    end
    foreach(x -> write(stream, pack(x)), value)
end

# Function pack with dispatch over map (Dict) value
function pack(
    stream::IOBuffer,
    value::Dict
    )
    l = length(value)
    if l < 0x10
        write(stream, UInt8(0xA0 + l))
    elseif  l < 0x100
        writehead(stream, markers[:Map8], UInt8(l))
    elseif l < 0x10000
        writehead(stream, markers[:Map16], UInt16(l))
    elseif l < 0x100000000
        writehead(stream, markers[:Map32], UInt32(l))
    else
        throw(ValueError("Dictionary too long to pack"))
    end
    foreach(x -> (write(stream, pack(x[1])); write(stream, pack(x[2]))), value)
end

# Function pack with dispatch over struct (Tuple) value
function pack(
    stream::IOBuffer,
    value::Tuple
    )
    signature, fields = value[1], value[2:end]
    l = length(fields)
    if l < 0x10
        write(stream, UInt8(0xB0 + l))
    elseif l < 0x100
        writehead(stream, markers[:Struct8], UInt8(l))
    elseif l < 0x10000
        writehead(stream, markers[:Struct16], UInt16(l))
    else
        throw(ValueError("Structure too big to pack"))
    end
    write(stream, UInt8(signature))
    foreach(x -> write(stream, pack(x)), fields)
end

"""
    pack(value::Packable)

Return a Byte Array (`Vector{UInt8}`) encoded through the PackageStream
specification.

The value type must be either of the following:
* The Nil value `nothing`
* A boolean value `true` or `false`
* An integer up to 64 bits `Int64`
* A float up to 64 bits `Float64`
* A string (up to 2^32 characters) `<:AbstractString`
* A list (up to 2^32 elements) `Vector`
* A map of keys and values (up to 2^32 elements) `Dict`
* A structure (signature, field 1, ..., field n) (up to 2^32 fields) `Tuple`
"""
function pack(
    value::Packable
    )
    stream = IOBuffer()
    pack(stream, value)
    return take!(stream)
end

# Read functions from big-endian to little-endian
function readn(
    stream::IOBuffer,
    t::DataType
    )
    return ntoh(read(stream, t))
end
function readi(
    stream::IOBuffer,
    t::DataType
    )
    return Int64(readn(stream, t))
end

# Sub-routine to unpack specific items: strings, lists, maps, structs
function unpackstring(
    stream::IOBuffer,
    length_values::Integer
    )
    return String(read(stream, length_values))
end
function unpacklist(
    stream::IOBuffer,
    length_values::Integer
    )
    return [unpack(stream) for i in 1:length_values]
end
function unpackmap(
    stream::IOBuffer,
    length_values::Integer
    )
    out = Dict()
    for i in 1:length_values
        k = unpack(stream)
        v = unpack(stream)
        out[k] = v
    end
    return out
end
function unpackstruct(
    stream::IOBuffer,
    length_values::Integer
    )
    return Tuple([unpack(stream) for i in 1:length_values+1])
end

# Unpack core dispatch (stream argument)
function unpack(
    stream::IOBuffer
    )
    headbyte = read(stream, UInt8)
    for (mark, dispatch) in values(markers)
        if mark[1] == headbyte
            return dispatch(stream)
        end
    end
    if -2^4 ≤ (hb = reinterpret(Int8, headbyte); hb) ≤ typemax(Int8)
        return Int64(hb)
    elseif 0x80 ≤ headbyte < 0x90
        unpackstring(stream, headbyte - 0x80)
    elseif 0x90 ≤ headbyte < 0xa0
        unpacklist(stream, headbyte - 0x90)
    elseif 0xa0 ≤ headbyte < 0xb0
        unpackmap(stream, headbyte - 0xa0)
    elseif 0xb0 ≤ headbyte < 0xc0
        unpackstruct(stream, headbyte - 0xb0)
    else
        throw(ValueError("Unknow marker byte $headbyte"))
    end
end

"""
    unpack(byte_array)

Argument is a Byte Array (`Vector{UInt8}`) following the encoding the
PackageStream specification for Neo4j.

The return value type can be either of the following:
* The Nil value `nothing`
* A boolean value `true` or `false`
* An integer up to 64 bits `Int64`
* A float up to 64 bits `Float64`
* A string (up to 2^32 characters) `<:AbstractString`
* A list (up to 2^32 elements) `Vector`
* A map of keys and values (up to 2^32 elements) `Dict`
* A structure (signature, field 1, ..., field n) (up to 2^32 fields) `Tuple`
"""
function unpack(
    stream::Vector{UInt8}
    )
    return unpack(IOBuffer(stream))
end

end # module


text:
using GlobalSensitivity, SafeTestsets
using Test

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

@time begin

if GROUP == "All" || GROUP == "GSA"
    @time @safetestset "Morris Method" begin include("morris_method.jl") end
    @time @safetestset "Sobol Method" begin include("sobol_method.jl") end
    @time @safetestset "DGSM Method" begin include("DGSM.jl") end
    @time @safetestset "eFAST Method" begin include("eFAST_method.jl") end
    @time @safetestset "RegressionGSA Method" begin include("regression_sensitivity.jl") end
    @time @safetestset "DeltaMoment Method" begin include("delta_method.jl") end
end

end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Logging
using JLD2
using BenchmarkTools
using Benchmarks

using Oceananigans
using Oceananigans.Models
using CUDA

Logging.global_logger(OceananigansLogger())

Nx = parse(Int, ARGS[1])
Ny = parse(Int, ARGS[2])

T = Threads.nthreads()

@info "Setting up threaded serial shallow water model with N=($Nx, $Ny) grid points and $T threads..."

topo = (Periodic, Periodic, Flat)   # Use Flat
grid = RectilinearGrid(topology=topo, size=(Nx, Ny), extent=(1, 1), halo=(3, 3))
model = ShallowWaterModel(architecture=CPU(), grid=grid, gravitational_acceleration=1.0)
set!(model, h=1.0)

@info "Warming up serial shallow water model..."

time_step!(model, 1) # warmup

@info "Benchmarking serial shallow water model..."

trial = @benchmark begin
    @sync_gpu time_step!($model, 1)
    #CUDA.@sync blocking=true time_step!($model, 1)
end samples=10 evals=1

t_median = BenchmarkTools.prettytime(median(trial).time)
@info "Done benchmarking. Median time: $t_median"

jldopen("distributed_shallow_water_model_threads$T.jld2", "w") do file
    file["trial"] = trial
end


text:
function report_routing_stats(m::Map{A,D}) where {A,D}
    # Get the link histogram first - this will make later analysis steps much
    # easier as we won't have to traverse through the architecture a bunch
    # of times.
    histogram = global_link_histogram(m)

    # Print the total number of communication links in the taskgraph
    num_links = sum(values(histogram))
    printstyled("Number of communication channels: "; color = :yellow)
    printstyled(num_links, "\n"; color = :white)

    # Report the total number of links used - this will be the sum of the
    # number of links with a given length multiplied by the length.
    total = sum(k * v for (k, v) in histogram)
    printstyled("Total global routing links used: "; color = :yellow)
    printstyled(total, "\n"; color = :white)

    # Report the average link length
    average_length = total / num_links
    printstyled("Average Link Length: "; color = :yellow)
    printstyled(average_length, "\n"; color = :white)

    # Print the maximum link distance
    printstyled("Maximum Link Distance: "; color = :yellow)
    printstyled(maximum(keys(histogram)), "\n"; color = :white)

    # Display the link histogram
    printstyled("Link Histogram: \n"; color = :yellow)
    display(histogram)
    println()

    return histogram
end

function total_global_links(m::Map{A,D}) where {A,D}
    count = 0
    for edge in m.mapping.edges
        count += count_global_links(edge)
    end
    return count
end

function global_link_histogram(m::Map{A,D}) where {A,D}
    histogram = SortedDict{Int64,Int64}()
    for edge in m.mapping.edges
        link_count = count_global_links(edge)
        add_to_dict(histogram, link_count)
    end
    return histogram
end

function count_global_links(path)
    link_count = 0
    for v in vertices(path)
        if isgloballink(v)
            link_count += 1
        end
    end
    return link_count
end


text:
# ------------------------------------------------------------------------------------------------------------------------------- #
# compute_log_return_array(data::DataFrame,map::Pair{Symbol,Symbol}; Δt::Float64 = 1.0) -> DataFrame

# Computes μ of historical price data.

# # Arguments:
# data_table  	DataFrame holding the historical price data
# map 	        Different data APIs return data with different field names. The map arg connects the time field to price field
# Δt 		    Time step size. Default is 1.0 in the units you want the μ to be calculated in
# ------------------------------------------------------------------------------------------------------------------------------ #
function compute_log_return_array(data_table::DataFrame, map::Pair{Symbol,Symbol}; Δt = (1.0 / 365.0))

    # initialize -
    (number_of_rows, _) = size(data_table)
    return_table = DataFrame(timestamp = Date[], P1 = Float64[], P2 = Float64[], μ = Float64[])

    # main loop -
    for row_index = 2:number_of_rows

        # grab the date -
        tmp_date = data_table[row_index, map.first]

        # grab the price data -
        yesterday_close_price = data_table[row_index-1, map.second]
        today_close_price = data_table[row_index, map.second]

        # compute the diff -
        μ = (1 / Δt) * log(today_close_price / yesterday_close_price)

        # push! -
        push!(return_table, (tmp_date, yesterday_close_price, today_close_price, μ))
    end

    # return -
    return return_table
end

function compute_log_return_array(ticker_symbol_array::Array{String,1}, data_tables::Dict{String,DataFrame},
    map::Pair{Symbol,Symbol}; Δt = (1.0 / 365.0))

    # initialize -
    number_of_ticker_symbols = length(ticker_symbol_array)
    data_dictionary = Dict{String,DataFrame}()

    # process -
    for ticker_symbol_index ∈ 1:number_of_ticker_symbols

        # get the ticker symbol for this index -
        ticker_symbol = ticker_symbol_array[ticker_symbol_index]

        # get the data frame for this index -
        df = data_tables[ticker_symbol]

        # compute the return -
        data_dictionary[ticker_symbol] = compute_log_return_array(df, map; Δt = Δt)
    end

    # return -
    return data_dictionary
end

function compute_fractional_return_array(data::Array{Float64,1};
    multiplier::Float64 = 1.0)

    # what is the length of the array?
    number_of_steps = length(data)
    Δ_array = Array{Float64,1}()

    for step_index ∈ 2:number_of_steps

        # compute -
        P₁ = data[step_index-1]
        P₂ = data[step_index]
        Δ = multiplier*(P₂ - P₁)/(P₁)

        # grab -
        push!(Δ_array,Δ)
    end

    # return -
    return Δ_array;
end

function compute_fractional_return_array(data::Array{Float64,2};
    multiplier::Float64 = 1.0)

    # what is the length of the array?
    (number_of_steps, number_of_cols) = size(data)
    Δ_array = Array{Float64,2}(undef, number_of_steps - 1, number_of_cols)

    for col_index ∈ 1:number_of_cols
        for step_index ∈ 2:number_of_steps

            # compute -
            P₁ = data[step_index-1,col_index]
            P₂ = data[step_index, col_index]
            Δ = multiplier*(P₂ - P₁)/(P₁)

            # grab -
            Δ_array[step_index - 1, col_index] = Δ
        end
    end

    # return -
    return Δ_array;
end

function compute_fractional_return_array(data_table::DataFrame, map::Pair{Symbol,Symbol};
    multiplier::Float64=1.0)

    # initialize -
    (number_of_rows, _) = size(data_table)
    return_table = DataFrame(timestamp = Date[], P1 = Float64[], P2 = Float64[], μ = Float64[])

    # main loop -
    for row_index = 2:number_of_rows

        # grab the date -
        tmp_date = data_table[row_index, map.first]

        # grab the price data -
        yesterday_close_price = data_table[row_index-1, map.second]
        today_close_price = data_table[row_index, map.second]

        # compute the diff -
        μ = ((today_close_price/yesterday_close_price) - 1.0)*multiplier;

        # push! -
        push!(return_table, (tmp_date, yesterday_close_price, today_close_price, μ))
    end

    # return -
    return return_table
end

function compute_fractional_return_array(ticker_symbol_array::Array{String,1}, data_tables::Dict{String,DataFrame}, map::Pair{Symbol,Symbol};
    multiplier::Float64=1.0)

    # initialize -
    number_of_ticker_symbols = length(ticker_symbol_array)
    data_dictionary = Dict{String,DataFrame}()

    # process -
    for ticker_symbol_index ∈ 1:number_of_ticker_symbols

        # get the ticker symbol for this index -
        ticker_symbol = ticker_symbol_array[ticker_symbol_index]

        # get the data frame for this index -
        df = data_tables[ticker_symbol]

        # compute the return -
        data_dictionary[ticker_symbol] = compute_fractional_return_array(df, map; multiplier=multiplier)
    end

    # return -
    return data_dictionary
end

function compute_discrete_geometric_brownian_motion_trajectory(model::GeometricBrownianMotionModel,
    initial_condition::Float64, number_of_steps::Int64; Δt::Float64 = 1.0, N::Int64 = 100)

    # initialize -
    state_array = Array{Float64,2}(undef, number_of_steps, N)
    state_array[1,:] .= initial_condition

    # get data from the model -
    μ = model.μ
    σ = model.σ

    # calculate the random term -
    d = Normal{Float64}(0.0, 1.0)

    # calculate the drift term -
    drift_term = μ*Δt

    # main: sample paths -
    for sample_path_index ∈ 1:N

        # generate new sequence of steps for this sample path -
        random_terms = σ*sqrt(Δt)*rand(d,number_of_steps)

        # solve the model -
        for time_step_index ∈ 2:number_of_steps

            # grab the current price P₁ -
            P₁ = state_array[time_step_index-1,sample_path_index]

            # compute the price difference -
            ΔP = P₁*(drift_term+random_terms[time_step_index])

            # capture -
            state_array[time_step_index,sample_path_index] = (P₁ + ΔP)
        end
    end

    # return -
    return state_array
end

function compute_rwm_cumulative_probabilty(compare::Function, price_array::Array{Float64,1})

    # initialize -
    number_of_samples = length(price_array)
    tmp_array = BitArray(undef, (number_of_samples, 1))

    # main -
    for sample_index = 1:number_of_samples

        # get the sample price -
        sample_price = price_array[sample_index]

        # check: which is larger, sample or target price?
        compare(sample_price) ? tmp_array[sample_index] = 1 : tmp_array[sample_index] = 0
    end

    # sum the tmp_array -
    number_of_larger_values = sum(tmp_array)

    # compute the probability -
    return (number_of_larger_values / number_of_samples)
end

function compute_rwm_cumulative_probabilty(price_array::Array{Float64,1}, target_price::Float64)

    # initialize -
    number_of_samples = length(price_array)
    tmp_array = Array{Int64,1}()

    # main -
    for sample_index = 1:number_of_samples

        # get the sample price -
        sample_price = price_array[sample_index]

        # check: which is larger, sample or target price?
        sample_price <= target_price ? push!(tmp_array, 1) : push!(tmp_array, 0)
    end

    # sum the tmp_array -
    number_of_larger_values = sum(tmp_array)

    # compute the probability -
    return (number_of_larger_values / number_of_samples)
end

function compute_minvar_portfolio_allocation(μ,Σ,target_return::Float64;
    w_lower::Float64 = 0.0, w_upper::Float64 = 1.0)

    # initialize -
    number_of_assets = length(μ)
    w = Variable(number_of_assets)
    risk = quadform(w,Σ)
    ret  = dot(w,μ)

    # setup problem -
    p = minimize(risk)
    p.constraints += [sum(w)==1.0, w_lower <= w, w <= w_upper, ret >= target_return]
    Convex.solve!(p, SCS.Optimizer(verbose = false))

    # return -
    return (p.status, evaluate(w), p.optval, evaluate(ret))
end

function compute_cybernetic_portfolio_allocation(μ,Σ)

    # how many assets do we have?
    𝒫 = length(μ)

    # initialize -
	term_array = Array{Float64,1}(undef, 𝒫)
	for term_index ∈ 1:𝒫
		term_array[term_index] = max(0.0, μ[term_index]/Σ[term_index,term_index])
	end

    # compute the u-variable -
	𝒵 = sum(term_array)

    # if we have no good options (all the returns are negative, then all allocations would be zero)
    if (𝒵 == 0)
        u_variable_array = zeros(𝒫)
    else
        u_variable_array = (1/𝒵)*term_array
    end

    # return -
    return u_variable_array
end

function compute_average_fractional_return_and_covariance(tickers::Array{String,1}, data::Dict{String,DataFrame},
    start::Date, stop::Date)

    # need to grab the number of time steps -
    test_df = Serenity.extract_data_block_for_date_range(data[tickers[1]], start, stop)
	number_of_time_steps = nrow(test_df)

    # initialize -
    μ_bar = Array{Float64,1}()
    R_array = Array{Float64,2}(undef, number_of_time_steps, length(tickers))

    # process each ticker -
    for (ticker_index, ticker) ∈ enumerate(tickers)

        # get the data for the data range -
        df = data[ticker]
        df_slice = extract_data_block_for_date_range(df,start,stop)

        # compute the return -
		avg_val = mean(df_slice[!,:μ])
		push!(μ_bar, avg_val)

        # build R_array -
		for step_index = 1:number_of_time_steps
			R_array[step_index, ticker_index] = df_slice[step_index,:μ]
		end
    end

    # compute the covariance -
    Σ = cov(R_array)

    # return -
    return (μ_bar,Σ)
end

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity