text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PCL"
version = v"1.11.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/PointCloudLibrary/pcl/releases/download/pcl-$version/source.tar.gz",
                  "19d1a0bee2bc153de47c05da54fc6feb23393f306ab2dea2e25419654000336e"),
    DirectorySource("./bundled")
]

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

# Patch to simplify CMake checks
atomic_patch -p1 ../patches/0001-Replace-run-checks-with-compile-checks.patch

if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ../patches/windows-cases.patch
    atomic_patch -p1 ../patches/pcl_io-link-ws2_32.patch
    atomic_patch -p1 ../patches/ssize_t-mingw.patch
fi

mkdir build && cd build

#see https://github.com/PointCloudLibrary/pcl/pull/4695 for -DPCL_WARNINGS_ARE_ERRORS flag

cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DWITH_VTK=OFF \
    -DWITH_LIBUSB=OFF \
    -DWITH_QT=OFF \
    -DWITH_CUDA=OFF \
    -DWITH_QHULL=OFF \
    -DWITH_OPENGL=OFF \
    -DWITH_PCAP=OFF \
    -DPCL_WARNINGS_ARE_ERRORS=OFF \
    -DCMAKE_BUILD_TYPE=Release

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

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [

    ExecutableProduct("pcl_add_gaussian_noise", :pcl_add_gaussian_noise),
    ExecutableProduct("pcl_boundary_estimation", :pcl_boundary_estimation),
    ExecutableProduct("pcl_cluster_extraction", :pcl_cluster_extraction),
    ExecutableProduct("pcl_compute_cloud_error", :pcl_compute_cloud_error),
    ExecutableProduct("pcl_compute_hausdorff", :pcl_compute_hausdorff),
    ExecutableProduct("pcl_concatenate_points_pcd", :pcl_concatenate_points_pcd),
    ExecutableProduct("pcl_convert_pcd_ascii_binary", :pcl_convert_pcd_ascii_binary),
    ExecutableProduct("pcl_crf_segmentation", :pcl_crf_segmentation),
    ExecutableProduct("pcl_demean_cloud", :pcl_demean_cloud),
    ExecutableProduct("pcl_elch", :pcl_elch),
    ExecutableProduct("pcl_extract_feature", :pcl_extract_feature),
    ExecutableProduct("pcl_fast_bilateral_filter", :pcl_fast_bilateral_filter),
    ExecutableProduct("pcl_fpfh_estimation", :pcl_fpfh_estimation),
    ExecutableProduct("pcl_generate", :pcl_generate),
    ExecutableProduct("pcl_gp3_surface", :pcl_gp3_surface),
    ExecutableProduct("pcl_grid_min", :pcl_grid_min),
    ExecutableProduct("pcl_hdl_grabber", :pcl_hdl_grabber),
    ExecutableProduct("pcl_icp", :pcl_icp),
    ExecutableProduct("pcl_icp2d", :pcl_icp2d),
    ExecutableProduct("pcl_linemod_detection", :pcl_linemod_detection),
    ExecutableProduct("pcl_local_max", :pcl_local_max),
    ExecutableProduct("pcl_lum", :pcl_lum),
    ExecutableProduct("pcl_marching_cubes_reconstruction", :pcl_marching_cubes_reconstruction),
    ExecutableProduct("pcl_match_linemod_template", :pcl_match_linemod_template),
    ExecutableProduct("pcl_mls_smoothing", :pcl_mls_smoothing),
    ExecutableProduct("pcl_morph", :pcl_morph),
    ExecutableProduct("pcl_ndt2d", :pcl_ndt2d),
    ExecutableProduct("pcl_ndt3d", :pcl_ndt3d),
    ExecutableProduct("pcl_normal_estimation", :pcl_normal_estimation),
    ExecutableProduct("pcl_outlier_removal", :pcl_outlier_removal),
    ExecutableProduct("pcl_passthrough_filter", :pcl_passthrough_filter),
    ExecutableProduct("pcl_pcd2ply", :pcl_pcd2ply),
    ExecutableProduct("pcl_pcd2vtk", :pcl_pcd2vtk),
    ExecutableProduct("pcl_pcd_change_viewpoint", :pcl_pcd_change_viewpoint),
    ExecutableProduct("pcl_pcd_convert_NaN_nan", :pcl_pcd_convert_NaN_nan),
    ExecutableProduct("pcl_pcd_introduce_nan", :pcl_pcd_introduce_nan),
    ExecutableProduct("pcl_pclzf2pcd", :pcl_pclzf2pcd),
    ExecutableProduct("pcl_plane_projection", :pcl_plane_projection),
    ExecutableProduct("pcl_ply2obj", :pcl_ply2obj),
    ExecutableProduct("pcl_ply2pcd", :pcl_ply2pcd),
    ExecutableProduct("pcl_ply2ply", :pcl_ply2ply),
    ExecutableProduct("pcl_ply2raw", :pcl_ply2raw),
    ExecutableProduct("pcl_plyheader", :pcl_plyheader),
    ExecutableProduct("pcl_poisson_reconstruction", :pcl_poisson_reconstruction),
    ExecutableProduct("pcl_progressive_morphological_filter", :pcl_progressive_morphological_filter),
    ExecutableProduct("pcl_radius_filter", :pcl_radius_filter),
    ExecutableProduct("pcl_sac_segmentation_plane", :pcl_sac_segmentation_plane),
    ExecutableProduct("pcl_spin_estimation", :pcl_spin_estimation),
    ExecutableProduct("pcl_train_linemod_template", :pcl_train_linemod_template),
    ExecutableProduct("pcl_train_unary_classifier", :pcl_train_unary_classifier),
    ExecutableProduct("pcl_transform_from_viewpoint", :pcl_transform_from_viewpoint),
    ExecutableProduct("pcl_transform_point_cloud", :pcl_transform_point_cloud),
    ExecutableProduct("pcl_unary_classifier_segment", :pcl_unary_classifier_segment),
    ExecutableProduct("pcl_uniform_sampling", :pcl_uniform_sampling),
    ExecutableProduct("pcl_vfh_estimation", :pcl_vfh_estimation),
    ExecutableProduct("pcl_voxel_grid", :pcl_voxel_grid),
    ExecutableProduct("pcl_xyz2pcd", :pcl_xyz2pcd),

    LibraryProduct("libpcl_common", :libpcl_common),
    LibraryProduct("libpcl_features", :libpcl_features),
    LibraryProduct("libpcl_filters", :libpcl_filters),
    LibraryProduct("libpcl_io_ply", :libpcl_io_ply),
    LibraryProduct("libpcl_io", :libpcl_io),
    LibraryProduct("libpcl_kdtree", :libpcl_kdtree),
    LibraryProduct("libpcl_keypoints", :libpcl_keypoints),
    LibraryProduct("libpcl_ml", :libpcl_ml),
    LibraryProduct("libpcl_octree", :libpcl_octree),
    LibraryProduct("libpcl_recognition", :libpcl_recognition),
    LibraryProduct("libpcl_registration", :libpcl_registration),
    LibraryProduct("libpcl_sample_consensus", :libpcl_sample_consensus),
    LibraryProduct("libpcl_search", :libpcl_search),
    LibraryProduct("libpcl_segmentation", :libpcl_segmentation),
    LibraryProduct("libpcl_stereo", :libpcl_stereo),
    LibraryProduct("libpcl_surface", :libpcl_surface),
    LibraryProduct("libpcl_tracking", :libpcl_tracking),

]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
    Dependency(PackageSpec(name="FLANN_jll", uuid="48b6455b-4cf5-590d-a543-2d733c79e793"))
    Dependency(PackageSpec(name="boost_jll", uuid="28df3c45-c428-5900-9ff8-a3135698ca75"))
    Dependency(PackageSpec(name="Eigen_jll", uuid="bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
]

# Build the tarballs, and possibly a `build.jl` as well.
#Unsure of actual lower bound on GCC version, PCL uses C++14 standard and their CI on github uses 7 so I went with that.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# Source: https://stackoverflow.com/questions/41636928/julia-find-the-indices-of-all-maxima
# returns max indices
function find_all_max_indices(arr)
    max_positions = Vector{Int}()
    min_val = typemin(eltype(arr))
    for i in eachindex(arr)
        if arr[i] > min_val
            min_val = arr[i]
            empty!(max_positions)
            push!(max_positions, i)
        elseif arr[i] == min_val
            push!(max_positions, i)
        end
    end
    max_positions
end

ext: jl


text:
"""
Set of functions for generating signals.
"""

function random_walk(process_noise,
                    measurement_noise,
                    mean_state0,
                    precision_state0;
                    time_horizon=100)
    "Generate data according to a random walk"

    # Dimensionality
    d = size(process_noise, 1)

    # Preallocate
    y = zeros(d, time_horizon)
    x = zeros(d, time_horizon+1)

    # Initialize state
    x[:, 1] = sqrt(inv(precision_state0))*randn(d) + mean_state0

    for t = 1:time_horizon

        # Evolve state
        x[:, t+1] = sqrt(process_noise)*randn(d) + x[:, t]

        # Observe
        y[:, t] = sqrt(measurement_noise)*randn(d) + x[:, t+1]

    end
    return y, x
end

function LGDS(transition_matrix,
              emission_matrix,
              process_noise,
              measurement_noise,
              prior_params;
              time_horizon=100)
    "Generate data according to a linear Gaussian dynamical system"

    # Dimensionality
    d = size(transition_matrix, 1)

    # Preallocate
    x = zeros(d, time_horizon+1)
    y = zeros(d, time_horizon)

    # Prior parameters
    mean0, var0 = prior_params

    # Initialize state
    x[:, 1] = sqrt(var0)*randn(d) .+ mean0

    for t = 1:time_horizon

        # Evolve state
        x[:, t+1] = sqrt(process_noise)*randn(d) .+ transition_matrix*x[:,t]

        # Observe
        y[:, t] = sqrt(measurement_noise)*randn(d) .+ emission_matrix*x[:,t+1]

    end
    return y, x
end

function NLGDS(transition_function,
               emission_function,
               process_noise,
               measurement_noise,
               state0_params;
               time_horizon=100)
    "Generate data according to a nonlinear Gaussian dynamical system"

    # Dimensionality
    d = size(process_noise, 1)

    # Preallocate
    x = zeros(d, time_horizon+1)
    y = zeros(d, time_horizon)

    # Prior parameters
    mean0, var0 = state0_params

    # Initialize state
    x[:, 1] = sqrt(var0)*randn(d) .+ mean0

    for t = 1:time_horizon

        # Evolve state
        x[:, t+1] = sqrt(process_noise)*randn(d) .+ transition_function(x[:,t])

        # Observe
        y[:, t] = sqrt(measurement_noise)*randn(d) .+ emission_function(x[:,t+1])

    end
    return y, x
end

max_stars_repo_name: wmkouw/kalman-filters

max_issues_repo_name: wmkouw/kalman-filters

max_forks_repo_name: wmkouw/kalman-filters


text:
function (p::Array{Float64, 1}, q::Array{Float64, 1})
    a = (p - q)^2
    b =  p + q
    a[abs(b) < 1e-12] = 0.0
    b[abs(b) < 1e-12] = 1.0

    r = a / b

    return 0.5*sum(r)

max_stars_repo_name: vladpopovici/WSItk

max_issues_repo_name: vladpopovici/WSItk

max_forks_repo_name: vladpopovici/WSItk


text:
# copied from CUDAdrv/deps/build.jl

using Libdl
use_artifacts = !parse(Bool, get(ENV, "JULIA_AMDGPU_DISABLE_ARTIFACTS", "false"))

function version_hsa(libpath)
    lib = Libdl.dlopen(libpath)
    sym = Libdl.dlsym(lib, "hsa_system_get_info")
    major_ref = Ref{Cushort}(typemax(Cushort))
    minor_ref = Ref{Cushort}(typemax(Cushort))
    status = ccall(sym, Cint, (Cint, Ptr{Cushort}), 0, major_ref)
    if status != 0
        @warn "HSA error: $status"
        return v"0"
    end
    status = ccall(sym, Cint, (Cint, Ptr{Cushort}), 1, minor_ref)
    if status != 0
        @warn "HSA error: $status"
        return v"0"
    end
    return VersionNumber(major_ref[], minor_ref[])
end

function init_hsa(libpath)
    lib = Libdl.dlopen(libpath)
    sym = Libdl.dlsym(lib, "hsa_init")
    ccall(sym, Cint, ())
end

function shutdown_hsa(libpath)
    lib = Libdl.dlopen(libpath)
    sym = Libdl.dlsym(lib, "hsa_shut_down")
    ccall(sym, Cint, ())
end

## auxiliary routines

status = 0
function build_warning(reason)
    println("$reason.")
    global status
    status = 1
end

function build_error(reason)
    println(reason)
    exit(1)
end

## library finding

function find_roc_paths()
    paths = split(get(ENV, "LD_LIBRARY_PATH", ""), ":")
    paths = filter(path->path != "", paths)
    paths = map(Base.Filesystem.abspath, paths)
    push!(paths, "/opt/rocm/hsa/lib") # shim for Ubuntu rocm packages...
    paths = filter(isdir, paths)
    @show paths
    return paths
end

function find_hsa_library(lib, dirs)
    path = Libdl.find_library(lib)
    if path != ""
        return Libdl.dlpath(path)
    end
    for dir in dirs
        files = readdir(dir)
        for file in files
            @info "$file: $(basename(file) == lib * ".so.1")"
            if basename(file) == lib * ".so.1"
                return joinpath(dir, file)
            end
        end
    end
end

function find_ld_lld()
    paths = split(get(ENV, "PATH", ""), ":")
    paths = filter(path->path != "", paths)
    paths = map(Base.Filesystem.abspath, paths)
    basedir = get(ENV, "ROCM_PATH", "/opt/rocm")
    ispath(joinpath(basedir, "llvm/bin/ld.lld")) &&
        push!(paths, joinpath(basedir, "llvm/bin/"))
    ispath(joinpath(basedir, "hcc/bin/ld.lld")) &&
        push!(paths, joinpath(basedir, "/hcc/bin/"))
    ispath(joinpath(basedir, "opencl/bin/x86_64/ld.lld")) &&
        push!(paths, joinpath(basedir, "opencl/bin/x86_64/"))
    for path in paths
        exp_ld_path = joinpath(path, "ld.lld")
        if ispath(exp_ld_path)
            try
                tmpfile = mktemp()
                run(pipeline(`$exp_ld_path -v`; stdout=tmpfile[1]))
                vstr = read(tmpfile[1], String)
                rm(tmpfile[1])
                vstr = replace(vstr, "AMD " => "")
                vstr_splits = split(vstr, ' ')
                if VersionNumber(vstr_splits[2]) >= v"6.0.0"
                    @info "Found useable ld.lld at $exp_ld_path"
                    return exp_ld_path
                end
            catch
                @warn "Failed running ld.lld in $exp_ld_path"
            end
        end
    end
    return ""
end

function find_roc_library(name::String)
    lib = Libdl.find_library(Symbol(name))
    lib == "" && return nothing
    return Libdl.dlpath(lib)
end

## main

const config_path = joinpath(@__DIR__, "ext.jl")
const previous_config_path = config_path * ".bak"

function write_ext(config, path)
    open(path, "w") do io
        println(io, "# autogenerated file, do not edit")
        for (key,val) in config
            println(io, "const $key = $(repr(val))")
        end
    end
end

function main()
    ispath(config_path) && mv(config_path, previous_config_path; force=true)
    config = Dict{Symbol,Any}(
        :configured => false,
        :build_reason => "unknown",
        :lld_configured => false,
        :lld_build_reason => "unknown",
        :hsa_configured => false,
        :hsa_build_reason => "unknown",
        :hip_configured => false,
        :hip_build_reason => "unknown",
        :device_libs_configured => false,
        :device_libs_build_reason => "unknown",
        :librocblas => nothing,
        :librocsparse => nothing,
        :librocalution => nothing,
        :librocfft => nothing,
        :librocrand => nothing,
        :libmiopen => nothing,
    )
    write_ext(config, config_path)

    # Skip build if running under AutoMerge
    if get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
        exit(0)
    end

    ## discover stuff

    # check that we're running Linux
    if !Sys.islinux()
        build_warning("Not running Linux, which is the only platform currently supported by the ROCm Runtime.")
        config[:build_reason] = "Unsupported OS: $(repr(Sys.KERNEL))"
        write_ext(config, config_path)
        return
    end

    # find some paths for library search
    roc_dirs = find_roc_paths()

    ### Find HSA
    libhsaruntime_path = nothing
    if use_artifacts
        try
            @eval using hsa_rocr_jll
        catch err
            iob = IOBuffer()
            println(iob, "`using hsa_rocr_jll` failed:")
            Base.showerror(iob, err)
            Base.show_backtrace(iob, catch_backtrace())
            config[:hsa_build_reason] = String(take!(iob))
            write_ext(config, config_path)
            return
        end
        libhsaruntime_path = hsa_rocr_jll.libhsa_runtime64
    else
        libhsaruntime_path = find_hsa_library("libhsa-runtime64.so.1", roc_dirs)
    end
    if libhsaruntime_path === nothing
        build_warning("Could not find HSA runtime library v1")
        config[:hsa_build_reason] = "HSA runtime library v1 not found"
        write_ext(config, config_path)
        return
    end

    # initializing the library isn't necessary, but flushes out errors that otherwise would
    # happen during `version` or, worse, at package load time.
    status = init_hsa(libhsaruntime_path)
    if status != 0
        build_warning("Initializing HSA runtime failed with code $status.")
        config[:hsa_build_reason] = "Failed to initialize HSA runtime, status code: $status"
        write_ext(config, config_path)
        return
    end

    libhsaruntime_version = version_hsa(libhsaruntime_path)

    # also shutdown just in case
    status = shutdown_hsa(libhsaruntime_path)
    if status != 0
        build_warning("Shutdown of HSA runtime failed with code $status.")
        config[:hsa_build_reason] = "Failed to shutdown HSA runtime, status code: $status"
        write_ext(config, config_path)
        return
    end
    config[:libhsaruntime_path] = libhsaruntime_path
    config[:libhsaruntime_version] = libhsaruntime_version
    config[:hsa_configured] = true

    ### Find HIP
    libhip_path = nothing
    if use_artifacts
        try
            @eval using HIP_jll
        catch err
            iob = IOBuffer()
            println(iob, "`using HIP_jll` failed:")
            Base.showerror(iob, err)
            Base.show_backtrace(iob, catch_backtrace())
            config[:hip_build_reason] = String(take!(iob))
            write_ext(config, config_path)
            return
        end
        libhip_path = HIP_jll.libamdhip64
    else
        libhip_path = Libdl.find_library(["libamdhip64", "libhip_hcc"])
    end
    if libhip_path === nothing
        build_warning("Could not find HIP runtime library")
        config[:hip_build_reason] = "HIP runtime library not found"
        write_ext(config, config_path)
        return
    end
    config[:libhip_path] = libhip_path
    config[:hip_configured] = true

    ### Find ld.lld
    ld_path = find_ld_lld()
    if ld_path == ""
        build_warning("Could not find ld.lld, please install it with your package manager")
        config[:lld_build_reason] = "ld.lld executable not found"
        write_ext(config, config_path)
        return
    end
    config[:ld_lld_path] = ld_path
    config[:lld_configured] = true

    ### Find/download device-libs
    device_libs_path = nothing
    device_libs_downloaded = nothing
    if use_artifacts
        try
            @eval using ROCmDeviceLibs_jll
        catch err
            iob = IOBuffer()
            println(iob, "`using ROCmDeviceLibs_jll` failed:")
            Base.showerror(iob, err)
            Base.show_backtrace(iob, catch_backtrace())
            config[:device_libs_build_reason] = String(take!(iob))
            write_ext(config, config_path)
            return
        end
        device_libs_path = ROCmDeviceLibs_jll.bitcode_path
        device_libs_downloaded = false
    else
        #include("download_device_libs.jl")
        device_libs_path = "/opt/rocm/amdgcn/bitcode"
        device_libs_downloaded = true
    end
    config[:device_libs_path] = device_libs_path
    config[:device_libs_downloaded] = device_libs_downloaded
    config[:device_libs_configured] = true

    ### Find external HIP-based libraries
    for name in ("rocblas", "rocsparse", "rocalution", "rocfft", "rocrand", "MIOpen")
        lib = Symbol("lib$(lowercase(name))")
        config[lib] = find_roc_library("lib$name")
        if config[lib] === nothing
            build_warning("Could not find library '$name'")
            # TODO: Save build reason?
        end
    end

    config[:configured] = true

    ## (re)generate ext.jl

    function globals(mod)
        all_names = names(mod, all=true)
        filter(name-> !any(name .== [nameof(mod), Symbol("#eval"), :eval]), all_names)
    end

    if isfile(previous_config_path)
        @eval module Previous; include($previous_config_path); end
        previous_config = Dict{Symbol,Any}(name => getfield(Previous, name)
                                           for name in globals(Previous))

        if config == previous_config
            mv(previous_config_path, config_path; force=true)
            return
        end
    end

    write_ext(config, config_path)

    if status != 0
        # we got here, so the status is non-fatal
        build_warning("""

            AMDGPU.jl has been built successfully, but there were warnings.
            Some functionality may be unavailable.""")
    end
end

# Load HSA, HIP, and friends, and ROCm external libraries
main()


text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

include("resource_nonlinear_20_30.jl")

const gurobi_env = Gurobi.Env()

const numVariables = numResources

# generate random samples of the demands
function generateRandomSamples(numSamples::Int64)

	lambda = rand(MvNormal(lambda_mean[1,:],lambda_covariance),numSamples)
	rho = Array{Float64}(numResources,numSamples)
	for i = 1:numResources
		rho[i,:] = rand(Normal(rho_mean[i],rho_stdev),numSamples)
		for samp = 1:numSamples
			rho[i,samp] = min(rho[i,samp],rho_max)
		end
	end

	return lambda, rho
end

srand(1234)

# generate samples to estimate actual risk level of solution
const numAnalyticalSamples = 100000
lambda_analytical, rho_analytical = generateRandomSamples(numAnalyticalSamples)

srand()

# compute conservative bound on the probability that
# constraints aren't satisfied
function boundRiskLevel(x::Array{Float64},reliabilityLevel::Float64=1E-06)

	const numSamples = size(lambda_analytical,2)
	simpleEst::Bool = false

	numViolated::Int64 = checkScenarioConstraints(x,lambda_analytical,rho_analytical)

	riskLevel::Float64 = Inf

	if(simpleEst)
		riskLevel = numViolated*1.0/numSamples
		return riskLevel
	else
		if(numViolated == 0)
			riskLevel = -log(reliabilityLevel)/numSamples
			return riskLevel
		end

		 # termination criterion for bisection
		gamma_tolerance::Float64 = 0.1/numSamples

		# use tail bounds on binomial distribution to estimate lower and upper bounds
		gamma_low::Float64 = (numViolated - sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_low < 0.0)
			gamma_low = 0.0
		end
		gamma_up::Float64 = (numViolated + sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_up > 1.0)
			gamma_up = 1.0
		end

		gamma::Float64 = 0.0
		func_value::Float64 = 0.0

		# use bisection to solve the nonlinear equation
		while (gamma_up - gamma_low > gamma_tolerance)

			gamma = (gamma_low + gamma_up)/2.0
			func_value = -reliabilityLevel
			for i = 1:numViolated
				tmp_func = i*log(gamma) + (numSamples-i)*log(1-gamma) + lfact(numSamples) - lfact(numSamples-i) - lfact(i)
				func_value += exp(tmp_func)
				if(func_value > 0)
					break
				end
			end
			if(func_value > 0)
				gamma_low = (gamma_low + gamma_up)/2.0
			else
				gamma_up = (gamma_low + gamma_up)/2.0
			end

		end

		riskLevel = gamma_up
		return riskLevel
	end

end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function checkScenarioConstraints(x::Array{Float64},lambda::Array{Float64},rho::Array{Float64})

	const numSamples = size(lambda,2)
	numViolated::Int64 = 0

	mod = Model(solver=GurobiSolver(gurobi_env,Presolve=0,OutputFlag=0))

	@variable(mod, y[1:numResources,1:numCustomers] >= 0)

	@objective(mod, Min, 0)
	@constraint(mod, con1[i=1:numResources], sum(y[i,j] for j = 1:numCustomers) <= 0)
	@constraint(mod, con2[j=1:numCustomers], sum(nu[i,j]*y[i,j] for i = 1:numResources) >= 0)

	for iter = 1:numSamples
		for i = 1:numResources
			JuMP.setRHS(con1[i],rho[i,iter]*(x[i])^2)
		end

		for j = 1:numCustomers
			JuMP.setRHS(con2[j],lambda[j,iter])
		end

		status = solve(mod,suppress_warnings=true)

		if(status != Symbol("Optimal"))
			numViolated += 1
		end
	end

	mod = 0

	return numViolated
end

numTrials = 10

# solution file name
const baseDirName = "C:/Users/rkannan/Desktop/SA for CCP/Stochastic approximation/experiments/resource_stochapprox/"

for trial = 1:numTrials

	const dirName = baseDirName * "rep" * string(trial) * "/"

	const fileName = dirName * "solution.txt"

	const outputFile = dirName * "trueRiskLevels.txt"
	open(outputFile, "w") do f
	end

	const timeFile = dirName * "trueRiskTime.txt"
	open(timeFile, "w") do f
	end

	soln = readdlm(fileName)

	const numPoints = size(soln,1)

	for iter = 1:numPoints
		tic()
		trueRiskLevel::Float64 = boundRiskLevel(soln[iter,:])
		riskTime = toq()
		println("Risk level #",iter," : ",trueRiskLevel,"  time: ",riskTime)
		open(outputFile, "a") do f
			write(f,"$trueRiskLevel \n")
		end
		open(timeFile, "a") do f
			write(f,"$riskTime \n")
		end
		gc()
	end
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
@testset "Charts tests:" begin

    umin   = 0.0
    umax   = 2.0
    unodes = 16

    ucoord = GaussLobatto("u", umin, umax, unodes)
    u, = Jecco.cheb(umin, umax, unodes)
    @test ucoord[:] ≈ u

    xmin   = -10.0
    xmax   =  10.0
    xnodes =  20

    xcoord = Cartesian{2}("x", xmin, xmax, xnodes, endpoint=false)
    @test xcoord[:] == collect(-10.0:1.0:9.0)

    ymin   = -20.0
    ymax   =  20.0
    ynodes =  20

    ycoord = Cartesian{3}("y", ymin, ymax, ynodes, endpoint=false)

    chart = Chart(ucoord, xcoord, ycoord)

    u, x, y = chart[:]

    @test x == collect(-10.0:1.0:9.0)
    @test y == collect(-20.0:2.0:19.0)

end


text:
using JLD2
using GLMakie

include("one_degree_artifacts.jl")
# bathymetry_path = download_bathymetry() # not needed because we uploaded to repo
boundary_conditions_path = download_boundary_conditions()

file = jldopen(bathymetry_path)
h = file["bathymetry"]
land = h .> 0
wet_Nx = sum(h -> h < 0, h, dims=1)
close(file)

ρ₀ = 1025
φ = -74.5:74.5
file = jldopen(boundary_conditions_path)

τˣ = file["τˣ"] ./ ρ₀
τʸ = file["τʸ"] ./ ρ₀
Tˢ = file["Tˢ"]
Sˢ = file["Sˢ"]

close(file)

τˣ[isnan.(τˣ)] .= 0
τʸ[isnan.(τʸ)] .= 0

τlim = 0.8 * max(maximum(abs, τˣ), maximum(abs, τʸ))

for n = 1:12
    view(τˣ, :, :, n)[land] .= NaN
    view(τʸ, :, :, n)[land] .= NaN
    view(Tˢ, :, :, n)[land] .= NaN
    view(Sˢ, :, :, n)[land] .= NaN
end

fig = Figure(resolution=(2000, 900))
axx = Axis(fig[2, 2], xlabel="Longitude", ylabel="Latitude", title="Zonal surface stress")
axy = Axis(fig[2, 3], xlabel="Longitude", ylabel="Latitude", title="Meridonal surface stress")
axT = Axis(fig[3, 2], xlabel="Longitude", ylabel="Latitude")
axS = Axis(fig[3, 3], xlabel="Longitude", ylabel="Latitude")
slider = Slider(fig[4, 2:3], range=1:12, startvalue=1)
n = slider.value

title = @lift string("Boundary conditions in month ", $n)
Label(fig[1, 2:3], title)

τˣn = @lift view(τˣ, :, :, $n)
τʸn = @lift view(τʸ, :, :, $n)
Tn = @lift view(Tˢ, :, :, $n)
Sn = @lift view(Sˢ, :, :, $n)

hmx = heatmap!(axx, τˣn, colorrange=(-τlim, τlim), colormap=:redblue)
hmy = heatmap!(axy, τʸn, colorrange=(-τlim, τlim), colormap=:redblue)
hmT = heatmap!(axT, Tn, colorrange=(-1, 31), colormap=:thermal)
hmS = heatmap!(axS, Sn, colorrange=(25, 38), colormap=:haline)

Colorbar(fig[2, 4], hmx, label="Stress (m² s⁻²)")
Colorbar(fig[3, 1], hmT, label="Surface temperature (ᵒC)", flipaxis=false)
Colorbar(fig[3, 4], hmS, label="Surface salinity (psu)")

display(fig)

record(fig, "one_degree_boundary_conditions.mp4", 1:12, framerate=1) do nn
    n[] = nn
end

#####
##### Zonal means
#####

filternan(t) = ifelse(isnan(t), 0.0, t)
zonal_mean_τˣ = sum(filternan, τˣ, dims=1) ./ wet_Nx
zonal_mean_τʸ = sum(filternan, τʸ, dims=1) ./ wet_Nx
zonal_mean_Tˢ = sum(filternan, Tˢ, dims=1) ./ wet_Nx
zonal_mean_Sˢ = sum(filternan, Sˢ, dims=1) ./ wet_Nx

zonal_mean_τˣ = dropdims(zonal_mean_τˣ, dims=1)'
zonal_mean_τʸ = dropdims(zonal_mean_τʸ, dims=1)'
zonal_mean_Tˢ = dropdims(zonal_mean_Tˢ, dims=1)'
zonal_mean_Sˢ = dropdims(zonal_mean_Sˢ, dims=1)'

fig = Figure(resolution=(2000, 900))

axx = Axis(fig[1, 2], xlabel="Month", ylabel="Latitude", title="Zonally-averaged zonal surface stress")
axy = Axis(fig[1, 3], xlabel="Month", ylabel="Latitude", title="Zonally-averaged meridonal surface stress")
axT = Axis(fig[2, 2], xlabel="Month", ylabel="Latitude")
axS = Axis(fig[2, 3], xlabel="Month", ylabel="Latitude")

τlim = 1.5e-4
hmx = heatmap!(axx, zonal_mean_τˣ, colorrange=(-τlim, τlim), colormap=:redblue)
hmy = heatmap!(axy, zonal_mean_τʸ, colorrange=(-τlim, τlim), colormap=:redblue)
hmT = heatmap!(axT, zonal_mean_Tˢ, colorrange=(-1, 31), colormap=:thermal)
hmS = heatmap!(axS, zonal_mean_Sˢ, colorrange=(25, 38), colormap=:haline)

Colorbar(fig[1, 4], hmx, label="Stress (m² s⁻²)")
Colorbar(fig[2, 1], hmT, label="Surface temperature (ᵒC)", flipaxis=false)
Colorbar(fig[2, 4], hmS, label="Surface salinity (psu)")

display(fig)

save("zonally_averaged_one_degree_boundary_conditions.png", fig)


text:
@testset "user_functions" begin

    @testset "generate_caos_rules" begin
        tree, character_labels, taxa_labels = generate_caos_rules("data/S10593.nex", "data/output")
        @test typeof(tree) == CAOS.Node
    end

    @testset "load_tree" begin
        tree, character_labels, taxa_labels = load_tree("data/output")
        @test typeof(tree) == CAOS.Node
    end

    @testset "classify_new_sequence" begin
        tree, character_labels, taxa_labels = load_tree("data/output")
        classification = classify_new_sequence(tree, character_labels, taxa_labels, "data/dna_seq.txt", "data/output")
        @test typeof(classification) == String
    end

end


text:
"""
$(TYPEDEF)
An abstract type that contains delay stochastic simulation algorithms:
- DelayDirect
- DelayRejection
- DelayMNRM
- DelayDirectCR
"""
abstract type AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)
Delay Rejection Method from Barrio, Manuel, Kevin Burrage, André Leier, and Tianhai Tian. "Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation." PLoS computational biology 2, no. 9 (2006): e117..
"""
struct DelayRejection <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)
Delay Direct Method from  Xiaodong Cai, "Exact stochastic simulation of coupled chemical reactions with delays", The Journal of Chemical Physics 126, 124108(2007).
"""
struct DelayDirect <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)
A modifed Composition-Rejection Direct Method with delays (DelayDirectCR), implementation combining features from the original article and from the code in  `DiffEqJump` package: DirectCR :
*A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks*,
by A. Slepoy, A.P. Thompson and S.J. Plimpton, J. Chem. Phys, 128, 205101 (2008).
and
*Efficient Formulations for Exact Stochastic Simulation of Chemical Systems*,
by S. Mauch and M. Stalzer, ACM Trans. Comp. Biol. and Bioinf., 8, No. 1, 27-35 (2010).
"""
struct DelayDirectCR <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)
A modifed version of the Delay Next Reaction Method from
David F. Anderson, "A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays", The Journal of Chemical Physics 128, 109903(2008).
"""
struct DelayMNRM <: AbstractDelayAggregatorAlgorithm end

is_spatial(aggregator::AbstractDelayAggregatorAlgorithm) = false
needs_vartojumps_map(aggregator::AbstractDelayAggregatorAlgorithm) = false
needs_depgraph(aggregator::AbstractDelayAggregatorAlgorithm) = false
needs_depgraph(aggregator::DelayMNRM) = true
needs_depgraph(aggregator::DelayDirectCR) = true


text:
include("../../src/PotentialLearning.jl")
using .PotentialLearning

params = get_conf_params("./")

# Get DFT data
dft_train_data, dft_val_data = generate_data("dft", params)

# Get reference data
ref_train_data, ref_val_data = generate_data("ref", params)

# Get potential learning problem (e.g. A β = b)
snap = learning_problem(dft_train_data, ref_train_data, params)

# Solve potential learning problem (e.g. β = A \ b)
learn(snap, params)

# Validate potentials, forces, and stresses
metrics = error_metrics(snap, dft_val_data - ref_val_data, params)
@show metrics


text:
# scripts for displaying R-trees using PlotlyJS
using PlotlyJS
using Printf: @sprintf

# webcolor palette for R-tree levels
const LevelsPalette = Dict(
    1 => "#228B22", # Forest Green
    2 => "#DAA520", # Goldenrod
    3 => "#FF6347", # Tomato
    4 => "#B22222", # Firebrick
    5 => "#800080", # Purple
    6 => "#4169E1", # Royal blue
    7 => "#008080",) # Teal

# create plotly trace for a single R-tree `node` (rectangle edges)
function node_trace(node::SI.Node{<:Any,2}, ix::Union{Integer, Nothing};
                    showlegend::Bool)
    nbr = SI.mbr(node)
    res = scatter(x = [nbr.low[1], nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1]],
                  y = [nbr.low[2], nbr.low[2], nbr.high[2], nbr.high[2], nbr.low[2]],
                  mode=:lines, line_color=get(LevelsPalette, SI.level(node), "#708090"),
                  line_width=SI.level(node),
                  hoverinfo=:text, hoveron=:points,
                  text="lev=$(SI.level(node)) ix=$(ix !== nothing ? ix : "none")
nchildren=$(length(node)) nelems=$(SI.nelements(node)) area=$(@sprintf("%.2f", SI.area(SI.mbr(node))))",
                  name="level $(SI.level(node))",
                  legendgroup="level $(SI.level(node))", showlegend=showlegend)
    return res
end

# create plotly trace for a single R-tree `node` (cube edges)
function node_trace(node::SI.Node, ix::Union{Integer, Nothing};
                    showlegend::Bool)
    nbr = SI.mbr(node)
    res = scatter3d(x = [nbr.low[1],  nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1],
                         nbr.low[1],  nbr.high[1], nbr.high[1], nbr.low[1], nbr.low[1],
                         NaN, nbr.low[1], nbr.low[1],
                         NaN, nbr.high[1], nbr.high[1],
                         NaN, nbr.high[1], nbr.high[1]],
                    y = [nbr.low[2],  nbr.low[2],  nbr.high[2], nbr.high[2], nbr.low[2],
                         nbr.low[2],  nbr.low[2],  nbr.high[2], nbr.high[2], nbr.low[2],
                         NaN, nbr.high[2], nbr.high[2],
                         NaN, nbr.high[2], nbr.high[2],
                         NaN, nbr.low[2], nbr.low[2]],
                    z = [nbr.low[3],  nbr.low[3],  nbr.low[3],  nbr.low[3], nbr.low[3],
                         nbr.high[3], nbr.high[3], nbr.high[3], nbr.high[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3],
                         NaN, nbr.low[3], nbr.high[3]],
                  mode=:lines, line_color=get(LevelsPalette, SI.level(node), "#708090"),
                  line_width=SI.level(node),
                  hoverinfo=:text, hoveron=:points,
                  text="lev=$(SI.level(node)) ix=$(ix !== nothing ? ix : "none")
nchildren=$(length(node)) nelems=$(SI.nelements(node)) area=$(@sprintf("%.2f", SI.area(SI.mbr(node))))",
                  name="level $(SI.level(node))",
                  legendgroup="level $(SI.level(node))", showlegend=showlegend)
    return res
end

# create plotly traces for the nodes in a subtree with the `node` root
# and append them to `node_traces`
function append_subtree_traces!(node_traces::Vector{PlotlyBase.AbstractTrace},
                                node::SI.Node, ix::Union{Integer, Nothing},
                                levels::Set{Int})
    push!(node_traces, node_trace(node, ix, showlegend = SI.level(node) ∉ levels))
    push!(levels, SI.level(node)) # show each level once in legend
    if node isa SI.Branch
        for (i, child) in enumerate(SI.children(node))
            append_subtree_traces!(node_traces, child, i, levels)
        end
    end
    return nothing
end

data_trace(tree::RTree{<:Any, 2}) =
    scatter(mode=:markers, name=:data, marker_color = "#333333", marker_size = 2,
            x=[SI.center(SI.mbr(x)).coord[1] for x in tree],
            y=[SI.center(SI.mbr(x)).coord[2] for x in tree],
            text=["id=$(SI.id(x))" for x in tree], hoverinfo=:text)

