text:
# Construct the phenotype network for a given parameter setting.
using HDF5, JLD
# Constructs the phenotype network matrix and outlist which is the counts
#   of genotypes that map to the corresponding phenotype before mutation
function construct_pheno_net_parallel( p::Parameters, nreps::Int64, numcircuits::Int64;
    use_lincircuit::Bool=false, csvfile::String="", pheno_file::String="")
  num_processes_per_processor = 4
  println("nreps: ",nreps,"  numcircuits: ",numcircuits)
  funcs = default_funcs(p.numinputs)
  if pheno_file == "" || !isfile(pheno_file)
    ph_net = create_empty_pheno_net( p )
    outlist = zeros(Int64,2^2^p.numinputs)
  elseif isfile(pheno_file)
    (outlist,ph_net)=jldopen(pheno_file,"r") do file
      (read( file, "outlist" ), read( file, "ph_net" ))
     end
  end
  #result_list = Tuple{Matrix{Int64},Vector{Int64}}[]
  pmap_list = nprocs() <= 1 ? [1] : collect(1:(num_processes_per_processor*(nprocs()-1)))
  denom = length(pmap_list)
  println("denom: ",denom,"  pmap_list: ",pmap_list)
  result_list = pmap(x->construct_pheno_net( p, div(nreps,denom), div(numcircuits,denom), use_lincircuit=use_lincircuit ), pmap_list )
  #println("result_list: ",result_list)
  ph_net = create_empty_pheno_net( p )
  outlist = zeros(Int64,2^2^p.numinputs)
  for (phn, outlst) in result_list
    ph_net .+= phn
    outlist .+= outlst
  end
  if pheno_file != ""
    save( pheno_file, "outlist", outlist, "ph_net", ph_net )
  end
  if csvfile != ""
    @assert csvfile[end-3:end] == ".csv"
    mdf = marginals_dataframe( ph_net, outlist, p )
    #println("mdf: ",mdf)
    pn_df = pheno_net_df( ph_net, p )
    ph_file = csvfile[1:end-4] * "_phnet" * ".csv"
    write_csv_file( ph_file, pn_df, funcs, num_processes_per_processor )
    rc_file = csvfile[1:end-4] * "_rowcol" * ".csv"
    write_csv_file( rc_file, mdf, funcs, num_processes_per_processor )
  end
  (ph_net, outlist)
end

# Constructs the phenotype network matrix and outlist which is the counts
#   of genotypes that map to the corresponding phenotype before mutation
# Note:  pmap paralleism failed badly on 9/28/21
function construct_pheno_net( p::Parameters, nreps::Int64, numcircuits::Int64;
    use_lincircuit::Bool=false)
  println("nreps: ",nreps,"  numcircuits: ",numcircuits)
  funcs = default_funcs(p.numinputs)
  ph_net = create_empty_pheno_net( p )
  (outlist, ch_ints_list) = count_outputs( nreps, p, numcircuits, use_lincircuit=use_lincircuit )
  #(new_outlist, ch_ints_list) = count_outputs_parallel( nreps, p, numcircuits, use_lincircuit=use_lincircuit )
  #outlist .+= new_outlist
  #println("length(outlist): ",length(outlist))
  println("length(ch_ints_list): ",length(ch_ints_list))
  #println("ch_ints_list: ",ch_ints_list)
  #new_ph_net_list = pmap(ch_ints->process_ch_ints( ch_ints, p, funcs, use_lincircuit=use_lincircuit ), ch_ints_list ) # pmap commented out 10/4/21
  new_ph_net_list = map(ch_ints->process_ch_ints( ch_ints, p, funcs, use_lincircuit=use_lincircuit ), ch_ints_list )
  for new_ph_net in new_ph_net_list
    ph_net .+= new_ph_net
  end
  (ph_net, outlist)
end

# Helper function for construct_pheno_net()
# Note:  pmap paralleism failed badly on 9/28/21
function process_ch_ints( ch_ints_list::Vector{Int128} , p::Parameters, funcs::Vector{Func};
      use_lincircuit::Bool=false )
  #println("process_ch_ints: ")
  ph_net = create_empty_pheno_net( p )
  for ch_int in ch_ints_list
    #println("ch_int: ",ch_int)
    if use_lincircuit
      #ch_vects = instruction_ints_to_instruction_vects( circuit_int_to_instruction_ints( ch_int, p, funcs ), p, funcs )
      #ch = LinCircuit( ch_vects, p )
      ch = circuit_int_to_circuit( ch_int, p, funcs )
    else
      ch = int_to_chromosome( ch_int, p, funcs )
      #print_circuit(ch)
    end
    src = output_values(ch)[1]   # assumes 1 output
    #println("ch_int: ",ch_int,"  src: ",src)
    dests = map(x->x[1], mutate_all( ch, funcs ))
    #println("dests: ",dests)
    for dest in dests
      ph_net[src+1,dest+1] += 1
      #println("dest: ",dest,"  ph_net[src+1,dest+1]: ",ph_net[src+1,dest+1])
    end
  end
  ph_net
end

function create_empty_pheno_net( p::Parameters; pheno_file::String="" )
  pheno_net = zeros( Int64, p.numoutputs*2^2^p.numinputs, p.numoutputs*2^2^p.numinputs )
  if pheno_file != ""
    save( pheno_file, "pheno_net", pheno_net )
  end
  pheno_net
end

function write_csv_file( csvfile::String, df::DataFrame, funcs::Vector{Func}, num_processes_per_processor::Int64 )
  open( csvfile, "w" ) do f
    hostname = chomp(open("/etc/hostname") do f read(f,String) end)
    println(f,"# date and time: ",Dates.now())
    println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
    print_parameters(f,p,comment=true)
    println(f,"# funcs: ",funcs)
    println(f,"# use_lincircuit: ",use_lincircuit)
    println(f,"# nreps: ",nreps)
    println(f,"# numcircuits: ",numcircuits)
    println(f,"# num_processes_per_processor: ",num_processes_per_processor)
    CSV.write( f, df, append=true, writeheader=true )
  end
end

function marginals_dataframe( ph_net::Array{Int64,2}, outlist::Vector{Int64}, p::Parameters )
  @assert p.numoutputs == 1
  mdf = DataFrame()
  mdf.indices = [ i for i = 0:2^2^p.numinputs-1 ]  # Decimal indices
  mdf.goal = [ @sprintf("0x%04x",i) for i = 0:2^2^p.numinputs-1 ]  # MyInt indices
  mdf.outlist = outlist
  mdf.row_sums = [ sum( ph_net[i,:] ) for i = 1:2^2^p.numinputs ]
  mdf.col_sums = [ sum( ph_net[:,i] ) for i = 1:2^2^p.numinputs ]
  result = evolvability_from_ph_net(ph_net,p)
  mdf.s_evolvability = result[1]
  mdf.d_evolvability = result[2]
  mdf.robustness = result[3]
  mdf
end

# Converts the ph_net matrix to a dataframe.
# The first column is the goals
function pheno_net_df( ph_net::Array{Int64,2}, p::Parameters )
  phdf = DataFrame()
  goals = [ @sprintf("0x%04x",i) for i = 0:2^2^p.numinputs-1 ]
  phdf.goal = goals
  for i = 1:2^2^p.numinputs
    colname = Symbol(goals[i])
    DataFrames.insertcols!(phdf, i+1, colname=>ph_net[:,i] )
  end
  phdf
end

function evolvability_from_ph_net( ph_net::Array{Int64,2}, p::Parameters )
  sigma( k::Int64 ) = ( k == 0 ? 0 : 1 )
  s_evolvability = [ sum( ph_net[i,:] ) - ph_net[i,i]  for i = 1:p.numoutputs*2^2^p.numinputs ]
  d_evolvability = [ sum( map(x->sigma(x), ph_net[i,:])) - sigma(ph_net[i,i])  for i = 1:p.numoutputs*2^2^p.numinputs ]
  robustness = [ ph_net[i,i]/sum( ph_net[i,:] )  for i = 1:p.numoutputs*2^2^p.numinputs ]
  (s_evolvability,d_evolvability,robustness)
end

function read_ph_net_from_csv( csvfile::String, p::Parameters )
  df = read_dataframe( csvfile )
  pheno_net = zeros( Int64, p.numoutputs*2^2^p.numinputs, p.numoutputs*2^2^p.numinputs )
  for i = 1:p.numoutputs*2^2^p.numinputs
    pheno_net[:,i] = df[:,i+1]
  end
  pheno_net
end

# returns the tuple (ph_net,outlist)
function read_jld_file( jld_file::String )
  (outlist,ph_net)=jldopen(jld_file,"r") do file
    (read( file, "ph_net" ), read( file, "outlist" ))
  end
end

# Approximates the Markov chain stationary distribution by post-multiplying init_distribution by the
#   transition matrix niters times.
# The Markov Chain transition matrix is computed by normalizing each row of the phenotype net matrix to sum to 1.
function markov_chain_stationary( niters::Int64, ph_net::Array{Int64,2}, init_distribution::AbstractArray )
  dim = size(ph_net)[1]
  println("dim: ",dim)
  @assert dim == size(ph_net)[2]
  if size(init_distribution)[1] == dim
    init_distribution = transpose(init_distribution)
  end
  @assert size(init_distribution)[1] == 1
  @assert size(init_distribution)[2] == dim
  nonzero_rows = transpose([ (sum(ph_net[i,:]) != 0 ? 1 : 0) for i = 1:dim ])
  init_distribution .*= nonzero_rows   # zero elements of init_distribution that correspond to zero rows of ph_net
  if sum(init_distribution) != 1.0
    init_distribution = init_distribution/sum(init_distribution)
  end
  println("init: ",init_distribution)
  T = zeros(Float64, dim, dim)  # Transition matrix
  for i = 1:dim
    T[i,:] =  sum(ph_net[i,:])>0 ? ph_net[i,:]/sum(ph_net[i,:]) : zeros(Float64,dim)
  end
  #println("T: ",T[1:20,1:20])
  for t = 1:niters
    init_distribution = init_distribution*T
    #println("init: ",init_distribution)
  end
  init_distribution
end


