text:
################################################################################
#
#  Basic field access
#
################################################################################

base_ring(A::AlgGrp{T}) where {T} = A.base_ring::parent_type(T)

Generic.dim(A::AlgGrp) = size(multiplication_table(A, copy = false), 1)

elem_type(::Type{AlgGrp{T, S, R}}) where {T, S, R} = AlgGrpElem{T, AlgGrp{T, S, R}}

group(A::AlgGrp) = A.group

has_one(A::AlgGrp) = true

function (A::AlgGrp{T, S, R})(c::Array{T, 1}) where {T, S, R}
  length(c) != dim(A) && error("Dimensions don't match.")
  return AlgGrpElem{T, typeof(A)}(A, c)
end

function multiplication_table(A::AlgGrp; copy::Bool = true)
  if copy
    return deepcopy(A.mult_table)
  else
    return A.mult_table
  end
end

################################################################################
#
#  Commutativity
#
################################################################################

iscommutative_known(A::AlgGrp) = (A.iscommutative != 0)

function iscommutative(A::AlgGrp)
  if iscommutative_known(A)
    return A.iscommutative == 1
  end
  for i in 1:dim(A)
    for j in 1:dim(A)
      if multiplication_table(A, copy = false)[i, j] != multiplication_table(A, copy = false)[j, i]
        A.iscommutative = 2
        return false
      end
    end
  end
  A.iscommutative = 1
  return true
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::AlgGrp)
  compact = get(io, :compact, false)
  if compact
    print(io, "Group algebra of dimension ", dim(A), " over ", base_ring(A))
  else
    print(io, "Group algebra of group\n", group(A), "\nover\n", base_ring(A))
  end
end

################################################################################
#
#  Deepcopy
#
################################################################################

# TODO: This is broken. I have to copy everything carefully by hand.
#function Base.deepcopy_internal(A::AlgGrp, dict::IdDict)
#  B = AlgGrp(base_ring(A), group(A))
#  for x in fieldnames(typeof(A))
#    if x != :base_ring && x != :group && isdefined(A, x)
#      setfield!(B, x, Base.deepcopy_internal(getfield(A, x), dict))
#    end
#  end
#  B.base_ring = A.base_ring
#  B.group = A.group
#  return B
#end

################################################################################
#
#  Equality
#
################################################################################

function ==(A::AlgGrp{T}, B::AlgGrp{T}) where {T}
  return base_ring(A) == base_ring(B) && group(A) == group(B)
end

###############################################################################
#
#  Trace Matrix
#
###############################################################################

function _assure_trace_basis(A::AlgGrp{T}) where {T}
  if !isdefined(A, :trace_basis_elem)
    A.trace_basis_elem = Vector{T}(undef, dim(A))
    A.trace_basis_elem[1] = base_ring(A)(dim(A))
    for i=2:length(A.trace_basis_elem)
      A.trace_basis_elem[i] = zero(base_ring(A))
    end
  end
  return nothing
end

function trace_matrix(A::AlgGrp)
  _assure_trace_basis(A)
  F = base_ring(A)
  n = dim(A)
  M = zero_matrix(F, n, n)
  for i = 1:n
    M[i,i] = tr(A[i]^2)
  end
  for i = 1:n
    for j = i+1:n
      x = tr(A[i]*A[j])
      M[i,j] = x
      M[j,i] = x
    end
  end
  return M
end

################################################################################
#
#  Center
#
################################################################################

@doc Markdown.doc"""
    center(A::AlgGrp{T}) where T

Returns the center C of A and the inclusion C \to A.
"""
function center(A::AlgGrp{T}) where {T}
  if iscommutative(A)
    B, mB = AlgAss(A)
    return B, mB
  end

  if isdefined(A, :center)
    return A.center::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
  end

  B, mB = AlgAss(A)
  C, mC = center(B)
  mD = compose_and_squash(mB, mC)
  A.center = C, mD
  return C, mD
end

################################################################################
#
#  Conversion to AlgAss
#
################################################################################

function AlgAss(A::AlgGrp{T, S, R}) where {T, S, R}
  K = base_ring(A)
  mult = Array{T, 3}(undef, dim(A), dim(A), dim(A))
  B = basis(A)
  d = dim(A)
  for i in 1:d
    for j in 1:d
      l = multiplication_table(A, copy = false)[i, j]
      for k in 1:d
        if k == l
          mult[i, j, k] = one(K)
        else
          mult[i, j, k] = zero(K)
        end
      end
    end
  end
  B = AlgAss(K, mult, one(A).coeffs)
  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
end

################################################################################
#
#  Misc
#
################################################################################

Base.enumerate(G::Generic.PermGroup) = enumerate(AllPerms(G.n))


text:
using Interpolations
using Polynomials

function ConcreteStressStrainModel(eps_co,eps_cu,f_c,f_t)

    #define concrete model strain range
    eps_c=range(0,eps_cu,length=100)

    #fill initial stress vector
    sig_c=zeros(Float64,length(eps_c))

    #define model for stress-strain relationship
    pre=findall(x->(x<=eps_co), eps_c)
    sig_c[pre]=f_c*(2*eps_c[pre]/eps_co-(eps_c[pre]./eps_co).^2)
    post=findall(x->(x>eps_co), eps_c)
    sig_c[post].=f_c.+((0.15.*f_c)./(eps_co.-eps_cu)).*(eps_c[post].-eps_co)

    #calculate initial elastic modulus
    E_ci=(sig_c[2]-sig_c[1])/(eps_c[2]-eps_c[1])

    #convert signs to represent compression
    eps_c=-eps_c
    sig_c=-sig_c

    #define concrete tension stress-strain curve
    f_cr=f_t  #tensile modulus of rupture
    eps_cr=f_cr/E_ci #strain associated with modulus of rupture

    eps_c=vcat(eps_cr,eps_c)
    # eps_c=vcat(eps_cr*1.000001,eps_c)
    # eps_c=vcat(0.003,eps_c)
    sig_c=vcat(f_cr,sig_c)
    # sig_c=vcat([0, 0],sig_c)

    return eps_c,sig_c,E_ci,eps_cr

end

function SteelStressStrainModel(E_si,sig_sy,sig_su,sig_sf,eps_sy,eps_sy1,eps_su,esp_sf)

    # define elastic part of steel stress-strain curve
    sig_s=[0 sig_sy]
    eps_s=[0 eps_sy]

    #define plateau and hardening part of steel stress-strain curve
    p = fit([eps_sy1, eps_su, eps_sf],[sig_sy, sig_su, sig_sf],2)

    sig_fit = zeros(11)
    eps_range = eps_sy1:(eps_sf-eps_sy1)./10:eps_sf

    for i=1:11
    sig_fit[i] = p(eps_range[i])
    end

    eps_s=vcat(vec(eps_s),collect(range(eps_sy1,length=11,stop=eps_sf)))
    sig_s=vcat(vec(sig_s), sig_fit)

    # add compression part of steel stress-strain curve
    sig_s=vcat(-reverse(sig_s[2:end]),sig_s)
    eps_s=vcat(-reverse(eps_s[2:end]),eps_s)

    return eps_s, sig_s

end

function Response5000(M,P,h,s,w,A_s,z_s,E_ci,E_si,eps_cr,SteelStrain,SteelStress,ConcreteStrain,ConcreteStress)

    # define initial concrete locations in cross-section (from bottom, z=0)
    A_c=ones(s,1)*h/s*w
    z_c=collect(range(h/s*0.5,h-h/s*0.5,length=s))

    # initialize elastic modulus
    E_c=ones(s,1)*E_ci
    E_s=ones(length(A_s),1)*E_si

    SteelFit = LinearInterpolation(SteelStrain, SteelStress)
    ConcreteFit = LinearInterpolation(reverse(ConcreteStrain), reverse(ConcreteStress))

    #initialize stress, strain, and curvature
    eps_P=0
    eps_M_c=zeros(length(A_c),1)
    eps_M_s=zeros(length(A_s),1)
    rho_inv_hist=zeros(length(M))
    rho_inv=0

    sig_c_hist=zeros(length(A_c),1)
    sig_s_hist=zeros(length(z_s),1)
    eps_c_hist=zeros(length(A_c),1)
    eps_s_hist=zeros(length(z_s),1)

    for i=2:length(M)

        # E_c
        # E_s
        # eps_P
        # rho_inv

        # define load increment
        dP=P[i]-P[i-1]
        dM=M[i]-M[i-1]

        # calculate cross-section centroid
        z_bar=(sum(A_c.*z_c.*E_c)+sum(A_s.*z_s.*E_s))./((sum(A_c.*E_c)+sum(A_s.*E_s)))

        # calculate concrete and steel strains from P
        eps_P=eps_P+dP./(sum(E_s.*A_s)+sum(E_c.*A_c))

        # calculate 1 over the radius of curvature from M
        rho_inv=rho_inv.+dM./(sum(A_s.*E_s.*(z_s.-z_bar).^2).+sum(A_c.*E_c.*(z_c.-z_bar).^2))

        rho_inv_hist[i]=rho_inv

        # calculate strains from M
        eps_M_c=-(z_c.-z_bar).*rho_inv   #strain in concrete
        eps_M_s=-(z_s.-z_bar).*rho_inv   #strain in steel

        # combine axial and flexural strain, calculate strain history
        eps_c_hist=eps_M_c.+eps_P  #total strain in concrete
        eps_s_hist=eps_M_s.+eps_P  #total strain in steel

        index=findall(x->x<eps_cr,eps_c_hist)
        antindex=findall(x->x>=eps_cr,eps_c_hist)

        if abs(minimum(eps_c_hist))>abs(minimum(ConcreteStrain))
            println("The concrete has crushed, P=",P[i]," M=",M[i])
            return eps_c_hist,eps_s_hist,sig_c_hist,sig_s_hist,z_bar,z_c,rho_inv_hist
        end
        steelfailure=0
        for i=1:length(z_s)
            if abs(eps_s_hist[i])>maximum(SteelStrain)
                steelfailure=1
            else
                steelfailure=0
            end
        end

        if steelfailure==1
            println("The steel has reached its ultimate strain, P=",P[i], " M=",M[i])
            return eps_c_hist,eps_s_hist,sig_c_hist,sig_s_hist,z_bar,z_c,rho_inv_hist
        end

        #store previous stress and strain history
        sig_s_hist_old = sig_s_hist
        sig_c_hist_old = sig_c_hist
        eps_s_hist_old = eps_s_hist
        eps_c_hist_old = eps_c_hist

        # calculate concrete and steel stresses
        sig_c_hist[index]=ConcreteFit(eps_c_hist[index])
        if isempty(antindex)==false
             sig_c_hist[antindex].=0
        end
        sig_s_hist=SteelFit(eps_s_hist)

        #tangent modulus
        E_c=(sig_c_hist .- sig_c_hist_old)./(eps_c_hist .- eps_c_hist_old)
        E_s=(sig_s_hist .- sig_s_hist_old)./(eps_s_hist .- eps_s_hist_old)
        # update elastic modulus
        # E_c=sig_c_hist./eps_c_hist
        # E_s=sig_s_hist./eps_s_hist

    end
    return eps_c_hist,eps_s_hist,sig_c_hist,sig_s_hist,z_bar,z_c,rho_inv_hist, E_c, E_s
end

# define concrete cross-section
w=50  #width
h=50  #depth/height
s=50*4  #number of cross-section divisions

# define rebar area and location (measured from bottom, z=0) in
# cross-section
     # number size location
     Rebar=[10 8 2.5]

RebarQuantity=Rebar[:,1]
RebarDiameter=Rebar[:,2]/8
RebarLocation=Rebar[:,3]

RebarSpacing=w./RebarQuantity

A_s=RebarQuantity.*pi.*(RebarDiameter./2).^2
z_s=RebarLocation

#define concrete compression engineering stress-strain model
#modified Hognestad
eps_co=0.002        #strain at f'c
eps_cu=0.0035       #ultimate strain at brittle failure
f_c=6.00            #concrete compressive strength
f_t=1/10*f_c        #concrete tensile strength

ConcreteStrain,ConcreteStress,E_ci,eps_cr=ConcreteStressStrainModel(eps_co,eps_cu,f_c,f_t)

E_si=29000                  #steel elastic modulus
sig_sy=60                   #steel yield stress
sig_su=80                   #steel ultimate stress
sig_sf=sig_sy*(1+0.10)      #steel fracture stress#
eps_sy=sig_sy/E_si          #steel yield strain
eps_sy1=sig_sy/E_si*(10)    #steel strain at end of yield plateau
eps_su=0.20                 #steel ultimate strain
eps_sf=0.30                 #steel fracture strain

SteelStrain,SteelStress=SteelStressStrainModel(E_si,sig_sy,sig_su,sig_sf,eps_sy,eps_sy1,eps_su,eps_sf)

M=range(0,38798.0,length=100)
P=range(0, 0,length=100)

eps_c_hist,eps_s_hist,sig_c_hist,sig_s_hist,z_bar,z_c,rho_inv_hist, E_c, E_s=Response5000(M,P,h,s,w,A_s,z_s,E_ci,E_si,eps_cr,SteelStrain,SteelStress,ConcreteStrain,ConcreteStress)

plot(eps_c_hist, z_c)

ext: jl


text:
function restart(cur_player::AbstractPlayer)
  end_round(cur_player, true)

  return
end


text:
using Printf

f(x)  = x^2 + sin(x)
df(x) = 2x + cos(x)
a = 1.0

for k = 1:10
    h = 10.0^(-k)
    central_diff = ( f(a+h) - f(a-h) ) / (2h)
    complex_diff = imag( f(complex(a, h)) ) / h
    deriv = df(a)
    @printf("\$10^{-%d}\$& %10.2e& %10.2e\\\\\n",
            k, central_diff-deriv, complex_diff-deriv)
end

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
function rotate(rot, ch::Char)
    for abc in [('a':'z'), ('A':'Z')]
        i = findfirst(x -> x == ch, abc)
        if !isnothing(i)
            return abc[1 + mod(i + rot - 1, length(abc))]
        end
    end
    return ch
end

rotate(rot, string::AbstractString) = map(x -> rotate(rot, x), string)

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# HellwarthTheory.jl

#####
"""
    HellwarthBScheme(LO)

LO an array assumed to be of [freq ; absolute ir activity ]

Multiple phonon mode reduction to a single effective frequency.
Hellwarth et al. 1999 PRB, 'B scheme'; the athermal method.  Averaging procedure is
constructed by considering the average effect of the action of multiple branches.

Follows Eqn. (58) in this paper, assuming typo on LHS, should actually be W_e.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
function HellwarthBScheme(LO)
    println("Hellwarth B Scheme... (athermal)")
    H58 = sum((LO[:, 2].^2) ./ LO[:, 1].^2)
    println("Hellwarth (58) summation: ", H58)

    H59 = sum(LO[:, 2].^2) # sum of total ir activity squarred
    println("Hellwarth (59) summation (total ir activity ^2): ", H59)
    println("Hellwarth (59) W_e (total ir activity ): ", sqrt(H59))

    omega = sqrt(H59 / H58)
    println("Hellwarth (61) Omega (freq): ", omega)

	return(omega)
end

# More complex scheme, involving thermodynamic Beta
# Hellwarth(50), RHS
"""
    HellwarthAScheme(phonon_modes; T = 295, convergence = 1e-6)

Multiple phonon mode reduction to a single effective frequency.
Temperature dependent, defaults to T = 295 K.

Solved iteratively by bisection until Δfreq < convergence.

Follows Hellwarth et al. 1999 PRB 'A' scheme, Eqn. (50) RHS.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
function HellwarthAScheme(phonon_modes; T = 295, convergence=1e-6)

	phonon_mode_freqs = phonon_modes[:, 1]
	ir_activities = phonon_modes[:, 2]

	condition(f) = coth(π * f * 1e12 * ħ / (kB * T)) / f
        - sum(ir_activities .*
              coth.(π .* phonon_mode_freqs .* 1e12 .* ħ ./ (kB * T)) ./ phonon_mode_freqs) / sum(ir_activities)

    # Solve by bisection
	minimum_frequency = minimum(phonon_mode_freqs)
	maximum_frequency = maximum(phonon_mode_freqs)
	middle_frequency = (maximum_frequency + minimum_frequency) / 2
	print("\n")

	while (maximum_frequency - minimum_frequency) / 2 > convergence

		if sign(condition(middle_frequency)) == sign(condition(minimum_frequency))
			minimum_frequency = middle_frequency
			middle_frequency = (maximum_frequency + minimum_frequency) / 2
		else
			maximum_frequency = middle_frequency
			middle_frequency = (maximum_frequency + minimum_frequency) / 2
		end
	end

	return middle_frequency
end


text:
import UUIDs

# This function is based off of a similar function here:
# https://github.com/JuliaRegistries/RegistryCI.jl/blob/master/src/RegistryCI.jl
function gather_stdlib_uuids()
    return Set{UUIDs.UUID}(x for x in keys(Pkg.Types.stdlib()))
end


text:
const KEYWORD_DOCS = """
# Arguments

## File layout options:

  * `header=1`: how column names should be determined; if given as an `Integer`, indicates the row to parse for column names; as an `AbstractVector{<:Integer}`, indicates a set of rows to be concatenated together as column names; `Vector{Symbol}` or `Vector{String}` give column names explicitly (should match # of columns in dataset); if a dataset doesn't have column names, either provide them as a `Vector`, or set `header=0` or `header=false` and column names will be auto-generated (`Column1`, `Column2`, etc.). Note that if a row number header and `comment` or `ignoreemptyrows` are provided, the header row will be the first non-commented/non-empty row _after_ the row number, meaning if the provided row number is a commented row, the header row will actually be the next non-commented row.
  * `normalizenames::Bool=false`: whether column names should be "normalized" into valid Julia identifier symbols; useful when using the `tbl.col1` `getproperty` syntax or iterating rows and accessing column values of a row via `getproperty` (e.g. `row.col1`)
  * `skipto::Integer`: specifies the row where the data starts in the csv file; by default, the next row after the `header` row(s) is used. If `header=0`, then the 1st row is assumed to be the start of data; providing a `skipto` argument does _not_ affect the `header` argument. Note that if a row number `skipto` and `comment` or `ignoreemptyrows` are provided, the data row will be the first non-commented/non-empty row _after_ the row number, meaning if the provided row number is a commented row, the data row will actually be the next non-commented row.
  * `footerskip::Integer`: number of rows at the end of a file to skip parsing.  Do note that commented rows (see the `comment` keyword argument) *do not* count towards the row number provided for `footerskip`, they are completely ignored by the parser
  * `transpose::Bool`: read a csv file "transposed", i.e. each column is parsed as a row
  * `comment::String`: string that will cause rows that begin with it to be skipped while parsing. Note that if a row number header or `skipto` and `comment` are provided, the header/data row will be the first non-commented/non-empty row _after_ the row number, meaning if the provided row number is a commented row, the header/data row will actually be the next non-commented row.
  * `ignoreemptyrows::Bool=true`: whether empty rows in a file should be ignored (if `false`, each column will be assigned `missing` for that empty row)
  * `select`: an `AbstractVector` of `Integer`, `Symbol`, `String`, or `Bool`, or a "selector" function of the form `(i, name) -> keep::Bool`; only columns in the collection or for which the selector function returns `true` will be parsed and accessible in the resulting `CSV.File`. Invalid values in `select` are ignored.
  * `drop`: inverse of `select`; an `AbstractVector` of `Integer`, `Symbol`, `String`, or `Bool`, or a "drop" function of the form `(i, name) -> drop::Bool`; columns in the collection or for which the drop function returns `true` will ignored in the resulting `CSV.File`. Invalid values in `drop` are ignored.
  * `limit`: an `Integer` to indicate a limited number of rows to parse in a csv file; use in combination with `skipto` to read a specific, contiguous chunk within a file; note for large files when multiple threads are used for parsing, the `limit` argument may not result in an exact # of rows parsed; use `threaded=false` to ensure an exact limit if necessary
  * `buffer_in_memory`: a `Bool`, default `false`, which controls whether a `Cmd`, `IO`, or gzipped source will be read/decompressed in memory vs. using a temporary file.
  * `ntasks::Integer=Threads.nthreads()`: [not applicable to `CSV.Rows`] for multithreaded parsed files, this controls the number of tasks spawned to read a file in concurrent chunks; defaults to the # of threads Julia was started with (i.e. `JULIA_NUM_THREADS` environment variable or `julia -t N`); setting `ntasks=1` will avoid any calls to `Threads.@spawn` and just read the file serially on the main thread; a single thread will also be used for smaller files by default (< 5_000 cells)
  * `rows_to_check::Integer=$DEFAULT_ROWS_TO_CHECK`: [not applicable to `CSV.Rows`] a multithreaded parsed file will be split up into `ntasks` # of equal chunks; `rows_to_check` controls the # of rows are checked to ensure parsing correctly found valid rows; for certain files with very large quoted text fields, `lines_to_check` may need to be higher (10, 30, etc.) to ensure parsing correctly finds these rows
  * `source`: [only applicable for vector of inputs to `CSV.File`] a `Symbol`, `String`, or `Pair` of `Symbol` or `String` to `Vector`. As a single `Symbol` or `String`, provides the column name that will be added to the parsed columns, the values of the column will be the input "name" (usually file name) of the input from whence the value was parsed. As a `Pair`, the 2nd part of the pair should be a `Vector` of values matching the length of the # of inputs, where each value will be used instead of the input name for that inputs values in the auto-added column.

## Parsing options:

  * `missingstring`: either a `nothing`, `String`, or `Vector{String}` to use as sentinel values that will be parsed as `missing`; if `nothing` is passed, no sentinel/missing values will be parsed; by default, `missingstring=""`, which means only an empty field (two consecutive delimiters) is considered `missing`
  * `delim=','`: a `Char` or `String` that indicates how columns are delimited in a file; if no argument is provided, parsing will try to detect the most consistent delimiter on the first 10 rows of the file
  * `ignorerepeated::Bool=false`: whether repeated (consecutive/sequential) delimiters should be ignored while parsing; useful for fixed-width files with delimiter padding between cells
  * `quoted::Bool=true`: whether parsing should check for `quotechar` at the start/end of cells
  * `quotechar='"'`, `openquotechar`, `closequotechar`: a `Char` (or different start and end characters) that indicate a quoted field which may contain textual delimiters or newline characters
  * `escapechar='"'`: the `Char` used to escape quote characters in a quoted field
  * `dateformat::Union{String, Dates.DateFormat, Nothing, AbstractDict}`: a date format string to indicate how Date/DateTime columns are formatted for the entire file; if given as an `AbstractDict`, date format strings to indicate how the Date/DateTime columns corresponding to the keys are formatted. The Dict can map column index `Int`, or name `Symbol` or `String` to the format string for that column.
  * `decimal='.'`: a `Char` indicating how decimals are separated in floats, i.e. `3.14` uses `'.'`, or `3,14` uses a comma `','`
  * `truestrings`, `falsestrings`: `Vector{String}`s that indicate how `true` or `false` values are represented; by default `"true", "True", "TRUE", "T", "1"` are used to detect `true` and `"false", "False", "FALSE", "F", "0"` are used to detect `false`; note that columns with only `1` and `0` values will default to `Int64` column type unless explicitly requested to be `Bool` via `types` keyword argument

## Column Type Options:

  * `types`: a single `Type`, `AbstractVector` or `AbstractDict` of types, or a function of the form `(i, name) -> Union{T, Nothing}` to be used for column types; if a single `Type` is provided, _all_ columns will be parsed with that single type; an `AbstractDict` can map column index `Integer`, or name `Symbol` or `String` to type for a column, i.e. `Dict(1=>Float64)` will set the first column as a `Float64`, `Dict(:column1=>Float64)` will set the column named `column1` to `Float64` and, `Dict("column1"=>Float64)` will set the `column1` to `Float64`; if a `Vector` is provided, it must match the # of columns provided or detected in `header`. If a function is provided, it takes a column index and name as arguments, and should return the desired column type for the column, or `nothing` to signal the column's type should be detected while parsing.
  * `typemap::Dict{Type, Type}`: a mapping of a type that should be replaced in every instance with another type, i.e. `Dict(Float64=>String)` would change every detected `Float64` column to be parsed as `String`; only "standard" types are allowed to be mapped to another type, i.e. `Int64`, `Float64`, `Date`, `DateTime`, `Time`, and `Bool`. If a column of one of those types is "detected", it will be mapped to the specified type.
  * `pool::Union{Bool, Real, AbstractVector, AbstractDict, Function}=$DEFAULT_POOL`: [not supported by `CSV.Rows`] controls whether columns will be built as `PooledArray`; if `true`, all columns detected as `String` will be pooled; alternatively, the proportion of unique values below which `String` columns should be pooled (by default $DEFAULT_POOL, meaning that if the # of unique strings in a column is under $(DEFAULT_POOL * 100)%, it will be pooled); if an `AbstractVector`, each element should be `Bool` or `Real` and the # of elements should match the # of columns in the dataset; if an `AbstractDict`, a `Bool` or `Real` value can be provided for individual columns where the dict key is given as column index `Integer`, or column name as `Symbol` or `String`. If a function is provided, it should take a column index and name as 2 arguments, and return a `Bool`, `Float64`, or `nothing` for each column.
  * `downcast::Bool=false`: controls whether columns detected as `Int64` will be "downcast" to the smallest possible integer type like `Int8`, `Int16`, `Int32`, etc.
  * `stringtype=$DEFAULT_STRINGTYPE`: controls how detected string columns will ultimately be returned; default is `InlineString`, which stores string data in a fixed-size primitive type that helps avoid excessive heap memory usage; if a column has values longer than 32 bytes, it will default to `String`. If `String` is passed, all string columns will just be normal `String` values. If `PosLenString` is passed, string columns will be returned as `PosLenStringVector`, which is a special "lazy" `AbstractVector` that acts as a "view" into the original file data. This can lead to the most efficient parsing times, but note that the "view" nature of `PosLenStringVector` makes it read-only, so operations like `push!`, `append!`, or `setindex!` are not supported. It also keeps a reference to the entire input dataset source, so trying to modify or delete the underlying file, for example, may fail
  * `strict::Bool=false`: whether invalid values should throw a parsing error or be replaced with `missing`
  * `silencewarnings::Bool=false`: if `strict=false`, whether invalid value warnings should be silenced
  * `maxwarnings::Int=$DEFAULT_MAX_WARNINGS`: if more than `maxwarnings` number of warnings are printed while parsing, further warnings will be silenced by default; for multithreaded parsing, each parsing task will print up to `maxwarnings`
  * `debug::Bool=false`: passing `true` will result in many informational prints while a dataset is parsed; can be useful when reporting issues or figuring out what is going on internally while a dataset is parsed
  * `validate::Bool=true`: whether or not to validate that columns specified in the `types`, `dateformat` and `pool` keywords are actually found in the data. If `false` no validation is done, meaning no error will be thrown if `types`/`dateformat`/`pool` specify settings for columns not actually found in the data.

## Iteration options:

  * `reusebuffer=false`: [only supported by `CSV.Rows`] while iterating, whether a single row buffer should be allocated and reused on each iteration; only use if each row will be iterated once and not re-used (e.g. it's not safe to use this option if doing `collect(CSV.Rows(file))` because only current iterated row is "valid")
"""


