text:
#using Distributions
#include("constants.jl")
#include("orbit.jl")
#include("planet.jl")

if !isdefined(:PlanetarySystemAbstract)
  @compat abstract type PlanetarySystemAbstract end

  immutable PlanetarySystem{StarT<:StarAbstract} <: PlanetarySystemAbstract
    star::StarT
    planet::Vector{Planet}         # TODO OPT: If want to optimize, consider using  something like ImmmutableArrays?
    orbit::Vector{Orbit}           # TODO OPT: If want to optimize, consider using  something like ImmmutableArrays?
    # TODO DETAIL: Setup inner constructor to enforce equal number of planets & orbits
    #PlanetarySystem{StarT}(s::StarAbstract, p::Vector{Planet}, o::Vector{Orbit})
    #  @assert(length(p)==length(o)) # else error(string("Number of planets must match number of orbits: Np= ",length(p)," No= ",length(o)))
    #  new(s,p,o)
    #end
  end
  PlanetarySystemSingleStar = PlanetarySystem{SingleStar}

end
PlanetarySystem{StarT<:StarAbstract}(s::StarT) = PlanetarySystem(s,Array{Planet}(0),Array{Orbit}(0))  # Constructor for a Planetary System with no planets

PlanetarySystem{StarT<:StarAbstract}(s::StarT, p::Planet, o::Orbit) = PlanetarySystem(s,[p],[o])  # Constructor for a single Planet System

function PlanetarySystem{StarT}(ps::PlanetarySystem{StarT}, keep::Vector{Int64}) # Why doesn't this work?
   PlanetarySystem{StarT}(ps.star,ps.planet[keep],ps.orbit[keep])
end

#function PlanetarySystemSingleStar(ps::PlanetarySystemSingleStar, keep::Vector{Int64})
#function PlanetarySystem{StarT<:StarAbstract}(ps::PlanetarySystem{StarT}, keep::Vector{Int64})
#   PlanetarySystem(ps.star,ps.planet[keep],ps.orbit[keep])
#end

flux(ps::PlanetarySystem{StarAbstract}) = flux(ps.star)
flux(ps::PlanetarySystem{Star}) = flux(ps.star)
flux(ps::PlanetarySystem{BinaryStar}) = flux(ps.star)
flux(ps::PlanetarySystem{MultipleStar}) = flux(ps.star)

function num_planets(s::PlanetarySystem{Star})
  @assert( length(s.planet) == length(s.orbit) )    # TODO OPT: Deactivate inner assert's like this for speed once tested
  return length(s.planet)
end

function generate_planet_mass_from_radius_powerlaw(r::Float64, sim_param::SimParam)
  const mr_power_index::Float64 = get_real(sim_param,"mr_power_index")
  const mr_const::Float64 = get_real(sim_param,"mr_const")
  const mr_max_mass::Float64 = get_real(sim_param,"mr_max_mass")
  m = mr_const*earth_mass*(r/earth_radius)^mr_power_index
  if m > mr_max_mass
     m = mr_max_mass
  end
  return m
end

function generate_planet_mass_from_radius_powerlaw(r::Float64, s::Star, o::Orbit, sim_param::SimParam)  # TODO SCI: IMPORTANT once have stability criteria, replace w/ better M-R relationship
  generate_planet_mass_from_radius_powerlaw(r,sim_param)
end

function generate_num_planets_poisson(lambda::Real, max_planets::Integer; min_planets::Integer = 0)
  if lambda < min_planets*1e-3
      return min_planets
  end
  bug_fixed = false # TODO OPT: true case should work, but Danley found bug in Distributions package.  Revert once fixed for speed.
  local n
  if bug_fixed
     d = Distributions.Truncated(Distributions.Poisson(lambda),min_planets,max_planets)
     n = rand(d)
  else
     if min_planets == 0
        min_planets = -1
     end
     d = Distributions.Truncated(Distributions.Poisson(lambda),min_planets,max_planets)
     n = rand(d)
     #=
     n = -1
     while !(min_planets<=n<=max_planets)
        n = rand(Distributions.Poisson(lambda))
     end
     =#
  end
  return n
end

function generate_num_planets_poisson(s::Star, sim_param::SimParam)
  const lambda::Float64 = exp(get_real(sim_param,"log_eta_pl"))
  const max_tranets_in_sys::Int64 = get_int(sim_param,"max_tranets_in_sys")
  generate_num_planets_poisson(lambda,max_tranets_in_sys)
end

function generate_period_and_sizes_log_normal(s::Star, sim_param::SimParam; num_pl::Integer = 1)  # TODO USER SCI:  User should make sure planetary properties are being drawn appropriately for their scientific purposes
    const mu_log_r::Float64 = get_real(sim_param,"mean_log_planet_radius")
    const sigma_log_r::Float64 = get_real(sim_param,"sigma_log_planet_radius")
    const mu_log_P::Float64 = get_real(sim_param,"mean_log_planet_period")
    const sigma_log_P::Float64 = get_real(sim_param,"sigma_log_planet_period")
    const min_period::Float64 = get_real(sim_param,"min_period")
    const max_period::Float64 = get_real(sim_param,"max_period")
    const min_radius::Float64 = get_real(sim_param,"min_radius")
    const max_radius::Float64 = get_real(sim_param,"max_radius")
    const max_draws::Int64 = 100

    if   sigma_log_r <= 0. || sigma_log_P<=0.
     println("# mu_log_r= ", mu_log_r, " sigma_log_r= ", sigma_log_r, " mu_log_P= ", mu_log_P, " sigma_log_P= ", sigma_log_P)
    end
    const rdist = LogNormal(mu_log_r,sigma_log_r)
    const Pdist = LogNormal(mu_log_P,sigma_log_P)
    #Rlist = rand(rdist,num_pl)
    #Plist = rand(Pdist,num_pl)
    #idx_keep = find(i->(min_radius<=Rlist[i]<=max_radius) && (min_period<=Plist[i]<=max_period), 1:num_pl )
    #return Plist[idx_keep], Rlist[idx_keep]  # replaced because want to return exactly num_pl.  Could use Truncated to restore above code.
    Rlist = zeros(num_pl)
    Plist = zeros(num_pl)
    for i in 1:num_pl
      j = 0
      while ! (min_radius<Rlist[i]<=max_radius) && j<max_draws
            Rlist[i] = rand(rdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw size for: ",mu_log_r, " ", sigma_log_r)
      end
      j = 0
      while ! (min_period<Plist[i]<=max_period) && j<max_draws
            Plist[i] = rand(Pdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw period for: ",mu_log_P, " ", sigma_log_P)
      end
    end
    return Plist, Rlist
end

function draw_power_law(n::Real, x0::Real, x1::Real, num_pl::Integer)
     ((x1^(n+1) - x0^(n+1)).*rand(num_pl) .+ x0^(n+1)).^(1/(n+1))
end

function draw_broken_power_law(n1::Real, n2::Real, x0::Real, x1::Real, xb::Real, num_pl::Integer)
    #x0 and x1 are the lower and upper truncation limits, and xb is the break point, i.e. x0 <= xb <= x1 (all must be positive)
    #n1 and n2 are the power law indices between x0 and xb, and xb and x1, respectively (can be positive or negative)
    @assert(x0 <= xb <= x1)
    @assert(num_pl >= 1)

    u_draws = rand(num_pl) #'num_pl' draws from the uniform distribution between 0 and 1
    x_draws = zeros(num_pl)

    if (n1 != -1) & (n2 != -1)
        C1 = 1./(((xb^(n1+1) - x0^(n1+1))/(n1+1)) + ((xb^(n1-n2)*(x1^(n2+1) - xb^(n2+1)))/(n2+1))) #normalization constant
        ub = (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = (((n1+1)*u)/C1 + x0^(n1+1))^(1/(n1+1))
            else #if u > ub
                x_draws[i] = (((n2+1)/(C1*xb^(n1-n2)))*(u - (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1)) + xb^(n2+1))^(1/(n2+1))
            end
        end
    elseif (n1 == -1) & (n2 != -1)
        C1 = 1./(log(xb/x0) + ((xb^(-1-n2))*(x1^(n2+1)) - 1)/(n2+1)) #normalization constant
        ub = C1*log(xb/x0) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = x0*exp(u/C1)
            else #if u > ub
                x_draws[i] = (((n2+1)/(C1*xb^(-1-n2)))*(u - C1*log(xb/x0)) + xb^(n2+1))^(1/(n2+1))
            end
        end
    elseif (n1 != -1) & (n2 == -1)
        C1 = 1./(((xb^(n1+1) - x0^(n1+1))/(n1+1)) + (xb^(n1+1))*log(x1/xb)) #normalization constant
        ub = (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1) #break point in u, between 0 and 1
        for (i,u) in enumerate(u_draws)
            if u <= ub
                x_draws[i] = (((n1+1)*u)/C1 + x0^(n1+1))^(1/(n1+1))
            else #if u > ub
                x_draws[i] = xb*exp((1/(C1*xb^(n1+1)))*(u - (C1*(xb^(n1+1) - x0^(n1+1)))/(n1+1)))
            end
        end
    else #if n1 == -1 and n2 == -1 (i.e. it is a single power-law with index of -1)
        for (i,u) in enumerate(u_draws)
            x_draws[i] = x0*exp(u*log(x1/x0))
        end
    end

    return x_draws
end

function generate_periods_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    const power_law_P::Float64 = get_real(sim_param,"power_law_P")
    const min_period::Float64 = get_real(sim_param,"min_period")
    const max_period::Float64 = get_real(sim_param,"max_period")
    Plist = power_law_P!=-1.0 ? draw_power_law(power_law_P,min_period,max_period, num_pl) : exp(log(min_period).+rand()*log(max_period/min_period))
    return Plist
end

function generate_sizes_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    const power_law_r::Float64 = get_real(sim_param,"power_law_r")
    const min_radius::Float64 = get_real(sim_param,"min_radius")
    const max_radius::Float64 = get_real(sim_param,"max_radius")
    Rlist = power_law_r!=-1.0 ? draw_power_law(power_law_r,min_radius,max_radius, num_pl) : exp(log(min_radius).+rand()*log(max_radius/min_radius))
    return Rlist
end

function generate_sizes_broken_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    const power_law_r1::Float64 = get_real(sim_param,"power_law_r1")
    const power_law_r2::Float64 = get_real(sim_param,"power_law_r2")
    const min_radius::Float64 = get_real(sim_param,"min_radius")
    const max_radius::Float64 = get_real(sim_param,"max_radius")
    const break_radius::Float64 = get_real(sim_param,"break_radius")
    Rlist = draw_broken_power_law(power_law_r1,power_law_r2,min_radius,max_radius,break_radius, num_pl)
    return Rlist
end

function generate_period_and_sizes_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    return (generate_periods_power_law(s, sim_param, num_pl=num_pl), generate_sizes_power_law(s, sim_param, num_pl=num_pl))
end

function TruncatedUpper(d::Distributions.UnivariateDistribution, u::Float64)
    zero(u) < u || error("lower bound should be less than upper bound.")
    lcdf = zero(u)
    ucdf = isinf(u) ? one(u) : cdf(d, u)
    tp = ucdf - lcdf
    Distributions.Truncated{typeof(d),Distributions.value_support(typeof(d))}(d, zero(u), u, lcdf, ucdf, tp, log(tp))
end

function generate_e_omega_rayleigh_direct(sigma_hk::Float64; max_e::Float64 = 1.0)
  @assert(0<max_e<=1.0)
  ecc::Float64 = rand( TruncatedUpper(Rayleigh(sigma_hk),max_e) )
  w::Float64 = 2pi*rand()
  return ecc, w
end

function generate_e_omega_rayleigh_two_gaussians(sigma_hk::Float64; max_e::Float64 = 1.0)
  @assert(0<max_e<=1.0)
  h = k = 1.0
  while h*h+k*k >= max_e*max_e
    h = sigma_hk*randn()
    k = sigma_hk*randn()
  end
  ecc::Float64 = sqrt(h*h+k*k)
  w::Float64 = atan2(k,h)
  return ecc, w
end

function generate_e_omega_rayleigh(sigma_hk::Float64; max_e::Float64 = 1.0)
   if max_e > sigma_hk
      return generate_e_omega_rayleigh_two_gaussians(sigma_hk,max_e=max_e)
   else
      return generate_e_omega_rayleigh_direct(sigma_hk,max_e=max_e)
   end
end

function generate_e_omega_rayleigh(sim_param::SimParam; max_e::Float64 = 1.0)
  sigma_hk::Float64 = get_real(sim_param,"sigma_hk")
  generate_e_omega_rayleigh(sigma_hk, max_e=max_e)
end

function generate_planetary_system_hardcoded_example(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # in this version we specify fixed functions that are known at compile time, allowing for additional optimizations (~0.6 second faster per Kepler catalog out of ~3.6 sec on my laptop w/ 1 core)
  const  generate_planet_mass_from_radius = generate_planet_mass_from_radius_powerlaw
  const  generate_num_planets = generate_num_planets_poisson
  const  generate_period_and_sizes = generate_period_and_sizes_log_normal

  const  generate_e_omega =  generate_e_omega_rayleigh

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl::Int64 = generate_num_planets(star, sim_param)

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
    (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(length(idx))
    orbit = Array{Orbit}(length(idx))
    a = map(i->semimajor_axis(Plist[i],star.mass),idx)
    max_e = ones(length(idx))
    max_e_factor = 0.999 # A factor just less than 1 to prevent numerical issues with near-crossing orbits
    if length(a)>=2
       for i in 1:length(a)
         if i==1
           max_e[i] = max_e_factor*(1-a[i]/a[i+1])/(1+a[i]/a[i+1])
         elseif i==length(a)
           max_e[i] = max_e_factor*(1-a[i-1]/a[i])/(1+a[i-1]/a[i])
         else
           max_e[i] = max_e_factor*min( (1-a[i]/a[i+1])/(1+a[i]/a[i+1]), (1-a[i-1]/a[i])/(1+a[i-1]/a[i]) )
         end
       end
    end
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param, max_e=max_e[i])
      incl::Float64 = acos(rand())
      orbit[i] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function generate_planetary_system_uncorrelated_incl(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
  const  generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
  const  generate_num_planets = get_function(sim_param,"generate_num_planets")
  const  generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
  const  generate_e_omega = get_function(sim_param,"generate_e_omega")

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl::Int64 = generate_num_planets(star, sim_param)::Int64
  const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)::PlanetarySystem
  else
    (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(length(idx))
    orbit = Array{Orbit}(length(idx))
    a = map(i->semimajor_axis(Plist[i],star.mass),idx)
    max_e = ones(length(idx))
    max_e_factor = 0.999 # A factor just less than 1 to prevent numerical issues with near-crossing orbits
    if length(a)>=2
       for i in 1:length(a)
          if i==1
            max_e[i] = max_e_factor*(1-a[i]/a[i+1])/(1+a[i]/a[i+1])
          elseif i==length(a)
            max_e[i] = max_e_factor*(1-a[i-1]/a[i])/(1+a[i-1]/a[i])
          else
            max_e[i] = max_e_factor*min( (1-a[i]/a[i+1])/(1+a[i]/a[i+1]), (1-a[i-1]/a[i])/(1+a[i-1]/a[i]) )
          end
        end
    end
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param, max_e=max_e[i])
      incl::Float64 = acos(rand())
      orbit[i] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      # set!(orbit[idx[i]],P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

# This version generates more systems roughly near a common plane, but until incorporate CORBITS data, ABC can't match input param
function generate_planetary_system_simple(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
  const  generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
  const  generate_num_planets = get_function(sim_param,"generate_num_planets")
  #const  generate_num_planets = generate_num_planets_christiansen
  #const  generate_period_and_sizes = generate_period_and_sizes_christiansen
  const  generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
  const  generate_e_omega = get_function(sim_param,"generate_e_omega")
  const  sigma_incl = deg2rad(get_real(sim_param,"sigma_incl"))

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl = generate_num_planets(star, sim_param)::Int64
  const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
    (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
    idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred
    incl_sys = acos(rand())

    min_a_in_rstar = 2.0
    min_P_orbit = day_in_year*sqrt((min_a_in_rstar*star.radius*rsol_in_au)^3 / star.mass) # minimum semi-major axis of two stellar radii
    idx = idx[Plist[idx] .> min_P_orbit]
    if( length(idx)==0 )
        return PlanetarySystem(star)
    end

    pl = Array{Planet}(length(idx))
    orbit = Array{Orbit}(length(idx))
    for i in 1:length(idx)
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      Rpl = Rlist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc,  omega) = generate_e_omega(sim_param)::Tuple{Float64,Float64}
      incl_mut = sigma_incl*sqrt(randn()^2+randn()^2) # rand(Distributions.Rayleigh(sigma_incl)) # sigma_incl*randn()
      asc_node = 2pi*rand()
      mean_anom = 2pi*rand()
      #incl = incl_sys + sigma_incl*randn()
      incl =  incl_mut!=zero(incl_mut) ? acos( cos(incl_sys)*cos(incl_mut) + sin(incl_sys)*sin(incl_mut)*cos(asc_node) ) : incl_sys
      orbit[i] = Orbit(P,ecc,incl,omega,asc_node,mean_anom)
      mass = generate_planet_mass_from_radius(Rpl, star, orbit[i], sim_param)::Float64
      pl[i] = Planet( Rpl,  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function test_planetary_system_constructors(sim_param::SimParam)
  generate_star = get_function(sim_param,"generate_star")
  star = generate_star(sim_param)
  empty_sys = PlanetarySystem(star)
  earth = Planet(earth_radius,earth_mass)
  earth_orbit = Orbit(365.2425,0.0167,0.5*pi,0.0,0.0,0.0)
  solar_sys = PlanetarySystem(star, earth,earth_orbit)
  m = generate_planet_mass_from_radius_powerlaw(0.02,star,earth_orbit,sim_param)/earth_mass
  generate_planetary_system_simple(star,sim_param,verbose=true)
end


text:
using IndEco
using Test

@testset "IndEco.jl" begin
    # Write your tests here.
end

ext: jl

max_stars_repo_name: MathAero/IndEco

max_issues_repo_name: MathAero/IndEco

max_forks_repo_name: MathAero/IndEco


text:
#=
 lhfs
 Description:
    Computes value and possibly source-derivatives of fundamental solution at
    a series of target points and energies, with the source at the origin.
 Parameters:
    x       (input) - array of x coordinates in the cartesian plane
    y       (input) - array of y coordinates in the cartesian plane
    energies(input) - parameter array to fundamental solution, one for each target pt
    derivs  (input) - 1: compute values and derivatives to fundamental solution.
                      0: only compute the fundamental solution values.
    n       (input) - size of x, y arrays
    u       (output) - array of size n of fundamental solution values.
            u[i] is the fundamental solution at the coordinate (x[i],y[i]).
    ux      (output) - array of the x derivative of each fundamental solution value
            (unused if derivs=0)
    uy      (output) - array of the y derivative of each fundamental solution value
            (unused if derivs=0)
    stdquad (input, int)     - global convergence params: # default quad pts per saddle pt
     h (input, double)        - PTR spacing relative to std Gaussian sigma
                           (if h<0, then -h is used and h is also used to
			   scale the parameters maxh and minsaddlen)
     meth (input) - int giving method for choosing # nodes (passed to quad_nodes)
     gamout (input) - 0: no file output (default), 1: diagnostic text to file nodes.dat
=#
function lhfs(x::Float64, y::Float64, energies::Float64, derivs::Bool, stdquad::Int, h::Float64, meth::Int)
    a = 0.25(x^2+y^2)
    b = 0.5y+energies
    if a == 0
        if derivs
            u = Inf+ZIM
            ux = Inf+ZIM
            uy = Inf+ZIM
            return u,ux,uy
        else
            u = Inf+ZIM
            return u
        end
    else
        nquad = quad_nodes!(ts, ws, a, b, x, y, derivs, stdquad, h, meth)
        contour!(gam, gamp, a, b, ts, nquad)
        integrand!(integ, integx, integy, gam, a, b, x, y, derivs, nquad)
        if derivs
            u,ux,uy = quadsum3(integ,integx,integy,gamp,ws,nquad)
            return u,ux,uy
        else
            u = quadsum(integ,gamp,ws,nquad)
            return u
        end
    end
end

function lhfs!(u::Vector{Complex{Float64}}, x::Vector{Float64}, y::Vector{Float64}, energies::Vector{Float64}, derivs::Bool, stdquad::Int, h::Float64, meth::Int, n::Int)
    for i=1:n
        a = 0.25(x[i]^2+y[i]^2)
        b = 0.5y[i]+energies[i]
        if a == 0
            u[i] = Inf+ZIM
        else
            nquad = quad_nodes!(ts, ws, a, b, x[i], y[i], derivs, stdquad, h, meth)
            contour!(gam, gamp, a, b, ts, nquad)
            integrand!(integ, integx, integy, gam, a, b, x[i], y[i], derivs, nquad)
            quadsum!(u,i,integ,gamp,ws,nquad)
        end
    end
end

function lhfs!(u::Vector{Complex{Float64}}, ux::Vector{Complex{Float64}}, uy::Vector{Complex{Float64}}, x::Vector{Float64}, y::Vector{Float64}, energies::Vector{Float64}, derivs::Bool, stdquad::Int, h::Float64, meth::Int, n::Int)
    for i=1:n
        a = 0.25(x[i]^2+y[i]^2)
        b = 0.5y[i]+energies[i]
        if a == 0
            u[i] = Inf+ZIM
            ux[i] = Inf+ZIM
            uy[i] = Inf+ZIM
        else
            nquad = quad_nodes!(ts, ws, a, b, x[i], y[i], derivs, stdquad, h, meth)
            contour!(gam, gamp, a, b, ts, nquad)
            integrand!(integ, integx, integy, gam, a, b, x[i], y[i], derivs, nquad)
            quadsum3!(u,ux,uy,i,integ,integx,integy,gamp,ws,nquad)
        end
    end
end

function quadsum!(u::Vector, i::Int, integ::Vector, gamp::Vector, ws::Vector, n::Int)
    @inbounds u[i] = integ[1]*gamp[1]*ws[1]
    @simd for k=2:n
        @inbounds u[i] += integ[k]*gamp[k]*ws[k]
    end
end

function quadsum3!(u::Vector, ux::Vector, uy::Vector, i::Int, integ::Vector, integx::Vector, integy::Vector, gamp::Vector, ws::Vector, n::Int)
    @inbounds temp = gamp[1]*ws[1]
    @inbounds u[i] = integ[1]*temp
    @inbounds ux[i] = integx[1]*temp
    @inbounds uy[i] = integy[1]*temp
    @simd for k=2:n
        @inbounds temp = gamp[k]*ws[k]
        @inbounds u[i] += integ[k]*temp
        @inbounds ux[i] += integx[k]*temp
        @inbounds uy[i] += integy[k]*temp
    end
end

function quadsum(integ::Vector, gamp::Vector, ws::Vector, n::Int)
    @inbounds ret = integ[1]*gamp[1]*ws[1]
    @simd for k=2:n
        @inbounds ret += integ[k]*gamp[k]*ws[k]
    end
    ret
end

function quadsum3(integ::Vector, integx::Vector, integy::Vector, gamp::Vector, ws::Vector, n::Int)
    @inbounds temp = gamp[1]*ws[1]
    @inbounds ret1 = integ[1]*temp
    @inbounds ret2 = integx[1]*temp
    @inbounds ret3 = integy[1]*temp
    @simd for k=2:n
        @inbounds temp = gamp[k]*ws[k]
        @inbounds ret1 += integ[k]*temp
        @inbounds ret2 += integx[k]*temp
        @inbounds ret3 += integy[k]*temp
    end
    ret1,ret2,ret3
end


text:
export monomialbasis, ideal

abstract type AbstractPolynomialIdeal end

struct EmptyPolynomialIdeal <: AbstractPolynomialIdeal end
ideal(p::FullSpace, algo=defaultgröbnerbasisalgorithm(p)) = EmptyPolynomialIdeal()
Base.rem(p::AbstractPolynomialLike, I::EmptyPolynomialIdeal) = p

mutable struct PolynomialIdeal{T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm} <: AbstractPolynomialIdeal
    p::Vector{PT}
    gröbnerbasis::Bool
    algo::A
end
function PolynomialIdeal{T, PT}(p::Vector{PT}, algo::A) where {T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm}
    PolynomialIdeal{T, PT, A}(p, false, algo)
end
function PolynomialIdeal(p::Vector{PT}, algo) where {T, PT<:APL{T}}
    S = Base.promote_op(/, T, T)
    PolynomialIdeal{S, polynomialtype(PT, S)}(AbstractVector{polynomialtype(PT, S)}(p), algo)
end
function PolynomialIdeal{T, PT}() where {T, PT<:APL{T}}
    PolynomialIdeal(PT[], defaultgröbnerbasisalgorithm(PT))
end

function Base.convert(::Type{PolynomialIdeal{T, PT, A}}, I::PolynomialIdeal{T, PT, A}) where {T, PT<:APL{T}, A<:AbstractGröbnerBasisAlgorithm}
    return I
end
function Base.convert(::Type{PolynomialIdeal{T, PT, A}}, I::PolynomialIdeal) where {T, PT, A}
    return PolynomialIdeal{T, PT, A}(I.p, I.gröbnerbasis, I.algo)
end

ideal(p, algo=defaultgröbnerbasisalgorithm(p)) = PolynomialIdeal(p, algo)

Base.:+(I::PolynomialIdeal, J::PolynomialIdeal) = PolynomialIdeal([I.p; J.p], I.algo)

MP.variables(I::PolynomialIdeal) = variables(I.p)

function Base.rem(p::AbstractPolynomialLike, I::PolynomialIdeal)
    computegröbnerbasis!(I)
    rem(p, I.p)
end

function computegröbnerbasis!(I::PolynomialIdeal)
    if !I.gröbnerbasis
        gröbnerbasis!(I.p, I.algo)
        I.gröbnerbasis = true
    end
end
function monomialbasis(I, vars=variables(I))
    computegröbnerbasis!(I)
    if isempty(I.p)
        return false, monomialtype(eltype(I.p))[]
    end
    mv = monovec(leadingmonomial.(I.p))
    # monovec makes sure all monomials have the same variables
    # if x_i^n is in the leading monomials, lv[i] <= n
    lv = -ones(Int, length(vars))
    for m in mv
        d = degree(m)
        if d == maximum(exponents(m))
            # univariate monomial
            if iszero(d)
                lv[:] .= 0
            else
                lv[argmax(exponents(m))] = d
            end
        end
    end
    if any(lv .< 0)
        return false, monomialtype(eltype(I.p))[]
    end
    return true, monomials(vars, 0:(sum(lv)-1), m -> !any(map(m2 -> divides(m2, m), mv)))
end


text:
export
    FoldSet,

    FOLD_TRAIN,
    FOLD_TEST,

    foldset_match,
    foldset_withhold,

    check_fold_match

#########################################

const FOLD_TRAIN = 1
const FOLD_TEST = 2

immutable FoldSet
    assignment::Vector{Int} # assignment[i] = j means the ith element is assigned to fold j
    fold::Int # target fold we either want to match or withhold, based on `match_fold`
    match_fold::Bool
end
foldset_match(assignment::Vector{Int}, fold::Int) = FoldSet(assignment, fold, true)
foldset_withhold(assignment::Vector{Int}, fold::Int) = FoldSet(assignment, fold, false)

function check_fold_match(fold::Integer, fold_assignment::Integer, match_fold::Bool)
    (fold != 0) && (fold_assignment != 0) && # NOTE(tim): zero never matches
        ((match_fold && fold_assignment == fold) || (!match_fold && fold_assignment != fold))
end
check_fold_match(value::Integer, foldset::FoldSet) = check_fold_match(foldset.fold, value, foldset.match_fold)

#########################################
# Iterator

function _find_next_valid_fold_match(foldset::FoldSet, state::Int)
    while state < length(foldset.assignment)
        state += 1
        if check_fold_match(foldset.assignment[state], foldset)
            return state
        end
    end
    state + 1 # returns length(foldset.assignment) + 1 on fail
end
Base.start(foldset::FoldSet) = _find_next_valid_fold_match(foldset, 0)
Base.done(foldset::FoldSet, state::Int) = state > length(foldset.assignment)
function Base.next(foldset::FoldSet, state::Int)
    @assert(check_fold_match(foldset.assignment[state], foldset))
    state, _find_next_valid_fold_match(foldset, state)
end
function Base.length(foldset::FoldSet)
    len = 0
    state = start(foldset)
    while !done(foldset, state)
        item, state = next(foldset, state)
        len += 1
    end
    len
end
function Base.collect(foldset::FoldSet)
    len = length(foldset)
    retval = Array{Int}(len)
    for (i,j) in enumerate(foldset)
        retval[i] = j
    end
    retval
end


text:
# Statistical computations for DischargeModels

nash_sutcliffe(M::DischargeModel) = 1-sum(abs2,residuals(M))/sum(abs2,model_response(M)-mean(model_response(M)))
nash_sutcliffe(M::DischargeModel,H,Q) = 1-sum(abs2,residuals(M,H,Q))/sum(abs2,Q-mean(Q))
function nash_sutcliffe(M::DischargeModel,p::Stage,dd::Discharge)
    Hm,Qm = preparedata(p,dd,1:length(dd),M.M)
    nash_sutcliffe(M,Hm,Qm)
end

function simpsonweights(N)
    w = zeros(Int,N)
    w[1] = 1
    w[end] = 1
    for j in 1:div(N-1,2)-1
        w[2j+1] = 2
        w[2j] = 4
    end
    w[N-1] = 4
    w
end

"""
Calculate the water balance for `Q` sampled at an interval of `dt` over the period `range`
"""
waterbalance(Q,dt,range) = dot(dt/3.*Q[r],simpsonweights(length(r)))

"""
Calculate the value of a Ljung-Box statistic for a model
"""
function ljung_box(M::DischargeModel,H,Q,p)
    rs = residuals(M,H,Q)
    n = length(rs)
    ρ = autocor(rs,1:p)
    n*(n+2)*sum((ρ.^2)./[n-k for k in 1:p])
end

ljung_box(M::DischargeModel,p) = ljung_box(M,M.H,M.Q,p)

function χ2test(Q,df)
    pdf(Chisq(df),Q)
end

"""
Calculate the p-value of a Ljung-Box test for a discharge model
"""
ljung_box_test(M::DischargeModel,H,Q,p) = χ2test(ljung_box(M,H,Q,p),p-M.M)

ljung_box_test(M::DischargeModel,p) = ljung_box_test(M,M.H,M.Q,p)

"""
Calculate the spectral flatness of a signal
"""
function flatness(x)
    psd = abs2.(fft(x))
    geomean(psd)/mean(psd)
end

flatness(M::DischargeModel,H,Q) = flatness(residuals(M,H,Q))
function flatness(M::DischargeModel,p::Stage,dd::Discharge)
    Hm,Qm = preparedata(p,dd,1:length(dd),M.M)
    flatness(M,Hm,Qm)
end

"""
Calculate the prediction error of a model
"""
pe(M::DischargeModel,H,Q) = var(residuals(M,H,Q))
pe(M::DischargeModel) = var(residuals(M))

"""
Calculate the coefficient of determination
"""
function r2(M::DischargeModel,H,Q)
    SStot = sum(abs2,Q-mean(Q))
    SSres = sum(abs2,residuals(M,H,Q))
    1-SSres/SStot
end

r2(M::DischargeModel) = r2(M,M.H,M.Q)


text:
"""
Julia script for kinematics analysis of four-bar mechanism.
Contains sections for position, velocity, and acceleration analyses.
"""

using Plots
plotlyjs()

function ex3()

# fixed parameters
AB=5
BC=6
CD=5
DA=6
DP=10
θ₁=π

#time limit
t=0:0.01:45

# position analysis

#independent coordinate
ζ₀=π/6 # angle subtended by the A-frame [rad]
#case 1
# β₀=dζ=90*π/(180*45*45)
# θ₂=ζ=ζ₀.+β₀.*t.^2
#case 2
β₀=dζ=90*π/(180*45)
θ₂=ζ=ζ₀.+β₀.*t

x=-AB.*cos.(θ₂).-DA.*cos.(θ₁)
y=-AB.*sin.(θ₂).-DA.*sin.(θ₁)

z=(x.^2 .+ y.^2 .+ CD^2 .- BC^2)./(2 .* CD)

θ₄=acos.(-z./sqrt.(x.^2 .+ y.^2)) .+ atan.(y,x)

y_C=CD.*sin.(θ₄)
y_B=AB.*sin.(θ₂)

x_C=CD.*cos.(θ₄).+DA
x_B=AB.*cos.(θ₂)

θ₃=atan.(y_C .- y_B,x_C .- x_B)

#plotting
p1=plot(t,θ₂,linewidth=2,xaxis="time (s)",yaxis="θ (rad)",
    label="θ₂",color="blue")
p2=plot(t,θ₃,linewidth=2,xaxis="time (s)",yaxis="θ (rad)",
    label="θ₃",color="purple",ylim=[-0.4,0.4])
p3=plot(t,θ₄,linewidth=2,xaxis="time (s)",yaxis="θ (rad)",
    label="θ₄",color="orange")
display(plot(p1,p2,p3,layout=(3,1)))
savefig("./plots/ex3-pb1.png")

#trajectory of point of suspension
xₚ=DA .+ DP.*cos.(θ₄)
yₚ=DP.*sin.(θ₄)
display(plot(xₚ,yₚ,label="trajectory",xaxis="x (m)",yaxis="y (m)",
    linewidth=2,xlim=[0,15],xticks=0:1:15,ylim=[0,15],yticks=0:1:15))
savefig("./plots/ex3-trajectory.png")

#velocity analysis

#independent coordinate
# case 1
# θ̇₂=dθ₂=2 .*β₀.*t
# case 2
θ̇₂=dθ₂=β₀.*ones(length(t))

X=ones(2,length(t))

#solving the equations
for i in 1:length(t)
    A=[BC.*sin.(θ₃[i])  -CD.*sin.(θ₄[i]); BC.*cos.(θ₃[i])  -CD.*cos.(θ₄[i])]
    # X=[θ̇₃;θ̇₄]
    B=[-AB.*θ̇₂[i].*sin.(θ₂[i]); -AB.*θ̇₂[i].*cos.(θ₂[i])]

    Aᴵ=inv(A)
    @assert size(Aᴵ,2)==size(B,1) "dimension mismatch"
    X[:,i]=Aᴵ*B
end

(θ̇₃,θ̇₄)=[X[i,:] for i in 1:size(X,1)]

#plotting
p4=plot(t,θ̇₂,linewidth=2,xaxis="time (s)",yaxis="ω (rad/s)",
    label="ω₂",color="blue",ylim=[0,0.1])
p5=plot(t,θ̇₃,linewidth=2,xaxis="time (s)",yaxis="ω (rad/s)",
    label="ω₃",color="purple",ylim=[-0.4,0.4])
p6=plot(t,θ̇₄,linewidth=2,xaxis="time (s)",yaxis="ω (rad/s)",
    label="ω₄",color="orange",ylim=[0,0.1])
display(plot(p4,p5,p6,layout=(3,1)))
savefig("./plots/ex3-pb2.png")

#acceleration analysis

#independent coordinate
#case 1
# θ̈₂=dθ̇₂=2*β₀.*ones(length(t))
#case 2
θ̈₂=dθ̇₂=zeros(length(t))

Y=ones(2,length(t))

#solving the equations
for i in 1:length(t)
    C=[BC*cos(θ₃[i])  -CD*cos(θ₄[i]); BC*sin.(θ₃[i])  -CD*sin.(θ₄[i])]
    # Y=[θ̈₃;θ̈₄]
    D=[(-AB*θ̈₂[i]*cos.(θ₂[i])+AB*(θ̇₂[i]^2)*sin(θ₂[i])+BC*(θ̇₃[i]^2)*sin.(θ₃[i])
            -CD*(θ̇₄[i]^2)*sin(θ₄[i]));
        (-AB*θ̈₂[i]*sin(θ₂[i])-AB*(θ̇₂[i]^2)*cos(θ₂[i])-BC*(θ̇₃[i]^2)*cos.(θ₃[i])
            +CD*(θ̇₄[i]^2)*cos(θ₄[i]))]

            Cᴵ=inv(C)
    @assert size(Cᴵ,2)==size(D,1) "dimension mismatch"
    Y[:,i]=Cᴵ*D
end

(θ̈₃,θ̈₄)=[Y[i,:] for i in 1:size(Y,1)]

#plotting
p7=plot(t,θ̈₂,linewidth=2,xaxis="time (s)",yaxis="α (rad/s^2)",
    label="α₂",color="blue",ylim=[0,0.0025])
p8=plot(t,θ̈₃,linewidth=2,xaxis="time (s)",yaxis="α (rad/s^2)",
    label="α₃",color="purple",ylim=[-0.4,0.4])
p9=plot(t,θ̈₄,linewidth=2,xaxis="time (s)",yaxis="α (rad/s^2)",
    label="α₄",color="orange",ylim=[0,0.0025])
display(plot(p7,p8,p9,layout=(3,1)))
savefig("./plots/ex3-pb3.png")

K=(θ₂,θ₃,θ₄,θ̇₂,θ̇₃,θ̇₄,θ̈₂,θ̈₃,θ̈₄)

return K

end

max_stars_repo_path: Dynamic analysis of the A-frame mechanism/scripts/ex3.jl


text:
"""
Type for correlated squared Euclidean
"""
mutable struct P_misfit_xcorr{T<:Real}
	pxcorr::Conv.P_xcorr{T}
	cy::Vector{Matrix{T}}
	dcg::Vector{Matrix{T}}
end

function P_misfit_xcorr(nt::Int, nr::Int,
			pxcorr=P_xcorr(nt, nr, norm_flag=false); y=nothing, cy=nothing)

	if(cy===nothing)
		cy=deepcopy(pxcorr.cg) # allocate cy
		(y===nothing) && error("need y or cy")
		(size(y)≠size(pxcorr.g)) && error("size y")

		Conv.mod!(pxcorr, cg=cy, g=y)
	end
	dcg=deepcopy(pxcorr.cg)
	fill!.(dcg, 0.0)
	return P_misfit_xcorr(pxcorr, cy, dcg)
end

function func_grad!(dfdx,  x,  pa::P_misfit_xcorr)
	nt=size(x,1)
	nr=size(x,2)
	cg=pa.pxcorr.cg
	cy=pa.cy
	dfdcg=pa.dcg

	# copy x
	copyto!(pa.pxcorr.g,x)

	# mod xcorr
	mod!(pa.pxcorr)

	J=0.0
	for ir in 1:length(cg)
		cgx=cg[ir]
		cyy=cy[ir]
		dfdcgx=dfdcg[ir]

		if(dfdx===nothing)
			JJ=Misfits.error_squared_euclidean!(nothing,  cgx,   cyy,   nothing, norm_flag=false)
		else
			JJ=Misfits.error_squared_euclidean!(dfdcgx,  cgx,   cyy,   nothing, norm_flag=false)
		end
		J+=JJ
	end
	if(!(dfdx===nothing))
		mod_grad!(dfdx, pa.pxcorr, dcg=dfdcg)
	end

	return J
end

mutable struct P_misfit_weighted_acorr{T<:Real}
	p_conv::Conv.P_conv{T,2,2,2}
	ds::Matrix{Float64}
end

function P_misfit_weighted_acorr(nt::Int, nr::Int)
	p_conv=P_conv(gsize=[nt,nr], dsize=[nt,nr], ssize=[2*nt-1,nr], slags=[nt-1, nt-1])
	ds=zero(p_conv.s)
	return P_misfit_weighted_acorr(p_conv, ds)
end

"""
Each column of the matrix x is correlated with itself,
the energy at non-zero las will be penalized and returned as J.
"""
function func_grad!(dfdx, x, pa::P_misfit_weighted_acorr)
	nt=size(x,1)
	nr=size(x,2)

	copyto!(pa.p_conv.g,x)
	copyto!(pa.p_conv.d,x)

	Conv.mod!(pa.p_conv, S())
	s=pa.p_conv.s
	J=0.0
	for ir in 1:nr
		for it in 1:size(s,1)
			J += abs2(s[it,ir]) * abs2((nt-it)/(nt-1))
		end
	end

	if(!(dfdx===nothing))
		fill!(pa.ds, 0.0)
		for ir in 1:nr
			for it in 1:size(s,1)
				pa.ds[it,ir] = 2.0 * (s[it,ir]) * abs2((nt-it)/(nt-1))
			end
		end

		Conv.mod!(pa.p_conv, G(), g=dfdx, s=pa.ds)

		rmul!(dfdx, 2.)
	end

	return J

end


text:
module Day20

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day20"), String)

export part1, part2

const pixel_map = Dict(
    "." => Int8(0),
    "#" => Int8(1)
)

read_algorithm(algorithm::AbstractString, hash_value::Int) = pixel_map[string(algorithm[hash_value+1])]

function parse_input(input::AbstractString)
    lines = split(input, "\n", keepempty = false)
    algorithm = popfirst!(lines)
    arrayed = [(x -> pixel_map[x]).(split(line, "")) for line in lines]
    image = Matrix(transpose(hcat(arrayed...))) # turn into a Matrix with this ugly code
    return algorithm, image
end

function out_of_bounds(m::Matrix{Int8}, row::Int, column::Int)
    nrows, ncols = size(m)
    row < 1 || row > nrows || column < 1 || column > ncols
end

function pixel(m::Matrix{Int8}, row::Int, column::Int, default::Int8)
    if out_of_bounds(m, row, column)
        return default
    end
    return m[row, column]
end

function neighbourhood(m::Matrix{Int8}, row::Int, column::Int, default::Int8)
    [pixel(m, i, j, default) for i = row-1:row+1 for j = column-1:column+1]
end

function neighbourhood_hash(m::Matrix{Int8}, row::Int, column::Int, default::Int8)
    neighbours = neighbourhood(m, row, column, default)
    neighbourhood_hash_string = join(string.(neighbours))
    parse(Int64, neighbourhood_hash_string, base = 2)
end

function add_border(image::Matrix{Int8}, pad_with::Int8 = Int8(0), border_size::Int = 2)
    pad_function = Bool(pad_with) ? ones : zeros
    nrows, ncols = size(image)
    vcat(
        pad_function(Int8, border_size, ncols + 2 * border_size),
        [pad_function(Int8, nrows, border_size) image pad_function(Int8, nrows, border_size)],
        pad_function(Int8, border_size, ncols + 2 * border_size)
    )
end

function enhance(algorithm::AbstractString, image::Matrix{Int8}, infinity::Int8 = Int8(0))
    padded = add_border(image, infinity)
    resulting_matrix = zeros(Int8, size(padded)...)
    nrows, ncols = size(padded)
    for row = 1:nrows
        for column = 1:ncols
            hash_value = neighbourhood_hash(padded, row, column, infinity)
            resulting_matrix[row, column] = read_algorithm(algorithm, hash_value)
        end
    end
    infinity_hash_value = parse(Int64, repeat(string(infinity), 9), base = 2)
    new_infinity = read_algorithm(algorithm, infinity_hash_value)
    return resulting_matrix, new_infinity
end

function part1(input = input)
    algorithm, image = parse_input(input)
    infinity = Int8(0)
    for i = 1:2
        image, infinity = enhance(algorithm, image, infinity)
    end
    sum(image)
end

function part2(input = input)
    algorithm, image = parse_input(input)
    infinity = Int8(0)
    for i = 1:50
        image, infinity = enhance(algorithm, image, infinity)
    end
    sum(image)
end

end #module


text:
using Rayden, Scallop

thick.behind_mirror = Inf*μm
ri.behind_mirror = 0.0
thick.water = 0μm

function format_radii{T <: Number}(x::Tuple{T})
    xx1, = x
    x1 = format(xx1)
    return Vec(x1, x1, x1)
end
function format_radii{T <: Number}(x::Tuple{T, T})
    xx1, xx2 = x
    x1 = format(xx1)
    x2 = format(xx2)
    return Vec(x1, x1, x2)
end
function format_radii{T <: Number}(xx::T)
    x = format(xx)
    return Vec(x, x, x)
end

_format(x::T) where T <: Number = uconvert(rad, x)
_format(x::T) where T <: Unitful.Length = uconvert(μm, x)
format(x::T) where T <: Number = Float64(ustrip(_format(x)))
format(x::Int) = Float64(x)
format(x::Float64) = x

morphing_factor = format(morphing_factor)
source_distance = format(source_distance)

volumes = Dict(String(k) => Volume(Float64(getfield(ri, k)), format(getfield(thick, k))) for k in fieldnames(ri))
planes = Dict(String(k) => Plane(format_radii(getfield(radii, k)), !r"distal$"(String(k))) for k in fieldnames(radii))
aperture = format(aperture)

ellipsoids = getmembranes(volumes, planes, morphing_factor, aperture)
opticunits = scallop(ellipsoids, volumes)
l = getLight(source_distance, aperture, ellipsoids["cornea_distal"], format(source_angle))

min_aperture = format(min_aperture)
max_aperture = format(max_aperture)
min_distance = format(min_distance)
max_distance = format(max_distance)

include("plottingFunctions.jl")


text:
function changes(amount::Int, coins::Array{Int})::Int128
    ways = zeros(Int128, amount + 1)
    ways[1] = 1
    for coin in coins, j in coin+1:amount+1
        ways[j] += ways[j - coin]
    end
    return ways[amount + 1]
end

@show changes(100, [1, 5, 10, 25])
@show changes(100000, [1, 5, 10, 25, 50, 100])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
$(DocStringExtensions.README)
"""
module MathML
using DocStringExtensions

using EzXML, Symbolics, Statistics, IfElse, AbstractTrees

include("utils.jl")
include("parse.jl")
include("maps.jl")

export extract_mathml, mathml_to_nums, @xml_str, @MathML_str
export parse_node, parse_apply, parse_file, parse_str

export parse_cn, parse_ci, parse_bvar, parse_lambda, parse_piecewise
end


text:
# ---
# title: 13. Roman to Integer
# id: problem13
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math, String
# link: <https://leetcode.com/problems/roman-to-integer/description/>
# hidden: true
# ---
#
# Roman numerals are represented by seven different symbols: `I`, `V`, `X`, `L`,
# `C`, `D` and `M`.
#
#
#
#     **Symbol**       **Value**
#     I             1
#     V             5
#     X             10
#     L             50
#     C             100
#     D             500
#     M             1000
#
# For example, `2` is written as `II` in Roman numeral, just two one's added
# together. `12` is written as `XII`, which is simply `X + II`. The number `27`
# is written as `XXVII`, which is `XX + V + II`.
#
# Roman numerals are usually written largest to smallest from left to right.
# However, the numeral for four is not `IIII`. Instead, the number four is
# written as `IV`. Because the one is before the five we subtract it making
# four. The same principle applies to the number nine, which is written as `IX`.
# There are six instances where subtraction is used:
#
#   * `I` can be placed before `V` (5) and `X` (10) to make 4 and 9.
#   * `X` can be placed before `L` (50) and `C` (100) to make 40 and 90.
#   * `C` can be placed before `D` (500) and `M` (1000) to make 400 and 900.
#
# Given a roman numeral, convert it to an integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "III"
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: s = "IV"
#     Output: 4
#
#
# **Example 3:**
#
#
#
#     Input: s = "IX"
#     Output: 9
#
#
# **Example 4:**
#
#
#
#     Input: s = "LVIII"
#     Output: 58
#     Explanation: L = 50, V= 5, III = 3.
#
#
# **Example 5:**
#
#
#
#     Input: s = "MCMXCIV"
#     Output: 1994
#     Explanation: M = 1000, CM = 900, XC = 90 and IV = 4.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 15`
#   * `s` contains only the characters `('I', 'V', 'X', 'L', 'C', 'D', 'M')`.
#   * It is **guaranteed**  that `s` is a valid roman numeral in the range `[1, 3999]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using OpenIDConnect
using Test
using Random
using HTTP

function test_state_store()
    ctx = OIDCCtx("https://accounts.google.com", "http://127.0.0.1:8888/auth/login", "test_client_id", "test_client_secret"; state_timeout_secs=5)
    state = randstring(10)
    OpenIDConnect.remember_state(ctx, state)
    @test length(ctx.states) == 1
    @test OpenIDConnect.validate_state(ctx, state)
    sleep(10)
    @test !OpenIDConnect.validate_state(ctx, state)
    @test length(ctx.states) == 0
    nothing
end

function test_oidc_flow()
    ctx = OIDCCtx("https://accounts.google.com", "http://127.0.0.1:8888/auth/login", "test_client_id", "test_client_secret"; state_timeout_secs=5)
    @test OpenIDConnect.authorization_endpoint(ctx) == "https://accounts.google.com/o/oauth2/v2/auth"
    @test OpenIDConnect.token_endpoint(ctx) == "https://oauth2.googleapis.com/token"

    # flow request authorization code
    uri_string = flow_request_authorization_code(ctx)
    uri = HTTP.URIs.URI(uri_string)
    @test uri.host == "accounts.google.com"

    query = HTTP.URIs.queryparams(uri)
    @test get(query, "client_id", "") == "test_client_id"
    @test get(query, "redirect_uri", "") == "http://127.0.0.1:8888/auth/login"
    @test get(query, "scope", "") == "openid profile email"
    @test get(query, "response_type", "") == "code"
    @test !isempty(get(query, "state", ""))

    uri_string = flow_request_authorization_code(ctx; nonce="test_nonce", display="test_display", prompt="test_prompt", max_age="12345", ui_locales="en", id_token_hint="test_id_tok_hint", login_hint="test_login_hint", acr_values="test_acr")
    uri = HTTP.URIs.URI(uri_string)
    @test uri.host == "accounts.google.com"

    query = HTTP.URIs.queryparams(uri)
    @test get(query, "client_id", "") == "test_client_id"
    @test get(query, "redirect_uri", "") == "http://127.0.0.1:8888/auth/login"
    @test get(query, "scope", "") == "openid profile email"
    @test get(query, "response_type", "") == "code"
    @test !isempty(get(query, "state", ""))
    @test get(query, "nonce", "") == "test_nonce"
    @test get(query, "display", "") == "test_display"
    @test get(query, "prompt", "") == "test_prompt"
    @test get(query, "max_age", "") == "12345"
    @test get(query, "ui_locales", "") == "en"
    @test get(query, "id_token_hint", "") == "test_id_tok_hint"
    @test get(query, "login_hint", "") == "test_login_hint"
    @test get(query, "acr_values", "") == "test_acr"

    # flow get authorization code
    @test isa(flow_get_authorization_code(ctx, Dict()), OpenIDConnect.APIError)
    @test isa(flow_get_authorization_code(ctx, Dict("state"=>"teststate")), OpenIDConnect.APIError)
    OpenIDConnect.remember_state(ctx, "teststate")
    @test isa(flow_get_authorization_code(ctx, Dict("state"=>"teststate")), OpenIDConnect.APIError)
    @test isa(flow_get_authorization_code(ctx, Dict("state"=>"teststate", "error"=>"testerror")), OpenIDConnect.AuthServerError)
    @test "testcode" == flow_get_authorization_code(ctx, Dict("state"=>"teststate", "code"=>"testcode"))
end

test_state_store()
test_oidc_flow()


text:
using GeometryTypes
import ImmutableArrays
using Benchmarks

macro mybench(expr)
	quote
		val = $(esc(expr))
		stats = Statistics(@benchmark for i=1:10
			val = $(esc(expr))
		end)
		println(stats.average_time, ": ", stats.interval[1], " - ", stats.interval[2])
		val
	end
end
function test()
	const a = Float64[1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]
	const b = Float64[1,2,3,4]
	const a2 = Mat{4,4,Float64}((
		(1,2,3,4),
		(1,2,3,4),
		(1,2,3,4),
		(1,2,3,4)
	))
	const b2 = Vec{4, Float64}(1,2,3,4)
	const a3 = ImmutableArrays.Matrix4x4{Float64}(
		ImmutableArrays.Vector4{Float64}(1,2,3,4),
		ImmutableArrays.Vector4{Float64}(1,2,3,4),
		ImmutableArrays.Vector4{Float64}(1,2,3,4),
		ImmutableArrays.Vector4{Float64}(1,2,3,4)
	)
	const b3 = ImmutableArrays.Vector4{Float64}(1,2,3,4)

	println("matmul 4x4 * 4: ")
	@mybench a*b
	@mybench a2*b2
	@mybench a3*b3

	println("matmul 4x4 * 4x4: ")
	@mybench a*a
	@mybench a2*a2
	@mybench a3*a3

	println("sum: ")
	@mybench sum(a)
	@mybench sum(a2)
	@mybench sum(a3)

	println("dot:")
	@mybench for i=1:100 dot(b,b) end
	@mybench for i=1:100 dot(b2,b2) end
	@mybench for i=1:100 dot(b3,b3) end

	println("column:")
	@mybench for i=1:1000 column(a2,1) end
	@mybench for i=1:1000 ImmutableArrays.column(a3,1) end

	println("row:")
	@mybench for i=1:1000 row(a2,1) end
	@mybench for i=1:1000 ImmutableArrays.row(a3,1) end
end
test()


text:
module TestTSML
using Test

include("test_tsclassifier.jl")
include("test_mlbase.jl")
include("test_valdate.jl")
include("test_statifier.jl")
include("test_monotonicer.jl")
include("test_outliernicer.jl")
include("test_ensemble.jl")
include("test_cliwrapper.jl")
include("test_normalizer.jl")
include("test_svm.jl")

# test if running windows
if !Base.Sys.iswindows()
  include("test_plotter.jl")
end

#include("test_timescaledb.jl")
#include("test_schemalizer.jl")

end


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using Test

include("apitest.jl")
include("testexamples.jl")
include("testshow.jl")
#include("test_mathprogbase.jl")
include("issues.jl")


text:
@testset "Rastrigin" begin

    function test_result(result::Vector, fitness::Float64, N::Int, tol::Float64)
        if round(fitness) == 0
            @test ≈(result, zeros(N), atol=tol)
            @test ≈(fitness,0.0, atol=tol)
        else
            # @warn("Found local minimum!!!")
            @test sum(round.(abs.(result))) < N
        end
    end

    # Objective function
    function rastrigin(x::AbstractVector{T}) where {T <: AbstractFloat}
        n = length(x)
        return 10n + sum([ x[i]^2 - 10cos(convert(T,2π*x[i])) for i in 1:n ])
    end

    # Parameters
    N = 3
    P = 100

    # Testing: (μ/μ_I,λ)-σ-Self-Adaptation-ES
    # with non-isotropic mutation operator y' := y + (σ_1 N_1(0, 1), ..., σ_N N_N(0, 1))
    result, fitness, cnt = optimize( rastrigin, ES(N=N;
            initStrategy = strategy(σ = .5ones(N), τ = 1/sqrt(2*N), τ0 = 1/sqrt(N)),
            recombination = average, srecombination = averageSigmaN,
            mutation = anisotropic, smutation = anisotropicSigma,
            selection=:comma,
            μ = 15, λ = P),iterations = 1000)
    println("(15/15,$(P))-σ-SA-ES-AS => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

    # Testing: CMA-ES
    result, fitness, cnt = optimize( rastrigin,CMAES(N=N; μ = 15, λ = P),tol = 1e-8)
    println("(15/15,$(P))-CMA-ES => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

    # Testing: GA
    selections = [roulette, sus, ranklinear(1.5)]
    crossovers = [discrete, intermediate(0.), intermediate(0.25), line(0.2)]
    mutations = [domainrange(fill(0.5,N)), domainrange(fill(1.0,N))]

    @testset "GA settings" for ss in selections, xovr in crossovers, ms in mutations
        result, fitness, cnt = optimize( rastrigin, GA(N=N;
                populationSize = P,
                ɛ = 0.1,
                selection = ss,
                crossover = xovr,
                mutation = ms),tol = 1e-5)
        println("GA(p=$(P),x=.8,μ=.1,ɛ=0.1) => F: $(fitness), C: $(cnt), OBJ: $(result)")
        test_result(result, fitness, N, 1e-1)
    end
end


text:
"""
Generate Julia wrapping code with type annotations,
by inspecting TensorFlow Python modules and docs
"""
module TFParser

export tfwritejulia

using PyCall

using ..PyInspector
using ..CoreTypes

@pyimport tensorflow as tf

const RESERVED = union(PyCall.reserved, ["begin","end"])

typealias TypeOrUnion Union{DataType,Union}
immutable TFArg
  name::Symbol
  default::Any
  doc::AbstractString
  doc1::AbstractString
  mandatory::Bool
  typesym::TypeOrUnion
end

"Avoid name collisions between TensorFlow functions and Julia reserved words/Base functions"
jlname(s::AbstractString) = (s in RESERVED || isdefined(symbol(s))) ? "$(s)_" : s
jlname(s::Symbol) = jlname(string(s))

# At code generation time, we can't yet import the TensorFlow package itself
# because it depends on the generated code.  So to get the types we need, we
# will import TensorFlow.CoreTypes as TFParser.CoreTypes.  Then we sometimes need to strip
# off the module prefix from the type names to use in the generated code.  Hence:
"Strips off module qualifications from type names"
relativename(dt::DataType) = rsplit(string(dt), '.', limit=2)[end]
relativename(u::Union) = "Union{$(join(sort(map(relativename, u.types)), ','))}"

"Format a TensorFlow function argument as a Julia function argument"
function jlformat(a::TFArg)
  s = jlname(a.name)
  if a.typesym != nothing
    s = "$(s)::$(relativename(a.typesym))"
  end
  if !a.mandatory
    if haskey(CoreTypes.DTYPES, a.default)
      "$s=$(CoreTypes.DTYPES[a.default])"
    else
      "$s=$(repr(a.default))"
    end
  else
    s
  end
end

"Pass Julia arguments into Python functios using a **kwargs dict to deal with reserved names"
function pyformat(a::TFArg)
  if a.typesym == Tensor
    # TODO extend/vary?
    ":$(a.name)=>$(jlname(a.name)).x"
  else
    ":$(a.name)=>$(jlname(a.name))"
  end
end

"Try to infer the type of a Tensorflow function argument from name, docstring and default value"
function guesstype(name::Symbol, doc1::AbstractString)
  if name == :dtype
    Dtype
  elseif name == :name
    AbstractString
  elseif name == :shape
    Union{AbstractTensor, DimsType}
  elseif name == :strides
    PyVectorType
  elseif name == :seed
    Int
  elseif name == :keep_dims
    Bool
  elseif name == :placeholder
    Placeholder
  elseif name == :Variable
    Variable
  elseif startswith(string(name), "use_")
    Bool
  elseif startswith(string(name), "num_")
    Int
  elseif endswith(string(name), "_size")
    Int
  else
    guesstype(doc1)
  end
end

function guesstyperet(name::Symbol, doc1::AbstractString)
  if name == :dtype
    Dtype
  elseif name == :shape
    Tensor
  elseif name == :placeholder
    Placeholder
  elseif name == :Variable
    Variable
  elseif name == :Session
    Session
  elseif name == :InteractiveSession
    Session
  elseif endswith(string(name), "Optimizer")
    Optimizer
  else
    guesstype(doc1)
  end
end
guesstyperet(name::Symbol, ::Void) = Any

"Try to infer the type of a Tensorflow function's argument or return value from its docstring"
function guesstype(docstring::AbstractString)
  if ismatch(r"[Tt]ensor", docstring) || ismatch(r"of (output )?type `dtype`", docstring)
    Tensor
  elseif ismatch(r"SparseTensor", docstring)
    SparseTensor
  elseif ismatch(r"d?type", docstring)
    Dtype
  elseif ismatch(r"[Aa]n (Python )?(integer|int|`int`)\b", docstring)
    Int
  elseif ismatch(r"[Aa](n optional) `bool`", docstring)
    Bool
  elseif ismatch(r"(^Boolean|^Whether |`?(True|False)`?)", docstring)
    Bool
  elseif ismatch(r"([Aa]|this) `string`", docstring)
    AbstractString
  elseif ismatch(r"\b([Pp]ath|[Dd]irectory|[Nn]ame)\b`", docstring)
    AbstractString
  elseif ismatch(r"[Aa] `tf.Dtype`", docstring)
    Dtype
  elseif ismatch(r": int,`", docstring)
    Int
  else
    # Give up
    Any
  end
end

"All the data we know about a TensorFlow function"
immutable TFFunction
  o::PyObject
  name::Symbol
  args::Vector{TFArg}
  returndoc::AbstractString
  doc::AbstractString
  typesym::DataType
end

"Return the TFArg of the given name in this TFFunction"
findarg(tff::TFFunction, name::Symbol) = tff.args[findfirst(a->a.name == name, tff.args)]

function showmatch(key::AbstractString, m::RegexMatch, ::AbstractString, leading_spaces::Int)
  join_str = "\n" * " "^(leading_spaces+2)
  return strip(replace(m.captures[1], join_str, ' '))
end

"Look in a TensorFlow docstring to find details about a single argument"
function greparg(docstring::AbstractString, key::AbstractString, num_leading_spaces=4)
  # We match everything non-greedily from the key, indented by num_leading_spaces
  # until we encounter either:
  # - a line with non-whitespace after at most num_leading_spaces
  # - or a line with at most some whitespace

  patt_str = "^ {$(num_leading_spaces)}$(key):(.*?)(?=\\n {0,$(num_leading_spaces)}\\w|^\\s*\\z)"
  patt = Regex(patt_str, "ms")
  m = match(patt, docstring)
  if m == nothing
    error("Couldn't find '$key' in docstring:\n$docstring")
  else
    showmatch(key, m, docstring, num_leading_spaces)
  end
end

"Look in a TensorFlow docstring to find details of the return value"
function grepreturn(docstring::AbstractString)
  m = match(r"  Returns:(.*)"ms, docstring)
  (m == nothing) ? "" : showmatch("Returns", m, docstring, 4)
end
grepreturn(::Void) = nothing

"Construct a TFArg from pieces of TensorFlow docstrings"
TFArg(arg::AbstractString, argdoc::AbstractString, default::Any, mandatory::Bool) = begin
  name = symbol(arg)
  doc1 = split(argdoc, ". ", limit=2)[1]
  typesym = guesstype(name, doc1)
  if typesym == Tensor
    typesym = AbstractTensor
  end
  if default == nothing
    typesym = Union{typesym, Void}
  end
  println("guessed type $typesym for $arg ($argdoc)")
  TFArg(name, default, argdoc, doc1, mandatory, typesym)
end

"Construct a TFFunction from a TensorFlow Python function"
TFFunction(f::PyObject, name::Symbol; skipself=false) = begin
  argspec = PyInspector.PyArgSpec(f)
  doc = PyInspector.pydoc(f)
  num_defaults = length(argspec.defaults)
  num_mandatory = length(argspec.args) - num_defaults
  tfargs = TFArg[]
  for (i, arg) in enumerate(argspec.args)
    if skipself && (symbol(arg) == :self) && (i==1)
      continue
    end
    j = i - num_mandatory
    if j<=0
      mandatory = true
      default = nothing
    else
      mandatory = false
      default = argspec.defaults[j]
    end

    argdoc = ""
    try
      argdoc = greparg(doc, arg)
    catch
      println("Failed to find doc for $(f[:__name__]) $(arg)")
    end

    push!(tfargs, TFArg(arg, argdoc, default, mandatory))
  end
  retdoc = grepreturn(doc)

  if doc == nothing
    doc = ""
    retdoc = ""
  else
    retdoc = grepreturn(doc)
  end

  TFFunction(f, name, tfargs, retdoc, doc, guesstyperet(name, retdoc))

end

TFFunction(f::PyObject) = TFFunction(f, symbol(f[:__name__]))

function tfwritejulia(mod::Module, pymodname::Symbol)
  src = IOBuffer()
  tfwritejulia(mod, pymodname, src)
  seekstart(src)
  readall(src)
end

tfwritejulia(pymod::PyObject, pymodname::Symbol, io::IO; withdoc=true) =
  tfwritejulia(pywrap(pymod), pymodname, io, withdoc=withdoc)

function tfwritejulia(mod::Module, pymodname::Symbol, io::IO; withdoc=true)
  tfwritejuliatypes(mod, pymodname, io; withdoc=withdoc)
    tfwritejuliafuncs(mod, pymodname, io; withdoc=withdoc)
end

"A function that takes a TensorFlow module and tries to output Julia code
  for all its top-level functions"
function tfwritejuliafuncs(mod::Module, pymodname::Symbol, io::IO; withdoc=true)
  for f in PyInspector.pyfuncs(mod)
    # we want to guess types for each of the arguments
    # identify wrapped PyTypes and annotate them, extracting their PyObjects in the call
    tff = TFFunction(f)
    jlarglist = join(map(jlformat, tff.args), ", ")
    pyarglist = join(map(pyformat, tff.args), ", ")
    definition = "$(pymodname).$(tff.name)(;Dict($pyarglist)...)"

    if tff.typesym != Any
      definition = "$(relativename(tff.typesym))($(definition))"
    end

    if withdoc
      write(io, "\n\n\"\"\"\n")
      write(io, replace(tff.doc, "\$", "\\\$"))
      write(io, "\"\"\"")
    end
    write(io, """
\n$(jlname(tff.name))($jlarglist) = $definition
export $(jlname(tff.name))
          """)
  end
end

"A function that takes a TensorFlow module and tries to output Julia code
  for all its top-level types"
function tfwritejuliatypes(mod::Module, pymodname::Symbol, io::IO; withdoc=true)
    for t in filter(t->haskey(t, :__init__),  PyInspector.pytypes(mod))
      name = symbol(t[:__name__])
      init = t[:__init__]
      doc = t[:__doc__]
      if name in names(CoreTypes)
        # Don't declare it again.
        # (Could still import it and add a new constructor?)
        continue
      end
      tff = nothing
      try
        tff = TFFunction(init, name; skipself=true) # TODO fix the type instead of guessing
      catch e
        println("Skipping $name")
        showerror(STDERR, e, catch_backtrace())
        continue
      end

    jlarglist = join(map(jlformat, tff.args), ", ")
    pyarglist = join(map(pyformat, tff.args), ", ")
    definition = "$(pymodname).$(tff.name)(;Dict($pyarglist)...)"

    if tff.typesym != Any
      definition = "$(relativename(tff.typesym))($(definition))"
    end

    if withdoc
      write(io, "\n\n\"\"\"\n")
      write(io, replace(tff.doc, "\$", "\\\$"))
      write(io, "\"\"\"")
    end
    write(io, """
\n$(jlname(tff.name))($jlarglist) = $definition
export $(jlname(tff.name))
          """)
  end
end

end # module


text:
"""
    agc!(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6)

Automatic gain control module for speech signals.
The input signal `wav` would be changed after calling this function.
# Arguments
- `wav     `: raw speech signal
- ` fs     `: sampling rate
- `maxvalue`: specified maximum value for speech signal, maxvalue ∈ (0, +Inf)
- `minstep `: minimum value to change the gain, minstep ∈ (-1, 0)
"""
function agc!(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6)
    minstep  = (-1.0 < minstep < 0.0) ? minstep : -0.6
    maxvalue = ( 0.0 < maxvalue ) ? maxvalue : 0.6
    gain   = 1.0
    winlen = floor(Int, 0.016 * fs)
    frames = div(length(wav)-winlen, winlen) + 1
    firstIds = (0:(frames-1)) .* winlen .+ 1       # 帧起始下标
    lasstIds = firstIds .+ (winlen - 1)            # 帧结束下标
    maxvalue = 1 / maxvalue
    for t = 1:frames
        index = firstIds[t]:lasstIds[t];
        frame = wav[index];
        fmax = maximum(abs.(frame));
        step = 1.0 - fmax*maxvalue;
        step = step * abs(step);
        step = max(minstep, min(step, 0.0));
        gain = 0.8*gain + 0.2*(1 + step);
        wav[index] .*= gain;
    end
end

"""
    agc(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6) -> Array

Automatic gain control module for speech signals.
# Arguments
- `wav     `: raw speech signal
- ` fs     `: sampling rate
- `maxvalue`: specified maximum value for speech signal, maxvalue ∈ (0, +Inf)
- `minstep `: minimum value to change the gain, minstep ∈ (-1, 0)
"""
function agc(wav::Array, fs::Real=16000.0; maxvalue::Real=0.6, minstep::Real=-0.6)
    data = deepcopy(wav)
    agc!(data, fs; maxvalue = maxvalue, minstep=minstep)
    return data
end

max_stars_repo_name: sonosole/AGC

max_issues_repo_name: sonosole/AutoGainCtrl

max_forks_repo_name: sonosole/AutoGainCtrl


text:
# We never need to permute indices on the CPU.
@inline   permute_index(solver_type, ::CPU, i, j, k, Nx, Ny, Nz) = i, j, k
@inline unpermute_index(solver_type, ::CPU, i, j, k, Nx, Ny, Nz) = i, j, k

"""
    _permute_index(i, N)

Perform the permutation [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]
on index `i` out of `N`.
"""
@inline function _permute_index(i, N)
    if (i & 1) == 1  # Same as isodd(i)
        return floor(Int, i/2) + 1
    else
        return N - floor(Int, (i-1)/2)
    end
end

"""
    _permute_index(i, N)

Undo the permutation [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]
on index `i` out of `N`.
"""
@inline function _unpermute_index(i, N)
    if i <= N/2
        return 2i-1
    else
        return 2(N-i+1)
    end
end

@inline   permute_index(::TriplyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, k
@inline unpermute_index(::TriplyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, k

@inline   permute_index(::HorizontallyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j,   _permute_index(k, Nz)
@inline unpermute_index(::HorizontallyPeriodic, ::GPU, i, j, k, Nx, Ny, Nz) = i, j, _unpermute_index(k, Nz)

@inline   permute_index(::Channel, ::GPU, i, j, k, Nx, Ny, Nz) = i,   _permute_index(j, Ny),   _permute_index(k, Nz)
@inline unpermute_index(::Channel, ::GPU, i, j, k, Nx, Ny, Nz) = i, _unpermute_index(j, Ny), _unpermute_index(k, Nz)


text:
"""This module wraps the rust pacosako dynamic library"""

# using Jtac

const DYNLIB_PATH = "../lib/target/debug/libpacosako.so"

mutable struct PacoSako # <: Game # We need to be using Jtac to subtype Game.
    ptr::Ptr{Nothing}
end

"""Initializer and general memory management"""
function PacoSako()::PacoSako
    ptr = ccall((:new, DYNLIB_PATH), Ptr{Nothing}, ())
    wrap_pacosako_ptr(ptr)
end

function wrap_pacosako_ptr(ptr::Ptr{Nothing})::PacoSako
    ps = PacoSako(ptr)
    finalizer(destroy!, ps)
    ps
end

"""Destructor. Never call that manually, or you'll double free / use after free."""
function destroy!(ps::PacoSako)
    # ccall to rust to tell it to release the object.
    ccall((:drop, DYNLIB_PATH), Nothing, (Ptr{Nothing},), ps.ptr)
end

################################################################################
## Implementation of the Game interface. #######################################
################################################################################

function Base.copy(ps::PacoSako)::PacoSako
    ptr = ccall((:clone, DYNLIB_PATH), Ptr{Nothing}, (Ptr{Nothing},), ps.ptr)
    wrap_pacosako_ptr(ptr)
end

# TODO: this one next.
# status(ps::PacoSako) = ccall((:status, DYNLIB_PATH), Int64, (Ptr{Nothing},), ps.ptr)

current_player(ps::PacoSako)::Int64 = ccall((:current_player, DYNLIB_PATH), Int64, (Ptr{Nothing},), ps.ptr)

function legal_actions(ps::PacoSako)
    out = zeros(UInt8, 128)
    ccall((:legal_actions, DYNLIB_PATH), Nothing, (Ptr{Nothing}, Ptr{UInt8}), ps.ptr, out)
    Int.(collect(Iterators.takewhile(x -> x > 0, out)))
end

function apply_action!(ps::PacoSako, action::Int)::PacoSako
    status_code = ccall((:apply_action_bang, DYNLIB_PATH), Int64, (Ptr{Nothing}, UInt8), ps.ptr, UInt8(action))
    ps
end

################################################################################
## Helpers #####################################################################
################################################################################

function println(ps::PacoSako)
    ccall((:print, DYNLIB_PATH), Nothing, (Ptr{Nothing},), ps.ptr)
end

max_stars_repo_name: sebalehnert/pacosako

max_issues_repo_name: sebalehnert/pacosako

max_forks_repo_name: sebalehnert/pacosako


text:
"Parameters for right-hand side function of diagonally implicit Runge-Kutta methods."
mutable struct ParametersDIRK{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    V::Vector{Vector{DT}}

    function ParametersDIRK{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt, zero(TT), zeros(DT,D), create_internal_stage_vector(DT, D, tab.s))
    end
end

"""
Diagonally implicit Runge-Kutta integrator cache.
"""
struct IntegratorCacheDIRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}
    y::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function IntegratorCacheDIRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y)
    end
end

function IntegratorCache{ST}(params::ParametersDIRK{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheDIRK{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersDIRK{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheDIRK{ST,D,S}

@doc raw"""
Diagonally implicit Runge-Kutta integrator solving the system
```math
\begin{aligned}
V_{n,i} &= v (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} .
\end{aligned}
```

"""
struct IntegratorDIRK{DT, TT, D, S, PT <: ParametersDIRK{DT,TT},
                                    ST,# <: NonlinearSolver{DT},
                                    IT <: InitialGuessODE{TT}} <: AbstractIntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorDIRK(params::ParametersDIRK{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorDIRK{DT,D}(equations::NamedTuple, tableau::Tableau{TT}, Δt::TT) where {DT, TT, D}
        # get number of stages
        S = tableau.s

        # check if tableau is diagonally implicit
        @assert !isfullyimplicit(tableau)

        if get_config(:verbosity) ≥ 1
            if isexplicit(tableau)
                @warn "Initializing IntegratorDIRK with explicit tableau $(tableau.name).\nYou might want to use IntegratorERK instead."
            end
        end

        # create params
        params = ParametersDIRK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solvers
        solvers = [create_nonlinear_solver(DT, D, params, caches, i) for i in 1:S]

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v], Δt)

        # create integrator
        IntegratorDIRK(params, solvers, iguess, caches)
    end

    function IntegratorDIRK{DT,D}(v::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorDIRK{DT,D}(NamedTuple{(:v,)}((v,)), tableau, Δt; kwargs...)
    end

    function IntegratorDIRK{DT,D}(v::Function, h::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorDIRK{DT,D}(NamedTuple{(:v,:h)}((v,h)), tableau, Δt; kwargs...)
    end

    function IntegratorDIRK(equation::ODE{DT}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorDIRK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorDIRK{DT,TT,D,S}) where {DT,TT,D,S} = D
@inline has_initial_guess(int::IntegratorDIRK) = true

function initialize!(int::IntegratorDIRK, cache::IntegratorCacheDIRK)
    # initialise initial guess
    cache.t̄ = cache.t - timestep(int)

    int.params.equs[:v](cache.t, cache.q, cache.v)

    initialize!(int.iguess, cache.t, cache.q, cache.v, cache.t̄, cache.q̄, cache.v̄)
end

function update_params!(int::IntegratorDIRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
end

function initial_guess!(int::IntegratorDIRK{DT}, sol::AtomicSolutionODE{DT},
                        cache::IntegratorCacheDIRK{DT}=int.caches[DT]) where {DT}

    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.q̃, cache.ṽ, int.params.tab.c[i])
        for k in eachindex(cache.V[i], cache.ṽ)
            cache.V[i][k] = cache.ṽ[k]
        end
    end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver[i].x[k] = 0
            for j in eachstage(int)
                int.solver[i].x[k] += int.params.tab.a[i,j] * cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{ST}, V::Vector{ST}, Y::Vector{ST},
                         params::ParametersDIRK{DT,TT,D,S}, i::Int) where {ST,DT,TT,D,S}

    local tᵢ::TT

    @assert D == length(Q) == length(V) == length(Y)

    # copy x to Y and compute Q = q + Δt Y
    for k in 1:D
        Y[k] = x[k]
        Q[k] = params.q[k] + params.Δt * Y[k]
    end

    # compute V = v(Q)
    tᵢ = params.t + params.Δt * params.tab.c[i]
    params.equs[:v](tᵢ, Q, V)
end

# Compute stages of fully implicit Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDIRK{DT,TT,D,S},
                          caches::CacheDict, i::Int) where {ST,DT,TT,D,S}

    # temporary variables
    local y1::ST
    local y2::ST

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q[i], cache.V[i], cache.Y[i], params, i)

    # compute b = - (Y-AV)
    for k in 1:D
        y1 = params.tab.a[i,i] * cache.V[i][k]
        y2 = params.tab.â[i,i] * cache.V[i][k]
        for j in 1:i-1
            y1 += params.tab.a[i,j] * params.V[j][k]
            y2 += params.tab.â[i,j] * params.V[j][k]
        end
        b[k] = - cache.Y[i][k] + (y1 + y2)
    end
end

function integrate_step!(int::IntegratorDIRK{DT,TT}, sol::AtomicSolutionODE{DT,TT},
                         cache::IntegratorCacheDIRK{DT}=int.caches[DT]) where {DT,TT}

     # update nonlinear solver parameters from atomic solution
     update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # consecutively solve for all stages
    for i in eachstage(int)
        # call nonlinear solver
        solve!(int.solver[i])

        # print solver status
        print_solver_status(int.solver[i].status, int.solver[i].params)

        # check if solution contains NaNs or error bounds are violated
        check_solver_status(int.solver[i].status, int.solver[i].params)

        # compute vector field at internal stages
        compute_stages!(int.solver[i].x, cache.Q[i], cache.V[i], cache.Y[i], int.params, i)

        # copy velocity field to parameters
        int.params.V[i] .= cache.V[i]
    end

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).b, tableau(int).b̂, timestep(int))

    # update vector field for initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end


text:
@doc raw"""
```
t = find_threshold(Balanced(), histogram, edges)
```
In balanced histogram thresholding, one interprets a  bin as a  physical weight
with a mass equal to its occupancy count. The balanced histogram method involves
iterating the following three steps: (1) choose the midpoint bin index as a
"pivot",  (2) compute the combined weight to the left and right of the pivot bin
and (3) remove the leftmost bin if the left side is the heaviest, and the
rightmost bin otherwise. The algorithm stops when only a single bin remains.
The last bin determines the sought-after threshold.

# Output

Returns a real number `t` in `edges`. The `edges` parameter represents an
`AbstractRange` which specifies the intervals associated with the histogram bins.

# Details
Let ``f_n`` (``n = 1 \ldots N``) denote the number of observations in the ``n``th
bin of the histogram. The balanced histogram method constructs a sequence
of nested intervals

```math
[1,N] \cap \mathbb{Z} \supset I_2 \supset I_3 \supset \ldots \supset I_{N-1},
```
where for ``k = 2 \ldots N-1``
```math
I_k = \begin{cases}
   I_{k-1} \setminus \{\min \left( I_{k-1} \right) \} &\text{if } \sum_{n = \min \left( I_{k-1} \right)}^{I_m}f_n \gt   \sum_{n =  I_m + 1}^{ \max \left( I_{k-1} \right)} f_n, \\
   I_{k-1} \setminus \{\max \left( I_{k-1} \right) \} &\text{otherwise},
\end{cases}
```
and ``I_m = \lfloor \frac{1}{2}\left(  \min \left( I_{k-1} \right) +  \max \left( I_{k-1} \right) \right) \rfloor ``.
The final interval ``I_{N-1}`` consists of a single element which is the bin index
corresponding to the desired threshold.

If one interprets a bin as a physical weight with a mass equal to its occupancy
count, then each step of the algorithm can be conceptualised as removing the
leftmost or rightmost bin to "balance" the resulting histogram on a pivot. The
pivot is defined to be the midpoint between the start and end points of the
interval under consideration.

If it turns out that the single element in ``I_{N-1}`` equals ``1`` or ``N`` then
the original histogram must have a single peak and the algorithm has failed to
find a suitable threshold. In this case the algorithm will fall back to using
the `UnimodalRosin` method to select the threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia

using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Balanced(), counts[1:end], edges)
```

# Reference

1. “BI-LEVEL IMAGE THRESHOLDING - A Fast Method”, Proceedings of the First International Conference on Bio-inspired Systems and Signal Processing, 2008. Available: [10.5220/0001064300700076](https://doi.org/10.5220/0001064300700076)
"""
function find_threshold(algorithm::Balanced, histogram::AbstractArray, edges::AbstractRange)
    # set initial start/middle/end points and weigths
    Iₛ = 1
    Iₑ = length(histogram)
    Iₘ = round(Int, (Iₛ + Iₑ) / 2)
    Wₗ = 0
    Wᵣ = 0
    for i = 1:Iₘ
        Wₗ += histogram[i]
    end
    for i = Iₘ:Iₑ
        Wᵣ += histogram[i]
    end
    while Iₛ < Iₑ
        if Wₗ < Wᵣ
            Wᵣ -= histogram[Iₑ]
            Iₑ -= 1
            if (Iₛ + Iₑ) / 2 < Iₘ
                Wₗ -= histogram[Iₘ]
                Wᵣ += histogram[Iₘ]
                Iₘ -= 1
            end
        else
            Wₗ -= histogram[Iₛ]
            Iₛ += 1
            if (Iₛ + Iₑ) / 2 > Iₘ
                Wₗ += histogram[Iₘ + 1]
                Wᵣ -= histogram[Iₘ + 1]
                Iₘ += 1
            end
        end
    end
    if Iₘ == 1 || Iₘ == length(histogram)
        @warn "Failed to threshold. Falling back to `UnimodalRosin` method."
        return find_threshold(UnimodalRosin(), histogram, edges)
    else
        return edges[Iₘ]
    end
end


text:
using Gtk
using Compose
using Cairo

function cm_in_context(cm)
	n = length(cm)
	r = Compose.rectangle(linspace(0,1-(1/n),n),0.1*ones(n),ones(n)*(1.1/(n)),0.8*ones(n));
	c = Compose.compose(context(),r,fill(cm));
	c
	end

function cm_in_window(cm)
	co = cm_in_context(cm);
	c = Gtk.@GtkCanvas(400,300);
	w = Gtk.@GtkWindow(c,"data win");
	show(c);

	Gtk.draw(c) do widget
    	Compose.draw(CAIROSURFACE(c.back),co)
    	end
	c
	end

function map_from_colormap(cm,z::Array{Float64,2})
     z1 = z .- minimum(z);
     z2 = round(Int,max(ceil(z1  * (length(cm) / maximum(z1))),1));

     display(minimum(z2));
     display(maximum(z2));

     c = similar(z,UInt32);

     for i1 in eachindex(z2)
     	c1 = convert(ColorTypes.RGB24,cm[z2[i1]]);
     	c[i1] = c1.color
     	end
    c
end

function cm_perlin_window(cm)
	z = perlin_noise(120,120);
	p = map_from_colormap(cm,z);
    imgsurf = Cairo.CairoRGBSurface(p)
    imgpatt = Cairo.CairoPattern(imgsurf)
    Cairo.pattern_set_filter(imgpatt, Cairo.FILTER_NEAREST)

	c = Gtk.@GtkCanvas(160,160);
	w = Gtk.@GtkWindow(c,"data win");
	show(c);

	Gtk.draw(c) do widget
		cr = Cairo.CairoContext(c.back)
		r = Gtk.allocation(c)
		iwidth = Float64(r.width)-20.0;
		iheight = Float64(r.height)-20.0;
		Cairo.rectangle(cr,20.0,20.0,iwidth,iheight)

    	imgmatrix = Cairo.CairoMatrix(120.0/iwidth, 0, 0, 120.0/iheight, 0, 0)
    	Cairo.set_matrix(imgpatt, imgmatrix)

	    Cairo.set_source(cr, imgpatt)
    	Cairo.fill(cr)
    	end
    c
    end

max_stars_repo_name: lobingera/ColorMaps


text:
let
    roadway,junction=gen_intersection(roadlength=30.0)
    connect_two_seg!(roadway.segments[2],roadway.segments[3],roadway)
    connect_two_seg!(roadway.segments[4],roadway.segments[5],roadway)
    connect_two_seg!(roadway.segments[6],roadway.segments[7],roadway)
    connect_two_seg!(roadway.segments[8],roadway.segments[1],roadway)
    doc,r = initialize_XML()
    convert_roadway!(r,roadway)
    junctions = [junction]
    handle_junctions(r,junctions,roadway)
    r1 = root(doc)
    r1_id = []
    for child1 in eachelement(r1)
        if haskey(child1, "id")
            push!(r1_id,child1["id"])
        end
        for child2 in eachelement(child1)
            if haskey(child2, "id")
                push!(r1_id,child2["id"])
            end
        end
    end
    doc2 = readxml("test_out.xodr")
    r2 = root(doc2)
    r2_id = []
    for child1 in eachelement(r2)
        if haskey(child1, "id")
            push!(r2_id,child1["id"])
        end
        for child2 in eachelement(child1)
            if haskey(child2, "id")
                push!(r2_id,child2["id"])
            end
        end
    end
    @test r1_id == r2_id
end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
using DescentMethods
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using NLPModels
using LinearOperators, LinearAlgebra
using Krylov

using Stopping
using LineSearch
using DescentMethods

function rosenbrock(x)
	n = 2; m = 2;
	f = []
	push!(f, 10 * (x[2]-x[1]^2))
	push!(f, (x[1]-1))
	return sum(f[i]^2 for i=1:m)
end
x₀ = [-1.2, 1.0]

solvers = [:CG_FR, :CG_HZ, :CG_HS, :CG_PR, :NewtonLDLT, :Newton, :Newlbfgs, :Shamanskii]

for solver in solvers
    nlp = ADNLPModel(rosenbrock, x₀)
    println("Testing $(String(solver))")
    nlpatx = NLPAtX(nlp.meta.x0)
    nlpstop = NLPStopping(nlp, Stopping.unconstrained, nlpatx)

    final_nlp_at_x, optimal = eval(solver)(nlp, nlpstop, verbose = true, linesearch = armijo_ls)
    println("optimal = $(string(optimal))")
    finalize(nlp)
end


text:
include("sampling.jl")
include("id_generator.jl")
include("span_processor.jl")
include("trace_provider.jl")


text:
using Flux.Tracker, NNlib
using Flux.Tracker: TrackedReal, gradcheck
using FluxExtensions

include("utils.jl")
include("layers/resdense.jl")
include("pdfs.jl")
include("triangularloss.jl")

@testset "sumnondiagonal" begin
	x = randn(4,4)
	@test abs(FluxExtensions.sumnondiagonal(x) .- (sum(x) - sum(diag(x)))) < 1e-14
	@test gradcheck(FluxExtensions.sumnondiagonal,x)
end

@testset "scatter" begin
	x = randn(4,4)
	@test gradcheck((x) -> sum(FluxExtensions.scatter(x,4)),x)
	@test gradcheck((x) -> sum(FluxExtensions.scatter(x,[1,2,3,4])),x)
end

@testset "conversion to Float32" begin
	m = Dense(3, 1)
	x = randn(Float32, 3,5)
	mm = to32(m)
	@test eltype(Flux.data(m(x))) == Float64
	@test eltype(Flux.data(mm(x))) == Float32
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
# Ordinary metrics

###########################################################
#
#   Metric types
#
###########################################################

type Euclidean <: Metric end
type SqEuclidean <: SemiMetric end
type Chebyshev <: Metric end
type Cityblock <: Metric end
type Jaccard <: Metric end
type RogersTanimoto <: Metric end

immutable Minkowski{T <: Real} <: Metric
    p::T
end

type Hamming <: Metric end

type CosineDist <: SemiMetric end
type CorrDist <: SemiMetric end

type ChiSqDist <: SemiMetric end
type KLDivergence <: PreMetric end

immutable RenyiDivergence{T <: Real} <: PreMetric
    p::T # order of power mean (order of divergence - 1)
    is_normal::Bool
    is_zero::Bool
    is_one::Bool
    is_inf::Bool
    function RenyiDivergence(q)
        # There are four different cases:
        #   simpler to separate them out now, not over and over in eval_op()
        is_zero = q ≈ zero(T)
        is_one = q ≈ one(T)
        is_inf = isinf(q)

        # Only positive Rényi divergences are defined
        !is_zero && q < zero(T) && throw(ArgumentError("Order of Rényi divergence not legal, $(q) < 0."))

        new(q - 1, !(is_zero || is_one || is_inf), is_zero, is_one, is_inf)
    end
end
RenyiDivergence{T}(q::T) = RenyiDivergence{T}(q)

type JSDivergence <: SemiMetric end

type SpanNormDist <: SemiMetric end

typealias UnionMetrics Union{Euclidean, SqEuclidean, Chebyshev, Cityblock, Minkowski, Hamming, Jaccard, RogersTanimoto, CosineDist, CorrDist, ChiSqDist, KLDivergence, RenyiDivergence, JSDivergence, SpanNormDist}

###########################################################
#
#  Define Evaluate
#
###########################################################

function evaluate(d::UnionMetrics, a::AbstractArray, b::AbstractArray)
    if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    s = eval_start(d, a, b)
    if size(a) == size(b)
        @simd for I in eachindex(a, b)
            @inbounds ai = a[I]
            @inbounds bi = b[I]
            s = eval_reduce(d, s, eval_op(d, ai, bi))
        end
    else
        for (Ia, Ib) in zip(eachindex(a), eachindex(b))
            @inbounds ai = a[Ia]
            @inbounds bi = b[Ib]
            s = eval_reduce(d, s, eval_op(d, ai, bi))
        end
    end
    return eval_end(d, s)
end
result_type{T1, T2}(dist::UnionMetrics, ::AbstractArray{T1}, ::AbstractArray{T2}) =
    typeof(eval_end(dist, eval_op(dist, one(T1), one(T2))))
eval_start(d::UnionMetrics, a::AbstractArray, b::AbstractArray) =
    zero(result_type(d, a, b))
eval_end(d::UnionMetrics, s) = s

evaluate{T <: Number}(dist::UnionMetrics, a::T, b::T) = eval_end(dist, eval_op(dist, a, b))

# SqEuclidean
@inline eval_op(::SqEuclidean, ai, bi) = abs2(ai - bi)
@inline eval_reduce(::SqEuclidean, s1, s2) = s1 + s2

sqeuclidean(a::AbstractArray, b::AbstractArray) = evaluate(SqEuclidean(), a, b)
sqeuclidean{T <: Number}(a::T, b::T) = evaluate(SqEuclidean(), a, b)

# Euclidean
@inline eval_op(::Euclidean, ai, bi) = abs2(ai - bi)
@inline eval_reduce(::Euclidean, s1, s2) = s1 + s2
eval_end(::Euclidean, s) = sqrt(s)
euclidean(a::AbstractArray, b::AbstractArray) = evaluate(Euclidean(), a, b)
euclidean(a::Number, b::Number) = evaluate(Euclidean(), a, b)

# Cityblock
@inline eval_op(::Cityblock, ai, bi) = abs(ai - bi)
@inline eval_reduce(::Cityblock, s1, s2) = s1 + s2
cityblock(a::AbstractArray, b::AbstractArray) = evaluate(Cityblock(), a, b)
cityblock{T <: Number}(a::T, b::T) = evaluate(Cityblock(), a, b)

# Chebyshev
@inline eval_op(::Chebyshev, ai, bi) = abs(ai - bi)
@inline eval_reduce(::Chebyshev, s1, s2) = max(s1, s2)
# if only NaN, will output NaN
@inline eval_start(::Chebyshev, a::AbstractArray, b::AbstractArray) = abs(a[1] - b[1])
chebyshev(a::AbstractArray, b::AbstractArray) = evaluate(Chebyshev(), a, b)
chebyshev{T <: Number}(a::T, b::T) = evaluate(Chebyshev(), a, b)

# Minkowski
@inline eval_op(dist::Minkowski, ai, bi) = abs(ai - bi) .^ dist.p
@inline eval_reduce(::Minkowski, s1, s2) = s1 + s2
eval_end(dist::Minkowski, s) = s .^ (1/dist.p)
minkowski(a::AbstractArray, b::AbstractArray, p::Real) = evaluate(Minkowski(p), a, b)
minkowski{T <: Number}(a::T, b::T, p::Real) = evaluate(Minkowski(p), a, b)

# Hamming
@inline eval_op(::Hamming, ai, bi) = ai != bi ? 1 : 0
@inline eval_reduce(::Hamming, s1, s2) = s1 + s2
hamming(a::AbstractArray, b::AbstractArray) = evaluate(Hamming(), a, b)
hamming{T <: Number}(a::T, b::T) = evaluate(Hamming(), a, b)

# Cosine dist
function eval_start{T<:AbstractFloat}(::CosineDist, a::AbstractArray{T}, b::AbstractArray{T})
    zero(T), zero(T), zero(T)
end
@inline eval_op(::CosineDist, ai, bi) = ai * bi, ai * ai, bi * bi
@inline function eval_reduce(::CosineDist, s1, s2)
    a1, b1, c1 = s1
    a2, b2, c2 = s2
    return a1 + a2, b1 + b2, c1 + c2
end
function eval_end(::CosineDist, s)
    ab, a2, b2 = s
    max(1 - ab / (sqrt(a2) * sqrt(b2)), zero(eltype(ab)))
end
cosine_dist(a::AbstractArray, b::AbstractArray) = evaluate(CosineDist(), a, b)

# Correlation Dist
_centralize(x::AbstractArray) = x .- mean(x)
evaluate(::CorrDist, a::AbstractArray, b::AbstractArray) = cosine_dist(_centralize(a), _centralize(b))
corr_dist(a::AbstractArray, b::AbstractArray) = evaluate(CorrDist(), a, b)
result_type(::CorrDist, a::AbstractArray, b::AbstractArray) = result_type(CosineDist(), a, b)

# ChiSqDist
@inline eval_op(::ChiSqDist, ai, bi) = abs2(ai - bi) / (ai + bi)
@inline eval_reduce(::ChiSqDist, s1, s2) = s1 + s2
chisq_dist(a::AbstractArray, b::AbstractArray) = evaluate(ChiSqDist(), a, b)

# KLDivergence
@inline eval_op(::KLDivergence, ai, bi) = ai > 0 ? ai * log(ai / bi) : zero(ai)
@inline eval_reduce(::KLDivergence, s1, s2) = s1 + s2
kl_divergence(a::AbstractArray, b::AbstractArray) = evaluate(KLDivergence(), a, b)

# RenyiDivergence
function eval_start{T<:AbstractFloat}(::RenyiDivergence, a::AbstractArray{T}, b::AbstractArray{T})
    zero(T), zero(T)
end

@inline function eval_op{T<:AbstractFloat}(dist::RenyiDivergence, ai::T, bi::T)
    if ai == zero(T)
        return zero(T), zero(T)
    elseif dist.is_normal
        return ai, ai .* ((ai ./ bi) .^ dist.p)
    elseif dist.is_zero
        return ai, bi
    elseif dist.is_one
        return ai, ai * log(ai / bi)
    else # otherwise q = ∞
        return ai, ai / bi
    end
end

@inline function eval_reduce{T<:AbstractFloat}(dist::RenyiDivergence,
                                               s1::Tuple{T, T},
                                               s2::Tuple{T, T})
    if dist.is_inf
        if s1[1] == zero(T)
            return s2
        elseif s2[1] == zero(T)
            return s1
        else
            return s1[2] > s2[2] ? s1 : s2
        end
    else
        return s1[1] + s2[1], s1[2] + s2[2]
    end
end

function eval_end(dist::RenyiDivergence, s)
    if dist.is_zero || dist.is_normal
        log(s[2] / s[1]) / dist.p
    elseif dist.is_one
        return s[2] / s[1]
    else # q = ∞
        log(s[2])
    end
end

renyi_divergence(a::AbstractArray, b::AbstractArray, q::Real) = evaluate(RenyiDivergence(q), a, b)

# JSDivergence
@inline function eval_op{T}(::JSDivergence, ai::T, bi::T)
    u = (ai + bi) / 2
    ta = ai > 0 ? ai * log(ai) / 2 : zero(log(one(T)))
    tb = bi > 0 ? bi * log(bi) / 2 : zero(log(one(T)))
    tu = u > 0 ? u * log(u) : zero(log(one(T)))
    ta + tb - tu
end
@inline eval_reduce(::JSDivergence, s1, s2) = s1 + s2
js_divergence(a::AbstractArray, b::AbstractArray) = evaluate(JSDivergence(), a, b)

# SpanNormDist
function eval_start(::SpanNormDist, a::AbstractArray, b::AbstractArray)
    a[1] - b[1], a[1] - b[1]
end
@inline eval_op(::SpanNormDist, ai, bi)  = ai - bi
@inline function eval_reduce(::SpanNormDist, s1, s2)
    min_d, max_d = s1
    if s2 > max_d
        max_d = s2
    elseif s2 < min_d
        min_d = s2
    end
    return min_d, max_d
end

eval_end(::SpanNormDist, s) = s[2] - s[1]
spannorm_dist(a::AbstractArray, b::AbstractArray) = evaluate(SpanNormDist(), a, b)
function result_type{T1, T2}(dist::SpanNormDist, ::AbstractArray{T1}, ::AbstractArray{T2})
    typeof(eval_op(dist, one(T1), one(T2)))
end

# Jaccard

@inline eval_start(::Jaccard, a::AbstractArray{Bool}, b::AbstractArray{Bool}) = 0, 0
@inline eval_start{T}(::Jaccard, a::AbstractArray{T}, b::AbstractArray{T}) = zero(T), zero(T)
@inline function eval_op(::Jaccard, s1, s2)
    abs_m = abs(s1 - s2)
    abs_p = abs(s1 + s2)
    abs_p - abs_m, abs_p + abs_m
end
@inline function eval_reduce(::Jaccard, s1, s2)
    @inbounds a = s1[1] + s2[1]
    @inbounds b = s1[2] + s2[2]
    a, b
end
@inline function eval_end(::Jaccard, a)
    @inbounds v = 1 - (a[1]/a[2])
    return v
end
jaccard(a::AbstractArray, b::AbstractArray) = evaluate(Jaccard(), a, b)

# Tanimoto

@inline eval_start(::RogersTanimoto, a::AbstractArray, b::AbstractArray) = 0, 0, 0, 0
@inline function eval_op(::RogersTanimoto, s1, s2)
  tt = s1 && s2
  tf = s1 && !s2
  ft = !s1 && s2
  ff = !s1 && !s2
  tt, tf, ft, ff
end
@inline function eval_reduce(::RogersTanimoto, s1, s2)
    @inbounds begin
        a = s1[1] + s2[1]
        b = s1[2] + s2[2]
        c = s1[3] + s2[3]
        d = s1[4] + s1[4]
    end
    a, b, c, d
end
@inline function eval_end(::RogersTanimoto, a)
    @inbounds numerator = 2(a[2] + a[3])
    @inbounds denominator = a[1] + a[4] + 2(a[2] + a[3])
    numerator / denominator
end
rogerstanimoto{T <: Bool}(a::AbstractArray{T}, b::AbstractArray{T}) = evaluate(RogersTanimoto(), a, b)

###########################################################
#
#  Special method
#
###########################################################

# SqEuclidean
function pairwise!(r::AbstractMatrix, dist::SqEuclidean, a::AbstractMatrix, b::AbstractMatrix)
    At_mul_B!(r, a, b)
    sa2 = sumabs2(a, 1)
    sb2 = sumabs2(b, 1)
    pdist!(r, sa2, sb2)
end
function pdist!(r, sa2, sb2)
    for j = 1 : size(r,2)
        sb = sb2[j]
        @simd for i = 1 : size(r,1)
            @inbounds r[i,j] = sa2[i] + sb - 2 * r[i,j]
        end
    end
    r
end
function pairwise!(r::AbstractMatrix, dist::SqEuclidean, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    sa2 = sumsq_percol(a)
    @inbounds for j = 1 : n
        for i = 1 : j-1
            r[i,j] = r[j,i]
        end
        r[j,j] = 0
        for i = j+1 : n
            r[i,j] = sa2[i] + sa2[j] - 2 * r[i,j]
        end
    end
    r
end

# Euclidean
function pairwise!(r::AbstractMatrix, dist::Euclidean, a::AbstractMatrix, b::AbstractMatrix)
    m, na, nb = get_pairwise_dims(r, a, b)
    At_mul_B!(r, a, b)
    sa2 = sumsq_percol(a)
    sb2 = sumsq_percol(b)
    @inbounds for j = 1 : nb
        for i = 1 : na
            v = sa2[i] + sb2[j] - 2 * r[i,j]
            r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
    end
    r
end

function pairwise!(r::AbstractMatrix, dist::Euclidean, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    sa2 = sumsq_percol(a)
    @inbounds for j = 1 : n
        for i = 1 : j-1
            r[i,j] = r[j,i]
        end
        @inbounds r[j,j] = 0
        for i = j+1 : n
            v = sa2[i] + sa2[j] - 2 * r[i,j]
            r[i,j] = isnan(v) ? NaN : sqrt(max(v, 0.))
        end
    end
    r
end

# CosineDist

function pairwise!(r::AbstractMatrix, dist::CosineDist, a::AbstractMatrix, b::AbstractMatrix)
    m, na, nb = get_pairwise_dims(r, a, b)
    At_mul_B!(r, a, b)
    ra = sqrt!(sumsq_percol(a))
    rb = sqrt!(sumsq_percol(b))
    for j = 1 : nb
        @simd for i = 1 : na
            @inbounds r[i,j] = max(1 - r[i,j] / (ra[i] * rb[j]), 0)
        end
    end
    r
end
function pairwise!(r::AbstractMatrix, dist::CosineDist, a::AbstractMatrix)
    m, n = get_pairwise_dims(r, a)
    At_mul_B!(r, a, a)
    ra = sqrt!(sumsq_percol(a))
    @inbounds for j = 1 : n
        @simd for i = j+1 : n
            r[i,j] = max(1 - r[i,j] / (ra[i] * ra[j]), 0)
        end
        r[j,j] = 0
        for i = 1 : j-1
            r[i,j] = r[j,i]
        end
    end
    r
end

# CorrDist
_centralize_colwise(x::AbstractVector) = x .- mean(x)
_centralize_colwise(x::AbstractMatrix) = x .- mean(x, 1)
function colwise!(r::AbstractVector, dist::CorrDist, a::AbstractMatrix, b::AbstractMatrix)
    colwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function colwise!(r::AbstractVector, dist::CorrDist, a::AbstractVector, b::AbstractMatrix)
    colwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function pairwise!(r::AbstractMatrix, dist::CorrDist, a::AbstractMatrix, b::AbstractMatrix)
    pairwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
end
function pairwise!(r::AbstractMatrix, dist::CorrDist, a::AbstractMatrix)
    pairwise!(r, CosineDist(), _centralize_colwise(a))
end


text:
#=
   From the Netica model Book Bags.neta
   """
   Book Bags                  Copyright 1998 Norsys Software Corp.

   There are two book bags each containing 10 poker chips.  In one
   bag there are 7 red and 3 blue.  In the other bag there are 3
   red and 7 blue. Five chips are drawn out of one of the bags and
   shown to the subject (one at a time then returned to the bag).
   The subject does not know which bag the chips came from.
   There is an equal chance that the draws are made from either
   bag. After each draw the subject reports which bag he believes
   the chips are coming from and provides a probability that the
   chips are being drawn from that bag.

   The problem comes from the early "revision of judgment" work
   that indicated that people were conservative with respect to Bayes.
   """

   Bag1 (7 blue, 3 red): 1 Bag2 (3 blue, 7 red): 2

   We've seen 1 blue chip(s).
   Distributions of variable bag
   bag1       =>   70262  (0.702620)
   bag2       =>   29738  (0.297380)

   We've seen 2 blue chip(s).
   Distributions of variable bag
   bag1       =>   84680  (0.846800)
   bag2       =>   15320  (0.153200)

   We've seen 3 blue chip(s).
   Distributions of variable bag
   bag1       =>   93152  (0.931520)
   bag2       =>    6848  (0.068480)

   We've seen 4 blue chip(s).
   Distributions of variable bag
   bag1       =>   90562  (0.905620)
   bag2       =>    9438  (0.094380)

   We've seen 5 blue chip(s).
   Distributions of variable bag
   bag1       =>   98388  (0.983880)
   bag2       =>    1612  (0.016120)

   Cf ~/blog/book_bags.blog
      ~/webppl/book_bags.wppl

=#

using Turing
# using Plots, StatsPlots
include("jl_utils.jl")

@model function book_bags(num_blue_chips=2)
    bag1 = 1 # 7 blue 3 red
    bag2 = 2 # 3 blue 7 red

    bag ~ Categorical([0.5,0.5])

    blue = 1
    red = 2

    n = 5
    draw = tzeros(num_blue_chips)
    for d in 1:num_blue_chips
        if bag == bag1
            draw[d] ~ Categorical([0.7,0.3])
        else
            draw[d] ~ Categorical([0.3,0.7])
        end
    end

    # We've seen $test number of blue chips
    for d in 1:num_blue_chips
        true ~ Dirac(draw[d] == blue)
    end

end

function run_book_bags(num_blue_chips)
    println("\nWe've seen $num_blue_chips blue chip(s).")
    model = book_bags(num_blue_chips)
    num_chains = 4
    # chs = sample(model, Prior(), 1000)
    chs = sample(model, MH(), 100_000)
    # chs = sample(model, PG(15), 10_000)
    # chs = sample(model, IS(), 10_000)
    # chs = sample(model, SMC(), 10_000)
    # chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chains)

    # chs = sample(model, SGLD(), 10_000)
    # chs = sample(model,NUTS(), 10_000)
    # chs = sample(model,HMC(0.01,5), 10_000)

    # display(chs)
    # display(plot(chs))

    show_var_dist_pct(chs,:bag,["bag1","bag2"])

end

println("Bag1 (7 blue, 3 red): 1 Bag2 (3 blue, 7 red): 2")
for test in 1:5
    run_book_bags(test)
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank