text:
using ESDL
using Test
using Dates
using IterTools
using WeightedOnlineStats

c=Cube()
d = subsetcube(c,variable=["air_temperature_2m", "gross_primary_productivity", "land_surface_temperature"],lon=(30,31),lat=(50,51),
              time=(Date("2002-01-01"),Date("2008-12-31")))
mytable = CubeTable(value=d,include_axes=("lon", "lat", "time"), expandaxes=("var",))

mytable2 = CubeTable(data=d,include_axes=("lon", "lat", "time", "variable"))

@testset "cubefittable and WeightedCovMatrix fittable" begin
    covmCube = cubefittable(mytable, WeightedCovMatrix, :value, weight=(x->cosd(x.lat)),showprog=false)
    covmVal = WeightedCovMatrix()
    for row in mytable
        obs = row.value
        obslat = row.lat
        if !any(ismissing.(obs))
            fit!(covmVal, obs, cosd(obslat))
        end
    end
    @test all(isapprox.(covmCube.data, value(covmVal)))
end

@testset "remaining functions in tablestats.jl" begin
    meanCube = cubefittable(mytable2, WeightedMean, :data, weight=(x->cosd(x.lat)),
                            by=(:variable,))
    weightedTableMeans = fittable(mytable2, WeightedMean, :data, by=(:variable,),
                                    weight=(x -> cosd(x.lat)))
    foreach(enumerate(["air_temperature_2m", "gross_primary_productivity", "land_surface_temperature"])) do (i,v)
      @test isapprox(value(weightedTableMeans)[(v,)], meanCube.data[i])
    end
end


text:
function plotFunc_timecourse(Sim::Module)
    rc("figure",figsize = (10.5,4));
    rc("font",family = "Arial");
    rc("font",size = 12);
    rc("axes",linewidth = 1);
    rc("lines",linewidth = 2);

    subplot(1,2,1);
    plot(Sim.t,Sim.p53_conc[:,1],"k",label="1μM etoposide");
    plot(Sim.t,Sim.p53_conc[:,2],"r",label="100μM etoposide");
    xlim(0,48);
    xticks([0,12,24,36,48]);
    ylim(0,8);
    yticks([1,3,5,7]);
    xlabel("Time (hours)");
    ylabel("[p53]");
    legend(loc="upper right",frameon=false);
    text(1,7.5,"U-2 OS");

    subplot(1,2,2);
    plot(Sim.t,Sim.p53_conc[:,3],"k",label="5μM etoposide");
    plot(Sim.t,Sim.p53_conc[:,4],"r",label="200μM etoposide");
    xlim(0,48);
    xticks([0,12,24,36,48]);
    ylim(0,5.5);
    yticks([0,1,2,3,4,5]);
    xlabel("Time (hours)");
    ylabel("[p53]");
    legend(loc="upper right",frameon=false);
    text(1,5.15,"MCF7")

    savefig("p53_dynamics.png",bbox_inches="tight");

end

ext: jl


text:
# indicator of the L2 norm ball with given radius

"""
  IndBallL2(r::Real=1.0)

Returns the function `g = ind{x : ||x|| ⩽ r}`, for a real parameter `r > 0`.
"""

immutable IndBallL2{R <: Real} <: IndicatorConvex
  r::R
  function IndBallL2(r::R)
    if r <= 0
      error("parameter r must be positive")
    else
      new(r)
    end
  end
end

IndBallL2{R <: Real}(r::R=1.0) = IndBallL2{R}(r)

function (f::IndBallL2){T <: RealOrComplex}(x::AbstractArray{T})
  if vecnorm(x) - f.r > 1e-14
    return +Inf
  end
  return 0.0
end

