text:
using Revise
using Plots; plotly();
using Statistics
using StatsPlots
using LinearAlgebra
using SparseArrays
using Distributions
using Random
using RandomNumbers.PCG
using DelimitedFiles
using BenchmarkTools

# using Distributed; addprocs(3);
# @everywhere include("hmc.jl")
includet("hmc.jl")
includet("convergence.jl")
# @everywhere include("models.jl")
includet("models.jl")

ndim = 64
Minv = Matrix(Diagonal(ones(ndim))); # ones(ndim);

# set up model
M = sparse([1.0 0.75; 0.75 1.0]);
S = Matrix(blockdiag([M for i in 1:ndim/2]...));
μ = zeros(ndim);

f = mvgaussian_(μ, Matrix(S));

normal_.(randn(10), 0, 10)

function Normal_(x::Vector{Float64}, μ::Vector{Float64}, Σ::Matrix{Float64})
    d = x - μ
    return 0.5 * (d' * (Σ \ d))
end

f(x) = Normal_(x, μ, S)

# single chain
samples, d = hmc(f, ndim;
                 M = cholesky(Symmetric(Minv)).L,
                 # M = Minv,
                 control = Dict(:skewsymmetric => true, :iterations => 2000));

samples2, d2 = hmc(f, ndim;
                 # M = cholesky(Symmetric(Minv)).L,
                 M = Minv,
                 control = Dict(:skewsymmetric => false, :iterations => 2000));

# multiple chains
# @benchmark
samples, d = stan(f, ndim;
                  M = Minv, control = Dict(:skewsymmetric => false, :iterations => 2000));

round.(mean(samples[:, :, :], dims=1), digits=1)
round.(std(samples[:, :, :], dims=1), digits=1)

round.(mean(samples2[:, :, :], dims=1), digits=1)
round.(std(samples2[:, :, :], dims=1), digits=1)

map(i -> round(ess_bulk(samples[:, :, 1])), 1:ndim)
map(i -> round(ess_tail(samples[:, :, 1])), 1:ndim)

map(i -> round(ess_std(samples[:, i, 1]), digits=2), 1:ndim)
map(i -> round(rhat(samples[:, i, 1]), digits=2), 1:ndim)

plot(1:10000, samples[:, 2, 2])

sum(d[:divergent], dims=1)
plot(d[:leapfrog][idx, :], seriestype=:density, title="Number leapfrog steps")
plot(d[:acceptstat][idx, :], seriestype=:density, title="Accept stats")
plot(1:15000, log10.(d[:stepsize][:, 2]), seriestype=:scatter, title="Stepsizes")
plot(d[:treedepth][2001:end, :], seriestype=:density, title="Tree depths")
plot(d[:massmatrix]', seriestype=:density, title="Mass matrix")

# TODO a future test
μ = zeros(2)
S = [1.0 0.99; 0.99 1.0];
N = MvNormal(μ, S)
w = WelfordState(zero(μ), zeros(length(μ)), 0)
W = WelfordState(zero(μ), zero(S), 0)
X = zeros(2, 10)

for n in 1:10
    x = rand(N)
    global X[:, n] = x
    global w = accmoments(w, x)
    global W = accmoments(W, x)
end

cov(X', corrected = true)
samplevariance(W, Dict(:skewsymmetric => true, :regularize => false))
samplevariance(w, Dict(:skewsymmetric => false, :regularize => false))

max_stars_repo_name: roualdes/yastan

max_issues_repo_name: roualdes/yastan

max_forks_repo_name: roualdes/yastan


text:
function get_FLay02H()

m = Model()

# ----- Variables ----- #
@variable(m, objvar)
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
@variable(m, x[x_Idx])
b_Idx = Any[43, 44, 45, 46]
@variable(m, b[b_Idx], Bin)
JuMP.set_lower_bound(x[36], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[16], 0.0)
JuMP.set_lower_bound(x[32], 0.0)
JuMP.set_lower_bound(x[27], 0.0)
JuMP.set_lower_bound(x[14], 0.0)
JuMP.set_lower_bound(x[17], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_lower_bound(x[25], 0.0)
JuMP.set_lower_bound(x[38], 0.0)
JuMP.set_lower_bound(x[30], 0.0)
JuMP.set_lower_bound(x[26], 0.0)
JuMP.set_lower_bound(x[23], 0.0)
JuMP.set_lower_bound(x[42], 0.0)
JuMP.set_lower_bound(x[34], 0.0)
JuMP.set_lower_bound(x[11], 0.0)
JuMP.set_lower_bound(x[29], 0.0)
JuMP.set_lower_bound(x[22], 0.0)
JuMP.set_lower_bound(x[12], 0.0)
JuMP.set_lower_bound(x[37], 0.0)
JuMP.set_lower_bound(x[19], 0.0)
JuMP.set_lower_bound(x[40], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_lower_bound(x[20], 0.0)
JuMP.set_lower_bound(x[24], 0.0)
JuMP.set_lower_bound(x[41], 0.0)
JuMP.set_lower_bound(x[39], 0.0)
JuMP.set_lower_bound(x[31], 0.0)
JuMP.set_lower_bound(x[18], 0.0)
JuMP.set_lower_bound(x[9], 0.0)
JuMP.set_lower_bound(x[15], 0.0)
JuMP.set_lower_bound(x[1], 0.0)
JuMP.set_lower_bound(x[33], 0.0)
JuMP.set_lower_bound(x[13], 0.0)
JuMP.set_lower_bound(x[21], 0.0)
JuMP.set_lower_bound(x[28], 0.0)
JuMP.set_lower_bound(x[35], 0.0)
JuMP.set_lower_bound(x[10], 0.0)
JuMP.set_upper_bound(x[1], 29.0)
JuMP.set_upper_bound(x[2], 29.0)
JuMP.set_upper_bound(x[3], 29.0)
JuMP.set_upper_bound(x[4], 29.0)
JuMP.set_lower_bound(x[5], 1.0)
JuMP.set_upper_bound(x[5], 40.0)
JuMP.set_lower_bound(x[6], 1.0)
JuMP.set_upper_bound(x[6], 50.0)
JuMP.set_lower_bound(x[7], 1.0)
JuMP.set_upper_bound(x[7], 40.0)
JuMP.set_lower_bound(x[8], 1.0)
JuMP.set_upper_bound(x[8], 50.0)
JuMP.set_upper_bound(x[9], 30.0)
JuMP.set_upper_bound(x[10], 30.0)

# ----- Constraints ----- #
@constraint(m, e1, -2*x[9]-2*x[10]+objvar == 0.0)
@constraint(m, e2, -x[1]-x[5]+x[9] >= 0.0)
@constraint(m, e3, -x[2]-x[6]+x[9] >= 0.0)
@constraint(m, e4, -x[3]-x[7]+x[10] >= 0.0)
@constraint(m, e5, -x[4]-x[8]+x[10] >= 0.0)
@NLconstraint(m, e6, 40/x[7]-x[5] <= 0.0)
@NLconstraint(m, e7, 50/x[8]-x[6] <= 0.0)
@constraint(m, e8, x[1]-x[11]-x[12]-x[13]-x[14] == 0.0)
@constraint(m, e9, x[2]-x[15]-x[16]-x[17]-x[18] == 0.0)
@constraint(m, e10, x[3]-x[19]-x[20]-x[21]-x[22] == 0.0)
@constraint(m, e11, x[4]-x[23]-x[24]-x[25]-x[26] == 0.0)
@constraint(m, e12, x[5]-x[27]-x[28]-x[29]-x[30] == 0.0)
@constraint(m, e13, x[6]-x[31]-x[32]-x[33]-x[34] == 0.0)
@constraint(m, e14, x[7]-x[35]-x[36]-x[37]-x[38] == 0.0)
@constraint(m, e15, x[8]-x[39]-x[40]-x[41]-x[42] == 0.0)
@constraint(m, e16, x[11]-29*b[43] <= 0.0)
@constraint(m, e17, x[12]-29*b[44] <= 0.0)
@constraint(m, e18, x[13]-29*b[45] <= 0.0)
@constraint(m, e19, x[14]-29*b[46] <= 0.0)
@constraint(m, e20, x[15]-29*b[43] <= 0.0)
@constraint(m, e21, x[16]-29*b[44] <= 0.0)
@constraint(m, e22, x[17]-29*b[45] <= 0.0)
@constraint(m, e23, x[18]-29*b[46] <= 0.0)
@constraint(m, e24, x[19]-29*b[43] <= 0.0)
@constraint(m, e25, x[20]-29*b[44] <= 0.0)
@constraint(m, e26, x[21]-29*b[45] <= 0.0)
@constraint(m, e27, x[22]-29*b[46] <= 0.0)
@constraint(m, e28, x[23]-29*b[43] <= 0.0)
@constraint(m, e29, x[24]-29*b[44] <= 0.0)
@constraint(m, e30, x[25]-29*b[45] <= 0.0)
@constraint(m, e31, x[26]-29*b[46] <= 0.0)
@constraint(m, e32, x[27]-40*b[43] <= 0.0)
@constraint(m, e33, x[28]-40*b[44] <= 0.0)
@constraint(m, e34, x[29]-40*b[45] <= 0.0)
@constraint(m, e35, x[30]-40*b[46] <= 0.0)
@constraint(m, e36, x[31]-40*b[43] <= 0.0)
@constraint(m, e37, x[32]-40*b[44] <= 0.0)
@constraint(m, e38, x[33]-40*b[45] <= 0.0)
@constraint(m, e39, x[34]-40*b[46] <= 0.0)
@constraint(m, e40, x[35]-40*b[43] <= 0.0)
@constraint(m, e41, x[36]-40*b[44] <= 0.0)
@constraint(m, e42, x[37]-40*b[45] <= 0.0)
@constraint(m, e43, x[38]-40*b[46] <= 0.0)
@constraint(m, e44, x[39]-40*b[43] <= 0.0)
@constraint(m, e45, x[40]-40*b[44] <= 0.0)
@constraint(m, e46, x[41]-40*b[45] <= 0.0)
@constraint(m, e47, x[42]-40*b[46] <= 0.0)
@constraint(m, e48, x[11]-x[15]+x[27] <= 0.0)
@constraint(m, e49, -x[12]+x[16]+x[32] <= 0.0)
@constraint(m, e50, x[21]-x[25]+x[37] <= 0.0)
@constraint(m, e51, -x[22]+x[26]+x[42] <= 0.0)
@constraint(m, e52, b[43]+b[44]+b[45]+b[46] == 1.0)

# ----- Objective ----- #
@objective(m, Min, objvar)

return m

end

ext: jl


text:
using AugmentedGaussianProcesses
using Test
using LinearAlgebra, Distributions
using Zygote
using PDMats
using MLDataUtils
using Random: seed!
seed!(42)

include("testingtools.jl")

# Global flags for the tests
@testset "AugmentedGaussianProcesses.jl tests" begin
    @info "Testing data"
    @testset "Data test" begin
        include(joinpath("data", "datacontainer.jl"))
        include(joinpath("data", "utils.jl"))
    end

    @info "Function tests"
    @testset "Functions" begin
        for f in readdir(joinpath(@__DIR__, "functions"))
            include(joinpath("functions", f))
        end
    end

    @info "Hyperparameter tests"
    @testset "Hyperparameters" begin
    end

    @info "Inference tests"
    @testset "Inference" begin
        for f in readdir(joinpath(@__DIR__, "inference"))
            include(joinpath("inference", f))
        end
    end

    # @info "Likelihood tests with Forward Diff"
    # AGP.setadbackend(:ForwardDiff)
    # @testset "Likelihoods (ForwardDiff)" begin
    #     for f in readdir(joinpath(@__DIR__, "likelihood"))
    #         include(joinpath("likelihood",f))
    #     end
    # end

    @info "Likelihood tests with Zygote"
    AGP.setadbackend(:Zygote)
    @testset "Likelihoods (Zygote)" begin
        for f in readdir(joinpath(@__DIR__, "likelihood"))
            include(joinpath("likelihood", f))
        end
    end

    @info "Model tests"
    @testset "Models" begin
        for f in readdir(joinpath(@__DIR__, "models"))
            include(joinpath("models", f))
        end
    end

    @info "Prior tests"
    @testset "Prior" begin
        for f in readdir(joinpath(@__DIR__, "prior"))
            include(joinpath("prior", f))
        end
    end
    include("training.jl")
    include("onlinetraining.jl")
    include("predictions.jl")
end


text:
# Convert the original example datasets to compressed CSV files

# See data/README.md for the sources of the input data files
# To regenerate the output files:
# 1) Have all input files ready in the data folder
# 2) Instantiate the package environment for data/src
# 3) Run this script with the root folder of the repository being the working directory

using CSV, CodecZlib, DataFrames, ReadStatTables

function adh()
    cols = [:czone, :year, :y, :x, :z, :wei, :l_sh_routine33, :t2, :Lsh_manuf]
    loc = DataFrame(readstat("data/location_level.dta", usecols=cols))
    loc.czone = convert(Vector{Int}, loc.czone)
    loc.year = convert(Vector{Int}, loc.year)
    open(GzipCompressorStream, "data/adh_location.csv.gz", "w") do stream
        CSV.write(stream, loc)
    end
    shares = DataFrame(readstat("data/Lshares.dta"))
    shares = shares[shares.ind_share.!=0,:]
    shares.czone = convert(Vector{Int}, shares.czone)
    shares.year = convert(Vector{Int}, shares.year)
    shares.sic87dd = convert(Vector{Int}, shares.sic87dd)
    open(GzipCompressorStream, "data/adh_share.csv.gz", "w") do stream
        CSV.write(stream, shares)
    end
    shocks = DataFrame(readstat("data/shocks.dta", usecols=[:year,:sic87dd,:g]))
    shocks.year = convert(Vector{Int}, shocks.year)
    shocks.sic87dd = convert(Vector{Int}, shocks.sic87dd)
    open(GzipCompressorStream, "data/adh_shock.csv.gz", "w") do stream
        CSV.write(stream, shocks)
    end
end

function main()
    adh()
end

main()


text:
"""
Methods to compute spin-weighted spheroidal harmonics
"""
module Spheroidal

include("CustomTypes.jl")

using .CustomTypes
using SparseArrays

import LinearAlgebra: eigen
using GenericSchur

"""
    compute_l_min(s,m)

Minimum l value
"""
function compute_l_min(s,m)
   return max(abs(s),abs(m))
end

"""
    compute_A(s::myI, l::myI, m::myI)

See the note.
"""
@inline function compute_A(s::myI, l::myI, m::myI)
   return sqrt(
      (l+tomyF(1)-s)*(l+tomyF(1)+s)*(l+tomyF(1)+m)*(l+tomyF(1)-m)
      /
      (((l+tomyF(1))^2)*(tomyF(2)*l+tomyF(1))*(tomyF(2)*l+tomyF(3)))
   )
end
"""
    compute_B(s::myI, l::myI, m::myI)

See the note.
"""
@inline function compute_B(s::myI, l::myI, m::myI)
   if (l==0)
      return tomyF(0)
   else
      return -(m*s)/(l*(l+tomyF(1)))
   end
end
"""
    compute_C(s::myI, l::myI, m::myI)

See the note.
"""
@inline function compute_C(s::myI, l::myI, m::myI)
   if (l==0)
      return tomyF(0)
   else
      return sqrt(
         tomyF(l-s)*tomyF(l+s)*tomyF(l+m)*tomyF(l-m)
         /
         ((l^2)*(tomyF(2)*l+tomyF(1))*(tomyF(2)*l-tomyF(1)))
      )
   end
end
"""
    mat_Y(
      nl::myI,
      s::myI,
      m::myI,
    )::SparseMatrixCSC

Returns the matrix for y in spectral space.
"""
function mat_Y(
      nl::myI,
      s::myI,
      m::myI,
   )::SparseMatrixCSC

   I = Vector{myI}(undef,0)
   J = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)

   lmin = compute_l_min(s,m)

   for i=1:nl

      l = (i-1) + lmin

      ## (l,l) component
      append!(I,i)
      append!(J,i)
      append!(V,compute_B(s,l,m))

      ## (l,l-1) component
      if (i>1)
         append!(I,i)
         append!(J,i-1)
         append!(V,compute_C(s,l,m))
      end
      ## (l,l+1) component
      if (i+1<=nl)
         append!(I,i)
         append!(J,i+1)
         append!(V,compute_A(s,l,m))
      end
   end
   return sparse(I,J,V)
end

"""
    mat_L(
      nl::myI,
      s::myI,
      m::myI,
    )::SparseMatrixCSC

Returns the matrix for spherical laplacian (s=0) in spectral space.
"""
function mat_L(
      nl::myI,
      s::myI,
      m::myI,
   )::SparseMatrixCSC

   I = Vector{myI}(undef,0)
   J = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)

   lmin = compute_l_min(s,m)

   for i=1:nl

      l = (i-1) + lmin

      ## (l,l) component
      append!(I,i)
      append!(J,i)
      append!(V,tomyF(l - s)*tomyF(l + s + 1.0))
   end
   return sparse(I,J,V)
end
"""
    compute_M_matrix(
      nl::myI,
      s::myI,
      m::myI,
      c::myC
    )

M matrix for computing spheroidal-spherical mixing and seperation coefficients.
"""
function compute_M_matrix(
      nl::myI,
      s::myI,
      m::myI,
      c::myC
   )
   Y1 = mat_Y(nl+2,s,m) # look at larger to capture last term in matrix mult
                        # at higher l
   Y2 = Y1*Y1
   Y1 = Y1[1:nl,1:nl]
   Y2 = Y2[1:nl,1:nl]
   L  = mat_L(nl,s,m)

   return sparse(
         transpose(
            -
            (c^2)*Y2
            +
            tomyF(2)*c*s*Y1
            +
            L
         )
      )

end

"""
    eig_vals_vecs(
      nl::myI,
      neig::myI,
      s::myI,
      m::myI,
      c::myC
    )

Compute eigenvectors and eigenvalues for the spheroidal equation.
Returns eigenvalues λs (smallest to largest), and eigenvectors as
an array (access nth eigenvector through v[:,n]); i.e. returns
(λs, vs)
"""
function eig_vals_vecs(
      nl::myI,
      neig::myI,
      s::myI,
      m::myI,
      c::myC
   )

   Mat = compute_M_matrix(nl, s, m, c)

   for i=1:nl
      Mat[i,i] += tomyF(1) # shift to avoid issues with a zero eigenvalue
   end
   # λs, vs
   t = eigen(Matrix(Mat), permute=true, scale=true, sortby=abs)
   return t.values .- tomyF(1), t.vectors
end

end # module


text:
using SafeTestsets

@safetestset "Testing" begin include("testing.jl") end
@safetestset "Layouts" begin include("layouts.jl") end
@safetestset "Indexing" begin include("indexing.jl") end
@safetestset "TriMatrix" begin include("trimatrix.jl") end
@safetestset "LinearAlgebra" begin include("linearalgebra.jl") end


text:
function _generate_precompile_execution_content(config::Config)
    pkgs = config.pkgs
    no_test = config.no_test
    exclude_packages_from_sysimage = config.exclude_packages_from_sysimage
    precompile_execution_env_vars = config.precompile_execution_env_vars
    pkg_names_to_import = Vector{String}(undef, 0)
    pkg_names_to_test = Vector{String}(undef, 0)
    for pkg in pkgs
        pkg_name = pkg[:name]
        if !(pkg_name in exclude_packages_from_sysimage)
            push!(pkg_names_to_import, pkg_name)
        end
        if ( !(pkg_name in no_test) ) && ( !(pkg_name in exclude_packages_from_sysimage) )
            push!(pkg_names_to_test, pkg_name)
        end
    end
    set_precompile_execution_env_vars_string = ""
    for (name, value) in precompile_execution_env_vars
        set_precompile_execution_env_vars_string *= "ENV[\"$(name)\"] = \"$(value)\"\n"
    end
    import_string = ""
    for pkg_name in pkg_names_to_import
        import_string *= "import $(pkg_name) # pkg_names_to_import\n"
    end
    lines = String[
        "$(import_string)",
        "import Pkg",
        "for (uuid, info) in Pkg.dependencies()",
        "if info.name in $(pkg_names_to_test) # pkg_names_to_test",
        "$(set_precompile_execution_env_vars_string)",
        "include(joinpath(info.source, \"test\", \"runtests.jl\"))",
        "end",
        "end",
    ]
    content = join(lines, "\n") * "\n"
    return content
end


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
#=
pokemon_csv_to_species_mjs: Converts the entries in pokemon.csv to an array in species.js.
- Julia version: 1.4.2
- Author: Wallace Watler
- Date: 2021-02-04
=#

function main()
    println("Reading entries from pokemon.csv...")
    lines = readlines("pokemon.csv")[2:end] # Skip header row

    lines = map(line -> begin
        values::Vector{Any} = split(line, ",")
        values[1] = values[1][2:end] # Remove '#' from dexNumber
        values[13] = lowercase(values[13])
        values
    end, lines)

    output =  "// All species, generated from pokemon.csv.\n"
    output *= "const SPECIES = [\n"
    output *= join(map(line -> begin
              "    {\n"*
              "        dexNumber: \"$(line[1])\",\n"*
              "        gen: $(line[2]),\n"*
              "        name: \"$(line[3])\",\n"*
              "        baseStats: [$(line[4]), $(line[5]), $(line[6]), $(line[7]), $(line[8]), $(line[9]), $(line[10])],\n"*
              "        type1: TYPE_NAMES.indexOf(\"$(line[11])\"),\n"*
              "        type2: " * (isempty(line[12]) ? "undefined" : "TYPE_NAMES.indexOf(\"$(line[12])\")") * ",\n"*
              "        isMega: $(line[13])\n"*
              "    }"
                       end, lines), ",\n")
    output *= "\n];\n"

    open("species.js", "w") do io
        println("Writing entries to species.js...")
        write(io, output)
    end
end

main()


text:
######################################
## comparing analytical (fictitious sources) and computed fields. minimum N
## routine runs on multiple cores and uses a lot of memory -- reduce max Nvec5,
## Nvecs if necessary
using PyPlot,Distributed,SpecialFunctions
import JLD
addprocs(4)
@everywhere using ParticleScattering
output_dir = homedir()

function unique_subset(N,v,rev = true)
	all(sort(v, rev = rev) .== v) || error("Only handles sorted vectors.")
	v_u = [v[1]]
	N_u = [N[1]]
	prev = v[1]
	for i = 2:length(N)
		if v[i] !== prev
			#encountered new val
			push!(v_u, v[i])
			push!(N_u, N[i])
			prev = v[i]
		end
	end
    N_u,v_u
end
# binary search here is suboptimal since we know that new P is "close" to old one,
# frequently P or P+1
function findMinP(N, errN, shapefun, N_points, k0, kin; P_last = 1, P_max = 100)
	E_multipole = Array{Complex{Float64}}(undef, N_points)
	errP = zeros(Float64, length(errN))
	Pmin = zeros(Int, length(errN))
	for iN in eachindex(errN)
		s = shapefun(N[iN])

		err_points = [2.0*s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)]
		#compute direct solution for comparison
		inner = get_potentialPW(k0, kin, s, 0.0)
		E_quadrature = scatteredfield(inner, k0, s, err_points)
		for P = P_last:P_max
			err = ParticleScattering.minimumP_helper(k0, kin, s, P, N_points,
					err_points, E_quadrature, E_multipole)
			if err <= errN[iN]
				errP[iN] = err
				Pmin[iN] = P
				P_last = P
				break
			elseif P == P_max
				@warn("Failed to find P for iN = $iN")
				return Pmin, errP
			end
		end
		display("findMinP: done with iN=$iN/$(length(errN)), matched P=$(Pmin[iN]) with $(errP[iN]) ≤ $(errN[iN])")
	end
	Pmin, errP
end

@everywhere begin
	k0 = 10.0
	kin = 1.5*k0
	l0 = 2π/k0
	a1 = 0.3l0
	a2 = 0.1l0
	N_points = 20_000
	myshapefun5(N) = rounded_star(a1,a2,5,N)
	myshapefun_squircle(N) = squircle(a1+0.5a2,N)
	Nvec5 = unique(round.(Int, 10 .^ range(log10(20), stop=log10(5000), length=200)))
	Nvecs = unique(round.(Int, 10 .^ range(log10(20), stop=log10(5000), length=200)))
end

dt_Nvec5 = @elapsed begin
	s = myshapefun5(400) #just for radius
	err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)]
	E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points)
	E_comp = Array{ComplexF64}(undef, length(E_ana))

	innerfunc1 = function(i)
		ParticleScattering.minimumN_helper(Nvec5[i], k0, kin,
						myshapefun5, err_points, E_comp, E_ana)
	end
	errNvec5 = pmap(innerfunc1, eachindex(Nvec5),
				on_error = ex->(isa(ex, OutOfMemoryError) ? NaN : rethrow(ex)))
	#or
	#errNvec5 = pmap(innerfunc1, eachindex(Nvec5))