data_trace(tree::RTree) =
    scatter3d(mode=:markers, name=:data, marker_color = "#333333", marker_size = 2,
              x=[SI.center(SI.mbr(x)).coord[1] for x in tree],
              y=[SI.center(SI.mbr(x)).coord[2] for x in tree],
              z=[SI.center(SI.mbr(x)).coord[3] for x in tree],
              text=["id=$(SI.id(x))" for x in tree], hoverinfo=:text)

# create Plotly plot of the given tree
function PlotlyJS.plot(tree::RTree)
    ndims(tree) == 1 && throw(ArgumentError("1-D R-trees not supported"))
    ndims(tree) > 3 && @warn("Only 1st-3rd dimensions would be shown for $(ndims(tree))-D trees")
    node_traces = Vector{PlotlyBase.AbstractTrace}()
    append_subtree_traces!(node_traces, tree.root, nothing, Set{Int}())
    PlotlyJS.plot([node_traces; [data_trace(tree)]], Layout(hovermode=:closest))
end


text:
immutable HessenbergQR{T}
    cosθ::Vector{T}
    sinθ::Vector{T}
end

## QR decomposition on an upper Hessenberg matrix
## H will be overwritten by Q'HQ=RQ
function hessenqr!{T}(H::Matrix{T})
    ## Size of matrix
    n = size(H, 1)
    if n != size(H, 2)
        error("matrix must be square")
    end

    ## Rotation factors
    ## Each pair (cosθ[i], sinθ[i]) forms a rotation matrix
    ## Gi = [cosθ[i]  -sinθ[i]]
    ##      [sinθ[i]   cosθ[i]]
    cosθ = ones(T, n)
    sinθ = zeros(T, n)

    for i = 1:(n - 1)
        ## Make sure H is upper Hessenberg
        ## Zero the elements below H[i+1, i]
        H[(i + 2):end, i] = zero(T)
        ## Calculate cosθ and sinθ
        x = H[i, i]
        y = H[i + 1, i]
        r = hypot(x, y)
        ## If r is too small, (cosθ, sinθ) stores the original values (1, 0)
        if r < eps(T)
            r = zero(T)
        else
            cosθ[i] = x / r
            sinθ[i] = -y / r
        end
        ## Apply the rotation on the left H -> Gi * H
        ## H[i, :]     <- cosθ[i] * H[i, :] - sinθ[i] * H[i +1, :]
        ## H[i + 1, :] <- sinθ[i] * H[i, :] + cosθ[i] * H[i +1, :]
        H[i, i] = r
        H[i + 1, i] = zero(T)
        c = cosθ[i]
        s = sinθ[i]
        for j = (i + 1):n
            tmp = H[i, j]
            H[i, j]     = c * tmp - s * H[i + 1, j]
            H[i + 1, j] = s * tmp + c * H[i + 1, j]
        end
    end

    ## Apply the rotations on the right H -> H * Gi'
    ## H[:, i]     <- cosθ[i] * H[:, i] - sinθ[i] * H[:, i + 1]
    ## H[:, i + 1] <- sinθ[i] * H[:, i] + cosθ[i] * H[:, i + 1]
    for i = 1:(n - 1)
        for j = 1:(i + 1)
            tmp = H[j, i]
            H[j, i]     = cosθ[i] * tmp - sinθ[i] * H[j, i + 1]
            H[j, i + 1] = sinθ[i] * tmp + cosθ[i] * H[j, i + 1]
        end
    end

    return HessenbergQR(cosθ, sinθ)
end

## QR decomposition on a tridiagonal matrix
## H will be overwritten by Q'HQ=RQ
function tridiagqr!{T}(H::Matrix{T})
    ## Size of matrix
    n = size(H, 1)
    if n != size(H, 2)
        error("matrix must be square")
    end

    ## Force H to be tridiagonal
    maindiag = diag(H)
    subdiag = diag(H, -1)
    fill!(H, zero(T))
    for i = 1:(n - 1)
        H[i, i] = maindiag[i]
        H[i + 1, i] = subdiag[i]
        H[i, i + 1] = subdiag[i]
    end
    H[n, n] = maindiag[n]

    ## Rotation factors
    ## Each pair (cosθ[i], sinθ[i]) forms a rotation matrix
    ## Gi = [cosθ[i]  -sinθ[i]]
    ##      [sinθ[i]   cosθ[i]]
    cosθ = ones(T, n)
    sinθ = zeros(T, n)

    for i = 1:(n - 1)
        ## Calculate cosθ and sinθ
        x = H[i, i]
        y = H[i + 1, i]
        r = hypot(x, y)
        ## If r is too small, (cosθ, sinθ) stores the original values (1, 0)
        if r < eps(T)
            r = zero(T)
        else
            cosθ[i] = x / r
            sinθ[i] = -y / r
        end
        ## Apply the rotation on the left H -> Gi * H
        ## Update H[i, i] and H[i + 1, i]
        H[i, i] = r
        H[i + 1, i] = zero(T)
        ## Update H[i, i + 1] and H[i + 1, i + 1]
        c = cosθ[i]
        s = sinθ[i]
        tmp = H[i, i + 1]
        H[i,     i + 1] = c * tmp - s * H[i + 1, i + 1]
        H[i + 1, i + 1] = s * tmp + c * H[i + 1, i + 1]
        ## Update H[i, i + 2] and H[i + 1, i + 2]
        if i < n - 1
            H[i,     i + 2] = -s * H[i + 1, i + 2]
            H[i + 1, i + 2] *= c
        end
    end

    ## Apply the rotations on the right H -> H * Gi'
    ## H[:, i]     <- cosθ[i] * H[:, i] - sinθ[i] * H[:, i + 1]
    ## H[:, i + 1] <- sinθ[i] * H[:, i] + cosθ[i] * H[:, i + 1]
    for i = 1:(n - 1)
        c = cosθ[i]
        s = sinθ[i]

        tmp = H[i + 1, i]
        H[i, i]         = c * H[i, i] - s * H[i, i + 1]
        H[i + 1, i]     = c * tmp     - s * H[i + 1, i + 1]
        H[i + 1, i + 1] = s * tmp     + c * H[i + 1, i + 1]

        H[i, i + 1] = H[i + 1, i]
        if i < n - 1
            H[i, i + 2] = zero(T)
        end
    end

    return HessenbergQR(cosθ, sinθ)
end

## Apply the QR factorization to the right of a matrix A
## A -> A * Q'
function applyright!{T}(qr::HessenbergQR{T}, A::Matrix{T})
    cosθ = qr.cosθ
    sinθ = qr.sinθ
    n = length(cosθ)
    ## A[:, i]     <- cosθ[i] * A[:, i] - sinθ[i] * A[:, i + 1]
    ## A[:, i + 1] <- sinθ[i] * A[:, i] + cosθ[i] * A[:, i + 1]
    for i = 1:(n - 1)
        for j = 1:size(A, 1)
            tmp = A[j, i]
            A[j, i]     = cosθ[i] * tmp - sinθ[i] * A[j, i + 1]
            A[j, i + 1] = sinθ[i] * tmp + cosθ[i] * A[j, i + 1]
        end
    end
end


text:
@testset "setup_test.jl" begin

    fov = LaserLab.Fov(1.0mm,1.0mm)
    @test fov.a ≈ π * (fov.d/2.)^2
    @test fov.v ≈ fov.a * fov.z
    @test fov.a ≈ π * (fov.d/2.)^2

    luv = LaserLab.CLaser(375.0nm, 1.0μW)
    @test luv.P/μW ≈ 1
    @test luv.λ/nm ≈ 375.0

    obj = LaserLab.Objective("Generic", 0.5, 40.0)
    @test obj.NA ≈ 0.5
    @test obj.M ≈ 40.0

    obj2 = LaserLab.Objective("Generic", 5.0mm, 5.0mm, 40.0)
    @test obj2.NA ≈ 0.5
    @test obj2.M ≈ 40.0
    @test obj2.f/mm ≈ 5.0
    @test obj2.d/mm ≈ 5.0

    epl375 = LaserLab.PulsedLaser(375.0nm, 140.0mW, 1.0MHz, 65.0ps)
    @test epl375.P == epl375.Pk * epl375.w * epl375.f

    lmu_40x_nuv = LaserLab.Objective("LMU-40X-NUV", 5.3mm, 5.1mm, 40.0)
    @test lmu_40x_nuv.NA  ≈ lmu_40x_nuv.d / (2.0*lmu_40x_nuv.f)
    @test LaserLab.lmu_40x_nuv_transmission(0.500μm) ==LaserLab.lmu_40x_nuv_transmission(500nm)

    @test isapprox(LaserLab.photon_energy(500.0nm)/eV, hnu_eV(500.0), rtol=0.01)
    @test uconvert(J, LaserLab.delivered_energy(luv, 1*s)) ≈ 1.0μJ
    @test LaserLab.n_photons(luv) ≈ LaserLab.n_photons(375.0nm, 1.0μW)
    @test LaserLab.n_photons(luv) ≈ LaserLab.n_photons(375.0nm, 1.0μW)
    @test LaserLab.n_photons_int(luv, 1.0s) ≈ LaserLab.n_photons(luv)/Hz
    @test LaserLab.photon_density(375.0nm, 1.0μW, fov.a) ≈ LaserLab.photon_density(luv, fov)
end

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
##
#
# Draw various histograms
#

using ChristensenSet

DATA_FILE = "ones_Float64.dat"

# set up images
images = (
  (RootsImage(-2.0-2.0im, 2.0+2.0im, 2000, 2000), "ones.png", "ones.tex"),
  (RootsImage(0.8-0.2im, 1.2+0.2im, 2000, 2000), "ones_1.png", "ones_1.tex"),
  (RootsImage(0.5-0.15im, 0.8+0.15im, 2000, 2000), "ones_2.png", "ones_2.tex"),
  (RootsImage(-0.2+0.65im, 0.2+0.8im, 5333, 2000), "ones_3.png", "ones_3.tex"),
  (RootsImage(0.3+0.4im, 0.6+0.7im, 2000, 2000), "ones_4.png", "ones_4.tex"),
)

# fill images with roots
fill_images!(collect(map(first, images)), DATA_FILE)

# save images to disk
for img in images
  save_image(img[1], img[2], latex_filename=img[3])
  run(`pdflatex $(img[3])`)
end

@info "Done."


text:
fldr="C:\\temp\\a\\"
cd(fldr)
cd("SomePackage")
using Pkg
Pkg.activate(".")
@show readdir(fldr)
fi=joinpath("src","file.jl")

absfi=abspath(fi)
@show isfile(absfi)
include(absfi) #works

@assert isfile(fi) #passes
@show pwd()
@show isfile(fi) #true
@show fi
include(fi) #fails
using SomePackage


text:
module FDTD1D
export make_advance_fields, YeeMesh1D

struct YeeMesh1D{T<:AbstractFloat, S<:Integer}
    ey::Vector{T}
    ez::Vector{T}
    hy::Vector{T}
    hz::Vector{T}
    N::S
    function YeeMesh1D{T}(N::S) where {T, S}
        new{T,S}(zeros(T, N+1), zeros(T, N+1), zeros(T, N), zeros(T, N), N)
    end
end

struct PML{T<:Integer, S<:AbstractFloat}
    N::T
    σ_max::S
    r₁::Vector{S}
    r₂::Vector{S}
    function PML(N::T, σ_max::S, Δx, Δt) where {T, S}
        σ = [σ_max*(i/2N)^3 for i = 1:2N]
        r₁ = exp.(-Δt.*σ)
        r₂ = (1.0 .- r₁)./(Δx.*σ)
        new{T,S}(N, σ_max, r₁, r₂)
    end
end

function make_advance_fields(f::YeeMesh1D{T,S}, cfl, pulse_shape, Δt, Δx, x_min, pml::PML = PML(10,1e3,Δt,Δx)) where {T,S}
    Nx = f.N

    function generate_fields_x_min!(t)
        f.ey[pml.N+2] -= cfl*pulse_shape.y(t, x_min + Δx)
        f.ez[pml.N+2] -= cfl*pulse_shape.z(t, x_min + Δx)

        f.hz[pml.N+2] -= cfl*pulse_shape.y(t + 0.5*Δt, x_min + 1.5*Δx)
        f.hy[pml.N+2] -= cfl*pulse_shape.z(t + 0.5*Δt, x_min + 1.5*Δx)
    end

    function update_ey!(jy)
        for i = 2:pml.N+1
            f.ey[i] = pml.r₁[1+2*(pml.N-i+1)]*f.ey[i] - pml.r₂[1+2*(pml.N-i+1)]*(f.hz[i] - f.hz[i-1])
        end
        for i = pml.N+2:Nx-pml.N
            f.ey[i] -= cfl*(f.hz[i] - f.hz[i-1])
        end
        for i = Nx-pml.N+1:Nx
            f.ey[i] = pml.r₁[1+2*(i-Nx+pml.N-1)]*f.ey[i] - pml.r₂[1+2*(i-Nx+pml.N-1)]*(f.hz[i] - f.hz[i-1])
        end
        for i = 1:Nx
            f.ey[i] += jy[i]
        end
    end

    function update_ez!(jz)
        for i = 2:pml.N+1
            f.ez[i] = pml.r₁[1+2*(pml.N-i+1)]*f.ez[i] + pml.r₂[1+2*(pml.N-i+1)]*(f.hy[i] - f.hy[i-1])
        end
        for i = pml.N+2:Nx-pml.N
            f.ez[i] += cfl*(f.hy[i] - f.hz[i-1])
        end
        for i = Nx-pml.N+1:Nx
            f.ez[i] = pml.r₁[1+2*(i-Nx+pml.N-1)]*f.ez[i] + pml.r₂[1+2*(i-Nx+pml.N-1)]*(f.hy[i] - f.hy[i-1])
        end
        for i = 1:Nx
            f.ez[i] += jz[i]
        end
    end

    function update_hy!()
        for i = 1:pml.N
            f.hy[i] = pml.r₁[2*(pml.N-i+1)]*f.hy[i] + pml.r₂[2*(pml.N-i+1)]*(f.ez[i+1] - f.ez[i])
        end
        for i = pml.N+1:Nx-pml.N
            f.hy[i] += cfl*(f.ez[i+1] - f.ez[i])
        end
        for i = Nx-pml.N+1:Nx
            f.hy[i] = pml.r₁[2*(i-Nx+pml.N)]*f.hy[i] + pml.r₂[2*(i-Nx+pml.N)]*(f.ez[i+1] - f.ez[i])
        end
    end

    function update_hz!()
        for i = 1:pml.N
            f.hz[i] = pml.r₁[2*(pml.N-i+1)]*f.hz[i] - pml.r₂[2*(pml.N-i+1)]*(f.ey[i+1] - f.ey[i])
        end
        for i = pml.N+1:Nx-pml.N
            f.hz[i] -= cfl*(f.ey[i+1] - f.ey[i])
        end
        for i = Nx-pml.N+1:Nx
            f.hz[i] = pml.r₁[2*(i-Nx+pml.N)]*f.hz[i] - pml.r₂[2*(i-Nx+pml.N)]*(f.ey[i+1] - f.ey[i])
        end
    end

    function make_step!(j)
        update_ey!(j.y)
        update_ez!(j.z)

        update_hz!()
        update_hy!()
    end

    function advance_fields!(t, j)
        generate_fields_x_min!(t)
        make_step!(j)
    end

    return advance_fields!
end

end


text:
# module Bukdu.Octo

import ..Bukdu

export Assoc

import Base: ==

type Assoc
    vector::Vector{Tuple{Symbol,Any}}
    Assoc(; kw...) = new(Vector(kw))
    Assoc(vector) = new(vector)
    function Assoc(dict::Dict{AbstractString,AbstractString})
        new([(Symbol(k),v) for (k,v) in dict])
    end
end

function Base.getindex(assoc::Assoc, key::Symbol)
    for (k,v) in assoc.vector
        k==key && return v
    end
    throw(KeyError("key $key not found"))
end

function Base.getindex(assoc::Assoc, key::String)
    getindex(assoc, Symbol(key))
end

Base.push!(assoc::Assoc, tup::Tuple{Symbol,Any}) = push!(assoc.vector, tup)

Base.keys(assoc::Assoc) = map(first, assoc.vector)
Base.values(assoc::Assoc) = map(last, assoc.vector)
Base.haskey(assoc::Assoc, key::Symbol) = key in keys(assoc)
Base.haskey(assoc::Assoc, key::String) = haskey(assoc, Symbol(key))

Base.length(assoc::Assoc) = length(assoc.vector)
Base.start(assoc::Assoc) = start(assoc.vector)
Base.next(assoc::Assoc, i::Int) = next(assoc.vector, i)
Base.done(assoc::Assoc, i::Int) = done(assoc.vector, i)

Base.isempty(assoc::Assoc) = isempty(assoc.vector)
Base.empty!(assoc::Assoc) = empty!(assoc.vector)

function Base.setindex!(assoc::Assoc, value::Any, key::Symbol)
    if haskey(assoc, key)
        assoc.vector = map((k,v)->k==key ? value : v, assoc.vector)
    else
        push!(assoc, (key,value))
    end
end

function Base.merge!(lhs::Assoc, rhs::Assoc)
    lkeys = keys(lhs)
    for (rk,rv) in rhs
        if rk in lkeys
            lhs.vector[findn(lkeys,rk)] = rv
        else
            push!(lhs.vector, (rk,rv))
        end
    end
    lhs
end

function ==(lhs::Assoc, rhs::Assoc)
    ==(lhs.vector, rhs.vector)
end

function Base.setdiff(lhs::Assoc, rhs::Assoc)::Assoc
    Assoc(setdiff(lhs.vector, rhs.vector))
end


text:
export ITensor,
       norm,
       delta,
       dims,
       δ,
       replaceindex!,
       inds,
       isNull,
       normalize!,
       multSiteOps,
       order,
       permute,
       randomITensor,
       scalar,
       store

mutable struct ITensor
  inds::IndexSet
  store::TensorStorage
  #TODO: check that the storage is consistent with the
  #total dimension of the indices (possibly only in debug mode);
  ITensor(is::IndexSet,st::T) where T = new(is,st)
end

ITensor() = ITensor(IndexSet(),Dense{Nothing}())
ITensor(is::IndexSet) = ITensor(Float64,is...)
ITensor(inds::Index...) = ITensor(IndexSet(inds...))

function ITensor(::Type{T},
                 inds::IndexSet) where {T<:Number}
  return ITensor(inds,Dense{float(T)}(zero(float(T)),dim(inds)))
end
ITensor(::Type{T},inds::Index...) where {T<:Number} = ITensor(T,IndexSet(inds...))

function ITensor(::UndefInitializer,
                 inds::IndexSet)
  return ITensor(inds,Dense{Float64}(Vector{Float64}(undef,dim(inds))))
end
ITensor(x::UndefInitializer,inds::Index...) = ITensor(x,IndexSet(inds...))

function ITensor(x::S,inds::IndexSet) where {S<:Number}
  return ITensor(inds,Dense{float(S)}(float(x),dim(inds)))
end
ITensor(x::S,inds::Index...) where {S<:Number} = ITensor(x,IndexSet(inds...))

#TODO: check that the size of the Array matches the Index dimensions
function ITensor(A::Array{S},inds::IndexSet) where {S<:Number}
    length(A) ≠ dim(inds) && throw(DimensionMismatch("In ITensor(Array,IndexSet), length of Array ($(length(A))) must match total dimension of IndexSet ($(dim(inds)))"))
  return ITensor(inds,Dense{float(S)}(float(vec(A))))
end
ITensor(A::Array{S},inds::Index...) where {S<:Number} = ITensor(A,IndexSet(inds...))

# Convert to complex
complex(T::ITensor) = ITensor(inds(T),storage_complex(store(T)))

inds(T::ITensor) = T.inds

# This constructor allows many IndexSet
# set operations to work with ITensors
IndexSet(T::ITensor) = inds(T)

store(T::ITensor) = T.store

eltype(T::ITensor) = eltype(store(T))

order(T::ITensor) = order(inds(T))

dims(T::ITensor) = dims(inds(T))

dim(T::ITensor) = dim(inds(T))

isNull(T::ITensor) = (store(T) isa Dense{Nothing})

copy(T::ITensor) = ITensor(copy(inds(T)),copy(store(T)))

Array(T::ITensor) = storage_convert(Array,store(T),inds(T))

