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

using ArgParse

function parse_cmd()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--snofile"
      help = "snoRNA fasta file"
      arg_type = ASCIIString
    "--regex"
      help = "snoRNA fasta header regex with capture [optional]"
      arg_type = ASCIIString
      default = "^>\\s*.*_.*_(\\S+)_.*"
  end
  return parse_args(s)
end

##################################################################
function lenientregex( s::ASCIIString )
   levels = ASCIIString[]
   nlevels = 6
   for i in 1:nlevels
      push!(levels, s)  # initialize levels
   end
   # this loop requires the first and last bases to match exactly
   for c in 2:(length(s)-1)
      pre = s[1:(c-1)]
      post = s[(c+1):end]
      levels[2] *= "|" * pre * "[atgcuATGCU]" * post
      levels[4] *= "|" * pre * "." * string(s[c]) * post
      levels[6] *= "|" * pre * "." * "[atgcuATGCU]" * post
   end
   # this one allows first and last bases to be lenient
   for c in 1:length(s)
      pre = s[1:(c-1)]
      post = s[(c+1):end]
      levels[3] *= "|" * pre * "[atgcuATGCU]" * post
      levels[5] *= "|" * pre * "." * string(s[c]) * post
   end
   ret = Regex[]
   for i in 1:nlevels
     push!(ret, Regex(levels[i]))
   end
   ret
end #-->Array{Regex,1}

function searchregex( regexarray::Array{Regex,1}, s::ASCIIString )
   for r in regexarray
      m = match(r, s)
      if m != nothing
        return (m.offset,m.match)
      end
   end
   (0,"")
end #-->Tuple{Int64,ASCIIString}

function annotate_cdbox( sno::ASCIIString )
   cboxreg = lenientregex("TGATGA")
   dboxreg = lenientregex("CTGA")

   cboxseg = sno[1:20]
   cbox = searchregex(cboxreg, cboxseg)
   #@assert(cbox[1] > 0, "Cannot find cbox in $cboxseg!")

   dboxseg = sno[(end-20):end]
   (dboxpos,dboxmatch) = searchregex(dboxreg, dboxseg)
   #@assert(dboxpos > 0, "Cannot find dbox in $dboxseg!")
   dboxpos += length(sno) - 20
   dbox = (dboxpos, dboxmatch)

   dprimeseg = sno[(cbox[1]+16):(dboxpos-14)]
   (dprimepos, dprimematch) = searchregex(dboxreg, dprimeseg)
   dprimepos += cbox[1] + 16
   dprime = (dprimepos, dprimematch)

   println(STDERR, "$(cbox[2]) ... $(dprime[2]) ... $(dbox[2])")

   (cbox, dprime, dbox)
end #-->Tuple{Tuple}

# this function reads a fasta file/io and returns a dict of {name,seq}
function readfasta( io; regex = r">\s*(\S+)" )
   rethash = Dict{ASCIIString,ASCIIString}()

   #internal funct for fasta header
   function checkname( head::ASCIIString, namereg::Regex )
      res = match(namereg, head)
      @assert(res != nothing, "$head looks to be an incorrectly formated fasta header!")
      chomp(res.captures[1])
   end

   curseq = ""
   head = readline(io)
   curname = checkname( head, regex )
   for line::ASCIIString in eachline( io )
      #finish up
      if line[1] == '>' #header line
         # push to dict
         rethash[curname] = curseq
         curname = checkname( line, regex )
         curseq = ""
      else #sequence line
        curseq *= chomp( uppercase(line) )
      end
   end
   rethash
end #--> Dict{ASCIIString,ASCIIString}

# find the distance from the 5' end of the transcript of the antisense region
function bind_distance( ind::Int64, ligstruct::String, startpos::Int64 )
   antiregex = r"([ATGCU]+(?:[\.\(\)]{1,5}[ATGCU]+)?)"
   m = match( antiregex, ligstruct ) # match antisense site
   cut = match( r"(\|)", ligstruct )
   m == nothing && return (1,"",1)
   cap = m.captures[1]
   len = length(cap)
   offset,lig = 0,0
   if ind <= 1
      offset = startpos + m.offsets[1]
      lig = startpos + (cut == nothing ? length(ligstruct) : cut.offsets[1])
   else
      barpos = cut == nothing ? 0 : cut.offsets[1]
      offset = startpos - barpos + m.offsets[1]
      lig = startpos
   end
   (offset, cap, lig)
end #--> Tuple{Int64,Int64,Int64}

function print_heatrow{I <: Integer}( io, rowname::ASCIIString, offset::I, antiseq::String, ligpos::I, anchorpos::I, totallen::I, boxes::Tuple; nais="NA", order="" )
   println(STDERR, "$io\t$rowname\t$offset\t$antiseq\t$ligpos\t$anchorpos")
   prestr = rowname * "\t"
   str = repeat( "0", totallen )
   for i = 1:3
     (off, seq) = boxes[i]
     off > 0 || continue
     asnum = replace( seq, r"[AUCGT]", string(i+2) )
     str = str[1:off-1] * asnum * str[off+length(asnum):end]
   end

   antinum = replace( antiseq, r"[AUCGT]", "1" ) |> x->replace( x, r"\.", "0" )
   str = str[1:offset-1] * antinum * str[offset+length(antinum):end]

   str = ligpos > totallen ? str : str[1:ligpos-1] * "2" * str[ligpos+1:end]
   rev = reverse(str)
   pad = totallen - anchorpos
   outiter = [ pad:-1:0; -1:-1:(-1*anchorpos)]
   for i = 1:length(rev)
      cur = rev[i]
      println( io, prestr * string(outiter[i]) * order * "\t$cur" )
   end
end

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

   fastareg = pargs["regex"] == nothing ? r"^>\s*(\S+)" : Regex(pargs["regex"])

   @assert(pargs["snofile"] != nothing, "you have to provide --snofile")
   snofile = pargs["snofile"]
   fh = open( snofile, "r" )
   snodict = readfasta( fh, regex=fastareg )
   close(fh)

   # try to annotate the C/D boxes
   cdboxhash = Dict{ASCIIString, Tuple}()
   for k in keys(snodict)
      cdboxhash[k] = annotate_cdbox( snodict[k] )
   end

   const genInd = 25

   output = Dict{ASCIIString, Tuple}()
   ligdist = Dict{ASCIIString, Int64}()
   capsize = Dict{ASCIIString, Int64}()
   annobox = Dict{ASCIIString, Tuple}()

   # now lets go through the pvl/lig data
   for l in eachline(STDIN)
      s = split( l, '\t' )
      @assert( ismatch(r"SNORD", s[genInd]), "$(s[genInd]) doesn't appear to have a snoRNA in it? (SNORD)\n" )
      ids = map(x->chomp(x), split( s[genInd], ':' ))
      ind = ismatch(r"SNORD", ids[1]) ? 1 : 2
      @assert( ismatch(r"SNORD", ids[ind]) )
      struct = s[ind+18]
      start = parse(Int, split( s[14], ',' )[ind])
      (antioffset, cap, lig) = bind_distance( ind, struct, start )
      @assert( haskey(cdboxhash, "SNORD9" ), "$(ids[ind]) doesn't have a key!" )
      cdboxes = cdboxhash[ ids[ind] ]
      snolen  = length( snodict[ ids[ind] ] )
      curname_key = ids[ind] * "_" * ids[ ind == 1 ? 2 : 1 ]
      if antioffset <= cdboxes[2][1] # D' box offset
         #curname_key = ids[ind] * "_DpBOX"
         outtup = (curname_key, antioffset, cap, lig, cdboxes[2][1], snolen )
         outdist = ( cdboxes[2][1] - (round(Int,length(cap) / 2) + antioffset) )
      elseif antioffset <= cdboxes[3][1] # D box offset
         #curname_key = ids[ind] * "_DBOX"
         outtup = (curname_key, antioffset, cap, lig, cdboxes[3][1], snolen )
         outdist = ( cdboxes[3][1] - (round(Int,length(cap) / 2) + antioffset) )
      else
         continue
      end

      # add to output dicts
      storedlen = get( capsize, curname_key, 0 )
      if length(cap) > storedlen
         output[curname_key]  = outtup # store anno tuple
         ligdist[curname_key] = outdist # store distance to nearest D box
         capsize[curname_key] = length(cap) # store length of antisense
         annobox[curname_key] = cdboxes #store tuple of c/d box annotations
      end
   end

   # print in order
   keysarr = collect( keys( ligdist ) )
   i = 1 # order iterator
   for ord in sortperm( collect( values( ligdist ) ) )
      outkey = keysarr[ord]
      print_heatrow( STDOUT, output[outkey]..., annobox[outkey], order="\t" * string(i) )
      i += 1
   end
end
###################################################################
main()

max_forks_repo_name: timbitz/aligater


text:
"""
ItoDiffusionProcess

Type for storing an Itô diffusion process.

It consists of a function μ(x) for the drift term, and a function σ(x) for the diffusion term.

See also https://en.wikipedia.org/wiki/It%C3%B4_diffusion
"""
struct ItoDiffusionProcess
	μ::Function
	σ::Function
	# interval::Tuple{T, T} where T<:Real
	# reflective::Tuple{Bool, Bool}

	# function ItoDiffusionProcess(mu::Function, sigma::Function, interval::Tuple{T, T} where T<:Real)

	# 	# If reflective nature of barriers is not specified, assume true
	# 	reflective = (True, True)

	# 	# Allocate object
	# 	new(mu, sigma, interval, reflective)
	# end

end

ext: jl


text:
using Unicode: graphemes

const TEST_GRAPHEMES = true

function myreverse(phrase::String)
    join(reverse(collect(graphemes(phrase))))
end


text:
function rhs!(du::AbstractArray{<:Any,4}, u, t,
              mesh::CurvedMesh, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, mesh, equations,
                                                                dg.volume_integral, dg, cache)

  # Calculate interface fluxes
  @timeit_debug timer() "interface flux" calc_interface_flux!(u, mesh, equations, dg, cache)

  # Calculate boundary fluxes
  @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, u, t, boundary_conditions, equations, mesh, dg)

  # Calculate surface integrals
  @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  @timeit_debug timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache)

  # Calculate source terms
  @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, mesh::CurvedMesh, equations,
                               volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis
  @unpack contravariant_vectors = cache.elements

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, element)

      flux1 = flux(u_node, 1, equations)
      flux2 = flux(u_node, 2, equations)

      # Compute the contravariant flux by taking the scalar product of the
      # first contravariant vector Ja^1 and the flux vector
      contravariant_flux1 = (contravariant_vectors[1, 1, i, j, element] * flux1 +
                             contravariant_vectors[1, 2, i, j, element] * flux2)

      for ii in eachnode(dg)
        integral_contribution = derivative_dhat[ii, i] * contravariant_flux1
        add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
      end

      # Compute the contravariant flux by taking the scalar product of the
      # second contravariant vector Ja^2 and the flux vector
      contravariant_flux2 = (contravariant_vectors[2, 1, i, j, element] * flux1 +
                             contravariant_vectors[2, 2, i, j, element] * flux2)

      for jj in eachnode(dg)
        integral_contribution = derivative_dhat[jj, j] * contravariant_flux2
        add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
      end
    end
  end

  return nothing
end

function calc_interface_flux!(u::AbstractArray{<:Any,4}, mesh::CurvedMesh{2},
                              equations, dg::DG, cache)
  @unpack elements = cache

  @threaded for element in eachelement(dg, cache)
    # Interfaces in negative directions
    # Faster version of "for orientation in (1, 2)"

    # Interfaces in x-direction (`orientation` = 1)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[1, element],
                         element, 1, u, equations, dg, cache)

    # Interfaces in y-direction (`orientation` = 2)
    calc_interface_flux!(elements.surface_flux_values,
                         elements.left_neighbors[2, element],
                         element, 2, u, equations, dg, cache)
  end

  return nothing
end

@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
                                      orientation, u, equations, dg::DG, cache)
  # This is slow for LSA, but for some reason faster for Euler (see #519)
  if left_element <= 0 # left_element = 0 at boundaries
    return surface_flux_values
  end

  @unpack surface_flux = dg
  @unpack contravariant_vectors = cache.elements

  right_direction = 2 * orientation
  left_direction = right_direction - 1

  for i in eachnode(dg)
    if orientation == 1
      u_ll = get_node_vars(u, equations, dg, nnodes(dg), i, left_element)
      u_rr = get_node_vars(u, equations, dg, 1,          i, right_element)

      # First contravariant vector Ja^1 as SVector
      normal_vector = get_contravariant_vector(1, contravariant_vectors, 1, i, right_element)
    else # orientation == 2
      u_ll = get_node_vars(u, equations, dg, i, nnodes(dg), left_element)
      u_rr = get_node_vars(u, equations, dg, i, 1,          right_element)

      # Second contravariant vector Ja^2 as SVector
      normal_vector = get_contravariant_vector(2, contravariant_vectors, i, 1, right_element)
    end

    flux = surface_flux(u_ll, u_rr, normal_vector, equations)

    for v in eachvariable(equations)
      surface_flux_values[v, i, right_direction, left_element] = flux[v]
      surface_flux_values[v, i, left_direction, right_element] = flux[v]
    end
  end

  return surface_flux_values
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
                             equations::AbstractEquations{2}, mesh::CurvedMesh{2}, dg::DG)
  @assert isperiodic(mesh)
end

function calc_boundary_flux!(cache, u, t, boundary_condition,
                             equations::AbstractEquations{2}, mesh::CurvedMesh{2}, dg::DG)
  @unpack surface_flux = dg
  @unpack surface_flux_values = cache.elements
  linear_indices = LinearIndices(size(mesh))

  for cell_y in axes(mesh, 2)
    # Negative x-direction
    direction = 1
    element = linear_indices[begin, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_condition, equations, mesh, dg, cache,
                                       direction, (1, j), (j,), element)
    end

    # Positive x-direction
    direction = 2
    element = linear_indices[end, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_condition, equations, mesh, dg, cache,
                                       direction, (nnodes(dg), j), (j,), element)
    end
  end

  for cell_x in axes(mesh, 1)
    # Negative y-direction
    direction = 3
    element = linear_indices[cell_x, begin]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_condition, equations, mesh, dg, cache,
                                       direction, (i, 1), (i,), element)
    end

    # Positive y-direction
    direction = 4
    element = linear_indices[cell_x, end]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_condition, equations, mesh, dg, cache,
                                       direction, (i, nnodes(dg)), (i,), element)
    end
  end
end

function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple},
                             equations::AbstractEquations{2}, mesh::CurvedMesh{2}, dg::DG)
  @unpack surface_flux = dg
  @unpack surface_flux_values = cache.elements
  linear_indices = LinearIndices(size(mesh))

  for cell_y in axes(mesh, 2)
    # Negative x-direction
    direction = 1
    element = linear_indices[begin, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction], equations, mesh, dg, cache,
                                       direction, (1, j), (j,), element)
    end

    # Positive x-direction
    direction = 2
    element = linear_indices[end, cell_y]

    for j in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1,
                                       boundary_conditions[direction], equations, mesh, dg, cache,
                                       direction, (nnodes(dg), j), (j,), element)
    end
  end

  for cell_x in axes(mesh, 1)
    # Negative y-direction
    direction = 3
    element = linear_indices[cell_x, begin]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction], equations, mesh, dg, cache,
                                       direction, (i, 1), (i,), element)
    end

    # Positive y-direction
    direction = 4
    element = linear_indices[cell_x, end]

    for i in eachnode(dg)
      calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2,
                                       boundary_conditions[direction], equations, mesh, dg, cache,
                                       direction, (i, nnodes(dg)), (i,), element)
    end
  end
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, u, t, orientation,
                                          boundary_condition::BoundaryConditionPeriodic, equations, mesh::CurvedMesh,
                                          dg::DG, cache, direction, node_indices, surface_node_indices, element)
  @assert isperiodic(mesh, orientation)
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, u, t, orientation,
                                          boundary_condition, equations, mesh::CurvedMesh, dg::DG, cache,
                                          direction, node_indices, surface_node_indices, element)
  @unpack node_coordinates, contravariant_vectors= cache.elements
  @unpack surface_flux = dg

  u_inner = get_node_vars(u, equations, dg, node_indices..., element)
  x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)

  # Contravariant vector Ja^i is the normal vector
  normal_vector = get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element)
  flux = boundary_condition(u_inner, normal_vector, direction, x, t, surface_flux, equations)

  for v in eachvariable(equations)
    surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
  end
end

function apply_jacobian!(du::AbstractArray{<:Any,4}, mesh::CurvedMesh, equations, dg::DG, cache)
  @unpack inverse_jacobian = cache.elements

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      factor = -inverse_jacobian[i, j, element]

      for v in eachvariable(equations)
        du[v, i, j, element] *= factor
      end
    end
  end

  return nothing
end


text:
# Implement number base conversion

function all_your_base(digits, base_in, base_out)

end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using DateUtils
using Compat.Test, Compat.Dates, Missings

@testset "DateUtils" begin
  include("utils.jl")
end


text:
using Documenter
using Literate
using RHEOS

"""
    docprepare()

Convert all source documentation (.jl and .md files)
into markdown, ready for documentation build, and
relevant .jl files into Jupyter notebooks ready for
deployment to notebooks branch.
"""
function docprepare()
    # remove remnants of previous build and recreate staging dir
    rm("docs/staging-docs", force=true, recursive=true)
    mkdir("docs/staging-docs")

    # copy readme to staging-docs, remove Logo image
    write("docs/staging-docs/index.md",
            open("README.md") do input
                readuntil(input, "<!-- delim -->", keep = true)
                read(input)
            end)

    # copy assets to staging directory
    cp("docs/src/assets", "docs/staging-docs/assets")

    # iterate through src and convert/copy as appropriate
    for file in readdir("docs/src")
        if endswith(file, "md")
            cp("docs/src/$file", "docs/staging-docs/$file")
        elseif endswith(file, "jl")
            Literate.markdown("docs/src/$file", "docs/staging-docs/"; documenter=true)
        end
    end
end

function notebookprepare()
    # create notebook staging dir
    mkdir("docs/staging-docs/notebooks")

    # copy assets to notebooks staging directory
    cp("docs/src/assets", "docs/staging-docs/notebooks/assets")

    for file in readdir("docs/src")
        if endswith(file, "jl")
            Literate.notebook("docs/src/$file", "docs/staging-docs/notebooks/")
        end
    end
end

function maindocbuilder()
    # prepare doc files from source
    docprepare()

    # prepare notebook files from source
    notebookprepare()

    # build docs from staging area
    makedocs(modules=[RHEOS],
            doctest = false, clean=true,
            format = Documenter.HTML(),
            sitename ="RHEOS.jl",
            source = "staging-docs",
            authors = "J Louis Kaplan, Alessandra Bonfati, Alexandre Kabla",
            pages = ["Home" => "index.md",
                     "Architecture" => "architecture.md",
                     "File I/O" => "fileio.md",
                     "Preprocessing" => "preprocessing.md",
                     "Generating Data" => "gendata.md",
                     "Fit and Predict" => ["Time data" => "fitpredictTime.md",
                                           "Frequency data" => "fitpredictFreq.md"],
                     "Numerical Aspects" => "numerical.md",
                     "Models" => ["Basic Elements" => "model_elements.md",
                                  "Maxwell" => "model_maxwell.md",
                                  "Kelvin-Voigt" => "model_kv.md",
                                  "Zener" => "model_zener.md",
                                  "Poynting-Thomson" => "model_pt.md",
                                  "Burgers" => "model_burgers.md",
                                  "Create Your Model" => "model_create.md"],
                     "Additional Examples" => "examples.md",
                     "API" => "API.md"])

    deploydocs(repo = "github.com/JuliaRheology/RHEOS.jl.git",
               deps = nothing,
               make = nothing,
               target = "build")

end

maindocbuilder()


text:
# Extensions to the C API
#
# These functions will only raise exceptions arising from julia. Python exceptions are propagated by return value as usual.
#
# Functions `CPy<type>_From(x)` convert `x` to a python `type`. If `x` is a `PyPtr`, this is equivalent to `<type>(x)` in python. Otherwise, it typically is a wrapper for one of the `CPy<type>_From<ctype>` functions in the C API. Returns `C_NULL` on error.
#
# Functions `CPy<type>_As(T, x)` convert `x` (which should be a `PyPtr` of python type `type`) to a `T`. Returns `nothing` on error.
#
# `CPyObject_From(x)` defines the default conversion from julia to python.

"""
	CPyDict_FromIterator([fk=identity, fv=identity,] x)

A python `dict` whose keys are `fk(k)` and values are `fv(v)` for `(k,v) ∈ x`.

The keys and values must be `PyPtr`. This steals references, so the elements must be new references.
"""
function CPyDict_FromIterator(fk, fv, x) :: PyPtr
	t = CPyDict_New()
	t == C_NULL && (return C_NULL)
	for (k, v) in x
		kk = fk(k) :: PyPtr
		kk == C_NULL && (CPy_DecRef(t); return C_NULL)
		vv = fv(v) :: PyPtr
		vv == C_NULL && (CPy_DecRef(kk); CPy_DecRef(t); return C_NULL)
		e = CPyDict_SetItem(t, kk, vv)
		CPy_DecRef(kk)
		CPy_DecRef(vv)
		e == -1 && (CPy_DecRef(t); return C_NULL)
	end
	t
end

CPyDict_FromIterator(x) =
	CPyDict_FromIterator(identity, identity, x)

function CPyDict_From(x::PyPtr) :: PyPtr
	CPyObject_CallFunction(CPyDict_Type[], x)
end

"""
    CPyTuple_FromIterator([f=identity,] x)

A python `tuple` whose elements come from `map(f, x)`.

The length of `x` must be known. The elements must be `PyPtr`. This steals references, so the elements must be new references.
"""
function CPyTuple_FromIterator(f, x) :: PyPtr
	t = CPyTuple_New(length(x))
	t == C_NULL && (return C_NULL)
	i = 0
	for y in x
	    z = f(y) :: PyPtr
	    z == C_NULL && (CPy_DecRef(t); return C_NULL)
	    e = CPyTuple_SetItem(t, i, z)
	    e == -1 && (CPy_DecRef(t); return C_NULL)
	    i += 1
	end
	t
end

CPyTuple_FromIterator(x) = CPyTuple_FromIterator(identity, x)

function CPyTuple_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyTuple_Type[], x)
end

"""
    CPyList_FromIterator([f=identity,] x)

A python `list` whose elements come from `map(f, x)`.

The elements must be `PyPtr`. This steals references, so the elements must be new references.
"""
function CPyList_FromIterator(f, x) :: PyPtr
	# potential optimization: check if the length of x is known in advance
	t = CPyList_New(0)
	t == C_NULL && (return C_NULL)
	for y in x
	    z = f(y) :: PyPtr
	    z == C_NULL && (CPy_DecRef(t); return C_NULL)
	    e = CPyList_Append(t, z)
	    e == -1 && (CPy_DecRef(t); return C_NULL)
	end
	t
end

CPyList_FromIterator(x) = CPyList_FromIterator(identity, x)

function CPyList_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyList_Type[], x)
end

function CPyUnicode_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyUnicode_Type[], x)
end

function CPyUnicode_From(x::String)::PyPtr
	CPyUnicode_DecodeUTF8(Base.unsafe_convert(Ptr{UInt8}, x), sizeof(x), C_NULL)
end

function CPyUnicode_From(x::AbstractString)::PyPtr
	CPyUnicode_From(convert(String, x))
end

function CPyUnicode_As(::Type{T}, x) where {T<:AbstractString}
	b = CPyUnicode_AsUTF8String(x)
	b == C_NULL && (return nothing)
	buf = Ref{Ptr{UInt8}}(C_NULL)
	len = Ref{Cssize_t}(0)
	z = CPyBytes_AsStringAndSize(b, buf, len)
	z == -1 && (CPy_DecRef(b); return nothing)
	s = unsafe_string(buf[], len[])
	CPy_DecRef(b)
	convert(T, s)
end

function CPyBool_From(x::PyPtr)::PyPtr
	r = CPyObject_IsTrue(x)
	r == -1 && (return C_NULL)
	CPyBool_From(r == 1)
end

function CPyBool_From(x::Bool)::PyPtr
	r = x ? CPy_True[] : CPy_False[]
	CPy_IncRef(r)
	r
end

function CPyLong_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyLong_Type[], x)
end

function CPyLong_From(x::AbstractString)::PyPtr
	y = CPyUnicode_From(x)
	y == C_NULL && (return C_NULL)
	z = CPyLong_From(y)
	CPy_DecRef(y)
	z
end

function CPyLong_From(x::Real)::PyPtr
	y = CPyFloat_From(x)
	y == C_NULL && (return C_NULL)
	z = CPyLong_From(y)
	CPy_DecRef(y)
	z
end

function CPyLong_From(x::Complex)::PyPtr
	y = CPyComplex_From(x)
	y == C_NULL && (return C_NULL)
	z = CPyLong_From(y)
	CPy_DecRef(y)
	z
end

function CPyLong_From(x::Number)::PyPtr
	CPyLong_From(convert(Real, x)::Real)
end

function CPyLong_From(x)::PyPtr
	CPyLong_From(convert(Number, x)::Number)
end

@static if pyversion < v"3"
	function CPyLong_From(x::T) where {T<:Integer}
		if isbitstype(T)
			if T <: Unsigned
				if sizeof(T) ≤ sizeof(Csize_t)
					return CPyInt_FromSize_t(x)
				end
			else
				if sizeof(T) ≤ sizeof(Cssize_t)
					return CPyInt_FromSsize_t(x)
				end
			end
		end
		CPyLong_From(string(x))
	end
else
	function CPyLong_From(x::T) where {T<:Integer}
		if isbitstype(T)
			if T <: Unsigned
				if sizeof(T) ≤ sizeof(Culonglong)
					return CPyLong_FromUnsignedLongLong(x)
				end
			else
				if sizeof(T) ≤ sizeof(Clonglong)
					return CPyLong_FromLongLong(x)
				end
			end
		end
		CPyLong_From(string(x))
	end
end

@static if pyversion < v"3"
    function CPyLong_As(::Type{T}, o) where {T<:Integer}
    	val = CPyInt_AsSsize_t(o)
    	if val == -1 && CPyErr_Occurred() != C_NULL
    		CPyErr_Clear()
    		convert(T, convert(BigInt, o))
    	else
    		convert(T, val)
    	end
    end
elseif pyversion < v"3.2"
	function CPyLong_As(::Type{T}, o) where {T<:Integer}
		val = CPyLong_AsLongLong(o)
		if val == -1 && CPyErr_Occurred() != C_NULL
			CPyErr_Clear()
			convert(T, convert(BigInt, o))
		else
	        convert(T, val)
	    end
    end
else
    function CPyLong_As(::Type{T}, o) where {T<:Integer}
        overflow = Ref{Cint}()
        val = CPyLong_AsLongLongAndOverflow(o, overflow)
        val == -1 && CPyErr_Occurred() != C_NULL && (return nothing)
        if iszero(overflow[])
        	convert(T, val)
        else
        	convert(T, convert(BigInt, o))
        end
    end
end

function CPyLong_As(::Type{BigInt}, o)
	s = CPyObject_Str(String, o)
	s===nothing && return nothing
	parse(BigInt, s)
end

function CPyFloat_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyFloat_Type[], x)
end

function CPyFloat_From(x::Float64)::PyPtr
	CPyFloat_FromDouble(x)
end

function CPyFloat_From(x::Real)::PyPtr
	CPyFloat_From(convert(Float64, x))
end

function CPyFloat_From(x::Number)::PyPtr
	CPyFloat_From(convert(Real, x))
end

function CPyFloat_From(x)::PyPtr
	CPyFloat_From(convert(Number, x))
end

function CPyFloat_From(x::AbstractString)::PyPtr
	y = CPyUnicode_From(x)
	y == C_NULL && (return C_NULL)
	z = CPyFloat_From(y)
	CPy_DecRef(y)
	z
end

function CPyFloat_As(::Type{T}, x) where {T<:Number}
	y = CPyFloat_AsDouble(x)
	y == -1 && CPyErr_Occurred() != C_NULL && (return nothing)
	convert(T, y)
end

function CPyComplex_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyComplex_Type[], x)
end

function CPyComplex_From(x::Complex)::PyPtr
	CPyComplex_FromDoubles(real(x), imag(x))
end

function CPyComplex_From(x::Number)::PyPtr
	y = CPyFloat_From(x)
	y == C_NULL && (return C_NULL)
	z = CPyComplex_From(y)
	CPy_DecRef(y)
	z
end

function CPyComplex_As(::Type{T}, x) where {T<:Number}
	y = CPyComplex_AsCComplex(x)
	CPyErr_Occurred() == C_NULL || (return nothing)
	convert(T, Complex(y.real, y.imag))
end

function CPyBytes_From(x::PyPtr)::PyPtr
	CPyObject_CallFunction(CPyBytes_Type[], x)
end

function CPyBytes_From(x::AbstractVector{UInt8})::PyPtr
	# try to avoid copying x
	try
		if stride(x, 1) == 1
			ptr = Base.unsafe_convert(Ptr{UInt8}, x)
			sz = sizeof(x)
			return CPyBytes_FromStringAndSize(ptr, sz)
		end
	catch
	end
	y = convert(Vector, x)
	ptr = Base.unsafe_convert(Ptr{UInt8}, y)
	sz = sizeof(y)
	return CPyBytes_FromStringAndSize(ptr, sz)
end

function CPyObject_Repr(::Type{T}, x) where {T<:AbstractString}
	y = CPyObject_Repr(x)
	y == C_NULL && (return nothing)
	s = CPyUnicode_As(T, y)
	CPy_DecRef(y)
	s
end

function CPyObject_ASCII(::Type{T}, x) where {T<:AbstractString}
	y = CPyObject_ASCII(x)
	y == C_NULL && (return nothing)
	s = CPyUnicode_As(T, y)
	CPy_DecRef(y)
	s
end

function CPyObject_Str(::Type{T}, x) where {T<:AbstractString}
	y = CPyObject_Str(x)
	y == C_NULL && (return nothing)
	s = CPyUnicode_As(T, y)
	CPy_DecRef(y)
	s
end

function CPyObject_As end

function CPyObject_As(::Type{Nothing}, x)
	nothing
end

function CPyObject_As(::Type{Missing}, x)
	missing
end

function CPyObject_As(::Type{Bool}, x)
	r = CPyObject_IsTrue(x)
	r == -1 ? nothing : r == 1
end

function CPyObject_As(::Type{T}, x) where {T<:AbstractString}
	CPyObject_Str(T, x)
end

function CPyObject_As(::Type{T}, x) where {T<:Real}
	CPyFloat_As(T, x)
end

function CPyObject_From(x::PyPtr)::PyPtr
	x == C_NULL || CPy_IncRef(x)
	return x
end

function CPyObject_From(x::Integer)::PyPtr
	CPyLong_From(x)
end

function CPyObject_From(x::Real)::PyPtr
	CPyFloat_From(x)
end

function CPyObject_From(x::Complex)::PyPtr
	CPyComplex_From(x)
end

function CPyObject_From(x::Nothing)::PyPtr
	CPy_None_NewRef()
end

function CPyObject_From(x::Bool)::PyPtr
	CPyBool_From(x)
end

function CPyObject_From(x::AbstractString)::PyPtr
	CPyUnicode_From(x)
end

function CPyObject_From(x::Symbol)::PyPtr
	CPyUnicode_From(string(x))
end

function CPyObject_From(x::Tuple)::PyPtr
	CPyTuple_FromIterator(CPyObject_From, x)
end

function CPyObject_From(x::AbstractDict)::PyPtr
	CPyDict_FromIterator(CPyObject_From, CPyObject_From, pairs(x))
end

function CPyObject_From(x::AbstractVector)::PyPtr
	CPyList_FromIterator(CPyObject_From, x)
end

function CPyObject_From(x::AbstractSet)::PyPtr
	CPySet_FromIterator(CPyObject_From, x)
end

function CPyObject_CallFunction(f, args...; kwargs...)::PyPtr
	if isempty(kwargs)
		_args = CPyTuple_FromIterator(CPyObject_From, args)
		_args == C_NULL && (return C_NULL)
		r = CPyObject_CallObject(f, _args)
		CPy_DecRef(_args)
		return r
	else
		_args = CPyTuple_FromIterator(CPyObject_From, args)
		_args == C_NULL && (return C_NULL)
		_kwargs = CPyDict_FromIterator(CPyObject_From, CPyObject_From, kwargs)
		_kwargs == C_NULL && (CPy_DecRef(_args); return C_NULL)
		r = CPyObject_Call(f, _args, _kwargs)
		CPy_DecRef(_args)
		CPy_DecRef(_kwargs)
		return r
	end
end

function cpyargdata(arg)
	# parse out the default
	if arg isa Expr && arg.head == :(=)
		lhs, dflt = arg.args
		dflt = Some(dflt)
	else
		lhs = arg
		dflt = nothing
	end
	# parse out the type
	if lhs isa Expr && lhs.head == :(::)
		argname, typ = lhs.args
	else
		argname = lhs
		typ = :PyPtr
	end
	# check the argname
	argname isa Symbol || error("invalid argument: $arg")
	argname = argname==:_ ? nothing : argname
	# done
	(name=argname, typ=typ, dflt=dflt)
end

function cpyargsdata(args)
	argsdata = map(cpyargdata, args)
	# number of required arguments
	numreq = 0
	for arg in argsdata
		if arg.dflt===nothing
			numreq += 1
		else
			break
		end
	end
	# (minimum) number of positional arguments
	numpos = 0
	for arg in argsdata
		if arg.name == nothing
			numpos += 1
		else
			break
		end
	end
	#
	(args=argsdata, nreq=numreq, npos=numpos)
end

cpyargparse(::Type{PyPtr}, x::PyPtr) = (CPy_IncRef(x); x)
cpyargparse(::Type{T}, x::PyPtr) where {T<:Integer} =
	CPyLong_As(T, x)
cpyargparse(::Type{T}, x::PyPtr) where {T<:Real} =
	CPyFloat_As(T, x)
cpyargparse(::Type{T}, x::PyPtr) where {T<:AbstractString} =
	CPyUnicode_As(T, x)

cpyargfree(x) = nothing
cpyargfree(x::PyPtr) = CPy_DecRef(x)