end
display("finished calculating minimum N for rounded star in $dt_Nvec5 seconds")
any(isnan.(errNvec5)) && @warn "encountered OutOfMemoryError, some values are NaN"

dt_Nvecs = @elapsed begin
	s = myshapefun_squircle(400) #just for radius
	err_points = [s.R*f(i*2*pi/N_points) for i=0:(N_points-1), f in (cos,sin)]
	E_ana = (0.25im*besselh(0,1,k0*s.R))*ones(Complex{Float64},N_points)
	E_comp = Array{ComplexF64}(undef, length(E_ana))

	innerfunc2 = function(i)
		ParticleScattering.minimumN_helper(Nvecs[i], k0, kin,
						myshapefun_squircle, err_points, E_comp, E_ana)
	end
	errNvecs = pmap(innerfunc2, eachindex(Nvecs),
				on_error = ex->(isa(ex, OutOfMemoryError) ? NaN : rethrow(ex)))
	#or
	#errNvecs = pmap(innerfunc2, eachindex(Nvecs))
end
display("finished calculating minimum N for squircle in $dt_Nvecs seconds")
any(isnan.(errNvecs)) && @warn "encountered OutOfMemoryError, some values are NaN"

N5, errN5 = unique_subset(Nvec5, errNvec5)
Ns, errNs = unique_subset(Nvecs, errNvecs)
inds = something(findfirst(errNs .< 5e-10), length(errNs))
ind5 = something(findfirst(errN5 .< 5e-10), length(errN5))
N5 = N5[1:ind5]
Ns = Ns[1:inds]
errN5 = errN5[1:ind5]
errNs = errNs[1:inds]

