text:
export Plane

"""
    Plane(𝗻̂::Vec{Dim,T}, d::T)
    Plane(P₁::Point{3}, P₂::Point{3}, P₃::Point{3})

A `Plane` in 3-dimensional space that satisfies X ⋅𝗻̂ = d, where
X is a `Dim`-dimensional point and 𝗻̂ is the unit normal to the plane.
"""
struct Plane{T}
    𝗻̂::Vec{3,T}
    d::T
end

function Plane(P₁::Point{3}, P₂::Point{3}, P₃::Point{3})
    𝗻̂ = normalize((P₂ - P₁) × (P₃ - P₁))
    return Plane(𝗻̂, coordinates(P₁) ⋅ 𝗻̂)
end

ext: jl

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
function print_result(io::IO, params::Array, obs::Array)
    if isa(params[1]["Model"], Union{Type{Ising}, Type{XY}})
        names = ["L", "T"]
    elseif isa(params[1]["Model"], Union{Type{Potts}, Type{Clock}})
        names = ["Q", "L", "T"]
    elseif isa(params[1]["Model"], Union{Type{QuantumXXZ}})
        names = ["S2", "L", "T", "Jz", "Jxy", "Gamma"]
    end
    print_result(io, params, obs, names)
end

function print_result(io::IO, params, obs::Array, param_names::Array)
    i  = 1
    for name in param_names
        println(io, "# $i : $name")
        i += 1
    end
    ks = collect(keys(obs[1]))
    sort!(ks)
    for name in ks
        println(io, "# $i, $(i+1) : $name")
        i += 2
    end

    for (p, o) in zip(params, obs)
        for name in param_names
            print(io, p[name], " ")
        end

        for name in ks
            print(io, mean(o[name]), " ")
            print(io, stderror(o[name]), " ")
        end
        println(io)
    end
end


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug, @sync, @context!

using GPUArrays

using LinearAlgebra

using BFloat16s

using CEnum

using Memoization

using DataStructures

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")
include("libcublas_deprecated.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

function math_mode!(handle, mode)
    flags = 0

    # https://github.com/facebookresearch/faiss/issues/1385
    if version() > v"11"
        flags = CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION
    end

    flags |= if mode == CUDA.PEDANTIC_MATH
        # prevent use of tensor cores
        if version() < v"11"
            CUBLAS_DEFAULT_MATH
        else
            CUBLAS_PEDANTIC_MATH
        end
    elseif mode == CUDA.DEFAULT_MATH
        # use tensor cores, but don't reduce precision
        if version() < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_DEFAULT_MATH
        end
    elseif mode == CUDA.FAST_MATH
        # we'll additionally select a compute-mode with reduced precision whenever possible
        if version() < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_TF32_TENSOR_OP_MATH
        end
    end

    cublasSetMathMode(handle, flags)

    return
end

# cache for created, but unused handles
const idle_handles = HandleCache{CuContext,cublasHandle_t}()
const idle_xt_handles = HandleCache{Any,cublasXtHandle_t}()

function handle()
    state = CUDA.active_state()
    handle, stream, math_mode = get!(task_local_storage(), (:CUBLAS, state.context)) do
        new_handle = pop!(idle_handles, state.context) do
            cublasCreate()
        end

        finalizer(current_task()) do task
            push!(idle_handles, state.context, new_handle) do
                @context! skip_destroyed=true state.context cublasDestroy_v2(handle)
            end
        end

        cublasSetStream_v2(new_handle, state.stream)

        math_mode!(new_handle, state.math_mode)

        new_handle, state.stream, state.math_mode
    end::Tuple{cublasHandle_t,CuStream,CUDA.MathMode}

    if stream != state.stream
        cublasSetStream_v2(handle, state.stream)
        task_local_storage((:CUBLAS, state.context), (handle, state.stream, math_mode))
        stream = state.stream
    end

    if math_mode != state.math_mode
        math_mode!(handle, state.math_mode)
        task_local_storage((:CUBLAS, state.context), (handle, stream, state.math_mode))
        math_mode = state.math_mode
    end

    return handle
end

function xt_handle()
    ctxs = Tuple(context(dev) for dev in devices())
    get!(task_local_storage(), (:CUBLASxt, ctxs)) do
        handle = pop!(idle_xt_handles, ctxs) do
            cublasXtCreate()
        end

        finalizer(current_task()) do task
            push!(idle_xt_handles, ctxs, handle) do
                # TODO: which context do we need to destroy this on?
                cublasXtDestroy(handle)
            end
        end

        # select all devices
        devs = convert.(Cint, devices())
        cublasXtDeviceSelect(handle, length(devs), devs)

        handle
    end::cublasXtHandle_t
end

## logging

const MAX_LOG_BUFLEN = UInt(1024*1024)
const log_buffer = Vector{UInt8}(undef, MAX_LOG_BUFLEN)
const log_cursor = Threads.Atomic{UInt}(0)
const log_cond = Ref{Base.AsyncCondition}()    # root

function log_message(ptr)
    # NOTE: this function may be called from unmanaged threads (by cublasXt),
    #       so we can't even allocate, let alone perform I/O.
    len = @ccall strlen(ptr::Cstring)::Csize_t
    cursor = Threads.atomic_add!(log_cursor, len+1)
    if cursor+len+1 <= MAX_LOG_BUFLEN
        @ccall memmove((pointer(log_buffer)+cursor)::Ptr{Nothing},
                       pointer(ptr)::Ptr{Nothing}, (len+1)::Csize_t)::Nothing
        @ccall uv_async_send(log_cond[]::Ptr{Nothing})::Cint
    end

    return
end

function _log_message(blob)
    # the message format isn't documented, but it looks like a message starts with a capital
    # and the severity (e.g. `I!`), and subsequent lines start with a lowercase mark (`!i`)
    #
    # lines are separated by a \0 if they came in separately, but there may also be multiple
    # actual lines separated by \n in each message.
    for message in split(blob, r"[\0\n]+(?=[A-Z]!)")
        code = message[1]
        lines = split(message[3:end], r"[\0\n]+[a-z]!")
        submessage = join(lines, '\n')
        if code == 'I'
            @debug submessage
        elseif code == 'W'
            @warn submessage
        elseif code == 'E'
            @error submessage
        elseif code == 'F'
            error(submessage)
        else
            @info "Unknown log message, please file an issue.\n$message"
        end
    end
    return
end

function __runtime_init__()
    # register a log callback
    log_cond[] = Base.AsyncCondition() do async_cond
        blob = ""
        while true
            blob = unsafe_string(pointer(log_buffer), log_cursor[])
            if Threads.atomic_cas!(log_cursor, UInt(length(blob)), UInt(0)) == length(blob)
                break
            end
        end
        _log_message(blob)
        return
    end
    if !Sys.iswindows() # NVIDIA bug #3321130
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end


text:
module Publish

export publish, publish_tex, setup, showit

import PGFPlots: Plots, Axis, save, TEX, PDF, SVG, Plottable, plot
export Plots, Axis, save

# Global variables
#  _output_buffer	Contains tex string
#  _num_plots		Number of plots in this document
#  _plot_list		List of the plot names
#  _publish_file	Name of file being pulished
#  _folder_name		Name of directory in which file to publish is housed
_output_buffer = nothing
_num_plots = 0
_plot_list = nothing
_publish_file = nothing
_folder_name = nothing

# Creates the tex file
# Called by publish
# filename includes the foldername...
function publish_tex(filename::String; runcode::Bool=true)
	global _output_buffer = IOBuffer()
	global _publish_file = filename
	global _plot_list = String[]

	# TODO: Here I should check for larger names (test.rar.jl)
	filename_arr = split(filename, ".")
	if filename_arr[end] == "jl"
		filename_short = filename_arr[1]
	end

	tex = open("$(filename_short).tex", "w")
	println(tex, "\\documentclass{article}")
	println(tex, "\\usepackage{caption}")

	# For quotation marks in lstlisting
	println(tex, "\\newcommand{\\dbllq}{``}")
	println(tex, "\\newcommand{\\dblrq}{''}")

	# Standalone stuff for figures
	println(tex, "\\usepackage{standalone}")
	println(tex, "\\usepackage[usenames,dvipsnames]{xcolor}")
	println(tex, "\\usepackage{pgfplots}")
	println(tex, "\\pgfplotsset{compat=newest}")

	# Change the Margins
	println(tex, "\\usepackage[margin=1.0in]{geometry}")

	# Headers, footers
	println(tex, "\\usepackage{fancyhdr}")
	println(tex, "\\pagestyle{fancy}")
	filename_sanitized = replace(filename, "_", "\\_")
	println(tex, "\\lhead{$filename_sanitized}")
	println(tex, "\\rhead{}")
	println(tex, "\\renewcommand{\\headrulewidth}{0.4pt}")
	println(tex, "\\renewcommand{\\footrulewidth}{0.4pt}")

	# listings stuff
	lst = readall(joinpath(Pkg.dir("Publish"), "src", "julia_listings.tex"))
	println(tex, lst)

	println(tex, "\\begin{document}")

	# Begin the code section
	println(tex, "\\begin{lstlisting}[style=JuliaStyle,numbers=left]")
	source_file = open(filename, "r")
	while (currentline = readline(source_file)) != ""
		print(tex, currentline)
	end
	close(source_file)
	println(tex, "\\end{lstlisting}")

	# Run the test file and print the output
	if runcode
		include(filename)
		output_string = takebuf_string(_output_buffer)
		if length(output_string) > 0
			println(tex, "\\begin{lstlisting}[style=Output]")
			println(tex, output_string)
			println(tex, "\\end{lstlisting}")
		end

		# We also have a bunch of tex plots
		global _num_plots
		global _plot_list
		for i = 1:_num_plots
			println(tex, "\\begin{figure}[!ht]")
			println(tex, "\\centering")
			println(tex, "\\input{$(_plot_list[i])}")
			println(tex, "\\end{figure}")
		end
		_num_plots = 0
	end
	_publish_file = nothing

	# End the document and close
	println(tex, "\\end{document}")
	close(tex)
end

# Does the same thing as normal println function,
#  but also prints to a buffer if it exists.
function Base.println(xs...)
	global _output_buffer
	if _output_buffer != nothing
		println(_output_buffer, xs...)
	end
	println(STDOUT, xs...)
end

# Does the same thing as normal print function,
#  but also prints to a buffer if it exists.
function Base.print(xs...)
	global _output_buffer
	if _output_buffer != nothing
		print(_output_buffer, xs...)
	end
	print(STDOUT, xs...)
end

function publish(filename::String; runcode::Bool=true)

	if !isfile(filename)
		error(filename, " does not exist.")
	end

	# Generate the .tex file and make pass along any possible errors
	IX = rsearch(filename,'/')              # Find the last slash
	foldername = filename[1:IX-1]           # Everything before that is the folder name
	foldername = (IX == 0 ? "." : foldername) # If there was no slash, folder is "."
	global _folder_name = foldername
	#save(TEX(f.filename * ".tex"), tp)        # Save the tex file in the directory that was given

	# This also sets the global _output_buffer to a value (as in, not nothing)
	# When it is done, set the global _output_buffer back to nothing
	publish_tex(filename; runcode=runcode)
	global _output_buffer = nothing

	# TODO: Here I should check for larger names (test.rar.jl)
	filename_arr = split(filename, ".")
	if filename_arr[end] == "jl"
		filename_short = filename_arr[1]
	end

	# From the .tex file, generate a pdf within the specified folder
	latexCommand = ``
	# TODO: I'm not sure what enableWrite18 does. I asssume it is some latex feature
	#if tp.enableWrite18
	#	latexCommand = `lualatex --enable-write18 --output-directory=$(foldername) $(f.filename)`
	#else
	#	latexCommand = `lualatex --output-directory=$(foldername) $(f.filename)`
	#end
	# Change to the folder...
	home = pwd()
	cd(foldername)
	latexCommand = `lualatex $(split(filename_short,"/")[end])`
	isPdfThere = success(latexCommand)
	cd(home)

	if !isPdfThere
		line1 = "ERROR: The pdf generation failed.\n"
		line2 = "   Be sure your latex libraries are fully up to date!\n"
		line3 = "   You tried: $(latexCommand)"
		error(line1 * line2 * line3)  # Throw an error
	end

	try
		# Shouldn't need to be try-catched anymore, but best to be safe
		# This failing is NOT critical either, so just make it a warning
		# Before, it checked for tikzDeleteIntermediate
		#if tikzDeleteIntermediate()
		#rm("$(filename).tex")
		rm("$(filename_short).aux")
		rm("$(filename_short).log")
	catch
		println("WARNING! Your intermediate files are not being deleted.")
	end
end

# Mostly a copy from the PGFPlots version
# However, here we also deal with some gobal stuff
function save(filename::String, o::Plottable)

	# Make sure you are referring to the global variables
	global _output_buffer
	global _num_plots
	global _plot_list
	global _publish_file
	global _folder_name

	# Determine the base and extension
    base_name, ext = splitext(filename)
    ext = lowercase(ext)

	# Additional functionality for Publish.jl
	# If this name matches the name they gave, do not accept
	# Check before the dot...
	# TODO: Ensure a clean break on error out
	if _publish_file != nothing
		file_base = splitext(_publish_file)[1]
		filename = "$(_folder_name)/$filename"
		if file_base == base_name
			_publish_file = nothing
			error("Plot name cannot equal filename")
		end
	end

	# Continue to the original save functionality
    if ext == ".pdf"
        save(PDF(filename), plot(o))
    elseif ext == ".svg"
        save(SVG(filename), plot(o))
    elseif ext == ".tex"
        save(TEX(filename), plot(o))
    elseif ext == "." || ext == ""
        error("You must specify a file extension.")
    else
        error("Unsupported file extensions: $ext")
    end

	# Additional functionality for Publish.jl
	# Save name and tex file so we can plot this later
	if _output_buffer != nothing
		_num_plots += 1
		temp_name = "$(base_name).tex"
		push!(_plot_list, temp_name)
		save(TEX("$(_folder_name)/$temp_name"), plot(o))
	end
end

end # module


text:
module GFit

using Printf, PrettyTables
using Statistics, Distributions
using DataStructures
using LsqFit
using MacroTools
using Dates
using ProgressMeter

import Base.show
import Base.ndims
import Base.size
import Base.length
import Base.haskey
import Base.keys
import Base.getindex
import Base.setindex!
import Base.reshape
import Base.propertynames
import Base.getproperty
import Base.setproperty!
import Base.iterate
import Base.push!

export Domain, CartesianDomain, coords, axis, roi, int_tabulated, Measures,
    Model, @expr, SumReducer, select_reducer!, domain,
    MultiModel, @patch!, @try_patch!, evaluate!, isfixed, thaw, freeze, fit!

include("domain.jl")
include("HashVector.jl")

# ====================================================================
mutable struct ExprFunction
    expr::Expr
    funct::Function
    args::Vector{Symbol}
end

macro exprfunc(_expr)
    @assert isexpr(longdef(_expr), :function)
    expr = prettify(_expr)
    args = convert(Vector{Symbol}, splitdef(expr)[:args])
    return esc(:(GFit.ExprFunction($(QuoteNode(expr)), $expr, $args)))
end

# ====================================================================
# Parameter and component identifiers
struct ParamID
    name::Symbol
    index::Int
    ParamID(pname::Symbol) = new(pname, 0)                  # scalar param
    ParamID(pname::Symbol, index::Int) = new(pname, index)  # vector of params
end

struct CompParamID
    cname::Symbol
    param::ParamID
end

mutable struct Parameter
    val::Float64
    low::Float64              # lower limit value
    high::Float64             # upper limit value
    step::Float64
    fixed::Bool
    patched::Float64
    unc::Float64
    Parameter(value::Number) = new(float(value), -Inf, +Inf, NaN, false, NaN, NaN)
end

# ====================================================================
# Components:
#
# A *component* is a generic implementation of a building block for a
# model. It must inherit `AbstractComponent` and implement the
# `evaluate!` method (optionally also `prepare!`).  The structure may
# contain zero or more field of type Parameter (see above)
abstract type AbstractComponent end

function getparams(comp::AbstractComponent)
    params = OrderedDict{ParamID, Parameter}()
    for pname in fieldnames(typeof(comp))
        par = getfield(comp, pname)
        if isa(par, Parameter)
            params[ParamID(pname)] = par
        elseif isa(par, Vector{Parameter})
            for i in 1:length(par)
                params[ParamID(pname, i)] = par[i]
            end
        end
    end
    return params
end

# CompEval: a wrapper for a component evaluated on a specific domain
mutable struct CompEval{TComp <: AbstractComponent, TDomain <: AbstractDomain}
    comp::TComp
    domain::TDomain
    params::OrderedDict{ParamID, Parameter}
    counter::Int
    lastvalues::Vector{Float64}
    buffer::Vector{Float64}
    cfixed::Bool

    function CompEval(_comp::AbstractComponent, domain::AbstractDomain)
        # Components internal state may be affected by `prepare!`
        # call.  Avoid overwriting input state with a deep copy.
        comp = deepcopy(_comp)

        params = getparams(comp)
        buffer = prepare!(comp, domain)
        return new{typeof(comp), typeof(domain)}(
            comp, domain, params, 0,
            fill(NaN, length(params)),
            buffer, false)
    end
end

# Fall back method
prepare!(comp::AbstractComponent, domain::AbstractDomain) = fill(NaN, length(domain))
#evaluate!(buffer::Vector{Float64}, comp::T, domain::AbstractDomain, pars...) where T =
#    error("No evaluate! method implemented for $T")

function evaluate_cached(c::CompEval, pvalues::Vector{Float64})
    @assert length(c.params) == length(pvalues)

    # Do we actually need a new evaluation?
    if (any(c.lastvalues .!= pvalues)  ||  (c.counter == 0))
        c.lastvalues .= pvalues
        c.counter += 1
        if !all(.!isnan.(pvalues))
            println("One or more parameter value(s) are NaN:")
            println(pvalues)
            @assert all(.!isnan.(pvalues))
        end
        evaluate!(c.buffer, c.comp, c.domain, pvalues...)
    end
    return c.buffer
end

# Built-in components
#
include("components/utils.jl")
include("components/CDomain.jl")
include("components/SimplePar.jl")
include("components/FuncWrap.jl")
include("components/OffsetSlope.jl")
include("components/Gaussian.jl")
include("components/Lorentzian.jl")

# ====================================================================
# Reducer
#

abstract type AbstractReducer end

prepare!(red::AbstractReducer, domain::AbstractDomain, args::OrderedDict{Symbol, Vector{Float64}}) =
    fill(NaN, length(domain))

struct ExprReducer <: AbstractReducer
    ef::ExprFunction
    function ExprReducer(ef::ExprFunction)
        @assert length(ef.args) == 1
        new(ef)
    end
end

prepare!(red::ExprReducer, domain::AbstractDomain, args::OrderedDict{Symbol, Vector{Float64}}) =
    fill(NaN, length(red.ef.funct(args)))

function evaluate!(buffer::Vector{Float64}, red::ExprReducer,
                   domain::AbstractDomain, args::OrderedDict{Symbol, Vector{Float64}})
    buffer .= red.ef.funct(args)
    nothing
end

macro expr(expr)
    return esc(:(GFit.ExprReducer(GFit.@exprfunc $expr)))
end

struct SumReducer <: AbstractReducer
    list::Vector{Symbol}
end

function evaluate!(buffer::Vector{Float64}, red::SumReducer,
                   domain::AbstractDomain, args::OrderedDict{Symbol, Vector{Float64}})
    buffer .= 0.
    for name in red.list
        buffer .+= args[name]
    end
    nothing
end

mutable struct ReducerEval{T <: AbstractReducer}
    red::T
    counter::Int
    buffer::Vector{Float64}
end

# ====================================================================
# Model
#
const PatchComp = HashVector{Float64}

struct ModelEval
    params::Vector{Parameter}
    ifree::Vector{Int}
    pvalues::Vector{Float64}
    patched::Vector{Float64}
    patchcomps::OrderedDict{Symbol, PatchComp}
    reducer_args::OrderedDict{Symbol, Vector{Float64}}
end

abstract type AbstractMultiModel end

# A model prediction suitable to be compared to a single empirical dataset
mutable struct Model
    parent::Union{Nothing, AbstractMultiModel}
    domain::AbstractDomain
    cevals::OrderedDict{Symbol, CompEval}
    revals::OrderedDict{Symbol, ReducerEval}
    rsel::Symbol
    patchfuncts::Vector{ExprFunction}
    meval::Union{Nothing, ModelEval}

    function Model(domain::AbstractDomain, args...)
        function parse_args(args::AbstractDict)
            out = OrderedDict{Symbol, Union{AbstractComponent, AbstractReducer}}()
            for (name, item) in args
                isa(item, Number)  &&  (item = SimplePar(item))
                @assert isa(name, Symbol)
                @assert isa(item, AbstractComponent) || isa(item, AbstractReducer)
                out[name] = item
            end
            return out
        end

        function parse_args(args::Vararg{Pair})
            out = OrderedDict{Symbol, Union{AbstractComponent, AbstractReducer}}()
            for arg in args
                out[arg[1]] = arg[2]
            end
            return parse_args(out)
        end

        model = new(nothing, domain,
                    OrderedDict{Symbol, CompEval}(),
                    OrderedDict{Symbol, ReducerEval}(),
                    Symbol(""),
                    Vector{ExprFunction}(),
                    nothing)
        evaluate!(model)  # populate meval
        for (name, item) in parse_args(args...)
            model[name] = item
        end
        if length(model.revals) == 0
            model[:default_sum] = SumReducer(collect(keys(model.cevals)))
        end
        return model
    end
end

function ModelEval(model::Model)
    params = Vector{Parameter}()
    ifree = Vector{Int}()
    pvalues = Vector{Float64}()
    patched = Vector{Float64}()
    patchcomps = OrderedDict{Symbol, PatchComp}()
    reducer_args = OrderedDict{Symbol, Vector{Float64}}()

    i = 1
    for (cname, ceval) in model.cevals
        patchcomps[cname] = HashVector{Float64}(patched)
        for (pid, par) in ceval.params
            if !(par.low <= par.val <= par.high)
                s = "Value outside limits for param [$(cname)].$(pid.name):\n" * string(par)
                error(s)
            end
            if (!par.fixed)  &&  (model.cevals[cname].cfixed == 0)
                push!(ifree, i)
            end
            push!(params, par)
            push!(pvalues, par.val)
            push!(patchcomps[cname], pid.name, NaN)
            i += 1
        end
    end
    @assert length(patched) == length(pvalues)

    for (cname, ceval) in model.cevals
        reducer_args[cname] = ceval.buffer
    end
    for (rname, reval) in model.revals
        reducer_args[rname] = reval.buffer
    end

    return ModelEval(params, ifree, pvalues, patched, patchcomps, reducer_args)
end

function evaluate!(model::Model)
    if !isnothing(model.parent)
        evaluate!(model.parent)
    else
        eval1!(model)
        eval2!(model)
        eval3!(model)
        eval4!(model)
    end
    return model
end

function eval1!(model::Model)
    # Update ModelEval structure
    model.meval = ModelEval(model)
end

function eval2!(model::Model; fromparent=false)
    # Evaluate patched parameters
    if !isnothing(model.parent)  &&  !fromparent
        eval2!(model.parent)
    else
        model.meval.patched .= model.meval.pvalues  # copy all values by default
        for pf in model.patchfuncts
            pf.funct(model.meval.patchcomps)
        end
    end
end

function eval3!(model::Model)
    # Evaluate model
    i1 = 1
    for (cname, ceval) in model.cevals
        if length(ceval.params) > 0
            i2 = i1 + length(ceval.params) - 1
            evaluate_cached(ceval, model.meval.patched[i1:i2])
            i1 += length(ceval.params)
        else
            evaluate_cached(ceval, Float64[])
        end
    end

    for (rname, reval) in model.revals
        (rname == model.rsel)  &&  continue  # this should be evaluated as last one
        reval.counter += 1
        evaluate!(reval.buffer, reval.red, model.domain, model.meval.reducer_args)
    end

    if length(model.revals) > 0
        reval = model.revals[model.rsel]
        reval.counter += 1
        evaluate!(reval.buffer, reval.red, model.domain, model.meval.reducer_args)
    end
end

function eval4!(model::Model, unc::Union{Nothing, Vector{Float64}}=nothing)
    # Update values, uncertainties and patched params from ModelEval
    # to the actual Model structure.
    i = 1
    j = 1
    for (cname, ceval) in model.cevals
        for (pid, par) in ceval.params
            if i in model.meval.ifree
                if !isnothing(unc)
                    par.unc = unc[j]
                elseif par.val != model.meval.pvalues[i]
                    # overwrite only if param. value has changed
                    # (i.e. avoid loosing uncertainty when invoking
                    # evaluate!(model)
                    par.unc = NaN
                end
                j += 1
            else
                par.unc = NaN
            end
            par.val = model.meval.pvalues[i]
            par.patched = model.meval.patched[i]
            i += 1
        end
    end
end

function setindex!(model::Model, comp::AbstractComponent, cname::Symbol)
    @assert !haskey(model.revals, cname)  "Name $cname already exists as a reducer name"
    ceval = CompEval(deepcopy(comp), model.domain)
    model.cevals[cname] = ceval
    evaluate!(model)
end

function setindex!(model::Model, reducer::T, rname::Symbol) where T <: AbstractReducer
    @assert !haskey(model.cevals, rname) "Name $rname already exists as a component name"
    model.revals[rname] = ReducerEval{T}(reducer, 1, prepare!(reducer, model.domain, model.meval.reducer_args))
    (model.rsel == Symbol(""))  &&  (model.rsel = rname)
    evaluate!(model)
    return model
end

function select_reducer!(model::Model, rname::Symbol)
    @assert haskey(model.revals, rname) "$rname is not a reducer name"
    model.rsel = rname
end

function geteval(model::Model, name::Symbol)
    if haskey(model.cevals, name)
        return model.cevals[name].buffer
    else
        if haskey(model.revals, name)
            return model.revals[name].buffer
        else
            return Float64[]
        end
    end
end

geteval(model::Model) = geteval(model, model.rsel)

function patch!(model::Model, exfunc::ExprFunction)
    push!(model.patchfuncts, exfunc)
    evaluate!(model)
    return model
end
include("macro_patch.jl")

function isfixed(model::Model, cname::Symbol)
    @assert cname in keys(model.cevals) "Component $cname is not defined"
    return model.cevals[cname].cfixed
end

function freeze(model::Model, cname::Symbol)
    @assert cname in keys(model.cevals) "Component $cname is not defined"
    model.cevals[cname].cfixed = true
    evaluate!(model)
    nothing
end

function thaw(model::Model, cname::Symbol)
    @assert cname in keys(model.cevals) "Component $cname is not defined"
    model.cevals[cname].cfixed = false
    evaluate!(model)
    nothing
end

Base.keys(p::Model) = collect(keys(p.cevals))
Base.haskey(p::Model, name::Symbol) = haskey(p.cevals, name)
function Base.getindex(model::Model, name::Symbol)
    if name in keys(model.cevals)
        return model.cevals[name].comp
    elseif name in keys(model.revals)
        return model.revals[name].red
    end
    error("Name $name not defined")
end
domain(model::Model) = model.domain
(model::Model)() = geteval(model)
(model::Model)(name::Symbol) = geteval(model, name)

include("multimodel.jl")
include("minimizers.jl")
include("fit.jl")
include("show.jl")

end


text:
module Lean
using LEAN_Community_jll

test = """
example (m n : nat) m + n = n + m :=
by simp [nat.add_comm]
"""
function check(test)
    lean() do lean
        run(`$lean --version`)
    end
end

process_exited(check(test))

end


text:
# -------------------------------------------------------
#         Measurements
# -------------------------------------------------------

# χ(qy,qx,iw) = C(qy, qx, iw) * Δτ/(N*M)
function measure_chi_dynamic(conf::AbstractArray{Float64, 3}; Δτ::Float64=0.1)
  N = size(conf, 2)
  M = size(conf, 3)
  Δτ/(N*M) * measure_phi_correlations(conf)
end

# don't sum over op components
function measure_chi_dynamic_components(conf::AbstractArray{Float64, 3}; Δτ::Float64=0.1)
  N = size(conf, 2)
  M = size(conf, 3)
  Δτ/(N*M) * measure_phi_correlations_components(conf)
end

# This is the rotated bosonic susceptibility with momenta relative to Q,
# i.e. χ_rotated(qy,qx,iw) = C(qy-pi, qx-pi, iw) * Δτ/(N*M)
#
# first element of C in w direction is 0 freq, last element of C in qx and qy direction is (-pi,-pi)
# (fftfreq(#elements, sampling rate = 1/a = 1 in our case))
# just reverse qx and qy dimensions to get q-Q dependence
function measure_rotated_chi_dynamic(conf::AbstractArray{Float64, 3}; Δτ::Float64=0.1)
  mapslices(x->rotr90(rotr90(x)), measure_chi_dynamic(conf; Δτ=Δτ), [1,2])
end

# χ_rotated(Q,0) = χ_rotated(pi,pi,0) = χ(0,0) = C(0,0) * Δτ/(N*M)
function measure_chi_static(conf::AbstractArray{Float64, 3}; Δτ::Float64=0.1)
  measure_chi_dynamic(conf; Δτ=Δτ)[1,1,1] # mean component == 0 component == idx 1 element in Julia
end

# Binder factors (m^2, m^4)
function measure_binder_factors(conf::AbstractArray{Float64, 3})
  m = mean(conf,dims=(2,3)) # average over sites and timeslices
  m2 = dot(m,m)
  return (m2, m2*m2)
end

function measure_op(conf::AbstractArray{Float64, 3})
  mean_abs_op = mean(abs.(conf))
  mean_op = vec(mean(conf,dims=(2,3)))
  return mean_abs_op, mean_op
end

# bosonic correlation function C(qy,qx,iw), not normalized
function measure_phi_correlations(conf::AbstractArray{Float64, 3})
  opdim, sites, slices = size(conf)
  L = Int(sqrt(sites))
  phiFT = rfft(reshape(conf,opdim,L,L,slices),[2,3,4])
  phiFT .*= conj(phiFT)
  n = floor(Int, L/2+1)
  nt = floor(Int, slices/2+1)
  return real(sum(phiFT,dims=1))[1,:,1:n,1:nt] # sum over op components
end

# measure correlations for every op component (don't sum over them)
function measure_phi_correlations_components(conf::AbstractArray{Float64, 3})
  opdim, sites, slices = size(conf)
  L = Int(sqrt(sites))
  phiFT = rfft(reshape(conf,opdim,L,L,slices),[2,3,4])
  phiFT .*= conj(phiFT)
  n = floor(Int, L/2+1)
  nt = floor(Int, slices/2+1)
  return real(phiFT)[:,:,1:n,1:nt] # do not yet sum over op components
end

# -------------------------------------------------------
#         Postprocessing/Analysis
# -------------------------------------------------------

"""
Returns qys, qxs and ωs corresponding to output of measure_phi_correlations etc.
"""
function get_momenta_and_frequencies(L::Int, M::Int, Δτ::Float64=0.1)
  qs = 2*pi*rfftfreq(L) # 2*pi because rfftfreq returns linear momenta
  ωs = 2*pi*rfftfreq(M, 1/Δτ)
  return qs, qs, ωs
end
get_momenta_and_frequencies(C::AbstractArray{Float64, 3}, Δτ::Float64=0.1) = get_momenta_and_frequencies(Int(sqrt(size(C,2))), size(C,3), Δτ)
get_momenta_and_frequencies(mc::AbstractDQMC) = get_momenta_and_frequencies(mc.l.L, mc.p.slices, mc.p.delta_tau)

# ------------------------------------------------------------------
#         "Macroscopic methods" - measure MC time series
# ------------------------------------------------------------------
function measure_chi_dynamics(confs::AbstractArray{Float64, 4})
    num_confs = size(confs,4)
    N = size(confs, 2)
    M = size(confs, 3)

    chi_dyn_symm = Observable(Array{Float64, 3}, "chi_dyn_symm"; alloc=num_confs)
    chi_dyn = Observable(Array{Float64, 3}, "chi_dyn"; alloc=num_confs)
    @inbounds @views for i in 1:num_confs
        chi = measure_chi_dynamic(confs[:,:,:,i])
        push!(chi_dyn, chi)
        chi = (permutedims(chi, [2,1,3]) + chi)/2 # C4 is basically flipping qx and qy (which only go from 0 to pi since we perform a real fft.)
        push!(chi_dyn_symm, chi)
    end

    return chi_dyn, chi_dyn_symm
end

function measure_binder(confs::AbstractArray{Float64, 4})
    num_confs = size(confs,4)
    N = size(confs, 2)
    M = size(confs, 3)

    m2s = Vector{Float64}(undef, num_confs)
    m4s = Vector{Float64}(undef, num_confs)
    @inbounds @views for i in 1:num_confs
        m = mean(confs[:,:,:,i],dims=(2,3))
        m2s[i] = dot(m, m)
        m4s[i] = m2s[i]*m2s[i]
    end

    m2ev2 = mean(m2s)^2
    m4ev = mean(m4s)

    binder = Observable(Float64, "binder")
    push!(binder, m4ev/m2ev2)

    return binder
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
# Copyright 2022 RelationalAI, Inc.
#
# 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.

"""
    RAI

The `RAI` module provides functionality for accessing the RelationalAI REST
APIs.
"""
module RAI

const PROJECT_VERSION = v"0.0.1"

export
    AccessToken,
    ClientCredentials,
    Config,
    Context,
    HTTPError,
    load_config

export
    create_database,
    delete_database,
    get_database,
    list_databases

export
    create_engine,
    delete_engine,
    get_engine,
    list_engines

export
    delete_model,
    get_model,
    list_models,
    load_model

export
    create_oauth_client,
    delete_oauth_client,
    get_oauth_client,
    list_oauth_clients

export
    create_user,
    disable_user,
    enable_user,
    get_user,
    list_users,
    update_user

export
    delete_model,
    get_model,
    list_models

export
    exec,
    exec_async,
    list_edbs,
    load_csv,
    load_json

export
    get_transaction,
    get_transaction_metadata,
    get_transaction_problems,
    get_transaction_results

export
    Relation,
    TransactionResult,
    getcolumn,
    getrow,
    schema,
    show_problems,
    show_result

include("creds.jl")
include("config.jl")
include("rest.jl")
include("api.jl")
include("results.jl")

end

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
# # GARTEUR SM-AG19 Testbed: Modal analysis

# ## Description

# This virtual test application is based on the test article
# used by the GARTEUR Structures & Materials Action Group 19
# which organized a Round Robin exercise where 12 European laboratories
# tested a single structure between 1995 and 1997. The benchmark structure
# was a laboratory structure built to simulate the dynamic behaviour
# of an aeroplane. The structure was initially built for a benchmark
# study on experimental modal analysis conducted by the
# Structures and Materials Action Group (SM-AG19) of the Group
# for Aeronautical Research and Technology in EURope (GARTEUR).
# The test-bed was designed and manufactured by ONERA, France.

# ![](IMAC97photo.png)

# ### References

# - [GARTEUR] Ground Vibration Test Techniques, compiled by A Gravelle, GARTEUR
#   Structures & Materials Action Group 19 Technical report TP-115, 1999.
# - [BW] Etienne Balmes, Jan R. Wright, GARTEUR GROUP ON GROUND VIBRATION
#   TESTING | RESULTS FROM THE TEST OF A SINGLE STRUCTURE BY 12 LABORATORIES IN
#   EUROPE, Proceedings of DETC'97, 1997 ASME Design Engineering Technical
#   Conferences, September 14-17, 1997, Sacramento, California.
# - [3M] 3M(TM) Viscoelastic Damping Polymer 112 Series,  Technical Data, May 2017.

# ## Goals

# - Show how to construct model from multiple connected beams.
# - Demonstrate the use of massless connectors.
# - Demonstrate the use of point masses.
# - Demonstrate the use of grounded springs.
# - Illustrate verification of the solution of the free vibration problem.

# ## Geometry of the testbed airplane.

# The aluminum testbed was a rather simple structure which was reasonably
# dynamically representative of a simple airplane structure [GARTEUR](@ref
# References). It was composed of several beams simulating a fuselage with wings
# and a tail. Wing tip drums allowed to adjust bending and torsion frequencies
# similarly to airplane ones, with some very close modal frequencies.

# ![](garteur-geom.png)

# The script included below defines the geometry of the structure, the
# cross-sectional properties, the connectivity, and the location of the nodes.
include("garteur_geometry_tut.jl")

# The geometry is visualized in the tutorial
# [`garteur_geometry_vis_tut.md`](garteur_geometry_vis_tut.jl).

# ## Material

# Material properties can be now used to create a material: isotropic elasticity
# model of the `FinEtoolsDeforLinear` package is instantiated.
using FinEtoolsDeforLinear
# The material of the structure is aluminum.
# The elastic modulus:
E = 70000.0 * phun("MPa")
nu = 0.31;
# The mass density:
rho = 2700 * phun("kg/m^3")
alu = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# The material of the constraining layer on top of the viscoelastic tape. It was
# aluminum.
layer = MatDeforElastIso(DeforModelRed3D, rho, alu.E, alu.nu, 0.0)

# Material for the massless connectors has the mass density set to zero;
# otherwise it has the same properties as the aluminum material  of the
# structure.

massless = MatDeforElastIso(DeforModelRed3D, 0.0, alu.E, alu.nu, 0.0)

# This simple function returns material based on the label of the beam elements.
getmaterial(labl) = begin
    if labl == 6
        return layer
    elseif labl == 7 || labl == 8
        return massless
    end
    return alu
end

# This is the assumed stifffness of the bungee cords (each one separately).
bungeecoefficient = 4000*phun("N/m");

# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to 0.
# This is that geometry field.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation field, three unknown rotations per node are represented
# with a rotation matrix, in total nine numbers. The utility function
# `initial_Rfield`
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node.
dchi = NodalField(zeros(size(fens.xyz, 1), 6))

# There are no support conditions.
applyebc!(dchi)
# The  the number of free(unknown) degrees of freedom is equal to the total
# number of degrees of freedom in the system.
numberdofs!(dchi);

# ## Identify support points and locations of sensors

# Suspension points
suspln = selectnode(fens; box = initbox!(Float64[], vec([0.0*L 0.0*L 0.805*L])), inflate = tolerance)
susprn = selectnode(fens; box = initbox!(Float64[], vec([0.0*L -0.0*L 0.805*L])), inflate = tolerance)
suspbn = selectnode(fens; box = initbox!(Float64[], vec([-2.0*L 0.0*L 0.0*L])), inflate = tolerance)

# The sensors at the tip of the left and right wing drum
sensor112n = selectnode(fens; box = initbox!(Float64[], vec([+1.8*L 9.8*L 0.96*L])), inflate = tolerance)
sensor12n = selectnode(fens; box = initbox!(Float64[], vec([+1.8*L -9.8*L .96*L])), inflate = tolerance)
sensor111n = selectnode(fens; box = initbox!(Float64[], vec([-1.8*L 9.8*L 0.96*L])), inflate = tolerance)
sensor11n = selectnode(fens; box = initbox!(Float64[], vec([-1.8*L -9.8*L .96*L])), inflate = tolerance)

# The joint between the horizontal and vertical tail parts
sensor202n = selectnode(fens; box = initbox!(Float64[], vec([-8*L 0 3.8*L])), inflate = tolerance)

# ## Assemble the global discrete system

# For disambiguation we will refer to the stiffness and mass functions by
# qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
using FinEtoolsFlexStructures.FEMMCorotBeamModule
CB = FEMMCorotBeamModule

# Note that we have an array of finite element sets. We compute the matrices for
# each set separately and accumulate them into the final overall matrix. Thus
# we can construct the stiffness and mass matrix as follows.
using  SparseArrays

K, M = let
    K = spzeros(dchi.nfreedofs, dchi.nfreedofs)
    M = spzeros(dchi.nfreedofs, dchi.nfreedofs)
    for fes in fesa
        labl  = fes.label[1]
        femm = CB.FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), getmaterial(labl));
        K += CB.stiffness(femm, geom0, u0, Rfield0, dchi);
        M += CB.mass(femm, geom0, u0, Rfield0, dchi);
    end
    K, M
end

# ## Additional concentrated masses.

using LinearAlgebra

using FinEtoolsFlexStructures.FEMMPointMassModule
PM = FEMMPointMassModule

# There is at the sensor on the tail.
femmcm1 =  PM.FEMMPointMass(IntegDomain(FESetP1(reshape([sensor202n;], 1, 1)), PointRule()), FFltMat(2*L*L/5*L/5*2*rho*LinearAlgebra.I(3)));

# These are the forward/interior locations on the wing drums.

mass1n = selectnode(fens; box = initbox!(Float64[], vec([1.8*L 9.2*L .96*L])), inflate = tolerance)
mass2n = selectnode(fens; box = initbox!(Float64[], vec([1.8*L -9.2*L .96*L])), inflate = tolerance)
femmcm2 =  PM.FEMMPointMass(IntegDomain(FESetP1(reshape([mass1n; mass2n;], 2, 1)), PointRule()), FFltMat(0.2*phun("kg")*LinearAlgebra.I(3)));

Mp = PM.mass(femmcm1, geom0, u0, Rfield0, dchi) + PM.mass(femmcm2, geom0, u0, Rfield0, dchi);

# ## Bungee supports

using LinearAlgebra

using FinEtoolsFlexStructures.FEMMPointGroundedSpringModule
BS = FEMMPointGroundedSpringModule

# There are three suspension points at the top of the fuselage. We assume that
# these bungee supports exert only reaction in the vertical direction.
femmbs =  BS.FEMMPointGroundedSpring(IntegDomain(FESetP1(reshape([suspln; susprn; suspbn;], 3, 1)), PointRule()),
FFltMat([bungeecoefficient*[0;0;1]*[0;0;1]' 0*LinearAlgebra.I(3); 0*LinearAlgebra.I(3) 0*LinearAlgebra.I(3)]));

Kb = BS.stiffness(femmbs, geom0, u0, Rfield0, dchi)

Kt = K + Kb
Mt = M + Mp

# We can compare the size of the stiffness matrix with the number of degrees of
# freedom that are unknown (20).
@show size(Kt)

# ## Solve the free-vibration problem

# Find this many natural frequencies:
neigvs = 20

# Since the structure is free-floating, mass shifting must be employed to obtain
# the solution with the singular stiffness matrix. We are simply guessing a
# frequency between zero and the first fundamental frequency.
oshift = (2*pi*0.5)^2;

# The Arnoldi algorithm implemented in the well-known `Arpack` package is used
# to solve the generalized eigenvalue problem with the sparse matrices. As is
# common in structural dynamics, we request the smallest eigenvalues in
# absolute value (`:SM`).
using Arpack
evals, evecs, nconv = eigs(Kt + oshift * Mt, Mt; nev=neigvs, which=:SM, explicittransform = :none);
# First  we should check that the requested eigenvalues actually converged:
@show nconv == neigvs

# The eigenvalues (i. e. the squares of the angular frequencies) are returned in
# the vector `evals`. The mode shapes constitute the columns of the matrix
# `evecs`.
@show size(evecs)
# The natural frequencies are obtained from the squares of the angular
# frequencies. We note the use of `sqrt.` which broadcast the square root over
# the array `evals`.
fs = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);

# ## Comparison of computed and analytical results

# Set of modes measured by participant C.

# | Mode | Description | Frequency |
# | ----- | ----- | ----- |
# | 1.-6. |  "Rigid body" modes | 0 |
# | 7. |  Two node bending | 6.37 Hz |
# | 8. |  Global fuselage rotation | 16.10 Hz |
# | 9. |  First antisymmetric wing torsion | 33.13 Hz |
# | 10. |  First symmetric wing torsion | 33.53 Hz |
# | 11. |  Three node bending | 35.65 Hz |

# The approximate and analytical frequencies are now reported.
sigdig(n) = round(n * 1000) / 1000

println("Frequencies 7 and higher")
println("Approximate: $(sigdig.(fs[7:end])) [Hz]")
println("Participant C experimental: $([6.37, 16.10, 33.13, 33.53, 35.65, 48.38, 49.43, 55.08]) [Hz]")

# ## Visualize vibration modes

# The animation will show one of the vibration modes overlaid on the undeformed
# geometry. The configuration during the animation needs to reflect rotations.
# The function `update_rotation_field!` will update the rotation field given a
# vibration mode.
using FinEtoolsFlexStructures.RotUtilModule: update_rotation_field!

# The visualization utilities take advantage of the PlotlyJS library.
using PlotlyJS
using VisualStructures: plot_space_box, plot_solid, render, react!, default_layout_3d, save_to_json

# The magnitude of the vibration modes (displacements  and rotations) will be
# amplified with this scale factor:
scale = 0.3

# The animation is encapsulated in a little function:
vis(mode) = begin
    tbox = plot_space_box(reshape(inflatebox!(boundingbox(fens.xyz), 4*L), 2, 3))
    tenv0 = tbox
    for fes in fesa
        t = plot_solid(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(125, 155, 125)", opacity=0.3);
        tenv0 = cat(tenv0, t; dims=1)
    end
    plots = tenv0
    layout = default_layout_3d(;title = "Mode $(mode), $(sigdig.(fs[mode])) [Hz]")
    layout[:scene][:aspectmode] = "data"
    pl = render(plots; layout=layout)
    sleep(0.115)
    for xscale in scale .* sin.(collect(0:1:89) .* (2 * pi / 21))
        scattersysvec!(dchi, xscale .* evecs[:, mode])
        u1 = deepcopy(u0)
        u1.values .= dchi.values[:, 1:3]
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        plots = tenv0
        for fes in fesa
            tenv1 = plot_solid(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 55, 125)");
            plots = cat(plots, tenv1; dims=1)
        end
        pl.plot.layout[:title] = "Mode $(mode), $(sigdig.(fs[mode])) [Hz]"
        react!(pl, plots, pl.plot.layout)
        sleep(0.08)
    end
end

# This is the mode that will be animated:

vis(7)

nothing


text:
# Constructing Gaussian quadrature weights for an
# arbitrary weight function and integration bounds.

# for numerical stability, we apply the usual Lanczos
# Gram–Schmidt procedure to the basis {T₀,T₁,T₂,…} of
# Chebyshev polynomials on [-1,1] rather than to the
# textbook monomial basis {1,x,x²,…}.

# evaluate Chebyshev polynomial p(x) with coefficients a[i]
# by a Clenshaw recurrence.
function chebeval(x, a)
    if length(a) ≤ 2
        length(a) == 1 && return a[1] + x * zero(a[1])
        return a[1]+x*a[2]
    end
    bₖ = a[end-1] + 2x*a[end]
    bₖ₊₁ = oftype(bₖ, a[end])
    for j = lastindex(a)-2:-1:2
        bⱼ = a[j] + 2x*bₖ - bₖ₊₁
        bₖ, bₖ₊₁ = bⱼ, bₖ
    end
    return a[1] + x*bₖ - bₖ₊₁
end

# if a[i] are coefficients of Chebyshev series, compute the coefficients xa (in-place)
# of the series multiplied by x, using recurrence xTₙ = 0.5 (Tₙ₊₁+Tₙ₋₁) for n > 0
function chebx!(xa, a)
    resize!(xa, length(a)+1)
    xa .= 0
    for n = 2:lastindex(a)
        c = 0.5*a[n]
        xa[n-1] += c
        xa[n+1] += c
    end
    if !isempty(a)
        xa[2] += a[1]
    end
    return xa
end

"""
    gauss(W, N, a, b; rtol=sqrt(eps), quad=quadgk)

Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to
integrate functions on the interval `(a, b)` multiplied by the weight function
`W(x)`.  That is, `sum(w .* f.(x))` approximates the integral `∫ W(x)f(x)dx`
from `a` to `b`.

This function performs `2N` numerical integrals of polynomials against `W(x)`
using the integration function `quad` (defaults to `quadgk`) with relative tolerance `rtol`
(which defaults to half of the precision `eps` of the endpoints).
This is followed by an O(N²) calculations. So, using a large order `N` is expensive.

If `W` has lots of singularities that make it hard to integrate numerically,
you may need to decrease `rtol`.   You can also pass in a specialized quadrature routine
via the `quad` keyword argument, which should accept arguments `quad(f,a,b,rtol=_,atol=_)`
similar to `quadgk`.  (This is useful if your weight function has discontinuities, in which
case you might want to break up the integration interval at the discontinuities.)
"""
function gauss(W, N, a::Real,b::Real; rtol::Real=sqrt(eps(typeof(float(b-a)))), quad=quadgk)
    (isfinite(a) && isfinite(b)) || throw(ArgumentError("a finite interval is required"))
    return _gauss(W, N, a, b, rtol, quad)
end

function _gauss(W, N, a, b, rtol, quad)
    # Uses the Lanczos recurrence described in Trefethen & Bau,
    # Numerical Linear Algebra, to find the `N`-point Gaussian quadrature
    # using O(N) integrals and O(N²) operations, applied to Chebyshev basis:
    xscale = 2.0/(b-a) # scaling from (a,b) to (-1,1)
    T = typeof(xscale)
    α = zeros(T, N)
    β = zeros(T, N)
    q₀ = sizehint!(T[0], N+1) # 0 polynomial
    wint = first(quad(W, a, b, rtol=rtol))
    (wint isa Real && wint > 0) ||
        throw(ArgumentError("weight W must be real and positive"))
    atol = rtol*wint
    q₁ = sizehint!(T[1/sqrt(wint)],N+1) # normalized constant polynomial
    v = copy(q₀)
    for n = 1:N
        chebx!(v, q₁) # v = x * q₁
        α[n] = first(quad(a, b, rtol=rtol, atol=atol) do x
            t = (x-a)*xscale - 1
            W(x) * chebeval(t, q₁) * chebeval(t, v)
        end)
        n == N && break
        for j = 1:length(q₀); v[j] -= β[n]*q₀[j]; end
        for j = 1:length(q₁); v[j] -= α[n]*q₁[j]; end
        β[n+1] = sqrt(first(quad(a, b, rtol=rtol, atol=atol) do x
            W(x) * chebeval((x-a)*xscale - 1, v)^2
        end))
        v .*= inv(β[n+1])
        q₀,q₁,v = q₁,v,q₀
    end

    # TODO: handle BigFloat etcetera — requires us to update eignewt() to
    #       support nonzero diagonal entries.
    E = eigen(SymTridiagonal(α, β[2:N]))

    w = E.vectors[1,:]
    w .= wint .* abs2.(w)
    return ((E.values .+ 1) ./ xscale .+ a, w)
end


text:
using Query
using DataFrames

struct Person
    Name::String
    Friends::Vector{String}
end

source = Array{Person}(undef, 0)
push!(source, Person("John", ["Sally", "Miles", "Frank"]))
push!(source, Person("Sally", ["Don", "Martin"]))

result = @from i in source begin
         @where length(i.Friends) > 2
         @select {i.Name, Friendcount=length(i.Friends)}
         @collect
end

println(result)

result = @from i in source begin
         @where length(i.Friends) > 2
         @select {i.Name, Friendcount=length(i.Friends)}
         @collect DataFrame
end

println(result)


text:
using JuliaInterpreter
using JuliaInterpreter: JuliaStackFrame, JuliaProgramCounter, @lookup
using JuliaInterpreter: finish_and_return!, @lookup, evaluate_call!, _step_expr!,
                        do_assignment!, getlhs, isassign, pc_expr, handle_err
using Base.Meta: isexpr
using Test, Random

# Execute a frame using Julia's regular compiled-code dispatch for any :call expressions
runframe(frame, pc=frame.pc[]) = Some{Any}(finish_and_return!(Compiled(), frame, pc))

# Execute a frame using the interpreter for all :call expressions (except builtins & intrinsics)
function runstack(frame::JuliaStackFrame, pc=frame.pc[])
    stack = JuliaStackFrame[]
    return Some{Any}(finish_and_return!(stack, frame, pc))
end

## For juliatests.jl

isdocexpr(ex) = isexpr(ex, :macrocall) && (a = ex.args[1]; a isa GlobalRef && a.mod == Core && a.name == Symbol("@doc"))

function read_and_parse(filename)
    src = read(filename, String)
    ex = Base.parse_input_line(src; filename=filename)
end

"""
    docexprs = lower_incrementally(f, mod::Module, ex::Expr)

Lower the sub-expressions of `ex` one-by-one, passing the `frame`s to `f(frame)`.
Any docstring-expressions encountered while processing module `M` (which might be `mod`
or one of its sub-modules) get returned in `docexprs[M]`.
"""
function lower_incrementally(@nospecialize(f), mod::Module, ex::Expr)
    docexprs = Dict{Module,Vector{Expr}}()
    aborts = LineNumberNode[]
    lower_incrementally!(f, docexprs, aborts, mod, ex)
    return docexprs, aborts
end

lower_incrementally!(@nospecialize(f), docexprs, aborts, mod::Module, ex::Expr) =
    lower_incrementally!(f, docexprs, aborts, Expr(:block), mod, ex)

function lower_incrementally!(@nospecialize(f), docexprs, aborts, lex:: Expr, mod::Module, ex::Expr)
    # lex is the expression we'll lower; it will accumulate LineNumberNodes and a
    # single top-level expression. We split blocks, module defs, etc.
    if ex.head == :toplevel || ex.head == :block
        lower_incrementally!(f, docexprs, aborts, lex, mod, ex.args)
    elseif ex.head == :module
        modname = ex.args[2]::Symbol
        newmod = isdefined(mod, modname) ? getfield(mod, modname) : Core.eval(mod, :(module $modname end))
        lower_incrementally!(f, docexprs, aborts, lex, newmod, ex.args[3])
    elseif isdocexpr(ex) && length(ex.args) >= 4
        docexs = get(docexprs, mod, nothing)
        if docexs === nothing
            docexs = docexprs[mod] = Expr[]
        end
        push!(docexs, ex)
        body = ex.args[4]
        if isa(body, Expr)
            lower_incrementally!(f, docexprs, aborts, lex, mod, body)
        end
    else
        # For map(x->x^2, a) we need to split out the anonymous function so that it
        # gets defined prior to using it (essentially a world-age issue for inference)
        if ex.head == :call || ex.head == :(=) && isa(ex.args[2], Expr) && ex.args[2].head == :call
            fas = split_anonymous!(ex)
            for fa in fas
                lex1 = copy(lex)
                push!(lex1.args, fa)
                lower!(f, docexprs, aborts, mod, lex1)
            end
        end
        push!(lex.args, ex)
        lower!(f, docexprs, aborts, mod, lex)
        empty!(lex.args)
    end
    return docexprs, aborts
end

function lower_incrementally!(@nospecialize(f), docexprs, aborts, lex, mod::Module, args::Vector{Any})
    for a in args
        if isa(a, Expr)
            lower_incrementally!(f, docexprs, aborts, lex, mod, a)
        else
            push!(lex.args, a)
        end
    end
end

function lower!(@nospecialize(f), docexprs, aborts, mod::Module, ex::Expr)
    lwr = Meta.lower(mod, ex)
    if isexpr(lwr, :thunk)
        frame = JuliaInterpreter.prepare_thunk(mod, lwr)
        cts = Test.get_testset()
        ret = Base.invokelatest(f, frame)  # if previous thunks define new methods, we need to update world age
        if isa(ret, Aborted)
            push!(aborts, abortline(ex))
        end
        tts = Test.get_testset()
        while tts != cts
            ttsinner = Test.pop_testset()
            tts = Test.get_testset()
            Test.record(tts, ttsinner)
        end
    elseif isa(lwr, Expr) && (lwr.head == :export || lwr.head == :using || lwr.head == :import)
    elseif isa(lwr, Symbol) || isa(lwr, Nothing)
    else
        @show mod ex
        error("lowering did not produce a :thunk Expr")
    end
    return docexprs, aborts
end

split_anonymous!(ex) = split_anonymous!(Expr[], ex)
function split_anonymous!(fs, ex)
    for (i,a) in enumerate(ex.args)
        if isa(a, Expr)
            if a.head == :->
                gs = gensym()
                push!(fs, Expr(:(=), gs, a))
                ex.args[i] = gs
            else
                split_anonymous!(fs, a)
            end
        end
    end
    return fs
end

struct Aborted end   # for signaling that some statement or test blocks were interrupted

function abortline(ex::Expr)
    if ex.head == :macrocall
        abortline(ex.args[2])
    elseif ex.head == :block
        i = findfirst(x->isa(x, LineNumberNode), ex.args)
        if i === nothing
            length(ex.args) == 1 && return abortline(ex.args[1])
            error("no LineNumberNodes in ", ex, "\nwith ", length(ex.args), " args")
        end
        abortline(ex.args[i])
    elseif ex.head == :let || ex.head == :for || ex.head == :if
        abortline(ex.args[2])
    elseif ex.head == :call
        error("aborted while running ", ex)
    else
        error("unhandled expr head ", ex.head, ":\n", ex)
    end
end
abortline(lnn::LineNumberNode) = lnn

"""
    ret = limited_finish_and_return!(stack, frame, nstmts, istoplevel::Bool)

Run `frame` until execution terminates or more than `nstmts` have been executed,
and pass back the computed return value. `stack` controls call evaluation; `stack = Compiled()`
evaluates :call expressions by normal dispatch, whereas a vector of `JuliaStackFrames`
will use recursive interpretation.
"""
function limited_finish_and_return!(stack, frame, nstmts::Int, pc::JuliaProgramCounter, istoplevel::Bool)
    refnstmts = Ref(nstmts)
    limexec!(s,f) = limited_exec!(s, f, refnstmts)
    # The following is like finish!, except we intercept :call expressions so that we can run them
    # with limexec! rather than the default finish_and_return!
    while nstmts > 0
        stmt = pc_expr(frame, pc)
        if isa(stmt, Expr)
            if stmt.head == :call
                refnstmts[] = nstmts
                try
                    rhs = evaluate_call!(stack, frame, stmt, pc; exec! = limexec!)
                    if isassign(frame, pc)
                        lhs = getlhs(pc)
                        do_assignment!(frame, lhs, rhs)
                    end
                catch err
                    return handle_err(frame, err), refnstmts[]
                end
                nstmts = refnstmts[]
                new_pc = pc + 1
            elseif stmt.head == :(=) && isexpr(stmt.args[2], :call)
                refnstmts[] = nstmts
                try
                    rhs = evaluate_call!(stack, frame, stmt.args[2], pc; exec! = limexec!)
                    do_assignment!(frame, stmt.args[1], rhs)
                catch err
                    return handle_err(frame, err), refnstmts[]
                end
                nstmts = refnstmts[]
                new_pc = pc + 1
            else
                new_pc = _step_expr!(stack, frame, stmt, pc, istoplevel)
                nstmts -= 1
            end
        else
            new_pc = _step_expr!(stack, frame, stmt, pc, istoplevel)
            nstmts -= 1
        end
        new_pc == nothing && break
        pc = new_pc
    end
    frame.pc[] = pc
    # Handle the return
    stmt = pc_expr(frame, pc)
    isexpr(stmt, :return) && return @lookup(frame, (stmt::Expr).args[1]), nstmts
    nstmts == 0 && return Aborted(), nstmts
    error("unexpected return statement ", stmt)
end
limited_finish_and_return!(stack, frame, nstmts::Int, istoplevel::Bool) =
    limited_finish_and_return!(stack, frame, nstmts, frame.pc[], istoplevel)

function limited_exec!(stack, newframe, refnstmts)
    ret, nleft = limited_finish_and_return!(stack, newframe, refnstmts[], newframe.pc[], false)
    refnstmts[] = nleft
    return ret
end

### Functions needed on workers for running tests

function configure_test()
    # To run tests efficiently, certain methods must be run in Compiled mode,
    # in particular those that are used by the Test infrastructure
    cm = JuliaInterpreter.compiled_methods
    empty!(cm)
    push!(cm, which(Test.eval_test, Tuple{Expr, Expr, LineNumberNode}))
    push!(cm, which(Test.get_testset, Tuple{}))
    push!(cm, which(Test.push_testset, Tuple{Test.AbstractTestSet}))
    push!(cm, which(Test.pop_testset, Tuple{}))
    for f in (Test.record, Test.finish)
        for m in methods(f)
            push!(cm, m)
        end
    end
    push!(cm, which(Random.seed!, Tuple{Union{Integer,Vector{UInt32}}}))
    push!(cm, which(copy!, Tuple{Random.MersenneTwister, Random.MersenneTwister}))
    push!(cm, which(copy, Tuple{Random.MersenneTwister}))
    push!(cm, which(Base.include, Tuple{Module, String}))
    push!(cm, which(Base.show_backtrace, Tuple{IO, Vector}))
    push!(cm, which(Base.show_backtrace, Tuple{IO, Vector{Any}}))
end

function runtest(frame, nstmts)
    stack = JuliaStackFrame[]
    # empty!(JuliaInterpreter.framedict)
    # empty!(JuliaInterpreter.genframedict)
    ret, nstmts = limited_finish_and_return!(stack, frame, nstmts, true)
    return ret
end

function dotest(test, fullpath, nstmts)
    println("Working on ", test, "...")
    mod = Core.eval(Main, :(
        module JuliaTests
        using Test, Random
        end
        ))
    ex = read_and_parse(fullpath)
    isexpr(ex, :error) && @error "error parsing $test: $ex"
    # so `include` works properly, we have to set up the relative path
    # oldpath = current_task().storage[:SOURCE_PATH]
    local ts, aborts
    try
        # current_task().storage[:SOURCE_PATH] = fullpath
        cd(dirname(fullpath)) do
            ts = Test.DefaultTestSet(test)
            Test.push_testset(ts)
            docexprs, aborts = lower_incrementally(frame->runtest(frame, nstmts), mod, ex)
            # Core.eval(JuliaTests, ex)
        end
    finally
        # current_task().storage[:SOURCE_PATH] = oldpath
    end
    println("Finished ", test)
    return ts, aborts
end

# Run a test in process id 1 (i.e., the main Julia process)
function dotest1(test, fullpath, nstmts)
    println("Working on ", test, "...")
    mod = Core.eval(Main, :(
        module JuliaTests
        using Test, Random
        end
        ))
    ex = read_and_parse(fullpath)
    isexpr(ex, :error) && @error "error parsing $test: $ex"
    # so `include` works properly, we have to set up the relative path
    oldpath = current_task().storage[:SOURCE_PATH]
    local ts, aborts
    try
        current_task().storage[:SOURCE_PATH] = fullpath
        cd(dirname(fullpath)) do
            ts = Test.DefaultTestSet(test)
            Test.push_testset(ts)
            docexprs, aborts = lower_incrementally(frame->runtest(frame, nstmts), mod, ex)
            # Core.eval(JuliaTests, ex)
        end
    finally
        current_task().storage[:SOURCE_PATH] = oldpath
    end
    println("Finished ", test)
    return ts, aborts
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CURL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CURL")
JLLWrappers.@generate_main_file("CURL", UUID("b21e61f3-bafc-59ac-ab14-4c5c62d6588d"))
end  # module CURL_jll


text:
#!/usr/bin/env julia-0.5

include("db.jl")

#
# File discovery
#

# Recursively scan a directory for matching files
function scan(dir, extensions)
    sources = Vector{String}()
    for entry in readdir(dir)
        path = joinpath(dir, entry)
        if isdir(path)
            append!(sources, scan(path, extensions))
        elseif isfile(path) && any(ex->endswith(path, ex), extensions)
            push!(sources, path)
        end
    end
    return sources
end

#
# Tree construction
#

using ProgressMeter

function create_functions!(parent::Node{FileNode}, repo::String)
    path = parent[:path]
    lines = parent[:lines]

    # parse functions and start line (end line = start next function)
    cmd = pipeline(`ctags -x $path`, `sort -k3 -n`)
    functions = Dict{String, Tuple{Int, Int}}()
    previous = nothing
    for line in readlines(cmd)
        ctags_re = r"^(\w+)\s+(\w+)\s+(\d+) [^ ]+ (.+)$"
        m = match(ctags_re, line)
        if m != nothing
            fn, ftype, line, sig = m.captures
            linenr = parse(Int, line)

            # fix end linenumber of previous def
            if previous != nothing
                functions[previous] = tuple(functions[previous][1], linenr-1)
            end

            if ftype == "function"
                functions[fn] = tuple(linenr, 0)
                previous = fn
            end
        end
    end

    # create function nodes and save range info
    nodes = Dict{String, Node}()
    for (fn,range) in functions
        nodes[fn] = Node{FunctionNode}()
        nodes[fn][:range] = range
    end

    # extract code snippets (extremely inefficient, I know)
    snippets = Dict{String,String}()
    for (id, f) in nodes
        code = open(readlines, path)
        range = f[:range]

        if range[2] == 0
            # ctags didn't manage to find the end, so do a quick scan for '}'
            endline = lines # worst case

            snippet = code[range[1]:end]
            for i in 1:length(snippet)
                if rstrip(snippet[i]) == "}"
                    snippet = snippet[1:i]
                    endline = range[1]+i
                    break
                end
            end

            range = tuple(range[1], min(lines, endline))
        end

        f[:range] = range
        f[:code] = join(code[range[1]:range[2]])
    end

    # gather last changed info
    # times = Vector{Int}()
    # for line in readlines(`git -C $repo blame $path --line-porcelain`)
    #     time_rt = r"^committer-time (\d+)$"
    #     m = match(time_rt, line)
    #     if m != nothing
    #         push!(times, parse(Int, m.captures[1]))
    #     end
    # end
    # @assert lines == length(times)
    # for (id, f) in nodes
    #     range = f[:range]
    #     f[:mtime] = maximum(times[range[1]:range[2]])
    # end

    for (id, f) in nodes
        parent.nodes[id] = f
    end
    return nodes
end

function create_file!(root::Node{DirectoryNode}, path::String)
    # make sure all branches towards the file exist
    parent, id = create_path!(root, path)

    # gather some properties
    size = stat(path).size
    lines = length(open(readlines, path))

    node = Node{FileNode}()
    node[:path] = path
    node[:size] = size
    node[:lines] = lines

    parent.nodes[id] = node
    return node
end

function propagate_info!(node::Node{DirectoryNode})
    node[:lines] = 0
    node[:size] = 0

    # process all director nodes
    dirs = filter((k,v) -> isa(v, Node{DirectoryNode}), node.nodes)
    for (id, dir) in dirs
        propagate_info!(dir)

        node[:lines] += dir[:lines]
        node[:size] += dir[:size]
    end

    # process all file nodes
    files = filter((k,v) -> isa(v, Node{FileNode}), node.nodes)
    if length(files) > 0
        node[:lines] += sum(pair->pair[2][:lines], files)
        node[:size] += sum(pair->pair[2][:size], files)
    end
end

function Tree(repo, sources)
    root = Node{DirectoryNode}()

    # add leaves to the tree
    p = Progress(length(sources), 1)
    for path in sources
        rel = relpath(path, repo)
        file = create_file!(root, rel)
        file.nodes = create_functions!(file, repo)
        next!(p)
    end

    # propagate information upwards in a DFS
    propagate_info!(root)

    return root
end

#
# Main
#

import Base: filter!

using JSON

const EXTENSIONS = [".c", ".cxx", ".cpp"]   # TODO: should also parse headers

# Recursively strip a node of certain properties
function filter!(cb, node::Node)
    filter!(cb, node.props)
    for (id,node) in node.nodes
        filter!(cb, node)
    end
end

function main(args)
    if length(args) != 1
        error("Usage: parse.jl SOURCEDIR")
    end
    repo = args[1]
    isdir(repo) || error("$repo isn't a directory")
    sources = scan(repo, EXTENSIONS)

    tree = cd(repo) do
        Tree(repo, sources)
    end

    # write to disk, only the info necessary to visualize the tree
    structure = deepcopy(tree)
    filter!((k,v) -> k!=:code, structure)
    open("sample_structure.json", "w") do io
        JSON.print(io, structure)
    end

    # write to disk, only the code
    code = deepcopy(tree)
    filter!((k,v) -> k==:code, code)
    open("sample_code.json", "w") do io
        JSON.print(io, code)
    end

    # HACK: get all valid function names, for easy matching in parse_callgrind.jl
    file_fns = Dict{String,Vector{String}}()
    function get_fns!(node)
        for (id, child) in node.nodes
            if isa(child, Node{FileNode})
                file_fns[id] = [k for (k,v) in child.nodes if isa(v, Node{FunctionNode})]
            else
                get_fns!(child)
            end
        end
    end
    get_fns!(tree)
    open("functions.dat", "w") do io
        serialize(io, file_fns)
    end
end
main(ARGS)

max_stars_repo_name: OhioAdam/Data-Tree-Modeler

max_issues_repo_name: OhioAdam/Data-Tree-Modeler

max_forks_repo_name: OhioAdam/Data-Tree-Modeler


text:
m = GenericModel()

# model paths. all this should work without errors
m.testing = true
addl_strings = ["foo=bar", "hat=head", "city=newyork"]
@testset "Check proper model paths" begin
    for fn in [:rawpath, :workpath, :tablespath, :figurespath]
        @eval $(fn)(m, "test")
        @eval $(fn)(m, "test", "temp")
        @eval $(fn)(m, "test", "temp", addl_strings)
    end
end

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors