text:
function zbl(xij, qi, qj, ti, tj, mu, eta, kappa)

    ng = length(ti)
    u = zeros(ng)
    u_xij = zeros(3,ng)
    for i = 1:ng
        xij1 = xij[1,i];
        xij2 = xij[2,i];
        xij3 = xij[3,i];
        t2 = xij1*xij1;
        t3 = xij2*xij2;
        t4 = xij3*xij3;
        t5 = t2+t3+t4;
        t6 = sqrt(t5);
        t7 = t6*1.0E2;
        t8 = t7-4.0E2;
        t9 = tanh(t8);
        t10 = t6-4.0;
        t11 = t10*t10;
        t12 = 1.0/sqrt(t5);
        t13 = t6*3.428291787990583E-2;
        t14 = t13-1.827782636402304E-1;
        t15 = t9*(1.0/2.0);
        t16 = t15+1.0/2.0;
        t21 = t6*2.308969885292557;
        t17 = exp(-t21);
        t23 = t6*4.614046060829141;
        t18 = exp(-t23);
        t25 = t6*3.664438963872198E1;
        t19 = exp(-t25);
        t27 = t6*1.079118754693148E1;
        t20 = exp(-t27);
        t22 = t17*2.16164490013485E3;
        t24 = t18*2.15028801532051E4;
        t26 = t19*1.394671496625875E4;
        t28 = t20*3.91244681854013E4;
        t29 = t22+t24+t26+t28;
        t30 = t9*t9;
        t31 = t30-1.0;
        t32 = t10*t11*t14;
        t33 = t32+1.625187736600612E-2;
        t34 = 1.0/t5^(3.0/2.0);
        u[i] = t9*8.125938683003058E-3+t12*t29-t16*t33-8.125938683003058E-3;
        u_xij[1,i] = -t12*(t12*t17*xij1*4.991172977107606E3+t12*t18*xij1*9.921527946737711E4+t12*t19*xij1*5.11068857403781E5+t12*t20*xij1*4.221994738626192E5)-t16*(t10*t11*t12*xij1*3.428291787990583E-2+t11*t12*t14*xij1*3.0)-t12*t31*xij1*8.125938683003058E-1-t29*t34*xij1+t12*t31*t33*xij1*5.0E1;
        u_xij[2,i] = -t12*(t12*t17*xij2*4.991172977107606E3+t12*t18*xij2*9.921527946737711E4+t12*t19*xij2*5.11068857403781E5+t12*t20*xij2*4.221994738626192E5)-t16*(t10*t11*t12*xij2*3.428291787990583E-2+t11*t12*t14*xij2*3.0)-t12*t31*xij2*8.125938683003058E-1-t29*t34*xij2+t12*t31*t33*xij2*5.0E1;
        u_xij[3,i] = -t12*(t12*t17*xij3*4.991172977107606E3+t12*t18*xij3*9.921527946737711E4+t12*t19*xij3*5.11068857403781E5+t12*t20*xij3*4.221994738626192E5)-t16*(t10*t11*t12*xij3*3.428291787990583E-2+t11*t12*t14*xij3*3.0)-t12*t31*xij3*8.125938683003058E-1-t29*t34*xij3+t12*t31*t33*xij3*5.0E1;
    end
    u = 0.5*u
    u_xij = 0.5*u_xij

    return u, u_xij
end

# function ezbl(r,zi,zj)
#     pzbl = 0.23;
#     a0 = 0.46850;
#     c1 = 0.02817;
#     c2 = 0.28022;
#     c3 = 0.50986;
#     c4 = 0.18175;
#     d1 = 0.20162;
#     d2 = 0.40290;
#     d3 = 0.94229;
#     d4 = 3.19980;
#     qqr2e = 14.399645;
#     qelectron = 1.0;

#     ainv = (sympy.Pow(zi,pzbl) + sympy.Pow(zj,pzbl))/(a0);
#     d1a = d1*ainv;
#     d2a = d2*ainv;
#     d3a = d3*ainv;
#     d4a = d4*ainv;
#     zze = zi*zj*qqr2e*qelectron*qelectron;

#     #[d1a d2a d3a d4a] #-> 2.3090    4.6140   10.7912   36.6444

#     e1 = exp(-d1a*r);
#     e2 = exp(-d2a*r);
#     e3 = exp(-d3a*r);
#     e4 = exp(-d4a*r);

#     rinv = 1/r;
#     sum = c1*e1;
#     sum = sum + c2*e2;
#     sum = sum + c3*e3;
#     sum = sum + c4*e4;

#     sum_p = -c1*d1a*e1;
#     sum_p = sum_p - c2*d2a*e2;
#     sum_p = sum_p - c3*d3a*e3;
#     sum_p = sum_p - c4*d4a*e4;

#     sum_pp = c1*e1*d1a*d1a;
#     sum_pp = sum_pp + c2*e2*d2a*d2a;
#     sum_pp = sum_pp + c3*e3*d3a*d3a;
#     sum_pp = sum_pp + c4*e4*d4a*d4a;

#     e = zze*sum*rinv;
#     ep = zze*(sum_p - sum*rinv)*rinv;
#     epp = zze*(sum_pp - 2.0*sum_p*rinv + 2.0*sum*rinv*rinv)*rinv;

#     return e, ep, epp
#   end

#   function zbl(xij, qi, qj, ti, tj, mu, eta, kappa)

#       r2 = xij[1,:].*xij[1,:] + xij[2,:]*xij[2,:] + xij[3,:]*xij[3,:];
#       r = sqrt.(r2);

#       cut_inner = 4;
#       cut_outer = 4.8;
#       zi = 73;
#       zj = 73;

#       # ZBL potential
#       u = ezbl(r,zi,zj)[1];

#       tc = cut_outer - cut_inner;
#       fc, fcp, fcpp = ezbl(cut_outer, zi, zj);
#       swa = (-3.0*fcp + tc*fcpp)/(tc*tc);
#       swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc);
#       swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp;

#       #sw1 = swa;
#       #sw2 = swb;
#       sw3 = swa/3.0;
#       sw4 = swb/4.0;
#       sw5 = swc;
#       #[sw3 sw4 sw5] # ->  [0.0456   -0.0343   -0.0163]

#       t = r - cut_inner;
#       sinner = sw5;
#       souter = sw5 + t*t*t * (sw3 + sw4*t);

#       # switching function
#       a = 0.5 + 0.5*tanh(1e2*t); # a = 0 when t<0, and a = 1 when t>0
#       s = a*souter + (1-a)*sinner;
#       u = u + s;
#       return u;
#   end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using RvSpectMLBase
using Test

@testset "ChunkList" begin
    using DataFrames
    num_pixels = 1000
    num_orders = 10
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    inst = TheoreticalInstrument2D(lambda2d)
    spectrum = TheoreticalInstrument.demo_generate_spectrum_line(inst)
    @test_nowarn make_orders_into_chunks(spectrum, inst)
    @test_nowarn make_orders_into_chunks(spectrum, orders_to_use=1:5, pixels_to_use=fill(20:80,5) )

    df = DataFrame(:lambda_lo=>[5501], :lambda_hi=>[5505])
    @test_nowarn make_chunk_list_around_lines(spectrum, df)
    cl = make_chunk_list_around_lines(spectrum, df)

    num_obs = 3
    times = sort(rand(num_obs))
    @test_skip ChunkListTimeseries(times, [cl, cl, cl])
end

ext: jl

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
# Autogenerated wrapper script for CUDA_loader_jll for x86_64-linux-gnu-cuda+11.1
export compute_sanitizer, libcublas, libcudadevrt, libcufft, libcupti, libcurand, libcusolver, libcusolverMg, libcusparse, libdevice, libnvtoolsext, libnvvm, nvdisasm, nvlink, ptxas