function prox!{T <: RealOrComplex}(f::IndBallL2, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  scal = f.r/vecnorm(x)
  if scal > 1
    y[:] = x
    return 0.0
  end
  for k in eachindex(x)
    y[k] = scal*x[k]
  end
  return 0.0
end

fun_name(f::IndBallL2) = "indicator of an L2 norm ball"
fun_dom(f::IndBallL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndBallL2) = "x ↦ 0 if ||x|| ⩽ r, +∞ otherwise"
fun_params(f::IndBallL2) = "r = $(f.r)"

function prox_naive{T <: RealOrComplex}(f::IndBallL2, x::AbstractArray{T}, gamma::Real=1.0)
  normx = vecnorm(x)
  if normx > f.r
    y = (f.r/normx)*x
  else
    y = x
  end
  return y, 0.0
end


text:
# ---
# title: 86. Partition List
# id: problem86
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Two Pointers
# link: <https://leetcode.com/problems/partition-list/description/>
# hidden: true
# ---
#
# Given a linked list and a value _x_ , partition it such that all nodes less
# than _x_ come before nodes greater than or equal to _x_.
#
# You should preserve the original relative order of the nodes in each of the
# two partitions.
#
# **Example:**
#
#
#
#     Input: head = 1->4->3->2->5->2, _x_ = 3
#     Output: 1->2->2->4->3->5
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""An `IntervalBox` is an `N`-dimensional rectangular box, given
by a Cartesian product of a vector of `N` `Interval`s.
"""
struct IntervalBox{N,T}
    v::SVector{N, Interval{T}}
end

# IntervalBox(x::Interval) = IntervalBox( SVector(x) )  # single interval treated as tuple with one element

IntervalBox(x::Interval...) = IntervalBox(SVector(x))
IntervalBox(x::SVector) = IntervalBox(interval.(x))
IntervalBox(x::Tuple) = IntervalBox(SVector(x))
IntervalBox(x::Real) = IntervalBox(interval.(x))
IntervalBox(x...) = IntervalBox(x)
IntervalBox(x) = IntervalBox(x...)
IntervalBox(X::IntervalBox, n) = foldl(×, Iterators.repeated(X, n))

Base.@propagate_inbounds Base.getindex(X::IntervalBox, i) = X.v[i]

setindex(X::IntervalBox, y, i) = IntervalBox( setindex(X.v, y, i) )

# iteration:

start(X::IntervalBox{N,T}) where {N,T} = 1

next(X::IntervalBox{N,T}, state) where {N,T} = (X[state], state+1)

done(X::IntervalBox{N,T}, state) where {N,T} = state > N

eltype(::Type{IntervalBox{N,T}}) where {N,T} = Interval{T} # Note that this is defined for the type

length(X::IntervalBox{N,T}) where {N,T} = N

# IntervalBox(xx) = IntervalBox(Interval.(xx))
# IntervalBox(xx::SVector) where {N,T} = IntervalBox(Interval.(xx))

## arithmetic operations
# Note that standard arithmetic operations are implemented automatically by FixedSizeArrays.jl
"""
    mid(X::IntervalBox, α=0.5)

Return a vector of the `mid` of each interval composing the `IntervalBox`.

See `mid(X::Interval, α=0.5)` for more informations.
"""
mid(X::IntervalBox) = mid.(X.v)
mid(X::IntervalBox, α) = mid.(X.v, α)

big(X::IntervalBox) = big.(X)

## set operations

# TODO: Update to use generator
⊆(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    all(X.v .⊆ Y.v)

∩(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    IntervalBox(X.v .∩ Y.v)
∪(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} =
    IntervalBox(X.v .∪ Y.v)

#=
On Julia 0.6 can now write
∩{N,T}(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) = IntervalBox(NTuple{N, Interval{Float64}}( (X[i] ∩ Y[i]) for i in 1:N))
=#

isempty(X::IntervalBox) = any(isempty, X.v)

diam(X::IntervalBox) = maximum(diam.(X.v))

emptyinterval(X::IntervalBox{N,T}) where {N,T} = IntervalBox(emptyinterval.(X.v))

isinf(X::IntervalBox) = any(isinf.(X))

isinterior(X::IntervalBox{N,T}, Y::IntervalBox{N,T}) where {N,T} = all(isinterior.(X, Y))

contains_zero(X::SVector) = all(contains_zero.(X))
contains_zero(X::IntervalBox) = all(contains_zero.(X))

×(a::Interval...) = IntervalBox(a...)
×(a::Interval, b::IntervalBox) = IntervalBox(a, b.v...)
×(a::IntervalBox, b::Interval) = IntervalBox(a.v..., b)
×(a::IntervalBox, b::IntervalBox) = IntervalBox(a.v..., b.v...)

IntervalBox(x::Interval, ::Type{Val{n}}) where {n} = IntervalBox(SVector(ntuple(i->x, Val{n})))

IntervalBox(x::Interval, n::Int) = IntervalBox(x, Val{n})

dot(x::IntervalBox, y::IntervalBox) = dot(x.v, y.v)


text:
export DouglasRachford
@doc raw"""
     DouglasRachford(M, F, proxMaps, x)
Computes the Douglas-Rachford algorithm on the manifold $\mathcal M$, initial
data $x_0$ and the (two) proximal maps `proxMaps`.

For $k>2$ proximal
maps the problem is reformulated using the parallelDouglasRachford: a vectorial
proximal map on the power manifold $\mathcal M^k$ and the proximal map of the
set that identifies all entries again, i.e. the Karcher mean. This henve also
boild down to two proximal maps, though each evauates proximal maps in parallel,
i.e. component wise in a vector.

# Input
* `M` – a Riemannian Manifold $\mathcal M$
* `F` – a cost function consisting of a sum of cost functions
* `proxes` – functions of the form `(λ,x)->...` performing a proximal map,
  where `λ` denotes the proximal parameter, for each of the summands of `F`.
* `x0` – initial data $x_0 ∈ \mathcal M$

# Optional values
the default parameter is given in brackets
* `λ` – (`(iter) -> 1.0`) function to provide the value for the proximal parameter
  during the calls
* `α` – (`(iter) -> 0.9`) relaxation of the step from old to new iterate, i.e.
  $t_{k+1} = g(α_k; t_k, s_k)$, where $s_k$ is the result
  of the double reflection involved in the DR algorithm
* `R` – ([`reflect`](@ref)) method employed in the iteration
  to perform the reflection of `x` at the prox `p`.
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(200),`[`StopWhenChangeLess`](@ref)`(10.0^-5))`) a [`StoppingCriterion`](@ref).
* `parallel` – (`false`) clarify that we are doing a parallel DR, i.e. on a
  `PowerManifold` manifold with two proxes. This can be used to trigger
  parallel Douglas–Rachford if you enter with two proxes. Keep in mind, that a
  parallel Douglas–Rachford implicitly works on a `PowerManifold` manifold and
  its first argument is the result then (assuming all are equal after the second
  prox.
* `return_options` – (`false`) – if actiavated, the extended result, i.e. the
    complete [`Options`](@ref) re returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `xOpt` if returned
...
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function DouglasRachford(
    M::MT,
    F::TF,
    proxes::Vector{<:Any},
    x;
    λ::Tλ=(iter) -> 1.0,
    α::Tα=(iter) -> 0.9,
    R::TR=reflect,
    parallel::Int=0,
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(200), StopWhenChangeLess(10.0^-5)
    ),
    return_options=false,
    kwargs..., #especially may contain decorator options
) where {MT<:Manifold,TF,Tλ,Tα,TR}
    if length(proxes) < 2
        throw(ErrorException("Less than two proximal maps provided, the (parallel) Douglas Rachford requires (at least) two proximal maps."))
    elseif length(proxes) == 2
        prox1 = proxes[1]
        prox2 = proxes[2]
        parallel = 0
    else # more than 2 -> parallelDouglasRachford
        parallel = length(proxes)
        prox1 = (λ, x) -> [proxes[i](λ, x[i]) for i in 1:parallel]
        prox2 = (λ, x) -> fill(mean(M.manifold, x), parallel)
    end
    if parallel > 0
        M = PowerManifold(M, NestedPowerRepresentation(), parallel)
        x = [copy(x) for i in 1:parallel]
        nF = x -> F(x[1])
    else
        nF = F
    end
    p = ProximalProblem(M, nF, (prox1, prox2))
    o = DouglasRachfordOptions(x, λ, α, reflect, stopping_criterion, parallel > 0)

    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(p::ProximalProblem, o::DouglasRachfordOptions) end
function step_solver!(p::ProximalProblem, o::DouglasRachfordOptions, iter)
    pP = get_proximal_map(p, o.λ(iter), o.s, 1)
    snew = o.R(p.M, pP, o.s)
    o.x = get_proximal_map(p, o.λ(iter), snew, 2)
    snew = o.R(p.M, o.x, snew)
    # relaxation
    return o.s = shortest_geodesic(p.M, o.s, snew, o.α(iter))
end
get_solver_result(o::DouglasRachfordOptions) = o.parallel ? o.x[1] : o.x


text:
Pkg.update()
Pkg.add("ScikitLearn")
Pkg.add("PyPlot")

using ScikitLearn
using PyPlot

@sk_import datasets: (make_circles, make_moons, make_blobs)
@sk_import cluster: (estimate_bandwidth, MeanShift, MiniBatchKMeans, AgglomerativeClustering, SpectralClustering)
@sk_import cluster: (DBSCAN, AffinityPropagation, Birch)
@sk_import preprocessing: StandardScaler
@sk_import neighbors: kneighbors_graph

srand(33)

# Generate datasets. We choose the size big enough to see the scalability
# of the algorithms, but not too big to avoid too long running times
n_samples = 1500
noisy_circles = make_circles(n_samples=n_samples, factor=.5, noise=.05)
noisy_moons = make_moons(n_samples=n_samples, noise=.05)
blobs = make_blobs(n_samples=n_samples, random_state=8)
no_structure = rand(n_samples, 2), nothing

colors0 = collect("bgrcmykbgrcmykbgrcmykbgrcmyk")
colors = vcat(fill(colors0, 20)...)

clustering_names = [
    "MiniBatchKMeans", "AffinityPropagation", "MeanShift",
    "SpectralClustering", "Ward", "AgglomerativeClustering",
    "DBSCAN", "Birch"];

figure(figsize=(length(clustering_names) * 2 + 3, 9.5))
subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,
                    hspace=.01)

plot_num = 1

datasets = [noisy_circles, noisy_moons, blobs, no_structure]

for (i_dataset, dataset) in enumerate(datasets)

    X, y = dataset
    # normalize dataset for easier parameter selection
    X = fit_transform!(StandardScaler(), X)

    # estimate bandwidth for mean shift
    bandwidth = estimate_bandwidth(X, quantile=0.3)

    # connectivity matrix for structured Ward
    connectivity = kneighbors_graph(X, n_neighbors=10,
                include_self=false)[:todense]()

    # PyCall does not support numpy sparse matrices
    # make connectivity symmetric
    connectivity = 0.5 * (connectivity + connectivity')

    # create clustering estimators
    ms = MeanShift(bandwidth=bandwidth, bin_seeding=true)
    two_means = MiniBatchKMeans(n_clusters=2)
    ward = AgglomerativeClustering(n_clusters=2, linkage="ward",
                                   connectivity=connectivity)
    spectral = SpectralClustering(n_clusters=2,
                                  eigen_solver="arpack",
                                  affinity="nearest_neighbors")
    dbscan = DBSCAN(eps=.2)
    affinity_propagation = AffinityPropagation(damping=.9, preference=-200)

    average_linkage = AgglomerativeClustering(
        linkage="average", affinity="cityblock", n_clusters=2,
        connectivity=connectivity)

    birch = Birch(n_clusters=2)
    clustering_algorithms = [
        two_means, affinity_propagation, ms, spectral, ward, average_linkage,
        dbscan, birch]

    for (name, algorithm) in zip(clustering_names, clustering_algorithms)
        fit!(algorithm, X)
        y_pred = nothing
        try
            y_pred = predict(algorithm, X)
        catch e
            if isa(e, KeyError)
                y_pred = map(Int, algorithm[:labels_])
                clamp!(y_pred, 0, 27) # not sure why some algorithms return -1
            else rethrow() end
        end
        subplot(4, length(clustering_algorithms), plot_num)
        if i_dataset == 1
            title(name, size=18)
        end

        for y_val in unique(y_pred)
            selected = y_pred.==y_val
            scatter(X[selected, 1], X[selected, 2], color=string(colors0[y_val+1]), s=10)
        end

        xlim(-2, 2)
        ylim(-2, 2)
        xticks(())
        yticks(())
        plot_num += 1
    end
end

max_stars_repo_path: Chapter 7/B05321_CH07_Code02.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 7/B05321_CH07_Code02.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 7/B05321_CH07_Code02.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
using Dolfin

import Dolfin:dot, Function, det
import PyCall: PyDict

parameters["form_compiler"]["cpp_optimize"] = true

ffc_options = PyDict()
ffc_options["optimize"]= true
ffc_options["eliminate_zeros"]= true
ffc_options["precompute_basis_const"]= true
ffc_options["precompute_ip_const"]=true

# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"))
r = Expression(("scale*0.0",
                "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
                "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
                scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3)

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, -0.5, 0.0))  # Body force per unit volume
T  = Constant((0.1,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = u["geometric_dimension"]()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F["T"]*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
# change ** to ^
psi = (mu/Constant(2))*(Ic - Constant(3)) - mu*ln(J) + (lmbda/Constant(2))*(ln(J)) * (ln(J))

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == Constant(0), u, bcs, J=J, form_compiler_parameters=ffc_options)

# Save solution in VTK format
# file = File("displacement.pvd");
# file << u;

# Plot and hold solution
plot(u, mode = "displacement", interactive = true)


text:
"""
    ONet
Application to upload the ONet data to the database.
"""
module ONet
    using HTTP: HTTP
    using JSON3: JSON3
    using CSV: CSV
    using Tables: columntable, rowtable, schema, transform, rows
    using LibPQ: Connection, execute, load!, CopyIn
    using ConfParser
    # Compat
    if !@isdefined(isnothing)
        isnothing(::Any) = false
        isnothing(::Nothing) = true
    end
    if !@isdefined(ismissing)
        isnothing(::Any) = false
        isnothing(::Missing) = true
    end
    foreach(include,
            file for file ∈ readdir(@__DIR__) if !isequal("ONet.jl", file))
    export sdad_setup!, onet!
    include_dependency(joinpath(dirname(@__FILE__), "..", "confs", "config.simple"))
end


text:
import Base: position

function position! end
function velocity! end
function state! end

function velocity end
function state end


text:
using BenchmarkTools, Symbolics

using LinearAlgebra

SUITE = BenchmarkGroup()

vars = @variables a,b,c,d,e,f,g,h,i

X = [0 b c;
     d e f;
     g h i]

z = det(X) - (-b * (d*i-f*g) + c * (d*h - e*g))

F = lu(X)
FbyX = F \ X

x = (f + ((((g*(c^2)*(e^2)) / d - e*h*(c^2)) / b + (-c*e*f*g) / d + c*e*i) /
          (i + ((c*e*g) / d - c*h) / b + (-f*g) / d) - c*e) / b +
     ((g*(f^2)) / d + ((-c*e*f*g) / d + c*f*h) / b - f*i) /
     (i + ((c*e*g) / d - c*h) / b + (-f*g) / d)) / d

o = (d + (e*((c*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d))) / b + (-f*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d)) / d

@show iszero(x)
@show isone(o)
SUITE["iszero/1"] = @benchmarkable iszero($((b*(h + (-e*g) / d)) / b + (e*g) / d - h))
SUITE["isone/1"] = @benchmarkable FbyX == I

SUITE["iszero/2"] = @benchmarkable iszero(x)
SUITE["isone/2"] = @benchmarkable isone(o)
SUITE["lu"] = @benchmarkable lu(X)

SUITE["_solve"] = @benchmarkable Symbolics._solve(X, X, true)


text:
module TestExclude

include("preamble.jl")
using Base: PkgId
using Aqua: getobj, normalize_exclude, normalize_and_check_exclude, rootmodule,
    reprexclude

@assert parentmodule(Tuple) === Core
@assert parentmodule(foldl) === Base
@assert parentmodule(Some) === Base
@assert parentmodule(Broadcast.Broadcasted) === Base.Broadcast
@assert rootmodule(Broadcast.Broadcasted) === Base

@testset "roundtrip" begin
    @testset for x in [
        foldl
        Some
        Tuple
        Broadcast.Broadcasted
    ]
        @test getobj(normalize_exclude(x)) == x
    end

    @testset "$(repr(last(spec)))" for spec in [
        (PkgId(Base) => "Base.foldl")
        (PkgId(Base) => "Base.Some")
        (PkgId(Core) => "Core.Tuple")
        (PkgId(Base) => "Base.Broadcast.Broadcasted")
    ]
        @test normalize_exclude(getobj(spec)) === spec
    end
end

@testset "normalize_and_check_exclude" begin
    @testset "$i" for (i, exclude) in enumerate([
        [foldl],
        [foldl, Some],
        [foldl, Tuple],
    ])
        local specs
        @test (specs = normalize_and_check_exclude(exclude)) isa Vector
        @test Base.include_string(@__MODULE__, reprexclude(specs)) == specs
    end
end

end  # module


text:
function pidigits(N::Int, printOut::Bool)
    """
    See http://shootout.alioth.debian.org/u64q/performance.php?test=pidigits#about

    Transliterated from Mario Pernici Python's program

    INPUT:

    - N -- a positive integer giving the number of digits of pi to be computed

    - printOut -- a boolean specifying if we want intermediate printouts of digits in packets of 10

    OUTPUT:

    - returns the last ten digits anyway

    - prints all the digits in packets of 10 iff printOut == true

    """

    i = k = ns = 0
    k1 = 1
    n,a,d,t,u = map(BigInt,(1,0,1,0,0))

    while true
        k += 1
        t = n << 1
        n *= k
        a += t
        k1 += 2
        a *= k1
        d *= k1

        if a >= n
            t,u = divrem(n*3 +a, d)
            u += n
            if d > u
                ns = ns*10 + t
                i += 1
                if mod(i,10) == 0
                    if printOut
                        print(ns)
                        @printf("\t:%d\n", i)
                    end
                    if i >= N
                        return ns
                    end
                    ns = 0
                end
                a -= d*t
                a *= 10
                n *= 10

            end
        end
    end
end

perf_pidigits(N::Int=1000) = pidigits(N,false)


text:
module APIGenerator

export make_api

"""
    make_api(pkg, dest; title = "\$pkg API", readme=true, clobber=true, toc=true)

Write a markdown file (`dest`) which includes the docstring of every exported name from `pkg`.

- `clobber = true` will overwrite `dest`
- `readme = true` will include the `pkg`'s readme as the first item
- `toc = true` adds a table of contents with section links
"""
function make_api(pkg::String, dest::String; title="$pkg API", readme=true, clobber=true,
                  toc=true)
    if dest[end-2:end] != ".md"
        dest *= ".md"
    end
    Pkg = Symbol(pkg)
    clobber && isfile(dest) && rm(dest)
    touch(dest)
    file = open(dest, "r+")
    write(file, "<!--- Generated at " * string(now()) * ".  Don't edit --->\n\n")
    write(file, "# $title\n\n")
    nms = setdiff(names(@eval(Main.$Pkg)), [Pkg])
    if toc
        write(file, "# Contents\n")
        readme && write(file, "- [$Pkg](#$Pkg)\n")
        for nm in nms
            write(file, "- [$nm](#$(lowercase(string(nm))))\n")
        end
        write(file, "---\n\n")
    end
    info("The following items are included in the output file:\n")
    if readme
        info("Writing $pkg README")
        write(file, Markdown.plain(Markdown.readme(@eval(Main.$Pkg))) * "\n")
    end
    for nm in nms
        try
            @eval obj = Main.$Pkg.$nm
            d = Docs.doc(obj)
            println("  > $nm")
            write(file, "## $(string(nm))\n$(Markdown.plain(d))\n")
            toc && write(file, "[top](#contents)\n")
        end
    end
    close(file)
end

end # module


text:
@testset "examples" begin

D = ParserCombinator.Parsers.DOT

# from http://graphs.grevian.org/example

d = parse_dot(open(p -> read(p, String), "dot/simple.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e"])
@test edges(d) == Set([("a","e"),("c","d"),("c","e"),("b","c"),("a","c"),("a","b")])

d = parse_dot(open(p -> read(p, String), "dot/k6.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e", "f"])
@test length(edges(d)) == 6*5/2

d = parse_dot(open(p -> read(p, String), "dot/simple-digraph.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "d"])
@test edges(d) == Set([("a","b"),("b","c"),("c","d"),("d","a")])

d = parse_dot(open(p -> read(p, String), "dot/full-digraph.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "e"])
@test edges(d) == Set([("a","b"),("a","c"),("c","b"),("e","b"),("e","e"),("c","e")])
@test length(d.stmts) == 6
@test length(d.stmts[1].nodes) == 2
# id.id are attribues of NodeID and StringID
@test d.stmts[1].nodes[1].id.id == "a"
@test isnull(d.stmts[1].nodes[1].port)
@test d.stmts[1].nodes[2].id.id == "b"
@test isnull(d.stmts[1].nodes[2].port)
@test length(d.stmts[1].attrs) == 2
@test d.stmts[1].attrs[1].name.id == "label"
@test d.stmts[1].attrs[1].value.id == "0.2"
@test d.stmts[1].attrs[2].name.id == "weight"
@test d.stmts[1].attrs[2].value.id == "0.2"

for f in ("path1.dot", "path2.dot")
    d = parse_dot(open(p -> read(p, String), "dot/$f"))[1]
    @test nodes(d) == Set(["a", "b", "c", "d", "e", "f"])
    @test edges(d) == Set([("a","b"),("a","d"),("b","c"),("b","d"),("c","d"),("c","f"),("e","f"),("d","e")])
end

# this example also has comments (ignored in parse)
d = parse_dot(open(p -> read(p, String), "dot/simple-subgraph.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "d", "f"])
@test edges(d) == Set([("a","b"),("a","f"),("b","c"),("c","d"),("f","c")])
@test length(d.stmts) == 2
@test isa(d.stmts[1], SubGraph)
@test !isnull(d.stmts[1].id)
@test get(d.stmts[1].id).id == "cluster_0"
@test d.stmts[1].stmts[1].name.id == "label"
@test d.stmts[1].stmts[1].value.id == "Subgraph A"

d = parse_dot(open(p -> read(p, String), "dot/bipartite-subgraph.dot"))[1]
@test nodes(d) == Set(["a", "b", "c", "d", "e"])
@test length(edges(d)) == 6

d1 = parse_dot(open(p -> read(p, String), "dot/large1.dot"))[1]
@test length(nodes(d1)) == 21
@test length(edges(d1)) == 42

d2 = parse_dot(open(p -> read(p, String), "dot/large2.dot"))[1]
@test length(nodes(d2)) == 21
@test length(edges(d2)) == 42

@test nodes(d1) == nodes(d2)
@test edges(d1) == edges(d2)

s = open(p -> read(p, String), "dot/tictactoe.dot")
@time d = parse_dot(s)[1]
@time n = nodes(d)
@time e = edges(d)
# these numbers are not really checks, just what the code returns
@test length(n) == 744
@test length(e) == 4361

d = parse_dot(open(p -> read(p, String), "dot/double.dot"))
d1, d2 = d
@test nodes(d1) == Set(["a", "b", "c", "d", "e"])
@test edges(d1) == Set([("a","e"),("c","d"),("c","e"),("b","c"),("a","c"),("a","b")])
@test nodes(d2) == Set(["a", "b", "c", "d", "e"])
@test edges(d2) == Set([("a","e"),("c","d"),("c","e"),("b","c"),("a","c"),("a","b")])

println("examples ok")

end


text:
module Shoot
    include("load2.jl")
    include("load2_module.jl")
    include("shoot_module.jl")
    using DifferentialEquations
    using .Load2
    using .Shoot_module

    const DELV = 1.0E-7
    const NVAR = 2
    const V0 = -1.588071022611375

    function construct(dx, data)
        load2_param, load2_val = Load2.construct()
        shoot_val = Shoot_module.Shoot_module_variables(
            dx,
            data.eps,
            data.grid_num,
            data.usethread,
            V0,
            Load2.make_vmax(data.xmax, load2_param, load2_val),
            data.xmin,
            data.x_matching_point,
            data.xmax,
            Shoot_module.Load2_module.Load2_module_param(load2_param.LAMBDA, load2_param.THRESHOLD, load2_param.K),
            Shoot_module.Load2_module.Load2_module_variables(load2_val.x, load2_val.yitp))

        return shoot_val
    end

    createresult(xarray1, xarray2, yarray1, yarray2) = let
        xarray = append!(xarray1[2:end], xarray2[2:end])
        yarray1[end] = (yarray1[end] + yarray2[1]) / 2.0
        yarray = append!(yarray1[2:end], yarray2[2:end])

        xar = [0.0]
        yar = [1.0]
        return append!(xar, xarray), append!(yar, yarray)
    end

    function funcx1!(dfdv, f1, f_vector, shoot_val)
        sav = shoot_val.v0
        shoot_val.v0 += DELV

        u0 = load1(shoot_val.xmin, shoot_val.v0)
        tspan = (shoot_val.xmin, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps)
        f = sol.u[end]

        # NVAR個の合致条件にある偏微分を数値的に計算
        for i in 1:NVAR
            dfdv[i, 1] = (f[i] - f1[i]) / DELV
        end

        # 境界におけるパラメータを格納
        shoot_val.v0 = sav
    end

    function funcx2!(dfdv, f2, f_vector, shoot_val)
        sav = shoot_val.vmax
        shoot_val.vmax += DELV

        u0 = Load2.load2(shoot_val.xmax, shoot_val.vmax, shoot_val.load2_param, shoot_val.load2_val)
        tspan = (shoot_val.xmax, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps)
        f = sol.u[end]

        for i in 1:NVAR
            dfdv[i, 2] = (f2[i] - f[i]) / DELV
        end

        shoot_val.vmax = sav
    end

    load1(xmin, v0) = let
        y = Vector{Float64}(undef, 2)

        # y[1] = 1.0 + v0[0] * xmin + 4.0 / 3.0 * xmin * sqrt(xmin) + 0.4 * v0[0] * xmin * xmin * sqrt(xmin) + 1.0 / 3.0 * xmin * xmin * xmin
        y[1] = (((1.0 / 3.0 * xmin + 0.4 * v0 * sqrt(xmin)) * xmin) + 4.0 / 3.0 * sqrt(xmin) + v0) * xmin + 1.0
        # y[2] = v0[0] + 2.0 * sqrt(xmin) + v0[0] * xmin * sqrt(xmin) + xmin * xmin + 0.15 * v0[0] * xmin * xmin * sqrt(xmin)
        y[2] = ((0.15 * v0 * sqrt(xmin) + 1.0) * xmin + v0 * sqrt(xmin)) * xmin + 2.0 * sqrt(xmin) + v0

        return y
    end

    function shootf!(shoot_val)
        f_vector(u, p, t) = [u[2], u[1] * sqrt(u[1] / t)]

        # 最良の仮の値v1_でx1からxfまで解いていく
        u0 = load1(shoot_val.xmin, shoot_val.v0)
        tspan = (shoot_val.xmin, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps)
        f1 = sol.u[end]

        # 最良の仮の値v2_でx2からxfまで解いていく
        u0 = Load2.load2(shoot_val.xmax, shoot_val.vmax, shoot_val.load2_param, shoot_val.load2_val)
        tspan = (shoot_val.xmax, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps)
        f2 = sol.u[end]

        dfdv = Matrix{Float64}(undef, NVAR, NVAR)

        if shoot_val.usethread
            # x1で用いる境界条件を変える
            fa = Threads.@spawn funcx1!(dfdv, f1, f_vector, shoot_val)

            # 次にx2で用いる境界条件を変える
            fb = Threads.@spawn funcx2!(dfdv, f2, f_vector, shoot_val)

            fetch(fa)
            fetch(fb)
        else
            # x1で用いる境界条件を変える
            funcx1!(dfdv, f1, f_vector, shoot_val)

            # 次にx2で用いる境界条件を変える
            funcx2!(dfdv, f2, f_vector, shoot_val)
        end

        f = f1 .- f2
        ff = -1.0 .* f

        solf = dfdv \ ff

        # x1の境界でのパラメータ値の増分
        shoot_val.v0 += solf[1]

        # x2の境界でのパラメータ値の増分
        shoot_val.vmax += solf[2]

        size1 = floor(Int, shoot_val.xf / shoot_val.dx) + 1
        size2 = floor(Int, (shoot_val.xmax - shoot_val.xf) / shoot_val.dx) + 1

        xarray1 = Vector{Float64}(undef, size1)
        yarray1 = Vector{Float64}(undef, size1)
        xarray2 = Vector{Float64}(undef, size2)
        yarray2 = Vector{Float64}(undef, size2)

        if shoot_val.usethread
            fa = Threads.@spawn solveode_xmintoxf(f_vector, shoot_val)
            fb = Threads.@spawn solveode_xmaxtoxf(f_vector, shoot_val)
            xarray1, yarray1 = fetch(fa)
            xarray2, yarray2 = fetch(fb)
        else
            xarray1, yarray1 = solveode_xmintoxf(f_vector, shoot_val)
            xarray2, yarray2 = solveode_xmaxtoxf(f_vector, shoot_val)
        end

        return createresult(xarray1, xarray2, yarray1, yarray2)
    end

    solveode_xmaxtoxf(f_vector, shoot_val) = let
        # 得られた条件でx2...xfまで微分方程式を解く
        u0 = Load2.load2(shoot_val.xmax, shoot_val.vmax, shoot_val.load2_param, shoot_val.load2_val)
        tspan = (shoot_val.xmax, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps, saveat = shoot_val.dx)
        a = vcat(sol.u...)
        yarray = deleteat!(a, 2:2:length(a))   # shoot_val.xmax...xfの結果を得る

        return reverse(sol.t), reverse(yarray)
    end

    solveode_xmintoxf(f_vector, shoot_val) = let
        # 得られた条件でx1...dxまで微分方程式を解く
        u0 = load1(shoot_val.xmin, shoot_val.v0)
        tspan = (shoot_val.xmin, shoot_val.dx)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps)
        y = sol.u[end]
        tarray = sol.t[1:1]

        # 得られた条件でdx...xfまで微分方程式を解く
        u0 = y
        tspan = (shoot_val.dx, shoot_val.xf)
        prob = ODEProblem(f_vector, u0, tspan)
        sol = solve(prob, VCABM(), abstol = shoot_val.eps, reltol = shoot_val.eps, saveat = shoot_val.dx)
        a = vcat(sol.u...)
        f = deleteat!(a, 2:2:length(a))

        y = load1(shoot_val.xmin, shoot_val.v0)
        yarray = append!(y[1:1], f)
        tarray = append!(tarray, sol.t)

        return tarray, yarray
    end
end


text:
# Acceleration index types for one or more columns

abstract type AbstractIndex{names}; end
abstract type AbstractUniqueIndex{names} <: AbstractIndex{names}; end

# By default, indices will be invalidated if some of the indexed columns are removed
@inline function (::Project{n})(i::AbstractIndex{names}) where {names, n}
    if _issubset(names, n)
        return i
    else
        return NoIndex()
    end
end

struct NoIndex <: AbstractIndex{()}
end
copy(index::NoIndex) = index

(::Rename{oldnames, newnames})(::NoIndex) where {oldnames, newnames} = NoIndex()

# This index states that the indicated columns have unique values according to `isequal`
struct UniqueIndex{names} <: AbstractIndex{names}
end
copy(index::UniqueIndex) = index

function (::Rename{oldnames, newnames})(::UniqueIndex{names}) where {oldnames, newnames, names}
    return UniqueIndex{_rename(Val(oldnames), Val(newnames), Val(names))}()
end

# Lexicographically ordered
struct SortIndex{names, V <: AbstractVector{Int}} <: AbstractIndex{names}
    order::V # The row indices in sort order
end
SortIndex{names}(v::V) where {names, V} = SortIndex{names, V}(v)
copy(index::SortIndex{names}) where {names} = SortIndex{names}(copy(index.order))

function (::Rename{oldnames, newnames})(index::SortIndex{names}) where {oldnames, newnames, names}
    return SortIndex{_rename(Val(oldnames), Val(newnames), Val(names))}(index.order)
end

@inline function (::Project{n})(i::SortIndex{names}) where {names, n}
    ns = _headsubset(Val(names), Val(n))
    if ns === ()
        return NoIndex()
    else # We retain lexicographical ordering
        return SortIndex{ns}(i.order)
    end
end

# Lexicographically ordered and unique
struct UniqueSortIndex{names, V <: AbstractVector{Int}} <: AbstractUniqueIndex{names}
    order::V # The row indices in sort order
end
UniqueSortIndex{names}(v::V) where {names, V} = UniqueSortIndex{names, V}(v)
copy(index::UniqueSortIndex{names}) where {names} = UniqueSortIndex{names}(copy(index.order))

function (::Rename{oldnames, newnames})(index::UniqueSortIndex{names}) where {oldnames, newnames, names}
    return UniqueSortIndex{_rename(Val(oldnames), Val(newnames), Val(names))}(index.order)
end

@inline function (::Project{n})(i::UniqueSortIndex{names}) where {names, n}
    ns = _headsubset(Val(names), Val(n))
    if ns === ()
        return NoIndex()
    elseif ns === names
        return i
    else # We drop uniqueness, but maintain lexicographical ordering
        return SortIndex{ns}(i.order)
    end
end

# Hash table acceleration index to unknown number of rows
struct HashIndex{names, D <: AbstractDict{<:Tuple, <:AbstractVector{Int}}} <: AbstractIndex{names}
    dict::D # Mapping from column values to list of matching indices
end
HashIndex{names}(d::D) where {names, D} = HashIndex{names, D}(d)
copy(index::HashIndex{names}) where {names} = HashIndex{names}(copy(index.d))

function (::Rename{oldnames, newnames})(index::HashIndex{names}) where {oldnames, newnames, names}
    return HashIndex{_rename(Val(oldnames), Val(newnames), Val(names))}(index.dict)
end

# Hash table acceleration index to unique rows
struct UniqueHashIndex{names, D <: AbstractDict{<:Tuple, Int}} <: AbstractUniqueIndex{names}
    dict::D # Mapping from column values to unique matching index
end
UniqueHashIndex{names}(d::D) where {names, D} = UniqueHashIndex{names, D}(d)
copy(index::UniqueHashIndex{names}) where {names} = UniqueHashIndex{names}(copy(index.d))

function (::Rename{oldnames, newnames})(index::UniqueHashIndex{names}) where {oldnames, newnames, names}
    return UniqueHashIndex{_rename(Val(oldnames), Val(newnames), Val(names))}(index.dict)
end


text:
using Test
using SymDiff

@test simplify(:((x+0)/(y+0))) == :(x/y)


text:
include("utils.jl")

@testset "Circle" begin
    M = Circle()
    @testset "Real Circle Basics" begin
        @test representation_size(M) == ()
        @test manifold_dimension(M) == 1
        @test !is_manifold_point(M, 9.)
        @test_throws DomainError is_manifold_point(M, 9., true)
        @test !is_tangent_vector(M, 9., 0.)
        @test_throws DomainError is_tangent_vector(M, 9., 0., true)
        @test flat(M,0.0, FVector(TangentSpace,1.0)) == FVector(CotangentSpace,1.0)
        @test sharp(M,0.0, FVector(CotangentSpace,1.0)) == FVector(TangentSpace,1.0)
        @test vector_transport_to(M,0.0,1.0,1.0, ParallelTransport()) == 1.0
        @test retract(M,0.0,1.0) == exp(M,0.0,1.0)
        @test injectivity_radius(M) ≈ π
        @test mean(M, [-π/2,0.,π]) ≈ π/2
        @test mean(M, [-π/2,0.,π],[1., 1., 1.]) == π/2
        v = MVector(0.0)
        x = SVector(0.0)
        log!(M,v,x,SVector(π/4))
        @test norm(M,x,v) ≈ π/4
        @test is_tangent_vector(M,x,v)
        @test project_point(M,1.0) == 1.0
        x = MVector(0.0)
        project_point!(M,x)
        @test x == MVector(0.0)
        x .+= 2*π
        project_point!(M,x)
        @test x == MVector(0.0)
    end
    types = [Float64, Float32]

    basis_types = (ArbitraryOrthonormalBasis(),)
    basis_types_real = (ArbitraryOrthonormalBasis(),
        DiagonalizingOrthonormalBasis([-1]),
        DiagonalizingOrthonormalBasis([1])
    )
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [-π/4,0.,π/4])
            test_manifold(
                M,
                pts,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vector_spaces = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                is_mutating = false
            )
            ptsS = SVector.(pts)
            test_manifold(
                M,
                ptsS,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                basis_types_vecs = basis_types_real,
                basis_types_to_from = basis_types_real
            )
        end
    end
    Mc = Circle(ℂ)
    @testset "Complex Circle Basics" begin
        @test representation_size(Mc) == ()
        @test manifold_dimension(Mc) == 1
        @test is_tangent_vector(Mc, 1im, 0.)
        @test is_manifold_point(Mc, 1im)
        @test !is_manifold_point(Mc, 1+1im)
        @test_throws DomainError is_manifold_point(Mc, 1+1im, true)
        @test !is_tangent_vector(Mc, 1+1im, 0.)
        @test_throws DomainError is_tangent_vector(Mc, 1+1im, 0., true)
        @test !is_tangent_vector(Mc, 1im, 2im)
        @test_throws DomainError is_tangent_vector(Mc, 1im, 2im, true)
        @test flat(Mc,0.0+0.0im, FVector(TangentSpace,1.0im)) == FVector(CotangentSpace,1.0im)
        @test sharp(Mc,0.0+0.0im, FVector(CotangentSpace,1.0im)) == FVector(TangentSpace,1.0im)
        @test norm(Mc,1.0,log(Mc, 1.0,-1.0)) ≈ π
        @test is_tangent_vector(Mc,1.0,log(Mc,1.0,-1.0))
        v = MVector(0.0+0.0im)
        x = SVector(1.0+0.0im)
        log!(Mc,v,x,SVector(-1.0+0.0im))
        @test norm(Mc,SVector(1.0),v) ≈ π
        @test is_tangent_vector(Mc,x,v)
        @test project_point(Mc,1.0) == 1.0
        project_point(Mc,1/sqrt(2.0) + 1/sqrt(2.0) * im) == 1/sqrt(2.0) + 1/sqrt(2.0) * im
        x = MVector(1.0+0.0im)
        project_point!(Mc,x)
        @test x == MVector(1.0+0.0im)
        x .*= 2
        project_point!(Mc,x)
        @test x == MVector(1.0+0.0im)

    end
    types = [Complex{Float64}, Complex{Float32}]
    for T in types
        @testset "Type $T" begin
            a = 1/sqrt(2.0)
            pts = convert.(Ref(T), [a-a*im, 1+0im, a+a*im])
            test_manifold(
                Mc,
                pts,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vector_spaces = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                is_mutating=false,
                exp_log_atol_multiplier = 2.0,
                is_tangent_atol_multiplier = 2.0
            )
            ptsS = SVector.(pts)
            test_manifold(
                Mc,
                ptsS,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                exp_log_atol_multiplier = 2.0,
                is_tangent_atol_multiplier = 2.0,
                basis_types_vecs = basis_types,
                basis_types_to_from = basis_types
            )
        end
    end
end


text:
export plaintext,
       nowrap,
       font,
       fontfamily,
       serif,
       sansserif,
       slabserif,
       monospace,
       normal,
       slanted,
       italic,
       uppercase,
       lowercase,
       xxsmall,
       xsmall,
       small,
       medium,
       large,
       xlarge,
       xxlarge,
       bold,
       bolder,
       lighter,
       textalign,
       raggedright,
       raggedleft,
       justify,
       centertext,
       fontweight,
       heading,
       h1, h2, h3, h4,
       paragraph,
       emph,
       codeblock,
       code,
       blockquote,
       caption,
       title,
       menufont

plaintext(x) = div(x)

# nowrap

@api nowrap => NoTextWrap <: Tile begin
    arg(tile::Tile)
end

# Font properties

abstract FontProperty

@api font => WithFont{T <: FontProperty} <: Tile begin
    typedarg(prop::T)
    curry(tile::TileList)
end
WithFont{T}(prop::T, t) = WithFont{T}(prop, t)

font(t::Union(Tile, String), args::FontProperty...) =
    foldr(font, t, props)

font(props::FontProperty...) =
    t -> foldr(font, t, props)

render(t::WithFont) =
    addclasses(wrapmany(t.tile, :span), classes(t))

type FontFamily <: FontProperty
    family::String
end
fontfamily(f) = FontFamily(f)

render(t::WithFont{FontFamily}) =
    render(t.tile) & [:style => [:fontFamily => t.prop.family]]

abstract FontType <: FontProperty

@terms FontType begin
    serif => Serif
    sansserif => SansSerif
    slabserif => SlabSerif
    monospace => Monospace
end
classes(::WithFont{Serif}) = "font-serif"
classes(::WithFont{SansSerif}) = "font-sansserif"
classes(::WithFont{SlabSerif}) = "font-serif font-slab"
classes(::WithFont{Monospace}) = "font-monospace"

abstract FontStyle <: FontProperty

@terms FontStyle begin
    normal => Normal
    slanted => Slanted
    italic => Italic
end
classes(::WithFont{Normal}) = "font-normal"
classes(::WithFont{Slanted}) = "font-slanted"
classes(::WithFont{Italic}) = "font-italic"

abstract FontCase <: FontProperty

@terms FontCase begin
    ucase => Uppercase
    lcase => Lowercase
end
classes(::WithFont{Uppercase}) = "font-uppercase"
classes(::WithFont{Lowercase}) = "font-lowercase"

abstract FontSize <: FontProperty
@terms FontSize begin
    xxsmall => XXSmall
    xsmall => XSmall
    small => Small
    medium => Medium
    large => Large
    xlarge => XLarge
    xxlarge => XXLarge
end
classes(::WithFont{XXSmall}) = "font-xx-small"
classes(::WithFont{XSmall}) = "font-x-small"
classes(::WithFont{Small}) = "font-small"
classes(::WithFont{Medium}) = "font-medium"
classes(::WithFont{Large}) = "font-large"
classes(::WithFont{XLarge}) = "font-x-large"
classes(::WithFont{XXLarge}) = "font-xx-large"

immutable AbsFontSize <: FontProperty
    size::Length
end
fontsize(size::Length) = AbsFontSize(size)

render(t::WithFont{AbsFontSize}) =
    render(t.tile) & [:style => [:fontSize => t.prop.size]]

abstract FontWeight <: FontProperty

@terms FontWeight begin
    bold => Bold
    bolder => Bolder
    lighter => Lighter
end
classes(::WithFont{Bold}) = "font-bold"
classes(::WithFont{Bolder}) = "font-bolder"
classes(::WithFont{Lighter}) = "font-lighter"

immutable NumericFontWeight{n} <: FontProperty
end

render{n}(t::WithFont{NumericFontWeight{n}}) =
    render(t.tile) & [:style => [:fontWeight => n]]

const allowed_font_weights = 100:100:900

function fontweight(w::Integer)
    if !(w in allowed_font_weights)
        error(string(w, " is not an allowed font weight"))
    end
    NumericFontWeight{w}()
end

abstract TextAlignment

@terms TextAlignment begin
    raggedright => RaggedRight
    raggedleft => RaggedLeft
    justifytext => JustifyText
    centertext => CenterText
end

@api textalign => AlignText{T <: TextAlignment} begin
    typedarg(alignment::T)
    curry(tile::Tile)
end

render(t::AlignText{RaggedRight}) =
    render(t.tile) & [:style => [:textAlign => :left]]
render(t::AlignText{RaggedLeft}) =
    render(t.tile) & [:style => [:textAlign => :right]]
render(t::AlignText{JustifyText}) =
    render(t.tile) & [:style => [:textAlign => :justify]]
render(t::AlignText{CenterText}) =
    render(t.tile) & [:style => [:textAlign => :center]]

# Themable fonts

heading(n::Int, txt) = class("heading-$n", txt, forcewrap=true, wrap=:h1)

h1(txt) = heading(1, txt)
h2(txt) = heading(2, txt)
h3(txt) = heading(3, txt)
h4(txt) = heading(4, txt)

paragraph(txt) = class("paragraph", txt, forcewrap=true, wrap=:p)
blockquote(txt) = class("blockquote", txt, forcewrap=true, wrap=:blockquote)
caption(txt) = class("caption", txt, wrap=:span)
emph(txt) = class("emph", txt, forcewrap=true, wrap=:em)

@api code => Code <: Tile begin
    arg(code::Any)
    kwarg(language::String="julia")
end

render(x::Code) = Elem(:code, x.code)

@api codeblock => CodeBlock <: Tile begin
    arg(code::Any)
    kwarg(language::String="julia")
end

render(x::CodeBlock) = Elem(:pre, x.code)

# colsize


text:
# -*- coding: utf-8 -*-
using Plots
using Random
using ProgressMeter
using Revise

# +
using Uchiyama

n = 50 # number of particles
ϵ = 0.02
c = trunc(Int, 200ϵ) # marker size

rng = MersenneTwister(1234)

particles = SquareParticles(rng, n, ϵ, option = :box)
pc = ParticleCollisions(particles)
bc = BoxCollisions(particles)

# +

steps = 1000
pbar = Progress(steps)

anim = @animate for _ in 1:steps

     dt = step!(particles, pc, bc)

     p = plot(size  = (200,200),
              xlims = (-0.5,0.5),
              ylims = (-0.5,0.5),
              grid  = false,
              axis  = nothing, framestyle= :box, legend=false, widen = false)

     scatter!( getindex.(particles.q,1) .- getindex.(particles.q,2),
               getindex.(particles.q,2) .+ getindex.(particles.q,1),
               markershape  = :square,
               markersize   = c,
               aspect_ratio = :equal)
     next!(pbar)

end every 10
# -

gif(anim, joinpath(@__DIR__, "event_driven.gif"), fps = 10)


text:
####################
## CIR model type ##
####################

type CIR
    alpha::Float64
    mu::Float64
    sigma::Float64
    scale::DatePeriod # which period is Δt = 1
end

## convenient constructors
##------------------------

function CIR(x::Array)
    return CIR(x[1], x[2], x[3], Dates.Year(1))
end

function CIR(x1, x2, x3)
    return CIR(x1, x2, x3, Dates.Year(1))
end

## customized display
##-------------------

import Base.display
function display(cirMod::CIR)
    println("α:        $(cirMod.alpha)")
    println("μ:        $(cirMod.mu)")
    println("σ:        $(cirMod.sigma)")
    println("scaling:  Δt = $(cirMod.scale)")
end

## information retrieval
##----------------------

function getParams(cirMod::CIR)
    return (cirMod.alpha, cirMod.mu, cirMod.sigma)
end

#################
## likelihoods ##
#################

function cirNllh(da::DataArray, params::Array{Float64, 1},
                 Δt::Float64 = 1)
    # Δt means time step of data relativ to time step represented by
    # parameters

    # calculate negative log-likelihood for given parameter values
    α, μ, σ = params

    # get constant values
    c = 2*α/(σ^2*(1-exp(-α*Δt)))
    q = 2*α*μ/σ^2 - 1

    # get changing values u and v
    uWithNA = c*da[1:(end-1)]*exp(-α*Δt)
    vWithNA = c*da[2:end]

    # eliminate observations with missing value for either u or v
    notElim = (!isna(uWithNA)) & (!isna(vWithNA))
    u = uWithNA[notElim]
    v = vWithNA[notElim]

    nObs = size(u, 1)

    return -(nObs)*log(c) +
    sum(u + v - 0.5*q*log(v./u) - log(besseli(q, 2*sqrt(u.*v))))
end

function cirNllh(da::DataArray, cirMod::CIR, Δt::Float64 = 1)
    return cirNllh(da, [getParams(cirMod)...], Δt)
end

function cirNllhx(da::DataArray, params::Array{Float64, 1},
                  Δt::Float64 = 1)
    # calculate negative log-likelihood for given parameter values
    # using scaled modified bessel function of first order
    α, μ, σ = params

    # get constant values
    c = 2*α/(σ^2*(1-exp(-α*Δt)))
    q = 2*α*μ/σ^2 - 1

    # get changing values u and v
    uWithNA = c*da[1:(end-1)]*exp(-α*Δt)
    vWithNA = c*da[2:end]

    # eliminate observations with missing value for either u or v
    notElim = (!isna(uWithNA)) & (!isna(vWithNA))
    u = uWithNA[notElim]
    v = vWithNA[notElim]

    nObs = size(u, 1)

    return -(nObs)*log(c) +
    sum(u + v - 0.5*q*log(v./u) -
        log(besselix(q, 2*sqrt(u.*v))) - 2*sqrt(u.*v))
end

function cirNllhx(da::DataArray, cirMod::CIR, Δt::Float64 = 1)
    return cirNllhx(da, [getParams(cirMod)...], Δt)
end

####################
## OLS estimation ##
####################

function cirOls(x::DataArray, Δt::Float64=1)
    y = (x[2:end] .- x[1:(end-1)])./sqrt(x[1:(end-1)])
    x1 = 1 ./ sqrt(x[1:(end-1)])
    x2 = sqrt(x[1:(end-1)])

    # conduct regression
    df = DataFrame(y = y, x1 = x1, x2 = x2)
    olsFit = lm(y ~ 0 + x1 + x2, df)

    a1, a2 = coef(olsFit)

    # get associated parameters
    α = -a2 / Δt
    μ = a1 / (Δt * α)
    σ = sqrt(var(residuals(olsFit)) /Δt)
    return (α, μ, σ)
end

##############################
## conditional distribution ##
##############################

function getCondDistr(cirMod::CIR, r0::Float64,  Δt::Float64=1.)
    # conditional cdf given r
    α, μ, σ = cirMod.alpha, cirMod.mu, cirMod.sigma

    df = 4*α*μ/σ^2
    c = 2*α/(σ^2*(1-exp(-α*Δt)))
    λ = 2*c*r0*exp(-α*Δt)
    return (df, λ, c)
end

import Distributions.cdf
function cdf(cirMod::CIR, x::Float64, r0::Float64, Δt::Float64=1.)
    df, λ, c = getCondDistr(cirMod, r0, Δt)

    return pnchisq(x*2*c, df, λ)
end

function cdf(cirMod::CIR, x::Array{Float64, 1}, Δt::Float64=1.)
    ## take first value as given input
    nObs = length(x)-1
    condDistrParams = [getCondDistr(cirMod, x[ii], Δt) for ii=1:nObs]

    u = zeros(nObs)
    for ii=1:nObs
        df, λ, c = condDistrParams[ii]
        u[ii] = pnchisq(x[ii+1]*2*c, df, λ)
    end

    return u
end

function quantile(cirMod::CIR, q::Float64, r0::Float64, Δt::Float64=1.)
    df, λ, c = getCondDistr(cirMod, r0, Δt)

    return qnchisq(q, df, λ)/(2*c)
end

####################
## simulating CIR ##
####################

function simulate(cirMod::CIR, r0::Float64, nObs::Int, Δt::Float64=1.)
    α, μ, σ = cirMod.alpha, cirMod.mu, cirMod.sigma

    df = 4*α*μ/σ^2 # constant
    c = 2*α/(σ^2*(1-exp(-α*Δt))) # constant for given Δt

    vals = [r0; Array(Float64, nObs)]
    for ii=1:nObs
        # get conditional non-centrality
        λ = 2*c*vals[ii]*exp(-α*Δt)
        vals[ii+1] = rnchisq(df, λ)/(2*c)
    end
    return vals
end


text:
# For numbers of a magnitude less than ~1000
function almostEqual(iterable1, iterable2, allowedError=0.0000000000001)
    if size(iterable1, 1) != size(iterable2, 1)
        s1 = size(iterable1, 1)
        s2  = size(iterable2, 1)
        println("Size mismatch: $s1 vs $s2")
        println(iterable1)
        println(iterable2)
        return false
    end

    if ndims(iterable1) >= 2
        return matrixAlmostEqual(iterable1, iterable2, allowedError)
    end

    for i in 1:size(iterable1, 1)
        if abs(iterable1[i] - iterable2[i]) > allowedError
            error = iterable1[i] - iterable2[i]
            println("Error: mismatch of $error at position $i")
            println(iterable1)
            println(iterable2)
            return false
        end
    end

    return true
end

function matrixAlmostEqual(m1, m2, allowedError=0.0000000000001)
    s1 = size(m1, 1)

    if size(m1) != size(m2)
        return false
    end

    for i in 1:s1
        if almostEqual(m1[i,:], m2[i,:], allowedError) != true
            return false
        end
    end

    return true
end

max_stars_repo_name: henrystoldt/fvCFD


text:
using TestCall
erfv(rand(3))
a = Vector{Float64}(undef, 3)
erfv2(a, rand(3))

max_stars_repo_name: el-oso/TestCall

max_issues_repo_name: el-oso/TestCall

max_forks_repo_name: el-oso/TestCall


text:
abstract RoaringContainer

# Constants for containers.

## Array containers.
const ARRAY_DEFAULT_MAX_SIZE = 4096
const ARRAY_DEFAULT_SIZE = 16

## Bitset containers.
const BITSET_CONTAINER_SIZE_IN_WORDS = (1 << 16) / 64

# Container type definitions and methods.
include("array_container.jl")
include("bitset_container.jl")


text:
module vehicle_platoon_10
using ReachabilityAnalysis
A = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1.702423734 3.929356551 -4.3607983776 -1.01374489 -1.6167727749 0.2653009364 -0.2375199245 -0.4793543458 0.06412815 -0.1079326841 -0.2463610381 0.0276872161 -0.0605561959 -0.1501445039 0.0151944922 -0.0374830081 -0.0986391305 0.009628751 -0.0242136837 -0.0665592904 0.0067836913 -0.015601062 -0.0442510048 0.0052325207 -0.0093924696 -0.0272127915 0.0043984935 -0.0044278796 -0.0129879863 0.0040303349;
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0.688678844 2.3125837761 0.2653009364 1.4649038095 3.4500022052 -4.2966702275 -1.1216775741 -1.863133813 0.2929881525 -0.2980761204 -0.6294988497 0.0793226422 -0.1454156921 -0.3450001686 0.0373159671 -0.0847698796 -0.2167037943 0.0219781835 -0.0530840701 -0.1428901352 0.0148612718 -0.0336061533 -0.0937720819 0.0111821848 -0.0200289416 -0.057238991 0.0092628557 -0.0093924696 -0.0272127915 0.0084288284;
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0.4511589195 1.8332294303 0.06412815 0.5807461599 2.066222738 0.2929881525 1.4043476136 3.2998577013 -4.2814757354 -1.1591605822 -1.9617729435 0.3026169036 -0.3222898041 -0.6960581401 0.0861063336 -0.1610167541 -0.3892511733 0.0425484878 -0.0941623492 -0.2439165858 0.026376677 -0.0575119497 -0.1558781215 0.0188916067 -0.0336061533 -0.0937720819 0.0152125197 -0.015601062 -0.0442510048 0.0136613491;
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0.3432262354 1.5868683922 0.0276872161 0.3906027236 1.6830849264 0.0793226422 0.5432631518 1.9675836075 0.3026169036 1.3801339299 3.2332984109 -4.274692044 -1.1747616442 -2.0060239482 0.3078494243 -0.3316822737 -0.7232709316 0.090504827 -0.1654446337 -0.4022391596 0.0465788228 -0.0941623492 -0.2439165858 0.0304070119 -0.0530840701 -0.1428901352 0.0232901001 -0.0242136837 -0.0665592904 0.0204450405;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0.2826700395 1.4367238883 0.0151944922 0.3057432273 1.4882292617 0.0373159671 0.3663890398 1.616525636 0.0861063336 0.5276620899 1.9233326028 0.3078494243 1.3707414603 3.2060856194 -4.2702935506 -1.1791895238 -2.0190119345 0.3118797592 -0.3316822737 -0.7232709316 0.094535162 -0.1610167541 -0.3892511733 0.0509773162 -0.0847698796 -0.2167037943 0.0356395326 -0.0374830081 -0.0986391305 0.0300737915;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0;
0.2451870315 1.3380847578 0.009628751 0.2584563558 1.3701645979 0.0219781835 0.2901421653 1.443978257 0.0425484878 0.3569965702 1.5893128445 0.090504827 0.5232342102 1.9103446165 0.3118797592 1.3707414603 3.2060856194 -4.2662632156 -1.1747616442 -2.0060239482 0.3162782527 -0.3222898041 -0.6960581401 0.0997676827 -0.1454156921 -0.3450001686 0.0577610076 -0.0605561959 -0.1501445039 0.0452682837;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0;
0.2209733477 1.2715254674 0.0067836913 0.2295859695 1.2938337531 0.0148612718 0.2490638862 1.3429518064 0.026376677 0.2857142857 1.4309902707 0.0465788228 0.3569965702 1.5893128445 0.094535162 0.5276620899 1.9233326028 0.3162782527 1.3801339299 3.2332984109 -4.2610306949 -1.1591605822 -1.9617729435 0.323061944 -0.2980761204 -0.6294988497 0.1093964337 -0.1079326841 -0.2463610381 0.0729554998;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0;
0.2053722857 1.2272744627 0.0052325207 0.2115808781 1.2443126759 0.0111821848 0.2251580898 1.2808457668 0.0188916067 0.2490638862 1.3429518064 0.0304070119 0.2901421653 1.443978257 0.0509773162 0.3663890398 1.616525636 0.0997676827 0.5432631518 1.9675836075 0.323061944 1.4043476136 3.2998577013 -4.2514019439 -1.1216775741 -1.863133813 0.3382564362 -0.2375199245 -0.4793543458 0.1370836498;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0;
0.1959798161 1.2000616712 0.0043984935 0.2009444061 1.2142864764 0.0092628557 0.2115808781 1.2443126759 0.0152125197 0.2295859695 1.2938337531 0.0232901001 0.2584563558 1.3701645979 0.0356395326 0.3057432273 1.4882292617 0.0577610076 0.3906027236 1.6830849264 0.1093964337 0.5807461599 2.066222738 0.3382564362 1.4649038095 3.4500022052 -4.2237147278 -1.01374489 -1.6167727749 0.4023845862;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1;
0.1915519365 1.1870736849 0.0040303349 0.1959798161 1.2000616712 0.0084288284 0.2053722857 1.2272744627 0.0136613491 0.2209733477 1.2715254674 0.0204450405 0.2451870315 1.3380847578 0.0300737915 0.2826700395 1.4367238883 0.0452682837 0.3432262354 1.5868683922 0.0729554998 0.4511589195 1.8332294303 0.1370836498 0.688678844 2.3125837761 0.4023845862 1.702423734 3.929356551 -3.9584137913;
]

function model(X0)
    S = @system(x' = Ax)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
NIScope_Vertical_Dict = Dict()
# Set trigger source, niScope channels start from 0, so subtract 1
NIScope_Vertical_Dict["1"]   = "0"
NIScope_Vertical_Dict["2"]   = "1"
NIScope_Vertical_Dict["3"]   = "2"
NIScope_Vertical_Dict["4"]   = "3"
NIScope_Vertical_Dict["5"]   = "4"
NIScope_Vertical_Dict["6"]   = "5"
NIScope_Vertical_Dict["7"]   = "6"
NIScope_Vertical_Dict["8"]   = "7"
# Channel enabled
NIScope_Vertical_Dict["on"]  = VI_TRUE
NIScope_Vertical_Dict["off"] = VI_FALSE
# Channel coupling
NIScope_Vertical_Dict["ac"]  = NISCOPE_VAL_AC
NIScope_Vertical_Dict["dc"]  = NISCOPE_VAL_DC
NIScope_Vertical_Dict["gnd"] = NISCOPE_VAL_GND
# Probe attenuation
NIScope_Vertical_Dict["1x"]   = 1
NIScope_Vertical_Dict["10x"]  = 10
NIScope_Vertical_Dict["100x"] = 100
# BW limit
NIScope_Vertical_Dict["full"]    = NISCOPE_VAL_BANDWIDTH_FULL
NIScope_Vertical_Dict["default"] = NISCOPE_VAL_BANDWIDTH_DEVICE_DEFAULT
NIScope_Vertical_Dict["20MHz"]   = NISCOPE_VAL_20MHZ_BANDWIDTH
NIScope_Vertical_Dict["100MHz"]  = NISCOPE_VAL_100MHZ_BANDWIDTH
# EnhanceRes - not available on niscope
NIScope_Vertical_Dict["0b"]   = 0
NIScope_Vertical_Dict["0.5b"] = 1
NIScope_Vertical_Dict["1b"]   = 2
NIScope_Vertical_Dict["1.5b"] = 3
NIScope_Vertical_Dict["2b"]   = 4
NIScope_Vertical_Dict["2.5b"] = 5
NIScope_Vertical_Dict["3b"]   = 6
# Input impedance
NIScope_Vertical_Dict["50"]   = NISCOPE_VAL_50_OHMS
NIScope_Vertical_Dict["1M"]   = NISCOPE_VAL_1_MEG_OHM

# Channel enabled
function set_state(obj::T;ch::String="1",st::String="on") where {T<:NIScope}
	#NISCOPE_ATTR_CHANNEL_ENABLED=1250005; # ((1000000 + 250000) + 5)
	ch = NIScope_Vertical_Dict[ch]
	st = NIScope_Vertical_Dict[st]
	status = niScope_SetAttributeViBoolean(obj, ch, NISCOPE_ATTR_CHANNEL_ENABLED, st)
end
# Channel input impedance
function set_input_impedance(obj::T;ch::String="1",val::String="1M") where {T<:NIScope}
	ch = NIScope_Vertical_Dict[ch]
	val = NIScope_Vertical_Dict[val]
	status = niScope_SetAttributeViReal64(obj, ch, NISCOPE_ATTR_INPUT_IMPEDANCE, Float64(val))
end
# Channel coupling
function set_coupling(obj::T;ch::String="1",cpl::String="dc") where {T<:NIScope}
	#NISCOPE_ATTR_VERTICAL_COUPLING = 1250003; # ((1000000 + 250000) + 3)
	ch  = NIScope_Vertical_Dict[ch]
	cpl = NIScope_Vertical_Dict[cpl]
	status = niScope_SetAttributeViInt32(obj, ch, NISCOPE_ATTR_VERTICAL_COUPLING, cpl)
end
# Probe attenuation
function set_atten(obj::T;ch::String="1",att::String="10x") where {T<:NIScope}
	#NISCOPE_ATTR_PROBE_ATTENUATION=1250004; # ((1000000 + 250000) + 4)
	ch  = NIScope_Vertical_Dict[ch]
	att = NIScope_Vertical_Dict[att]
	status = niScope_SetAttributeViReal64(obj, ch, NISCOPE_ATTR_PROBE_ATTENUATION, Float64(att))
end
# BW limit
function set_bw_lim(obj::T;ch::String="1",bw::String="full") where {T<:NIScope}
	#NISCOPE_ATTR_MAX_INPUT_FREQUENCY=1250006; # ((1000000 + 250000) + 6)
	ch = NIScope_Vertical_Dict[ch]
	bw = NIScope_Vertical_Dict[bw]
	status = niScope_SetAttributeViReal64(obj, ch, NISCOPE_ATTR_MAX_INPUT_FREQUENCY, Float64(bw))
end
# Nr of averages
function set_nr_avg(obj::T;ch::String="1",navg::Int=1) where {T<:NIScope}
	error("not implemented yet")
end
# EnhanceRes
function set_eres(obj::T;ch::String="1",bits::String="0b") where {T<:NIScope}
	error("not available for instruments of type ::NIScope")
end
# Degauss
function degauss(obj::T;ch::String="1",probe::String="CP030A") where {T<:NIScope}
	error("not available for instruments of type ::NIScope")
end
# Vertical range
function set_vrange(obj::T;ch::String="1",vrang=8) where {T<:NIScope}
	#NISCOPE_ATTR_VERTICAL_RANGE=1250001; # ((1000000 + 250000) + 1)
	ch = NIScope_Vertical_Dict[ch]
	status = niScope_SetAttributeViReal64(obj, ch, NISCOPE_ATTR_VERTICAL_RANGE, Float64(vrang))
end
# Vertical scale (V/div)
function set_vscale(obj::T;ch::String="1",volt=8) where {T<:NIScope}
	error("not implemented yet")
end
# Vertical offset
function set_voffs(obj::T;ch::String="1",voffs=0) where {T<:NIScope}
	#NISCOPE_ATTR_VERTICAL_OFFSET=1250002; # ((1000000 + 250000) + 2)
	ch = NIScope_Vertical_Dict[ch]
	status = niScope_SetAttributeViReal64(obj, ch, NISCOPE_ATTR_VERTICAL_OFFSET, Float64(voffs))
end


text:
module PkgTest
using Test
@testset "Generating Random Deviates" begin
    include("unit_test/generate_random_deviates/rand_continuous.jl")
    include("unit_test/generate_random_deviates/rand_discrete.jl")
    include("unit_test/generate_random_deviates/rand_multivariate.jl")
end

## VCM
@testset "VCM Covariance" begin
    include("unit_test/VCM/VCM_model_interface.jl")
    include("unit_test/VCM/singlerun_nbVCM.jl")
    include("unit_test/VCM/singlerun_bernoulliVCM.jl")
    include("unit_test/VCM/singlerun_normalVCM.jl")
    include("unit_test/VCM/singlerun_poissonVCM.jl")
    include("unit_test/VCM/singlerun_poisson_bernoulli_mixedVCM.jl")
end

## AR
@testset "AR Covariance" begin
    include("unit_test/AR/AR_model_interface.jl")
    include("unit_test/AR/singlerun_nbAR.jl")
    include("unit_test/AR/singlerun_bernoulliAR.jl")
    include("unit_test/AR/singlerun_normalAR.jl")
    include("unit_test/AR/singlerun_poissonAR.jl")
end

## CS
@testset "CS Covariance" begin
    include("unit_test/CS/CS_model_interface.jl")
    include("unit_test/CS/singlerun_nbCS.jl")
    include("unit_test/CS/singlerun_bernoulliCS.jl")
    include("unit_test/CS/singlerun_normalCS.jl")
    include("unit_test/CS/singlerun_poissonCS.jl")
end

end


text:
# T: coeff type
# Degree: degree
# Dim1: dimension of the support + 1
mutable struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T,NF} end

numfunctions(s::LagrangeRefSpace{T,D,2}) where {T,D} = D+1
numfunctions(s::LagrangeRefSpace{T,0,3}) where {T} = 1
numfunctions(s::LagrangeRefSpace{T,1,3}) where {T} = 3

valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
        SVector{numfunctions(ref), Tuple{T,T}}

# Evaluate constant lagrange elements on anything
(ϕ::LagrangeRefSpace{T,0})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),))

# Evaluate linear Lagrange elements on a segment
function (f::LagrangeRefSpace{T,1,2})(mp) where T
    u = mp.bary[1]
    j = jacobian(mp)
    SVector(
        (value=  u, derivative=-1/j),
        (value=1-u, derivative= 1/j))
end

# Evaluete linear lagrange elements on a triangle
function (f::LagrangeRefSpace{T,1,3})(t) where T
    u,v,w, = barycentric(t)
    SVector(
        (value=u,),
        (value=v,),
        (value=w,))
end

"""
    f(tangent_space, Val{:withcurl})

Compute the values of the shape functions together with their curl.
"""
function (f::LagrangeRefSpace{T,1,3})(t, ::Type{Val{:withcurl}}) where T
    # Evaluete linear Lagrange elements on a triange, together with their curl
    j = jacobian(t)
    u,v,w, = barycentric(t)
    p = t.patch
    SVector(
        (value=u, curl=(p[3]-p[2])/j),
        (value=v, curl=(p[1]-p[3])/j),
        (value=w, curl=(p[2]-p[1])/j)
    )
end

# Evaluate constant Lagrange elements on a triangle, with their curls
function (f::LagrangeRefSpace{T,0,3})(t, ::Type{Val{:withcurl}}) where T
    i = one(T)
    z = zero(cartesian(t))
    (
        (value=i, curl=z,),
    )
end

function curl(ref::LagrangeRefSpace, sh, el)
    sh1 = Shape(sh.cellid, mod1(sh.refid+1,3), -sh.coeff)
    sh2 = Shape(sh.cellid, mod1(sh.refid+2,3), +sh.coeff)
    return [sh1, sh2]
end

function strace(x::LagrangeRefSpace, cell, localid, face)

    Q = zeros(scalartype(x),2,3)

    p1 = neighborhood(face, 1)
    p2 = neighborhood(face, 0)

    u1 = carttobary(cell, cartesian(p1))
    u2 = carttobary(cell, cartesian(p2))

    P1 = neighborhood(cell, u1)
    P2 = neighborhood(cell, u2)

    vals1 = x(P1)
    vals2 = x(P2)

    for j in 1:numfunctions(x)
        Q[1,j] = vals1[j].value
        Q[2,j] = vals2[j].value
    end

    Q
end

function restrict(refs::LagrangeRefSpace{T,0}, dom1, dom2) where T
    #Q = eye(T, numfunctions(refs))
    Q = Matrix{T}(I, numfunctions(refs), numfunctions(refs))
end

function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T

    D = numfunctions(f)
    Q = zeros(T, D, D)

    # for each point of the new domain
    for i in 1:D
        v = dom2.vertices[i]

        # find the barycentric coordinates in dom1
        uvn = carttobary(dom1, v)

        # evaluate the shape functions in this point
        x = neighborhood(dom1, uvn)
        fx = f(x)

        for j in 1:D
            Q[j,i] = fx[j][1]
        end
    end

    return Q
end


text:
module EmpiricalCDFs
import Statistics
import Printf

"""
    module EmpiricalCDFs

Build and use empirical cumulative distribution functions.

Types/constructors: `AbstractEmpiricalCDF`, `EmpiricalCDF`, `EmpiricalCDFHi`

Functions/methods: Many of the methods defined for `AbstractEmpiricalCDF` are forwarded to the underlying
array.

`logprint`, `linprint`,
`length`, `counts`, `sort!`, `issorted`,
`empty!`, `mean`, `median`, `middle`,
`minimum`, `maximum`, `extrema`,
`std`, `stdm`, `var`, `varm`, `quantile`,
`iterate`

The statistical methods extend functions in `Statistics`. One way to use them with
`EmpiricalCDFs` is `using Statistics`.

Binary IO is supported by the submodule `EmpiricalCDFs.IOcdf`.
"""
EmpiricalCDFs

# cdfs.jl
export AbstractEmpiricalCDF, EmpiricalCDF, EmpiricalCDFHi,  linprint, logprint,
       getcdfindex, data, counts
export finv

include("cdfs.jl")

include("IOcdf.jl")
using .IOcdf

end # module


text:
@recipe function f(sim::T, mode=:ψ) where T <: Union{Calc, Sim}
    # Set Defaults
    res_x --> 500
    res_t --> 500
    n_lines --> 10
    skip --> 1

    # Extract Values
    res_x = plotattributes[:res_x]
    res_t = plotattributes[:res_t]
    skip = plotattributes[:skip]
    n_lines = plotattributes[:n_lines]

    # Compute sampling interval
    xₛ = Int(ceil(sim.box.Nₓ/res_x))
    tₛ = Int(ceil(sim.box.Nₜ/res_t))

    if mode==:ψ
        # Downsample
        x_ax = sim.box.t[1:tₛ:end]
        y_ax =  sim.box.x[1:xₛ:end]
        ψ = sim.ψ[1:tₛ:end, 1:xₛ:end]

        xguide --> L"t"
        yguide --> L"x"
        colorbar_title --> L"|\psi|"
        seriescolor --> :viridis
        fontfamily --> "Computer Modern"
        tick_direction --> :out
        s = get(plotattributes, :seriestype, :auto)
        if s == :surface
            zguide --> L"|\psi|"
            colorbar --> false
            camera --> (35, 75)
            grid --> false
        end
        @series begin
            # return series data
            x_ax, y_ax, abs.(ψ)'
        end
    elseif mode==:ψ̃
        s = get(plotattributes, :seriestype, :auto)
        if s == :auto
            # Downsample
            y_ax =  sim.box.x[1:xₛ:end]
            #ψ̃ = sim.ψ̃[1:tₛ:end, 1:xₛ:end]
            # Extract Lines
            start = sim.box.Nₜ÷2 + 1
            ψ̃ = sim.ψ̃[start:skip:start+skip*n_lines-1, 1:xₛ:end]
            labs = reshape([L"\tilde{\psi}_{%$(i-1)}" for i in 1:skip:skip*n_lines], 1, :)
            tick_direction --> :out
            linewidth --> 1.5
            xguide --> L"x"
            yguide --> L"\log|\tilde{\psi}|"
            fontfamily --> "Computer Modern"
            tick_direction --> :out
            #margin --> 4mm
            legend --> :outertopright
            for i=1:size(ψ̃)[1]
                @series begin
                    label --> labs[i]
                    y_ax,  log.(abs.(ψ̃[i, :]))
                end
            end
        else
            # Downsample
            x_ax = sim.box.ω[1:tₛ:end]
            y_ax =  sim.box.x[1:xₛ:end]
            ψ̃ = sim.ψ̃[1:tₛ:end, 1:xₛ:end]

            xguide --> L"\omega"
            yguide --> L"x"
            colorbar_title --> L"\log|\tilde{\psi}|"
            seriescolor --> :viridis
            fontfamily --> "Computer Modern"
            tick_direction --> :out
            @series begin
                seriestype := s
                # return series data
                x_ax, y_ax, log.(abs.(ψ̃'))
            end
        end
    elseif mode==:IoM
        # Downscale
        x = sim.box.x[1:xₛ:end]
        E = sim.E[1:xₛ:end]
        KE = sim.KE[1:xₛ:end]
        PE = sim.PE[1:xₛ:end]
        δE = sim.E[1:xₛ:end] .- sim.E[1]
        δP = sim.P[1:xₛ:end] .- sim.P[1]
        δN = sim.N[1:xₛ:end] .- sim.N[1]

        layout --> (2, 2)
        fontfamily --> "Computer Modern"
        tick_direction --> :out
        @series begin
            linewidth --> 1.5
            subplot := 1
            label --> L"E"
            x, E
        end
        @series begin
            linewidth --> 1.5
            subplot := 1
            label --> L"K"
            x, KE
        end
        @series begin
            linewidth --> 1.5
            subplot := 1
            xguide := L"x"
            yguide := L"E"
            label --> L"V"
            x, PE
        end
        @series begin
            linewidth --> 1.5
            subplot := 2
            xguide := L"x"
            yguide := L"\delta E"
            label --> ""
            x, δE
        end
        @series begin
            linewidth --> 1.5
            subplot := 3
            xguide := L"x"
            yguide := L"\delta P"
            label --> ""
            x, δP
        end
        @series begin
            linewidth --> 1.5
            subplot := 4
            xguide := L"x"
            yguide := L"\delta N"
            label --> ""
            x, δN
        end
    else
        @error "Unknown Plotting Mode for NonlinearSchrodinger.jl Objects $mode"
    end

    return nothing
end


text:
# whileの中身はiを参照できないからglobalをつける
i = 0
while i < 10
  global i += 1
  println(i)
end

# iterableならなんでもいい
fruits = ["apple","banana","orange"]

i = 1
while i<= length(fruits)
    println(fruits[i])
    global i+=1
end

# for

for i = 1:10
  println(i)
end

for i in fruits
    println(i)
end

fruits = ("apple","banana","orange")
for i in fruits
    println(i)
end

fruits = Dict("apple"=>1,"banana"=>2,"orange"=>3)
# かっこをつけないとエラー
for (k,v) in fruits
    println(k,v)
end

apple = "apple"
for i in apple
    println(i)
end

fruits = Set(["apple","banana","orange"])
for i in fruits
    println(i)
end

## 少し書き方が特殊なループ

## 二重forループがあっても一行でかけてしまう．
for i = 1:2, j=2:3
    println((i,j))
end

for i =1:5, j=2:20
    println((i,j))
end

# 同時にループを回したいときはzip関数をつかう
for ( i,j ) in zip(1:3,2:4)
    println(i,j)
end

# 内包表記
# 内包表記を使って，行列がすぐ作れる

mat = [i+j for i=1:10,j=1:10 ]
println(mat)
println(size(mat))


text:
function defines_CM_field(F::FieldsTower)
  K = F.field
  if isodd(degree(K)) || !istotally_complex(K)
    return false, id_hom(K)
  end
  autsK = F.generators_of_automorphisms
  permGC = GAP.julia_to_gap([_perm_to_gap_perm(x) for x in permutations(autsK)])
  perm_group = GAP.Globals.GroupByGenerators(permGC)
  Z = GAP.Globals.Center(perm_group)
  if isodd(GAP.Globals.Size(Z))
    return false, id_hom(K)
  end
  els = GAP.Globals.Elements(Z)
  for i = 1:length(els)
    @show i
    el = els[i]
    if GAP.Globals.Order(els[i]) != 2
      continue
    end
    word = GAP.Globals.Factorization(perm_group, el)
    rep_of_word = GAP.Globals.ExtRepOfObj(word)
    aut = autsK[rep_of_word[1]]^rep_of_word[2]
    for j = 2:div(length(rep_of_word), 2)
      aut *= autsK[rep_of_word[2*j-1]]^rep_of_word[2*j]
    end
    if iscomplex_conjugation(aut)
      return true, aut
    end
  end
  return false, id_hom(K)
end


text:
# Show AbstractAxis types
Base.show(io::IO, ::MIME"text/plain", ::Axis{IdxMap}) where IdxMap = print(io, "Axis$IdxMap")
Base.show(io::IO, ::Axis{IdxMap}) where IdxMap = print(io, "Axis$IdxMap")

Base.show(io::IO, ::FlatAxis) = print(io, "FlatAxis()")

Base.show(io::IO, ::NullAxis) = print(io, "NullAxis()")

Base.show(io::IO, ::MIME"text/plain", ::PartitionedAxis{PartSz, IdxMap, Ax}) where {PartSz, IdxMap, Ax} =
    print(io, "PartitionedAxis($PartSz, $(Ax()))")
Base.show(io::IO, ::PartitionedAxis{PartSz, IdxMap, Ax}) where {PartSz, IdxMap, Ax} =
    print(io, "PartitionedAxis($PartSz, $(Ax()))")

Base.show(io::IO, ::ShapedAxis{Shape, IdxMap}) where {Shape, IdxMap} =
    print(io, "ShapedAxis($Shape, $IdxMap)")

Base.show(io::IO, ::MIME"text/plain", ::ViewAxis{Inds, IdxMap, Ax}) where {Inds, IdxMap, Ax} =
    print(io, "ViewAxis($Inds, $(Ax()))")
Base.show(io::IO, ::ViewAxis{Inds, IdxMap, <:Ax}) where {Inds, IdxMap, Ax} =
    print(io, "ViewAxis($Inds, $(Ax()))")
Base.show(io::IO, ::ViewAxis{Inds, IdxMap, <:NullorFlatAxis}) where {Inds, IdxMap} =
    print(io, Inds)

Base.show(io::IO, ci::ComponentIndex) = print(io, "ComponentIndex($(ci.idx), $(ci.ax))")

# Show ComponentArrays
Base.show(io::IO, ::MIME"text/plain", ::Type{ComponentArray{T,N,A,Ax}}) where {T,N,A,Ax} = print(io, "ComponentArray{$T,$N,$A,$Ax}") # make `typeof(u)` show the full type
Base.show(io::IO, ::Type{<:ComponentArray{T,N,<:Array}}) where {T,N} = print(io, "ComponentArray{$T,$N}") # do not pollute the stacktrace with verbose type printing
Base.show(io::IO, ::Type{<:ComponentArray{T,1,<:Array}}) where T = print(io, "ComponentVector{$T}")
Base.show(io::IO, ::Type{<:ComponentArray{T,2,<:Array}}) where T = print(io, "ComponentMatrix{$T}")
Base.show(io::IO, ::Type{<:ComponentArray{T,N,A}}) where {T,N,A} = print(io, "ComponentArray{$T,$N,$(A.name)...}")
Base.show(io::IO, ::Type{<:ComponentArray{T,1,A}}) where {T,A} = print(io, "ComponentVector{$T,$(A.name)...}")
Base.show(io::IO, ::Type{<:ComponentArray{T,2,A}}) where {T,A} = print(io, "ComponentMatrix{$T,$(A.name)...}")

function Base.show(io::IO, x::ComponentVector)
    print(io, "(")
    for (i,key) in enumerate(keys(x))
        if i==1
            print(io, "$key = ")
        else
            print(io, ", $key = ")
        end
        show(io, x[key])
    end
    print(io, ")")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", x::ComponentVector)
    show(io, typeof(x))
    show(io, x)
    return nothing
end

function Base.show(io::IO, a::AbstractVector{<:T}) where T<:ComponentVector
    elem = a[1]
    print(io, "[$elem")
    for idx in 2:length(a)
        print(io, ", $(a[idx])")
    end
    print(io, "]")
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", x::ComponentMatrix{T,A,Axes}) where {T,A,Axes}
    if !haskey(io, :compact) && length(axes(x, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    axs = getaxes(x)
    sz = size(x)
    println(io, "$(sz[1])×$(sz[2]) $(typeof(x)) with axes $(axs[1]) × $(axs[2])")
    Base.print_matrix(io, getdata(x))
    return nothing
end


text:
export mat2CDS

"""
Convert sparse and square matrix A to compressed diagonal format (CDS)
Currently not multithreaded
"""
function mat2CDS(A::SparseMatrixCSC{TF,TI}) where {TF<:Real,TI<:Integer}

  # Find all nonzero diagonals
  (i,j,dummy) = findnz(A) #[i,j] = find(A);
  d           = sort(j-i)
  #julia 0.6: offset      = d[find(diff([-Inf;d]))] #matlab:d = d(find(diff([-inf; d])));
  offset      = d[findall(!iszero, diff([-Inf;d]))]
  offset      = convert(Array{TI,1},offset)

  # initialize
  (m,n) = size(A)
  p     = length(offset)
  R     = zeros(TF,m,p)

  #construct tall dense matrix with diagonals of A
  for i=1:length(offset)
    dA = diag(A,offset[i])
    if offset[i]>=0
      R[1:length(dA),i]=dA;
    elseif offset[i]<0
      R[end-length(dA)+1:end,i]=dA;
    end
  end

  return R, offset
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
struct Planar
end

function calcArea(pg::Planar, r::Float64)::Float64
    return 1
end

function calcVolume(pg::Planar, r::Float64)::Float64
    return r
end

struct Cylindrical
end

function calcArea(pg::Cylindrical, r::Float64)::Float64
    return 2*pi*r
end

function calcVolume(pg::Cylindrical, r::Float64)::Float64
    return pi*r^2
end

struct Spherical
end

function calcArea(pg::Spherical, r::Float64)::Float64
    return 4*pi*r^2
end

function calcVolume(pg::Spherical, r::Float64)::Float64
    return 4*pi*r^3/3
end

struct SphericalComplamentary
end

function calcSource(source::SphericalComplamentary,
                    state::HydroSnapshot,
                    t::Float64,
                    dt::Float64)::Array{Conserved,1}
    cumulative_volumes = [4*pi*r^3/3 for r in state.grid]
    volumes = diff(cumulative_volumes)
    r_list = midGrid(state.grid)
    return [Conserved(0,2*volume*c.pressure/r,0)
            for (c,r,volume) in zip(state.cells,r_list,volumes)]
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
module TestAstTools

using FLoops: unbound_rhs
using Test

function test_unbound_rhs()
    @test unbound_rhs(:a) == [:a]
    @test unbound_rhs(:(let x = 1; (x, y); end)) == [:y]
    @test unbound_rhs(quote end) == []
    @test unbound_rhs(quote () end) == []
    @test unbound_rhs(quote [] end) == []
    @test unbound_rhs(nothing) == []
end

end  # module


text:
export @memoize

"""
    @memoize [arg::T]... begin
        # expensive computation
    end::T

Low-level, no-frills memoization macro that stores values in a thread-local, typed Dict. The
types of the dictionary are derived from the syntactical type assertions.

When there are no arguments to key the cache with, instead of a dictionary a simple array
with per-thread elements is used. This further improves performance to 2ns per access.
"""
macro memoize(ex...)
    code = ex[end]
    args = ex[1:end-1]

    # decode the code body
    @assert Meta.isexpr(code, :(::))
    rettyp = code.args[2]
    code = code.args[1]

    # decode the arguments
    argtyps = []
    argvars = []
    for arg in args
        @assert Meta.isexpr(arg, :(::))
        push!(argvars, arg.args[1])
        push!(argtyps, arg.args[2])
    end

    # the global cache is an array with one entry per thread. if we don't have to key on
    # anything, that entry will be the memoized new_value, or else a dictionary of values.
    @gensym global_cache

    # generate code to access memoized values
    # (assuming the global_cache can be indexed with the thread ID)
    if isempty(args)
        # if we don't have to key on anything, use the global cache directly
        global_cache_eltyp = :(Union{Nothing,$rettyp})
        ex = quote
            cache = get!($(esc(global_cache))) do
                [nothing for _ in 1:Threads.nthreads()]
            end
            cached_value = @inbounds cache[Threads.threadid()]
            if cached_value !== nothing
                cached_value
            else
                new_value = $(esc(code))::$rettyp
                @inbounds cache[Threads.threadid()] = new_value
                new_value
            end
        end
    else
        if length(args) == 1
            global_cache_eltyp = :(Dict{$(argtyps[1]),$rettyp})
            global_init = :(Dict{$(argtyps[1]),$rettyp}())
            key = :($(esc(argvars[1])))
        else
            global_cache_eltyp = :(Dict{Tuple{$(argtyps...)},$rettyp})
            global_init = :(Dict{Tuple{$(argtyps...)},$rettyp}())
            key = :(tuple($(map(esc, argvars)...)))
        end
        ex = quote
            cache = get!($(esc(global_cache))) do
                [$global_init for _ in 1:Threads.nthreads()]
            end
            local_cache = @inbounds cache[Threads.threadid()]
            cached_value = get(local_cache, $key, nothing)
            if cached_value !== nothing
                cached_value
            else
                new_value = $(esc(code))::$rettyp
                local_cache[$key] = new_value
                new_value
            end
        end
    end

    # define the per-thread cache
    @eval __module__ begin
        const $global_cache = LazyInitialized{Vector{$(global_cache_eltyp)}}()
    end

    quote
        $ex
    end
end


text:
"""
    $(TYPEDEF)

`DiagGaussianPolicy` policy represents a stochastic control policy, represented as a
multivariate Gaussian distribution of the form:

```math
\\pi_{\\theta}(a | o) = \\mathcal{N}(\\mu_{\\theta_1}(o), \\Sigma_{\\theta_2})
```

where ``\\mu_{\\theta_1}`` is a neural network, parameterized by ``\\theta_1``, that maps
an observation to a mean action and ``\\Sigma_{\\theta_2}`` is a diagonal
covariance matrix parameterized by ``\\theta_2``, the diagonal entries of the matrix.
Rather than tracking ``\\Sigma_{\\theta_2}`` directly, we track the log standard deviations,
which are easier to learn. Note that ``\\mu_{\\theta_1}`` is a _state-dependent_ mean
while ``\\Sigma_{\\theta_2}`` is a _global_ covariance.
"""
struct DiagGaussianPolicy{Mean,Logstd<:AbstractVector}
    meanNN::Mean
    logstd::Logstd
    fixedlogstd::Bool
end

"""
    $(SIGNATURES)

Construct a `DiagGaussianPolicy` with a state-dependent mean `meanNN` and initial
log-standard deviation `logstd`. If `fixedlogstd` is true, `logstd` will be treated as
a constant. `meanNN` should be object that is compatible with Flux.jl and have the
following signatures:
- `meanNN(obs::AbstractVector)` --> `action::AbstractVector`
- `meanNN(obs::AbstractMatrix)` --> `action::AbstractMatrix`
"""
function DiagGaussianPolicy(meanNN, logstd::AbstractVector; fixedlogstd::Bool = false)
    DiagGaussianPolicy(meanNN, logstd, fixedlogstd)
end

Flux.@functor DiagGaussianPolicy

function Flux.trainable(policy::DiagGaussianPolicy)
    if policy.fixedlogstd
        (meanNN = policy.meanNN, ) \propagate_inbounds
    else
        (meanNN = policy.meanNN, logstd = policy.logstd)
    end
end

"""
    sample!([rng = GLOBAL_RNG, ]action, policy, feature)

Treating `policy` as a stochastic policy, sample an action from `policy`, conditioned on
`feature`, and store it in `action`.
"""
function sample!(rng::AbstractRNG, action::AbsVec, policy::DiagGaussianPolicy, feature::AbsVec)
    randn!(rng, action)
    action .= action .* exp.(policy.logstd) .+ policy(feature)
end
function sample!(action::AbsVec, policy::DiagGaussianPolicy, feature::AbsVec)
    sample!(default_rng(), action, policy, feature)
end

(policy::DiagGaussianPolicy)(features) = policy.meanNN(features)

"""
    $(SIGNATURES)

Treating `policy` as a deterministic policy, compute the mean action of `policy`,
conditioned on `feature`, and store it in `action`.
"""
function LyceumBase.getaction!(action::AbsVec, policy::DiagGaussianPolicy, feature::AbsVec)
    action .= action .* exp.(policy.logstd) .+ policy(feature)
end

"""
    $(SIGNATURES)

Return loglikelihood of `action` conditioned on `feature` for `policy`.
"""
function loglikelihood(policy::DiagGaussianPolicy, action::AbsVec, feature::AbsVec)
    meanact = policy(feature)
    ls = policy.logstd
    ll = -length(ls) * log(2pi) / 2
    for i = 1:length(action)
        ll -= ((meanact[i] - action[i]) / exp(ls[i]))^2 / 2
        ll -= ls[i]
    end
    ll
end

"""
    $(SIGNATURES)

Treating each column of `actions` and `features` as a single action/feature, return
a vector of the loglikelihoods of `actions` conditioned on `features` for `policy`.
"""
function loglikelihood(policy::DiagGaussianPolicy, actions::AbsMat, features::AbsMat)
    constterm = -length(policy.logstd) * log(2pi) / 2 - sum(policy.logstd)
    meanactions = policy(features)

    #zs2 = ((meanactions .- actions) ./ exp.(m.logstd)) .^ 2
    # NOTE: this is about 2x faster than the above line
    zs = ((meanactions .- actions) ./ exp.(policy.logstd))
    zs = zs .^ 2

    #ll = DT(-0.5) * sum(zs, dims=1) .+ constterm
    ll = constterm .- sum(zs, dims = 1) / 2
    dropdims(ll, dims = 1)
end

"""
    $(SIGNATURES)

Return the gradient of the loglikelihood of `action` conditioned on `feature` with respect
to `policy`'s parameters.
"""
function grad_loglikelihood!(
    gradll::AbsVec,
    policy::DiagGaussianPolicy,
    action::AbsVec,
    feature::AbsVec
)
    ps = params(policy)
    gs = gradient(() -> loglikelihood(policy, action, feature), ps)
    copygrad!(gradll, ps, gs)
    gradll
end

"""
    $(SIGNATURES)

Treating each column of `actions` and `features` as a single action/feature,
return a vector of gradients of the loglikelihood of `action` conditioned on `feature` with respect
to `policy`'s parameters.

This computation is done in parallel using `nthreads` threads.
"""
@propagate_inbounds function grad_loglikelihood!(
    gradlls::AbsMat,
    policy::DiagGaussianPolicy,
    actions::AbsMat,
    features::AbsMat,
    nthreads::Int = Threads.nthreads(),
)
    @boundscheck begin
        if !(axes(gradlls, 2) == axes(actions, 2) == axes(features, 2))
            throw(ArgumentError("gradlls, actions, and features must have same 2nd axis"))
        end
        require_one_based_indexing(gradlls, actions, features)
    end

    bthreads = max(1, div(Threads.nthreads(), nthreads))
    ranges = splitrange(size(gradlls, 2), nthreads)
    @with_blasthreads bthreads begin
        @sync for r in ranges
            Threads.@spawn _threaded_gradll!(gradlls, policy, actions, features, r)
        end
    end
    gradlls
end

function _threaded_gradll!(
    gradlls::AbsMat,
    m::DiagGaussianPolicy,
    acts::AbsMat,
    feats::AbsMat,
    range::UnitRange{Int},
)
    for i in range
        feat = uview(feats, :, i)
        act = uview(acts, :, i)
        gradll = uview(gradlls, :, i)
        grad_loglikelihood!(gradll, m, act, feat)
    end
    gradlls
end


text:
# FWI on the 2D Overthrust model using spectral projected gradient descent
# Author: mlouboutin3@gatech.edu
# Date: February 2021
#
using DrWatson
@quickactivate :TimeProbeSeismic
import TimeProbeSeismic: qr_data

# Setup a 2 layers model
n = (101, 101)
d = (10., 10.)
o = (0., 0.)

m = 1.5f0^(-2) * ones(Float32, n);
m[:, 51:end] .= .25f0;
# collect(m[i, i÷2+30:end] .= .35f0 for i=1:101);

model = Model(n, d, o, m; nb=80)
model0 = smooth(model; sigma=5)
dm = model0.m - model.m

# Simple geometry
# Src/rec sampling and recording time
timeD = 1000f0   # receiver recording time [ms]
dtD = get_dt(model0)    # receiver sampling interval [ms]
nsrc = 1

nxrec = n[1]
xrec = range(0f0, stop=(n[1] - 1)*d[1], length=nxrec)
yrec = 0f0
zrec = range(2*d[2], stop=2*d[2], length=nxrec)

# xsrc = convertToCell(range(0f0, stop=(n[1] - 1)*d[1], length=nsrc))
# ysrc =  convertToCell(range(0f0, stop=0f0, length=nsrc))
# zsrc =  convertToCell(range(2*d[2], stop=2*d[2], length=nsrc))

xsrc = 500f0
ysrc = 0f0
zsrc = 2*d[2]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtD, t=timeD, nsrc=nsrc)
# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtD, t=timeD)

# setup wavelet
f0 = 0.01f0     # kHz
wavelet = ricker_wavelet(timeD, dtD, f0)
q = judiVector(srcGeometry, wavelet)

# Forward operator
opt = Options(space_order=16, sum_padding=true, isic=true)
F = judiModeling(model, srcGeometry, recGeometry; options=opt)
F0 = judiModeling(model0, srcGeometry, recGeometry; options=opt)
J = judiJacobian(F0, q)

# data
# Nonlinear modeling
dobs = F*q
d0 = F0*q
residual = d0 - dobs
δd = J*dm

# gradient
g = J'*residual

# Probe
ge = Array{Any}(undef, 8)
for ps=1:8
    ge[ps] = judiJacobian(F0, q, 2^ps, dobs)'*residual
end

clip = maximum(g)/10

figure()
subplot(331)
imshow(g', cmap="seismic", vmin=-clip, vmax=clip, aspect=.5)
title("Reference")
for ps=1:8
    subplot(3,3,ps+1)
    imshow(ge[ps]', cmap="seismic", vmin=-clip, vmax=clip, aspect=.5)
    title("ps=$(2^ps)")
end
tight_layout()

clip = maximum(g)/10

figure()
subplot(331)
imshow(g', cmap="seismic", vmin=-clip, vmax=clip, aspect=.5)
title("Reference")
for ps=1:8
    subplot(3,3,ps+1)
    imshow(g' - ge[ps]', cmap="seismic", vmin=-clip, vmax=clip, aspect=.5)
    title("Difference ps=$(2^ps)")
end
tight_layout()

# Similarities
similar = [simil(g, ge[i]) for i=1:8]
# Similarities with muted wated
mw(x::PhysicalParameter) = mw(x.data)
mw(x::Array) = x[:, 30:end]
similar2 = [simil(mw(g), mw(ge[i])) for i=1:8]

figure()
semilogx([2^i for i=1:8], similar, "-or", label="Full gradient", basex=2)
semilogx([2^i for i=1:8], similar2, "-ob", label="Muted water layer", basex=2)
title(L"Similarity $\frac{<g, g_e>}{||g|| ||g_e||}$")
xlabel("Number of probing vectors")
legend()

# Adjoint test
at = Array{Any}(undef, 8)
J0 = dot(δd, residual)
for ps=1:8
    at[ps] = dot(ge[ps], dm)
    @printf("ps=%d,  <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n",
            2^ps, J0, at[ps], (J0 - at[ps])/(J0 + at[ps]))
end

att = [abs(1 - J0/at[i]) for i=1:8]
figure()
semilogx([2^i for i=1:8], att, basex=2, label="Probed adjoint test")
semilogx([2^i for i=1:8], [0 for i=1:8], label="target")
title(L"$1 - \frac{<J \delta m, g_e>}{<J' \delta d, \delta m>}$")

# Probing vectors
figure();imshow(dobs.data[1]*dobs.data[1]', cmap="seismic", vmin=-10, vmax=10)
title(L"$d_{obs} d_{obs}^\top$")

figure()
for ps=1:9
    subplot(3,3,ps)
    Q = qr_data(dobs.data[1]*dobs.data[1]', 2^ps)
    imshow(Q, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
    title("Probing vectors ps=$(2^ps)")
end
tight_layout()

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 359a6c00-c9a6-441d-b258-55bfb5deb4b5
begin
    import Pkg
    Pkg.develop("PoreMatMod")
end

# ╔═╡ a948f8b3-4ec5-40b9-b2c1-fcf5b8ad67fa
# load required packages (Pluto.jl will automatically install them)
using PoreMatMod, PlutoUI

# ╔═╡ 77c63e60-2708-4fef-a6ea-cb394f114b88
using PoreMatMod.ExampleHelpers

# ╔═╡ 8d523993-6e85-443a-9949-12030552b457
md"""
## example: generating hypothetical MOF structures with different replacement styles
"""

# ╔═╡ 74b45651-21d5-4332-a4a5-866ea1bb02b8
input_file_message()

# ╔═╡ 3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
xtal_folder_message()

# ╔═╡ ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
rc[:paths][:crystals]

# ╔═╡ b66dd0d5-5bac-4b23-a33f-17305b0d4728
moiety_folder_message()

# ╔═╡ a76e79ee-c82c-4771-818d-5380a5fb4c18
rc[:paths][:moieties]

# ╔═╡ be86f07e-0669-46c2-8c79-80e65dfcc6f2
md"""
!!! example \"the task\"
    we have the IRMOF-1 crystal structure, and wish to explore appending the acetamido functional group its BDC linkers using different replacement styles.

**Parent crystal structure**: first, we read in the `.cif` file describing the parent IRMOF-1 crystal structure.
"""

# ╔═╡ b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
begin
    # read in the parent xtal
    parent = Crystal("IRMOF-1.cif") # load .cif file
    infer_bonds!(parent, true)      # infer bonds
    view_structure(parent)          # view structure
end

# ╔═╡ b402e79a-784b-4d8b-82f1-df4fe1cedad1
fragment_construction_note()

# ╔═╡ 9b2d534a-4f78-4950-add1-9ba645669bb9
md"""
**Query fragment**: next, we construct a query fragment as a _p_-phenylene fragment to match that on the BCD linker of the IRMOF-1 parent structure. we mark one hydrogen atom on the query fragment as "masked" (shown in pink) by tagging its species label with `!` in the input. we need to mask this hydrogen atom because it will eventually be replaced by the acetamido functional group.
"""

# ╔═╡ 4be03110-61ab-4cd6-b3fe-7d51ac5ee771
query = moiety("2-!-p-phenylene.xyz");

# ╔═╡ e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
view_query_or_replacement("2-!-p-phenylene.xyz")

# ╔═╡ 559ef3c3-a176-4d65-8958-810c9b0b32c5
with_terminal() do
    display_query_or_replacement_file("2-!-p-phenylene.xyz")
end

# ╔═╡ d1aa8a19-3702-40f6-b73e-b9ebfc1a7a71
md"""
**Replacement fragment**: next, we construct a replacement fragment as a modified version of the query fragment (with acetamido group in place of one hydrogen atom).
"""

# ╔═╡ 44e7e665-ae68-4cd4-b45f-138a0fb8910e
replacement = moiety("2-nitro-p-phenylene.xyz");

# ╔═╡ af606f02-fbd5-4033-a5f3-56a8f740b1a1
view_query_or_replacement("2-nitro-p-phenylene.xyz")

# ╔═╡ c8f8dbd3-4191-460d-944f-f5e456ce8b83
with_terminal() do
    display_query_or_replacement_file("2-nitro-p-phenylene.xyz")
end

# ╔═╡ 127a2bef-3903-4801-bc75-00a6dde2bc6e
md"""
### Replacement Modes

There are several options when performing a replacement that control how the replacement fragments will be positioned in the parent structure.

With all three file inputs loaded (IRMOF-1 as `parent`, 2-!-*p*-phenylene as `query`, and 2-nitro-*p*-phenylene as `replacement`) and a `search` performed, replacements can be made.

`PoreMatMod.jl` has several replacement modes, one of which must be specified.
"""

# ╔═╡ 0e3f7334-9e7f-483b-a0ac-70777902bf51
md"""
!!! note \"Note\"
    Multiple replacements can be done with a single search.
"""

# ╔═╡ b7b46022-aee0-4d51-8a5e-0c8c005f341a
search = query ∈ parent

# ╔═╡ 68557426-3204-45a1-8801-84c459e5ac66
md"""
#### Default: optimal orientation at all locations

Optimal configurations will be chosen for each location in `search.isomorphisms`, so that each occurrence of the `query` in the `parent` is replaced with minimal perturbation of the conserved atoms from the parent structure.
"""

# ╔═╡ 74aa19d2-b1a4-4333-9ff9-e6ea74e7d989
begin
    local child = substructure_replace(search, replacement)
    view_structure(child)
end

# ╔═╡ bed9d62f-d61d-4fb1-9a12-188a864fff21
md"""
#### Optimal Orientation at *n* random locations

If only some of the linkers should be replaced, the `nb_loc` argument lets us specify how many.
"""

# ╔═╡ 9dbf7859-762d-4edc-8300-ac3bba151a8a
begin
    local child = substructure_replace(search, replacement, nb_loc=8)
    view_structure(child)
end

# ╔═╡ 95215fe9-a7e6-4563-94b6-5ed1d5dde03b
md"""
#### Optimal orientation at specific locations

Specific locations are chosen by providing the `loc` argument.
"""

# ╔═╡ 8c22f94c-b235-40b2-8fc8-6052c31a6b6e
begin
    local child = substructure_replace(search, replacement, loc=[17, 18, 19, 20])
    view_structure(child)
end

# ╔═╡ 6297833a-ea26-4958-ad56-fc76aeabee69
md"""
#### Specific replacements

Providing both the `loc` and `ori` arguments allows specifying the exact configuration used in each replacement.  A zero value for any element of `ori` means to use the optimal replacement at the corresponding location.
"""

# ╔═╡ 183f9ff1-810a-4062-a3c6-cdb82dbd8a7a
begin
    local child = substructure_replace(
        search, replacement, loc=[1, 2, 3, 13], ori=[0, 1, 2, 3])
    view_structure(child)
end

# ╔═╡ 0347c45c-899b-4c4d-9b31-f09a205636f0
md"""
#### Random orientations

By using the `random` keyword argument, the search for optimal alignment can be skipped, and an arbitrary alignment option will be used.
"""

# ╔═╡ b7064b05-6e57-4e81-8cbc-b2075166a1af
begin
    local child = substructure_replace(search, replacement, nb_loc=24, random=true)
    view_structure(child)
end

# ╔═╡ Cell order:
# ╟─359a6c00-c9a6-441d-b258-55bfb5deb4b5
# ╟─8d523993-6e85-443a-9949-12030552b457
# ╠═a948f8b3-4ec5-40b9-b2c1-fcf5b8ad67fa
# ╠═77c63e60-2708-4fef-a6ea-cb394f114b88
# ╟─74b45651-21d5-4332-a4a5-866ea1bb02b8
# ╟─3ff0cd63-8bdf-4ba5-92b6-e2a52f7573c2
# ╠═ccc7fc86-5f48-4cc9-bd5c-349fc1d55b0c
# ╟─b66dd0d5-5bac-4b23-a33f-17305b0d4728
# ╠═a76e79ee-c82c-4771-818d-5380a5fb4c18
# ╟─be86f07e-0669-46c2-8c79-80e65dfcc6f2
# ╠═b53e7c38-d8f5-4f28-a9dc-f7902a86fdb2
# ╟─b402e79a-784b-4d8b-82f1-df4fe1cedad1
# ╟─9b2d534a-4f78-4950-add1-9ba645669bb9
# ╠═4be03110-61ab-4cd6-b3fe-7d51ac5ee771
# ╟─e5eaaef4-13a0-48bd-9f2b-5040b2d10ac1
# ╟─559ef3c3-a176-4d65-8958-810c9b0b32c5
# ╟─d1aa8a19-3702-40f6-b73e-b9ebfc1a7a71
# ╠═44e7e665-ae68-4cd4-b45f-138a0fb8910e
# ╟─af606f02-fbd5-4033-a5f3-56a8f740b1a1
# ╟─c8f8dbd3-4191-460d-944f-f5e456ce8b83
# ╟─127a2bef-3903-4801-bc75-00a6dde2bc6e
# ╟─0e3f7334-9e7f-483b-a0ac-70777902bf51
# ╠═b7b46022-aee0-4d51-8a5e-0c8c005f341a
# ╟─68557426-3204-45a1-8801-84c459e5ac66
# ╠═74aa19d2-b1a4-4333-9ff9-e6ea74e7d989
# ╟─bed9d62f-d61d-4fb1-9a12-188a864fff21
# ╠═9dbf7859-762d-4edc-8300-ac3bba151a8a
# ╟─95215fe9-a7e6-4563-94b6-5ed1d5dde03b
# ╠═8c22f94c-b235-40b2-8fc8-6052c31a6b6e
# ╟─6297833a-ea26-4958-ad56-fc76aeabee69
# ╠═183f9ff1-810a-4062-a3c6-cdb82dbd8a7a
# ╟─0347c45c-899b-4c4d-9b31-f09a205636f0
# ╠═b7064b05-6e57-4e81-8cbc-b2075166a1af


text:
export transform_iter

"""
    transform_iter(f, it) = Base.Iterators.accumulate((a,b)->f(b), it)

Transform iterator that applies `f` to each element of `it`.
"""
transform_iter(f, it) = (f(v) for v in it)
# transform_iter(f, it) = Base.Iterators.accumulate((a,b)->f(b), it)