text:
"""
    point_variables(x, icol)
    point_variables(x, icol, prescribed_conditions)

Extract u, θ, F, M for the point described by the point state variables at `icol` in
x after incorporating the prescribed conditions in `prescribed_conditions`

Note that the degrees of freedom that are not specified in `prescribed_conditions`
are used as state variables (e.g. prescribing F[2] would mean u[2] = x[icol+1])
"""
@inline function point_variables(x, icol)
    u = SVector(x[icol], x[icol+1], x[icol+2])
    θ = SVector(x[icol+3], x[icol+4], x[icol+5])
    F = zero(u)
    M = zero(θ)

    return u, θ, F, M
end

@inline function point_variables(x, icol, prescribed_conditions, force_scaling)

    value = prescribed_conditions.value
    follower = prescribed_conditions.follower

    # get degrees of freedom to fix/set
    force_dof = prescribed_conditions.force

    # get the displacement and rotations of the point
    u = SVector(ifelse(force_dof[1], x[icol  ], value[1]),
                ifelse(force_dof[2], x[icol+1], value[2]),
                ifelse(force_dof[3], x[icol+2], value[3]))
    θ = SVector(ifelse(force_dof[4], x[icol+3], value[4]),
                ifelse(force_dof[5], x[icol+4], value[5]),
                ifelse(force_dof[6], x[icol+5], value[6]))

    # get the transformation matrix for the point
    C = get_C(θ)

    # solve for the force applied at the point due to the prescribed loads
    Fp = zero(u)
    for i = 1:3
        if force_dof[i]
            # add dead force_dof
            Fp += SVector(I3[i,1], I3[i,2], I3[i,3])*value[i]
            # add follower force_dof
            Fp += SVector(C[i,1], C[i,2], C[i,3])*follower[i]
        end
    end

    # solve for the moment applied at the point due to the prescribed loads
    Mp = zero(θ)
    for i = 4:6
        if force_dof[i]
            # add dead moment
            Mp += SVector(I3[i-3,1], I3[i-3,2], I3[i-3,3])*value[i]
            # add follower moment
            Mp += SVector(C[i-3,1], C[i-3,2], C[i-3,3])*follower[i]
        end
    end

    # overwrite external forces/moments with solved for forces/moments
    F = SVector(ifelse(force_dof[1], Fp[1], x[icol  ] * force_scaling),
                ifelse(force_dof[2], Fp[2], x[icol+1] * force_scaling),
                ifelse(force_dof[3], Fp[3], x[icol+2] * force_scaling))
    M = SVector(ifelse(force_dof[4], Mp[1], x[icol+3] * force_scaling),
                ifelse(force_dof[5], Mp[2], x[icol+4] * force_scaling),
                ifelse(force_dof[6], Mp[3], x[icol+5] * force_scaling))

    return u, θ, F, M
end

"""
    insert_point_residual!(resid, irow_elem, irow_point, u, θ, F, M, side)

Modify the equilibrium and constitutive equations to account for the point
variables given by u, θ, F, M

If `irow_elem != irow_point`, assume that the equilibrium equations have already
been modified

# Arguments
 - resid: Residual vector
 - irow_elem: Row index of the first equilibrium/compatability equation for one
       side of the beam element
 - irow_point: Row index of the first equilibrium equation for the point
 - u: Displacement of the point
 - θ: Rotation of the point
 - F: External forces imposed on the point
 - M: External moments imposed on the point
 - side: Side of beam element (-1 (left) or 1 (right))
"""
@inline function insert_point_residual!(resid, irow_elem, irow_point, u, θ, F, M, side, force_scaling)

    if irow_elem == irow_point
        # add to equilibrium and compatability equations
        for i = 1:3
            resid[irow_elem+i-1] -= F[i] ./ force_scaling
            resid[irow_elem+i+2] -= M[i] ./ force_scaling
            resid[irow_elem+i+5] += side*u[i]
            resid[irow_elem+i+8] += side*θ[i]
        end
    else
        # add to compatability equations
        for i = 1:3
            resid[irow_elem+i-1] += side*u[i]
            resid[irow_elem+i+2] += side*θ[i]
        end
    end

    return resid
end

"""
    point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
        force_scaling, icol, irow_point, irow_elem1, irow_elem2)

Adds a points contributions to the residual vector

# Arguments
 - `resid`: residual vector
 - `x`: current state vector
 - `ipoint`: index of point
 - `assembly`: assembly of interconnected beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `force_scaling`: scaling parameter for forces
 - `icol`: starting index for the point's state variables
 - `irow_point`: Row index of the first equilibrium equation for the point
 - `irow_elem1`: Row index of first equation for the left side of each beam element
 - `irow_elem2`: Row index of first equation for the right side of each beam element
"""
@inline function point_residual!(resid, x, ipoint, assembly, prescribed_conditions,
    force_scaling, icol, irow_point, irow_elem1, irow_elem2)

    nelem = length(assembly.elements)

    # incorporate prescribed condition if applicable
    prescribed = haskey(prescribed_conditions, ipoint)
    if prescribed
        u, θ, F, M = point_variables(x, icol, prescribed_conditions[ipoint], force_scaling)
    else
        u, θ, F, M = point_variables(x, icol)
    end

    # search for beam elements that are connected to the specified point
    for ielem = 1:nelem
        # check left side of beam element
        if ipoint == assembly.start[ielem]
            # add to residual equations for the beam element endpoint
            side = -1
            irow_elem = irow_elem1[ielem]
            insert_point_residual!(resid, irow_elem, irow_point, u, θ, F, M, side, force_scaling)
        end
        # check right side of beam element
        if ipoint == assembly.stop[ielem]
            # add to residual equations for the beam element endpoint
            side = 1
            irow_elem = irow_elem2[ielem]
            insert_point_residual!(resid, irow_elem, irow_point, u, θ, F, M, side, force_scaling)
        end
    end

    return resid
end