text:
using Logging
using IterTools
using LinearAlgebra

include("../../phantoms/sperm_phantom.jl")
include("../snake_forward.jl")
include("../snake.jl")
include("../curve_utils.jl")

cwd = @__DIR__
savepath = normpath(joinpath(@__DIR__, "result"))
!isdir(savepath) && mkdir(savepath)

#object to store center and radius of circle in complex plane
struct circle
    center::Complex{Float64}
    radius::Float64
end

#return the parametric equation for a circle in the complex plane
function get_parametric_circle_equation(c::circle)
    return t -> c.radius*ℯ^(t*im) + c.center
end

H,W = 609, 609
detmin, detmax = -38.0, 38.0
grid = collect(detmin:0.1:detmax)
r(s) = 1.0

images, tracks = get_sperm_phantom(301,r,grid)

grid = collect(detmin:0.1:detmax)
bins = collect(detmin:0.125:detmax)
angles = collect(range(0,π,length=9))[1:end-1]

@info  "using " length(angles) " angles"
iterations = 1000
α = 0.01
β = 0.3
step_size = 0.2
r(s) = 1.0
tol = 0.001

num_points = 62
t = range(0, 2π, length=num_points)
c = circle((0.0-0.0*im), 15.0)
γ = get_parametric_circle_equation(c)
compl_circ = γ.(t)
outline_templ = cat(real(compl_circ), imag(compl_circ), dims=2)

frames2reconstruct = collect(1:10:301)
reconstructions = zeros(num_points,2,length(frames2reconstruct))

best_residual = ones(length(frames2reconstruct),4)

for step_size in collect(1.0:-0.1:0.01)
    for α in collect(1.0:-0.1:0.01)
        for β in collect(1.0:-0.1:0.01)
            for (iter, frame_nr) in Base.Iterators.reverse(enumerate(frames2reconstruct))
                @info iter frame_nr

                #Get projection
                @info "making forward projection for frame: " frame_nr
                outline, normals = get_outline(tracks[frame_nr], r)
                sinogram = parallel_forward(outline, angles, bins)

                #Add noise
                @info "adding gaussian noise at level 0.01"
                rho = 0.01
                e = randn(size(sinogram));
                e = rho*norm(sinogram)*e/norm(e);
                sinogram = sinogram + e;
                #plot(Gray.(sinogram./maximum(sinogram)))
                #cd(savepath)
                #savefig(@sprintf "sinogram_%d" frame_nr)

                @info "reconstructing"
                try
                    reconstruction, residuals = vedrana(outline_templ,angles,bins,sinogram, iterations, α, β, step_size, tol)

                if minimum(residuals) > 0.01 || minimum(residuals) > best_residual[iter,1]
                    @info "combination not better:" frame_nr step_size α β
                    continue
                else
                    best_residual[iter,:] =  [minimum(residuals), α, β, step_size]
                end

                reconstructions[:,:,iter] = reconstruction

                @info "plotting"
                heatmap(grid, grid, Gray.(images[:,:,frame_nr]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=true)
                plot!(reconstruction[:,1], reconstruction[:,2], aspect_ratio=:equal, linewidth=2, label=frame_nr)
                plot!(outline_templ[:,1], outline_templ[:,2], aspect_ratio=:equal, linewidth=2, label="template")
                cd(savepath)
                savefig(@sprintf "result_original_vedrana_%d" frame_nr)

                @info "plotting"
                plot(collect(1:1:length(residuals)), residuals)
                cd(savepath)
                savefig(@sprintf "residuals_%d" frame_nr)
                catch
                    @info "combination not error:" frame_nr step_size α β
                    continue
                end
            end
        end
    end
end

global ps = AbstractPlot[]
for (iter, frame_nr) in enumerate(frames2reconstruct)
    best_recon = reconstructions[:,:,iter]

    l = @layout [a b c]
    p = heatmap(grid, grid, Gray.(images[:,:,frame_nr]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    plot!(best_recon[:,1],best_recon[:,2], aspect_ratio=:equal, linewidth=5)
    push!(ps,p)
    if length(ps) == 3
        plot(ps[1], ps[2], ps[3], layout = l, size=(2000,600), linewidth=5)
        savefig(@sprintf "result_all_original_vedrana_%d_%d" length(angles) frame_nr)
        global ps = AbstractPlot[]
    end
end
cd(cwd)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
function assembleInferenceAlgorithm!(algo::InferenceAlgorithm)
    # Generate structures for fast lookup
    algo.interface_to_schedule_entry = interfaceToScheduleEntry(algo)
    algo.target_to_marginal_entry = targetToMarginalEntry(algo)

    for (id, pf) in algo.posterior_factorization
        pf.algorithm_id = algo.id
        assemblePosteriorFactor!(pf)
    end

    return algo
end

function assemblePosteriorFactor!(pf::PosteriorFactor)
    assembleSchedule!(pf)
    assembleInitialization!(pf)
    assembleMarginalTable!(pf)

    return pf
end

function assembleSchedule!(pf::PosteriorFactor)
    # Collect inbounds and assign message index per schedule entry
    for (msg_idx, schedule_entry) in enumerate(pf.schedule)
        schedule_entry.initialize = false # Preset initialization to false
        schedule_entry.inbounds = collectInbounds(schedule_entry, schedule_entry.message_update_rule)
        schedule_entry.schedule_index = msg_idx
    end

    return pf
end

function assembleInitialization!(pf::PosteriorFactor)
    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry

    # Collect outbound types from schedule
    outbound_types = Dict{Interface, Type}()
    for entry in pf.schedule
        outbound_types[entry.interface] = outboundType(entry.message_update_rule)
    end

    # Find breaker types and dimensions
    pf_update_clamp_flag = false # Flag that tracks whether the update of a clamped variable is required
    pf_initialize_flag = false # Indicate need for an initialization block
    for entry in pf.schedule
        partner = ultimatePartner(entry.interface)

        # Assemble breakers
        if entry.interface in pf.breaker_interfaces
            (_, dims) = breakerParameters(entry.interface)
            assembleBreaker!(entry, family(outbound_types[entry.interface]), dims)
            pf_initialize_flag = true # Signifies the need for an initialization block
        elseif isClamped(partner)
            assembleBreaker!(entry, family(outbound_types[entry.interface]), size(partner.node.value))
            pf_initialize_flag = true
            pf_update_clamp_flag = true # Signifies the need for creating a custom `step!` function for optimizing clamped variables
        end
    end

    pf.optimize = pf_update_clamp_flag
    pf.initialize = pf_initialize_flag

    return pf
end

function assembleMarginalTable!(pf::PosteriorFactor)
    interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry

    for entry in pf.marginal_table
        if entry.marginal_update_rule == Nothing
            iface = entry.interfaces[1]
            inbounds = [interface_to_schedule_entry[iface]]
        elseif entry.marginal_update_rule == Product
            iface1 = entry.interfaces[1]
            iface2 = entry.interfaces[2]
            inbounds = [interface_to_schedule_entry[iface1],
                        interface_to_schedule_entry[iface2]]
        else
            inbounds = collectInbounds(entry)
        end

        entry.marginal_id = entry.target.id
        entry.inbounds = inbounds
    end

    return pf
end

function assembleBreaker!(breaker_entry::ScheduleEntry, family::Type, dimensionality::Tuple)
    breaker_entry.initialize = true
    breaker_entry.family = family
    breaker_entry.dimensionality = dimensionality

    return breaker_entry
end

"""
Depending on the origin of the Clamp node message, contruct a message or marginal inbound
"""
function assembleClamp!(inbound::Clamp, dist_or_msg::Type)
    inbound.dist_or_msg = dist_or_msg
    if inbound in keys(ForneyLab.current_graph.placeholders)
        # Message comes from data buffer
        (buffer, idx) = ForneyLab.current_graph.placeholders[inbound]
        inbound.buffer_id = buffer
        inbound.buffer_index = idx # Can also be 0, in which case index is ignored
    end

    return inbound
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 175c73f6-f2e5-11ea-3458-eddf1465c5c4
using Pkg, DrWatson

# ╔═╡ 175d464e-f2e5-11ea-318a-dd4ad9573a95
using StatisticalRethinking

# ╔═╡ 0b40c1ca-f2e4-11ea-1658-0d38b13854c0
md"## Clip-03-02-05t.jl"

# ╔═╡ 175ca6a0-f2e5-11ea-066c-49e37107aea2
@quickactivate "StatisticalRethinkingTuring"

# ╔═╡ 176c899a-f2e5-11ea-3f17-59b92f0a2483
md"### snippet 3.2"

# ╔═╡ 176d23ca-f2e5-11ea-3a5c-93353273f6eb
begin
	p_grid = range(0, stop=1, length=1000)
	prior = ones(length(p_grid))
	likelihood = pdf.(Binomial.(9, p_grid), 6)
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
end

# ╔═╡ 177885bc-f2e5-11ea-0d11-f9ae48c08d8f
md"### snippet 3.3"

# ╔═╡ 17798f98-f2e5-11ea-299e-f3424848225c
md"##### Draw 10000 samples from this posterior distribution."

# ╔═╡ 1784a202-f2e5-11ea-051b-09f8abb69875
begin
	N = 10000
	samples = sample(p_grid, Weights(posterior), N)
end;

# ╔═╡ 178bd8ba-f2e5-11ea-21dd-61c16076a938
md"##### Create an MCMCChains.Chains object."

# ╔═╡ 178d0d36-f2e5-11ea-2593-bf4c46edb720
chn = MCMCChains.Chains(reshape(samples, N, 1, 1), [:p]);

# ╔═╡ df229f44-7dd0-11eb-2ebe-9b12eff28aa1
typeof(chn)

# ╔═╡ 179c95ce-f2e5-11ea-2358-714343a64afa
md"##### Describe the chain."

# ╔═╡ 17a5374c-f2e5-11ea-1bfe-3fb1cce48d48
Text(sprint(show, "text/plain", summarize(chn)))

# ╔═╡ 17a70c16-f2e5-11ea-2e3f-2dc8006c1ee7
md"##### Plot the chain."

# ╔═╡ 17aece92-f2e5-11ea-311a-1b4591b1172a
p1 = plot(chn; seriestype=:traceplot)

# ╔═╡ 17bdf552-f2e5-11ea-28e5-1f5f22b43146
md"### snippet 3.4"

# ╔═╡ 17c6af64-f2e5-11ea-0b28-c55bf548b3ea
p2 = plot(chn; seriestype=:density)

# ╔═╡ 17cf1094-f2e5-11ea-1a80-858534e81348
md"### snippet 3.5"

# ╔═╡ 17d80aa0-f2e5-11ea-1c8f-01770152e142
md"##### Compare with analytical (conjugate) solution."

# ╔═╡ 17e08cca-f2e5-11ea-31a6-8931f98059ba
begin
	w = 6
	n = 9
	x = 0:0.01:1
	plot( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
	density!(samples, lab="Sample density")
end

# ╔═╡ 17f4a430-f2e5-11ea-048a-73ae33e192ed
begin
	density(samples, lab="Sample2 density")
	vline!(hpdi(samples), lab="hpdi samples2")
	vline!(quantile(samples, [0.25, 0.75]), lab="quantiles [0.25, 0.75]")
end

# ╔═╡ 17fecdb6-f2e5-11ea-2c5b-7bc56e8b1360
md"## End of clip-03-02-05t.jl"

# ╔═╡ Cell order:
# ╟─0b40c1ca-f2e4-11ea-1658-0d38b13854c0
# ╠═175c73f6-f2e5-11ea-3458-eddf1465c5c4
# ╠═175ca6a0-f2e5-11ea-066c-49e37107aea2
# ╠═175d464e-f2e5-11ea-318a-dd4ad9573a95
# ╟─176c899a-f2e5-11ea-3f17-59b92f0a2483
# ╠═176d23ca-f2e5-11ea-3a5c-93353273f6eb
# ╟─177885bc-f2e5-11ea-0d11-f9ae48c08d8f
# ╟─17798f98-f2e5-11ea-299e-f3424848225c
# ╠═1784a202-f2e5-11ea-051b-09f8abb69875
# ╟─178bd8ba-f2e5-11ea-21dd-61c16076a938
# ╠═178d0d36-f2e5-11ea-2593-bf4c46edb720
# ╠═df229f44-7dd0-11eb-2ebe-9b12eff28aa1
# ╟─179c95ce-f2e5-11ea-2358-714343a64afa
# ╠═17a5374c-f2e5-11ea-1bfe-3fb1cce48d48
# ╟─17a70c16-f2e5-11ea-2e3f-2dc8006c1ee7
# ╠═17aece92-f2e5-11ea-311a-1b4591b1172a
# ╟─17bdf552-f2e5-11ea-28e5-1f5f22b43146
# ╠═17c6af64-f2e5-11ea-0b28-c55bf548b3ea
# ╟─17cf1094-f2e5-11ea-1a80-858534e81348
# ╟─17d80aa0-f2e5-11ea-1c8f-01770152e142
# ╠═17e08cca-f2e5-11ea-31a6-8931f98059ba
# ╠═17f4a430-f2e5-11ea-048a-73ae33e192ed
# ╟─17fecdb6-f2e5-11ea-2c5b-7bc56e8b1360


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

using FEMQuad, Test
using FEMQuad: get_rule, integrate_1d, integrate_2d, integrate_3d

f(i) = x -> sum(sum(x)^j for j=0:i)

@testset "Gauss-Legendre quadratures in one dimension" begin
    # Gaussian quadrature using N points can provide the exact integral if
    # polynomial degree of f is 2N-1 or less => N = (deg(f)+1)/2

    rules = (:GLSEG1, :GLSEG2, :GLSEG3, :GLSEG4, :GLSEG5)
    r(i::Int) = get_rule(i, rules...)

    @test isapprox(integrate_1d(f(0), r(0)), 2.0)
    @test isapprox(integrate_1d(f(1), r(1)), 2.0)
    @test isapprox(integrate_1d(f(2), r(2)), 8/3)
    @test isapprox(integrate_1d(f(3), r(3)), 8/3)
    @test isapprox(integrate_1d(f(4), r(4)), 46/15)
    @test isapprox(integrate_1d(f(5), r(5)), 46/15)
    @test isapprox(integrate_1d(f(6), r(6)), 352/105)
    @test isapprox(integrate_1d(f(7), r(7)), 352/105)
    @test isapprox(integrate_1d(f(8), r(8)), 1126/315)
    @test isapprox(integrate_1d(f(9), r(9)), 1126/315)
end

@testset "Gauss-Legendre quadratures in quadrilaterals" begin
    rules = (:GLQUAD1, :GLQUAD4, :GLQUAD9, :GLQUAD16, :GLQUAD25)
    r(i::Int) = get_rule(i, rules...)

    @test isapprox(integrate_2d(f(0), r(0)), 4.0)
    @test isapprox(integrate_2d(f(1), r(1)), 4.0)
    @test isapprox(integrate_2d(f(2), r(2)), 20/3)
    @test isapprox(integrate_2d(f(3), r(3)), 20/3)
    @test isapprox(integrate_2d(f(4), r(4)), 164/15)
    @test isapprox(integrate_2d(f(5), r(5)), 164/15)
    @test isapprox(integrate_2d(f(6), r(6)), 2108/105)
    @test isapprox(integrate_2d(f(7), r(7)), 2108/105)
    @test isapprox(integrate_2d(f(8), r(8)), 13492/315)
    @test isapprox(integrate_2d(f(9), r(9)), 13492/315)
end

@testset "Gauss-Legendre quadratures in hexahedrons" begin
    rules = (:GLHEX1, :GLHEX8, :GLHEX27, :GLHEX64, :GLHEX125)
    r(i::Int) = get_rule(i, rules...)

    @test isapprox(integrate_3d(f(0), r(0)), 8.0)
    @test isapprox(integrate_3d(f(1), r(1)), 8.0)
    @test isapprox(integrate_3d(f(2), r(2)), 16.0)
    @test isapprox(integrate_3d(f(3), r(3)), 16.0)
    @test isapprox(integrate_3d(f(4), r(4)), 184/5)
    @test isapprox(integrate_3d(f(5), r(5)), 184/5)
    @test isapprox(integrate_3d(f(6), r(6)), 12064/105)
    @test isapprox(integrate_3d(f(7), r(7)), 12064/105)
    @test isapprox(integrate_3d(f(8), r(8)), 9928/21)
    @test isapprox(integrate_3d(f(9), r(9)), 9928/21)
end


text:
module WeakRefStrings

export WeakRefString, WeakRefStringArray, StringArray, StringVector

########################################################################
# WeakRefString
########################################################################

"""
A custom "weakref" string type that only points to external string data.
Allows for the creation of a "string" instance without copying data,
which allows for more efficient string parsing/movement in certain data processing tasks.

**Please note that no original reference is kept to the parent string/memory, so `WeakRefString` becomes unsafe
once the parent object goes out of scope (i.e. loses a reference to it)**

Internally, a `WeakRefString{T}` holds:

  * `ptr::Ptr{T}`: a pointer to the string data (code unit size is parameterized on `T`)
  * `len::Int`: the number of code units in the string data

See also [`StringArray`](@ref)
"""
struct WeakRefString{T} <: AbstractString
    ptr::Ptr{T}
    len::Int # of code units
end

WeakRefString(ptr::Ptr{T}, len) where {T} = WeakRefString(ptr, Int(len))
WeakRefString(t::Tuple{Ptr{T}, Int}) where {T} = WeakRefString(t...)
WeakRefString(a::Vector{T}) where T <: Union{UInt8,UInt16,UInt32} = WeakRefString(pointer(a), length(a))

const NULLSTRING = WeakRefString(Ptr{UInt8}(0), 0)
const NULLSTRING16 = WeakRefString(Ptr{UInt16}(0), 0)
const NULLSTRING32 = WeakRefString(Ptr{UInt32}(0), 0)
Base.zero(::Type{WeakRefString{T}}) where {T} = WeakRefString(Ptr{T}(0), 0)

import Base: ==
function ==(x::WeakRefString{T}, y::WeakRefString{T}) where {T}
    x.len == y.len && (x.ptr == y.ptr || ccall(:memcmp, Cint, (Ptr{T}, Ptr{T}, Csize_t),
                                           x.ptr, y.ptr, x.len) == 0)
end
function ==(x::String, y::WeakRefString{T}) where {T}
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{T}, Ptr{T}, Csize_t),
                                 pointer(x), y.ptr, y.len) == 0
end
==(y::WeakRefString, x::String) = x == y

function Base.cmp(a::WeakRefString{T}, b::WeakRefString{T}) where T
    al, bl = a.len, b.len
    c = ccall(:memcmp, Int32, (Ptr{T}, Ptr{T}, Csize_t),
              a.ptr, b.ptr, min(al,bl))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(al,bl)
end

function Base.hash(s::WeakRefString{T}, h::UInt) where {T}
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{T}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

function Base.show(io::IO, x::WeakRefString{T}) where {T}
    print(io, '"')
    print(io, string(x))
    print(io, '"')
    return
end
Base.print(io::IO, s::WeakRefString) = print(io, string(s))
Base.textwidth(s::WeakRefString) = textwidth(string(s))

chompnull(x::WeakRefString{T}) where {T} = unsafe_load(x.ptr, x.len) == T(0) ? x.len - 1 : x.len

Base.string(x::WeakRefString) = x == NULLSTRING ? "" : unsafe_string(x.ptr, x.len)
Base.string(x::WeakRefString{UInt16}) = x == NULLSTRING16 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))
Base.string(x::WeakRefString{UInt32}) = x == NULLSTRING32 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))