P5,errP5 = findMinP(N5, errN5, myshapefun5, N_points, k0, kin)
Ps,errPs = findMinP(Ns, errNs, myshapefun_squircle, N_points, k0, kin)

JLD.@save joinpath(output_dir, "mindata.jld") k0 kin l0 a1 a2 N5 errN5 Ns errNs Ps errPs P5 errP5

##################################################
#plot
fig = figure(figsize=[3.5, 2.5])
loglog(errN5, N5, "b-.", linewidth=1,
		label="\$N_{\\mathrm{min}} \\, \\mathrm{(star)}\$")
loglog(errNs, Ns, "r:", linewidth=1,
		label="\$N_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$")
loglog(errN5, P5, "tab:green", linewidth=1,
		label="\$P_{\\mathrm{min}} \\, \\mathrm{(star)}\$")
loglog(errNs, Ps, "k--", linewidth=1,
		label="\$P_{\\mathrm{min}} \\, \\mathrm{(squircle)}\$")
xlim([5e-10;1e-1])
xlabel("\$\\Delta \\, u\$")
tick_params(which="both", direction="in")
legend(fontsize=8)
grid()


text:
module ISA

include("components/analyticComponents.jl")
export AMFMtriplet, AMFMcomp

include("components/numericalComponents.jl")
export numComp, numTriplet

include("models/analyticModels.jl")
export compSet, AMFMmodel

include("models/fourierSeries.jl")
export fourierSeries

include("models/numericalModels.jl")
export numSet, numModel

include("demod/demodulation.jl")
export derivApprox, AMFMdemod

include("other/specialFunctions.jl")
export u, 𝒩ᵤ, 𝒩, δ, δn

include("plots/GRrecipes.jl")

end # module


text:
"""
    A = JopDifference(R::JetSpace[, dim=1])

`A*m` is similar to `diff(a,dims=dim)` where `a` is a vector in the space `R`.
In other words, it is the one-sided difference of `a` along the dimension `dim`.
"""
function JopDifference(spc::JetAbstractSpace, dim::Integer=1)
    JopLn(df! = JopDifference_df!, df′! = JopDifference_df′!, dom = spc, rng = spc, s = (dim=dim,))
end
export JopDifference

function JopDifference_df!(d, m; dim, kwargs...)
    d .= 0
    n,Rpre,Rpost = JopDifference_indices(m, dim)
    JopDifference_df_dim!(d, m, n, Rpre, Rpost)
end

function JopDifference_df′!(m, d; dim, kwargs...)
    m .= 0
    n,Rpre,Rpost = JopDifference_indices(m, dim)
    JopDifference_df′_dim!(m, d, n, Rpre, Rpost)
end

function JopDifference_indices(m, dim)
    Rpre = CartesianIndices(size(m)[1:dim-1])
    Rpost = CartesianIndices(size(m)[dim+1:end])
    n = size(m, dim)
    n,Rpre, Rpost
end

function JopDifference_df_dim!(d, m, n, Rpre, Rpost)
    for Ipost in Rpost
        for i = 1:n-1
            for Ipre in Rpre
                d[Ipre, i, Ipost] += m[Ipre, i+1, Ipost]
                d[Ipre, i, Ipost] -= m[Ipre, i, Ipost]
            end
        end
    end
    d
end

function JopDifference_df′_dim!(m, d, n, Rpre, Rpost)
    for Ipost in Rpost
        for i = 1:n-1
            for Ipre in Rpre
                m[Ipre, i+1, Ipost] += d[Ipre, i, Ipost]
                m[Ipre, i, Ipost] -= d[Ipre, i, Ipost]
            end
        end
    end
    m
end


text:
# rotate2d.jl
# Copyright 2019-03-05, Jeff Fessler, University of Michigan

using Test: @test

"""
`(xr,yr) = rotate2d(x, y, theta)`
2D rotation
"""
function rotate2d(x, y, theta)
	xr = cos(theta) * x + sin(theta) * y
	yr = -sin(theta) * x + cos(theta) * y
	return (xr, yr)
end

"""
`rotate2d(:test)`
self test
"""
function rotate2d(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	@test isapprox([rotate2d(1, 1, pi/2)...], [(1.,-1.)...])
	true
end

rotate2d(:test)


text:
"""
    Collision

    abstract type defining interaction between two bodies
"""
abstract type Collision{T,O,I,OI} end

# contact point origin
function contact_point_origin(x, q, k)
    x + vector_rotate(k, q)
end

∂contact_point_origin∂x(x, q, k) = 1.0 * I(3)
∂contact_point_origin∂q(x, q, k) = ∂vector_rotate∂q(k, q)
∂contact_point_origin∂k(x, q, k) = rotation_matrix(q)

"""
    contact_normal(collision, xp, qp, xc, qc)

    the contact normal (from child to parent) between two contact points

    collision: Collision
    xp: parent body position
    qp: parent body orientation
    xc: child body position
    qc: child body orientation
"""
function contact_normal(collision::Collision, xp, qp, xc, qc)
    # contact points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # normalized direction
    if dis >= 0.0
        return normalize(dir)'
    else
        return -1.0 * normalize(dir)'
    end
end

function ∂contact_normal_transpose∂x(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    # Jacobians
    X = ∂normalize∂x(dir) * (∂contact_point∂x(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂x(:child, jacobian, collision, xp, qp, xc, qc))

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # normalized direction
    if dis >= 0.0
        return X
    else
        return -1.0 * X
    end
end

function ∂contact_normal_transpose∂q(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    Q = ∂normalize∂x(dir) * (∂contact_point∂q(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂q(:child, jacobian, collision, xp, qp, xc, qc))

    # Jacobians
    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_normal(collision, xp, Quaternion(q..., false), xc, qc)', vector(qp))
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_normal(collision, xp, qp, xc, Quaternion(q..., false))', vector(qc))
    # end

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # @assert norm((dis >= 0.0 ? 1.0 : -1.0) * Q - FD, Inf) < 1.0e-2

    # normalized direction
    if dis >= 0.0
        return Q
    else
        return -1.0 * Q
    end
end

function contact_tangent_one(collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)

    # tangents
    w = [1.0; 0.0; 0.0]
    t1 = skew(w) * n' # tangent
    if !(norm(t1) > 1.0e-6)
        # @warn "edge case!"
        w = [0.0; 1.0; 0.0]
        t1 = skew(w) * n' # tangent
    end

    return t1'
end

function contact_tangent_two(collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)

    # tangent one
    t1 = contact_tangent_one(collision, xp, qp, xc, qc)

    # tangent two
    t2 = skew(t1') * n'

    return t2'
end

"""
    contact_tangent(collision, xp, qp, xc, qc)

    contact tangents between two contact points

    collision: Collision
    xp: parent body position
    qp: parent body orientation
    xc: child body position
    qc: child body orientation
"""
function contact_tangent(collision::Collision, xp, qp, xc, qc)
   [
       contact_tangent_one(collision, xp, qp, xc, qc);
       contact_tangent_two(collision, xp, qp, xc, qc);
   ]
end

function ∂contact_tangent_one_transpose∂x(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)

    # tangents
    w = [1.0; 0.0; 0.0] # candidate
    t1 = skew(w) * n' # tangent
    if !(norm(t1) > 1.0e-5)
        w = [0.0; 1.0; 0.0] # candidate
        t1 = skew(w) * n' # tangent
    end

    ∂t1∂x = skew(w) * ∂contact_normal_transpose∂x(jacobian, collision, xp, qp, xc, qc)

    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, x, qp, xc, qc)[1, :]', xp)
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, xp, qp, x, qc)[1, :]', xc)
    # end

    # @assert norm(FD - ∂t1∂x, Inf) < 1.0e-2

    return ∂t1∂x
end

function ∂contact_tangent_two_transpose∂x(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)
    t1 = contact_tangent_one(collision, xp, qp, xc, qc)

    ∂t1∂x = ∂contact_tangent_one_transpose∂x(jacobian, collision, xp, qp, xc, qc)
    ∂t2∂x = skew(t1') * ∂contact_normal_transpose∂x(jacobian, collision, xp, qp, xc, qc) + ∂skew∂p(n') * ∂t1∂x

    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, x, qp, xc, qc)[2, :]', xp)
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, xp, qp, x, qc)[2, :]', xc)
    # end

    # @assert norm(FD - ∂t2∂x, Inf) < 1.0e-2

    return ∂t2∂x
end

function ∂contact_tangent_one_transpose∂q(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)

    # tangents
    w = [1.0; 0.0; 0.0] # candidate
    t1 = skew(w) * n' # tangent
    if !(norm(t1) > 1.0e-5)
        w = [0.0; 1.0; 0.0] # candidate
        t1 = skew(w) * n' # tangent
    end

    ∂t1∂q = skew(t1) * ∂contact_normal_transpose∂q(jacobian, collision, xp, qp, xc, qc)

    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xp, Quaternion(q..., false), xc, qc)[1, :]', vector(qp))
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xp, qp, xc, Quaternion(q..., false))[1, :]', vector(qc))
    # end

    # return FD
    # @assert norm(FD - ∂t1∂q, Inf) < 1.0e-2

    return ∂t1∂q
end

function ∂contact_tangent_two_transpose∂q(jacobian::Symbol, collision::Collision, xp, qp, xc, qc)
    # normal
    n = contact_normal(collision, xp, qp, xc, qc)

    # tangents
    w = [1.0; 0.0; 0.0] # candidate
    t1 = skew(w) * n' # tangent
    if !(norm(t1) > 1.0e-5)
        w = [0.0; 1.0; 0.0] # candidate
        t1 = skew(w) * n' # tangent
    end
    t2 = skew(t1) * n' # tangent

    ∂t1∂q = skew(w) * ∂contact_normal_transpose∂q(jacobian, collision, xp, qp, xc, qc)

    ∂t2∂q = skew(t1) * ∂contact_normal_transpose∂q(jacobian, collision, xp, qp, xc, qc) + ∂skew∂p(n') * ∂t1∂q

    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xp, Quaternion(q..., false), xc, qc)[2, :]', vector(qp))
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xp, qp, xc, Quaternion(q..., false))[2, :]', vector(qc))
    # end

    # @assert norm(FD - ∂t2∂q, Inf) < 1.0e-2

    return ∂t2∂q
end


text:
# using DataFrames
# using Gadfly
# using Proj4
# using Compose
# using ChoroplethMaps

function graticule(;longs::Union{Array, Range}=[],
                   lats::Union{Array, Range}=[],
                   src_srs="epsg:4326",
                   dst_srs="epsg:3857",
                   showlabels=true,
                   latoflonglabels=Void(),
                   longoflatlabels=Void())
  xs = sort(collect(longs))
  ys = reverse(sort(collect(lats))) # reverse is only to make debugging easier
  src_grid = [[x, y] for y in ys, x in xs]
  dst_grid = mapgrid(src_grid, src_srs, dst_srs)

  vlineary = [line(dst_grid[:,col]) for col in 1:size(dst_grid,2)]
  hlineary = [line(vec(dst_grid[row,:])) for row in 1:size(dst_grid,1)]
  lines = compose(context(), vlineary..., hlineary...)

  layers = [lines]

  if showlabels
    src_proj = Projection(projstr(src_srs))
    dst_proj = Projection(projstr(dst_srs))

    # Place labels at outside left and bottom edges if label pos was unspecified
    labelx = isa(longoflatlabels, Void) ? xs[1] : longoflatlabels
    labely = isa(latoflonglabels, Void) ? ys[end] : latoflonglabels

    texts = vcat(longlabels(xs, labely, src_proj, dst_proj),
                 latlabels(ys, labelx, src_proj, dst_proj))
    labels = compose(context(), texts..., fill(colorant"gray"), fontsize(9pt))
    push!(layers, labels)
  end

  return compose(layers..., stroke(colorant"gray"), linewidth(0.1mm), strokedash([0.5mm, 0.5mm]))
end

function graticule(df::DataFrame; src_srs="epsg:4326", other...)
  # Notice src_srs refers to the *grid*, not the DataFrame!
  dst_srs = df[1, :CM_P]
  lons, lats = ticksfrombbox(
                  mapbbox(
                    bbox(df),
                    dst_srs,
                    src_srs))
  labelx = lons[Integer(floor(size(lons,1) / 4))] # Left edge of original bbox
  labely = lats[Integer(ceil(size(lats,1) / 4))] # Bottom edge of original bbox
  return graticule(longs=lons, lats=lats, src_srs=src_srs, dst_srs=dst_srs,
                   latoflonglabels=labely, longoflatlabels=labelx; other...)
end

function longlabels(xs, y, src_proj, dst_proj)
  map(xs) do x
    px, py = transform(src_proj, dst_proj, [x, y])
    text(px, py, "$x", hcenter, vtop)
  end
end

function latlabels(ys, x, src_proj, dst_proj)
  map(ys) do y
    px, py = transform(src_proj, dst_proj, [x, y])
    text(px, py, "$y", hleft, vcenter)
  end
end

# function testgrid(prj)
#   Gadfly.set_default_graphic_size(25cm, 19cm)
#   df = DataFrame(NAME=["Arkansas", "Florida", "Mississippi", "Alabama", "Louisiana", "Tennessee", "Georgia", "South Carolina"], feature=[8.7,13.2,36.7,30.0,25.6,23.7,24.6,14.6])
#   mp = mapify(df, Provider.STATESUMMARY(), key=:NAME, projection=prj)

#   guido = graticule(mp)
#   guida = graticule(longs=-91.0:0.5:-89.0,
#                     lats=23.0:0.5:25,
#                     src_srs="epsg:4236", dst_srs=prj,
#                     latoflonglabels=22.5, longoflatlabels=-88.5)

#   choroplethmap(mp, group=:NAME, color=:feature, Guide.annotation(guido), #Guide.annotation(guida),
#                 Guide.xticks( label=false), Guide.yticks( label=false),  #removing ticks also removes zoom
#                 Guide.XLabel(""), Guide.YLabel(""), Theme(grid_color=colorant"rgba(0,0,0,0.0)"))
# end

function grid(df::DataFrame, gridcrs::AbstractString)
  datacrs = df[1, :CM_P]
  unprojgrid = gridfrombbox(
                 mapbbox(
                   bbox(df),
                   datacrs,
                   gridcrs))
  return (unprojgrid, mapgrid( unprojgrid, gridcrs, datacrs))
end

function bbox(df::DataFrame)
  return Array[[extrema(df[:CM_X])...], [extrema(df[:CM_Y])...]]
end

function mapgrid(grid, from::AbstractString, to::AbstractString)
  projfrom = Projection(projstr(from))
  projto = Projection(projstr(to))
  return map((el) -> (transform(projfrom, projto, el)...), grid)
end

function gridfrombbox(bbox)
  xs = gridondim(bbox, 1, 180)
  ys = gridondim(bbox, 2, 90)
  return [[x, y] for y in reverse(ys), x in xs]
end

function ticksfrombbox(bbox)
  xs = gridondim(bbox, 1, 180)
  ys = gridondim(bbox, 2, 90)
  return (xs, ys)
end

function gridondim(bbox, dim, bmax)
  rng = getrange(bbox, dim)
  allowed_steps = Any[0.1, 0.2, 0.5, 1, 2, 5, 10, 15, 20, 30]
  step = allowed_steps[findfirst((s) -> rng / s <= 12, allowed_steps)]

  # expand the range by 50% in each direction to ensure graticule covers entire map
  c = mean(bbox[dim])
  s = lock(c - 1.5 * rng / 2, bmax)
  t = lock(c + 1.5 * rng / 2, bmax)

  low  = nudge(s, step, floor)
  high = nudge(t, step, ceil)
  return collect(low:step:high)
end

function lock(x, bmax)
  abs(x) > bmax && return sign(x) * bmax
  return x
end

function getrange(bbox, dim)
  return bbox[dim][2] - bbox[dim][1]
end

function nudge(x, step, func)
  return step * Integer(func(x / step))
end

# Map each of the 4 bbox corners into the new coordinate space and return
# a new bounding box in that cs.
function mapbbox(bbox, from::AbstractString, to::AbstractString)
  tfp = transform(Projection(projstr(from)), Projection(projstr(to)),
                  [[bbox[1][1] bbox[2][1]]; [bbox[1][2] bbox[2][1]];
                   [bbox[1][2] bbox[2][2]]; [bbox[1][1] bbox[2][2]]])

  return Array[[extrema(tfp[:, 1])...], [extrema(tfp[:, 2])...]]
end


text:
module BasisModule
export GaussianBasis, computeBasis, computeBasisShells

using ..BaseModule
using ..BasisFunctionsModule
using ..BasisSetModule
using ..GeometryModule
using ..AtomModule
using ..IntegralsModule
using ..ShellModule

import Base.display
import Base.convert

abstract type Basis end

"""
A GaussianBasis is a Basis constructed purely from Gaussian type functions (see PrimitiveGaussianBasisFunction).
The list of basis functions defines the basis. It is most easily constructed from a BasisSet and a Geometry by computeBasis().
"""
immutable GaussianBasis <: Basis	# doesn't change much to the user as contractedBFs is mutable anyways
  contractedBFs::Vector{ContractedGaussianBasisFunction}
end

function computeBasis_primitive(basSet::BasisSet,geo::Geometry)
  bas=GaussianBasis([])
  for atom in geo.atoms
    contrDefinitions = basSet.definitions[atom.element]
    for cdef in contrDefinitions
      for mqn in MQuantumNumbers(cdef.lQuantumNumber)
	coeffs = [pdef.prefactor for pdef in cdef.primitives]
	pgbfs  = [PrimitiveGaussianBasisFunction(atom.position,pdef.exponent,mqn) for pdef in cdef.primitives]
	push!(bas.contractedBFs,ContractedGaussianBasisFunction(coeffs,pgbfs))
      end
    end
  end
  return bas
end

function renormprimitives!(basis::GaussianBasis)
  for cgbf in basis.contractedBFs
    for p in eachindex(cgbf.primitiveBFs)
      pgbf = cgbf.primitiveBFs[p]
      primnorm = computeIntegralOverlap(pgbf,pgbf)
      cgbf.coefficients[p] /= sqrt(primnorm)
    end
  end
end

function computeBasis(basSet::BasisSet,geo::Geometry,normalize::Bool=true)
  bas = computeBasis_primitive(basSet,geo)
  if normalize
    renormprimitives!(bas)
    normalize!.(bas.contractedBFs)
  end
  return bas
end

function display(basis::GaussianBasis)
  println(typeof(basis))
  indent = "  "

  println(indent * "contractedBFs: ")
  indent = indent * "  "

  for cgbf in basis.contractedBFs
    display(cgbf,indent)
  end
end

function computeBasisShells(basSet::BasisSet,geo::Geometry)
  shells = Shell[]
  for atom in geo.atoms
    for contractedDefinition in basSet.definitions[atom.element]
      exponents = [prim.exponent for prim in contractedDefinition.primitives]
      coefficients = [prim.prefactor for prim in contractedDefinition.primitives]
      sh = Shell(atom.position,contractedDefinition.lQuantumNumber,exponents,coefficients)
      push!(shells,sh)
    end
  end
  return shells
end

convert(::Type{GaussianBasis},sh::Shell) = GaussianBasis(expandShell(sh::Shell))

function convert(::Type{GaussianBasis},shells::Vector{Shell})
  bas = GaussianBasis([])
  for sh in shells
    append!(bas.contractedBFs,expandShell(sh))
  end
  return bas
end

end # module


text:
"""
    isCartesian(sim)
"""
CoordinateSystems.isCartesian(sim::Simulation) = isCartesian(sim.dis)
CoordinateSystems.isCartesian(dis::Discretization) = isCartesian(dis.coordinate_system)

"""
    isPolar(sim)
"""
CoordinateSystems.isPolar(sim::Simulation) = isPolar(sim.dis)
CoordinateSystems.isPolar(dis::Discretization) = isPolar(dis.coordinate_system)


text:
@ch_struct struct ClientInfo
    query_kind::UInt8
    initial_user::String
    initial_query_id::String
    initial_address_string::String
    read_interface::UInt8
    os_user::String
    client_hostname::String
    client_name::String
    client_ver_major::VarUInt
    client_ver_minor::VarUInt
    client_rev::VarUInt
    @has_quota_key quota_key::String = ""
    @has_version_patch client_ver_patch::VarUInt = VarUInt(0)
end
ClientInfo() = ClientInfo(
    0x01,
    "",
    "",
    "0.0.0.0:0",
    0x01,
    "",
    "",
    CLIENT_NAME,
    DBMS_VER_MAJOR,
    DBMS_VER_MINOR,
    DBMS_VER_REV,
    "",
    2,
)

@ch_struct struct ClientHello
    client_name::String
    client_dbms_ver_major::VarUInt
    client_dbms_ver_minor::VarUInt
    client_dbms_ver_rev::VarUInt
    database::String
    username::String
    password::String
end

@ch_struct struct ClientPing
end

@ch_struct struct ClientQuery
    query_id::String
    @has_client_info client_info::ClientInfo = ClientInfo()
    settings::String
    query_stage::VarUInt
    compression::VarUInt
    query::String
end


text:
using Test

include("../../src/io.jl")

@testset "ChunkedInputBuffer Tests" begin

    @testset "test_should_start_empty" begin
        buffer = ChunkedInputBuffer()
        @test length(view(buffer)) == 0
    end

    @testset "test_should_be_able_to_set_capacity" begin
        buffer = ChunkedInputBuffer(capacity=10)
        @test capacity(buffer) == 10
    end

    @testset "test_should_be_able_to_load_data" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello")
        @test view(buffer) == b"\x00\x05hello"
    end

    @testset "test_should_be_able_to_load_multiple_times" begin
        buffer = ChunkedInputBuffer()

        load(buffer, b"\x00\x05hello")
        load(buffer, b"\x00\x05world")

        @test view(buffer) == b"\x00\x05hello\x00\x05world"
    end

    @testset "test_should_be_able_to_load_after_discard" begin
        buffer = ChunkedInputBuffer()

        load(buffer, b"\x00\x05hello\x00\x00")
        frame_message(buffer)
        frame = buffer.frame

        @test frame.panes == [(2,7)]

        discard_message(buffer)
        load(buffer, b"\x00\x07bonjour\x00\x00")
        frame_message(buffer)
        frame = buffer.frame

        @test frame.panes == [(2, 9)]
    end

    @testset "test_should_auto_extend_on_load" begin
        buffer = ChunkedInputBuffer(capacity=10)
        load(buffer, b"\x00\x05hello")

        load(buffer, b"\x00\x07bonjour")

        @test capacity(buffer) == 16
    end

    @testset "test_should_start_with_no_frame" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")

        @test buffer.frame === nothing
    end

    @testset "test_should_be_able_to_frame_message" begin
        buffer = ChunkedInputBuffer(capacity=10)
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")

        framed = frame_message(buffer)

        @test framed
        @test buffer.frame.panes == [(2, 7)]
    end

    @testset "test_should_be_able_to_frame_empty_message" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x00")

        framed = frame_message(buffer)

        @test framed
        @test buffer.frame.panes == []
    end

    @testset "test_should_not_be_able_to_frame_empty_buffer" begin
        buffer = ChunkedInputBuffer()

        framed = frame_message(buffer)

        @test !framed
    end

    @testset "test_should_not_be_able_to_frame_partial_message" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello")

        framed = frame_message(buffer)

        @test !framed
    end

    @testset "test_should_be_able_to_discard_message" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")
        frame_message(buffer)

        discard_message(buffer)

        @test buffer.frame === nothing
    end

    @testset "test_should_be_able_to_discard_empty_message" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x00")
        frame_message(buffer)

        discard_message(buffer)

        @test buffer.frame === nothing
    end

    @testset "test_discarding_message_should_move_read_pointer" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")
        frame_message(buffer)
        discard_message(buffer)

        framed = frame_message(buffer)

        @test framed
        @test buffer.frame.panes == [(2,9)]
    end

    @testset "test_should_be_able_to_frame_successive_messages_without_discarding" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")
        frame_message(buffer)

        framed = frame_message(buffer)

        @test framed
        @test buffer.frame.panes == [(2,9)]
    end

    @testset "test_discarding_message_should_not_recycle_buffer" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")
        load(buffer, b"\x00\x07bonjour\x00\x00")
        frame_message(buffer)

        discard_message(buffer)

        @test view(buffer) == b"\x00\x05hello\x00\x00\x00\x07bonjour\x00\x00"
    end

    @testset "test_should_not_be_able_to_frame_message_if_empty" begin
        buffer = ChunkedInputBuffer()

        @test !frame_message(buffer)
    end

    @testset "test_should_not_be_able_to_frame_message_if_incomplete" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello")

        @test !frame_message(buffer)
    end

    @testset "test_should_be_able_to_frame_message_if_complete" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00")

        @test frame_message(buffer)
    end

    @testset "test_should_be_able_to_frame_message_if_complete_with_more" begin
        buffer = ChunkedInputBuffer()
        load(buffer, b"\x00\x05hello\x00\x00\x00\x05world\x00\x00")

        @test frame_message(buffer)
    end
end

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
#!/usr/bin/env julia

using ArgParse
using GZip
using MIToS.Utils # get_n_words, check_file

function parse_commandline()
    s = ArgParseSettings(description = "Splits a file with multiple sequence alignments in Stockholm format, creating one compressed file per MSA in Stockholm format: accessionumber.gz",
                         version = "MIToS $(Pkg.installed("MIToS"))",
                         add_version = true)

    @add_arg_table s begin
        "file"
            help = "Input file"
            required = true
        "--path", "-p"
        		help = "Path for the output files [default: execution directory]"
		        arg_type = ASCIIString
    		    default = ""
    end

    s.epilog = """

    \n
    MIToS $(Pkg.installed("MIToS"))\n
    \n
    Bioinformatics Unit\n
    Leloir Institute Foundation\n
    Av. Patricias Argentinas 435, CP C1405BWE, Buenos Aires, Argentina
    """

    return parse_args(s)
end

const Args = parse_commandline()

function main(input)
	infh = GZip.open(input)
	lines = []
	id = "no_accessionumber"
	for line in eachline(infh)
		if length(line) > 7 && line[1:7] == "#=GF AC"
			id = get_n_words(line, 3)[3]
		end
		push!(lines, line)
		if line == "//\n"
			filename = joinpath(Args["path"], string(id, ".gz"))
			outfh = GZip.open(filename, "w")
			for l in lines
				write(outfh, string(l))
			end
			close(outfh)
			id = "no_accessionumber"
			empty!(lines)
		end
	end
	close(infh)
end

main(check_file(Args["file"]))


text:
# using Revise
using Test
using JDF, DataFrames

@testset "JDF.jl Char" begin
    ac = ['a', 'b', 'a', 'c', 'd']
    metadata_ac = open("io.jdf", "w") do io
        compress_then_write(ac, io)
    end

    a = ['a', 'b', 'a', missing, 'd']
    metadata = open("io.jdf", "w") do io
        compress_then_write(a, io)
    end

    aa = open("io.jdf", "r") do io
        column_loader(Union{Missing,Symbol}, io, metadata)
    end

    df = DataFrame(a = a, ac = ac)
    savejdf("pls_del.jdf", df)
    df2 = loadjdf("pls_del.jdf")

    @test size(df) == size(df2)
    @test all(isequal.(df.a, df2.a))
    @test all(isequal.(df.ac, df2.ac))

    rm("io.jdf", force = true)
    rm("pls_del.jdf", force = true, recursive = true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using Test

using Persa
using Surprise
using DatasetsCF

dataset = DatasetsCF.MovieLens()

@testset "Memory Based Algorithms" begin
    model = Surprise.KNNBasic(dataset)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0

    model = Surprise.KNNBaseline(dataset)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0

    model = Surprise.KNNWithMeans(dataset)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0

    model = Surprise.SlopeOne(dataset)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0
end

@testset "Model Based Algorithms" begin
    model = Surprise.IRSVD(dataset, n_epochs = 2)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0

    model = Surprise.RSVD(dataset, n_epochs = 2)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0

    model = Surprise.NMF(dataset, n_epochs = 2)
    Persa.train!(model, dataset)
    @test model[1,1] >= 0
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# All types and constructors are defined here.

## Axis configuration
mutable struct Axes
    axesconf::Dict{Symbol, Any}
end

Axes(pairs::Pair...) = Axes(Dict(pairs))
Axes(;args...) = Axes(args)
keys(a::Axes) = keys(a.axesconf)
getindex(a::Axes, args... ; kwargs...) = getindex(a.axesconf, args... ; kwargs...)
merge(a::Axes, b::Axes) = Axes(merge(a.axesconf, b.axesconf))

## Real coordinates
const NCoord = Union{AbstractRange{T}, AbstractArray{T}} where T # actual coordinates
const Coord = Union{NCoord, Nothing}  # nullable coordinates

## Complex coordinates
const ComplexCoord = AbstractArray{T} where T <: Complex

## Curves

# A curve is a set of coordinates and a plot string.
struct Curve{X, Y, Z, S}
    x::X
    y::Y
    z::Z
    supp::S
    conf::String
end

## Plots
# A plot has a datafile, dims, a configuration, and a set of curves.
Base.@kwdef mutable struct Plot
    datafile::String = tempname()   # file to store plot data
    dims::Int        = 2            # 2-D or 3-D plot
    axesconf::String = ""           # axes configuration
    curves::Vector{Curve} = Curve[] # a vector of curves
end
getindex(p::Plot, args... ; kwargs...) = getindex(p.curves, args... ; kwargs...)
isempty(p::Plot) = isempty(p.curves)

# A subplot can be a Plot, or can be nothing.
const SubPlot = Union{Plot, Nothing}
const SubPlots = Vector{SubPlot}

# A layout is a tuple of integers.
const Layout = Tuple{Int, Int}

## Figure
# At the top level, a figure is a handle, a layout and a set of Axes.
Base.@kwdef mutable struct Figure
    handle::Int                   # each figure has a unique handle
    layout::Layout     = (0,0)    # multiplot layout
    subplots::SubPlots = [Plot()] # a figure contains a number of plots
end
Figure(handle::Int) = Figure(handle=handle)

getindex(f::Figure, args... ; kwargs...) = getindex(f.subplots, args... ; kwargs...)

isempty(f::Figure; subplot = 1) = isempty(f[subplot])

length(f::Figure) = length(f.subplots)

function push!(f::Figure, c::Curve ; subplot = 1)
    if isempty(f, subplot=subplot)
        f[subplot].curves = [c]
    else
        push!(f[subplot].curves, c)
    end
end

function push!(f::Figure, sp::SubPlot)
    if isempty(f)
        f.subplots = [sp]
    else
        push!(f.subplots, sp)
    end
end

function push!(dest::Figure, src::Figure)
    for sp in src.subplots
        push!(dest, sp)
    end
end

# Types needed for multiplot
const FigArray = Union{Array{Figure}, Array{Union{Figure, Nothing}}}

# We need a global variable to keep track of gnuplot's state
mutable struct GnuplotState
    current                      # current figure -- "pointer" to figs
    gnuplot_available::Bool      # is gnuplot installed on this system?
    gp_stderr::String            # last data read from  gnuplot's stderr
    gp_lasterror::String         # gnuplot's last error output
    gp_error::Bool               # true if last command resulted in gp error
    figs::Vector{Figure}         # storage for all figures
end

GnuplotState() = GnuplotState(nothing,false,"","",false,Figure[])


text:
# Amino Acid
# ==========
#
# Amino acid type.
#
# This file is a part of the BioJulia ecosystem.
# License is MIT: https://github.com/BioJulia/NucleicAcids.jl/blob/master/LICENSE.md

"""
An amino acid type.
"""
primitive type AminoAcid <: BioSymbol 8 end

prefix(::AminoAcid) = "AA"
type_text(::AminoAcid) = "Amino Acid"
isterm(symbol::AminoAcid) = symbol == AA_Term
bytemask(symbol::AminoAcid) = 0b11111

# Conversion from/to integers
# ---------------------------

Base.convert(::Type{AminoAcid}, aa::UInt8) = reinterpret(AminoAcid, aa)
Base.convert(::Type{UInt8}, aa::AminoAcid) = reinterpret(UInt8, aa)
Base.convert(::Type{T}, aa::AminoAcid) where T <: Number = convert(T, convert(UInt8, aa))
Base.convert(::Type{AminoAcid}, aa::T) where T <: Number = convert(AminoAcid, convert(UInt8, aa))
AminoAcid(nt::Integer) = convert(AminoAcid, nt)

# Conversion from/to Char
# -----------------------

function Base.convert(::Type{AminoAcid}, c::Char)
    aa = tryparse(AminoAcid, c)
    if aa == nothing
        throw(InexactError(:convert, AminoAcid, c))
    end
    return aa
end
AminoAcid(c::Char) = convert(AminoAcid, c)

Base.convert(::Type{Char}, aa::AminoAcid) = aa_to_char[convert(UInt8, aa) + 1]
Char(aa::AminoAcid) = convert(Char, aa)

# Amino acid encoding definition
# ------------------------------

"Invalid Amino Acid"
const AA_INVALID = convert(AminoAcid, 0x1c)  # Used during conversion from strings

# lookup table for characters
const char_to_aa = [AA_INVALID for _ in 0x00:0x7f]
const aa_to_char = Vector{Char}(undef, 0x1c)

# compatibility bits
const compatbits_aa = Vector{UInt32}(undef, 28)

# This set of amino acids is defined by IUPAC-IUB Joint Commission on Biochemical Nomenclature.
# Reference: http://www.insdc.org/documents/feature_table.html#7.4.3

for (aa, doc, code) in [
        ('A', "Alanine",                           0x00),
        ('R', "Arginine",                          0x01),
        ('N', "Asparagine",                        0x02),
        ('D', "Aspartic Acid",                     0x03),
        ('C', "Cysteine",                          0x04),
        ('Q', "Glutamine",                         0x05),
        ('E', "Glutamic Acid",                     0x06),
        ('G', "Glycine",                           0x07),
        ('H', "Histidine",                         0x08),
        ('I', "Isoleucine",                        0x09),
        ('L', "Leucine",                           0x0a),
        ('K', "Lysine",                            0x0b),
        ('M', "Methionine",                        0x0c),
        ('F', "Phenylalanine",                     0x0d),
        ('P', "Proline",                           0x0e),
        ('S', "Serine",                            0x0f),
        ('T', "Threonine",                         0x10),
        ('W', "Tryptophan",                        0x11),
        ('Y', "Tyrosine",                          0x12),
        ('V', "Valine",                            0x13),
        ('O', "Pyrrolysine",                       0x14),  # non-standard
        ('U', "Selenocysteine",                    0x15),  # non-standard
        ('B', "Aspartic Acid or Asparagine",       0x16),  # ambiguous
        ('J', "Leucine or Isoleucine",             0x17),  # ambiguous
        ('Z', "Glutamine or Glutamic Acid",        0x18),  # ambiguous
        ('X', "Unspecified or Unknown Amino Acid", 0x19)]  # ambiguous
    var = Symbol("AA_", aa)
    @eval begin
        @doc $doc const $var = convert(AminoAcid, $code)
        char_to_aa[$(Int(aa)+1)] = char_to_aa[$(Int(lowercase(aa))+1)] = $var
        aa_to_char[$(code)+1] = $aa
    end
    if code ≤ 0x15
        compatbits_aa[code+1] = 1 << code
    elseif code == 0x16
        compatbits_aa[code+1] = 1 << 0x02 | 1 << 0x03
    elseif code == 0x17
        compatbits_aa[code+1] = 1 << 0x09 | 1 << 0x0a
    elseif code == 0x18
        compatbits_aa[code+1] = 1 << 0x05 | 1 << 0x06
    elseif code == 0x19
        compatbits_aa[code+1] = (1 << 0x16) - 1
    end
end

"Terminal"
const AA_Term = convert(AminoAcid, 0x1a)
char_to_aa[Int('*')+1] = AA_Term
aa_to_char[0x1a+1] = '*'
compatbits_aa[0x1a+1] = 1 << 0x1a

"Amino Acid Gap"
const AA_Gap = convert(AminoAcid, 0x1b)
char_to_aa[Int('-') + 1] = AA_Gap
aa_to_char[0x1b+1] = '-'
compatbits_aa[0x1b+1] = 0

@eval alphabet(::Type{AminoAcid}) = $(tuple([reinterpret(AminoAcid, x) for x in 0x00:0x1b]...))

"""
    alphabet(AminoAcid)

Get all symbols of `AminoAcid` in sorted order.

Examples
--------

```jldoctest
julia> alphabet(AminoAcid)
(AA_A, AA_R, AA_N, AA_D, AA_C, AA_Q, AA_E, AA_G, AA_H, AA_I, AA_L, AA_K, AA_M, AA_F, AA_P, AA_S, AA_T, AA_W, AA_Y, AA_V, AA_O, AA_U, AA_B, AA_J, AA_Z, AA_X, AA_Term, AA_Gap)

julia> issorted(alphabet(AminoAcid))
true

```
"""
alphabet(::Type{AminoAcid})

# lookup table of 20 standard amino acids
const threeletter_to_aa = Dict(
    "ALA" => AA_A, "ARG" => AA_R, "ASN" => AA_N, "ASP" => AA_D, "CYS" => AA_C,
    "GLN" => AA_Q, "GLU" => AA_E, "GLY" => AA_G, "HIS" => AA_H, "ILE" => AA_I,
    "LEU" => AA_L, "LYS" => AA_K, "MET" => AA_M, "PHE" => AA_F, "PRO" => AA_P,
    "SER" => AA_S, "THR" => AA_T, "TRP" => AA_W, "TYR" => AA_Y, "VAL" => AA_V,
    "ASX" => AA_B, "XLE" => AA_J, "GLX" => AA_Z, "XAA" => AA_X,
    "PYL" => AA_O, "SEC" => AA_U,
)

# Generate an amino acid parser.
let
    re = Automa.RegExp
    function aapat(three, aa)
        one = convert(Char, aa)
        pat = re.alt(
            re.cat([re.alt(lowercase(x), uppercase(x)) for x in three]...),
            lowercase(one),
            uppercase(one))
        pat.actions[:exit] = [Symbol(three)]
        return pat
    end
    aminoacids = [
        ("ALA", AA_A),
        ("ARG", AA_R),
        ("ASN", AA_N),
        ("ASP", AA_D),
        ("CYS", AA_C),
        ("GLN", AA_Q),
        ("GLU", AA_E),
        ("GLY", AA_G),
        ("HIS", AA_H),
        ("ILE", AA_I),
        ("LEU", AA_L),
        ("LYS", AA_K),
        ("MET", AA_M),
        ("PHE", AA_F),
        ("PRO", AA_P),
        ("SER", AA_S),
        ("THR", AA_T),
        ("TRP", AA_W),
        ("TYR", AA_Y),
        ("VAL", AA_V),
        ("ASX", AA_B),
        ("XLE", AA_J),
        ("GLX", AA_Z),
        ("XAA", AA_X),
        ("PYL", AA_O),
        ("SEC", AA_U),
    ]
    aa_term = re"\*"
    aa_term.actions[:exit] = [:Term]
    aa_gap = re"-"
    aa_gap.actions[:exit] = [:Gap]
    whitespaces = re"[ \t\r\n]*"
    aminoacid = re.cat(
        whitespaces,
        re.alt(vcat((aapat(three, aa) for (three, aa) in aminoacids)..., aa_term, aa_gap)...),
        whitespaces)
    machine = Automa.compile(aminoacid)
    actions = Dict(Symbol(three) => :(aa = $(aa)) for (three, aa) in aminoacids)
    actions[:Term] = :(aa = AA_Term)
    actions[:Gap]  = :(aa = AA_Gap)

    ctx = Automa.CodeGenContext(checkbounds=false)
    @eval function Base.tryparse(::Type{AminoAcid}, data::AbstractString)
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = sizeof(data)
        aa = AA_INVALID
        $(Automa.generate_exec_code(ctx, machine, actions))
        if cs != 0
            return nothing
        end
        return aa
    end
end

function Base.tryparse(::Type{AminoAcid}, c::Char)
    @inbounds aa = c <= '\x7f' ? char_to_aa[Int(c)+1] : AA_INVALID
    if aa == AA_INVALID
        return nothing
    else
        return aa
    end
end

function Base.parse(::Type{AminoAcid}, c::Union{AbstractString,Char})
    aa = tryparse(AminoAcid, c)
    if aa == nothing
        throw(ArgumentError("invalid amino acid"))
    end
    return aa
end

# Arithmetic and Order
# --------------------

# These methods are necessary when deriving some algorithims
# like iteration, sort, comparison, and so on.

# 0x1c is the size of the amino acid alphabet
Base.:-(x::AminoAcid, y::Integer) = x + mod(-y, 0x1c)
Base.:+(x::AminoAcid, y::Integer) = reinterpret(AminoAcid, mod((convert(UInt8, x) + y) % UInt8, 0x1c))

Base.isvalid(::Type{AminoAcid}, x::Integer) = 0 ≤ x ≤ 0x1b
Base.isvalid(aa::AminoAcid) = aa ≤ AA_Gap

"""
    isambiguous(aa::AminoAcid)

Test if `aa` is an ambiguous amino acid.
"""
isambiguous(aa::AminoAcid) = AA_B ≤ aa ≤ AA_X

"""
    iscertain(aa::AminoAcid)

Test if `aa` is a non-ambiguous amino acid.
"""
function iscertain(aa::AminoAcid)
    return AA_A ≤ aa ≤ AA_U || aa == AA_Term
end

"""
    gap(AminoAcid)

Return `AA_Gap`.
"""
gap(::Type{AminoAcid}) = AA_Gap

"""
    compatbits(aa::AminoAcid)

Return the compatibility bits of `aa` as `UInt32`.

Examples
--------

```jldoctest
julia> compatbits(AA_A)
0x00000001

julia> compatbits(AA_J)
0x00000600

```
"""
compatbits(aa::AminoAcid) = @inbounds compatbits_aa[reinterpret(UInt8, aa) + 1]


text:
export
    pckcls,
    pckcov!,
    pckcov,
    pckfrm!,
    pckfrm,
    pcklof,
    pckopn,
    pckuof,
    pckw02,
    pcpool,
    pdpool,
    pgrrec,
    phaseq,
    pipool,
    pjelpl,
    pl2nvc,
    pl2nvp,
    pl2psv,
    pltar,
    pltexp,
    pltnp,
    pltnrm,
    pltvol,
    polyds,
    pos,
    posr,
    prop2b,
    prsdp,
    prsint,
    psv2pl,
    pxform,
    pxfrm2

"""
    pckcls(handle)

Close an open PCK file.

### Arguments ###

- `handle`: Handle of the PCK file to be closed

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckcls_c.html)
"""
function pckcls(handle)
    ccall((:pckcls_c, libcspice), Cvoid, (SpiceInt,), handle)
    handleerror()
end

"""
    pckcov!(cover, pck, idcode)

Find the coverage window for a specified reference frame in a specified binary PCK file.

### Arguments ###

- `cover`: An initalized window `SpiceDoubleCell`
- `pck`: Path of PCK file
- `idcode`: Class ID code of PCK reference frame

### Output ###

Returns `cover` containing coverage in `pck` for `idcode`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckcov_c.html)
"""
function pckcov!(cover, pck, idcode)
    ccall((:pckcov_c, libcspice), Cvoid, (Cstring, SpiceInt, Ref{Cell{SpiceDouble}}),
          pck, idcode, cover.cell)
    handleerror()
    cover
end

@deprecate pckcov pckcov!

"""
    pckfrm!(ids, pck)

Find the set of reference frame class ID codes of all frames in a specified binary PCK file.

### Arguments ###

- `ids`: An initalized `SpiceIntCell`
- `pck`: Path of PCK file

### Output ###

Returns `ids` containing a set of frame class ID codes of frames in PCK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckfrm_c.html)
"""
function pckfrm!(ids, pck)
    ccall((:pckfrm_c, libcspice), Cvoid, (Cstring, Ref{Cell{SpiceInt}}),
          pck, ids.cell)
    handleerror()
    ids
end

@deprecate pckfrm pckfrm!

"""
    pcklof(filename)

Load a binary PCK file for use by the readers. Return the handle of the loaded
file which is used by other PCK routines to refer to the file.

### Arguments ###

- `filename`: Path of the PCK file

### Output ###

Returns an integer handle.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pcklof_c.html)
"""
function pcklof(filename)
    handle = Ref{SpiceInt}()
    ccall((:pcklof_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}), filename, handle)
    handleerror()
    handle[]
end

"""
    pckopn(name, ifname, ncomch)

Create a new PCK file, returning the handle of the opened file.

### Arguments ###

- `name`: The name of the PCK file to be opened
- `ifname`: The internal filename for the PCK
- `ncomch`: The number of characters to reserve for comments

### Output ###

Returns the handle of the opened PCK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckopn_c.html)
"""
function pckopn(name, ifname, ncomch)
    handle = Ref{SpiceInt}()
    ccall((:pckopn_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceInt, Ref{SpiceInt}),
          name, ifname, ncomch, handle)
    handleerror()
    handle[]
end

"""
    pckuof(handle)

Unload a binary PCK file so that it will no longer be searched by the readers.

### Arguments ###

- `handle`: Integer handle of a PCK file

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckuof_c.html)
"""
function pckuof(handle)
    ccall((:pckuof_c, libcspice), Cvoid, (SpiceInt,), handle)
end

"""
    pckw02(handle, clssid, frame, first, last, segid, intlen, cdata, btime)

Write a type 2 segment to a PCK binary file given the file handle, frame class ID, base frame, time
range covered by the segment, and the Chebyshev polynomial coefficients.

### Arguments ###

- `handle`: Handle of binary PCK file open for writing.
- `clssid`: Frame class ID of body-fixed frame.
- `frame`: Name of base reference frame.
- `first`: Start time of interval covered by segment.
- `last`: End time of interval covered by segment.
- `segid`: Segment identifier.
- `intlen`: Length of time covered by logical record.
- `cdata`: Array of Chebyshev coefficients.
- `btime`: Begin time of first logical record.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pckw02_c.html)
"""
function pckw02(handle, clssid, frame, first, last, segid, intlen, cdata, btime)
    n = length(cdata)
    polydg = length(cdata[1]) - 1
    cdata_ = array_to_cmatrix(cdata, n=polydg + 1)
    ccall((:pckw02_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Cstring, SpiceDouble, SpiceDouble, Cstring, SpiceDouble, SpiceInt,
           SpiceInt, Ref{SpiceDouble}, SpiceDouble),
          handle, clssid, frame, first, last, segid, intlen, n, polydg, cdata_, btime)
    handleerror()
end