JLLWrappers.@generate_wrapper_header("CUDA_loader")
JLLWrappers.@declare_executable_product(compute_sanitizer)
JLLWrappers.@declare_library_product(libcublas, "libcublas.so.11")
JLLWrappers.@declare_file_product(libcudadevrt)
JLLWrappers.@declare_library_product(libcufft, "libcufft.so.10")
JLLWrappers.@declare_library_product(libcupti, "libcupti.so.11.1")
JLLWrappers.@declare_library_product(libcurand, "libcurand.so.10")
JLLWrappers.@declare_library_product(libcusolver, "libcusolver.so.11")
JLLWrappers.@declare_library_product(libcusolverMg, "libcusolverMg.so.11")
JLLWrappers.@declare_library_product(libcusparse, "libcusparse.so.11")
JLLWrappers.@declare_file_product(libdevice)
JLLWrappers.@declare_library_product(libnvtoolsext, "libnvToolsExt.so.1")
JLLWrappers.@declare_library_product(libnvvm, "libnvvm.so.3")
JLLWrappers.@declare_executable_product(nvdisasm)
JLLWrappers.@declare_executable_product(nvlink)
JLLWrappers.@declare_executable_product(ptxas)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        compute_sanitizer,
        "bin/compute-sanitizer",
    )

    JLLWrappers.@init_library_product(
        libcublas,
        "lib/libcublas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libcudadevrt,
        "lib/libcudadevrt.a",
    )

    JLLWrappers.@init_library_product(
        libcufft,
        "lib/libcufft.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcupti,
        "lib/libcupti.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcurand,
        "lib/libcurand.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolver,
        "lib/libcusolver.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolverMg,
        "lib/libcusolverMg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusparse,
        "lib/libcusparse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libdevice,
        "share/libdevice/libdevice.10.bc",
    )

    JLLWrappers.@init_library_product(
        libnvtoolsext,
        "lib/libnvToolsExt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvvm,
        "lib/libnvvm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nvdisasm,
        "bin/nvdisasm",
    )

    JLLWrappers.@init_executable_product(
        nvlink,
        "bin/nvlink",
    )

    JLLWrappers.@init_executable_product(
        ptxas,
        "bin/ptxas",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module JLLStats

import ..JLLExceptions: InputException

export nanmean;

"""
    nanmean(x; dims=nothing, count_nans=false)

Average an array (`x`) excluding NaNs. By default, the average is calculated as
the sum of all non-NaN values in `x` divided by the number of non-NaN values in `x`.
Setting `count_nans=true` will change it to be divided by the total number of
elements in `x`, include NaNs.

If dims is given, then `x` will be averaged along that dimension. Like `sum()`,
that dimension will not be removed in the output (it will have length = 1).

While using `missing` is generally preferred in Julia to mark values that are
actually missing and not the result of an undefined operation (e.g. `0. / 0.`),
currently `missing`s are not compatible with Unitful quantities.
"""
function nanmean(x; dims=nothing, count_nans=false)
    count_nans::Bool;

    notnans = .!isnan.(x);
    if dims == nothing
        # if no dim given, average over all elements
        n = count_nans ? length(x) : sum(notnans);
        return sum(x[notnans]) / n;
    else
        # require dims is an Integer
        dims :: Integer;
        if dims < 1 || dims > ndims(x)
            throw(InputException("x must be between 1 and ndims(x)"))
        end
        # if dim given, average along that element.
        n = count_nans ? size(x, dims) : sum(notnans, dims=dims);
        # this creates an anonymous function in the do block that's passed as the
        # first argument to mapslices, which iterates over x such that if e.g. x
        # is 3D and dims=2, slices x[1,:,1], x[1,:,2], x[2,:,1], etc. are taken
        # in turn. So this function gets each slice and adds up the non-NaN
        # values.
        s = mapslices(x; dims=(dims,)) do v
            vclean = v[.!isnan.(v)];
            sum(vclean);
        end
        return s ./ n;
    end
end

"""
    combinations(options::AbstractArray{<:Any,1}...)

Given multiple 1D arrays, returns an array of arrays with
all possible combinations of the input options. Example:

    julia> combinations([1.,2.,3.],["a","b","c"])
    9-element Array{Any,1}:
     Any[1.0, "a"]
     Any[2.0, "a"]
     Any[3.0, "a"]
     Any[1.0, "b"]
     Any[2.0, "b"]
     Any[3.0, "b"]
     Any[1.0, "c"]
     Any[2.0, "c"]
     Any[3.0, "c"]
"""
function combinations(options::AbstractArray{<:Any,1}...)
    dims = Tuple([length(o) for o in options]);
    combos = Array{Any,1}();
    for inds in CartesianIndices(dims)
        # CartesianIndices(dims) produces an array with dimensions dims
        # where each element gives its cartesian indices. The list
        # comprehension essentially takes each options vector as the
        # possible values along one dimension
        this_combo = [options[i][inds[i]] for i = 1:length(inds)];
        push!(combos, this_combo);
    end
    return combos;
end

"""
    normalize_to(data::AbstractArray{<:Integer}, range=(0.0, 1.0))
    normalize_to(data::AbstractArray{<:AbstractFloat}, range=(0.0, 1.0))

Return a copy of the input array, `data`, normalized to the specified range.
If `data` is an array of integers, then it is converted to floats internally.
Does not affect the original `data`.

Use the `range` keyword to specify a range to normalize to other that 0 to 1.
The value may be any object where `range[1]` returns the bottom of the range
and `range[2]` the top of it.

Examples:

    julia> x = collect(1:5)
    5-element Array{Int64,1}:
     1
     2
     3
     4
     5

    julia> normalize_to(x)
    5-element Array{Float64,1}:
     0.0
     0.25
     0.5
     0.75
     1.0

    julia> normalize_to(x, (50, 100))
    5-element Array{Float64,1}:
     50.0
     62.5
     75.0
     87.5
     100.0
"""
function normalize_to(data::AbstractArray{<:Integer}, range=(0.0, 1.0))
    data = convert.(Float64, data);
    return normalize_to!(data, range);
end # normalize_to(AbstractArray{<:Integer})

function normalize_to(data::AbstractArray{<:AbstractFloat}, range=(0.0, 1.0))
    data = copy(data);
    normalize_to!(data, range);
    return data;
end # normalize_to(AbstractArray{<:AbstractFloat})

"""
    normalize_to!(data::AbstractArray{<:AbstractFloat}, range=(0.0, 1.0))

An in-place version of `normalize_to`, this overwrites `data` with the normalized
version, but may be faster since it does not have to create a new array. Because
of that, this function requires you to give it an array containing floats. An integer
array will not be accepted, and trying to call this with an integer array will
result in a `MethodError`.

Example:
    julia> x = collect(1.0:5.0)
    5-element Array{Float64,1}:
     1.0
     2.0
     3.0
     4.0
     5.0

    julia> normalize_to!(x);
    julia> x
    5-element Array{Float64,1}:
     0.0
     0.25
     0.5
     0.75
     1.0
"""
function normalize_to!(data::AbstractArray{<:AbstractFloat}, range=(0.0, 1.0))
    # first make all values in data between 0 and 1. don't use -= to avoid
    data .-= minimum(data);
    data ./= maximum(data);

    # Now rescale to the requested range
    range_diff = range[2] - range[1];
    data .*= range_diff;
    data .+= range[1];

end # normalize_to!

end # module

max_stars_repo_name: joshua-laughner/Julia-JLLUtils

max_issues_repo_name: joshua-laughner/Julia-JLLUtils

max_forks_repo_name: joshua-laughner/Julia-JLLUtils


text:
@doc raw"""
    CircleGroup <: GroupManifold{Circle{ℂ},MultiplicationOperation}

The circle group is the complex circle ([`Circle(ℂ)`](@ref)) equipped with
the group operation of complex multiplication ([`MultiplicationOperation`](@ref)).
"""
const CircleGroup = GroupManifold{ℂ,Circle{ℂ},MultiplicationOperation}

CircleGroup() = GroupManifold(Circle{ℂ}(), MultiplicationOperation())

@inline function active_traits(f, M::CircleGroup, args...)
    return merge_traits(
        IsGroupManifold(M.op),
        IsDefaultMetric(EuclideanMetric()),
        HasBiinvariantMetric(),
        active_traits(f, M.manifold, args...),
        IsExplicitDecorator(),
    )
end

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

adjoint_action(::CircleGroup, p, X) = X

adjoint_action!(::CircleGroup, Y, p, X) = copyto!(Y, X)

function compose(
    ::MultiplicationGroupTrait,
    G::CircleGroup,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
)
    return map((pp, qq) -> compose(G, pp, qq), p, q)
end

function compose!(
    ::MultiplicationGroupTrait,
    G::CircleGroup,
    x,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
)
    return copyto!(x, compose(G, p, q))
end

identity_element(G::CircleGroup) = 1.0
identity_element(::CircleGroup, p::Number) = one(p)

Base.inv(G::CircleGroup, p::AbstractArray{<:Any,0}) = map(pp -> inv(G, pp), p)

function inverse_translate(
    ::CircleGroup,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
    ::LeftAction,
)
    return map(/, q, p)
end
function inverse_translate(
    ::CircleGroup,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
    ::RightAction,
)
    return map(/, q, p)
end

lie_bracket(::CircleGroup, X, Y) = zero(X)

lie_bracket!(::CircleGroup, Z, X, Y) = fill!(Z, 0)

translate_diff(::GT, p, q, X, ::ActionDirection) where {GT<:CircleGroup} = map(*, p, X)
function translate_diff(
    ::CircleGroup,
    ::Identity{MultiplicationOperation},
    q,
    X,
    ::ActionDirection,
)
    return X
end

function translate_diff!(G::CircleGroup, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, translate_diff(G, p, q, X, conv))
end

function exp_lie(::CircleGroup, X)
    return map(X) do imθ
        θ = imag(imθ)
        sinθ, cosθ = sincos(θ)
        return Complex(cosθ, sinθ)
    end
end

exp_lie!(G::CircleGroup, q, X) = (q .= exp_lie(G, X))

function _log_lie(::CircleGroup, q)
    return map(q) do z
        cosθ, sinθ = reim(z)
        θ = atan(sinθ, cosθ)
        return θ * im
    end
end

_log_lie!(G::CircleGroup, X, q) = (X .= _log_lie(G, q))

function number_of_coordinates(G::CircleGroup, B::AbstractBasis)
    return number_of_coordinates(base_manifold(G), B)
end

@doc raw"""
    RealCircleGroup <: GroupManifold{Circle{ℝ},AdditionOperation}

The real circle group is the real circle ([`Circle(ℝ)`](@ref)) equipped with
the group operation of addition ([`AdditionOperation`](@ref)).
"""
const RealCircleGroup = GroupManifold{ℝ,Circle{ℝ},AdditionOperation}

RealCircleGroup() = GroupManifold(Circle{ℝ}(), AdditionOperation())

@inline function active_traits(f, M::RealCircleGroup, args...)
    return merge_traits(
        IsGroupManifold(M.op),
        IsDefaultMetric(EuclideanMetric()),
        HasBiinvariantMetric(),
        active_traits(f, M.manifold, args...),
        IsExplicitDecorator(),
    )
end

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

is_default_metric(::RealCircleGroup, ::EuclideanMetric) = true

# Lazy overwrite since this is a rare case of nonmutating foo.
compose(::RealCircleGroup, p, q) = sym_rem(p + q)
compose(::RealCircleGroup, ::Identity{AdditionOperation}, q) = sym_rem(q)
compose(::RealCircleGroup, p, ::Identity{AdditionOperation}) = sym_rem(p)
function compose(
    ::RealCircleGroup,
    e::Identity{AdditionOperation},
    ::Identity{AdditionOperation},
)
    return e
end

compose!(::RealCircleGroup, x, p, q) = copyto!(x, sym_rem(p + q))
compose!(::RealCircleGroup, x, ::Identity{AdditionOperation}, q) = copyto!(x, sym_rem(q))
compose!(::RealCircleGroup, x, p, ::Identity{AdditionOperation}) = copyto!(x, sym_rem(p))
function compose!(
    ::RealCircleGroup,
    ::Identity{AdditionOperation},
    e::Identity{AdditionOperation},
    ::Identity{AdditionOperation},
)
    return e
end

identity_element(G::RealCircleGroup) = 0.0
identity_element(::RealCircleGroup, p::AbstractArray) = map(i -> zero(eltype(p)), p)

Base.inv(G::RealCircleGroup, p::AbstractArray{<:Any,0}) = map(pp -> inv(G, pp), p)

function inverse_translate(
    ::RealCircleGroup,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
    ::LeftAction,
)
    return map((x, y) -> sym_rem(x - y), q, p)
end
function inverse_translate(
    ::RealCircleGroup,
    p::AbstractArray{<:Any,0},
    q::AbstractArray{<:Any,0},
    ::RightAction,
)
    return map((x, y) -> sym_rem(x - y), q, p)
end

function exp_lie(::RealCircleGroup, X)
    return sym_rem(X)
end

exp_lie!(G::RealCircleGroup, q, X) = (q .= sym_rem(X))

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
#Aim: have map operate on FB

#Probably to be moved to AbstractAlgebra
function orbits(v::Vector{Perm{Int}})
  orbs = Vector{Vector{Int}}()
  n = length(v[1].d)
  found = falses(n)
  for i = 1:n
    if found[i]
      continue
    end
    found[i] = true
    orbit = Vector{Int}()
    push!(orbit, i)
    ind = 1
    while ind <= length(orbit)
      for j in v
        new = j[orbit[ind]]
        if !found[new]
          found[new] = true
          push!(orbit, new)
        end
      end
      ind += 1
    end
    push!(orbs, orbit)
  end
  return orbs
end

@doc Markdown.doc"""
    induce_action(primes::Vector{NfOrdIdl}, A::Map) -> Perm{Int}

Given a set of prime ideals invariant under the action of $A$, this function
returns the corresponding permutation induced by $A$.
"""
function induce_action(primes::Vector{NfOrdIdl}, A::Map)
  K = domain(A)
  f = A(gen(K)) # essentially a polynomial in the primitive element

  O = order(primes[1])
  prm = Vector{Tuple{Int, Int}}()

  G = SymmetricGroup(length(primes))
  if f == gen(K)
    return one(G)
  end

  primes_underneath = Set{fmpz}([minimum(x, copy = false) for x in primes])
  for p in primes_underneath
    indices = [i for i in 1:length(primes) if minimum(primes[i], copy = false) == p]
    lp = NfOrdIdl[primes[i] for i in indices]
    prm_p = _induce_action_p(lp, A)
    for (i, j) in prm_p
      push!(prm, (indices[i], indices[j]))
    end
  end
  sort!(prm, by = a -> a[1])
  return G([x[2] for x = prm])
end

#As above, but assumes that the prime ideals are lying over the same prime number p.
function _induce_action_p(lp::Vector{NfOrdIdl}, A::Map)
  O = order(lp[1])
  K = nf(O)
  prm = Tuple{Int, Int}[]
  p = minimum(lp[1], copy = false)

  if length(lp) < 3 || is_index_divisor(O, p) || !fits(Int, p)
    #TODO: Put some more thought. At least, do not check ideals that have already been found!
    found = falses(length(lp))
    for (i, P) in enumerate(lp)
      Q = induce_image(A, P)
      id = 0
      for j = 1:length(lp)
        if found[j]
          continue
        end
        if lp[j] == Q
          found[j] = true
          id = j
          break
        end
      end
      @assert !iszero(id)
      push!(prm, (i, id))
    end
  else
    px = PolynomialRing(GF(Int(p), cached=false), "x", cached=false)[1]
    fpx = px(A(gen(K)))
    gpx = px(K.pol)
    #idea/ reason
    # if p is no index divisor, then a second generator is just
    #   an irreducible factor of gpx (Kummer/ Dedekind)
    # an ideal is divisible by P iff the canonical 2nd generator of the prime ideal
    # divides the 2nd generator of the target (CRT)
    # so
    lpols = gfp_poly[gcd(px(K(P.gen_two)), gpx) for P in lp]
    # this makes lp canonical (should be doing nothing actually)

    for (i, P) in enumerate(lp)
      hp = px(K(P.gen_two))
      if degree(hp)==1
        im = fpx + coeff(hp, 0)
      else
        im = compose_mod(hp, fpx, gpx)
        # the image, directly mod p...
      end
      im = Hecke.gcd!(im, gpx, im)
      # canonical
      push!(prm, (i, findfirst(isequal(im), lpols)))
      # and find it!
    end
  end
  return prm
end

function induce(FB::Hecke.NfFactorBase, A::Map)
  K = domain(A)
  f = A(gen(K)) # essentially a polynomial in the primitive element

  O = order(FB.ideals[1])
  prm = Vector{Tuple{Int, Int}}()

  G = SymmetricGroup(length(FB.ideals))
  if f == gen(K)
    return one(G)
  end

  for p in FB.fb_int.base
    FP = FB.fb[p]
    lp = NfOrdIdl[x[2] for x = FP.lp]
    prm_p = _induce_action_p(lp, A)
    for (a, b) in prm_p
      push!(prm, (FP.lp[a][1], FP.lp[b][1]))
    end
  end
  sort!(prm, by = a -> a[1])
  return G([x[2] for x = prm])
end

#= implementation from Butler's Fundamntal Algorithm for Permutation Groups
  Algo 4: Dimino
  Tested for cyclic groups - unfortunately only.
  I still need to generate other input
=#
#function orbit_in_FB(op::Vector{Tuple{Map, Generic.Perm}}, a::nf_elem, s::SRow)
function orbit_in_FB(op::Array, a::nf_elem, s::SRow)
  function op_smat(n::SRow, p::Generic.Perm)
    r = [(p[i], v) for (i,v) = n]
    sort!(r, lt = (a,b)->a[1]<b[1])
    return typeof(n)(r)
  end

  Ss = Dict{nf_elem, typeof(s)}()
#  Ss = Dict{typeof(s), nf_elem}()
  Ss[a] = s
  # start with the cyclic group be op[1]
  n = s
  b = op[1][1](a)
  while b != a
    n = op_smat(n, op[1][2])
    Ss[b] = n
    b = op[1][1](b)
  end

  for i=2:length(op)
    bb = op[i][1](a)
    if haskey(Ss, bb)
      continue
    end
    old = collect(Ss)
    for (bb, n) in old # one redundant - step
      Ss[op[i][1](bb)] = op_smat(n, op[i][2])
    end
    while true
      done = true
      for j = 1:length(op)
        bb = op[j][1](b)
        if haskey(Ss, bb)
          continue;
        end
        done = false
        b = bb
        for (bbb, n) in old
          # TODO: What is going on with the b's?
          Ss[op[j][1](bbb)] = op_smat(n, op[j][2])
        end
      end
      if done
        break
      end
    end
  end
  return Ss
end

function generated_subgroup(op::Array) #pairs: permutations and Map
  elt = Vector{Any}()
  push!(elt, (x->x, parent(op[1][2])()))
  ord = 1
  g = op[1]
  while g[2] != parent(op[1][2])()
    let c_g = g
      push!(elt, c_g)
#      g = (x->op[1][1](c_g[1](x)), op[1][2]*c_g[2])
      g = (x->op[1][1](c_g[1](x)), c_g[2]*op[1][2])
    end
  end
  ord = length(elt)

  for i=2:length(op)
    c_i = i
    if op[i][2] in [x[2] for x=elt]
      continue
    end
    pord = ord
    push!(elt, op[i])
    for j=2:pord
      c_j = j
#      push!(elt, (x->elt[c_j][1](op[c_i][1](x)), elt[c_j][2]*op[c_i][2]))
      push!(elt, (x->elt[c_j][1](op[c_i][1](x)), op[c_i][2]*elt[c_j][2]))
    end
    ord = length(elt)
    rpos = pord + 1
    while true
      for s in op
        let c_rpos = rpos, c_s = s
#          g = (x->elt[c_rpos][1](c_s[1](x)), elt[c_rpos][2]*c_s[2])
          g = (x->elt[c_rpos][1](c_s[1](x)), c_s[2]*elt[c_rpos][2])
          if g[2] in [x[2] for x=elt]
            continue
          end
        end
        let c_g = g
          push!(elt, c_g)
          for j = 2:pord
            c_j = j
#            push!(elt, (x->elt[c_j][1](c_g[1](x)), elt[c_j][2]*c_g[2]))
            push!(elt, (x->elt[c_j][1](c_g[1](x)), c_g[2]*elt[c_j][2]))
          end
        end
        ord = length(elt)
      end
      rpos += pord
      if rpos > length(elt)
        break
      end
    end
  end
  return elt
end

function class_group_add_auto(ctx::ClassGrpCtx, auts::Vector{NfToNfMor})
  #I am assuming that auts contains all the automorphisms of K
  K = domain(auts[1])
  p = 11
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  while degree(fmod) != degree(K) || !is_squarefree(fmod)
    p = next_prime(p)
    R = GF(p, cached = false)
    Rx, x = PolynomialRing(R, "x", cached = false)
    fmod = Rx(K.pol)
  end
  S = small_generating_set(auts)

  Dpols = Dict{gfp_poly, NfToNfMor}()
  for i = 1:length(auts)
    Dpols[Rx(image_primitive_element(auts[i]))] = auts[i]
  end
  Gperm = SymmetricGroup(length(ctx.FB.ideals))

  elements = Vector{Tuple{NfToNfMor, Generic.Perm{Int}}}(undef, length(auts))
  elements[1] = (id_hom(K), one(Gperm))
  if length(auts) == 1
    return elements
  end
  ind_elem = 3
  pols = gfp_poly[x, Rx(image_primitive_element(S[1]))]
  perms = Generic.Perm{Int}[one(Gperm), induce(ctx.FB, S[1])]
  elements[2] = (S[1], perms[2])
  gperm = perms[2]*perms[2]
  gpol = compose_mod(pols[2], pols[end], fmod)
  while gpol != pols[1]
    elements[ind_elem] = (Dpols[gpol], gperm)
    @hassert :ClassGroup 1 induce(ctx.FB, elements[ind_elem][1]) == elements[ind_elem][2]
    ind_elem += 1
    push!(pols, gpol)
    push!(perms, gperm)
    gperm = perms[2]*gperm
    gpol = compose_mod(pols[2], pols[end], fmod)
  end
  if length(pols) == length(auts)
    return elements
  end

  for i in 2:length(S)
    pi = Rx(image_primitive_element(S[i]))
    if !(pi in pols)
      permi = induce(ctx.FB, S[i])
      previous_order = length(pols)
      elements[ind_elem] = (S[i], permi)
      push!(pols, pi)
      push!(perms, permi)
      ind_elem += 1
      for j in 2:previous_order
        push!(pols, compose_mod(pols[j], pi, fmod))
        push!(perms, perms[j]*permi)
        elements[ind_elem] = (Dpols[pols[end]], perms[end])
        @hassert :ClassGroup 1 induce(ctx.FB, elements[ind_elem][1]) == elements[ind_elem][2]
        ind_elem += 1
      end
      if length(pols) == length(auts)
        return elements
      end
      rep_pos = previous_order + 1
      while rep_pos <= length(pols)
        for k in 1:i
          s = S[k]
          po = Rx(image_primitive_element(s))
          permo = elements[1][2]
          for l = 1:length(elements)
            if s == elements[l][1]
              permo = elements[l][2]
              break
            end
          end
          permatt = perms[rep_pos]*permo
          att = compose_mod(pols[rep_pos], po, fmod)
          if !(att in pols)
            push!(perms, permatt)
            push!(pols, att)
            elements[ind_elem] = (Dpols[pols[end]], perms[end])
            @hassert :ClassGroup 1 induce(ctx.FB, elements[ind_elem][1]) == elements[ind_elem][2]
            ind_elem += 1
            for j in 2:previous_order
              push!(pols, compose_mod(pols[j], att, fmod))
              push!(perms, perms[j]*permatt)
              elements[ind_elem] = (Dpols[pols[end]], perms[end])
              @hassert :ClassGroup 1 induce(ctx.FB, elements[ind_elem][1]) == elements[ind_elem][2]
              ind_elem += 1
            end
            if length(pols) == length(elements)
              return elements
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  return elements
end


text:
#!/usr/bin/env julia

using PortAudio
using Test

@testset "PortAudio Tests" begin
    @testset "Reports version" begin
        io = IOBuffer()
        PortAudio.versioninfo(io)
        result = split(String(take!((io))), "\n")
        # make sure this is the same version I tested with
        @test startswith(result[1], "PortAudio V19")
    end

    @testset "Can list devices without crashing" begin
        PortAudio.devices()
    end
end


text:
# NOTE: The code in this file comes from SNOW.jl, eventually it might be removed
# from here and imported from there.

"""
    SparsePattern(A::SparseMatrixCSC)
construct sparse pattern from representative sparse matrix
# Arguments
- `A::SparseMatrixCSC`: sparse jacobian
"""
function SparsePattern(A::SparseMatrixCSC)
    rows, cols, _ = findnz(A)
    return SparsePattern(rows, cols)
end

"""
    SparsePattern(A::Matrix)
construct sparse pattern from representative matrix
# Arguments
- `A::Matrix`: sparse jacobian
"""
function SparsePattern(A::Matrix)
    return SparsePattern(sparse(A))
end

"""
    SparsePattern(::ForwardAD, func!, ng, x1, x2, x3)
detect sparsity pattern by computing derivatives (using forward AD)
at three different locations. Entries that are zero at all three
spots are assumed to always be zero.
# Arguments
- `func!::Function`: function of form f = func!(g, x)
- `ng::Int`: number of constraints
- `x1,x2,x3::Vector{Float}`:: three input vectors.
"""
function SparsePattern(::ForwardAD, func!, ng, x1, x2, x3)

    g = zeros(ng)
    config = ForwardDiff.JacobianConfig(func!, g, x1)
    J1 = ForwardDiff.jacobian(func!, g, x1, config)
    J2 = ForwardDiff.jacobian(func!, g, x2, config)
    J3 = ForwardDiff.jacobian(func!, g, x3, config)
    @. J1 = abs(J1) + abs(J2) + abs(J3)
    Jsp = sparse(J1)

    return SparsePattern(Jsp)
end

"""
    SparsePattern(::AbstractFD, func!, ng, x1, x2, x3)
detect sparsity pattern by computing derivatives (using finite differencing)
at three different locations. Entries that are zero at all three
spots are assumed to always be zero.
# Arguments
- `func!::Function`: function of form f = func!(g, x)
- `ng::Int`: number of constraints
- `x1,x2,x3::Vector{Float}`:: three input vectors.
"""
function SparsePattern(dtype::AbstractFD, func!, ng,
    x1, x2, x3)

    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.JacobianCache(x1, zeros(ng), fdtype)

    nx = length(x1)
    J1 = zeros(ng, nx)
    J2 = zeros(ng, nx)
    J3 = zeros(ng, nx)
    FiniteDiff.finite_difference_jacobian!(J1, func!, x1, cache)
    FiniteDiff.finite_difference_jacobian!(J2, func!, x2, cache)
    FiniteDiff.finite_difference_jacobian!(J3, func!, x3, cache)

    @. J1 = abs(J1) + abs(J2) + abs(J3)
    Jsp = sparse(J1)

    return SparsePattern(Jsp)
end

"""
    SparsePattern(diffmethod, func!, ng, lx, ux)
detect sparsity pattern by computing derivatives
at three randomly generating locations w/in bounds.
Entries that are zero at all three spots are assumed to always be zero.
# Arguments
- `diffmethod::AbstractDiffMethod`: method to compute derivatives
- `func!::Function`: function of form f = func!(g, x)
- `ng::Int`: number of constraints
- `lx::Vector{Float}`: lower bounds on x
- `ux::Vector{Float}`: upper bounds on x
"""
function SparsePattern(diffmethod, func!, ng, lx, ux)
    r = rand(3)
    x1 = @. (1-r[1])*lx + r[1]*ux
    x2 = @. (1-r[2])*lx + r[2]*ux
    x3 = @. (1-r[3])*lx + r[3]*ux
    return SparsePattern(diffmethod, func!, ng, x1, x2, x3)
end

"""
    getsparsity(::DensePattern, nx, nf)

Get rows and cols for dense jacobian
"""
function getsparsity(::DensePattern, nx, nf)
    len = nf*nx
    rows = [i for i = 1:nf, j = 1:nx][:]
    cols = [j for i = 1:nf, j = 1:nx][:]
    return rows, cols
end

"""
    getsparsity(sp::SparsePattern, nx, nf)

Get rows and cols for sparse jacobian
"""
getsparsity(sp::SparsePattern, nx, nf) = sp.rows, sp.cols


text:
module Cos

#using ..Common:ispdfdigit,ispdfspace,PERCENT,MINUS_SIGN,DIGIT_ZERO,PERIOD,
#               TRAILER
#using ..Common.Parser

using ..Common

include("CosObject.jl")
include("CosObjectHelpers.jl")
include("CosStream.jl")
include("CosReader.jl")
include("CosObjStream.jl")
include("CosDoc.jl")
include("CosCrypt.jl")
include("StdSecHandler.jl")
include("PKISecHandler.jl")
end


text:
module CacheVariables

using BSON
import Logging: @info

export @cache

function _cachevars(ex::Expr)
    (ex.head === :(=))   && return Symbol[ex.args[1]]
    (ex.head === :block) && return collect(
        Iterators.flatten([_cachevars(exi) for exi in ex.args if isa(exi,Expr)])
    )
    return Vector{Symbol}(undef,0)
end

macro cache(path, ex::Expr, overwrite=false)
    vars = _cachevars(ex)
    vardesc = join(string.(vars), "\n")
    varkws  = [:($(var) = $(var)) for var in vars]
    varlist = :($(varkws...),)
    vartuple = :($(vars...),)

    return quote
        if !isfile($(esc(path))) || $(esc(overwrite))
            _ans = $(esc(ex))
            _msg = isfile($(esc(path))) ? "Overwriting " : "Saving to "
            @info(string(_msg, $(esc(path)), "\n", $(vardesc)))
            bson($(esc(path)); $(esc(varlist))...,ans=_ans)
            _ans
        else
            @info(string("Loading from ", $(esc(path)), "\n", $(vardesc)))
            data = BSON.load($(esc(path)))
            $(esc(vartuple)) = getindex.(Ref(data), $vars)
            data[:ans]
        end
    end
end

end # module


text:
"""
    sample_observables(x₀, sampler, observables; options=MDOptions())

Run the `sampler` starting at `x₀`, evaluating the trajectory on a tuple of
`observables` scalar functions.  Number of iterations and interval between saves
are set using the `options` argument.  Only the computed observables are returned.

### Fields
* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
* `observables` - Observables on which to evaluate the trajectory
### Optional Fields
* `TO`- Observable data type, if needed, should be entered as the first argument
* `options`   - Sampling options, including number of iteration
"""
function sample_observables(x₀::Tx, sampler::S, observables::Tuple{Vararg{<:Function,NO}}; options=MDOptions()) where {Tx,  S<:AbstractSampler, NO}

    state = InitState(x₀, sampler);
    observable_samples = zeros(NO, options.n_save);
    save_index = 1;
    for i = 1:options.n_iters
        UpdateState!(state, sampler);
        if(mod(i,options.n_save_iters)==0)
            ntuple(k-> observable_samples[k,save_index] =(observables[k])(state.x), NO);
            # Base.Cartesian.@nexprs $NO k -> observable_samples[k,save_index] = (observables[k])(state.x);
            save_index+=1;
        end
    end
    return observable_samples

end

function sample_observables(TO::Type, x₀::Tx, sampler::S, observables::Tuple{Vararg{<:Function,NO}}; options=MDOptions()) where {Tx,  S<:AbstractSampler, NO}

    state = InitState(x₀, sampler);
    observable_samples = zeros(TO, NO, options.n_save);
    save_index = 1;
    for i = 1:options.n_iters
        UpdateState!(state, sampler);
        if(mod(i,options.n_save_iters)==0)
            ntuple(k-> observable_samples[k,save_index] =(observables[k])(state.x), NO);
            # Base.Cartesian.@nexprs $NO k -> observable_samples[k,save_index] = convert(TO,(observables[k])(state.x));
            save_index+=1;
        end
    end
    return observable_samples
end

"""
    sample_observables(x₀, sampler, observables; options=MDOptions())

Run the `sampler` starting at `x₀`, evaluating the trajectory on a tuple of
`observables` scalar functions.  Number of iterations and interval between saves
are set using the `options` argument.  Only the computed observables are returned.

### Fields
* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
* `recycler`  - Recycling structure for A→B transitions
* `observables` - Observables on which to evaluate the trajectory
### Optional Fields
* `TO`- Observable data type, if needed, should be entered as the first argument
* `options`   - Sampling options, including number of iteration
"""
function sample_observables(x₀::Tx, sampler::S,recycler::R, observables::Tuple{Vararg{<:Function,NO}}; options=MDOptions()) where {Tx,  S<:AbstractSampler,R<:AbstractRecycler, NO}

    state = InitState(x₀, sampler);
    observable_samples = zeros(NO, options.n_save);
    save_index = 1;
    for i = 1:options.n_iters
        if(mod(i-1,recycler.n_recycle_iters)==0)
            if (recycler.inB(state))
                recycler.restartA!(state)
            end
        end
        UpdateState!(state, sampler);
        if(mod(i,options.n_save_iters)==0)
            ntuple(k-> observable_samples[k,save_index] =(observables[k])(state.x), NO);
            # Base.Cartesian.@nexprs $NO k -> observable_samples[k,save_index] = (observables[k])(state.x);
            save_index+=1;
        end
    end
    return observable_samples
end

function sample_observables(TO::Type, x₀::Tx, sampler::S,recycler::R, observables::Tuple{Vararg{<:Function,NO}}; options=MDOptions()) where {Tx,  S<:AbstractSampler,R<:AbstractRecycler, NO}

    state = InitState(x₀, sampler);
    observable_samples = zeros(TO, NO, options.n_save);
    save_index = 1;
    for i = 1:options.n_iters
        if(mod(i-1,recycler.n_recycle_iters)==0)
            if (recycler.inB(state))
                recycler.restartA!(state)
            end
        end
        UpdateState!(state, sampler);
        if(mod(i,options.n_save_iters)==0)
            ntuple(k-> observable_samples[k,save_index] =(observables[k])(state.x), NO);
            # Base.Cartesian.@nexprs $NO k -> observable_samples[k,save_index] = convert(TO,(observables[k])(state.x));
            save_index+=1;
        end
    end
    return observable_samples
end

max_stars_repo_name: gideonsimpson/JuBasicMD

max_issues_repo_name: gideonsimpson/JuBasicMD

max_forks_repo_name: gideonsimpson/JuBasicMD


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

using Markdown
using InteractiveUtils

# ╔═╡ d02ecf8a-fdd0-11ea-0bec-b561da2f2a02
using DataFrames,CSV,Test,PovertyAndInequalityMeasures,Plots,Test,StatsBase

# ╔═╡ 97c87a0a-fde0-11ea-1419-2f9939bfa6e0
md"""

# Building A MicroSimulation Tax-Benefit Model of Scotland (in 25 minutes)

### Graham Stark
[graham.stark@virtual-worlds-research.com](mailto:graham.stark@virtual-worlds-research.com)

Things we need:

* data - Family Resources Survey(FRS)/LFS/Living Costs and Food Survey(LCF)/Understanding Society... This example uses FRS
* A programming language:
   - Spreadsheets, R, Fortran, Python ..
   - We use [Julia](https://julialang.org)
* Structures: we need to model:
   - people, families and households;
   - the fiscal system (taxes, benefits);
   - outcomes: incomes net of taxes and benefits, revenues raised, inequalities, poverty and so on.

Let's build one from scratch, in real time (with some 'here's some I prepared earlier' where needs be).

"""

# ╔═╡ d1e09cc6-fde1-11ea-07c7-0fce2411a060
md"""
## A quick intro to julia.

This web page is created automatically by Julia.

"""

# ╔═╡ 00f12bde-fde2-11ea-1026-c11c72cea702
2+2

# ╔═╡ 0651656c-fde2-11ea-251f-5747f4746430
x = 24

# ╔═╡ 0a7f7e3a-fde2-11ea-3825-37b7be64f5e2
x

# ╔═╡ 0c3b3004-fde2-11ea-0567-4bd9576443af
function f(x)
	x*x-1
end

# ╔═╡ 2b19845e-fde2-11ea-2353-7ddeb1182fc9
f(10)

# ╔═╡ 53b4ca1c-fde2-11ea-2955-f3c19026972b
struct Dog
	age
	sex
	barks
end

# ╔═╡ 708ab7c8-fde2-11ea-32ac-b56e9df87ebe
hope = Dog(4,"Female",true)

# ╔═╡ 83d83a9c-fde2-11ea-3512-9f140251a4e9
hope.barks

# ╔═╡ 940fd034-fde2-11ea-1fe2-97533b9dc40f
hope.age == 4

# ╔═╡ a8ff38e0-45c7-11eb-39ac-119d98bad749
	struct Pers
		age :: Int
		income :: Float64
		sex :: Int
	end

# ╔═╡ e8aa0fae-45c7-11eb-03b2-cf5021684e7e
g = Pers(62,25_000.0,1)

# ╔═╡ d871fa52-fde2-11ea-1e58-a9ddbea34508
md"""
## Data

this is the 'here's one I prepared earlier' bit... This is a large extract from the Family Resources Survey (pooled 2015-2018 surveys, Scotland only households).

Julia provides a structure called a [DataFrame](https://juliadata.github.io/DataFrames.jl/stable/) to hold large datasets like this.

"""

# ╔═╡ e8d7b7a4-fdd0-11ea-3e4e-d3cfcd789742
sp=CSV.File( "/home/graham_s/julia/vw/ScottishTaxBenefitModel/data/model_people_scotland.tab")|>DataFrame

# ╔═╡ 28ce666e-fdd1-11ea-2ae6-395356fc23e5
begin
    sp[((sp.age.>60).&(sp.sex.==1).&(sp.income_self_employment_income.>0)),[:age,:usual_hours_worked,:income_wages,:income_self_employment_income]]
	# sp = tsp[tsp.from_child_record.!=1,:]
end

# ╔═╡ d9eb5ebc-fdd2-11ea-28ab-2fa7c0cfaba9
names(sp)

# ╔═╡ 4f4ef062-fde3-11ea-1a36-1964d61271ca
md"""
## Structures

As I mentioned, we need structures for households, people,the tax and benefit system, and the results of our calculations. A good model is mostly structures and tests - get those right and the actual calculations are usually easy.

We'll skip most things for now; just:
* people (just age and income)
* simple Income Tax
* a results structure.

"""

# ╔═╡ 5ce6e028-fdd3-11ea-2924-897de20ca1c2
begin
	struct Person
		age  :: Int
		wage :: Real
	end

	struct TaxSystem
		rates :: Vector
		bands :: Vector
		allow :: Real
	end

	struct Result
		taxes :: Real
		net_income :: Real
	end

end

# ╔═╡ ba302158-fde3-11ea-2bae-939a6cab3511
md"""

### Initialise some tax systems

The current Scottish one and a flat-rate one.

"""

# ╔═╡ e118d6a8-fdd3-11ea-2d61-ad17c31e0d1a
begin
	sys1=TaxSystem(
		[0.190,0.20,0.21,0.41,0.46], # rates
		[2_049.0, 12_444.0, 30_930.0, 150_000.0], # thresholds
		12_500.00) # allowance
	# flat - no allowance, 1 rate
	sys2=TaxSystem([0.1475],[],0)
end

# ╔═╡ eff4a08e-fde3-11ea-2e58-9b4e4969819e

# ╔═╡ 055e7cc4-fde4-11ea-145b-c5261a647e10
md"""

### A very simple tax system

"""

# ╔═╡ e89d49e0-fdd3-11ea-199c-2d4bf8a5563b
begin

	"""
Tax due on `taxable` income, given rates and thresholds
rates can be one more than thresholds, in which case the last band is assumed infinite.
Rates should be (e.g.) 0.12 for 12%.
"""
function calctaxdue(
   taxable    :: Number,
   rates      :: Vector,
   thresholds :: Vector ) :: Real
   nthresholds = length(thresholds)[1]
   nrates = length(rates)[1]

   @assert (nrates >= 1) && ((nrates - nthresholds) in 0:1 ) # allow thresholds to be 1 less & just fill in the top if we need it
   due = 0.0
   remaining = taxable
   i = 0
   if nthresholds > 0
      maxv = typemax( typeof( thresholds[1] ))
      gap = thresholds[1]
   else
      maxv = typemax( typeof( taxable ))
      gap = maxv
   end
   while remaining > 0.0
      i += 1
      if i > 1
         if i < nrates
            gap = thresholds[i]-thresholds[i-1]
         else
            gap = maxv
         end
      end
      t = min( remaining, gap )
      due += t*rates[i]
      remaining -= gap
   end
   due
end

	function dotax( pers::Person, sys :: TaxSystem )::Result
		taxable = max(0.0, pers.wage-sys.allow)
		tax = calctaxdue( taxable, sys.rates, sys.bands )
		net = pers.wage - tax
		return Result(tax,net)
	end
end

# ╔═╡ 44f853e2-fdd4-11ea-1a13-d5740807f2d0
begin
	pers = Person(40,20_000.0)
	res = dotax(pers,sys1)
end

# ╔═╡ 6a698960-fdd5-11ea-3b8d-81aaf65a4d22
tt = (20_000-12500)*0.19

# ╔═╡ a5f30538-fde9-11ea-39c1-1f7f412d2c04
md"""
### Sample Weights

`target population/people (rows)` in the dataset.

"""

# ╔═╡ 7c9954a6-fdd5-11ea-32cc-819552e44bd2
begin
 	ss = size(sp)[1] # count of rows
	popn = 5_500_000 # approx popn of Scotland
	w = popn/ss
end

# ╔═╡ b265ebfa-fdd5-11ea-10fd-2f02021f86a0

function do_all(data::DataFrame, sys :: TaxSystem)::Matrix
	tot_tax = 0.0
	ss = size(data)[1]
	popn = 5_500_000
	w = popn/ss
	out = zeros(ss,3)
	inflation = 1.1
	i = 0
	for p in eachrow(data)
		i+=1
		pers = Person( p.age,p.income_wages*52*inflation)
		res = dotax(pers,sys)
		out[i,1]=w
		out[i,2]=res.net_income
		out[i,3]=res.taxes
	end
	out
end

# ╔═╡ a1f92cf4-fdd6-11ea-37ba-5b4efab11252
begin
	allr1=do_all(sp,sys1)
	allr2=do_all(sp,sys2)
end

# ╔═╡ d16e31fa-fdd6-11ea-3476-71dbfa596605
begin
	ineq1 = PovertyAndInequalityMeasures.make_inequality( allr1, 1, 2 )
	ineq2 = PovertyAndInequalityMeasures.make_inequality( allr2, 1, 2 )
end

# ╔═╡ 8fd8dc52-fdd8-11ea-0d46-db0fbd96cd90
begin
	dc1 = ineq1[:deciles]
	dc2 = ineq2[:deciles]
	m = hcat(dc1[:,1:2],dc2[:,2])
	plot(m,labels=["equality" "pre" "post"],title="Gini")
	#m
end

# ╔═╡ 8f1e6616-ffb2-11ea-0e84-f35d2f98221b
begin
	decs = (dc2[:,3] - dc1[:,3])./52
	bar(decs,labels="avg. gain",title="deciles")
end

# ╔═╡ 5667cce2-fdda-11ea-1eb2-130e04fd8cb3
begin
	r1=Int(trunc(sum(allr1[:,1].*allr1[:,3])/1_000_000))
	r2=Int(trunc(sum(allr2[:,1].*allr2[:,3])/1_000_000))
	cost = r2-r1
	(cost=cost,r1=r1,r2=r2)
end

# ╔═╡ 8f3045e0-fdda-11ea-2f94-0dec0b980c14
(ineq1[:gini],ineq2[:gini],ineq1[:gini]-ineq2[:gini])

# ╔═╡ f4a96ad2-fde9-11ea-28c3-f9376758d82f
@testset "X" begin
	@test 1 ≈ 1.0000000000000000000008
end

# ╔═╡ c985e5d6-0084-11eb-2569-c3d8af344837
combine(sp,:income_wages=>mean)[1,:income_wages_mean]*52|>trunc|>Integer

# ╔═╡ 20950768-0093-11eb-3e65-391d2057f808
last(select(sp, :income_wages => tiedrank,:income_wages => ordinalrank,:income_wages),9)

# ╔═╡ Cell order:
# ╠═d02ecf8a-fdd0-11ea-0bec-b561da2f2a02
# ╟─97c87a0a-fde0-11ea-1419-2f9939bfa6e0
# ╟─d1e09cc6-fde1-11ea-07c7-0fce2411a060
# ╠═00f12bde-fde2-11ea-1026-c11c72cea702
# ╠═0651656c-fde2-11ea-251f-5747f4746430
# ╠═0a7f7e3a-fde2-11ea-3825-37b7be64f5e2
# ╠═0c3b3004-fde2-11ea-0567-4bd9576443af
# ╠═2b19845e-fde2-11ea-2353-7ddeb1182fc9
# ╠═53b4ca1c-fde2-11ea-2955-f3c19026972b
# ╠═708ab7c8-fde2-11ea-32ac-b56e9df87ebe
# ╠═83d83a9c-fde2-11ea-3512-9f140251a4e9
# ╠═940fd034-fde2-11ea-1fe2-97533b9dc40f
# ╠═a8ff38e0-45c7-11eb-39ac-119d98bad749
# ╠═e8aa0fae-45c7-11eb-03b2-cf5021684e7e
# ╟─d871fa52-fde2-11ea-1e58-a9ddbea34508
# ╠═e8d7b7a4-fdd0-11ea-3e4e-d3cfcd789742
# ╠═28ce666e-fdd1-11ea-2ae6-395356fc23e5
# ╠═d9eb5ebc-fdd2-11ea-28ab-2fa7c0cfaba9
# ╟─4f4ef062-fde3-11ea-1a36-1964d61271ca
# ╠═5ce6e028-fdd3-11ea-2924-897de20ca1c2
# ╠═ba302158-fde3-11ea-2bae-939a6cab3511
# ╠═e118d6a8-fdd3-11ea-2d61-ad17c31e0d1a
# ╠═eff4a08e-fde3-11ea-2e58-9b4e4969819e
# ╟─055e7cc4-fde4-11ea-145b-c5261a647e10
# ╠═e89d49e0-fdd3-11ea-199c-2d4bf8a5563b
# ╠═44f853e2-fdd4-11ea-1a13-d5740807f2d0
# ╠═6a698960-fdd5-11ea-3b8d-81aaf65a4d22
# ╠═a5f30538-fde9-11ea-39c1-1f7f412d2c04
# ╠═7c9954a6-fdd5-11ea-32cc-819552e44bd2
# ╠═b265ebfa-fdd5-11ea-10fd-2f02021f86a0
# ╠═a1f92cf4-fdd6-11ea-37ba-5b4efab11252
# ╠═d16e31fa-fdd6-11ea-3476-71dbfa596605
# ╠═8fd8dc52-fdd8-11ea-0d46-db0fbd96cd90
# ╠═8f1e6616-ffb2-11ea-0e84-f35d2f98221b
# ╠═5667cce2-fdda-11ea-1eb2-130e04fd8cb3
# ╠═8f3045e0-fdda-11ea-2f94-0dec0b980c14
# ╠═f4a96ad2-fde9-11ea-28c3-f9376758d82f
# ╠═c985e5d6-0084-11eb-2569-c3d8af344837
# ╠═20950768-0093-11eb-3e65-391d2057f808

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Rankers.Storages
using LOLTools: SummonerV4, LeagueV4, http_action
using Rankers: Chart

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"

read_from_file = true # false
queue = "RANKED_SOLO_5x5"
league = LeagueV4.challengerleagues(api_key, region, queue; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
@info map(x -> (x.summonerName, x.leaguePoints), Chart.entries(league)[1:10])

using Millboard
table(Chart.entries(league)[1:20])

summonerName = "DRX Deft"
summonerdto = SummonerV4.summoner_by_name(api_key, region, summonerName; action=store_with(region, http_action, debug=true))
encryptedSummonerId = summonerdto.id

entries = LeagueV4.entries_by_summoner_id(api_key, region, encryptedSummonerId; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
@info entries

queue, tier, division = "RANKED_SOLO_5x5", "DIAMOND", "I"
entries = LeagueV4.entries_by_queue_tier_division(api_key, region, queue, tier, division; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
@info first(entries)


text:
# Question 1: Modifisere binærsøk
## Binærsøk er en enkel, men viktig algoritme med mange bruksområder. Les om den i appendiks C («Binærsøk») i pensumheftet.
## Under vises julia-koden for funksjonen bisect_left som er en variant av binærsøk. Denne har noen forskjeller sammenliknet med Bisect-algoritmen fra pensumheftet. Dersom det finnes flere forekomster av v i tabellen, returneres en indeks slik at om v settes inn på denne indeksen, havner elementet direkte til venstre for de andre forekomstene. Dersom v ikke finnes i A er indeksen som returneres slik at v kan settes inn på denne indeksen i sortert rekkefølge.

function bisect_left(A, p, r, v)
    i = p
    if p < r
       q = floor(Int, (p+r)/2)  # floor() er en innebygd funksjon som runder ned. ceil() runder opp.
       if v <= A[q]
           i = bisect_left(A, p, q, v)
       else
           i = bisect_left(A, q+1, r, v)
       end
    end
    return i
end

#Implementer funksjonen => bisect_right(A, p, r, v)

function bisect_right(A, p, r, v)
    i = p
    if p < r
        q = floor(Int, (p+r)/2)
        if v < A[q]
            i = bisect_right(A, p, q, v)
       	else
            i = bisect_right(A, q+1, r, v)
        end
    end
    return i
end

A = [0,1,1,1,2,3,4,5,6,7,8,9]

V = 1
#@info "V" V
#@info "index where we want to insert V"  bisect_left(A,1,length(A)+1,V)
#@info "The list is still" A

@info "Bisect_left" bisect_left(A,1,length(A)+1,V) "should be 2"
@info "Bisect_right" bisect_right(A,1,length(A)+1,V) "Should be 5"


text:
@testset "Measure" begin
    Mod.@polyvar x y
    @test_throws ArgumentError measure([1, 2], [x, x*y, y])
    @test_throws ArgumentError measure([1, 2, 3, 4], [x, x*y, y])
    μ = measure([1, 0, 2, 3], [x^2*y^2, y*x^2, x*y*x^2, x*y^2])
    @test base_functions(μ) == MonomialBasis(monovec([x^3*y, x^2*y^2, x^2*y, x*y^2]))
    @test base_function.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
    @test MultivariateMoments.moment_value.(moments(μ)) == [2, 1, 0, 3]
    @test all(nvariables.(moments(μ)) .== 2)
    @test degree.(moments(μ)) == [4, 4, 3, 3]
    @test μ.values == [2, 1, 0, 3]
    @test (-μ).values == [-2, -1, 0, -3]
    @test (2 * μ).values == [4, 2, 0, 6]
    @test (μ * 3).values == [6, 3, 0, 9]
    #@test_throws ArgumentError measure([1], [x]) + measure([1], [y])

    @testset "Measure MB" begin
        Mod.@polyvar x y
        @test_throws AssertionError measure([1, 0, 2, 3, 0], maxdegree_basis(ChebyshevBasis, [x, y], 2))
        @test_throws AssertionError measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 3))
        μ = measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 2))
        @test base_functions(μ) ==  maxdegree_basis(ChebyshevBasis, [x, y], 2)
        @test base_function.(moments(μ)) == [2.0x^2 - 1.0, x*y, 2.0y^2 - 1.0, x, y, 1.0]
        @test MultivariateMoments.moment_value.(moments(μ)) == [1, 0, 2, 3, 0, 4]
        @test all(nvariables.(moments(μ)) .== 2)
        @test degree.(moments(μ)) == [2, 2, 2, 1, 1, 0]
        @test μ.values == [1, 0, 2, 3, 0, 4]
        @test (-μ).values == -[1, 0, 2, 3, 0, 4]
        @test (2 * μ).values == [1, 0, 2, 3, 0, 4].*2
        @test (μ * 3).values == [1, 0, 2, 3, 0, 4].*3

    end
end


text:
# Test the simple 3D BEM implementation and asymptotic compression
# See popsomer/bempp/meshes for additional meshes
using PyCall
@pyimport matplotlib.pyplot as plt
nrMesh = 0
k = 15

include("bea.jl")
include("window.jl")

(xyz, elmt) = readMesh(string("sphere", nrMesh, ".msh") )
# Prints to ignore first enty which is origin:
xyz = xyz[2:end,:]

A = @time(getSystemMatrix(xyz,elmt,k))
b = -exp(1im*k*xyz[:,1])

c = A\b
f = open(string("ck", k, "sph", nrMesh), "w") # save results
writedlm(f,[real(c) imag(c)])
close(f)
resid = A*c-b

print("Potential at -1 0 0 is ", potential(xyz,elmt,k,[-1 0 0],c), k, " =k, nrMesh = ", nrMesh, "\n\n")
checkSolution(xyz,elmt,k,c)

# Compare with compressed system matrix
Acom = @time(getComprSystemMatrix(xyz,elmt,k,cuto))
ccom = Acom\b
f = open(string("windowck", k, "sph", nrMesh), "w")
writedlm(f,[real(ccom) imag(ccom)])
close(f)
resid = A*ccom-b
residcom = Acom*ccom-b

print(resid[1:10], "=resid, residcom = ", residcom[1:10], '\n')
print("Potential at -1 0 0 is ", potential(xyz,elmt,k,[-1 0 0],ccom), "\n\n")
checkSolution(xyz,elmt,k,ccom)

