text:
#!/usr/bin/env julia
#=
   Author: Tim Sterne-Weiler, 3/16/2015
   e-mail: tim.sterne.weiler@utoronto.ca

   Calculate enrichment of chimeras over expression level background
   using multinomial draws of n=2 to calculate a second multinomial
   distribution of gene pairs which is used to compare observed frequency
   to expected where each pair is binomially distributed
   k occurences of chimera (i,j) is k[(i,j)]~ Bin( n, p[(i,j)] )
=#

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ## INITIALIZATION  ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

global head = "[aligater stats]:"
println(STDERR, "$head Loading Packages..")

using ArgParse
using StatsBase
using Distributions

function parse_cmd()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--back", "-b"
      help = "lig file to build background from"
      arg_type = ASCIIString
    "--fore", "-f"
      help = "lig file to evaluate"
      arg_type = ASCIIString
    "--gi"
      help = "index for geneid:geneid field in foreground"
      arg_type = Int
      default = 25
    "--ci"
      help = "index for counting field if applicable (non integer fields increment by 1 instead)"
      arg_type = Int
    "--bi"
      help = "index for background calibration (default is [end] of collapsed lig without any post)"
      arg_type = Int
      default = 19
    "--nd"
      help = "comma delmited set of strings to substitute for wildcard.  The ratio of each pair is output"
      arg_type = ASCIIString
    "--filt"
      help = "pattern filter [syntax= column:regex]"
      arg_type = ASCIIString
      default = "1:."
    "--vs"
      help = "extra variable string etc.., 'col:type,col:type', where type = [pcfdn]"
      arg_type = ASCIIString
      default = ""
    "--alpha"
      help = "maximum p-value for output, only applied to the numerator (first element) of --nd if supplied"
      arg_type = Float64
      default = 1.0
    "--nc"
      help = "normalization constants, applied to each of --nd as a comma separated list."
      arg_type = ASCIIString
    "--rpm"
      help = "function to combine the background expression values for RPM column"
      arg_type = ASCIIString
      default = "min"
  end
  return parse_args(s)
end

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ###  FUNCTIONS ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

include("dictext.jl")

# this accepts a --filt string argument in the form 'column:regexstring'
function parse_colfilt(toparse::ASCIIString)
  cap = match(r"(\d+)\:(.*)", toparse).captures
  @assert( length(cap) == 2 )
  col = parse(Int, cap[1])
  reg = Regex(cap[2])
  col, reg
end #--> (Integer, Regex)