"""
    point_follower_jacobians(x, icol, prescribed_conditions)

Calculate the jacobians of the follower forces/moments with respect to θ

# Arguments
 - x: Current state variable vector
 - icol: Row/Column index of the first state variable for the point
 - prescribed_conditions: Prescribed conditions for the point
"""
@inline function point_follower_jacobians(x, icol, prescribed_conditions)

    TF = eltype(x)

    # Apply time functions
    value = prescribed_conditions.value
    follower = prescribed_conditions.follower

    force_dof = prescribed_conditions.force

    θ = SVector(ifelse(force_dof[4], x[icol+3], value[4]),
                ifelse(force_dof[5], x[icol+4], value[5]),
                ifelse(force_dof[6], x[icol+5], value[6]))

    C = get_C(θ)
    C_θ1, C_θ2, C_θ3 = get_C_θ(C, θ)

    # solve for the jacobian wrt theta of the follower forces
    Fp_θ = @SMatrix zeros(TF, 3, 3)
    for i = 1:3
        if force_dof[i]
            rot_θ = @SMatrix [
                C_θ1[i,1] C_θ2[i,1] C_θ3[i,1];
                C_θ1[i,2] C_θ2[i,2] C_θ3[i,2];
                C_θ1[i,3] C_θ2[i,3] C_θ3[i,3]
                ]
            Fp_θ += rot_θ*follower[i]
        end
    end

    # solve for the jacobian wrt theta of the follower moments
    Mp_θ = @SMatrix zeros(TF, 3, 3)
    for i = 1:3
        if force_dof[i+3]
            rot_θ = @SMatrix [
                C_θ1[i,1] C_θ2[i,1] C_θ3[i,1];
                C_θ1[i,2] C_θ2[i,2] C_θ3[i,2];
                C_θ1[i,3] C_θ2[i,3] C_θ3[i,3]
                ]
            Mp_θ += rot_θ*follower[i+3]
        end
    end

    # if displacement is specified, corresponding component of follower jacobian is zero
    F1_θ = ifelse(force_dof[1], SVector(Fp_θ[1,1], Fp_θ[1,2], Fp_θ[1,3]), @SVector zeros(TF, 3))
    F2_θ = ifelse(force_dof[2], SVector(Fp_θ[2,1], Fp_θ[2,2], Fp_θ[2,3]), @SVector zeros(TF, 3))
    F3_θ = ifelse(force_dof[3], SVector(Fp_θ[3,1], Fp_θ[3,2], Fp_θ[3,3]), @SVector zeros(TF, 3))
    F_θ = vcat(F1_θ', F2_θ', F3_θ')

    M1_θ = ifelse(force_dof[4], SVector(Mp_θ[1,1], Mp_θ[1,2], Mp_θ[1,3]), @SVector zeros(TF, 3))
    M2_θ = ifelse(force_dof[5], SVector(Mp_θ[2,1], Mp_θ[2,2], Mp_θ[2,3]), @SVector zeros(TF, 3))
    M3_θ = ifelse(force_dof[6], SVector(Mp_θ[3,1], Mp_θ[3,2], Mp_θ[3,3]), @SVector zeros(TF, 3))
    M_θ = vcat(M1_θ', M2_θ', M3_θ')

    return F_θ, M_θ
end

"""
    insert_point_jacobian!(jacob, irow_elem, irow_point, icol,
        prescribed_conditions, side, F_θ, M_θ)

Modify the jacobian entries for the equilibrium and constitutive equations to
account for the point variables at icol

If irow_elem != irow_point, assume that the equilibrium equations have already been modified

# Arguments
 - jacob: Jacobian of residual vector with respect to state vectors
 - irow_elem: Row index of the first equilibrium/compatability equation for one side of the beam element
 - irow_point: Row index of the first equilibrium equation for the point
 - icol: Row/Column index of the first state variable for the point
 - prescribed_conditions: Prescribed force/displacement and moment/rotation on point
 - side: Side of beam element (-1 (left) or 1 (right))
"""
@inline function insert_point_jacobian!(jacob, irow_elem, irow_point, icol, side, prescribed_force, F_θ, M_θ, force_scaling)

    if irow_elem == irow_point
        # add jacobian entries for equilibrium and compatability equations
        for i = 1:3 # forces/displacements
            if prescribed_force[i]
                # F[i] is prescribed, u[i] is a state variable
                for j = 4:6
                    # check if θ[j-3] is a state variable
                    if prescribed_force[j]
                        # add equilibrium equation jacobian component
                        jacob[irow_elem+i-1, icol+j-1] = -F_θ[i,j-3] ./ force_scaling
                    end
                end
                # add compatability equation jacobian component
                jacob[irow_elem+i+5, icol+i-1] = side # u_u = I
            else
                # u[i] is prescribed, F[i] is a state variable

                # add equilibrium equation jacobian component
                jacob[irow_elem+i-1, icol+i-1] = -1 # F_F = I
            end
        end
        for i = 4:6 # moments/rotations
            if prescribed_force[i]
                # M[i-3] is prescribed, θ[i-3] is a state variable
                for j = 4:6
                    # check if θ[j-3] is a state variable
                    if prescribed_force[j]
                        # add equilibrium equation jacobian component
                        jacob[irow_elem+i-1, icol+j-1] = -M_θ[i-3,j-3] ./ force_scaling
                    end
                end
                # add compatability equation jacobian component
                jacob[irow_elem+i+5, icol+i-1] = side # θ_θ = I
            else
                # θ[i-3] is prescribed, M[i-3] is a state variable

                # add equilibrium equation jacobian component
                jacob[irow_elem+i-1, icol+i-1] = -1 # M_M = I
            end
        end
    else
        # add jacobian entries for compatability equations only
        for i = 1:3
            if prescribed_force[i]
                # F[i] is prescribed, u[i] is a state variable
                jacob[irow_elem+i-1, icol+i-1] = side # u_u = I
            end
        end
        for i = 4:6
            if prescribed_force[i]
                # M[i-3] is prescribed, θ[i-3] is a state variable
                jacob[irow_elem+i-1, icol+i-1] = side # θ_θ = I
            end
        end
    end

    return jacob
end

"""
    point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions, icol,
        irow_point, irow_elem1, irow_elem2)

Adds a points contributions to the residual vector

# Arguments
 - `jacob`: residual vector
 - `x`: current state vector
 - `ipoint`: index of point
 - `assembly`: assembly of interconnected beam elements
 - `prescribed_conditions`: dictionary of prescribed conditions
 - `force_scaling`: scaling parameter for forces
 - `icol`: starting index for the point's state variables
 - `irow_point`: Row index of the first equilibrium equation for the point
 - `irow_elem1`: Row index of first equation for the left side of each beam element
 - `irow_elem2`: Row index of first equation for the right side of each beam element
"""
@inline function point_jacobian!(jacob, x, ipoint, assembly, prescribed_conditions,
    force_scaling, icol, irow_point, irow_elem1, irow_elem2)

    nelem = length(assembly.elements)

    # incorporate prescribed condition if applicable
    prescribed = haskey(prescribed_conditions, ipoint)
    if prescribed
        F_θ, M_θ = point_follower_jacobians(x, icol, prescribed_conditions[ipoint])
        force_dof = prescribed_conditions[ipoint].force
    else
        F_θ = @SMatrix zeros(3,3)
        M_θ = @SMatrix zeros(3,3)
        force_dof = @SVector ones(Bool, 6)
    end

    # search for beam elements that are connected to the specified point
    for ielem = 1:nelem
        # check left side of beam element
        if ipoint == assembly.start[ielem]
            # add to residual equations for the beam element endpoint
            side = -1
            irow_elem = irow_elem1[ielem]
            jacob = insert_point_jacobian!(jacob, irow_elem, irow_point, icol, side, force_dof, F_θ, M_θ, force_scaling)
        end
        # check right side of beam element
        if ipoint == assembly.stop[ielem]
            # add to residual equations for the beam element endpoint
            side = 1
            irow_elem = irow_elem2[ielem]
            jacob = insert_point_jacobian!(jacob, irow_elem, irow_point, icol, side, force_dof, F_θ, M_θ, force_scaling)
        end
    end

    return jacob
end


text:
using BinaryBuilder, Pkg

function configure(version_offset, min_julia_version, proj_jll_version)
    # The version of this JLL is decoupled from the upstream version.
    # Whenever we package a new upstream release, we initially map its
    # version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
    # So for example version 2.6.3 would become 200.600.300.

    name = "GDAL"
    upstream_version = v"3.2.1"
    version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                            upstream_version.minor * 100 + version_offset.minor,
                            upstream_version.patch * 100 + version_offset.patch)

    # Collection of sources required to build GDAL
    sources = [
        ArchiveSource("https://github.com/OSGeo/gdal/releases/download/v$upstream_version/gdal-$upstream_version.tar.gz",
            "43d40ba940e3927e38f9e98062ff62f9fa993ceade82f26f16fab7e73edb572e"),
        DirectorySource("../bundled"),
    ]

    # Bash recipe for building across all platforms
    script = raw"""
    cd $WORKSPACE/srcdir/gdal-*/

    if [[ ${target} == *mingw* ]]; then
        export LDFLAGS="-L${libdir}"
        # Apply patch to customise PROJ library
        atomic_patch -p1 "$WORKSPACE/srcdir/patches/configure_ac_proj_libs.patch"
        autoreconf -vi
        export PROJ_LIBS="proj_7_2"
    elif [[ "${target}" == *-linux-* ]]; then
        # Hint to find libstdc++, required to link against C++ libs when using C compiler
        if [[ "${nbits}" == 32 ]]; then
            export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib"
        else
            export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
        fi
        # Use same flags also for GEOS
        atomic_patch -p1 "$WORKSPACE/srcdir/patches/geos-m4-extra-cflags.patch"
        export EXTRA_GEOS_CFLAGS="${CFLAGS}"
        if [[ "${target}" == powerpc64le-* ]]; then
            atomic_patch -p1 "$WORKSPACE/srcdir/patches/sqlite3-m4-extra-libs.patch"
            export EXTRA_GEOS_LIBS="${EXTRA_GEOS_LIBS} -lm"
            export EXTRA_SQLITE3_LIBS="-lm"
            # libpthread and libldl are needed for libgdal, so let's always use them
            export LDFLAGS="$LDFLAGS -lpthread -ldl"
        fi
        autoreconf -vi
    fi

    # Clear out `.la` files since they're often wrong and screw us up
    rm -f ${prefix}/lib/*.la

    ./configure --help
    ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
        --with-geos=${bindir}/geos-config \
        --with-proj=$prefix \
        --with-tiff=$prefix \
        --with-geotiff=$prefix \
        --with-libz=$prefix \
        --with-expat=$prefix \
        --with-zstd=$prefix \
        --with-sqlite3=$prefix \
        --with-curl=${bindir}/curl-config \
        --with-openjpeg \
        --with-python=no \
        --enable-shared \
        --disable-static
    # Make sure that some important libraries are found
    grep "HAVE_GEOS='yes'" config.log
    grep "HAVE_SQLITE='yes'" config.log
    grep "CURL_SETTING='yes'" config.log
    grep "ZSTD_SETTING='yes'" config.log
    grep "HAVE_EXPAT='yes'" config.log

    make -j${nproc}
    make install
    """

    platforms = expand_cxxstring_abis(supported_platforms())

    # The products that we will ensure are always built
    products = [
        LibraryProduct("libgdal", :libgdal),
        ExecutableProduct("gdal_contour", :gdal_contour_path),
        ExecutableProduct("gdal_grid", :gdal_grid_path),
        ExecutableProduct("gdal_rasterize", :gdal_rasterize_path),
        ExecutableProduct("gdal_translate", :gdal_translate_path),
        ExecutableProduct("gdaladdo", :gdaladdo_path),
        ExecutableProduct("gdalbuildvrt", :gdalbuildvrt_path),
        ExecutableProduct("gdaldem", :gdaldem_path),
        ExecutableProduct("gdalinfo", :gdalinfo_path),
        ExecutableProduct("gdallocationinfo", :gdallocationinfo_path),
        ExecutableProduct("gdalmanage", :gdalmanage_path),
        ExecutableProduct("gdalsrsinfo", :gdalsrsinfo_path),
        ExecutableProduct("gdaltindex", :gdaltindex_path),
        ExecutableProduct("gdaltransform", :gdaltransform_path),
        ExecutableProduct("gdalwarp", :gdalwarp_path),
        ExecutableProduct("nearblack", :nearblack_path),
        ExecutableProduct("ogr2ogr", :ogr2ogr_path),
        ExecutableProduct("ogrinfo", :ogrinfo_path),
        ExecutableProduct("ogrlineref", :ogrlineref_path),
        ExecutableProduct("ogrtindex", :ogrtindex_path),
    ]

    # Dependencies that must be installed before this package can be built
    dependencies = [
        Dependency("GEOS_jll"),
        Dependency(PackageSpec(name="PROJ_jll", version=proj_jll_version)),
        Dependency("Zlib_jll"),
        Dependency("SQLite_jll"),
        Dependency("OpenJpeg_jll"),
        Dependency("Expat_jll", v"2.2.7"; compat="~2.2.7"),
        Dependency("Zstd_jll"),
        # TODO: v4.3.0 is available, use that next time
        Dependency("Libtiff_jll"; compat="4.1.0"),
        Dependency("libgeotiff_jll"),
    ]

    jll_stdlibs = Dict(
        v"1.3" => [
            Dependency("LibCURL_jll", v"7.71.1"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll", v"1.9.0"),
            Dependency("MbedTLS_jll", v"2.16.8"),
            Dependency("nghttp2_jll", v"1.40.0"),
        ],
        v"1.6" => [
            Dependency("LibCURL_jll"),
            # The following libraries are dependencies of LibCURL_jll which is now a
            # stdlib, but the stdlib doesn't explicitly list its dependencies
            Dependency("LibSSH2_jll"),
            Dependency("MbedTLS_jll", v"2.24.0"),
            Dependency("nghttp2_jll"),
        ]
    )

    append!(dependencies, jll_stdlibs[min_julia_version])

    return name, version, sources, script, platforms, products, dependencies
end

ext: jl

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
# These are some functions to allow static type-checking of Julia programs

module TypeCheck
export checkreturntypes, checklooptypes, checkmethodcalls,
  methodswithdescendants

## Modifying functions from Base

# return the type-inferred AST for each method of a generic function
function Base.code_typed(f::Function)
  Expr[code_typed(m) for m in f.env]
end

# return the type-inferred AST for one method of a generic function
function Base.code_typed(m::Method)
 linfo = m.func.code
 (tree,ty) = Base.typeinf(linfo,m.sig,())
 if !isa(tree,Expr)
     ccall(:jl_uncompress_ast, Any, (Any,Any), linfo, tree)
 else
    tree
 end
end

# specify a method via a method, function, or function+argument-type-tuple
# prints out the local variables of that method, and their types
function Base.whos(f,args...)
  for e in code_typed(f,args...)
    display(MethodSignature(e))
    for x in sort(e.args[2][2];by=x->x[1])
      println("\t",x[1],"\t",x[2])
    end
    println("")
  end
end

## Basic Operations on Function Exprs

# given an Expr representing a method, return its inferred return type
returntype(e::Expr) =  e.args[3].typ

# given an Expr representing a method, return an Array of Exprs representing its body
body(e::Expr) = e.args[3].args

# given an Expr representing a method, return all of the return statement in its body
returns(e::Expr) = filter(x-> typeof(x) == Expr && x.head==:return,body(e))

# given an Expr representing a method,
# return all function all Exprs contained in return statements in the method body
function extractcallsfromreturns(e::Expr)
  rs = returns(e)
  rs_with_calls = filter(x->typeof(x.args[1]) == Expr && x.args[1].head == :call,rs)
  Expr[expr.args[1] for expr in rs_with_calls]
end

# A type that covers both Types and TypeVars
AType = Union(Type,TypeVar)

# given an Expr representing a method, get the type of each of the arguments in the signature
function argumenttypes(e::Expr)
  argnames = Symbol[typeof(x) == Symbol ? x : x.args[1] for x in e.args[1]]
  argtuples = filter(x->x[1] in argnames, e.args[2][2]) #only arguments, no local vars
  AType[t[2] for t in argtuples]
end

# given an Expr, determine if it is calling a TopNode
# (this affects how we should handle resolving the callee name)
istop(e) = Base.is_expr(e,:call) && typeof(e.args[1]) == TopNode

# given a call Expr (:call, :call1, :new), determine its return type
function returntype(e::Expr,context::Expr) #must be :call,:new,:call1
  if Base.is_expr(e,:new); return e.typ; end
  if Base.is_expr(e,:call1) && isa(e.args[1], TopNode); return e.typ; end
  if !Base.is_expr(e,:call); error("Expected :call Expr"); end

  if istop(e)
    return e.typ
  end

  callee = e.args[1]
  if istop(callee)
    return returntype(callee,context)
  elseif isa(callee,SymbolNode) # only seen (func::F), so non-generic function
    return Any
  elseif is(callee,Symbol)
    if e.typ != Any || any([isa(x,LambdaStaticData) for x in e.args[2:end]])
      return e.typ
    end

    if isdefined(Base,callee)
      f = eval(Base,callee)
      if !isa(f,Function) || !isgeneric(f)
        return e.typ
      end
      fargtypes = tuple([argumenttype(ea,context) for ea in e.args[2:end]])
      return Union([returntype(ef) for ef in code_typed(f,fargtypes)]...)
    else
      return @show e.typ
    end
  end

  return e.typ
end

# for an Expr `e` used as an argument to a call in function Expr `context`,
# determine the type of `e`
function argumenttype(e::Expr,context::Expr)
 if Base.is_expr(e,:call) || Base.is_expr(e,:new) || Base.is_expr(e,:call1)
   return returntype(e,context)
 end

 @show e
 return Any
end

# for a Symbol `s` used as an argument to a call in a function Expr `e`,
# determine the type of `s`.
function argumenttype(s::Symbol,e::Expr)
  vartypes = [x[1] => x[2] for x in e.args[2][2]]
  s in vartypes ? (vartypes[@show s]) : Any
end

# as above, but for different call argument types
argumenttype(s::SymbolNode,e::Expr) = s.typ
argumenttype(t::TopNode,e::Expr) = Any
argumenttype(l::LambdaStaticData,e::Expr) = Function
argumenttype(q::QuoteNode,e::Expr) = argumenttype(q.value,e)

# as above, but for various literal values
argumenttype(n::Number,e::Expr) = typeof(n)
argumenttype(c::Char,e::Expr) = typeof(c)
argumenttype(s::String,e::Expr) = typeof(s)
argumenttype(i,e::Expr) = typeof(i) #catch all, hopefully for more literals

# start, next, and done are the functions for-loops use to iterate
# having implementations of them makes a type iterable
# this defines iterating over a DataType to mean iterating over its descendants
# (breadth-first-search ordering)
Base.start(t::DataType) = [t]
function Base.next(t::DataType,arr::Vector{DataType})
  c = pop!(arr)
  append!(arr,[x for x in subtypes(c)])
  (c,arr)
end
Base.done(t::DataType,arr::Vector{DataType}) = length(arr) == 0

# this is like the `methodswith` function from Base,
# but it considers all descendants of the type
# rather than the type itself (or the type + super types)
# it produces a list of function names and statistics
# about how many descendants implement that function
# this is good for discovering implicit interfaces
function methodswithdescendants(t::DataType;onlyleaves::Bool=false,lim::Int=10)
  d = Dict{Symbol,Int}()
  count = 0
  for s in t
    if !onlyleaves || (onlyleaves && isleaftype(s))
      count += 1
      fs = Set{Symbol}()
      for m in methodswith(s)
        push!(fs,m.func.code.name)
      end
      for sym in fs
        d[sym] = get(d,sym,0) + 1
      end
    end
  end
  l = [(k,v/count) for (k,v) in d]
  sort!(l,by=(x->x[2]),rev=true)
  l[1:min(lim,end)]
end

# check all the generic functions in a module
function checkallmodule(m::Module;test=checkreturntypes,kwargs...)
  score = 0
  for n in names(m)
    f = eval(m,n)
    if isgeneric(f) && typeof(f) == Function
      fm = test(f;mod=m,kwargs...)
      score += length(fm.methods)
      display(fm)
    end
  end
  println("The total number of failed methods in $m is $score")
end

# use checkallmodule to implement the Module version of other checks
checkreturntypes(m::Module;kwargs...) = checkallmodule(m;test=checkreturntypes,kwargs...)
checklooptypes(m::Module) = checkallmodule(m;test=checklooptypes)
checkmethodcalls(m::Module) = checkallmodule(m;test=checkmethodcalls)

## Checking that return values are base only on input *types*, not values.

type MethodSignature
  typs::Vector{AType}
  returntype::Union(Type,TypeVar) # v0.2 has TypeVars as returntypes; v0.3 does not
end
MethodSignature(e::Expr) = MethodSignature(argumenttypes(e),returntype(e))
function Base.writemime(io, ::MIME"text/plain", x::MethodSignature)
  println(io,"(",join([string(t) for t in x.typs],","),")::",x.returntype)
end

type FunctionSignature
  methods::Vector{MethodSignature}
  name::Symbol
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionSignature)
  for m in x.methods
    print(io,string(x.name))
    display(m)
  end
end

# given a function, run checkreturntypes on each method
function checkreturntypes(f::Function;kwargs...)
  results = MethodSignature[]
  for e in code_typed(f)
    (ms,b) = checkreturntype(e;kwargs...)
    if b push!(results,ms) end
  end
  FunctionSignature(results,f.env.name)
end

# given an Expr representing a Method,
# determine whether its return type is based
# only on the arugment types or whether it is
# also influenced by argument values
# (the Method fails the check if the return type depends on values)
function checkreturntype(e::Expr;kwargs...)
  (typ,b) = isreturnbasedonvalues(e;kwargs...)
  (MethodSignature(argumenttypes(e),typ),b)
end

# Determine whether this method's return type might change based on input values rather than input types
function isreturnbasedonvalues(e::Expr;mod=Base)
  rt = returntype(e)
  ts = argumenttypes(e)
  if isleaftype(rt) || rt == None return (rt,false) end

  for t in ts
    if !isleaftype(t)
      return (rt,false)
    end
  end

  cs = [returntype(c,e) for c in extractcallsfromreturns(e)]
  for c in cs
    if rt == c
       return (rt,false)
    end
  end

  return (rt,true) # return is not concrete type; all args are concrete types
end

## Checking that variables in loops have concrete types

type LoopResult
  msig::MethodSignature
  lines::Vector{(Symbol,Type)} #TODO should this be a specialized type? SymbolNode?
  LoopResult(ms::MethodSignature,ls::Vector{(Symbol,Type)}) = new(ms,unique(ls))
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResult)
  display(x.msig)
  for (s,t) in x.lines
    println(io,"\t",string(s),"::",string(t))
  end