function _cpyargsparse(t, k, args)
	data = cpyargsdata(args)
	tnames = [Symbol(:t, i-1) for i in 1:length(args)]
	anames = [Symbol(:x, i-1) for i in 1:length(args)]
	body = quote
		# check the length
		len = CPyTuple_Size(t)
		len == -1 && @goto(ERROR)
		len ≤ $(length(args)) || @goto(ERROR) # TODO: set an error
		# default PyPtr values (these are borrowed references)
		$([:($n :: PyPtr = C_NULL) for n in tnames]...)
		# parse the tuple
		$([:(len ≥ $i && ($n = CPyTuple_GetItem(t, $(i-1)); $n == C_NULL && @goto(ERROR))) for (i,(a,n)) in enumerate(zip(data.args, tnames))]...)
		# parse the kwargs
		$(k===nothing ? nothing : error("parsing keywords not implemented"))
		# extra parsing
		$([quote
			$n =
				if $t == C_NULL
					$(a.dflt===nothing ? :(#=TODO: set an error=# nothing) : :(Some($(esc(something(a.dflt))))))
				else
					cpyargparse($(esc(a.typ)), $t)
				end
			$n === nothing && @goto($(Symbol(:ERROR_,n)))
		end for (a,n,t) in zip(data.args, anames, tnames)]...)
		# return the tuple
		return ($([:(something($n)) for n in anames]...),)
		# errors
		$([:(cpyargfree($n); @label($(Symbol(:ERROR_,n)))) for n in reverse(anames)]...)
		@label ERROR
		return nothing
	end
	:((function (t, k); $body; end)($(esc(t)), $(k===nothing ? C_NULL : esc(k))))
end

"""
	@CPyArg_ParseTuple(args, NAME::TYPE=DEFAULT)

Similar to `@CPyArg_ParseTupleAndKeywords` but only takes an `args` tuple.

Note that the `NAME` is ignored in this case.
"""
macro CPyArg_ParseTuple(t, args...)
	_cpyargsparse(t, nothing, args)
end

"""
	@CPyArg_ParseTupleAndKeywords(args, kwargs, NAME::TYPE=DEFAULT, ...)

Parse `args` (a `PyPtr` to a tuple) and `kwargs` (a `PyPtr` to a dict) according to the argument specifiers of the form `NAME::TYPE=DEFAULT`. Return a tuple whose entries are of types `TYPE`, or `nothing` on error.

The `NAME` may be `_` to signify a positional-only argument. Otherwise, it is the name of a possibly-keyword argument.

The `TYPE` is optional, and defaults to `PyPtr`. Currently it can be one of `PyPtr`, `<:Integer`, `<:Real`, `<:AbstractString`.

The `DEFAULT` is optional. When not specified, the argument is required.
"""
macro CPyArg_ParseTupleAndKeywords(t, k, args...)
	_cpyargsparse(t, k, args...)
end


text:
using Test, FFTW, BasisFunctions, LinearAlgebra, SparseArrays

# Verify types of FFT and DCT plans by FFTW
# If anything changes here, the aliases in fouriertransforms.jl have to change as well
d1 = FFTW.plan_fft!(zeros(Complex{Float64}, 10), 1:1)
@test typeof(d1) == FFTW.cFFTWPlan{Complex{Float64},-1,true,1,UnitRange{Int64}}
d2 = FFTW.plan_fft!(zeros(Complex{Float64}, 10, 10), 1:2)
@test typeof(d2) == FFTW.cFFTWPlan{Complex{Float64},-1,true,2,UnitRange{Int64}}
d3 = FFTW.plan_bfft!(zeros(Complex{Float64}, 10), 1:1)
@test typeof(d3) == FFTW.cFFTWPlan{Complex{Float64},1,true,1,UnitRange{Int64}}
d4 = FFTW.plan_bfft!(zeros(Complex{Float64}, 10, 10), 1:2)
@test typeof(d4) == FFTW.cFFTWPlan{Complex{Float64},1,true,2,UnitRange{Int64}}

d5 = FFTW.plan_dct!(zeros(10), 1:1)
@test typeof(d5) == FFTW.DCTPlan{Float64,5,true}
d6 = FFTW.plan_idct!(zeros(10), 1:1)
@test typeof(d6) == FFTW.DCTPlan{Float64,4,true}

D = Diagonal(1:400)
R = BasisFunctions.RestrictionIndexMatrix{Float64}((400,),2:3)
E = BasisFunctions.ExtensionIndexMatrix{Float64}((400,),2:3)
@test R' == E
@test E' == R
@test R*D == Matrix(R)*Matrix(D)
@test D*E == Matrix(D)*Matrix(E)
@test R*D isa AbstractSparseArray
@test D*E isa AbstractSparseArray

V = BasisFunctions.VerticalBandedMatrix(400,200, collect(1:3))
H = BasisFunctions.HorizontalBandedMatrix(200,400, collect(1:3))
@test V' == H
@test H' == V
@test R*V == Matrix(R)*Matrix(V)
@test H*E == Matrix(H)*Matrix(E)
@test R*V isa AbstractSparseArray
@test H*E isa AbstractSparseArray

B = rand(400,2)
@test R*B == Matrix(R)*Matrix(B)
B = rand(2,400)
@test B*E == Matrix(B)*Matrix(E)

include("test_normalization.jl")
include("test_space.jl")
include("test_discrete_measure.jl")
include("test_measure.jl")


text:
@testset "AbstractAxis" begin
    @testset "keys" begin
        axis = Axis(2:3 => 1:2)

        @test keytype(typeof(Axis(1.0:10.0))) <: Float64
        @test keytype(typeof(parent(axis))) <: Int
        @test haskey(axis, 3)
        @test !haskey(axis, 4)
        @test keys.(axes(axis)) == (2:3,)

        A = AxisArray(ones(3,2), [:one, :two, :three], nothing)
        @testset "reverse" begin
            x = [1, 2, 3]
            y = AxisArray(x)
            z = AxisArray(x, Axis([:one, :two, :three]))

            revx = reverse(x)
            revy = @inferred(reverse(y))
            revz = @inferred(reverse(z))

            @testset "reverse vectors values properly" begin
                @test revx == revz == revy
            end

            @testset "reverse vectors keys" begin
                @test keys(axes(revy, 1)) == [3, 2, 1]
                @test keys(axes(revz, 1)) == [:three, :two, :one]
            end

            @testset "reverse arrays" begin
                b = [1 2; 3 4]
                x = AxisArray(b, [:one, :two], ["a", "b"])

                xrev1 = reverse(x, dims=1)
                xrev2 = reverse(x, dims=2)
                @test keys.(axes(xrev1)) == ([:two, :one], ["a", "b"])
                @test keys.(axes(xrev2)) == ([:one, :two], ["b", "a"])
            end
        end
    end

    @testset "step(r)" begin
        for (r,b) in ((SimpleAxis(OneToMRange(10)), OneToMRange(10)),)
            @test @inferred(step(r)) === step(b)
            @test @inferred(Base.step_hp(r)) === Base.step_hp(b)
            if b isa StepRangeLen
                @test @inferred(stephi(r)) == stephi(b)
                @test @inferred(steplo(r)) == steplo(b)
            end
        end

        @test @inferred(step(SimpleAxis(2))) == 1
        @test @inferred(firstindex(Axis(1:10))) == firstindex(1:10)
    end

    @testset "to_index" begin
        a = Axis(2:10)
        @test @inferred(to_index(a, 1)) == 1
        @test @inferred(to_index(a, 1:2)) == 1:2
        x = Axis([:one, :two])
        @test @inferred(to_index(x, :one)) == 1
        @test @inferred(to_index(x, [:one, :two])) == [1, 2]

        x = Axis(0.1:0.1:0.5)
        @test @inferred(to_index(x, 0.3)) == 3

        x = Axis(["a", "b"])
        @test @inferred(to_index(x, "a")) == 1
        @test @inferred(to_index(x, ==("b"))) == 2
        @test @inferred(to_index(x, "b")) == 2

        @test @inferred(to_index(x, 2)) == 2
        @test @inferred(to_index(x, ["a", "b"])) == [1, 2]
        @test @inferred(to_index(x, in(["a", "b"]))) == [1, 2]
        @test @inferred(to_index(x, 1:2)) == [1, 2]
        @test @inferred(to_index(x, [false, true])) == [2]
        @test @inferred(to_index(x, [true, true])) == [1, 2]
        @test @inferred(to_index(x, true)) == 1
        @test @inferred(to_index(x, :)) == Base.Slice(parent(x))

        @test_throws BoundsError to_index(x, ==("c"))
        @test_throws BoundsError to_index(x, "c")
        @test_throws BoundsError to_index(x, 3)
        @test_throws BoundsError to_index(x, ["a", "b", "c"])
        @test_throws BoundsError to_index(x, 1:3)
        @test_throws BoundsError to_index(x, [true, true, true])
        @test_throws BoundsError to_index(x, false)
    end

    @testset "to_indices" begin
        A = AxisArray(ones(2,2),  (Axis(1:2), Axis(1.0:2.0)));
        V = AxisArray(ones(2), ["a", "b"]);

        @testset "linear indexing" begin
            @test @inferred(AxisIndices.to_indices(A, (1,))) == (1,)
            @test @inferred(AxisIndices.to_indices(A, (1:2,))) == (1:2,)

            @testset "Linear indexing doesn't ruin vector indexing" begin
                @test @inferred(AxisIndices.to_indices(V, (1:2,))) == (1:2,)
                @test @inferred(AxisIndices.to_indices(V, (1,))) == (1,)
                @test @inferred(AxisIndices.to_indices(V, ("a",))) == (1,)
            end
        end

        @test @inferred(AxisIndices.to_indices(A, (1, 1))) == (1, 1)
        @test @inferred(AxisIndices.to_indices(A, (1, 1:2))) == (1, 1:2)
        @test @inferred(AxisIndices.to_indices(A, (1:2, 1))) == (1:2, 1)
        @test @inferred(AxisIndices.to_indices(A, (1, :))) == (1, Base.Slice(Axis(1.0:2.0)))
        @test @inferred(AxisIndices.to_indices(A, (:, 1))) == (Base.Slice(Axis(1:2)), 1)
        @test @inferred(AxisIndices.to_indices(A, ([true, true], :))) == (Base.LogicalIndex(Bool[1, 1]), Base.Slice(Axis(1.0:2.0)))
        @test @inferred(AxisIndices.to_indices(A, (CartesianIndices((1,)), 1))) == (Axis(1:1 => 1:1), 1)
        @test @inferred(AxisIndices.to_indices(A, (1, 1.0))) == (1,1)

        A = AxisArray(reshape(1:27, 3, 3, 3));
        @test @inferred AxisIndices.to_indices(A, ([CartesianIndex(1,1,1), CartesianIndex(1,2,1)],)) == (CartesianIndex{3}[CartesianIndex(1, 1, 1), CartesianIndex(1, 2, 1)],)
        @test @inferred(A[[CartesianIndex(1,1,1), CartesianIndex(1,2,1)]]) == [1, 4]
    end

    @testset "to_axes" begin
        A = ones(2,2)
        args = (:,:)
        new_indices = axes(A)

        @test @inferred(to_axes(AxisArray(A), new_indices)) ==
              (SimpleAxis(Base.OneTo(2)), SimpleAxis(Base.OneTo(2)))

        A = ones(2)
        args = (:,)
        new_indices = axes(A)
        @test @inferred(to_axes(AxisArray(A), new_indices)) == (SimpleAxis(Base.OneTo(2)),)

        A = ones(2,2)
        args = (:,1)
        new_indices = (Base.OneTo(2),)
        @test @inferred(to_axes(AxisArray(A), (Base.OneTo(2), 1))) == (SimpleAxis(Base.OneTo(2)),)
    end

    @testset "checkindex" begin
        @test @inferred(checkindex(Bool, Axis([:a, :b]), :a))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), :c))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), 1))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), 3))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), true))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), [:a, :b]))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), [:a, :c]))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), [1, 2]))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), [1, 3]))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), [true, true]))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), [true, true, true]))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), 1..2))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), in([:a, :b])))
        #@test !@inferred(checkindex(Bool, Axis([:a, :b]), in([:a, :c])))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), ==(:a)))
        @test !@inferred(checkindex(Bool, Axis([:a, :b]), ==(:c)))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), <(:b)))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), <(2)))
        @test @inferred(checkindex(Bool, Axis([:a, :b]), :))

    @testset "checkbounds" begin
        x = Axis(1:10)
        @test Base.checkindex(Bool, x, Base.Slice(1:10))
        @test Base.checkindex(Bool, x, [1,2,3])
        @test !Base.checkindex(Bool, x, [0, 1,2,3])

        @test checkbounds(Bool, Axis(1:2), CartesianIndex(1))
        @test !checkbounds(Bool, Axis(1:2), CartesianIndex(3))

        x2 = Axis(1:2)

        @test checkbounds(Bool, x2, 2)
        @test !checkbounds(Bool, x2, 3)

        @test Base.checkindex(Bool, x2, x2 .> 3)

        @test Base.checkindex(Bool, x2, 1:2)
        @test !Base.checkindex(Bool, x2, 1:3)

        @test Base.checkindex(Bool, x2, [true, true])
        @test !Base.checkindex(Bool, x2, [true, true, true])

        @test Base.checkindex(Bool, x2, first(to_indices(1:2, ([true, true],))))
        @test !Base.checkindex(Bool, x2, first(to_indices(1:2, ([true, true, true],))))

        @test Base.checkindex(Bool, x2, 1:1:2)
        @test !Base.checkindex(Bool, x2, 1:1:3)

        # trigger errors when functions return bad indices
        @test_throws BoundsError AxisIndices.to_index(Axis(1:10), ==(11))

    end
        #= TODO I don't think this should ever happen
            @testset "CartesianElement" begin
                @test @inferred(checkindex(Bool, Axis([:a, :b]), CartesianIndex(1)))
                #@test !@inferred(checkindex(Bool, Axis([:a, :b]), CartesianIndex(3)))
            end
        =#
    end

    @testset "values/indices" begin
        @test valtype(typeof(Axis(1.0:10.0))) <: Int
        a1 = Axis(2:3 => 1:2)

        @test allunique(a1)
        @test in(2, a1)
        @test !in(3, a1)
        @test eachindex(a1) == 1:2

        @testset "Floats as keys #13" begin
            A = AxisArray(collect(1:5), 0.1:0.1:0.5)
            @test @inferred(A[0.3]) == 3
        end
    end

    @testset "isequal" begin
        x = 1:10
        y = 1:1:10
        z = StaticRanges.GapRange(1:5,6:10)
        axis = SimpleAxis(x)
        @test ==(x, axis)
        @test ==(axis, x)
        @test ==(y, axis)
        @test ==(axis, y)
        @test ==(z, axis)
        @test ==(axis, z)
        @test ==(axis, axis)
    end
end

@testset "pop" begin
    for x in (Axis(UnitMRange(1,10),UnitMRange(1,10)),
              SimpleAxis(UnitMRange(1,10)))
        y = collect(x)
        #@test pop(x) == pop(y)
        @test pop!(x) == pop!(y)
        @test x == y
    end

    r = UnitMRange(1, 1)
    y = collect(r)
    @test pop!(r) == pop!(y)
    @test isempty(r) == true
end

#= FIXME popfirst
@testset "popfirst" begin
    for x in (Axis(UnitMRange(1,10),UnitMRange(1,10)),
              SimpleAxis(UnitMRange(1,10)))
        y = collect(x)
        #@test popfirst(x) == popfirst(y)
        @test popfirst!(x) == popfirst!(y)
        @test x == y
    end
    r = UnitMRange(1, 1)
    y = collect(r)
    @test popfirst!(r) == popfirst!(y)
    @test isempty(r) == true
end
=#

@testset "last" begin
   @testset "can_set_last" begin
       @test @inferred(can_set_last(typeof(Axis(UnitMRange(1:2))))) == true
       @test @inferred(can_set_last(typeof(Axis(UnitSRange(1:2))))) == false
    end

    for (r1,b,v,r2) in ((SimpleAxis(UnitMRange(1,3)), true, 2, SimpleAxis(UnitMRange(1,2))),
                        (Axis(UnitMRange(1,3),UnitMRange(1,3)), true, 2, Axis(UnitMRange(1,2),UnitMRange(1,2))))
        @testset "set_last-$(r1)" begin
            x = @inferred(can_set_last(typeof(r1)))
            @test x == b
            if x
                @test @inferred(set_last!(r1, v)) == r2
            end
            if x
                @test @inferred(set_last(r1, v)) == r2
            end
        end
    end

    @test last(Axis(2:3)) == 2
    @test last(Axis(2:3, 2:3)) == 3
end

# FIXME? should we be able to change the first index of an axis?
#=
@testset "first" begin
    @testset "can_set_first" begin
        @test @inferred(!StaticRanges.can_set_first(Axis{Int,Int,UnitRange{Int},Base.OneTo{Int}}))
        @test @inferred(StaticRanges.can_set_first(Axis{Int,Int,UnitMRange{Int},UnitMRange{Int}}))
    end

    for (r1,b,v,r2) in ((SimpleAxis(UnitMRange(1,3)), true, 2, SimpleAxis(UnitMRange(2,3))),
                        (Axis(UnitMRange(1,3),UnitMRange(1,3)), true, 2, Axis(UnitMRange(2,3),UnitMRange(2,3))))
        @testset "set_first-$(r1)" begin
            x = @inferred(can_set_first(r1))
            @test x == b
            if x
                set_first!(r1, v)
                @test r1 == r2
            end
            @test set_first(r1, v) == r2
        end
    end

    @test first(Axis(2:3)) == 1
    @test first(Axis(2:3, 2:3)) == 2
end
=#

@testset "length - tests" begin
    @testset "length(r)" begin
        for (r,b) in ((SimpleAxis(UnitMRange(1,3)), 1:3),
                      (Axis(UnitMRange(1,3),UnitMRange(1,3)), 1:3)
                     )
            @test @inferred(length(r)) == length(b)
            @test @inferred(length(r)) == length(b)
            if b isa StepRangeLen
                @test @inferred(stephi(r)) == stephi(b)
                @test @inferred(steplo(r)) == steplo(b)
            end
        end
    end

    @testset "can_set_length" begin
        @test @inferred(!StaticRanges.can_set_length(Axis{Int,Int,UnitRange{Int},Base.OneTo{Int}}))
        @test @inferred(StaticRanges.can_set_length(Axis{Int,Int,UnitMRange{Int},OneToMRange{Int}}))
    end

    @testset "set_length!" begin
        @test @inferred(set_length!(SimpleAxis(OneToMRange(10)), UInt32(11))) == SimpleAxis(OneToMRange(11))
        @test @inferred(set_length!(Axis(OneToMRange(10), OneToMRange(10)), UInt32(11))) == Axis(OneToMRange(11), OneToMRange(11))
    end

    @testset "set_length" begin
        @test @inferred(set_length(SimpleAxis(OneToMRange(10)), UInt32(11))) == SimpleAxis(OneToMRange(11))
        @test @inferred(set_length(Axis(OneToMRange(10), OneToMRange(10)), UInt32(11))) == Axis(OneToMRange(11), OneToMRange(11))
    end

    @testset "empty length" begin
        @test length(empty!(Axis(UnitMRange(1, 10)))) == 0
        @test length(empty!(SimpleAxis(UnitMRange(1, 10)))) == 0
    end
end


text:
"""
Kinetic.jl: A Portable Framework for Computational Fluid Dynamics and Scientific Machine Learning

Copyright (c) 2021 Tianbai Xiao (tianbaixiao@gmail.com)
"""
module Kinetic

export 転

"""
Portable Framework for Computational Fluid Dynamics and Scientific Machine Learning

轻量化的计算流体力学建模和计算框架

"転" means "rolling" in Japannese
"""
const 転 = Kinetic

using Reexport
@reexport using KitBase
@reexport using KitML

end


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 3.1 Basic Operations

p=1000;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("1"); n=BigInt(201); iN=convert(Int64,n); dx=BigFloat("2")*L/(n-1); x=-L:dx:L;
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n)
inv_V_comp = V\I; det_V_comp=det(V);
inv_V_analyt,det_V_analyt=TaylorBigF.calc_Vandermonde_inv_formula(x,n);
diff_V=maximum(abs.(inv_V_comp.-inv_V_analyt));@sprintf("diff_V=%1.3e",diff_V)
diff_det=det_V_analyt-det_V_comp;@sprintf("det_V_comp=%1.3e det_V_analyt=%1.3e diff_det=%1.3e",det_V_comp,det_V_analyt,diff_det)

f=sin.(x); a=V\f;
# plot(a)
r_root=Array{Float64}(undef,0); for i=1:Int64(n) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
# ticks = round.(range(minimum(r_rootF),stop=maximum(r_rootF),length = 4),digits=2);
# plot(1:Int64(n),r_rootF,color=:black,legend=false,xlabel="n",ylabel="1/r",yticks=ticks)
# scatter!(1:Int64(n),r_rootF,color=:black)
# savefig(string("r_sin",p,".pdf"))
# r_ratio=Array{Float64}(undef,0); for i=1:Int64(n)-1 push!(r_ratio,abs(a[i+1]/a[i])) end; plot(r_ratio[2:2:end])

t = Taylor1(Float64, Int64(n)-1); ts=sin(t)
a_analyt=Array{Float64}(undef,0); for i=1:Int64(n) push!(a_analyt,getcoeff(ts,i-1)) end
diff_a=maximum(abs.(a_analyt-a));@sprintf("diff_a=%1.3e",diff_a)

f_comp=V*a
diff_f=maximum(abs.(f-f_comp));@sprintf("diff_f=%1.3e",diff_f)
@sprintf("Rn=%1.3e",1/factorial(n))


text:
import HoTRG_lemon.SpinModule:
	initialize!

type Classical2dFractalLattice{T} <: Classical2dLattice{T}
	spinmodel::ClassicalSpinModel{T}
	legextension::Int
	tensorT::Array{T,4} # ; x,xp,y,yp
	tensorTtilde::Array{T,4}
	tensorP::Array{Array{T,3},1}# leg tensor; x,xp,s
	tensorQ::Array{T,2} # corner tensor; x,y
	function Classical2dFractalLattice{T}(spinmodel::ClassicalSpinModel{T}, legextension::Int)
		this = new{T}()
 		this.spinmodel = spinmodel
		this.legextension = legextension
 		initialize!(this)
 		return this
 	end
end

Classical2dFractalLattice{T}(spinmodel::ClassicalSpinModel{T}, legextension::Int) =
Classical2dFractalLattice{T}(spinmodel, legextension)
Classical2dFractalLattice{T}(spinmodel::ClassicalSpinModel{T}) =
Classical2dFractalLattice{T}(spinmodel, 1)

#---
# initialize functions for constructor
function initialize!{T}(lattice::Classical2dFractalLattice{T})
	lattice.tensorT = initializeTensorT2D_fractal(lattice.spinmodel)
	lattice.tensorTtilde = initializeTensorT2D_fractal(lattice.spinmodel; tilde = true)
	factorW = getFactorW(lattice)
	tensorP = makeTensorP(factorW)
	lattice.tensorP = fill(tensorP, getSpaceDimension(lattice))
	setTensorQ!(lattice, makeTensorQ(factorW))
end

function initializeTensorT2D_fractal{T}(spinmodel::ClassicalSpinModel{T}; tilde::Bool = false)
	if tilde
		measureOperator = getMeasureOperator(spinmodel)
	else
		measureOperator = ones(T,getStates(spinmodel))
	end
	factorW = getFactorW(spinmodel)
	tensorT = makeTensorT(factorW, measureOperator)
	return tensorT
end

#---
# functions for Classical2dFractalLattice
function getCoarserate(lattice::Classical2dFractalLattice)
	L = getLegextension(lattice)
	return 4*(L+1)
end

#---
# functions for test
#####
function testClassical2dFractalLattice()
	spinmodel = testIsingModel()
	latticecode = "classical_2d_fractal_1"
	lattice = buildLattice(latticecode, spinmodel)
	return lattice
end


text:
using DataStructures

const data = split.(readlines(), Ref("-"))
const graph = begin
    g = DefaultDict{String, Vector{String}}([])
    for (u, v) in data
        push!(g[u], v)
        push!(g[v], u)
    end
    g
end

issmall(node) = all(islowercase.(collect(node)))

function dfs(node, small_seen; double_visit=nothing)
    node == "end" && return 1
    s = 0
    for nei in graph[node]
        if !issmall(nei) || nei ∉ small_seen
            small_seen_cp = copy(small_seen)
            issmall(nei) && push!(small_seen_cp, nei)
            s += dfs(nei, small_seen_cp; double_visit)
        elseif isnothing(double_visit) && nei != "start"
            small_seen_cp = copy(small_seen)
            s += dfs(nei, small_seen_cp; double_visit=nei)
        end
    end
    s
end

# Part One
dfs("start", Set(["start"]); double_visit="start") |> println

# Part Two
dfs("start", Set(["start"])) |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
struct ScalarField{T<:Union{Float64,Float32},N,N2,L} <: AbstractPaddedArray{T,N}
    field::PaddedArray{T,N,N2,L}
    fplan::FFTW.rFFTWPlan{T,-1,true,N}
    bplan::FFTW.rFFTWPlan{Complex{T},1,true,N}
    k::VecArray{T,3,HomogeneousArray{T,3,RKvec{T},1},HomogeneousArray{T,3,Kvec{T},2},HomogeneousArray{T,3,Kvec{T},3}}
    realspace::Base.RefValue{Bool}

    function ScalarField(field::PaddedArray{T,N,N2,L},l) where {T,N,N2,L}
        fplan = plan_rfft!(field,1:N,flags=FFTW.MEASURE)
        bplan = plan_brfft!(field,1:N,flags=FFTW.MEASURE)
        s = size(field)
        k = VecArray(HomogeneousArray{1}(RKvec{T}(size(field.r,1),l[1]),s),
                     HomogeneousArray{2}(Kvec{T}(size(field.r,2),l[2]),s),
                     HomogeneousArray{3}(Kvec{T}(size(field.r,3),l[3]),s))
        return new{T,N,N2,L}(field,fplan,bplan,k,Ref(true))
    end

end

@inline function Base.getproperty(a::S,s::Symbol) where {S<:ScalarField}
    if s === :rr
        return getfield(getfield(a,:field),:data)
    elseif s === :r
        return getfield(getfield(a,:field),:r)
    elseif s === :c
        return getfield(a,:field)
    elseif s === :kx
        return getfield(getfield(getfield(a,:k),:x),:data)
    elseif s === :ky
        return getfield(getfield(getfield(a,:k),:y),:data)
    elseif s === :kz
        return getfield(getfield(getfield(a,:k),:z),:data)
    else
        return getfield(a,s)
    end
end

ScalarField{T}(dims::NTuple{3,Int},l::NTuple{3,Real}) where {T} = ScalarField(PaddedArray{T}(dims...),l)
ScalarField(dims::NTuple{3,Int},l::NTuple{3,Real}) = ScalarField{Float64}(dims,l)

function ScalarField(file,dim::NTuple{N,Int},l::Vararg{real,3}) where N
    r = ScalarField(PaddedArray(dim...),l)
    read!(file,r.field.data)
    return r
end

Base.similar(a::ScalarField) = ScalarField(similar(a.field),(a.kx.l,a.ky.l,a.kz.l))

@inline InplaceRealFFT.data(a::ScalarField) = InplaceRealFFT.data(a.field)
@inline Base.real(a::ScalarField) = real(a.field)
@inline InplaceRealFFT.complex_view(a::ScalarField) = InplaceRealFFT.complex_view(a.field)

function Base.:*(p::FFTW.rFFTWPlan{T,FFTW.FORWARD,true,N},f::ScalarField{T,N,N2,L}) where {T<:Float3264,N,N2,L}
    isrealspace(f) || error("Field already is in Fourier space")
    mul!(InplaceRealFFT.complex_view(f),p,real(f))
    f.realspace[] = false
    f
end

InplaceRealFFT.rfft!(a::ScalarField) =
    a.fplan * a

function Base.:*(p::FFTW.rFFTWPlan{Complex{T},FFTW.BACKWARD,true,N},f::ScalarField{T,N,N2,L}) where {T<:Float3264,N,N2,L}
    isrealspace(f) && error("Field already is in real space")
    mul!(real(f),p,InplaceRealFFT.complex_view(f))
    f.realspace[] = true
    real(f)
end

InplaceRealFFT.brfft!(a::ScalarField) =
    a.bplan * a

function InplaceRealFFT.irfft!(a::ScalarField)
    brfft!(a)
    rmul!(InplaceRealFFT.data(a),inv(prod(size(real(a)))))
    real(a)
end

isrealspace(a::ScalarField) =
    a.realspace[]

Base.show(io::IO,m::MIME"text/plain",a::ScalarField) =
    a.realspace[] ? show(io,m,real(a)) : show(io,m,a.field)

setreal!(a::ScalarField) = (a.realspace[] = true)
setfourier!(a::ScalarField) = (a.realspace[] = false)

function fourier!(a::ScalarField)
    rfft!(a)
    rmul!(InplaceRealFFT.data(a),inv(prod(size(real(a)))))
    return a
end

real!(a::ScalarField) = brfft!(a)


text:
using Gencode

function findexec(filename, version)

status = Sys.which(filename);
if status != nothing
    return filename;
end
status = Sys.which("/usr/bin/" * filename);
if status != nothing
    filename = "/usr/bin/" * filename;
    return filename;
end
status = Sys.which("/usr/local/bin/" * filename);
if status != nothing
    filename = "/usr/local/bin/" * filename;
    return filename;
end
status = Sys.which("/opt/local/bin/" * filename);
if status != nothing
    filename = "/opt/local/bin/" * filename;
    return filename;
end
if status == nothing
    print("Exasim could not find " * filename * " in /usr/bin, /usr/local/bin, /opt/local/bin\n");
    if Sys.isapple()
        print("Exasim tries to find " * filename * " in /Applications. It may take a while.\n");
        a = read(Gencode.string2cmd("find /Applications -name " * filename * " -type f"), String);
        if length(a)>0
            ii = findlast("/MacOS/" * filename, a);
            newfilename = a[1:ii[end]];
            print("Exasim found " * filename * " at " * newfilename * "\n");
            print("Please open initializepde.jl in the folder Exasim/" * version * "/Julia/Preprocessing\n");
            qstr1 = """ "$filename" """;
            qstr2 = """ "$newfilename" """;
            print("Replace pde." * filename * " =" * qstr1 * " with pde." * filename * " =" * qstr2 * "\n");
            print("Doing so will prevent Exasim from searching " * filename * " again to save time.\n");
            print("Read the above instructions carefully and press any key to continue...\n");
            sleep(20);
            filename = newfilename;
            return filename;
        end
    # elseif Sys.isunix()
    # elseif Sys.windows()
    end
    mystr = "Exasim could not find " * filename * " on your system.\n";
    mystr = mystr * "Please see the documentation to install " * filename * ".\n"
    mystr = mystr * "After installing " * filename * ", open initializepde.jl in the folder Exasim/" * version * "/Julia/Preprocessing\n";
    qstr1 = """ "$filename" """;
    newfilename = "/path/to/executable/" * filename;
    qstr2 = """ "$newfilename" """;
    mystr = mystr * "and replace pde." * filename * " =" * qstr1 * "with pde." * filename * " =" * qstr2 * "\n";
    print(mystr);
    error("Exit Exasim!")
end

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# models.jl
using Flux
using Flux: glorot_uniform
struct GCNLayer
    W
    b
    adj
    act
end

struct GCNModel
    layer1::GCNLayer
    layer2::GCNLayer
    drop
end

function GCNLayer(in::Integer, out::Integer, adj, act=relu; initW=glorot_uniform,
                        initb=zeros)
    return GCNLayer(param(initW(in, out)), param(initb(out)), adj, act)
end

function (a::GCNLayer)(x::AbstractArray)
    support = x * a.W
    return a.act.(a.adj * support .+ a.b')
end

Flux.@treelike GCNLayer

function GCNModel(nfeat, nhid, nclass, drop, adj)
    return GCNModel(GCNLayer(nfeat, nhid, adj), GCNLayer(nhid, nclass, adj, identity), drop)
end

function (a::GCNModel)(x::AbstractArray)
    z = Chain(a.layer1, Dropout(a.drop), a.layer2)
    return softmax(z(x)')'
end

Flux.@treelike GCNModel

max_stars_repo_name: mcbaron/fluxGCN

max_issues_repo_name: mcbaron/fluxGCN

max_forks_repo_name: mcbaron/fluxGCN


text:
#!/usr/bin/env julia

# not the correct approach

function solve(case, N, C)
    N - 1 <= C <= sum(1:N) - Int(iseven(N)) || (println("Case #$case: IMPOSSIBLE"); return)
    S = collect(1:N)

    i = lastindex(S) - 1
    while C > 0
        c = lastindex(S) - i + 1
        if c <= C - i
            #@info C i c
            @views S[i:end] .= Iterators.reverse(S[i:end])
            C -= c
        else
            C -= 1
        end
        i -= 1
    end
    println("Case #$case: $(join(S, ' '))")
end

function main()
    T = parse(Int, readline())
    for case in 1:T
        N, C = map(x -> parse(Int, x), split(readline()))
        solve(case, N, C)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Stheno: CondCache, ConditionalMean, ConditionalKernel, ConditionalCrossKernel,
    ConstantMean, ZeroMean, ZeroKernel, ConstantKernel, pairwise, chol

@testset "conditional" begin

    # Tests for conditioning with independent processes.
    let
        rng, N, N′, D = MersenneTwister(123456), 10, 13, 2
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        μf, μg, μh = ConstantMean(randn(rng)), ZeroMean{Float64}(), ConstantMean(randn(rng))
        kff, kgg, khh = EQ(), EQ() + Noise(1e-6), Noise(1e-6)
        kfg, kfh, kgh = EQ(), ZeroKernel{Float64}(), ZeroKernel{Float64}()

        # Run all tests for both scalar and vector input domains.
        for (X0, X1, X2) in [(X0, X1, X2), (x0, x1, x2)]

            # Construct conditioned objects.
            y = map(μf, X0) + chol(LazyPDMat(pairwise(kff, X0)))' * randn(rng, N)
            cache = CondCache(kff, μf, X0, y)
            kffX = LhsFiniteCrossKernel(kff, X0)
            kfgX = LhsFiniteCrossKernel(kfg, X0)
            kfhX = LhsFiniteCrossKernel(kfh, X0)
            μ′f = ConditionalMean(cache, μf, kffX)
            μ′g = ConditionalMean(cache, μg, kfgX)
            μ′h = ConditionalMean(cache, μh, kfhX)
            k′ff = ConditionalKernel(cache, kffX, kff)
            k′gg = ConditionalKernel(cache, kfgX, kgg)
            k′hh = ConditionalKernel(cache, kfhX, khh)
            k′fg = ConditionalCrossKernel(cache, kffX, kfgX, kfg)
            k′fh = ConditionalCrossKernel(cache, kffX, kfhX, kfh)
            k′gh = ConditionalCrossKernel(cache, kfgX, kfhX, kgh)

            # Run standard consistency tests on all of the above.
            for (n, μ′) in enumerate([μ′f, μ′g, μ′h])
                mean_function_tests(μ′, X0)
                mean_function_tests(μ′, X2)
                mean_function_tests(μ′, BlockData([X0, X2]))
            end
            for (n, k′) in enumerate([k′ff, k′gg, k′hh])
                kernel_tests(k′, X0, X1, X2, 1e-6)
                kernel_tests(k′, BlockData([X0]), BlockData([X1]), BlockData([X2]), 1e-6)
                kernel_tests(k′, BlockData([X0]), BlockData([X1]), BlockData([X2, X1]), 1e-6)
            end
            for k′ in [k′fg, k′fh, k′gh]
                cross_kernel_tests(k′, X0, X1, X2)
                cross_kernel_tests(k′, BlockData([X0]), BlockData([X1]), BlockData([X2]))
                cross_kernel_tests(k′, BlockData([X0]), BlockData([X1]), BlockData([X2, X1]))
            end

            # Test that observing the mean function shrinks the posterior covariance
            # appropriately, but leaves the posterior mean at the prior mean (modulo noise).
            cache = CondCache(kff, μf, X0, map(μf, X0))
            @test cache.α ≈ zeros(Float64, N)

            # Posterior covariance at the data should be _exactly_ zero.
            @test maximum(abs.(pairwise(k′ff, X0))) < 1e-6

            # Posterior for indep. process should be _exactly_ the same as the prior.
            @test map(μ′h, X0) == map(μh, X0)
            @test map(μ′h, X1) == map(μh, X1)
            @test map(μ′h, X2) == map(μh, X2)
            @test map(k′hh, X0) == map(khh, X0)
            @test map(k′hh, X1) == map(khh, X1)
            @test map(k′hh, X2) == map(khh, X2)
            @test map(k′hh, X0, X1) == map(khh, X0, X1)
            @test pairwise(k′hh, X0) == pairwise(khh, X0)
            @test pairwise(k′hh, X0, X2) == pairwise(khh, X0, X2)
            @test pairwise(k′fh, X0, X2) == pairwise(kfh, X0, X2)
            @test pairwise(k′fh, X1, X2) == pairwise(kfh, X0, X2)

            # Ensure that printing doesn't error.
            show(IOBuffer(), μ′f)
            show(IOBuffer(), k′ff)
            show(IOBuffer(), k′fg)

            # Ensure lengths are propagated correctly.
            @test length(μ′f) == Inf
            @test size(k′ff, 1) == Inf && size(k′ff, 2) == Inf
            @test size(k′fg, 1) == Inf && size(k′fg, 2) == Inf
        end
    end

    #     # Should be able to get similar results using a ConditionalCrossKernel.
    #     @test xcov(ConditionalCrossKernel(cache, kff, kff, kff), X, X) ≈
    #         Matrix(cov(ConditionalKernel(cache, kff, kff), X))
end


text:
using ZonalFlow
using Test

domain  = Domain(extent=(2π,2π),res=(5,5));
coeffs  = Coefficients(Ω=2π,θ=0.0,μ=0.0,ν=0.0,ν₄=1.0,linear=true);
forcing = Stochastic(kf=3,dk=1,ε=0.0);
prob    = BetaPlane(domain,coeffs,forcing);

tspan   = (0.0,1000.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            saveat=500
           );

eqs     = [NL(),CE2()];
eqs     = append!(eqs,[GQL(Λ=l) for l=0:prob.d.nx-1])
eqs     = append!(eqs,[GCE2(Λ=l) for l=0:prob.d.nx-1])
sols    = integrate(prob,eqs,tspan;tsargs...);

@testset "Linear Equations" begin
        for sol in sols
                E,Z = energy(length(prob.d)...,size(prob.d)...,sol.t,sol.u);
                @test E[1] == E[end]
                @test Z[1] == Z[end]
        end
end

coeffs  = Coefficients(Ω=2π,θ=0.0,μ=0.0,ν=0.0,ν₄=1.0,linear=false);
prob    = BetaPlane(domain,coeffs,forcing);
sols    = integrate(prob,eqs,tspan;tsargs...);

@testset "Nonlinear Equations" begin
        for sol in sols
                E,Z = energy(length(prob.d)...,size(prob.d)...,sol.t,sol.u);
                @test E[1] == E[end]
                @test Z[1] == Z[end]
        end
end


text:
# from FastTransforms
function backend_fft_pow2!{T}(x::Vector{T})
    n,big2,bigpi=length(x),2one(T),T(pi)
    nn,j=n÷2,1
    @inbounds for i=1:2:n-1
        if j>i
            x[j], x[i] = x[i], x[j]
            x[j+1], x[i+1] = x[i+1], x[j+1]
        end
        m = nn
        while m ≥ 2 && j > m
            j -= m
            m = m÷2
        end
        j += m
    end
    logn = 2
    while logn < n
        piθ=-bigpi*big2/logn
        wtemp = sin(piθ/2)
        wpr, wpi = -2wtemp^2, sin(piθ)
        wr, wi = one(T), zero(T)
        for m=1:2:logn-1
            @inbounds for i=m:2logn:n
                j=i+logn
                mixr, mixi = wr*x[j]-wi*x[j+1], wr*x[j+1]+wi*x[j]
                x[j], x[j+1] = x[i]-mixr, x[i+1]-mixi
                x[i], x[i+1] = x[i]+mixr, x[i+1]+mixi
            end
            wr = (wtemp=wr)*wpr-wi*wpi+wr
            wi = wi*wpr+wtemp*wpi+wi
        end
        logn = logn << 1
    end
    return x
end

function fft_pow2{T}(x::Vector{T})
    y = FastTransforms.interlace(real(x),imag(x))
    backend_fft_pow2!(y)
@compat    return complex.(y[1:2:end],y[2:2:end])
end

## TAYLOR

taylorinterpolationerror(r,Ninterp::Int,M,R) = M*(r/R)^Ninterp/(1-r/R)

function taylortransform{T}(vals::Vector{Complex{T}},r::T,M::T,R::T)
  Ninterp = length(vals)
  cfs = fft_pow2(vals)/Ninterp
  err= erroradjustment(Complex{T},taylorinterpolationerror(r,Ninterp,M,R))
  powr = one(T)
  @inbounds for i = 1:Ninterp
    cfs[i] *= powr
    cfs[i] += err
    err /= R
    powr /= r
  end
  cfs
end

## CHEBYSHEV

function chebyshevtransform{T}(a::AbstractArray{Complex{T}})
  N = length(a) #big(length(a))
    ispow2(N) || error("N not a power of 2")
    c = fft_pow2([a; flipdim(a,1)])
    d = c[1:N]
 @compat    d .*= exp.((-im*convert(T,pi)).*(0:N-1)./(2*N))
     d[1] = d[1] / 2 #sqrt(big(2))
     scale!(inv(convert(T,N)), d)
end
@compat chebyshevtransform{T<:Real}(a::AbstractArray{T}) = real(chebyshevtransform(complex(a)))

function cos_vec!(A::Vector)
  for i in eachindex(A)
    A[i] = cos(A[i])
  end
  A
end
chebyshevpoints(N_interp,T) = cos_vec!(collect(1:2:2N_interp-1)*(T(pi)/2N_interp));

chebyshev_shifted_sq(x,k) = (-1)^k * cos(2k*asin(x))
chebyshev_shifted(x,k) = chebyshev_shifted_sq(sqrt(x),k)

function chebyshev_shifted_by_squaring{T}(x::T,k)
  tot=k;
  s0sq = 4x*(1-x)
  cn = 2x-1; sn_div = one(T)
  ctot = one(T); stot_div = zero(T)
  while tot>0
    if isodd(tot)
      ctot,stot_div = ctot*cn - stot_div*sn_div*s0sq, stot_div*cn - ctot*sn_div
    end
    sn_div *= 2cn; cn *= 2cn; cn-=1;
    tot ÷=2
  end
  ctot
end

function chebyshev_shifted_multiel{T}(x::T,K)
  Tk = Array{T}(K)
  mul = 4x-2
  K > 0 && (Tk[1] = one(typeof(x)))
  K > 1 && (Tk[2] = 2x-1)
  @inbounds for kk = 1:div(K-1,2)
    Tk[2kk+1] = 2Tk[kk+1]^2-1
    Tk[2kk+2] = mul*Tk[2kk+1]-Tk[2kk]
  end
  isodd(K) && (Tk[K] = 2Tk[div(K,2)+1]^2-1)
  Tk
end

function fejerweights(n,T::DataType=Float64)
  v = Array{Complex{T}}(n)
  v[1] = 2
  # twid = exp(-im*T(pi)/n) # conjugating here to get ifft
  # twidp = copy(twid)
  @inbounds for k = 1:div(n-1,2)
    twidp = exp(-im*k*T(pi)/n)
    v[k+1] = 2twidp/(1-4k^2)
    # twidp *= twid
  end
  mod(n,2)==0 && (v[div(n,2)+1]=0)
  @inbounds for k = 1:div(n-1,2)
    v[n-k+1] = conj(v[k+1])
  end
   real(fft_pow2(v))/n
end
 # fejerweights{T<:Real}(n,::Type{Complex{T}}) = fejerweights(n,T)

fejererror(N,M,ρ) = 2M/(ρ-1)*ρ^(-N)


text:
using RecipesBase
using ColorTypes
using LaTeXStrings
import PlotUtils: cgrad

#using Compat
#using Compat: range

const mygreen = RGBA{Float64}(151/255,180/255,118/255,1)
const mygreen2 = RGBA{Float64}(113/255,161/255,103/255,1)
const myblue = RGBA{Float64}(74/255,144/255,226/255,1)

@recipe function plot(b::Body)
    x = [b.x; b.x[1]]
    y = [b.y; b.y[1]]
    linecolor --> mygreen
    fillrange --> 0
    fillcolor --> mygreen
    aspect_ratio := 1
    legend := :none
    grid := false
    x := x
    y := y
    ()
end

@recipe function f(m::RigidBodyMotion;tmax=10)

    t = 0.0:0.01:tmax
    ux = map(ti -> real(m(ti)[2]),t)
    uy = map(ti -> imag(m(ti)[2]),t)
    adot = map(ti -> m(ti)[5],t)
    xlims --> (0,tmax)
  layout := 3
  grid --> :none
  linewidth --> 1
  legend --> :none
  framestyle --> :frame
  xguide --> L"t"
  ulim = min(minimum(ux),minimum(uy)),max(maximum(ux),maximum(uy))
  alim = extrema(adot)
  @series begin
      subplot := 1
      ylims --> ulim
      yguide --> L"u"
      t, ux
    end
  @series begin
      subplot := 2
      ylims --> ulim
      yguide --> L"v"
      t, uy
    end
   @series begin
      subplot := 3
      ylims --> alim
      yguide --> L"\dot{\alpha}"
      t, adot
    end
end

function RecipesBase.RecipesBase.apply_recipe(plotattributes::Dict{Symbol, Any}, bl::BodyList)
    series_list = RecipesBase.RecipeData[]
    for b in bl
        append!(series_list, RecipesBase.RecipesBase.apply_recipe(copy(plotattributes), b) )
    end
    series_list
end


text:
module GarchTest

using Distributions
using TimeData
using Base.Test
using NLopt

import Econometrics

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

############################
## constructors / display ##
############################

## GARCH(1,1), normal innovations
##-------------------------------

## test construction of GARCH(1,1) with normal
params = [0.02, 0.000004, 0.0906, 0.8959]

gMod = Econometrics.GARCH_1_1(params..., Normal())
gMod = Econometrics.GARCH_1_1(params, Normal())

@test isa(gMod, Econometrics.GARCH_1_1{Normal})

params = [0.02, -0.000004, 0.0906, 0.8959]
@test_throws Exception Econometrics.GARCH_1_1(params, Normal)

params = [0.02, -0.000004, 0.0906, 0.8959]
@test_throws Exception Econometrics.GARCH_1_1(params, Normal())

params = [0.02, 0.000004, -0.0906, 0.8959]
@test_throws Exception Econometrics.GARCH_1_1(params, Normal())

params = [0.02, 0.000004, 0.0906, -0.8959]
@test_throws Exception Econometrics.GARCH_1_1(params, Normal())

params = [0.02, 0.000004, 0.0906, 0.9959]
@test_throws Exception Econometrics.GARCH_1_1(params, Normal())

## GARCH(1,1), t innovations
##--------------------------

## test construction of GARCH(1,1) with Student's t
params = [0.02, 0.000004, 0.0906, 0.8959]

gMod = Econometrics.GARCH_1_1(0.02, 0.000004, 0.0906, 0.8959,
                              TDist(3.2))
@test isa(gMod, Econometrics.GARCH_1_1{TDist})

##############
## simulate ##
##############

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

nSim = 1000
params = [0.02, 0.000004, 0.0906, 0.8959]

## with normal innovations
gMod = Econometrics.GARCH_1_1(params, Normal())
simVals = Econometrics.simulate(gMod, nSim)

## with t innovations
gMod = Econometrics.GARCH_1_1(params, Econometrics.TDist(3.2))
simVals = Econometrics.simulate(gMod, nSim)

## with normal innovations and initial sigma
gMod = Econometrics.GARCH_1_1(params, Normal())
simVals = Econometrics.simulate(gMod, nSim, 0.01)

## with t innovations and initial sigma
gMod = Econometrics.GARCH_1_1(params, Econometrics.TDist(3.2))
simVals = Econometrics.simulate(gMod, nSim, 0.01)

###############
## getSigmas ##
###############

params = [0.02, 0.000004, 0.0906, 0.8959]
gMod = Econometrics.GARCH_1_1(params, Normal())

sigmas = Econometrics.getSigmas(gMod, simVals, 0.02)

##############
## estimate ##
##############

## GARCH(1,1), normal innovations
##-------------------------------

nSim = 1000
params = [0.02, 0.000004, 0.0906, 0.8959]
gMod = Econometrics.GARCH_1_1(params, Normal())
simVals = Econometrics.simulate(gMod, nSim)
gModHat = Econometrics.estimate(Econometrics.GARCH_1_1{Normal},
                           simVals)

## t innovations for normal data
gModHat = Econometrics.estimate(Econometrics.GARCH_1_1{TDist},
                           simVals)

## t innovations, t data
gMod = Econometrics.GARCH_1_1(params, TDist(4.3))
simVals = Econometrics.simulate(gMod, nSim)
gModHat = Econometrics.estimate(Econometrics.GARCH_1_1{TDist},
                           simVals)

gMod

#########
## fit ##
#########

gModFit = Econometrics.fit(Econometrics.GARCH_1_1{Normal},
                           simVals)

gModFit = Econometrics.fit(Econometrics.GARCH_1_1{TDist},
                           simVals)

end


text:
## tree
struct Tree{newicktype, Teltype, Pareltype, idseltype, nameseltype}
    newick::newicktype   #   (:id, :data, :parent, :children)
    T::Vector{Teltype}
    Par::Vector{Pareltype}
    ids::Vector{idseltype}
    names::Vector{nameseltype}
    lids::Vector{idseltype}
    lnames::Vector{nameseltype}
    lastone::Vector{Bool} # true if it is the last child when running forwards
    n::Int64
end

function Tree(S::String)
    nwtree = readnw(S)

    ids = map(x->x.id, prewalk(nwtree))
    names = map(x->x.data.name, prewalk(nwtree))
    n = length(ids)
    # leaves
    l = getleaves(nwtree)
    lids = map(x -> Int(x.id), l)
    lnames = map(x->x.data.name, l)
    #@assert issubset(tipslist[:,1], lnames)

    # extract parents and time spans (on segments)
    Par = zeros(Int, n)  # parent ids
    T = zeros(n)  # time spans
    for node in prewalk(nwtree)
        isroot(node) && continue
        p = node.parent.id
        Par[node.id] = p
        T[node.id] = T[p] + node.data.distance
    end
    @assert all(Par .< 1:n)

    lastone = zeros(Bool,n)
    parentdone = zeros(Bool,n)
    for i ∈ n:-1:2
        if !parentdone[Par[i]]
            lastone[i] = true
            parentdone[Par[i]] = true
        end
    end
    Tree{typeof(nwtree), eltype(T), eltype(Par), eltype(ids), eltype(names)}(nwtree,  T, Par, ids, names, lids, lnames,lastone, n)
end

## read tree (very simple example)
S = "(A:0.6,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
S20 = "(((t20:0.1377639426,t8:0.3216105914):0.461882598,(((t14:0.09650261351,t12:0.5344695284):0.1058950284,(t1:0.7008816295,t13:0.8587384557):0.7271672823):0.8026451592,(t6:0.5240857233,(t16:0.451590274,t10:0.4654036982):0.8569358122):0.140687139):0.7207935331):0.06727131596,((((t17:0.03115223325,t18:0.8930883685):0.6862687438,t9:0.4896655723):0.3056753092,t3:0.2549535078):0.07862622454,(((t4:0.03973490326,t5:0.4096179632):0.3020580804,t15:0.06550096255):0.2660011558,(t7:0.3895205685,((t19:0.165174433,t2:0.6469203965):0.2730371242,t11:0.2700119715):0.8393620371):0.02117527928):0.206376995):0.6012785432);"
S50 = "((((((t23:0.5766009684,t33:0.9917489749):0.3572315329,(t40:0.4535627225,t25:0.1296473157):0.9540394924):0.4740875189,(t6:0.946340692,t47:0.386682119):0.1535291784):0.07034701738,((((t49:0.6612032412,t27:0.1606025242):0.1527153626,t15:0.1854858669):0.5782065149,t50:0.3087045723):0.0759275395,((t17:0.258283512,t35:0.7651973176):0.2003453802,t2:0.4033703087):0.2418845498):0.6834833182):0.4268306631,(((t26:0.1744227773,(t16:0.6275201498,t10:0.09013780276):0.1857208994):0.8804157388,(t28:0.626011228,t46:0.8026940082):0.7066663445):0.6334459039,(t45:0.3075582262,t42:0.2759609232):0.05279792193):0.4003514573):0.5333229457,((((t30:0.05955170072,t41:0.2500256761):0.04172409023,((t19:0.3332289669,t31:0.6650944015):0.2450174955,t9:0.3372890623):0.2041792406):0.3337947708,t1:0.9396123393):0.0896910606,(((((t12:0.1379979164,t43:0.9007850592):0.703310129,(t11:0.3347885907,(t22:0.7190679635,t14:0.8988074451):0.1574828231):0.9644375315):0.137505424,((t37:0.7556577872,t8:0.3481238757):0.7712922962,(t7:0.7946650374,t48:0.8638419164):0.338934626):0.7507952554):0.5193078597,(t4:0.2408366969,(((t29:0.7862726359,t38:0.4319011718):0.8655222375,((t20:0.09156441619,t3:0.5975058537):0.3240908289,(t34:0.7324223334,t32:0.1348467385):0.4636785307):0.7708284601):0.7991992962,t39:0.2996040138):0.4087550903):0.1688272723):0.7225795928,(((t44:0.02794037666,t5:0.09033886855):0.5778779227,(t24:0.8744501146,((t21:0.2940180032,t18:0.03399693617):0.8231791619,t13:0.5192477896):0.5669933448):0.9547001948):0.8561708115,t36:0.9225172002):0.02867545444):0.1595393196):0.1254606678);"

# S20coal = "(((((t5:0.01408886814,t6:0.01408886814):0.03926824924,t3:0.05335711738):0.3769529348,t17:0.4303100522):0.6904821712,((t2:0.05400246327,(t19:0.04515631227,t7:0.04515631227):0.008846151006):0.02843109882,t15:0.08243356209):1.038358661):0.6759374458,(t12:0.5094458799,((t11:0.09067183363,(t4:0.08830667144,t1:0.08830667144):0.002365162193):0.2392304005,((t13:0.1499179245,((t10:0.01410801788,t20:0.01410801788):0.01360477778,t8:0.02771279566):0.1222051288):0.1061152139,(t16:0.1612987658,((t18:0.03933683292,t9:0.03933683292):0.008828557137,t14:0.04816539006):0.1131333758):0.09473437257):0.07386909574):0.1795436458):1.287283789);"

S20coal =  "(((t19:1.201139146,(t16:0.2402278292,t11:0.2402278292):0.9609113167):2.642506121,(t12:3.603417438,((t14:1.921822633,t17:1.921822633):0.9609113167,(t2:2.642506121,t7:2.642506121):0.2402278292):0.7206834875):0.2402278292):0.7206834875,((t20:3.122961779,((t6:0.4804556584,t4:0.4804556584):1.201139146,((t10:0.7206834875,t3:0.7206834875):0.2402278292,t9:0.9609113167):0.7206834875):1.441366975):1.201139146,(((t18:2.162050463,(t5:1.441366975,t8:1.441366975):0.7206834875):1.201139146,t13:3.363189608):0.7206834875,(t15:2.402278292,t1:2.402278292):1.681594804):0.2402278292):0.2402278292);"


text:
function rotate(	q0::Float64,
										q1::Float64,
										q2::Float64,
										q3::Float64,
										sigma_rotation::Float64)

	ux_rot::Float64 = randn()
	uy_rot::Float64 = randn()
	uz_rot::Float64 = randn()
	u_norm::Float64 = sqrt( ux_rot^2 + uy_rot^2 + uz_rot^2 )
	ux_rot /= u_norm
	uy_rot /= u_norm
	uz_rot /= u_norm

	theta_rot::Float64 = sigma_rotation * randn()

	q0_rot::Float64 = cos(0.5 * theta_rot)
	q1_rot::Float64 = ux_rot * sin(0.5 * theta_rot)
	q2_rot::Float64 = uy_rot * sin(0.5 * theta_rot)
	q3_rot::Float64 = uz_rot * sin(0.5 * theta_rot)

	(q0_star::Float64, q1_star::Float64, q2_star::Float64, q3_star::Float64) = quaternion_mult(q0_rot, q1_rot, q2_rot, q3_rot, q0, q1, q2, q3)

	# Normalize proposal quaternion to avoid numerical problems.
	q_norm::Float64 = sqrt( q0_star^2 + q1_star^2 + q2_star^2 + q3_star^2 )
	q0_star /= q_norm
	q1_star /= q_norm
	q2_star /= q_norm
	q3_star /= q_norm

	return (q0_star, q1_star, q2_star, q3_star)
end


text:
"""
    eom_logistic2(dx, x, p, n) -> function

Equations of motions for a system consisting of two coupled logistic maps where
X unidirectionally influences Y.

The parameter `c` controls how strong the dynamical forcing is. Parameters `r₁`
and `r₂` are set to the chaotic regime by default

# References
1. Diego, D., Agasøster Haaga, K., & Hannisdal, B. (2018, November 1).
Transfer entropy computation using the Perron-Frobenius operator.
Eprint ArXiv:1811.01677. Retrieved from
https://ui.adsabs.harvard.edu/#abs/2018arXiv181101677D
"""
function eom_logistic2(dx, x, p, n)
    c, r₁, r₂, σ = (p...,)
    ξ = rand() # random number from flat distribution on [0, 1]
    x, y = x[1], x[2]
    f_xy = (y +  (c*(x + σ*ξ)/2) ) / (1 + (c/2)*(1+σ))

    dx[1] = r₁ * x * (1 - x)
    dx[2] = r₂ * (f_xy) * (1 - f_xy)
    return
end

function logistic2(u₀, c, r₁, r₂, σ)
    p = [c, r₁, r₂, σ]
    DiscreteDynamicalSystem(eom_logistic2, u₀, p)
end

"""
    logistic2(;u₀ = rand(2), c = 0.1, σ = 0.05,
        r₁ = 3.78, r₂ = 3.66) -> DiscreteDynamicalSystem

Initialise a system consisting of two coupled logistic maps where X
unidirectionally influences Y. By default, the parameters `r₁` and `r₂` are set
to values yielding chaotic behaviour.

The parameter `c` controls how strong the dynamical forcing is. If `σ > 0`,
dynamical noise masking the influence of  `x` on `y` equivalent to
``\\sigma \\cdot \\xi`` is added at each iteration. Here,``\\xi`` is a draw from a
flat distribution on ``[0, 1]``. Thus, setting `σ = 0.05` is equivalent to
add dynamical noise corresponding to a maximum of ``5 \\%`` of the possible
range of values of the logistic map.

The equations of motion are

```math
\\begin{aligned}
dx &= r_1 x(1 - x) \\
dy &= r_2 f(x,y)(1 - f(x,y)),
\\end{aligned}
```
with
```math
\\begin{aligned}
f(x,y) = \\dfrac{y + \\frac{c(x \\xi )}{2}}{1 + \\frac{c}{2}(1+ \\sigma )}
\\end{aligned}
```

where

# References

1. Diego, D., Agasøster Haaga, K., & Hannisdal, B. (2018, November 1).
Transfer entropy computation using the Perron-Frobenius operator.
Eprint ArXiv:1811.01677. Retrieved from
https://ui.adsabs.harvard.edu/#abs/2018arXiv181101677D

"""
logistic2(;u₀ = rand(2), c = 0.1, r₁ = 3.78, r₂ = 3.66, σ = 0.05) =
    logistic2(u₀, c, r₁, r₂, σ)


text:
function evaluate(net::Network,test::Vector{Tuple{Vector{Float64},Vector{Float64}}},rule::Function)::Float64
    c = sum([1 for t in test if rule(feedforward(net,t[0])) == rule(t[1])])
    return "$(c) / $(length(test)) = $(100*(c/length(test)))%"
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using BinningAnalysis

mutable struct Observables
    energy::ErrorPropagator{Float64,32}
    magnetization::LogBinner{Float64,32,BinningAnalysis.Variance{Float64}}
    magnetizationVector::LogBinner{Vector{Float64},32,BinningAnalysis.Variance{Vector{Float64}}}
    correlation::LogBinner{Vector{Float64},32,BinningAnalysis.Variance{Vector{Float64}}}
end

function Observables(lattice::T) where T<:Lattice
    return Observables(ErrorPropagator(Float64), LogBinner(Float64), LogBinner(zeros(Float64,3)), LogBinner(zeros(Float64,lattice.length)))
end

function performMeasurements!(observables::Observables, lattice::T, energy::Float64) where T<:Lattice
    #measure energy and energy^2
    push!(observables.energy, energy / length(lattice), energy * energy / (length(lattice) * length(lattice)))

    #measure magnetization
    m = getMagnetization(lattice)
    push!(observables.magnetization, norm(m))
    push!(observables.magnetizationVector, m)

    #measure spin correlations
    push!(observables.correlation, getCorrelation(lattice))
end


text:
using OffsetArrays
using Printf # for @sprintf, etc.
#import Base.close

#=========================================================#
function out1(meqn::Int, mbc::Int, mx::Int, xlower::Float64, dx::Float64,
              q::OffsetArray{Float64,2}, t::Float64, iframe::Int,
              aux::OffsetArray{Float64}, maux::Int, outaux::Bool)
#=========================================================#
#
#     # Output the results for a general system of conservation laws
#     # in 1 dimension
#
#     # Write the results to the file fort.q<iframe>
#     # Use format required by matlab script  plotclaw1.m
#     # The same format is used by the amrclaw package.
#     # Here it's adapted to output just the single grid.
#
#     # set outaux = .true. to also output the aux arrays to fort.a<iframe>
#
#
##       implicit double precision (a-h,o-z)
##       dimension q(meqn,1-mbc:mx+mbc)
##       dimension aux(maux,1-mbc:mx+mbc)
##       character*10 fname1, fname2, fname3
##       logical outaux

#
#     # Write the results to the file fort.q<iframe>
#     # Use format required by matlab script  plotclaw1.m
#     # The same format is used by the amrclaw package.
#     # Here it's adapted to output just the single grid.
#
#     # first create the file name and open file
#
##      fname1 = 'fort.qxxxx'
##      fname2 = 'fort.txxxx'
##      fname3 = 'fort.axxxx'
##      nstp = iframe
##      do 55 ipos = 10, 7, -1
##         idigit = mod(nstp,10)
##         fname1(ipos:ipos) = char(ichar('0') + idigit)
##         fname2(ipos:ipos) = char(ichar('0') + idigit)
##         fname3(ipos:ipos) = char(ichar('0') + idigit)
##         nstp = nstp / 10
## 55   continue
##

    fname1 = Printf.@sprintf("fort.q%04d", iframe)
    fname2 = Printf.@sprintf("fort.t%04d", iframe)
    fname3 = Printf.@sprintf("fort.a%04d", iframe)

##      open(unit=50,file=fname1,status='unknown',form='formatted')
##      open(unit=60,file=fname2,status='unknown',form='formatted')

    unit50 = open(fname1, "w")
    unit60 = open(fname2, "w")

    # the following parameters are used in amrclaw where there are
    # multiple grids.  Here they are all set to 1:
    ngrids = 1
    mptr = 1
    level = 1

    # write(50,1001) mptr,level,mx

## 1001 format(i5,'                 grid_number',/,
##     &     i5,'                 AMR_level',/,
##     &     i5,'                 mx')
##

    @printf(unit50,"%5d                 grid_number\n",mptr)
    @printf(unit50,"%5d                 AMR_level\n",level)
    @printf(unit50,"%5d                 mx\n",mx)

##      write(50,1002) xlower,dx
## 1002 format(e18.8,'    xlow', /,
##     &     e18.8,'    dx', /)
###

    @printf(unit50,"%18.8e    xlow\n",xlower)
    @printf(unit50,"%18.8e    dx\n",dx)

##      do 10 i=1,mx
##         do m=1,meqn
###     # exponents with more than 2 digits cause problems reading
###     # into matlab... reset tiny values to zero:
##            if (dabs(q(m,i)) .lt. 1d-99) q(m,i) = 0.d0
##         enddo
###
##         write(50,1005) (q(m,i), m=1,meqn)
## 1005    format(4e16.8)
###
## 10   continue
##      write(50,*) ' '
##      write(50,*) ' '
##

    for i=1:mx
        for m=1:meqn
            if abs(q[m,i]) < 1e-99
                q[m,i] = 0.0
            end
            @printf(unit50,"%16.8e",q[m,i])
            # FIXME
            ## if i % 4 == 1 && i > 1
            ##     @printf(unit50,"\n")
            ## end
        end
        @printf(unit50,"\n")
    end
    @printf(unit50," \n \n")

##
##      if (outaux) then
###     # also output the aux arrays:
##         open(unit=70,file=fname3,status='unknown',form='formatted')
##         write(70,1001) mptr,level,mx
##         write(70,1002) xlower,dx
##         do 110 i=1,mx
##            do m=1,maux
###     # exponents with more than 2 digits cause problems reading
###     # into matlab... reset tiny values to zero:
##               if (dabs(aux(m,i)) .lt. 1d-99) aux(m,i) = 0.d0
##            enddo
###
##            write(70,1005) (aux(m,i), m=1,maux)
###
## 110     continue
##         write(70,*) ' '
##         close(unit=70)
##      endif
##
##      write(60,1000) t,meqn,ngrids,maux,1
##
## 1000 format(e26.16,'    time', /,
##     &     i5,'                 meqn'/,
##     &     i5,'                 ngrids'/,
##     &     i5,'                 maux'/,
##     &     i5,'                 ndim'/,/)
###
##
    Base.close(unit50)
    Base.close(unit60)

end


text:
# Molecular Dynamics Result Abstract Interface

"""
    MolecularDynamicsResult
Abstract type to be extended by all concrete structs representing the result of a molecular dynamics simulation.
"""
abstract type MolecularDynamicsResult end

"""
    get_system(result::MolecularDynamicsResult, t::Integer = 0)::AbstractSystem

Extract the underlying system at a particular timestep from the simulation result.

An implementer of this API should implement a method of this function for their custom result type.
"""
function get_system(result::MolecularDynamicsResult, t::Integer = 0)::AbstractSystem
    throw(UnimplementedError(:get_system, result))
end
"""
    get_time_range(result::MolecularDynamicsResult)::Vector{<:Real}

Extract the time range from the simulation result in atomic units.

An implementer of this API should implement a method of this function for their custom result type.
"""
function get_time_range(result::MolecularDynamicsResult)::Vector{<:Real}
    throw(UnimplementedError(:get_time_range, result))
end

"""
    temperature(result::MolecularDynamicsResult, t::Integer = 0)::Real

Extract the temperature of the simulation at a particular timestep from the simulation result in atomic units.

An implementer of this API should implement a method of this function for their custom result type.
"""
function temperature(result::MolecularDynamicsResult, t::Integer = 0)::Real
    throw(UnimplementedError(:temperature, result))
end
"""
    reference_temperature(result::MolecularDynamicsResult)::Union{Real,Missing}

Extract the reference temperature of the simulation from the simulation result.
If there is no thermostat with a reference temperature in this simulation, return missing.

An implementer of this API should implement a method of this function for their custom result type if it supports thermostats.
If not implmented, the default implemenation just returns missing.
"""
function reference_temperature(result::MolecularDynamicsResult)::Union{Real,Missing}
    missing
end

"""
    kinetic_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real

Extract the kinetic energy of the simulation at a particular timestep from the simulation result in atomic units.

An implementer of this API should implement a method of this function for their custom result type.
"""
function kinetic_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real
    throw(UnimplementedError(:kinetic_energy, result))
end
"""
    potential_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real

Extract the potential energy of the simulation at a particular timestep from the simulation result in atomic units.

An implementer of this API should implement a method of this function for their custom result type.
"""
function potential_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real
    throw(UnimplementedError(:potential_energy, result))
end
"""
    total_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real

Extract the total energy of the simulation at a particular timestep from the simulation result in atomic units.

The default implementation simply sums the kinetic and potential energy at the timestep.
An implementer of this API could implement a method of this function for their custom result type if it supports a more efficient way to calculate this quantity.
"""
function total_energy(result::MolecularDynamicsResult, t::Integer = 0)::Real
    kinetic_energy(result, t) + potential_energy(result, t)
end

"""
    rdf(result::MolecularDynamicsResult, sample_fraction::Float64 = 1.0)::Tuple{Vector{<:Real},Vector{<:Real}}

Calculate the radial distribution function from the simulation result.

To include only a trailing portion of the timesteps for reduced noise and faster computation, set sample_fraction to be less than 1.0.
The result is a tuple of vectors which represent the radial distances (in atomic units) and the value of the rdf at each distance respectively.
"""
function rdf(result::MolecularDynamicsResult, sample_fraction::Float64 = 1.0)::Tuple{Vector{<:Real},Vector{<:Real}}
    throw(UnimplementedError(:rdf, result))
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
using ModelObjectsLH
using Test

include("test_model.jl");

@testset "All" begin
    include("single_id_test.jl");
    include("object_id_test.jl");
    include("m_objects_test.jl");
end

# ----------


text:
# ---
# title: 1275. Find Winner on a Tic Tac Toe Game
# id: problem1275
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/find-winner-on-a-tic-tac-toe-game/description/>
# hidden: true
# ---
#
# Tic-tac-toe is played by two players _A_ and _B_ on a  _3_  x  _3_  grid.
#
# Here are the rules of Tic-Tac-Toe:
#
#   * Players take turns placing characters into empty squares (" ").
#   * The first player _A_ always places "X" characters, while the second player _B_  always places "O" characters.
#   * "X" and "O" characters are always placed into empty squares, never on filled ones.
#   * The game ends when there are 3 of the same (non-empty) character filling any row, column, or diagonal.
#   * The game also ends if all squares are non-empty.
#   * No more moves can be played if the game is over.
#
# Given an array `moves` where each element is another array of size 2
# corresponding to the row and column of the grid where they mark their
# respective character in the order in which _A_ and _B_ play.
#
# Return the winner of the game if it exists ( _A_ or _B_ ), in case the game
# ends in a draw return "Draw", if there are still movements to play return
# "Pending".
#
# You can assume that `moves` is  **valid** (It follows the rules of Tic-Tac-
# Toe), the grid is initially empty and _A_ will play **first**.
#
#
#
# **Example 1:**
#
#
#
#     Input: moves = [[0,0],[2,0],[1,1],[2,1],[2,2]]
#     Output: "A"
#     Explanation: "A" wins, he always plays first.
#     "X  "    "X  "    "X  "    "X  "    " **X**  "
#     "   " -> "   " -> " X " -> " X " -> " **X** "
#     "   "    "O  "    "O  "    "OO "    "OO **X** "
#
#
# **Example 2:**
#
#
#
#     Input: moves = [[0,0],[1,1],[0,1],[0,2],[1,0],[2,0]]
#     Output: "B"
#     Explanation: "B" wins.
#     "X  "    "X  "    "XX "    "XXO"    "XXO"    "XX **O** "
#     "   " -> " O " -> " O " -> " O " -> "XO " -> "X **O** "
#     "   "    "   "    "   "    "   "    "   "    " **O**  "
#
#
# **Example 3:**
#
#
#
#     Input: moves = [[0,0],[1,1],[2,0],[1,0],[1,2],[2,1],[0,1],[0,2],[2,2]]
#     Output: "Draw"
#     Explanation: The game ends in a draw since there are no moves to make.
#     "XXO"
#     "OOX"
#     "XOX"
#
#
# **Example 4:**
#
#
#
#     Input: moves = [[0,0],[1,1]]
#     Output: "Pending"
#     Explanation: The game has not finished yet.
#     "X  "
#     " O "
#     "   "
#
#
#
#
# **Constraints:**
#
#   * `1 <= moves.length <= 9`
#   * `moves[i].length == 2`
#   * `0 <= moves[i][j] <= 2`
#   * There are no repeated elements on `moves`.
#   * `moves` follow the rules of tic tac toe.
#
#
## @lc code=start
using LeetCode

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


text:
module TGH
# this is woefully incomplete

using LinearAlgebra
using LinearAlgebra: RealHermSymComplexHerm
using GenericLinearAlgebra: symtri!
using Test
using Quadmath
using GenericMRRR

# Note: we use Frobenius norms here

function runtest(A::RealHermSymComplexHerm{T,<:StridedMatrix}) where T
    λ, V = geigen!(copy(A))
    n = size(A,1)
    m = size(V,2)
    ntol = min(m,n)
    myeps = eps(real(T))
    if m == n
        resnorm = norm(V*Diagonal(λ)*V' - A)
    else
        resnorm = norm(Diagonal(λ) - V' * A * V)
    end
    Tnorm = max(opnorm(A,1), floatmin(real(T)))
    if resnorm < Tnorm
        d_err =  (resnorm / Tnorm) / (myeps * ntol)
    elseif Tnorm < 1
        d_err = (min(resnorm, Tnorm * ntol) / Tnorm) / (myeps * ntol)
    else
        d_err = min(resnorm / Tnorm, ntol) / (myeps * ntol)
    end
    if m == n
        v_err = norm(V * V' - I) / (myeps * ntol)
    else
        v_err = norm(V' * V - I) / (myeps * ntol)
    end
    d_err, v_err
end

for T in [Float128, Complex{Float128}, Float64, ComplexF64]
    @testset "dense $T" begin
        n = 16
        A = Hermitian(rand(T,n,n))
        de, ve = runtest(A)
        @test de < 30
        @test ve < 30
    end
end

end # module