text:
# This code shows the functions to calculate the rotor losses in high speed PM machines
# author: Jaime Renedo Anglada (renedo.jaime@ieee.org)

function Calc_rotor_loss( n_mech, space_order, p, k_time, B_r_ext,L,R_1,R_2,R_3,R_4,R_5,sigma_1,sigma_2,sigma_3,sigma_4,sigma_5,mu_r,r_wave)
# Calculation of the rotor losses in one function

    # define geometry:
    mu_0=4*pi*10.0^-7; # [m kg s^-2 A^-2]
    q=p*space_order; # harmonic order
    k_time=1*6; # time order
    f_1=n_mech/60*p; # [Hz]
    f=f_1;

    # Permeability:
    mu_1=750*mu_0; # [m kg s^-2 A^-2]
    mu_2=mu_r*mu_0; # [m kg s^-2 A^-2]
    mu_3=1*mu_0; # [m kg s^-2 A^-2]
    mu_4=1*mu_0; # [m kg s^-2 A^-2]
    mu_5=5000*mu_0; # [m kg s^-2 A^-2]

    B_r_ext=B_r_ext*(R_3/r_wave)^(q-1); # where you measured the amplitude of the harmonics

    ### No eddy currents:
    f_no_eddy=0.01;
    omega=2*pi*f_no_eddy; # [rad/s] frequency
    J_kq=1; # current sheet density
    sigma_air=3*10.0^-15; # [S/m]

    # Propagation constants:
    k_p_1=sqrt(im*omega*mu_1*sigma_air);
    k_p_2=sqrt(im*omega*mu_2*sigma_air);
    k_p_3=sqrt(im*omega*mu_3*sigma_air);
    k_p_4=sqrt(im*omega*mu_4*sigma_air);
    k_p_5=sqrt(im*omega*mu_5*sigma_air);

    # Region 1:
    #beta_1=mu_0*im*k_p_1*R_1*besseli_d(q,k_p_1*R_1)/besseli(q,k_p_1*R_1)/mu_1/q;
    beta_1=mu_0*im/mu_1/q;

    # Region 2:
    #T_2=TransferMatrix_nc( R_1, R_2, mu_2, sigma_air, omega, q);
    T_2=TransferMatrix_simp( R_1, R_2, mu_2, sigma_air, omega, q, N_div);
    #         T_2_pm=TransferMatrix_medium_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);
    #         T_2_pm=TransferMatrix_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);
    #         T_2_pm=TransferMatrix_super_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);
    println("$T_2")

    # Region 3:
    #T_3=TransferMatrix_nc( R_2, R_3, mu_3, sigma_air, omega, q);
    T_3=TransferMatrix_simp( R_2, R_3, mu_3, sigma_air, omega, q, N_div);
    println("$T_3")

    # Region 4:
    #T_4=TransferMatrix_nc( R_3, R_4, mu_4, sigma_air, omega, q);
    T_4=TransferMatrix_simp( R_3, R_4, mu_4, sigma_air, omega, q, N_div);
    println("$T_4")

    # Region 5:
    #T_5=TransferMatrix_nc( R_4, R_5, mu_5, sigma_air, omega, q);
    T_5=TransferMatrix_simp( R_4, R_5, mu_5, sigma_air, omega, q, N_div);
    println("$T_5")

    # Estimation of the coefficients:

    Mat_D=[0.0 0.0; 0.0 1.0]-T_5*T_4*T_3*T_2*[1.0 0.0; beta_1 0.0];
    println("Mat_D = $Mat_D")
    #         temp_mat=inv(Mat_D);
    vec_sol=inv(Mat_D)*T_5*[0; mu_0*J_kq];
    println("vec sol = $vec_sol")
    # Get the fields:
    B_1=vec_sol[1];
    H_1=-beta_1*B_1/mu_0;
    println("B_1 = $B_1")
    println("H_1 = $H_1 [A/m]")

    H_5=vec_sol[2]/mu_0;

    temp_2=T_2*[B_1; mu_0*(H_1)];
    B_2=temp_2[1];
    H_2=temp_2[2]/mu_0;
    println("B_2 = $B_2 [T]")
    println("H_2 = $H_2 [A/m]")

    temp_3=T_3*[B_2; mu_0*(H_2)];
    B_3=temp_3[1];
    H_3=temp_3[2]/mu_0;
    println("B_3 = $B_3 [T]")
    println("H_3 = $H_3 [A/m]")

    temp_4=T_4*[B_3; mu_0*(H_3)];
    B_4=temp_4[1];
    H_4=temp_4[2]/mu_0;
    println("B_4 = $B_4 [T]")
    println("H_4 = $H_4 [A/m]")

    # Field at the rotor hub:
    vec_hub_1=real(B_1+mu_1*H_1);
    vec_hub_1=abs(B_1);

    vec_hub_2=imag(B_1+mu_1*H_1);
    vec_hub_2=abs(mu_1*H_1);

    B_r_hub_ne=vec_hub_1;
    B_theta_hub_ne=vec_hub_2;

    # Field at the outer part of the PM:
    vec_out_1=real(B_2+mu_2*H_2);
    vec_out_1=abs(B_2);

    vec_out_2=imag(B_2+mu_2*H_2);
    vec_out_2=abs(mu_2*H_2);

    B_r_pm_ne=vec_out_1;
    B_theta_pm_ne=vec_out_2;

    # Field at the outer part of the sleeve:
    vec_out_1=real(B_3+mu_3*H_3);
    vec_out_1=abs(B_3);

    vec_out_2=imag(B_3+mu_3*H_3);
    vec_out_2=abs(mu_3*H_3);

    B_r_sleeve_ne=vec_out_1;
    B_theta_sleeve_ne=vec_out_2;

    ### with eddy currents:
    f_2=f
    q=q;

    omega=2*pi*f_2*k_time; # [rad/s] frequency
    println("omega = $omega")

    # Propagation constants:
    k_p_1=sqrt(im*omega*mu_1*sigma_1);
    k_p_2=sqrt(im*omega*mu_2*sigma_2);
    k_p_3=sqrt(im*omega*mu_3*sigma_3);
    k_p_4=sqrt(im*omega*mu_4*sigma_4);
    k_p_5=sqrt(im*omega*mu_5*sigma_5);

    # Region 1:
    beta_1=mu_0*im*k_p_1*R_1*besseli_d(q,k_p_1*R_1)/besseli(q,k_p_1*R_1)/mu_1/q;
    beta_1=mu_0*im/mu_1/q;

    # Region 2:
    T_2=TransferMatrix_comp( R_1, R_2, mu_2, sigma_2, omega, q);
    #         T_2_pm=TransferMatrix_medium_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);
    #         T_2_pm=TransferMatrix_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);
    #         T_2_pm=TransferMatrix_super_simp( R_1, R_2, mu_2, sigma_2, omega, q, N_div);

    # Region 3:
    T_3=TransferMatrix_comp( R_2, R_3, mu_3, sigma_3, omega, q);

    # Region 4:
    T_4=TransferMatrix_comp( R_3, R_4, mu_4, sigma_4, omega, q);

    # Region 5:
    T_5=TransferMatrix_comp( R_4, R_5, mu_5, sigma_5, omega, q);

    # Estimation of the coefficients:

    Mat_D=[0.0 0.0; 0.0 1.0]-T_5*T_4*T_3*T_2*[1.0 0.0; beta_1 0.0];
    #         temp_mat=inv(Mat_D);
    vec_sol=inv(Mat_D)*T_5*[0; mu_0*J_kq];
    # Get the fields:
    B_1=vec_sol[1];
    H_1=-beta_1*B_1/mu_0;
    println("B_1 = $B_1 [T]")
    println("H_1 = $H_1 [A/m]")

    H_5=vec_sol[2]/mu_0;
    println("H_5 = $H_5 [A/m]")

    temp_2=T_2*[B_1; mu_0*(H_1)];
    B_2=temp_2[1];
    H_2=temp_2[2]/mu_0;
    println("B_2 = $B_2 [T]")
    println("H_2 = $H_2 [A/m]")

    temp_3=T_3*[B_2; mu_0*(H_2)];
    B_3=temp_3[1];
    H_3=temp_3[2]/mu_0;
    println("B_3 = $B_3 [T]")
    println("H_3 = $H_3 [A/m]")

    temp_4=T_4*[B_3; mu_0*(H_3)];
    B_4=temp_4[1];
    H_4=temp_4[2]/mu_0;
    println("B_4 = $B_4 [T]")
    println("H_4 = $H_4 [A/m]")

    # Field at the rotor hub:
    vec_hub_1=real(B_1+mu_1*H_1);
    vec_hub_1=abs(B_1);

    vec_hub_2=imag(B_1+mu_1*H_1);
    vec_hub_2=abs(mu_1*H_1);

    B_r_hub=vec_hub_1;
    B_theta_hub=vec_hub_2;

    # Field at the outer part of the PM:
    vec_out_1=real(B_2+mu_2*H_2);
    vec_out_1=abs(B_2);

    vec_out_2=imag(B_2+mu_2*H_2);
    vec_out_2=abs(mu_2*H_2);

    B_r_pm=vec_out_1;
    B_theta_pm=vec_out_2;

    # Field at the outer part of the sleeve:
    vec_out_1=real(B_3+mu_3*H_3);
    vec_out_1=abs(B_3);

    vec_out_2=imag(B_3+mu_3*H_3);
    vec_out_2=abs(mu_3*H_3);

    B_r_sleeve=vec_out_1;
    B_theta_sleeve=vec_out_2;

    ### Calculate the coefficients and power losses:

    S_1=2*pi*R_1*L; # area surf rotor steel
    S_2=2*pi*R_2*L; # area surf PM
    S_3=2*pi*R_3*L; # area surf sleeve
    S_4=2*pi*R_4*L; # area surf bore

    # region 2 (PM):
    mu=mu_2
    sigma=sigma_2
    q=q
    k_p=k_p_2;

    F_2=mu_0*im*k_p*q/(mu*R_1)*(besselk(q,k_p*R_1)*besseli_d(q,k_p*R_1)-besseli(q,k_p*R_1)*besselk_d(q,k_p*R_1));
    M_2=[-mu_0*k_p/mu*besselk_d(q,k_p*R_1) -im*q/R_1*besselk(q,k_p*R_1); mu_0*k_p/mu*besseli_d(q,k_p*R_1) im*q/R_1*besseli(q,k_p*R_1)];
    N_2=[im*q/R_2*besseli(q,k_p*R_2) im*q/R_2*besselk(q,k_p*R_2); -mu_0*k_p/mu*besseli_d(q,k_p*R_2) -mu_0*k_p/mu*besselk_d(q,k_p*R_2)];
    T_mat=[1.0 0.0; 0.0 -1.0]*N_2*M_2/F_2*[1.0 0.0; 0.0 -1.0];
    println("F_2 = $F_2")
    println("M_2 = $M_2")

    temp=M_2*[B_1; mu_0*H_1]/F_2;

    C_2=temp[1];
    D_2=temp[2];
    println("C_2 = $C_2")
    println("D_2 = $D_2")

    E_2=-im*omega*(C_2*besseli(q,k_p_2*R_2)+D_2*besselk(q,k_p_2*R_2)); # electric field intensity
    H_theta_2=-k_p_2*(C_2*besseli_d(q,k_p_2*R_2)+D_2*besselk_d(q,k_p_2*R_2))/mu_0; # H_theta complex value

    E_2=omega*R_2/q*B_2; # electric field intensity
    H_theta_2=H_2; # H_theta complex value

    println("E_2 = $E_2")
    println("H_theta_2 = $H_theta_2")

    # region 3 (sleeve):
    mu=mu_3
    sigma=sigma_3
    q=q
    k_p=k_p_3;

    F_3=mu_0*im*k_p*q/(mu*R_2)*(besselk(q,k_p*R_2)*besseli_d(q,k_p*R_2)-besseli(q,k_p*R_2)*besselk_d(q,k_p*R_2));
    M_3=[-mu_0*k_p/mu*besselk_d(q,k_p*R_2) -im*q/R_2*besselk(q,k_p*R_2); mu_0*k_p/mu*besseli_d(q,k_p*R_2) im*q/R_2*besseli(q,k_p*R_2)];
    N_3=[im*q/R_3*besseli(q,k_p*R_3) im*q/R_3*besselk(q,k_p*R_3); -mu_0*k_p/mu*besseli_d(q,k_p*R_3) -mu_0*k_p/mu*besselk_d(q,k_p*R_3)];
    T_mat=[1.0 0.0; 0.0 -1.0]*N_3*M_3/F_3*[1.0 0.0; 0.0 -1.0];
    println("F_3 = $F_3")
    println("M_3 = $M_3")

    temp=M_3*[B_2; mu_0*H_2]/F_3;

    C_3=temp[1];
    D_3=temp[2];
    println("C_3 = $C_3")
    println("D_3 = $D_3")

    E_3=-im*omega*(C_3*besseli(q,k_p_3*R_3)+D_3*besselk(q,k_p_3*R_3)); # electric field intensity
    H_theta_3=-k_p_3*(C_3*besseli_d(q,k_p_3*R_3)+D_2*besselk_d(q,k_p_3*R_3))/mu_0; # H_theta complex value

    E_3=omega*R_3/q*B_3; # electric field intensity
    H_theta_3=H_3; # H_theta complex value

    println("E_3 = $E_3")
    println("H_theta_3 = $H_theta_3")

    # Test of matrix:
    T_2_comp=T_3*T_2
    println("T_2_comp = $T_2_comp")

    # region 4 (air-gap):
    mu=mu_4
    sigma=sigma_4
    q=q
    k_p=k_p_4;

    F_4=mu_0*im*k_p*q/(mu*R_3)*(besselk(q,k_p*R_3)*besseli_d(q,k_p*R_3)-besseli(q,k_p*R_3)*besselk_d(q,k_p*R_3));
    M_4=[-mu_0*k_p/mu*besselk_d(q,k_p*R_3) -im*q/R_3*besselk(q,k_p*R_3); mu_0*k_p/mu*besseli_d(q,k_p*R_3) im*q/R_3*besseli(q,k_p*R_3)];
    N_4=[im*q/R_4*besseli(q,k_p*R_4) im*q/R_4*besselk(q,k_p*R_4); -mu_0*k_p/mu*besseli_d(q,k_p*R_4) -mu_0*k_p/mu*besselk_d(q,k_p*R_4)];
    T_mat=[1.0 0.0; 0.0 -1.0]*N_3*M_3/F_3*[1.0 0.0; 0.0 -1.0];

    temp=M_4*[B_3; mu_0*H_3]/F_4;

    C_4=temp[1];
    D_4=temp[2];
    #println("C_4 = $C_4")
    #println("D_4 = $D_4")

    # Scale coefficient:
    K_B=B_r_ext/B_r_sleeve_ne;
    println("K_B = $K_B")

    # Power going towards the rotor sleeve:
    P3_temp=1/2*real(E_3*conj(H_theta_3))*S_3*(K_B)^2;
    Loss_factor_temp=P3_temp/B_r_ext^2/omega/S_3;

    # Power going towards the rotor PM:
    P2_temp=1/2*real(E_2*conj(H_theta_2))*S_2*(K_B)^2;

    # Power going towards the rotor hub:
    E_1=-im*omega*(C_2*besseli(q,k_p_2*R_1)+D_2*besselk(q,k_p_2*R_1)); # electric field intensity
    H_theta_1=-k_p_2*(C_2*besseli_d(q,k_p_2*R_1)+D_2*besselk_d(q,k_p_2*R_1))/mu_0; # H_theta complex value

    E_1=omega*R_1/q*B_1; # electric field intensity
    H_theta_1=H_1; # H_theta complex value

    P1_temp=1/2*real(E_1*conj(H_theta_1))*S_1*(K_B)^2;

    P_mag=P2_temp-P1_temp;
    P_hub=P1_temp;
    P_sleeve=P3_temp-P2_temp;
    Loss_factor=Loss_factor_temp;

    println("Electrical frequency f = $f_2 [Hz]")
    println("Rotor hub radius R_1= $(R_1*1000) [mm]")
    println("Thickness of the sleeve t_sl= $(t_sl*1000) [mm]")

    println("External field at the sleeve B_r_ext = $(B_r_ext*1000) [mT]")
    println("wavelength of the field (q=$q), lambda= $(2*pi*R_3/q*1000) [mm]")

    println("Total rotor loss P_sleeve = $P_sleeve [W]")
    println("Total magnet loss P_mag = $P_mag [W]")
    println("P_2 = $P2_temp [W]")
    println("P_3 = $P3_temp [W]")