"""
    pcpool(name, vals)

Insert character data into the kernel pool.

### Arguments ###

- `name`: The kernel pool name to associate with `vals`
- `vals`: An array of values to insert into the kernel pool

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pcpool_c.html)
"""
function pcpool(name, vals)
    vals_, n, lenvals = chararray(vals)
    ccall((:pcpool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
          name, n, lenvals, vals_)
    handleerror()
end

"""
    pdpool(name, vals)

Insert double precision data into the kernel pool.

### Arguments ###

- `name`: The kernel pool name to associate with `vals`
- `vals`: An array of values to insert into the kernel pool

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pdpool_c.html)
"""
function pdpool(name, vals)
    n = length(vals)
    ccall((:pdpool_c, libcspice), Cvoid, (Cstring, SpiceInt, Ref{SpiceDouble}), name, n, vals)
    handleerror()
end

"""
    pgrrec(body, lon, lat, alt, re, f)

Convert planetographic coordinates to rectangular coordinates.

### Arguments ###

- `body`: Body with which coordinate system is associated.
- `lon`: Planetographic longitude of a point (radians).
- `lat`: Planetographic latitude of a point (radians).
- `alt`: Altitude of a point above reference spheroid.
- `re`: Equatorial radius of the reference spheroid.
- `f`: Flattening coefficient.

### Output ###

Returns the rectangular coordinates of the point.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pgrrec_c.html)
"""
function pgrrec(body, lon, lat, alt, re, f)
    rectan = Array{SpiceDouble}(undef, 3)
    ccall((:pgrrec_c, libcspice), Cvoid,
          (Cstring, SpiceDouble, SpiceDouble, SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          body, lon, lat, alt, re, f, rectan)
    handleerror()
    rectan
end

"""
    phaseq(et, target, illmn, obsrvr, abcorr)

Compute the apparent phase angle for a target, observer, illuminator set of ephemeris objects.

### Arguments ###

- `et`: Ephemeris seconds past J2000 TDB
- `target`: Target body name
- `illmn`: Illuminating body name
- `obsrvr`: Observer body
- `abcorr`: Aberration correction flag

### Output ###

Returns the value of the phase angle.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/phaseq_c.html)
"""
function phaseq(et, target, illmn, obsrvr, abcorr)
    retval = ccall((:phaseq_c, libcspice), SpiceDouble,
                   (SpiceDouble, Cstring, Cstring, Cstring, Cstring),
                   et, target, illmn, obsrvr, abcorr)
    handleerror()
    retval
end

# Deprecated
_pi() = ccall((:pi_c, libcspice), SpiceDouble, ())

"""
    pipool(name, ivals)

Insert integer data into the kernel pool.

### Arguments ###

- `name`: The kernel pool name to associate with the values
- `ivals`: An array of integers to insert into the pool

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pipool_c.html)
"""
function pipool(name, ivals)
    n = length(ivals)
    ivals_ = SpiceInt.(ivals)
    ccall((:pipool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}),
          name, n, ivals_)
    handleerror()
end

"""
    pjelpl(elin, plane)

Project an ellipse onto a plane, orthogonally.

### Arguments ###

- `elin`: An ellipse to be projected
- `plane`: A plane onto which `elin` is to be projected

### Output ###

Returns the ellipse resulting from the projection.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pjelpl_c.html)
"""
function pjelpl(elin, plane)
    elout = Ref{Ellipse}()
    ccall((:pjelpl_c, libcspice), Cvoid,
          (Ref{Ellipse}, Ref{Plane}, Ref{Ellipse}),
          elin, plane, elout)
    handleerror()
    elout[]
end

"""
    pl2nvc(plane)

Return a unit normal vector and constant that define a specified plane.

### Arguments ###

- `plane`: A plane

### Output ###

Returns a tuple consisting of

- `normal`: A normal vector and...
- `constant`: ... constant defining the geometric plane represented by `plane`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pl2nvc_c.html)
"""
function pl2nvc(plane)
    normal = Array{SpiceDouble}(undef, 3)
    constant = Ref{SpiceDouble}()
    ccall((:pl2nvc_c, libcspice), Cvoid,
          (Ref{Plane}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          plane, normal, constant)
    normal, constant[]
end

"""
    pl2nvp(plane)

Return a unit normal vector and point that define a specified plane.

### Arguments ###

- `plane`: A plane

### Output ###

Returns a tuple consisting of

- `normal`: A normal vector and...
- `point`: ... point defining the geometric plane represented by `plane`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pl2nvp_c.html)
"""
function pl2nvp(plane)
    normal = Array{SpiceDouble}(undef, 3)
    point = Array{SpiceDouble}(undef, 3)
    ccall((:pl2nvp_c, libcspice), Cvoid,
          (Ref{Plane}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          plane, normal, point)
    normal, point
end

"""
    pl2psv(plane)

Return a point and two orthogonal spanning vectors that define a specified plane.

### Arguments ###

- `plane`: A plane

### Output ###

Returns a tuple consisting of a point in the `plane` and two vectors spanning
the input plane.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pl2psv_c.html)
"""
function pl2psv(plane)
    point = Array{SpiceDouble}(undef, 3)
    span1 = Array{SpiceDouble}(undef, 3)
    span2 = Array{SpiceDouble}(undef, 3)
    ccall((:pl2psv_c, libcspice), Cvoid,
          (Ref{Plane}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          plane, point, span1, span2)
    point, span1, span2
end

"""
    pltar(vrtces, plates)

Compute the total area of a collection of triangular plates.

### Arguments ###

- `vrtces`: Array of vertices
- `plates`: Array of plates

### Output ###

Returns the area.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pltar_c.html)
"""
function pltar(vrtces, plates)
    nv = length(vrtces)
    np = length(plates)
    vrtces_ = array_to_cmatrix(vrtces)
    plates_ = array_to_cmatrix(plates)
    res = ccall((:pltar_c, libcspice), SpiceDouble,
                (SpiceInt, Ref{SpiceDouble}, SpiceInt, Ref{SpiceInt}),
                nv, vrtces_, np, plates_)
    handleerror()
    res
end

"""
    pltexp(iverts, delta)

Expand a triangular plate by a specified amount. The expanded plate is co-planar with,
and has the same orientation as, the original. The centroids of the two plates coincide.

### Arguments ###

- `iverts`: Vertices of the plate to be expanded
- `delta`: Fraction by which the plate is to be expanded

### Output ###

Returns the vertices of the expanded plate.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pltexp_c.html)
"""
function pltexp(iverts, delta)
    iverts_ = array_to_cmatrix(iverts)
    overts = Array{SpiceDouble}(undef, 3, 3)
    ccall((:pltexp_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble}),
          iverts_, delta, overts)
    cmatrix_to_array(overts)
end

"""
    pltnp(point, v1, v2, v3)

Find the nearest point on a triangular plate to a given point.

### Arguments ###

- `point`: A point in 3-dimensional space.
- `v1`, `v2`, `v3`: Vertices of a triangular plate

### Output ###

Returns a tuple consisting of

- `pnear`: Nearest point on the plate to `point`
- `dist`: Distance between `pnear` and `point`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pltnp_c.html)
"""
function pltnp(point, v1, v2, v3)
    @checkdims 3 point v1 v2 v3
    pnear = Array{SpiceDouble}(undef, 3)
    dist = Ref{SpiceDouble}()
    ccall((:pltnp_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceDouble}),
          point, v1, v2, v3, pnear, dist)
    pnear, dist[]
end

"""
    pltnrm(v1, v2, v3)

Compute an outward normal vector of a triangular plate.  The vector does not
necessarily have unit length.

### Arguments ###

- `v1`, `v2`, `v3`: Vertices of a plate

### Output ###

Returns the plate's outward normal vector.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pltnrm_c.html)
"""
function pltnrm(v1, v2, v3)
    @checkdims 3 v1 v2 v3
    normal = Array{SpiceDouble}(undef, 3)
    ccall((:pltnrm_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          v1, v2, v3, normal)
    normal
end

"""
    pltvol(vrtces, plates)

Compute the volume of a three-dimensional region bounded by a collection of triangular plates.

### Arguments ###

- `vrtces`: Array of vertices
- `plates`: Array of plates

### Output ###

Returns the volume of the spatial region bounded by the plates.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pltvol_c.html)
"""
function pltvol(vrtces, plates)
    nv = length(vrtces)
    np = length(plates)
    vrtces_ = array_to_cmatrix(vrtces, n=3)
    plates_ = array_to_cmatrix(plates, n=3)
    res = ccall((:pltvol_c, libcspice), SpiceDouble,
                (SpiceInt, Ref{SpiceDouble}, SpiceInt, Ref{SpiceInt}),
                nv, vrtces_, np, plates_)
    handleerror()
    res
end

"""
    polyds(coeffs, nderiv, t)

Compute the value of a polynomial and it's first `nderiv` derivatives at the value `t`.

### Arguments ###

- `coeffs`: Coefficients of the polynomial to be evaluated
- `nderiv`: Number of derivatives to compute
- `t`: Point to evaluate the polynomial and derivatives

### Output ###

Returns the value of the polynomial and the derivatives as an array.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/polyds_c.html)
"""
function polyds(coeffs, nderiv, t)
    deg = length(coeffs) - 1
    p = Array{SpiceDouble}(undef, nderiv + 1)
    ccall((:polyds_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceDouble, Ref{SpiceDouble}),
          coeffs, deg, nderiv, t, p)
    p
end

function _pos(str, substr, start)
    res = ccall((:pos_c, libcspice), SpiceInt,
                (Cstring, Cstring, SpiceInt),
                str, substr, start - 1)
    handleerror()
    res == -1 ? res : res + 1
end

@deprecate pos(str, substr, start) first(findnext(substr, str, start))

"""
    pos(str, substr, start)

!!! warning "Deprecated"
    Use `first(findnext(substr, str, start))` instead.
"""
pos

function _posr(str, substr, start)
    res = ccall((:posr_c, libcspice), SpiceInt,
                (Cstring, Cstring, SpiceInt),
                str, substr, start - 1)
    handleerror()
    res == -1 ? res : res + 1
end

@deprecate posr(str, substr, start) first(findprev(substr, str, start))

"""
    posr(str, substr, start)

!!! warning "Deprecated"
    Use `first(findprev(substr, str, start))` instead.
"""
posr

"""
    prop2b(gm, pvinit, dt)

Given a central mass and the state of massless body at time `t_0`, this routine determines
the state as predicted by a two-body force model at time `t_0 + dt`.

### Arguments ###

- `gm`: Gravity of the central mass.
- `pvinit`: Initial state from which to propagate a state.
- `dt`: Time offset from initial state to propagate to.

### Output ###

Returns the propagated state.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/prop2b_c.html)
"""
function prop2b(gm, pvinit, dt)
    @checkdims 6 pvinit
    pvprop = Array{SpiceDouble}(undef, 6)
    ccall((:prop2b_c, libcspice), Cvoid,
          (SpiceDouble, Ref{SpiceDouble}, SpiceDouble, Ref{SpiceDouble}),
          gm, pvinit, dt, pvprop)
    handleerror()
    pvprop
end

function _prsdp(str)
    dp = Ref{SpiceDouble}()
    ccall((:prsdp_c, libcspice), Cvoid, (Cstring, Ref{SpiceDouble}), str, dp)
    dp[]
end

@deprecate prsdp(str) parse(Float64, str)

"""
    prsdp(str)

!!! warning "Deprecated"
    Use `parse(Float64, str)` instead.
"""
prsdp

function _prsint(str)
    int = Ref{SpiceInt}()
    ccall((:prsint_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}), str, int)
    int[]
end

@deprecate prsint(str) parse(Int, str)

"""
    prsint(str)

!!! warning "Deprecated"
    Use `parse(Int, str)` instead.
"""
prsint

"""
    psv2pl(point, span1, span2)

Make a plane from a point and two spanning vectors.

### Arguments ###

- `point`, `span1`, `span2`: A point and two spanning vectors defining a plane

### Output ###

Returns the plane.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/psv2pl_c.html)
"""
function psv2pl(point, span1, span2)
    @checkdims 3 point span1 span2
    plane = Ref{Plane}()
    ccall((:psv2pl_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{Plane}),
          point, span1, span2, plane)
    handleerror()
    plane[]
end

"""
    pxform(from, to, et)

Return the matrix that transforms position vectors from one specified frame to
another at a specified epoch.

### Arguments ###

- `from`: Name of the frame to transform from
- `to`: Name of the frame to transform to
- `et`: Epoch of the rotation matrix

### Output ###

Returns the rotation matrix.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pxform_c.html)
"""
function pxform(from, to, et)
    rot = Array{SpiceDouble}(undef, 3, 3)
    ccall((:pxform_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceDouble, Ref{SpiceDouble}),
          from, to, et, rot)
    handleerror()
    permutedims(rot)
end

"""
    pxfrm2(from, to, etfrom, etto)

Return the 3x3 matrix that transforms position vectors from one specified frame at a
specified epoch to another specified frame at another specified epoch.

### Arguments ###

- `from`: Name of the frame to transform from
- `to`: Name of the frame to transform to
- `etfrom`: Evaluation time of `from` frame
- `etto`: Evaluation time of `to` frame

### Output ###

Returns a position transformation matrix from frame `from` to frame `to`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/pxfrm2_c.html)
"""
function pxfrm2(from, to, etfrom, etto)
    rot = Array{SpiceDouble}(undef, 3, 3)
    ccall((:pxfrm2_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          from, to, etfrom, etto, rot)
    handleerror()
    permutedims(rot)
end


text:
##########################################################
# Rafał Nowak
#
# Instytut Informatyki
# Uniwersytet Wrocławski
#
#
# Pomocniczne materiały do wykładu z analizy numerycznej
#
##########################################################

# Arytmetyka pojedynczej precyzji
# Float32  (float w języku C)
x = 1.0f0
y = 0.5f-8
typeof(x)
typeof(y)

# Domyślnie używana jest arytmetyka podwójnej precyzji,
# Float64  (double w języku C)
x = 1.0;
y = 2.0e-7;
typeof(x)
typeof(y)

# Konwersja to wybranego typu
x = Float32( 1.0 )
y = Float64( 1.0 )
typeof(x)
typeof(y)

function test(x, n = 10^9)
  y = zero(typeof(x)); # funkcja zero( .. ) zwraca zero w odpowiedniej arytmetyce
  for i=1:n
    y = y+x;
  end
  y
end

tic()
x=test( Int32(1) , 10^8); # ok. 6 sekund
toc()
tic()
x=test( Int64(1) , 10^8); # ok. 6 microsekund
toc()
tic()
x=test( 1.0f-5 , 10^9 );  # ok. 1 sek.
toc()
tic()
x=test( 1.0e-5 , 10^9 );  # ok. 1 sek., ale zupełnie inny wynik
toc()

# Reprezentacja zmiennopozycyjna
(s,c,m) = ( 1.0,
           -5,
           1.0 + 2.0^(-3) + 2.0^(-5) ) # s-znak, c-cecha, m-mantysa

x = s*2.0^c*m;
typeof(x)

repr = bits(x); # bits zwraca reprezentację binarną
rs = repr[1];
rc = repr[2:2+11-1];
rm = repr[2+11:end]
@printf("Reprezentacja liczby x = %f  to\n%c|%s|%s\n", x, rs,rc,rm);

tmp = bits( c+1023 );
@printf("  %s  - reprezentacja binarna liczby całkowitej c+1023  (c=%d)\n", tmp[end-10:end], c);

#####
# (s,c,m) = (1.0, 1024, 1.5)  =>  NaN, np. x = 1.0 / 0.0
# (s,c,m) = (1.0, 1024, 1.0)  =>  Inf
# println(bits(NaN))
# println(bits(Inf))
# println(bits(e64))
##########################################################
x = Float16(1.3); # Half precision (16 bits)
@printf("Wartość 1.3 w arytmetyce 16 bitowej = %c.8f  = %.8f\n",  '%', x);  # 1.29980469
@printf("Wartość 1.3 w arytmetyce 16 bitowej = %c.16f = %.16f\n", '%', x);  # 1.2998046875000000
##########################################################
# Obliczanie *precyzji arytmetyki*
# one - wartość 1.0 w danej arytmetyce
# Przykład:  myeps( float16 )
function myeps(ftype = float64)
  one = ftype(1.0);
  two = ftype(2.0);
  x   = one;
  while x+one>one
    x=x/two;
  end;
  return x;
end

e16 = myeps( Float16 );
e32 = myeps( Float32 );
e64 = myeps( Float64 );

@printf("Precyzja arytmetyki float16 = %.4e\n", e16);
@printf("Precyzja arytmetyki float32 = %.4e\n", e32);
@printf("Precyzja arytmetyki float64 = %.4e\n", e64);

# Julia ma wbudowaną funkcję eps(...)
# help("Base.eps");
# Wniosek:  2*myeps = eps
2.0*e32 == eps(Float32(1.0))

t = 10; # +1 oznacza jeden dodatkowy bit mantysy w standardzie IEEE 754
e16 == 0.5 * 2.0^(-t)

for (ftype,t) in [ (Float16,10), (Float32,23), (Float64,52) ]
  assert( myeps(ftype) == 0.5 * 2.0 ^ (-t) );
end;

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
# HalidePerovskites
# Uses the newly forked off FeynmanKadanoffOsakaHellwarth.jl package
# Codes by Jarvist Moore Frost, 2017
# These codes were developed with Julia 0.5.0, and requires the Optim and Plots packages.

# This file, when run under Julia, should regenerate all data (and plots) associated with Arxiv paper:
# https://arxiv.org/abs/1704.05404
# Polaron mobility in halide perovskites
# Jarvist Moore Frost
# (Submitted on 18 Apr 2017 [v1])

push!(LOAD_PATH,"../src/") # load module from local directory

using PolaronMobility

##### load in library routines... #####
# Plot figures with Plots, which defaults to Pyplot backend
using Plots
#pyplot()
gr()
#default(grid=false) # No silly dotted grid lines
#default(size=(400,300)) # A good small size for two-column EPS output
#gr()
#default(size=(800,600)) # Nice size for small-ish PNGs for slides

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space

#####
# Call simulation

# CsSnX3 X={Cl,Br,I}
# L. Huang, W.Lambrecht - PRB 88, 165203 (2013)
# Dielectric consts, from TABLE VII
# Effective masses from TABLE VI, mh*
const cm1=2.997e10 # cm-1 to Herz

#Ts,Kμs, Hμs, FHIPμs, ks, Ms, As, Bs, Cs, Fs, Taus
#effectivemass=0.12 # the bare-electron band effective-mass.
# --> 0.12 for electrons and 0.15 for holes, in MAPI. See 2014 PRB.
# MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
println("OK, solving Polaron problem...")
MAPIe = polaronmobility("MAPI-electron", [150,300], 4.5, 24.1, 2.25E12, 0.12; verbose=true, figures=false)
MAPIh = polaronmobility("MAPI-hole",     [150,300], 4.5, 24.1, 2.25E12, 0.15; verbose=true, figures=false)

MAPIe = polaronmobility("MAPI-electron", [150,300], 4.5, 24.1, 1.25E12, 0.12; verbose=true, figures=false)

s=ImX(0.1:0.1:6,MAPIe.v[2],MAPIe.w[2],MAPIe.βred[2], MAPIe.α[1],MAPIe.ω[1],MAPIe.mb[1])
plot( s.nu,s.ImX,label="ImX",
    markersize=3,marker=:downtriangle, xlab="nu (units Omega)",ylab="ImX")
savefig("MAPIe-ImX.png")
plot( s.nu,s.μ,label="mu",
    markersize=3,marker=:uptriangle, xlab="nu (units Omega)",ylab="Mob")
savefig("MAPIe-mu.png")

println("That's me!")


text:
export PIBasis

# ---------------------- Implementation of the PIBasisSpec

"""
`struct PIBasisSpec`
"""
struct PIBasisSpec
   orders::Vector{Int}     # order (length) of ith basis function
   iAA2iA::Matrix{Int}     # where in A can we find the ith basis function
end

==(B1::PIBasisSpec, B2::PIBasisSpec) = (
         (B1.b2iA == B2.b2iA) &&
         (B1.iAA2iA == B2.iAA2iA) )

Base.length(spec::PIBasisSpec) = length(spec.orders)

function _get_pibfcn(spec0, Aspec, vv)
   vv1 = vv[2:end]
   vvnz = vv1[findall(vv1 .!= 0)]
   return (spec0[vv[1]], Aspec[vvnz])
end

function _get_pibfcn(Aspec, vv)
   vvnz = vv[findall(vv .!= 0)]
   return Aspec[vvnz]
end

function PIBasisSpec( basis1p::OneParticleBasis,
                      maxν::Integer, maxdeg::Real;
                      Deg = NaiveTotalDegree(),
                      property = nothing,
                      filterfun = _->true,
                      constant = false )
   # would make sense to construct the basis1p spec here?

   # get the basis spec of the one-particle basis
   #  Aspec[i] described the basis function that will get written into A[i]
   Aspec = get_spec(basis1p)

   # we assume that `Aspec` is sorted by degree, but best to double-check this
   # since the notion of degree used to construct `Aspec` might be different
   # from the one used to construct AAspec.
   if !issorted(Aspec; by = b -> degree(b, Deg, basis1p))
      error("""PIBasisSpec : AAspec construction failed because Aspec is not
               sorted by degree. This could e.g. happen if an incompatible
               notion of degree was used to construct the 1-p basis spec.""")
   end
   # An AA basis function is given by a tuple 𝒗 = vv. Each index 𝒗ᵢ = vv[i]
   # corresponds to the basis function Aspec[𝒗ᵢ] and the tuple
   # 𝒗 = (𝒗₁, ...) to a product basis function
   #   ∏ A_{vₐ}
   tup2b = vv -> _get_pibfcn(Aspec, vv)

   #  degree of a basis function ↦ is it admissible?
   admissible = b -> (degree(b, Deg, basis1p) <= maxdeg)

   if property != nothing
      filter1 = b -> filterfun(b) && filter(property, b)
   else
      filter1 = filterfun
   end

   # we can now construct the basis specification; the `ordered = true`
   # keyword signifies that this is a permutation-invariant basis
   AAspec = gensparse(; NU = maxν,
                        tup2b = tup2b,
                        admissible = admissible,
                        ordered = true,
                        maxvv = [length(Aspec) for _=1:maxν],
                        filter = filter1,
                        constant = constant )

   return PIBasisSpec(AAspec)
end

function PIBasisSpec(AAspec)
   orders = zeros(Int, length(AAspec))
   iAA2iA = zeros(Int, (length(AAspec), length(AAspec[1])))
   for (iAA, vv) in enumerate(AAspec)
      # we use reverse because gensparse constructs the indices in
      # ascending order, but we want descending here.
      # (I don't remember why though)
      iAA2iA[iAA, :] .= reverse(vv)
      orders[iAA] = length( findall( vv .!= 0 ) )
   end
   return PIBasisSpec(orders, iAA2iA)
end

get_spec(AAspec::PIBasisSpec, i::Integer) = AAspec.iAA2iA[i, 1:AAspec.orders[i]]

# --------------------------------- PIBasis implementation

"""
`mutable struct PIBasis:` implementation of a permutation-invariant
basis based on the density projection trick.

The standard constructor is
```
PIBasis(basis1p, N, D, maxdeg)
```
* `basis1p` : a one-particle basis
* `N` : maximum interaction order
* `D` : an abstract degee specification, e.g., SparsePSHDegree
* `maxdeg` : the maximum polynomial degree as measured by `D`
"""
mutable struct PIBasis{BOP, REAL} <: ACEBasis
   basis1p::BOP             # a one-particle basis
   spec::PIBasisSpec
   real::REAL     # could be `real` or `identity` to keep AA complex
   # evaluator    # classic vs graph
end

cutoff(basis::PIBasis) = cutoff(basis.basis1p)

==(B1::PIBasis, B2::PIBasis) = ACE._allfieldsequal(B1, B2)

# TODO: allow the option of converting to real part?
fltype(basis::PIBasis) = basis.real( fltype(basis.basis1p) )
rfltype(basis::PIBasis) = real( fltype(basis) )

Base.length(basis::PIBasis) = length(basis.spec)

PIBasis(basis1p, args...; isreal = true, kwargs...) =
   PIBasis(basis1p, PIBasisSpec(basis1p, args...; kwargs...),
           isreal ? Base.real : Base.identity )

get_spec(pibasis::PIBasis) =
   [ get_spec(pibasis, i) for i = 1:length(pibasis) ]

get_spec(pibasis::PIBasis, i::Integer) =
      get_spec.( Ref(pibasis.basis1p), get_spec(pibasis.spec, i) )

setreal(basis::PIBasis, isreal::Bool) =
   PIBasis(basis.basis1p, basis.spec, isreal)

# function scaling(pibasis::PIBasis, p)
#    ww = zeros(Float64, length(pibasis))
#    for iz0 = 1:numz(pibasis)
#       wwin = @view ww[pibasis.inner[iz0].AAindices]
#       for i = 1:length(pibasis.inner[iz0])
#          bspec = get_basis_spec(pibasis, iz0, i)
#          wwin[i] = scaling(bspec, p)
#       end
#    end
#    return ww
# end

# graphevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, DAGEvaluator())
#
# standardevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, StandardEvaluator())
#
#
#
# -------------------------------------------------
# FIO codes

write_dict(basis::PIBasis) =
   Dict(  "__id__" => "ACE_PIBasis",
         "basis1p" => write_dict(basis.basis1p),
            "spec" => write_dict(basis.spec),
            "real" => basis.real )

read_dict(::Val{:ACE_PIBasis}, D::Dict) =
   PIBasis( read_dict(D["basis1p"]),
            read_dict(D["spec"]),
            D["real"] )

write_dict(spec::PIBasisSpec) =
   Dict( "__id__" => "ACE_PIBasisSpec",
         "orders" => spec.orders,
         "iAA2iA" => write_dict(spec.iAA2iA) )

read_dict(::Val{:ACE_PIBasisSpec}, D::Dict) =
   PIBasisSpec( D["orders"], read_dict(D["iAA2iA"]) )

# -------------------------------------------------
# Evaluation codes

alloc_B(basis::PIBasis, args...) = zeros( fltype(basis), length(basis) )

alloc_temp(basis::PIBasis, args...) =
      ( A = alloc_B(basis.basis1p, args...),
        tmp1p = alloc_temp(basis.basis1p, args...),
      )

function evaluate!(AA, tmp, basis::PIBasis, config::AbstractConfiguration)
   A = evaluate!(tmp.A, tmp.tmp1p, basis.basis1p, config)
   fill!(AA, 1)
   for iAA = 1:length(basis)
      aa = one(eltype(A))
      for t = 1:basis.spec.orders[iAA]
         aa *= A[ basis.spec.iAA2iA[ iAA, t ] ]
      end
      AA[iAA] = basis.real(aa)
   end
   return AA
end

# -------------------------------------------------
# gradients

gradtype(basis::PIBasis) = gradtype(basis.basis1p)

alloc_dB(basis::PIBasis, nmax::Integer) =
      zeros(gradtype(basis), (length(basis), nmax))

alloc_temp_d(basis::PIBasis, nmax::Integer) =
      (
        A = alloc_B(basis.basis1p),
        dA = alloc_dB(basis.basis1p, nmax),
        tmp_basis1p = alloc_temp(basis.basis1p),
        tmpd_basis1p = alloc_temp_d(basis.basis1p),
        # ---- adjoint stuff
        dAAdA = zeros(fltype(basis.basis1p),
                      maximum(basis.spec.orders))
      )

# TODO: This is a naive forwardmode implementation; we need
#       to switch this to reverse mode differentiation

function evaluate_ed!(AA, dAA, tmpd, basis::PIBasis,
                      cfg::AbstractConfiguration)
   evaluate_ed!(tmpd.A, tmpd.dA, tmpd.tmpd_basis1p, basis.basis1p, cfg)
   evaluate_ed!(AA, dAA, tmpd, basis, tmpd.A, tmpd.dA)
end

function evaluate_ed!(AA, dAA, tmpd, basis::PIBasis,
                      A::AbstractVector, dA::AbstractMatrix)
   dAAdA = tmpd.dAAdA
   orders = basis.spec.orders
   iAA2iA = basis.spec.iAA2iA

   for iAA = 1:length(basis)
      ord = orders[iAA]
      # ----- compute the local adjoints dAA / dA
      # dAAdA[a] ← ∏_{t ≂̸ a} A_{v_t}
      # TODO - optimize?
      # TODO: maybe this code could be shared with the evaluator
      # Forward pass:
      dAAdA[1] = 1
      AAfwd = A[iAA2iA[iAA, 1]]
      for a = 2:ord
         dAAdA[a] = AAfwd
         AAfwd *= A[iAA2iA[iAA, a]]
      end
      AA[iAA] = basis.real(AAfwd)
      # backward pass
      AAbwd = A[iAA2iA[iAA, ord]]
      for a = ord-1:-1:1
         dAAdA[a] *= AAbwd
         AAbwd *= A[iAA2iA[iAA, a]]
      end

      # ----- now convert them into dAA / dX
      for j = 1:size(dA, 2)
         dAA[iAA, j] = sum(dAAdA[a] * dA[iAA2iA[iAA, a], j]
                           for a = 1:ord) |> basis.real
      end
   end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis: subs
using Test

expression = :(1 + u + v + u*v^2)
data = (:u => 1.0, :v => 2.0)
result = subs(expression, data)
@debug "subs(expression, data) where" expression data
@test isapprox(result, 8.0)


text:
# Kernels for elementwise Array,Array,Array->Array ops with equal sized
# arrays used for activation backward functions (with input x,y,dy output dx)

fp = open("cuda111.cu","w")
#using Knet: actback_ops

function cuda111src(f, j=f, ex="$f(xi,yi)"; BLK=256, THR=256)
  sprint() do s
    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,
"""
__global__ void _$(F)_111(int n, $T *x_, $T *y_, $T *dy_, $T *dx_) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T x=x_[i];
    $T y=y_[i];
    $T dy=dy_[i];
    dx_[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void $(F)_111(int n, $T *x, $T *y, $T *dy, $T *dx) {
    _$(F)_111<<<$BLK,$THR>>>(n,x,y,dy,dx);
  }
  $DLLEXPORT void $(F)_111_stream(int n, $T *x, $T *y, $T *dy, $T *dx, cudaStream_t STR) {
    _$(F)_111<<<$BLK,$THR,0,STR>>>(n,x,y,dy,dx);
  }
}
""")
    end
  end
end

for a in actback_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda111src(a...))
end

close(fp)


text:
@testset "sNES" begin

function assign_weights_wrapper(candi_ixs::Vector{Int})
    candidates = BlackBoxOptim.Candidate{Float64}[BlackBoxOptim.Candidate{Float64}([0.0], i, NaN) for i in candi_ixs]
    u = BlackBoxOptim.fitness_shaping_utilities_linear(length(candi_ixs))
    BlackBoxOptim.assign_weights!(similar(u), candidates, u)
end

@testset "assign_weights!()" begin
    @testset "when indices are already ordered" begin
        u = assign_weights_wrapper([1, 2])
        @test length(u) == 2
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test u[1] == 0.5
        @test u[2] == -0.5

        u = assign_weights_wrapper([1, 2, 3])
        @test length(u) == 3
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test u[1] == 1.0 - 1/3
        @test u[2] == -1/3
        @test u[3] == -1/3

        u = assign_weights_wrapper([1, 2, 3, 4])
        @test length(u) == 4
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test isapprox(u[1], 2/3-1/4)
        @test isapprox(u[2], 1/3-1/4)
        @test u[3] == -1/4
        @test u[4] == -1/4

        u = assign_weights_wrapper([1, 2, 3, 4, 5])
        @test length(u) == 5
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test isapprox(u[1], 2/3-1/5)
        @test isapprox(u[2], 1/3-1/5)
        @test u[3] == -1/5
        @test u[4] == -1/5
        @test u[5] == -1/5

        u = assign_weights_wrapper([1, 2, 3, 4, 5, 6])
        @test length(u) == 6
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test isapprox(u[1], (3/3)/(3/3+2/3+1/3)-1/6)
        @test isapprox(u[2], (2/3)/(3/3+2/3+1/3)-1/6)
        @test isapprox(u[3], (1/3)/(3/3+2/3+1/3)-1/6)
        @test u[4] == -1/6
        @test u[5] == -1/6
        @test u[6] == -1/6
    end

    @testset "when indices are not ordered" begin
        u = assign_weights_wrapper([4, 1, 2, 3])
        @test isapprox(sum(u), 0.0; atol=1E-10)
        @test isapprox(u[4], 2/3-1/4)
        @test isapprox(u[1], 1/3-1/4)
        @test u[2] == -1/4
        @test u[3] == -1/4
    end
end

end