text:
facts("Utilities") do
  context("Units") do
    nb = 1.
    h = 1.m
    v = 1.m^3
    vv = (v, v)
    q = 1.m^3/s
    p = 1.W

    @fact h::Height --> h
    @fact v::Volume --> v
    @fact q::Discharge --> q
    @fact p::Power --> p

    @fact ReservoirManagement._from_unitful(h)::Float64 --> ReservoirManagement._from_unitful(h)
    @fact ReservoirManagement._from_unitful(v)::Float64 --> ReservoirManagement._from_unitful(v)
    @fact ReservoirManagement._from_unitful(q)::Float64 --> ReservoirManagement._from_unitful(q)
    @fact ReservoirManagement._from_unitful(p)::Float64 --> ReservoirManagement._from_unitful(p)
    @fact ReservoirManagement._from_unitful(nb)::Float64 --> ReservoirManagement._from_unitful(nb)
    @fact ReservoirManagement._from_unitful(vv)::Tuple{Float64, Float64} --> ReservoirManagement._from_unitful(vv)
  end

  context("Season separation: SeasonsSeparator") do
    f = y -> DateTime(y, October, 1)
    @fact_throws SeasonsSeparator(Dict(:a => f), Dict(:b => f)) # IDs mismatch
    @fact typeof(SeasonsSeparator(Dict(:a => f), Dict(:a => f))) <: SeasonsSeparator --> true # No error.
  end

  context("Season separation: getters") do
    @fact length(seasons(HydrologicalWesternEurope)) --> 2
    @fact in(:Dry, seasons(HydrologicalWesternEurope)) --> true
    @fact in(:Wet, seasons(HydrologicalWesternEurope)) --> true

    @fact isSeasonKnown(HydrologicalWesternEurope, :Strange) --> false
    @fact isSeasonKnown(HydrologicalWesternEurope, :Dry) --> true
    @fact isSeasonKnown(HydrologicalWesternEurope, :Wet) --> true

    @fact_throws ensureSeasonKnown(HydrologicalWesternEurope, :Strange)
    @fact ensureSeasonKnown(HydrologicalWesternEurope, :Dry)::Void --> ensureSeasonKnown(HydrologicalWesternEurope, :Dry)
    @fact ensureSeasonKnown(HydrologicalWesternEurope, :Wet)::Void --> ensureSeasonKnown(HydrologicalWesternEurope, :Wet)
  end

  context("Season separation: HydrologicalWesternEurope and getters") do
    @fact begins(HydrologicalWesternEurope, :Wet, 2000) --> Date(2000, October, 1)
    @fact ends(HydrologicalWesternEurope, :Wet, 2000) --> Date(2001, April, 30)
    @fact begins(HydrologicalWesternEurope, :Dry, 2000) --> Date(2001, May, 1)
    @fact ends(HydrologicalWesternEurope, :Dry, 2000) --> Date(2001, September, 30)

    @fact length(HydrologicalWesternEurope, :Wet, 2000) --> Day(212)
    @fact length(HydrologicalWesternEurope, :Dry, 2000) --> Day(153)
    @fact length(HydrologicalWesternEurope, :Wet, 2000, Week(1)) --> 30
    @fact length(HydrologicalWesternEurope, :Dry, 2000, Week(1)) --> 22

    bs, es = seasonBoundariesGuarantee(HydrologicalWesternEurope, Week(1), 2000)
    @fact length(bs[:Wet] : Week(1) : es[:Wet]) --> 30
    @fact length(bs[:Dry] : Week(1) : es[:Dry]) --> 22
    @fact length(bs[:Wet] : Week(1) : es[:Dry]) --> 52
  end

  context("Season separation: separation for HydrologicalWesternEurope, no special guarantee") do
    ts = gen_ts(25 * 52, 2000)
    new_dates = [DateTime(2000, January, 1) + Week(i) for i in 0:length(ts.timestamp) - 1]
    ts = TimeArray(new_dates, values(ts), colnames(ts), meta(ts))

    ts_sep = separate(ts, HydrologicalWesternEurope, 2000)
    ts1w = ts_sep[:Wet]
    ts1w_ts = timestamp(ts1w)
    ts1d = ts_sep[:Dry]
    ts1d_ts = timestamp(ts1d)

    @fact ts1w --> not(ts1d)
    for t in ts1w_ts; @fact in(t, ts1d_ts) --> false; end
    for t in ts1d_ts; @fact in(t, ts1w_ts) --> false; end
    for t in ts1w_ts; @fact t >= begins(HydrologicalWesternEurope, :Wet, 2000) --> true; end
    for t in ts1w_ts; @fact t <=   ends(HydrologicalWesternEurope, :Wet, 2000) --> true; end
    for t in ts1d_ts; @fact t >= begins(HydrologicalWesternEurope, :Dry, 2000) --> true; end
    for t in ts1d_ts; @fact t <=   ends(HydrologicalWesternEurope, :Dry, 2000) --> true; end

    # This separator does not have the property that all years are 52-week long.
    for y in 2000:2024
      ts_sep = separate(ts, HydrologicalWesternEurope, 2000)
      @fact length(ts_sep[:Wet]) + length(ts_sep[:Dry]) --> less_than_or_equal(53)
    end
  end

  context("Season separation: separation for HydrologicalWesternEurope, with 52-week year guarantee") do
    ts = gen_ts(25 * 52, 2000)
    new_dates = [DateTime(2000, January, 1) + Week(i) for i in 0:length(ts.timestamp) - 1]
    ts = TimeArray(new_dates, values(ts), colnames(ts), meta(ts))

    ts_sep = separate(ts, HydrologicalWesternEurope, 2000, guaranteeTimeStep=true)
    ts1w = ts_sep[:Wet]
    ts1w_ts = timestamp(ts1w)
    ts1d = ts_sep[:Dry]
    ts1d_ts = timestamp(ts1d)

    @fact ts1w --> not(ts1d)
    for t in ts1w_ts; @fact in(t, ts1d_ts) --> false; end
    for t in ts1d_ts; @fact in(t, ts1w_ts) --> false; end
    for t in ts1w_ts; @fact t >= begins(HydrologicalWesternEurope, :Wet, 2000) - Week(1) --> true; end
    for t in ts1w_ts; @fact t <=   ends(HydrologicalWesternEurope, :Wet, 2000) + Week(1) --> true; end
    for t in ts1d_ts; @fact t >= begins(HydrologicalWesternEurope, :Dry, 2000) - Week(1) --> true; end
    for t in ts1d_ts; @fact t <=   ends(HydrologicalWesternEurope, :Dry, 2000) + Week(1) --> true; end

    # Specific property of HydrologicalWesternEurope: each year is 52-week long, exactly.
    for y in 2000:2024
      ts_sep = separate(ts, HydrologicalWesternEurope, 2000)
      @fact length(ts_sep[:Wet]) + length(ts_sep[:Dry]) --> 52
    end
  end

  context("Non leap years iterator") do
    for y in collect(NonLeapYearIterator(1900, 500)) # Strange rules every 400 years, hence go for 500.
      @fact y --> not(isleapyear)
    end
  end

  context("Making scenarios from read time series") do
    # Generate some hourly data for a few years.
    hours_per_year = [24 * daysinyear(y) for y in 2000:2005]
    vals = [2 + sin(2π * t / 52) + rand() for t in 1:sum(hours_per_year)]
    dats = [DateTime(2000, January, 1) + Hour(i) for i in 0:sum(hours_per_year) - 1]
    TS = TimeArray(dats, vals, ["Discharge (m³/s)"], Hour(1)) # m^3/s, recorded every hour for six calendar years.

    # Check these randomly generated inputs are consistent.
    @fact timestamp(TS) --> unique(timestamp(TS))
    @fact sum(values(TS) .< 0.) --> 0

    r = makeScenarios(TS, scenarioDuration=Year(1), sampleDuration=Day(1), resamplingRate=1,
                          yearBeginning=begins(HydrologicalWesternEurope, :Wet, 2000))

    # Check the results for r. Complex logic for the values, so quite limited for now...
    @fact length(r) --> 5

    first_date = begins(HydrologicalWesternEurope, :Wet, 2000)
    first_value_origs = TS[first_date : Hour(1) : first_date + Hour(23) + Minute(59)]
    @fact values(r[1])[1] --> roughly(sum(values(first_value_origs) * 3600 / 1.e6))
    for i in 1:5
      y = 2000 + i - 1
      beg_date = begins(HydrologicalWesternEurope, :Wet, y)
      end_date = ends(HydrologicalWesternEurope, :Dry, y)
      if isleapyear(y + 1) # The leap day is in the calendar year after y (always in February,
        # while the hydrological year starts in October).
        end_date -= Day(1)
      end

      @fact meta(r[i]) --> Day(1)
      @fact timestamp(r[i])[1] --> beg_date
      @fact timestamp(r[i])[end] --> end_date
    end
  end

  context("Ensuring an array has the right size, with integer copies allowed to reach that size") do
    array = collect(1.:5.)

    @fact toLength(array, 5) --> array
    @fact toLength(array, 4) --> array[1:4]
    @fact_throws toLength(array, 6) # To reach 6, would need to copy parts of the array.
    @fact toLength(array, 10) --> vcat(array, array) # Exactly two copies of the array.
  end
end

ext: jl


text:
@doc """
```math
LFST = \\frac{\\sum^q_{j=1} X^2_{j,T}}{\\sum^q_{j=1}}
```
""" ->
function lfst_test(X, g)
  q = size(X, 1)
  tmp = (collect(1:1:q)*π).^2
  cov_rw = 1./tmp
  cov_ha = ones(q).+(g*g*cov_rw)
  cov_ha_inv = 1./cov_ha
  chol_ha = sqrt(cov_ha_inv)
  Xa = X.*chol_ha
  lfst = sum(X.^2)/sum(Xa.^2)
  nsim = 50000
  Xsim = randn(q, nsim)
  Xasim = Xsim.*repmat(chol_ha, 1, nsim)
  nx = sum(Xsim.^2, 1)
  dx = sum(Xasim.^2, 1)
  lfst_sim = (nx./dx)'
  lfst_pvalue = mean(lfst_sim.>=lfst)
  return lfst, lfst_pvalue
end
function lfst_pv(X, g)
  q = size(X, 1)
  tmp = (collect(1:1:q)*π).^2
  cov_rw = 1./tmp
  cov_ha = ones(q).+(g*g*cov_rw)
  cov_ha_inv = 1./cov_ha
  chol_ha = sqrt(cov_ha_inv)
  Xa = X.*chol_ha
  lfst = sum(X.^2)/sum(Xa.^2)
  nsim = 50000
  Xsim = randn(q, nsim)
  Xasim = Xsim.*repmat(chol_ha, 1, nsim)
  nx = sum(Xsim.^2, 1)
  dx = sum(Xasim.^2, 1)
  lfst_sim = (nx./dx)'
  lfst_pvalue = mean(lfst_sim.>=lfst)
  return lfst, lfst_sim
end

function  mw_lfst_pv(x_t, itrend, param_m, param_t, nsim; stat=false)
# mwlfst_pv - rewrite of Muller and Watson (2015) GAUSS code. Calculates LFST
#   test pvalues.
# Syntax:  pv::Float64 = mwlfst_pv(x_t, itrend, param_m, param_t, nsim)
#
# Inputs:
#   x_t::Vector{Float64} - Discrete cosine transforms of data.
#   itrend::Int64 - 1 if x_t is demeaned, 2 if x_t is detrended.
#   param_m::Float64 - "g" parameter for demeaned alternative.
#   param_t::Float64 - "g" parameter for detrended alternative.
#   nsim::Int64 - Number of simulations to compute p-value
# Outputs:
#   pv::Float64 - p-value
#------------- BEGIN CODE --------------
t = 1000 # Sample Size Used for Asymptotic Approximation
q = size(x_t, 1)
# Compute Null and Alternative Covariance Matrices
if itrend == 1
  psi = psi_compute(t, q)
  cov_h0 = psi'psi
  g_t = param_m/t
  gt2 = g_t.^2
  tmp = mw_rw_cov(t)
  tmp = gt2*tmp
  tmp1 = diag(tmp) + ones(t, 1)
  cov_u = diagrv(tmp, tmp1)
  cov_ha = psi'cov_u*psi
end
# Compute Test Statistic @
  chol_cov_h0 = chol(cov_h0)
  a = inv(chol_cov_h0')
  x_t_trans = a*x_t # @ X transformed under Ho @
  cov_ha_trans = a*cov_ha*a' # @ Covariance under Ha Transformed @
#   @ Note ... LF under Ho is unity @
  lr = mw_marg_lf(x_t_trans, cov_ha_trans)
# @ Simulations for Pvalue @
  xsim = rand(Normal(), size(x_t,1), nsim)
  lr_sim = mw_marg_lf(xsim, cov_ha_trans)
  pv = mean(lr_sim .>(lr - 0.0001))
return lr, lr_sim
end

function diagrv(x, v)
  y = x - diagm(diag(x)) + diagm(v[:])
  return y
end

function mw_marg_lf(x::Vector{Float64}, cov::Array{Float64,2})
# This is a scaled version of the density in King and Hillier (1985) Equation (2) @
# -- Note x is ncov x n ... compute n marginal likelihood values @
# Compute Marginal Likelihood value for x with covariance matrix cov
  x = x./sqrt(sum(x.^2, 1))'
  m = size(x, 1)
  d = det(cov)
  # (gauss code) x = chol(invpd(cov))
  c = chol(Hermitian(inv(x'x)))
  #x = c*x
  tmp = sum(x, 1)
  tmp = d*(tmp.^m)
  mlf = sqrt(1.0./tmp)
  return mlf
end
function mw_marg_lf(x::Array{Float64,2}, cov::Array{Float64,2})
  nsim = size(x, 2)
  mlf_vec = Vector{Float64}(nsim)
  for isim ∈ 1:nsim
    mlf_vec[isim, :] = mw_marg_lf(x[:, isim], cov)
  end
return mlf_vec
end

function mw_rw_cov(t)
  sum_mat = zeros(t,t)
  for i ∈ 1:t
    sum_mat[i, 1:i] = ones(1, i)
  end
  cov = sum_mat*sum_mat'
  return cov
end


text:
# ---
# title: 22. Generate Parentheses
# id: problem22
# author: Qling
# date: 2021-03-08
# difficulty: Medium
# categories: String, Backtracking
# link: <https://leetcode.com/problems/generate-parentheses/description/>
# hidden: true
# ---
#
# Given `n` pairs of parentheses, write a function to _generate all combinations
# of well-formed parentheses_.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3
#     Output: ["((()))","(()())","(())()","()(())","()()()"]
#
#
# **Example 2:**
#
#
#
#     Input: n = 1
#     Output: ["()"]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 8`
#
#
## @lc code=start
using LeetCode

function _generate!(p::String, left::Int, right::Int, res::Vector{String})
    if right == 0
        push!(res, p)
    else
        if left > 0
            _generate!(string(p, "("), left - 1, right, res)
        end

        if right > left
            _generate!(string(p, ")"), left, right - 1, res)
        end
    end
end

function generate_parenthesis(n::Int)::Vector{String}
    res = String[]
    _generate!("", n, n, res)
    return res
end

## @lc code=end


text:
"""
    convex_quadratic(Q, b)

Returns a convex quadratic function `fdf` with signature for `DescentMethods.strong_backtracking!`

    fdf(x, α, d) = f(x + α * d), ∇f(x + α * d)

where

      f(x) = 1/2 xᵀQx - bᵀx
     ∇f(x) = Qx - b
    ∇²f(x) = Q (hessian)

For f(x) to be convex, `Q` must be positive semi-definite matrix.
"""
function convex_quadratic(Q::AbstractMatrix, b::AbstractVector)
    @assert isposdef(Q)

    function f(x::AbstractVector, α::Number, d::AbstractVector)
        y = x + α * d
        y = convert(Vector{Float64}, y)
        val = 0.5 * y' * Q * y - b' * y
        grad = Q * y - b
        return (val, grad)
    end
    return f
end

function random_convex_quadratic(dim::Int)
    Q = rand(test_rng, Float64, (dim, dim))
    Q = Q * Q'  # QQᵀ is always positive semi-definite
    b = rand(test_rng, Float64, (dim, ))
    fdf = convex_quadratic(Q, b)
    return fdf, Q, b
end

"""
    convex_test(fdf, x0, Q[, d])

Perform linear search on a convex quadratic function `fdf` with hessian `Q`.
The search starts at `x0` with descent direction `d`.
Default `d` is opposite to gradient of `fdf`.

Return a tuple `(α, α₀)` where `α` is line search minimizer and `α₀` is analytical minimizer.

# Todo
1. Add randomness to default `d` vector.
"""
function convex_test(fdf, x0, Q, d=nothing)
    y0, grad0 = fdf(x0, 0, x0)  # f(x0 + α x0) = f(x0)
    if isnothing(d)
        d = - 0.01 * grad0  # for greater values of d ⟹ α ≉ α₀
    end

    α = Downhill.strong_backtracking!(fdf, x0, d, y0, grad0)

    "Exact minimizer of convex quadratic function. See Nocedal p. 56, eq. 3.55."
    α₀ = - (grad0' * d) / (d' * Q * d)

    return (α, α₀)
end

@testset "Convex quadratic" begin
    dim = 10
    fdf, Q, _ = random_convex_quadratic(dim)
    for i in 1:20
        x0 = rand(test_rng, Float64, (dim, ))
        α, α₀ = convex_test(fdf, x0, Q)
        println(α, '\t', α₀)
        @test isapprox(α, α₀, rtol=0.05)
    end
end


text:
### ESC LINK
@testset "esclink" begin
    s = Markdown.htmlesc("[hello]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test isnothing(c)
    s = Markdown.htmlesc("[hello][]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test isempty(c)
    s = Markdown.htmlesc("[hello][id]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test isnothing(a)
    @test b == "hello"
    @test c == "id"
    s = Markdown.htmlesc("![hello][id]")
    a,b,c = match(F.ESC_LINK_PAT, s).captures
    @test a == Markdown.htmlesc("!")
    @test b == "hello"
    @test c == "id"
    s = Markdown.htmlesc("[hello]:")
    @test isnothing(match(F.ESC_LINK_PAT, s))
end

### HBLOCK REGEXES

@testset "hb-if" begin
    for s in (
            "{{if var1}}",
            "{{if  var1 }}",
            "{{ if  var1 }}",
        )
        m = match(F.HBLOCK_IF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{if var1 var2}}",
            "{{ifvar1}}"
        )
        m = match(F.HBLOCK_IF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-else" begin
    for s in (
            "{{else}}",
            "{{ else}}",
            "{{  else   }}",
        )
        m = match(F.HBLOCK_ELSE_PAT, s)
        @test !isnothing(m)
    end
end
@testset "hb-elseif" begin
    for s in (
            "{{elseif var1}}",
            "{{else if  var1 }}",
            "{{ elseif  var1 }}",
        )
        m = match(F.HBLOCK_ELSEIF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{else if var1 var2}}",
            "{{elif var1}}"
        )
        m = match(F.HBLOCK_ELSEIF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-end" begin
    for s in (
            "{{end}}",
            "{{ end}}",
            "{{  end   }}",
        )
        m = match(F.HBLOCK_END_PAT, s)
        @test !isnothing(m)
    end
end
@testset "hb-isdef" begin
    for s in (
            "{{isdef var1}}",
            "{{ isdef  var1 }}",
            "{{ isdef  var1 }}",
            "{{ ifdef  var1 }}",
        )
        m = match(F.HBLOCK_ISDEF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{isdef var1 var2}}",
            "{{is def var1}}",
            "{{if def var1}}"
        )
        m = match(F.HBLOCK_ISDEF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-isndef" begin
    for s in (
            "{{isnotdef var1}}",
            "{{ isndef  var1 }}",
            "{{ ifndef  var1 }}",
            "{{ ifnotdef  var1 }}",
        )
        m = match(F.HBLOCK_ISNOTDEF_PAT, s)
        @test m.captures[1] == "var1"
    end
    for s in (
            "{{isnotdef var1 var2}}",
            "{{isnot def var1}}",
            "{{ifn def var1}}"
        )
        m = match(F.HBLOCK_ISNOTDEF_PAT, s)
        @test isnothing(m)
    end
end
@testset "hb-ispage" begin
    for s in (
            "{{ispage var1 var2}}",
        )
        m = match(F.HBLOCK_ISPAGE_PAT, s)
        @test m.captures[1] == "var1 var2"
    end
end
@testset "hb-isnotpage" begin
    for s in (
            "{{isnotpage var1 var2}}",
        )
        m = match(F.HBLOCK_ISNOTPAGE_PAT, s)
        @test m.captures[1] == "var1 var2"
    end
end
@testset "hb-for" begin
    for s in (
            "{{for (v1,v2,v3) in iterate}}",
            "{{for (v1, v2,v3) in iterate}}",
            "{{for ( v1, v2, v3) in iterate}}",
            "{{for ( v1  , v2 , v3 ) in iterate}}"
        )
        m = match(F.HBLOCK_FOR_PAT, s)
        @test isapproxstr(m.captures[1], "(v1, v2, v3)")
    end
    s = "{{for v1 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1")

    # WARNING: NOT RECOMMENDED / NEEDS CARE
    s = "{{for v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1,v2")
    s = "{{for (v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "(v1,v2")
    s = "{{for v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s)
    @test isapproxstr(m.captures[1], "v1,v2)")
end
@testset "hb-toc" begin
    for s in (
            "{{toc}}",
            "{{ toc }}"
            )
        m = match(F.HBLOCK_TOC_PAT, s)
        @test !isnothing(m)
    end
end

# ========
# Checkers
# ========
@testset "ch-for" begin
    s = "{{for v in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test isnothing(F.check_for_pat(m))
    s = "{{for (v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test isnothing(F.check_for_pat(m))
    s = "{{for (v in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
    s = "{{for v1,v2) in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
    s = "{{for v1,v2 in iterate}}"
    m = match(F.HBLOCK_FOR_PAT, s).captures[1]
    @test_throws F.HTMLBlockError F.check_for_pat(m)
end


text:
module Fenkar

using StaticArrays

export fenkar!, fenkar

# define fenkar system
function H(x;k=100.0)
	return 0.5*(1.0 + tanh(k*x))
end
function fenkar!(dx, x, p, t) # in-place, no SA, fast

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi = p

	# fenkar dynamics
	dx[1] = -((x[1]/to)*H(uc-x[1]) + H(x[1]-uc)/tr - (x[2]/td)*H(x[1]-uc)*(1.0-x[1])*(x[1]-uc) - (x[3]/tsi)*H(x[1]-ucsi;k=xk))
	dx[2] = H(uc-x[1])*(1.0-x[2])/(tv1m*H(uv-x[1]) + tv2m*H(x[1]-uv)) - H(x[1]-uc)*x[2]/tvp
	dx[3] = H(uc-x[1])*(1.0-x[3])/twm - H(x[1]-uc)*x[3]/twp

	return nothing
end
function fenkar(x,p,t)		# out-of-place, SA, faster

	# parameters
	@views tsi,tv1m,tv2m,tvp,twm,twp,td,to,tr,xk,uc,uv,ucsi = p

	# fenkar dynamics
	du = -((x[1]/to)*H(uc-x[1]) + H(x[1]-uc)/tr - (x[2]/td)*H(x[1]-uc)*(1.0-x[1])*(x[1]-uc) - (x[3]/tsi)*H(x[1]-ucsi;k=xk))
	dv = H(uc-x[1])*(1.0-x[2])/(tv1m*H(uv-x[1]) + tv2m*H(x[1]-uv)) - H(x[1]-uc)*x[2]/tvp
	dw = H(uc-x[1])*(1.0-x[3])/twm - H(x[1]-uc)*x[3]/twp

	return SA[du,dv,dw]
end

# fenkar system parameters in p[1:13] (BR; from https://doi.org/10.1063/1.166311)
#         tsi,  tv1m,   tv2m,  tvp,  twm,    twp,   td,   to,   tr,   xk,  uc,     uv,  ucsi
#   p = [29.0,  19.6, 1250.0, 3.33, 41.0,  870.0, 0.25, 12.5, 33.3, 10.0, 0.13,  0.04,  0.85]
   p = SA[ 22.0, 333.0,   40.0, 10.0, 65.0, 1000.0, 0.12, 12.5, 25.0, 10.0, 0.13, 0.025,  0.85]
   u0= SA[0.2, 0.5, 0.8]
   tspan = (0.0,1000.0)
end


text:
using Revise
using VCFTools
using SnpArrays
using MendelImpute
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLSO
using BenchmarkTools
using LinearAlgebra
using CSV
using DataFrames

function filter_and_mask(chr::Int)
    # filter chromosome data for unique snps
#         data = "../beagle_raw/chr$chr.1kg.phase3.v5a.vcf.gz"
#         full_record_index = .!find_duplicate_marker(data)
#         VCFTools.filter(data, full_record_index, 1:nsamples(data),
#             des = "chr$chr.uniqueSNPs.vcf.gz", allow_multiallelic=false)

    # import VCF data with only unique SNPs
    _, vcf_sampleID, _, vcf_record_pos, _, _, _ = convert_ht(Bool,
        "chr$chr.uniqueSNPs.vcf.gz", save_snp_info=true, msg="importing")
    total_snps = length(vcf_record_pos)
    samples = length(vcf_sampleID)

    # generate target panel with all snps
    n = 1000
    sample_idx = falses(samples)
    sample_idx[1:n] .= true
    shuffle!(sample_idx)
    VCFTools.filter("chr$chr.uniqueSNPs.vcf.gz", 1:total_snps,
        sample_idx, des = "target.chr$chr.full.vcf.gz", allow_multiallelic=false)

    # generate reference panel without target samples
    VCFTools.filter("chr$chr.uniqueSNPs.vcf.gz", 1:total_snps,
        .!sample_idx, des = "ref.chr$chr.excludeTarget.vcf.gz", allow_multiallelic=false)

    # generate target file with 1000 samples whose snps are in the omni 5.4 chip
    bed = CSV.read("/home/biona001/omni_chips/InfiniumOmni5-4v1-2_A1.bed",
        skipto = 2, header=false, delim='\t', DataFrame)
    typed_snppos = bed[findall(x -> x == "chr$chr", bed[!, 1]), 3]
    match_idx = indexin(vcf_record_pos, typed_snppos)
    record_idx = falses(total_snps)
    record_idx[findall(!isnothing, match_idx)] .= true
    VCFTools.filter("chr$chr.uniqueSNPs.vcf.gz", record_idx, sample_idx,
        des = "target.chr$chr.typedOnly.vcf.gz", allow_multiallelic=false)

    # unphase and mask 0.1% entries in target file
    p = nrecords("target.chr$chr.typedOnly.vcf.gz")
    masks = falses(p, n)
    missingprop = 0.001
    for j in 1:n, i in 1:p
        rand() < missingprop && (masks[i, j] = true)
    end
    mask_gt("target.chr$chr.typedOnly.vcf.gz", masks,
        des="target.chr$chr.typedOnly.masked.vcf.gz", unphase=true)
end
Random.seed!(2020)
@time filter_and_mask(10)

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using InteractionNetworkModels, Distributions, BenchmarkTools
using StatsBase, Plots, StatsPlots

# Testing truncated geom difference
p = 0.7
n = 30
d = TrGeometric(p,n)
a = Geometric(p)
data = [n-rand(d)+rand(a) for i in 1:1000]
bar(counts(data, 0:2*n))

p = 0.8
p_del = TrGeometric(p,0,length(curr))
p_ins = Geometric(p)
η = 0.3
V = 1:10

ind_del = [1,2]
ind_add_1 = [1]
ind_add_2 = [1,6]
vals_1 = [10]
vals_2 = [9,9]
curr = [1,1,1,1,1,1]
prop1 = copy(curr)
prop2 = copy(curr)

# In place (with two proposals given)
@time split_noisy!(
    curr,
    prop1, prop2,
    ind_del,
    ind_add_1, ind_add_2,
    vals_1, vals_2,
    η, V
)
curr,prop1,prop2

# When given current value and storage
curr = [1,2,1,2,1,2]
prop1 = copy(curr) # Current val
prop2 = Int[] # Storage for new val
@time split_noisy!(
    prop1,
    prop2,
    ind_del,
    ind_add_1, ind_add_2,
    vals_1, vals_2,
    η, V
)
prop1, prop2

# Standard
split_noisy(
    curr,
    ind_del,
    ind_add_1, ind_add_2,
    vals_1, vals_2,
    η, V
)

# With auxiliary sampling
curr = [1,2,1,2,1,2]
p = 0.7
p_del = TrGeometric(p,0,length(curr))
p_ins = Geometric(p)
η = 0.7
V = 1:10
prop1, prop2 = (copy(curr), copy(curr))
# In place
rand_split_noisy!(
    curr,
    prop1, prop2,
    p_del, p_ins,
    η, V
)
curr, prop1, prop2

# With alloc
rand_split_noisy(
    curr,
    p_del, p_ins,
    η, V
)

# ========
# Merging
# =======

# With keep indexing
curr1 = [1,1,1,1]
curr2 = [2,2,2,2]
ind1_keep = [1,2,3]
ind2_keep = [1,3,4]
ind_add = [1,4]
vals = [10,10]
prop = zeros(Int, length(ind1_keep))
@time merge_noisy!(
    curr1, curr2,
    prop,
    ind1_keep,ind2_keep,
    ind_add, vals,
)
curr1, curr2, prop

@time merge_noisy(
    curr1, curr2,
    ind1_keep, ind2_keep,
    ind_add, vals
)

# With deletion indexing
ind1_del = [1,4]
ind2_del = [1]
ind_add = [1,4]
vals = [10,10]
curr1 = [2,1,2,1]
curr2 = [2,2,2]
@time merge_noisy_ins_dels!(
    curr1, curr2,
    ind1_del, ind2_del,
    ind_add, vals
)
curr1, curr2

# Random
p = 0.8
V = 1:10
η = 0.9
curr1 = [1,1,1]
curr2 = [2,2,2,2]
p_ins = Geometric(p)
p_del = TrGeometric(p,0,min(length(curr1), length(curr2)))
@time rand_merge_noisy_ins_dels!(
    curr1, curr2,
    p_del, p_ins,
    V,
    η
)
curr1, curr2

# =====================
# Multiple split/merges
# =====================

# Multiple splits
ind_split = [1,3]
p = 0.5
p_ins = Geometric(p)
η = 0.5
V = 1:10
S_curr = [[1,1,1], [2,2,2], [3,3,3]]
S_store = [Int[] for i in 1:4]
@time multiple_adj_split_noisy!(
    S_curr,
    S_store,
    ind_split,
    p_ins, η, V
)
S_curr
S_store

# Multiple merges
ind_merge = [1,3]
p = 0.9
p_ins = Geometric(p)
η = 0.1
V = 1:10
S_curr = [[1,1,1], [1,1,2], [3,3,3], [2,3,2], [3,3,3]]
S_store = [Int[] for i in 1:2]
@time multiple_adj_merge_noisy!(
    S_curr,
    S_store,
    ind_merge,
    p_ins, η, V
)
S_curr
S_store

# With loc sampling
# ------------------

# Split
p = 0.4
p_ins = Geometric(p)
η = 0.7
V = 1:10
ν_td = 10
S_curr = [[1,1,1], [2,2,2], [3,3,3]]
S_store = [Int[] for i in 1:4]
@time rand_multiple_adj_split_noisy!(
    S_curr,
    S_store,
    ν_td,
    p_ins, η, V
)
S_curr
S_store

# Merge
p = 0.99
p_ins = Geometric(p)
η = 0.1
V = 1:10
ν_td = 10
S_curr = [[1,1,1], [2,2,2], [3,3,3]]
S_store = [Int[] for i in 1:4]
@time rand_multiple_adj_merge_noisy!(
    S_curr,
    S_store,
    ν_td,
    p_ins, η, V
)
S_curr
S_store

# Subsequence sampling
n = 10
k = 5
ind = zeros(Int,k)
@btime StatsBase.seqsample_a!(1:n, ind)
@btime StatsBase.seqsample_c!(1:n, ind)
@btime StatsBase.seqsample_d!(1:n, ind)

StatsBase.seqsample_a!(1:n, ind)
ind
StatsBase.seqsample_c!(1:n, ind)
ind
StatsBase.seqsample_d!(1:n, ind)
ind


text:
using Base.Broadcast: BroadcastFunction, broadcasted, materialize

ChainRulesCore.@non_differentiable Base.getindex(m::AbstractAttenMask, I::Integer...)
ChainRulesCore.@non_differentiable Base.getindex(m::MaskIndexer, I::Integer...)
ChainRulesCore.@non_differentiable Base.getindex(m::AbstractAttenMask, I::Tuple)
ChainRulesCore.@non_differentiable Base.getindex(m::MaskIndexer, I::Tuple)
ChainRulesCore.@non_differentiable (::Type{<:AbstractAttenMask})(args...)
ChainRulesCore.@non_differentiable (::Type{<:AbstractAttenMaskOp})(args...)
ChainRulesCore.@non_differentiable getmask(arg...)

function ChainRulesCore.rrule(::typeof(apply_mask), op::NaiveAttenMaskOp, mask, score)
    m = as_bool(randomness(mask)) ? getmask(mask, score) : mask
    naive_apply_mask_pullback(Ȳ) = (NoTangent(), NoTangent(), NoTangent(), @thunk unthunk(Ȳ) .* m)
    return score .* m, naive_apply_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), f, mask, score, scale)
    tmp = getmask(mask, score, scale)
    Y, back = rrule_via_ad(config, broadcast, f, tmp, score)
    function fallback_apply_broadcast_mask_pullback(Ȳ)
        Ȳs = (back∘unthunk)(Ȳ)
        return (NoTangent(), Ȳs[2], NoTangent(), Ȳs[4], NoTangent())
    end
    Y, fallback_apply_broadcast_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), ::typeof(+), mask, score, scale)
    apply_broadcast_mask_pullback(Ȳ) = (NoTangent(), NoTangent(), NoTangent(), Ȳ, NoTangent())
    apply_broadcast_mask(+, mask, score, scale), apply_broadcast_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), ::typeof(*), mask, score, scale)
    rnm = as_bool(randomness(mask))
    m = rnm ? getmask(mask, score, scale) : mask
    function apply_broadcast_mask_pullback(Ȳ)
        thk = @thunk begin
            tmp = unthunk(Ȳ)
            rnm ? tmp .* m : @. tmp * (mask * scale)
        end
        return (NoTangent(), NoTangent(), NoTangent(), thk, NoTangent())
    end
    fwd = rnm ? score .* m : apply_broadcast_mask(*, mask, score, scale)
    fwd, apply_broadcast_mask_pullback
end

@init @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    Zygote.unbroadcast(x::AbstractAttenMask, _) = nothing

    function Zygote._pullback(ctx::Zygote.AContext, ::typeof(Broadcast.broadcasted), ::MaskStyle{M}, f, args...) where M
        return y, ∇broadcasted = Zygote._pullback(ctx, Broadcast.broadcasted, M(), f, args...)
    end
end


text:
function dynamics(ev::CPEGWorkspace, x::SVector{7,T}, u::SVector{1,W}) where {T,W}

    # scaled variables
    r_scaled = x[SA[1,2,3]]
    v_scaled = x[SA[4,5,6]]
    σ = x[7]

    # unscale
    r, v = unscale_rv(ev.scale,r_scaled,v_scaled)

    # altitude
    h = norm(r) - ev.params.gravity.R

    # density
    ρ = density(ev.params.density, h)

    # lift and drag magnitudes
    L, D = LD_mags(ev.params.aero,ρ,r,v)

    # basis for e frame
    e1, e2 = e_frame(r,v)

    # drag and lift accelerations
    D_a = -(D/norm(v))*v
    L_a = L*sin(σ)*e1 + L*cos(σ)*e2

    # gravity
    g = gravity(ev.params.gravity,r)

    # acceleration
    ω = ev.planet.ω
    a = D_a + L_a + g - 2*cross(ω,v) - cross(ω,cross(ω,r))

    # rescale units
    v,a = scale_va(ev.scale,v,a)

    return SA[v[1],v[2],v[3],a[1],a[2],a[3],u[1]*ev.scale.uscale]
end

function rk4(
    ev::CPEGWorkspace,
    x_n::SVector{7,T},
    u::SVector{1,W},
    dt_s::Float64) where {T,W}

    k1 = dt_s*dynamics(ev,x_n,u)
    k2 = dt_s*dynamics(ev,x_n+k1/2,u)
    k3 = dt_s*dynamics(ev,x_n+k2/2,u)
    k4 = dt_s*dynamics(ev,x_n+k3,u)
    return (x_n + (1/6)*(k1 + 2*k2 + 2*k3 + k4))
end

function rollout(ev::CPEGWorkspace,x0::SVector{7,T},U::Vector{SVector{1,T}}) where T
    """everything in and out of the function is scaled"""

    # scaled dt
    dt_s = ev.dt/ev.scale.tscale

    # input U
    U_in = copy(ev.U)

    N = 1000
    X = [@SVector zeros(length(x0)) for i = 1:N]
    U = [@SVector zeros(length(U_in[1])) for i = 1:N]
    X[1] = x0
    end_idx = NaN
    for i = 1:N-1
        U[i] = (i>length(U_in)) ? SA[0.0] : U_in[i]

        X[i+1] = rk4(ev,X[i],U[i],dt_s)

        # for debugging purposes
        # @show (norm(X[i+1][1:3])*ev.scale.dscale - ev.params.gravity.R)/1e3

        if (norm(X[i+1][1:3])*ev.scale.dscale - ev.params.gravity.R) < 10e3
            # @info "hit alt"
            # @show i
            end_idx = i+1
            break
        end
    end

    # trim relavent
    if isnan(end_idx)
        error("didn't hit the altitude during the rollout")
    end

    X = X[1:end_idx]
    # ev.U = U[1:(end_idx-1)]

    return X, U[1:(end_idx-1)]
end

function get_jacobians(
    ev::CPEGWorkspace,
    X::Vector{SVector{7,T}}) where T

    dt_s = ev.dt/ev.scale.tscale

    N = length(X)
    A = [ForwardDiff.jacobian(_x->rk4(ev,_x,ev.U[i],dt_s),X[i]) for i = 1:N-1]
    B = [ForwardDiff.jacobian(_u->rk4(ev,X[i],_u,dt_s),ev.U[i]) for i = 1:N-1]
    return A,B
end

# let
#
#     ev = CPEGWorkspace()
#
#     # # @show ev.scale.uscale
#     #
#     # x = @SVector randn(7)
#     # u = SA[1.45]
#     # #
#     # @btime dynamics($ev,$x,$u)
#     #
#     # @btime rk4($ev,$x,$u,0.1)
#     # #
#     # A = ForwardDiff.jacobian(_x -> dynamics(ev,_x,u),x)
#     # # @btime ForwardDiff.jacobian(_x -> dynamics($ev,_x,$u),$x)
#     # @btime ForwardDiff.jacobian(_x -> rk4($ev,_x,$u,0.1),$x)
#     # ForwardDiff.jacobian(_x -> rk4(ev,_x,u,0.1),x)
#     Rm = ev.params.gravity.R
#     r0 = SA[Rm+125e3, 0.0, 0.0] #Atmospheric interface at 125 km altitude
#     V0 = 5.845e3 #Mars-relative velocity at interface of 5.845 km/sec
#     γ0 = -15.474*(pi/180.0) #Flight path angle at interface
#     v0 = V0*SA[sin(γ0), cos(γ0), 0.0]
#     σ0 = deg2rad(90)
#
#     r0sc,v0sc = scale_rv(ev.scale,r0,v0)
#
#     dt = 1.0/ev.scale.tscale
#
#     x0 = SA[r0sc[1],r0sc[2],r0sc[3],v0sc[1],v0sc[2],v0sc[3],σ0]
#
#     # @show norm(r0)
#     #
#     # @show norm(r0)*ev.scale.dscale - ev.params.gravity.R
#     N = 100
#     U = [SA[0.0] for i = 1:N-1]
#     #
#     X,U = rollout(ev, x0, U, dt)
#     @show length(X)
#     # N = 100
#     # X = [@SVector zeros(7) for i = 1:N]
#
#     # @show typeof([r0;v0;σ0])
#     # X[1] = SA[r0[1],r0[2],r0[3],v0[1],v0[2],v0[3],σ0]
#     # @btime $X[1] = SA[$r0[1],$r0[2],$r0[3],$v0[1],$v0[2],$v0[3],$σ0]
#
#     # @show X[1]
#     # @show typeof(X[1])
#
#     # i = 7
#     # @show ForwardDiff.jacobian(_x->rk4(ev,_x,U[i],dt),X[i])
#     # @show ForwardDiff.jacobian(_u->rk4(ev,X[i],_u,dt),U[i])
#     # A,B= get_jacobians(ev,X,U,dt)
#     #
#     # N = length(X)
#     # Xm = zeros(7,N)
#     # for i = 1:N
#     #     Xm[:,i] = X[i]
#     # end
#     #
#     # nr = [norm(X[i][1:3]) for i = 1:N]
#     # nv = [norm(X[i][4:6]) for i = 1:N]
#     # mat"
#     # figure
#     # hold on
#     # plot($Xm(1:3,:)')
#     # hold off
#     # "
#     #
#     # mat"
#     # figure
#     # hold on
#     # plot($Xm(4:6,:)')
#     # hold off
#     # "
#     # mat"
#     # figure
#     # hold on
#     # plot($nr)
#     # plot($nv)
#     # hold off"
#
#     # X = unscale_X(ev.scale,X)
#     # # @show typeof(X[1][SA[1,2,3]])
#     #
#     # alt, dr, cr = postprocess(ev,X,[r0;v0])
#     #
#     #
#     # mat"
#     # figure
#     # hold on
#     # plot($alt/1e3)
#     # hold off "
#     #
#     # mat"
#     # figure
#     # hold on
#     # plot($dr/1e3,$cr/1e3)
#     # hold off"
#
# end


text:
const GraphMatrix = Union{Matrix, BitMatrix}

function clique_size(seq::Union{Vector{Bool}, BitArray}, A::GraphMatrix)
  set = find(seq)
  clique = true
  for i in set, j in set
    if (i != j) && !(A[i,j] > 0)
      clique = false
      break
    end
  end
  return clique*length(set)
end

function clique_loss(seq::Union{Vector{Bool}, BitArray}, A::GraphMatrix, kappa::Float)
  set = find(seq)
  clique = true
  edges = Float(0)
  n = length(set)
  for i in set, j in set
    if i != j
      aij = (A[i,j] > 0)
      edges += aij
      clique = clique && aij
    end
  end

  return (Float(1) - Float(edges)/Float(max(n*(n-1+kappa), 1)), Float(-n*clique))
end

function clique_loss(seq::Vector{Int}, A::GraphMatrix, kappa::Float)
  n = size(A, 1)
  if length(seq) < n
    idx = falses(n)
    idx[seq] = true
  else
    idx = seq .> 1
  end
  return clique_loss(idx, A, kappa)
end

"""
### Clique finding using the cakewalk method
Minimize the soft-clique-loss function with cakewalk. Usage: \n
`opt_seq, opt_val = clique_cw(A, kappa)` \n
where `kappa` is a size booster parameter beteen 0.0 and 1.0.
For details about the loss, kappa, and the reasoning behind them see our paper.\n\n

When `opt_seq[i]` is 2, it means vertex `i` is selected, and 1 means it is not. \n
`opt_val` is a tuple of two floats. The first is the value being optimized, i.e.,
the soft-clique-loss, and the second is the actual loss that we are interested in,
the negative clique size.
"""
clique_cw(A::GraphMatrix, args...; kwargs...) = cakewalk(clique_loss, 2, size(A, 1), A, args...; kwargs...)

inc_clique_size(ind, A::GraphMatrix) = maximum(map(v -> (w = copy(ind); w[v] = true; clique_size(w, A)), 1:size(A,1)))

"""
### Validate that a clique is locally maximal
First, find some clique, \n
`opt_seq, opt_val = clique_cw(A, kappa)` \n
Then, validate the clique is locally maximal, \n
`is_local_max_clique(opt_seq .== 2, A)`\n
The returned value is binary.
"""
is_local_max_clique(ind, A::GraphMatrix) = (cs = clique_size(ind,A); (cs > 0) && (cs >= inc_clique_size(ind, A)))


text:
"""
This script serves the Pollen.jl documentation on a local file server
so that it can be loaded by the frontend in development mode.
files should be stored.

    > julia docs/serve.jl

Use `./make.jl` to export the generated documents to disk.

There are two modes for interactive development: Lazy and Regular.
In lazy mode, each document will be built only if it is requested in
the frontend, while for Regular mode, each document will be built
once before serving.
"""

using Pollen

project = include("project.jl")

Pollen.serve(
    project;
    lazy = get(ENV, "POLLEN_LAZY", "false") == "true",
    port = Base.parse(Int, get(ENV, "POLLEN_PORT", "8000")),
    format = JSONFormat()
)


text:
using Pkg

Pkg.develop(PackageSpec(path=pwd()))
Pkg.instantiate()

using Documenter, GenomicFeatures

makedocs(
    format = Documenter.HTML(
        edit_link = :commit
    ),
    modules = [GenomicFeatures],
    sitename = "GenomicFeatures.jl",
    pages = [
        "Home" => "index.md",
        "Intervals" => "man/intervals.md",
        "API Reference" => "api/api.md"
    ],
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ) * ", The BioJulia Organisation, and other contributors."
)
deploydocs(
    repo = "github.com/BioJulia/GenomicFeatures.jl.git",
    devbranch = "develop",
    push_preview = true
)


text:
export iauTttai
"""
Time scale transformation:  Terrestrial Time, TT, to International
Atomic Time, TAI.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tt1,tt2    double    TT as a 2-part Julian Date

Returned:
   tai1,tai2  double    TAI as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Note:

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tai1,tai2 follow
   suit.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTttai(double tt1, double tt2, double *tai1, double *tai2)

function iauTttai(tt1::Real, tt2::Real)
   ref_tai1 = Ref{Float64}(0.0)
   ref_tai2 = Ref{Float64}(0.0)

   status = ccall((:iauTttai, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tt1), convert(Float64, tt2), ref_tai1, ref_tai2)

   return status, ref_tai1[], ref_tai2[]
end


text:
#!/usr/bin/env julia
using CSV, DataFrames, PowerModels, Tables, ArgParse
using JuMP, Ipopt, Mosek, MosekTools
using Statistics, Distributions, Measures, StatsBase
using SparseArrays, LinearAlgebra
using CDDLib, Polyhedra

# cd to directory
cd(dirname(@__FILE__))
PowerModels.silence()

# load functions
include("scr/aux_fun.jl")
include("scr/pwr_fun.jl")
include("scr/gas_fun.jl")

# simulation settings
set = Dict(
# gas and power networks
:gas_case => "case_48", :pwr_case => "pglib_opf_case118_ieee.m",
# control horizon data
:T => 5, :k_t => [1 4 7 10 13],
# wind data
:w̅ => 100, :σ => args["σ²"],
# number of uncertainty samples
:Ns => args["Ns"],
# constraint violation tolerance
:ε̅_p => args["ε̅_p"], :ε̅_g => args["ε̅_g"],
# variability parameters
:α_ϱ => args["α_ϱ"], :α_ϑ => args["α_ϑ"], :α_ψ => args["α_ψ"],
# netowrk topology opt data
:b_valve_edge => [21 30], :M̅ => 10e3, :M̲ => 10e3,
)

## solve CC OPF
net = load_pwr_data(set)
# wind power forecast
ξ = wind_forecast(set)
# solve OPF models
sol_opf = solve_chance_constrained_OPF(net,set,ξ)
# save gas consumption at gas-power coupling points
CSV.write("data/gas_uncertainty.csv",  Tables.table(net[:Λ]*net[:M_p]*sol_opf[:P][set[:T]]), writeheader=false)

## solve CC gas network optimization
net             = load_gas_data(set)
# solve non-convex gas network optimization
sol_non_convex  = non_convex_opt(net)
# linearize Weymouth equation
lin_res         = linearization(net,sol_non_convex[:model])
# solve the linearized model
sol_lin         = linearized_opt(net,lin_res)
# optimize chance-constrained gas network optimization
sol_cc  = chance_constrained_gas_opt(net,lin_res,set)
ofs_res_sto = out_of_sample(sol_cc,net,set,ξ)

## display some and export all results
println("expected cost:                \$$(round(sol_cc[:cost]/1000,digits=1))×10^3")
println("pressure variability measure:   $(round(sol_cc[:ϱ_std_tot]/1000,digits=1))×10^3")
open("$(outdir)/summary.txt","a") do io
   println(io,"simulation settings:")
   println(io,"α_ϱ = $(args["α_ϱ"]), α_ϑ = $(args["α_ϑ"]), α_ψ = $(args["α_ψ"]), σ = $(args["σ²"])")
   println(io,"Ns = $(args["Ns"]), ε̅_p = $(args["ε̅_p"]), ε̅_g = $(args["ε̅_g"])")
   println(io,"results:")
   println(io,"exp_cost = ", sol_cc[:cost])
   println(io,"variability = ", sol_cc[:ϱ_std_tot])
   println(io,"mean_comp_deployment = ", ofs_res_sto[:mean_compr_deployment])
   println(io,"mean_valve_deployment = ", ofs_res_sto[:mean_valve_deployment])
   println(io,"mass_con_violation_exp = ", ofs_res_sto[:mass_violation_exp])
   println(io,"press_con_violation_exp = ", ofs_res_sto[:press_violation_exp])
end


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Read in data

delim = ';'
df = CSV.read(joinpath("..", "..", "data", "rugged.csv"), DataFrame; delim)
df = filter(row -> !(ismissing(row[:rgdppc_2000])), df)
df.log_gdp = log.(df.rgdppc_2000)
df.cont_africa = Array{Float64}(convert(Array{Int}, df.cont_africa))

Base.@kwdef mutable struct RuggedModel{Ty <: AbstractVector,
  Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::RuggedModel)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :rugged], df[:, :cont_africa],
  df[:, :rugged] .* df[:, :cont_africa]);
model = RuggedModel(;y=df[:, :log_gdp], x=x)

# Model callable with *a single argument*.

function (problem::RuggedModel)(θ)
    @unpack y, x = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(0, 100), x[1])
    ll += logpdf(Normal(0, 10), x[2])
    ll += logpdf(Normal(0, 10), x[3])
    ll += logpdf(Normal(0, 10), x[4])
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0, 1.0, 2.0], σ = 1.0))
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

p = as_particles(posterior)
display(p)

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.22 0.14  9.00  9.46   282    1
bR    -0.21 0.08 -0.33 -0.08   275    1
bA    -1.94 0.24 -2.33 -1.59   268    1
bAR    0.40 0.14  0.18  0.62   271    1
sigma  0.96 0.05  0.87  1.04   339    1
"

# End of `08/m8.1s.jl`


text:
using LinearAlgebra
using Plots

########################
## Obstacle Avoidance ##
########################
model = Dynamics.dubinscar[1]
n,m = model.n, model.m

x0 = [0.0;0.0;0.]
xf = [10.0;10.0;0.]
tf =  5.0
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

u_bnd = 4.
u_min = [-u_bnd; -u_bnd]
u_max = [u_bnd; u_bnd]
x_min = [-Inf; -Inf; -Inf]
x_max = [Inf; Inf; Inf]

r = 1.0
circles = ((1.,2.,r),(5.,5.,r),(8.,9,r),(1.5,5,r),(7.5,5,r))
n_circles = length(circles)

function cI(c,x,u)
    for i = 1:n_circles
        c[i] = circle_constraint(x,circles[i][1],circles[i][2],circles[i][3])
    end
    c
end

obj_con = ConstrainedObjective(obj,u_min=u_min, u_max=u_max,x_min=x_min,x_max=x_max,cI=cI)

opts = SolverOptions()
opts.verbose = true
opts.use_gradient_aula = true
opts.constraint_tolerance = 1e-7
opts.cost_tolerance = 1e-7
opts.penalty_initial = 0.01
opts.penalty_scaling = 10
opts.cost_tolerance_intermediate = 1e-2
opts.outer_loop_update_type = :default

solver = Solver(model, obj_con, integration=:rk3, N=301, opts=opts)
U0 = rand(m,solver.N)
X0 = line_trajectory(solver)
X0_rollout = rollout(solver,U0)

solver.opts.use_nesterov = true
res, stats = solve(solver,U0)
stats["iterations"]
maximum(norm.(res.d,Inf))
plot(res.U)

solver2 = Solver(solver)
solver2.opts.square_root = true
solver2.opts.use_gradient_aula = false
solver2.opts.use_nesterov = false
res2, stats2 = solve(solver2,U0)

solver2.opts.use_nesterov = true
res3, stats3 = solve(solver2,U0)

plot(stats["c_max"],yscale=:log10)
plot!(stats2["c_max"],yscale=:log10)
plot!(stats3["c_max"])

# Obstacle Avoidance
plt = plot(title="Obstacle Avoidance",aspect_ratio=:equal)
plot_obstacles(circles)
plot_trajectory!(to_array(res.X),width=2,color=:green,label="Constrained")
plot_trajectory!(to_array(res2.X),width=2,color=:blue,label="Constrained")
plot_trajectory!(to_array(res3.X),width=2,color=:yellow,label="Constrained")

display(plt)


text:
module NelderMead

include("Vertexes.jl")
include("Simplexes.jl")

"""
    optimise(f, initial_vertex_positions; kwargs...)

Find minimum of function, `f`, first creating a Simplex from vertices at
`initial_vertex_positions`, and options passed in via kwargs.
"""
function optimise(f::T, initial_vertex_positions::U; kwargs...
   ) where {T<:Function, W<:Number, V<:AbstractVector{W}, U<:AbstractVector{V}}
  return optimise!(Simplex(f, initial_vertex_positions), f; kwargs...)
end

"""
    optimise(f, initial_positions, initial_step; kwargs...)

Find minimum of function, `f`, first creating a Simplex using a starting
vertex position, `initial_position`, and other vertices `initial_step` away from
that point in all directions, and options passed in via kwargs.
"""
function optimise(f::T, initial_position::AbstractVector{U},
               initial_step::AbstractVector{V}; kwargs...
               ) where {T<:Function, U<:Number, V<:Number}
  return optimise!(Simplex(f, initial_position, initial_step), f; kwargs...)
end

"""
    optimise(f, lower, upper, gridsize; kwargs...)