end

type LoopResults
  name::Symbol
  methods::Vector{LoopResult}
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResults)
  for lr in x.methods
    print(io,string(x.name))
    display(lr)
  end
end

# for a given Function, run checklooptypes on each Method
function checklooptypes(f::Function;kwargs...)
  lrs = LoopResult[]
  for e in code_typed(f)
    lr = checklooptypes(e)
    if length(lr.lines) > 0 push!(lrs,lr) end
  end
  LoopResults(f.env.name,lrs)
end

# for an Expr representing a Method,
# check that the type of each variable used in a loop
# has a concrete type
checklooptypes(e::Expr;kwargs...) = LoopResult(MethodSignature(e),loosetypes(loopcontents(e)))

# This is a function for trying to detect loops in the body of a Method
# Returns lines that are inside one or more loops
function loopcontents(e::Expr)
  b = body(e)
  loops = Int[]
  nesting = 0
  lines = {}
  for i in 1:length(b)
    if typeof(b[i]) == LabelNode
      l = b[i].label
      jumpback = findnext(
        x-> (typeof(x) == GotoNode && x.label == l) || (Base.is_expr(x,:gotoifnot) && x.args[end] == l),
        b, i)
      if jumpback != 0
        push!(loops,jumpback)
        nesting += 1
      end
    end
    if nesting > 0
      push!(lines,(i,b[i]))
    end

    if typeof(b[i]) == GotoNode && in(i,loops)
      splice!(loops,findfirst(loops,i))
      nesting -= 1
    end
  end
  lines
end

# given `lr`, a Vector of expressions (Expr + literals, etc)
# try to find all occurances of a variables in `lr`
# and determine their types
# `method` is only used as part of constructing the return type
function loosetypes(lr::Vector)
  lines = (Symbol,Type)[]
  for (i,e) in lr
    if typeof(e) == Expr
      es = copy(e.args)
      while !isempty(es)
        e1 = pop!(es)
        if typeof(e1) == Expr
          append!(es,e1.args)
        elseif typeof(e1) == SymbolNode && !isleaftype(e1.typ) && typeof(e1.typ) == UnionType
          push!(lines,(e1.name,e1.typ))
        end
      end
    end
  end
  return lines
end

## Check method calls

type CallSignature
  name::Symbol
  argumenttypes::Vector{AType}
end
function Base.writemime(io, ::MIME"text/plain", x::CallSignature)
  println(io,string(x.name),"(",join([string(t) for t in x.argumenttypes],","),")")
end

type MethodCalls
  m::MethodSignature
  calls::Vector{CallSignature}
end

function Base.writemime(io, ::MIME"text/plain", x::MethodCalls)
  display(x.m)
  for c in x.calls
    print(io,"\t")
    display(c)
  end
end

type FunctionCalls
  name::Symbol
  methods::Vector{MethodCalls}
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionCalls)
  for mc in x.methods
    print(io,string(x.name))
    display(mc)
  end
end

# given a Function, run `checkmethodcalls` on each Method
# and collect the results into a FunctionCalls
function checkmethodcalls(f::Function;kwargs...)
  calls = MethodCalls[]
  for m in f.env
    e = code_typed(m)
    mc = checkmethodcalls(e,m;kwargs...)
    if !isempty(mc.calls)
      push!(calls, mc)
    end
  end
  FunctionCalls(f.env.name,calls)
end

# given an Expr representing a Method,
# and the Method it represents,
# check the Method body for calls to non-existant Methods
function checkmethodcalls(e::Expr,m::Method;kwargs...)
  if Base.arg_decl_parts(m)[3] == symbol("deprecated.jl")
    CallSignature[]
  end
  nomethoderrors(e,methodcalls(e);kwargs...)
end

# Find any methods that match the given CallSignature
function hasmatches(mod::Module,cs::CallSignature)
  if isdefined(mod,cs.name)
    f = eval(mod,cs.name)
    if isgeneric(f)
      opts = methods(f,tuple(cs.argumenttypes...))
      if isempty(opts)
        return false
      end
    end
  else
    #println("$mod.$(cs.name) is undefined")
  end
  return true
end

# Find any CallSignatures that indicate potential NoMethodErrors
function nomethoderrors(e::Expr,cs::Vector{CallSignature};mod=Base)
  output = CallSignature[]
  for callsig in cs
    if !hasmatches(mod,callsig)
      push!(output,callsig)
    end
  end
  MethodCalls(MethodSignature(e),output)
end

# Look through the body of the function for `:call`s
function methodcalls(e::Expr)
  b = body(e)
  lines = CallSignature[]
  for s in b
    if typeof(s) == Expr
      if s.head == :return
        append!(b, s.args)
      elseif s.head == :call
        if typeof(s.args[1]) == Symbol
          push!(lines,CallSignature(s.args[1], [argumenttype(e1,e) for e1 in s.args[2:end]]))
        end
      end
    end
  end
  lines
end

end  #end module


text:
using StringUtils

function main()
    product = 0
    for a = 1:998
        for b = (a+1):999
            c = 1000 - a - b
            if c > a && c > b && c^2 == a^2 + b^2
                @printf("(a,b,c) = (%d, %d, %d)\n", a, b, c)
                product = a * b * c
            end
        end
    end

    @printf("The Pythagorean triplet product is %d\n", product)
end

main()


text:
#========================================================================================#
#	Laboratory 2
#
# Programs, types and functions
#
# Author: Niall Palfreyman, 04/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we look at the very heart of Julia: How it uses types to dispatch
		functions - that is, how Julia decides which particular implementation method
		it will use to calculate the value of a function call.

		First let's start with the structure of a Julia program. Every Julia program
		starts off as a string:

		prog = "8 / (2 + 3.0)"

		Enter this mini-program, then use Meta.parse() to analyse the program into an expression
		'expr' and tell me the typeof() expr:
		""",
		"expr = Meta.parse(prog)",
		x -> (x==Expr)
	),
	Activity(
		"""
		An Expr is a structure with two parts: a head (expr.head) and a vector of arguments
		(expr.args). What is the head of expr?
		""",
		"expr.head",
		x -> x==:call
	),
	Activity(
		"""
		Use the function dump() to view the entire structure of expr. As you can see,
		the .args field of expr is a Vector. Use your experience of accessing vector elements
		to find out the typeof the second argument of the third argument of expr:
		""",
		"dump(expr) or expr.args[3].args[2] - Try both! :-)",
		x -> x==Int64
	),
	Activity(
		"""
		And of course we can execute (i.e., evaluate) our expression using the function eval().
		What is the value of expr?
		""",
		"eval(expr)",
		x -> x==1.6
	),
	Activity(
		"""
		An expression is a tree that can be evaluated. We call such an evaluation tree a
		METHOD: a single, particular implementation of some function we wish to evaluate.
		Let's look at a more complicated evaluation tree. Create the following method for
		calculating the factorial of an Integer value:

		fact(n::Integer) = Meta.parse("(\$n ≤ 1) ? 1 : \$n * eval(fact(\$(n-1)))")

		Notice how we're using string interpolation here to insert the value of the
		argument n into an Expr tree structure that we can study. Use dump, .head and .args
		to study the tree structure of fact(5), and then evaluate that tree ...
		""",
		"eval(fact(5))",
		x -> x==120
	),
	Activity(
		"""
		Our fact() method for Integers is working fine, however we also have a problem:
		What if we wished to calculate the factorial of the floating-point value 5.0? Try
		doing this now to see that it does not work...

		This problem actually makes perfect sense. For example, if we allow float values,
		how would we go about calculating fact(5.01)? To calculate the factorial of a float
		value, we need to use the Gamma function gamma(). First load the SpecialFunctions library:

		using SpecialFunctions

		Now experiment with the gamma function to find out its relationship to factorial.
		For example, what value of n yields the result gamma(n) == fact(5)?
		""",
		"n == 6",
		x -> x==6
	),
	Activity(
		"""
		Now write a new factorial method for Real numbers:

		fact(n::Real) = Meta.parse("gamma(\$(n+1))")

		Now use dump, .head and .args to study the tree structure of fact(5.0), and then
		evaluate that tree ...
		""",
		"eval(fact(5.0))",
		x -> x==120.0
	),
	Activity(
		"""
		What is the value of fact(5.01)?
		""",
		"eval(fact(5.01))",
		x -> round(x,digits=2) == 122.07
	),
	Activity(
		"""
		Looking at the method tree for fact(5) and for fact(5.0), you can see that
		Julia is using the type of the function's arguments to decide which method to
		use for calculating the function fact(). This decision-process is called
		DISPATCHING, and it is extremely useful!

		Our first fact() method is quicker and simpler in cases where the argument name
		is an Integer, whereas the second method is more effective in cases where n is
		a Real number. The Dispatcher automatically chooses which method is needed for
		the various argument types.

		Use methods() to see a list of the different methods you have written for the
		function fact(). If you are unsure how to use the methods() function, type '?'
		at the Julia prompt and then enter "methods".

		What is the typeof the return result of methods()?
		""",
		"methods(fact)",
		x -> x == Base.MethodList
	),
	Activity(
		"""
		Since types are so important for dispatching functions, let's try creating our
		own user-defined types. First create a few ABSTRACT types to represent various
		biological organisms:

		abstract type Organism end
		abstract type Animal <: Organism end

		Now use the supertype() function to find the supertype of Animal. What type is
		returned if you ask for the subtypes() of Organism?
		""",
		"subtypes(Organism)",
		x -> (x <: Vector)
	),
	Activity(
		"""
		Next create a CONCRETE subtype of our abstract Animal type:

		struct Weasel <: Animal
			name::String
			weight::Integer
			female::Bool
		end

		Use the function fieldnames() to inspect the individual fields of the Weasel
		type, then give me the descriptor of the third field
		""",
		"fieldnames(Weasel)[3]",
		x -> x == :female
	),
	Activity(
		"""
		We use struct types to instantiate concrete OBJECTs in computer memory by
		entering specific values for the individual fields of the Weasel struct:

		wendy = Weasel( "Wendy", 101, true)
		willy = Weasel( "Willy", 115, false)

		Notice that types start with an UPpercase letter, whereas objects start with a
		lowercase letter. By default, Julia creates structs as IMMUTABLE - that is, we
		cannot modify the value of their fields. This means Julia can always know exactly
		how much memory an object needs, which provides major performance advantages. Try
		changing the value of Wendy's gender to see that this is not allowed. What word
		does the resulting exception message use to describe the struct you tried to modify?
		""",
		"wendy.female = false",
		x -> lowercase(strip(x)) == "immutable"
	),
	Activity(
		"""
		Immutability might at first seem awkward, but very often we don't want to change
		an object's fields, but instead want to replace one object by another. This way
		of using objects is far faster and less error-prone. If, however, we do want to
		change the fields of a type, we must define it as MUTABLE, like this:

		mutable struct Rabbit <: Animal
			name::String
			length::Integer
		end

		rabia = Rabbit( "Rabia", 27)

		Change Rabia's length to 29 cm. What is the supertype of her type?
		""",
		"rabia.length = 29",
		x -> (x == Main.Animal)
	),
	Activity(
		"""
		Unlike object-oriented languages, which only dispatch on the first argument of a
		function, a major design feature of Julia is that it uses MULTIPLE DISPATCHING,
		which is particularly important when we use Julia to perform biological simulations.
		To see multiple dispatching in action, let's use our above type definitions.

		When Animals meet each other, they react in different ways according to their type:
		Weasels challenge each other, but they attack Rabbits. We could implement these
		different interactions using if-else conditionals, but it is easier to use
		multiple dispatching. Enter the following definitions at the Julia prompt:

		meet( meeter::Weasel, meetee::Rabbit) = "attacks"
		meet( meeter::Weasel, meetee::Weasel) = "challenges"
		meet( meeter::Rabbit, meetee::Rabbit) = "sniffs"
		meet( meeter::Rabbit, meetee::Weasel) = "hides"
		meet( meeter::Organism, meetee::Organism) = "ignores"

		Test these definitions by finding out what happens when Rabia meets Wendy:
		""",
		"meet(rabia,wendy)",
		x -> occursin( "hide", lowercase(x))
	),
	Activity(
		"""
		The dispatcher must make these decisions during the execution time of our
		simulation program. Enter the following function definition:

		function encounter( meeter::Organism, meetee::Organism)
			println( meeter.name, " meets ", meetee.name, " and ", meet(meeter,meetee), ".")
		end

		Now test this function by calling encounter() with various combinations of
		Wendy, Willy and Rabia.

		What happens if you create a new Rabbit called Robby, and Rabia encounters him?
		""",
		"meet(robby,rabia)",
		x -> occursin( "sniff", lowercase(x))
	),
	Activity(
		"""
		Now, to see the full power of multiple dispatching, add your own new concrete
		type and then check how an encounter between your type and Rabia works out. Do
		it something like this:

		struct Tree <: Organism; name::String end
		tilly = Tree( "Tilly")

		How does Rabia react to Tilly?
		""",
		"encounter(rabia,tilly)",
		x -> occursin( "ignores", lowercase(x))
	),
	Activity(
		"""
		And now one final Activity for you: Can you add a new type of Organism called
		Grass, and arrange for Rabia to eat it? It is possible to do this in just 3-4
		new lines of code.
		""",
		"encounter(meeter::Rabbit,meetee::Grass) = \"eats\"",
		x -> occursin( "eats", lowercase(x))
	),
]


text:
export initialize

function initialize(model;groups=[],expr="",parallel=false)
  var=match(r"^.*(?==)",expr).match
  value=match(r"(?<==).*",expr).match
  if !contains(model.initialization,var)
    N=model.groups[:N]
    model.initialization*="$var=SharedArray(Float64,N)\n"
  end
  if groups==[] || groups=="all"
    model.initialization*="for i=1:N $var[i]=$value end\n"
  else
    for name in [groups;]
    i=model.groups[name][:localindex]
    ii="(i-$(i[1]-1))"
    value=replace(value,"i",ii)
    model.initialization*="for i=$i $var[i]=$value end\n"
    end
  end
end


text:
"""
    rasterize_points(points, dx::AbstractFloat)

Rasterize points in 2D by a cell size `dx`.
Returns a dictionary containing the indices points that are in a cell.
"""
immutable Raster{T<:SVector,U<:Integer}
    pixels::Dict{Tuple{U,U}, Vector{U}}
    r_min::T
    r_max::T
    cellsize::Float64
end

# from: https://github.com/andyferris/HeightOrderedGrid.jl/blob/master/src/HeightOrderedGrid.jl
# as opposed to: SVector{N, Float64}(minimum(map(x->x[1], points)), minimum(map(x->x[2], points)), 0)
function bounds{T}(points::Vector{T})
    xmin = typemax(eltype(eltype(T)))
    xmax = typemin(eltype(eltype(T)))
    ymin = typemax(eltype(eltype(T)))
    ymax = typemin(eltype(eltype(T)))

    for p ∈ points
        @inbounds x = p[1]
        if x < xmin
            xmin = x
        end
        if x > xmax
            xmax = x
        end

        @inbounds y = p[2]
        if y < ymin
            ymin = y
        end
        if y > ymax
            ymax = y
        end
    end

    return (xmin, xmax, ymin, ymax)
end

"""
`rasterize_points{T <: AbstractVector}(points::Vector{T}, dx::AbstractFloat) -> Raster`

Returns a `Raster` of `points` with quadratic cellsize `dx`.
"""
function rasterize_points{T <: AbstractVector}(points::Vector{T}, dx::AbstractFloat)
    (xmin, xmax, ymin, ymax) = bounds(points)

    if size(eltype(points))[1] == 2
        min_xy = eltype(points)(xmin, ymin)
        max_xy = eltype(points)(xmax, ymax)
    elseif size(eltype(points))[1] == 3
        min_xy = eltype(points)(xmin, ymin, 0)
        max_xy = eltype(points)(xmax, ymax, 0)
    else
        throw(DimensionMismatch("Unsupported input point dimensions"))
    end

    pixels = Dict{Tuple{UInt32, UInt32}, Vector{UInt32}}()
    inv_dx = 1.0/dx
    for i = 1:length(points)
        @inbounds p = points[i] - min_xy
        key = (floor(UInt32, p[1]*inv_dx), floor(UInt32, p[2]*inv_dx))
        if haskey(pixels, key)
            push!(pixels[key], i)
        else
            pixels[key] = Vector{UInt32}()
            push!(pixels[key], i)
        end
    end
    return Raster(pixels,min_xy,max_xy,dx)
end

function rasterize_points{T <: Number}(points::Matrix{T}, dx::AbstractFloat)
    @assert isbits(T)
    ndim = size(points, 1)
    npoints = size(points, 2)
    new_data = reinterpret(SVector{ndim, T}, points, (npoints,))
    rasterize_points(new_data, dx)
end

Base.keys(r::Raster) = keys(r.pixels)
Base.values(r::Raster) = values(r.pixels)
Base.getindex(r::Raster, ind) = r.pixels[ind]

function Base.show{T,U}(io::IO, raster::Raster{T,U})
    println(io, typeof(raster))
    println(io, "  Number of pixels: ", length(raster.pixels))
    print(io, "  Cellsize: ", raster.cellsize)
end


text:
using Test, LinearMaps, LinearAlgebra

function left_tester(L::LinearMap{T}) where {T}
    M, N = size(L)
    A = Matrix(L)

    x = rand(T, N)
    y = rand(T, M)

    # *
    @test y' * A ≈ y' * L
    @test (y' * A) * x ≈ y' * (L * x) # associative
    @test transpose(y) * A ≈ transpose(y) * L

    # mul!
    α = rand(T); β = rand(T)
    b1 = y' * A
    b2 = similar(b1)
    bt = copy(b1')'
    bm = Matrix(bt)
    @test mul!(b2, y', L) ≈ mul!(bt, y', L) ≈ mul!(bm, y', L)
    @test mul!(b2, y', L) === b2
    @test b1 ≈ b2 ≈ bt
    @test mul!(copy(b2), y', L, α, β) ≈ b2*β + y'A*α

    b1 = transpose(y) * A
    b2 = similar(b1)
    bt = transpose(copy(transpose(b1)))
    @test mul!(b2, transpose(y), L) ≈ mul!(bt, transpose(y), L)
    @test b1 ≈ b2 ≈ bt
    @test mul!(copy(b2), transpose(y), L, α, β) ≈ b2*β + transpose(y)*A*α

    Y = rand(T, M, 3)
    X = similar(Y, 3, N)
    Xt = copy(X')'
    @test Y'L isa LinearMap
    @test Matrix(Y'L) ≈ Y'A
    @test mul!(X, Y', L) ≈ mul!(Xt, Y', L) ≈ Y'A
    @test mul!(Xt, Y', L) === Xt
    @test mul!(copy(X), Y', L, α, β) ≈ X*β + Y'A*α
    @test mul!(X, Y', L) === X
    @test mul!(X, Y', L, α, β) === X

    @test transpose(Y)*L isa LinearMap
    @test Matrix(transpose(Y)*L) ≈ transpose(Y)*A
    @test mul!(X, transpose(Y), L) ≈ mul!(Xt, transpose(Y), L) ≈ transpose(Y)*A
    @test mul!(Xt, transpose(Y), L) === Xt
    @test mul!(copy(X), transpose(Y), L, α, β) ≈ X*β + transpose(Y)*A*α
    @test mul!(X, transpose(Y), L) === X
    @test mul!(X, transpose(Y), L, α, β) === X

    return true
end

@testset "left multiplication" begin
    T = ComplexF32
    N = 5
    L = LinearMap{T}(cumsum, reverse ∘ cumsum ∘ reverse, N)

    @test left_tester(L) # FunctionMap
    @test left_tester(adjoint(L)) # AdjointMap
    @test left_tester(transpose(L)) # TransposeMap
    @test left_tester(L'*L) # CompositeMap
    @test left_tester(2rand(T)*L) # ScaledMap
    @test left_tester(kron(L, L')) # KroneckerMap
    @test left_tester(2rand(T)*L + 3rand(T)*L') # LinearCombination
    @test left_tester([L L]) # BlockMap

    W = LinearMap(randn(T,5,4)) # WrappedMap
    @test left_tester(W)
    @test left_tester(W')
    @test left_tester(transpose(W))

    J = LinearMap(1.0I, 5) # UniformScalingMap
    @test left_tester(J)
end


text:
using Random

import Base.Order: ForwardOrdering, ReverseOrdering, Forward, Reverse
import BioSimulator: enqueue!, dequeue!, heapify!

@testset "PQBinaryHeap" begin
  # Test dequeing in sorted order.
  function test_issorted!(pq::PQBinaryHeap, rev=false)
    (_, last_val) = dequeue!(pq)
    while !isempty(pq)
      (_, next_val) = dequeue!(pq)
      if !rev
        @test last_val <= next_val
      else
        @test next_val <= last_val
      end
      last_val = next_val
    end
  end

  # TODO...
  function test_isrequested!(pq::PQBinaryHeap, keys)
    i = 0
    while !isempty(pq)
      krqst    = keys[i+=1]
      krcvd, _ = dequeue!(pq, krqst)
      @test krcvd == krqst
    end
  end

  # create the PQBinaryHeap but do not build the heap
  # function init_maxpq(keys, values)
  #   K = eltype(keys)
  #   V = eltype(values)

  #   keymap = collect(1:length(keys))
  #   valmap = copy(keymap)

  #   return PQBinaryHeap{K,V,ForwardOrdering}(keys, values, keymap, valmap)
  # end

  # function init_minpq(keys, values)
  #   K = eltype(keys)
  #   V = eltype(values)

  #   keymap = collect(1:length(keys))
  #   valmap = copy(keymap)

  #   return PQBinaryHeap{K,V,ReverseOrdering}(keys, values, keymap, valmap)
  # end

  # create the PQBinaryHeap by building the heap
  function build_maxpq(keys, values)
    K = eltype(keys)
    V = eltype(values)

    return PQBinaryHeap{K,V,ReverseOrdering}(keys, values)
  end

  function build_minpq(keys, values)
    K = eltype(keys)
    V = eltype(values)

    return PQBinaryHeap{K,V,ForwardOrdering}(keys, values)
  end

  @testset "Low-Level Heap Operations" begin
    # pmax = 1000
    # n = 10000

    # keys, values = collect(1:n), rand(1:pmax, n)

    # @testset "heapify!" begin
    #   maxpq = init_maxpq(keys, values)
    #   minpq = init_minpq(keys, values)

    #   heapify!(maxpq)
    #   heapify!(minpq)

    #   @test issorted([dequeue!(maxpq)[2] for _ in length(values)])
    #   @test issorted([dequeue!(minpq)[2] for _ in length(values)])
    # end

    # TODO: What else should be in here?
  end

  @testset "Constructors" begin
    pmax = 1000
    n = 10000

    keys, values = collect(1:n), rand(1:pmax, n)

    @testset "build from parallel arrays" begin
      test_issorted!(build_minpq(keys, values), false)
      test_issorted!(build_maxpq(keys, values), true)
    end

    @testset "different array sizes => ArgumentError" begin
      @test_throws ArgumentError PQBinaryHeap(['a'], [1, 2])
      @test_throws ArgumentError PQBinaryHeap(['a', 'b'], [1])
    end

    @testset "duplicate key => ArgumentError" begin
      @test_throws ArgumentError PQBinaryHeap(['a', 'a'], [1, 2])
    end
  end

  @testset "Methods" begin
    pmax = 1000
    n = 10000

    easy_sorted_keys,   easy_sorted_values   = ['a', 'b'],   [1, 2]
    easy_unsorted_keys, easy_unsorted_values = ['b', 'a'],   [2, 1]

    easy_cases = [
      "sorted/sorted"     => (easy_sorted_keys,   easy_sorted_values),
      "unsorted/sorted"   => (easy_unsorted_keys, easy_sorted_values),
      "sorted/unsorted"   => (easy_sorted_keys,   easy_unsorted_values),
      "unsorted/unsorted" => (easy_unsorted_keys, easy_unsorted_values)
    ]

    random_unsorted_keys, random_unsorted_values = shuffle(1:n),               rand(1:pmax, n)
    random_sorted_keys,   random_sorted_values   = sort(random_unsorted_keys), sort(random_unsorted_values)

    random_cases = [
      "sorted/sorted"     => (random_sorted_keys,   random_sorted_values),
      "unsorted/sorted"   => (random_unsorted_keys, random_sorted_values),
      "sorted/unsorted"   => (random_sorted_keys,   random_unsorted_values),
      "unsorted/unsorted" => (random_unsorted_keys, random_unsorted_values)
    ]

    examples = ["easy" => easy_cases, "random" => random_cases]
    # examples = ["easy" => easy_cases]

    @testset "$(exampletype) cases" for (exampletype, example) in examples
      @testset "$(case)" for (case, problem) in example
        (keys, values) = problem

        minpq = build_minpq(keys, values)
        maxpq = build_maxpq(keys, values)

        # we can only check for the correct value, key mappings will be stored in a non-deterministic fashion
        @testset "peektop" begin
          minval, minidx = findmin(values)
          maxval, maxidx = findmax(values)

          @test peektop(minpq)[2] == minval
          @test peektop(maxpq)[2] == maxval
        end

        @testset "get/getindex" begin
          if exampletype == "easy"
            # key does not exist
            @test get(minpq, 'c', 0) == 0
            @test get(maxpq, 'c', 0) == 0

            @test_throws KeyError minpq['c']
            @test_throws KeyError maxpq['c']

            # key exists
            a_value = values[something(findfirst(isequal('a'), keys), 0)]

            @test get(minpq, 'a', 0) == a_value
            @test get(maxpq, 'a', 0) == a_value
            @test minpq['a'] == a_value
            @test maxpq['a'] == a_value
          else
            # key does not exist
            @test get(minpq, n + 1, 0) == 0
            @test get(maxpq, n + 1, 0) == 0

            @test_throws KeyError minpq[n + 1]
            @test_throws KeyError maxpq[n + 1]

            # key exists
            random_key   = rand(keys)
            random_value = values[something(findfirst(isequal(random_key), keys), 0)]

            @test get(minpq, random_key, 0) == random_value
            @test get(maxpq, random_key, 0) == random_value
            @test minpq[random_key] == random_value
            @test maxpq[random_key] == random_value
          end
        end

        @testset "enqueue!" begin
          K, V = eltype(keys), eltype(values)

          minpq_copy = PQBinaryHeap{K,V}(Forward)
          maxpq_copy = PQBinaryHeap{K,V}(Reverse)

          for i in eachindex(keys)
            key   = keys[i]
            value = values[i]

            enqueue!(minpq_copy, value, key)
            enqueue!(maxpq_copy, value, key)
          end

          test_issorted!(minpq_copy, false)
          test_issorted!(maxpq_copy, true)
        end

        # again, we can only test the value returned by dequeue!
        @testset "dequeue!" begin
          minpq_copy = deepcopy(minpq)
          maxpq_copy = deepcopy(maxpq)

          sorted_order = sortperm(values)

          for i in eachindex(sorted_order)
            ifwd = sorted_order[i]
            irev = sorted_order[end - i + 1]

            nextmin = values[ifwd]
            nextmax = values[irev]

            @test dequeue!(minpq_copy)[2] == nextmin
            @test dequeue!(maxpq_copy)[2] == nextmax
          end
        end

        @testset "setindex!" begin
          @testset "changing priorities" begin
            if exampletype == "easy"
              minpq_copy = deepcopy(minpq)
              maxpq_copy = deepcopy(maxpq)

              minpq_copy['a'] = 100
              maxpq_copy['a'] = 100

              # make sure we do not add a duplicate key
              @test length(minpq_copy) == length(minpq)
              @test length(maxpq_copy) == length(maxpq)

              # was the value updated correctly?
              @test minpq_copy['a'] == 100
              @test maxpq_copy['a'] == 100
            end
          end

          @testset "adding a new item" begin

          end
        end

        # remember, key mappings are non-deterministic
        @testset "Iteration" begin
          K, V = eltype(keys), eltype(values)

          minpq_copy = PQBinaryHeap{K,V}(Forward)
          maxpq_copy = PQBinaryHeap{K,V}(Reverse)

          for (key, value) in minpq
            enqueue!(minpq_copy, value, key)
          end

          @test minpq.keys   == minpq_copy.keys
          @test minpq.values[minpq.keymap] == minpq_copy.values[minpq_copy.keymap]

          for (key, value) in maxpq
            enqueue!(maxpq_copy, value, key)
          end

          @test maxpq.keys == maxpq_copy.keys
          @test maxpq.values[maxpq.keymap] == maxpq_copy.values[maxpq_copy.keymap]
        end
      end
    end
  end
end


text:
module ColorButtons

Base.Experimental.@optlevel 1

using Redux
using CImGui
using ..Buttons
import ..Buttons: AbstractButtonAction, AbstractButtonState, get_label
import CImGui: ImVec4

# actions
const Rename = Buttons.Rename
const SetTriggeredTo = Buttons.SetTriggeredTo
const Resize = Buttons.Resize
const ChangeWidth = Buttons.ChangeWidth
const ChangeHeight = Buttons.ChangeHeight

"""
    SetButtonColorTo(label, h, s, v)
    SetButtonColorTo(label, color::CImGui.ImVec4)
Change button's color to `color`.
"""
struct SetButtonColorTo <: AbstractButtonAction
    label::String
    color::ImVec4
end
SetButtonColorTo(label, h, s, v) = SetButtonColorTo(label, CImGui.HSV(h, s, v))

"""
    SetHoveredColorTo(label, h, s, v)
    SetHoveredColorTo(label, color::ImVec4)
Change hovered color to `color`.
"""
struct SetHoveredColorTo <: AbstractButtonAction
    label::String
    color::ImVec4
end
SetHoveredColorTo(label, h, s, v) = SetHoveredColorTo(label, CImGui.HSV(h, s, v))

"""
    SetActiveColorTo(label, h, s, v)
    SetActiveColorTo(label, color::ImVec4)
Change active color to `color`.
"""
struct SetActiveColorTo <: AbstractButtonAction
    label::String
    color::ImVec4
end
SetActiveColorTo(label, h, s, v) = SetActiveColorTo(label, CImGui.HSV(h, s, v))

"""
    AddButton(label::AbstractString, size = (0,0))
Add a [`ColorButton`](@ref).
"""
struct AddButton <: AbstractButtonAction
    button::Buttons.State
    button_color::ImVec4
    hovered_color::ImVec4
    active_color::ImVec4
end
AddButton(label::AbstractString, size) = AddButton(
    Buttons.State(label, size),
    CImGui.HSV(4 / 7.0, 0.6, 0.6),
    CImGui.HSV(4 / 7.0, 0.7, 0.7),
    CImGui.HSV(4 / 7.0, 0.8, 0.8),
)
AddButton(label::AbstractString) = AddButton(label, (0, 0))

get_label(a::AddButton) = a.button.label

"""
    DeleteButton(label::AbstractString)
Delete the [`ColorButton`](@ref).
"""
struct DeleteButton <: AbstractButtonAction
    label::String
end

# state
"""
    ColorButtons.State(label::AbstractString)
    ColorButtons.State(label::AbstractString, size)
    ColorButtons.State(label::AbstractString, size, button_color, hovered_color, active_color)
A color button state which contains a common button plus three colors.
"""
struct State <: AbstractButtonState
    button::Buttons.State
    button_color::ImVec4
    hovered_color::ImVec4
    active_color::ImVec4
end
State(label::AbstractString, size) = State(
    Buttons.State(label, size),
    CImGui.HSV(4 / 7.0, 0.6, 0.6),
    CImGui.HSV(4 / 7.0, 0.7, 0.7),
    CImGui.HSV(4 / 7.0, 0.8, 0.8),
)
State(label::AbstractString) = State(label, (0, 0))

get_label(s::State) = s.button.label
get_size(s::State) = s.button.size
get_color(s::State) = s.button_color
get_hovered_color(s::State) = s.hovered_color
get_active_color(s::State) = s.active_color
is_triggered(s::State) = s.button.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::State, a::AbstractButtonAction) =
    State(Buttons.reducer(s.button, a), s.button_color, s.hovered_color, s.active_color)
reducer(s::State, a::SetButtonColorTo) = State(s.button, a.color, s.hovered_color, s.active_color)
reducer(s::State, a::SetHoveredColorTo) = State(s.button, s.button_color, a.color, s.active_color)
reducer(s::State, a::SetActiveColorTo) = State(s.button, s.button_color, s.hovered_color, a.color)

reducer(s::Dict{String,<:AbstractButtonState}, a::AbstractButtonAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractButtonState}, a::AbstractButtonAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end
reducer(s::Vector{<:AbstractButtonState}, a::AddButton) =
    [s..., State(a.button, a.button_color, a.hovered_color, a.active_color)]
reducer(s::Vector{<:AbstractButtonState}, a::DeleteButton) =
    filter(s -> get_label(s) !== get_label(a), s)

# helper
"""
    ColorButton(store::AbstractStore, get_state=Redux.get_state) -> Bool
A wrapper over [`Buttons.Button`](@ref) with additional color states.
"""
function ColorButton(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    CImGui.PushStyleColor(CImGui.ImGuiCol_Button, get_color(s))
    CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonHovered, get_hovered_color(s))
    CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonActive, get_active_color(s))
    is_triggered = CImGui.Button(get_label(s), get_size(s))
    CImGui.PopStyleColor(3)
    dispatch!(store, Buttons.SetTriggeredTo(get_label(s), is_triggered))
    return is_triggered
end

end # module


text:
module Autotrack

using Dates, LinearAlgebra, Statistics
using VideoIO, ColorVectorSpace, Dierckx, ImageTransformations, ImageDraw, StaticArrays, PaddedViews, OffsetArrays, ImageFiltering, ColorTypes
SV = SVector{2, Float64}

export track

# include("video.jl")

window_radius = 4
nframes = 25 # data points
σ = 0.85
smoothing_factor = 200
# sz = (108, 144)
scale = 10

w = CartesianIndex(window_radius, window_radius)
window = -w:w

# _parse_time(::Nothing, T) = 0.0
# _parse_time(x, T) = Millisecond(T(parse(Int, x)))/Millisecond(Second(1))
# _parse_time(xT) = _parse_time(xT...)
#
# function _get_seconds(txt)
#   m = match(r"(?:(\d+):)*(\d+):(\d+)(?:,(\d+))*", txt)
#   sum(_parse_time, zip(m.captures, (Hour, Minute, Second, Millisecond)))
# end

function seekread!(img, vid, t)
  seek(vid, t)
  read!(vid, img)
end

# function get_times(start::Time, stop::Time)
#   t1 = _get_seconds(start)
#   t2 = _get_seconds(stop)
#   range(t1, t2, length = nframes)
# end

function get_next(guess, img, bkgd)
  centered_img = OffsetArrays.centered(img, Tuple(guess))
  centered_bkgd = OffsetArrays.centered(bkgd, Tuple(guess))
  x = centered_bkgd[window] .- centered_img[window]
  imfilter!(x, x, Kernel.DoG(σ))
  _, i = findmax(x)
  guess + window[i]
end

function get_imgs(vid, ts)
  img = read(vid)
  # t₀ = gettime(vid)
  h, w = size(img)
  width_ind = 1:scale:w
  height_ind = 1:scale:h
  sz = (length(height_ind), length(width_ind))

  unpadded_imgs = [similar(img, sz) for _ in 1:nframes]
  for (i, t) in enumerate(ts)
    seekread!(img, vid, t)
    # seekread!(img, vid, t + t₀)
    unpadded_imgs[i] .= img[height_ind, width_ind]
  end

  height, width = sz
  padded_axes = (1-window_radius:height+window_radius, 1-window_radius:width+window_radius)
  imgs = PaddedView.(zero(eltype(img)), unpadded_imgs, Ref(padded_axes))
  bkgd = PaddedView(zero(eltype(img)), mean(unpadded_imgs), padded_axes)

  return sz, imgs, bkgd
end

function get_spline(sz, imgs, bkgd, ts)
  guess = CartesianIndex(sz .÷ 2)
  coords = accumulate((guess, img) -> get_next(guess, img, bkgd), imgs, init = guess)
  ParametricSpline(ts, hcat(SV.(Tuple.(coords))...); s = smoothing_factor, k = 2)
end

function track(start_file, start_time, stop_file, stop_time; debug = nothing)

  if start_file == stop_file
    vid = VideoIO.openvideo(start_file, target_format=VideoIO.AV_PIX_FMT_GRAY8)
    ts = range(start_time, stop_time, length = nframes)
    sz, imgs, bkgd = get_imgs(vid, ts)
    spl = get_spline(sz, imgs, bkgd, ts)
  else
    @error "waaaat"
  end

  ar = VideoIO.aspect_ratio(vid)
  if !isnothing(debug)
    save(debug, ts, imgs, spl, ar)
  end

  return ts[1], ts[end], spl, scale*[1, ar]
end

# function get_data(videofile, start, stop)
#
#   vid = VideoIO.openvideo(videofile, target_format=VideoIO.AV_PIX_FMT_GRAY8)
#   ts = get_times(start, stop)
#   sz, imgs, bkgd = get_imgs(vid, ts)
#   spl = get_spline(sz, imgs, bkgd, ts)
#
#   return ts, imgs, spl, VideoIO.aspect_ratio(vid)
# end

function save(result_file, ts, imgs, spl, aspect_ratio)
  # ar = SV(1, aspect_ratio)
  # txy = (NamedTuple{(:t, :x, :y)}((t, ar .* spl(t)...)) for t in ts[1] : 1/10 : ts[end])
  # # txy = (NamedTuple{(:t, :x, :y)}((t, spl(t)...)) for t in ts[1] : 1/10 : ts[end])
  # CSV.write("$result_file.csv", txy)

  encoder_options = (color_range=2, crf=23, preset="medium")
  time = 2
  framerate = round(Int, length(ts)/time)
  open_video_out("$result_file.mp4", parent(imgs[1]), framerate=framerate, encoder_options=encoder_options) do writer
    for (i, t) in enumerate(ts)
      path = Path([CartesianIndex(Tuple(round.(Int, spl(i)))) for i in ts[1] : 1/10 : t])
      img = parent(imgs[i])
      draw!(img, path)
      write(writer, img)
    end
  end

end

# function process_track(video_file, start, stop, result_file)
#   data = get_data(video_file, start, stop)
#   save(result_file, data...)
# end
#
# function _process_tracks(rows, f, video_dir, result_dir)
#   @showprogress for (row, (file, start, stop)) in enumerate(f)
#     process_track(joinpath(video_dir, file), start, stop, joinpath(result_dir, string(row)))
#   end
# end
#
# function process_tracks(csvfile; video_dir = dirname(csvfile), result_dir = dirname(csvfile))
#   f = CSV.File(csvfile)
#   rows = 1:length(f)
#   _process_tracks(rows, f, video_dir, result_dir)
# end
#
# function process_tracks(csvfile, row::Int; video_dir = dirname(csvfile), result_dir = dirname(csvfile))
#   f = CSV.File(csvfile)
#   file, start, stop = f[row]
#   process_track(joinpath(video_dir, file), start, stop, joinpath(result_dir, string(row)))
# end
#
# function process_tracks(csvfile, rows::AbstractVector{Int}; video_dir = dirname(csvfile), result_dir = dirname(csvfile))
#   f = CSV.File(csvfile)
#   _process_tracks(rows, f, video_dir, result_dir)
# end

end


text:
abstract type IndicesArray{T,N,P,I} <: AbstractArray{T,N} end

###
### Interface
###
parent_type(::T) where {T<:AbstractArray} = axes_type(T)
parent_type(::Type{<:IndicesArray{T,N,P,I}}) where {T,N,P,I} = P

axes_type(::T) where {T<:AbstractArray} = axes_type(T)
axes_type(::Type{<:IndicesArray{T,N,P,I}}) where {T,N,P,I} = I

Base.parentindices(x::IndicesArray) = axes(parent(x))

Base.length(a::IndicesArray) = axes2length(axes(a))

Base.size(a::IndicesArray) = axes2size(axes(a))

@propagate_inbounds Base.size(a::IndicesArray, i) = length(axes(a, i))

Base.axes(a::IndicesArray, i) = unsafe_axes(axes(a), to_dims(a, i))
unsafe_axes(axs::Tuple, idx::Int) = axs[idx]

IndexCore.indices(x::IndicesArray) = axes(x)

"""
    IndicesMatrix{T,I1,I2}

