text:
module Continuum
 using RvSpectMLBase
 using Statistics, NaNMath # For mean, median, etc
 using DSP                 # For smoothing flux
 using DataInterpolations  # For predicting continuum between anchors
 using RollingFunctions    # For rolling maximum.  Used in setting radius of rolling-pin.
 using SortFilters         # For rolling quantile for fitting polynomial while avoiding lines + computing IQR for sigma clipping.  https://github.com/sairus7/SortFilters.jl
 using FITSIO              # For reading SED
 using Polynomials         # For fitting out low-order polynomial

# constants
speed_of_light_mks = 3e8  # m/s  TODO: Update number
fwhm_sol = 7.9e3          # m/s  TODO: Update number

function read_master_sed_neid(;filename::String = "neidMaster_HR_SmoothLampSED_20210101.fits", path::String = "." )
  f = FITS(joinpath(path,filename))
  sed = read(f[2])
  close(f)
  return sed
end

function normalize_by_sed(λ::AA1, flux::AA2, var::AA3, sed::AA4; poly_order::Integer = 4, orders_to_use::AV5 = 1:size(flux,2),
      half_width::Integer = 50, quantile::Real = 0.9, verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Integer, AA1<:AbstractArray{T1,2}, AA2<:AbstractArray{T2,2}, AA3<:AbstractArray{T3,2}, AA4<:AbstractArray{T4,2}, AV5<:AbstractVector{T5}  }
  @assert size(λ) == size(flux) == size(var) == size(sed)
  f_norm = fill(NaN,size(flux,1),size(flux,2))
  var_norm = fill(NaN,size(var,1),size(var,2))
  min_usable_pixels_in_order = 7
  for order_idx in orders_to_use
      pix_mask = .!(isnan.(λ[:,order_idx]) .| isnan.(flux[:,order_idx]) .| (flux[:,order_idx].<0) .| isnan.(var[:,order_idx])  .| (var[:,order_idx].<=0.0) .| (sed[:,order_idx] .<= 0 ))
      if sum(pix_mask) < min_usable_pixels_in_order continue   end
      lambda = view(λ,pix_mask,order_idx)
      f_obs = view(flux,pix_mask,order_idx)
      sed_view = view(sed,pix_mask,order_idx)
      var_obs = view(var,pix_mask,order_idx)
      if poly_order >= 1 && half_width >=1
        upper_envelope = calc_rolling_quantile(lambda,f_obs./sed_view, half_width=half_width, quantile=quantile)
        try
           coeff = Polynomials.fit(lambda,upper_envelope,poly_order,weights=sed_view.^2 ./var_obs)
           polyfit = coeff.(view(λ,:,order_idx))
        catch
           println("# Warning:  Polynomial fit to continuum for order index ", order_idx, " failed (order = ", poly_order, ".")
           #println("#           size(lambda) = ", size(lambda), " size(upper_envelope) = ", size(upper_envelope), " size(weights) = ", size(sed_view.^2 ./var_obs))
           #println("#           ", sum(isnan.(lambda)), ",  ", sum(isnan.(upper_envelope)), ",  ", sum(isnan.(sed_view)), ",  ", sum(isnan.(var_obs)), ",  ", sum(isinf.(sed_view.^2 ./var_obs)))
           polyfit = 1
        end
      else
        if verbose  println("# Skipping polynomial fit to flux/sed for order index ", order_idx)  end
        polyfit = 1
      end
      f_norm[:,order_idx]   .= view(flux,:,order_idx)./(view(sed,:,order_idx).*polyfit)
      var_norm[:,order_idx] .=  view(var,:,order_idx)./(view(sed,:,order_idx).*polyfit).^2
      idx_invalid_norm = view(sed,:,order_idx).*polyfit .<= 0.0
      f_norm[idx_invalid_norm,order_idx] .= NaN
      var_norm[idx_invalid_norm,order_idx] .= NaN
    end
  return (flux=f_norm, var=var_norm)
end

function calc_mean_snr( flux::AV1, var::AV2 ) where { T1<:Real, T2<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2} }
   idx_bad = isnan.(flux) .| isnan.(var) .| (var .<=0.0)
   return mean(flux[.!idx_bad]./sqrt.(var[.!idx_bad]))
end

function smooth(f::AV2; half_width::Integer = 6 ) where {  T2<:Real, AV2<:AbstractVector{T2} }
  #println("# size(f) to smooth: ", size(f))
  w=DSP.hanning(1+2*half_width)
  if 1+2*half_width<3
    return f
  elseif length(f)<1+2*half_width
    return f
  else
    return y_smooth = conv(f, w/sum(w))[1+half_width:end-half_width]
  end
end

function calc_rolling_max(f::AV2; width::Integer = 13 ) where {  T2<:Real, AV2<:AbstractVector{T2} }
 #shift = max(1,floor(Int64,width//2))
 @assert width%2 == 1
 shift = max(1,convert(Int64,(width+1)//2))
 #println(" size(f) = ",size(f), " width= ", width, " shift = ", shift, " rolling_max = ", size(rolling(NaNMath.maximum,f,width)), " target = ", size((shift+1):(length(f)-shift)))
 z = similar(f)
 #z[(shift):(length(f)-shift+1)] .= rollmax(f,width)
 z[(shift):(end-shift+1)] .= rolling(NaNMath.maximum,f,width)
 z[(end-shift+2):end] .= z[end-shift+1]
 z[1:shift-1] .= z[shift]
 return z
end

function find_local_maxima(f::AV2; half_width::Integer = 6 ) where {  T2<:Real, AV2<:AbstractVector{T2} }
  rollingmax = calc_rolling_max(f,width=half_width*2+1)
  findall(f .== rollingmax)
end

function calc_rolling_median(λ::AV1, f::AV2; half_width::Integer = 4, verbose::Bool = false ) where { T1<:Real, T2<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2} }
  @assert length(λ) == length(f)
  @assert half_width >= 1
  #shift = half_width # floor(Int64,width//2)
  f_filtered = similar(f)
  #if any(isnan.(f))
  #   f_filtered .= running(NaNMath.median,f,width)
  #else
     f_filtered[(1+half_width):(length(f)-half_width)] .= movsort(f, 2*half_width+1, 0.5)[1+2*half_width:end]
     f_filtered[1:half_width] .= f_filtered[1+half_width]
     f_filtered[(length(f)-half_width+1):end] .= f_filtered[length(f)-half_width]
  #end
  return f_filtered
end

function calc_rolling_quantile(λ::AV1, f::AV2; half_width::Integer = 4, quantile::Real, verbose::Bool = false ) where { T1<:Real, T2<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2} }
  @assert length(λ) == length(f)
  @assert half_width >= 1
  output = similar(λ)
  moving_quartiles = movsort(f, 2*half_width+1, quantile)
  output[(1+half_width):(end-half_width)] .= moving_quartiles[1+2*half_width:end]
  output[1:half_width] .= output[1+half_width]
  output[end-half_width+1:end] .= output[end-half_width]
  return output
end