return P3_temp
end


text:
function relative_prime(a, b)
  gcd(a, b) == 1
end

function find_an_e_that_is_relative_prime(p, q)
  p1q1 = (p-1) * (q-1)

  candidate = -1
  found = false

  while !found
    # Select a random candidate and then loop until relative_prime or reached p1q1
    candidate = round( 1 + rand() * p1q1 )

    while candidate < p1q1 && !relative_prime(candidate, p1q1)
      candidate += 1
    end

    if relative_prime(candidate, p1q1)
      found = true
    end
  end

  return candidate
end

function find_a_d(p, q, e)
  p1q1 = (p-1) * (q-1)
  sgd = gcd(e, p1q1)
  d = e + p1q1 / sgd
  return d
end

function rand_big_int(keylength = 4096)
  bigint = BigInt(0)
  bits_left = keylength
  max_bits = 30
  while bits_left > 0
    numbits = min(bits_left, max_bits)
    bigint = (bigint << numbits) + (rand(1:2^numbits) - 1)
    bits_left -= numbits
  end
  bigint
end

type RSAKey
  # Public key = (n, e) and Private key = (n, d)
  n::Integer
  e::Integer
  d::Integer

  RSAKey(p = rand_big_int(4096)) = begin

    q = next_prime(p)
    n = p * q
    e = find_an_e_that_is_relative_prime(p, q)
    d = find_a_d(p, q, e)

    new(n, e, d)
  end
end

type RSAPublicKey
  n::Integer
  e::Integer
  RSAPublicKey(n::Integer, e::Integer) = new(n, e)
end

function public_key(rsa::RSAKey)
  RSAPublicKey(rsa.n, rsa.e)
end

function encrypt(rsa::Union(RSAKey, RSAPublicKey), P::Integer)
  fast_exponentiation(P, rsa.e, rsa.n)
end

function uint8_array_to_integer(arr::Array{Uint8,1})
  integer = BigInt(0)
  for i = 1:length(arr)
    integer = (integer << 8) + arr[i]
  end
  integer
end

function integer_to_uint8_array(x::Integer)
  bytes = Uint8[]

  while x > 0
    b = uint8(x & 0xff)
    push!(bytes, b)
    x = x >>> 8
  end

  return reverse(bytes)
end

# Convert a string to an integer for RSA encryption.
function string_to_integer(s::String)
  bytes = convert(Vector{Uint8}, utf8(s))
  uint8_array_to_integer(bytes)
end

using Codecs

# Convert a non-negative integer i into an octet string and then base64
# encode it so we can pass it around.
function integer_to_base64_string(x::Integer, len = nothing)
  if typeof(len) <: Integer && (x >= 256^len)
    throw("integer is too large")
  end

  if x < 0
    throw("integer is negative")
  end

  bytes = Uint8[]

  while x > 0
    b = uint8(x & 0xff)
    push!(bytes, b)
    x = x >>> 8
  end

  bytes = reverse(bytes)

  if typeof(len) <: Integer && (length(bytes) < len)
    bytes = [zeros(Uint8, len - str), bytes]
  end

  encoded = encode(Base64, bytes)

  return bytestring(encoded)
end

function base64_string_to_integer(s::String)
  decoded = decode(Base64, s)
  uint8_array_to_integer(decoded)
end

function encrypt(rsa::Union(RSAKey, RSAPublicKey), P::String)
  C = fast_exponentiation(string_to_integer(P), rsa.e, rsa.n)
end

function decrypt(rsa::RSAKey, C::String)
  P = fast_exponentiation(string_to_integer(C), rsa.d, rsa.n)
  integer_to_string(P)
end

ext: jl


text:
tryconvert(::Type{T}, x::T) where {T} = Ok(x)  # TODO: should it be `Ok{T}(x)`?

const MightHaveSize = Union{AbstractArray,AbstractDict,AbstractSet,AbstractString,Number}

trygetlength(xs::MightHaveSize)::Result =
    if IteratorSize(xs) isa Union{HasLength,HasShape}
        return Ok(length(xs))
    else
        return Causes.notimplemented(trygetlength, (xs,))
    end

trygeteltype(xs) = trygeteltype(typeof(xs))
trygeteltype(T::Type) = Causes.notimplemented(trygeteltype, (T,))
trygeteltype(::Type{Union{}}) = Causes.notimplemented(trygeteltype, (Union{},))
trygeteltype(::Type{<:AbstractArray{T}}) where {T} = Ok(T)
trygeteltype(::Type{AbstractSet{T}}) where {T} = Ok(T)

trygeteltype(::Type{Dict}) where {K,V,Dict<:AbstractDict{K,V}} = eltype_impl(Dict)
trygeteltype(::Type{Num}) where {Num<:Number} = eltype_impl(Num)
trygeteltype(::Type{Str}) where {Str<:AbstractString} = eltype_impl(Str)

eltype_impl(::Type{T}) where {T} =
    if IteratorEltype(T) isa HasEltype
        return Ok(eltype(T))
    else
        return Causes.notimplemented(trygeteltype, (T,))
    end

@inline function trygetindex(a::AbstractArray, i...)::Result
    checkbounds(Bool, a, i...) || return Err(BoundsError(a, i))
    return Ok(@inbounds a[i...])
end

@inline function trysetindex!(a::AbstractArray, v, i...)::Result
    checkbounds(Bool, a, i...) || return Err(BoundsError(a, i))
    @inbounds a[i...] = v
    return Ok(v)
end

@inline function trygetindex(xs::Tuple, i::Integer)::Result
    i < 1 && return Err(BoundsError(xs, i))
    i > length(xs) && return Err(BoundsError(xs, i))
    return Ok(xs[i])
end

struct NotFound end

function trygetindex(dict::AbstractDict, key)::Result
    value = get(dict, key, NotFound())
    value isa NotFound && return Err(KeyError(key))
    return Ok(value)
end

function trysetindex!(dict::AbstractDict, value, key)::Result
    dict[key] = value
    return Ok(value)
end

trygetindex(dict::AbstractDict, k1, k2, ks...) = trygetindex(dict, (k1, k2, ks...))
trysetindex!(dict::AbstractDict, v, k1, k2, ks...) = trysetindex!(dict, v, (k1, k2, ks...))

trygetfirst(xs) = trygetindex(xs, 1)
trygetlast(xs) = trygetindex(xs, lastindex(xs))

trypop!(a::Vector)::Result = isempty(a) ? Causes.empty(a) : Ok(pop!(a))
trypopfirst!(a::Vector)::Result = isempty(a) ? Causes.empty(a) : Ok(popfirst!(a))

function trypush!(a::Vector, x)::Result
    y = @? tryconvert(eltype(a), x)
    push!(a, y)
    return Ok(a)
end

function trypushfirst!(a::Vector, x)::Result
    y = @? tryconvert(eltype(a), x)
    pushfirst!(a, y)
    return Ok(a)
end

function trytake!(ch::Channel)::Result
    y = iterate(ch)
    y === nothing && return Causes.empty(ch)
    return Ok(first(y))
end

function tryput!(ch::Channel, x)::Result
    isopen(ch) || return Causes.closed(ch)
    y = @? tryconvert(eltype(ch), x)
    try
        put!(ch, x)
    catch err
        err isa InvalidStateException && !isopen(ch) && return Causes.closed(ch)
        rethrow()
    end
    return Ok(y)
end


text:
using BioSeq
importall Base, BioSeq

immutable Alphabet{ T <: BioUnit }
  alphabet::Array{T,1}
  indice::Array{UInt8,1}
  case_sensitive::Bool
end

function alphabet{T <: BioUnit}(seq::Array{T,1}, case_sensitive::Bool)
  indices = zeros(UInt8, 256 )
  if case_sensitive
    seq = unique(seq)
    indices[seq] = 1:length(seq)
  else
    seq = unique(lowercase(seq))
    indices[seq] = 1:length(seq)
    seq = uppercase!(seq)
    indices[seq] = 1:length(seq)
  end

  Alphabet{T}(seq,indices,case_sensitive)
end

length{T <: BioUnit}(ab::Alphabet{T}) = length(ab.alphabet)

start{T <: BioUnit}(ab::Alphabet{T}) = 1
next{T <: BioUnit}(ab::Alphabet{T},i::Int) = (ab.alphabet[i],i+1)
done{T <: BioUnit}(ab::Alphabet{T},i::Int) = (i > length(ab))

join{T <: BioUnit}(ab::Alphabet{T}) = join(ab.alphabet)

function show{T <: BioUnit}(io::IO,ab::Alphabet{T})
  println(io,"Case $(ab.case_sensitive ? "S":"Ins")ensitive Alphabet{$(T)} of $(length(ab)) elements:\n")
  println(io," indice   : $(length(ab.indice))-element $UInt8 Array")
  println(io," alphabet : $(length(ab.alphabet))-element $T Array\n")
  println(" alphabet\t\t\tindice[alphabet]")
  println(io," $T ($Int)\t\tUInt8 ($Int)\n")
  for letter in ab
    println(io," $letter ($(int(letter))) \t\t\t$(ab.indice[letter]) ($(int(ab.indice[letter])))")
  end
end

getindex{T <: BioUnit, I <: Union{UInt8,Int}}(ab::Alphabet{T}, i::I ) = ab.alphabet[ i ]
getindex{T <: BioUnit, I <: Union{UInt8,Int}}(ab::Alphabet{T}, i::AbstractArray{I} ) = ab.alphabet[ i ]
getindex{T <: BioUnit}(ab::Alphabet{T}, i::T ) = ab.indice[ i ]
getindex{T <: BioUnit}(ab::Alphabet{T}, i::AbstractArray{T} ) = ab.indice[ i ]

in{T <: BioUnit}(element::T, ab::Alphabet{T}) = ab.indice[element] != 0

## Is In Alphabet ? ##

function isin{T,A<:Union{IntSet,Set,Alphabet}}(s::AbstractArray{T},y::A)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = in(s[i],y)
  end
  res
end

function isin{T,A<:Associative}(s::AbstractArray{T},y::A)
  res = BitArray(size(s))
  for i in 1:length(s)
    res[i] = haskey(y,s[i])
  end
  res
end

isin(s::ASCIIString,y)= isin(s.data,y)

## Check Alphabet ##

function check{T,A<:Union{IntSet,Set,Alphabet}}(s::AbstractArray{T},y::A)
  for x in s
    if !in(x,y)
      println("$(Char(x)) is not in this Alphabet")
      return(false)
    end
  end
  true
end

function check{T,A<:Associative}(s::AbstractArray{T},y::A)
  for x in s
    if !haskey(y,x)
      println("$(Char(x)) is not in this Alphabet")
      return(false)
    end
  end
  true
end

check(s::ASCIIString,y)= check(s.data,y)

## Interchange between Dic and Array values ##

function swap!{T<:Number}(s::AbstractArray{T},y::Dict{T,T})
  for i in 1:length(s)
    s[i] = get(y,s[i],s[i])
  end
  s
end

function swap{T,Tk,Tv}(s::AbstractArray{T},y::Dict{Tk,Tv})
  res = similar(s,Tv)
  for i in 1:length(s)
    res[i] = get(y,s[i],s[i])
  end
  res
end

swap!{T<:Integer,Tv}(s::AbstractArray{T},y::Array{Tv}) = (for i in length(s) s[i] = y[s[i]] end)
swap{T<:Integer,Tv}(s::AbstractArray{T},y::Array{Tv}) = [y[b] for b in s]

swap(s::ASCIIString,y)= ASCIIString(convert(Vector{UInt8},swap(s.data,y)))

## TODO: rand


text:
module CommutativeRings

