text:
# import Particles package
using Particles

# parameters
Δt = 0.1    # time step size
nt = 100    # number of time steps
ns = 5      # number of samples

# random initial conditions
x₀ = [ [rand() .* 2 .- 1, rand() .* 2 .- 1, 0, rand() .- 0.5, rand() .- 0.5, 1] for i in 1:ns ]

# electric field
E(x::Vector{DT}) where {DT} = DT[0, 0, cos(2π*x[3])]

# magnetic field
B(x::Vector{DT}) where {DT} = DT[0, 0, 1]

# vector field
function lorentz_force!(ż, z)
    for i in eachindex(ż, z)
        x = z[i][1:3]
        v = z[i][4:6]

        e = E(x)
        b = B(x)

        ż[i][1] = v[1]
        ż[i][2] = v[2]
        ż[i][3] = v[3]

        ż[i][4] = e[1] + v[2] * b[3] - v[3] * b[2]
        ż[i][5] = e[2] + v[3] * b[1] - v[1] * b[3]
        ż[i][6] = e[3] + v[1] * b[2] - v[2] * b[1]
    end
end

# create an Equation instance
equ = Equation(lorentz_force!, x₀)

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt)

# run simulation
x = sim()

# load Plots package
# using Plots

# set plot ranges
# xlim = (-2, +2)
# ylim = (-2, +2)
# zlim = ( 0, 20)

# create empty 3d plot with ns empty series
# plt = plot3d(
#     ns,
#     xlim = xlim,
#     ylim = ylim,
#     zlim = zlim,
#     title = "Charged Particles",
#     legend = false,
#     marker = 5,
#     size = (800, 600)
# )

# plot solutions
# for j in axes(x,3)
#     for i in axes(x,2)
#         push!(plt, i, x[1,i,j], x[2,i,j], x[3,i,j])
#     end
# end

# save figure to file
# savefig("charged_particles.png")

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
module AbsTypes

using DataFrames

export fit!, transform!, fit_transform!
export Machine, Learner, Transformer, Workflow, Computer

abstract type Machine end
abstract type Computer <: Machine end # does computation: learner and transformer
abstract type Workflow <: Machine end # types: Linear vs Combine
abstract type Learner <: Computer end
abstract type Transformer <: Computer end

"""
    fit!(mc::Machine, input::DataFrame, output::Vector)

Generic trait to be overloaded by different subtypes of Machine.
Multiple dispatch for fit!.
"""
function fit!(mc::Machine, input::DataFrame, output::Vector)
	error(typeof(mc),"not implemented")
end

"""
    transform!(mc::Machine, input::DataFrame)

Generic trait to be overloaded by different subtypes of Machine.
Multiple dispatch for transform!.
"""
function transform!(mc::Machine, input::DataFrame)
	error(typeof(mc),"not implemented")
end

"""
    fit_transform!(mc::Machine, input::DataFrame, output::Vector)

Dynamic dispatch that calls in sequence `fit!` and `transform!` functions.
"""
function fit_transform!(mc::Machine, input::DataFrame=DataFrame(), output::Vector=Vector())
	fit!(mc,input,output)
	transform!(mc,input)
end

end

ext: jl


text:
@testset "gp" begin

    # Ensure that GP implements the AbstractGP API consistently.
    @testset "GP" begin
        rng, N, N′ = MersenneTwister(123456), 5, 6
        m, k = CustomMean(sin), Matern32Kernel()
        f = GP(m, k)
        x = collect(range(-1.0, 1.0; length=N))
        x′ = collect(range(-1.0, 1.0; length=N′))

        @test mean(f, x) == AbstractGPs._map_meanfunction(m, x)
        @test cov(f, x) == kernelmatrix(k, x)
        TestUtils.test_internal_abstractgps_interface(rng, f, x, x′)
    end

    # Check that mean-function specialisations work as expected.
    @testset "sugar" begin
        @test GP(5, Matern32Kernel()).mean isa ConstMean
        @test GP(Matern32Kernel()).mean isa ZeroMean
    end
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
oggm_config()

## Retrieving gdirs and climate for the following glaciers
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351", "RGI60-01.02170",
"RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051",
"RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316"]

gdirs = init_gdirs(rgi_ids, force=false)
tspan = (0.0,5.0) # period in years for simulation
gdirs_climate = get_gdirs_with_climate(gdirs, tspan, overwrite=false, plot=false)

# Load input data for the SIA PDE simulation
# gdirs_climate = JLD.load(joinpath(ODINN.root_dir, "test/data/gdirs_climate.jld"), "gdirs_climate")
# Load reference values for the simulation
PDE_refs = load(joinpath(ODINN.root_dir, "data/PDE_refs.jld2"))

# Run the forward PDE simulation
H_refs, V̄x_refs, V̄y_refs = @time generate_ref_dataset(gdirs_climate, tspan)

@test H_refs ≈ PDE_refs["H_refs"]
@test V̄x_refs ≈ PDE_refs["V̄x_refs"]
@test V̄y_refs ≈ PDE_refs["V̄y_refs"]


text:
# TODO (cg 2022/05/29 14:42):
#
# 5. Longer-term, this constructor occasionally redundantly computes the
#    cross-covariance for individual chunks (j,k), or at least has the potential
#    to. It would be nice to offer users the option to allocate even more, but in
#    doing so create a cache that computes all the pairwise cross-covariances once
#    and only once, and then uses that cache to fill in chunks of the potentially
#    blocked cross-covs that are used when your condix has more than a single past
#    chunk in the conditioning sets.
#
#    This will probably not be earth-shattering in terms of perf unless the
#    kernel is really expensive...but some kernels are pretty expensive. I
#    frequently end up working with expensive kernels, even. So worth doing.
#

# Just allocates all the memory for the struct. This does NOT fill in values.
function RCholesky_alloc(V::AbstractVecchiaConfig{H,D,F}, T) where{H,D,F}
  szs    = map(length, V.pts)
  diags  = map(sz->UpperTriangular(Matrix{T}(I(sz))), szs)
  odiags = map(enumerate(V.condix)) do (j,cj)
    isempty(cj) && return zeros(T, 0, 0)
    c_size = sum(k->length(V.pts[k]), cj)
    zeros(T, c_size, szs[j])
  end
  RCholesky{T}(diags, odiags, V.condix, globalidxs(V.pts), false)
end