max_stars_repo_name: popsomer/asyBEM

max_issues_repo_name: popsomer/asyBEM

max_forks_repo_name: popsomer/asyBEM


text:
## Docstring Templates

using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) = """
                                         $(DOCSTRING)

                                         ---
                                         # Signatures
                                         $(TYPEDSIGNATURES)
                                         ---
                                         ## Methods
                                         $(METHODLIST)
                                         """

@template (TYPES) = """
                    $(TYPEDEF)
                    $(DOCSTRING)

                    ---
                    ## Fields
                    $(TYPEDFIELDS)
                    """

@template MODULES = """
                    $(DOCSTRING)

                    ---
                    ## Imports
                    $(IMPORTS)
                    """


text:
mutable struct Run{F <: AbstractFloat, I <: Integer, C <: Configuration}
    cost::Function
    cost_arguments::Dict{Symbol, Any}
    cost_evaluations::I
    cost_values::Array{F, 1}
    starting_point::C
    starting_cost::F
    report_after::I
    reporting_criterion::Function
    measurement_method::Function
    stopping_criterion::Function
    duration::I
    methods::Array{Any, 2}
    channel::RemoteChannel
end

function Run(;cost::Function                    = (c(x) = 0),
              cost_arguments::Dict{Symbol, Any} = Dict{Symbol, Any}(),
              cost_evaluations::I               = 1,
              cost_values::Array{F, 1}          = [0.0],
              starting_point::C                 = Configuration("empty"),
              starting_cost::F                  = 0.0,
              report_after::I                   = 100,
              reporting_criterion::Function     = iterations_reporting_criterion,
              measurement_method::Function      = measure_mean!,
              stopping_criterion::Function      = iterations_criterion,
              duration::I                       = 1_000,
              methods::Array{Any, 2}            = [[:simulated_annealing 2];],
              channel::RemoteChannel            = RemoteChannel(()->Channel{Any}(128))
             ) where {F <: AbstractFloat, I <: Integer, C <: Configuration}

    Run{F, I, C}(cost,
                 cost_arguments,
                 cost_evaluations,
                 cost_values,
                 starting_point,
                 starting_cost,
                 report_after,
                 reporting_criterion,
                 measurement_method,
                 stopping_criterion,
                 duration,
                 methods,
                 channel)
end


text:
module TestLeafVararg

using Comonicon.AST
using Comonicon.JuliaExpr
using Comonicon.JuliaExpr: emit, emit_body, emit_norm_body, emit_dash_body
using Test

const test_args = Ref{Vector{Any}}()

function foo(a, b = 2, c...)
    test_args[] = [a, b, c...]
end

cmd = Entry(;
    version = v"1.1.0",
    root = LeafCommand(;
        fn = foo,
        name = "leaf",
        args = [
            Argument(; name = "a", type = Int),
            Argument(; name = "b", type = Int, require = false),
        ],
        vararg = Argument(; name = "c", type = Float64, vararg = true),
    ),
)

eval(emit(cmd))

@testset "test leaf optional argument" begin
    @test command_main(["3", "2", "5.5", "6", "7"]) == 0
    @test test_args[] == Any[3, 2, 5.5, 6.0, 7.0]
    @test command_main(["3", "2"]) == 0
    @test test_args[] == Any[3, 2]
end

end


text:
# This file was generated, do not modify it. # hide
#hideall
colors=(:pink, :lightpink, :hotpink, :deeppink, :mediumvioletred, :palevioletred, :coral, :tomato, :orangered, :darkorange, :orange, :gold, :yellow)
print("@@ccols ")
print.("\\circle{$c}" for c in colors)
println("@@")

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
#module KERNEL
    include("atom.jl")
    include("atomContainer.jl")
    include("bond.jl")
    include("chain.jl")
    include("residue.jl")
    include("molecule.jl")
    include("PTE.jl")
    include("secondaryStructure.jl")
    include("system.jl")

    const CompositeType = Union{Atom,AtomContainer,Bond ,Nothing, Residue, Chain}

    #export AtomContainerMod.AtomContainer#, Atom, AtomContainer, Bond, Chain, Fragment, Molecule,
    #Element, SecondaryStructure, System

#export System, SecondaryStructure, Symbol, Element, Molecule,
#Fragment, Chain, Bond, AtomContainer, Atom

#end
#push!(LOAD_PATH,joinpath(pwd(),"src","KERNEL","common.jl"))
#=
common:
- Julia version:
- Author: Dan
- Date: 2021-06-07
=#


text:
module BlueMesh7

using Agents
using LightGraphs
using GraphPlot: gplot
using Random
using StatsBase: sample, Weights

using DataStructures: CircularBuffer
using Distributions
using LoopVectorization
using SpecialFunctions

include("actors.jl")
include("model.jl")
include("pomdp.jl")

end


text:
# vandermonde.jl tests

using SpecialMatrices: Vandermonde
import SpecialMatrices # dvand!, pvand!
using Test: @test, @testset, @test_throws, @inferred

# abstract vector construction
V = @inferred Vandermonde(5:7)
@test V == (5:7) .^ (0:2)'

a = [1,2,3+1im,8,5]
V = @inferred Vandermonde(a)
M = @inferred Matrix(V)

# Test element
@testset "basics" begin
    @test V[4,5] == a[4]^4

    @test isassigned(V, 2)
    @test !isassigned(V, 26)
    @test size(V) == (5, 5)
end

# Test full matrix conversions
@testset "convert" begin
    @test Matrix(V') == M'
    @test Matrix(transpose(V)) == transpose(M)
end

@testset "solve" begin
    # Test solving with vector and matrix rhs
    y = [1im,1,5,0,2]
    Y = [y 2*y]
    for rhs in [y, Y]
        # Test solution
        @test isapprox(V \ rhs, M \ rhs)
        @test isapprox(V' \ rhs, M' \ rhs)
        @test isapprox(rhs' / V, rhs' / M)
        @test isapprox(transpose(V) \ rhs, transpose(M) \ rhs)

        # Check that overloading works
        x = zero(M \ rhs)
        copyto!(x, rhs)
        SpecialMatrices.dvand!(a, x)
#       @which V \ rhs
        @test V \ rhs == x

        copyto!(x, rhs)
        SpecialMatrices.pvand!(a', x)
#       @which V' \ rhs
        @test V' \ rhs == x

#       @which rhs' / V
        @test rhs' / V == x'

        copyto!(x, rhs)
        SpecialMatrices.pvand!(a, x)
#       @which transpose(V) \ rhs
        @test transpose(V) \ rhs == x
    end
end

# Test dimension errors
@testset "dims" begin
    rhs = zeros(2,2)
    @test_throws DimensionMismatch V \ rhs
    @test_throws DimensionMismatch V' \ rhs
end


text:
#############################
##  Numerical experiments associated with section 5.1.1
#############################

using ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, SpecialFunctions, Plots
using SparseVolterraExamples

#############################
## Problem in Equation (11)
#############################
####
## This function sets up g in the correct basis as well as kernel K.
function solveSec511A(n)
    gf(x) = 1
    gF = Fun(x->gf(x),Jacobi(2,3, 0..1))
    Kfun(x,y) = x-y
        V = triVolterraFullKernelOpP01(Kfun,n,true)
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V
        V = Conversion(Jacobi(0,0,0..1),Jacobi(2,3,0..1))[1:n,1:n]*V[1:n,1:n]
        V = Derivative(Jacobi(0,1,0..1),2)[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);DirectEvalLHSP10atPrime0(n)*Derivative(Jacobi(0,1,0..1),1)[1:n,1:n];V[1:n,1:n]] \ [1;0;pad(gF.coefficients,n)]
end

####
## Compute the solution
coeff = solveSec511A(50)
u = Fun(Jacobi(0,1,0..1), coeff)
####
## Plot the computed and analytic solution
plot(u,grid=false,ylabel="u(x)",xlabel="x",label="sparse method",legend=:topleft)
plot!(x->cosh(x),0,1,grid=false,ylabel="u(x)",xlabel="x",label="analytic cosh(x)",legend=:topleft)
####
## Plot the error of the computed compared to analytic solution.
plot(x->abs(cosh(x)-u(x)),0,1,grid=false,ylabel="error",xlabel="x")
####
## Looping over the the degree of polynomials gives a convergence plot like in Figure 1.
function convSec511A(max,stepsize)
    errorvec = []
    for n=1:stepsize:max
        diff = maximum(abs.([maximum(Fun(Jacobi(0,1,0..1),solveSec511A(n))-Fun(x->cosh(x),0..1)),minimum(Fun(Jacobi(0,1,0..1),solveSec511A(n))-Fun(x->cosh(x),0..1))]))
        errorvec = append!(errorvec,diff)
    end
    return errorvec
end
plot(Array(1:1:14),convSec511A(14,1),grid=false,ylabel="error",xlabel="n",yaxis=:log10,label=false)

#############################
## Problem in Equation (12)
#############################
####
## This function sets up g in the correct basis as well as kernel K.
function solveSec511B(n)
    gf(x) = -1+x
    gF = Fun(x->gf(x),Jacobi(4,5, 0..1))
    Kfun(x,y) = y-x
    V = triVolterraFullKernelOpP01(Kfun,n,true)
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V
        V = Conversion(Jacobi(0,0,0..1),Jacobi(4,5,0..1))[1:n,1:n]*V[1:n,1:n]
        V = Derivative(Jacobi(0,1,0..1),4)[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);DirectEvalLHSP10atPrime0(n)*Derivative(Jacobi(0,1,0..1),1)[1:n,1:n];DirectEvalLHSP10atPPrime0(n)*Derivative(Jacobi(0,1,0..1),2)[1:n,1:n];DirectEvalLHSP10atPPPrime0(n)*Derivative(Jacobi(0,1,0..1),3)[1:n,1:n];V[1:n,1:n]] \ [-1;1;1;-1;pad(gF.coefficients,n)]
end
####
## Compute the solution
coeff = solveSec511B(15)
u = Fun(Jacobi(0,1,0..1), coeff)
####
## Plot the computed and analytic solution
plot(u,grid=false,ylabel="u(x)",xlabel="x",label="sparse method",legend=:topleft)
plot!(x->sin(x)-cos(x),0,1,grid=false,ylabel="u(x)",xlabel="x",label="analytic sin(x)-cos(x)",legend=:topleft)
####
## Plot the error of the computed compared to analytic solution.
plot(x->abs(sin(x)-cos(x)-u(x)),0,1,grid=false,ylabel="error",xlabel="x")
####
## Looping over the the degree of polynomials gives a convergence plot like in Figure 1.
function convSec511B(max,stepsize)
    errorvec = []
    for n=1:stepsize:max
        diff = maximum(abs.([maximum(Fun(Jacobi(0,1,0..1),solveSec511B(n))-Fun(x->sin(x)-cos(x),0..1)),minimum(Fun(Jacobi(0,1,0..1),solveSec511B(n))-Fun(x->sin(x)-cos(x),0..1))]))
        errorvec = append!(errorvec,diff)
    end
    return errorvec
end
plot(Array(1:4:40),convSec511B(40,4),grid=false,ylabel="error",xlabel="n",yaxis=:log10,label=false)

#############################
## Problem in Equation (13)
#############################
####
## This function sets up g in the correct basis as well as kernel K.
function solveSec511C(n)
    gf(x) = 1+x+(x^2)/2-(x^4)/(4*3*2)
    gF = Fun(x->gf(x),Jacobi(3,4, 0..1))
    Kfun(x,y) = (x-y)^2/2
    V = triVolterraFullKernelOpP01(Kfun,n,true)
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V
        V = Conversion(Jacobi(0,0,0..1),Jacobi(3,4,0..1))[1:n,1:n]*V[1:n,1:n]
        V = Derivative(Jacobi(0,1,0..1),3)[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);DirectEvalLHSP10atPrime0(n)*Derivative(Jacobi(0,1,0..1),1)[1:n,1:n];DirectEvalLHSP10atPPrime0(n)*Derivative(Jacobi(0,1,0..1),2)[1:n,1:n];V[1:n,1:n]] \ [1;2;1;pad(gF.coefficients,n)]
end
####
## Compute the solution
coeff = solveSec511C(15)
u = Fun(Jacobi(0,1,0..1), coeff)
####
## Plot the computed and analytic solution
plot(u,grid=false,ylabel="u(x)",xlabel="x",label="sparse method",legend=:topleft)
plot!(x->x+exp(x),0,1,grid=false,ylabel="u(x)",xlabel="x",label="analytic x+exp(x)",legend=:topleft)
####
## Plot the error of the computed compared to analytic solution.
plot(x->abs(x+exp(x)-u(x)),0,1,grid=false,ylabel="error",xlabel="x")
####
## Looping over the the degree of polynomials gives a convergence plot like in Figure 1.
function convSec511C(max,stepsize)
    errorvec = []
    for n=1:stepsize:max
        diff = maximum(abs.([maximum(Fun(Jacobi(0,1,0..1),solveSec511C(n))-Fun(x->x+exp(x),0..1)),minimum(Fun(Jacobi(0,1,0..1),solveSec511C(n))-Fun(x->x+exp(x),0..1))]))
        errorvec = append!(errorvec,diff)
    end
    return errorvec
end
plot(Array(1:1:20),convSec511C(20,1),grid=false,ylabel="error",xlabel="n",yaxis=:log10,label=false)

max_stars_repo_path: examples/Section 5-1-1 - Basic convolution kernels, second kind VIDEs.jl


text:
"""
    AbstractLattice

Contains information about a crystal lattice.
"""
abstract type AbstractLattice
end

"""
    RealLattice{N}

Describes an N-dimensional real space crystal lattice, with an associated primitive and
conventional lattice. These lattices may be the same depending on the lattice type.
"""
struct RealLattice{N} <: AbstractLattice
    # Primitive and conventional lattice vectors are stored together
    prim::SVector{N,SVector{N,Float64}}
    conv::SVector{N,SVector{N,Float64}}
end

# TODO: Implement lattice checking and conversion
# Some of this functionality might already be present in Crystalline.jl
# However, we need to be able to recognize primitive lattices with their own unique conditions

"""
    nuniqlen(basis::AbstractVector{<:AbstractVector{<:Real}}) -> Int

Returns the number of vectors that are the same length.
"""
function nuniqlen(basis::AbstractVector{<:AbstractVector{<:Real}})
    # unique!() does better than unique() here (per @btime)
    return length(unique!(norm.(basis)))
end

# Same thing for matrices
nuniqlen(basis::AbstractMatrix{<:Real}) = nuniqlen(_tovectors(basis))

"""
    northog(basis::AbstractVector{<:AbstractVector{<:Real}}) -> Int

Returns the number of pairs of orthogonal vectors.
"""
function northog(basis::AbstractVector{<:AbstractVector{<:Real}})
    # FIXME: This doesn't work for arbitrary dimensions
    # It needs to use the upper portion of the dot product matrix
    return count(iszero, dot.(basis, basis[circshift(1:length(basis), 1)]))
end

northog(basis::AbstractMatrix{<:Real}) = northog(_tovectors(basis))

# TODO: Write a number of tests for checking if lattices meet certain criteria

#=
Criteria for primitive unit cells in 3 dimensions:
    aP: 6 parameters (a, b, c, α, β, γ)
    mP: 4 parameters (a, b, c, β), α = γ = 90°
    mS: 4 parameters (a, c, α, γ), a = b, α = β
    oP: 3 parameters (a, b, c), α = β = γ = 90°
    oS: 3 parameters (a, c, γ), a = b, α = β = 90°
    oI: 3 parameters (a, α, γ), a = b = c, |cos(α)| + |cos(β)| + |cos(γ)| = 1
    oF: 3 parameters (a, b, c), sin(α)/a = sin(β)/b = sin(γ)/c, α + β + γ = 180°
        3 parameters (a, c, γ), [forgot the rest of this alternate test]
    tP: 2 parameters (a, c), a = b, α = β = γ = 90°
    tI: 2 parameters (a, α), a = b = c, |cos(α)| + |cos(β)| + |cos(γ)| = 1
    hP: 2 parameters (a, c), a = b, α = β = 90°, γ = 120°
    hR: 2 parameters (a, α), a = b = c, α = β = γ
    cP: 1 parameter (a), a = b = c, α = β = γ = 90°
    cI: 1 parameter (a), a = b = c, α = β = γ = arccos(-1/3)°
    cF: 1 parameter (a), a = b = c, α = β = γ = 60°

Rules that need to be checked:
    Number of unique lengths
    Number of orthogonal vectors
    Number of unique angles
    Cosine rule: |cos(α)| + |cos(β)| + |cos(γ)| = 1
    Sine rule: sin(α)/a = sin(β)/b = sin(γ)/c, α + β + γ = 180°

    The last three are probably going to be trickier to implement, or at the very least there
    might be some unintuitive linear algebra tricks to make those checks happen faster.
=#

# TODO: what are the crtieria for a "good" lattice basis?
# Potentially worth looking at the Niggli reduction algorithm to do this

# TODO: Tools for working with reciprocal space lattices


text:
# a simple CUDA host array type for testing purposes.
# use CuArrays.jl for real applications!

mutable struct CuHostArray{T,N}
    ptr::CuPtr{T}
    dims::Dims{N}
end

Base.pointer(xs::CuHostArray) = xs.ptr

Base.size(xs::CuHostArray) = xs.dims

Base.length(xs::CuHostArray) = prod(xs.dims)

Base.unsafe_convert(::Type{<:CuPtr}, x::CuHostArray) = pointer(x)

## ctors & finalizer

function CuHostArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
    len = prod(dims)
    sz = len*sizeof(T)
    buf = try
        Mem.alloc(Mem.Device, sz)
    catch err
        (isa(err, CuError) && err.code == CUDAdrv.ERROR_OUT_OF_MEMORY) || rethrow()
        GC.gc(true)
        Mem.alloc(Mem.Device, sz)
    end
    ptr = convert(CuPtr{T}, buf)

    obj = CuHostArray{T,N}(ptr, dims)
    finalizer(obj) do a
        CUDAdrv.isvalid(buf.ctx) && Mem.free(buf)
    end
    return obj
end

# underspecified
CuHostArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = CuHostArray{T,N}(undef, dims)
CuHostArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = CuHostArray{T,N}(undef, dims)
CuHostArray{T}(::UndefInitializer, dims::Integer...) where {T} =
  CuHostArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

Base.similar(a::CuHostArray{T,N}) where {T,N} = CuHostArray{T,N}(undef, size(a))

## memory copy operations

function CuHostArray(src::Array{T,N}) where {T,N}
    dst = CuHostArray{T,N}(undef, size(src))
    unsafe_copyto!(pointer(dst), pointer(src), length(src))
    return dst
end

function Base.Array(src::CuHostArray{T,N}) where {T,N}
    dst = Array{T,N}(undef, size(src))
    unsafe_copyto!(pointer(dst), pointer(src), length(src))
    return dst
end

## conversions

using Adapt
function Adapt.adapt_storage(::CUDAnative.Adaptor, a::CuHostArray{T,N}) where {T,N}
    ptr = pointer(a)
    devptr = CUDAnative.DevicePtr{T,AS.Global}(ptr)
    CuDeviceArray{T,N,AS.Global}(a.dims, devptr)
end


text:
"""
    A = JopGradient(dom, δ)

Gradient for a 2D or 3D arrays. The range will have one more
dimension than the domain and contain the components of the gradient
in each dimension.

* `dom::JetSpace{T,N}` is the domain of the operator.
* `δ::NTuple{T,N}` is the grid spacing in each dimension.

# Examples:

## 2D
```julia
nz,nx = 11,12
dom = JetSpace(Float32, nz, nx)
A = JopGradient(dom, (1.0,1.0))
size(domain(A)) # (11,12)
size(range(A))  # (11,12,2)
```

## 3D
```julia
nz,ny,nx = 11,12,13
dom = JetSpace(Float32, nz, ny, nx)
A = JopGradient(dom,(1.0,1.0,1.0))
size(domain(A)) # (11,12,13)
size(range(A))  # (11,12,13,3)
```
"""
function JopGradient(dom::JetSpace{T,N}, δ::NTuple{N,T}) where {T,N}
    rng = JetSpace(T, (size(dom)...,length(size(dom))))
    if length(size(dom)) != length(δ)
        error("dimensionality of length(size(dom)) ($(length(size(dom)))) and length(δ) ($(length(δ))) do not agree!")
    end
    JopLn(dom = dom, rng = rng, df! = JopGradient_df!, df′! = JopGradient_df′!, s = (δ=δ,))
end

export JopGradient

# ....................................................................................
# 2D
# ....................................................................................
function JopGradient_df!(d::AbstractArray{T,3}, m::AbstractArray{T,2}; δ, kwargs...) where {T}
    d .= 0
    n1,n2 = size(m)
    for k1=1:n1-1, k2=1:n2
        d[k1,k2,1] += m[k1+1,k2] / δ[1]
        d[k1,k2,1] -= m[k1+0,k2] / δ[1]
    end
    for k1=1:n1, k2=1:n2-1
        d[k1,k2,2] += m[k1,k2+1] / δ[2]
        d[k1,k2,2] -= m[k1,k2+0] / δ[2]
    end
    d
end

function JopGradient_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,3}; δ, kwargs...) where {T}
    m .= 0
    n1,n2 = size(m)
    for k1=1:n1-1, k2=1:n2
        m[k1+1,k2] += d[k1,k2,1] / δ[1]
        m[k1+0,k2] -= d[k1,k2,1] / δ[1]
    end
    for k1=1:n1, k2=1:n2-1
        m[k1,k2+1] += d[k1,k2,2] / δ[2]
        m[k1,k2+0] -= d[k1,k2,2] / δ[2]
    end
    m
end

# ....................................................................................
# 3D
# ....................................................................................
function JopGradient_df!(d::AbstractArray{T,4}, m::AbstractArray{T,3}; δ, kwargs...) where {T}
    d .= 0
    n1,n2,n3 = size(d)
    for k1=1:n1-1, k2=1:n2, k3=1:n3
        d[k1,k2,k3,1] += m[k1+1,k2,k3] / δ[1]
        d[k1,k2,k3,1] -= m[k1+0,k2,k3] / δ[1]
    end
    for k1=1:n1, k2=1:n2-1, k3=1:n3
        d[k1,k2,k3,2] += m[k1,k2+1,k3] / δ[2]
        d[k1,k2,k3,2] -= m[k1,k2+0,k3] / δ[2]
    end
    for k1=1:n1, k2=1:n2, k3=1:n3-1
        d[k1,k2,k3,3] += m[k1,k2,k3+1] / δ[3]
        d[k1,k2,k3,3] -= m[k1,k2,k3+0] / δ[3]
    end
    d
end

function JopGradient_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,4}; δ, kwargs...) where {T}
    m .= 0
    n1,n2,n3 = size(m)
    for k1=1:n1-1, k2=1:n2, k3=1:n3
        m[k1+1,k2,k3] += d[k1,k2,k3,1] / δ[1]
        m[k1+0,k2,k3] -= d[k1,k2,k3,1] / δ[1]
    end
    for k1=1:n1, k2=1:n2-1, k3=1:n3
        m[k1,k2+1,k3] += d[k1,k2,k3,2] / δ[2]
        m[k1,k2+0,k3] -= d[k1,k2,k3,2] / δ[2]
    end
    for k1=1:n1, k2=1:n2, k3=1:n3-1
        m[k1,k2,k3+1] += d[k1,k2,k3,3] / δ[3]
        m[k1,k2,k3+0] -= d[k1,k2,k3,3] / δ[3]
    end
    m
end


text:
#utilities for QUtilities
#some of these are candidates for Finometrics
const CRSP_MISSINGS = [
  "A", "B", "C", "D", "E", "S", "T", "P",
  -55.0, -66.0, -77.0, -88.0, -99.0,
  -55, -66, -77, -88, -99,
  "-55.0", "-66.0", "-77.0", "-88.0", "-99.0"]

const CRSP_PARSED_MISSINGS = (-55.0, -66.0, -77.0, -88.0, -99.0)

const USED_COMPA_VALUE_FIELDS = [:at]
const USED_COMPQ_VALUE_FIELDS = [:atq, :ltq, :ibq, :txditcq, :pstkrq, :seqq, :ceqq, :pstkq]

const RETAINED_COLS_COMP_Q = [:gvkeyq, :dateq, :fyearq]

#NOTE: uncomment when we have the missings const COMP_PARSED_MISSINGS = (nothing)

#lags a variable in the time series
#=function lagwithin!(df::DataFrame, targets::Vector{Symbol};
  period::Symbol = :fisyr,
  lags::Int=1,
  laggednames::Vector{Symbol} = (lags ≠ 1 ?
    (s::Symbol->Symbol(:L, lags, s)).(targets) : (s::Symbol->Symbol(:L, s)).(targets)),
  group=:gvkey,
  sorted::Bool=false)::Nothing
end=#

#WARNING quick and dirty hack, could create issues if multiple versions of Int8
const QQuarter = Int8
const QYear = Int16

#converts gzip to lz4
function gz2lz4OLD(infile::String; outfile::String = infile[1:(length(infile)-4)] * lowercase(
  replace(uppercase(infile[(length(infile)-3):end]),".GZ"=>".LZ4")),
  deleteold::Bool = false)

  @assert isfile(infile)

  #Below doesnt work
  indata = open(infile, "r")
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4CompressorStream(outdata))
  write(s, indata)
  close(s)

  #run(``)

end

#converts gzip to lz4
function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  @assert isfile(infile)

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4CompressorStream(outdata))
  write(s, indata)
  close(s)

  #run(``)

end

FILES_TO_RECOMPRESS = ["COMP-A", "COMP-Q", "CCM", "CRSP-M", "CRSP-D"]
function recompressqfactorfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2lz4(f, location=datapath)
  end

  return nothing
end

recompressqfactorfile(file; datapath=DATA_PATH)::Nothing = recompressqfactorfiles(
  [file]; datapath=datapath)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
############################################################################################
"""
$(TYPEDEF)

SMC tuning container.

# Fields
$(TYPEDFIELDS)
"""
struct SMCTune{
    T<:Tagged,
    F<:Function,
    B<:BaytesCore.UpdateBool,
    R<:BaytesCore.ResamplingMethod,
    U<:BaytesCore.UpdateBool,
    S<:BaytesCore.UpdateBool
    }
    "Tagged Model parameter."
    tagged::T
    "Iterator that contains current data iteration ~ starts with Ndata."
    iter::Iterator
    "ParticleTune struct to possibly adaptively change number of chains."
    chains::BaytesCore.ChainsTune
    "Resampling Method of parameter particle indices."
    resample::BaytesCore.ResampleTune{R}
    "Function that is applied to correlation of parameter to determine if jittering can be stopped."
    jitter::BaytesCore.JitterTune{B}
    "Determines number of jittering steps in jitter kernel."
    jitterfun::F
    "UpdateBool if parameter can be captured after initial jitter step."
    capture::U
    "Number of warmup Tuning steps for kernels."
    Ntuning::Int64
    "Batchsize for parallel processing of jittering steps."
    batchsize::Int64
    "Boolean if generated quantities should be generated while sampling"
    generated::S
    function SMCTune(
        tagged::T,
        Nchains::Int64,
        Ndata::Int64,
        Ntuning::Int64,
        resamplingthreshold::Float64,
        resampling::D,
        capture::U,
        jitterfun::F,
        jitteradaption::B,
        jitterthreshold::Float64,
        Njitter_min::Integer,
        Njitter_max::Integer,
        generated::S
    ) where {
        T<:Tagged,
        F<:Function,
        B<:BaytesCore.UpdateBool,
        D<:BaytesCore.ResamplingMethod,
        U<:BaytesCore.UpdateBool,
        S<:BaytesCore.UpdateBool
    }
        ## Check input for correctness
        ArgCheck.@argcheck 0 < Ntuning "Number of tuning iterations has to be positive"
        ArgCheck.@argcheck 0.0 <= resamplingthreshold "threshold has to be positive"
        ArgCheck.@argcheck 0.0 <= jitterthreshold "jitterthreshold has to be positive"
        ArgCheck.@argcheck 0 < Njitter_min <= Njitter_max
        ## Assign tuning structs
        jitter_maybe = Updater(false)
        #!NOTE: Start with Ndata, so first propose!() step has correct iteration after training.
        iter = Iterator(Ndata)
        particletune = BaytesCore.ChainsTune(
            Nchains / Ndata, resamplingthreshold, Nchains, Ndata
        )
        jitter = BaytesCore.JitterTune(jitteradaption, jitterthreshold, Njitter_min, Njitter_max)
        resample = BaytesCore.ResampleTune(resampling, BaytesCore.Updater(false))
        ## Compute batchsize for parallel processing in resampling step
        batchsize = compute_batchsize(Nchains)
        ## Return SMC tune
        return new{T,F,B,D,U,S}(
            tagged, iter, particletune, resample, jitter, jitterfun, capture, Ntuning, batchsize, generated
        )
    end
end

############################################################################################
function update!(tune::SMCTune)
    ##  Update iteration counter
    update!(tune.iter)
    ## Set jitter count back to 0
    update!(tune.jitter)
    return nothing
end

############################################################################################
#export
export SMCTune


text:
using LinearAlgebra

function recursion_partition_function(num_elec, eff_spec_SpinUp, eff_spec_SpinDn)
    Z_up, Z_dn = zeros(num_elec), zeros(num_elec)
    zₖ_up, zₖ_dn = zeros(num_elec), zeros(num_elec)
    for k = 1 : num_elec
        zₖ_up[k] = sum(eff_spec_SpinUp .^ k)
        zₖ_dn[k] = sum(eff_spec_SpinDn .^ k)
    end
    for N = 1 : num_elec
        for k = 1 : N
            sgn_k = 2 * isodd(k) - 1
            if k == N
                Z_up[N] += sgn_k * zₖ_up[k]
                Z_dn[N] += sgn_k * zₖ_dn[k]
            else
                Z_up[N] += sgn_k * Z_up[N-k] * zₖ_up[k]
                Z_dn[N] += sgn_k * Z_dn[N-k] * zₖ_dn[k]
            end
        end
        Z_up[N] = Z_up[N] / N
        Z_dn[N] = Z_dn[N] / N
    end
    return Z_up, Z_dn
end

function recursion_OccNum(Z_up, eff_spec_SpinUp, Z_dn, eff_spec_SpinDn, num_elec, num_sites) # Calculating <nᵢ> * Z_N
    n_up_mean = zeros(num_elec, num_sites)
    n_dn_mean = zeros(num_elec, num_sites)
    n_up_mean[1, :], n_dn_mean[1, :] = eff_spec_SpinUp, eff_spec_SpinDn
    for N = 2 : num_elec
        n_up_mean[N, :] = eff_spec_SpinUp .* (Z_up[N-1] .- n_up_mean[N-1, :])
        n_dn_mean[N, :] = eff_spec_SpinDn .* (Z_dn[N-1] .- n_dn_mean[N-1, :])
    end
    return n_up_mean ./ Z_up, n_dn_mean ./ Z_dn
end

function measure_tot_energy(num_sites, t, U, U_array_SpinUp, U_array_SpinDn, n_up_mean, n_dn_mean, num_SpinUp, num_SpinDn, Z_up, Z_dn)
    E_onebody, E_twobody = 0, 0
    for i = 1 : num_sites
        for j = 1 : num_sites
            if abs(i - j) == 1 || abs(i - j) == num_sites - 1
                D_ij_up = sum(U_array_SpinUp[i, j, :] .* n_up_mean[num_SpinUp, :])
                D_ij_dn = sum(U_array_SpinDn[i, j, :] .* n_dn_mean[num_SpinDn, :])
                E_onebody += -t * ( D_ij_up + D_ij_dn )
            end
        end
    end
    for i = 1 : num_sites
        U_matrix = U_array_SpinUp[i, i, :] * U_array_SpinDn[i, i, :]'
        n_matrix = n_up_mean[num_SpinUp, :] * n_dn_mean[num_SpinDn, :]'
        D_iiii = sum(U_matrix .* n_matrix)
        E_twobody += U * D_iiii
    end
    return E_onebody, E_twobody
end

function recursion(self::Lattice, expϵup::Array{Float64,1}, expϵdn::Array{Float64,1})
    nup_mean = zeros(Float64, self.nsites, self.nsites)
    ndn_mean = zeros(Float64, self.nsites, self.nsites)
    Znup_ratio = zeros(Float64, self.ntot)
    Zndn_ratio = zeros(Float64, self.ntot)
    expϵup = expϵup ./ expϵup[self.nsites]
    expϵdn = expϵdn ./ expϵdn[self.nsites]
    Znup_ratio[1], Zndn_ratio[1] = sum(expϵup), sum(expϵdn)
    nup_mean[:, 1], ndn_mean[:, 1] = expϵup / Znup_ratio[1], expϵdn / Zndn_ratio[1]
    for N = 2 : max(self.nup, self.ndn)
        Znup_ratio[N] = sum(expϵup .* (1 .- nup_mean[:, N-1])) / N
        Zndn_ratio[N] = sum(expϵdn .* (1 .- ndn_mean[:, N-1])) / N
        nup_mean[:, N] = expϵup .* (1 .- nup_mean[:, N - 1]) / Znup_ratio[N]
        ndn_mean[:, N] = expϵdn .* (1 .- ndn_mean[:, N - 1]) / Zndn_ratio[N]
    end
    return nup_mean, ndn_mean, Znup_ratio, Zndn_ratio
end

function measure_energy(
    self::Lattice, Uup_array::Array{Float64,3}, Udn_array::Array{Float64,3},
    nup_mean::Array{Float64,2}, ndn_mean::Array{Float64,2}
)
    Eonebody, Etwobody = 0, 0
    for i = 1 : self.nsites
        for j = 1 : self.nsites
            if abs(i - j) == 1 || abs(i - j) == self.nsites - 1
                Dij_up = sum(Uup_array[:, i, j] .* nup_mean[:, self.nup])
                Dij_dn = sum(Udn_array[:, i, j] .* ndn_mean[:, self.ndn])
                Eonebody += -self.t * ( Dij_up + Dij_dn )
            end
        end
    end
    for i = 1 : self.nsites
        Umatrix = Uup_array[:, i, i] * Udn_array[:, i, i]'
        nmatrix = nup_mean[:, self.nup] * ndn_mean[:, self.ndn]'
        Diiii = sum(Umatrix .* nmatrix)
        Etwobody += self.U * Diiii
    end
    return Eonebody, Etwobody
end

# rank-one perturbation calculator
function rankone_perturb!(
    γ::Float64, σil::Int64, Ns::Int64,
    MatProdUp_new::Array{Float64,2}, MatProdDn_new::Array{Float64,2}
)
    Δup, Δdn = Matrix{Float64}(I, size(MatProdUp_new)), Matrix{Float64}(I, size(MatProdDn_new))
    Δup[Ns, Ns] = exp(-4 * γ * σil)
    Δdn[Ns, Ns] = exp(4 * γ * σil)
    MatProdUp_new = Δup * MatProdUp_new
    MatProdDn_new = Δdn * MatProdDn_new
    return MatProdUp_new, MatProdDn_new
end

# determine the sign of the accept_ratio
@inline function Rsign(R::Float64)
    Rsgn = 2 * (R > 0) - 1
end

function MCIntegration(lattice::Lattice, qmc::QMC, opt::Opt)
    # initialize vectors that store MC data
    Ek_array = Vector{Float64}() # kinetic energy
    Ep_array = Vector{Float64}() # potential energy
    E_array = Vector{Float64}() # total energy
    sgn_array = Vector{Int64}() # sign
    NiMax_array = Vector{Float64}()
    # initialize a random configuration
    σlist = 2 * (rand(Float64, (lattice.nsites, qmc.L)) .< 0.5) .- 1
    MatProdUp, MatProdDn, EigensUp, EigensDn = matrix_product_QR(σlist, lattice, qmc, opt.stab_interval)
    nup_mean, ndn_mean, Znup_ratio, Zndn_ratio = recursion(lattice, real(EigensUp.values), real(EigensDn.values))
    sgn = Rsign(prod(Znup_ratio[1 : lattice.nup]) * prod(Zndn_ratio[1 : lattice.nup]))
    for n = 1 : qmc.samples
        σlist_new = copy(σlist)
        # random flip a spin
        #flipL, flipNs = rand(1 : qmc.L), rand(1 : lattice.nsites)
        flip = rand(1 : qmc.L * lattice.nsites)
        σlist_new[flip] *= -1
        #rankone_perturb!(qmc.γ, σlist[flipNs, flipL], flipNs, MatProdUp_new, MatProdDn_new)
        MatProdUp_new, MatProdDn_new, EigensUp_new, EigensDn_new = matrix_product_QR(σlist_new, lattice, qmc, opt.stab_interval)
        nup_mean_new, ndn_mean_new, Znup_ratio_new, Zndn_ratio_new = recursion(lattice, real(EigensUp_new.values), real(EigensDn_new.values))
        Rup = prod(Znup_ratio_new[1 : lattice.nup] ./ Znup_ratio[1 : lattice.nup])
        Rdn = prod(Zndn_ratio_new[1 : lattice.ndn] ./ Zndn_ratio[1 : lattice.ndn])
        sgn_new = Rsign(Rup * Rdn) * sgn
        R = abs(Rup * Rdn)
        accept_ratio = R / (1 + R)
        if rand() <= accept_ratio   # accept
            σlist = σlist_new
            Znup_ratio, Zndn_ratio = Znup_ratio_new, Zndn_ratio_new
            EigensUp, EigensDn = EigensUp_new, EigensDn_new
            nup_mean, ndn_mean = nup_mean_new, ndn_mean_new
            sgn = sgn_new
            Uup_array, Udn_array = overlap_coefficient(lattice.nsites, EigensUp.vectors, EigensDn.vectors)
            Eonebody, Etwobody = measure_energy(lattice, Uup_array, Udn_array, nup_mean, ndn_mean)
            push!(sgn_array, sgn)
            push!(Ek_array, Eonebody * sgn)
            push!(Ep_array, Etwobody * sgn)
            push!(E_array, Eonebody * sgn + Etwobody * sgn)
            push!(NiMax_array, maximum(nup_mean[:, lattice.nup]))
        else    # reject
            Uup_array, Udn_array = overlap_coefficient(lattice.nsites, EigensUp.vectors, EigensDn.vectors)
            Eonebody, Etwobody = measure_energy(lattice, Uup_array, Udn_array, nup_mean, ndn_mean)
            push!(sgn_array, sgn)
            push!(Ek_array, Eonebody * sgn)
            push!(Ep_array, Etwobody * sgn)
            push!(E_array, Eonebody * sgn + Etwobody * sgn)
            push!(NiMax_array, maximum(nup_mean[:, lattice.nup]))
        end
    end
    E_mean = mean(E_array[opt.nburn_in : qmc.samples])
    E_error = std(E_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
    Ek_mean = mean(Ek_array[opt.nburn_in : qmc.samples])
    Ek_error = std(Ek_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
    Ep_mean = mean(Ep_array[opt.nburn_in : qmc.samples])
    Ep_error = std(Ep_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
    sgn_avg = mean(sgn_array[opt.nburn_in : qmc.samples])
    NiMax_avg = mean(NiMax_array[opt.nburn_in : qmc.samples])
    return E_mean, E_error, Ek_mean, Ek_error, Ep_mean, Ep_error, sgn_avg, NiMax_avg
end

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
@testset "Z4Group" begin
    for (n, x) in zip(1:4, (im, -1, -im, 1))
        @test z4group(n) == Z4Group(x)
    end
end

@testset "Z4Group0" begin
    for x in (im, -1, -im, 1, 0)
        @test Z4Group0(x) == x
    end
    for (x, y, z) in ((0, 0, 0), (0, 1, 0), (0, -1, 0), (0, im, 0), (0, -im, 0),
                      (1, 0, 0), (1, 1, 1), (1, -1, -1), (1, im, im), (1, -im, -im),
                      (-1, 0, 0), (-1, 1, -1), (-1, -1, 1), (-1, im, -im), (-1, -im, im),
                      (im, 0, 0), (im, 1, im), (im, -1, -im), (im, im, -1), (im, -im, 1),
                      (-im, 0, 0), (-im, 1, -im), (-im, -1, im), (-im, im, 1), (-im, -im, -1))
        @test Z4Group0(x) * Z4Group0(y) == Z4Group0(z)
    end
end


text:
using AlgThys
using AlgThys.Examples
using Test

@testset "AlgThys" begin
  @test check_presentation(UnlabeledModule)

  @present Module begin
    A::SORT
    S::SORT
    m::((A,A) -> A)
    ap::((A,S) -> S)
    assoc::(ap(m(x,y),z) == ap(x,ap(y,z)) ⊣ (x::A,y::A,z::S))
  end

  @test check_presentation(unlabel(Module))
end


text:
## input the tolerance for mismatch
tol = 1e-2
function check_difference(a,b,tol)
	return abs((a-b)/(a+b+1e-7)) <= tol
end

## Parse the network data
network = _WM.parse_file("data/epanet/van_zyl.inp")
network_mn = _WM.make_multinetwork(network)

network_ids = sort([parse(Int, nw) for (nw, nw_data) in network_mn["nw"]])

# parse a feasible solution
result = JSON.parsefile("data/solution/van_zyl_PWLRD_feasible_sol.json"; dicttype=Dict,
                        inttype=Int64, use_mmap=true)

solution = result["solution"]
wntr_network = initialize_wntr_network(network_mn)
update_wntr_controls(wntr_network, network_mn, solution, float(network_mn["time_step"]))
wntr_result = simulate_wntr(wntr_network)


text:
using Plots
using Random
Random.seed!(9)

function square(c,side)
  x = [ c[1]-side/2, c[1]+side/2, c[1]+side/2, c[1]-side/2, c[1]-side/2]
  y = [ c[2]-side/2, c[2]-side/2, c[2]+side/2, c[2]+side/2, c[2]-side/2]
  return x, y
end

plot()

x,y=square([1,1],2)
plot!(x,y,seriestype=[:shape],
      linewidth=2,fillalpha=0.2,color="green",label="")

x,y=square([1,1],6)
plot!(x,y,seriestype=[:shape],
      linewidth=2,fillalpha=0.1,color="orange",label="")

lines = collect(-2:2:6)
vline!(lines,color="gray",label="")
hline!(lines,color="gray",label="")

px = [ 0.1 + 9.8*rand() for i in 1:100 ]
py = [ 0.1 + 9.8*rand() for i in 1:100 ]
scatter!(px,py,label="",alpha=0.2,color="blue")

for i in 1:5
  for j in 1:5
    annotate!(2*j-1,2*i-1,text("($i,$j)"))
  end
end

color=:blue
annotate!(-1,-1,text("(3,3)",color))
for i in 1:2:5
  annotate!(-1,i,text("($((i-1)÷2+1),3)",color))
  annotate!(i,-1,text("(3,$((i-1)÷2+1))",color))
end

plot!(size=(400,400),
      xlim=(-2,6),
      ylim=(-2,6),
      framestyle=:box,
      xlabel="x",ylabel="y",grid=false)

savefig("lcells3.pdf")


text:
using Flux
println("forward_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
println("Starting Soon!")
using SciMLBase
using DomainSets
import ModelingToolkit: Interval

@testset "ODE" begin
    @parameters x
    @variables u(..)

    Dx = Differential(x)
    eq = Dx(u(x)) ~ 0.
    bcs = [u(0.) ~ u(0.)]
    domains = [x ∈ Interval(0.0,1.0)]
    chain = FastChain((x,p) -> x.^2)

    chain([1],Float64[])
    strategy_ = NeuralPDE.GridTraining(0.1)
    discretization = NeuralPDE.PhysicsInformedNN(chain,strategy_;init_params = Float64[])
    @named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
    prob = NeuralPDE.discretize(pde_system,discretization)
    discretized_functions = NeuralPDE.discretize_inner_functions(pde_system,discretization)

    eqs = pde_system.eqs
    bcs = pde_system.bcs
    domains = pde_system.domain
    dx = strategy_.dx
    eltypeθ = eltype(discretized_functions.flat_initθ)
    depvars,indvars,dict_indvars,dict_depvars, dict_depvar_input = NeuralPDE.get_vars(pde_system.ivs,
                                                        pde_system.dvs)

    train_sets = generate_training_sets(domains,dx,eqs,bcs,eltypeθ,
                                        dict_indvars,dict_depvars)

    pde_train_sets, bcs_train_sets = train_sets
    pde_train_sets = NeuralPDE.adapt(eltypeθ,pde_train_sets)[1]

    train_data =pde_train_sets
    pde_loss_function = discretized_functions.inner_pde_loss_functions[1]

    dudx(x) = @. 2*x
    @test pde_loss_function(train_data, Float64[]) ≈ dudx(train_data) rtol = 1e-8
end

@testset "derivatives" begin
    chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
    initθ = Float64.(DiffEqFlux.initial_params(chain))

    eltypeθ = eltype(initθ)
    parameterless_type_θ = DiffEqBase.parameterless_type(initθ)
    phi = NeuralPDE.get_phi(chain,parameterless_type_θ)
    derivative = NeuralPDE.get_numeric_derivative()

    u_ = (cord, θ, phi)->sum(phi(cord, θ))

    phi([1,2], initθ)

    phi_ = (p) -> phi(p, initθ)[1]
    dphi = Zygote.gradient(phi_,[1.,2.])

    function get_ε(dim, der_num,eltypeθ)
        epsilon = cbrt(eps(eltypeθ))
        ε = zeros(eltypeθ, dim)
        ε[der_num] = epsilon
        ε
    end

    eps_x = get_ε(2, 1,Float64)
    eps_y = get_ε(2, 2,Float64)

    dphi_x = derivative(phi,u_,[1.,2.],[eps_x],1,initθ)
    dphi_y = derivative(phi,u_,[1.,2.],[eps_y],1,initθ)

    #first order derivatives
    @test isapprox(dphi[1][1], dphi_x, atol=1e-8)
    @test isapprox(dphi[1][2], dphi_y, atol=1e-8)

    dphi_x = derivative(phi,u_,[1.,2.],[[ 0.0049215667, 0.0]],1,initθ)
    dphi_y = derivative(phi,u_,[1.,2.],[[0.0,  0.0049215667]],1,initθ)

    hess_phi = Zygote.hessian(phi_,[1,2])

    dphi_xx = derivative(phi,u_,[1.,2.],[eps_x,eps_x],2,initθ)
    dphi_xy = derivative(phi,u_,[1.,2.],[eps_x,eps_y],2,initθ)
    dphi_yy = derivative(phi,u_,[1.,2.],[eps_y,eps_y],2,initθ)

    #second order derivatives
    @test isapprox(hess_phi[1], dphi_xx, atol=1e-5)
    @test isapprox(hess_phi[2], dphi_xy, atol=1e-5)
    @test isapprox(hess_phi[4], dphi_yy, atol=1e-5)
end

@testset "Integral" begin
    #semi-infinite intervals
    @parameters x
    @variables u(..)
    I = Integral(x in ClosedInterval(0, Inf))
    eq = I(u(x)) ~ 0
    bcs = [u(1.) ~ exp(1)/(exp(2) + 3)]
    domains = [x ∈ Interval(1.0, 2.0)]
    chain = FastChain((x,p) -> exp.(x) ./ (exp.(2 .*x) .+ 3))
    chain([1],Float64[])
    strategy_ = NeuralPDE.GridTraining(0.1)
    discretization = NeuralPDE.PhysicsInformedNN(chain,strategy_;init_params = Float64[])
    @named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
    sym_prob = SciMLBase.symbolic_discretize(pde_system, discretization)
    prob = NeuralPDE.discretize(pde_system,discretization)
    discretized_functions = NeuralPDE.discretize_inner_functions(pde_system,discretization)
    inner_loss =discretized_functions.inner_pde_loss_functions[1]
    exact_u = π/(3*sqrt(3))
    @test  inner_loss(ones(1,1), Float64[])[1] ≈ exact_u  rtol = 1e-5

    #infinite intervals
    @parameters x
    @variables u(..)
    I = Integral(x in ClosedInterval(-Inf, Inf))
    eqs = I(u(x)) ~ 0
    domains = [x ∈ Interval(1.0, 2.0)]
    bcs = [u(1) ~ u(1)]
    chain = FastChain((x,p) -> x .* exp.(-x .^2))
    chain([1],Float64[])

    discretization = NeuralPDE.PhysicsInformedNN(chain,strategy_;init_params = Float64[])
    @named pde_system = PDESystem(eqs, bcs, domains, [x], [u(x)])
    sym_prob = SciMLBase.symbolic_discretize(pde_system, discretization)
    prob = SciMLBase.discretize(pde_system, discretization)
    discretized_functions = NeuralPDE.discretize_inner_functions(pde_system,discretization)
    inner_loss =discretized_functions.inner_pde_loss_functions[1]
    exact_u = 0
    @test  inner_loss(ones(1,1), Float64[])[1] ≈ exact_u  rtol = 1e-9
end


text:
# Arithmetic Operators

x = 1;
y = 2;

+x;
-x;
x + y;
x - y;
x * y;
x / y; # Floating point division
x ÷ y; # Integer division
x \ y; # Inverse division
y ^ 2;
x % y;

# Bitwise Operators

b = Int8(10);
println(bitstring(b));

~b; # bitwise NOT
b & b; # bitwise and
b | b; # bitwise or
b ⊻ b; # bitwise xor
b >>> 2; # logical right shift
b >> 2; # bitwise right shift
b << 2; # bitwise left shift

# Logical Operators

a = true;
b = false;

a && b;
a || b;
!a;

# Assignment Operators

x = 1;
x += 1;
x -= 1;
x *= 1;
x /= 1;

# Vectorized/dot Operators
# - tells compiler to do operator to all elements in array

A = [1, 2, 3, 4];
println(A);
A = A .^ 2;
println(A);
A .^= 2;
println(A);

# Relational Operators

x = 2;
y = 3;

x == y;
x < y;
x > y;
x <= y;
x >= y;
x != y;

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
function FieldDenseOutputModel()
    return FieldDenseOutputModel(())
end

function append(obj::FieldDenseOutputModel, arg0::FieldDenseOutputModel)
    return jcall(obj, "append", void, (FieldDenseOutputModel,), arg0)
end

function get_final_time(obj::FieldDenseOutputModel)
    return jcall(obj, "getFinalTime", RealFieldElement, ())
end

function get_initial_time(obj::FieldDenseOutputModel)
    return jcall(obj, "getInitialTime", RealFieldElement, ())
end

function get_interpolated_state(obj::FieldDenseOutputModel, arg0::RealFieldElement)
    return jcall(obj, "getInterpolatedState", FieldODEStateAndDerivative, (RealFieldElement,), arg0)
end

function handle_step(obj::FieldDenseOutputModel, arg0::FieldODEStateInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (FieldODEStateInterpolator, jboolean), arg0, arg1)
end

function init(obj::FieldDenseOutputModel, arg0::FieldODEStateAndDerivative, arg1::RealFieldElement)
    return jcall(obj, "init", void, (FieldODEStateAndDerivative, RealFieldElement), arg0, arg1)
end


text:
type Counter
  X::Int64
end

GLOBAL_COUNTER = Counter(0)
inc(c::Counter) = c.X +=1
genint() = (inc(GLOBAL_COUNTER);GLOBAL_COUNTER.X-1)
genvar(prefix="x") = "$prefix$(genint())"
restart_counter!() = GLOBAL_COUNTER.X = 0