Base.convert(::Type{WeakRefString{UInt8}}, x::String) = WeakRefString(pointer(x), sizeof(x))
Base.convert(::Type{String}, x::WeakRefString) = convert(String, string(x))
Base.String(x::WeakRefString) = string(x)
Base.Symbol(x::WeakRefString{UInt8}) = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), x.ptr, x.len)

Base.pointer(s::WeakRefString) = s.ptr
Base.pointer(s::WeakRefString, i::Integer) = s.ptr + i - 1

Base.ncodeunits(s::WeakRefString{T}) where {T} = s.len
Base.codeunit(s::WeakRefString{T}) where {T} = T

@inline function Base.codeunit(s::WeakRefString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(pointer(s, i))
end

Base.thisind(s::WeakRefString, i::Int) = Base._thisind_str(s, i)
Base.nextind(s::WeakRefString, i::Int) = Base._nextind_str(s, i)
Base.isvalid(s::WeakRefString, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i

Base.@propagate_inbounds function Base.iterate(s::WeakRefString, i::Int=firstindex(s))
    i > ncodeunits(s) && return nothing
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return iterate_continued(s, i, u)
end

if isdefined(Base, :next_continued)
    const iterate_continued = Base.next_continued
else
    const iterate_continued = Base.iterate_continued
end

function iterate_continued(s::WeakRefString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

########################################################################
# StringArray
########################################################################

const STR = Union{Missing, <:AbstractString}

"""
`StringArray{T,N}`

Efficient storage for N dimensional array of strings.

`StringArray` stores underlying string data for all elements of the array
in a single contiguous buffer. It maintains offset and length for each
element.

`T` can be `String`, `WeakRefString`, `EscapedString`, `Union{Missing, String}`,
`Union{Missing, WeakRefString}`, or `Union{Missing, EscapedString}`.
`getindex` will return this type except `EscapedString` returns the
unescaped `String`.

You can use `convert(StringArray{U}, ::StringArray{T})` to change the
element type (e.g. to `WeakRefString` for efficiency) without copying
the data.

# Example construction

```
julia> sa = StringArray(["x", "y"]) # from Array{String}
2-element WeakRefStrings.StringArray{String,1}:
 "x"
 "y"

julia> sa = StringArray{WeakRefString}(["x", "y"])
2-element WeakRefStrings.StringArray{WeakRefStrings.WeakRefString,1}:
 "x"
 "y"

julia> sa = StringArray{Union{Missing, String}}(["x", "y"]) # with Missing
2-element WeakRefStrings.StringArray{Union{Missing, String},1}:
 "x"
 "y"

julia> sa = StringArray{Union{Missing, String}}(2,2) # undef
2×2 WeakRefStrings.StringArray{Union{Missing, String},2}:
 #undef  #undef
 #undef  #undef
```
"""
struct StringArray{T, N} <: AbstractArray{T, N}
    buffer::Vector{UInt8}
    offsets::Array{UInt64, N}
    lengths::Array{UInt32, N}
end

const StringVector{T} = StringArray{T, 1}

const UNDEF_OFFSET = typemax(UInt64)
const MISSING_OFFSET = typemax(UInt64)-1

Base.size(a::StringArray) = size(a.offsets)
Base.IndexStyle(::Type{<:StringVector}) = IndexLinear()

# no-copy convert between eltypes
"""
`convert(StringArray{U}, A::StringArray{T})`

convert `A` to StringArray of another element type (`U`) without
copying the underlying data.
"""
function Base.convert(::Type{<:StringArray{T}}, x::StringArray{<:STR,N}) where {T, N}
    StringArray{T, ndims(x)}(x.buffer, x.offsets, x.lengths)
end

function StringArray{T, N}(::UndefInitializer, dims::Tuple{Vararg{Integer}}) where {T,N}
    StringArray{T,N}(similar(Vector{UInt8}, 0), fill(UNDEF_OFFSET, dims), fill(zero(UInt32), dims))
end

function StringArray{T}(::UndefInitializer, dims::Tuple{Vararg{Integer}}) where {T}
    StringArray{T,length(dims)}(undef, dims)
end

function StringArray(::UndefInitializer, dims::Tuple{Vararg{Integer}})
    StringArray{String}(undef, dims)
end

StringArray{T, N}(::UndefInitializer, dims::Vararg{Integer,N}) where {T,N} = StringArray{T,N}(undef, dims)
StringArray{T}(::UndefInitializer, dims::Integer...) where {T} = StringArray{T,length(dims)}(undef, dims)
StringArray(::UndefInitializer, dims::Integer...) = StringArray{String,length(dims)}(undef, dims)

function Base.convert(::Type{<:StringArray{T}}, arr::AbstractArray{<:STR, N}) where {T,N}
    s = StringArray{T, N}(undef, size(arr))
    @inbounds for i in eachindex(arr)
        if _isassigned(arr, i)
            s[i] = arr[i]
        else
            s.offsets[i] = UNDEF_OFFSET
            s.lengths[i] = 0
        end
    end
    s
end
Base.convert(::Type{StringArray}, arr::AbstractArray{T}) where {T<:STR} = StringArray{T}(arr)
Base.convert(::Type{StringArray{T, N} where T}, arr::AbstractArray{S}) where {S<:STR, N} = StringVector{S}(arr)
StringVector{T}() where {T} = StringVector{T}(UInt8[], UInt64[], UInt32[])
StringVector() = StringVector{String}()
StringVector{T}(::UndefInitializer, len::Int) where {T} = StringArray{T}(undef, len)
StringVector(::UndefInitializer, len::Int) = StringArray{String}(undef, len)
# special constructor where a full byte buffer is provided and offsets/lens will be filled in later
StringVector{T}(buffer::Vector{UInt8}, len::Int) where {T} = StringVector{T}(buffer, fill(UNDEF_OFFSET, len), fill(zero(UInt32), len))

(T::Type{<:StringArray})(arr::AbstractArray{<:STR}) = convert(T, arr)

_isassigned(arr, i...) = isassigned(arr, i...)
_isassigned(arr, i::CartesianIndex) = isassigned(arr, i.I...)
@inline Base.@propagate_inbounds function Base.getindex(a::StringArray{T}, i::Integer...) where T
    offset = a.offsets[i...]
    if offset == UNDEF_OFFSET
        throw(UndefRefError())
    end

    if Missing <: T && offset === MISSING_OFFSET
        return missing
    end

    convert(T, WeakRefString(pointer(a.buffer) + offset, a.lengths[i...]))
end

function Base.similar(a::StringArray, T::Type{<:STR}, dims::Tuple{Vararg{Int64, N}}) where N
    StringArray{T, N}(undef, dims)
end

function Base.empty!(a::StringVector)
    empty!(a.buffer)
    empty!(a.offsets)
    empty!(a.lengths)
    a
end

Base.copy(a::StringArray{T, N}) where {T,N} = StringArray{T, N}(copy(a.buffer), copy(a.offsets), copy(a.lengths))

@inline function Base.setindex!(arr::StringArray, val::WeakRefString, idx::Integer...)
    p = pointer(arr.buffer)
    if val.ptr <= p + sizeof(arr.buffer)-1 && val.ptr >= p
        # this WeakRefString points to data entirely within arr's buffer
        # don't add anything to the buffer in this case.
        # this optimization helps `permute!`
        arr.offsets[idx...] = val.ptr - p
        arr.lengths[idx...] = val.len
        val
    else
        _setindex!(arr, val, idx...)
    end
end

@inline function Base.setindex!(arr::StringArray, val::STR, idx::Integer...)
    _setindex!(arr, val, idx...)
end

@inline function Base.setindex!(arr::StringArray, val::STR, idx::Integer)
    _setindex!(arr, val, idx)
end

function _setindex!(arr::StringArray, val::AbstractString, idx...)
    buffer = arr.buffer
    l = length(arr.buffer)
    resize!(buffer, l + sizeof(val))
    unsafe_copyto!(pointer(buffer, l+1), pointer(val,1), sizeof(val))
    arr.lengths[idx...] = sizeof(val)
    arr.offsets[idx...] = l
    val
end

function _setindex!(arr::StringArray, val::AbstractString, idx)
    buffer = arr.buffer
    l = length(arr.buffer)
    resize!(buffer, l + sizeof(val))
    unsafe_copyto!(pointer(buffer, l+1), pointer(val,1), sizeof(val))
    arr.lengths[idx] = sizeof(val)
    arr.offsets[idx] = l
    val
end

function _setindex!(arr::StringArray{Union{T, Missing}, N}, val::Missing, idx) where {T, N}
    arr.lengths[idx] = 0
    arr.offsets[idx] = MISSING_OFFSET
    val
end

function _setindex!(arr::StringArray{Union{T, Missing}, N}, val::Missing, idx...) where {T, N}
    arr.lengths[idx...] = 0
    arr.offsets[idx...] = MISSING_OFFSET
    val
end

function Base.resize!(arr::StringVector, len)
    l = length(arr)
    resize!(arr.offsets, len)
    resize!(arr.lengths, len)
    if l < len
        arr.offsets[l+1:len] .= UNDEF_OFFSET # undef
        arr.lengths[l+1:len] .= 0
    end
    arr
end

function Base.push!(arr::StringVector, val::AbstractString)
    l = length(arr.buffer)
    resize!(arr.buffer, l + sizeof(val))
    unsafe_copyto!(pointer(arr.buffer, l + 1), pointer(val,1), sizeof(val))
    push!(arr.offsets, l)
    push!(arr.lengths, sizeof(val))
    arr
end

function Base.push!(arr::StringVector{Union{T, Missing}}, val::Missing) where {T}
    push!(arr.offsets, MISSING_OFFSET)
    push!(arr.lengths, 0)
    arr
end

function Base.deleteat!(arr::StringVector, idx)
    deleteat!(arr.lengths, idx)
    deleteat!(arr.offsets, idx)
    arr
end

function Base.insert!(arr::StringVector, idx::Integer, item::AbstractString)
    l = length(arr.buffer)
    resize!(arr.buffer, l + sizeof(item))
    unsafe_copyto!(pointer(arr.buffer, l + 1), pointer(item), sizeof(item))
    insert!(arr.offsets, idx, l)
    insert!(arr.lengths, idx, sizeof(item))
    arr
end

function Base.insert!(arr::StringVector{Union{T, Missing}}, idx::Integer, item::Missing) where {T}
    insert!(arr.offsets, idx, MISSING_OFFSET)
    insert!(arr.lengths, idx, 0)
    arr
end

function Base.permute!(arr::StringArray{String}, p::AbstractVector)
    permute!(convert(StringArray{WeakRefString{UInt8}}, arr), p)
    arr
end

function Base.sortperm(arr::StringArray{String})
    sortperm(convert(StringArray{WeakRefString{UInt8}}, arr))
end

function Base.sort!(arr::StringArray{String})
    permute!(arr, sortperm(arr))
    arr
end

function Base.vcat(a::StringVector{T}, b::StringVector{T}) where T
    StringVector{T}(vcat(a.buffer, b.buffer), vcat(a.offsets, b.offsets .+ length(a.buffer)), vcat(a.lengths, b.lengths))
end

function Base.append!(a::StringVector{T}, b::StringVector) where T
    append!(a.offsets, b.offsets .+ length(a.buffer))
    append!(a.buffer, b.buffer)
    append!(a.lengths, b.lengths)
    a
end

function Base.append!(a::StringVector{T}, b::AbstractVector) where T
    for x in b
        push!(a, x)
    end
    a
end

function _growat!(a::StringVector, i, len)
    Base._growat!(a.offsets, i, len)
    Base._growat!(a.lengths, i, len)
    return
end

function _deleteat!(a::StringVector, i, len)
    Base._deleteat!(a.offsets, i, len)
    Base._deleteat!(a.lengths, i, len)
    return
end

const _default_splice = []

function Base.splice!(a::StringVector, i::Integer, ins=_default_splice)
    v = a[i]
    m = length(ins)
    if m == 0
        deleteat!(a, i)
    elseif m == 1
        a[i] = ins[1]
    else
        _growat!(a, i, m-1)
        k = 1
        for x in ins
            a[i+k-1] = x
            k += 1
        end
    end
    return v
end

function Base.splice!(a::StringVector, r::UnitRange{<:Integer}, ins=_default_splice)
    v = a[r]
    m = length(ins)
    if m == 0
        deleteat!(a, r)
        return v
    end

    n = length(a)
    f = first(r)
    l = last(r)
    d = length(r)

    if m < d
        delta = d - m
        _deleteat!(a, (f - 1 < n - l) ? f : (l - delta + 1), delta)
    elseif m > d
        _growat!(a, (f - 1 < n - l) ? f : (l + 1), m - d)
    end

    k = 1
    for x in ins
        a[f+k-1] = x
        k += 1
    end
    return v
end

function _growbeg!(a::StringVector, n)
    Base._growbeg!(a.offsets, n)
    Base._growbeg!(a.lengths, n)
    return
end

function Base.prepend!(a::StringVector, items::AbstractVector)
    itemindices = eachindex(items)
    n = length(itemindices)
    _growbeg!(a, n)
    if a === items
        copyto!(a, 1, items, n+1, n)
    else
        copyto!(a, 1, items, first(itemindices), n)
    end
    return a
end

Base.prepend!(a::StringVector, iter) = _prepend!(a, Base.IteratorSize(iter), iter)
Base.pushfirst!(a::StringVector, iter...) = prepend!(a, iter)

function _prepend!(a, ::Union{Base.HasLength,Base.HasShape}, iter)
    n = length(iter)
    _growbeg!(a, n)
    i = 0
    for item in iter
        @inbounds a[i += 1] = item
    end
    a
end

function _prepend!(a, ::Base.IteratorSize, iter)
    n = 0
    for item in iter
        n += 1
        pushfirst!(a, item)
    end
    reverse!(a, 1, n)
    a
end

function Base.pop!(a::StringVector)
    if isempty(a)
        throw(ArgumentError("array must be non-empty"))
    end
    item = a[end]
    deleteat!(a, length(a))
    return item
end

function Base.pushfirst!(a::StringVector, item)
    _growbeg!(a, 1)
    a[1] = item
    return a
end

function Base.popfirst!(a::StringVector)
    if isempty(a)
        throw(ArgumentError("array must be non-empty"))
    end
    item = a[1]
    deleteat!(a, 1)
    return item
end

end # module


text:
struct MultivariateNormal <: Distribution{Vector{Float64}} end

"""
    mvnormal(mu::AbstractVector{T}, cov::AbstractMatrix{U}} where {T<:Real,U<:Real}

Samples a `Vector{Float64}` value from a multivariate normal distribution.
"""
const mvnormal = MultivariateNormal()

function logpdf(::MultivariateNormal, x::AbstractVector{T}, mu::AbstractVector{U},
                cov::AbstractMatrix{V}) where {T,U,V}
    dist = Distributions.MvNormal(mu, cov)
    Distributions.logpdf(dist, x)
end

function logpdf_grad(::MultivariateNormal, x::AbstractVector{T}, mu::AbstractVector{U},
                cov::AbstractMatrix{V}) where {T,U,V}
    dist = Distributions.MvNormal(mu, cov)
    inv_cov = Distributions.invcov(dist)

    x_deriv = Distributions.gradlogpdf(dist, x)
    mu_deriv = -x_deriv
    cov_deriv = -0.5 * (inv_cov - (mu_deriv * transpose(mu_deriv)))

    (x_deriv, mu_deriv, cov_deriv)
end

function random(::MultivariateNormal, mu::AbstractVector{U},
                cov::AbstractMatrix{V}) where {T,U,V}
    rand(Distributions.MvNormal(mu, cov))
end

(::MultivariateNormal)(mu, cov) = random(MultivariateNormal(), mu, cov)

has_output_grad(::MultivariateNormal) = true
has_argument_grads(::MultivariateNormal) = (true, true)

export mvnormal

max_stars_repo_name: cameronraysmith/Gen

max_issues_repo_name: cameronraysmith/Gen

max_forks_repo_name: cameronraysmith/Gen


text:
using Gen
import Random

include("static_model.jl")
include("dataset.jl")

@gen (static) function slope_proposal(trace)
    slope = trace[:slope]
    @trace(normal(slope, 0.5), :slope)
end

@gen (static) function intercept_proposal(trace)
    intercept = trace[:intercept]
    @trace(normal(intercept, 0.5), :intercept)
end

@gen (static) function inlier_std_proposal(trace)
    log_inlier_std = trace[:log_inlier_std]
    @trace(normal(log_inlier_std, 0.5), :log_inlier_std)
end

@gen (static) function outlier_std_proposal(trace)
    log_outlier_std = trace[:log_outlier_std]
    @trace(normal(log_outlier_std, 0.5), :log_outlier_std)
end

@gen (static) function flip_z(z::Bool)
    @trace(bernoulli(z ? 0.0 : 1.0), :z)
end

@gen (static) function is_outlier_proposal(trace, i::Int)
    prev_z = trace[:data => i => :z]
    @trace(bernoulli(prev_z ? 0.0 : 1.0), :data => i => :z)
end

@gen (static) function is_outlier_proposal(trace, i::Int)
    prev_z = trace[:data => i => :z]
    @trace(bernoulli(prev_z ? 0.0 : 1.0), :data => i => :z)
end

Gen.load_generated_functions()

function do_inference(xs, ys, num_iters)
    observations = choicemap()
    for (i, y) in enumerate(ys)
        observations[:data => i => :y] = y
    end

    # initial trace
    (trace, _) = generate(model, (xs,), observations)

    scores = Vector{Float64}(undef, num_iters)
    for i=1:num_iters

        # steps on the parameters
        for j=1:5
            (trace, _) = metropolis_hastings(trace, slope_proposal, ())
            (trace, _) = metropolis_hastings(trace, intercept_proposal, ())
            (trace, _) = metropolis_hastings(trace, inlier_std_proposal, ())
            (trace, _) = metropolis_hastings(trace, outlier_std_proposal, ())
        end

        # step on the outliers
        for j=1:length(xs)
            (trace, _) = metropolis_hastings(trace, is_outlier_proposal, (j,))
        end

        score = get_score(trace)
        scores[i] = score

        # print
        slope = trace[:slope]
        intercept = trace[:intercept]
        inlier_std = exp(trace[:log_inlier_std])
        outlier_std = exp(trace[:log_outlier_std])
        println("score: $score, slope: $slope, intercept: $intercept, inlier_std: $inlier_std, outlier_std: $outlier_std")
    end
    return scores
end

(xs, ys) = make_data_set(200)
do_inference(xs, ys, 10)
@time do_inference(xs, ys, 50)

# handcoded version

mutable struct State
    xs::Vector{Float64}
    log_inlier_std::Float64
    log_outlier_std::Float64
    slope::Float64
    intercept::Float64
    is_outlier::Vector{Bool}
    ys::Vector{Float64}
end

function init_state(xs, ys)
    log_inlier_std = random(normal, 0, 2)
    log_outlier_std = random(normal, 0, 2)
    slope = random(normal, 0, 2)
    intercept = random(normal, 0, 2)
    is_outlier = Vector{Bool}(undef, length(xs))
    for i=1:length(xs)
        is_outlier[i] = (rand() < 0.5)
    end
    State(xs, log_inlier_std, log_outlier_std, slope, intercept, is_outlier, ys)
end

function log_likelihood(x, y, is_outlier, log_inlier_std, log_outlier_std, slope, intercept)
    std = is_outlier ? exp(log_outlier_std) : exp(log_inlier_std)
    logpdf(normal, y, slope * x + intercept, std)
end

function log_likelihood(state)
    ll = 0.
    for (x, is_outlier, y) in zip(state.xs, state.is_outlier, state.ys)
        ll += log_likelihood(x, y, is_outlier,
                    state.log_inlier_std, state.log_outlier_std, state.slope, state.intercept)
    end
    ll
end

macro parameter_mh_move!(state, addr)
    quote
        prev_ll = log_likelihood($(esc(state)))
        prev_val = $(esc(state)).$addr
        new_val = normal(prev_val, 0.5)
        prev_prior = logpdf(normal, prev_val, 0, 2)
        new_prior = logpdf(normal, new_val, 0, 2)
        $(esc(state)).$addr = new_val
        new_ll = log_likelihood($(esc(state)))
        if log(rand()) >= new_ll - prev_ll + new_prior - prev_prior
            $(esc(state)).$addr = prev_val # reject
        end
    end
end

function is_outlier_mh_move!(state, i::Int)
    prev_val = state.is_outlier[i]
    new_val = !prev_val
    x = state.xs[i]
    y = state.ys[i]
    prev_ll = log_likelihood(x, y, prev_val,
        state.log_inlier_std, state.log_outlier_std, state.slope, state.intercept)
    new_ll = log_likelihood(x, y, new_val,
        state.log_inlier_std, state.log_outlier_std, state.slope, state.intercept)
    if log(rand()) < new_ll - prev_ll
        state.is_outlier[i] = new_val
    end
end

function handcoded_inference(xs, ys, iters)
    state = init_state(xs, ys)
    for iter=1:iters
        for j=1:5
            @parameter_mh_move!(state, slope)
            @parameter_mh_move!(state, intercept)
            @parameter_mh_move!(state, log_inlier_std)
            @parameter_mh_move!(state, log_outlier_std)
        end
        for j=1:length(xs)
            is_outlier_mh_move!(state, j)
        end
    end
    println("slope: $(state.slope)")
    println("intercept: $(state.intercept)")
end

handcoded_inference(xs, ys, 50)
@time handcoded_inference(xs, ys, 50)

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using MPI, Test
include(joinpath("..", "..", "testhelpers.jl"))

@testset "DGMethods" begin
    runmpi(joinpath(@__DIR__, "courant.jl"), ntasks = 2)
    runmpi(joinpath(@__DIR__, "grad_test.jl"))
    runmpi(joinpath(@__DIR__, "grad_test_sphere.jl"))
    runmpi(joinpath(@__DIR__, "horizontal_integral_test.jl"))
    runmpi(joinpath(@__DIR__, "integral_test.jl"))
    runmpi(joinpath(@__DIR__, "integral_test_sphere.jl"))
    runmpi(joinpath(@__DIR__, "vars_test.jl"))
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
include("roadway/main.jl")
include("definitions.jl")
include("vehicles/main.jl")
include("actions/main.jl")
include("utils/main.jl")
include("features/main.jl")
include("behaviors/main.jl")
# include("evaluation/main.jl")

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
#=
Helper function for score2pvalue. Find the lowest
score s such that P-value(M_ϵ,s) = P-value(M_ϵ,s-E).

Input:
    Q_m: a sorted dict
    E: the max column-wise error between matrix M and M_ϵ
Output:
=#
function find_s(Q_m::SortedDict{T,T}, E::T) where T <: Real
    keys_ = collect(keys(Q_m)); ℓ = length(keys_);
    if ℓ > 1
        for i = 1:(ℓ-1)
            @inbounds if keys_[i]+E < keys_[i+1]
                return keys_[i];
            end
        end
    end
    return nothing
end

"""
    score2pvalue(pwm, α, ϵ=1e-1, k=100, bg=[.25,.25,.25,.25])
Returns P-value(M,α) of a `pwm` with a given threshold `α`.

Input:
* `pwm`: a 4 x m matrix
* `α`: the score
* `ϵ`: initial granularity  (optional)
* `k`: Refinement parameter (optional)
* `bg`: multinomial background (optional)

Output:
* `pval`: p-value
"""
function score2pvalue(pwm::Matrix{T}, α::Real, ϵ=1e-2, k=10, bg=[.25,.25,.25,.25]) where T <: Real
    @assert size(pwm,1) == 4 "The input matrix must have 4 and only 4 rows"
    pwm_double64 = Double64.(pwm);
    α_double64 = Double64(α);
    bg_double64 = Double64.(bg);
    mpwm = min_score_range(pwm_double64);
    m = size(mpwm,2);
    β = best_score(mpwm)+1;
    pval = Double64(0);
    s = 0;
    i = 1;
    @inbounds while !(α_double64 ≈ s)
        # println(α_double64-s)
        ϵ = i == 1 ? ϵ : ϵ/k; i+=1;
        # isinf(ϵ) || isnan(ϵ) && (println(i); break);
        pwm_ϵ = round_pwm(mpwm, ϵ);
        # any(isnan.(pwm_ϵ)) && (println(i); break);
        E = calc_E(mpwm, pwm_ϵ);
        Q = create_Q(m);
        Q = score_distribution(pwm_ϵ, α-E, β, bg_double64);
        #=
        note:
        Sometimes the score range is small
        (i.e. β-(α-E) is very  small) and hence we won't be
        able to find such s. In such case, we return the
        calculated p-value.

        Note that this p-value is an underestimate. To get a
        more accurate result, set k to be a larger number,
        e.g. k=100.
        =#
        s = find_s(Q[m],E);
        if !isnothing(s)
            for (k,v) in Q[m]
                pval += k ≥ s ? v : 0;
            end
            β = s;
        else
            break;
        end
    end
    return pval
end


text:
# julia
using ArgParse
using BioStructures
# input parameters function
function parse_commandline()

    s = ArgParseSettings(description = "input a pdb ids from PDB and subset by selecting the model and chain")
    @add_arg_table s begin
        "--id"
            help = "input pdb ids"
        "--model"
            help = "model to make the pdb file for"
            arg_type = Int
            default = 1
            required = false
        "--chain"
            help = "chain in model to make the pdb file for"
        "--out"
            help = "output pdb file"
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println(parsed_args)
# main
    downloadpdb(parsed_args["id"]) do fp
        s = read(fp, PDB)
        subpdb = s[parsed_args["model"]][parsed_args["chain"]]
        writepdb(parsed_args["out"], subpdb)
    end
end

main()


text:
"""
    differentiate(f::Expr, x::Symbol)

Differentiate f with respect to x.
"""
function differentiate(f::Expr, x::Symbol)
    @assert f.head == :call
    return differentiate(Val{f.args[1]}, f, x)
end

"""
    differentiate(::Number, ::Symbol)

Differentiate a constant.
"""
function differentiate(::Number, ::Symbol)
    return 0
end

"""
    differentiate(f::Symbol, x::Symbol)

Differentiate a symbol.
"""
function differentiate(f::Symbol, x::Symbol)
    return f == x ? 1 : 0
end

"""
    differentiate(Val{:*}, f::Expr, x::Symbol)

Differentiate (fg)' using chain rule: ```(fgh)' = f'gh + fg'h + fgh'```.
"""
function differentiate(::Type{Val{:*}}, f::Expr, x::Symbol)
    op = first(f.args)
    @assert op == :*
    res_args = Any[:+]
    for i in 2:length(f.args)
        new_args = copy(f.args)
        new_args[i] = differentiate(f.args[i], x)
        push!(res_args, Expr(:call, new_args...))
    end
    return Expr(:call, res_args...)
end

"""
    differentiate(Val{:^}, f::Expr, x::Symbol)

Differentiate d/dx f^a = a * f ^ (a - 1) * diff(f, x)
"""
function differentiate(::Type{Val{:^}}, ex::Expr, x::Symbol)
    op, f, a = ex.args
    @assert op == :^
    df = differentiate(f, x)
    return :($a * $f ^ ($a - 1) * $df)
end

"""
    differentiate(Val{:/}, ex::Expr, x::Symbol)

Differentiate d/dx f/g.
"""
function differentiate(::Type{Val{:/}}, ex::Expr, x::Symbol)
    op, g, h = ex.args
    @assert op == :/
    dg = differentiate(g, x)
    dh = differentiate(h, x)
    return :(($dg * $h + $g * $dh) / $h^2)
end

"""
    differentiate(Val{:+}, ex::Expr, x::Symbol)

Differentiate d/dx (f + g)
"""
function differentiate(::Type{Val{:+}}, ex::Expr, x::Symbol)
    op = first(ex.args)
    @assert op == :+
    args = differentiate.(ex.args[2:end], x)
    return Expr(:call, :+, args...)
end

"""
    differentiate(Val{:-}, ex::Expr, x::Symbol)

Differentiate d/dx (f - g)
"""
function differentiate(::Type{Val{:-}}, ex::Expr, x::Symbol)
    op = first(ex.args)
    @assert op == :-
    args = differentiate.(ex.args[2:end], x)
    return Expr(:call, :-, args...)
end


text:
module NetworkViz
using ThreeJS

    export NodeProperty, EdgeProperty

    type NodeProperty
      color
      size
      shape
    end

    type EdgeProperty
      color
      width
    end

    NodeProperty(c) = NodeProperty(color=c,size=0.2,shape=0)
    NodeProperty(c,s) = NodeProperty(color=c,size=s,shape=0)

    include("spring.jl")
    include("utils.jl")

end


text:
# -------------------------------------------------------
# simple ann template
# -------------------------------------------------------

m = Dense(5, 1)
function loss(x, y)
    sum(y .- m(x))^2
end
opt = Descent(0.01)
x, y = rand(5, 5), fill(0.5, 1, 5)
Flux.train!(loss, params(m), [(x, y)], opt)

# view results
params(m)
loss(x, y)


text:
using Dates

tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_model(config)
@unpack network = model

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    row = csv_first_row(model.writer.csv_path)

    @test row.time == DateTime("2000-01-02T00:00:00")
    @test row.Q ≈ 8.1396354649197f0
    @test row.volume ≈ 2.7535003939625636f7
    @test row.temp_bycoord ≈ 2.390000104904175f0
    @test row.temp_byindex ≈ 2.390000104904175f0
    @test row.Q_6336050 ≈ 0.006160282939850074f0
    @test row.Q_6336510 ≈ 0.029177309556712334f0
    @test row.Q_6836100 ≈ 0.19621282450614713f0
    @test row.Q_6336500 ≈ 0.006089112638001381f0
    @test row.Q_6836190 ≈ 0.0031262850749354237f0
    @test row.Q_6336800 ≈ 0.007770868657277307f0
    @test row.Q_6336900 ≈ 0.006403194169947582f0
    @test row.Q_6336930 ≈ 0.08888787154163148f0
    @test row.Q_6336910 ≈ 0.007071851236520184f0
    @test row.Q_6136500 ≈ 0.0016367337487926633f0
    @test row.Q_6136520 ≈ 0.002084670434294102f0
    @test row.Q_6136150 ≈ 0.006095549758915344f0
    @test row.Q_6136151 ≈ 0.007643634432992056f0
    @test row.Q_6136160 ≈ 3.9199531493174726f0
    @test row.Q_6136202 ≈ 1.4125847550988493f0
    @test row.recharge_1 ≈ -0.05653226176238641f0
end

@testset "NetCDF scalar output" begin
    ds = model.writer.dataset_scalar
    @test ds["time"][1] == DateTime("2000-01-02T00:00:00")
    @test ds["Q"][:][1:20] ≈ [
        0.73993874f0,
        1.4125848f0,
        1.4389194f0,
        1.4036233f0,
        5.7276225f0,
        2.7586424f0,
        2.1080604f0,
        4.1026044f0,
        0.008365918f0,
        3.919953f0,
        4.0615087f0,
        0.006050462f0,
        0.0076436345f0,
        0.00780229f0,
        0.00346879f0,
        0.7066179f0,
        0.0022446192f0,
        1.3323202f0,
        3.8271446f0,
        1.6729931f0,
    ]
    @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
    @test ds["temp_index"][:] ≈ [2.39f0]
    @test ds["temp_coord"][:] ≈ [2.39f0]
    @test keys(ds.dim) == ["time", "Q_gauges", "temp_bycoord", "temp_byindex"]
end

@testset "first timestep" begin
    sbm = model.vertical

    @test sbm.tt[50063] ≈ 0.0f0

    @test model.clock.iteration == 2

    @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
    @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
    @test sbm.runoff[50063] == 0.0
    @test sbm.soilevap[50063] == 0.0
    @test sbm.snow[5] ≈ 3.592840840467347f0
end

# run the second timestep
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "second timestep" begin
    sbm = model.vertical
    @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
    @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
    @test sbm.runoff[50063] == 0.0
    @test sbm.soilevap[50063] ≈ 0.006358004660566856f0
    @test sbm.snow[5] ≈ 3.667748983774868f0
end

@testset "subsurface flow" begin
    ssf = model.lateral.subsurface.ssf
    @test sum(ssf) ≈ 6.370399148012509f7
    @test ssf[network.land.order[1]] ≈ 7.169036749244327f2
    @test ssf[network.land.order[end-100]] ≈ 2335.2465707069578f0
    @test ssf[network.land.order[end]] ≈ 288.19428729403944f0
end

@testset "overland flow" begin
    q = model.lateral.land.q_av
    @test sum(q) ≈ 290.5520014030802f0
    @test q[26625] ≈ 0.0
    @test q[39308] ≈ 0.0
    @test q[network.land.order[end]] ≈ 1.0f-30
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 3615.2084127447865f0
    @test q[1622] ≈ 0.0005986272622329333f0
    @test q[43] ≈ 12.036342425160155f0
    @test q[network.river.order[end]] ≈ 0.03835913312643948f0
end

@testset "reservoir simple" begin
    res = model.lateral.river.reservoir
    @test res.outflow[1] ≈ 0.21750000119148086f0
    @test res.inflow[1] ≈ 43.18479982574888f0
    @test res.volume[1] ≈ 2.751299001489657f7
    @test res.precipitation[1] ≈ 0.17999997735023499f0
    @test res.evaporation[1] ≈ 0.5400000810623169f0
end

# set these variables for comparison in "changed dynamic parameters"
precip = copy(model.vertical.precipitation)
evap = copy(model.vertical.potential_evaporation)
lai = copy(model.vertical.leaf_area_index)
res_evap = copy(model.lateral.river.reservoir.evaporation)

Wflow.close_files(model, delete_output = false)

# test for setting a pit and multithreading multiple basins (by setting 2 extra pits
# resulting in 3 basins)
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)
config["model"]["pits"] = true
config["input"]["pits"] = "wflow_pits"
config.endtime = DateTime(2000, 1, 9)
config.loglevel = "info"

model = Wflow.run(config)

@testset "timing" begin
    # clock has been reset
    calendar = get(config, "calendar", "standard")::String
    @test model.clock.time == Wflow.cftime(config.starttime, calendar)
    @test model.clock.iteration == 1
end

@testset "river flow at basin outlets and downstream of one pit" begin
    q = model.lateral.river.q_av
    @test q[4009] ≈ 8.60480399680283f0 # pit/ outlet, CartesianIndex(141, 228)
    @test q[4020] ≈ 0.006779014715290862f0 # downstream of pit 4009, CartesianIndex(141, 229)
    @test q[2508] ≈ 150.15321264134985f0 # pit/ outlet
    @test q[5808] ≈ 0.12625654862968252f0 # pit/ outlet
end

# test changing forcing and cyclic LAI parameter
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

config.input.vertical.precipitation =
    Dict("scale" => 2.0, "netcdf" => Dict("variable" => Dict("name" => "precip")))
config.input.vertical.potential_evaporation = Dict(
    "scale" => 3.0,
    "offset" => 1.50,
    "netcdf" => Dict("variable" => Dict("name" => "pet")),
)
config.input.vertical.leaf_area_index =
    Dict("scale" => 1.6, "netcdf" => Dict("variable" => Dict("name" => "LAI")))

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "changed dynamic parameters" begin
    res = model.lateral.river.reservoir
    vertical = model.vertical
    @test vertical.precipitation[2] / precip[2] ≈ 2.0f0
    @test (vertical.potential_evaporation[100] - 1.50) / evap[100] ≈ 3.0f0
    @test vertical.leaf_area_index[100] / lai[100] ≈ 1.6f0
    @test (res.evaporation[2] - 1.50) / res_evap[2] ≈ 3.0000012203408635f0
end

# test cyclic river inflow
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)

config.input.cyclic = ["vertical.leaf_area_index", "lateral.river.inflow"]
config.input.lateral.river.inflow = "inflow"

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "river inflow (cyclic)" begin
    @test model.lateral.river.inflow[44] ≈ 0.75
    @test model.lateral.river.q_av[44] ≈ 10.698591283662008
end

# test fixed forcing (precipitation = 2.5)
config = Wflow.Config(tomlpath)
config.input.vertical.precipitation = Dict("value" => 2.5)
model = Wflow.initialize_sbm_model(config)
Wflow.load_fixed_forcing(model)

@testset "fixed precipitation forcing (initialize)" begin
    @test maximum(model.vertical.precipitation) ≈ 2.5
    @test minimum(model.vertical.precipitation) ≈ 0.0
    @test all(isapprox.(model.lateral.river.reservoir.precipitation, 2.5))
end

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "fixed precipitation forcing (first timestep)" begin
    @test maximum(model.vertical.precipitation) ≈ 2.5
    @test minimum(model.vertical.precipitation) ≈ 0.0
    @test all(isapprox.(model.lateral.river.reservoir.precipitation, 2.5))
end

Wflow.close_files(model, delete_output = false)

# test local-inertial option for river flow river_routing
tomlpath = joinpath(@__DIR__, "sbm_config.toml")
config = Wflow.Config(tomlpath)
config.model.river_routing = "local-inertial"

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "river flow and depth (local inertial)" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 3906.99107445248f0
    @test q[1622] ≈ 6.0094181857060604f-5
    @test q[43] ≈ 11.900372477232786f0
    @test q[501] ≈ 3.536628093804679f0
    h = model.lateral.river.h_av
    @test h[1622] ≈ 0.0018099650639472893f0
    @test h[43] ≈ 0.4362704420866992f0
    @test h[501] ≈ 0.05669956233680629f0
end
Wflow.close_files(model, delete_output = false)

# test local-inertial option for river flow river_routing
tomlpath = joinpath(@__DIR__, "sbm_swf_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_model(config)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "river and overland flow and depth (local inertial)" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 2369.4745120995767f0
    @test q[1622] ≈ 6.027872269041313f-5
    @test q[43] ≈ 5.3952531643382144f0
    @test q[501] ≈ 1.6023076297110879f0
    h = model.lateral.river.h_av
    @test h[1622] ≈ 0.001812085885794477f0
    @test h[43] ≈ 0.30122213937420667f0
    @test h[501] ≈ 0.03172037092237322f0
    qx = model.lateral.land.qx
    qy = model.lateral.land.qy
    @test qx[[26, 35, 631]] ≈
          [0.13689085916734353f0, 0.02571727530678835f0, 0.06163019045555399f0]
    @test qy[[26, 35, 631]] ≈
          [0.189187686873958f0, 0.017683557812671172f0, 0.09182928055862745f0]
    h = model.lateral.land.h
    @test h[[26, 35, 631]] ≈
          [0.07359964438114004f0, 0.0070632177001078555f0, 0.0009883334841610835f0]
end
Wflow.close_files(model, delete_output = false)


text:
# using Revise, Plots
using Test
	using BifurcationKit, Parameters, Setfield, LinearAlgebra, ForwardDiff, SparseArrays
	const BK = BifurcationKit
##################################################################
# The goal of these tests is to test all combinaisons of options
##################################################################

norminf = x -> norm(x, Inf)

function Fsl!(f, u, p, t)
	@unpack r, μ, ν, c3 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2)
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1)
	return f