Find minimum of function, `f`, first creating a `2 .* prod(gridsize)` Simplices
within a hyper-rectangle spanning from `lower` to `upper`, and
options passed in via kwargs.
"""
function optimise(f::F, lower::AbstractVector{T}, upper::AbstractVector{T},
    gridsize::AbstractVector{<:Integer}; kwargs...) where {F<:Function, T<:Number}

  all(gridsize .> 0) || throw(ArgumentError("gridsize, $gridsize, must .> 0"))
  if !(length(lower) == length(upper) == length(gridsize))
    throw(ArgumentError("The lengths of lower $lower, upper $upper, and
                        gridsize $gridsize must be the same"))
  end
  if !all(lower .< upper)
    throw(ArgumentError("lower $lower must .< upper $upper"))
  end

  dim = length(lower)
  index2position(i) = (i .- 1) ./ gridsize .* (upper .- lower) .+ lower
  index2values = Dict()
  totaltime = @elapsed for ii ∈ CartesianIndices(Tuple(gridsize .+ 1))
    index = collect(Tuple(ii))
    x = index2position(index)
    index2values[index] = f(x)
  end
  U = typeof(first(index2values)[2])
  simplices = Set{Simplex{T, U}}()
  function generatesimplices!(simplices, direction)
    for ii ∈ CartesianIndices(Tuple(((gridsize .+ 1) .* ones(Int, dim))))
      vertices = Vector{Vertex{T, U}}()
      index = collect(Tuple(ii))
      for i ∈ 1:dim + 1
        vertexindex = [index[j] + ((j == i) ? direction : 0) for j ∈ 1:dim]
        all(1 .<= vertexindex .<= gridsize .+ 1) || continue
        vertex = Vertex{T, U}(index2position(vertexindex),
                              index2values[vertexindex])
        push!(vertices, vertex)
      end
      length(vertices) != dim + 1 && continue
      push!(simplices, Simplex(vertices))
    end
  end

  totaltime += @elapsed generatesimplices!(simplices, 1)
  totaltime += @elapsed generatesimplices!(simplices, -1)
  @assert length(simplices) == 2 * prod(gridsize)
  totaltime += @elapsed solutions = map(s->optimise!(s, f; kwargs...), collect(simplices))
  return solutions
end

"""
    optimise!(s, f; kwargs...)

Find minimum of function, `f`, starting from Simplex, `s`, with options
passed in via kwargs.

# Keyword Arguments
-  stopval (default sqrt(eps())): stopping criterion when function evaluates
equal to or less than stopval
-  xtol_abs (default zeros(T)) .* ones(Bool, dimensionality(s)): stop if
the vertices of simplex get within this absolute tolerance
-  xtol_rel (default eps(T)) .* ones(Bool, dimensionality(s)): stop if
the vertices of simplex get within this relative tolerance
-  ftol_abs (default zero(real(U))): stop if function evaluations at the
vertices are close to one another by this absolute tolerance
-  ftol_rel (default 1000eps(real(U))): stop if function evaluations at the
vertices are close to one another by this relative tolerance
-  maxiters (default 1000): maximum number of iterations of the Nelder Mead
algorithm
-  timelimit (default Inf): stop if it takes longer than this in seconds
-  α (default 1): Reflection factor
-  β (default 0.5): Contraction factor
-  γ (default 2): Expansion factor
-  δ (default 0.5): Shrinkage factor
"""
function optimise!(s::Simplex{T,U}, f::F; kwargs...) where {F<:Function, T<:Real, U}
  kwargs = Dict(kwargs)
  stopval = get(kwargs, :stopval, sqrt(eps()))
  xtol_abs = get(kwargs, :xtol_abs, zeros(T)) .* ones(Bool, dimensionality(s))
  xtol_rel = get(kwargs, :xtol_rel, eps(T)) .* ones(Bool, dimensionality(s))
  ftol_abs = get(kwargs, :ftol_abs, zero(real(U)))
  ftol_rel = get(kwargs, :ftol_rel, 1000eps(real(U)))
  maxiters = get(kwargs, :maxiters, 1000)
  timelimit = get(kwargs, :timelimit, Inf)
  α = get(kwargs, :α, 1)
  β = get(kwargs, :β, 0.5)
  γ = get(kwargs, :γ, 2)
  δ = get(kwargs, :δ, 0.5)

  @assert α >= 0 "$α >= 0"
  @assert 0 <= β < 1 "0 <= $β < 1"
  @assert γ > 1 "$γ > 1"
  @assert γ > α "$γ > $α"

  reflect(this, other) = Vertex(newposition(this, α, other), f)
  expand(this, other) = Vertex(newposition(this, -γ, other), f)
  contract(this, other) = Vertex(newposition(this, -β, other), f)
  shrink(this, other) = Vertex(newposition(this, δ, other), f)

  function shrink!(s::Simplex)
    lengthbefore = length(s)
    best = bestvertex(s)
    newvertices = [shrink(best, v) for v ∈ s if !isequal(v, best)]
    remove!(s, findall(v->!isequal(v, best), s.vertices))
    push!(s, newvertices...)
    sort!(s, by=v->value(v))
    @assert length(s) == lengthbefore
    return nothing
  end

  iters, totaltime = 0, 0.0
  returncode = assessconvergence(s,
    xtol_abs, xtol_rel, ftol_abs, ftol_rel, stopval)
  history = deepcopy(s.vertices)
  while returncode == :CONTINUE && iters < maxiters && totaltime < timelimit
    totaltime += @elapsed begin
      iters += 1
      best = bestvertex(s)
      worst = worstvertex(s)
      secondworst = secondworstvertex(s)
      reflected = reflect(centroidposition(s), worst)

      if any(h->isequal(h, reflected), history)
        returncode = :ENDLESS_LOOP
        break
      end
      history .= circshift(history, 1)
      history[1] = reflected

      if best <= reflected < secondworst
        swapworst!(s, reflected)
      else
        centroid = Vertex(centroidposition(s), f)
        if reflected < best
          expanded = expand(centroid, reflected)
          expanded < reflected && swapworst!(s, expanded)
          expanded >= reflected && swapworst!(s, reflected)
        elseif secondworst <= reflected < worst
          contracted = contract(centroid, reflected)
          contracted <= reflected ? swapworst!(s, contracted) : shrink!(s)
        elseif reflected >= worst
          contracted = contract(centroid, worst)
          contracted < worst ? swapworst!(s, contracted) : shrink!(s)
        end
      end
      returncode = assessconvergence(s,
        xtol_abs, xtol_rel, ftol_abs, ftol_rel, stopval)
    end
  end

  iters == maxiters && (returncode = :MAXITERS_REACHED)
  best = bestvertex(s)
  return value(best), position(best), returncode, iters
end # optimise

end


text:
#!/usr/bin/env julia

"""
alignment_overlaps.jl - An alternative implemention of the BioAlignments.jl
BAM overlap.jl code to convert Record types into fragment-based and read-based Interval types

Original code at https://github.com/BioJulia/BioAlignments.jl/blob/master/src/bam/overlap.jl

By: Shaun Adkins (sadkins@som.umaryland.edu)
"""

import BGZFStreams

### OverlapIterator - Extended version of original BAM OverlapIterator

struct OverlapIterator{T}
    reader::BAM.Reader{T}
    refname::String
    interval::UnitRange{Int}
    strand_type::String
    max_frag_size::UInt
end

function GenomicFeatures.eachoverlap(reader::BAM.Reader, interval::Interval, strand_type::String="no", max_frag_size::UInt=1000)
    return GenomicFeatures.eachoverlap(reader, interval.seqname, interval.first:interval.last, strand_type, max_frag_size)
end

function GenomicFeatures.eachoverlap(reader::BAM.Reader, interval, strand_type::String="no", max_frag_size::UInt=1000)
    return GenomicFeatures.eachoverlap(reader, convert(Interval, interval), strand_type, max_frag_size)
end

function GenomicFeatures.eachoverlap(reader::BAM.Reader, refname::AbstractString, interval::UnitRange, strand_type::String="no", max_frag_size::UInt=1000)
    return OverlapIterator(reader, String(refname), interval, strand_type, max_frag_size)
end

### SuperBAMRecord - BAM Record with additional information

struct SuperBAMRecord
    record::BAM.Record
    interval::Interval{Bool}
    strand::Char
    #alignmenttype::T where {T<:AbstractAlignment}
end

#function SuperBAMRecord(record::BAM.Record, interval::Interval, strand::Char)
#    return SuperBAMRecord(record, interval, strand, getalignmenttype(metadata(interval)))
#end

record(s::SuperBAMRecord) = s.record
interval(s::SuperBAMRecord) = s.interval
strand(s::SuperBAMRecord) = s.strand
alignment_type(s::SuperBAMRecord) = getalignmenttype(metadata(interval(s)))

# Iterator
# --------

function Base.iterate(iter::OverlapIterator)
    refindex = findfirst(isequal(iter.refname), iter.reader.refseqnames)
    if refindex == 0
        throw(ArgumentError("sequence name $(iter.refname) is not found in the header"))
    end
    @assert iter.reader.index !== nothing
    chunks = Indexes.overlapchunks(iter.reader.index.index, refindex, iter.interval)
    if !isempty(chunks)
        seek(iter.reader, first(chunks).start)
    end
    state = BAM.OverlapIteratorState(refindex, chunks, 1, BAM.Record())
    return iterate(iter, state)
end

function Base.iterate(iter::OverlapIterator, state)
    while state.chunkid ≤ lastindex(state.chunks)
        chunk = state.chunks[state.chunkid]

        # Ran into issue where the BAM.Reader processed all records but kept going which
        # lead to the stream block_index to exceed the number of blocks.  This solves that.
        if eof(iter.reader.stream)
            break
        end

        while BGZFStreams.virtualoffset(iter.reader.stream) < chunk.stop
            read!(iter.reader, state.record)
            # Determine if fragment or read alignment.  Fragment alignments may overlap the interval whereas the single read may be to the right
            alignmentinterval = get_alignment_interval(state.record, iter.max_frag_size, iter.strand_type)
            alignmentinterval === nothing && continue

            c = GenomicFeatures.compare_overlap(alignmentinterval, Interval(iter.refname, iter.interval), isless)
            if c == 0
                alignmentstrand = getstrand(alignmentinterval, isstranded(iter.strand_type))
                return SuperBAMRecord(copy(state.record), alignmentinterval, alignmentstrand), state
            elseif c > 0
                # no more overlapping records in this chunk since records are sorted
                break
            end

        end
        state.chunkid += 1
        if state.chunkid ≤ lastindex(state.chunks)
            seek(iter.reader, state.chunks[state.chunkid].start)
        end
    end
    return nothing
end

max_stars_repo_name: IGS/FADU

max_issues_repo_name: adkinsrs/FADU

max_forks_repo_name: adkinsrs/FADU


text:
module Requestor

using Artifacts: load_artifacts_toml
using DocStringExtensions: METHODLIST
using HTTP: get
using JSON: parse
using PkgUtility: terror

export request_LUT

"""
This function requests data from our server for supported datasets (only latest datasets supported). Supported methods are

$(METHODLIST)
"""
function request_LUT end;

"""
    request_LUT(artname::String, lat::Number, lon::Number, cyc::Int = 0; user::String="Anonymous", interpolation::Bool = false, server::String = "tofu.gps.caltech.edu", port::Int = 5055)

Request data from the server, given
- `artname` Artifact full name such as `LAI_MODIS_2X_8D_2017_V1`
- `lat` Latitude
- `lon` Longitude
- `cyc` Cycle index, default is 0 (read entire time series)
- `user` User name (non-registered users need to wait for 5 seconds before the server processes the request)
- `interpolation` If true, interpolate the data linearly
- `server` Server address
- `port` Port number for the GriddingMachine server

---
# Examples
```julia
request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5);
request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5; interpolation=true);
request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5, 8);
request_LUT("LAI_MODIS_2X_8D_2017_V1", 30.5, 115.5, 8; interpolation=true);
```
"""
request_LUT(artname::String, lat::Number, lon::Number, cyc::Int = 0; user::String="Anonymous", interpolation::Bool = false, server::String = "tofu.gps.caltech.edu", port::Int = 5055) = (
    # make sure the artifact is within our collection
    _metas = load_artifacts_toml(joinpath(@__DIR__, "../Artifacts.toml"));
    _artns = [_name for (_name,_) in _metas];
    @assert artname in _artns;

    # send a request to our webserver at tofu.gps.caltech.edu:5055 and translate it back to Dictionary
    _response = get("http://$(server):$(port)/request.json?user=$(user)&artifact=$(artname)&lat=$(lat)&lon=$(lon)&cyc=$(cyc)&interpolate=$(interpolation)");
    _json_str = String(_response.body);
    _results  = parse(_json_str);

    # if there is no result item in the Dictionary
    if !haskey(_results, "Result")
        @show _results;
        @error terror("There is something wrong with the request, please check the details about it!");
    end;

    # replace -9999 to NaN
    _data = _results["Result"];
    _std = _results["Error"];
    if typeof(_data) <: Array
        _data[_data .== -9999] .= NaN;
        _std[_std .== -9999] .= NaN;
    else
        if _data == -9999 _data = NaN; end;
        if _std == -9999 _std = NaN; end;
    end;

    # if length of data is 1 return a Number, else return a Vector
    if length(_data) > 1
        return Float32.(_data), Float32.(_data)
    end;

    return Float32(_data[1]), Float32(_std[1])
);

end


text:
module Ops21

import Knet.Ops20
using AutoGrad: AutoGrad, @primitive
using SpecialFunctions: erf

const dropout = Ops20.dropout; export dropout
const bmm = Ops20.bmm; export bmm
include("activation.jl"); export elu, gelu, relu, selu, sigm
include("mmul.jl"); export mmul

end


text:
runge(x) = 1 / (1 + 25 * x^2)

function test_fh_runge(tol=1e-8)
    xx = [-1.0:0.05:1.0;]
    yy = runge.(xx)
    fh = FHInterp(xx, yy, order=3)
    xxx = [-1.0:0.01:1.0;]
    vl = fh.(xxx)
    rf = runge.(xxx)
    return norm(vl - rf) < 5e-5
end

# Ok, so this is kinda brittle..
function test_abs_x(order=3)
    pts = [10, 20, 40, 80, 160, 320, 640]
    tol = [0.5, 0.2, 0.09, 0.05, 0.03, 0.01, 0.005]
    pass = true
    for (p, t) in zip(pts, tol)
        xx = collect(range(-5.0, 5.0, length=2p-1))
        xi = xx[1:2:end]
        xt = xx[2:2:end]
        yy = abs.(xi)
        fh = FHInterp(xi, yy, order=order, grid=true)
        err = maximum(abs.(fh.(xt) .- abs.(xt)))
        pass *= err < t
    end
    return pass
end


text:
using ..TopOpt: GPU, @init_cuda
import ..TopOpt: whichdevice
using ..CUDASupport: CUDASupport
using ..GPUUtils

@init_cuda

@define_cu(FilterMetadata, :cell_neighbouring_nodes, :cell_node_weights)
whichdevice(m::FilterMetadata) = whichdevice(m.cell_neighbouring_nodes)

@define_cu(SensFilter, :metadata, :nodal_grad, :last_grad, :cell_weights)
whichdevice(c::SensFilter) = whichdevice(c.nodal_grad)

function update_nodal_grad!(nodal_grad::CuVector, node_cells, args...)
    T = eltype(nodal_grad)
    allargs = (nodal_grad, node_cells.offsets, node_cells.values, args...)
    callkernel(dev, cheq_kernel1, allargs)
    CUDAdrv.synchronize(ctx)
    return nothing
end

function cheq_kernel1(
    nodal_grad,
    node_cells_offsets,
    node_cells_values,
    cell_weights,
    cells,
    cellvolumes,
    black,
    white,
    varind,
    grad,
)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        nodal_grad[n] = zero(T)
        cell_weights[n] = zero(T)
        r = node_cells_offsets[n]:(node_cells_offsets[n + 1] - 1)
        for i in r
            c = node_cells_values[i][1]
            if black[c] || white[c]
                continue
            end
            ind = varind[c]
            w = cellvolumes[c]
            cell_weights[n] += w
            nodal_grad[n] += w * grad[ind]
        end
        n += offset
    end
end

function normalize_grad!(nodal_grad::CuVector, cell_weights)
    T = eltype(nodal_grad)
    args = (nodal_grad, cell_weights)
    callkernel(dev, cheq_kernel2, args)
    CUDAdrv.synchronize(ctx)
    return nothing
end
function cheq_kernel2(nodal_grad, cell_weights)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        w = cell_weights[n]
        w = ifelse(w > 0, w, one(T))
        nodal_grad[n] /= w
        n += offset
    end
    return nothing
end

function update_grad!(
    grad::CuVector,
    black,
    white,
    varind,
    cell_neighbouring_nodes,
    cell_node_weights,
    nodal_grad,
)
    T = eltype(grad)
    allargs = (
        grad,
        black,
        white,
        varind,
        cell_neighbouring_nodes.offsets,
        cell_neighbouring_nodes.values,
        cell_node_weights.values,
        nodal_grad,
    )
    callkernel(dev, cheq_kernel3, allargs)
    CUDAdrv.synchronize(ctx)
    return nothing
end
function cheq_kernel3(
    grad,
    black,
    white,
    varind,
    cell_neighbouring_nodes_offsets,
    cell_neighbouring_nodes_values,
    cell_node_weights_values,
    nodal_grad,
)
    T = eltype(nodal_grad)
    i = @thread_global_index()
    offset = @total_threads()
    while i <= length(black)
        if black[i] || white[i]
            continue
        end
        ind = varind[i]
        r =
            cell_neighbouring_nodes_offsets[ind]:(cell_neighbouring_nodes_offsets[ind + 1] - 1)
        length(r) == 0 && continue
        grad[ind] = zero(T)
        sum_weights = zero(T)
        for linear_ind in r
            w = cell_node_weights_values[linear_ind]
            sum_weights += w
            node_ind = cell_neighbouring_nodes_values[linear_ind]
            grad[ind] += nodal_grad[node_ind] * w
        end
        grad[ind] /= sum_weights
        i += offset
    end

    return nothing
end


text:
using Random
using StatsBase
using BenchmarkTools

# Initializations
T = 1000;                           # sample size
s = zeros(Int, 2);                  # anchors
c = [1 + (i - 1) % 4 for i in 1:T]; # vector of cluster indicators
s̄_clu = zeros(Int, T);              # s̄ (cluster indicators)
s̄_idx = zeros(Int, T);              # s̄ (indices)
σ = collect(1:T);

# Update s
function update_s!(s, T)
    StatsBase.sample!(1:T, s, replace=false)
end
update_s!(s, T);

# Update s̄
function update_s̄!(s̄_clu, s̄_idx, s, c)
    card_s̄ = 0
    T = length(c)
    resize!(s̄_clu, T)
    resize!(s̄_idx, T)
    c1, c2 = c[s[1]], c[s[2]]
    for i ∈ 1:T
        i ∈ s && continue
        (c[i] == c1) && (card_s̄ += 1; s̄_clu[card_s̄] = 1; s̄_idx[card_s̄] = i)
        (c[i] == c2) && (card_s̄ += 1; s̄_clu[card_s̄] = 2; s̄_idx[card_s̄] = i)
    end
    card_s̄ += 2
    s̄_clu[card_s̄ - 1] = c1; s̄_idx[card_s̄ - 1] = s[1]
    s̄_clu[card_s̄ - 0] = c2; s̄_idx[card_s̄ - 0] = s[2]
    resize!(s̄_clu, card_s̄)
    resize!(s̄_idx, card_s̄)
    return card_s̄
end
card_s̄ = update_s̄!(s̄_clu, s̄_idx, s, c);

# Update σ
function update_σ!(σ, card_s̄)
    resize!(σ, card_s̄ - 2);
    randperm!(σ);
    resize!(σ, card_s̄);
    σ[card_s̄ - 1] = card_s̄ - 1;
    σ[card_s̄ - 0] = card_s̄ - 0;
    return σ
end
update_σ!(σ, card_s̄);

# Tests
@btime update_s!($s, $T);
@btime update_s̄!($s̄_clu, $s̄_idx, $s, $c);
@btime update_σ!($σ, $card_s̄);


text:
import .DayCounts: yearfraction, Actual365, DateRoll, adjust
import .TermStructure: YieldTermStructure, ConstantContinuousYieldCurve, discount, startdate, frequency, compounding, daycount

struct SingleStock <: Contract
end

abstract type AbstractCoreModel <: AbstractModel

end

"""
    CoreModel(startdate, startprice, interestrate, carryrate)

Contains the core parameters for the model of `SingleStock`:
* `startprice`: initial price at `startdate`
* `yieldcurve`: risk free rate of return for the numeraire.
* `carrycurve`: the carry rate, i.e. the net income for holding the asset:
  - for stocks this is typically positive (i.e. dividends)
  - for commodities this is typically negative (i.e. cost-of-carry)

The `yieldcurve` and `carryrate` are all specified on a continously compounded, Actual/365 basis.

The exact dynamics of the model are unspecified, but has the property that

```math
E[S_t] = S_0 \\times e^{(r-q)t}
```

"""
struct CoreModel{T, S<:YieldTermStructure, R<:YieldTermStructure} <: AbstractCoreModel
    startprice::T
    yieldcurve::S
    carrycurve::R
end

CoreModel(startprice::T, yieldcurve::S, carryrate::Float64=0.0) where {T, S} =
    CoreModel(startprice, yieldcurve,
              ConstantContinuousYieldCurve(daycount(yieldcurve), carryrate, startdate(yieldcurve)))

CoreModel(startdate::Date, startprice::T, yieldrate::V, carryrate::Float64=0.0) where {T,V} =
    CoreModel(startprice,
              ConstantContinuousYieldCurve(Actual365(), yieldrate, startdate),
              ConstantContinuousYieldCurve(Actual365(), carryrate, startdate))

numeraire(m::CoreModel) = unit(m.startprice)
startdate(m::CoreModel) = startdate(m.yieldcurve)

yearfractionto(m::CoreModel, dt::Date) =
    yearfraction(daycount(m.yieldcurve),  startdate(m.yieldcurve),  dt)

function value(m::CoreModel, c::WhenAt{Receive{T}}) where T
    value(m, c.c) * discount(m.yieldcurve, maturitydate(c))
end

value(m::CoreModel, c::SingleStock) = m.startprice
function value(m::CoreModel, c::WhenAt{SingleStock})
    m.startprice * discount(m.carrycurve, maturitydate(c))
end

function forwardprice(m::CoreModel, ::SingleStock, dt::Date)
    m.startprice * discount(m.carrycurve, dt) / discount(m.yieldcurve, dt)
end

"""A model for the time value of money
"""
struct YieldModel{S<:YieldTermStructure, T<:DateRoll, U<:HolidayCalendar} <: AbstractModel
    yieldcurve::S
    dateroll::T
    holidaycalendar::U
end

value(m::YieldModel, c::When{At, Receive{T}}) where {T} = value(m, c.c) * discount(m.yieldcurve, adjust(m.dateroll, m.holidaycalendar, maturitydate(c)))


text:
include("../src/Dataset.jl")
import .Dataset:
  get_MINST,
  TrainTest,
  split_X_by_batches

include("../src/Model.jl")
import .Model:
  split_encoder_result,
  create_vae,
  model_sample,
  conv_MINST_model

include("../src/Utils.jl")
import .Utils:
  gen_images,
  save_model,
  load_model

using Test
using Flux
using Flux.Tracker: TrackedReal
using Printf
using Distributions: Uniform

temp_dir = tempdir()

@testset "Model save/load" begin
  n_sample = Int(floor(rand(Uniform(10,100))))
  n_latent = Int(floor(rand(Uniform(1,50)))) * 2

  ps, loss_fn, f, g = create_vae(n_sample, n_latent)
  file = tempname()
  save_model(f, g, file)

  fp, gp = load_model
  @test true
end

@testset "Image Utilities" begin
  n_sample = Int(floor(rand(Uniform(10,100))))
  n_latent = Int(floor(rand(Uniform(1,50)))) * 2
  dataset = get_MINST(n_sample)
  ps, loss_fn, f, g = create_vae(n_latent, n_sample)

  # zero initialized μ̂, logσ̂
  outfile_zero_latent = joinpath(temp_dir, "zero_latent.png")
  @test typeof(model_sample(f)) == BitArray{4}
  gen_images(outfile_zero_latent, f)
  @test isfile(outfile_zero_latent)
  @sprintf("saved file to %s", outfile_zero_latent)

  # f(g(x)) test
  outfile_X_latent = joinpath(temp_dir, "X_latent.png")
  gen_images(outfile_X_latent, g, f, dataset.test_x)
  @test isfile(outfile_X_latent)
  @sprintf("saved file to %s", outfile_X_latent)
end

@testset "ADAM optimization can run" begin
  n_sample = Int(floor(rand(Uniform(10,100))))
  n_latent = Int(floor(rand(Uniform(1,50)))) * 2
  dataset = get_MINST(n_sample)
  X = dataset.train_x
  ps, loss_fn, f, g = create_vae(n_latent, n_sample)
  opt = ADAM()

  @test typeof(loss_fn(X)) == TrackedReal{Float32}
  X = float.(X .> 0.5)
  Flux.train!(loss_fn, ps, zip([X]), opt)
  @test true == true
end

@testset "Convolution and transpose is isomorphic" begin
  dataset = get_MINST()

  n_sample = 100
  n_latent = 10
  X = dataset.train_x[:,:,:,1:n_sample]
  f, g = conv_MINST_model(n_latent)

  X_transformed = g(reshape(X, 28, 28, 1, n_sample))

  x_mean, x_std = split_encoder_result(X_transformed, n_latent)

  Xp = f(x_mean)
  @test size(Xp) == size(X)
end

@testset "MINST dataset: size/shape okay" begin
  ds = get_MINST()
  img_shape = (28, 28, 1)
  n_train = 60000
  n_test = 10000
  @test size(ds.train_x) == (img_shape..., n_train)
  @test size(ds.train_y) == (n_train, )
  @test size(ds.test_x) ==  (img_shape..., n_test)
  @test size(ds.test_y) == (n_test, )
end

max_stars_repo_name: aniketjivani/variational-autoencoders

max_issues_repo_name: aniketjivani/variational-autoencoders

max_forks_repo_name: aniketjivani/variational-autoencoders


text:
filename = "audio"
extension = ".wav"
filepath = joinpath("..", "data", filename * extension)
@assert(isfile(filepath),
        "The test audio wave file does not exist in path $path.")
channel = 1 # the test audio file is mono
signal = Signal(filepath, channel = channel)

## get properties

framerate = getframerate(signal)
@test isapprox(framerate, 44100)
data = getdata(signal)
@test typeof(data) <: Array{Complex,1}

## get chunks
seconds = countseconds(signal)
@test seconds == 4

chunkspersecond = rand(24:30)
chunksnumber = countchunks(signal, chunkspersecond)
@test chunksnumber > 1

i = rand(1:chunksnumber)
chunk = getchunk(signal, i, chunkspersecond = chunkspersecond)
@test typeof(chunk) <: Signal
@test length(getdata(chunk)) > 1

## analyse and synthesize
frequencysignal = analyse(chunk)
@test typeof(getdata(frequencysignal)) <: Array{Complex,1}
timesignal = synthesize(frequencysignal)
@test typeof(getdata(timesignal)) <: Array{Complex,1}

## perform DFT and FFT
dft = performdft(chunk)
@test typeof(getdata(dft)) <: Array{Complex,1}
minisignal = Signal(getdata(signal)[1:128], framerate)
fft = performfft(minisignal)
@test typeof(getdata(fft)) <: Array{Complex,1}

samplesnumber = rand(3:5)
frequencies = getfftchunk(signal, i, chunkspersecond, samplesnumber)
@test typeof(frequencies) <: Signal
@test length(getdata(frequencies)) > 1

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
export IISO,
       SLF,
       IWB

abstract FDTDScheme

immutable IISO <: FDTDScheme
	λ::Float64              # Courant number
	a::Float64              # scheme parameter a
	b::Float64              # scheme parameter b
	c::Float64              # scheme parameter c
	name::AbstractString    # name of scheme
	IISO() = new(sqrt(3/4), 1/6, 0., 0., "Interp. Isotropic (IISO)")
end

immutable SLF <: FDTDScheme
	λ::Float64              # Courant number
	a::Float64              # scheme parameter a
	b::Float64              # scheme parameter b
	c::Float64              # scheme parameter c
	name::AbstractString    # name of scheme
	SLF() = new(sqrt(1/3), 0., 0., 0., "Standard Leapfrog Scheme (SLF)")
end

immutable IWB <: FDTDScheme
	λ::Float64              # Courant number
	a::Float64              # scheme parameter a
	b::Float64              # scheme parameter b
	c::Float64              # scheme parameter c
	name::AbstractString    # name of scheme
	IWB() = new(1., 1/4, 1/16, 0., "Interp. Wideband (IWB)")
end

function Base.show(io::IO, scm::FDTDScheme)
	println(io, "FDTD Scheme    : $(scm.name)")
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
module FindBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools, Statistics

const SUITE = BenchmarkGroup()

########################
# findall              #
########################

g = addgroup!(SUITE, "findall")

const VEC_LENGTH = 1000

for (name, x) in (("50-50", Vector{Bool}(samerand(Bool, VEC_LENGTH))),
                  ("10-90", Vector{Bool}(samerand(VEC_LENGTH) .> .9)),
                  ("90-10", Vector{Bool}(samerand(VEC_LENGTH) .> .1)))
    g[string(typeof(x)), name] = @benchmarkable findall($x)

    bx = BitArray(x)
    g[string(typeof(bx)), name] = @benchmarkable findall($bx)
end

ispos(x) = x > 0

for T in (Bool, Int8, Int, UInt8, UInt, Float32, Float64)
    y = samerand(T, VEC_LENGTH)
    g["ispos", string(typeof(y))] = @benchmarkable findall($ispos, $y)

    # vary the proportion of returned elements (i.e size of returned array)
    ps = [0.5, 0.8, 0.95, 0.99]
    qs = quantile(y, ps)
    for (p, q) in zip(ps, qs)
        g["> q$p", string(typeof(y))] = @benchmarkable findall($(x -> x > q), $y)
    end
end

########################
# findnext/findprev    #
########################

gn = addgroup!(SUITE, "findnext")
gp = addgroup!(SUITE, "findprev")

const VEC_LENGTH = 1000

function perf_findnext(x)
    s = findfirst(x)
    while s !== nothing
        s = findnext(x, nextind(x, s))
    end
    s
end

function perf_findprev(x)
    s = findlast(x)
    while s !== nothing
        s = findprev(x, prevind(x, s))
    end
    s
end

function perf_findnext(pred, x)
    s = findfirst(pred, x)
    while s !== nothing
        s = findnext(pred, x, nextind(x, s))
    end
    s
end

function perf_findprev(pred, x)
    s = findlast(pred, x)
    while s !== nothing
        s = findprev(pred, x, prevind(x, s))
    end
    s
end

for (name, x) in (("50-50", samerand(Bool, VEC_LENGTH)),
                  ("10-90", samerand(VEC_LENGTH) .> .9),
                  ("90-10", samerand(VEC_LENGTH) .> .1))
    bx = BitArray(x)

    gn[string(typeof(x)), name] = @benchmarkable perf_findnext($x)
    gn[string(typeof(bx)), name] = @benchmarkable perf_findnext($bx)

    gp[string(typeof(x)), name] = @benchmarkable perf_findprev($x)
    gp[string(typeof(bx)), name] = @benchmarkable perf_findprev($bx)
end

for T in (Bool, Int8, Int, UInt8, UInt, Float32, Float64)
    y = samerand(T, VEC_LENGTH)
    gn["ispos", string(typeof(y))] = @benchmarkable perf_findnext($ispos, $y)
    gp["ispos", string(typeof(y))] = @benchmarkable perf_findprev($ispos, $y)
end

end


text:
module TestExamples1DBurgers

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")

@testset "Inviscid Burgers" begin
  @testset "elixir_burgers_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_basic.jl"),
      l2   = [2.967470209082194e-5],
      linf = [0.00016152468882624227])
  end

  @testset "elixir_burgers_linear_stability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_linear_stability.jl"),
      l2   = [0.5660569881106876],
      linf = [1.9352238038313998])
  end
end

end # module


text:
using Gtk.GConstants.GtkAlign: GTK_ALIGN_START, GTK_ALIGN_END, GTK_ALIGN_FILL

function change_red(col::CLim{T}, chan::CLim{T2}) where {T<:AbstractRGB, T2<:GrayLike}
    cmin = T(chan.min, green(col.min), blue(col.min))
    cmax = T(chan.max, green(col.max), blue(col.max))
    return CLim(cmin, cmax)
end
function change_green(col::CLim{T}, chan::CLim{T2}) where {T<:AbstractRGB, T2<:GrayLike}
    cmin = T(red(col.min), chan.min, blue(col.min))
    cmax = T(red(col.max), chan.max, blue(col.max))
    return CLim(cmin, cmax)
end
function change_blue(col::CLim{T}, chan::CLim{T2}) where {T<:AbstractRGB, T2<:GrayLike}
    cmin = T(red(col.min), green(col.min), chan.min)
    cmax = T(red(col.max), green(col.max), chan.max)
    return CLim(cmin, cmax)
end
change_red(col::Signal, chan::CLim) = change_red(value(col), chan)
change_green(col::Signal, chan::CLim) = change_green(value(col), chan)
change_blue(col::Signal, chan::CLim) = change_blue(value(col), chan)

function contrast_gui(enabled::Signal{Bool}, hists::Vector, clim::Signal{CLim{T}}) where {T<:AbstractRGB}
    @assert length(hists) == 3 #one signal per color channel
    chanlims = channel_clims(value(clim))
    rsig = Signal(chanlims[1])
    gsig = Signal(chanlims[2])
    bsig = Signal(chanlims[3])
    #make sure that changes to individual channels update the color clim signal and vice versa
    bindmap!(clim, x->change_red(clim, x), rsig, x->channel_clim(red, x); initial = false)
    bindmap!(clim, x->change_green(clim, x), gsig, x->channel_clim(green, x); initial = false)
    bindmap!(clim, x->change_blue(clim, x), bsig, x->channel_clim(blue, x); initial = false)
    names = ["Red Contrast"; "Green Contrast"; "Blue Contrast"]
    csigs = [rsig; gsig; bsig]
    cguis = []
    for i=1:length(hists)
        push!(cguis, contrast_gui(enabled, hists[i], csigs[i]; wname = names[i]))
    end
    return cguis
end

contrast_gui(enabled, hist::Vector, clim) = contrast_gui(enabled, hist[1], clim)

function contrast_gui(enabled::Signal{Bool}, hist::Signal, clim::Signal; wname="Contrast")
    vhist, vclim = value(hist), value(clim)
    T = eltype(vclim)
    Δ = T <: Integer ? T(1) : eps(T)
    rng = vhist.edges[1]
    cmin, cmax = vclim.min, vclim.max
    if !(cmin < cmax)
        cmin, cmax = first(rng), last(rng)
        if !(cmin < cmax)
            cmin, cmax = zero(cmin), one(cmax)
        end
    end
    smin = Signal(convert(eltype(rng), cmin))
    smax = Signal(convert(eltype(rng), cmax))
    cgui = contrast_gui_layout(smin, smax, rng; wname=wname)
    signal_connect(cgui["window"], :destroy) do widget
        push!(enabled, false)
    end
    updateclim = map(smin, smax) do cmin, cmax
        # if min/max is outside the current range, update the sliders
        adj = Gtk.Adjustment(widget(cgui["slider_min"]))
        rmin, rmax = Gtk.G_.lower(adj), Gtk.G_.upper(adj)
        if cmin < rmin || cmax > rmax || cmax-cmin < Δ
            # Also, don't cross the sliders
            bigmax = max(cmin,cmax,rmin,rmax)
            bigmin = min(cmin,cmax,rmin,rmax)
            thismax = min(typemax(T), max(cmin, cmax, rmax))
            thismin = max(typemin(T), min(cmin, cmax, rmin))
            rng = range(thismin, stop=thismax, length=255)
            cminT, cmaxT = T(min(cmin, cmax)), T(max(cmin, cmax))
            if cminT == cmaxT
                cminT = min(cminT, cminT-Δ)
                cmaxT = max(cmaxT, cmaxT+Δ)
            end
            mn, mx = minimum(rng), maximum(rng)
            cmin, cmax = clamp(cminT, mn, mx), clamp(cmaxT, mn, mx)
            push!(cgui["slider_min"], rng, cmin)
            push!(cgui["slider_max"], rng, cmax)
        end
        # Update the image contrast
        push!(clim, CLim(cmin, cmax))
    end
    # TODO: we might want to throttle this?
    redraw = draw(cgui["canvas"], hist) do cnvs, hst
        if get_gtk_property(cgui["window"], :visible, Bool) # protects against window destruction
            rng, cl = hst.edges[1], value(clim)
            mn, mx = minimum(rng), maximum(rng)
            push!(cgui["slider_min"], rng, clamp(cl.min, mn, mx))
            push!(cgui["slider_max"], rng, clamp(cl.max, mn, mx))
            drawhist(cnvs, hst)
        end
    end
    GtkReactive.gc_preserve(cgui["window"], (cgui, redraw, updateclim))
    cgui
end

function contrast_gui_layout(smin::Signal, smax::Signal, rng; wname="Contrast")
    win = Window(wname) |> (g = Grid())
    slmax = slider(rng; signal=smax)
    slmin = slider(rng; signal=smin)
    for sl in (slmax, slmin)
        set_gtk_property!(sl, :draw_value, false)
    end
    g[1,1] = widget(slmax)
    g[1,3] = widget(slmin)
    cnvs = canvas(UserUnit)
    g[1,2] = widget(cnvs)
    set_gtk_property!(cnvs, :expand, true)
    emax_w = Entry(; width_chars=5, hexpand=false, halign=GTK_ALIGN_END, valign=GTK_ALIGN_START)
    emin_w = Entry(; width_chars=5, hexpand=false, halign=GTK_ALIGN_END, valign=GTK_ALIGN_END)
    g[2,1] = emax_w
    g[2,3] = emin_w
    # By not specifying the range on the textbox, we let the user
    # enter something out-of-range, which can be handy in some
    # circumstances.
    emax = textbox(eltype(smax); widget=emax_w, signal=smax) # , range=rng)
    emin = textbox(eltype(smin); widget=emin_w, signal=smin) #, range=rng)

    Gtk.showall(win)
    Dict("window"=>win, "canvas"=>cnvs, "slider_min"=>slmin, "slider_max"=>slmax, "textbox_min"=>emin, "textbox_max"=>emax)
end

# We could use one of the plotting toolkits, but most are pretty slow
# to load and/or produce the first plot. So let's just do it manually.
@guarded function drawhist(canvas, hist)
    ctx = getgc(canvas)
    fill!(canvas, colorant"white")
    edges, counts = hist.edges[1], hist.weights
    xmin, xmax = first(edges), last(edges)
    cmax = maximum(counts)
    if cmax <= 0 || !(xmin < xmax)
        return nothing
    end
    set_coordinates(ctx, BoundingBox(xmin, xmax, log10(cmax+1), 0))
    set_source(ctx, colorant"black")
    move_to(ctx, xmax, 0)
    line_to(ctx, xmin, 0)
    for (i, c) in enumerate(counts)
        line_to(ctx, edges[i], log10(c+1))
        line_to(ctx, edges[i+1], log10(c+1))
    end
    line_to(ctx, xmax, 0)
    fill(ctx)
    nothing
end


text:
module Kaggler

include("common/common.jl"); reexport dim

end


text:
using Docile

config = [
    :output => [plain, html, helpdb]
]

Docile.build("Blending", config)