#=
function calc_rolling_median_and_max_to_clip(λ::AV1, f::AV2; width::Integer = 8, cr_pad_width::Integer = 1,
          num_clippings::Integer = 3, verbose::Bool = false ) where { T1<:Real, T2<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2} }
  @assert length(λ) == length(f)
  @assert width >= 4
  @assert width%2 == 0
  @assert 0<= cr_pad_width <= 7
  @assert 1 <= num_clippings <= 10
  shift = floor(Int64,width//2)
  high_quantile = (length(λ) > 200) ? 0.99 : 0.95
  p = (0.25, 0.5, 0.75, high_quantile)
  function iqr(mq::AVNT, i::Integer) where { T<:Real, NT<:NTuple{4,T}, AVNT<:AbstractVector{NT} }
      (mq[i][3]-mq[i][1])    # standard iqr
  end
  function iqr_upper(mq::AVNT, i::Integer) where { T<:Real, NT<:NTuple{4,T}, AVNT<:AbstractVector{NT} }
      2*(mq[i][3]-mq[i][2])  # double Q3-Q2
  end
  q75(mq::AVNT, i::Integer) where { T<:Real, NT<:NTuple{4,T}, AVNT<:AbstractVector{NT} } = mq[i][3]
  qhi(mq::AVNT, i::Integer) where { T<:Real, NT<:NTuple{4,T}, AVNT<:AbstractVector{NT} } = mq[i][4]

  output = similar(λ)
  mask = trues(length(λ))
  f_median_filtered = similar(f)
  f_clean = deepcopy(f)

  if verbose flush(stdout)  end
  for i in 1:num_clippings
     moving_quartiles = movsort(f_clean, width, p)
     f_median_filtered[(1+shift):(length(f)-shift)] .= map(i->moving_quartiles[i][2], (1+2*shift):length(moving_quartiles) )
     f_median_filtered[1:shift] .= f_median_filtered[1+shift]
     f_median_filtered[(length(f)-shift+1):end] .= f_median_filtered[length(f)-shift]
     max_iqr = map(i->q75(moving_quartiles,i)+1.5*iqr_upper(moving_quartiles,i),1:length(λ))
     idx = (1+2*shift):length(moving_quartiles)
     output[(1+shift):(length(f)-shift)] .= max_iqr[idx]
     #=
     # Rassine uses a second pass with simple high quantile with a wider window (to avoid being affected by lines).
     # For now we use the same window
     width_pass2 = 100  # TODO: Automate for orders.  FOr now arbtrary rather than 100A, while testing on small chunks
     shift_pass2 = floor(Int64,width_pass2//2)
     max_99 = movsort(f, width_pass2, 0.99)
     output[(1+shift):(length(f)-shift)] .= max.(output[idx]),max_99[idx])
     =#
     #output[(1+shift):(length(f)-shift)] .= max.(output[idx],map(i->qhi(moving_quartiles,i),idx))
     output[(1+shift):(length(f)-shift)] .= map(i->max(output[i],qhi(moving_quartiles,i)),idx)

     # Fill in edges
     output[1:shift] .= output[1+shift]
     output[(length(f)-shift+1):end] .= output[length(f)-shift]

     # Identify & replace high outliers, plus neighboring pixels
     for i in 1:length(mask)
       if f[i] <= output[i] continue end
       for j in max(1,i-cr_pad_width):min(i+cr_pad_width,length(mask))
          mask[j] = false
       end
     end
     #mask .&= (f.<=output)
     if verbose println("# pass ", i, ", removing ", sum((f.>output)), " hi pixels for a total of ", sum(.!mask), ".") end
     f_clean[.!mask] .= f_median_filtered[.!mask]
  end
  return (f_filtered=f_median_filtered, f_threshold=output, f_clean=f_clean)
end
=#

function longest_cluster_same_sign(y::AV) where { T<:Real, AV<:AbstractVector{T} }
  longest = 0
  current = 1
  this_start = 1
  start = 1
  stop = 1
  sign_last = sign(first(y))
  for i in 2:length(y)
    sign_this = sign(y[i])
    if sign_this == sign_last
      current += 1
    else
      if current > longest
        longest = current
        start = this_start
        stop = i-1
      end
      current = 1
      this_start = i
    end
    sign_last = sign_this
  end
  if current > longest
    longest = current
    start = this_start
    stop = length(y)
  end
  return (len=longest, start=start, stop=stop)
end

function calc_rollingpin_radius(λ::AV1, f::AV2; fwhm::Real = fwhm_sol,
            S1_width_factor::Real = 40, S2_width_factor::Real = 10,
            min_R_factor::Real = 100, ν::Real = 1, verbose::Bool = false ) where { T1<:Real, AV1<:AbstractVector{T1}, T2<:Real, AV2<:AbstractVector{T2}  }
   λlo = minimum(λ)
   S1_width = S1_width_factor*fwhm/speed_of_light_mks*λlo
   S1_width = min(2*floor(Int64,S1_width/2),2)+1
   S2_width = S2_width_factor*S1_width
   S2_width = 2*floor(Int64,S2_width/2)+1
   f_max1 = calc_rolling_max(f,width=S1_width)
   f_max2 = calc_rolling_max(f,width=S2_width)
   (longest_len, longest_start, longest_stop) = longest_cluster_same_sign(f_max2.-f_max1)
   longest_Δλoverλ = (λ[longest_stop]-λ[longest_start]) / ((λ[longest_start]+λ[longest_stop])/2)
   #penalty = (f_max2.-f_max1)./f_max2
   #conversion_fwhm_sig = (10*min_λ/(sqrt(8*log(2))*speed_of_light_kmps))
   #min_radius = fwhm * conversion_fwhm_sig # par_R
   #radius = min_radius .* λ./min_λ
   rollingpin_R_min = min_R_factor*λlo*(fwhm/speed_of_light_mks)/log(8*log(2))
   #rollingpin_R_max = max_R_factor*rollingpin_R_min  # Arbitrary, not same as Rassine's auto setting, see https://github.com/MichaelCretignier/Rassine_public/blob/master/Rassine.py#L718
   rollingpin_R_max = max(λlo*longest_Δλoverλ,rollingpin_R_min*2)
   #rollingpin_R_min = 2*floor(Int64,rollingpin_R_min/2)
   #rollingpin_R_max = 2*floor(Int64,rollingpin_R_max/2)
   if verbose
     println("longest_len = ", longest_len)
     println("longest_Δλoverλ = ", longest_Δλoverλ)
     println("rollingpin_R_min = ", rollingpin_R_min)
     println("rollingpin_R_max = ", rollingpin_R_max)
   end
   rollingpin_radius = (λ./λlo).*(rollingpin_R_min .+ (rollingpin_R_max-rollingpin_R_min).*((f_max2.-f_max1)./f_max2).^ν)
   #plot!(lambda,f_max1,color=:red)
   #plot!(lambda,f_max2,color=:magenta)
   #plot!(lambda,penalty)
   return rollingpin_radius
end

function calc_continuum_anchors(λ::AV1, f::AV2; radius::AV3,
            stretch_factor::Real = 1.0, fwhm::Real = 7.3 #= km/s=#,
                        verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3} }
  @assert length(λ) == length(f)
  #speed_of_light_kmps = 3e5
  #par_stretching = 2# 0.5# 0.5 #11.385
  extrema_λ = extrema(λ)
  min_λ = first(extrema_λ)
  extrema_f = extrema(f)
  #normalization = (extrema_f[2]-extrema_f[1])/(extrema_λ[2]-extrema_λ[1])
  quantile_for_extrema = max(0.5,(length(f) >= 4000) ? 0.995 : 1-20/length(f))
  (min_f, max_f) = quantile(Iterators.filter(!isnan, f),[1-quantile_for_extrema,quantile_for_extrema])
  min_f = 0.0
  normalization = (max_f-min_f)/(extrema_λ[2]-extrema_λ[1])
  normalization *= stretch_factor
  #= v1
  waves = λ.-λ'
  mask = sign.(waves).>0
  #distance = sign.(waves).*sqrt.(waves.^2 .+((f.-f')./normalization).^2)
  distance[mask] = sqrt.(waves[mask].^2 .+((f.-f')[mask]./normalization).^2)
  =#
  #= v2
  distance = zeros(length(λ),length(λ))
  mask = sign.(λ.-λ').>0
  distance[mask] .= sqrt.((λ.-λ')[mask].^2 .+((f.-f')[mask]./normalization).^2)
  =#
  #= v3
  distance = zeros(length(λ),length(λ))
  for i in 1:length(λ)
    for j in 1:length(λ)
      if λ[i]<=λ[j] continue end
      distance[i,j] = sqrt((λ[i]-λ[j])^2 +((f[i]-f[j])/normalization)^2)
    end
  end
  #return distance
  =#
  # v4
  function calc_dist(i::Integer, j::Integer)  # TODO OPT: eliminate sqrt
    if λ[i]<=λ[j]
      return 0.0
    else
      sqrt((λ[i]-λ[j])^2 +((f[i]-f[j])/normalization)^2)
    end
  end
  numero = range(1,stop=length(λ))
  mask = falses(length(λ))
  dist = zeros(length(λ))
  #conversion_fwhm_sig = (10*min_λ/(sqrt(8*log(2))*speed_of_light_kmps))
  #min_radius = fwhm * conversion_fwhm_sig # par_R
  #println("# old min_radius = ", min_radius)
  if verbose println("# passed radius = ", extrema(radius)) end
  #radius = min_radius .* λ./min_λ
  keep = Vector{Int64}()
  sizehint!(keep,max(32,min(256,2^floor(Int64,log2(length(λ))-2))))
  push!(keep,1)
  j = 1
  while length(λ)-j>2
    par_R = radius[j]
    if j==1 par_R /= 1.5 end
    map!(i->calc_dist(i,j),dist,1:length(λ))
    #mask .= (0.0 .< distance[:,j].<2.0*par_R)
    mask .= (0.0 .< dist.<2.0*par_R)
    while sum(mask)==0
      par_R *= 1.5
      #mask .= (0.0 .< distance[:,j].<2.0*par_R)
      mask .= (0.0 .< dist.<2.0*par_R)
    end
    p1 = [ λ[j] f[j]/normalization ]
    p2 = hcat(λ[mask], f[mask]./normalization)
    delta = p2 .- p1
    #c = sqrt.(delta[:,1].^2 .+delta[:,2].^2)
    c = sqrt.(vec(sum(delta.^2,dims=2)))
    harg = (par_R.^2 .-0.25.*c.^2)
    #=
    if any(harg.<0)
      println("# j = ", j, " dist = ", distance[j:j+10,j])
      println("# R^2 = ", par_R^2)
      println("# c^2/4 = ", c.^2/4)
      println("# harg = ", harg)
    end
    h = zeros(length(harg))
    hmask = harg.>0
    h[hmask] .= sqrt.(harg[hmask])
    =#
    h = sqrt.(harg)
    cx = p1[1] .+ 0.5.*delta[:,1] .-h./c.*delta[:,2]
    cy = p1[2] .+ 0.5.*delta[:,2] .+h./c.*delta[:,1]
    #return (cx,cy)
    #cond1 = (cy.-p1[2]).>=0
    #theta = zeros(length(cy))
    mintheta = Inf
    argmintheta = 0
    for i in 1:length(cy)
      if cy[i]-p1[2] >0
        acos_arg = (cx[i]-p1[1])/par_R
        if acos_arg > 1.0 acos_arg = 1.0
        elseif  acos_arg < -1.0 acos_arg = -1.0  end
        thetai = -acos(acos_arg)+π
      else
        asin_arg = (cy[i]-p1[2])/par_R
        if asin_arg > 1.0 asin_arg = 1.0
        elseif  asin_arg < -1.0 asin_arg = -1.0  end
        thetai = -asin(asin_arg)+π
      end
      if thetai < mintheta
        mintheta = thetai
        argmintheta = i
      end
    end
    #theta[cond1] .= -acos.((cx[cond1].-p1[1])./par_R).+π
    #theta[.!cond1] .= -asin.((cy[.!cond1].-p1[2])./par_R).+π
    #j2 = argmin(theta)
    j2 = argmintheta
    j = numero[mask][j2]
    push!(keep,j)
  end
  if verbose println("# using ", length(keep), " anchor points") end
  return keep
end

function calc_continuum_from_anchors( λ::AV1, f::AV2, anchors::AV3; λout::AV4 = λ,
              verbose::Bool = false ) where {
                T1<:Real, T2<:Real, T3<:Integer, T4<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}  }
  calc_continuum_from_anchors_hybrid( λ, f, anchors, λout=λout,verbose=verbose)
end

function calc_continuum_from_anchors_linear( λ::AV1, f::AV2, anchors::AV3; λout::AV4 = λ,
              verbose::Bool = false ) where {
                T1<:Real, T2<:Real, T3<:Integer, T4<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}  }
  @assert length(anchors) >= 2
  #if length(anchors) <= 1   return ones(size(λout))  end
  interp_linear = LinearInterpolation(f[anchors],λ[anchors])
  function extrap(x::Real)
    if x<first(interp_linear.t)
      slope = (interp_linear.u[2]-interp_linear.u[1])/(interp_linear.t[2]-interp_linear.t[1])
      return first(interp_linear.u)+(x-interp_linear.t[1])*slope
    elseif x>last(interp_linear.t)
      slope = (interp_linear.u[end]-interp_linear.u[end-1])/(interp_linear.t[end]-interp_linear.t[end-1])
      return last(interp_linear.u)+(x-interp_linear.t[end])*slope
    else
      return interp_linear(x)
    end
  end
  return extrap.(λout)
end

function calc_continuum_from_anchors_cubic( λ::AV1, f::AV2, anchors::AV3; λout::AV4 = λ,
              verbose::Bool = false ) where {
                T1<:Real, T2<:Real, T3<:Integer, T4<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}  }
  interp_cubic = CubicSpline(f[anchors],λ[anchors])
  interp_cubic.(λout)
end

function calc_continuum_from_anchors_hybrid( λ::AV1, f::AV2, anchors::AV3; λout::AV4 = λ,
              verbose::Bool = false ) where {
                T1<:Real, T2<:Real, T3<:Integer, T4<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}  }
  output = similar(λout)
  if length(anchors) < 2
    output .= ones(size(λout))
  elseif length(anchors) < 9
    output .= calc_continuum_from_anchors_linear( λ, f, anchors, λout=λout,verbose=verbose)
  else
    λmin = λ[anchors[4]]
    λmax = λ[anchors[end-3]]
    idx_pix_cubic = searchsortedfirst(λout,λmin):searchsortedlast(λout,λmax)
    output[idx_pix_cubic] .= calc_continuum_from_anchors_cubic( λ, f, anchors, λout=λout[idx_pix_cubic],verbose=verbose)
    idx_pix_linear1 = 1:first(idx_pix_cubic)
    idx_pix_linear2 = last(idx_pix_cubic):length(λout)
    idx_pix_linear = vcat(idx_pix_linear1,idx_pix_linear2)
    if verbose println("# idx_pix_linear = ", idx_pix_linear, " \n # λ_linear = ", λout[idx_pix_linear]) end
    output[idx_pix_linear] .= calc_continuum_from_anchors_linear( λ, f, anchors, λout=λout[idx_pix_linear],verbose=verbose)
  end
  return output