#= make joint interaction distribution
     P(x,y) = [ x * y  when x != y
              |
              [   0    when x == y
=#
## Use this function to find the denominator for x != y
#  By dividing each x * y pair by 1-Sum(x == y)
function sumSamePairs(pset::Dict)
  sum = 1.0
  for i in keys(pset)
    sum -= 2(pset[i] ^ 2)
  end
  sum
end #--> Float64

# make a pseudoCount
function pseudoCnt(pset::Dict)
  pseudo = minimum( collect( values( pset ) ) )^2
  pseudo * 0.5
end #--> Float

# get probability of pairs
function jointInteracProb{K,V}(pset::Dict{K,V}, key1::K, key2::K, denom = 1.0, pseudo = 0.0)
  if(!haskey(pset, key1) || !haskey(pset, key2))
    return( pseudo )
  end
  2(pset[key1] * pset[key2]) / denom
end #--> Float

#= This function takes a foreground count set for gene pairs
   and a background probability set and derives a probability
   set for those pairs and computes binomial p-values for the
   foreground counts =#
function calculateBinomialStats{T <: Tuple, S <: AbstractString}(forecnt::Dict{T,Float64}, backprob::Dict{S,Float64})
  const n = sum( collect( values(forecnt) ) )
  const den = sumSamePairs(backprob)
  const pseudo = pseudoCnt(backprob)
  const bonfCor = length(keys(forecnt))

  foreprob = Dict{Tuple,Float64}()
  retval = Dict{Tuple{ASCIIString,ASCIIString},Any}()

  obsDen = 0.0

  for (k1,k2) in keys(forecnt)
    prob = jointInteracProb(backprob, k1, k2, den, pseudo)
    obsDen += prob
    foreprob[(k1,k2)] = prob
  end

  for (k1,k2) in keys(forecnt)
    prob = foreprob[(k1,k2)] / obsDen
    bin = Binomial(Int(n), prob)
    k = forecnt[(k1,k2)]
    pval = ccdf(bin, k-1)
    pval *= bonfCor
    pval = (pval > 1.0) ? one(pval) : pval
    expval = prob * n
    retval[(k1,k2)] = (k, expval, pval, k/expval)
  end
  retval
end #--> Dict{(String,String), Any}

# if the pairs --TODODOC
function loadInteractionFile(io::IOStream, gInd::Int, cntInd::Int, col, reg)
  cset = Dict{Tuple,Float64}()
  sizehint!(cset, 100000) # allocate 100K pairs
  for l::ASCIIString in eachline( io )
    s = split(chomp(l), '\t')
    if isa(col, Integer) && isa(reg, Regex) && length(s) >= col
      ismatch(reg, s[col]) || continue
    end
    # make sure we don't get out of bounds error
    (ismatch(r":", s[gInd])) || continue
    k1,k2 = sort( split(s[gInd], ':') )
    cnum = parse( s[cntInd], raise=false )
    dinc!(cset, (k1,k2), isa(cnum, Number) ? convert(Float64, cnum) : 1.0)
  end
  cset
end #--> Dict{Tuple,Float64}

# this does io and creates a cset structure of type Dict
function loadBackgroundFile(io::IOStream, ind::Int; keyType=AbstractString, valType=Float64)
  cset = Dict{keyType,valType}()
  sizehint!(cset, 30000) # approx number of genes
  for i::ASCIIString in eachline(io)
    s = split(chomp(i), '\t')
    #=genes = replace(s[ind], ":*NA:*", "") |> x->split(x, ':')
    genes = split(s[ind], ':')
    length(genes) > 0 && dinc!(cset, genes[1], one(valType)) =# # 04/15/15 deprecate for performance
    dinc!(cset, s[ind], one(valType))
  end
  cset
end #--> Dict{keyType,valType}

# This takes a cmd line string that contains pairs=> "column:type,column:type"
# where column is an int and type is a Char of [sfdn]
function varStatSummary( openwhat, refInd::Int, varstr::ASCIIString, col, reg )
  # helper functions:
  function parseVarStr( str::ASCIIString )
    sets = split(str, ',')
    ret  = Tuple{Int,Char}[]
    for i in sets
      vcol,typ = map(parse, split(i, ':'))
      tchar = string(typ)[1]
      @assert( ismatch( r"[psfdn]", string(tchar) ) )
      push!(ret, (vcol, tchar) )
    end
    ret
  end #--> Array{(Int,Char),1}

  function letterToType( letter::Char )
    letter == 'p' && return typeof(("",""))
    letter == 's' && return ASCIIString
    letter == 'f' && return Float64
    letter == 'd' && return Int64
    letter == 'n' && return Number
    Any
  end #--> Type{T}

  reOrder( t::Tuple ) = (length(t) <= 1 || t[1] <= t[2]) ? t : (t[2],t[1])
  toStrings( arr ) = map( x->convert(ASCIIString,x), arr )

  # varStatSummary code:
  varSet = parseVarStr( varstr )
  summary = Dict{Int,Dict}()
  open( openwhat , "r" ) do fh
    for l::ASCIIString in eachline(fh)
      s = split(chomp(l), "\t")
      # follow the same filter rules as interactionFile:
      if isa(col, Integer) && isa(reg, Regex) && length(s) >= col
        ismatch(reg, s[col]) || continue
      end
      # get reference key.
      refKey = sort( split(s[refInd], ':') )
      k = tuple(refKey...)[1:2] # main refKey k
      # iterate through column & type pairs
      for (i,c) in varSet
        cType = letterToType( c ) # convert char to type
        parSi = cType <: Number ? parse(cType, string(s[i]) ) : string( s[i] )
        cVal  = cType <: Tuple ? begin # if tuple check if properly ordered
                                    a,b = split(parSi, ':') |> toStrings
                                    reOrder( (a,b) )
                                 end : convert(cType, parSi)
                                  #otherwise initial parse was fine
        @assert( isa(cVal, cType) )
        if !haskey(summary, i)
          vType = cType <: ASCIIString || cType <: Tuple ? cType : Array{cType,1}
          summary[i] = Dict{typeof(k), vType}()
        end
        if cType <: ASCIIString || cType <: Tuple # set to single value unless already set
          !haskey( summary[i], k ) && ( summary[i][k] = cVal )
        else # numeric, so push the current element
          dush!( summary[i], k, cVal, arrType=cType )
        end
      end
    end
  end
  summary
end #--> Dict{Int,Dict{Tuple, ? }}

# central function for processing a foreground background filename pair
#=function loadFilesAndCalculate(forefile::ASCIIString, backfile::ASCIIString, pargs)
  # fetch options.
  const backInd = pargs["bi"]
  const geneInd = pargs["gi"]
  const cntInd  = pargs["ci"] == nothing ? geneInd : pargs["ci"]

  c,r = parse_colfilt(pargs["filt"])

  print(STDERR, "Loading Background $backfile..\n")
  backhndl = open(backfile, "r")
  cset = loadBackgroundFile(backhndl, backInd, valType=Float64)
  close(backhndl)

  print(STDERR, "Processing probability distribution..\n")
  # process the probability distribution
  pset = deepcopy(cset) # multinomial probabilities
  dnorm!( pset ) # normalize dict

  # if --rpm then print abundance levels also
  @assert( ismatch(r"^[A-Z|a-z|0-9]+$", pargs["rpm"]) )
  rpmFunc = eval(parse(pargs["rpm"]))
  @assert( typeof(rpmFunc) <: Function )
  rset = deepcopy(pset)
  dnorm!(rset, 1/1000000) # reads per M

  # iterate through foreground
  print(STDERR, "Loading Foreground $forefile..\n")
  forehndl = open(forefile, "r")
  forecnt  = loadInteractionFile(forehndl, geneInd, cntInd, c, r)
  close(forehndl)

  print(STDERR, "Calculating binomial stats..\n")
  # calculate p-values and return hash
  calculateBinomialStats( forecnt, pset, rset, rpmFunc ) # return stat hash
end =#  ##Deprecated 7/16/2015
#--> Dict{(String,String), Any}

# final print function for stat array
function printStats( io, statarr::Array, rpmarr::Array; normarr=[1,1], alpha=1.0, vardict=Dict(), rpmFunc=min )
  aHsh = statarr[1]
  aRpm = rpmarr[1]
  @assert( isa(aHsh, Dict) )
  @assert( isa(aRpm, Dict) )

  # internal print function for variable keys
  function printVardict( io, dict, refkey )
    for col in keys(dict)
      if haskey( dict[col], refkey )
        val = dict[col][refkey]
        if isa( val, Array )
          asfloat = typeof(val) <: Array{Int64,1} ? convert(Array{Float64,1}, val) : val
          aMed = median( asfloat )
          aMad = mad( asfloat )
          @printf( io, "\t%d\t%.2f±%.2f", col, aMed, aMad)
        elseif isa( val, Tuple )
          @printf( io, "\t%d\t%s", col, join( val, ',' ) )
        else
          @printf( io, "\t%d\t%s", col, val )
        end
      else
        @printf( io, "\t%d\tNA", col )
      end
    end
    @printf( io, "\n" )
  end #--> nothing

  if length(statarr) == 1
    # just a single set print output normally
    for (k1,k2) in keys( aHsh )
      k, kExp, pval, obsExp = aHsh[(k1,k2)]
      aPval <= alpha || continue  #filter p-val
      @printf( io, "%s,%s\t%d\t%.1f\t%.2e\t%.2e", k1, k2, k, k/kExp, pval, obsExp )
      printVardict( io, vardict, (k1,k2) )
    end
  elseif length(statarr) == 2
    bHsh = statarr[2]
    bRpm = rpmarr[2]
    @assert( isa(bHsh, Dict) )
    @assert( isa(bRpm, Dict) )
    ksone  = collect( keys( aHsh ) )
    kstwo  = collect( keys( bHsh ) )
    both = union( ksone, kstwo )
    # function to check existance of and reduce rpm values.
    function apply_rpmfunc{K <: AbstractString,V <: Float64}( hsh1::Dict{K,V}, hsh2::Dict{K,V}, key1::K, key2::K, func::Function)
      begin
        val1 = get(hsh1, key1, zero(V))
        val2 = get(hsh2, key2, zero(V))
        func(val1, val2)
      end
    end #--> V
    # iterat ehrough shared set.
    for (k1,k2) in both
      aK, aExp, aPval, aObsExp = haskey(aHsh,(k1,k2)) ? aHsh[(k1,k2)] : (0.5, 0.5, 1.0, 1.0)
      bK, bExp, bPval, bObsExp = haskey(bHsh,(k1,k2)) ? bHsh[(k1,k2)] : (0.5, 0.5, 1.0, 1.0)
      aRpmVal = apply_rpmfunc( aRpm, aRpm, k1, k2, rpmFunc )
      bRpmVal = apply_rpmfunc( bRpm, bRpm, k1, k2, rpmFunc )
      a_b = (aK / normarr[1]) / (bK / normarr[2])
      exp_a_b = (aExp / normarr[1]) / (bExp / normarr[2])
      full_a_b = a_b / exp_a_b
      aPval <= alpha || continue # filter if p-value is above cutoff
      @printf( io, "%s,%s\t%.1f\t%.1f\t%.1f\t%d\t%.1f\t%d\t%.1f\t%.2e\t%.2e\t%.1f\t%.1f", k1, k2, full_a_b, a_b, exp_a_b, aK, aK/aExp, bK, bK/bExp, aPval, bPval, aRpmVal, bRpmVal )
      printVardict( io, vardict, (k1,k2) )
    end #endfor
  end #endif
end #--> nothing

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ### MAIN  ### ### ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #
function main()
  pargs = parse_cmd()

  const backInd = pargs["bi"]
  const geneInd = pargs["gi"]
  const cntInd  = pargs["ci"] == nothing ? geneInd : pargs["ci"]

  c,r = parse_colfilt(pargs["filt"])

  # check options.
  ndArray = pargs["nd"] == nothing ? [""] : split(pargs["nd"], ",")
  normArray = pargs["nc"] == nothing ? [1,1] : map(parse, split(pargs["nc"], ","))

  @assert( pargs["fore"] != nothing, "$head --fore and --back must be given!")
  @assert( pargs["back"] != nothing, "$head --fore and --back must be given!")

  ndParsed = ASCIIString[]
  stats = Dict[]
  rpms  = Dict[]

  varstat = nothing  #scope

  # parse --rpm flag to Function
  @assert( ismatch(r"^[A-Z|a-z|0-9]+$", pargs["rpm"]) )
  rpmFunc = eval(parse(pargs["rpm"]))
  @assert( typeof(rpmFunc) <: Function, "typeof(--rpm) must be Julia Function!" )

  # iterate through delimited names list
  for (ndIter,nd) in enumerate(ndArray)

    # if wilcard present, substitute with current pattern
    backfile = replace(pargs["back"], "%", nd)
    forefile = replace(pargs["fore"], "%", nd)

    push!(ndParsed, forefile)

    # calculate p-values
    #stathsh = loadFilesAndCalculate( forefile, backfile, pargs ) ##Deprecated 7/16/2015

    print(STDERR, "Loading Background $backfile..\n")
    backhndl = open(backfile, "r")
    cset = loadBackgroundFile(backhndl, backInd, valType=Float64)
    close(backhndl)

    print(STDERR, "Processing probability distribution..\n")
    # process the probability distribution
    pset = deepcopy(cset) # multinomial probabilities
    dnorm!( pset ) # normalize dict

    # if --rpm then print abundance levels also
    rset = deepcopy(pset)
    dnorm!(rset, 1/1000000) # reads per M

    # iterate through foreground
    print(STDERR, "Loading Foreground $forefile..\n")
    forehndl = open(forefile, "r")
    forecnt  = loadInteractionFile(forehndl, geneInd, cntInd, c, r)
    close(forehndl)

    print(STDERR, "Calculating binomial stats..\n")
    # calculate p-values and return hash
    stathsh = calculateBinomialStats( forecnt, pset ) # return stat hash

    push!(stats, stathsh)
    push!(rpms, rset)
  end

  if length(pargs["vs"]) > 0
    println(STDERR, "Loading --vs variables")
    c,r = parse_colfilt(pargs["filt"])
    varstat = varStatSummary( `cat $ndParsed` , pargs["gi"], pargs["vs"], c, r )
  end

  @assert(0 < length(stats) <= 2, "--nd does not contain properly formatted arguments!")
  printStats( STDOUT, stats, rpms, alpha=pargs["alpha"], normarr=normArray, vardict=varstat, rpmFunc=rpmFunc )

end
#####################################################
# main execution here
main()
# eof

max_forks_repo_name: timbitz/aligater


text:
######################################################################
# ergodic_manager/r3.jl
#
# handles stuff needed for ergodicity
######################################################################

export ErgodicManagerR3

"""
`ErgodicManagerR3(d::Domain, phi::Matrix{Float64}, K::Int)`
"""
type ErgodicManagerR3 <: ErgodicManager
	domain::Domain				# spatial domain
	K::Int						# number of Fourier coefficients
	phi::Array{Float64,3}		# spatial distribution
	phik::Array{Float64,3}		# distribution's Fourier coefficients

	# constant regardless of phi (depend on k1,k2)
	Lambda::Array{Float64,3}
	hk::Array{Float64,3}

	# to speed up computation
	kpixl::Matrix{Float64}
	kpiyl::Matrix{Float64}
	kpizl::Matrix{Float64}

	function ErgodicManagerR3(d::Domain, phi::Array{Float64,3}, K::Int=5)
		em = new()
		em.domain = deepcopy(d)
		em.K = K
		em.hk = zeros(K+1,K+1,K+1)
		em.phi = deepcopy(phi)
		em.phik = zeros(K+1,K+1,K+1)
		em.Lambda = zeros(K+1,K+1,K+1)
		em.kpixl = zeros(K+1, d.cells[1])
		em.kpiyl = zeros(K+1, d.cells[2])
		em.kpizl = zeros(K+1, d.cells[3])

		Lambda!(em)
		kpixl!(em)
		hk!(em)
		decompose!(em)
		return em
	end
end

# fills each entry Lambda[k1,k2] in the Lambda matrix
function Lambda!(em::ErgodicManagerR3)
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				den = (1.0 + k1*k1 + k2*k2 + k3*k3) ^ 2.0
				em.Lambda[k1+1, k2+1, k3+1] = 1.0 / den
			end
		end
	end
end

# TODO: check if the row/col ordering of this is julia-efficient
# TODO: also, I could probably speed up by iterating over k first?
function kpixl!(em::ErgodicManagerR3)
	Lx = em.domain.lengths[1]
	xmin = x_min(em)
	for xi = 1:x_cells(em)
		# TODO: seems like I just cancel xmin out. Ignore it maybe?
		x = xmin + (xi-0.5)*x_size(em)
		for k = 0:em.K
			em.kpixl[k+1,xi] = cos(k*pi*(x-xmin) / Lx)
		end
	end

	Ly = em.domain.lengths[2]
	ymin = y_min(em)
	for yi = 1:y_cells(em)
		y = ymin + (yi-0.5)*y_size(em)
		for k = 0:em.K
			em.kpiyl[k+1,yi] = cos(k*pi*(y-ymin) / Ly)
		end
	end

	Lz = em.domain.lengths[3]
	zmin = z_min(em)
	for zi = 1:z_cells(em)
		z = zmin + (zi-0.5)*z_size(em)
		for k = 0:em.K
			em.kpizl[k+1,zi] = cos(k*pi*(z-zmin) / Lz)
		end
	end
end

# generates the hk coefficients for the ergodic manager
# these coefficients only need to be computed once
function hk!(em::ErgodicManagerR3)
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				em.hk[k1+1,k2+1,k3+1] = hk_ij(em, k1, k2, k3)
			end
		end
	end
end

# computes the coefficients for a specific value of k1 and k2
# called by hk!
# TODO: this is wrecked by domain
function hk_ij(em::ErgodicManagerR3, k1::Int, k2::Int, k3::Int)
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		cx2 = cx * cx
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			cy2 = cy * cy
			for zi = 1:z_cells(em)
				cz = em.kpizl[k3+1,zi]
				cz2 = cz * cz
				val += cx2 * cy2 * cz2 * em.domain.cell_size
			end
		end
	end

	return sqrt(val)
end

######################################################################
# Computing Fourier coefficients
######################################################################
# update the Fourier coefficients based on some distribution
# Here, I assume it is discrete, but I should change this...
# TODO: maybe I should do some bounds checking?
function decompose!(em::ErgodicManagerR3, d::Array{Float64,3})
	for k1 = 0:em.K
		for k2 = 0:em.K
			for k3 = 0:em.K
				em.phik[k1+1,k2+1,k3+1] = phi_ij(em, k1, k2, k3, d)
			end
		end
	end
	em.phi = d
end

# iterate over the state space
function phi_ij(em::ErgodicManagerR3, k1::Int, k2::Int, k3::Int, d::Array{Float64,3})
	val = 0.0
	for xi = 1:x_cells(em)
		cx = em.kpixl[k1+1,xi]
		for yi = 1:y_cells(em)
			cy = em.kpiyl[k2+1,yi]
			for zi = 1:z_cells(em)
				cz = em.kpizl[k3+1,zi]
				val += d[xi,yi,zi] * cx * cy * cz * em.domain.cell_size
			end
		end
	end
	return val / em.hk[k1+1,k2+1,k3+1]
end

function decompose(em::ErgodicManagerR3, traj::VVF)
	N = length(traj)-1

	# Array to hold trajectory's Fourier coefficients
	ck = zeros(em.K+1, em.K+1, em.K+1)

	# length of each dimension
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	Lz = em.domain.lengths[3]

	# minimum values in each dimension
	xmin = x_min(em)
	ymin = y_min(em)
	zmin = z_min(em)

	for k1 = 0:em.K
		kpiL1 = k1 * pi / Lx
		for k2 = 0:em.K
			kpiL2 = k2 * pi / Ly
			for k3 = 0:em.K
				kpiL3 = k3 * pi / Lz
				hk = em.hk[k1+1, k2+1, k3+1]
				fk_sum = 0.0
				# now loop over time
				for n = 0:N
					xn = traj[n + 1]
					c1 = cos(kpiL1 * (xn[1] - xmin))
					c2 = cos(kpiL2 * (xn[2] - ymin))
					c3 = cos(kpiL3 * (xn[3] - zmin))
					fk_sum += c1*c2*c3
				end
				ck[k1+1, k2+1, k3+1] = fk_sum / (hk * (N+1))
			end
		end
	end
	return ck
end

ext: jl


text:
@testset "Select" begin
    @testset "Dense" begin
        #dense first
        X = GBMatrix(rand(-10.0:0.05:10.0, 10, 10))
        test_frule(select, diag, X)
        test_rrule(select, diag, X)
        test_frule(select, offdiag, X)
        test_rrule(select, offdiag, X)
        test_frule(select, tril, X)
        test_rrule(select, tril, X)
        test_frule(select, triu, X)
        test_rrule(select, triu, X)
        test_frule(select, nonzeros, X)
        test_rrule(select, nonzeros, X)
        test_frule(select, >, X, 0.)
        test_rrule(select, >, X, 0.)
        test_frule(select, <=, X, 0.)
        test_rrule(select, <=, X, 0.)
    end
    @testset "Sparse" begin
        X = GBMatrix(sprand(50, 50, 0.15))
        test_frule(select, diag, X)
        test_rrule(select, diag, X)
        test_frule(select, offdiag, X)
        test_rrule(select, offdiag, X)
        test_frule(select, tril, X)
        test_rrule(select, tril, X)
        test_frule(select, triu, X)
        test_rrule(select, triu, X)
        test_frule(select, nonzeros, X)
        test_rrule(select, nonzeros, X)
        test_frule(select, >, X, 0.)
        test_rrule(select, >, X, 0.)
        test_frule(select, <=, X, 0.)
        test_rrule(select, <=, X, 0.)
    end
end


text:
function a(input)
    input =  string(input, input[1])
    checkSum = 0

    for i in 1:length(input)-1
        if input[i] == input[i+1]
            checkSum += parse(Int32, input[i])
        end
    end

    return checkSum
end

function b(input)

end

let
    input = "237369991482346124663395286354672985457326865748533412179778188397835279584149971999798512279429268727171755461418974558538246429986747532417846157526523238931351898548279549456694488433438982744782258279173323381571985454236569393975735715331438256795579514159946537868358735936832487422938678194757687698143224139243151222475131337135843793611742383267186158665726927967655583875485515512626142935357421852953775733748941926983377725386196187486131337458574829848723711355929684625223564489485597564768317432893836629255273452776232319265422533449549956244791565573727762687439221862632722277129613329167189874939414298584616496839223239197277563641853746193232543222813298195169345186499866147586559781523834595683496151581546829112745533347796213673814995849156321674379644323159259131925444961296821167483628812395391533572555624159939279125341335147234653572977345582135728994395631685618135563662689854691976843435785879952751266627645653981281891643823717528757341136747881518611439246877373935758151119185587921332175189332436522732144278613486716525897262879287772969529445511736924962777262394961547579248731343245241963914775991292177151554446695134653596633433171866618541957233463548142173235821168156636824233487983766612338498874251672993917446366865832618475491341253973267556113323245113845148121546526396995991171739837147479978645166417988918289287844384513974369397974378819848552153961651881528134624869454563488858625261356763562723261767873542683796675797124322382732437235544965647934514871672522777378931524994784845817584793564974285139867972185887185987353468488155283698464226415951583138352839943621294117262483559867661596299753986347244786339543174594266422815794658477629829383461829261994591318851587963554829459353892825847978971823347219468516784857348649693185172199398234123745415271222891161175788713733444497592853221743138324235934216658323717267715318744537689459113188549896737581637879552568829548365738314593851221113932919767844137362623398623853789938824592"

    println(a(input))
end


text:
export some

import Base: show

"""
    some()

Creates a some operator, which filters out `nothing` items by the source Observable by emitting only
those that not equal to `nothing`.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream `<: Subscribable{Union{L, Nothing}}`

# Examples
```jldoctest
using Rocket

source = from([ 1, nothing, 3 ])
subscribe!(source |> some(), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 3
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`max`](@ref), [`min`](@ref), [`logger`](@ref)
"""
some() = SomeOperator()

struct SomeOperator <: InferableOperator end

function on_call!(::Type{Union{L, Nothing}}, ::Type{L}, ::SomeOperator, source) where L
    return proxy(L, source, SomeProxy())
end

operator_right(::SomeOperator, ::Type{Union{L, Nothing}}) where L = L
operator_right(::SomeOperator, ::Type{L})                 where L = error("some() operator can operate on streams with data type '<: Union{Nothing, $L}', but '<: $L' was found.")

struct SomeProxy <: ActorProxy end

actor_proxy!(::Type{L}, proxy::SomeProxy, actor::A) where { L, A } = SomeActor{L, A}(actor)

struct SomeActor{L, A} <: Actor{Union{L, Nothing}}
    actor :: A
end

function on_next!(actor::SomeActor{L}, data::Union{L, Nothing}) where L
    if data !== nothing
        next!(actor.actor, data)
    end
end

on_error!(actor::SomeActor, err) = error!(actor.actor, err)
on_complete!(actor::SomeActor)   = complete!(actor.actor)

Base.show(io::IO, ::SomeOperator)         = print(io, "SomeOperator()")
Base.show(io::IO, ::SomeProxy)            = print(io, "SomeProxy()")
Base.show(io::IO, ::SomeActor{L}) where L = print(io, "SomeActor($L)")


text:
module EllipticalSliceSampling

using AbstractMCMC: AbstractMCMC
using ArrayInterfaceCore: ArrayInterfaceCore
using Distributions: Distributions

using Random: Random
using Statistics: Statistics

export ESSModel, ESS

# reexports
using AbstractMCMC: sample, MCMCThreads, MCMCDistributed, MCMCSerial
export sample, MCMCThreads, MCMCDistributed, MCMCSerial

include("abstractmcmc.jl")
include("model.jl")
include("interface.jl")

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# mutable struct of FreeMotion is defined in Modia3D/Composition/object3d.jl,
# since Object3D references FreeMotion and FreeMotion references Object3D

function driveJoint!(joint::FreeMotion)
   joint.isDriven = true
   joint.r.numericType  = ModiaMath.WR
   joint.v.numericType  = ModiaMath.WR
   joint.a.numericType  = ModiaMath.WR
   joint.q.numericType    = ModiaMath.WR
   joint.derq.numericType = ModiaMath.WR
   joint.w.numericType  = ModiaMath.WR
   joint.z.numericType  = ModiaMath.WR
end

function computeKinematics!(joint::FreeMotion, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D = obj.parent
   absolute::Bool   = parent.r_abs ≡ ModiaMath.ZeroVector3D && parent.R_abs ≡ ModiaMath.NullRotation

   obj.r_rel = joint.r.value
   obj.R_rel = ModiaMath.from_q( joint.q.value )

   if absolute
      obj.r_abs = obj.r_rel
      obj.R_abs = obj.R_rel
   else
      obj.r_abs = parent.r_abs + parent.R_abs'*obj.r_rel
      obj.R_abs = obj.R_rel*parent.R_abs
   end

   q = joint.q.value
   joint.residue_q.value = dot(q, q) - 1.0

   if analysis == ModiaMath.DynamicAnalysis
      dynamics::Object3Ddynamics       = obj.dynamics
      parentDynamics::Object3Ddynamics = parent.dynamics

      joint.residue_w.value = joint.w.value - 2.0*( @SMatrix( [ q[4]  q[3] -q[2] -q[1];
                                                               -q[3]  q[4]  q[1] -q[2];
                                                                q[2] -q[1]  q[4] -q[3]] )*joint.derq.value)

      if absolute
         dynamics.v0 = joint.v.value
         dynamics.a0 = joint.a.value
         dynamics.w  = joint.w.value
         dynamics.z  = joint.z.value
      else
         dynamics.v0 = parentDynamics.v0 + parent.R_abs'*(joint.v.value + cross(parentDynamics.w, obj.r_rel))
         dynamics.a0 = parentDynamics.a0 + parent.R_abs'*(joint.a.value + cross(parentDynamics.z, obj.r_rel) +
                                                          cross(parentDynamics.w, cross(parentDynamics.w, obj.r_rel)))
         dynamics.w  = parent.R_abs*parentDynamics.w + joint.w.value
         dynamics.z  = parent.R_abs*parentDynamics.z + joint.z.value
      end
   end
   return nothing
end

function computeForceTorqueAndResidue!(joint::FreeMotion, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   joint.residue_f.value = obj.dynamics.f
   joint.residue_t.value = obj.dynamics.t
   return nothing
end

# Sets the relative position vector
function set_r!(obj::Object3D, r::ModiaMath.Vector3D)::NOTHING
   if isNotFree(obj)
      error("Not allowed to call set_r!(..) on Object3D ", ModiaMath.fullName(obj), " because not freely moving")
   end
   obj.joint.r.value = r
   obj.r_rel         = r
   return nothing
end
set_r!(obj::Object3D, r::AbstractVector)::NOTHING = set_r!(obj, ModiaMath.Vector3D(r))

# Sets the absolute position vector
function set_r_abs!(obj::Object3D, r::ModiaMath.Vector3D)::NOTHING
   if isNotFree(obj)
      error("Not allowed to call set_r!(..) on Object3D ", ModiaMath.fullName(obj), " because not freely moving")
   end
   obj.joint.r.value = r
   obj.r_abs         = r
   return nothing
end
set_r_abs!(obj::Object3D, r::AbstractVector)::NOTHING = set_r_abs!(obj, ModiaMath.Vector3D(r))

# Sets the relative quaternion
function set_q!(obj::Object3D, q::ModiaMath.Quaternion)::NOTHING
   if isNotFree(obj)
      error("Not allowed to call set_q!(..) on Object3D ", ModiaMath.fullName(obj), " because not freely moving")
   end
   obj.joint.q.value = q
   obj.R_rel         = ModiaMath.from_q(q)
   return nothing
end
set_q!(obj::Object3D, q::AbstractVector)::NOTHING = set_q!(obj, ModiaMath.Quaternion(q))


text:
"""
    index(a, x)

Locate the leftmost value exactly equal to `x` in `a`.
"""
function index(a, x)
    i = bisect_left(a, x)
    return i != length(a) + 1 && a[i] == x ? i : nothing
end

"""
    find_lt(a, x)

Find rightmost value less than `x` in `a`.
"""
function find_lt(a, x)
    i = bisect_left(a, x)
    return i > 1 ? a[i-1] : nothing
end

"""
    find_le(a, x)

Find rightmost value less than or equal to `x` in `a`.
"""
function find_le(a, x)
    i = bisect_right(a, x)
    return i > 1 ? a[i-1] : nothing
end

"""
    find_gt(a, x)

Find leftmost value greater than `x` in `a`.
"""
function find_gt(a, x)
    i = bisect_right(a, x)
    return i != length(a) + 1 ? a[i] : nothing
end

"""
    find_ge(a, x)

Find leftmost item greater than or equal to `x` in `a`.
"""
function find_ge(a, x)
    i = bisect_left(a, x)
    return i != length(a) + 1 ? a[i] : nothing
end


text:
# Deal with atoms

export
    inferelement,
    vdw_radius,
    hard_sphere_radius,
    areatomssame,
    findcommonatoms,
    atommap,
    atomid,
    calphaindices,
    atomcoords,
    atomcoords!,
    mod_atom_info,
    ismodulator

"""
Determines the element of an atom from the atom name.
Returns the element as a string.
"""
function inferelement(atom_name::AbstractString)
    if occursin("C", atom_name)
        element = "C"
    elseif occursin("N", atom_name)
        element = "N"
    elseif occursin("O", atom_name)
        element = "O"
    elseif occursin("S", atom_name)
        element = "S"
    elseif occursin("H", atom_name)
        element = "H"
    elseif occursin("P", atom_name)
        element = "P"
    else
        println("Could not find element for atom with name ", atom_name)
        element = "-"
    end
    return element
end

"Dictionary of atomic VDW radii in Angstroms."
const vdw_radius = Dict{String, Float64}(
    "C"=> 1.70,
    "N"=> 1.55,
    "O"=> 1.52,
    "S"=> 1.80,
    "H"=> 1.20,
    "P"=> 1.80,
    "-"=> 1.70, # Say unknown elements are like carbon
)

"Dictionary of atomic hard sphere radii in Angstroms."
const hard_sphere_radius = Dict{String, Float64}(
    "C"=> 0.67,
    "N"=> 0.56,
    "O"=> 0.48,
    "S"=> 0.88,
    "H"=> 0.53,
    "P"=> 0.98,
    "-"=> 0.67, # Say unknown elements are like carbon
)

"""
Determines if two atoms are the same.
Atoms are the same if they have the same `atom_name`, `chain_id`, `res_n` and `element`.
`res_name` is not considered.
"""
function areatomssame(atom_one::Atom, atom_two::Atom)
    if atom_one.atom_name == atom_two.atom_name && atom_one.chain_id == atom_two.chain_id && atom_one.res_n == atom_two.res_n && atom_one.element == atom_two.element
        answer = true
    else
        answer = false
    end
    return answer
end

"""
Checks two atom arrays and finds common atoms.
Returns copies of input atom arrays with only common atoms remaining.
The returned arrays are therefore of equal length.
"""
function findcommonatoms(atoms_one::Array{Atom,1}, atoms_two::Array{Atom,1})
    common_atoms_one = Atom[]
    common_atoms_two = Atom[]
    for atom_one in atoms_one
        for atom_two in atoms_two
            if areatomssame(atom_one, atom_two)
                push!(common_atoms_one, atom_one)
                push!(common_atoms_two, atom_two)
                break
            end
        end
    end
    return common_atoms_one, common_atoms_two
end

"""
Finds the mapping of the atoms in the first list to those in the second.
Returns an array of length the same as the second list where the number
is the atom index in the first list.
Atoms missing from the first list are given 0.
"""
function atommap(atoms_one::Array{Atom,1}, atoms_two::Array{Atom,1})
    atom_map = Int[]
    for atom_two in atoms_two
        found = false
        for (i, atom_one) in enumerate(atoms_one)
            if areatomssame(atom_one, atom_two)
                push!(atom_map, i)
                found = true
                break
            end
        end
        if !found
            push!(atom_map, 0)
        end
    end
    @assert length(atom_map) == length(atoms_two)
    return atom_map
end

"A string describing an `Atom` in the form `atom_name`/`res_name`/`res_n`/`chain_id`."
atomid(atom::Atom) = "$(atom.atom_name)/$(atom.res_name)/$(atom.res_n)/$(atom.chain_id)"

"Return as an array the indices in an array of atoms that correspond to C-alpha atoms."
function calphaindices(atoms::Array{Atom,1})
    ca_inds = Int[]
    for (i, atom) in enumerate(atoms)
        if atom.atom_name == "CA"
            push!(ca_inds, i)
        end
    end
    return ca_inds
end

"Get or set the list of coordinates for a list of atoms."
function atomcoords(atoms::Array{Atom,1})
    coords = zeros(3, length(atoms))
    for (i, atom) in enumerate(atoms)
        coords[:,i] = atom.coords
    end
    return coords
end

function atomcoords(atoms::Array{Atom,1}, coords::Array{Float64})
    new_atoms = deepcopy(atoms)
    atomcoords!(new_atoms, coords)
    return new_atoms
end

"Set the coordinates of a list of atoms from a list of coordinates."
function atomcoords!(atoms::Array{Atom,1}, coords::Array{Float64})
    n_atoms = length(atoms)
    @assert n_atoms == size(coords, 2) "The number of atoms is not the same as the number of coordinates"
    for (i, atom) in enumerate(atoms)
        atom.coords = coords[:,i]
    end
end

"Dictionary of atomic information for fake modulator atoms."
const mod_atom_info = Dict{String, Any}(
    "atom_name"=> "X",
    "res_name" => "MOD",
    "chain_id" => 'X',
    "res_n"    => 0,
    "element"  => "X",
)

"Determines whether an `Atom` is a fake modulator atom."
function ismodulator(atom::Atom)
    mod_atom = Atom(mod_atom_info["atom_name"], mod_atom_info["res_name"], mod_atom_info["chain_id"], mod_atom_info["res_n"], [0.0, 0.0, 0.0], mod_atom_info["element"])
    if areatomssame(atom, mod_atom)
        answer = true
    else
        answer = false
    end
    return answer
end


text:
avxt() = false

macro threaded(code)
    return esc(:(
        if $(@__MODULE__).avxt()
            @avxt($code)
        else
            @threads($code)
        end
    ))
end

# multithreaded kernels
function copy_threaded(C, A)
    @assert length(C) == length(A)
    @threaded for i in eachindex(C, A)
        @inbounds C[i] = A[i]
    end
    return nothing
end

function scale_threaded(B, C, s)
    @assert length(C) == length(B)
    @threaded for i in eachindex(C)
        @inbounds B[i] = s * C[i]
    end
    return nothing
end

function add_threaded(C, A, B)
    @assert length(C) == length(B) == length(A)
    @threaded for i in eachindex(C)
        @inbounds C[i] = A[i] + B[i]
    end
    return nothing
end

function triad_threaded(A, B, C, s)
    @assert length(C) == length(B) == length(A)
    @threaded for i in eachindex(C)
        @inbounds A[i] = B[i] + s * C[i]
    end
    return nothing
end

# kernels
function copy(C, A)
    @assert length(C) == length(A)
    for i in eachindex(C, A)
        @inbounds C[i] = A[i]
    end
    return nothing
end

function scale(B, C, s)
    @assert length(C) == length(B)
    for i in eachindex(C)
        @inbounds B[i] = s * C[i]
    end
    return nothing
end

function add(C, A, B)
    @assert length(C) == length(B) == length(A)
    for i in eachindex(C)
        @inbounds C[i] = A[i] + B[i]
    end
    return nothing
end

function triad(A, B, C, s)
    @assert length(C) == length(B) == length(A)
    for i in eachindex(C)
        @inbounds A[i] = B[i] + s * C[i]
    end
    return nothing
end


text:
mutable struct EmpiricalMean{T<:Real,V<:AbstractVector{<:Real}} <: PriorMean{T}
    C::V
    opt::Optimizer
end

"""
**EmpiricalMean**
```julia`
function EmpiricalMean(c::V=1.0;opt::Optimizer=Adam(α=0.01)) where {V<:AbstractVector{<:Real}}
```
Construct a constant mean with values `c`
Optionally give an optimizer `opt` (`Adam(α=0.01)` by default)
"""
function EmpiricalMean(c::V=1.0;opt::Optimizer=Adam(α=0.01)) where {V<:AbstractVector{<:Real}}
    EmpiricalMean{eltype(c),V}(c,opt)
end

function update!(μ::EmpiricalMean{T},grad::AbstractVector{T}) where {T<:Real}
    μ.C .+= update!(μ.opt,grad)
end

function (μ::EmpiricalMean{T})(x::AbstractMatrix) where {T<:Real}
    @assert size(x,1)==length(μ)
    return μ.C
end


text:
function eftAdd{T<:StdFloat}(a::T, b::T)
  hi = a + b
  t = hi - a
  lo = (a - (hi - t)) + (b - t)
  hi,lo
end

function eftAddGTE{T<:StdFloat}(a::T, b::T)
  hi = a + b
  lo = b - (hi - a)
  hi,lo
end

function eftAddLTE{T<:StdFloat}(a::T, b::T)
  hi = a + b
  lo = a - (hi - b)
  hi,lo
end

function eftSub{T<:StdFloat}(a::T, b::T)
  hi = a - b
  t = hi - a
  lo = (a - (hi - t)) - (b + t)
  hi,lo
end

function eftSubGTE{T<:StdFloat}(a::T, b::T)
  hi = a - b
  lo = (a - hi) - b
  hi,lo
end

function eftMul{T<:StdFloat}(a::T, b::T)
    hi = a * b
    lo = fma(a, b, -hi)
    hi,lo
end

function accDiv{T<:StdFloat}(a::T, b::T)
    hi = a/b
    remainder = fma(-b,hi,a)
    lo = remainder/b
    remainder = fma(-b,lo,remainder)
    lo += remainder/b
    hi,lo
end

# for use with directed rounding, when only the sign of the lo part is needed
function accDivForSign{T<:StdFloat}(a::T, b::T)
    hi = a/b
    loApprox = -fma(b,hi,-a) # sign is correct, value is actual*b
    hi,loApprox
end

for fn in (:eftAdd, :eftAddGTE, :eftSub, :eftSubGTE, :eftMul, :accDiv, :accDivForSign)
  @eval begin
    ($fn){T1<:Real,T2<:StdFloat}(a::T1, b::T2) = ($fn)(convert(Float64,a), b)
    ($fn){T1<:Real,T2<:StdFloat}(a::T2, b::T1) = ($fn)(a, convert(Float64,b))
    ($fn){T1<:Real,T2<:Real}(a::T1, b::T2) = ($fn)(convert(Float64,a), convert(Float64,b))
  end
end


text:
mutable struct GebModel
    size::Int8

    organisms::Vector{GebOrganism}
    grid::Array{Union{GebOrganism, Nothing}, 2}

    next_key::UInt64

    time::UInt64

    logger::Logger

    function GebModel(; size=20)
        GebModel(DoNothingLogger(), size=size)
    end
    function GebModel(logger::Logger; size=20)
        initial_organisms = [
            GebOrganism(size*(x-1) + y, "0", (x - 0.5, y - 0.5))
            for x in 1:size for y in 1:size
        ]
        grid = [initial_organisms[size*(x-1) + y] for x in 1:size, y in 1:size]

        return new(size, initial_organisms, grid, size*size + 1, 0, logger)
    end
end

function _get_grid_coordinates(continuous_coordinates)
    (
        @fastmath max(1, Int(ceil(continuous_coordinates[1]))),
        @fastmath max(1, Int(ceil(continuous_coordinates[2])))
    )
end

function get_organism(model::GebModel, key::UInt64)
    for organism in model.organisms
        if organism.key == key return organism end
    end
end

function add_organism!(model::GebModel, organism::GebOrganism)
    x, y = _get_grid_coordinates(organism.coordinates)
    model.grid[x, y] = organism

    push!(model.organisms, organism)

    organism.key = model.next_key
    model.next_key += 1

    return organism
end

function kill!(model::GebModel, organism::GebOrganism)
    if length(model.organisms) <= MIN_ORGANISMS return end

    x, y = _get_grid_coordinates(organism.coordinates)
    model.grid[x, y] = nothing

    @assert organism in model.organisms
    delete!(model.organisms, organism)

    log_death(model.logger, model, organism)
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function fixup_blocks!(ir, n)
  last = length(blocks(ir))
  for bl in blocks(ir), i = 1:length(branches(bl))
    br = branches(bl)[i]
    if br.block > last
      branches(bl)[i] = Branch(br, block = br.block + n - 1)
    end
  end
end

function inlinehere!(ir, source, args...)
  source = merge_returns!(copy(source)) # TODO preserve type info
  offset = length(blocks(ir.to))-1
  env = Dict()
  retvalue = nothing
  rename(x::Variable) = env[x]
  rename(x::Expr) = Expr(x.head, rename.(x.args)...)
  rename(x::Statement) = stmt(x, expr = rename(x.expr))
  rename(x) = x
  for (name, arg) in zip(arguments(source), args)
    env[name] = arg
  end
  for bl in blocks(source)
    if bl.id != 1
      block!(ir)
      for (arg, T) in zip(arguments(bl), argtypes(bl))
        env[arg] = blockargument!(ir, T)
      end
    end
    for (v, st) in bl
      env[v] = push!(ir, rename(st))
    end
    for br in branches(bl)
      if isreturn(br)
        retvalue = rename(returnvalue(br))
      else
        branch!(ir, br.block+offset, rename.(br.args)..., unless = rename(br.condition))
      end
    end
  end
  return retvalue
end

"""
    inline(ir, location, source)

Replace the function call at `ir[location]` with the IR `source`. The inlined IR
will use the function arguments at `ir[location]` as its input.

```
julia> foo(x, y) = max(x, y)+1

julia> ir = @code_ir foo(1, 1)
1: (%1, %2, %3)
  %4 = Main.max(%2, %3)
  %5 = %4 + 1
  return %5

julia> inline(ir, var(4), @code_ir(max(1,1)))
1: (%1, %2, %3)
  %4 = %3 < %2
  %5 = Base.ifelse(%4, %2, %3)
  %6 = %5 + 1
  return %6
```
"""
function inline(ir::IR, loc::Variable, source::IR)
  pr = Pipe(ir)
  startblock = 0
  branches(pr) do br
    if startblock != 0 && br.block > startblock
      Branch(br, block = br.block + length(blocks(source)) - 1)
    else
      return br
    end
  end
  for (v, st) in pr
    if v === loc
      startblock = length(blocks(pr.to))
      fixup_blocks!(pr.to, length(blocks(source)))
      ex = ir[loc].expr
      delete!(pr, v)
      v′ = inlinehere!(pr, source, ex.args...)
      substitute!(pr, v, substitute(pr, v′))
    end
  end
  return finish(pr)
end


text:
"""
Convert a matrix A which is the output of a real FFT to a real vector, keeping
only unqiue real/imaginary entries of A
"""
function rfft2vec(A::AbstractMatrix)
    m,n = size(A)
    ireal,iimag = fftsyms(Val(m),Val(n))
    [real(A[ireal]); imag(A[iimag])]
end

"""
Convert a vector produced by rfft2vec back into a complex matrix.
"""
function vec2rfft(v::AbstractVector{<:Real})
    n = round(Int,sqrt(length(v)))
    m = n÷2+1
    nreal = (n^2)÷2 + (iseven(n) ? 2 : 1)
    ireal,iimag,inegks = fftsyms(Val(m),Val(n))
    A = fill(NaN+im*NaN,m,n)
    A[ireal] = v[1:nreal]
    A[iimag] += im*v[nreal+1:end]
    for i=1:m, j=1:n
        if isnan(A[i,j])
            A[i,j] = A[inegks[i,j]...]'
        end
    end
    A
end

"""
Convert an M×N matrix (with M=N÷2+1) which is the output a real FFT to a full
N×N one via symmetries.
"""
unfold(Tls::AbstractArray{<:Complex,3}, Ny) = mapslices(X -> unfold(X, Ny), Array(Tls), dims=(1,2))
unfold(Tl::AbstractMatrix{<:Complex}, Ny) = unfold(Array(Tl), Ny)
function unfold(Tl::Matrix{<:Complex}, Ny::Int)
    m,n = size(Tl)
    @assert m==Ny÷2+1
    m2 = iseven(Ny) ? 2m : 2m+1
    n2 = iseven(n) ?  n+2 : n+3
    Tlu = similar(Tl,Ny,n)
    Tlu[1:m,1:n] = Tl
    @inbounds for i=m+1:Ny
        Tlu[i,1] = Tl[m2-i, 1]'
        @simd for j=2:n
            Tlu[i,j] = Tl[m2-i, n2-j]'
        end
    end
    Tlu
end

"""
Arguments `m` and `n` refer to the sizes of an `m`×`n` matrix (call it `A`) that is the output of a
real FFT (thus `m=n÷2+1`)

Returns a tuple of (ireal, iimag, negks) where these are

* `ireal` — `m`×`n` mask corrsponding to unique real entries of `A`
* `iimag` — `m`×`n` mask corrsponding to unique imaginary entries of `A`
* `negks` — `m`×`n` matrix of giving the index into A where the negative k-vector
            is, s.t. `A[i,j] = A[negks[i,j]]'`
"""
@generated function fftsyms(::Val{m},::Val{n}) where {m,n}
    k = ifftshift(-n÷2:(n-1)÷2)
    ks = tuple.(k',k)[1:n÷2+1,:]
    wrapk(k) = mod(k+n÷2,n) - n÷2
    negk(k) = @. wrapk(-k)
    k_in_ks(k) = -n÷2<=k[1]<=(n-1)÷2 && (0<=k[2]<=(n-1)÷2 || k[2]==-n÷2)
    ireal = map(k->(negk(k)==k || !k_in_ks(negk(k)) || (k[1]>0 || k[2]>0)), ks)
    iimag = map(k->(negk(k)!=k && !k_in_ks(negk(k)) || (k[1]>0 || k[2]>0)), ks)
    indexof(k) = (mod(k[2],n)+1, mod(k[1],n)+1)
    inegks = indexof.(negk.(ks))
    inegks[.!k_in_ks.(negk.(ks))] .= Ref((0,0))
    ireal,iimag,inegks#,ks,negk.(ks)#,k_in_ks.(negk.(ks)),map(k->k_in_ks(negk(k)),ks)
end


text:
# PiecewiseLinearContinuousFitting
# Piecewise Linear Continouous Interpolation Constraints
# Sparse L2 Optimal Fitting subject Continuity Constraints
# Integrated Square
using EllZeroTrendFiltering: minimize_wrt_x2

using Polynomials
using IterTools

N = 240

t = range(0, stop=4π, length=50)

#g = Poly([0,0,1,-1])
ran = 0.1*randn(N)
z(t) = ran[floor(Int64, t/4pi*(N-1))+1]
g = sin
g2(t) = g(t) + z(t)

@time ℓ = compute_transition_costs(g, t);

K = 7
#@time I, Y, f = brute_force_search(ℓ, K-1);

Λ_0 = [create_new_pwq(minimize_wrt_x2(ℓ[i, N])) for i in 1:N-1];

@time Λ = pwq_dp_constrained(Λ_0, ℓ, 7);
tot

using Plots
plot(t, g.(t), lab="g(t)")
for k = 7:7
    @time I2, y2, f2 = recover_solution(Λ[k, 1], 1, N)
    Y2, _ = find_optimal_y_values(ℓ, I2)

    println("Comparison: ", sqrt(f), " --- ", sqrt(f2))

    #find_optimal_y_values

    l = zeros(size(Λ))
    for i=1:size(l,1)
        for j=1:size(l,2)
            if isassigned(Λ,i,j)
                l[i,j] = length(Λ[i, j])
            end
        end
    end

    #x = range(0, stop=1, length=100)

    #closefig()
    plot!(t[I2],Y2, m=:circle, lab="k=$k")
end
plot!()
#using PyPlot
# #close()
# figure(1)
# plot(t, g.(t), "g", t[I2], Y2, "bo-")
#
# figure(2)
# plot(1:length(t)-1, l')


text:
module Wanderer
using MLStyle
export @select, @where, @groupby, @having, @limit, @orderby

include("ConstantNames.jl")
include("DynamicInfer.jl")
include("Interfaces.jl")
include("MacroProcessors.jl")
include("Implementation.jl")

end # module


text:
module arr

end

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
function next_variable_index!(model::Optimizer)
    return model.last_variable_index_added += 1
end

function next_parameter_index!(model::Optimizer)
    return model.last_parameter_index_added += 1
end

function update_number_of_parameters!(model::Optimizer)
    return model.number_of_parameters_in_model += 1
end

function is_parameter_in_model(model::Optimizer, v::MOI.VariableIndex)
    return PARAMETER_INDEX_THRESHOLD <
           v.value <=
           model.last_parameter_index_added
end

function is_variable_in_model(model::Optimizer, v::MOI.VariableIndex)
    return 0 < v.value <= model.last_variable_index_added
end

function has_quadratic_constraint_caches(model::Optimizer)
    return !isempty(model.quadratic_added_cache)
end

function function_has_parameters(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
) where {T}
    for term in f.terms
        if is_parameter_in_model(model, term.variable_index)
            return true
        end
    end
    return false
end

function function_has_parameters(model::Optimizer, f::MOI.VectorOfVariables)
    for variable in f.variables
        if is_parameter_in_model(model, variable)
            return true
        end
    end
    return false
end

function function_has_parameters(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
) where {T}
    for term in f.terms
        if is_parameter_in_model(model, term.scalar_term.variable_index)
            return true
        end
    end
    return false
end

function function_has_parameters(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    return function_affine_terms_has_parameters(model, f.affine_terms) ||
           function_quadratic_terms_has_parameters(model, f.quadratic_terms)
end

function function_affine_terms_has_parameters(
    model::Optimizer,
    affine_terms::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    for term in affine_terms
        if is_parameter_in_model(model, term.variable_index)
            return true
        end
    end
    return false
end

function function_quadratic_terms_has_parameters(
    model::Optimizer,
    quadratic_terms::Vector{MOI.ScalarQuadraticTerm{T}},
) where {T}
    for term in quadratic_terms
        if is_parameter_in_model(model, term.variable_index_1) ||
           is_parameter_in_model(model, term.variable_index_2)
            return true
        end
    end
    return false
end

function separate_possible_terms_and_calculate_parameter_constant(
    model::Optimizer,
    terms::Vector{MOI.ScalarAffineTerm{T}},
) where {T}
    vars = MOI.ScalarAffineTerm{T}[]
    params = MOI.ScalarAffineTerm{T}[]
    param_constant = zero(T)

    for term in terms
        if is_variable_in_model(model, term.variable_index)
            push!(vars, term)
        elseif is_parameter_in_model(model, term.variable_index)
            push!(params, term)
            param_constant +=
                term.coefficient * model.parameters[term.variable_index]
        else
            error("Constraint uses a variable that is not in the model")
        end
    end
    return vars, params, param_constant
end

# This version is used on SQFs
function separate_possible_terms_and_calculate_parameter_constant(
    model::Optimizer,
    terms::Vector{MOI.ScalarAffineTerm{T}},
    variables_associated_to_parameters::Vector{MOI.VariableIndex},
) where {T}
    vars = MOI.ScalarAffineTerm{T}[]
    params = MOI.ScalarAffineTerm{T}[]
    terms_with_variables_associated_to_parameters = MOI.ScalarAffineTerm{T}[]
    param_constant = zero(T)

    for term in terms
        if is_variable_in_model(model, term.variable_index)
            push!(vars, term)
            if term.variable_index in variables_associated_to_parameters
                push!(terms_with_variables_associated_to_parameters, term)
            end
        elseif is_parameter_in_model(model, term.variable_index)
            push!(params, term)
            param_constant +=
                term.coefficient * model.parameters[term.variable_index]
        else
            error("Constraint uses a variable that is not in the model")
        end
    end
    return vars,
    params,
    terms_with_variables_associated_to_parameters,
    param_constant
end

function separate_possible_terms_and_calculate_parameter_constant(
    model::Optimizer,
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
) where {T}
    quad_params = MOI.ScalarQuadraticTerm{T}[] # parameter x parameter
    quad_aff_vars = MOI.ScalarQuadraticTerm{T}[] # parameter x variable
    quad_vars = MOI.ScalarQuadraticTerm{T}[] # variable x variable
    aff_terms = MOI.ScalarAffineTerm{T}[]
    variables_associated_to_parameters = MOI.VariableIndex[]
    quad_param_constant = zero(T)

    # TODO (guilherme) do we really need the convention ?
    # When we have a parameter x variable or a variable x parameter the convention is to rewrite
    # the SQT with parameter as variable_index_1 and variable as variable_index_2
    for term in terms
        if (
            is_variable_in_model(model, term.variable_index_1) &&
            is_variable_in_model(model, term.variable_index_2)
        )
            push!(quad_vars, term) # if there are only variables, it remains a quadratic term

        elseif (
            is_parameter_in_model(model, term.variable_index_1) &&
            is_variable_in_model(model, term.variable_index_2)
        )
            push!(quad_aff_vars, term)
            push!(variables_associated_to_parameters, term.variable_index_2)
            push!(model.multiplicative_parameters, term.variable_index_1.value)
            push!(
                aff_terms,
                MOI.ScalarAffineTerm(
                    term.coefficient * model.parameters[term.variable_index_1],
                    term.variable_index_2,
                ),
            )

        elseif (
            is_variable_in_model(model, term.variable_index_1) &&
            is_parameter_in_model(model, term.variable_index_2)
        )
            # Check convention defined above
            push!(
                quad_aff_vars,
                MOI.ScalarQuadraticTerm(
                    term.coefficient,
                    term.variable_index_2,
                    term.variable_index_1,
                ),
            )
            push!(variables_associated_to_parameters, term.variable_index_2)
            push!(
                aff_terms,
                MOI.ScalarAffineTerm(
                    term.coefficient * model.parameters[term.variable_index_2],
                    term.variable_index_1,
                ),
            )
            push!(model.multiplicative_parameters, term.variable_index_2.value)

        elseif (
            is_parameter_in_model(model, term.variable_index_1) &&
            is_parameter_in_model(model, term.variable_index_2)
        )
            push!(quad_params, term)
            push!(model.multiplicative_parameters, term.variable_index_1.value)
            push!(model.multiplicative_parameters, term.variable_index_2.value)
            quad_param_constant +=
                term.coefficient *
                model.parameters[term.variable_index_1] *
                model.parameters[term.variable_index_2]
        else
            throw(
                ErrorException(
                    "Constraint uses a variable or parameter that is not in the model",
                ),
            )
        end
    end
    return quad_vars,
    quad_aff_vars,
    quad_params,
    aff_terms,
    variables_associated_to_parameters,
    quad_param_constant
end

function fill_quadratic_constraint_caches!(
    model::Optimizer,
    new_ci::MOI.ConstraintIndex,
    quad_aff_vars::Vector{MOI.ScalarQuadraticTerm{T}},
    quad_params::Vector{MOI.ScalarQuadraticTerm{T}},
    aff_params::Vector{MOI.ScalarAffineTerm{T}},
    terms_with_variables_associated_to_parameters::Vector{
        MOI.ScalarAffineTerm{T},
    },
) where {T}
    if !isempty(quad_aff_vars)
        model.quadratic_constraint_cache_pv[new_ci] = quad_aff_vars
    end
    if !isempty(quad_params)
        model.quadratic_constraint_cache_pp[new_ci] = quad_params
    end
    if !isempty(aff_params)
        model.quadratic_constraint_cache_pc[new_ci] = aff_params
    end
    if !isempty(terms_with_variables_associated_to_parameters)
        model.quadratic_constraint_variables_associated_to_parameters_cache[new_ci] =
            terms_with_variables_associated_to_parameters
    end
    return nothing
end

function quadratic_constraint_cache_map_check(
    model::Optimizer,
    idxs::MOI.ConstraintIndex{F,S},
) where {F,S}
    cached_constraints = values(model.quadratic_added_cache)
    # Using this becuase some custom brodcast method throws errors if
    # inner_idex .∈ cached_constraints is used
    return [i ∈ cached_constraints for i in idxs]
end

# Vector Affine
function separate_possible_terms_and_calculate_parameter_constant(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    set::S,
) where {T,S<:MOI.AbstractVectorSet}
    vars = MOI.VectorAffineTerm{T}[]
    params = MOI.VectorAffineTerm{T}[]
    n_dims = length(f.constants)
    param_constants = zeros(T, n_dims)
    for term in f.terms
        oi = term.output_index

        if is_variable_in_model(model, term.scalar_term.variable_index)
            push!(vars, term)
        elseif is_parameter_in_model(model, term.scalar_term.variable_index)
            push!(params, term)
            param_constants[oi] +=
                term.scalar_term.coefficient *
                model.parameters[term.scalar_term.variable_index]
        else
            error("Constraint uses a variable that is not in the model")
        end
    end
    return vars, params, param_constants
end


text:
"""
    delete!(dict, k)

Deletes the entry for `k`, if it exists. Returns -1 if `k` was found, 0
if not found.
"""
function Base.delete!(dict::MuttsDict{K,V}, k::K) where {K,V}
    @dassert1 is_mutable(dict)
    _delete!(dict, hash(k), k)
    return dict
end

function _delete!(dict::MuttsDict{K,V}, hashcode::UInt64, k::K) where {K,V}
    root = dict.root
    pathmask = UInt64(0)
    pathhash = UInt64(0)
    if isa(root, Leaf{K,V})
        # Branch specialized for isa(root,Leaf{K,V})
        (!is_mutable(root)) && error("delete!($(typeof(dict)), ..): dict is immutable")
        (new_root,size_change) = _delete!(root, hashcode, k, pathmask, pathhash)
        dict.root = new_root
    else
        # Branch specialized for isa(root,INode{K,V})
        (!is_mutable(root)) && error("delete!($(typeof(dict)), ..): dict is immutable")
        config = dict.config
        height = length(config.inode_capacity)
        size_change = _delete!(root, config, height, hashcode, k,
            pathmask, pathhash)
    end
    dict.n += size_change
    return size_change
end

# delete!() is not designed to be a frequent operation, this method is
# going to always clone the leaf node. This is so that get() and insert!()
# can stop probing when they encounter an empty slot.
function _delete!(leaf::Leaf{K,V},
    hashcode::UInt64,
    k::K,
    pathmask::UInt64,
    pathhash::UInt64
) where {K,V}
    num_slots = length(leaf.entries)
    new_leaf = Leaf{K,V}(num_slots)

    found = false
    for i=1:num_slots
        slot = reinterpret(UInt64, i)
        if _slot_occupied(leaf, slot, num_slots)
            (ki,vi) = @inbounds leaf.entries[i]
            if (ki === k) || isequal(ki,k)
                found = true
                continue
            end
            h = hash(ki)
            if (h & pathmask) == pathhash
                (new_leaf,_) = _insert!(new_leaf, h, ki, vi)
            end
        end
    end

    if found
        return (new_leaf,-1)
    else
        return (leaf,0)
    end
end

function _delete!(inode::INode{K,V},
    config::NodeConfig,
    height::Int,
    hashcode::UInt64,
    k::K,
    pathmask=UInt64(0),
    pathhash=UInt64(0)
) where {K,V}
    @dassert1 inode.is_mutable
    mask = reinterpret(UInt64, length(inode.entries)-1)
    shift = reinterpret(UInt64, height+1)*UInt64(8)
    slot = 1 + ((hashcode >> shift) & mask)
    pathmask = pathmask | (mask << shift)
    pathhash = pathhash | ((slot-1) << shift)
    entries = inode.entries
    if isa(entries, Vector{Leaf{K,V}})
        @dassert1 height==1
        leaf = @inbounds entries[slot]
        (new_leaf,size_change) = _delete!(leaf, hashcode, k, pathmask, pathhash)
        if new_leaf != leaf
            entries[slot] = new_leaf
        end
        return size_change
    else
        child_inode = @inbounds entries[slot]
        child_inode = _get_mutable_version(child_inode, config, height-1)
        size_change = _delete!(child_inode, config, height-1, hashcode, k, pathmask, pathhash)
        @inbounds inode.entries[slot] = child_inode
        return size_change
    end
end


text:
plan_fejer1(μ) = FFTW.plan_r2r!(μ, FFTW.REDFT01)

"""
Compute nodes of Fejer's first quadrature rule.
"""
fejernodes1(::Type{T}, N::Int) where T = chebyshevpoints(T, N, Val(1))

"""
Compute weights of Fejer's first quadrature rule with modified Chebyshev moments of the first kind ``\\mu``.
"""
fejerweights1(μ::Vector) = fejerweights1!(copy(μ))
fejerweights1!(μ::Vector) = fejerweights1!(μ, plan_fejer1(μ))
function fejerweights1!(μ::Vector{T}, plan) where T
    N = length(μ)
    rmul!(μ, inv(T(N)))
    return plan*μ
end

plan_fejer2(μ) = FFTW.plan_r2r!(μ, FFTW.RODFT00)

"""
Compute nodes of Fejer's second quadrature rule.
"""
fejernodes2(::Type{T}, N::Int) where T = T[sinpi((N-2k-one(T))/(2N+two(T))) for k=0:N-1]

"""
Compute weights of Fejer's second quadrature rule with modified Chebyshev moments of the second kind ``\\mu``.
"""
fejerweights2(μ::Vector) = fejerweights2!(copy(μ))
fejerweights2!(μ::Vector) = fejerweights2!(μ, plan_fejer2(μ))
function fejerweights2!(μ::Vector{T}, plan) where T
    N = length(μ)
    Np1 = N+one(T)
    rmul!(μ, inv(Np1))
    plan*μ
    @inbounds for i=1:N μ[i] = sinpi(i/Np1)*μ[i] end
    return μ
end


text:
"""
# abstract.jl

A collection of geometrical data types
==============================================================
Author: Michael Kuegeler
Email: mkuegeler@mac.com
Date: 17.5.2020
==============================================================

## About
This script comprises a collection of data types, related to generic two-dimensional geometries.
The intend is to provide a set of data types which facilitate the creation of parametric graphics.

The script is structured into these sections:

## Section 1: Essential data types
## Section 2: Composite data types

"""

"""
==============================================================

# Section 1: Essential data types
==============================================================
"""

# default digits
global d = 3

"""
Abstract point definition.
fields: x,y,z coordinates
"""

mutable struct Point
   x::Float64
   y::Float64
   z::Float64 # z is optional
    function Point(x,y,z::Float64=0.0)
        new(Float64(x, RoundUp) , Float64(y, RoundUp) , Float64(z, RoundUp))
    end
end

"""
Abstract line definition.
fields: start and end point (type of point is previously defined Point)
"""

mutable struct Line
   p1::Point
   p2::Point
    function Line(p1,p2)
        new(p1, p2)
    end
end

"""
Raw rectangle
"""

# Rectangle support function: Calculation of coordinates, depending on inion point
function get_rectangle_points(p0::Point,w::Float64,h::Float64)
    # clock wise direction
    # p0 = top left (default)
    p1 = Point((p0.x+(w/2)),p0.y) # mid top
    p2 = Point((p0.x+w),p0.y) # top right
    p3 = Point((p0.x+w),(p0.y+(h/2))) # mid right
    p4 = Point((p0.x+w),(p0.y+h)) # buttom right
    p5 = Point((p0.x+(w/2)),(p0.y+h)) # mid buttom
    p6 = Point(p0.x,(p0.y+h)) # buttom left
    p7 = Point(p0.x,(p0.y+(h/2))) # mid left
    p8 = Point((p0.x+(w/2)),(p0.y+(h/2))) # center

    Dict("tl"=>p0,"mt"=>p1,"tr"=>p2,"mr"=>p3,"br"=>p4,"mb"=>p5,"bl"=>p6,"ml"=>p7,"cp"=>p8)
end

mutable struct Rectangle
    in::Point # in point of rectangle
    w::Float64 # w
    h::Float64 # h
    start::Int64 # indicates the location of inion point
    # 0 = top left (default) 1 = top right, 2 = buttom right, 3 = buttom left, 4 = center
    points::Dict # Dict with all four edge points of the rectangle plus center point
    r::Float64 # maximum value of a radius within the rectangle, depends on w and h
    function Rectangle(in::Point=Point(0.0,0.0,0.0),w::Float64=100.0,h::Float64=100.0,start::Int64=0)
        r = w >= h ? (h/2) : (w/2)
        if start == 1 # inion: top right
            p = Point((in.x-w),in.y)
        elseif start == 2 # inion: buttom right
            p = Point( (in.x-w),(in.y-h) )
        elseif start == 3 # inion: buttom left
            p = Point(in.x, (in.y+h) )
        elseif start == 4 # inion: center
            p = Point( (in.x-(w/2) ) ,(in.y-(h/2)) )
        else # default
            p = in
        end
        points = get_rectangle_points(p,w,h)

        new(in, w, h, start, points, r)
    end
end

# Wrapper function for struct Rectangle: Convert string values to Float64 and returns a tuple

# Call with default variables

function abstract_rectangle()
        rect = Rectangle()
        points  = Dict()
        for (key,value) in rect.points
            push!(points, key => (x=string(value.x), y=string(value.y), z=string(value.z)))
        end
        (w=string(rect.w), h=string(rect.h), points=points, r=string(rect.r))
end
# x,y,w,h,start
function abstract_rectangle(at::Dict)
        x,y,w,h,start = at["x"],at["y"],at["w"],at["h"],at["start"]
        # We may get the values as different types, so let's convert them first into strings
        # Next, convert them into types required by the Rectangle composite
        p = map(v -> parse(Float64, v), map(v -> string(v), [x,y,w,h]))
        rect = Rectangle(Point(p[1],p[2]),p[3],p[4],parse(Int64, string(start)))

        points  = Dict()
        for (key,value) in rect.points
            push!(points, key => (x=string(value.x), y=string(value.y), z=string(value.z)))
        end
        (w=string(rect.w), h=string(rect.h), points=points, r=string(rect.r))
end

"""
Abstract ScaleneTriangle
A scalene triangle has all its sides of different lengths

Other types are:
By lengths of sides: Equilateral, Isosceles

By internal angles: Right, Obtuse, Acute

"""
mutable struct ScaleneTriangle
    pa::Point
    pb::Point
    pc::Point
    function ScaleneTriangle(pa,pb,pc)

        new(pa, pb, pc)
    end

end

"""
==============================================================

# Section 2: Composite data types
==============================================================
"""

"""
# Abstract Grid
"""

mutable struct Grid
    in::Point
    nx::UInt16
    ny::UInt16
    w::Float64
    h::Float64
    f::UInt16
    grid::Array{Line,1}
    function Grid(in=Point(0.0,0.0,0.0), nx=5, ny=5,  w=100.0, h=100.0, f=1)
        # if factor f is larger or equal to either nx or ny, set f to 1
        f = (f >= nx || f >= ny) ? 1 : f

        grid = Line[]

        for i= 1:(nx-f)

            x1 = round( (in.x+((w/nx)*i)), digits=d)
            y1 = round(in.y, digits=d)

            x2 = round((in.x+((w/nx)*i)),digits=d)
            y2 = round((in.y+h), digits=d)

            push!(grid, Line(Point(x1,y1),Point(x2,y2)))

        end

        for j= 1:(ny-f)

            x1 = round(in.x, digits=d)
            y1 = round((in.y+((h/ny)*j)), digits=d)

            x2 = round((in.x+w),digits=d)
            y2 = round((in.y+((h/ny)*j)), digits=d)

            push!(grid, Line(Point(x1,y1),Point(x2,y2)))

        end

        new(in, nx, ny, w, h, f, grid)
    end

end

# Wrapper for Grid
function abstract_grid()
    g = Grid()
    lines  = Array{Any,1}()
    for l in g.grid
         push!(lines, (x1=string(l.p1.x), y1=string(l.p1.y), x2=string(l.p2.x), y2=string(l.p2.y)))
    end
    (x=string(g.in.x), y=string(g.in.y),nx=string(g.nx), ny=string(g.ny), w=string(g.w), h=string(g.h), lines=lines,f=string(g.f) )
end

function abstract_grid(at::Dict)
    x,y,nx,ny,w,h,f = at["x"],at["y"],at["nx"],at["ny"],at["w"],at["h"],at["f"]
    p = map(v -> parse(Float64, v), map(v -> string(v), [x,y,w,h]))
    nx,ny = parse(Int16,string(nx)),parse(Int16,string(ny))
    g = Grid(Point(p[1],p[2]),nx,ny,p[3],p[4],parse(Int16,string(f)))

    lines  = Array{Any,1}()
    for l in g.grid
         push!(lines, Dict("x1"=>string(l.p1.x), "y1"=>string(l.p1.y),"x2"=>string(l.p2.x), "y2"=>string(l.p2.y)))
    end
    (x=string(g.in.x), y=string(g.in.y),nx=string(g.nx), ny=string(g.ny), w=string(g.w), h=string(g.h), lines=lines,f=string(g.f) )
end

"""
Abstract NodeGrid
"""

mutable struct NodeGrid
    in::Point # inion: top left
    nx::UInt16
    ny::UInt16
    w::Float64
    h::Float64
    f::UInt16
    nodes::Array{Point,1}
    function NodeGrid(in=Point(0.0,0.0,0.0), nx=5, ny=5,  w=100.0, h=100.0, f=1)
        # if factor f is larger or equal to either nx or ny, set f to 1
        f = (f >= nx || f >= ny) ? 1 : f
        nodes = Point[]
        for j= 1:(ny-f)
           for i= 1:(nx-f)
             x = round((in.x+((w/nx)*i)), digits=d)
             y = round((in.y+((h/ny)*j)), digits=d)
             push!(nodes, Point(x,y))
           end
        end
        new(in, nx, ny, w, h, f, nodes)
    end

end

# Wrapper for NodeGrid
function abstract_nodegrid()
    n = NodeGrid()
    points  = Array{Any,1}()
    for e in n.nodes
        push!(points, Dict("x"=>string(e.x), "y"=>string(e.y)))

    end
    (x=string(n.in.x), y=string(n.in.y), w=string(n.w), h=string(n.h), points=points, f=string(n.f))
end

function abstract_nodegrid(at::Dict)
    x,y,nx,ny,w,h,f = at["x"],at["y"],at["nx"],at["ny"],at["w"],at["h"],at["f"]
    p = map(v -> parse(Float64, v), map(v -> string(v), [x,y,w,h]))
    nx,ny = parse(Int16,string(nx)),parse(Int16,string(ny))
    n = NodeGrid(Point(p[1],p[2]),nx,ny,p[3],p[4],parse(Int16,string(f)))
    points  = Array{Any,1}()
    for e in n.nodes
        push!(points, Dict("x"=>string(e.x), "y"=>string(e.y), "z"=>string(e.z)))

    end
    (x=string(n.in.x), y=string(n.in.y), w=string(n.w), h=string(n.h), points=points, f=string(n.f))
end

"""
Abstract PanelGrid
nx,ny,in::Point,w,h,of
in, nx, ny,  w, h, f
"""

mutable struct PanelGrid
    in::Point # inion: top left
    nx::UInt16
    ny::UInt16
    w::Float64
    h::Float64
    of::Float64
    panels::Array{Rectangle,1}
    function PanelGrid(in=Point(0.0,0.0,0.0), nx=5, ny=5,  w=100.0, h=100.0, of=1.0)
        # if of is larger or equal to either panel w or panel h, set of to 0
        of = (of >= (w/nx) || of >= (h/ny)) ? 0 : of

        panels = Rectangle[]

        x = round((in.x+(of/2)), digits=d)
        y = round((in.y+(of/2)), digits=d)

        x1 = x
        y1 = y

        for j= 1:ny
            for i= 1:nx
               push!(panels, Rectangle(Point(x1,y1),round(((w/nx)-of), digits=d),round(((h/ny)-of), digits=d)))
               x1 = round(((in.x+(of/2))+((w/nx)*i)), digits=d)
            end
            x1 = round((in.x+(of/2)),digits=d)
            y1= round(((in.y+(of/2))+((h/ny)*j)), digits=d)
        end

        new(in, nx, ny, w, h, of, panels)
    end

end

# Wrapper function PanelGrid
function abstract_panelgrid()
    g = PanelGrid()
    panels  = Array{Any,1}()
    for r in g.panels
        push!(panels, Dict("x"=>string(r.in.x), "y"=>string(r.in.y), "w"=>string(r.w), "h"=>string(r.h), "points"=>r.points))
    end
    (x=string(g.in.x), y=string(g.in.y), w=string(g.w), h=string(g.h), panels=panels, of=string(g.of))
end

function abstract_panelgrid(at::Dict)
    x,y,nx,ny,w,h,of = at["x"],at["y"],at["nx"],at["ny"],at["w"],at["h"],at["of"]
    p = map(v -> parse(Float64, v), map(v -> string(v), [x,y,w,h,of]))
    nx,ny = parse(Int16,string(nx)),parse(Int16,string(ny))

    g = PanelGrid(Point(p[1],p[2]),nx,ny,p[3],p[4],p[5])
    panels  = Array{Any,1}()
    for r in g.panels
        push!(panels, Dict("x"=>string(r.in.x), "y"=>string(r.in.y), "w"=>string(r.w), "h"=>string(r.h), "r"=>string(r.r), "points"=>r.points))
    end
    (x=string(g.in.x), y=string(g.in.y), w=string(g.w), h=string(g.h), panels=panels, of=string(g.of))
end


text:
using DiscreteEvents, Distributions

function boy(sim::Simulation, simlog, nr::Int64, bowls)
  roll = logvar("Boy"*string(nr), 0)
  logvar2log(simlog, roll)
  while true
    yield(Timeout(sim, 1))
    mystock = bowls[nr]
    roll.value = rand(DiscreteUniform(1, 6))
    matches = min(roll.value, mystock)
    bowls[nr]   -= matches
    bowls[nr+1] += matches
  end
end

# initalization
srand(1234)       # seed random number generator for reproducibility
simlog = newlog()
bowls = Dict(1=>999, 2=>0, 3=>0, 4=>0, 5=>0, 6=>0) # initialize bowls with matches
dict2log(simlog, bowls, redict=Dict(1=>"Bowl1",2=>"Bowl2",3=>"Bowl3",
                                    4=>"Bowl4",5=>"Bowl5",6=>"Output"))
sim = Simulation()
for i = 1:5
  @process boy(sim, simlog, i, bowls)  # create 5 players
end
@process logtick(sim, simlog, 1)         # add logging after each run
run(sim, 101)
d = log2df(simlog);
tail(d)


text:
using ThreeBodyDecay
using Test

ms = ThreeBodyMasses(1.0,2.0,3.0; m0=120.0)

@testset "Three Body Masses structure" begin
    #
    @test ms.m1 == ms[1] == 1
    @test ms.m2 == ms[2] == 2
    @test ms.m3 == ms[3] == 3.0
    @test ms.m0 == ms[0] == ms[4] == 120
    #
    @test_throws ErrorException ms[5]
    @test_throws ErrorException ThreeBodyMasses(0,1,1; m0=1)
    @test_throws ErrorException ThreeBodyMasses(0,1,1)
end

@testset "operations and interate" begin
    @test sum(ms^2) == sum(abs2, ms)
end

@testset "Limits for ThreeBodyMasses" begin
    m1 = 0.938; m2 = 0.49367; m3 = 0.13957; m0 = 2.46867
    ms = ThreeBodyMasses(m1=m1,m2=m2,m3=m3,m0=m0)
    #
    @test lims1(ms)[1] == (m2+m3)^2
    @test lims2(ms)[1] == (m3+m1)^2
    @test lims3(ms)[1] == (m1+m2)^2
    #
    @test lims1(ms)[2] == (m0-m1)^2
    @test lims2(ms)[2] == (m0-m2)^2
    @test lims3(ms)[2] == (m0-m3)^2
    #
    @test lims1(ms) == lims(1, ms)
    @test lims2(ms) == lims(2, ms)
    @test lims3(ms) == lims(3, ms)
    #
    σs = randomPoint(ms)
    @test ThreeBodyDecay.inrange(σs[1], lims1(ms))
    @test ThreeBodyDecay.inrange(σs[2], lims2(ms))
    @test ThreeBodyDecay.inrange(σs[3], lims3(ms))
    @test inphrange(σs, ms)
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text: x = 2


text:
# Particle Gibbs sampler

immutable PG <: InferenceAlgorithm
  n_particles :: Int
  n_iterations :: Int
  resampler :: Function
  resampler_threshold :: Float64
  PG(n1::Int,n2::Int) = new(n1,n2,resampleSystematic,0.5)
end

assume(spl :: Sampler{PG}, distr :: Distribution)   = randr(current_trace(), distr)

function Base.run(spl::Sampler{PG})
  chain = Chain()
  logevidence = Vector{Float64}(spl.alg.n_iterations)

  ## custom resampling function for pgibbs
  ## re-inserts reteined particle after each resampling step
  ref_particle = nothing
  for tt = 1:spl.alg.n_iterations
    dprintln(-1, "[PG]: Iter $(tt) out of $(spl.alg.n_iterations)")
    spl.particles = ParticleContainer{TraceR}(spl.model)
    if ref_particle == nothing
      push!(spl.particles, spl.alg.n_particles)
    else
      push!(spl.particles, spl.alg.n_particles-1)
      push!(spl.particles, ref_particle)
    end

    while consume(spl.particles) != Val{:done}
      ess = effectiveSampleSize(spl.particles)
      if ess <= spl.alg.resampler_threshold * length(spl.particles)
        resample!(spl.particles, spl.alg.resampler, ref_particle)
      end
    end

    logevidence[tt] = spl.particles.logE
    ## pick a particle to be retained.
    Ws, _ = weights(spl.particles)
    indx = rand(Categorical(Ws))
    ref_particle = fork2(spl.particles[indx])

    s = getsample(spl.particles, indx)
    push!(chain, Sample(1/spl.alg.n_iterations, s.value))
  end
  chain.weight = exp(mean(logevidence))
  chain
end


text:
using Documenter, MetidaStats, Weave, PrettyTables, CSV, DataFrames
#using DocumenterLaTeX

makedocs(
        modules = [MetidaStats],
        sitename = "MetidaStats.jl",
        authors = "Vladimir Arnautov",
        pages = [
            "Home" => "index.md",
            ],
        )

deploydocs(repo = "github.com/PharmCat/MetidaStats.jl.git", devbranch = "main", forcepush = true
)


text:
using Benchmark

function krons()
  for i=1:3000
    kron(eye(5), eye(5))
  end
end

function muls()
  for i=1:3000
    eye(5) * eye(5)
  end
end

walltimes = Array(Float64, 10)
for i=1:10
  d = benchmark(krons, "Kronecker Products of two 5x5 eyes", 1)
  walltimes[i] = float(d[4])[1]
end
output = "E:\\Cherokee\\simulator\\kron_output.txt"
f = open(output, "w")
print(f, walltimes)
close(f)
println("DONE")

max_stars_repo_name: younisshah/cherokee

max_issues_repo_name: younisshah/cherokee

max_forks_repo_name: younisshah/cherokee


text:
# Helpers for creating data containers and iterators over
# encoded obseravations.

# ## Data container

"""
    taskdataset(data, task, context)

Transform data container `data` of samples into a data container
of encoded samples.
Maps `encodesample(task, context, sample)` over the observations in
`data`. Also handles in-place `MLUtils.getobs!` through `encodesample!`.
"""
struct TaskDataset{TData, TTask<:LearningTask, TContext<:Context}
    data::TData
    task::TTask
    context::TContext
end

Base.length(ds::TaskDataset) = numobs(ds.data)

function Base.getindex(ds::TaskDataset, idx)
    return encodesample(ds.task, ds.context, getobs(ds.data, idx))
end

function MLUtils.getobs!(buf, ds::TaskDataset, idx)
    return encodesample!(buf, ds.task, ds.context, getobs(ds.data, idx))
end

"""
    taskdataset(data, task, context)

Transform data container `data` of samples into a data container of `(x, y)`-pairs.
Maps `encodesample(task, context, sample)` over the observations in `data`.
"""
const taskdataset = TaskDataset

# ## Data iterator

"""
    taskdataloaders(data, task[, batchsize])
    taskdataloaders(traindata, validdata, task[, batchsize; shuffle = true, dlkwargs...])

Create training and validation `DataLoader`s from two data containers
`(traindata, valdata)`. If only one container `data` is passed, splits
it into two, with `pctgvalid`% of the data going into the validation split.

## Arguments

Positional:
- `batchsize = 16`

Keyword:
- `shuffle = true`: Whether to shuffle the training data container
- `validbsfactor = 2`: Factor to multiply batchsize for validation data loader with (validation
    batches can be larger since no GPU memory is needed for the backward pass)

All remaining keyword arguments are passed to [`DataLoader`](#).

## Examples

Basic usage:

```julia
traindl, validdl = taskdataloaders(data, task, 128)
```

Explicit validation data container and no shuffling of training container:

```julia
traindl, validdl = taskdataloaders(traindata, validdata, task, shuffle=false)
```

Customizing the [`DataLoader`](#)

```julia
traindl, validdl = taskdataloaders(data, task, parallel=false, buffered=false)
```
"""
function taskdataloaders(
        traindata,
        validdata,
        task::LearningTask,
        batchsize = 16;
        shuffle = true,
        validbsfactor = 2,
        parallel = true,
        collate = true,
        kwargs...)
    return (
        DataLoader(taskdataset(traindata, task, Training()); batchsize, shuffle, collate, parallel, kwargs...),
        DataLoader(taskdataset(validdata, task, Validation());
                   batchsize = validbsfactor * batchsize, collate, parallel, kwargs...),
    )
end

function taskdataloaders(
        data,
        task::LearningTask,
        batchsize = 16;
        pctgval = 0.2,
        kwargs...)
    traindata, validdata = splitobs(shuffleobs(data), at = 1-pctgval)
    taskdataloaders(traindata, validdata, task, batchsize; kwargs...)
end


text:
"""
    tempcorr(w::CTASensor, E, temp, T0)

Corrects the anemometer output due to differences in flow temperature.

 * `w::CTASensor` object
 * `E` Anemometer output
 * `temp` Temperature during measurement
 * `T0` Reference calibration temperature. If not provided, use ref. temp of Rw

This correction assumes that only the temperature changed and not fluid properties. Another function that takes that into account will be developped later.

The correction works by using the definition of convective heat transfer:

``Q = h A (Tw - Too)``

where ``Tw`` is the operating temperature and ``Too`` is the flow temperature. Imagine that for the *same* flow speed the flow temperature changes to ``Too'``. Now, the assuming the the convection coefficient (``h``) remains the same, the heat output will change:

``Q' = h A (Tw - Too')``

Since the anemometer output is proportial to the voltage accross the sensor element, dividing the above equations and setting

``Q = α²E²/Rw``

we get the following equation for output correction:

``E' = E × ( (Tw - Too') / (Tw - Too) )^(1/2)``

If between the resistance element and the flow is an insulator (very common in thermistors and hot film probes), the correction equation remains the same event though this insulator should be taken into account:

``Q = h A (Tw - Too - Q⋅β)``

where ``Q⋅β`` corrects the sensor temperature (``Tw``) to the outer surface mean temperature.

"""
function tempcorr(w::CTASensor, E, temp, T0)
    Tw = optemperature(w)
    return E * sqrt( (Tw - T0) / (Tw - temp) )
end

tempcorr(w::CTASensor, E, temp) = tempcorr(w, E, temp, temperature(w.R))

"""
    tempcorr(w::CCASensor, E, temp, T0)

Corrects the anemometer output due to differences in flow temperature.

 * `w::CCASensor` Constant current sensor
 * `E` Anemometer output
 * `temp` Temperature during measurement
 * `T0` Reference calibration temperature. If not provided, use ref. temp of R.

This correction assumes that only the temperature changed and not fluid properties. Another function that takes that into account will be developped later.

The correction works by using the definition of convective heat transfer:

``Q = h A (Tw - Too)``

where ``Tw`` is the operating temperature and ``Too`` is the flow temperature. Imagine that for the *same* flow speed the flow temperature changes to ``Too'``. Now, the assuming the the convection coefficient (``h``) remains the same, the heat output will change:

``Q' = h A (Tw - Too')``

Since the anemometer output is proportial to the voltage accross the sensor element, dividing the above equations and setting

``Q = α²E²/Rw``

we get the following equation for output correction:

``E' = E × ( (Tw - Too') / (Tw - Too) )^(1/2)``

If between the resistance element and the flow is an insulator (very common in thermistors and hot film probes), the correction equation remains the same event though this insulator should be taken into account:

``Q = h A (Tw - Too - Q⋅β)``

where ``Q⋅β`` corrects the sensor temperature (``Tw``) to the outer surface mean temperature.

"""
function tempcorr(w::CCASensor{Resistor}, E, temp, T0)

    α = w.R.α
    return E / (1.0 + α*(temp - T0))
end

tempcorr(w::CCASensor{Resistor}, E, temp) = tempcorr(w, E, temp, temperature(w.R))

function funroot(f, x1, x2, eps=1e-7, maxiter=100)
    y1 = f(x1)
    y2 = f(x2)
    if abs(y1) < eps return x1 end
    if abs(y2) < eps return x2 end

    err = 0.0

    for i in 1:maxiter
        a = (y2-y1) / (x2 - x1)
        b = y1 - a*x1
        xn = -b/a
        err = abs(xn - x2)
        if err < eps
            return xn
        end
        x1 = x2
        x2 = xn
        y1 = y2
        y2 = f(x2)
        if abs(y2) < eps
            return x2
        end

    end

    error("funroot failed to converge after $maxiter iterations and with $err error!")
    return x2
end

function tempcorr(w::CCASensor, E, temp, T0; eps=1e-7, maxiter=100)

    R = w.R
    g = gain(w)
    I = current(w)
    Tw = temperature(R, E/(g*I))
    To = T0
    ΔT = Tw - temp

    f = Eo -> E*(temperature(R, Eo/(g*I)) - To) - Eo * ΔT
    return funroot(f, E, 1.01*E, eps, maxiter)
    #error("tempcorr(w::CCASensor, $E, $temp failed to converge in $maxiter iterations with $err error!")
    return Eo
end

tempcorr(w::CCASensor, E, temp; eps=1e-7, maxiter=100) =
    tempcorr(w, E, temp, temperature(w.R), eps=eps, maxiter=maxiter)


text:
using NFFT, LinearAlgebra#, CuNFFT

include("../Wrappers/NFFT3.jl")
include("../Wrappers/FINUFFT.jl")

NFFT.FFTW.set_num_threads(Threads.nthreads())
ccall(("omp_set_num_threads",NFFT3.lib_path_nfft),Nothing,(Int64,),convert(Int64,Threads.nthreads()))
@info ccall(("nfft_get_num_threads",NFFT3.lib_path_nfft),Int64,())
NFFT._use_threads[] = (Threads.nthreads() > 1)

function nfft_performance_simple(;N = (1024,1024), M = prod(N), m = 4,
  σ = 2.0, threading=false, pre=NFFT.LINEAR, T=Float64, blocking=true,
  storeDeconvolutionIdx=false, fftflags=NFFT.FFTW.MEASURE, ctor=NFFTPlan)

  timing = TimingStats()
  x = T.(rand(length(N),M) .- 0.5)
  x .= sortslices(x, dims=2) # sort nodes to gain cache locality

  fHat = randn(Complex{T}, M)
  f = randn(Complex{T}, N)

  #=if ctor == CuNFFT.CuNFFTPlan
    fHat = CuNFFT.CuArray(fHat)
    f = CuNFFT.CuArray(f)
  end=#

  NFFT._use_threads[] = threading
  NFFT.FFTW.set_num_threads( threading ? Threads.nthreads() : 1)

  tpre = @elapsed p = ctor(x, N; m, σ, window=:kaiser_bessel, precompute=pre,
                            fftflags, storeDeconvolutionIdx, blocking)

  tadjoint = @elapsed mul!(f, adjoint(p), fHat; timing)
  ttrafo = @elapsed mul!(fHat, p, f; timing)

  @info tpre, ttrafo, tadjoint

  println(timing)

end


text:
# This file was generated, do not modify it. # hide
Xtrain = X[train,:]
ytrain = y[train];


text:
using Libdl
print(abspath(Libdl.dlpath((ccall(:jl_is_debugbuild, Cint, ()) != 0) ? "libjulia-debug" : "libjulia")))

max_stars_repo_name: dmikushin/libcxxwrap-julia

max_issues_repo_name: dmikushin/libcxxwrap-julia

max_forks_repo_name: dmikushin/libcxxwrap-julia


text:
function getsolution(filename,dmd,npe)

nproc = length(dmd);
if nproc==1
    fn = filename * "_np0.bin";
    UDG = reinterpret(Float64,read(fn));
    ne = length(dmd[1].elempart[:]);
    nc = Int64(round(length(UDG)/(npe*ne)));
    UDG = reshape(UDG,(npe,nc,ne));
else
    nei = zeros(Int,nproc);
    for i = 1:nproc
        nei[i] = sum(dmd[i].elempartpts[1:2]);
    end
    ne = sum(nei);

    fn = filename * "_np0.bin";
    tmp = reinterpret(Float64,read(fn));
    nc = Int64(round(length(tmp)/(npe*nei[1])));

    UDG = zeros(npe,nc,ne);
    for i = 1:nproc
        elempart = dmd[i].elempart[1:nei[i]];
        fn = filename * "_np" * string(i-1) * ".bin"
        UDG[:,:,elempart]  = reshape(reinterpret(Float64,read(fn)),(npe,nc,nei[i]));
    end
end

return UDG

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
############################################################
# joAbstractLinearOperatorInplace ##########################
############################################################

# includes
include("joAbstractLinearOperatorInplace/constructors.jl")
include("joAbstractLinearOperatorInplace/base_functions.jl")
include("joAbstractLinearOperatorInplace/extra_functions.jl")


text:
"""
    @kwdef typedef
This is a helper macro that automatically defines a keyword-based constructor for the type
declared in the expression `typedef`, which must be a `struct` or `mutable struct`
expression. The default argument is supplied by declaring fields of the form `field::T =
default`. If no default is provided then the default is provided by the `kwdef_val(T)`
function.
```julia
@kwdef struct Foo
    a::Cint            # implied default Cint(0)
    b::Cint = 1        # specified default
    z::Cstring         # implied default Cstring(C_NULL)
    y::Bar             # implied default Bar()
end
```
"""
macro kwdef(expr)
    expr = macroexpand(__module__, expr) # to expand @static
    expr isa Expr && expr.head === :struct || error("Invalid usage of @kwdef")
    expr = expr::Expr
    T = expr.args[2]
    if T isa Expr && T.head === :<:
        T = T.args[1]
    end

    params_ex = Expr(:parameters)
    call_args = Any[]

    _kwdef!(expr.args[3], params_ex.args, call_args)
    # Only define a constructor if the type has fields, otherwise we'll get a stack
    # overflow on construction
    if !isempty(params_ex.args)
        if T isa Symbol
            kwdefs = :(($(esc(T)))($params_ex) = ($(esc(T)))($(call_args...)))
        elseif T isa Expr && T.head === :curly
            T = T::Expr
            # if T == S{A<:AA,B<:BB}, define two methods
            #   S(...) = ...
            #   S{A,B}(...) where {A<:AA,B<:BB} = ...
            S = T.args[1]
            P = T.args[2:end]
            Q = Any[U isa Expr && U.head === :<: ? U.args[1] : U for U in P]
            SQ = :($S{$(Q...)})
            kwdefs = quote
                ($(esc(S)))($params_ex) =($(esc(S)))($(call_args...))
                ($(esc(SQ)))($params_ex) where {$(esc.(P)...)} =
                    ($(esc(SQ)))($(call_args...))
            end
        else
            error("Invalid usage of @kwdef")
        end
    else
        kwdefs = nothing
    end
    quote
        Base.@__doc__($(esc(expr)))
        $kwdefs
    end
end

# @kwdef helper function
# mutates arguments inplace
function _kwdef!(blk, params_args, call_args)
    for i in eachindex(blk.args)
        ei = blk.args[i]
        if ei isa Symbol
            #  var
            push!(params_args, ei)
            push!(call_args, ei)
        elseif ei isa Expr
            if ei.head === :(=)
                lhs = ei.args[1]
                if lhs isa Symbol
                    #  var = defexpr
                    var = lhs
                elseif lhs isa Expr && lhs.head === :(::) && lhs.args[1] isa Symbol
                    #  var::T = defexpr
                    var = lhs.args[1]
                else
                    # something else, e.g. inline inner constructor
                    #   F(...) = ...
                    continue
                end
                defexpr = ei.args[2]  # defexpr
                push!(params_args, Expr(:kw, var, esc(defexpr)))
                push!(call_args, var)
                blk.args[i] = lhs
            elseif ei.head === :(::) && ei.args[1] isa Symbol
                # var::Typ
                dec = ei # var::Typ
                var = dec.args[1] # var
                def = :(kwdef_val($(ei.args[2])))
                push!(params_args, Expr(:kw, var, def))
                push!(call_args, dec.args[1])
            elseif ei.head === :block
                # can arise with use of @static inside type decl
                _kwdef!(ei, params_args, call_args)
            end
        end
    end
    blk
end

"""
    kwdef_val(T)
The default value for a type for use with the `@kwdef` macro. Returns:
 - null pointer for pointer types (`Ptr{T}`, `Cstring`, `Cwstring`)
 - zero for integer types
 - no-argument constructor calls (e.g. `T()`) for all other types
"""
function kwdef_val end

kwdef_val(::Type{Ptr{T}}) where T = Ptr{T}(C_NULL)
kwdef_val(::Type{Cstring}) = Cstring(C_NULL)
kwdef_val(::Type{Cwstring}) = Cwstring(C_NULL)
kwdef_val(::Type{T}) where T<: Real = zero(T)
kwdef_val(::Type{T}) where T = T()
kwdef_val(::Type{T}) where T<: String = ""
kwdef_val(::Type{T}) where T<: Symbol = :NULL
kwdef_val(::Type{Array{T,N}}) where {T,N} = Array{T}(undef, zeros(Int,N)...)

# @kwdef type Foo
#     a::Cint            # implied default Cint(0)
#     b::Cint = 1        # specified default
#     z::Cstring         # implied default Cstring(C_NULL)
#     y::Float64             # implied default Bar()
#     w::Bool             # implied default Bar()
#     ss::String
# end

#############################
struct heuristicOptions
    constraint_per_iteration::Int
end

@kwdef mutable struct heuristicModel
    model::Model = Model()
    functions::Array = []
    sets::Array = []
    violation::Array{Float64,1} = []
    functionValue::Array = []

end


text:
# see #1305
@test isfinite(mean(Chi(1000)))


text:
# This is based on the COCO (COmparing Continuous Optimizers) source code
# and tries to stay close to its methods/functions in order to allow the same
# type of comparisons in Julia.
module COCO

# BBOB test functions are optimization problems
abstract type BBOBFunction <: OptimizationProblem end

# Number of dimensions.
ndims(f::BBOBFunction) = f.ndims

instantiated(f::BBOBFunction) = ndims(f) != Nothing

# Returns the best function value of this instance of the function.
getfopt(f::BBOBFunction) = instantiated(f) ? f.fopt : Exception("This function class has not been instantiated yet.")

# Set/get the argument of the optimum of the function.
setxopt(f::BBOBFunction, xopt) = f.xopt = xopt
getxopt(f::BBOBFunction) = f.xopt

# Returns the objective function value.
evaluate(f::BBOBFunction, x) = evalfull(f, x)[0]

# Returns a float penalty for being outside of boundaries [-5, 5]
function defaultboundaryhandling(x, fac = 1.0)
    xoutside = maximum(vcat(zeros(size(x)), abs(x) - 5), 1) .* sign(x)
    return fac * sum(xoutside.^2)
end

# Most functions do not have a penalty outside the boundaries though so return
# a zero penalty by default.
boundaryhandling(f::BBOBFunction, x) = 0

function evalfull(f::BBOBFunction, x)
    fadd = getfopt(f)
    curshape, dim = size(x)
    # it is assumed x are row vectors

    if f.lastshape != curshape
      initwithsize(curshape, dim)
    end

    # BOUNDARY HANDLING
    fadd += boundaryhandling(f, x)

    # TRANSFORMATION IN SEARCH SPACE
    x = x - arrxopt(f) # cannot be replaced with x -= arrxopt! RF: Why???

    # COMPUTATION core
    ftrue = compute_core(f, x)
    fval = noise(f.noisefunc, ftrue)

    # FINALIZE
    ftrue += fadd
    fval += fadd
    return (fval, ftrue)
end

# All BBOBFunction functions must implement evalfull, which return noisy and
# noise-free value, the latter for recording purpose.
# evalfull(f::BBOBFunction, x) = ...

# Boundary handling
# boundaryhandling(f::BBOBFunction, x) = ...

# We depart from the BBOB type/class hierarchy they use in Python since
# we can use Julia's type system to model the noise functions in a better
# way than how it is done in Python. Thus a noise function need not be a
# BBOBFunction.
abstract type BBOBNoiseFunction end

# Different types depending on the noise function used.
struct BBOBNfreeFunction <: BBOBNoiseFunction
end
noise(f::BBOBNfreeFunction, ftrue) = copy(ftrue) # no noise added

struct BBOBGaussFunction <: BBOBNoiseFunction
    gaussbeta::Float64
end
noise(f::BBOBGaussFunction, ftrue) = fGauss(ftrue, f.gaussbeta)

struct BBOBUniformFunction <: BBOBNoiseFunction
    unifalphafac::Float64
    unifbeta::Float64
end
noise(f::BBOBUniformFunction, ftrue) =
    fUniform(ftrue, f.unifalphafac * (0.49 + 1. / f.dim), f.unifbeta)

struct BBOBCauchyFunction <: BBOBNoiseFunction
    cauchyalpha::Float64
    cauchyp::Float64
end
noise(f::BBOBCauchyFunction, ftrue) = fCauchy(ftrue, f.cauchyalpha, f.cauchyp)

abstract type FSphere{NoiseFunc} <: BBOBFunction end
compute_core(f::FSphere, x) = sum(x.^2)

# Sphere without noise
struct F1 <: FSphere{BBOBNfreeFunction}
    funId::Int
    noisefunc::BBOBNoiseFunction
    F1() = new(1, BBOBNfreeFunction())
end
boundaryhandling(f::F1, x) = 0

# Sphere with Gaussian noise
struct F101 <: FSphere{BBOBGaussFunction}
    funId::Int
    noisefunc::BBOBNoiseFunction
    F101 = new(101, BBOBGaussFunction(0.01))
end

# Run an optimizer like in the COCO (COmparing Continuous Optimizers) sw,
# i.e. with a given number of dimensions, specific target value and max number
# of function evaluations.
# COCO method: MY_OPTIMIZER in MY_OPTIMIZER.m
# COCO differences:
#  - We send in an optimizer instead of changing the code each time we want to evaluate an optimizer.
#function coco_run_optimizer(optimizer, problem, numDimensions,
#  funcTarget, maxFunevals; populationSize = 200)
#  maxfunevals = min(1e8 * numDimensions, maxFunevals)
#  popsize = min(maxFunevals, populationSize)
#  fbest = inf
#  for iter in 1:ceil(maxfunevals/popsize)
#
#    xpop = 10 * rand(numDimensions, popsize) - 5;      % new solutions
#    [fvalues, idx] = sort(feval(FUN, xpop)); % evaluate
#    if fbest > fvalues(1)                    % keep best
#      fbest = fvalues(1);
#      xbest = xpop(:,1);
#    end
#    if feval(FUN, 'fbest') < ftarget         % COCO-task achieved
#      break;                                 % (works also for noisy functions)
#    end
#  end
#end
#
#function coco


text:
using SafeTestsets
using Test
using TraceEstimators

const METHOD = get(ENV, "METHOD", "all")

# SLQ
if METHOD in ("all", "slq")
    println("Running tests for SLQ")
    @safetestset "slq.jl" begin
        include("slq.jl")
    end
end

# cheby
if METHOD in ("all", "chebyhutch")
    println("Running tests for ChebyHutch")
    @safetestset "chebyhutch.jl" begin
        include("chebyhutch.jl")
    end
end

# diagonalapprox
if METHOD in ("all", "diagonalapprox")
    println("Running tests for DiagonalApprox")
    @safetestset "diagonalapprox.jl" begin
        include("diagonalapprox.jl")
    end
end

# Hutchinson
if METHOD in ("all", "hutchinson")
    println("Running tests for Hutchinson")
    @safetestset "hutchinson.jl" begin
        include("hutchinson.jl")
    end
end

#Diagonal Approximation
if METHOD in ("all", "diagapp")
    println("Running tests for diagapp")
    @safetestset "diagapp.jl" begin
        include("diagapp.jl")
    end
end


text:
using StructArrays

# A basic chains constructor that works with the Transition struct we defined.
function bundle_samples(
    rng::AbstractRNG,
    model::DensityModel,
    s::Metropolis,
    N::Integer,
    ts::Vector{T},
    chain_type::Type{StructArray};
    param_names=missing,
    kwargs...
) where {ModelType<:AbstractModel, T<:AbstractTransition}
    return StructArray(bundle_samples(rng, model, s, N, ts, Vector{NamedTuple}; param_names=param_names, kwargs...))
end

AbstractMCMC.chainscat(cs::StructArray...) = [cs[i] for i in 1:length(cs)]


text:
export MeasureResult, gensample, QPEPSRunTime, QPEPSMachine
export chbasis!

struct QPEPSRunTime
    circuit
    rotblocks
    mbasis
    mblocks
end

function chbasis!(rt::QPEPSRunTime, basis)
    for m in rt.mbasis
        setcontent!(m, chcontent(m.content, basis_rotor(basis)))
    end
    return rt
end

struct QPEPSMachine{CT<:AbstractQPEPSConfig}
    config::CT
    runtime::QPEPSRunTime
    reg0::AbstractRegister
end

function QPEPSMachine(config::AbstractQPEPSConfig, reg0::AbstractRegister)
    c = get_circuit(config)
    rt = QPEPSRunTime(c, collect_blocks(RotationGate, c), collect_blocks(Bag, c), collect_blocks(Measure, c))
    QPEPSMachine(config, rt, reg0)
end

struct MeasureResult{AT<:AbstractMatrix}
    nm::Int
    nbatch::Int
    nx::Int
    data::AT
end
MeasureResult(nm::Int, data::AbstractMatrix) = MeasureResult(nm, size(data)..., data)

Base.getindex(m::MeasureResult, i::Int,j::Int) = readbit.(view(m.data,:,i), j)
Base.getindex(m::MeasureResult, ij) = getindex(m, (ij-1)%m.nx+1, (ij-1)÷m.nx+1)
Base.getindex(m::MeasureResult, ci::CartesianIndex) = getindex(m, ci.I...)
Base.size(m::MeasureResult) = (m.nx, m.nm)
Base.size(m::MeasureResult,i::Int) = i==1 ? m.nx : m.nm

"""
generate samples
"""
function gensample(qpeps::QPEPSMachine, basis)
    c, rt = qpeps.config, qpeps.runtime
    nx = c.nrepeat+nbath(c)
    chbasis!(rt, basis)
    copy(qpeps.reg0) |> rt.circuit
    res = zeros(Int, nbatch(qpeps.reg0), nx)

    for j=1:nx
        res[:,j] = Vector(rt.mblocks[j].results)
    end
    return MeasureResult(nm(c), res)
end