Alias for two-dimensional `IndicesArray`  with elements of type `T`, parent of
type `P`, and index of type `I1` and `I2`. Alias for [`IndicesArray{T,1}`](@ref).
"""
const IndicesMatrix{T,P,I1,I2} = IndicesArray{T,2,P,Tuple{I1,I2}}

"""
    IndicesVector{T,I}

Alias for one-dimensional `IndicesArray`  with elements of type `T`, parent of
type `P`, and index of type `I`. Alias for [`IndicesArray{T,1}`](@ref).
"""
const IndicesVector{T,P,I1} = IndicesArray{T,1,P,Tuple{I1}}

"""
    IVecOrMat{T}

Union type of [`IndicesVector{T}`](@ref) and [`IndicesMatrix{T}`](@ref).
"""
const IndicesVecOrMat{T} = Union{IndicesMatrix{T},IndicesVector{T}}

const IndicesAdjoint{T,P,I1,I2} = IndicesMatrix{T,Adjoint{T,P},I1,I2}

const IndicesTranspose{T,P,I1,I2} = IndicesMatrix{T,Transpose{T,P},I1,I2}

const IndicesDiagonal{T,V,I1,I2} = IndicesMatrix{T,Diagonal{T,V},I1,I2}

const IndicesQRCompactWY{T,P,I1,I2} = LinearAlgebra.QRCompactWY{T,IndicesMatrix{T,P,I1,I2}}

const IndicesQRPivoted{T,P,I1,I2} = QRPivoted{T,IndicesMatrix{T,P,I1,I2}}

const IndicesQR{T,P,I1,I2} = QR{T,IndicesMatrix{T,P,I1,I2}}

const IndicesQRUnion{T,P,I1,I2} = Union{IndicesQRCompactWY{T,P,I1,I2},
                                        IndicesQRPivoted{T,P,I1,I2},
                                        IndicesQR{T,P,I1,I2}}

_maybe_indices_array(x, inds::Tuple) = IndicesArray(x, inds, AllUnique, LengthChecked)
_maybe_indices_array(x, inds::Tuple{}) = x

function Base.empty!(a::IndicesArray)
    empty!(axes(a, 1))
    empty!(parent(a))
    return a
end
###
### Utilities
###
_catch_empty(x::Tuple) = x
_catch_empty(x::NamedTuple) = x
_catch_empty(::Tuple{}) = nothing
_catch_empty(::NamedTuple{(),Tuple{}}) = nothing

StaticRanges.Size(::Type{<:IndicesArray{T,N,P,I}}) where {T,N,P,I} = StaticRanges._Size(I)


text:
function sp_value_from_primitive_values(x...)
    #=
    assume the values provided are x1,x2...
    generates the list of values for e1,e2...
    For example, x=1,y=2,z=3
    Then e1=x+y+z = 6,e2=xy+yz+xz = 11, e3=xyz = 6
    =#
    length(x) < 1 && return []
    p = [1,x[1]]
    for i = 2:length(x)
        tmp = [0,x[i]*p...]
        tmp[1:length(p)] += p
        p=tmp
    end
    return p[2:end]
end

function evaluate_monomial(sp_values,m::semi_elementary_monomial)
    if !is_elementary(m)
        throw(ArgumentError("currently only support evaluating elementary monomials"))
    elseif length(sp_values) != dim(m)
        throw(DimensionMismatch("first argument must have same length"))
    end
    return prod([sp_values[i]^m.powers[i] for i=1:length(sp_values)])
end

function evaluate_polynomial(sp_values,m::semi_elementary_polynomial)
    return sum([evaluate_monomial(sp_values,k)*v for (k,v) in m.terms])
end

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
using NamedDims
using NamedDims: names
using Test
using Tracker
using Tracker: data, TrackedReal

@testset "Tracker Compat" begin

    tw = TrackedArray(2ones(3,5))
    tb = TrackedArray(5ones(3))

    nx = NamedDimsArray{(:features, :obs)}(ones(5,10))
    ny = NamedDimsArray{(:targets, :obs)}(ones(3,10))

    y = tw*nx .+ tb

    @test data(y) == 15ones(3,10)
    @test dimnames(y) == (:_, :obs)
    @test dimnames(data(y)) == (:_, :obs)

    loss = sum((y .- ny).^2)
    @test loss isa TrackedReal

    Tracker.back!(loss)
    @test size(Tracker.grad(tw)) == (3,5)
    @test all(!isequal(0.0), Tracker.grad(tw))
    @test size(Tracker.grad(tb)) == (3,)
    @test all(!isequal(0.0), Tracker.grad(tb))
end


text:
using HypothesisTests, Base.Test

@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:both), 1.0)
@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:both), 1.0)
@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:left), 0.05)
@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:right), 0.05)

srand(12345)
@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:both), 1.0)
@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:both), 1.0)
@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:left), 0.08)
@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:right), 0.07)


text:
# Univariate bcdModel type and constructors.
type bcdPar
  val::Float64
  grid::Vector{Float64}
  lik::Vector{Float64}
  lik_max_index::Int64
  function bcdPar(val::Float64, grid::Vector{Float64})
    self = new()
    self.val = val
    self.grid = grid
    lik = NaN*zeros(length(self.grid))
    lik_max_index = 0
    return self
  end
end

type bcdModel
  # Class Attribute
    name::String # Series name
    vals::Vector{Float64}#Raw series
    dte::Vector{Date} # series dates in Date format
    q::Int64 # number of cosine transforms
    T::Int64
    avg::Float64
    lrstd::Float64
    std::Float64
    cost::Vector{Float64,} # low frequency cosine transforms
    scost::Array{Float64,} # standardized low frequency cosine transforms
    proj::Vector{Float64} #
    Omega::Array{Float64,} # Sample covariance matrix in (qxq) form
    #parameterization
    b::bcdPar
    c::bcdPar
    d::bcdPar
    Sigma::Array{Float64,} # qΣbcd model-specific covariance matrix. Note this depends only on parameters b, c, d, and on cosine transform weights cost.
    stuf # dict to dump other shit
  # Class methods
    getCov::Function
    mle_d::Function
  # bcdModel Constructor
  function bcdModel(vals::Array{Float64,1}; q=12::Int64, b=0.0::Float64, c=5.0::Float64, d=0.3::Float64, name="Series"::String, start_dte=1900, nc = 5, nd = 50, cov_calc=false)
    #####################
    # Class constructor #
    #####################
    self = new()
    self.name = name
    self.vals = vals
    self.dte = collect(Date(start_dte):Dates.Year(1):Date(start_dte+length(vals)-1))
    self.q = q
    self.T = length(self.vals)
    self.avg = mean(self.vals)
    self.std = std(self.vals)
    #self.cost = NaN*zeros(q)
      Ψ = psi_compute(self.T, self.q)
    self.cost = Ψ'*(self.vals - self.avg)
    self.lrstd = sqrt(mean(self.cost.^2)*length(self.vals))
    self.scost = (self.cost)./sqrt(self.cost'*self.cost)
    #self.scost = NaN*similar(self.cost)
    #self.proj = NaN*similar(self.vals)
    self.proj = self.avg + (Ψ/(Ψ'Ψ))*self.cost
    self.Omega = NaN*zeros(q, q)
    self.b = bcdPar(b, [0.])
      (nc <= 1) && (nc = nc + 1)
      c_grid = zeros(nc)
      for ic in 2:nc
        c_grid[ic] = 2.*(200)^(2.^(ic./(nc .- 1.0)))
      end
    self.c = bcdPar(c, c_grid)
      (nd <= 1) && (nd = nd + 1)
      d_grid = collect(linspace(-0.49, 1.49, nd))
    self.d = bcdPar(d, d_grid)
    self.stuf = Dict()
    #################
    # Class methods #
    #################
    self.getCov = function()
      Sigma_grid = q_Sigma_cd_grid(self.q, self.c.grid, self.d.grid)
      Sigma = interp_q_Sigma_cd_grid(Sigma_grid, self.c.val, self.d.val, self.c.grid, self.d.grid)
      self.Sigma = Sigma
      return self
    end
    #
    self.mle_d = function()
      Sigma_d_mat = NaN*zeros(self.q, self.q, length(self.d.grid))
      dcount = 0.
      self.d.lik = NaN*similar(self.d.grid)
      for dlik in 1:length(self.d.grid)
        tic
        Sigma = Sigma_Compute2(0.0, 0.0, self.d.grid[dlik], self.q)
        Sigma_inv = inv(Sigma)
        self.d.lik[dlik] = den_invariant(self.scost, Sigma_inv)[1]
#        Sigma_d_mat[:, :, dlik] = Sigma
#        Sigma_inv_d_mat[:, :, dlik] = Sigma_inv
        tic()
      end
      tmp, dmax_index = findmax(self.d.lik)
      self.d.val = self.d.grid[dmax_index]
      self.d.lik_max_index = dmax_index
      println("d_MLE: $(self.d.val) (estimated from an equally-spaced grid of $(length(self.d.grid)) points between -0.49 and 1.49 for the fractional difference model, i.e. bcd(0, 0, 1))")
      return self
    end
    return self
  end
end


text:
"""
$(SIGNATURES)