end

function replace_edge_anchor_vals!(f::AV1, n::Integer = 3 ) where { T1<:Real, AV1<:AbstractVector{T1} }
  @assert length(f) >= 2*n+1
  f[1:(n+1)] .= f[n+1]
  f[end-n:end] .= f[length(f)-n-1]
  return f
end

function find_clean_anchors_by_slope(anchors::AV1, f::AV2; threshold::Real = 0.995, verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3} }
  nanchors = length(anchors)
  @assert nanchors == length(f)
  @assert nanchors >= 8
  Δabsslope = zeros(length(anchors))
  for i in 2:nanchors-1
    slope_hi = (f[i+1]-f[i])/(anchors[i+1]-anchors[i])
    slope_lo = (f[i]-f[i-1])/(anchors[i]-anchors[i-1])
    Δabsslope[i] = abs(slope_hi-slope_lo)
  end
  threshold = quantile(Δabsslope[2:end-1],threshold)
  mask = Δabsslope.<threshold
  return mask
end

# Currently, not used
function filter_anchor_groups(anchors::AV1; min_num::Integer = 3, verbose::Bool = false ) where { T1<:Real, AV1<:AbstractVector{T1} }
  #Δpix = anchors[2:end] .- anchors[1:end-1]
  mask = trues(length(anchors))
  return mask
end