end

Fsl(x, p) = Fsl!(similar(x), x, p, 0.)
jet = BK.getJet(Fsl; matrixfree=false)
####################################################################################################
par_sl = (r = 0.1, μ = 0., ν = 1.0, c3 = 1.0)
u0 = [.001, .001]
par_hopf = (@set par_sl.r = 0.1)
####################################################################################################
# continuation, Hopf bifurcation point detection
optconteq = ContinuationPar(ds = -0.01, detectBifurcation = 3, pMin = -0.5, nInversion = 8)
br, = continuation(jet[1], jet[2], u0, par_hopf, (@lens _.r), optconteq)
####################################################################################################
Ntst = 3
m = 4
N = 3
coll_cache = BK.POOrthogonalCollocationCache(Ntst, m)
const Mf = rand(N,N)
prob_col = BK.PeriodicOrbitOCollProblem(F = (x,p) -> Mf * x.^2, N = N, coll_cache = coll_cache, ϕ = ones(N*( 1 + m * Ntst)), xπ = zeros(N*( 1 + m * Ntst)))
size(prob_col)
length(prob_col)
BK.getTimes(prob_col)
BK.getMaxTimeStep(prob_col)
size(coll_cache)
BK.updateMesh!(prob_col, prob_col.coll_cache.mesh)
PeriodicOrbitOCollProblem(10,2) |> BK.getMeshSize
BK.getLs(prob_col)
show(prob_col)