using LinearAlgebra
using Base.Checked
using Primes

export category_trait, isfield
export Ring, RingInt, FractionRing, QuotientRing, Polynomial
export ZZ, QQ, ZZmod, Frac, Quotient, UnivariatePolynomial, MultivariatePolynomial
export GaloisField

export Hom, Ideal

export isunit, deg, content, primpart, isnegative
export LC, LM, LT, lcunit, multideg, modulus, value
export isdiv, pdivrem, pgcd, pgcdx, resultant, discriminant
export basetype, basetypes, depth, iszerodiv
export monom, ismonom, ismonic, issimpler, iscoprime, evaluate, derive
export isirreducible, irreducible, irreducibles, monic, factorise
export num_irreducibles, isreducible, reducible, reducibles
export characteristic, dimension, order
export ofindex, index
export log_zech, generator
export GF, homomorphism
export num_primitives, isprimitive

export VectorSpace, complement, sum, intersect, isequal, issubset
export groebnerbase, SPOL, lextend
export generators, varnames, varname, factors

export characteristic_polynomial, adjugate, companion

export coeffbounds

import Base: +, -, *, /, inv, ^, \, //, getindex, sign, log
import Base: iszero, isone, zero, one, div, rem, divrem, ==, hash, gcd, gcdx, lcm
import Base: copy, show, promote_rule, convert, abs, isless
import Primes: factor
import Base: numerator, denominator
import LinearAlgebra: checksquare, det

# RingClass subtypes describe the different categories
abstract type RingClass end
struct ZZClass{T<:Integer} <: RingClass end
abstract type FractionRingClass <:RingClass end
struct FractionClass{P} <: FractionRingClass end
struct QQClass{S<:Integer} <:FractionRingClass end
abstract type QuotientRingClass <:RingClass end
struct QuotientClass{Id,M} <: QuotientRingClass
    modulus::M
end
struct ZZmodClass{m,T<:Integer} <: QuotientRingClass
    modulus::T
end
abstract type PolyRingClass <: RingClass end
struct UniPolyRingClass{X,R} <: PolyRingClass end
struct MultiPolyRingClass{X,R,N} <: PolyRingClass
    varnames::Vector{Symbol}
end

struct GaloisFieldClass{Id,T,Q} <: QuotientRingClass
    exptable::Vector{T}
    zechtable::Vector{T}
end

const NCT = Val{:nocheck}
const NOCHECK = Val(:nocheck)

# Ring subtypes describe the ring elements
# They are connected to a corresponding RingClass by type parameters
# and accessible by method `owner`.
"""
    Ring{<:RingClass}

Abstract superclass of all commutative ring element types.
All subtypes support the arithmetic operators `+, -, *, and ^(power)`.
The operators `inv, /` may be defined on a subset.
"""
abstract type Ring{T<:RingClass} end

"""
    RingInt

Union of system `Integer` types and any `Ring` subtype.
"""
const RingInt = Union{Ring,Integer}

"""
    FractionRing{S<:RingInt,T<:FractionRingClass}

Rational extension of elements of ring `S`. Example the field of rational numbers.
"""
abstract type FractionRing{S<:RingInt,T<:FractionRingClass} <: Ring{T} end

"""
    QuotientRing{S<:Union{Integer,Ring},T<:QuotientRingClass}

Quotient ring of ring `S` by some ideal. If S = Z, the ring of integer numbers, and p
a positive number Z/p is calculation modulo p.
"""
abstract type QuotientRing{S<:RingInt,T<:QuotientRingClass} <: Ring{T} end

"""
    Polynomial{S<:Ring,T<:PolyRingClass}
"""
abstract type Polynomial{S<:Ring,T<:PolyRingClass} <: Ring{T} end

"""
    ZZ{S<:Signed}

The ring of integer numbers. There may be restrictions on the set of representable
elemets according to the chosen `S`.
"""
struct ZZ{T<:Signed} <: Ring{ZZClass{T}}
    val::T
    ZZ{T}(val::Integer) where T = new{T}(val)
    ZZ(val::T) where T<:Signed = ZZ{T}(val)
end

"""
    ZZmod{m,S<:Integer}

Quotient ring modulo integer `m > 0`, also noted as `ZZ/m`.
If `p` is a prime number, `ZZmod{p}` is the field `ZZ/p`.
"""
struct ZZmod{m,S<:Integer} <: QuotientRing{S,ZZmodClass{m,S}}
    val::S
    ZZmod{m,T}(a::Integer, ::NCT) where {m,T} = new{m,T}(T(a))
end

"""
    Frac{R}

The ring of fractions of `R`. The elements consist of pairs `num::R,den::R`.
During creation the values may be canceled to achieve `gcd(num, den) == one(R)`.
The special case of `R<:Integer` is handled by `QQ{R}`.
"""
struct Frac{P<:Union{Polynomial,ZZ}} <: FractionRing{P,FractionClass{P}}
    num::P
    den::P
    Frac{P}(num::P, den::P, ::NCT) where P = new{P}(num, den)
end

"""
    Quotient{R,I,m}

The quotient ring of `R` modulo `m` of type `I`, also written as `R / m`.
`m` may be an ideal of `R` or a (list of) element(s) of `R` generating the ideal.
Typically `m` is replaced by a symbolic `X`, and the actual `m` is given as argument(s)
to the type constructor like  `new_class(Quotient{ZZ,`X`}, m...)`.
If the ``X``is omitted, an anonymous symbol is used.

The preferred way of construction is via `Zm = Z/m`.
"""
struct Quotient{R<:Ring,I,X,Id} <: QuotientRing{R,QuotientClass{X,I}}
    val::R
    Quotient{R,I,X,Id}(v::R, ::NCT) where {I,X,R<:Ring,Id} = new{R,I,X,Id}(v)
end

"""
    QQ{S<:Integer}

"""
struct QQ{S<:Integer} <: FractionRing{S,QQClass{S}}
    num::S
    den::S
    QQ{T}(num::Integer, den::Integer, ::NCT) where T = new{T}(num, den)
end

"""
    UnivariatePolynomial{S<:RingInt,Id}

Polynomials of ring elemets `S` in one variable `Id` (by default `:x`).
The variable name is specified as a `Symbol`.

A convenience constructor `S[:x]` is the preferred way to construct this class.
"""
struct UnivariatePolynomial{S<:Ring,X} <: Polynomial{S,UniPolyRingClass{X,S}}
    coeff::Vector{S}
    UnivariatePolynomial{S,X}(v::Vector{S}, ::NCT) where {X,S<:Ring} = new{S,X}(v)
end

"""
    MultivariatePolynomial{S<:RingInt,N,Id}

Polynomials of ring elemets `S` in `N` variables.
The `Id` identifies on object of type `MultiPolyRingClass` which is needed to store
the variable names and properties.

A convenience constructor `S[:x,:y...]` is the preferred way to construct this class.
"""
struct MultivariatePolynomial{S<:Ring,N,Id,T,B} <: Polynomial{S,MultiPolyRingClass{Id,S,N}}
    ind::Vector{T}
    coeff::Vector{S}
end

struct GaloisField{Id,T,Q} <: QuotientRing{ZZmod{T},GaloisFieldClass{Id,T,Q}}
    val::T
    GaloisField{Id,T,Q}(v::Integer, ::NCT) where {Id,T,Q} = new{Id,T,Q}(T(v))
end

# Categorial traits specify algebraic properties of ring types
# (cf. https://en.wikipedia.org/wiki/Integral_domain)
# rings ⊃ commutative rings ⊃ integral domains ⊃ integrally closed domains ⊃ GCD domains ⊃
# unique factorization domains ⊃ principal ideal domains ⊃ Euclidean domains ⊃ fields ⊃ algebraically closed fields
abstract type RingTrait end
abstract type CommutativeRingTrait <: RingTrait end
abstract type IntegralDomainTrait <: CommutativeRingTrait end
abstract type IntegrallyClosedDomainTrait <: IntegralDomainTrait end
abstract type GCDDomainTrait <: IntegrallyClosedDomainTrait end
abstract type UniqueFactorizationDomainTrait <: GCDDomainTrait end
abstract type PrincipalIdealDomainTrait <: UniqueFactorizationDomainTrait end
abstract type EuclidianDomainTrait <: PrincipalIdealDomainTrait end
abstract type FieldTrait <: EuclidianDomainTrait end
abstract type AlgebraicallyClosedFieldTrait <: FieldTrait end

## End of Ring classes

"""
    Ideal{R}

Respresent in Ideal of Ring `R`. Only Ideals with a finite (Groebner) basis.
"""
struct Ideal{R<:Ring}
    base::Vector{R}
    Ideal{R}(a::Vector{R}) where R = new{R}(a)
end

"""
    Hom{function,R,R'}

Represent a ring homomorphism `function:R -> R'`.
"""
struct Hom{F,R<:RingInt,S<:RingInt}
    f::F
    Hom{R,S}(f::F) where {F<:Union{Function,Type},R,S} = new{F,R,S}(f)
end

struct VectorSpace{T}
    base::T
    pivr::Vector{Int} # row permutation vector
end

# implementation

include("typevars.jl")
include("generic.jl")
include("zz.jl")
include("qq.jl")
include("zzmod.jl")
include("quotient.jl")
include("fraction.jl")
include("univarpolynom.jl")
include("multivarpolynom.jl")
include("ideal.jl")
include("enumerations.jl")
include("factorization.jl")
include("intfactorization.jl")
include("numbertheoretical.jl")
include("galoisfields.jl")
include("linearalgebra.jl")

end # module


text:
@safetestset "JS rendering" begin

  @safetestset "Plain JS rendering" begin
    using Genie
    using Genie.Renderer.Js

    script = raw"var app = new Vue({el: '#app', data: { message: 'Hello Vue!' }})"

    r = js(script)

    @test String(r.body) == "var app = new Vue({el: '#app', data: { message: 'Hello Vue!' }})"
    @test r.headers[1]["Content-Type"] == "application/javascript; charset=utf-8"
  end;

  @safetestset "Vars JS rendering" begin
    using Genie
    using Genie.Renderer.Js
    using Genie.Renderer.Json.JSONParser

    data = JSON.json(("message" => "Hi Vue"))

    script = raw"var app = new Vue({el: '#app', data: $data})"

    r = js(script, data = data)

    @test String(r.body) == "var app = new Vue({el: '#app', data: {\"message\":\"Hi Vue\"}})\n"
    @test r.headers[1]["Content-Type"] == "application/javascript; charset=utf-8"
  end;

end;


text:
function readBaseNodes( mpSim::MPsim, xf::XF, catXF::XF )

    if !XLSX.hassheet( xf, "Base Nodes" )
        @warn "Configuration file is missing a 'Base Nodes' sheet, cannot configure base nodes."
        return
    end  # if !XLSX.hassheet( xf, "Base Nodes" )

    ws = xf["Base Nodes"]
    nNodes = ws["B4"] isa Integer ? ws["B4"] : 0
    lineNrs = (1:nNodes) .+ 6

    if XLSX.hassheet( catXF, "States" )
        catalogue = catXF["States"]
        nodes = readBaseNode.( Ref( mpSim ), Ref( ws ), Ref( catalogue ),
            lineNrs )
    else
        @warn "Catalogue file is missing a 'States' sheet, cannot configure attributes defining base nodes."
        nodes = readBaseNode.( Ref( ws ), lineNrs )
    end  # if XLSX.hassheet( catXF, "States" )

    setSimulationBaseNodes!( mpSim, nodes )

    # Set the preferred node order.
    nodes = ws[CR( 7, 1, 6 + nNodes, 1 )]
    order = ws[CR( 7, 4, 6 + nNodes, 4 )]
    isValid = .!isa.( nodes, Missing ) .& isa.( order, Integer )

    if !any( isValid )
        return
    end  # if !any( isValid )

    nodes = string.( nodes[isValid] )
    order = Vector{Int}( order[isValid] )
    setSimulationBaseNodeOrder!( mpSim, nodes, order )

end  # readBaseNodes( mpSim, wf, catWF )

function readBaseNode( mpSim::MPsim, ws::WS, catalogue::WS, sLine::Int )

    # Read the base node's target.
    node = readBaseNode( ws, sLine )

    # If it's in the catalogue, read the base node's attribute requirements and
    #   associated attrition scheme.
    cLine = ws[sLine, 3]

    if !( cLine isa Integer )
        return node
    end  # if !( cLine isa Integer )

    cLine += 1

    # Set the node's attrition scheme.
    isAttritionFixed = catalogue[cLine, 3] isa Missing ? "YES" :
        string( catalogue[cLine, 3] )
    isAttritionFixed = isAttritionFixed != "NO"

    if isAttritionFixed
        attrition = Attrition( string( node.name, " attrition" ) )
        period = catalogue[cLine, 4] isa Real ? catalogue[cLine, 4] : 12.0
        rate = catalogue[cLine, 5] isa Real ? catalogue[cLine, 5] : 0.0
        setAttritionPeriod!( attrition, period )
        setAttritionRate!( attrition, rate )
        addSimulationAttrition!( mpSim, attrition )
        setNodeAttritionScheme!( node, attrition.name )
    else
        nodeAttrition = catalogue[cLine, 6] isa Missing ? "default" :
            string( catalogue[cLine, 6] )
        setNodeAttritionScheme!( node, nodeAttrition )
    end  # if isAttritionFixed

    # Set the node's attributes.
    nRequirements = catalogue[cLine, 7] isa Real ?
        floor( Int, catalogue[cLine, 7] ) : 0
    requirements = Dict{String, String}(
        string( catalogue[cLine, 6 + 2 * ii] ) =>
        string( catalogue[cLine, 7 + 2 * ii] ) for ii in 1:nRequirements )
    setNodeRequirements!( node, requirements )

    return node

end  # readBaseNode( ws, catalogue, sLine )

function readBaseNode( ws::WS, sLine::Int )

    node = BaseNode( string( ws[sLine, 1] ) )
    nodeTarget = ws[sLine, 2]

    if nodeTarget isa Integer
        setNodeTarget!( node, nodeTarget )
    end  # if nodeTarget isa Integer

    return node

end  # function readBaseNode( ws, sLine )