function merge_nearby_anchors(anchors::AV1, λ::AV2; threshold::Real = 1, verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3} }
  Δλ = λ[anchors[2:end]] .- λ[anchors[1:end-1]]
  close_anchor_pair_mask = Δλ .<= threshold
  nclose = sum(close_anchor_pair_mask)
  if verbose println("# Found ", nclose, " close anchor pairs out of ", length(anchors), ".")  end
  if nclose == 0
    return anchors
  else
    if verbose
      println("# Anchors = ", anchors)
      println("Close anchor pairs: ", findall(close_anchor_pair_mask), " => ", anchors[findall(close_anchor_pair_mask)])
    end
  end
  old_anchor_mask = falses(length(anchors))
  if first(Δλ) > threshold
    old_anchor_mask[1] = true
  end
  for i in 2:length(Δλ)
    if (Δλ[i-1] > threshold) && (Δλ[i] > threshold)
      old_anchor_mask[i] = true
    end
  end
  if last(Δλ) > threshold
    old_anchor_mask[length(anchors)] = true
  end
  old_anchors = anchors[old_anchor_mask]
  if verbose   println("# Old anchors = ", old_anchors)  end
  close_anchor_idx = (1:(length(anchors)-1))[close_anchor_pair_mask]
  new_anchors = zeros(Int64,nclose)
  resize!(new_anchors,0)
  idx_start = 1
  for i in 2:nclose
    if close_anchor_idx[i]-close_anchor_idx[i-1] == 1
      continue
    else
      idx_stop = i-1

      if idx_start == idx_stop
        merged_anchor_idx = idx_start
      elseif (idx_stop-idx_start)%2 == 0
        merged_anchor_idx = idx_start + convert(Int64,(idx_stop-idx_start)//2)
      else
        merged_anchor_idx = idx_start + convert(Int64,(idx_stop-idx_start+1)//2)
      end
      #=
      num_in_group = idx_stop-idx_start+1
      min_equidist = Inf
      merged_anchor_idx = 0
      for j in 0:(num_in_group-1)
        λ_anchor_prev = λ[anchors[max(close_anchor_idx[idx_start]-1,1)]]
        λ_anchor_next = λ[anchors[min(close_anchor_idx[idx_stop]+1,length(anchors))]]
        λthis = λ[anchors[close_anchor_idx[idx_start+j]]]
        equidist = abs((λthis-λ_anchor_prev)-(λ_anchor_next-λthis))
        #equidist = abs(log(λthis^2/(λ_anchor_prev*λ_anchor_next)))
        if equidist<min_equidist
          min_equidist = equidist
          merged_anchor_idx = idx_start+j
        end
      end
      =#
      merged_anchor = anchors[close_anchor_idx[merged_anchor_idx]]
      push!(new_anchors,merged_anchor)
      idx_start = i
    end
  end
  if idx_start == nclose == 1
    push!(new_anchors,anchors[close_anchor_idx[idx_start]])
  elseif idx_start != nclose
    idx_stop = nclose
    if idx_start == idx_stop
      merged_anchor_idx = idx_start
    elseif (idx_stop-idx_start)%2 == 0
      merged_anchor_idx = idx_start + convert(Int64,(idx_stop-idx_start)//2)
    else
      merged_anchor_idx = idx_start + convert(Int64,(idx_stop-idx_start+1)//2)
    end
    #=
    num_in_group = idx_stop-idx_start+1
    min_equidist = Inf
    merged_anchor_idx = 0
    for j in 0:(num_in_group-1)
      λ_anchor_prev = λ[anchors[max(close_anchor_idx[idx_start]-1,1)]]
      λ_anchor_next = λ[anchors[min(close_anchor_idx[idx_stop]+1,length(anchors))]]
      λthis = λ[anchors[close_anchor_idx[idx_start+j]]]
      equidist = abs(log(λthis^2/(λ_anchor_prev*λ_anchor_next)))
      if equidist<min_equidist
        min_equidist = equidist
        merged_anchor_idx = idx_start+j
      end
    end
    =#
    merged_anchor = anchors[close_anchor_idx[merged_anchor_idx]]
    push!(new_anchors,merged_anchor)
  end
  if verbose println("# Keeping ", length(new_anchors), " new anchors, ", length(old_anchors), " old anchors.")  end
  anchors_out = mergesorted(old_anchors,new_anchors)
  return anchors_out
end

function calc_continuum(λ::AV1, f_obs::AV2, var_obs::AV3; λout::AV4 = λ, fwhm::Real = fwhm_sol, ν::Real =1.0,
  smoothing_half_width::Integer = 6, local_maximum_half_width::Integer = smoothing_half_width+1,
   stretch_factor::Real = 5.0, merging_threshold::Real = 0.25,  min_R_factor::Real = 100.0, verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}  }
 #clip_width_A = 1.0
 #clip_width_pix = clip_width_A/(λ[2]-λ[1])
 #@assert 1 <= clip_width_pix < Inf
 #clip_width_pix = 2*floor(Int64,clip_width_pix/2)
 mean_snr_per_pix = calc_mean_snr(f_obs,var_obs)
 #smoothing_half_width = (mean_snr_per_pix >= 30) ? smoothing_half_width : 40
 if mean_snr_per_pix < 20
    smoothing_half_width = min(100, ceil(Int64,6*(20/mean_snr_per_pix)^2))
 end
 f_smooth = Continuum.smooth(f_obs, half_width=smoothing_half_width)
 idx_local_maxima = find_local_maxima(f_smooth, half_width=local_maximum_half_width)
 if verbose println("# Found ", length(idx_local_maxima), " local maxima." )  end
 if length(idx_local_maxima) < 7
   println("# Warning only ",  length(idx_local_maxima), " local maxima, aborting order.")
   half_width = min_R_factor*(fwhm/speed_of_light_mks)/log(8*log(2))*minimum(λ)/(λ[2]-λ[1])
   f_alt_continuum = Continuum.smooth(f_obs, half_width=floor(Int64,half_width/2)*2 )
   return (anchors=Int64[], continuum=f_alt_continuum, f_filtered=f_smooth)
 end
 #(f_median_filtered, f_clip_threshold, f_clean) = calc_rolling_median_and_max_to_clip(λ,f_obs, width=clip_width_pix, verbose=verbose)
 rollingpin_radius = calc_rollingpin_radius(λ, f_smooth, fwhm=fwhm, min_R_factor=min_R_factor, verbose=verbose, ν=ν)

 anch_orig = calc_continuum_anchors(λ[idx_local_maxima],f_smooth[idx_local_maxima],radius=rollingpin_radius[idx_local_maxima], stretch_factor=stretch_factor, verbose=verbose)
 if verbose println("# Found ", length(anch_orig), " potential anchors." )  end
 anch_orig = idx_local_maxima[anch_orig]
 anchor_vals = f_smooth[anch_orig]
 if length(anch_orig) >= 8  # 2n+1, for replace_edge_anchor_vals!
   replace_edge_anchor_vals!(anchor_vals)
   anch_mask = find_clean_anchors_by_slope(anch_orig,anchor_vals, threshold = 0.95, verbose=verbose)
   if verbose println("# After rejected high-slope anchors ", sum(anch_mask), " anchors left." )  end
 else
   anch_mask = trues(length(anch_orig))
 end
 #anchors_merged = anch_orig[anch_mask]
 if merging_threshold > 0  && length(anch_orig) >= 4
     anchors_merged = merge_nearby_anchors(anch_orig[anch_mask],λ, threshold=merging_threshold, verbose=verbose)
     if verbose println("# After merging ", length(anchors_merged), " anchors left." )  end
 else
    anchors_merged = anch_orig[anch_mask]
 end
 if length(anchors_merged) < 4
   println("# Warning only ",  length(anchors_merged), " anchors found, aborting order.")
   half_width = min_R_factor*(fwhm/speed_of_light_mks)/log(8*log(2))*minimum(λ)/(λ[2]-λ[1])
   f_alt_continuum = Continuum.smooth(f_obs, half_width=floor(Int64,half_width/2)*2 )
   return (anchors=Int64[], continuum=f_alt_continuum, f_filtered=f_smooth)
 end
  #anchors_merged = anch_orig[anch_mask]
 continuum = calc_continuum_from_anchors_hybrid(λ,f_smooth,anchors_merged, λout=λout) # , verbose=verbose)
 return (anchors=anchors_merged, continuum=continuum, f_filtered=f_smooth)
end

function calc_continuum(λ::AV1, f_obs::AV2, var_obs::AV3, anchors::AV5; λout::AV4 = λ, smoothing_half_width::Integer = 6, verbose::Bool = false ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Integer, AV1<:AbstractVector{T1}, AV2<:AbstractVector{T2}, AV3<:AbstractVector{T3}, AV4<:AbstractVector{T4}, AV5<:AbstractVector{T5}  }
 #clip_width_A = 1.0
 #clip_width_pix = clip_width_A/(λ[2]-λ[1])
 #@assert 1 <= clip_width_pix < Inf
 #clip_width_pix = 2*floor(Int64,clip_width_pix/2)
 mean_snr_per_pix = calc_mean_snr(f_obs,var_obs)
 #smoothing_half_width = (mean_snr_per_pix >= 30) ? smoothing_half_width : 40
 if mean_snr_per_pix < 30
    smoothing_half_width = min(100, ceil(Int64,6*(30/mean_snr_per_pix)^2))
 end
 f_smooth = Continuum.smooth(f_obs, half_width=smoothing_half_width)
 continuum = calc_continuum_from_anchors_hybrid(λ,f_smooth,anchors, λout=λout) # , verbose=verbose)
 return (anchors=anchors, continuum=continuum, f_filtered=f_smooth)
end

abstract type AbstractGetPixelRangeFunctor end

struct GetPixelRangeFromWavelengthGrid{AA} <: AbstractGetPixelRangeFunctor where { T<:Real, AA<:AbstractArray{T,2} }
      data::AA
end

function (gpr::GetPixelRangeFromWavelengthGrid{AA1})(ord_idx::Integer) where { T1<:Real, AA1<:AbstractArray{T1,2} }
    return 1:size(gpr.data,1)
 end

#=
function default_get_pixel_range(λ::AA1, ord_idx::Integer) where { T1<:Real, AA1<:AbstractArray{T1,2} }
  return 1:size(λ,1)
end
=#

struct GetPixelRangeFromInstrument{InstT<:AbstractInstrument} <:  AbstractGetPixelRangeFunctor
      inst::InstT
end
function (gpr::GetPixelRangeFromInstrument{InstT})(ord_idx::Integer) where {  InstT<:AbstractInstrument }
   #inst = get_inst(spec)
   pix = min_pixel_in_order(gpr.inst):max_pixel_in_order(gpr.inst)
end

function calc_continuum(λ::AA1, f_obs::AA2, var_obs::AA3; λout::AA4 = λ, fwhm::Real = fwhm_sol, ν::Real =1.0,
  stretch_factor::Real = 5.0, merging_threshold::Real = 0.25, smoothing_half_width::Integer = 6, min_R_factor::Real = 100.0,
        orders_to_use::AbstractVector{<:Integer} = 1:size(λ,2), get_pixel_range::GPRT = GetPixelRangeFromWavelengthGrid(λ), verbose::Bool = false ) where {
        T1<:Real, T2<:Real, T3<:Real, T4<:Real, AA1<:AbstractArray{T1,2}, AA2<:AbstractArray{T2,2}, AA3<:AbstractArray{T3,2} , AA4<:AbstractArray{T4,2}, GPRT<:AbstractGetPixelRangeFunctor }
  @assert size(λ) == size(f_obs) == size(var_obs)
  @assert size(λout,2) == size(λout,2)
  num_orders = size(λout,2)
  anchors_2d = fill(Int64[],num_orders)
  continuum_2d = fill(NaN,size(λout,1),size(λout,2))
  f_filtered_2d = fill(NaN,size(λout,1),size(λout,2))
  #Threads.@threads for ord_idx in orders_to_use
  for ord_idx in orders_to_use
    if verbose
      println("# Order index = ", ord_idx)
      flush(stdout)
    end
    pix = get_pixel_range(ord_idx)
    λ_use = view(λ,pix,ord_idx)
    f_obs_use = convert.(Float64,view(f_obs,pix,ord_idx))
    var_obs_use = convert.(Float64,view(var_obs,pix,ord_idx))
    if all(isnan.(λ_use)) || all(isnan.(f_obs_use)) || all(isnan.(var_obs_use))   continue   end
    λout_use = view(λout,pix,ord_idx)
    (anchors_1d, continuum_order_1d, f_filtered_1d) = Continuum.calc_continuum(λ_use,f_obs_use,var_obs_use, λout=λout_use,
         stretch_factor=stretch_factor, merging_threshold=merging_threshold, smoothing_half_width=smoothing_half_width, min_R_factor=min_R_factor, verbose=verbose)
    anchors_2d[ord_idx] = anchors_1d
    continuum_2d[pix,ord_idx] .= continuum_order_1d
    f_filtered_2d[pix,ord_idx] .= f_filtered_1d
  end
  return (anchors=anchors_2d, continuum=continuum_2d, f_filtered=f_filtered_2d)
end

function calc_continuum(λ::AA1, f_obs::AA2, var_obs::AA3, anchors::AVV; λout::AA4 = λ,
   smoothing_half_width::Integer = 6,
   #fwhm::Real = fwhm_sol, ν::Real =1.0,
  #stretch_factor::Real = 5.0, merging_threshold::Real = 0.25, min_R_factor::Real = 100.0,
        orders_to_use::AbstractVector{<:Integer} = 1:size(λ,2), get_pixel_range::GPRT = GetPixelRangeFromWavelengthGrid(λ), verbose::Bool = false ) where {
        T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Integer, AA1<:AbstractArray{T1,2}, AA2<:AbstractArray{T2,2}, AA3<:AbstractArray{T3,2} , AA4<:AbstractArray{T4,2}, AVV<:AbstractVector{Vector{T5}}, GPRT<:AbstractGetPixelRangeFunctor }
  @assert size(λ) == size(f_obs) == size(var_obs)
  @assert size(λout,2) == size(λout,2)
  num_orders = size(λout,2)
  @assert length(anchors) == num_orders
  #anchors_2d = fill(Int64[],num_orders)
  continuum_2d = fill(NaN,size(λout,1),size(λout,2))
  f_filtered_2d = fill(NaN,size(λout,1),size(λout,2))
  Threads.@threads for ord_idx in orders_to_use
  #for ord_idx in orders_to_use
    if verbose
      println("# Order index = ", ord_idx)
      flush(stdout)
    end
    pix = get_pixel_range(ord_idx)
    λ_use = view(λ,pix,ord_idx)
    λout_use = view(λout,pix,ord_idx)
    f_obs_use = convert.(Float64,view(f_obs,pix,ord_idx))
    var_obs_use = convert.(Float64,view(var_obs,pix,ord_idx))
    if all(isnan.(λ_use)) || all(isnan.(f_obs_use)) || all(isnan.(var_obs_use))   continue   end

    (anchors_1d, continuum_order_1d, f_filtered_1d) = Continuum.calc_continuum(λ_use,f_obs_use,var_obs_use,
          anchors[ord_idx], λout=λout_use, verbose=verbose)
    #anchors_2d[ord_idx] = anchors_1d
    continuum_2d[pix,ord_idx] .= continuum_order_1d
    f_filtered_2d[pix,ord_idx] .= f_filtered_1d
  end
  return (anchors=anchors, continuum=continuum_2d, f_filtered=f_filtered_2d)
end

# Code for mergesorted from
# https://discourse.julialang.org/t/looking-for-merge-algorithm/47473/2
import Base.Order: Ordering, Forward, ord, lt

indices(x,i)=eachindex(x)

# merge sorted vectors vl and vr into v
# from indices lo to hi in v
function mergesorted!(v::AbstractVector,
                      lo::Int, hi::Int,
                      vl::AbstractVector,
                      lol::Int, hil::Int,
                      vr::AbstractVector,
                      lor::Int, hir::Int,
                      order::Ordering)
    c = lol
    p = lor
    nl = hil
    nr = hir
    i = lo
    @inbounds while c <= nl && p <= nr && i <= hi
        if lt(order, vr[p], vl[c])
            v[i] = vr[p]
            p = p+1
            i = i+1
        else
            v[i] = vl[c]
            c = c+1
            i = i+1
        end
    end
    @inbounds while p <= nr && i <= hi
        v[i] = vr[p]
        i = i+1
        p = p+1
    end
    @inbounds while c <= nl && i <= hi
        v[i] = vl[c]
        i = i+1
        c = c+1
    end
    v
end

function mergesorted!(v::AbstractVector, vl::AbstractVector,
                      vr::AbstractVector, order::Ordering)
    inds = indices(v,1)
    indsl = indices(vl,1)
    indsr = indices(vr,1)
    mergesorted!(v,first(inds),last(inds),vl,first(indsl),last(indsl),
                 vr,first(indsr),last(indsr),order)
end

"""
    mergesorted!(v, vl, vr; lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)
Merge sorted vectors `vl` and `vr`, overwriting vector `v`. Assumes
that `vl` and `vr` are sorted and does not check whether `vl` or `vr`
are sorted. You could used `issorted` to check if `vl` and `vr` are sorted.
If length of `v` is less than the sum of the lengths of
`vl` and `vr`, this simply stops when all indices in `v` are filled.
The `by` keyword lets you provide a function that will be applied to
each element before comparison; the `lt` keyword allows providing a
custom "less than" function; use `rev=true` to reverse the sorting
order. These options are independent and can be used together in all
possible combinations: if both `by` and `lt` are specified, the `lt`
function is applied to the result of the `by` function; `rev=true`
reverses whatever ordering specified via the `by` and `lt` keywords.
"""
function mergesorted!(v::AbstractVector,
                      vl::AbstractVector,
                      vr::AbstractVector;
                      lt=isless,
                      by=identity,
                      rev::Bool=false,
                      order::Ordering=Forward)
    ordr = ord(lt,by,rev,order)
    mergesorted!(v, vl, vr, ordr)
end

"""
    mergesorted(vl, vr; lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)
Merge sorted vectors `vl` and `vr`. See [`mergesorted!`](@ref).
"""
function mergesorted(vl::AbstractVector,
                      vr::AbstractVector;
                      lt=isless,
                      by=identity,
                      rev::Bool=false,
                      order::Ordering=Forward)
    v = similar(promote_type(typeof(vl),typeof(vr)), length(vl)+length(vr))
    ordr = ord(lt,by,rev,order)
    mergesorted!(v, vl, vr, ordr)
end

end # module Continuum


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -13.266290625856113 <= q <= 11.094199160852137)

                     add_NL_constraint(m, :(softplus(-0.2977454197977769 + 0.7722469146312885*softplus(0.9100824262656433 + 0.0704650714228956*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.5815923034237733*softplus(-0.15598180224250147 + 0.7627892034632455*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + -0.731486602090937*softplus(-0.9577813343303001 + 0.660774759467007*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + softplus(0.23088458854123672 + 0.9025271253252471*softplus(0.9100824262656433 + 0.0704650714228956*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.9612479731086991*softplus(-0.15598180224250147 + 0.7627892034632455*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.9633061422344626*softplus(-0.9577813343303001 + 0.660774759467007*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + softplus(-0.8299198997441475 + 0.7709117851978275*softplus(0.9100824262656433 + 0.0704650714228956*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.5338511667423695*softplus(-0.15598180224250147 + 0.7627892034632455*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.3695524647159538*softplus(-0.9577813343303001 + 0.660774759467007*softplus(-0.7110372864748489 + -0.3189000958940382*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*softplus(0.7299768289984225 + 0.7045346620386517*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*softplus(-0.5099367345236407 + -0.3810722083004521*softplus(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*softplus(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*softplus(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function Base.intersect(l::LineSegment2D{T}, poly::QuadraticPolygon{N, 2, T}
                       ) where {N, T <:Union{Float32, Float64}}
    # Create the quadratic segments that make up the polygon and intersect each one
    points = zeros(MVector{N, Point2D{T}})
    npoints = 0x0000
    M = N ÷ 2
    for i ∈ 1:M-1
        hits, ipoints = l ∩ QuadraticSegment2D(poly[i], poly[i + 1], poly[i + M])
        for j in 1:hits
            npoints += 0x0001
            points[npoints] = ipoints[j]
        end
    end
    hits, ipoints = l ∩ QuadraticSegment2D(poly[M], poly[1], poly[N])
    for j in 1:hits
        npoints += 0x0001
        points[npoints] = ipoints[j]
    end
    return npoints, SVector(points)
end

# Cannot mutate BigFloats in an MVector, so we use a regular Vector
function Base.intersect(l::LineSegment2D{BigFloat}, poly::QuadraticPolygon{N, 2, BigFloat}
                       ) where {N}
    # Create the quadratic segments that make up the polygon and intersect each one
    points = zeros(Point2D{BigFloat}, N)
    npoints = 0x0000
    M = N ÷ 2
    for i ∈ 1:M-1
        hits, ipoints = l ∩ QuadraticSegment2D(poly[i], poly[i + 1], poly[i + M])
        for j in 1:hits
            npoints += 0x0001
            points[npoints] = ipoints[j]
        end
    end
    hits, ipoints = l ∩ QuadraticSegment2D(poly[M], poly[1], poly[N])
    for j in 1:hits
        npoints += 0x0001
        points[npoints] = ipoints[j]
    end
    return npoints, SVector{N, Point2D{BigFloat}}(points)
end

# A quadratic triangle, defined in 3D.

# Summary of intersection methods:
#       Triangulation
#           - Speed varies dramatically with precision
#               - The faster method when precision at or below 2 decimal places is desired
#           - Can provide false positives and false negatives due to approximation of the surface
#           - Constant time for a given number of triangles
#               - The only way to determine if two intersections exist is to test until two unique
#                 points are found. Since this happens infrequently, all points are tested, so time
#                 isn't wasted on control logic.
#       Newton-Raphson (iterative)
#           - Speed varies slightly based upon line segment length and orientation
#               - A shorter line segment will converge faster
#               - Converges based upon Jacobian matrix.
#                   - If derivatives are small, the iteration can become slow
#               - The faster method when precision beyond 2 decimal places is desired
#           - Less accurate
#               - May falsely give one intersection instead of two, especially for longer segments.
#               - This is due to the point of convergence being dependent on the initial guess
#                 point. The two starting points are placed close to the line segment start/stop
#                 to try to mitigate this.
#           - Precision to 6+ decimal places.
#       Overall
#           - Triangulation is predictable in speed, slow for high precision, and largely accurate
#           - Newton-Raphson is unpredictable in speed, fast for high precision, and
#               can be inaccurate for 2 intersections
#           - Consider the difference in timing between intersections for porting to GPU
#               - Newton-Raphson may cause thread divergence

# Triangulate, then intersect
function intersect(l::LineSegment_3D{T}, tri6::Triangle6_3D{T};
        N::Int64 = 25) where {T <: AbstractFloat}
    triangles = triangulate(tri6, N)
    npoints = 0
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    intersections = l .∩ triangles
    bools = map(x->x[1], intersections)
    points = map(x->x[2], intersections)
    npoints = count(bools)
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    if npoints === 0
        return false, 0, p₁, p₂
    elseif npoints === 1
        p₁ = points[argmax(bools)]
        return true, 1, p₁, p₂
    elseif npoints === 2
        indices = findall(bools)
        p₁ = points[indices[1]]
        p₂ = points[indices[2]]
        # Check uniqueness
        if p₁ ≈ p₂
            return true, 1, p₁, p₂
        else
            return true, 2, p₁, p₂
        end
    else
        # Account for 3 points and 4 points?
        # If intersection is on edge shared by two triangles on entrance and/or exit 3/4 intersections
        # can be detected
        return true, -1, p₁, p₂
    end
end

# A more exact intersection algorithm that triangulation, uses Newton-Raphson.
function intersect_iterative(l::LineSegment_3D{T}, tri6::Triangle6_3D{T};
        N::Int64=30) where {T <: AbstractFloat}
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    npoints = 0
    u⃗ = l.points[2] - l.points[1]
    ray_start = real_to_parametric(l(0), tri6; N=10) # closest r,s to the ray start
    ray_stop  = real_to_parametric(l(1), tri6; N=10) # closest r,s to the ray stop
    # The parametric coordinates corresponding to the start of the line segment
    r₁ = ray_start[1]
    s₁ = ray_start[2]
    t₁ = T(0)
    # The parametric coordinates corresponding to the start of the line segment
    r₂ = ray_stop[1]
    s₂ = ray_stop[2]
    t₂ = T(1)
    # Iteration for point 1
    err₁ = tri6(r₁, s₁) - l(t₁)
    for i = 1:N
        ∂r₁, ∂s₁ = derivatives(tri6, r₁, s₁)
        J₁ = hcat(∂r₁.x, ∂s₁.x, -u⃗.x)
        # If matrix is singular, it's probably bad luck. Just perturb it a bit.
        if abs(det(J₁)) < 1e-5
            r₁, s₁, t₁ = [r₁, s₁, t₁] + rand(3)/10
        else
            r₁, s₁, t₁ = [r₁, s₁, t₁] - inv(J₁) * err₁.x
            errₙ = tri6(r₁, s₁) - l(t₁)
            if norm(errₙ - err₁) < 5e-6
                break
            end
            err₁ = errₙ
        end
    end
    # Iteration for point 2
    err₂ = tri6(r₂, s₂) - l(t₂)
    for j = 1:N
        ∂r₂, ∂s₂ = derivatives(tri6, r₂, s₂)
        J₂ = hcat(∂r₂.x, ∂s₂.x, -u⃗.x)
        if abs(det(J₂)) < 1e-5
            r₂, s₂, t₂ = [r₂, s₂, t₂] + rand(3)/10
        else
            r₂, s₂, t₂ = [r₂, s₂, t₂] - inv(J₂) * err₂.x
            errₘ = tri6(r₂, s₂) - l(t₂)
            if norm(errₘ - err₂) < 5e-6
                break
            end
            err₂ = errₘ
        end
    end

    p₁ = l(t₁)
    if (0 ≤ r₁ ≤ 1) && (0 ≤ s₁ ≤ 1) && (0 ≤ t₁ ≤ 1) && (p₁ ≈ tri6(r₁, s₁))
        npoints += 1
    end

    p₂ = l(t₂)
    if (0 ≤ r₂ ≤ 1) && (0 ≤ s₂ ≤ 1) && (0 ≤ t₂ ≤ 1) && (p₂ ≈ tri6(r₂, s₂))
        npoints += 1
        # If only point 2 is valid, return it as p₁
        # If points are duplicate, reduce npoints
        if npoints === 2 && p₁ ≈ p₂
            npoints -= 1
        elseif npoints === 1
            p₁ = p₂
        end
    end
    return npoints > 0, npoints, p₁, p₂
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Distributed, SharedArrays, Random
include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))

addprocs_with_testenv(4)
@test nprocs() == 5

@everywhere using Test, SharedArrays

id_me = myid()
id_other = filter(x -> x != id_me, procs())[rand(1:(nprocs()-1))]

dims = (20,20,20)

if Sys.islinux()
    S = SharedArray{Int64,3}(dims)
    @test startswith(S.segname, "/jl")
    @test !ispath("/dev/shm" * S.segname)

    S = SharedArray{Int64,3}(dims; pids=[id_other])
    @test startswith(S.segname, "/jl")
    @test !ispath("/dev/shm" * S.segname)
end

# TODO : Need a similar test of shmem cleanup for OSX

##### SharedArray tests

function check_pids_all(S::SharedArray)
    pidtested = falses(size(S))
    for p in procs(S)
        idxes_in_p = remotecall_fetch(p, S) do D
            parentindices(D.loc_subarr_1d)[1]
        end
        @test all(sdata(S)[idxes_in_p] .== p)
        pidtested[idxes_in_p] = true
    end
    @test all(pidtested)
end

d = SharedArrays.shmem_rand(1:100, dims)
a = convert(Array, d)

partsums = Vector{Int}(uninitialized, length(procs(d)))
@sync begin
    for (i, p) in enumerate(procs(d))
        @async partsums[i] = remotecall_fetch(p, d) do D
            sum(D.loc_subarr_1d)
        end
    end
end
@test sum(a) == sum(partsums)

d = SharedArrays.shmem_rand(dims)
for p in procs(d)
    idxes_in_p = remotecall_fetch(p, d) do D
        parentindices(D.loc_subarr_1d)[1]
    end
    idxf = first(idxes_in_p)
    idxl = last(idxes_in_p)
    d[idxf] = Float64(idxf)
    rv = remotecall_fetch(p, d,idxf,idxl) do D,idxf,idxl
        @assert D[idxf] == Float64(idxf)
        D[idxl] = Float64(idxl)
        D[idxl]
    end
    @test d[idxl] == rv
end

@test fill(1., 10, 10, 10) == SharedArrays.shmem_fill(1.0, (10,10,10))
@test zeros(Int32, 10, 10, 10) == SharedArrays.shmem_fill(0, (10,10,10))

d = SharedArrays.shmem_rand(dims)
s = SharedArrays.shmem_rand(dims)
copyto!(s, d)
@test s == d
s = SharedArrays.shmem_rand(dims)
copyto!(s, sdata(d))
@test s == d
a = rand(Float64, dims)
@test sdata(a) == a

d = SharedArray{Int}(dims, init = D->fill!(D.loc_subarr_1d, myid()))
for p in procs(d)
    idxes_in_p = remotecall_fetch(p, d) do D
        parentindices(D.loc_subarr_1d)[1]
    end
    idxf = first(idxes_in_p)
    idxl = last(idxes_in_p)
    @test d[idxf] == p
    @test d[idxl] == p
end

d = @inferred(SharedArray{Float64,2}((2,3)))
@test isa(d[:,2], Vector{Float64})

### SharedArrays from a file

# Mapping an existing file
fn = tempname()
write(fn, 1:30)
sz = (6,5)
Atrue = reshape(1:30, sz)

S = @inferred(SharedArray{Int,2}(fn, sz))
@test S == Atrue
@test length(procs(S)) > 1
@everywhere procs(S) begin
    $fill!($S.loc_subarr_1d, $myid())
end
check_pids_all(S)

filedata = similar(Atrue)
read!(fn, filedata)
@test filedata == sdata(S)
finalize(S)

# Error for write-only files
@test_throws ArgumentError SharedArray{Int,2}(fn, sz, mode="w")

# Error for file doesn't exist, but not allowed to create
@test_throws ArgumentError SharedArray{Int,2}(joinpath(tempdir(),randstring()), sz, mode="r")

# Creating a new file
fn2 = tempname()
S = SharedArray{Int,2}(fn2, sz, init=D->D[localindices(D)] = myid())
@test S == filedata
filedata2 = similar(Atrue)
read!(fn2, filedata2)
@test filedata == filedata2
finalize(S)

# Appending to a file
fn3 = tempname()
write(fn3, fill(0x1, 4))
S = SharedArray{UInt8}(fn3, sz, 4, mode="a+", init=D->D[localindices(D)]=0x02)
len = prod(sz)+4
@test filesize(fn3) == len
filedata = Vector{UInt8}(uninitialized, len)
read!(fn3, filedata)
@test all(filedata[1:4] .== 0x01)
@test all(filedata[5:end] .== 0x02)
finalize(S)

# call gc 3 times to avoid unlink: operation not permitted (EPERM) on Windows
S = nothing
@everywhere GC.gc()
@everywhere GC.gc()
@everywhere GC.gc()
rm(fn); rm(fn2); rm(fn3)

### Utility functions

# construct PR #13514
S = @inferred(SharedArray{Int}((1,2,3)))
@test size(S) == (1,2,3)
@test typeof(S) <: SharedArray{Int}
S = @inferred(SharedArray{Int}(2))
@test size(S) == (2,)
@test typeof(S) <: SharedArray{Int}
S = @inferred(SharedArray{Int}(1,2))
@test size(S) == (1,2)
@test typeof(S) <: SharedArray{Int}
S = @inferred(SharedArray{Int}(1,2,3))
@test size(S) == (1,2,3)
@test typeof(S) <: SharedArray{Int}

# reshape

d = SharedArrays.shmem_fill(1.0, (10,10,10))
@test fill(1., 100, 10) == reshape(d,(100,10))
d = SharedArrays.shmem_fill(1.0, (10,10,10))
@test_throws DimensionMismatch reshape(d,(50,))

# rand, randn
d = SharedArrays.shmem_rand(dims)
@test size(rand!(d)) == dims
d = SharedArrays.shmem_fill(1.0, dims)
@test size(randn!(d)) == dims

# similar
d = SharedArrays.shmem_rand(dims)
@test size(similar(d, ComplexF64)) == dims
@test size(similar(d, dims)) == dims

# issue #6362
d = SharedArrays.shmem_rand(dims)
s = copy(sdata(d))
ds = deepcopy(d)
@test ds == d
pids_d = procs(d)
remotecall_fetch(setindex!, pids_d[findfirst(id->(id != myid()), pids_d)::Int], d, 1.0, 1:10)
@test ds != d
@test s != d
copyto!(d, s)
@everywhere setid!(A) = A[localindices(A)] = myid()
@everywhere procs(ds) setid!($ds)
@test d == s
@test ds != s
@test first(ds) == first(procs(ds))
@test last(ds)  ==  last(procs(ds))

# SharedArray as an array
# Since the data in d will depend on the nprocs, just test that these operations work
a = d[1:5]
@test_throws BoundsError d[-1:5]
a = d[1,1,1:3:end]
d[2:4] = 7
d[5,1:2:4,8] = 19

AA = rand(4,2)
A = @inferred(convert(SharedArray, AA))
B = @inferred(convert(SharedArray, copy(AA')))
@test B*A == AA'*AA

d=SharedArray{Int64,2}((10,10); init = D->fill!(D.loc_subarr_1d, myid()), pids=[id_me, id_other])
d2 = map(x->1, d)
@test reduce(+, d2) == 100

@test reduce(+, d) == ((50*id_me) + (50*id_other))
map!(x->1, d, d)
@test reduce(+, d) == 100

@test fill!(d, 1) == fill(1., 10, 10)
@test fill!(d, 2.) == fill(2, 10, 10)
@test d[:] == fill(2, 100)
@test d[:,1] == fill(2, 10)
@test d[1,:] == fill(2, 10)

# Boundary cases where length(S) <= length(pids)
@test 2.0 == remotecall_fetch(D->D[2], id_other, SharedArrays.shmem_fill(2.0, 2; pids=[id_me, id_other]))
@test 3.0 == remotecall_fetch(D->D[1], id_other, SharedArrays.shmem_fill(3.0, 1; pids=[id_me, id_other]))

# Shared arrays of singleton immutables
@everywhere struct ShmemFoo end
for T in [Nothing, ShmemFoo]
    local s = @inferred(SharedArray{T}(10))
    @test T() === remotecall_fetch(x->x[3], workers()[1], s)
end

# Issue #14664
d = SharedArray{Int}(10)
@sync @distributed for i=1:10
    d[i] = i
end

for (x,i) in enumerate(d)
    @test x == i
end

# complex
sd = SharedArray{Int}(10)
se = SharedArray{Int}(10)
@sync @distributed for i=1:10
    sd[i] = i
    se[i] = i
end
sc = convert(SharedArray, complex.(sd,se))
for (x,i) in enumerate(sc)
    @test i == complex(x,x)
end

# Once finalized accessing remote references and shared arrays should result in exceptions.
function finalize_and_test(r)
    finalize(r)
    @test_throws ErrorException fetch(r)
end

for id in [id_me, id_other]
    local id
    finalize_and_test(Future(id))
    finalize_and_test((r=Future(id); put!(r, 1); r))
    finalize_and_test(RemoteChannel(id))
    finalize_and_test((r=RemoteChannel(id); put!(r, 1); r))
end

d = SharedArray{Int}(10)
finalize(d)
@test_throws BoundsError d[1]

# Issue 22139
let
    aorig = a1 = SharedArray{Float64}((3, 3))
    a1 = remotecall_fetch(fill!, id_other, a1, 1.0)
    @test objectid(aorig) == objectid(a1)
    id = a1.id
    aorig = nothing
    a1 = remotecall_fetch(fill!, id_other, a1, 1.0)
    GC.gc(); GC.gc()
    a1 = remotecall_fetch(fill!, id_other, a1, 1.0)
    @test haskey(SharedArrays.sa_refs, id)
    finalize(a1)
    @test !haskey(SharedArrays.sa_refs, id)
end

#14399
let s = convert(SharedArray, [1,2,3,4])
    @test pmap(i->length(s), 1:2) == [4,4]
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
@testset "Sampling" begin
  @testset "RegularSampling" begin
    b = Box(P2(0, 0), P2(2, 2))
    ps = sample(b, RegularSampling(3))
    @test collect(ps) == P2[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)]
    ps = sample(b, RegularSampling(2, 3))
    @test collect(ps) == P2[(0,0),(2,0),(0,1),(2,1),(0,2),(2,2)]

    s = Sphere(P2(0, 0), T(2))
    ps = sample(s, RegularSampling(4))
    ts = P2[(2,0),(0,2),(-2,0),(0,-2)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    s = Sphere(P3(0, 0, 0), T(2))
    ps = sample(s, RegularSampling(2,2))
    ts = P3[(1.7320508075688772, 0.0, 1.0),
            (1.7320508075688772, 0.0, -1.0),
            (-1.7320508075688772, 0.0, 1.0),
            (-1.7320508075688772, 0.0, -1.0)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    b = Ball(P2(0, 0), T(2))
    ps = sample(b, RegularSampling(4,3))
    ts = P2[(1.0, 0.0), (0.0, 1.0), (-1.0, 0.0), (0.0, -1.0),
            (1.5, 0.0), (0.0, 1.5), (-1.5, 0.0), (0.0, -1.5),
            (2.0, 0.0), (0.0, 2.0), (-2.0, 0.0), (0.0, -2.0)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end

    b = Ball(P3(0, 0, 0), T(2))
    ps = sample(b, RegularSampling(3,2,3))
    ts = P3[(0.7071067811865475, 0.0, 0.7071067811865476),
            (1.0, 0.0, 6.123233995736766e-17),
            (0.7071067811865476, 0.0, -0.7071067811865475),
            (-0.7071067811865475, 8.659560562354932e-17, 0.7071067811865476),
            (-1.0, 1.2246467991473532e-16, 6.123233995736766e-17),
            (-0.7071067811865476, 8.659560562354934e-17, -0.7071067811865475),
            (1.0606601717798212, 0.0, 1.0606601717798214),
            (1.5, 0.0, 9.184850993605148e-17),
            (1.0606601717798214, 0.0, -1.0606601717798212),
            (-1.0606601717798212, 1.2989340843532398e-16, 1.0606601717798214),
            (-1.5, 1.8369701987210297e-16, 9.184850993605148e-17),
            (-1.0606601717798214, 1.29893408435324e-16, -1.0606601717798212),
            (1.414213562373095, 0.0, 1.4142135623730951),
            (2.0, 0.0, 1.2246467991473532e-16),
            (1.4142135623730951, 0.0, -1.414213562373095),
            (-1.414213562373095, 1.7319121124709863e-16, 1.4142135623730951),
            (-2.0, 2.4492935982947064e-16, 1.2246467991473532e-16),
            (-1.4142135623730951, 1.7319121124709868e-16, -1.414213562373095)]
    for (p, t) in zip(ps, ts)
      @test p ≈ t
    end
  end

  @testset "UniformSampling" begin
    Random.seed!(2021)
    d = CartesianGrid{T}(100,100)
    s = sample(d, UniformSampling(100))
    μ = mean(coordinates(s, 1:nelements(s)), dims=2)
    @test nelements(s) == 100
    @test isapprox(μ, [50.,50.], atol=T(10))
  end

  @testset "WeightedSampling" begin
    # uniform weights => uniform sampler
    Random.seed!(2020)
    d = CartesianGrid{T}(100,100)
    s = sample(d, WeightedSampling(100))
    μ = mean(coordinates(s, 1:nelements(s)), dims=2)
    @test nelements(s) == 100
    @test isapprox(μ, [50.,50.], atol=T(10))
  end

  @testset "BallSampling" begin
    d = CartesianGrid{T}(100,100)
    s = sample(d, BallSampling(T(10)))
    n = nelements(s)
    X = coordinates(s, sample(1:n, 2, replace=false))
    x, y = X[:,1], X[:,2]
    @test n < 100
    @test sqrt(sum((x - y).^2)) ≥ T(10)

    d = CartesianGrid{T}(100,100)
    s = sample(d, BallSampling(T(20)))
    n = nelements(s)
    X = coordinates(s, sample(1:n, 2, replace=false))
    x, y = X[:,1], X[:,2]
    @test n < 50
    @test sqrt(sum((x - y).^2)) ≥ T(20)
  end

  @testset "Utilities" begin
    # uniform sampling
    d = CartesianGrid{T}(10,10)
    s = sample(d, 50)
    @test nelements(s) == 50
    @test s[1] isa Quadrangle

    # weighted sampling
    d = CartesianGrid{T}(10,10,10)
    s = sample(d, 100, rand([1,2], 1000))
    @test nelements(s) == 100
    @test s[1] isa Hexahedron
  end
end


text:
@inline function accumarray(a,b)
    # TODO: consider blocking for larger sizes (if exceeding 1-2 ms)
    v = [fill(0.0,maximum(a)) for i in 1:Threads.nthreads()]
    Threads.@threads for i in axes(b,1)
         @inbounds v[Threads.threadid()][a[i]] += b[i]
    end
    Threads.nthreads() > 1 ? sum(v) : v[1]
end

@inline function _angular_distances!(VCOORD_th,ang_dist)

    @inbounds begin
        ang_dist[1,:].= @. abs(@views (VCOORD_th[1,:] - VCOORD_th[2,:]))
        ang_dist[2,:].= @. abs(@views (VCOORD_th[2,:] - VCOORD_th[3,:]))
        ang_dist[3,:].= @. abs(@views (VCOORD_th[3,:] - VCOORD_th[1,:]))
    end

    idx = sum(ang_dist .> pi,dims=1).>1
    c   = sum(idx)
    if c > 0
        # idx = findall(x->x==true,idx)
        @inbounds for i in 1:length(idx)
            if idx[i] == true
                j = view(VCOORD_th,:,i) .< pi
                @views VCOORD_th[j,i] .+=  2*pi
            end
        end
    end

end

function elementcoordinate(GlobC,e2n)
    m,n = size(e2n)
    x = Matrix{Float64}(undef,m,n)
    z = Matrix{Float64}(undef,m,n)
    @inbounds for j in 1:n
        idx = view(e2n,:,j)
        for i in 1:m
            ii = idx[i]
            x[i,j] = GlobC[ii].x
            z[i,j] = GlobC[ii].z
        end
    end
    return x,z

end

function fixangles!(x)
    @inbounds for j in axes(x,2)
        # v1,v2,v3 = x[1,j],x[2,j],x[3,j]
        v = ntuple(i->x[i,j],3)
        a1 = abs(v[1] - v[2])
        a2 = abs(v[2] - v[3])
        a3 = abs(v[3] - v[1])

        if (a1+a2+a3) > π
            for k in 1:3
                if v[k] < π
                    x[k,j] += 2π
                end
            end
        end
    end
end

function fixangles6!(x)
    @inbounds for j in axes(x,2)
        v = ntuple(i->x[i,j],6)
        a1 = abs(v[1] - v[2])
        a2 = abs(v[2] - v[3])
        a3 = abs(v[3] - v[1])

        if (a1+a2+a3) > π
            for k in 1:6
                if v[k] < π
                    x[k,j] += 2π
                end
            end
        end
    end
end

getspeed(A::Point2D{Cartesian}) = √(A.x^2 + A.z^2)

getspeed(x::Real, z::Real) = √(x^2 + z^2)

# function calculate_Δt(Ucartesian, nθ, min_inradius)
#     # C = 0.1 # Courant number
#     # r = 1.22 # inner radius
#     # velocity = maximum(getspeed.(Ucartesian))
#     # dx_limit = 2π*r/nθ/4
#     # dt_adv = dx_limit / velocity * C
#     # min(dt_adv, 1.0)

#     min_inradius/maximum(getspeed.(Ucartesian))
# end

calculate_Δt(Ucartesian, nθ, min_inradius) = min_inradius/maximum(getspeed.(Ucartesian))

function getpoints(P)
    n = length(P)
    x = Vector{Float64}(undef,n)
    z = similar(x)

    Threads.@threads for i in eachindex(P)
        @inbounds x[i] = P[i].x
        @inbounds z[i] = P[i].z
    end
    return x,z
end

getpointsX(P) = [@inbounds P[i].x for i in eachindex(P)]
getpointsZ(P) = [@inbounds P[i].z for i in eachindex(P)]

function add_bubble_node(EL2NOD, thr)

    nel = size(EL2NOD,2)
    nnod = maximum(EL2NOD)
    θ, r = deepcopy(thr[1,:]), deepcopy(thr[2,:])
    θel = @views θ[EL2NOD[1:3,:]]
    rel = @views r[EL2NOD[1:3,:]]
    fixangles!(θel)
    θ_new = mean(θel, dims=1)
    r_new = mean(rel, dims=1)
    θ = vcat(θ, vec(θ_new))
    r = vcat(r, vec(r_new))

    return vcat(EL2NOD, transpose(nnod.+(1:nel))), hcat(thr,[θ_new; r_new]), θ, r

end

function getips(EL2NOD, θ, r)
    ndim = 2
    nvert = 3
    nnodel = 6
    nel = size(EL2NOD,2)
    nip = 6
    ni, nn  = Val(nip), Val(nnodel)
    N, dNds, _, w_ip = _get_SF(ni,nn)

    θip = Array{Float64,2}(undef, nel, nip)
    rip = similar(θip)

    @inbounds for iel in 1:nel

        VCOORD_th = SVector{6}(view(θ, :, iel))
        VCOORD_r = SVector{6}(view(r, :, iel))

        for ip in 1:nip
            θip[iel, ip] = (N[ip]*VCOORD_th)[1]
            rip[iel, ip] = (N[ip]*VCOORD_r)[1]
        end
    end

    return θip, rip
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
export OutputContext, @rocprint, @rocprintln

"Internal representation of a static string."
struct DeviceStaticString{N} end
Base.sizeof(dss::DeviceStaticString{N}) where N = N

function Base.unsafe_load(ptr::DevicePtr{DeviceStaticString{N},AS.Global}) where N
    vec_ptr = convert(Ptr{UInt8}, ptr)
    vec_raw = Base.unsafe_wrap(Vector{UInt8}, vec_ptr, (N,))
    idx = findfirst(x->x==0, vec_raw)
    idx = idx === nothing ? N : idx
    return vec_raw[1:idx-1]
end
Base.unsafe_store!(ptr::DevicePtr{<:DeviceStaticString}, x) = nothing

struct OutputContext{HC}
    hostcall::HC
end
function OutputContext(io::IO=stdout; agent=get_default_agent(), buf_len=2^16, kwargs...)
    hc = HostCall(Int64, Tuple{DeviceStaticString{buf_len}}; agent=agent, continuous=true, kwargs...) do bytes
        print(io, String(bytes))
        Int64(length(bytes))
    end
    OutputContext(hc)
end

const GLOBAL_OUTPUT_CONTEXT_TYPE = OutputContext{HostCall{UInt64,Int64,Tuple{DeviceStaticString{2^16}}}}

### macros

macro rocprint(oc, str)
    rocprint(oc, str)
end
macro rocprintln(oc, str)
    rocprint(oc, str, true)
end

macro rocprint(str)
    @gensym oc_ptr oc
    ex = quote
        $(esc(oc_ptr)) = AMDGPU.get_global_pointer(Val(:__global_output_context),
                                                         $GLOBAL_OUTPUT_CONTEXT_TYPE)
        $(esc(oc)) = Base.unsafe_load($(esc(oc_ptr)))
    end
    push!(ex.args, rocprint(oc, str))
    ex
end
macro rocprintln(str)
    @gensym oc_ptr oc
    ex = quote
        $(esc(oc_ptr)) = AMDGPU.get_global_pointer(Val(:__global_output_context),
                                                         $GLOBAL_OUTPUT_CONTEXT_TYPE)
        $(esc(oc)) = Base.unsafe_load($(esc(oc_ptr)))
    end
    push!(ex.args, rocprint(oc, str, true))
    ex
end

### parse-time helpers

function rocprint(oc, str, nl::Bool=false)
    ex = Expr(:block)
    if !(str isa Expr)
        str = Expr(:string, str)
    end
    @assert str.head == :string
    for (idx,arg) in enumerate(str.args)
        if nl && idx == length(str.args)
            arg *= '\n'
        end
        N = rocprint!(ex, 1, oc, arg)
        N = rocprint!(ex, N, oc, '\0')
        dstr = DeviceStaticString{N}()
        push!(ex.args, :(hostcall!($(esc(oc)).hostcall, $dstr)))
    end
    push!(ex.args, :(nothing))
    return ex
end
function rocprint!(ex, N, oc, str::String)
    @gensym str_ptr
    push!(ex.args, :($str_ptr = AMDGPU.alloc_string($(Val(Symbol(str))))))
    push!(ex.args, :(AMDGPU.memcpy!($(esc(oc)).hostcall.buf_ptr+$(N-1), $str_ptr, $(length(str)))))
    return N+length(str)
end
function rocprint!(ex, N, oc, char::Char)
    @assert char == '\0' "Non-null chars not yet implemented"
    byte = UInt8(char)
    ptr = :(Base.unsafe_convert(DevicePtr{UInt8,AS.Global}, $(esc(oc)).hostcall.buf_ptr))
    push!(ex.args, :(Base.unsafe_store!($ptr, $byte, $N)))
    return N+1
end
function rocprint!(ex, N, oc, iex::Expr)
    for arg in iex.args
        N = rocprint!(ex, N, oc, arg)
    end
    return N
end
#= TODO: Support printing arbitrary values?
function rocprint!(ex, N, oc, sym::Symbol)
    push!(ex.args, :($rocprint!($(esc(oc)), $(QuoteNode(sym)))))
    return N+4
end
=#

### runtime helpers


text:
# Struct for allocating the vectors necessary for interpolation.
struct Interpolator{U, V}
    A ::Vector{U} # coefficient for `u`
    Sf::Vector{V} # pseudo-structure of the final state
    ΔS::Vector{V} # difference between initial structure and `Sf`
end

function Interpolator(S₀::T, S₁::T) where {T}
    @assert S₀.grid == S₁.grid

    D₀ = LiquidsDynamics.diffusion_coeff(S₀.f.liquid)
    K  = nodes(S₀.grid)
    Sf = project.(S₁.f.(K))
    ΔS = project.(S₀.f.(K)) .- Sf
    A  = LiquidsDynamics.bsfactors(-2D₀, K, Sf) .* inv.(Sf)

    U = eltype(A)
    V = eltype(Sf)

    return Interpolator{U, V}(A, Sf, ΔS)
end

function interpolate!(S, interpolator, u)
    @unpack A, Sf, ΔS = interpolator
    @. S = muladd(ΔS, exp(reduce_dof(A * u)), Sf)
    return S
end


text:
# MIT License

# Copyright (c) 2021 Fran Urbano
# Copyright (c) 2021 Val Lyashov
# Copyright (c) 2020 Bryan Palmer

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Henderson symmetric weights

# Derive an n-term array of symmetric 'Henderson Moving Average' weights
# formula from ABS (2003), 'A Guide to Interpreting Time Series', page 41.
# returns a numpy array of symmetric Henderson weights indexed from 0 to n-1
function hmaSymmetricWeights(n::Integer)

    m = (n-1)÷2
    m1 = (m+1)^2
    m2 = (m+2)^2
    m3 = (m+3)^2
    d = 8*(m+2)*(m2-1)*(4*m2-1)*(4*m2-9)*(4*m2-25)

    w = map(x -> 315*(m1-x^2)*(m2-x^2)*(m3-x^2)*(3*m2-11*x^2-16)/d, 0:m+1)
    u = vcat(w[end-1:-1:1],w[2:end-1])

    return mod(n, 2) != 0 ? u : vcat(u, missing)

end

# Calculate the asymmetric end-weights
# w --> an array of symmetrical henderson weights (from above function)
# m --> the number of asymmetric weights sought; where m < len(w);
# returns a numpy array of asymmetrical weights, indexed from 0 to m-1;
# formula from Mike Doherty (2001), 'The Surrogate Henderson Filters in X-11',
# Aust, NZ J of Stat. 43(4), 2001, pp901-999; see formula (1) on page 903
function hmaAsymmetricWeights(m::Integer, w::AbstractArray{<:Real})

    n = length(w)

    @assert m <= n "The m argument must be less than w"
    @assert m >= (n-1)÷2 "The m argument must be greater than (n-1)/2"

    sumResidual = sum(w[range(m + 1, n, step = 1)]) / m
    sumEnd =  sum(map(x -> (x-(m+1.0)/2.0)*w[x], m+1:n))

    ic = n < 13 ? 1.0 : (13 <= n < 15 ? 3.5 : 4.5)

    b2s2 = 4/pi/ic^2
    denominator = 1.0 + ((m*(m-1.0)*(m+1.0) / 12.0 ) * b2s2)
    map(r -> w[r] + sumResidual + (((r-1+1.0) - (m+1.0)/2.0) * b2s2) / denominator * sumEnd, 1:m)

end

"""
Package: Forecast

    hma(s, n)

Applies the Henderson moving average filter to dataset `s` with `n`-term.

"Henderson moving averages are filters which were derived by Robert Henderson in 1916
for use in actuarial applications. They are trend filters, commonly used in time series
analysis to smooth seasonally adjusted estimates in order to generate a trend estimate.

They are used in preference to simpler moving averages because they can reproduce
polynomials of up to degree 3, thereby capturing trend turning points.

The ABS uses Henderson moving averages to produce trend estimates from a seasonally
adjusted series. The trend estimates published by the ABS are typically derived using
a 13 term Henderson filter for monthly series, and a 7 term Henderson filter for quarterly series.

Henderson filters can be either symmetric or asymmetric. Symmetric moving averages can be applied
at points which are sufficiently far away from the ends of a time series. In this case, the
smoothed value for a given point in the time series is calculated from an equal number of values
on either side of the data point." - Australian Bureau of Statistics (www.abs.gov.au)

# Arguments
- `s`: Observations' support.
- `n`: Observation values. Tipically 13 or 7 for quarterly data, larger values may need a BigInt type.

# Returns
An array of Henderson filter smoothed values provided in `s`.

# Examples
```julia-repl
julia> hma(rand(1000), BigInt(303)))
1000-element Vector{BigFloat}:
[...]
```
"""
function hma(s::AbstractArray{<:Real}, n::Integer)

    @assert isodd(n) "n must be odd"
    @assert n >= 5 "n must be >= 5"
    @assert length(s) >= n "dataset must be >= than n"

    w = hmaSymmetricWeights(n)
    m = (n-1) ÷ 2
    ls = length(s)

    function hmai(i)
        if i - 1 < m
            u = hmaAsymmetricWeights(m + i, w)[end:-1:1]
            sum(s[1:i + m] .* u)
        elseif i - 1 + m >= ls
            u = hmaAsymmetricWeights(m + ls - i + 1, w)
            sum(s[i-m:ls] .* u)
        else
            sum(s[i-m:i+m] .* w)
        end
    end

    map(x -> hmai(x), 1:ls)
end


text:
module adapters_sqlite_repo_test

using Test # @test
using Octo.Adapters.SQLite # Repo Schema from SELECT FROM WHERE

import SQLite # pathof(SQLite)
dbfile = joinpath(dirname(pathof(SQLite)), "..", "test", "Chinook_Sqlite.sqlite")
dbfile2 = joinpath(@__DIR__, "test.sqlite")
cp(dbfile, dbfile2; force=true)
chmod(dbfile2, 0o666)

Repo.debug_sql()

repo = Repo.connect(
    adapter = Octo.Adapters.SQLite,
    dbfile = dbfile2
)

struct Employee
end
Schema.model(Employee,
    table_name = "Employee",
    primary_key = "EmployeeId"
)

df = Repo.query(Employee)
@test df isa Vector{<:NamedTuple}
@test size(df) == (8,)

df = Repo.get(Employee, 2)
@test size(df) == (1,)

em = from(Employee)
df = Repo.query([SELECT * FROM em WHERE em.EmployeeId == 2])
@test size(df) == (1,)

# using Octo.Adapters.SQLite # CREATE TABLE IF NOT EXISTS
Repo.query([CREATE TABLE IF NOT EXISTS :temp AS SELECT * FROM :Album])

struct Temp
end
Schema.model(Temp, table_name="temp", primary_key="AlbumId")

df = Repo.query(Temp)
@test size(df) == (347,)

changes = (AlbumId=0, Title="Test Album", ArtistId=0)
Repo.insert!(Temp, changes)
df = Repo.query(Temp)
@test size(df) == (348,)

df = Repo.get(Temp, 6)
@test df[1].Title == "Jagged Little Pill"

df = Repo.get(Temp, (Title="Jagged Little Pill",))
@test df[1].Title == "Jagged Little Pill"

changes = (AlbumId=6, Title="Texas")
Repo.update!(Temp, changes)
df = Repo.get(Temp, 6)
@test df[1].Title == "Texas"

Repo.delete!(Temp, changes)
df = Repo.get(Temp, 6)
@test size(df) == (0,)

❔  = Octo.PlaceHolder

a = from(Temp)
df = Repo.query([SELECT * FROM a WHERE a.Title == ❔  AND a.ArtistId == ❔ ], ["The Doors", 140])
@test size(df) == (1,)

Repo.disconnect()

end # module adapters_sqlite_repo_test


text:
## Methods to be used with the Stellar Cluster associated with a Graph.
## Using Lightgraphs

## Zagreb-Randic indices
##

## first Zagreb index
## Normalized here with vertex nmber
function zagreb_first(g, α=2)
    nvertex= nv(g)
    deg= degree(g)
    M1= sum(deg .^ α) / nvertex

    return(M1)
end

## second Zagreb index
##
function zagreb_second(g, α=1 )
    nvertex= nv(g)
    nedge= ne(g)
    edg= edges(g)

    M2= 0
    for e in edg
        d1= degree(g,src(e))
        d2= degree(g,dst(e))
        M2 += (d1*d2)^α
    end

    M2 /= nedge

    return(M2)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
module ACP

include("ACP_MC.jl")
include("BackwardKernel.jl")

using .ACP_MC, .BackwardKernel

#export generate, log_marginal, alloc_dist
export Particle, smc_weight
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
export Exp

"""
`Exp([domainType=Float64::Type,] dim_in::Tuple)`

Creates the exponential non-linear operator with input dimensions `dim_in`:
```math
e^{ \\mathbf{x} }.
```

"""
struct Exp{T,N} <: NonLinearOperator
	dim::NTuple{N,Int}
end

function Exp(DomainType::Type, DomainDim::NTuple{N,Int}) where {N}
	Exp{DomainType,N}(DomainDim)
end

Exp(DomainDim::NTuple{N,Int}) where {N} = Exp{Float64,N}(DomainDim)
Exp(DomainDim::Vararg{Int}) = Exp{Float64,length(DomainDim)}(DomainDim)

function mul!(y::AbstractArray{T,N}, L::Exp{T,N}, x::AbstractArray{T,N}) where {T,N}
	y .= exp.(x)
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{A,TT}},
              b::AbstractArray) where {T,N, A<: Exp{T,N}, TT <: AbstractArray{T,N} }
    L = J.A
    y .= conj.(exp.(L.x)).*b
end

fun_name(L::Exp) = "e"

size(L::Exp) = (L.dim, L.dim)

domainType(L::Exp{T,N}) where {T,N} = T
codomainType(L::Exp{T,N}) where {T,N} = T


text:
function perf_rev_load_slow!(x)
    n = size(x, 1)
    for i = (n - 1):-1:1
        x[i] = x[i + 1]
    end
    return x
end

function perf_rev_load_fast!(x)
    n = size(x, 1)
    x_ = x[n]
    for i = (n - 1):-1:1
        x_ = x[i] = x_
    end
    return x
end

function perf_rev_loadmul_slow!(x, c)
    n = size(x, 1)
    for i = (n - 1):-1:1
        x[i] = c[i + 1] * x[i + 1]
    end
    return x
end

function perf_rev_loadmul_fast!(x, c)
    n = size(x, 1)
    x_ = x[n]
    for i = (n - 1):-1:1
        x_ = x[i] = c[i + 1] * x_
    end
    return x
end


text:
@testset "UmDesign" begin

  y = [[1 0]; [0 1]]
  sc = DataFrame(a=[1,2], b=[3,4])
  oc = DataFrame(c=[5,6,7,8], d=[7,8,9,10])
  dat = UmData(y, sc, oc)

  umd = Unmarked.UmDesign(:test, @formula(ψ~a+b),
                          Unmarked.LogitLink(), dat.site_covs)

  #Test existing dataframes are unchanged
  @test isequal(sc, DataFrame(a=[1,2], b=[3,4]))

  #Check slots
  @test isequal(umd.name, :test)
  @test isequal(umd.formula, @formula(ψ~a+b))
  @test isequal(umd.link, Unmarked.LogitLink())
  @test isequal(umd.data, dat.site_covs)
  @test isequal(umd.coefnames, ["(Intercept)", "a", "b"])
  @test isequal(umd.mat, [[1,1] [1,2] [3,4]])
  @test isnothing(umd.idx)

  #Test getting coefficient indices
  umd2 = Unmarked.UmDesign(:test2, @formula(p~1),
                           Unmarked.LogitLink(), dat.obs_covs)
  Unmarked.add_idx!([umd, umd2])
  @test isequal(umd.idx, 1:3)
  @test isequal(umd2.idx, 4:4)
  @test isequal(Unmarked.get_np([umd, umd2]), 4)

  #Test transform
  tr = logistic.(umd.mat * [0,0.3,0.5])
  @test isequal(Unmarked.transform(umd, [0,0.3,0.5]), tr)

  #Test bad formula
  @test_throws ErrorException Unmarked.UmDesign(:test3, @formula(p~fake),
                                Unmarked.LogitLink(), dat.obs_covs)

end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
module VisualTests

using Gadfly, Interpolations
const nx = 10
xg = 1:nx
xf = -1.2nx:.1:3.3nx

f1(x) = sin(2pi/nx * (x-1))

y1 = map(f1,xg)

p = plot()

if true

for (IT, EB) in (
        (Constant{OnCell}, ExtrapConstant),
        (Constant{OnGrid}, ExtrapConstant),
        (Linear{OnCell}, ExtrapNaN),
        (Linear{OnGrid}, ExtrapNaN),
        (Quadratic{Flat,OnCell}, ExtrapConstant),
        (Quadratic{Flat,OnGrid}, ExtrapConstant),
        (Quadratic{Line,OnCell}, ExtrapLinear),
        (Quadratic{Line,OnGrid}, ExtrapLinear),
        (Quadratic{Line,OnGrid}, ExtrapReflect),
        (Quadratic{Flat,OnCell}, ExtrapReflect),
        (Quadratic{Free,OnCell}, ExtrapReflect),
        (Quadratic{Free,OnGrid}, ExtrapReflect),
        (Quadratic{Periodic,OnGrid}, ExtrapPeriodic),
        (Quadratic{Periodic,OnCell}, ExtrapPeriodic),
    )
    itp = Interpolation(y1, IT(), EB())

    stuff = Any[]

    push!(stuff, layer(x=xg,y=y1,Geom.point,Theme(default_color=color("green"))))
    push!(stuff, layer(x=xf,y=[itp[x] for x in xf],Geom.path,Theme(default_point_size=2px)))
    title = "$(IT.name.name){$(join(map(t -> t.name.name, IT.parameters), ','))}, $(EB.name.name)"
    push!(stuff, Guide.title(title))
    display(plot(stuff...))
end

nx = 10
xg = 1:nx
ny = 6
yg = 1:ny

f2(x,y) = sin(2pi/nx * (x-1)) * cos(2pi/ny * (y-1))
zg = Float64[f2(x,y) for x in xg, y in yg]

xf = -1:.1:nx+1
yf = -1:.1:ny+1

itp2 = Interpolation(zg, Quadratic(Flat(),OnCell()),ExtrapLinear())

display(plot(
    layer(x=xf,y=yf,z=[itp2[x,y] for x in xf, y in yf], Geom.contour),
    Guide.title("Quadratic{Flat,Oncell}, ExtrapLinear")
))

end

end