# TODO (cg 2022/05/30 11:05): Continually look to squeeze allocations out of
# here. Maybe I can pre-allocate things for the BLAS calls, even?
function rchol_instantiate!(strbuf::RCholesky{T}, V::VecchiaConfig{H,D,F},
                            params::AbstractVector{T},
                            execmode=ThreadedEx()) where{H,D,F,T}
  @assert !strbuf.is_instantiated "This instantiation function makes extensive use of in-place algebraic operations and makes certain assumptions about the values of those buffers coming in. Please make a new struct to pass in here, or manually reset your current one."
  strbuf.is_instantiated=true
  kernel  = V.kernel
  Z       = promote_type(H,T)
  cpts_sz = V.chunksize*V.blockrank
  pts_sz  = V.chunksize
  # allocate three buffers:
  @floop execmode for j in 1:length(V.condix)
    # allocate work buffers in the thread-safe way:
    @init buf_pp  = Array{Z}(undef,  pts_sz,  pts_sz)
    @init buf_cp  = Array{Z}(undef, cpts_sz,  pts_sz)
    @init buf_cc  = Array{Z}(undef, cpts_sz, cpts_sz)
    @init buf_pts = Array{SVector{D,Float64}}(undef, cpts_sz)
    pts = V.pts[j]
    dat = V.data[j]
    cov_pp = view(buf_pp, 1:length(pts), 1:length(pts))
    if isone(j)
      # In this special case, I actually can skip the lower triangle.
      updatebuf!(cov_pp, pts, pts, kernel, params, skipltri=true)
      empty_dict = Dict{Tuple{Int64,Int64},Matrix{Z}}()
      cov_pp_f = cholesky!(Symmetric(cov_pp))
      buf      = strbuf.diagonals[1]
      ldiv!(cov_pp_f.U, buf)
      continue
    end
    # If j != 1, then I'm going to use an in-place mul! on this buffer, so I
    # need to fill it all in.
    updatebuf!(cov_pp, pts, pts, kernel, params, skipltri=false)
    # prepare conditioning points:
    idxs = V.condix[j]
    cpts = updateptsbuf!(buf_pts, V.pts, idxs)
    # prepare and fill in the matrix buffers pertaining to the cond.  points:
    cov_cp = view(buf_cp, 1:length(cpts), 1:length(pts))
    cov_cc = view(buf_cc, 1:length(cpts), 1:length(cpts))
    updatebuf!(cov_cc, cpts, cpts, kernel, params, skipltri=false)
    updatebuf!(cov_cp, cpts,  pts, kernel, params, skipltri=false)
    # pre-factorize the cpts:cpts marginal matrix, then get the two remaining
    # pieces, like the conditional covaraince of the prediction points. I
    # acknowledge that this is a little hard to read, but it really nicely cuts
    # out all the unnecessary allocations. If you do a manual check, you can
    # confirm that cov_pp becomes the conditional covariance of pts | cpts, etc.
    cov_cc_f = cholesky!(Symmetric(cov_cc))
    ldiv!(cov_cc_f.U', cov_cp)
    mul!(cov_pp, adjoint(cov_cp), cov_cp, -one(Z), one(Z))
    ldiv!(cov_cc_f.U, cov_cp)
    Djf = cholesky!(Symmetric(cov_pp))
    # Update the struct buffers. Note that the diagonal elements are actually
    # UpperTriangular, and I am not supposed to mutate those. But we do the ugly
    # hack of directly working with the data buffer backing the UpperTriangular
    # wrapper, which is probably not really recommended, but it works.
    #strbuf.odiagonals[j] .= Bt_chunks
    strbuf.odiagonals[j] .= cov_cp
    strbuf_Djf = strbuf.diagonals[j].data
    ldiv!(Djf.U, strbuf_Djf)
    # now update the block with the rmul!, being careful to now use the
    # UpperTriangular matrix, NOT the raw buffer!
    strbuf_Bt = strbuf.odiagonals[j]
    strbuf_Bt .*= -one(Z)
    rmul!(strbuf_Bt, strbuf.diagonals[j])
  end
  nothing
end

function rchol(V::VecchiaConfig{H,D,F}, params::AbstractVector{T};
               execmode=ThreadedEx()) where{H,D,F,T}
  out = RCholesky_alloc(V, T)
  rchol_instantiate!(out, V, params, execmode)
  out
end

function nll(U::RCholesky{T}, data::Matrix{Float64}) where{T}
  -logdet(U) + sum(x->x^2, U'*data)/2
end

# This is really just for debugging.
function nll_rchol(V::VecchiaConfig{H,D,F}, params::AbstractVector{T};
                   execmode=ThreadedEx()) where{H,D,F,T}
  U    = rchol(V, params; execmode=execmode)
  data = reduce(vcat, V.data)
  nll(U, data)
end

# So much manual indexing here. I'm sorry to anybody who tries to read this. I
# just really wanted to keep allocations down, so I really did some grinding on
# the logic and a lot of manual book-keeping of the indices so I could use
# setindex! instead of push!.
#
# TODO (cg 2022/05/30 15:51): This is mucho serial. Which is probably best
# because of how fast it should run for even enormous matrices. But at some
# point should consider a parallel constructor that uses the fancy tools in the
# Transducers ecosystem (@floop, append!!, push!!, etc) to see if a parallel
# constructor works faster, even though it will allocate much more.
function SparseArrays.sparse(U::RCholesky{T}) where{T}
  master_len = rchol_nnz(U)
  Iv = Vector{Int64}(undef, master_len)
  Jv = Vector{Int64}(undef, master_len)
  Vv = Vector{T}(undef, master_len)
  master_ix = 1
  # Diagonal blocks. This section makes no allocations and probably can't really
  # be improved upon.
  for (j,ix) in enumerate(U.idxs)
    Dj = U.diagonals[j]
    sz = size(Dj, 1)
    offset = ix[1]-1
    for col_ix in 1:sz
      for row_ix in 1:col_ix
        Iv[master_ix] = row_ix+offset
        Jv[master_ix] = col_ix+offset
        Vv[master_ix] = Dj[row_ix, col_ix]
        master_ix    += 1
      end
    end
  end
  # Off-diagonal blocks. There are a couple enumerates in here that will
  # inevitably make a few allocations, but I think for the moment this is about
  # as optimized as is sensible.
  for (j,ix_c) in enumerate(U.condix)
    isempty(ix_c) && continue
    Bj = U.odiagonals[j]
    c_offset = 0
    for ix_c_k in ix_c
      idxs_c = U.idxs[ix_c_k]
      for (_i1, i1) in enumerate(idxs_c)
        for (_i2, i2) in enumerate(U.idxs[j])
          Iv[master_ix] = i1
          Jv[master_ix] = i2
          Vv[master_ix] = Bj[_i1+c_offset,_i2]
          master_ix    += 1
        end
      end
      c_offset += length(idxs_c)
    end
  end
  sparse(Iv, Jv, Vv)
end


text:
using FerriteMeshParser
using Documenter

const is_ci = get(ENV, "CI", "false") == "true"

include("generate.jl")
GENERATEDEXAMPLES = [joinpath("examples", f) for f in (
    "compact_tension.md",
    "user_element.md")]

DocMeta.setdocmeta!(FerriteMeshParser, :DocTestSetup, :(using FerriteMeshParser); recursive=true)

makedocs(;
    authors="Knut Andreas Meyer and contributors",
    repo="https://github.com/Ferrite-FEM/FerriteMeshParser.jl/blob/{commit}{path}#{line}",
    sitename="FerriteMeshParser.jl",
    format=Documenter.HTML(;
        prettyurls=is_ci,
        canonical="https://Ferrite-FEM.github.io/FerriteMeshParser.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => GENERATEDEXAMPLES,
    ],
)

deploydocs(;
    repo="github.com/Ferrite-FEM/FerriteMeshParser.jl",
    devbranch="main",
    push_preview=true,
)


text:
# would be good to submit a PR to ODBC.jl...
import Base.convert
convert(::Type{DateTime}, dt::ODBC.API.SQLTimestamp) =
    DateTime(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second,
        round(dt.fraction/1e6))

"""
```
newjob(dsn; cryostat="", username="", device="",
    nmeas=1, jobstart="", jobstop="", jobstatus=0, dataserver="")
```

Create a new job in the `jobs` table. This function will return a `DataFrame`
containing the columns `job_id` and `jobsubmit` with the inserted job id
and job submission time.
"""
function newjob(dsn; cryostat="", username="", device="",
    nmeas=1, jobstart="", jobstop="", jobstatus=0, dataserver="")

    dformat(x) = x == "" ? :NULL : "'$x'"
    df = ODBC.query(dsn, """
        INSERT INTO jobs (cryostat, username, device, nmeas, jobsubmit, jobstart,
            jobstop, jobstatus, dataserver) VALUES ('$cryostat', '$username',
            '$device', $nmeas, $(dformat(now())), $(dformat(jobstart)),
            $(dformat(jobstop)), $jobstatus, '$dataserver')
            RETURNING job_id, jobsubmit;
    """)
    get(df[:job_id][1]), convert(DateTime, get(df[:jobsubmit][1]))
end

"""
```
updatejob(dsn, job_id; kwargs...)
```

Update an existing job in the `jobs` table based on its `job_id`. Specify
the fields to update with keyword arguments specified in
[`ICDataServer.newjob`](@ref).
"""
function updatejob(dsn, job_id; kwargs...)
    isempty(kwargs) && return
    pstr = reduce((a,b)->a*","*b, "$k = '$v'" for (k,v) in kwargs)
    for (k,v) in kwargs
        ODBC.execute!(dsn, "UPDATE jobs SET "*pstr*" WHERE job_id = $job_id;")
    end
end

precompile(newjob, (ODBC.DSN, ))
precompile(updatejob, (ODBC.DSN, Integer))


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the methods in `jlcodestyle/cli.jl`.
"""

# --- Imports

# Standard library
using Test

# External packages
using JuliaFormatter
using Suppressor

# Local modules
using TestTools.jltest: EnhancedTestSet
using TestTools.jlcodestyle: cli

# --- Tests

@testset EnhancedTestSet "jlcodestyle.cli.parse_args()" begin

    # --- Default arguments

    raw_args = Vector{String}()
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # --- overwrite

    # Case: raw_args = "--overwrite"
    raw_args = ["--overwrite"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => true,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-o"
    raw_args = ["-o"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => true,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # --- style

    # Case: raw_args = "--style"
    raw_args = ["--style", "yas"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => YASStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-s"
    raw_args = ["-s", "default"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => DefaultStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = invalid style
    invalid_style = "InvalidStyle"
    raw_args = ["--style", invalid_style]
    output = @capture_err begin
        args = cli.parse_args(; raw_args=raw_args)
    end
    @test startswith(output, "┌ Warning: Invalid style: $invalid_style. Using BlueStyle")
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # --- verbose

    # Case: raw_args = "--verbose"
    raw_args = ["--verbose"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => true,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-v"
    raw_args = ["-v"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => true,
        "version" => false,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # --- version

    # Case: raw_args = "--version"
    raw_args = ["--version"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => true,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-V"
    raw_args = ["-V"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => true,
        "paths" => Vector{String}(),
    )
    @test args == expected_args

    # --- paths

    raw_args = ["path/to/file-1.jl", "file-2.jl"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "overwrite" => false,
        "style" => BlueStyle(),
        "verbose" => false,
        "version" => false,
        "paths" => raw_args,
    )
    @test args == expected_args
end

@testset EnhancedTestSet "jlcodestyle.cli.run()" begin
    # --- Preparations

    # Get current directory
    cwd = pwd()

    # Construct path to test data directory
    test_file_dir = abspath(joinpath(@__DIR__, "data"))

    # --- Exercise functionality and check results

    # Case: `paths` is empty
    cd(test_file_dir)

    error = @capture_err begin
        cli.run([])
    end

    expected_error = """
                     Style errors found. Files not modified.
                     """
    @test error == expected_error

    cd(cwd)  # Restore current directory

    # Case: `paths` contains only source files without style errors
    cd(test_file_dir)

    output = @capture_out begin
        cli.run([joinpath(test_file_dir, "blue-style.jl")])
    end

    expected_output = """
                      No style errors found.
                      """
    @test output == ""

    # Case: verbose = true
    cd(test_file_dir)

    local output
    error = @capture_err begin
        output = @capture_out begin
            cli.run([joinpath(test_file_dir, "blue-style.jl")]; verbose=true)
        end
    end

    expected_error = """
                     [ Info: Style = BlueStyle
                     [ Info: Overwrite = false
                     """
    @test error == expected_error

    expected_output = """
        Formatting $(joinpath(pwd(), "blue-style.jl"))

        No style errors found.
        """
    @test output == expected_output

    # Case: `paths` contains only source files without style errors, overwrite = true
    cd(test_file_dir)

    local output
    error = @capture_err begin
        output = @capture_out begin
            cli.run(
                [joinpath(test_file_dir, "blue-style.jl")]; overwrite=true, verbose=true
            )
        end
    end

    expected_error = """
                     [ Info: Style = BlueStyle
                     [ Info: Overwrite = true
                     """
    @test error == expected_error

    expected_output = """
        Formatting $(joinpath(pwd(), "blue-style.jl"))

        No style errors found.
        """
    @test output == expected_output

    # Case: `paths` contains only source files with style errors, overwrite = true
    cd(test_file_dir)
    bluestyle_pass_file = joinpath(test_file_dir, "blue-style.jl")
    yasstyle_fail_file = joinpath(test_file_dir, "yasstyle-fail.jl")
    cp(bluestyle_pass_file, yasstyle_fail_file; force=true)

    local output
    error = @capture_err begin
        output = @capture_out begin
            cli.run(
                [yasstyle_fail_file];
                style=JuliaFormatter.YASStyle(),
                overwrite=true,
                verbose=true,
            )
        end
    end

    expected_error = """
                     [ Info: Style = YASStyle
                     [ Info: Overwrite = true

                     Style errors found. Files modified to correct errors.
                     """
    @test error == expected_error

    expected_output = """
        Formatting $(yasstyle_fail_file)
        """
    @test output == expected_output

    rm(yasstyle_fail_file)
end

@testset EnhancedTestSet "jlcodestyle.cli.run(): invalid arguments" begin

    # Case: invalid `paths` arg
    @test_throws MethodError cli.run([1, 2, 3])
end


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

module I18n

export locale

LOCALE = nothing
CALLBACKS = Function[]

function locale()
    if LOCALE === nothing
        # XXX:TBD return default locale
        return ""
    end
    LOCALE
end

function locale(s::String)
    global LOCALE = s
    # XXX:TBD call setlocale
    for cb in CALLBACKS
        cb()
    end
end

end # module

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
M=Model(:[
    W*η+R*b[-1]-b-c
    (b)-β*Expect(R[+1]*λ[+1]*(b))/λ
],:[
    b   = (0.01,100,7)
    η   = (1,0.9,0.1,1)
    ηz  := (1,0.9,0.001,1)
],:[
    b   = (0,200,b,0.9)
    c   = (0,5,0.4)
],:[
    λ 	= c^-σc
    R   = 1-δ+(α)*ηz*A[-1]^(α-1)*H^(1-α)
    W   = (1-α)*ηz*A[-1]^α*H^(-α)
    A = ∫(b,60)
    H = ∫(η,0.9)
],:[
    α   = 0.35
    δ   = 0.025
    β   = 0.99
    σc  = 2.5
    ϕh  = 1.6
    σh  = 2.0
])

@profile solve(M,400,0.5)

function solveks(M)
    solve(M,400,0.5)
    i= 0
    while i<100 && maximum(abs(M.temp.error))>5e-5
        i+=1
        solve(M,100,0.5,crit = 1e-4)
        updateA(M,0.05)
        println(round(∫(M,:A),3),"  ",round(log10(maximum(abs(M.temp.error),1)),2))
    end

    i=0
    while i<200 && maximum(abs(M.temp.error))>1e-7
        i+=1
    solve(M,1,disp = 50)
    updateA(M)
    mod(i,10)==0 && println(round(∫(M,:A),3),"  ",round(log10(maximum(abs(M.temp.error),1)),2))
    end

    println(maximum(abs(M.temp.error)))
    println("∫b : ",round(∫(M,:b),3))
end

T = 5000
nh = 50000
id = ones(Int,nh)*2
state = ones(nh)*[40 1 1]
State = zeros(T,3)
zid = 2
srand(1)
for t=-500:T
	zid = M[:ηz]([zid])[1]
	state=[M(:b,state) M[:η].x[M[:η](id)] M[:ηz].x[zid]*ones(nh)]
    t>0 && (State[t,:]=mean(state,1))
	t>0 && mod(t,25)==0 && println(State[t,1])
end


text:
## Pretty printing of Junet types ##

typestr(g::Graph) = string(
    nodecount(g), "-node ",
    edgecount(g), "-edge ",
    isdirected(g) ? "directed " : "undirected ",
    ismultigraph(g) ? "multigraph" : "graph")

memstr(m::Integer) =
    if m < 1024
        "< 1 KiB"
elseif m < 1024^2
        string(round(m / 1024, digits=1), " KiB")
elseif m < 1024^3
    string(round(m / 1024^2, digits=2), " MiB")
    else
    string(round(m / 1024^3, digits=3), " GiB")
    end

summary(d::AttributeDict) = join(
    (string(k, " (", eltype(d[k]), ")") for k = keys(d)), ", ")

function edgestr(g::Graph, count=20)
    es = edges(g)
    e = map(string, Base.Iterators.take(es, count))
    res = join(e, ", ")
    if length(es) > count
        res *= "…"
    end
    return res
end

function summary(g::Graph, long=true)
    typ = typestr(g) * " occupying " * memstr(memsize(g))
    nc = length(g.nodeattrs)
    natt = string(nc, " node attributes", nc > 0 ? (": " * summary(g.nodeattrs)) : "")
    ec = length(g.edgeattrs)
    eatt = string(ec, " edge attributes", ec > 0 ? (": " * summary(g.edgeattrs)) : "")
    println(join([typ, natt, eatt], "\n"))
    es = edgecount(g)
    if es > 0
        println("Sample edges: ", edgestr(g))
    end
end

show(io::IO, g::Graph) = print(io, typestr(g))

show(io::IO, n::Node) = print(io,
    "Node with ", ptr_length(n, Both), " adjacent edges")

function show(io::IO, e::NodePtr)
    print(io, "→ ", e.node)
    if e.id != nothing
        print(io, " (id ", e.id, ")")
    else
        print(io, " (no id)")
    end
end

show(io::IO, e::Edge) = print(io,
    e.source, e.isdir ? " → " : " ←→ ", e.target)

# TODO: make this work with EdgeIter and potential other iterators
function show(io::IO, x::PtrView)
    print(io, length(x), "-element Junet.PtrView with ", eltype(x), " values")
    if !get(io, :limit, false)
        screenheight = 20
    else
        sz = displaysize(io)
        screenheight = sz[1] - 4
    end
    screenheight < 5 || length(x) == 0 && return
    println(io)
    if length(x) < screenheight
        for i = 1:length(x)
            if i < length(x)
                println(io, ' ', x[i])
            else
                print(io, ' ', x[i])
            end
        end
    else
        halfheight = div(screenheight, 2) - 1
        for i = 1:halfheight
            println(io, ' ', x[i])
        end
        println(io, " ⋮")
        for i = length(x) - halfheight:length(x)
            if i < length(x)
                println(io, ' ', x[i])
            else
                print(io, ' ', x[i])
            end
        end
    end
end


text:
using BSplineFit, Plots

knots = [0.0,1.0,2.0,4.0,7.0]
data_x = sort(rand(100)*7)
data_y = rand(100)*0.1 .+ 0.5 .+ sin.(data_x)./2

plts = []
for order ∈ 1:3
    p = scatter(data_x, data_y, label="", title="Fit with BSpline of order $(order).")
    for id ∈ 1:length(knots)+order
        b = BSplineFit.BSpline1DBasis(order, knots, id)
        plot!(t->b(t), 0, 7, lw=2, ls=:dash, label="")
    end

    bf = fit(data_x, data_y, order, knots)
    plot!(t->bf(t), 0, 7, lc=:black, lw=2, label="", ylims=(-0.5,1.5))

    push!(plts, p)
end

p = plot(plts..., layout=grid(3,1), size=(800,800))

savefig(p, "examples/fit.svg")


text:
struct Cglsr
    B::Matrix{Float64}
    g::Vector{Float64}
    xmeans::Vector{Float64}
    ymeans::Vector{Float64}
    F::Union{Array{Float64}, Nothing}
end

"""
    cglsr(X, y; nlv, reorth = true, filt = false)
Conjugate gradient algorithm for the normal equations (CGLS) (Björck 1996).
* `X` : X-data.
* `y` : Univariate Y-data.
* `nlv` : Nb. CG iterations.
* `reorth` : If `true`, a Gram-Schmidt reorthogonalization of the normal equation
    residual vectors is done.
* `filt` : If `true`, the CG filter factors are computed (output `F`).

CGLS algorithm "7.4.1" Bjorck 1996, p.289

The code computing the re-orthogonalization (Hansen 1998) and filter factors (Vogel 1987, Hansen 1998)
is a transcription (with few adaptations) of the Matlab function `cgls`
(Saunders et al. https://web.stanford.edu/group/SOL/software/cgls/; Hansen 2008).

`X` and `y` are internally centered. The model is computed with an intercept.

## References

Björck, A., 1996. Numerical Methods for Least Squares Problems, Other Titles in Applied Mathematics.
Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611971484

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation.
Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3.
Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Manne R. Analysis of two partial-least-squares algorithms for multivariate calibration. Chemometrics Intell.
Lab. Syst. 1987; 2: 187–197.

Phatak A, De Hoog F. Exploiting the connection between
PLS, Lanczos methods and conjugate gradients: alternative proofs of some properties of PLS. J. Chemometrics
2002; 16: 361–367.

Vogel, C. R.,  "Solving ill-conditioned linear systems using the conjugate gradient method",
Report, Dept. of Mathematical Sciences, Montana State University, 1987.
"""
function cglsr(X, y; nlv, reorth = true, filt = false)
    cglsr!(copy(X), copy(y); nlv = nlv, reorth = reorth, filt = filt)
end
function cglsr!(X, y; nlv, reorth = true, filt = false)
    X = ensure_mat(X)
    y = vec(y)
    n = size(X, 1)
    p = size(X, 2)
    xmeans = colmeans(X)
    ymeans = colmeans(y)
    center!(X, xmeans)
    center!(y, ymeans)
    # Pre-allocation and initialization
    B = similar(X, p, nlv)
    b = zeros(p)
    r = y       # r = y - X * b, with b = 0
    s = X' * r
    zp = copy(s)
    q = similar(X, n)
    g = dot(s, s)
    gnew = similar(X, nlv)
    if(reorth)
        A = similar(X, p, nlv + 1)
        A[:, 1] .= s ./ sqrt(g)
    end
    F = nothing
    if(filt)
        eig = svd(X).S.^2
        if(n < p)
            eig = [eig ; zeros(p - n)]
        end
        F = similar(X, p, nlv)
        Fd = similar(X, p)
    end
    fudge = 1e-4
    # End
    @inbounds for j in 1:nlv
        mul!(q, X, zp)
        alpha = g / dot(q, q)
        b .+= alpha * zp
        r .-= alpha * q
        mul!(s, X', r)
        # Reorthogonalize s to previous s-vectors
        if reorth
            @inbounds for i in 1:j
                v = vcol(A, i)
                s .-= dot(v, s) * v
            end
            A[:, j + 1] .= s ./ sqrt(dot(s, s))
        end
        # End
        gnew[j] = dot(s, s)
        beta = gnew[j] / g
        g = copy(gnew[j])
        zp .= s .+ beta * zp
        B[:, j] .= b
        # Filter factors
        # fudge threshold is used to prevent filter factors from exploding
        if filt
            if j == 1
                F[:, 1] .= alpha * eig
                Fd .= eig .- eig .* F[:, 1] .+ beta * eig
            else
                F[:, j] .= F[:, j - 1] .+ alpha * Fd
                Fd .= eig - eig .* F[:, j] .+ beta * Fd
            end
            if j > 2
                u = (abs.(F[:, j - 1] .- 1) .< fudge) .* (abs.(F[:, j - 2] .- 1) .< fudge)
                u = findall(u)
                if length(u) > 0
                    F[u, j] .= ones(length(u))
                end
            end
        end
        # End
    end
    Cglsr(B, gnew, xmeans, ymeans, F)
end

"""
    coef(object::Cglsr)
Compute the b-coefficients of a fitted model.
* `object` : The fitted model.
* `nlv` : Nb. iterations to consider. If nothing, it is the maximum nb. iterations.
"""
function coef(object::Cglsr; nlv = nothing)
    a = size(object.B, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    B = object.B[:, nlv:nlv]
    int = object.ymeans' .- object.xmeans' * B
    (B = B, int = int)
end

"""
    predict(object::Cglsr, X; nlv = nothing)
Compute Y-predictions from a fitted model.
* `object` : The maximal fitted model.
* `X` : X-data for which predictions are computed.
* `nlv` : Nb. iterations, or collection of nb. iterations, to consider.
If nothing, it is the maximum nb. iterations.
"""
function predict(object::Cglsr, X; nlv = nothing)
    a = size(object.B, 2)
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    le_nlv = length(nlv)
    pred = list(le_nlv, Matrix{Float64})
    @inbounds for i = 1:le_nlv
        z = coef(object; nlv = nlv[i])
        pred[i] = z.int .+ X * z.B
    end
    le_nlv == 1 ? pred = pred[1] : nothing
    (pred = pred,)
end


text:
using ProgressMeter

function md(data::Data,MDinputs::MDinput)

  @unpack mass, N   = data
  @unpack dt, nsteps, total_time, iprint = MDinputs

  box            = initial_point(data);
  nc, fatm, natm = linkedlist(box,data);
  force_vec      = fvec(data)
  vel            = initial_velocity(data.N);
  times = 0.

  #implementar arquivo para salvar trajetórias

  file = open("traj.xyz","w")

  println(file,"$N")
  println(file,"Frame 1")
  for i in 1:length(box)
    println(file,"H","   ", box[i][1],"   ",box[i][2] ,"   " ,0.0)
  end

  k = 1

@showprogress 1 "Running - MD " for i in 1:nsteps
    ut = force!(box, data, force_vec, fatm, natm, nc)
    times = times + dt

   # println(i)

    for ip in 1:data.N
      # positions update @. box[ip] = box[ip] + vel[ip]*dt + force_vec[ip] * (dt*dt)/(2*mass)
      @. box[ip] = box[ip] + vel[ip]*dt + force_vec[ip] * (dt*dt)/(2*mass)
      # velocity update
      @. vel[ip] = vel[ip] + force_vec[ip] * dt/mass
    end

    if i%iprint ==0
      println("Total energy at step $i = $ut")

      k += 1
      println(file,"$N")
      println(file,"Frame $k")
      for i in 1:length(box)
        println(file,"H","   ", box[i][1],"   ",box[i][2], "   " ,0.0)
      end
    end

  end
  close(file)
end


text:
using BenchmarkTools
using OffsetArrays

const dim = 1000

x = Array{Float64}(undef, 2*dim);
y = OffsetArray{Float64}(undef, -dim + 1 : dim);
x2d = Array{Float64}(undef, 2*dim, 2*dim);
y2d = OffsetArray{Float64}(undef, -dim + 1 : dim, -dim + 1 : dim);

s = OffsetVector(1:dim, 0);
sur = 1:dim;
sior = OffsetArrays.IdOffsetRange(parent(s));

fill1d(x) = for i in axes(x,1); x[i] = i; end
fill2d(x) = for j in axes(x,2); for i in axes(x,1); x[i,j] = i + j; end; end
update(x) = for i in axes(x,1); x[i] = x[i] + i; end
update2d(x) = for j in axes(x,2); for i in axes(x,1); x[i,j] = x[i,j] + i + j; end; end
update_eachindex(x) = for i in eachindex(x); x[i] = x[i] + i; end

unsafe_fill(x) = @inbounds(for i in axes(x,1); x[i] = i; end)
unsafe_fill2d(x) = @inbounds(for j in axes(x,2); for i in axes(x,1); x[i,j] = i + j; end; end)
unsafe_update(x) = @inbounds(for i in axes(x,1); x[i] = x[i] + i; end)
unsafe_update2d(x) = @inbounds(for j in axes(x,2); for i in axes(x,1); x[i,j] = x[i,j] + i + j; end; end)
unsafe_update_eachindex(x) = @inbounds(for i in eachindex(x); x[i] = x[i] + i; end)

vectorlinearindexing(a, ax) = a[ax]
vectorCartesianindexing(a, ax1, ax2) = a[ax1, ax2]
nestedvectorlinearindexing(a, ax1, ax2) = a[ax1[ax2]]

macro showbenchmark(ex)
	print(ex, " : ")
	quote
		println(@benchmark esc($ex))
	end
end

@showbenchmark fill1d(x)
@showbenchmark fill1d(y)
@showbenchmark fill2d(x2d)
@showbenchmark fill2d(y2d)
@showbenchmark update(x)
@showbenchmark update(y)
@showbenchmark update2d(x2d)
@showbenchmark update2d(y2d)
@showbenchmark update_eachindex(x)
@showbenchmark update_eachindex(y)

@showbenchmark unsafe_fill(x)
@showbenchmark unsafe_fill(y)
@showbenchmark unsafe_fill2d(x2d)
@showbenchmark unsafe_fill2d(y2d)
@showbenchmark unsafe_update(x)
@showbenchmark unsafe_update(y)
@showbenchmark unsafe_update2d(x2d)
@showbenchmark unsafe_update2d(y2d)
@showbenchmark unsafe_update_eachindex(x)
@showbenchmark unsafe_update_eachindex(y)

# Benchmarks of vector indexing using OffsetRanges as axes
@showbenchmark vectorlinearindexing(x, s)
@showbenchmark vectorlinearindexing(x, sur)
@showbenchmark vectorlinearindexing(x, sior)
@showbenchmark vectorlinearindexing(y, s)
@showbenchmark vectorlinearindexing(y, sur)
@showbenchmark vectorlinearindexing(y, sior)

@showbenchmark vectorlinearindexing(sur, s)
@showbenchmark vectorlinearindexing(sur, sur)
@showbenchmark vectorlinearindexing(sur, sior)

@showbenchmark vectorCartesianindexing(x2d, s, s)
@showbenchmark vectorCartesianindexing(x2d, sur, sur)
@showbenchmark vectorCartesianindexing(x2d, sior, sior)

@showbenchmark nestedvectorlinearindexing(x, s, s)
@showbenchmark nestedvectorlinearindexing(x, sur, sur)
@showbenchmark nestedvectorlinearindexing(x, s, sior)
@showbenchmark nestedvectorlinearindexing(x, sur, sior)
@showbenchmark nestedvectorlinearindexing(x, sior, sior)
@showbenchmark vectorlinearindexing(x, sior[sior])
@showbenchmark nestedvectorlinearindexing(x, sur, sior)
@showbenchmark vectorlinearindexing(x, sur[sior])
@showbenchmark nestedvectorlinearindexing(x, sior, sur)
@showbenchmark vectorlinearindexing(x, sior[sur])
@showbenchmark nestedvectorlinearindexing(x, sur, sur)


text:
"SPARK tableau for Gauss-Legendre Runge-Kutta method with s stages."
function TableauSPARKGLRK(s; name=Symbol("GLRK($s)"))
    g = TableauGauss(s)
    ω = get_GLRK_ω_matrix(s)
    δ = zeros(0, s)

    return TableauSPARK(name, g.o,
                        g.a, g.a, g.a, g.a,
                        g.a, g.a, g.a, g.a,
                        g.b, g.b, g.b, g.b,
                        g.c, g.c, g.c, g.b,
                        ω, δ)
end

"SPARK tableau for Gauss-Lobatto methods."
function TableauSPARKLobatto(A, B; name=Symbol("Lobatto($s)"))
    @assert A.s == B.s

    s = A.s
    o = min(A.o, B.o)
    δ = zeros(0, s)

    return TableauSPARK(name, o,
                        A.a, B.a, A.a, B.a,
                        A.a, B.a, A.a, B.a,
                        A.b, B.b, A.b, B.b,
                        A.c, B.c, A.c, A.b,
                        get_lobatto_ω_matrix(s), δ,
                        get_lobatto_nullvector(s))
end

"SPARK tableau for Gauss-Lobatto IIIA-IIIB method with s stages."
function TableauSPARKLobattoIIIAIIIB(s)
    TableauSPARKLobatto(TableauLobattoIIIA(s), TableauLobattoIIIB(s); name=Symbol("LobattoIIIAIIIB($s)"))
end

"SPARK tableau for Gauss-Lobatto IIIB-IIIA method with s stages."
function TableauSPARKLobattoIIIBIIIA(s)
    TableauSPARKLobatto(TableauLobattoIIIB(s), TableauLobattoIIIA(s); name=Symbol("LobattoIIIBIIIA($s)"))
end

"SPARK tableau for Gauss-Legendre/Gauss-Lobatto-IIIA-IIIB methods."
function TableauSPARKGLRKLobattoIIIAIIIB(s, σ=s+1)
    g = TableauGauss(s)
    A = TableauLobattoIIIA(σ)
    B = TableauLobattoIIIB(σ)

    α = A.a[2:σ, 1:σ]
    ã = B.a[1:σ, 1:s]

    o = min(g.o, A.o, B.o)
    δ = zeros(0, σ)

    return TableauSPARK(Symbol("GLRK($s)LobattoIIIAIIIB($σ)"), o,
                        g.a, g.a, α,   α,
                        ã,   ã,   A.a, B.a,
                        g.b, g.b, A.b, B.b,
                        g.c, g.c, A.c, A.b,
                        get_lobatto_ω_matrix(σ), δ)
end

"SPARK tableau for Gauss-Legendre/Gauss-Lobatto-IIIB-IIIA methods."
function TableauSPARKGLRKLobattoIIIBIIIA(s, σ=s+1)
    g = TableauGauss(s)
    A = TableauLobattoIIIA(σ)
    B = TableauLobattoIIIB(σ)

    α = B.a[2:σ, 1:σ]
    ã = A.a[1:σ, 1:s]

    o = min(g.o, B.o, A.o)
    δ = zeros(0, σ)

    return TableauSPARK(Symbol("GLRK($s)LobattoIIIAIIIB($σ)"), o,
                        g.a, g.a, α,   α,
                        ã,   ã,   B.a, A.a,
                        g.b, g.b, B.b, A.b,
                        g.c, g.c, A.c, A.b,
                        get_lobatto_ω_matrix(σ), δ)
end

function TableauSPARKLobatto(name, l1::Tableau{T}, l2::Tableau{T},
                                   l3::Tableau{T}, l4::Tableau{T}=l3, d=[]; R∞=1) where {T}

    @assert l1.s == l2.s == l3.s == l4.s

    ω = get_lobatto_ω_matrix(l3.s)
    δ = zeros(T, 0, 1)

    if length(d) == 0
        return TableauSPARK(name, min(l1.o, l2.o, l3.o, l4.o),
                            l3.a, l1.a, l4.a, l2.a,
                            l3.a, l1.a, l4.a, l2.a,
                            l3.b, l1.b, l4.b, l2.b,
                            l3.c, l1.c, l4.c, l4.b,
                            ω, δ)
    else
        @assert length(d) == q.s == p.s

        return TableauSPARK(name, min(l1.o, l2.o, l3.o, l4.o),
                            l3.a, l1.a, l4.a, l2.a,
                            l3.a, l1.a, l4.a, l2.a,
                            l3.b, l1.b, l4.b, l2.b,
                            l3.c, l1.c, l4.c, l4.b,
                            ω, δ, d)
    end
end

"Tableau for Gauss-Lobatto IIIA-IIIB-IIIC method with s stages."
function TableauSPARKLobABC(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobc = TableauLobattoIIIC(s)
    TableauSPARKLobatto(Symbol("LobattoIIIABC($s)"), loba, lobc, lobb; R∞=(-1)^(s+1))
end

"Tableau for Gauss-Lobatto IIIA-IIIB-IIID method with s stages."
function TableauSPARKLobABD(s)
    loba = TableauLobattoIIIA(s)
    lobb = TableauLobattoIIIB(s)
    lobd = TableauLobattoIIID(s)
    TableauSPARKLobatto(Symbol("LobattoIIIABD($s)"), loba, lobd, lobb; R∞=(-1)^(s+1))
end

"SPARK Tableau for Variational Partitioned Runge-Kutta Methods."
function TableauSPARKVPRK(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    ω = zeros(T, q.s, q.s+1)
    δ = zeros(T, 0, 1)

    for i in 1:q.s
        ω[i,i] = 1
    end

    if length(d) == 0
        return TableauSPARK(name, min(q.o, p.o),
                            q.a, p.a, q.a, p.a,
                            q.a, p.a, q.a, p.a,
                            q.b, p.b, q.b, p.b,
                            q.c, p.c, q.c, q.b,
                            ω, δ)
    else
        @assert length(d) == q.s == p.s

        return TableauSPARK(name, min(q.o, p.o),
                            q.a, p.a, q.a, p.a,
                            q.a, p.a, q.a, p.a,
                            q.b, p.b, q.b, p.b,
                            q.c, p.c, q.c, q.b,
                            ω, δ, d)
    end
end

"Tableau for Variational Gauss-Legendre method with s stages."
function TableauSPARKGLVPRK(s)
    glrk = TableauGauss(s)
    TableauSPARKVPRK(Symbol("GLVPRK($s)"), glrk, glrk; R∞=(-1)^s)
end


text:
module TimeSeriesA1P2

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTIO, CTReg, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]

#does common tasks on a Yahoo dataframe
function processYahooDF!(df::DataFrame)::DataFrame

  #read and pre-process the data
  dateFormat::String = YAHOO_DATE_FORMAT
  if :Date ∈ names(df)
    rename!(df, :Date, :DATE)
  elseif :_Date ∈ names(df)
    rename!(df, :_Date, :DATE)
  end

  df[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(df[:,:DATE])::DataVector{Date}
  sort!(df, cols = [:DATE])
  df[:,:returns] = similar(df[:,:Adj_Close])
  df[2:end,:returns] = log.(df[2:end,:Adj_Close] ./  df[1:(end-1),:Adj_Close])

  return df
end

#pre-process the s&p500 data
function preProcessSP500()::Void

  #read and pre-process the data
  sp500DF::DataFrame = readtable("$DATA_PATH\\$SP500_NAME.csv")
  sp500DF = processYahooDF!(sp500DF)

  #write to JLS binary file
  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls", "w")
  serialize(stream, sp500DF)
  close(stream)

  return nothing
end

function preProcessDBV()::Void

  #read and pre-process the data
  dbvDF::DataFrame = readtable("$DATA_PATH\\$DBV_NAME.csv")
  dbvDF = processYahooDF!(dbvDF)

  #write to JLS binary file
  stream::IOStream = open("$DATA_PATH\\$DBV_NAME.jls", "w")
  serialize(stream, dbvDF)
  close(stream)

  return nothing
end

function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

#plots the descriptive graphs
function visualizeA1P2(df::DataFrame)

  #convenience containers
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  outDF::DataFrame = meltdf(df, :DATE) #put data in long form
  completecases!(outDF)

  push!(plotNames, "Histogram")
  push!(plots, plot(outDF, x=:value, color=:variable,
      Guide.title("Histograms"),
      Guide.xlabel("Log Return"),Guide.ylabel("Density"),
      Coord.Cartesian(xmin=-0.12, xmax=0.12),
      Geom.histogram(bincount=20, density=true, position=:dodge)))

  push!(plotNames, "Time Series")
  push!(plots, vstack(
        plot(outDF[outDF[:,:variable].==:sp500,:], y=:value, x=:DATE, Geom.line,
          Guide.xlabel("Year"),Guide.ylabel("Log Return"), Guide.title("SP500")),
        plot(outDF[outDF[:,:variable].==:dbv,:], y=:value, x=:DATE, Geom.line,
          Guide.xlabel("Year"),Guide.ylabel("Log Return"), Guide.title("DBV"))))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
  end

end

#calculates the test statistic and T value of a data vector
skewnessStat(data::Vector{Float64})::Float64 = skewness(data)/√(6./length(data))
kurtosisStat(data::Vector{Float64})::Float64 = kurtosis(data)/√(24./length(data))
JBStat(data::Vector{Float64})::Float64 =
  skewness(data)^2/(6./length(data)) + kurtosis(data)^2/(24./length(data))

mutable struct UnivariateTest
  statistic::Function
  distribution::Distribution
  α::Float64
  H₀::Float64
  testType::Symbol
end

#convenience constructor for univariate test
UnivariateTest(statistic::Function, distribution::Distribution, α::Float64;
  H₀::Float64=0.0, testType::Symbol=:upper)::UnivariateTest =
  UnivariateTest(statistic, distribution, α, H₀, testType)

formNum(x::T where T<:Real)::String = format(NUM_FORMAT, x)

function testUni(U::UnivariateTest, data::Vector{Float64}; verbose::Bool = true)::NTuple{2,Float64}
  θ::Float64 = U.statistic(data .- U.H₀)

  p::Float64 = cdf(U.distribution, θ)

  #adjust the p value for the specific test type
  if U.testType == :upper
    p = 1.0 - p
  elseif  U.testType == :lower
  elseif U.testType == :twoTailed
    p = p > 0.5 ? (1-p)*2 : 2*p
  else
    warn("Test type not found in univariate test")
  end

  return (θ, p)
end

#script for doing the hypothesis test
function hypothesisA1P2(df::DataFrame)
  α::Float64 = 0.05

  #label the table
  colNames::Vector{Vector{String}} =
    [["SP500", "DBV"],
    ["statistic", "p-value", "reject \$H_0\$?", "statistic", "p-value", "reject \$H_0\$?"]]
  numCols = length(colNames[end])

  rowNames::Vector{String} = ["Skewness", "Kurtosis", "Normality"]
  numRows::Int = length(rowNames)

  #form the data
  sp500::Vector{Float64} = dropna(df[:,:dbv])
  dbv::Vector{Float64} = dropna(df[:,:dbv])

  NSP500::Int = length(sp500)
  NDBV::Int = length(dbv)

  #set the test parameters
  sp500Tests::Vector{UnivariateTest} =
    [UnivariateTest(skewnessStat, TDist(NSP500-1), α, testType=:twoTailed),
    UnivariateTest(kurtosisStat, TDist(NSP500-1), α, testType=:twoTailed),
    UnivariateTest(JBStat, TDist(NSP500), α, testType=:upper)]

  dbvTests::Vector{UnivariateTest} =
    [UnivariateTest(skewnessStat, TDist(NDBV-1), α, testType=:twoTailed),
    UnivariateTest(kurtosisStat, TDist(NDBV-1), α, testType=:twoTailed),
    UnivariateTest(JBStat, Chisq(2), α, testType=:upper)]

  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #run the tests and record the results
  for r::Int ∈ 1:numRows
    sp500θ::Float64, sp500P::Float64 = testUni(sp500Tests[r], sp500)
    content[r][1] = formNum(sp500θ)
    content[r][2] = formNum(sp500P)
    content[r][3] = sp500P < α ? "X" : ""

    dbvθ::Float64, dbvP::Float64 = testUni(dbvTests[r], dbv)
    content[r][4] = formNum(dbvθ)
    content[r][5] = formNum(dbvP)
    content[r][6] = dbvP < α ? "X" : ""
  end

  #write the table
  hypothesisTable::String = texTable( "Summary of ST and LT Interest Yields",
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
    widthColNames = [[3, 3], ones(Int,6)],
    alignmentColNames = [["c", "c"], ["r" for i::Int ∈ 1:length(6)]]
  )

  return hypothesisTable
end

#script for doing the regressions
function regressionA1P2(df::DataFrame, SE::Function; caption::String="TBD", title::String="TBD")

  #model specification
  XSpec::Symbol = parse("sp500")
  XName::Vector{Symbol} = [:intercept; :sp500]
  YSpec::Symbol = :dbv
  colNames::Vector{String} = ["DBV (1)"]

  #run the regression
  reg::Vector{CTLM} = [CTLM(df, XSpec, YSpec, XNames=XName, YName = :dbv)]

  #get some descriptive stats
  descRowNames::Vector{String} = ["N", "\$R^{2}\$"]
  descContent::Vector{Vector{String}} = [["$(reg[1].N)"], ["$(round(getR(reg[1])^2.0,2))"]]

  TableText::String = texTable(reg, SE, [:sp500; :intercept],
    titleCaption = title,
    #colNames = [["($i)" for i∈1:length(models)]],
    colNames = [colNames],
    contentRowNames = ["SP500", "Intercept"],
    descRowNames = descRowNames,
    descContent = descContent,
    stars=true,
    decimalDigits = 4,
    caption=caption,
    starStrings = OVERRIDE_STAR_STRINGS)

end

function A1P2Script(; refreshData::Bool = false)
  sp500DF::DataFrame = getDF("$DATA_PATH\\$SP500_NAME.jls", refreshData, preProcessSP500)
  dbvDF::DataFrame = getDF("$DATA_PATH\\$DBV_NAME.jls", refreshData, preProcessDBV)

  #merge the dataframes
  rename!(sp500DF, :returns, :sp500)
  rename!(dbvDF, :returns, :dbv)
  df::DataFrame = join(sp500DF[:,[:DATE, :sp500]],dbvDF[:,[:DATE, :dbv]], on=[:DATE], kind=:inner)

  #make the graphs
  visualizeA1P2(df)

  #do the hypothesis tast
  hypothesisTable::String = hypothesisA1P2(df)
  homoskedasticTable::String = regressionA1P2(df, getHomoskedΣ!, title="DBV vs S\\\&P500",
    caption="Homoskedastic Standard Errors")
  modWhiteTable::String = regressionA1P2(df, getModWhiteΣ!, title="DBV vs S\\\&P500",
    caption="Modified White Standard Errors")
  whiteTable::String = regressionA1P2(df, getWhiteΣ!, title="DBV vs S\\\&P500",
    caption="White Standard Errors")
  #=whiteTableSlow::String = regressionA1P2(df, getWhiteΣSlow, title="DBV vs S\\\&P500",
    caption="White Standard Errors (Test)")=#

  #record the tables
  writeTables2File([hypothesisTable, homoskedasticTable, modWhiteTable, whiteTable],
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "OutputTables.tex")

end

@time begin
  A1P2Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using MeshCat, GeometryBasics, Colors

vis = Visualizer()
render(vis) # open(vis)

# circle
θ = range(0, stop = 2 * π, length = 100)
posy = 0.2 * cos.(θ)
posz = 0.8 .+ 0.2 .* sin.(θ)

ee_positions = hcat([[0.0; posy[t]; posz[t]] for t = 1:100]...)
points = collect(eachcol(ee_positions))
material = LineBasicMaterial(color=colorant"orange", linewidth=10.0)
setobject!(vis["ee path"], Object(PointCloud(points), material, "Line"))


text:
using SafeTestsets

@safetestset "Test my_f" begin include("test_my_f.jl") end


text:
"""
Internal structure that remembers that `quantity` was requested by `request`.
"""
struct ResourceRequest{Q<:Real}
    quantity::Q
    request::Function
end

"""
Resource type for holding numeric values (like amount of liquid).
It stores current `quantity` of matter and its allowed `lo` and `hi` amounts.
Servers can get matter from the resource with optional maximum number of
requests pending for fulfillment.

Fields:
* `quantity`       current quantity in resource
* `lo`            minimum quantity of resource
* `hi`            maximum quantity of resource
* `fifo_requests` if `true` `requests` is FIFO, otherwise LIFO
* `max_requests`  maximum `requests` size
* `requests`      vector of request and requested quantity

Functions in `requests` must accept one argument `Scheduler`, so they should
know the amount they requested. When resource arrives to a queue there is a try
to immediately dispatch it to pending requests. When new request arrives there
is a try to immediately fulfill it.

Initially an empty `SimResource` with no requests is constructed.
Initial `quantity`, `lo` and `hi` may be provided. By default `SimResource` is
empty, and has minimum quantity of zero and unbounded maximum.
"""
mutable struct SimResource{Q<:Real} <: AbstractReservoir
    quantity::Q
    lo::Q
    hi::Q
    fifo_requests::Bool
    max_requests::Int
    requests::Vector{ResourceRequest{Q}}

    function SimResource{Q}(;quantity::Q=zero(Q), lo::Q=zero(Q), hi::Q=typemax(Q),
                      fifo_requests::Bool=true,
                      max_requests::Int=typemax(Int)) where Q
        lo <= quantity <= hi || error("wrong quantity/lo/hi combination")
        max_requests > 0 || error("max_requests must be positive")
        new(quantity, lo, hi, fifo_requests, max_requests,
            Vector{ResourceRequest{Q}}())
    end
end

function dispatch!(s::Scheduler, r::SimResource)
    # while is needed as large provide! can fulfill many requests
    while (!isempty(r.requests) &&
           r.lo <= r.quantity - r.requests[end].quantity <= r.hi)
        req = pop!(r.requests)
        r.quantity -= req.quantity
        register!(s, req.request)
    end
end

function request!(s::Scheduler, r::SimResource{Q}, quantity::Q,
                  request::Function) where Q
    rr = ResourceRequest{Q}(quantity, request)
    length(r.requests) < r.max_requests || return false, rr
    qend = r.fifo_requests ? pushfirst! : push!
    qend(r.requests, rr)
    dispatch!(s, r)
    return true, rr
end

function waive!(r::SimResource{Q}, res_request::ResourceRequest{Q}) where Q
    idx = findfirst(isequal(res_request), r.requests)
    isa(idx, Nothing) && return false
    deleteat!(r.requests, idx)
    return true
end

function provide!(s::Scheduler, r::SimResource{Q}, quantity::Q) where Q
    startq = r.quantity
    # we do not want to oveflow the container, but allow large provisions
    r.quantity = clamp(startq+quantity, r.lo, r.hi)
    added = r.quantity - startq
    dispatch!(s, r)
    return added # return how much was added
end


text:
abstract type TimeKind end
struct Discrete <: TimeKind end
struct Continuous <: TimeKind end

timekind(::T) where T = TimeKind(T)

TimeKind(::Type{<: AbstractContinuationCache{PC}}) where PC = TimeKind(PC)
TimeKind(::Type{<: AbstractProblemCache{P}}) where P = TimeKind(P)


text:
"""
    lyapunov_exponents(phase_dynamics!, u0, tspan; <keyword arguments>)

Calculate Lyapunov exponents of a dynamical system.
"""
function lyapunov_exponents(phase_dynamics!,
                            u0,
                            tspan,
                            t_attr::Integer;
                            discrete=false,
                            progress=-1,
                            kwargs...)
    if discrete
        prob = ContinuousLEProblem(phase_dynamics!, u0, tspan, kwargs...)
    else
        prob = DiscreteLEProblem(phase_dynamics!, u0, tspan, kwargs...)
    end
    solver = init(prob; progress=progress)
    solve!(solver, t_attr; progress=progress)
    lyapunov_exponents(solver)
end

de_prob(stage::Union{PhaseRelaxer, AbstractRenormalizer}) =
    de_prob(stage.integrator)
de_prob(integrator::DiscreteIterator) = integrator.prob
de_prob(integrator::ODEIntegrator) = integrator.sol.prob


text:
"""
    struct UniformCartesianMesh{N,T} <: AbstractMesh{N,T}

An `N`-dimensional cartesian grid given as the tensor-product of `N`
one-dimensional `LinRange{T}` grids.

Iterating over a `UniformCartesianMesh` generates the elements which compose the mesh;
i.e. the `HyperRectangle` cells.
"""
struct UniformCartesianMesh{N,T} <: AbstractMesh{N,T}
    grids::NTuple{N,LinRange{T}}
end

grids(g::UniformCartesianMesh)     = g.grids
grids(g::UniformCartesianMesh,dim) = g.grids[dim]

Base.keys(m::UniformCartesianMesh{N,T}) where {N,T} = (HyperRectangle{N,T},)

Base.step(m::UniformCartesianMesh) = step.(grids(m))

"""
    UniformCartesianMesh(domain::HyperRectangle,sz::NTuple)
    UniformCartesianMesh(domain::HyperRectangle;step::NTuple)

Construct a uniform `UniformCartesianMesh` with `sz[d]` elements along dimension
`d`. If the kwarg `step` is passed, construct a `UniformCartesianMesh` with
elements of approximate size `step`.
"""
function UniformCartesianMesh(domain::HyperRectangle{N,T},sz::NTuple{N}) where {N,T}
    lc = low_corner(domain)
    uc = high_corner(domain)
    grids1d = ntuple(N) do n
        # to have sz elements, need sz+1 points in LinRange.
        npts = sz[n] + 1
        LinRange{float(T)}(lc[n], uc[n], npts) # use float(T) in case T<:Integer
    end
    UniformCartesianMesh(grids1d)
end
UniformCartesianMesh(domain::HyperRectangle{N,T},sz::Int) where {N,T} = UniformCartesianMesh(domain,ntuple(i->sz,N))

# in case you pass arguments like UniformCartesianMesh(xgrid,ygrid), convert
# them to a Tuple
UniformCartesianMesh(grids::Vararg{LinRange{T}}) where {T} = UniformCartesianMesh(Tuple(grids))

function UniformCartesianMesh(domain::HyperRectangle{N};step::NTuple{N}) where {N}
    lc = low_corner(domain)
    hc = high_corner(domain)
    sz = ntuple(N) do i
        (hc[i] - lc[i]) / step[i] |> ceil |> Int
    end
    UniformCartesianMesh(domain,sz)
end

ambient_dimension(g::UniformCartesianMesh{N}) where {N} = N

xgrid(g::UniformCartesianMesh) = g.grids[1]

ygrid(g::UniformCartesianMesh) = g.grids[2]

zgrid(g::UniformCartesianMesh) = g.grids[3]

# implement ElementIterator interface to UniformCartesianMesh

function Base.size(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh})
    g = mesh(iter)
    length.(g.grids) .- 1
end

Base.length(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh}) = prod(size(iter))

function Base.CartesianIndices(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh})
    CartesianIndices(size(iter))
end

function Base.getindex(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh}, I::CartesianIndex)
    m = mesh(iter)
    N = ambient_dimension(m)
    @assert N == length(I)
    lc = ntuple(N) do dim
        i = I[dim]
        m.grids[dim][i]
    end
    hc = ntuple(N) do dim
        i = I[dim] + 1
        m.grids[dim][i]
    end
    HyperRectangle(lc, hc)
end
function Base.getindex(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh},I...)
    iter[CartesianIndex(I)]
end
function Base.getindex(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh}, i::Int)
    I = CartesianIndices(iter)
    iter[I[i]]
end

function Base.iterate(iter::ElementIterator{<:HyperRectangle,<:UniformCartesianMesh},state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

# since UniformCartesianMesh has only one elment type, calling ElementIterator without
# specifying the has clear sense
function ElementIterator(m::UniformCartesianMesh)
    E = keys(m) |> first
    ElementIterator(m,E)
end

# NodeIterator over UniformCartesianMesh

function Base.IteratorSize(iter::Type{NodeIterator{UniformCartesianMesh{N,T}}}) where {N,T}
    Base.HasShape{N}()
end

function Base.size(iter::NodeIterator{<:UniformCartesianMesh})
    g = mesh(iter)
    length.(g.grids)
end

Base.length(iter::NodeIterator{<:UniformCartesianMesh}) = prod(size(iter))

function Base.CartesianIndices(iter::NodeIterator{<:UniformCartesianMesh})
    CartesianIndices(size(iter))
end

function Base.getindex(iter::NodeIterator{<:UniformCartesianMesh}, I::CartesianIndex)
    m = mesh(iter)
    N = ambient_dimension(m)
    @assert N == length(I)
    svector(N) do dim
        i = I[dim]
        m.grids[dim][i]
    end
end
function Base.getindex(iter::NodeIterator{<:UniformCartesianMesh},I...)
    iter[CartesianIndex(I)]
end
function Base.getindex(iter::NodeIterator{<:UniformCartesianMesh}, i::Int)
    I = CartesianIndices(iter)
    iter[I[i]]
end

function Base.iterate(iter::NodeIterator{<:UniformCartesianMesh},state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
import Core: Expr

# Terminals
function julia_normalization_map(c::Int32, x::Ptr{Nothing})::Int32
    return c == 0x00B5 ? 0x03BC : # micro sign -> greek small letter mu
           c == 0x025B ? 0x03B5 : # latin small letter open e -> greek small letter
           c
end

# Note: This code should be in julia base
function utf8proc_map_custom(str::String, options)
    norm_func = @cfunction julia_normalization_map Int32 (Int32, Ptr{Nothing})
    nwords = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Nothing}, Ptr{Nothing}),
                   str, sizeof(str), C_NULL, 0, options, norm_func, C_NULL)
    nwords < 0 && Base.Unicode.utf8proc_error(nwords)
    buffer = Base.StringVector(nwords * 4)
    nwords = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Nothing}, Ptr{Nothing}),
                   str, sizeof(str), buffer, nwords, options, norm_func, C_NULL)
    nwords < 0 && Base.Unicode.utf8proc_error(nwords)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && Base.Unicode.utf8proc_error(nbytes)
    return String(resize!(buffer, nbytes))
end

function normalize_julia_identifier(str::AbstractString)
    options = Base.Unicode.UTF8PROC_STABLE | Base.Unicode.UTF8PROC_COMPOSE
    utf8proc_map_custom(String(str), options)
end

function sized_uint_literal(s::AbstractString, b::Integer)
    # We know integers are all ASCII, so we can use sizeof to compute
    # the length of ths string more quickly
    l = (sizeof(s) - 2) * b
    l <= 8   && return Base.parse(UInt8,   s)
    l <= 16  && return Base.parse(UInt16,  s)
    l <= 32  && return Base.parse(UInt32,  s)
    l <= 64  && return Base.parse(UInt64,  s)
    # l <= 128 && return Base.parse(UInt128, s)
    l <= 128 && return Expr(:macrocall, GlobalRef(Core, Symbol("@uint128_str")), nothing, s)
    return Base.parse(BigInt, s)
end

function sized_uint_oct_literal(s::AbstractString)
    s[3] == 0 && return sized_uint_literal(s, 3)
    len = sizeof(s)
    (len < 5  || (len == 5  && s <= "0o377")) && return Base.parse(UInt8, s)
    (len < 8  || (len == 8  && s <= "0o177777")) && return Base.parse(UInt16, s)
    (len < 13 || (len == 13 && s <= "0o37777777777")) && return Base.parse(UInt32, s)
    (len < 24 || (len == 24 && s <= "0o1777777777777777777777")) && return Base.parse(UInt64, s)
    # (len < 45 || (len == 45 && s <= "0o3777777777777777777777777777777777777777777")) && return Base.parse(UInt128, s)
    # return Base.parse(BigInt, s)
    (len < 45 || (len == 45 && s <= "0o3777777777777777777777777777777777777777777")) && return Expr(:macrocall, GlobalRef(Core, Symbol("@uint128_str")), nothing, s)
    return Meta.parse(s)
end

function _literal_expr(x)
    if headof(x) === :TRUE
        return true
    elseif headof(x) === :FALSE
        return false
    elseif is_nothing(x)
        return nothing
    elseif headof(x) === :INTEGER || headof(x) === :BININT || headof(x) === :HEXINT || headof(x) === :OCTINT
        return Expr_int(x)
    elseif isfloat(x)
        return Expr_float(x)
    elseif ischar(x)
        return Expr_char(x)
    elseif headof(x) === :MACRO
        return Symbol(valof(x))
    elseif headof(x) === :STRING || headof(x) === :TRIPLESTRING
        return valof(x)
    elseif headof(x) === :CMD
        return Expr_cmd(x)
    elseif headof(x) === :TRIPLECMD
        return Expr_tcmd(x)
    end
end

const TYPEMAX_INT64_STR = string(typemax(Int))
const TYPEMAX_INT128_STR = string(typemax(Int128))
function Expr_int(x)
    is_hex = is_oct = is_bin = false
    val = replace(valof(x), "_" => "")
    if sizeof(val) > 2 && val[1] == '0'
        c = val[2]
        c == 'x' && (is_hex = true)
        c == 'o' && (is_oct = true)
        c == 'b' && (is_bin = true)
    end
    is_hex && return sized_uint_literal(val, 4)
    is_oct && return sized_uint_oct_literal(val)
    is_bin && return sized_uint_literal(val, 1)
    # sizeof(val) <= sizeof(TYPEMAX_INT64_STR) && return Base.parse(Int64, val)
    return Meta.parse(val)
    # # val < TYPEMAX_INT64_STR && return Base.parse(Int64, val)
    # sizeof(val) <= sizeof(TYPEMAX_INTval < TYPEMAX_INT128_STR128_STR) && return Base.parse(Int128, val)
    # # val < TYPEMAX_INT128_STR && return Base.parse(Int128, val)
    # Base.parse(BigInt, val)
end

function Expr_float(x)
    if !startswith(valof(x), "0x") && 'f' in valof(x)
        return Base.parse(Float32, replace(replace(valof(x), 'f' => 'e'), '_' => ""))
    end
    Base.parse(Float64, replace(valof(x), "_" => ""))
end
function Expr_char(x)
    val = _unescape_string(valof(x)[2:prevind(valof(x), lastindex(valof(x)))])
    # one byte e.g. '\xff' maybe not valid UTF-8
    # but we want to use the raw value as a codepoint in this case
    sizeof(val) == 1 && return Char(codeunit(val, 1))
    length(val) == 1 || error("Invalid character literal: $(Vector{UInt8}(valof(x)))")
    val[1]
end

# Expressions

function Expr(x::EXPR)
    if isidentifier(x)
        if headof(x) === :NONSTDIDENTIFIER
            if startswith(valof(x.args[1]), "@")
                Symbol("@", normalize_julia_identifier(valof(x.args[2])))
            else
                Symbol(normalize_julia_identifier(valof(x.args[2])))
            end
        else
            return Symbol(normalize_julia_identifier(valof(x)))
        end
    elseif iskeyword(x)
        if headof(x) === :BREAK
            return Expr(:break)
        elseif headof(x) === :CONTINUE
            return Expr(:continue)
        else
            return Symbol(lowercase(string(headof(x))))
        end
    elseif isoperator(x)
        return Symbol(valof(x))
    elseif ispunctuation(x)
        if headof(x) === :DOT
            return :(.)
        else
            # We only reach this if we have a malformed expression.
            Expr(:error)
        end
    elseif isliteral(x)
        return _literal_expr(x)
    elseif isbracketed(x)
        return Expr(x.args[1])
    elseif x.head isa EXPR
        Expr(Expr(x.head), Expr.(x.args)...)
    elseif x.head === :quotenode
        QuoteNode(Expr(x.args[1]))
    elseif x.head === :globalrefdoc
        GlobalRef(Core, Symbol("@doc"))
    elseif x.head === :globalrefcmd
        if VERSION > v"1.1"
            GlobalRef(Core, Symbol("@cmd"))
        else
            Symbol("@cmd")
        end
    elseif x.head === :macrocall && is_getfield_w_quotenode(x.args[1]) && !ismacroname(x.args[1].args[2].args[1])
        # Shift '@' to the right
        valofrhs = valof(x.args[1].args[2].args[1])
        valofrhs = valofrhs === nothing ? "" : valofrhs
        new_name = Expr(:., remove_at(x.args[1].args[1]), QuoteNode(Symbol("@", valofrhs)))
        Expr(:macrocall, new_name, Expr.(x.args[2:end])...)
    elseif x.head === :macrocall && isidentifier(x.args[1]) && valof(x.args[1]) == "@."
        Expr(:macrocall, Symbol("@__dot__"), Expr.(x.args[2:end])...)
    elseif x.head === :macrocall && length(x.args) == 3 && x.args[1].head === :globalrefcmd && x.args[3].head == :string
        Expr(:macrocall, Expr(x.args[1]), Expr(x.args[2]), join(valof.(x.args[3])))
    elseif x.head === :string && length(x.args) > 0 && (x.args[1].head === :STRING || x.args[1].head === :TRIPLESTRING) && isempty(valof(x.args[1]))
        # Special conversion needed - the initial text section is treated as empty for the represented string following lowest-common-prefix adjustments, but exists in the source.
        Expr(:string, Expr.(x.args[2:end])...)
    elseif x.args === nothing
        Expr(Symbol(lowercase(String(x.head))))
    elseif x.head === :errortoken
        Expr(:error)
    else
        Expr(Symbol(lowercase(String(x.head))), Expr.(x.args)...)
    end
end

function remove_at(x)
    if isidentifier(x) && valof(x) !== nothing && first(valof(x)) == '@'
        return Symbol(valof(x)[2:end])
    elseif is_getfield_w_quotenode(x)
        Expr(:., remove_at(x.args[1]), QuoteNode(remove_at(x.args[2].args[1])))
    else
        Expr(x)
    end
end

# cross compatability for line number insertion in macrocalls
if VERSION > v"1.1-"
    Expr_cmd(x) = Expr(:macrocall, GlobalRef(Core, Symbol("@cmd")), nothing, valof(x))
    Expr_tcmd(x) = Expr(:macrocall, GlobalRef(Core, Symbol("@cmd")), nothing, valof(x))
else
    Expr_cmd(x) = Expr(:macrocall, Symbol("@cmd"), nothing, valof(x))
    Expr_tcmd(x) = Expr(:macrocall, Symbol("@cmd"), nothing, valof(x))
end

function clear_at!(x)
    if x isa Expr && x.head == :.
        if x.args[2] isa QuoteNode && string(x.args[2].value)[1] == '@'
            x.args[2].value = Symbol(string(x.args[2].value)[2:end])
        end
        if x.args[1] isa Symbol && string(x.args[1])[1] == '@'
            x.args[1] = Symbol(string(x.args[1])[2:end])
        else
            clear_at!(x.args[1])
        end
    end
end

"""
    remlineinfo!(x)
Removes line info expressions. (i.e. Expr(:line, 1))
"""
function remlineinfo!(x)
    if isa(x, Expr)
        if x.head == :macrocall && x.args[2] !== nothing
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (@isdefined(LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
            insert!(x.args, 2, nothing)
        else
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (@isdefined(LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
        end
        if x.head == :elseif && x.args[1] isa Expr && x.args[1].head == :block && length(x.args[1].args) == 1
            x.args[1] = x.args[1].args[1]
        end
    end
    x
end


text:
using DelimitedFiles

function savepagedata(model::Model, component::Symbol, variable::Symbol, filepath::AbstractString)
    # Handle relative paths
    if filepath[1] ∉ ['.', '/'] && !isfile(filepath)
        filepath = joinpath(@__DIR__, "..", "..", filepath)
    end

    dims = variable_dimensions(model, component, variable)
    open(filepath, "w") do io
        write(io, "# Index: " * join(dims, ", ") * "\n\n")

        if dims == [:region]
            write(io, join(dim_keys(model, :region), ",") * "\n")
            writedlm(io, model[component, variable], ',')
        elseif dims == [:time]
            write(io, join(dim_keys(model, :time), ",") * "\n")
            writedlm(io, model[component, variable], ',')
        elseif dims == [:time, :region]
            writedlm(io, [["year" dim_keys(model, :region)...]; [["$x" for x in dim_keys(model, :time)] model[component, variable]]], ',')
        else
            error("Unknown dimensions in variable $component.$variable.")
        end
    end
end


text:
# with @eval, we test that include=:static works even for a non-toplevel
# @testset, what matters is that a testset is defined after including the file
@eval @testset "include_static_included1 $i" for i=1:2
    @test true
    @testset "nested include_static_included1" begin
        @test true
        push!(Hijack.RUN, 2)
        # test that `include` kwarg is forwarded
        include("include_static_included2.jl")
    end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Transforms

using StaticArrays: SVector, SMatrix
using LinearAlgebra: cross, normalize

# export Vec3,
#     unitX3,
#     unitY3,
#     unitZ3,
#     Vec4,
#     unitX4,
#     unitY4,
#     unitZ4,
#     unitW4,
#     get_orthogonal_vectors,
#     Transform,
#     identitytransform,
#     right,
#     up,
#     forward,
#     # it was imported from OpticSim, not required?? IDK
#     #   origin TODO origin(t::Transform{<:Real}) = Vec3(t[1, 4], t[2, 4], t[3, 4])
#     rotationX,
#     rotationY,
#     rotationZ,
#     translation,
#     rotation,
#     rotationd,
#     translation,
#     scale,
#     local2world,
#     world2local,
#     decomposeRTS,
#     rotmatd,
#     rotmat

"""
`Vec3{T}` provides an immutable vector of fixed length 3 and type `T`.

`Vec3` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3)` or `Vec3([1.0, 2.0, 3.0])`.
It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec3(rand(3))`.
"""
Vec3{T} = SVector{3, T}

# empty constructor - initialized with zeros
Vec3(::Type{T} = Float64) where {T <: Real} = zeros(Vec3{T})

"""
returns the unit vector `[1, 0, 0]`
"""
function unitX3(::Type{T} = Float64) where {T <: Real}
    return Vec3{T}(one(T), zero(T), zero(T))
end
"""
returns the unit vector `[0, 1, 0]`
"""
function unitY3(::Type{T} = Float64) where {T <: Real}
    return Vec3{T}(zero(T), one(T), zero(T))
end
"""
returns the unit vector `[0, 0, 1]`
"""
function unitZ3(::Type{T} = Float64) where {T <: Real}
    return Vec3{T}(zero(T), zero(T), one(T))
end

"""
`Vec4{T}` provides an immutable vector of fixed length 4 and type `T`.

`Vec4` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3, 4)` or `Vec3([1.0, 2.0, 3.0, 4.0])`.
It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec4(rand(4))`.
"""
Vec4{T} = SVector{4, T}

# empty constructor - initialized with zeros
Vec4(::Type{T} = Float64) where {T <: Real} = zeros(Vec4{T})

"""
    Vec4(v::SVector{3, T}) where {T<:Real} -> Vec4{T}

Accept `SVector` and create a `Vec4` type [v[1], v[2], v[3], 1]
"""
function Vec4(v::SVector{3, T}) where {T <: Real}
    return Vec4{T}(v[1], v[2], v[3], one(T))
end

"""
    Vec4(m::SMatrix{3,N,T} where{N,T<:Real} -> SMatrix{3,N,T})
Input is matrix of 3d points, each column is one point. Returns matrix of 3d points with 1 appended in the last row.
"""
function Vec4(m::SMatrix{3, N, T}) where {N, T <: Real}
    return vcat(m, ones(SMatrix{1, N, T}))::SMatrix{4, N, T}
end
"""
returns the unit vector `[1, 0, 0, 0]`
"""
function unitX4(::Type{T} = Float64) where {T <: Real}
    return Vec4{T}(one(T), zero(T), zero(T), zero(T))
end
"""
returns the unit vector `[0, 1, 0, 0]`
"""
function unitY4(::Type{T} = Float64) where {T <: Real}
    return Vec4{T}(zero(T), one(T), zero(T), zero(T))
end
"""
returns the unit vector `[0, 0, 1, 0]`
"""
function unitZ4(::Type{T} = Float64) where {T <: Real}
    return Vec4{T}(zero(T), zero(T), one(T), zero(T))
end
"""
returns the unit vector `[0, 0, 0, 1]`
"""
function unitW4(::Type{T} = Float64) where {T <: Real}
    return Vec4{T}(zero(T), zero(T), zero(T), one(T))
end

# utility function that given a vector, return 2 orthogonal vectors to that one
function get_orthogonal_vectors(direction::Vec3{T}) where {T <: Real}
    axis1 = normalize(direction)

    dp = dot(unitX3(T), axis1)
    # check the special case where the input vector is parallel to the X axis
    if (dp >= one(T) - eps(T))
        # axis1 = unitX(T);
        axis2 = unitY3(T)
        axis3 = unitZ3(T)
        return (axis2, axis3)
    elseif (dp <= -one(T) + eps(T))
        # axis1 = -unitX(T);
        axis2 = -unitY3(T)
        axis3 = -unitZ3(T)
        return (axis2, axis3)
    end

    axis3 = normalize(cross(axis1, unitX3(T)))
    axis2 = normalize(cross(axis3, axis1))
    return (axis2, axis3)
end

#---------------------------------------
# 3D Transform / Local Frame
#---------------------------------------
"""
    Transform{S<:Real}

Transform encapsulating rotation, translation and scale in 3D space. Translation happens **after** rotation.

```julia
Transform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T)
Transform(rotation::SMatrix{3,3,S}, translation::SVector{3,S})
Transform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1})
```
`θ`, `ϕ` and `ψ` in first constructor are in **radians**.
"""
struct Transform{T}
    matrix::SMatrix{4, 4, T, 16}

    """ This is a private internal function. In general don't want to allow users to populate Transform matrices with arbitrary elements. Not to be called by code outside of the Transform module. Don't use Transform{Float64}(...) for example. Instead use Transform(..)"""
    function Transform{T}(
        a11::T,
        a21::T,
        a31::T,
        a41::T,
        a12::T,
        a22::T,
        a32::T,
        a42::T,
        a13::T,
        a23::T,
        a33::T,
        a43::T,
        a14::T,
        a24::T,
        a34::T,
        a44::T,
    ) where {T <: Real}
        return new{T}(
            SMatrix{4, 4, T, 16}(
                a11,
                a21,
                a31,
                a41,
                a12,
                a22,
                a32,
                a42,
                a13,
                a23,
                a33,
                a43,
                a14,
                a24,
                a34,
                a44,
            ),
        )
    end

    Transform{T}(mat::SMatrix{4, 4, T, 16}) where {T <: Real} = new{T}(mat)
end

#functions to make Transform compatible with base matrix API
matrix(a::Transform) = a.matrix

# Base.length(a::Transform) = length(matrix(a))
function Base.getindex(a::Transform, indices::Vararg{Int, N}) where {N}
    return getindex(matrix(a), indices...)
end
# Base.iterate(a::Transform) = iterate(matrix(a))
# Base.iterate(a::Transform{Float64}, b::Tuple{StaticArrays.SOneTo{16}, Int64}) = iterate(matrix(a),b)

Base.collect(a::Transform) = collect(matrix(a))

function Base.:*(transa::Transform{T}, transb::Transform{T}) where {T <: Real}
    return Transform{T}(matrix(transa) * matrix(transb))
end

function Base.:*(
    transform::Transform{T},
    v::SVector{3, S},
) where {T <: Real, S <: Number}
    t = matrix(transform)

    res = t * Vec4(v)
    if (t[4, 4] == one(T))
        return SVector{3, S}(res[1], res[2], res[3])
    else
        return SVector{3, S}(res[1] / res[4], res[2] / res[4], res[3] / res[4])
    end
end

#Transform is not necessarily constrained to be a rigid body transformation so use general invers.
Base.inv(a::Transform{T}) where {T <: Real} = Transform{T}(inv(matrix(a)))

""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
```
id = identitytransform()
m = fill(1mm,3,4)
id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
```
"""
function Base.:*(
    transform::Transform{T},
    m::SMatrix{3, N, S},
) where {N, T <: Real, S <: Number}
    res = MMatrix{3, N, T}(undef)
    t = matrix(transform)

    for outcol in 1:N
        for row in 1:3
            sum = T(0)
            for incol in 1:3
                sum += t[row, incol] * m[incol, outcol]
            end
            #implicit 1 w coordinate value
            sum += t[row, 4]
            res[row, outcol] = sum
        end
        if t[4, 4] != 1
            res[:, outcol] /= t[4, 4]
        end
    end
    return SMatrix{3, N, S}(res)
end

""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
WARNING: this doesn't work. The translation component of the transform matrix has to be in Unitful units but the rotation part has to be in unitless units for this to work. Only works if one assumes that the translation part of the transform implicitly has the same units as the Unitful vectors being transformed. Brittle and likely to cause obscure bugs.
```
id = identitytransform()
m = fill(1mm,3,4)
id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
```
"""
function Base.:*(
    transform::Transform{T},
    m::SMatrix{4, N, S},
) where {N, T <: Real, S <: Number}
    return matrix(transform) * m
end

Base.transpose(a::Transform{T}) where {T <: Real} = Transform{T}(a.matrix')

# END of functions for compatibility with base matrix API

# for compatability ith the "old" RigidBodyTransform
"""
identitytransform([S::Type]) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform.
"""
function identitytransform(::Type{T} = Float64) where {T <: Real}
    return Transform{T}(
        one(T),
        zero(T),
        zero(T),
        zero(T),
        zero(T),
        one(T),
        zero(T),
        zero(T),
        zero(T),
        zero(T),
        one(T),
        zero(T),
        zero(T),
        zero(T),
        zero(T),
        one(T),
    )
end

"""
    Transform([S::Type]) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform.
"""
function Transform(::Type{T} = Float64) where {T <: Real}
    return identitytransform(T)
end

"""
    Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real}

Costruct a transform from the input columns.
"""
function Transform(
    colx::Vec3{T},
    coly::Vec3{T},
    colz::Vec3{T},
    colw::Vec3{T} = zero(Vec3{T}),
) where {T <: Real}
    return Transform{T}(
        vcat(
            hcat(colx, coly, colz, colw),
            SMatrix{1, 4, T}(zero(T), zero(T), zero(T), one(T)),
        ),
    )
end

"""
    Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real}

Costruct a transform from the input columns.
"""
function Transform(
    colx::Vec4{T},
    coly::Vec4{T},
    colz::Vec4{T},
    colw::Vec4{T},
) where {T <: Real}
    return Transform{T}(hcat(colx, coly, colz, colw))
end

"""
    Transform(origin, forward) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the local frame with origin and forward direction. the other 2 axes are computed automaticlly.
"""
function Transform(
    origin::Vec3{T},
    forward::Vec3{T} = unitZ3(),
) where {T <: Real}
    forward = normalize(forward)
    right, up = get_orthogonal_vectors(forward)
    return Transform(right, up, forward, origin)
end

function Transform(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) where {T <: Number}
    return Transform(rotmat(T, θ, ϕ, ψ), Vec3{T}(x, y, z))
end

"""
    Transform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real} -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix and translation vector.
"""
function Transform(
    rotation::SMatrix{3, 3, T},
    translation::SVector{3, T},
) where {T <: Real}
    return Transform{T}(
        rotation[1, 1],
        rotation[2, 1],
        rotation[3, 1],
        zero(T),
        rotation[1, 2],
        rotation[2, 2],
        rotation[3, 2],
        zero(T),
        rotation[1, 3],
        rotation[2, 3],
        rotation[3, 3],
        zero(T),
        translation[1],
        translation[2],
        translation[3],
        one(T),
    )
end

"""
    Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real} -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix (3x3) and translation vector of length 3.
"""
function Transform(
    rotation::AbstractArray{T, 2},
    translation::AbstractArray{T, 1},
) where {T <: Real}
    @assert size(rotation)[1] == size(rotation)[2] == length(translation) == 3
    return Transform{T}(
        rotation[1, 1],
        rotation[2, 1],
        rotation[3, 1],
        zero(T),
        rotation[1, 2],
        rotation[2, 2],
        rotation[3, 2],
        zero(T),
        rotation[1, 3],
        rotation[2, 3],
        rotation[3, 3],
        zero(T),
        translation[1],
        translation[2],
        translation[3],
        one(T),
    )
end

# define some utility functions

"""
    right(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the first column, representing the "X" axis.
"""
right(t::Transform{<:Real}) = normalize(Vec3(t[1, 1], t[2, 1], t[3, 1]))

"""
    up(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the second column, representing the "Y" axis.
"""
up(t::Transform{<:Real}) = normalize(Vec3(t[1, 2], t[2, 2], t[3, 2]))

"""
    forward(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the third column, representing the "Z" axis.
"""
forward(t::Transform{<:Real}) = normalize(Vec3(t[1, 3], t[2, 3], t[3, 3]))

# """
#     origin(t::Transform{<:Real}) -> Vec3
#
# Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the fourth column, containing the translation part of the transform in 3D.
# """
# OpticSim.origin(t::Transform{<:Real}) = Vec3(t[1, 4], t[2, 4], t[3, 4])

"""
    rotationX(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the x-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationX(angle::T) where {T <: Real}
    c = cos(angle)
    s = sin(angle)

    row1 = unitX4()
    row2 = Vec4(zero(T), c, s, zero(T))
    row3 = Vec4(zero(T), -s, c, zero(T))
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end

"""
    rotationY(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the y-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationY(angle::T) where {T <: Real}
    c = cos(angle)
    s = sin(angle)

    row1 = Vec4(c, zero(T), -s, zero(T))
    row2 = unitY4()
    row3 = Vec4(s, zero(T), c, zero(T))
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end

"""
    rotationZ(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the z-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationZ(angle::T) where {T <: Real}
    c = cos(angle)
    s = sin(angle)

    row1 = Vec4(c, s, zero(T), zero(T))
    row2 = Vec4(-s, c, zero(T), zero(T))
    row3 = unitZ4()
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end

"""
    rotation(t::Transform{T}) where {T<:Real} -> SMatrix{3,3,T}

returns the rotation part of the transform `t` - a 3x3 matrix.
"""
function rotation(t::Transform{T}) where {T <: Real}
    rot = SMatrix{3, 3, T}(
        t[1, 1],
        t[2, 1],
        t[3, 1],
        t[1, 2],
        t[2, 2],
        t[3, 2],
        t[1, 3],
        t[2, 3],
        t[3, 3],
    )
    return Transform(rot, zero(Vec3{T}))
end

"""
    rotate(a::Transform{T}, vector::Union{Vec3{T}, SVector{3,T}}) where {T<:Real} -> Vec3{T}

apply the rotation part of the transform `a` to the vector `vector` - this operation is usually used to rotate direction vectors.
"""
function rotate(
    a::Transform{T},
    vector::Union{Vec3{T}, SVector{3, T}},
) where {T <: Real}
    return rotation(a) * vector
end

"""
    rotation([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**.
"""
rotation(θ::T, ϕ::T, ψ::T) where {T <: Number} = rotation(Float64, θ, ϕ, ψ)
function rotation(::Type{S}, θ::T, ϕ::T, ψ::T) where {T <: Number, S <: Real}
    return Transform(rotmat(S, θ, ϕ, ψ), zeros(SVector{3, S}))
end
"""
    rotationd([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**.
"""
rotationd(θ::T, ϕ::T, ψ::T) where {T <: Number} = rotationd(Float64, θ, ϕ, ψ)
function rotationd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T <: Number, S <: Real}
    return Transform(rotmatd(S, θ, ϕ, ψ), zeros(SVector{3, S}))
end
export translation, rotation, rotationd
"""
    translation(x::T, y::T, z::T) where {T<:Real}

Creates a translation transform
"""
function translation(::Type{S}, x::T, y::T, z::T) where {T <: Number, S <: Real}
    return convert(Transform{S}, translation(x, y, z))
end

function translation(x::T, y::T, z::T) where {T <: Real}
    col1 = unitX4(T)
    col2 = unitY4(T)
    col3 = unitZ4(T)
    col4 = Vec4(x, y, z, one(T))

    return Transform(col1, col2, col3, col4)
end

"""
    translation(x::T, y::T, z::T) where {T<:Real}

Creates a translation transform
"""
function translation(t::Vec3{T}) where {T <: Real}
    return translation(t[1], t[2], t[3])
end

"""
    scale(x::T, y::T, z::T) where {T<:Real}

Creates a scaling transform
"""
function scale(x::T, y::T, z::T) where {T <: Real}
    col1 = unitX4(T) * x
    col2 = unitY4(T) * y
    col3 = unitZ4(T) * z
    col4 = unitW4(T)

    return Transform(col1, col2, col3, col4)
end

"""
    scale(s::T) where {T<:Real}

Creates a uniform scaling transform
"""
function scale(s::T) where {T <: Real}
    return scale(s, s, s)
end

"""
    scale(t::Vec3{T}) where {T<:Real}

Creates a scaling transform
"""
function scale(t::Vec3{T}) where {T <: Real}
    return scale(t[1], [2], [3])
end

"""
    local2world(t::Transform{T}) where {T<:Real}

return the transform matrix that takes a point in the local coordinate system to the global one
"""
function local2world(t::Transform{T}) where {T <: Real}
    return t
end

"""
    world2local(t::Transform{T}) where {T<:Real}

return the transform matrix that takes a point in the global coordinate system to the local one
"""
function world2local(t::Transform{T}) where {T <: Real}
    return inv(t)
end

"""
    decomposeRTS(tr::Transform{T}) where {T<:Real}

return a touple containing the rotation matrix, the translation vector and the scale vecto represnting the transform.
"""
function decomposeRTS(tr::Transform{T}) where {T <: Real}
    t = Vec3(tr[1, 4], tr[2, 4], tr[3, 4])
    sx = norm(Vec3(tr[1, 1], tr[2, 1], tr[3, 1]))
    sy = norm(Vec3(tr[1, 2], tr[2, 2], tr[3, 2]))
    sz = norm(Vec3(tr[1, 3], tr[2, 3], tr[3, 3]))
    s = Vec3(sx, sy, sz)
    rot = SMatrix{4, 4, T}(
        tr[1, 1] / sx,
        tr[2, 1] / sx,
        tr[3, 1] / sx,
        0,
        tr[1, 2] / sy,
        tr[2, 2] / sy,
        tr[3, 2] / sy,
        0,
        tr[1, 3] / sz,
        tr[2, 3] / sz,
        tr[3, 3] / sz,
        0,
        0,
        0,
        0,
        1,
    )

    return rot, t, s
end

"""
    rotmatbetween([S::Type], a::SVector{3,T}, b::SVector{3,T}) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation between vetors `a` and `b`, i.e. rotation(a,b) * a = b.
"""
function rotmatbetween(a::SVector{3, T}, b::SVector{3, T}) where {T <: Real}
    return rotmatbetween(Float64, a, b)
end
function rotmatbetween(
    ::Type{S},
    a::SVector{3, T},
    b::SVector{3, T},
) where {T <: Real, S <: Real}
    # TODO: Brian, is there a hidden assumption that a and b are normalized?
    v = cross(a, b)
    c = dot(a, b)
    V = SMatrix{3, 3, S, 9}(0, v[3], -v[2], -v[3], 0, v[1], v[2], -v[1], 0)
    R = I + V + V^2 * one(T) / (one(T) + c)
    return SMatrix{3, 3, S, 9}(R)
end

"""
    rotmatd([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**.
"""
function rotmatd(θ::T, ϕ::T, ψ::T) where {T <: Number}
    return rotmat(Float64, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ))
end
function rotmatd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T <: Number, S <: Real}
    return rotmat(S, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ))
end
"""
    rotmat([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**.
"""
rotmat(θ::T, ϕ::T, ψ::T) where {T <: Number} = rotmat(Float64, θ, ϕ, ψ)
function rotmat(::Type{S}, θ::T, ϕ::T, ψ::T) where {T <: Number, S <: Real}
    sinψ = sin(ψ)
    sinϕ = sin(ϕ)
    sinθ = sin(θ)
    cosψ = cos(ψ)
    cosϕ = cos(ϕ)
    cosθ = cos(θ)
    return SMatrix{3, 3, S, 9}(
        cosψ * cosϕ,
        sinψ * cosϕ,
        -sinϕ,
        cosψ * sinϕ * sinθ - sinψ * cosθ,
        sinψ * sinϕ * sinθ + cosψ * cosθ,
        cosϕ * sinθ,
        cosψ * sinϕ * cosθ + sinψ * sinθ,
        sinψ * sinϕ * cosθ - cosψ * sinθ,
        cosϕ * cosθ,
    )
end

end # module