function readBaseNodes( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    nodePars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Base Node'" ) ) )

    if !isempty( nodePars )
        nodes = map( eachindex( nodePars[:, :parName] ) ) do ii
            node = BaseNode( nodePars[ii, :parName] )
            pars = split( nodePars[ii, :parValue], ";" )

            if tryparse( Int, pars[1] ) isa Int
                setNodeTarget!( node, parse( Int, pars[1] ) )
            end  # if tryparse( Int, pars[1] )

            setNodeAttritionScheme!( node, string( pars[2] ) )

            if length( pars[3] ) > 2
                nodeReqs = split.( split( pars[3][2:(end-1)], "," ), ":" )
                nodeReqDict = Dict{String, String}()

                for nodeReq in nodeReqs
                    nodeReqDict[nodeReq[1]] = nodeReq[2]
                end  # for nodeReq in nodeReqs

                setNodeRequirements!( node, nodeReqDict )
            end  # if length( pars[3] ) > 2

            return node
        end  # map( eachindex( nodePars[:parName] ) ) do ii

        setSimulationBaseNodes!( mpSim, nodes )
    end  # if !isempty( nodePars )

    nodeOrder = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Base Node Order'" ) ) )

    if !isempty( nodeOrder ) && ( length( nodeOrder[1, :parValue] ) > 2 )
        order = split.( split( nodeOrder[1, :parValue], ";" ), ":" )
        orderDict = Dict{String, Int}()

        for orderPair in order
            orderDict[orderPair[1]] = parse( Int, orderPair[2] )
        end  # for orderPair in order

        setSimulationBaseNodeOrder!( mpSim, orderDict )
    end  # if !isempty( nodeOrder ) && ...

end  # readBaseNodes( mpSim, configDB, configName )


text:
## This shows the simplest use of BAMF with a Gaussian PSF

using Revise
using BayesMultiFit
BAMF=BayesMultiFit
using ReversibleJumpMCMC
const RJMCMC = ReversibleJumpMCMC
using Plots
using Distributions
using MicroscopePSFs
const PSF=MicroscopePSFs

plot_flag = true

## simulation config
n=6             # number of emitters
μ=1000          # mean photons per emitter
sz=16           # ROI size in pixels
iterations=5000 # RJMCMC iterations
burnin=5000     # RJMCMC iterations for burn-in

## PSF config
σ=1.3 # Gaussian PSF Sigma in Pixels
pixelsize=1.0
psf=PSF.Gauss2D(σ,pixelsize)

## setup prior distribution on intensity
α=4.0
θ=μ/α
prior_photons=Gamma(α,θ)
pdf_x=range(1,quantile(prior_photons,.999),length=1024)
mypdf=pdf(prior_photons,pdf_x)
if plot_flag
    plt=plot(pdf_x,mypdf)
    display(plt)
end

## setup prior on background
prior_background=Gamma(1.0,1.0)

# set emitter positions and intensity
x=sz*rand(Float32,n)
y=sz*rand(Float32,n)

photons=Float32.(rand(prior_photons,n))
bg=1f-6
datastate=BAMF.StateFlatBg(n,x,y,photons,bg)

## Create synthetic data
data=BAMF.ArrayDD(sz)
BAMF.genmodel!(datastate,psf,data)
BAMF.poissrnd!(data.data)
# imshow(data.data)  #look at data

## run chain. This is the call to the main algorithm
# this example doesn't update background
mychain=BAMF.bamf_roi(data,psf,prior_photons,prior_background; jumpprobability=[0.5,0,.1,.1,.1,.1])

## Display
if plot_flag
    plotly()
    zm=4
    plt=BAMF.histogram2D(mychain.states,sz,zm,datastate)
    display(plt)
end

map_n,posterior_n,traj_n=BAMF.getn(mychain.states)
if plot_flag
    plt2=plot(traj_n)
    display(plt2)
    out=BAMF.showoverlay(mychain.states,data,psf)
    plt=heatmap(out[:,:,end]')
    display(plt)
end

## MAPN Results
states_mapn,n=BAMF.getmapnstates(mychain.states)
if plot_flag
    plt=BAMF.histogram2D(states_mapn,sz,zm,datastate)
    display(plt)
end

Results_mapn=BAMF.getmapn(mychain.states)
if plot_flag
    gr()
    plt=BAMF.plotstate(datastate,Results_mapn;fig=plot())
    display(plt)
end


text:
module Node2Vec
    using LightGraphs,SimpleWeightedGraphs,Random,Word2Vec,DelimitedFiles

    export node2vec_walk, simulate_walks, learn_embeddings

    include("rw.jl")
end # module


text:
using Test
import ArchGDAL as AG

@testset "test_featurelayer.jl" begin
    @testset "Testing FeatureLayer Methods" begin
        AG.read("data/point.geojson") do dataset
            AG.copy(dataset) do tmpcopy
                @test AG.nlayer(tmpcopy) == 1
                tmplayer = AG.getlayer(tmpcopy, 0)
                @test sprint(print, AG.getspatialref(tmplayer)) ==
                      "NULL Spatial Reference System"
                AG.getspatialref(tmplayer) do spref
                    @test sprint(print, spref) ==
                          "NULL Spatial Reference System"
                end
                @test AG.isignored(AG.getgeomdefn(AG.layerdefn(tmplayer), 0)) ==
                      false
                AG.setignoredfields!(tmplayer, ["OGR_GEOMETRY"])
                @test AG.isignored(AG.getgeomdefn(AG.layerdefn(tmplayer), 0)) ==
                      true
            end

            AG.copy(dataset, driver = AG.getdriver("Memory")) do tmpcopy
                tmplayer = AG.getlayer(dataset, 0)
                AG.copy(tmplayer) do copylayer
                    @test AG.nfeature(tmplayer) == AG.nfeature(copylayer)
                end
                @test sprint(print, AG.getspatialref(tmplayer)) ==
                      "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
                AG.getspatialref(tmplayer) do spref
                    @test sprint(print, spref) ==
                          "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
                end
                newlayer = AG.createlayer(
                    name = "new layer",
                    dataset = tmpcopy,
                    spatialref = AG.getspatialref(tmplayer),
                    geom = AG.wkbPoint,
                )
                @test AG.ngeom(AG.layerdefn(newlayer)) == 1
                @test sprint(print, newlayer) == """
                    Layer: new layer
                      Geometry 0 (): [wkbPoint]
                    """
                AG.writegeomdefn(newlayer, "new geom", AG.wkbLineString) do gfd
                    @test AG.getname(gfd) == "new geom"
                    @test AG.gettype(gfd) == AG.wkbLineString
                end
                @test sprint(print, newlayer) == """
                    Layer: new layer
                      Geometry 0 (): [wkbPoint]
                      Geometry 1 (new geom): [wkbLineString]
                    """
                @test AG.ngeom(AG.layerdefn(newlayer)) == 2
                @test AG.nfeature(newlayer) == 0
                AG.createfeature(newlayer) do newfeature
                    AG.setgeom!(newfeature, 0, AG.createpoint())
                    AG.setgeom!(newfeature, 1, AG.createlinestring())
                    return nothing
                end
                @test AG.nfeature(newlayer) == 1
                @test sprint(print, newlayer) == """
                    Layer: new layer
                      Geometry 0 (): [wkbPoint]
                      Geometry 1 (new geom): [wkbLineString]
                    """
                AG.deletefeature!(newlayer, 0)
                @test AG.nfeature(newlayer) == 0
                @test sprint(print, newlayer) == """
                    Layer: new layer
                      Geometry 0 (): [wkbPoint]
                      Geometry 1 (new geom): [wkbLineString]
                    """
                AG.writegeomdefn!(newlayer, "new poly", AG.wkbPolygon)
                @test AG.ngeom(AG.layerdefn(newlayer)) == 3
                @test sprint(print, newlayer) == """
                    Layer: new layer
                      Geometry 0 (): [wkbPoint]
                      Geometry 1 (new geom): [wkbLineString]
                      Geometry 2 (new poly): [wkbPolygon]
                    """

                AG.addfeature(newlayer) do newfeature
                    AG.setgeom!(newfeature, 0, AG.createpoint())
                    AG.setgeom!(newfeature, 1, AG.createlinestring())
                    AG.setgeom!(
                        newfeature,
                        2,
                        AG.createpolygon([[
                            [0.0, 0.0],
                            [1.0, 1.0],
                            [0.0, 1.0],
                        ]]),
                    )

                    @test sprint(print, AG.getgeom(newfeature)) ==
                          "Geometry: POINT EMPTY"
                    @test sprint(print, AG.getgeom(newfeature, 0)) ==
                          "Geometry: POINT EMPTY"
                    @test sprint(print, AG.getgeom(newfeature, 1)) ==
                          "Geometry: LINESTRING EMPTY"
                    @test sprint(print, AG.getgeom(newfeature, 2)) ==
                          "Geometry: POLYGON ((0 0,1 1,0 1))"
                    AG.getgeom(newfeature) do g
                        @test sprint(print, g) == "Geometry: POINT EMPTY"
                    end
                    AG.getgeom(newfeature, 1) do g
                        @test sprint(print, g) == "Geometry: LINESTRING EMPTY"
                    end
                    AG.getgeom(newfeature, 2) do g
                        @test sprint(print, g) ==
                              "Geometry: POLYGON ((0 0,1 1,0 1))"
                    end
                end
            end

            layer = AG.getlayer(dataset, 0)
            @test sprint(print, AG.getspatialref(layer)) ==
                  "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
            @test AG.fidcolumnname(layer) == ""
            @test AG.geomcolumnname(layer) == ""
            @test AG.nreference(layer) == 0
            AG.reference(layer)
            @test AG.nreference(layer) == 1
            AG.dereference(layer)
            @test AG.nreference(layer) == 0

            @test AG.nfeature(layer) == 4
            @test AG.getfield.(layer, 1) == ["point-a", "point-b", "a", "b"]
            AG.setspatialfilter!(layer, 100, -1, 100.1, 1)
            @test AG.toWKT(AG.getspatialfilter(layer)) ==
                  "POLYGON ((100 -1,100 1,100.1 1.0,100.1 -1,100 -1))"
            @test AG.nfeature(layer) == -1
            @test AG.getfield.(layer, 1) == ["point-a", "a"]
            AG.clearspatialfilter!(layer)

            @test AG.nfeature(layer) == 4
            @test AG.getfield.(layer, 1) == ["point-a", "point-b", "a", "b"]
            @test AG.findgeomindex(AG.layerdefn(layer)) == 0
            AG.setspatialfilter!(layer, 0, 100, -1, 100.1, 1)
            @test AG.toWKT(AG.getspatialfilter(layer)) ==
                  "POLYGON ((100 -1,100 1,100.1 1.0,100.1 -1,100 -1))"
            @test AG.nfeature(layer) == -1
            @test AG.getfield.(layer, 1) == ["point-a", "a"]

            AG.clone(AG.getspatialfilter(layer)) do poly
                n = 0
                for feature in layer
                    n += 1
                end
                @test n == 2
                AG.clearspatialfilter!(layer)
                @test sprint(print, AG.getspatialfilter(layer)) ==
                      "NULL Geometry"
                n = 0
                for feature in layer
                    n += 1
                end
                @test n == 4

                @testset "Test with setting to index of geomfield" begin
                    AG.setspatialfilter!(layer, 0, poly)
                    n = 0
                    for feature in layer
                        n += 1
                    end
                    @test n == 2
                    AG.clearspatialfilter!(layer, 0)
                    n = 0
                    for feature in layer
                        n += 1
                    end
                    @test n == 4

                    AG.setattributefilter!(layer, "FID = 2")
                    n = 0
                    for feature in layer
                        n += 1
                    end
                    @test n == 1
                    AG.setattributefilter!(layer, "FID = 3")
                    n = 0
                    for feature in layer
                        n += 1
                    end
                    @test n == 2
                    AG.clearattributefilter!(layer)
                    n = 0
                    for feature in layer
                        n += 1
                    end
                    @test n == 4
                    AG.nextfeature(layer) do feature
                        @test AG.getfield(feature, 1) == "point-a"
                    end
                    AG.setnextbyindex!(layer, 2)
                    AG.nextfeature(layer) do feature
                        @test AG.getfield(feature, 1) == "a"
                    end
                    @test AG.testcapability(layer, "OLCRandomWrite") == false
                end
            end
            @test AG.findfieldindex(layer, "FID", true) == 0
            @test AG.findfieldindex(layer, :FID, false) == 0
            @test AG.findfieldindex(layer, "pointname", true) == 1
            @test AG.findfieldindex(layer, "pointname", false) == 1
            @test AG.findfieldindex(layer, "geom", true) == -1
            @test AG.findfieldindex(layer, "geom", true) == -1
            @test AG.findfieldindex(layer, "rubbish", true) == -1
            @test AG.findfieldindex(layer, "rubbish", false) == -1
            @test sprint(print, AG.envelope(layer, 0, true)) ==
                  "GDAL.OGREnvelope(100.0, 100.2785, 0.0, 0.0893)"
            @test sprint(print, AG.envelope(layer, true)) ==
                  "GDAL.OGREnvelope(100.0, 100.2785, 0.0, 0.0893)"
        end
    end

    # Not implemented yet:

    # intersection(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)
    # union(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressdata = C_NULL, progressfunc::Function = GDAL.gdaldummyprogress)
    # symdifference(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)
    # identity(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)
    # update(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)
    # clip(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)
    # erase(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.gdaldummyprogress, progressdata = C_NULL)

    # deletefielddefn!(layer::FeatureLayer, i::Integer)
    # reorderfielddefn!(layer::FeatureLayer, indices::Vector{Cint})
    # reorderfielddefn!(layer::FeatureLayer, oldpos::Integer, newpos::Integer)
    # updatefielddefn!(layer::FeatureLayer, i::Integer, newfielddefn::FieldDefn, flags::UInt8)

    # starttransaction(layer)
    # committransaction(layer)
    # rollbacktransaction(layer)

end


text:
module MyUtils

include("defunpack.jl")
include("my_threads.jl")
include("my_distributed.jl")
include("printf_functions.jl")
include("remove_one_from.jl")
include("showimg.jl")
include("savevar.jl")

using .PrintfFunctions
export printf, sprintf

using .DefUnPack
export @defunpack

end


text:
using Logging
using ArgParse

@info "# Parsing arguments..."
 function parse_commandline()
     s = ArgParseSettings( description = "Combined NEID solar daily reports into one csv file.")
     add_arg_group!(s, "Files", :argg_files)
     @add_arg_table! s begin
         "input_path"
             help = "Path for input files"
             arg_type = String
             default = pwd()
         "output"
             help = "Filename for output (csv)"
             arg_type = String
             default = "summary.csv"
         "--input_filename"
             help = "Filename for daily tables (csv)"
             arg_type = String
             default = "daily_summary.toml"
         "--exclude_filename"
             help = "Filename with dates to exclude (csv)"
             arg_type = String
         "--overwrite"
            help = "Specify it's ok to overwrite the output file."
            #default = true
            action = :store_true
         "--log_info"
            help = "Print info-level messages."
            action = :store_true
         "--log_debug"
            help = "Print debug-level messages."
            action = :store_true
      end

     return parse_args(s)
 end
 args = parse_commandline()
 if args["log_debug"]
    global_logger(ConsoleLogger(stderr, Logging.Debug))
 elseif args["log_info"]
    global_logger(ConsoleLogger(stderr, Logging.Info))
 end

 input_fn = args["input_filename"]
 input_path = args["input_path"]
 output_fn = joinpath(input_path,args["output"])

if !(args["overwrite"] || !isfile(output_fn) || (filesize(output_fn)==0))
   @error "Can't overwrite " output_filename=output_fn
   exit(1)
else
   #touch(output_fn)  # Should we create empty file as a lock?
end

@info "# Loading packages"
using Dates, TOML
using CSV, DataFrames, Query, Glob

@info "# Finding files named $input_fn"
files1 = glob([r"\d{2}",args["input_filename"]],args["input_path"])
files2 = glob([r"\d{2}",r"\d{2}",args["input_filename"]],args["input_path"])
files3 = glob([r"\d{4}",r"\d{2}",r"\d{2}",args["input_filename"]],args["input_path"])
files = vcat(files1,files2,files3)
@assert length(files) >= 1

if !isnothing(args["exclude_filename"]) && isfile(args["exclude_filename"]) && (filesize(args["exclude_filename"])>0)
   @info "# Reading days to exclude."
   df_exclude = CSV.read(args["exclude_filename"],DataFrame)
   @info "# Found $(size(df_exclude,1)) dates to exclude"
else
   df_exclude = DataFrame()
end

@info "# Parsing daily csv files."  num_files=length(files)
daily = Vector{DataFrame}(undef, size(files,1) )
println("files = ", files)
flush(stdout)
flush(stderr)
j = 0
for file in files
    if filesize(file) >0
       #=
       fn_toml = replace(file,"csv"=>"toml")
       if filesize(fn_toml)>0
          @info "# Processing $fn_toml"
          t = TOML.parsefile(fn_toml)

          if haskey(t,"obs_date") && haskey(t["obs_date"],"string") && (size(df_exclude,1)>=1) && ( in(t["obs_date"]["string"], df_exclude.date_to_exclude) )
             continue
          end
       end
       =#
       @info "# Processing $file"
       d = CSV.read(file,DataFrame)
       if hasproperty(df_exclude,:date_to_exclude) && in(Date(julian2datetime(first(d.jd_drp))), df_exclude.date_to_exclude) continue end
       global j += 1
       daily[j] = d
    end
end
num_days_with_usable_obs = j
resize!(daily,num_days_with_usable_obs)

@info "# Making dataframe"
df = DataFrame()
for j in 1:num_days_with_usable_obs
  append!(df,daily[j],cols=:union)
end

df_sorted = sort!(df,"jd_drp")

#=
daily_flat = flatten_dict.(daily)
first_cols = ["obs_date.string", "obs_date.mean_bjd","num_rvs.good"]
cols_to_exclude = ["report.hostname", "title", "report.input_md5", "report.input_file"]
df =  array_of_dict_to_dataframe( daily_flat, first_cols=first_cols, exclude_cols=cols_to_exclude )

@info "# Sorting dataframe"
df_sorted = sort!(df,Symbol("obs_date.mean_bjd") )
=#

#=
df.obs_date = map(day->day["obs_date"]["string"], daily)
df.obs_date = map(day->day["obs_date"]["string"], daily)
df.mean_bjd = map(day->day["obs_date"]["mean_bjd"], daily)
df.num_rvs_usable = map(day->day["num_rvs"]["usable"], daily)
df.num_rvs_good = map(day->day["num_rvs"]["good"], daily)
df.mean_rv_drp = map(day->day["rv"]["drp"]["mean_rv"], daily)
df.median_rv_drp = map(day->day["rv"]["drp"]["mean_rv"], daily)
df.median_σ_rv = map(day->day["rv"]["drp"]["median_σ_rv"], daily)
df.rms_rv_drp = map(day->day["rv"]["drp"]["rms_rvs"], daily)
df.winsor_mean_rv_drp = map(day->day["rv"]["drp"]["winsor_mean_rv"], daily)
df.winsor_rms_rv_drp = map(day->day["rv"]["drp"]["winsor_rms_rv"], daily)
@info "# Sorting dataframe"
df_sorted = df |> @orderby( _.mean_bjd ) |> DataFrame
=#

@info "# Writing CSV file"
CSV.write(output_fn,df_sorted )


text:
# utility functions for Gibbs sampling

mutable struct GibbsSamplerControls
  """
    Structure to store the parameters of the Gibbs sampler
  """
  n_iter::Int64 # total number of Gibbs iterations to run
  n_burn::Int64  # burn-in
  n_thin::Int64 #  collect every `n_thin` samples
  n_print::Int64 # print status updates with timings every `n_print` iterations

  # set which components to update
  gibbs_psi::Bool # NTL Ψ paramters
  gibbs_alpha::Bool  # NTL alpha parameter
  gibbs_T::Bool # arrival times
  gibbs_ia_params::Bool # arrival time distribution parameters
  gibbs_sigma::Bool # order of blocks in partition/vertices in graph

  # set arrival model
  arrivals::String
  ia_distribution::Function
  # prior parameters
  ia_prior_params::Vector{Float64}
  n_ia_params::Int64

  # model-specific functions used within the sampler
  initialize_arrival_params::Function # for initializing arrival distn parameters
  update_arrival_params::Function # for updating arrival distn parameters
  ntl_alpha_log_prior::Function # for calculating the prior density of alpha
  ntl_alpha_prior_dist::UnivariateDistribution # for initializing alpha
  lpdf_ia_param_prior::Function # log pdf for interarrival parameter priors

  # slice-sampling parameters
  ss_w_alpha::Float64 # slice-sampling window parameter for alpha

  alpha_fixed::Float64 # when alpha is fixed (i.e., not to be updated)
  arrival_params_fixed::Vector{Float64} # when arrival distn parameters are fixed

end

struct GibbsSamplerOutput
 """
  Structure to store posterior samples
 """
  psi::Array{Float64,2}
  T::Array{Int64,2}
  alpha::Array{Float64,1}
  ia_params::Array{Float64}
  sigma::Array{Int64,2}
  log_joint::Vector{Float64}
  t_elapsed::Float64

end

function GibbsSampler(G::GibbsSamplerControls,PP::Vector{Int},T_data::Vector{Int})
  """
    Main function for running the Gibbs sampler.

    Inputs:
      - `G`: a GibbsSamplerControls structure
      - `PP`: a partition or degree vector
      - `T_data`: a vector of arrival times (pass an empty vector if unknown)
  """

  n_iter = G.n_iter
  n_burn = G.n_burn
  n_thin = G.n_thin

  gibbs_psi = G.gibbs_psi
  gibbs_T = G.gibbs_T
  gibbs_alpha = G.gibbs_alpha
  gibbs_ia_params = G.gibbs_ia_params
  gibbs_sigma = G.gibbs_sigma

  arrivals = G.arrivals
  ia_dist = G.ia_distribution
  ia_prior_params = G.ia_prior_params
  n_ia_params = G.n_ia_params

  w_alpha = G.ss_w_alpha

  initialize_arrival_params = G.initialize_arrival_params
  update_arrival_params! = G.update_arrival_params
  ntl_alpha_log_prior = G.ntl_alpha_log_prior
  ntl_alpha_prior_dist = G.ntl_alpha_prior_dist

  lpdf_ia_param_prior = G.lpdf_ia_param_prior

  alpha_fixed = G.alpha_fixed
  arrival_params_fixed = G.arrival_params_fixed

  K = size(PP,1)
  N = sum(PP)

  # println("Initializing sampler.")
  gibbs_psi ? psi_gibbs = zeros(Float64,K,Int(ceil((n_iter-n_burn)/n_thin))) : psi_gibbs = []
  gibbs_T ? T_gibbs = zeros(Int,K,Int(ceil((n_iter-n_burn)/n_thin))) : T_gibbs = []
  gibbs_alpha ? alpha_gibbs = zeros(Float64,Int(ceil((n_iter-n_burn)/n_thin))) : alpha_gibbs = []
  gibbs_ia_params ? ia_params_gibbs = zeros(Float64,n_ia_params,Int(ceil((n_iter-n_burn)/n_thin))) : ia_params_gibbs = []
  gibbs_sigma ? perm_gibbs = zeros(Int,K,Int(ceil((n_iter-n_burn)/n_thin))) : perm_gibbs = []
  log_joint_gibbs = zeros(Float64,Int(ceil((n_iter-n_burn)/n_thin)))

  # initialize
  gibbs_psi ? psi_current = 0.5*ones(Float64,K) : nothing
  gibbs_alpha ? alpha_current = initialize_alpha(ntl_alpha_prior_dist) : alpha_current = [alpha_fixed]
  if gibbs_ia_params
    arrival_params_current = initialize_arrival_params(ia_prior_params)
  else
    arrival_params_current = arrival_params_fixed
  end
  if gibbs_T
    T_current = initialize_arrival_times(PP,alpha_current[1],Geometric((K-1)/(N-1)))
  elseif isempty(T_data)
    error("Empty `T_data` is not compatible with not inferring arrival times.")
  else
    T_current = T_data
  end
  perm_current = collect(1:K) # initial permutation order

  # Gibbs sampler
  n_print = G.n_print
  ct_gibbs = 0 # counts when to store state of Markov Chain

  t_elapsed = 0.
  # println("Finished initializing sampler.")

  # println("Running Gibbs sampler.")
  tic();
  for s in 1:n_iter
    gibbs_psi ? update_psi_parameters_partition!(psi_current,PP[perm_current],alpha_current[1]) : nothing ;
    gibbs_T ? update_arrival_times!(T_current,PP[perm_current],alpha_current[1],ia_dist(arrival_params_current)) : nothing ;
    gibbs_sigma ? update_block_order!(perm_current,PP[perm_current],T_current,alpha_current[1]) : nothing ;
    gibbs_ia_params ? update_arrival_params!(arrival_params_current,PP[perm_current],T_current,N,ia_prior_params) : nothing ;
    if arrivals=="crp-coupled"
      alpha_current[1] = arrival_params_current[2]
    else
      gibbs_alpha ? update_ntl_alpha!(alpha_current,PP[perm_current],T_current,ntl_alpha_log_prior,w_alpha) : nothing ;
    end

    if (s > n_burn) && mod(s - n_burn,n_thin)==0
      ct_gibbs += 1 ;
      gibbs_psi ? psi_gibbs[:,ct_gibbs] = psi_current : nothing ;
      gibbs_T ? T_gibbs[:,ct_gibbs] = T_current : nothing ;
      gibbs_alpha ? alpha_gibbs[ct_gibbs] = alpha_current[1] : nothing ;
      gibbs_ia_params ? ia_params_gibbs[:,ct_gibbs] = arrival_params_current : nothing ;
      gibbs_sigma ? perm_gibbs[:,ct_gibbs] = perm_current : nothing ;
      log_joint_gibbs[ct_gibbs] = logp_partition(PP[perm_current],T_current,alpha_current[1],ia_dist(arrival_params_current),true) +
                      ntl_alpha_log_prior(1 - alpha_current[1]) + lpdf_ia_param_prior(arrival_params_current)
    end
    if mod(s,n_print)==0
      t_elapsed += toq();
      println("  >> Finished with ",s," / ",n_iter," samples. Elapsed time is ",t_elapsed," seconds.")
      tic();
    end
  end

  println("Finished running Gibbs sampler.")

  return GibbsSamplerOutput(psi_gibbs,T_gibbs,alpha_gibbs,ia_params_gibbs,perm_gibbs,log_joint_gibbs,t_elapsed)

end


text:
function PositionBuilder(arg0::CorrelatedRandomVectorGenerator, arg1::jdouble, arg2::jdouble, arg3::ObservableSatellite)
    return PositionBuilder((CorrelatedRandomVectorGenerator, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3)
end

function build(obj::PositionBuilder, arg0::Vector{SpacecraftState})
    return jcall(obj, "build", Position, (Vector{SpacecraftState},), arg0)
end


text:
# Implementation of Franz and Scholkopf
# A unifying view of Wiener and Volterra theory and polynomial
# kernel regression

abstract type VolterraKernel end

export AdaptiveKernel,
InhomogeneousPolyKernel,
GaussianKernel,
VovkKernel

immutable AdaptiveKernel <: VolterraKernel
    a::Vector{Float64}
end

(k::AdaptiveKernel)(h) = sum(k.a.^2.*[h^n for n in 0:length(k.a)-1])
(k::AdaptiveKernel)(x1,x2) = k(dot(x1,x2))
weights(k::AdaptiveKernel,n) = k.a[n+1]
weights(k::AdaptiveKernel) = k.a

immutable InhomogeneousPolyKernel <: VolterraKernel
    p
end

(k::InhomogeneousPolyKernel)(h) = (1+h)^k.p
(k::InhomogeneousPolyKernel)(x1,x2) = k(dot(x1,x2))
weights(k::InhomogeneousPolyKernel,n) = binomial(k.p,n)
weights(k::InhomogeneousPolyKernel) = [weights(k,n) for n in 0:p]

immutable GaussianKernel <: VolterraKernel
end

(k::GaussianKernel)(h) = exp(h)
(k::GaussianKernel)(x1,x2) = k(dot(x1,x2))
weights(k::GaussianKernel,n) = 1/factorial(n)

immutable VovkKernel <: VolterraKernel
    α
end

(k::VovkKernel)(h) = (1-h)^(-k.α)
(k::VovkKernel)(x1,x2) = k(dot(x1,x2))
weights(k::VovkKernel,n) = binomial(-k.α,n)*(-1)^n

## Now we need to do kernel regression

"""
Form the kernel matrix out of a data matrix.
"""
kernel_matrix(k::VolterraKernel,X1,X2) = map(k,X1'X2)
kernel_matrix(k::VolterraKernel,X) = kernel_matrix(k,X,X)

type VolterraModel <: DischargeModel
    M::Int
    k::VolterraKernel
    β::Vector{Float64}
    H::Matrix{Float64}
    Q::Vector{Float64}
end

VolterraModel(M,N,k=1) = VolterraModel(M,k,zeros(N),zeros(M,N),zeros(N))

"""
Fitting a Volterra Series Model simply requires inversion of the kernel matrix.
"""
function StatsBase.fit(M::VolterraModel,H::Matrix{Float64},Q::Vector{Float64})
    VolterraModel(M.M,M.k,kernel_matrix(M.k,H)\Q,H,Q)
end

function StatsBase.fit!(M::VolterraModel,H::Matrix{Float64},Q::Vector{Float64})
    M.β = kernel_matrix(M.k,H)\Q
    M.H = H
    M.Q = Q
    M
end

StatsBase.coef(M::VolterraModel) = M.β
StatsBase.nobs(M::VolterraModel) = length(M.Q)
StatsBase.model_response(M::VolterraModel) = M.Q
function StatsBase.predict(M::VolterraModel,H)
    vec(M.β'kernel_matrix(M.k,M.H,H))
end
StatsBase.predict(M::VolterraModel) = predict(M,M.H)
StatsBase.residuals(M::VolterraModel) = M.Q-predict(M)
StatsBase.residuals(M::VolterraModel,H,Q) = Q-predict(M,H)

type RegularizedVolterraModel<: DischargeModel
    M::Int
    k::VolterraKernel
    λ::Float64
    R::AbstractMatrix
    β::Vector{Float64}
    H::Matrix{Float64}
    Q::Vector{Float64}
end

RegularizedVolterraModel(M,N,k=1,λ=0.0) = RegularizedVolterraModel(M,InhomogeneousPolyKernel(k),λ,eye(N-M+1,N-M+1),zeros(N),zeros(M,N),zeros(N))

function StatsBase.fit(M::RegularizedVolterraModel,H::Matrix{Float64},Q::Vector{Float64})
    RegularizedVolterraModel(M.M,
                             M.k,
                             M.λ,
                             M.R,
                             (kernel_matrix(M.k,H)+M.λ*M.R)\Q,
                             H,
                             Q)
end

function StatsBase.fit!(M::RegularizedVolterraModel,H::Matrix{Float64},Q::Vector{Float64})
    M.β = (kernel_matrix(M.k,H)+M.λ*M.R)\Q
    M.H = H
    M.Q = Q
    M
end

StatsBase.coef(M::RegularizedVolterraModel) = M.β
StatsBase.nobs(M::RegularizedVolterraModel) = length(M.Q)
StatsBase.model_response(M::RegularizedVolterraModel) = M.Q
function StatsBase.predict(M::RegularizedVolterraModel,H)
    vec(M.β'kernel_matrix(M.k,M.H,H))
end
StatsBase.predict(M::RegularizedVolterraModel) = predict(M,M.H)
StatsBase.residuals(M::RegularizedVolterraModel) = M.Q-predict(M)
StatsBase.residuals(M::RegularizedVolterraModel,H,Q) = Q-predict(M,H)

"""
To get the Volterra operator, we have to calculate a design matrix.

This is Φn = [ϕn(x_1),ϕn(x_2),...]

where ϕn(x) collects all monomials of order n of the entries of x
"""
function ϕ(x,n)
    out = ones(ntuple(y->length(x),n))
    for I in CartesianRange(size(out))
        for i in 1:length(I)
            out[I]*=x[I[i]]
        end
    end
    vec(out)
end

function designmatrix(H,n)
    a,b = size(H)
    Φ = zeros(a^n,b)
    for i in 1:b
        Φ[:,i] = ϕ(H[:,i],n)
    end
    Φ
end

"""
Derive the Volterra operator of order n from a model
"""
operator(M::VolterraModel,n) = weights(M.k,n).*designmatrix(M.H,n)*M.β
operator(M::RegularizedVolterraModel,n) = weights(M.k,n).*designmatrix(M.H,n)*M.β


text:
######################################################################
# bearing.jl
#
# bearing only sensing modality
######################################################################

struct BearingOnly <: Sensor
    noise_sigma::Float64		# noise std deviation (degrees)

    BearingOnly(ns::Real) = new(ns)
    BearingOnly() = new(10.0)
end

function observe(theta::LocTuple, s::BearingOnly, p::Pose)
    truth = true_bearing(p, theta)
    noise = s.noise_sigma * randn()
    return mod(truth + noise, 360.0)
end

function O(bo::BearingOnly, theta::TargetTuple, p::Pose, o)

    # Calculate true bearing, and find distance to bin edges
    ang_deg = true_bearing(p, theta)
    o_diff = fit_180(o - ang_deg)

    # now look at probability
    d = Normal(0, bo.noise_sigma)
    p = pdf(d, o_diff)
    return p
end


text:
using StatsBase

"""
    interaural_coherence(x::SampleBuf, lags::Unitful.Time)

Compute the interaural coherence of a two channel sound signals.

Interaural coherence (IAC) is commonly defined as the
peak of the cross-correlation coefficient of the signals at the two ears [1, 2].
It is commonly computed over a restricted range of `lags` of the cross-correlation function.

Inputs
------
* `x` data in the form of SampledSignals.SampleBuf. Must be two channels of audio.
* `lags` time range of lags to be used for finding maximum in cross correlation function.
  If lags=0, then the entire function will be used, effecively same as lags=Inf.

References
----------

1. Chait, M., Poeppel, D., de Cheveigne, A., and Simon, J.Z. (2005). Human auditory cortical processing of changes in interaural correlation. J Neurosci 25, 8518-8527.
2. Aaronson, N.L., and Hartmann, W.M. (2010). Interaural coherence for noise bands: waveforms and envelopes. J Acoust Soc Am 127, 1367-1372.

Example
-------
```julia
correlation = 0.6
source = CorrelatedNoiseSource(Float64, 48000, 2, 0.1, correlation)
a = read(source, 3u"s")
@test interaural_coherence(a.data) ≈ correlation atol = 0.025
```

"""
function interaural_coherence(x::T; lags::AbstractQuantity=0u"s")::Float64 where {T<:SampledSignals.SampleBuf}

    lags_seconds = lags |> u"s" |> ustrip
    lags_samples = Int(lags_seconds * x.samplerate)
    interaural_coherence(x.data; lags=lags_samples)
end

function interaural_coherence(x::T; lags::AbstractQuantity=0u"s")::Float64 where {T<: DummySampleSink}

    lags_seconds = lags |> u"s" |> ustrip
    lags_samples = Int(lags_seconds * x.samplerate)
    interaural_coherence(x.buf; lags=lags_samples)
end

function interaural_coherence(x::Array{T, 2}; lags::Int=0)::Float64 where {T<:Number}

    if lags == 0
        lags = size(x, 1) - 1
    end

    cc = crosscor(x[:, 1], x[:, 2], -lags:lags)
    iac = maximum(cc)
end


text:
code_0 = 0
code_1 = 1
l_a = [code_0, code_1]
code_a = "a"
code_b = "b"
l_b = [code_a, code_b]
if abspath(PROGRAM_FILE) == @__FILE__
    for i in l_a
        println(i)
    end
    for j in l_b
        println(j)
    end
    if "a" ∈ ["a", "b"]
        println("OK")
    end
end


text:
export cycle, cycleir, @cycle, @cycleir

# zt - Fixme this is type unstable
"`cycle(f, args...)` yields `xs_` such `f(xs_...) == f(args...)``"
cycle(θ, f, args...) =
  invertinvoke(f, Base.typesof(args...), f(args...), θ)

cycleir(f, args...) =
  invertir(typeof(f), Base.typesof(args...))

cycle(f, args...) =
  cycle(defθ(), f, args...)

function cyclem(ex)
  if IRTools.isexpr(ex, :call)
    f, args = ex.args[1], ex.args[2:end]
  elseif IRTools.isexpr(ex, :do)
    f, args = ex.args[1].args[1], vcat(ex.args[2], ex.args[1].args[2:end])
  else
    error("@code_ir f(args...)")
  end
  esc(:(cycle(ParametricInversion.defθ(), $f, $args...)))
end

function cycleirm(ex)
  if IRTools.isexpr(ex, :call)
    f, args = ex.args[1], ex.args[2:end]
  elseif IRTools.isexpr(ex, :do)
    f, args = ex.args[1].args[1], vcat(ex.args[2], ex.args[1].args[2:end])
  else
    error("@code_ir f(args...)")
  end
  esc(:(cycleir($f, $args...)))
end

macro cycle(ex)
  cyclem(ex)
end

macro cycleir(ex)
  cycleirm(ex)
end


text:
using Test
using MPI

MPI.Init()

info = MPI.Info()
@test typeof(info) == MPI.Info
@test info.val != MPI.Info
MPI.Info_set(info, "foo", "fast")
MPI.Info_set(info, "bar", "evenfaster")
@test MPI.Info_get(info, "foo" ) == "fast"
MPI.Info_delete(info, "bar")
@test MPI.Info_get(info, "bar" ) == ""
finalize(info)
@test info.val == MPI.MPI_INFO_NULL

MPI.Finalize()


text:
function twodOseen(x,eConn,innerNodes,outerNodes,advectionVelocity,
                   μ::Float64=1.0, ω::Float64=1.0)
#
#  Solves the twod-dimensional steady-state Oseen equations
#       (v⋅∇)u + (u⋅∇)v = -∇p + ∇⋅(μ(∇u+∇u')) + f,
#                     0 = ∇⋅u
#  Ω is the domain between two concentric circles.  v is a prescribed
#  advection velocity and f is used to generate a closed-form solution.
#
#  We prescribe zero Dirichlet conditions on the inner boundary,
#  and a uniform rotation of the outer boundary.
#
#  The parameter ϵ is used to implement a penalty method (computed using μ)
#
#  As a unit test, choose a value of f such that the solution u when
#  v = (-ωy,ωx) becomes u = (r-1)*(-ωy,ωx)/r and p = x^2+y^2.
#
#  v = (u,v) as the advection velocity, then u = (u,v) as the solution
#---------------------------------------------------------------------------78--

#  include("twodQuadratureRule.jl")
#  include("twodShape.jl")
#  include("twodBilinear.jl")
#  include("twodLinForm.jl")

  ϵ = 1e-8/μ   # penalty parameter
  rule  = 7    # points in quadrature formula

  function f(x::Array{Float64,2},μ,ω)
    d = size(x,1)
    ansX = zeros(Float64,d,1)
    ansY = zeros(Float64,d,1)
    for i=1:d
      rad = sqrt(x[i,1]^2+x[i,2]^2)
      ansX[i] = (2*x[i,1]*rad^5-3μ*ω*x[i,2]^3+ω^2*x[i,1]*(rad^2-rad^5) -
                  3*μ*ω*x[i,1]^2*x[i,2] + 4*μ*ω*x[i,2]*rad^2)/rad^5
      ansY[i] = (2*x[i,2]*rad^5+3μ*ω*x[i,1]^3+ω^2*x[i,2]*(rad^2-rad^5) +
                  3*μ*ω*x[i,1]*x[i,2]^2 - 4*μ*ω*x[i,1]*rad^2)/rad^5
    end
    return ansX,ansY
  end

  function uEx(x::Array{Float64,1},μ=1.0,ω=-1.0)
    rad = sqrt(x[1]^2+x[2]^2)
    factor = (rad-1.0)/rad

    ansX = -x[2]*ω*factor
    ansY =  x[1]*ω*factor
    ansP = rad^2

    return ansX, ansY, ansP
  end

  #  Get problem dimensions
  #-----------------------------------------------------------------------------
  nNodes = size(x,1)
  nElements = size(eConn,1)

#  innerNodes = sort!(innerNodes)  # assume they come in sorted
#  outerNodes = sort!(outerNodes)
  nInnerNodes = length(innerNodes)
  nOuterNodes = length(outerNodes)
#  nDirichlet = nInnerNodes+nOuterNodes   # use for Dirichlet BCs
  nDirichlet = nOuterNodes

  #  Set the index into equation numbers
  #-----------------------------------------------------------------------------
  ideU = zeros(Int64,nNodes,2)

  global nUnk = 0
  for i=1:nNodes
    global nUnk = nUnk+1
    ideU[i,1] = nUnk

    global nUnk = nUnk+1
    ideU[i,2] = nUnk
  end

  ideP = zeros(Int64,nNodes)
  for n_el=1:nElements
    vertex = eConn[n_el,1:3]
    if ( ideP[vertex[1]]==0 )
      global nUnk = nUnk + 1
      ideP[vertex[1]] = nUnk
    end

    if ( ideP[vertex[2]]==0 )
      global nUnk = nUnk + 1
      ideP[vertex[2]] = nUnk
    end

    if ( ideP[vertex[3]]==0 )
      global nUnk = nUnk + 1
      ideP[vertex[3]] = nUnk
    end
  end
  nP = nUnk - 2*nNodes

  #  Integrate system matrices, element-by-element (each independently)
  #-------------------------------------------------------------------------78--
  r,s,w = twodQuadratureRule(rule)
  μ_g   = μ*ones(rule)
  ϵ_g   = ϵ*ones(rule)
  one   = ones(rule)

  nElDOF   = size(eConn,2)
  nElDOFsq = nElDOF*nElDOF
  nEntries = nElements*(2*nElDOF+3)*(2*nElDOF+3)

  II  = Array{Int64,1}(undef, nEntries)
  JJ  = Array{Int64,1}(undef, nEntries)
  AA  = Array{Float64,1}(undef, nEntries)
  b   = zeros(Float64,nUnk,1)

  for k=1:nElements
    nLocal = eConn[k,:][:]
    xLocal = x[nLocal,:]

    xg,wg,ϕ,ϕ_x,ϕ_y = twodShape( xLocal       , r, s, w )
    xg,wg,ψ,ψ_x,ψ_y = twodShape( xLocal[1:3,:], r, s, w )

    # evaluate forcing function at quadrature points
    fx_g, fy_g = f(xg,μ,ω)

    u_g   = ϕ*advectionVelocity[nLocal,1]
    u_xg  = ϕ_x*advectionVelocity[nLocal,1]
    u_yg  = ϕ_y*advectionVelocity[nLocal,1]
    v_g   = ϕ*advectionVelocity[nLocal,2]
    v_xg  = ϕ_x*advectionVelocity[nLocal,2]
    v_yg  = ϕ_y*advectionVelocity[nLocal,2]

    A11Loc = -twodBilinear( 2*μ_g, ϕ_x, ϕ_x, wg) -
              twodBilinear(   μ_g, ϕ_y, ϕ_y, wg) -
              twodBilinear(   u_g, ϕ_x, ϕ  , wg) -
              twodBilinear(   v_g, ϕ_y, ϕ  , wg) -
              twodBilinear(  u_xg, ϕ  , ϕ  , wg)

    A12Loc = -twodBilinear(   μ_g, ϕ_x, ϕ_y, wg) -
              twodBilinear(  u_yg, ϕ  , ϕ  , wg)

    A21Loc = -twodBilinear(   μ_g, ϕ_y, ϕ_x, wg) -
              twodBilinear(  v_xg, ϕ  , ϕ  , wg)

    A22Loc = -twodBilinear(   μ_g, ϕ_x, ϕ_x, wg) -
              twodBilinear( 2*μ_g, ϕ_y, ϕ_y, wg) -
              twodBilinear(   u_g, ϕ_x, ϕ  , wg) -
              twodBilinear(   v_g, ϕ_y, ϕ  , wg) -
              twodBilinear(  v_yg, ϕ  , ϕ  , wg)

    B1Loc  =  twodBilinear(   one, ϕ_x, ψ  , wg)
    B2Loc  =  twodBilinear(   one, ϕ_y, ψ  , wg)

    MLoc   = -twodBilinear(   ϵ_g, ψ  , ψ  , wg)

    F1Loc  = twodLinForm( fx_g, ϕ, wg)
    F2Loc  = twodLinForm( fy_g, ϕ, wg)

    index = (k-1)*(2*nElDOF+3)^2  # compute base index (since k could be in any order)
    lDOF = ideU[nLocal,1][:]

    # u-momentum equations
    for nt = 1:nElDOF
      nTest = ideU[nLocal[nt],1]
      for nu = 1:nElDOF
        nUnkU = ideU[nLocal[nu],1]
        nUnkV = ideU[nLocal[nu],2]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkU
        AA[index] = A11Loc[nt,nu]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkV
        AA[index] = A12Loc[nt,nu]
      end
      for np = 1:3  # special to this linear element
        nUnkP = ideP[nLocal[np]]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkP
        AA[index] = B1Loc[np,nt]
      end

      b[nTest] = b[nTest] + F1Loc[nt]
    end

    # v-momentum equations
    for nt = 1:nElDOF
      nTest = ideU[nLocal[nt],2]
      for nu = 1:nElDOF
        nUnkU = ideU[nLocal[nu],1]
        nUnkV = ideU[nLocal[nu],2]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkU
        AA[index] = A21Loc[nt,nu]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkV
        AA[index] = A22Loc[nt,nu]
      end
      for np = 1:3  # special to this linear element
        nUnkP = ideP[nLocal[np]]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkP
        AA[index] = B2Loc[np,nt]
      end

      b[nTest] = b[nTest] + F2Loc[nt]
    end

    # continuity equations
    for nt = 1:3
      nTest = ideP[nLocal[nt]]
      for nu = 1:nElDOF
        nUnkU = ideU[nLocal[nu],1]
        nUnkV = ideU[nLocal[nu],2]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkU
        AA[index] = B1Loc[nt,nu]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkV
        AA[index] = B2Loc[nt,nu]
      end
      for np = 1:3  # special to this linear element
        nUnkP = ideP[nLocal[np]]

        index = index + 1
        II[index] = nTest
        JJ[index] = nUnkP
        AA[index] = MLoc[nt,np]
      end

    end
  end

  A = sparse(II,JJ,AA)

  interiorNodes = convert(Array{Int64,1},1:nNodes)
  innerN = convert(Array{Int64,1},innerNodes)
  outerN = convert(Array{Int64,1},outerNodes)
  setdiff!(interiorNodes,innerN)
  setdiff!(interiorNodes,outerN)

  iB = vcat(innerN,outerN)
  nDirichlet = length(iB)
# knownIndexU = Array{Int64,1}(undef,nDirichlet)
# knownIndexV = Array{Int64,1}(undef,nDirichlet)
  knownIndexU = zeros(Int64,nDirichlet,1)
  knownIndexV = zeros(Int64,nDirichlet,1)

  dirichletU = zeros(Float64,nDirichlet,1)
  dirichletV = zeros(Float64,nDirichlet,1)
  for i=1:length(innerN)
    knownIndexU[i] = ideU[innerN[i],1]
    knownIndexV[i] = ideU[innerN[i],2]
    uE,vE,pE = uEx(x[innerN[i],:],μ,ω)
    dirichletU[i] = uE
    dirichletV[i] = vE
  end

  ni = length(innerN)
  for i=1:length(outerN)
    idx = ni+i
    knownIndexU[idx] = ideU[outerN[i],1]
    knownIndexV[idx] = ideU[outerN[i],2]
    uE,vE,pE = uEx(x[outerN[i],:],μ,ω)
    dirichletU[idx] = uE
    dirichletV[idx] = vE
  end

  nUnknowns = 2*length(interiorNodes) + nP
  unknownIndex = Array{UInt64,1}(undef,nUnknowns)

  global nUnknownCounter = 0
  for i ∈ interiorNodes
    global nUnknownCounter = nUnknownCounter+1
    unknownIndex[nUnknownCounter] = ideU[i,1]

    global nUnknownCounter = nUnknownCounter+1
    unknownIndex[nUnknownCounter] = ideU[i,2]
  end

  getP = zeros(Int64,3*nNodes,1)
  for nNode=1:nNodes
    if ideP[nNode]>0
      global nUnknownCounter = nUnknownCounter + 1
      unknownIndex[nUnknownCounter] = ideP[nNode]
      getP[ideP[nNode]] = nUnknownCounter
    end
  end

  rhs = -b[unknownIndex]-A[unknownIndex,knownIndexU[:]]*dirichletU-A[unknownIndex,knownIndexV[:]]*dirichletV

  uvp = A[unknownIndex,unknownIndex]\rhs

  velocity = zeros(Float64,nNodes,2)

  global nUnknownCounter = 0
  for i ∈ interiorNodes
    global nUnknownCounter = nUnknownCounter+1
    velocity[i,1] = uvp[nUnknownCounter]

    global nUnknownCounter = nUnknownCounter + 1
    velocity[i,2] = uvp[nUnknownCounter]
  end

  for i=1:nDirichlet
    velocity[iB[i],1] = dirichletU[i]
    velocity[iB[i],2] = dirichletV[i]
  end

  pressure = zeros(Float64,nNodes,1)

  vv = eConn[:,1:3]
  vvv = vv[:]
  vertices = unique(sort(vvv))
  for i ∈ vertices
    pressure[i] = uvp[getP[ideP[i]]]
  end
  pressure = TriMesh_PromoteL2Q(pressure,eConn)

#  return x,eConn,velocity,pressure

  #  Compute the exact solution for the unit test
  velocityExact = zeros(Float64,nNodes,2)
  pressureExact = zeros(Float64,nNodes,1)
  for n=1:nNodes
    uE,vE,pE = uEx(x[n,:],μ,ω)
    velocityExact[n,1] = uE
    velocityExact[n,2] = vE
    pressureExact[n] = pE
  end

  # Compute the norm of the approximation error for the unit test
  M = twodMassMatrix(x,eConn)
#  vDiff = velocity-velocityExact
  C     = pressureExact[1]-pressure[1]
#  pDiff = pressure-pressureExact.+C
  pressureExact = pressureExact.-C
  return velocity,pressure,velocityExact,pressureExact
#  return sqrt(dot(vDiff[:,1],M,vDiff[:,1])) +
#         sqrt(dot(vDiff[:,2],M,vDiff[:,2])) +
#         sqrt(dot(pDiff,M,pDiff))
end


text:
using Documenter

makedocs()

deploydocs(
   #options
   deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-cinder"),
   repo = "github.com/pranavtbhat/JuliaCookbook.git",
   julia = "0.5"
)

max_stars_repo_name: pranavtbhat/JuliaCookbook

max_issues_repo_name: pranavtbhat/JuliaCookbook

max_forks_repo_name: pranavtbhat/JuliaCookbook


text:
#!/usr/bin/env julia
using Test
using Random
using Parameters
using ResumableFunctions

no = nothing

want(x) = begin println("# ",x);
                include("./"*x*".jl") end

want("it")
it=It()

Random.seed!(it.seed)

want("lib")
want("col")
want("tests")

go() = want("nb")
# say(Some(w=2))

max_stars_repo_name: timm/jut

max_issues_repo_name: timm/jut

max_forks_repo_name: timm/jut


text:
using ServerMetrics
using HTTP
using Dates

Base.@kwdef struct UptimeMetrics <: AbstractMetricCollection
    server_uptime_seconds = Gauge()
    server_heartbeats_total = Counter()
end

metrics = UptimeMetrics()
ServerMetrics.publish_metrics_from(metrics; overwrite=true)
startup_timestamp = Dates.datetime2unix(Dates.now())
ServerMetrics.@spawn_sticky_periodic_task "UptimeTracker" Dates.Second(1) begin
    time_delta = Dates.datetime2unix(Dates.now()) - startup_timestamp
    set!(metrics.server_uptime_seconds, time_delta)
    inc!(metrics.server_heartbeats_total)
end

HTTP.serve() do http_request
    ServerMetrics.handle_metrics()
end


text:
"""
    sampling_codebooks(X, k, m)

Build `m` codebooks for input matrix `X` by sampling `k` prototype vectors.
"""
function sampling_codebooks(X::AbstractMatrix{T}, k::Int, m::Int) where {T}
    nrows, ncols = size(X)
    cbooks = Vector{Matrix{T}}(undef, m)
    @inbounds for i in 1:m
        rr = rowrange(nrows, m, i)
        vectors = Matrix{T}(X[rr, sample(1:ncols, k, replace=false)])
        cbooks[i] = vectors
    end
    return cbooks
end

"""
    pq_codebooks(X, k, m [;distance=DEFAULT_DISTANCE, maxiter=DEFAULT_PQ_MAXITER])

Build `m` codebooks for input matrix `X` using `k` sub-space centers obtained
using k-means clustering, the distance `distance` and `maxiter` iterations.

# References:
  * [Jègou et al. 2011](https://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf)
"""
function pq_codebooks(X::AbstractMatrix{T}, k::Int, m::Int;
                      distance::Distances.PreMetric=DEFAULT_DISTANCE,
                      maxiter::Int=DEFAULT_PQ_MAXITER) where {T}
    nrows = size(X, 1)
    cbooks = Vector{Matrix{T}}(undef, m)
    @inbounds for i in 1:m
        rr = rowrange(nrows, m, i)
        model = kmeans(X[rr, :], k,
                       maxiter=maxiter,
                       distance=distance,
                       init=:kmpp,
                       display=:none)
        cbooks[i] = model.centers
    end
    return cbooks
end

"""
    opq_codebooks(X, k, m [;distance=DEFAULT_DISTANCE, maxiter=DEFAULT_PQ_MAXITER])

Build `m` codebooks for input matrix `X` using `k` sub-space centers obtained
using 'cartesian' k-means clustering, the distance `distance` and `maxiter` iterations.

# References:
  * [Ge et al. 2014](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/opq_tr.pdf)
"""
function opq_codebooks(X::AbstractMatrix{T}, k::Int, m::Int;
                       distance::Distances.PreMetric=DEFAULT_DISTANCE,
                       maxiter::Int=DEFAULT_OPQ_MAXITER) where {T}
    # Initialize R, rotated data
    nrows, ncols = size(X)
    X̂ = zeros(T, size(X))
    R = diagm(0 => ones(T, nrows))
    Xr = R * X

    # Initialize codebooks, codes
    cweights   = zeros(Int, k)
    costs      = zeros(ncols)
    counts     = zeros(Int, k)
    unused     = Int[]
    to_update_z = zeros(Bool, k)
    codes = fill(zeros(Int, ncols), m)
    cbooks = sampling_codebooks(X, k, m)
    @inbounds for i in 1:m
        rr = rowrange(nrows, m, i)
        dists = Distances.pairwise(distance, cbooks[i], Xr[rr, :], dims=2)
        Clustering.update_assignments!(dists, true, codes[i], costs, counts, to_update_z, unused)
        X̂[rr,:] .= cbooks[i][:, codes[i]]
    end

    # Run optimization
    to_update = fill(ones(Bool, k), m)
    for it = 1:maxiter
        # Update R using orthogonal Procustes closed form solution
        # and update rotated data matrix X̂
        U, _, V = svd(X * X̂', full=false)
        R = V * U'
        Xr = R * X
        @inbounds for i in 1:m
            rr = rowrange(nrows, m, i)
            # Update subspace cluster centers
            Clustering.update_centers!(Xr[rr, :], nothing, codes[i], to_update[i], cbooks[i], cweights)
            # Update subspace data assignments
            dists = Distances.pairwise(distance, cbooks[i], Xr[rr, :], dims=2)
            Clustering.update_assignments!(dists, false, codes[i], costs, counts, to_update[i], unused)
            X̂[rr, :] .= cbooks[i][:, codes[i]]
        end
        @debug "OPQ: Iteration $it, error = $(sum((R'*X̂ .- X).^2)./ncols)"
    end
    return cbooks, R
end

"""
    rvq_codebooks(X, k, m [;distance=DEFAULT_DISTANCE, maxiter=DEFAULT_RVQ_MAXITER])

Build `m` codebooks for input matrix `X` using `k` layers of residual vector cluster centers,
obtained using k-means clustering, the distance `distance` and `maxiter` iterations.

# References:
  * [Chen et al. 2010](https://www.mdpi.com/1424-8220/10/12/11259/htm)
"""
function rvq_codebooks(X::AbstractMatrix{T}, k::Int, m::Int;
                       distance::Distances.PreMetric=DEFAULT_DISTANCE,
                       maxiter::Int=DEFAULT_RVQ_MAXITER) where {T}
    nrows = size(X, 1)
    cbooks = Vector{Matrix{T}}(undef, m)
    Xr = copy(X)
    @inbounds for i in 1:m
        model = kmeans(Xr, k,
                       maxiter=maxiter,
                       distance=distance,
                       init=:kmpp,
                       display=:none)
        cbooks[i] = model.centers
        codes = model.assignments
        Xr .-= cbooks[i][:, codes]
    end
    return cbooks
end


text:
using DashBootstrapComponents, DashHtmlComponents

standalone_radio_check = html_div([
    dbc_formgroup(
        [
            dbc_checkbox(id="standalone-checkbox", className="form-check-input"),
            dbc_label(
                "This is a checkbox",
                html_for="standalone-checkbox",
                className="form-check-label",
            ),
        ],
        check=true,
    ),
    dbc_formgroup(
        [
            dbc_radiobutton(id="standalone-radio", className="form-check-input"),
            dbc_label(
                "This is a radio button",
                html_for="standalone-radio",
                className="form-check-label",
            ),
        ],
        check=true,
    ),
    html_br(),
    html_p(id="standalone-radio-check-output"),
]);

callback!(
    app,
    Output("standalone-radio-check-output", "children"),
    Input("standalone-checkbox", "checked"),
    Input("standalone-radio", "checked"),
) do checkbox_checked, radio_checked
    if checkbox_checked == 1 && radio_checked == 1
        return "Both checked."
    elseif checkbox_checked == 1 || radio_checked == 1
        return "One checked."
    else
        return "None checked."
    end
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediastore_data
using AWS.Compat
using AWS.UUIDs

"""
    delete_object(path)
    delete_object(path, params::Dict{String,<:Any})

Deletes an object at the specified path.

# Arguments
- `path`: The path (including the file name) where the object is stored in the container.
  Format: <folder name>/<folder name>/<file name>

"""
delete_object(Path; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("DELETE", "/$(Path)"; aws_config=aws_config)
delete_object(Path, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("DELETE", "/$(Path)", params; aws_config=aws_config)

"""
    describe_object(path)
    describe_object(path, params::Dict{String,<:Any})

Gets the headers for an object at the specified path.

# Arguments
- `path`: The path (including the file name) where the object is stored in the container.
  Format: <folder name>/<folder name>/<file name>

"""
describe_object(Path; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("HEAD", "/$(Path)"; aws_config=aws_config)
describe_object(Path, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("HEAD", "/$(Path)", params; aws_config=aws_config)

"""
    get_object(path)
    get_object(path, params::Dict{String,<:Any})

Downloads the object at the specified path. If the object’s upload availability is set to
streaming, AWS Elemental MediaStore downloads the object even if it’s still uploading the
object.

# Arguments
- `path`: The path (including the file name) where the object is stored in the container.
  Format: <folder name>/<folder name>/<file name> For example, to upload
  the file mlaw.avi to the folder path premiumcanada in the container movies, enter the path
  premium/canada/mlaw.avi. Do not include the container name in this path. If the path
  includes any folders that don't exist yet, the service creates them. For example, suppose
  you have an existing premium/usa subfolder. If you specify premium/canada, the service
  creates a canada subfolder in the premium folder. You then have two subfolders, usa and
  canada, in the premium folder.  There is no correlation between the path to the source and
  the path (folders) in the container in AWS Elemental MediaStore. For more information about
  folders and how they exist in a container, see the AWS Elemental MediaStore User Guide. The
  file name is the name that is assigned to the file that you upload. The file can have the
  same name inside and outside of AWS Elemental MediaStore, or it can have the same name. The
  file name can include or omit an extension.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Range"`: The range bytes of an object to retrieve. For more information about the Range
  header, see http://www.w3.org/Protocols/rfc2616/rfc2616-sec14.html#sec14.35. AWS Elemental
  MediaStore ignores this header for partially uploaded objects that have streaming upload
  availability.
"""
get_object(Path; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("GET", "/$(Path)"; aws_config=aws_config)
get_object(Path, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("GET", "/$(Path)", params; aws_config=aws_config)

"""
    list_items()
    list_items(params::Dict{String,<:Any})

Provides a list of metadata entries about folders and objects in the specified folder.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return per API request. For example, you
  submit a ListItems request with MaxResults set at 500. Although 2,000 items match your
  request, the service returns no more than the first 500 items. (The service also returns a
  NextToken value that you can use to fetch the next batch of results.) The service might
  return fewer results than the MaxResults value. If MaxResults is not included in the
  request, the service defaults to pagination with a maximum of 1,000 results per page.
- `"NextToken"`: The token that identifies which batch of results that you want to see. For
  example, you submit a ListItems request with MaxResults set at 500. The service returns the
  first batch of results (up to 500) and a NextToken value. To see the next batch of results,
  you can submit the ListItems request a second time and specify the NextToken value. Tokens
  expire after 15 minutes.
- `"Path"`: The path in the container from which to retrieve items. Format: <folder
  name>/<folder name>/<file name>
"""
list_items(; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("GET", "/"; aws_config=aws_config)
list_items(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("GET", "/", params; aws_config=aws_config)

"""
    put_object(body, path)
    put_object(body, path, params::Dict{String,<:Any})

Uploads an object to the specified path. Object sizes are limited to 25 MB for standard
upload availability and 10 MB for streaming upload availability.

# Arguments
- `body`: The bytes to be stored.
- `path`: The path (including the file name) where the object is stored in the container.
  Format: <folder name>/<folder name>/<file name> For example, to upload
  the file mlaw.avi to the folder path premiumcanada in the container movies, enter the path
  premium/canada/mlaw.avi. Do not include the container name in this path. If the path
  includes any folders that don't exist yet, the service creates them. For example, suppose
  you have an existing premium/usa subfolder. If you specify premium/canada, the service
  creates a canada subfolder in the premium folder. You then have two subfolders, usa and
  canada, in the premium folder.  There is no correlation between the path to the source and
  the path (folders) in the container in AWS Elemental MediaStore. For more information about
  folders and how they exist in a container, see the AWS Elemental MediaStore User Guide. The
  file name is the name that is assigned to the file that you upload. The file can have the
  same name inside and outside of AWS Elemental MediaStore, or it can have the same name. The
  file name can include or omit an extension.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Cache-Control"`: An optional CacheControl header that allows the caller to control the
  object's cache behavior. Headers can be passed in as specified in the HTTP at
  https://www.w3.org/Protocols/rfc2616/rfc2616-sec14.html#sec14.9. Headers with a custom
  user-defined value are also accepted.
- `"Content-Type"`: The content type of the object.
- `"x-amz-storage-class"`: Indicates the storage class of a Put request. Defaults to
  high-performance temporal storage class, and objects are persisted into durable storage
  shortly after being received.
- `"x-amz-upload-availability"`: Indicates the availability of an object while it is still
  uploading. If the value is set to streaming, the object is available for downloading after
  some initial buffering but before the object is uploaded completely. If the value is set to
  standard, the object is available for downloading only when it is uploaded completely. The
  default value for this header is standard. To use this header, you must also set the HTTP
  Transfer-Encoding header to chunked.
"""
put_object(Body, Path; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("PUT", "/$(Path)", Dict{String, Any}("Body"=>Body); aws_config=aws_config)
put_object(Body, Path, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = mediastore_data("PUT", "/$(Path)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Body"=>Body), params)); aws_config=aws_config)


text:
#!/usr/bin/julia
#
# Implementation of simulation variations common across multiple simulation methods.
#
export AdaptiveΔt;

"Adaptive timestep parameters."
struct AdaptiveΔt <: SimVariation
  "Maximum error before reducing Δt."
  ϵmax::Float64;
  "Minimum error before increasing Δt."
  ϵmin::Float64;

  "Factor to reduce Δt by."
  Δtdown::Float64;
  "Factor to increase Δt by."
  Δtup::Float64;
  "Maximum allowable Δt."
  Δtmax::Float64;
  "Minimum allowable Δt."
  Δtmin::Float64;

  "Maximum number of repeated evaluations per iteration."
  maxeval::UInt;

  """
  Maximum number successive iterations where the evaluation limit was reached
  before a different parameter should be altered.

  This might be undefined. As an example, in the Lanczos method, when this
  limit is reached the basis is regenerated.
  """
  maxevallimit::Union{UInt, Nothing};

  # TODO Implement inner constructor to check validity of passed arguments.
end
AdaptiveΔt(;ϵmax::Real=1e-8, ϵmin::Real=1e-12,
            Δtdown::Real=1.2, Δtup::Real=1.5, Δtmax::Real=0.5, Δtmin::Real=1e-5,
            maxeval::Real=5, maxevallimit::Real=3) =
  AdaptiveΔt(ϵmax, ϵmin, Δtdown, Δtup, Δtmax, Δtmin, maxeval, maxevallimit);

"""
    adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number; debug::Bool=true)

Adapt the time step `Δt` using the error in the last propagation step, `ϵ`
according to the settings passed in `par`.
"""
function adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number; debug::Bool=true)
  # Check passed arguments.
  if (ϵ <= 0)
    throw(DomainError(ϵ, "error must be positive."));
  elseif (imag(Δt) == 0 && real(Δt) <= 0)
    throw(DomainError(Δt, "timestep must be positive."));
  elseif (real(Δt) == 0 && imag(Δt) <= 0)
    throw(DomainError(Δt, "timestep must be positive."));
  elseif (real(Δt) > 0 && imag(Δt) > 0)
    throw(DomainError(Δt, "timestep must be either real or imaginary."));
  end

  # Increase timestep if less than the minimum error.
  if (ϵ < par.ϵmin)
    # Check that increasing the timestep is allowed.
    if (abs(Δt) > par.Δtmax)
      println("Reached maximum allowable timestep");
    else
      Δt *= par.Δtup;
      if (debug) println("Increasing timestep to ", Δt); end
    end
  elseif (ϵ > par.ϵmax)
    # Check that decreasing timestep is allowed.
    if (abs(Δt) < par.Δtmin)
      println("Reached minimum allowable timestep");
    # Decrease timestep.
    else
      Δt /= par.Δtdown;
      if (debug) println("Decreasing timestep to ", Δt); end
    end
  end

  return Δt;
end

"""
    adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, evals::Real; debug::Bool=true)

Return adapted time step and decision on if previous evaluation should be
repeated using the number of evaluations already completed.

This returns a `Tuple` with the first value being the new time step to use and
the second a `Boolean` indicating if re-evaluation is required.
"""
function adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, evals::Real;
                  debug::Bool=true)
  newΔt = adapt_Δt(par, ϵ, Δt, debug=debug);

  # Repeat evaluation if time step was reduced and maximum evaluations not
  # reached.
  if (newΔt < Δt && evals < par.maxeval) return (newΔt, true);
  else return (newΔt, false); end
end
"""
    adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, evals::Vector{T};
             debug::Bool=true) where T<:Union{Int,UInt}

Return adapted time step, re-evaluation decision and if the `maxevallimit` was
reached.

The `evals` vector contains the number of evaluations for all completed
propagation steps (including the current one) with zero values for steps not yet
completed.
"""
function adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, evals::Vector{T};
                  debug::Bool=true) where T<:Union{Int,UInt}
  # Make copy of `evals` to avoid changing its values.
  evalscopy = copy(evals);
  # Remove evaluations which haven't occurred yet (evaluation counter is zero).
  filter!(v->(v!=0), evalscopy);

  (Δt, reeval) = adapt_Δt(par, ϵ, Δt, evalscopy[end], debug=debug);

  # Forward result if re-evaluation was not required.
  if !reeval return (Δt, false, false); end

  maxevallimitreached = true;
  # Convert to Int, otherwise result will be cast to UInt and will overflow if result is negative.
  min = max(1, length(evalscopy)-convert(Int,par.maxevallimit));
  for i ∈ evalscopy[min]:evalscopy[end]
    if i < par.maxeval
      maxevallimitreached = false;
      break;
    end
  end

  return (Δt, reeval, maxevallimitreached);
end
"""
    adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, curreval::Real,
             evals::Vector{T}; debug::Bool=true) where T<:Union{Int,UInt}

Here `curreval` contains the index storing the number of evaluations for the
current propagation step within `evals`.
"""
function adapt_Δt(par::AdaptiveΔt, ϵ::Real, Δt::Number, curreval::Real,
                  evals::Vector{T}; debug::Bool=true) where T<:Union{Int,UInt}
  if par.maxevallimit === nothing
    error("maxium evalution limit was not initalised.");
  end

  (Δt, reeval) = adapt_Δt(par, ϵ, Δt, evals[curreval], debug=debug);

  # Forward result if re-evaluation was not required.
  if !reeval return (Δt, false, false); end

  # TODO: Check UInt overflow characteristics here.
  maxevallimitreached = true;
  for i ∈ max(curreval-par.maxevallimit, 1):curreval
    if evals[i] < par.maxeval
      maxevallimitreached = false;
      break;
    end
  end

  return (Δt, reeval, maxevallimitreached);
end