Find markers for indented lines (i.e. a line return followed by a tab or 4
spaces).
"""
function find_indented_blocks!(tokens::Vector{Token}, st::String)::Nothing
    # index of the line return tokens
    lr_idx = [j for j in eachindex(tokens) if tokens[j].name == :LINE_RETURN]
    remove = Int[]
    # go over all line return tokens; if they are followed by either four spaces
    # or by a tab, then check if the line is empty or looks like a list, otherwise
    # change the token for a LR_INDENT token which will be captured as part of code
    # blocks.
    for i in 1:length(lr_idx)-1
        # capture start and finish of the line (from line return to line return)
        start  = from(tokens[lr_idx[i]])   # first :LINE_RETURN
        finish = from(tokens[lr_idx[i+1]]) # next  :LINE_RETURN
        line   = subs(st, start, finish)
        indent = ""
        if startswith(line, "\n    ")
            indent = "    "
        elseif startswith(line, "\n\t")
            indent = "\t"
        else
            continue
        end
        # is there something on that line? if so, does it start with a list indicator
        # like `*`, `-`, `+` or [0-9](.|\)) ? in which case this takes precedence (commonmark)
        code_line = subs(st, nextind(st, start+length(indent)), prevind(st, finish))
        scl       = strip(code_line)
        isempty(scl) && continue
        # if here, it looks like an indented block
        tokens[lr_idx[i]] = Token(:LR_INDENT, subs(st, start, start+length(indent)), i+1)
    end
    return nothing
end

"""
$(SIGNATURES)

In initial phase, discard all `:LR_INDENT` that don't meet the requirements for
an indented code block. I.e.: where the first line is not preceded by a blank
line and then subsequently followed by indented lines.
"""
function filter_lr_indent!(vt::Vector{Token}, s::String)::Nothing
    tind_idx = [i for i in eachindex(vt) if vt[i].name == :LR_INDENT]
    isempty(tind_idx) && return nothing
    disable = zeros(Bool, length(tind_idx))
    open = false
    i    = 1
    while i <= length(tind_idx)
        t = vt[tind_idx[i]]
        if !open
            # check if previous line is blank
            prev_idx = prevind(s, from(t))
            if prev_idx < 1 || s[prev_idx] != '\n'
                disable[i] = true
            else
                open = true
            end
            i += 1
        else
            # if touching previous line, keep open and go to next
            prev = vt[tind_idx[i-1]]
            if prev.lno + 1 == t.lno
                i += 1
            else
                # otherwise close, and continue (will go to !open)
                open = false
                continue
            end
        end
    end
    for i in eachindex(disable)
        disable[i] || continue
        vt[tind_idx[i]] = Token(:LINE_RETURN, subs(vt[tind_idx[i]].ss, 1))
    end
    return nothing
end

"""
$SIGNATURES

When two indented code blocks follow each other and there's nothing in between
(empty line(s)), merge them into a super block.
"""
function merge_indented_blocks!(blocks::Vector{OCBlock}, mds::AS)::Nothing
    # indices of CODE_BLOCK_IND
    idx = [i for i in eachindex(blocks) if blocks[i].name == :CODE_BLOCK_IND]
    isempty(idx) && return
    # check if they're separated by something or nothing
    inter_space = [(subs(mds, to(blocks[idx[i]]), from(blocks[idx[i+1]])) |> strip |> length) > 0
                    for i in 1:length(idx)-1]

    curseq     = Int[] # to keep track of current list of blocks to merge
    del_blocks = Int[] # to keep track of blocks that will be removed afterwards

    # if there's no inter_space, add to the list, if there is, close and merge
    for i in eachindex(inter_space)
        if inter_space[i] && !isempty(curseq)
            # close and merge all in curseq and empty curseq
            form_super_block!(blocks, idx, curseq, del_blocks)
        elseif !inter_space[i]
            push!(curseq, i)
        end
    end
    !isempty(curseq) && form_super_block!(blocks, idx, curseq, del_blocks)
    # remove the blocks that have been merged
    deleteat!(blocks, del_blocks)
    return
end

"""
$SIGNATURES

Helper function to [`merge_indented_blocks`](@ref).
"""
function form_super_block!(blocks::Vector{OCBlock}, idx::Vector{Int},
                           curseq::Vector{Int}, del_blocks::Vector{Int})::Nothing
    push!(curseq, curseq[end]+1)
    first_block = blocks[idx[curseq[1]]]
    last_block  = blocks[idx[curseq[end]]]
    # replace the first block with the super block
    blocks[idx[curseq[1]]] = OCBlock(:CODE_BLOCK_IND, (otok(first_block) => ctok(last_block)))
    # append all blocks but the first to the delete list
    append!(del_blocks, curseq[2:end])
    empty!(curseq)
    return
end

"""
$SIGNATURES

Discard any indented block that is within a larger block to avoid ambiguities
(see #285).
"""
function filter_indented_blocks!(blocks::Vector{OCBlock})::Nothing
    # retrieve the indices of the indented blocks
    idx = [i for i in eachindex(blocks) if blocks[i].name == :CODE_BLOCK_IND]
    isempty(idx) && return nothing
    # keep track of the ones that are active
    active    = ones(Bool, length(idx))
    # retrieve their span
    indblocks = blocks[idx]
    froms     = indblocks .|> from
    tos       = indblocks .|> to
    # retrieve the max/min span for faster processing
    minfrom   = froms |> minimum
    maxto     = tos |> maximum
    update    = false
    # go over all blocks and check if they contain an indented block
    for block in blocks
        # discard if the block is before or after all indented blocks
        from_, to_ = from(block), to(block)
        (to_ < minfrom || from_ > maxto) && continue
        # otherwise check and deactivate if it's contained
        for k in eachindex(active)
            active[k] || continue
            indblock = blocks[idx[k]]
            if from_ < from(indblock) && to_ > to(indblock)
                active[k] = false
                update = true
            end
        end
        if update
            froms[.!active] .= typemax(Int)
            tos[.!active]   .= 0
            minfrom          = froms |> minimum
            maxto            = tos |> maximum
            update           = false
        end
    end
    deleteat!(blocks, idx[.!active])
    return nothing
end


text:
using ExoplanetsSysSim

start_dir = pwd()
data_dir = joinpath(dirname(pathof(ExoplanetsSysSim)),"..","data")
pkg_dev_home = joinpath(dirname(pathof(ExoplanetsSysSim)),"..")

cd(pkg_dev_home)
if true
   if !isfile(".gitmodules")
      println("# Adding git@github.com:ExoJulia/SysSimData.git a submodule in the data directory...")
      flush(stdout)
      run(`git submodule add git@github.com:ExoJulia/SysSimData.git data`)
   end
   println("# Initializing the data submodule...")
   flush(stdout)
   run(`git submodule init`)
   run(`git submodule update`)
   cd("data")
   println("# Initializing data's submodules...")
   flush(stdout)
   run(`git submodule init`)
   println("# git lfs pull just in case binary files not downloaded already...")
   flush(stdout)
   run(`git lfs pull`)
   cd("..")
end
println("# Recursively updating submodules...")
flush(stdout)
run(`git submodule update --recursive`)

if data_dir != joinpath(pkg_dev_home,"data")
  cd(pkg_dev_home)
  symlink(data_dir,"data")
end

cd(start_dir)


text:
"""
    process_session(filepath::String)

Given a raw data file location it computes pokes, bouts and streaks dataframe adding informations about daily variables:
[:Session, :Day, :Turn, :MouseID, :Gen]
"""
function process_session(filepath::String)
    filename = splitpath(filepath)[end]
    Session = replace(filename, ".csv" => "")
    Day = Date("20" * match(r"\d{6}",filename).match, "yyyymmdd")
    Turn = Session[end]
    MouseID = match(r"[a-zA-Z]{1,2}\d{1,2}",filename).match
    Gen = gen(MouseID)
    p = process_pokes(filepath)
    b = process_bouts(p)
    s = process_streaks(p)
    for df in (p,b,s)
        for (n,v) in zip([:Session, :Day, :Turn, :MouseID, :Gen],[Session, Day, Turn, MouseID, Gen])
            df[!,n] .= v
        end
    end
    return p,b,s
end


text:
using BenchmarkTools
using Distributions
import Cubature
import HCubature
import QuadGK

p(x::Float64) = exp(-0.5 * x^2) / √(2π)
p1(x) = p(x[1])
p2(xy) = p(xy[1]) * p(xy[2])
p3(xyz) = p(xyz[1]) * p(xyz[2]) * p(xyz[3])
p4(xyzw) = p(xyzw[1]) * p(xyzw[2]) * p(xyzw[3]) * p(xyzw[4])

function mybench(m::Module)
    @info "start bench $m"
    # 1 次元用のルーチン
    @btime $m.hquadrature(p, -10, 10)[1]
    # 多次元用のルーチンで試す
    @btime $m.hcubature(p1, [-10], [10])[1]
    @btime $m.hcubature(p2, [-10, -10], [10, 10])[1]
    @btime $m.hcubature(p3, [-10, -10, -10], [10, 10, 10])[1]
    @btime $m.hcubature(p4, [-10, -10, -10, -10], [10, 10, 10, 10])[1]
    @info "end bench $m"
end

mybench(Cubature)
mybench(HCubature)
@info "Appendix QuadGK"
# one-dimensional numerical integration
@btime QuadGK.quadgk(p, -10, 10)[1]

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
math : keeps track of a few mathy functions.
"""

module math

using StatsBase

function mylogsumexp(b)
	B = maximum(b, dims=2)
	return log.(sum(exp.(broadcast(-, b, B)), dims=2)) + B

end

function abs_tol(older, newer, thresh=0.0002)
	if abs(older - newer) < thresh
		return true
	else
		return false
	end
end

function rel_tol(older, newer, thresh=0.0004)
	fun = abs(1 - (older / newer))
	if fun < thresh
		return true
	else
		return false
	end
end

function softmax(vec, num_samples::Int=1, temperature::Float64=1.0)
	vec = vec ./ temperature
	max_value = maximum(vec)
	log_partition = log.(sum(exp.(vec .- max_value))) .+ max_value

	prob_vec = exp.(vec .- log_partition)

	wv = Weights(prob_vec)
	result = sample(1:length(prob_vec), wv, num_samples)
	return result
end

end


text:
using BenchmarkTools
using SFGTools

println("Listing Spectra:")
df = @btime list_spectra(inexact="twinscan", group=true)

println("Loading 7 Spectra (tiff):")
data = @btime load_spectra($df[:id][1])
println("Sizes: $(size(data[1]))\n")

println("Loading single Spectrum (tiff)")
data = @btime load_spectra($df[:id][1][1])
println("Size: $(size(data))\n")

println("Event Removal:")
@btime rm_events!.($data)

println("Get Metadata:")
@btime get_metadata($data[1])