_orbit(t) = [cos(2pi*t)] * sqrt(par_sl.r/par_sl.c3)
_ci = BK.generateSolution(prob_col, _orbit, 1.)
BK.getPeriodicOrbit(prob_col, _ci, par_sl)
BK.getMaximum(prob_col, _ci, par_sl)
BK.∂(sin, 2)(0.)
prob_col(_ci, par_sl) #|> scatter
BK.getTimeSlices(prob_col, _ci)
# interpolate solution
sol = BK.POOcollSolution(prob_col, _ci)
sol(rand())

# using ForwardDiff
# J(x,p) = ForwardDiff.jacobian(u -> prob_col(u,  p), x)
# _J = J(vcat(vec(_ci), 1),  par_sl)
# 	heatmap(_J .!= 0, yflip = true)
####################################################################################################
Ntst = 50
m = 4
N = 2
coll_cache = BK.POOrthogonalCollocationCache(Ntst, m)
prob_col = BK.PeriodicOrbitOCollProblem(F = Fsl, N = 2, coll_cache = coll_cache, ϕ = zeros(N*( 1 + m * Ntst)), xπ = zeros(N*( 1 + m * Ntst)))
prob_col.ϕ[2] = 1

_orbit(t) = [cos(t), sin(t)] * sqrt(par_sl.r/par_sl.c3)
_ci = BK.generateSolution(prob_col, _orbit, 2pi)
prob_col(_ci, par_sl)
@test prob_col(_ci, par_sl) |> norminf < 1e-7

# using ForwardDiff
# J(x,p) = ForwardDiff.jacobian(u -> prob_col(u,  p), x)
# _J = J(vcat(vec(_ci), 2pi),  par_sl)
# 	heatmap(_J .!=0, yflip = true)

args = (tangentAlgo = BorderedPred(),
	# jacobianPO = :autodifSparse,
	# linearAlgo = MatrixBLS(),
	recordFromSolution = (x,p) -> (norminf(x[1:end-1])),
	plotSolution = (x,p; k...) -> begin
		outt = getPeriodicOrbit(prob_col, x, p)
		plot!(vec(outt.t), outt.u[1,:]; k...)
	end,
	finaliseSolution = (z, tau, step, contResult; k...) -> begin
		# BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals)
		return true
	end,)

optcontpo = setproperties(optconteq; detectBifurcation = 2, precisionStability = 1e-7)
@set! optcontpo.ds = -0.01
@set! optcontpo.newtonOptions.verbose = false

newton(prob_col, _ci, par_sl, optcontpo.newtonOptions)

br_po, = @time continuation(prob_col, _ci, par_sl, (@lens _.r), optcontpo;
	verbosity = 0, plot = false,
	args...,
	updateSectionEveryStep = 1,
	)
####################################################################################################
# test  Hopf aBS
br_po, = continuation(jet..., br, 1, ContinuationPar(optcontpo; ds = 0.01),
	PeriodicOrbitOCollProblem(20, 5);
	jacobianPO = :autodiffDense,
	δp = 0.1,
	usedeflation = true,
	# regular continuation options
	verbosity = 0,	plot = false,
	updateSectionEveryStep = 1,
	)


text:
# TODO: this file is pretty sparse...

using InteractiveUtils: subtypes
using Core.Compiler: SSAValue, Const, GotoNode

function const_bool_retval(f, sig)
    ctyped = code_typed(f, sig)[1].first
    retval = ctyped.code[end]
    retval = isa(retval, GotoNode) ? ctyped.code[retval.label] : retval
    retval = Meta.isexpr(retval, :return) ? retval.args[1] : retval
    retval = isa(retval, SSAValue) ? ctyped.ssavaluetypes[retval.id] : retval
    retval = isa(retval, Const) ? retval.val : retval
    isa(retval, Bool) && return retval
    error("did not infer constant boolean return value for ", f, sig)
end

@context canrecurseTestCtx
ctx = canrecurseTestCtx()
# Test that the result value of `canrecurse` is inferred exactly for `T <: Core.Builtin`.
# Note that these aren't necessarily valid full call signatures, this is just
# testing that Cassette's `canrecurse` implementation has an inferrable
# short-circuiting `Core.Builtin` check.
for T in subtypes(Core.Builtin)
    if !(T <: typeof(Core._apply))
        @test !const_bool_retval(canrecurse, (typeof(ctx), T))
        @test !const_bool_retval(canrecurse, (typeof(ctx), typeof(Core._apply), T))
        @test !const_bool_retval(canrecurse, (typeof(ctx), typeof(Core.invoke), T))
    end
end

@test canrecurse(ctx, hypot, 1, 2)
@test canrecurse(ctx, Core.invoke, hypot, Tuple{Int,Int}, 1, 2)
@test canrecurse(ctx, Core._apply, hypot, (1, 2))
@test canrecurse(ctx, Core._apply, Core.invoke, (hypot, Tuple{Int,Int}, 1, 2))

###########

# Declaring a `@context` with the same name twice, should be a No-Op, not an error
# Code below will throw an error if this is not true

@context FooBar
@context FooBar

# Test keyword function detection
kwtest_1(; x) = x  # Normal function
const kwtest_2 = (; x) -> x  # Anonymous function
const kwtest_3  = let; (; x) -> x end  # Closure
for n in 1:3
    F = typeof(getfield(Main, Symbol(:kwtest_, n)))
    @test iskwftype(Core.kwftype(F))
    @test !iskwftype(F)
end


text:
"""
```julia
p3p(points, pdn_pixels::AbstractVector{SVector{3, T}}, K)
```

Recover pose using P3P algorithm.

# Arguments:

- `points: Vector of 3D points in `(x, y, z)` format.
- `pdn_pixels::AbstractVector{SVector{3, Float64}}`:
    Corresponding projections of `points` onto image plane,
    predivided by `K` intrinsics and normalized.
    E.g.: `Ki = inv(K); p = Ki [x, y, 1]; p /= norm(p)`.
- `K::SMatrix{3, 3, Float64}`: Camera intrinsics.

# Returns:

`Vector{SMatrix{3, 4, Float64}}` vector of up to 4 possible solutions.
Each element is a projection matrix `P = K * [R|t]`.
To get pure transformation matrix, multiply `P` by `inv(K)`.

# References:

```
Link: https://cmp.felk.cvut.cz/~pajdla/gvg/GVG-2016-Lecture.pdf
chapter: 7.3 Calibrated camera pose computation.
pages: 51-59
```
"""
function p3p(
    points, pdn_pixels::AbstractVector{SVector{3, T}}, K,
) where T <: Real
    models = SMatrix{3, 4, Float64}[]
    ϵ = 1e-4

    d12 = norm(points[1] .- points[2])
    d13 = norm(points[1] .- points[3])
    d23 = norm(points[2] .- points[3])
    (d12 < ϵ || d13 < ϵ || d23 < ϵ) && return models

    sd12, sd13, sd23 = d12^2, d13^2, d23^2

    α12 = pdn_pixels[1] ⋅ pdn_pixels[2]
    α13 = pdn_pixels[1] ⋅ pdn_pixels[3]
    α23 = pdn_pixels[2] ⋅ pdn_pixels[3]

    # Compute coefficients of a polynomial [7.80].
    m1, m2 = sd12, sd13 - sd23
    p1 = PPolynomial{Float64, :x}([0.0, -2.0 * α23 * sd12])
    p2 = PPolynomial{Float64, :x}([2.0 * α13 * sd23, - 2.0 * α23 * sd13])
    q1 = sd23 * PPolynomial{Float64, :x}([1.0, - 2.0 * α12, 1.0]) +
        PPolynomial{Float64, :x}([0.0, 0.0, -sd12])
    q2 = PPolynomial{Float64, :x}([sd23, 0.0, -sd13])

    P::PPolynomial{Float64, :x} =
        (m1 * q2 - m2 * q1)^2 - (m1 * p2 - m2 * p1) * (q1 * p2 - q2 * p1)
    P_roots::Vector{ComplexF64} = P |> Polynomials.roots

    for η12r in P_roots
        isreal(η12r) || continue
        η12 = real(η12r)
        η12 ≤ 0 && continue
        # Evaluate [7.89] at η12 to get η13.
        η13 = (m1 * q2 - m2 * q1)(η12) / (m1 * p2 - m2 * p1)(η12)
        # [7.91 - 7.93].
        denom = 1.0 + η12^2 - 2.0 * η12 * α12
        denom ≤ 0 && continue
        η1 = d12 / √denom
        η2 = η1 * η12
        η3 = η1 * η13
        (η1 ≤ 0 || η2 ≤ 0 || η3 ≤ 0) && continue
        # Compute and check errors.
        (
            abs((√(η1^2 + η2^2 - 2 * η1 * η2 * α12) - d12) / d12) > ϵ ||
            abs((√(η1^2 + η3^2 - 2 * η1 * η3 * α13) - d13) / d13) > ϵ ||
            abs((√(η2^2 + η3^2 - 2 * η2 * η3 * α23) - d23) / d23) > ϵ
        ) && continue

        nx1 = η1 * pdn_pixels[1]
        z2 = η2 * pdn_pixels[2] - nx1; z2 /= norm(z2)
        z3 = η3 * pdn_pixels[3] - nx1; z3 /= norm(z3)
        z1 = z2 × z3; z1 /= norm(z1)
        z3z1 = z3 × z1

        zw2 = @. (points[2] - points[1]) / d12
        zw3 = @. (points[3] - points[1]) / d13
        zw1 = zw2 × zw3; zw1 /= norm(zw1)
        zw3zw1 = zw3 × zw1

        # Recover rotation R [7.130 - 7.134].
        Z = SMatrix{3, 3, Float64}(z1..., z2..., z3z1...)
        ZW = SMatrix{3, 3, Float64}(zw1..., zw2..., zw3zw1...)
        R = Z * inv(ZW)
        KR = K * R
        Kt = -KR * (points[1] - R' * nx1)
        push!(models, SMatrix{3, 4, Float64}(KR..., Kt...))
    end
    models
end

"""
```julia
p3p_select_model(models, points, pixels; threshold = 1.0)
```

Select best pose from `models`.

# Arguments:

- `models`: Projection matrices from which to select best one.
- `points`: 3D points in `(x, y, z)`
- `pixels`: Corresponding projections onto image plane in `(x, y)` format.
- `threshold`: Maximum distance in pixels between projected point
    and its target pixel, for the point to be considered inlier.
    Default value is `1.0`.

# Returns:

`n_inliers, (projection, inliers, error)`.

- `n_inliers`: Number of inliers for the `projection`.
- `projection`: `K * P` projection matrix that projects `points`
    onto image plane.
- `inliers`: Boolean vector indicating which point is inlier.
- `error`: Average reprojection error for the `pose`.
"""
function p3p_select_model(models, points, pixels; threshold = 1.0)
    best_n_inliers = 0
    best_error = maxintfloat()
    best_projection = nothing
    best_inliers = nothing

    for P in models
        avg_error = 0.0
        inliers = fill(false, length(points))
        n_inliers = 0
        for (i, (pixel, point)) in enumerate(zip(pixels, points))
            projected = P * SVector{4}(point..., 1.0)
            projected = projected[[1, 2]] ./ projected[3]
            error = norm(pixel .- projected)

            avg_error += error
            error ≥ threshold && continue

            inliers[i] = true
            n_inliers += 1
        end
        avg_error /= length(points)
        n_inliers < 3 && continue

        if avg_error < best_error
            best_projection = P
            best_inliers = inliers
            best_n_inliers = n_inliers
            best_error = avg_error
        end
    end

    best_n_inliers, (best_projection, best_inliers, best_error)
end

"""
```julia
pre_divide_normalize(pixels, K)
```

Divide `pixels` by `K` intrinsic matrix and normalize them.
`pixels` in `(x, y)` format.
"""
function pre_divide_normalize(pixels, K)
    res = Vector{SVector{3, Float64}}(undef, length(pixels))
    @inbounds for (i, px) in enumerate(pixels)
        p = SVector{3, Float64}(
            (px[1] - K[1, 3]) / K[1, 1], (px[2] - K[2, 3]) / K[2, 2], 1)
        res[i] = normalize(p)
    end
    res
end

function p3p(points, pixels::AbstractVector{SVector{2, T}}, K) where T <: Real
    p3p(points, pre_divide_normalize(pixels, K), K)
end

"""
```julia
p3p_ransac(points, pixels, pdn_pixels, K; threshold = 1.0, ransac_kwargs...)
```

Recover pose `K*[R|t]` using P3P Ransac algorithm.

# Arguments:

- `points`: 3D points in `(x, y, z)`
- `pixels`: Corresponding projections onto image plane in `(x, y)` format.
- `pdn_pixels`: Corresponding projections onto image plane,
    predivided by `K` intrinsics and normalized.
    E.g.: `Ki = inv(K); p = Ki [x, y, 1]; p /= norm(p)`.
- `K`: Camera intrinsics.
- `threshold`: Maximum distance in pixels between projected point
    and its target pixel, for the point to be considered inlier.
    Default value is `1.0`.
- `ransac_kwargs...`: Keyword arguments passed to [`ransac`](@ref).

# Returns:

`n_inliers, (projection, inliers, error)`.

- `n_inliers`: Number of inliers for the `projection`.
- `projection`: `K * P` projection matrix that projects `points`
    onto image plane.
- `inliers`: Boolean vector indicating which point is inlier.
- `error`: Average reprojection error for the `pose`.

# References:

```
Link: https://cmp.felk.cvut.cz/~pajdla/gvg/GVG-2016-Lecture.pdf
chapter: 7.3 Calibrated camera pose computation.
pages: 51-59
```
"""
function p3p_ransac(
    points, pixels, pdn_pixels, K; threshold = 1.0, ransac_kwargs...,
)
    sample_selection(sample_ids) = (points[sample_ids], pdn_pixels[sample_ids], K)
    rank(models; sample_ids) = p3p_select_model(models, points, pixels; threshold)
    ransac(sample_selection, p3p, rank, length(points), 3; ransac_kwargs...)
end

"""
```julia
p3p_ransac(points, pixels, K; threshold = 1.0, ransac_kwargs...)
```

Recover pose `K*[R|t]` using P3P Ransac algorithm.

# Arguments:

- `points`: 3D points in `(x, y, z)`
- `pixels`: Corresponding projections onto image plane in `(x, y)` format.
    These values, **WILL** be predivided and normalized by the intrinsic matrix.
    E.g.: `Ki = inv(K); p = Ki [x, y, 1]; p /= norm(p)`.
- `K`: Camera intrinsics.
- `threshold::Real`:
    Maximum distance in pixels between projected point and its target pixel,
    for the point to be considered inlier. Default value is `1.0`.
- `ransac_kwargs...`: Keyword arguments passed to [`ransac`](@ref).

# Returns:

`Vector{SMatrix{3, 4, Float64}}` vector of up to 4 possible solutions.
Each element is a projection matrix `P = K * [R|t]`.
To get pure transformation matrix, multiply `P` by `inv(K)`.

# References:

```
Link: https://cmp.felk.cvut.cz/~pajdla/gvg/GVG-2016-Lecture.pdf
chapter: 7.3 Calibrated camera pose computation.
pages: 51-59
```
"""
function p3p_ransac(points, pixels, K; threshold = 1.0, ransac_kwargs...)
    pdn_pixels = pre_divide_normalize(pixels, K)
    p3p_ransac(points, pixels, pdn_pixels, K; threshold, ransac_kwargs...)
end