function getindex(T::ITensor,vals::Int...)
  if order(T) ≠ length(vals)
    error("In getindex(::ITensor,::Int..), number of \\
           values provided ($(length(vals))) must equal \\
           order of ITensor ($(order(T)))")
  end
  storage_getindex(store(T),inds(T),vals...)
end

function getindex(T::ITensor,ivs::IndexVal...)
  p = calculate_permutation(inds(T),ivs)
  vals = val.(ivs)[p]
  return getindex(T,vals...)
end

function getindex(T::ITensor,ivs::Union{IndexVal, AbstractVector{IndexVal}}...)
  p = calculate_permutation(inds(T),map(x->x isa IndexVal ? x : x[1], ivs))
  vals = map(x->x isa IndexVal ? val(x) : val.(x), ivs[p])
  return storage_getindex(store(T),inds(T),vals...)
end

getindex(T::ITensor) = scalar(T)

setindex!(T::ITensor,x::Number,vals::Int...) = storage_setindex!(store(T),inds(T),x,vals...)

function setindex!(T::ITensor,x::Number,ivs::IndexVal...)
  p = calculate_permutation(inds(T),ivs)
  vals = val.(ivs)[p]
  return setindex!(T,x,vals...)
end

function setindex!(T::ITensor,
                   x::Union{<:Number, AbstractArray{<:Number}},
                   ivs::Union{IndexVal, AbstractVector{IndexVal}}...)
  remap_ivs = map(x->x isa IndexVal ? x : x[1], ivs)
  p = calculate_permutation(inds(T),remap_ivs)
  vals = map(x->x isa IndexVal ? val(x) : val.(x), ivs[p])
  storage_setindex!(store(T),inds(T),x,vals...)
  return T
end

function fill!(T::ITensor,
               x::Number)
  # TODO: automatically switch storage type if needed
  storage_fill!(store(T),x)
  return T
end

function replaceindex!(A::ITensor,i::Index,j::Index)
  pos = indexpositions(A,i)
  inds(A)[pos[1]] = j
  return A
end

prime(A::ITensor,vargs...)= ITensor(prime(inds(A),vargs...),store(A))
prime!(A::ITensor,vargs...)= ( prime!(inds(A),vargs...); return A )

adjoint(A::ITensor) = prime(A)

setprime(A::ITensor,vargs...) = ITensor(setprime(inds(A),vargs...),store(A))

noprime(A::ITensor,vargs...) = ITensor(noprime(inds(A),vargs...),store(A))

mapprime(A::ITensor,vargs...) = ITensor(mapprime(inds(A),vargs...),store(A))

swapprime(A::ITensor,vargs...) = ITensor(swapprime(inds(A),vargs...),store(A))

addtags(A::ITensor,vargs...) = ITensor(addtags(inds(A),vargs...),store(A))

removetags(A::ITensor,vargs...) = ITensor(removetags(inds(A),vargs...),store(A))

replacetags(A::ITensor,vargs...) = ITensor(replacetags(inds(A),vargs...),store(A))

settags(A::ITensor,vargs...) = ITensor(settags(inds(A),vargs...),store(A))

swaptags(A::ITensor,vargs...) = ITensor(swaptags(inds(A),vargs...),store(A))

function ==(A::ITensor,B::ITensor)
  !hassameinds(A,B) && return false
  p = calculate_permutation(inds(B),inds(A))
  for i ∈ CartesianIndices(dims(A))
    A[Tuple(i)...] ≠ B[Tuple(i)[p]...] && return false
  end
  return true
end

function isapprox(A::ITensor,
                  B::ITensor;
                  atol::Real=0.0,
                  rtol::Real=Base.rtoldefault(eltype(A),eltype(B),atol))
    return norm(A-B) <= atol + rtol*max(norm(A),norm(B))
end

function scalar(T::ITensor)
  !(order(T)==0 || dim(T)==1) && throw(ArgumentError("ITensor with inds $(inds(T)) is not a scalar"))
  return storage_scalar(store(T))
end

function randn!(T::ITensor)
  storage_randn!(store(T))
  return T
end

const Indices = Union{IndexSet,Tuple{Vararg{Index}}}

function randomITensor(::Type{S},inds::Indices) where {S<:Number}
  T = ITensor(S,IndexSet(inds))
  randn!(T)
  return T
end
randomITensor(::Type{S},inds::Index...) where {S<:Number} = randomITensor(S,IndexSet(inds...))
randomITensor(inds::Indices) = randomITensor(Float64,IndexSet(inds))
randomITensor(inds::Index...) = randomITensor(Float64,IndexSet(inds...))

norm(T::ITensor) = storage_norm(store(T))
dag(T::ITensor) = ITensor(storage_dag(store(T),inds(T))...)

function permute(T::ITensor,permTinds)
  permTis = IndexSet(permTinds)
  permTstore = typeof(store(T))(dim(T))
  storage_permute!(permTstore,permTis,store(T),inds(T))
  return ITensor(permTis,permTstore)
end
permute(T::ITensor,inds::Index...) = permute(T,IndexSet(inds...))

function *(A::ITensor,x::Number)
    storeB = storage_mult(store(A), x)
    return ITensor(inds(A),storeB)
end
*(x::Number,A::ITensor) = A*x
#TODO: make a proper element-wise division
/(A::ITensor,x::Number) = A*(1.0/x)

-(A::ITensor) = -one(eltype(A))*A
function +(A::ITensor,B::ITensor)
  C = copy(A)
  add!(C,B)
  return C
end
function -(A::ITensor,B::ITensor)
  C = copy(A)
  add!(C,-1,B)
  return C
end

#TODO: Add special case of A==B
#A==B && return ITensor(norm(A)^2)
#TODO: Add more of the contraction logic here?
#We can move the logic of getting the integer labels,
#etc. since they are generic for all storage types
function *(A::ITensor,B::ITensor)
  (Cis,Cstore) = storage_contract(store(A),inds(A),store(B),inds(B))
  C = ITensor(Cis,Cstore)
  return C
end

dot(A::ITensor,B::ITensor) = scalar(dag(A)*B)

#######################################################################
#
# In-place operations
#

"""
normalize!(T::ITensor)

Normalize an ITensor in-place, such that norm(T)==1.
"""
normalize!(T::ITensor) = scale!(T,1/norm(T))

"""
copyto!(B::ITensor, A::ITensor)

Copy the contents of ITensor A into ITensor B.

B .= A
"""
function copyto!(A::ITensor,B::ITensor)
  storage_copyto!(store(A),inds(A),store(B),inds(B))
  return A
end

"""
add!(B::ITensor, A::ITensor)

Add ITensors B and A and store the result in B.

B .+= A
"""
function add!(B::ITensor,A::ITensor)
  B.store = storage_add!(store(B),inds(B),store(A),inds(A))
  return B
end

"""
add!(B::ITensor,α::Number,A::ITensor)

Add ITensors B and α*A and store the result in B.

B .+= α .* A
"""
function add!(A::ITensor,x::Number,B::ITensor)
  A.store = storage_add!(store(A),inds(A),store(B),inds(B),x)
  return A
end

"""
add!(A::ITensor, α::Number, β::Number, B::ITensor)

Add ITensors α*A and β*B and store the result in A.

A .= α .* A .+ β .* B
"""
function add!(A::ITensor,y::Number,x::Number,B::ITensor)
  A.store = storage_add!(y*store(A),inds(A),store(B),inds(B),x)
  return A
end

"""
axpy!(a::Number,v::ITensor,w::ITensor)

w .+= a .* v
"""
axpy!(a::Number,v::ITensor,w::ITensor) = add!(w,a,v)

# This is not implemented correctly
#"""
#w .= a .* v + b .* w
#"""
#axpby!(a::Number,v::ITensor,b::Number,w::ITensor) = add!(w,b,w,a,v)

"""
scale!(A::ITensor,x::Number)

Scale the ITensor A by x in-place.

A .*= x
"""
function scale!(A::ITensor,x::Number)
  storage_mult!(store(A), x)
  return A
end

"""
mul!(A::ITensor,x::Number,B::ITensor)

Multiply ITensor B with x and store the result in A.

A .= x .* B
"""
function mul!(A::ITensor,x::Number,B::ITensor)
  storage_copyto!(store(A),inds(A),store(B),inds(B),x)
  return A
end
mul!(R::ITensor,T::ITensor,fac::Number) = mul!(R,fac,T)

rmul!(T::ITensor,fac::Number) = scale!(T,fac)

#TODO: This is just a stand-in for a proper delta/diag storage type
"""
delta(::Type{T},inds::Index...)

Make a diagonal ITensor with all diagonal elements 1.

WARNING: This is just a stand-in for a proper delta/diag storage type
"""
function delta(::Type{T},inds::Index...) where {T}
  d = ITensor(zero(T),inds...)
  minm = min(dims(d)...)
  for i ∈ 1:minm
    d[IndexVal.(inds,i)...] = one(T)
  end
  return d
end
delta(inds::Index...) = delta(Float64,inds...)
const δ = delta

function show_info(io::IO,
                   T::ITensor)
  print(io,"ITensor ord=$(order(T))")
  for i = 1:order(T)
    print(io," ",inds(T)[i])
  end
  print(io,"\n",typeof(store(T)))
end

function show(io::IO,T::ITensor)
  show_info(io,T)
  print(io,"\n")
  if !isNull(T)
    Base.print_array(io,reshape(data(store(T)),dims(T)))
  end
end

function show(io::IO,
              mime::MIME"text/plain",
              T::ITensor)
  show_info(io,T)
end

function similar(T::ITensor,
                 element_type=eltype(T))::ITensor
  if element_type != eltype(T)
    error("similar(::ITensor) currently only defined for same element type")
  end
  return copy(T)
end

function multSiteOps(A::ITensor,
                     B::ITensor)::ITensor
  R = copy(A)
  prime!(R,"Site")
  R *= B
  return mapprime(R,2,1)
end


text:
using Dates, WeakRefStrings, CategoricalArrays, Tables

@testset "CSV.write" begin

    io = IOBuffer()
    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io)
    @test String(take!(io)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; delim='\t')
    @test String(take!(io)) == "col1\tcol2\tcol3\n1\t4\t7\n2\t5\t8\n3\t6\t9\n"

    (col1=[1,2,3], col2=["hey", "the::re", "::sailor"], col3=[7,8,9]) |> CSV.write(io; delim="::")
    @test String(take!(io)) == "col1::col2::col3\n1::hey::7\n2::\"the::re\"::8\n3::\"::sailor\"::9\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; header=[:Col1, :Col2, :Col3])
    @test String(take!(io)) == "Col1,Col2,Col3\n1,4,7\n2,5,8\n3,6,9\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; header=["Col1", "Col2", "Col3"])
    @test String(take!(io)) == "Col1,Col2,Col3\n1,4,7\n2,5,8\n3,6,9\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; writeheader=false)
    @test String(take!(io)) == "1,4,7\n2,5,8\n3,6,9\n"

    # various types
    data = codeunits("hey there sailor")
    weakrefs = WeakRefStringArray(data, [WeakRefString(pointer(data), 3), WeakRefString(pointer(data), 3), WeakRefString(pointer(data), 3)])
    cats = CategoricalVector{String, UInt32}(["b", "a", "b"])

    (col1=[true, false, true],
     col2=[4.1,5.2,4e10],
     col3=[NaN, Inf, -Inf],
     col4=[Date(2017, 1, 1), Date(2018, 1, 1), Date(2019, 1, 1)],
     col5=[DateTime(2017, 1, 1, 4, 5, 6, 7), DateTime(2018, 1, 1, 4, 5, 6, 7), DateTime(2019, 1, 1, 4, 5, 6, 7)],
     col6=["hey", "there", "sailor"],
     col7=[WeakRefString(pointer(data), 3), WeakRefString(pointer(data), 3), WeakRefString(pointer(data), 3)],
     col8=weakrefs,
     col9=cats,
    ) |> CSV.write(io)
    @test String(take!(io)) == "col1,col2,col3,col4,col5,col6,col7,col8,col9\ntrue,4.1,NaN,2017-01-01,2017-01-01T04:05:06.007,hey,hey,hey,b\nfalse,5.2,Inf,2018-01-01,2018-01-01T04:05:06.007,there,hey,hey,a\ntrue,4.0e10,-Inf,2019-01-01,2019-01-01T04:05:06.007,sailor,hey,hey,b\n"

    (col4=[Date(2017, 1, 1), Date(2018, 1, 1), Date(2019, 1, 1)],
     col5=[DateTime(2017, 1, 1, 4, 5, 6, 7), DateTime(2018, 1, 1, 4, 5, 6, 7), DateTime(2019, 1, 1, 4, 5, 6, 7)],
    ) |> CSV.write(io; dateformat="mm/dd/yyyy")
    @test String(take!(io)) == "col4,col5\n01/01/2017,01/01/2017\n01/01/2018,01/01/2018\n01/01/2019,01/01/2019\n"

    (col1=[1,missing,3], col2=[missing, missing, missing], col3=[7,8,9]) |> CSV.write(io)
    @test String(take!(io)) == "col1,col2,col3\n1,,7\n,,8\n3,,9\n"

    (col1=[1,missing,3], col2=[missing, missing, missing], col3=[7,8,9]) |> CSV.write(io; missingstring="NA")
    @test String(take!(io)) == "col1,col2,col3\n1,NA,7\nNA,NA,8\n3,NA,9\n"

    (col1=["hey, there, sailor", "this, also, has, commas", "this\n has\n newlines\n", "no quoting", "just a random \" quote character", ],) |> CSV.write(io)
    @test String(take!(io)) == "col1\n\"hey, there, sailor\"\n\"this, also, has, commas\"\n\"this\n has\n newlines\n\"\nno quoting\njust a random \" quote character\n"

    (col1=["\"hey there sailor\""],) |> CSV.write(io)
    @test String(take!(io)) == "col1\n\"\\\"hey there sailor\\\"\"\n"

    (col1=["{\"key\": \"value\"}", "{\"key\": null}"],) |> CSV.write(io; openquotechar='{', closequotechar='}')
    @test String(take!(io)) == "col1\n{\\{\"key\": \"value\"\\}}\n{\\{\"key\": null\\}}\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io)
    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; append=false) # this is the default
    @test String(take!(io)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"

    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io)
    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io; append=true)
    @test String(take!(io)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n1,4,7\n2,5,8\n3,6,9\n"

    file = "test.csv"
    (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(file)
    @test String(read(file)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"
    rm(file)

    open(file, "w") do io
        (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]) |> CSV.write(io)
    end
    @test String(read(file)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"
    rm(file)

    # #247
    open(file, "w") do io
        write(io, "or5a2ztZo\n")
        (A=1:3, B=[17, 17, 19], C=["Wg5", "SJ4", "w48"]) |> CSV.write(io; append=true, writeheader=true)
    end
    @test String(read(file)) == "or5a2ztZo\nA,B,C\n1,17,Wg5\n2,17,SJ4\n3,19,w48\n"
    rm(file)

    # unknown schema case
    io = IOBuffer()
    CSV.write(nothing, Tables.rows((col1=[1,2,3], col2=[4,5,6], col3=[7,8,9])), io)
    @test String(take!(io)) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"

    rt = [(a=1, b=4.0, c=7), (a=2.0, b=missing, c="8"), (a=3, b=6.0, c="9")]
    CSV.write(nothing, rt, io)
    @test String(take!(io)) == "a,b,c\n1,4.0,7\n2.0,,8\n3,6.0,9\n"

    CSV.write(nothing, Tables.rows((col1=Int[], col2=Float64[])), io)
    @test String(take!(io)) == ""

    CSV.write(nothing, Tables.rows((col1=Int[], col2=Float64[])), io; header=["col1", "col2"])
    @test String(take!(io)) == "col1,col2\n"

    # 280
    io = IOBuffer()
    CSV.write(io, (x=[',','\n', ','],))
    @test String(take!(io)) == "x\n\",\"\n\"\n\"\n\",\"\n"

    CSV.write(io, (x=['-'],y=['-']), delim='-')
    @test String(take!(io)) == "x-y\n\"-\"-\"-\"\n"

    CSV.write(io, (x= [[1 2; 3 4]],y=[[5 6; 7 8]]), delim=';')
    @test String(take!(io)) == "x;y\n\"[1 2; 3 4]\";\"[5 6; 7 8]\"\n"

    try
        io = open("$file.gz", "w")
        open(`gzip`, "w", io) do f
            CSV.write(f, (col1=[1,2,3], col2=[4,5,6], col3=[7,8,9]))
        end
        run(`gunzip $file.gz`)
        @test String(read("$file")) == "col1,col2,col3\n1,4,7\n2,5,8\n3,6,9\n"
        rm(file)
    catch e
        @error "error running test" exception=(e, stacktrace(catch_backtrace()))
    end
end


text:
using Test, ODBC, DBInterface, Tables, Dates, DecFP, MariaDB_Connector_ODBC_jll, MariaDB_Connector_C_jll

tracefile = abspath(joinpath(@__DIR__, "odbc.log"))
ODBC.setdebug(true, tracefile)
@show ODBC.drivers()
@show ODBC.dsns()
ODBC.setdebug(false)
@test filesize(tracefile) > 0
rm(tracefile)

PLUGIN_DIR = joinpath(MariaDB_Connector_C_jll.artifact_dir, "lib", "mariadb", "plugin")
if Sys.islinux()
    if Int == Int32
        libpath = joinpath(ENV["TRAVIS_BUILD_DIR"], "mariadb32/lib/libmaodbc.so")
    else
        libpath = joinpath(ENV["TRAVIS_BUILD_DIR"], "mariadb64/lib/libmaodbc.so")
    end
elseif Sys.iswindows()
    if Int == Int32
        libpath = joinpath(ENV["TRAVIS_BUILD_DIR"], "mariadb-connector-odbc-3.1.7-win32", "maodbc.dll")
    else
        @show readdir(joinpath(ENV["TRAVIS_BUILD_DIR"], "mariadb-connector-odbc-3.1.7-win64", "SourceDir", "MariaDB", "MariaDB ODBC Driver 64-bit"))
        libpath = joinpath(ENV["TRAVIS_BUILD_DIR"], "mariadb-connector-odbc-3.1.7-win64", "SourceDir", "MariaDB", "MariaDB ODBC Driver 64-bit", "maodbc.dll")
    end
else
    libpath = MariaDB_Connector_ODBC_jll.libmaodbc_path
end
ODBC.adddriver("ODBC_Test_MariaDB", libpath)
ODBC.adddsn("ODBC_Test_DSN_MariaDB", "ODBC_Test_MariaDB"; SERVER="127.0.0.1", UID="root", PLUGIN_DIR=PLUGIN_DIR, Option=67108864, CHARSET="utf8mb4")

conn = DBInterface.connect(ODBC.Connection, "ODBC_Test_DSN_MariaDB")
DBInterface.close!(conn)
conn = DBInterface.connect(ODBC.Connection, "Driver={ODBC_Test_MariaDB};SERVER=127.0.0.1;PLUGIN_DIR=$PLUGIN_DIR;Option=67108864;CHARSET=utf8mb4;USER=root")

DBInterface.execute(conn, "DROP DATABASE if exists mysqltest")
DBInterface.execute(conn, "CREATE DATABASE mysqltest")
DBInterface.execute(conn, "use mysqltest")
DBInterface.execute(conn, """CREATE TABLE Employee
                 (
                     ID INT NOT NULL AUTO_INCREMENT,
                     OfficeNo TINYINT,
                     DeptNo SMALLINT,
                     EmpNo BIGINT UNSIGNED,
                     Wage FLOAT(7,2),
                     Salary DOUBLE,
                     Rate DECIMAL(5, 3),
                     LunchTime TIME,
                     JoinDate DATE,
                     LastLogin DATETIME,
                     LastLogin2 TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP,
                     Initial CHAR(1),
                     Name VARCHAR(255) CHARACTER SET utf8mb4,
                     Photo BLOB,
                     JobType ENUM('HR', 'Management', 'Accounts'),
                     Senior BIT(1),
                     PRIMARY KEY (ID)
                 );""")

DBInterface.execute(conn, """INSERT INTO Employee (OfficeNo, DeptNo, EmpNo, Wage, Salary, Rate, LunchTime, JoinDate, LastLogin, LastLogin2, Initial, Name, Photo, JobType, Senior)
                 VALUES
                 (1, 2, 1301, 3.14, 10000.50, 1.001, '12:00:00', '2015-8-3', '2015-9-5 12:31:30', '2015-9-5 12:31:30', 'A', 'John', 'abc', 'HR', b'1'),
                 (1, 2, 1422, 3.14, 20000.25, 2.002, '13:00:00', '2015-8-4', '2015-10-12 13:12:14', '2015-10-12 13:12:14', 'B', 'Tom', 'def', 'HR', b'1'),
                 (1, 2, 1567, 3.14, 30000.00, 3.003, '12:30:00', '2015-6-2', '2015-9-5 10:05:10', '2015-9-5 10:05:10', 'C', 'Jim', 'ghi', 'Management', b'0'),
                 (1, 2, 3200, 3.14, 15000.50, 2.5, '12:30:00', '2015-7-25', '2015-10-10 12:12:25', '2015-10-10 12:12:25', 'D', '望研測来白制父委供情治当認米注。規', 'jkl', 'Accounts', b'1');
              """)

expected = (
  ID         = Union{Missing, Int32}[1, 2, 3, 4],
  OfficeNo   = Union{Missing, Int8}[1, 1, 1, 1],
  DeptNo     = Union{Missing, Int16}[2, 2, 2, 2],
  EmpNo      = Union{Missing, Int64}[1301, 1422, 1567, 3200],
  Wage       = Union{Missing, Float32}[3.14, 3.14, 3.14, 3.14],
  Salary     = Union{Missing, Float64}[10000.5, 20000.25, 30000.0, 15000.5],
  Rate       = Union{Missing, Dec64}[d64"1.001", d64"2.002", d64"3.003", d64"2.5"],
  LunchTime  = Union{Missing, Dates.Time}[Dates.Time(12,00,00), Dates.Time(13,00,00), Dates.Time(12,30,00), Dates.Time(12,30,00)],
  JoinDate   = Union{Missing, Dates.Date}[Date("2015-08-03"), Date("2015-08-04"), Date("2015-06-02"), Date("2015-07-25")],
  LastLogin  = Union{Missing, Dates.DateTime}[DateTime("2015-09-05T12:31:30"), DateTime("2015-10-12T13:12:14"), DateTime("2015-09-05T10:05:10"), DateTime("2015-10-10T12:12:25")],
  LastLogin2 = Union{Missing, Dates.DateTime}[DateTime("2015-09-05T12:31:30"), DateTime("2015-10-12T13:12:14"), DateTime("2015-09-05T10:05:10"), DateTime("2015-10-10T12:12:25")],
  Initial    = Union{Missing, String}["A", "B", "C", "D"],
  Name       = Union{Missing, String}["John", "Tom", "Jim", "望研測来白制父委供情治当認米注。規"],
  Photo      = Union{Missing, Vector{UInt8}}[b"abc", b"def", b"ghi", b"jkl"],
  JobType    = Union{Missing, String}["HR", "HR", "Management", "Accounts"],
  Senior     = Union{Missing, Bool}[true, true, false, true],
)

cursor = DBInterface.execute(conn, "select * from Employee")
@test eltype(cursor) == ODBC.Row
@test Tables.istable(cursor)
@test Tables.rowaccess(cursor)
@test Tables.rows(cursor) === cursor
@test Tables.schema(cursor) == Tables.Schema(propertynames(expected), eltype.(collect(expected)))
@test Base.IteratorSize(typeof(cursor)) == Base.HasLength()
@test length(cursor) == 4

row = first(cursor)
@test Base.IndexStyle(typeof(row)) == Base.IndexLinear()
@test length(row) == length(expected)
@test propertynames(row) == collect(propertynames(expected))
for (i, prop) in enumerate(propertynames(row))
    @test getproperty(row, prop) == row[prop] == row[i] == expected[prop][1]
end

res = DBInterface.execute(conn, "select * from Employee") |> columntable
@test length(res) == 16
@test length(res[1]) == 4
@test res == expected

# as a prepared statement
stmt = DBInterface.prepare(conn, "select * from Employee")
cursor = DBInterface.execute(stmt)
@test eltype(cursor) == ODBC.Row
@test Tables.istable(cursor)
@test Tables.rowaccess(cursor)
@test Tables.rows(cursor) === cursor
@test Tables.schema(cursor) == Tables.Schema(propertynames(expected), eltype.(collect(expected)))
@test Base.IteratorSize(typeof(cursor)) == Base.HasLength()
@test length(cursor) == 4

row = first(cursor)
@test Base.IndexStyle(typeof(row)) == Base.IndexLinear()
@test length(row) == length(expected)
@test propertynames(row) == collect(propertynames(expected))
for (i, prop) in enumerate(propertynames(row))
    @test getproperty(row, prop) == row[prop] == row[i] == expected[prop][1]
end

res = DBInterface.execute(stmt) |> columntable
@test length(res) == 16
@test length(res[1]) == 4
@test res == expected

@test DBInterface.close!(stmt) === nothing
@test_throws ErrorException DBInterface.execute(stmt)

# insert null row
DBInterface.execute(conn, "INSERT INTO Employee () VALUES ();")
for i = 1:length(expected)
    if i == 1
        push!(expected[i], 5)
    elseif i == 11
    else
        push!(expected[i], missing)
    end
end

res = DBInterface.execute(conn, "select * from Employee") |> columntable
@test length(res) == 16
@test length(res[1]) == 5
for i = 1:length(expected)
    if i != 11
        @test isequal(res[i], expected[i])
    end
end

stmt = DBInterface.prepare(conn, "select * from Employee")
res = DBInterface.execute(stmt) |> columntable
DBInterface.close!(stmt)
@test length(res) == 16
@test length(res[1]) == 5
for i = 1:length(expected)
    if i != 11
        @test isequal(res[i], expected[i])
    end
end

# ODBC.load
ODBC.load(Base.structdiff(expected, NamedTuple{(:LastLogin2, :Wage,)}), conn, "Employee_copy"; limit=4, columnsuffix=Dict(:Name=>"CHARACTER SET utf8mb4"))
res = DBInterface.execute(conn, "select * from Employee_copy") |> columntable
@test length(res) == 14
@test length(res[1]) == 4
for nm in keys(res)
    @test isequal(res[nm], expected[nm][1:4])
end

# now test insert/parameter binding
DBInterface.execute(conn, "DELETE FROM Employee")
for i = 1:length(expected)
    if i != 11
        pop!(expected[i])
    end
end

stmt = DBInterface.prepare(conn,
    "INSERT INTO Employee (OfficeNo, DeptNo, EmpNo, Wage, Salary, Rate, LunchTime, JoinDate, LastLogin, LastLogin2, Initial, Name, Photo, JobType, Senior)
     VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)")

DBInterface.executemany(stmt, Base.structdiff(expected, NamedTuple{(:ID,)}))

stmt2 = DBInterface.prepare(conn, "select * from Employee")
res = DBInterface.execute(stmt2) |> columntable
DBInterface.close!(stmt2)
@test length(res) == 16
@test length(res[1]) == 4
for i = 1:length(expected)
    if i != 11 && i != 1
        @test isequal(res[i], expected[i])
    end
end

DBInterface.execute(stmt, [missing, missing, missing, missing, missing, missing, missing, missing, missing, DateTime("2015-09-05T12:31:30"), missing, missing, missing, missing, missing])
DBInterface.close!(stmt)

stmt = DBInterface.prepare(conn, "select * from Employee")
res = DBInterface.execute(stmt; ignore_driver_row_count=true) |> columntable
DBInterface.close!(stmt)
for i = 1:length(expected)
    if i != 11 && i != 1
        @test res[i][end] === missing
    end
end

DBInterface.execute(conn, """
CREATE PROCEDURE get_employee()
BEGIN
   select * from Employee;
END
""")
res = DBInterface.execute(conn, "call get_employee()") |> columntable
@test length(res) > 0
@test length(res[1]) == 5
res = DBInterface.execute(conn, "call get_employee()") |> columntable
@test length(res) > 0
@test length(res[1]) == 5
# test that we can call multiple stored procedures in a row w/o collecting results (they get cleaned up properly internally)
res = DBInterface.execute(conn, "call get_employee()")
res = DBInterface.execute(conn, "call get_employee()")

# and for prepared statements
stmt = DBInterface.prepare(conn, "call get_employee()")
res = DBInterface.execute(stmt) |> columntable
@test length(res) > 0
@test length(res[1]) == 5
res = DBInterface.execute(stmt) |> columntable
@test length(res) > 0
@test length(res[1]) == 5
res = DBInterface.execute(stmt)
res = DBInterface.execute(stmt)

results = DBInterface.executemultiple(conn, """
select ID from Employee;
select DeptNo, OfficeNo from Employee where OfficeNo IS NOT NULL
""")
state = iterate(results)
@test state !== nothing
res, st = state
@test !st
# @test_broken length(res) == 5
ret = columntable(res)
# @test_broken length(ret[1]) == 5
state = iterate(results, st)
@test state !== nothing
res, st = state
@test !st
# @test_broken length(res) == 4
ret = columntable(res)
# @test_broken length(ret[1]) == 4

DBInterface.execute(conn, """CREATE TABLE Employee2
                 (
                     ID INT NOT NULL AUTO_INCREMENT,
                     望研 VARCHAR(255) CHARACTER SET utf8mb4,
                     PRIMARY KEY (ID)
                 )""")
DBInterface.execute(conn, "INSERT INTO Employee2 (望研) VALUES ('hey'), ('ho')")
ret = DBInterface.execute(conn, "select * from Employee2") |> columntable
@test length(ret.望研) == 2

DBInterface.execute(conn, """CREATE TABLE big_decimal
                 (
                     ID INT NOT NULL AUTO_INCREMENT,
                     `dec` DECIMAL(20, 2),
                     PRIMARY KEY (ID)
                 )""")
DBInterface.execute(conn, "INSERT INTO big_decimal (`dec`) VALUES (123456789012345678.91)")
ret = DBInterface.execute(conn, "select * from big_decimal") |> columntable
@test ret.dec[1] == d128"1.2345678901234567891e17"


text:
module StatisticalRethinking

using Reexport

@reexport using Distributions, RDatasets, DataFrames
@reexport using StatsBase, StatsPlots, StatsFuns
@reexport using CSV, DelimitedFiles, Serialization
@reexport using MCMCChains, KernelDensity
@reexport using LogDensityProblems
@reexport using TransformVariables, MCMCDiagnostics
@reexport using Parameters, Random

using DataStructures
import StatsBase: sample
import MCMCChains: describe

const src_path = @__DIR__

"""

# rel_path

Relative path using the StatisticalRethinking src/ directory. Copied from
[DynamicHMCExamples.jl](https://github.com/tpapp/DynamicHMCExamples.jl)

### Example to get access to the data subdirectory
```julia
rel_path("..", "data")
```
"""
rel_path(parts...) = normpath(joinpath(src_path, parts...))

include("link.jl")
include("scriptentry.jl")
include("generate.jl")
include("quap.jl")
include("sampling.jl")
include("scale.jl")
include("shading.jl")
include("HMC.jl")
include("generate_n_samples.jl")

export
  rel_path,
  ScriptEntry,
  scriptentry,
  script_dict,
  generate,
  quap,
  link,
  describe,
  #shade,
  scale,
  standardize

end # module


text:
import DataStructures

function (f::DiffableFunction)(vars...)
    f_gradfield = f.grad_field
    f_gradfield.inputs = collect(vars)

    xs = get_values.(vars)
    ys = as_tuple(forward(f, xs...))

    f_gradfield.generation = minimum((x -> x.generation), f_gradfield.inputs)
    f_gradfield.outputs = [Variable(y, creator=f, grad=nothing, generation=f_gradfield.generation - 1) for y in ys]

    return length(f_gradfield.outputs)  == 1 ? f_gradfield.outputs[1] : f_gradfield.outputs
end

"""
    backward(y::Variable; retain_grad=false)
Compute the back propagation, `y` as the end of the computational graph.
For a more detailed explanation, see the documentation (here).

```julia-repl
julia> x
name: nothing
values: 1
creator: User-Defined(nothing)

julia> x = Variable(1)
name: nothing
values: 1
creator: User-Defined(nothing)

julia> y = x +  2
name: nothing
values: 3
creator: JITrench.Add

julia> backward!(y)

julia> x.grad
name: nothing
values: 1
creator: User-Defined(nothing)
```
"""
function backward!(y::Variable; retain_grad=false)
    funcs = DataStructures.PriorityQueue{DiffableFunction,Int}()
    seen_set = Set{DiffableFunction}()
    if y.grad isa Nothing
        y.grad = Variable(ones_like(y.values))
    end
    DataStructures.enqueue!(funcs, y.creator, 1)
    push!(seen_set, y.creator)
    while !(isempty(funcs))
        f = DataStructures.dequeue!(funcs)
        gy = [output.grad for output in f.grad_field.outputs]
        gxs = backward(f, gy...)
        gxs = as_tuple(gxs)
        for (x, gx) in zip(f.grad_field.inputs, gxs)
            if x.grad isa Nothing
                x.grad =  gx
            else
                x.grad = x.grad + gx
            end
            if (!(isnothing(x.creator))) && (!(x.creator in seen_set))
                push!(seen_set, x.creator)
                DataStructures.enqueue!(funcs, x.creator, x.creator.grad_field.generation)
            end
        end
        if !(retain_grad)
            for y in f.grad_field.outputs
                y.grad = nothing
            end
        end
    end
    return nothing
end


text:
function featherwrite(filename::AbstractString, columns, colnames; description::AbstractString="", metadata::AbstractString="")
    ncol = length(columns)
    nrows = length(columns[1])
    cols = ArrowVector[arrowformat(_first_col_convert_pass(col)) for col in columns]

    open(filename, "w+") do io
        writepadded(io, FEATHER_MAGIC_BYTES)
        colmetadata = Metadata.Column[writecolumn(io, string(colnames[i]), cols[i]) for i in 1:ncol]
        ctable = Metadata.CTable(description, nrows, colmetadata, FEATHER_VERSION, metadata)
        len = writemetadata(io, ctable)
        write(io, Int32(len))  # these two writes combined are properly aligned
        write(io, FEATHER_MAGIC_BYTES)
    end
    return nothing
end

# NOTE: the below is very inefficient, but we are forced to do it by the Feather format
# Feather requires us to encode any vector that doesn't have a missing value
# as a normal list.
_first_col_convert_pass(col) = col
function _first_col_convert_pass(col::AbstractVector{Union{T,Missing}}) where T
    hasmissing = findfirst(ismissing, col)
    return hasmissing == nothing ? convert(AbstractVector{T}, col) : col
end

function Metadata.PrimitiveArray(A::ArrowVector{J}, off::Integer, nbytes::Integer) where J
    Metadata.PrimitiveArray(feathertype(J), Metadata.PLAIN, off, length(A), nullcount(A), nbytes)
end
function Metadata.PrimitiveArray(A::DictEncoding, off::Integer, nbytes::Integer)
    Metadata.PrimitiveArray(feathertype(eltype(references(A))), Metadata.PLAIN, off, length(A),
                            nullcount(A), nbytes)
end

writecontents(io::IO, A::Primitive) = writepadded(io, A)
writecontents(io::IO, A::NullablePrimitive) = writepadded(io, A, bitmask, values)
writecontents(io::IO, A::List) = writepadded(io, A, offsets, values)
writecontents(io::IO, A::NullableList) = writepadded(io, A, bitmask, offsets, values)
writecontents(io::IO, A::BitPrimitive) = writepadded(io, A, values)
writecontents(io::IO, A::NullableBitPrimitive) = writepadded(io, A, bitmask, values)
writecontents(io::IO, A::DictEncoding) = writecontents(io, references(A))
function writecontents(::Type{Metadata.PrimitiveArray}, io::IO, A::ArrowVector)
    a = position(io)
    writecontents(io, A)
    b = position(io)
    Metadata.PrimitiveArray(A, a, b-a)
end

function writecolumn(io::IO, name::AbstractString, A::ArrowVector{J}) where J
    vals = writecontents(Metadata.PrimitiveArray, io, A)
    Metadata.Column(String(name), vals, getmetadata(io, J, A), "")
end

function writemetadata(io::IO, ctable::Metadata.CTable)
    meta = FlatBuffers.build!(ctable)
    rng = (meta.head+1):length(meta.bytes)
    writepadded(io, view(meta.bytes, rng))
    Int32(length(rng))
end


text:
"""
    plsrosa(X, Y, weights = ones(size(X, 1)); nlv)
Partial Least Squares Regression (PLSR) with the ROSA algorithm (Liland et al. 2016).
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to consider.

The function has the following differences with the original
algorithm of Liland et al. (2016):
* Scores T are not normed to 1.
* Multivariate Y is allowed. In such a case, the entered
    Y-columns should have the same scale (for finding the winning blocks,
    the squared residuals are summed over the columns)

Vector `weights` (row-weighting) is internally normalized to sum to 1.
See the help of `plskern` for details.

`X` and `Y` are internally centered.

## References

Liland, K.H., Næs, T., Indahl, U.G., 2016. ROSA—a fast extension of partial least
squares regression for multiblock data analysis. Journal of Chemometrics 30,
651–662. https://doi.org/10.1002/cem.2824

"""
function plsrosa(X, Y, weights = ones(size(X, 1)); nlv)
    plsrosa!(copy(X), copy(Y), weights; nlv = nlv)
end

function plsrosa!(X, Y, weights = ones(size(X, 1)); nlv)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    n = size(X, 1)
    p = size(X, 2)
    q = size(Y, 2)
    nlv = min(nlv, n, p)
    weights = mweights(weights)
    D = Diagonal(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    # Pre-allocation
    XtY = similar(X, p, q)
    T = similar(X, n, nlv)
    P = similar(X, p, nlv)
    W = copy(P)
    C = similar(X, q, nlv)
    TT = similar(X, nlv)
    t = similar(X, n)
    dt = similar(X, n)
    zp = similar(X, p)
    w = similar(X, p)
    c = similar(X, q)
    # End
    @inbounds for a = 1:nlv
        XtY .= X' * (D * Y)
        if q == 1
            w .= vec(XtY)
            w ./= sqrt(dot(w, w))
        else
            w .= svd!(XtY).U[:, 1]
        end
        mul!(t, X, w)
        if a > 1
            z = vcol(T, 1:(a - 1))
            t .= t .- z * inv(z' * D * z) * z' * (D * t)
            z = vcol(W, 1:(a - 1))
            w = w .- z * (z' * w)
            w ./= sqrt(dot(w, w))
        end
        dt .= weights .* t
        tt = dot(t, dt)
        mul!(c, Y', dt)
        c ./= tt
        mul!(zp, X', dt)
        zp ./= tt
        Y .-= t * c'
        P[:, a] .= zp
        T[:, a] .= t
        W[:, a] .= w
        C[:, a] .= c
        TT[a] = tt
     end
     R = W * inv(P' * W)
     Plsr(T, P, R, W, C, TT, xmeans, ymeans, weights, nothing)
end


text:
function sum_2(x::A) where {T, N, A<:AbstractArray{T,N}}
    n = length(x)
    if n < 2
        iszero(n) && return zero(T)
        isone(n)  && return x[1]
    end
    return summation_(n, x)
end

@inline sum_(x::A) where {T, N, A<:AbstractArray{T,N}} = sum_2(x)

function summation_(n::Int, x::A) where {T, N, A<:AbstractArray{T,N}}
   hi = x[1]
   lo = zero(T)
   for i in 2:n
       hi, low = add_(hi, x[i])
       lo += low
   end
   hi, lo = add_hilo_(hi, lo)
   return hi, lo
end


text:
using StaticArrays, HybridArrays, Test, LinearAlgebra

@testset "AbstractArray interface" begin
    @testset "size and length" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2 3; 4 5 6])

        @test length(M) == 6
        @test size(M) == (2, 3)
        @test Base.isassigned(M, 2, 2) == true
        @test (@inferred Size(M)) == Size(Tuple{2, StaticArrays.Dynamic()})
        @test (@inferred Size(typeof(M))) == Size(Tuple{2, StaticArrays.Dynamic()})
    end

    @testset "reshape" begin
        # TODO?
    end

    @testset "convert" begin
        MM = [1 2 3; 4 5 6]
        M = HybridMatrix{2, StaticArrays.Dynamic()}(MM)
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()}}, MM)).data == M
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Int}, MM)).data == M
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64}, MM)).data == M
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64,2}, MM)).data == M
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64,2,2}, MM)).data == M
        @test (@inferred convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64,2,2,Matrix{Float64}}, MM)).data == M
        @test convert(typeof(M), M) === M
        if VERSION >= v"1.1"
            @test convert(HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64}, M) == M
        end
        @test convert(Array, M) === M.data
        @test convert(Array{Int}, M) === M.data
        @test convert(Matrix, M) === M.data
        @test convert(Matrix{Int}, M) === M.data

        @test Array(M) == M
        @test Array(M) !== M.data
        @test Matrix(M) == M
        @test Matrix(M) !== M.data
        @test Matrix{Int}(M) == M
        @test Matrix{Int}(M) !== M.data
        @test_throws MethodError Vector(M)
    end

    @testset "copy" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2; 3 4])
        @test @inferred(copy(M))::HybridMatrix == M
        @test copy(M).data !== M.data
    end

    @testset "similar" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])

        @test isa(@inferred(similar(M)), HybridMatrix{2, StaticArrays.Dynamic(), Int})
        @test isa(@inferred(similar(M, Float64)), HybridMatrix{2, StaticArrays.Dynamic(), Float64})
        @test isa(@inferred(similar(M, Float64, Size(Tuple{3, 3}))), HybridMatrix{3, 3, Float64})
    end

    @testset "IndexStyle" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])
        MT = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4]')
        @test (@inferred IndexStyle(M)) === IndexLinear()
        @test (@inferred IndexStyle(MT)) === IndexCartesian()
        @test (@inferred IndexStyle(typeof(M))) === IndexLinear()
        @test (@inferred IndexStyle(typeof(MT))) === IndexCartesian()
    end

    @testset "vec" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])
        Mv = vec(M)
        @test Mv == [1, 3, 2, 4]
        @test Mv isa Vector{Int}
        Mv[2] = 100
        @test M[2, 1] == 100
    end

    @testset "errors" begin
        @test_throws TypeError HybridArrays.new_out_size_nongen(Size{Tuple{1,2}}, 'a')
    end

    @testset "strides" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])

        @test strides(M) == strides(M.data)
    end

    @testset "pointer" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])
        MT = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4]')

        @test pointer(M) == pointer(M.data)
        if VERSION >= v"1.5"
            # pointer on Adjoint is not available on earilier versions of Julia
            @test pointer(MT) == pointer(MT.data)
        end
    end

    @testset "unsafe_convert" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Int}([1 2; 3 4])
        @test Base.unsafe_convert(Ptr{Int}, M) === pointer(M.data)
    end

    @test HybridArrays._totally_linear() === true
end


text:
include("4.3.1.1 (a+b tan)^n.jl")
include("4.3.1.2 (d sec)^m (a+b tan)^n.jl")
include("4.3.1.3 (d sin)^m (a+b tan)^n.jl")
include("4.3.10 (c+d x)^m (a+b tan)^n.jl")
include("4.3.11 (e x)^m (a+b tan(c+d x^n))^p.jl")
include("4.3.12 (d+e x)^m tan(a+b x+c x^2)^n.jl")
include("4.3.2.1 (a+b tan)^m (c+d tan)^n.jl")
include("4.3.2.3 (g tan)^p (a+b tan)^m (c+d tan)^n.jl")
include("4.3.3.1 (a+b tan)^m (c+d tan)^n (A+B tan).jl")
include("4.3.4.1 (a+b tan)^m (A+B tan+C tan^2).jl")
include("4.3.4.2 (a+b tan)^m (c+d tan)^n (A+B tan+C tan^2).jl")
include("4.3.7 (d trig)^m (a+b (c tan)^n)^p.jl")
include("4.3.9 trig^m (a+b tan^n+c tan^(2 n))^p.jl")

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/.jl