text:
module HeatEquationTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

import Gridap: ∇
import GridapODEs.TransientFETools: ∂t

θ = 0.0

# Analytical functions
# u(x,t) = (x[1]+x[2])*t
# u(x,t) = (2*x[1]+x[2])*t
u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
v(x) = t -> u(x,t)
∂tu(t) = x -> ForwardDiff.derivative(v(x),t)
∂tu(x,t) = ∂tu(t)(x)
∂t(::typeof(u)) = ∂tu
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

V0 = FESpace(
  reffe=lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="boundary")
U = TransientTrialFESpace(V0,u)

trian = Triangulation(model)
degree = 2*order
quad = CellQuadrature(trian,degree)

#
a(u,v) = ∇(v)⋅∇(u)
b(v,t) = v*f(t)

res(t,(u,ut),v) = a(u,v) + ut*v - b(v,t)
jac(t,(u,ut),du,v) = a(du,v)
jac_t(t,(u,ut),dut,v) = dut*v

t_Ω = FETerm(res,jac,jac_t,trian,quad)
op = TransientFEOperator(U,V0,t_Ω)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)

ls = LUSolver()
using Gridap.Algebra: NewtonRaphsonSolver
nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
odes = ThetaMethod(ls,dt,θ)
solver = TransientFESolver(odes)

sol_t = solve(solver,op,uh0,t0,tF)

l2(w) = w*w

tol = 1.0e-4
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  e = u(tn) - uh_tn
  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
  @test el2 < tol
end

end #module


text:
# ---
# title: 600. Non-negative Integers without Consecutive Ones
# id: problem600
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/non-negative-integers-without-consecutive-ones/description/>
# hidden: true
# ---
#
# Given a positive integer n, find the number of **non-negative** integers less
# than or equal to n, whose binary representations do NOT contain **consecutive
# ones**.
#
# **Example 1:**
#
#
#
#     Input: 5
#     Output: 5
#     Explanation:
#     Here are the non-negative integers <= 5 with their corresponding binary representations:
#     0 : 0
#     1 : 1
#     2 : 10
#     3 : 11
#     4 : 100
#     5 : 101
#     Among them, only integer 3 disobeys the rule (two consecutive ones) and the other 5 satisfy the rule.
#
#
# **Note:** 1 <= n <= 109
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
# consider using the project of this directory for tested versions
using PGFPlotsX, KernelDensity, HighestDensityRegions, Distributions, Contour, LinearAlgebra

d = MixtureModel([MultivariateNormal([-1, -1], Diagonal([0.5, 0.5])),
                  MultivariateNormal([1, 1], Diagonal([0.7, 1.5]))], [0.6, 0.4])
x = rand(d, 2000)
x1, x2 = x[1, :], x[2, :]
k = kde((x1, x2))
ik = InterpKDE(k)
ps = pdf.(Ref(ik), x1, x2)
qs = 0.05:0.1:0.95
thresholds = hdr_thresholds(qs, ps)

xgrid = range(minimum(x); stop = maximum(x), length = 100)
pgrid = pdf.(Ref(ik), xgrid, xgrid')

p = @pgf Axis({ xlabel = raw"$x_1$", ylabel = raw"$x_2$", width = "10cm", height = "10cm" },
              Plot({ thick }, Table(contours(xgrid, xgrid, pgrid, thresholds))))
pgfsave("2d-density.svg", p)


text:
"""
    getdriver(i::Integer)

Fetch driver by index.
"""
getdriver(i::Integer) = Driver(GDAL.gdalgetdriver(i))

"""
    getdriver(name::AbstractString)

Fetch a driver based on the short name (such as `GTiff`).
"""
getdriver(name::AbstractString) = Driver(GDAL.gdalgetdriverbyname(name))

"""
    destroy(drv::Driver)

Destroy a `GDALDriver`.

This is roughly equivalent to deleting the driver, but is guaranteed to take
place in the GDAL heap. It is important this that function not be called on a
driver that is registered with the `GDALDriverManager`.
"""
function destroy(drv::Driver)
    GDAL.gdaldestroydriver(drv.ptr)
    drv.ptr = C_NULL
end

"""
    register(drv::Driver)

Register a driver for use.
"""
register(drv::Driver) = GDAL.gdalregisterdriver(drv.ptr)

"""
    deregister(drv::Driver)

Deregister the passed driver.
"""
deregister(drv::Driver) = GDAL.gdalderegisterdriver(drv.ptr)

"""
    options(drv::Driver)

Return the list of creation options of the driver [an XML string].
"""
options(drv::Driver) = GDAL.gdalgetdrivercreationoptionlist(drv.ptr)

driveroptions(name::AbstractString) = options(getdriver(name))

"""
    shortname(drv::Driver)

Return the short name of a driver (e.g. `GTiff`).
"""
shortname(drv::Driver) = GDAL.gdalgetdrivershortname(drv.ptr)

"""
    longname(drv::Driver)

Return the long name of a driver (e.g. `GeoTIFF`), or empty string.
"""
longname(drv::Driver) = GDAL.gdalgetdriverlongname(drv.ptr)

"""
    ndriver()

Fetch the number of registered drivers.
"""
ndriver() = GDAL.gdalgetdrivercount()

"""
    listdrivers()

Returns a listing of all registered drivers.
"""
listdrivers() = Dict{String,String}([
    shortname(getdriver(i)) => longname(getdriver(i)) for i in 0:(ndriver()-1)
])

"""
    identifydriver(filename::AbstractString)

Identify the driver that can open a raster file.

This function will try to identify the driver that can open the passed filename
by invoking the Identify method of each registered `GDALDriver` in turn. The
first driver that successful identifies the file name will be returned. If all
drivers fail then `NULL` is returned.
"""
identifydriver(filename::AbstractString) =
    Driver(GDAL.gdalidentifydriver(filename, C_NULL))

"""
    validate(drv::Driver, options::Vector{<:AbstractString})

Validate the list of creation options that are handled by a drv.

This is a helper method primarily used by `create()` and `copy()` to
validate that the passed in list of creation options is compatible with the
`GDAL_DMD_CREATIONOPTIONLIST` metadata item defined by some drivers.

### Parameters
* `drv`     the handle of the driver with whom the lists of creation option
            must be validated
* `options` the list of creation options. An array of strings, whose last
            element is a `NULL` pointer

### Returns
`true` if the list of creation options is compatible with the `create()` and
`createcopy()` method of the driver, `false` otherwise.

### Additional Remarks
See also: `options(drv::Driver)`

If the `GDAL_DMD_CREATIONOPTIONLIST` metadata item is not defined, this
function will return ``true``. Otherwise it will check that the keys and values
in the list of creation options are compatible with the capabilities declared
by the `GDAL_DMD_CREATIONOPTIONLIST` metadata item. In case of incompatibility
a (non fatal) warning will be emited and ``false`` will be returned.
"""
validate(drv::Driver, options::Vector{T}) where {T <: AbstractString} =
    Bool(GDAL.gdalvalidatecreationoptions(drv.ptr, options))

"""
    copyfiles(drv::Driver, new::AbstractString, old::AbstractString)
    copyfiles(drvname::AbstractString, new::AbstractString, old::AbstractString)

Copy all the files associated with a dataset.
"""
function copyfiles end

function copyfiles(drv::Driver, new::AbstractString, old::AbstractString)
    result = GDAL.gdalcopydatasetfiles(drv.ptr, new, old)
    @cplerr result "Failed to copy dataset files"
end

copyfiles(drvname::AbstractString, new::AbstractString, old::AbstractString) =
    copyfiles(getdriver(drvname), new, old)


text:
let
    s = CBSEnv.State()
    @test states_match(s,s)
    a = CBSEnv.Action()
    env = CBSEnv.LowLevelEnv()

    a = CRCBS.wait(s)
    sp = get_next_state(s,a)
    sp = get_next_state(env,s,a)
    @test states_match(s,sp)
    @test is_goal(env,s)
end
let
    s = CBSEnv.State()
    a = CBSEnv.Action()
    env = CBSEnv.LowLevelEnv()
    get_next_state(s,a)
    get_transition_cost(env,s,a,s)
end
let
    G = Graph(3)
    add_edge!(G,1,2)
    add_edge!(G,2,2)
    add_edge!(G,2,3)
    env = CBSEnv.LowLevelEnv(graph=G)
    for a in get_possible_actions(env,CBSEnv.State(2,0))
        @test a.e.src == 2
        @test a.e.dst ∈ [1,2,3]
    end
end
let
    CBSEnv.LowLevelEnv(cost_model=SumOfTravelTime())
end
let
    vtx_grid = initialize_dense_vtx_grid(4,4)
    # 1  5   9  13
    # 2  6  10  14
    # 3  7  11  15
    # 4  8  12  16
    G = initialize_grid_graph_from_vtx_grid(vtx_grid)
    T = 10
    n_agents = 2
    env = CBSEnv.LowLevelEnv(graph=G,agent_idx=1,cost_model=HardConflictCost(G,T,n_agents))
    set_path!(get_cost_model(env),2,[2,6,10,14])
    @test get_transition_cost(env,CBSEnv.State(5,0),CBSEnv.Action(Edge(5,6),1),CBSEnv.State(6,1)) == 1
end
# solve!
let
    S = CBSEnv.State
    A = CBSEnv.Action
    vtx_grid = initialize_regular_vtx_grid(;n_obstacles_x=1,n_obstacles_y=1)
    #  1   2   3   4   5   6
    #  7   8   9  10  11  12
    # 13  14   0   0  15  16
    # 17  18   0   0  19  20
    # 21  22  23  24  25  26
    # 27  28  29  30  31  32
    G = initialize_grid_graph_from_vtx_grid(vtx_grid)
    starts = [CBSEnv.State(1,0),CBSEnv.State(2,0)]
    goals = [CBSEnv.State(vtx=6),CBSEnv.State(vtx=5)]
    heuristic = PerfectHeuristic(G,map(s->s.vtx,starts),map(s->s.vtx,goals))
    env = CBSEnv.LowLevelEnv(graph=G,heuristic=heuristic)
    mapf = MAPF(env,starts,goals)
    solver = CBSSolver()
    initialize_root_node(solver,mapf)
    default_solution(mapf)
    # println("Testing CBS")
    # low_level_search
    let
        solver = CBSSolver()
        node = CBSEnv.initialize_root_node(mapf)
        initialize_child_search_node(solver,mapf,node)
        low_level_search!(solver,mapf,node)
    end
    # high_level_search
    let
        solver = CBSSolver()
        set_iteration_limit!(solver,0)
        # @test_throws SolverException CRCBS.solve!(solver,mapf)
    end
    let
        solver = CBSSolver()
        set_verbosity!(solver,global_verbosity())
        set_iteration_limit!(solver,100)
        solution, cost = CRCBS.solve!(solver,mapf)
        @test cost == 10
    end
    # with CompositeCost and CompositeHeuristic
    let
        cost_model = construct_composite_cost_model(
            FullCostModel(sum,NullCost()),
            FullCostModel(sum,TravelTime())
        )
        heuristic = construct_composite_heuristic(
            NullHeuristic(),
            PerfectHeuristic(G,map(s->s.vtx,starts),map(s->s.vtx,goals))
        )
        env = CBSEnv.LowLevelEnv(graph=G,cost_model=cost_model,heuristic=heuristic)
        mapf = MAPF(env,starts,goals)
        solver = CBSSolver(AStar{cost_type(mapf)}())
        set_verbosity!(solver,global_verbosity())
        solution, cost = CRCBS.solve!(solver,mapf)
    end
    let
        deadline = ne(G)+1.0
        cost_model = construct_composite_cost_model(
            FullDeadlineCost(DeadlineCost(deadline)),
            FullCostModel(sum,NullCost()),
            SumOfTravelTime()
        )
        heuristic_model = construct_composite_heuristic(
            PerfectHeuristic(G,map(s->s.vtx,starts),map(s->s.vtx,goals)),
            HardConflictHeuristic(G,ne(G),num_agents(mapf)),
            PerfectHeuristic(G,map(s->s.vtx,starts),map(s->s.vtx,goals)),
        )
        env = CBSEnv.LowLevelEnv(graph=G,cost_model=cost_model,heuristic=heuristic_model)
        mapf = MAPF(env,starts,goals)
        default_solution(mapf)
        solver = CBSSolver(AStar{cost_type(mapf)}())
        set_verbosity!(solver,global_verbosity())
        solution, cost = CRCBS.solve!(solver,mapf)
    end
end
# problem instances
let
    solver = CBSSolver()
    set_iteration_limit!(solver,1000)
    for f in [init_mapf_1,init_mapf_2,init_mapf_4,CRCBS.init_mapf_8]
        # @info "$(f)"
        prob = f()
        reset_solver!(solver)
        solution, cost = solve!(solver,prob)
    end

end
# test copying
let
    solver = CBSSolver()
    prob = init_mapf_1()
    solution, cost = solve!(solver,prob)
    solution2 = copy(solution)
    # should not affect solution2
    set_solution_path!(solution,get_paths(solution)[1],2)
    @test convert_to_vertex_lists(solution)[2] != convert_to_vertex_lists(solution2)[2]
    n = get_paths(solution)[1].path_nodes[end]
    # should affect both solutions
    push!(get_paths(solution)[1],n)
    @test length(get_paths(solution)[1]) == length(get_paths(solution2)[1])

end


text:
using Test
import ClimaCore: Domains, Meshes, Topologies
import ClimaCore.Geometry: Geometry
using StaticArrays
using IntervalSets

function rectangular_grid(
    n1,
    n2,
    x1periodic,
    x2periodic;
    x1min = 0.0,
    x1max = 1.0,
    x2min = 0.0,
    x2max = 1.0,
)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(x1min)..Geometry.XPoint(x1max),
        Geometry.YPoint(x2min)..Geometry.YPoint(x2max),
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )
    mesh = Meshes.EquispacedRectangleMesh(domain, n1, n2)
    grid_topology = Topologies.GridTopology(mesh)
    return (domain, mesh, grid_topology)
end

function regular_tensorproduct_grid(
    n1,
    n2,
    x1periodic,
    x2periodic;
    x1min = 0.0,
    x1max = 1.0,
    x2min = 0.0,
    x2max = 1.0,
)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(x1min)..Geometry.XPoint(x1max),
        Geometry.YPoint(x2min)..Geometry.YPoint(x2max),
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )
    mesh = Meshes.TensorProductMesh(domain, n1, n2)
    r_ts_topology = Topologies.GridTopology(mesh)
    return (domain, mesh, r_ts_topology)
end

function irregular_tensorproduct_grid(
    n1,
    n2,
    x1periodic,
    x2periodic;
    x1min = 0.0,
    x1max = 1.0,
    x2min = 0.0,
    x2max = 1.0,
)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(x1min)..Geometry.XPoint(x1max),
        Geometry.YPoint(x2min)..Geometry.YPoint(x2max),
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )

    mesh = Meshes.EquispacedRectangleMesh(domain, n1, n2)
    nv2 = x2periodic ? n2 : n2 + 1
    range1 = mesh.range1
    range2 = mesh.range2
    elemheight = range2[2] - range2[1]

    coordinates =
        Vector{Geometry.XYPoint{typeof(x1min)}}(undef, (n1 + 1) * (n2 + 1))

    # Iterate mesh vertices
    grid_topo = Topologies.GridTopology(mesh)
    v_iter = Topologies.VertexIterator(grid_topo)

    for v in v_iter
        (i, j) = v.num
        # Shift 0-based indices
        i += 1
        j += 1
        vcoords = Geometry.XYPoint(range1[i], range2[j])
        # If I am on a boundary, don't apply any warping
        if j == 1 || j == nv2
            coordinates[(i - 1) * (n2 + 1) + j] = vcoords
        else # Interior vertex case, shift x2-coord of a fixed amount
            coordinates[(i - 1) * (n2 + 1) + j] =
                Geometry.XYPoint(vcoords.x, vcoords.y + 0.5 * elemheight)
        end
    end
    ts_mesh = Meshes.TensorProductMesh(domain, n1, n2, coordinates)
    i_ts_topology = Topologies.GridTopology(ts_mesh)
    return (domain, ts_mesh, i_ts_topology)
end

@testset "opposing face tests" begin

    _, _, grid_topology = rectangular_grid(1, 1, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, true)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "assert correct element numbering" begin
            @test_throws AssertionError Topologies.opposing_face(topology, 0, 1)
            @test_throws AssertionError Topologies.opposing_face(topology, 2, 1)
        end

        @testset "assert correct face numbering" begin
            @test_throws AssertionError Topologies.opposing_face(topology, 1, 5)
            @test_throws AssertionError Topologies.opposing_face(topology, 1, 0)
        end

        @testset "1×1 element quad mesh with all periodic boundries" begin
            @test Topologies.opposing_face(topology, 1, 1) == (1, 2, false)
            @test Topologies.opposing_face(topology, 1, 2) == (1, 1, false)
            @test Topologies.opposing_face(topology, 1, 3) == (1, 4, false)
            @test Topologies.opposing_face(topology, 1, 4) == (1, 3, false)
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, true, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with 1 periodic boundary" begin
            _, _, topology = rectangular_grid(1, 1, true, false)
            @test Topologies.opposing_face(topology, 1, 1) == (1, 2, false)
            @test Topologies.opposing_face(topology, 1, 2) == (1, 1, false)
            @test Topologies.opposing_face(topology, 1, 3) == (0, 3, false)
            @test Topologies.opposing_face(topology, 1, 4) == (0, 4, false)
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with non-periodic boundaries" begin
            @test Topologies.opposing_face(topology, 1, 1) == (0, 1, false)
            @test Topologies.opposing_face(topology, 1, 2) == (0, 2, false)
            @test Topologies.opposing_face(topology, 1, 3) == (0, 3, false)
            @test Topologies.opposing_face(topology, 1, 4) == (0, 4, false)
        end
    end

    _, _, grid_topology = rectangular_grid(2, 2, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 2, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "2×2 element quad mesh with non-periodic boundaries" begin
            @test Topologies.opposing_face(topology, 1, 1) == (0, 1, false)
            @test Topologies.opposing_face(topology, 1, 2) == (2, 1, false)
            @test Topologies.opposing_face(topology, 1, 3) == (0, 3, false)
            @test Topologies.opposing_face(topology, 1, 4) == (3, 3, false)
            @test Topologies.opposing_face(topology, 2, 1) == (1, 2, false)
            @test Topologies.opposing_face(topology, 2, 2) == (0, 2, false)
            @test Topologies.opposing_face(topology, 2, 3) == (0, 3, false)
            @test Topologies.opposing_face(topology, 2, 4) == (4, 3, false)
        end
    end
end

@testset "simple rectangular mesh interior faces iterator" begin

    _, _, grid_topology = rectangular_grid(1, 1, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, true)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with all periodic boundries" begin
            @test length(Topologies.interior_faces(topology)) == 2
            faces = collect(Topologies.interior_faces(topology))
            @test faces[1] == (1, 1, 1, 2, false)
            @test faces[2] == (1, 3, 1, 4, false)
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, true, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with 1 periodic boundary" begin
            @test length(Topologies.interior_faces(topology)) == 1
            faces = collect(Topologies.interior_faces(topology))
            @test faces[1] == (1, 1, 1, 2, false)
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.interior_faces(topology)) == 0
            faces = collect(Topologies.interior_faces(topology))
            @test isempty(faces)
        end
    end

    _, _, grid_topology = rectangular_grid(2, 2, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 2, false, false)
    _, _, i_ts_topology = irregular_tensorproduct_grid(2, 2, false, false)
    topologies = (grid_topology, r_ts_topology, i_ts_topology)
    for topology in topologies
        @testset "2×2 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.interior_faces(topology)) == 4
            faces = collect(Topologies.interior_faces(topology))
            @test faces[1] == (2, 1, 1, 2, false)
            @test faces[2] == (3, 3, 1, 4, false)
            @test faces[3] == (4, 1, 3, 2, false)
            @test faces[4] == (4, 3, 2, 4, false)
        end
    end
end

@testset "simple rectangular mesh boundry faces iterator" begin
    _, _, grid_topology = rectangular_grid(1, 1, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, true)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with all periodic boundries" begin
            @test length(Topologies.boundary_faces(topology, 1)) == 0
            @test length(Topologies.boundary_faces(topology, 2)) == 0
            @test length(Topologies.boundary_faces(topology, 3)) == 0
            @test length(Topologies.boundary_faces(topology, 4)) == 0
            @test isempty(collect(Topologies.boundary_faces(topology, 1)))
            @test isempty(collect(Topologies.boundary_faces(topology, 2)))
            @test isempty(collect(Topologies.boundary_faces(topology, 3)))
            @test isempty(collect(Topologies.boundary_faces(topology, 4)))
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, true, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with 1 periodic boundary" begin
            @test length(Topologies.boundary_faces(topology, 1)) == 0
            @test length(Topologies.boundary_faces(topology, 2)) == 0
            @test length(Topologies.boundary_faces(topology, 3)) == 1
            @test length(Topologies.boundary_faces(topology, 4)) == 1
            @test isempty(collect(Topologies.boundary_faces(topology, 1)))
            @test isempty(collect(Topologies.boundary_faces(topology, 2)))
            @test collect(Topologies.boundary_faces(topology, 3)) == [(1, 3)]
            @test collect(Topologies.boundary_faces(topology, 4)) == [(1, 4)]
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.boundary_faces(topology, 1)) == 1
            @test length(Topologies.boundary_faces(topology, 2)) == 1
            @test length(Topologies.boundary_faces(topology, 3)) == 1
            @test length(Topologies.boundary_faces(topology, 4)) == 1
            @test collect(Topologies.boundary_faces(topology, 1)) == [(1, 1)]
            @test collect(Topologies.boundary_faces(topology, 2)) == [(1, 2)]
            @test collect(Topologies.boundary_faces(topology, 3)) == [(1, 3)]
            @test collect(Topologies.boundary_faces(topology, 4)) == [(1, 4)]
            @test Topologies.boundary_tag(topology, :west) == 1
            @test Topologies.boundary_tag(topology, :east) == 2
            @test Topologies.boundary_tag(topology, :south) == 3
            @test Topologies.boundary_tag(topology, :north) == 4
        end
    end

    _, _, grid_topology = rectangular_grid(2, 3, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 3, false, false)
    _, _, i_ts_topology = irregular_tensorproduct_grid(2, 3, false, false)
    topologies = (grid_topology, r_ts_topology, i_ts_topology)
    for topology in topologies
        @testset "2×3 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.boundary_faces(topology, 1)) == 3
            @test length(Topologies.boundary_faces(topology, 2)) == 3
            @test length(Topologies.boundary_faces(topology, 3)) == 2
            @test length(Topologies.boundary_faces(topology, 4)) == 2

            @test collect(Topologies.boundary_faces(topology, 1)) ==
                  [(1, 1), (3, 1), (5, 1)]
            @test collect(Topologies.boundary_faces(topology, 2)) ==
                  [(2, 2), (4, 2), (6, 2)]
            @test collect(Topologies.boundary_faces(topology, 3)) ==
                  [(1, 3), (2, 3)]
            @test collect(Topologies.boundary_faces(topology, 4)) ==
                  [(5, 4), (6, 4)]
        end
    end
end

@testset "simple rectangular mesh vertex iterator" begin

    _, _, grid_topology = rectangular_grid(1, 1, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, true)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with all periodic boundries" begin
            @test length(Topologies.vertices(topology)) == 1
            V = collect(Topologies.vertices(topology))
            @test V[1] isa Topologies.Vertex
            @test length(V[1]) == 4
            @test collect(V[1]) == [(1, 1), (1, 2), (1, 3), (1, 4)]
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, true, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with 1 periodic boundary" begin
            @test length(Topologies.vertices(topology)) == 2
            V = collect(Topologies.vertices(topology))
            @test V[1] isa Topologies.Vertex
            @test length(V[1]) == 2
            @test collect(V[1]) == [(1, 1), (1, 2)]
            @test length(V[2]) == 2
            @test collect(V[2]) == [(1, 3), (1, 4)]
        end
    end

    _, _, grid_topology = rectangular_grid(1, 1, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.vertices(topology)) == 4
            V = collect(Topologies.vertices(topology))
            @test V[1] isa Topologies.Vertex
            @test length(V[1]) == 1
            @test collect(V[1]) == [(1, 1)]
            @test length(V[2]) == 1
            @test collect(V[2]) == [(1, 2)]
            @test length(V[3]) == 1
            @test collect(V[3]) == [(1, 3)]
            @test length(V[4]) == 1
            @test collect(V[4]) == [(1, 4)]
        end
    end

    _, _, grid_topology = rectangular_grid(2, 3, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 3, false, false)
    _, _, i_ts_topology = irregular_tensorproduct_grid(2, 3, false, false)
    topologies = (grid_topology, r_ts_topology, i_ts_topology)
    for topology in topologies
        @testset "2×3 element quad mesh with non-periodic boundaries" begin
            @test length(Topologies.vertices(topology)) == 3 * 4
            V = collect(Topologies.vertices(topology))
            @test length(V[1]) == 1
            @test collect(V[1]) == [(1, 1)]
        end
    end

    _, _, grid_topology = rectangular_grid(2, 3, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 3, true, true)
    _, _, i_ts_topology = irregular_tensorproduct_grid(2, 3, true, true)
    topologies = (grid_topology, r_ts_topology, i_ts_topology)
    for topology in topologies
        @testset "2×3 element quad mesh with periodic boundaries" begin
            @test length(Topologies.vertices(topology)) == 2 * 3
            V = collect(Topologies.vertices(topology))
            @test length(V) == 6
            @test collect(V[1]) == [(1, 1), (2, 2), (5, 3), (6, 4)]
            @test collect(V[2]) == [(2, 1), (1, 2), (6, 3), (5, 4)]
            @test collect(V[3]) == [(3, 1), (4, 2), (1, 3), (2, 4)]
            @test collect(V[4]) == [(4, 1), (3, 2), (2, 3), (1, 4)]
            @test collect(V[5]) == [(5, 1), (6, 2), (3, 3), (4, 4)]
            @test collect(V[6]) == [(6, 1), (5, 2), (4, 3), (3, 4)]
        end
    end
end

@testset "simple rectangular mesh coordinates" begin

    _, _, grid_topology = rectangular_grid(1, 1, true, true)
    _, _, r_ts_topology = regular_tensorproduct_grid(1, 1, true, true)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with all periodic boundries" begin
            c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
            @test c1 == Geometry.XYPoint(0.0, 0.0)
            @test c2 == Geometry.XYPoint(1.0, 0.0)
            @test c3 == Geometry.XYPoint(0.0, 1.0)
            @test c4 == Geometry.XYPoint(1.0, 1.0)

        end
    end

    _, _, grid_topology = rectangular_grid(
        1,
        1,
        false,
        false;
        x1min = -1.0,
        x1max = 1.0,
        x2min = -1.0,
        x2max = 1.0,
    )
    _, _, r_ts_topology = regular_tensorproduct_grid(
        1,
        1,
        false,
        false;
        x1min = -1.0,
        x1max = 1.0,
        x2min = -1.0,
        x2max = 1.0,
    )
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "1×1 element quad mesh with non-periodic boundries" begin
            c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
            @test c1 == Geometry.XYPoint(-1.0, -1.0)
            @test c2 == Geometry.XYPoint(1.0, -1.0)
            @test c3 == Geometry.XYPoint(-1.0, 1.0)
            @test c4 == Geometry.XYPoint(1.0, 1.0)
        end
    end

    _, _, grid_topology = rectangular_grid(2, 4, false, false)
    _, _, r_ts_topology = regular_tensorproduct_grid(2, 4, false, false)
    topologies = (grid_topology, r_ts_topology)
    for topology in topologies
        @testset "2×4 element quad mesh with non-periodic boundaries" begin
            c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
            @test c1 == Geometry.XYPoint(0.0, 0.0)
            @test c2 == Geometry.XYPoint(0.5, 0.0)
            @test c3 == Geometry.XYPoint(0.0, 0.25)
            @test c4 == Geometry.XYPoint(0.5, 0.25)

            c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 8)
            @test c1 == Geometry.XYPoint(0.5, 0.75)
            @test c2 == Geometry.XYPoint(1.0, 0.75)
            @test c3 == Geometry.XYPoint(0.5, 1.0)
            @test c4 == Geometry.XYPoint(1.0, 1.0)
        end
    end
    _, _, i_ts_topology = irregular_tensorproduct_grid(2, 4, false, false)
    @testset "2×4 element quad mesh with non-periodic boundaries" begin
        c1, c2, c3, c4 = Topologies.vertex_coordinates(i_ts_topology, 1)
        @test c1 == Geometry.XYPoint(0.0, 0.0)
        @test c2 == Geometry.XYPoint(0.5, 0.0)
        @test c3 == Geometry.XYPoint(0.0, 0.375)
        @test c4 == Geometry.XYPoint(0.5, 0.375)

        c1, c2, c3, c4 = Topologies.vertex_coordinates(i_ts_topology, 8)
        @test c1 == Geometry.XYPoint(0.5, 0.875)
        @test c2 == Geometry.XYPoint(1.0, 0.875)
        @test c3 == Geometry.XYPoint(0.5, 1.0)
        @test c4 == Geometry.XYPoint(1.0, 1.0)
    end

    _, _, grid_topology = rectangular_grid(
        3,
        1,
        false,
        false,
        x1min = big(0.0),
        x1max = big(1.0),
        x2min = big(0.0),
        x2max = big(1.0),
    )
    _, _, r_ts_topology = regular_tensorproduct_grid(
        3,
        1,
        false,
        false,
        x1min = big(0.0),
        x1max = big(1.0),
        x2min = big(0.0),
        x2max = big(1.0),
    )

    _, _, i_ts_topology = irregular_tensorproduct_grid(
        3,
        1,
        false,
        false,
        x1min = big(0.0),
        x1max = big(1.0),
        x2min = big(0.0),
        x2max = big(1.0),
    )
    topologies = (grid_topology, r_ts_topology, i_ts_topology)
    for topology in topologies
        @testset "check coordinate type accuracy" begin
            c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
            @test eltype(c2) == BigFloat
            @test getfield(c2, 1) ≈ big(1.0) / big(3.0) rtol = eps(BigFloat)
            @test getfield(c2, 2) == 0.0
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpdateTrustedIdProviderProperties <: SwaggerModel
    idProvider::Any # spec type: Union{ Nothing, String } # spec name: idProvider

    function UpdateTrustedIdProviderProperties(;idProvider=nothing)
        o = new()
        validate_property(UpdateTrustedIdProviderProperties, Symbol("idProvider"), idProvider)
        setfield!(o, Symbol("idProvider"), idProvider)
        o
    end
end # type UpdateTrustedIdProviderProperties

const _property_map_UpdateTrustedIdProviderProperties = Dict{Symbol,Symbol}(Symbol("idProvider")=>Symbol("idProvider"))
const _property_types_UpdateTrustedIdProviderProperties = Dict{Symbol,String}(Symbol("idProvider")=>"String")
Base.propertynames(::Type{ UpdateTrustedIdProviderProperties }) = collect(keys(_property_map_UpdateTrustedIdProviderProperties))
Swagger.property_type(::Type{ UpdateTrustedIdProviderProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_UpdateTrustedIdProviderProperties[name]))}
Swagger.field_name(::Type{ UpdateTrustedIdProviderProperties }, property_name::Symbol) =  _property_map_UpdateTrustedIdProviderProperties[property_name]

function check_required(o::UpdateTrustedIdProviderProperties)
    true
end

function validate_property(::Type{ UpdateTrustedIdProviderProperties }, name::Symbol, val)
end


text:
using StatsBase

"""
    function mutation!(pop::populations)

Mutate a base pair in a random species, chosen depending on subpopulation size.
If the species is extinct after mutation, it is removed. If mutation creates a new species,
it is added to the population, otherwise the subpopulation size of species the individual mutates to
is simply increased by one.
"""
function mutation!(pop::populations)
    # Choose random sequence depending on subpopulation size
    mutation_sequence = sample(collect(1:length(pop.seqs)), Weights(pop.freqs ./ pop.N))
    # Add sequence to mutate
    push!(pop.seqs, copy(pop.seqs[mutation_sequence]))
    push!(pop.freqs, 1)
    pop.freqs[mutation_sequence] -= 1
    # Choose random base to mutate
    mutation_base = rand(1:pop.l)
    # Exclude self mutations
    poss_bases = filter(x->x != pop.seqs[end][mutation_base], collect(1:pop.n))
    # Choose random new base
    pop.seqs[end][mutation_base] = rand(poss_bases)
    copied = findfirst(x->~any(.~(x.==pop.seqs[end])), pop.seqs[1:end-1])
    # Look if mutated sequence already exists
    if ~isnothing(copied)
        pop.freqs[copied] += 1
        deleteat!(pop.freqs, length(pop.seqs))
        deleteat!(pop.seqs, length(pop.seqs))
    end
    # Remove extinct species
    remove_empty!(pop.freqs, pop.seqs)
end

"""
    function mutation!(pop::driver_trailer_l)

Mutate a base pair in a random species, chosen depending on subpopulation size.
If the species is extinct after mutation, it is removed. If mutation creates a new species,
it is added to the population, otherwise the subpopulation size of species the individual mutates to
is simply increased by one.
"""
function mutation!(pop::driver_trailer_l)
    # Choose random sequence depending on subpopulation size
    mutation_sequence = sample(collect(1:length(pop.seqs)), Weights(pop.freqs ./ pop.N))
    # Add sequence to mutate
    push!(pop.seqs, deepcopy(pop.seqs[mutation_sequence]))
    push!(pop.freqs, 1)
    push!(pop.l, pop.l[mutation_sequence])
    pop.freqs[mutation_sequence] -= 1
    # Choose random base to mutate
    mutation_base = rand(1:pop.L)
    # Exclude self mutations
    poss_bases = filter(x->x != pop.seqs[end][mutation_base], collect(1:pop.n))
    # Choose random new base
    pop.seqs[end][mutation_base] = rand(poss_bases)
    copied = findfirst(x->~any(.~(x.==pop.seqs[end])), pop.seqs[1:end-1])
    # Look if mutated sequence already exists
    if ~isnothing(copied)
        pop.freqs[copied] += 1
        deleteat!(pop.freqs, length(pop.seqs))
        deleteat!(pop.seqs, length(pop.seqs))
        deleteat!(pop.l, length(pop.l))
    end
    # Remove extinct species
    remove_empty!(pop.freqs, pop.seqs, pop.l)
end

"""
    function driver_mutation!(pop::DT_population)

Mutate a random base in the driver sequence.
"""
function driver_mutation!(pop::DT_population)
    # Choose random base to mutate
    mutation_base = rand(1:length(pop.driver))
    # Exclude self mutations
    poss_bases = filter(x->x != pop.driver[mutation_base], collect(1:pop.m))
    # Choose random new base
    pop.driver[mutation_base] = rand(poss_bases)
end

"""
    function driver_mutation!(pop::driver_trailer)

Mutate the length of species and create a new species
"""
function length_mutation!(pop::driver_trailer_l)
    # Choose random sequence depending on subpopulation size
    mutation_sequence = sample( collect(1:length(pop.seqs)),
                                Weights(pop.freqs ./ pop.N))
    # Add sequence to mutate

    if rand() < 0.5 && pop.l[mutation_sequence] < pop.L
        var = 1
    elseif pop.l[mutation_sequence] > 2
        var = -1
    else
        return
    end

    push!(pop.seqs, copy(pop.seqs[mutation_sequence]))
    push!(pop.freqs, 1)
    push!(pop.l, pop.l[mutation_sequence])

    pop.freqs[mutation_sequence] -= 1
    pop.l[end] += var
    # Remove extinct species
    remove_empty!(pop.freqs, pop.seqs, pop.l)
end

"""
    bp_substitution!(pop::mono_pop, emat::Array{T, 2}, fitness_function::fitness_functions) where {T<:Real}

Attempt a base pair substitution. The population has to consist of one species
only. Mutation is accepted with probability given by the Kimura fixation
probability.
"""
function bp_substitution!(pop::mono_pop, emat::Array{T, 2}, fitness_function::fitness_functions) where {T<:Real}

    # Choose random base to mutate
    mutation_base = rand(eachindex(pop.seqs))
    # Exclude self mutations
    possible_bases = filter(x->x != pop.seqs[mutation_base], collect(1:pop.n))
    # Choose random new base
    temp_pop = deepcopy(pop)
    temp_pop.seqs[mutation_base] = rand(possible_bases)
    # Compute energies
    E = get_energy(pop, emat)
    E_mutant = get_energy(temp_pop, emat)
    # Compute fitness and selection coefficients
    F = fitness(E, length(pop.seqs), fitness_function)
    F_mutant = fitness(E_mutant, length(pop.seqs), fitness_function)
    s = F_mutant - F
    # Accept mutatant depending on Kimura probability
    if rand() < kimura_prob(s, pop.N)
        pop.seqs = temp_pop.seqs
    end
end

"""
    bp_substitution!(pop::populations, emat::Array{T, 2}, fitness_function::fitness_functions) where {T<:Real}

Attempt a base pair substitution. The population has to consist of one species
only. Mutation is accepted with probability given by the Kimura fixation
probability.
"""
function bp_substitution!(pop::populations, emat::Array{T, 2}, fitness_function::fitness_functions) where {T<:Real}
    if length(pop.seqs) > 1
        throw(ArgumentError("Input population has more than one species."))
    end
    # Choose random base to mutate
    mutation_base = rand(1:pop.L)
    # Exclude self mutations
    poss_bases = filter(x->x != pop.seqs[1][mutation_base], collect(1:pop.n))
    # Choose random new base
    temp_pop = deepcopy(pop)
    temp_pop.seqs[1][mutation_base] = rand(poss_bases)
    # Compute energies
    E = get_energy(pop, emat)
    E_mutant = get_energy(temp_pop, emat)
    # Compute fitness and selection coefficients
    F = fitness(E[1], pop.l[1], fitness_function)
    F_mutant = fitness(E_mutant[1], pop.l[1], fitness_function)
    s = F_mutant - F
    # Accept mutatant depending on Kimura probability
    if rand() < kimura_prob(s, pop.N)
        pop.seqs = temp_pop.seqs
    end
end

"""
    l_substitution!(
        pop::driver_trailer_l,
        emat::Array{T, 2},
        fitness_function::fitness_functions) where {T<:Real}

Attempt a substitution of a length mutation. Population has to be consisting
of one species only.
"""
function l_substitution!(
    pop::driver_trailer_l,
    emat::Array{T, 2},
    fitness_function::fitness_functions) where {T<:Real}

    if length(pop.seqs) > 1
        throw(ArgumentError("Input population has more than one species."))
    end

    if rand() < 0.5 && pop.l[1] < pop.L
        var = 1
    elseif pop.l[1] > 2
        var = -1
    else
        return
    end
    # Choose random new base
    temp_pop = deepcopy(pop)
    temp_pop.l[1] += var
    # Compute energies
    E = get_energy(pop, emat)
    E_mutant = get_energy(temp_pop, emat)
    # Compute fitness and selection coefficients
    F = fitness(E[1], pop.l[1], fitness_function)
    F_mutant = fitness(E_mutant[1], temp_pop.l[1], fitness_function)
    s = F_mutant - F
    # Accept mutatant depending on Kimura probability
    if rand() < kimura_prob(s, pop.N)
        pop.l = temp_pop.l
    end
end

"""
    l_substitution!(
        pop::mono_pop,
        emat::Array{T, 2},
        fitness_function::fitness_functions) where {T<:Real}

Attempt a substitution of a length mutation. Population has to be consisting
of one species only.
"""
function l_substitution!(
    pop::mono_pop,
    emat::Array{T, 2},
    fitness_function::fitness_functions) where {T<:Real}

    temp_pop = deepcopy(pop)
    if rand() < 0.5
        push!(temp_pop.seqs, rand(1:pop.n))
        push!(temp_pop.driver, rand(1:pop.m))
    elseif length(temp_pop.seqs) > 1
        pop!(temp_pop.seqs)
        pop!(temp_pop.driver)
    end
    # Compute energies
    E = get_energy(pop, emat)
    E_mutant = get_energy(temp_pop, emat)
    # Compute fitness and selection coefficients
    F = fitness(E, length(pop.seqs), fitness_function)
    F_mutant = fitness(E_mutant, length(temp_pop.seqs), fitness_function)
    s = F_mutant - F
    # Accept mutatant depending on Kimura probability
    if rand() < kimura_prob(s, pop.N)
        pop.seqs = temp_pop.seqs
        pop.driver = temp_pop.driver
    end
end

"""
    kimura_prob(s, N)

Probability of fixation of a mutation with selection coefficient
`s` in a monomorphic population of size `N`.
"""
function kimura_prob(s::Real, N::Int)
    if abs(s) <= 10^-8
        return 1 / N
    elseif s < -10^3
        return 0
    else
        return (1 - exp(-2s)) / (1 - exp(-2N * s))
    end
end


text:
# Autogenerated wrapper script for ASL_jll for i686-w64-mingw32
export libasl

JLLWrappers.@generate_wrapper_header("ASL")
JLLWrappers.@declare_library_product(libasl, "libasl.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libasl,
        "bin\\libasl.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# We shall say that an n-digit number is pandigital if it makes use of all the
# digits 1 to n exactly once; for example, the 5-digit number, 15234, is 1
# through 5 pandigital.
#
# The product 7254 is unusual, as the identity, 39 × 186 = 7254, containing
# multiplicand, multiplier, and product is 1 through 9 pandigital.
#
# Find the sum of all products whose multiplicand/multiplier/product identity
# can be written as a 1 through 9 pandigital.
#
# HINT: Some products can be obtained in more than one way so be sure to only
# include it once in your sum.

using ProjectEulerSolutions

# Uses factorization to find sets of digits from the products.  Shortcuts full
# comparisons when the total number of digits is incorrect.
function p032solution_factors(n::Integer=3)::Integer

    startval = 10^(fld(n, 2) - 1)
    endval = 10^fld(n, 2)

    desired_digits = Set{Integer}(1:n)
    products = Integer[]
    for k in startval:endval
        facs = factors(k)
        factorlen = length(facs)
        numdigits = map(ndigits, facs[2:end])
        numdigits = numdigits .+ reverse(numdigits) .+ ndigits(k)
        # Only look at half the factors
        for i in 2:cld(factorlen, 2)
            if numdigits[i-1] == n
                alldigits = Set(vcat(digits(facs[i]),
                                     digits(facs[factorlen - i + 2]),
                                     digits(k)))
                if (desired_digits == alldigits)
                    push!(products, k)
                end
            end
        end
    end
    return sum(Set(products))
end

# Uses a digits from multipliers of the correct lengths.  If a, b, c are the
# number of digits in the multiplicands and product, then a + b + c = n.  We
# know that a + b < c < a + b + 1, so 2(a + b) < n < 2(a + b) + 1.  So we pick
# a: 1 < a < 0.5*n and b: 0.5*n - 0.5 - a < b < 0.5*n - a
function p032solution_multipliers(n::Integer=5)::Integer
    products = Integer[]
    desired_digits = Set{Integer}(1:n)
    r = n / 2
    as = 2:floor(Integer, 10^r)
    for a in as
        na = ndigits(a)
        minb = max(10, floor(Integer, 10^(r - 0.5 - na)))
        maxb = max(minb, floor(Integer, 10^(r - na)))

        for b in minb:maxb
            c = a * b
            if na + ndigits(b) + ndigits(c) == n
                alldigits = Set(vcat(digits(a), digits(b), digits(c)))
                if (desired_digits == alldigits)
                    push!(products, c)
                end
            end
        end
    end
    return sum(Set(products))
end

p032 = Problems.Problem(Dict("Factors" => p032solution_factors,
                             "Multipliers" => p032solution_multipliers))

Problems.benchmark(p032, 9)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

using Base.Test

@testset "Some tests" begin
	@test sum([1,2,3]) == 6
	@test sum([1,2,3,4]) == 10
end


text:
using GeometricFlux
using Flux
using Flux: throttle
using Flux.Losses: logitbinarycrossentropy
using Flux: @epochs
using JLD2
using Statistics: mean
using SparseArrays
using Graphs.SimpleGraphs
using CUDA

@load "data/cora_features.jld2" features
@load "data/cora_graph.jld2" g

num_nodes = 2708
num_features = 1433
hidden1 = 32
hidden2 = 16
target_catg = 7
epochs = 200

## Preprocessing data
fg = FeaturedGraph(g) |> gpu
train_X = Matrix{Float32}(features) |> gpu  # dim: num_features * num_nodes
train_y = fg  # dim: num_nodes * num_nodes

## Model
encoder = Chain(GCNConv(fg, num_features=>hidden1, relu),
                GCNConv(fg, hidden1=>hidden2))
model = Chain(GAE(encoder, σ)) |> gpu

## Loss
loss(x, y) = logitbinarycrossentropy(model(x), y)

## Training
ps = Flux.params(model)
train_data = [(train_X, train_y)]
opt = ADAM(0.01)
evalcb() = @show(loss(train_X, train_y))

@epochs epochs Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))


text:
## pav.jl  Pool Adjacent Violators algorithm for isotonic regression
## (c) 2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md
## This code is largely based on code from the Bosaris toolkit

## pav(y) returns the isotonic regression of the predictor y.
function pav{T<:Real}(y::Vector{T})
    yy = similar(y, Float64)
    n = length(y)
    i = similar(y, Int)         # index
    l = similar(y, Int)         # length
    ci = 1
    i[ci] = l[ci] = 1
    yy[ci] = y[1]
    for j=2:n
        ci += 1
        i[ci] = j
        l[ci] = 1
        yy[ci] = y[j]
        while ci ≥ 2 && yy[max(ci-1,1)] ≥ yy[ci]
            nl = l[ci-1] + l[ci]
            yy[ci-1] += l[ci] * (yy[ci] - yy[ci-1]) / nl
            ci -= 1
            l[ci] = nl
        end
    end
    while n > 0
        for j=i[ci]:n
            yy[j] = yy[ci]
        end
        n = i[ci] - 1
        ci -= 1
    end
    yy
end

## return the optimal log-likelihood ratios for target and non-target scores
function optllr{T<:Real}(tar::Vector{T}, non::Vector{T}; laplace=true)
    ntar = length(tar)
    nnon = length(non)
    ntar > 0 && nnon > 0 || error("Lenghts must be nonzero")
    o = sortperm([tar; non])
    pideal = zeros(ntar + nnon + 4laplace)
    ## set pideal=1 for target scores
    for i=1:ntar+nnon
        if o[i] ≤ ntar         # tar sorted before non, a target score
            pideal[i+2laplace] = 1.0
        end
    end
    if laplace
        pideal[1] = pideal[end-1] = 1.0
    end
    ## optimal posterior, given monotonicity constraints
    popt = pav(pideal)          # optimal posterior
    postlo = logit(popt)[1+2laplace:end-2laplace] # posterior log odds
    priorlo = log(length(tar) / length(non))
    llrs = Array(T, ntar+nnon)
    llrs[o] = postlo - priorlo
    llrs[1:ntar], llrs[ntar+1:end]
end


text:
module TightBinding

using JuLIP, StaticArrays
using JuLIP: @protofun
using JuLIP.Potentials: @pot, SitePotential

import JuLIP: energy, forces, cutoff
import JuLIP.Potentials: evaluate, evaluate_d, site_energy, site_energy_d,
                        partial_energy, partial_energy_d

export hamiltonian, densitymatrix, TBModel, TightBindingModel

# ============================================================================

# abstractions
include("types.jl")

# chemical potentials
# define how to go from spectrum (energy-levels) to potential energy
include("potentials.jl")

# BZ integration: MPGrid and GammaPoint, iterators, etc
include("bzintegration.jl")

# Construction of contours for PEXSI
include("FermiContour.jl")

# ============= SLATER KOSTER TYPE MODELS ================

# basics for slater-koster type hamiltonians
include("sk_core.jl")

# assembling hamiltonian and hamiltonian derivatives
include("matrixassembly.jl")

# the TB toy model for quick tests (a slater-koster s-orbital model)
include("toymodel.jl")

# the NRLTB model
include("NRLTB.jl")

# The Kwon model - a simple orthogonal TB model for Silicon
include("kwon.jl")

# generic code, such as computing spectral decoposition, energy, forces
include("calculators.jl")

# pole expansion (contour integration) based calculator for TBModel
include("pexsi.jl")

# TODO: perturbation theory module
# include("perturbation.jl")

end    # module


text:
@testset "cora" begin
    graph = graphdata(Cora())
    @test typeof(graph) == SimpleDiGraph{UInt32}
    @test nv(graph) == 19793
    @test ne(graph) == 65311

    all_X, all_y = alldata(Cora())
    @test typeof(all_X) == SparseMatrixCSC{Float32,Int64}
    @test size(all_X) == (8710, 19793)
    @test typeof(all_y) == Matrix{UInt16}
    @test size(all_y) == (1, 19793)

    raw = rawdata(Cora())
    @test typeof(raw[:graph]) == SparseMatrixCSC{Float32,Int64}
    @test typeof(raw[:all_X]) == SparseMatrixCSC{Float32,Int64}
    @test typeof(raw[:all_y]) == Vector{Int64}

    meta = metadata(Cora())
    @test meta.graph.num_V == nv(graph)
    @test meta.graph.num_E == ne(graph)
end


text:
function get_dimension(obj::Euclidean2D)
    return jcall(obj, "getDimension", jint, ())
end

function get_instance(::Type{Euclidean2D})
    return jcall(Euclidean2D, "getInstance", Euclidean2D, ())
end

function get_sub_space(obj::Euclidean2D)
    return jcall(obj, "getSubSpace", Euclidean1D, ())
end


text:
abstract type AbstractCorrelator{T} end

"""
$(SIGNATURES)

EarlyPromptLateCorrelator for the three correlators: Early, Prompt and Late
"""
struct EarlyPromptLateCorrelator{T} <: AbstractCorrelator{T}
    early::T
    prompt::T
    late::T
end

"""
$(SIGNATURES)

EarlyPromptLateCorrelator constructor without parameters assumes a single antenna.
"""
function EarlyPromptLateCorrelator()
    EarlyPromptLateCorrelator(NumAnts(1))
end

function EarlyPromptLateCorrelator(num_ants::NumAnts{1})
    EarlyPromptLateCorrelator(
        complex(0.0, 0.0),
        complex(0.0, 0.0),
        complex(0.0, 0.0)
    )
end

"""
$(SIGNATURES)

EarlyPromptLateCorrelator constructor that considers multiple antennas. The number of
antennas has to be specified by `num_ants::NumAnts{N}` where N is the number of antenna
elements.
"""
function EarlyPromptLateCorrelator(num_ants::NumAnts{N}) where N
    EarlyPromptLateCorrelator(
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64})
    )
end

"""
$(SIGNATURES)

Get the early correlator
"""
@inline get_early(correlator::EarlyPromptLateCorrelator) = correlator.early

"""
$(SIGNATURES)

Get the prompt correlator
"""
@inline get_prompt(correlator::EarlyPromptLateCorrelator) = correlator.prompt

"""
$(SIGNATURES)

Get the late correlator
"""
@inline get_late(correlator::EarlyPromptLateCorrelator) = correlator.late

"""
$(SIGNATURES)

Reset the correlator
"""
function zero(correlator::EarlyPromptLateCorrelator{T}) where T
    EarlyPromptLateCorrelator(zero(T), zero(T), zero(T))
end

"""
$(SIGNATURES)

Filter the correlator by the function `post_corr_filter`.
"""
function filter(post_corr_filter, correlator::EarlyPromptLateCorrelator)
    EarlyPromptLateCorrelator(
        post_corr_filter(get_early(correlator)),
        post_corr_filter(get_prompt(correlator)),
        post_corr_filter(get_late(correlator))
    )
end

"""
$(SIGNATURES)

Calculate the shift between the early and late in samples.
"""
function get_early_late_sample_shift(
    ::Type{S},
    correlator::EarlyPromptLateCorrelator,
    sample_frequency,
    preferred_code_shift
) where S <: AbstractGNSSSystem
    round(Int, preferred_code_shift * sample_frequency / get_code_frequency(S))
end

"""
$(SIGNATURES)

Normalize the correlator.
"""
function normalize(correlator::EarlyPromptLateCorrelator, integrated_samples)
    EarlyPromptLateCorrelator(
        get_early(correlator) / integrated_samples,
        get_prompt(correlator) / integrated_samples,
        get_late(correlator) / integrated_samples
    )
end

"""
$(SIGNATURES)

Perform a correlation iteration.
"""
@inline function correlate_iteration(
    ::Type{S},
    correlator::EarlyPromptLateCorrelator,
    current_signal,
    carrier,
    code_register::Integer,
    early_late_sample_shift::Integer,
    early_code
) where S <: AbstractGNSSSystem
    prompt_code = ((code_register & 1 << early_late_sample_shift) > 0) * 2 - 1
    late_code = ((code_register & 1 << (2 * early_late_sample_shift)) > 0) * 2 - 1
    early = get_early(correlator) + current_signal * carrier * early_code
    prompt = get_prompt(correlator) + current_signal * carrier * prompt_code
    late = get_late(correlator) + current_signal * carrier * late_code
    EarlyPromptLateCorrelator(early, prompt, late)
end

struct VeryEarlyPromptLateCorrelator{T} <: AbstractCorrelator{T}
    very_early::T
    early::T
    prompt::T
    late::T
    very_late::T
end

function VeryEarlyPromptLateCorrelator(num_ants::NumAnts{1})
    VeryEarlyPromptLateCorrelator(
        complex(0.0, 0.0),
        complex(0.0, 0.0),
        complex(0.0, 0.0),
        complex(0.0, 0.0),
        complex(0.0, 0.0)
    )
end

function VeryEarlyPromptLateCorrelator(num_ants::NumAnts{N}) where N
    VeryEarlyPromptLateCorrelator(
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64}),
        zero(SVector{N, ComplexF64})
    )
end

@inline get_very_early(correlator::VeryEarlyPromptLateCorrelator) = correlator.early
@inline get_early(correlator::VeryEarlyPromptLateCorrelator) = correlator.early
@inline get_prompt(correlator::VeryEarlyPromptLateCorrelator) = correlator.prompt
@inline get_late(correlator::VeryEarlyPromptLateCorrelator) = correlator.late
@inline get_very_late(correlator::VeryEarlyPromptLateCorrelator) = correlator.late

function zero(correlator::VeryEarlyPromptLateCorrelator{T}) where T
    EarlyPromptLateCorrelator(zero(T), zero(T), zero(T), zero(T), zero(T))
end

function filter(post_corr_filter, correlator::VeryEarlyPromptLateCorrelator)
    EarlyPromptLateCorrelator(
        post_corr_filter(get_very_early(correlator)),
        post_corr_filter(get_early(correlator)),
        post_corr_filter(get_prompt(correlator)),
        post_corr_filter(get_late(correlator)),
        post_corr_filter(get_very_late(correlator))
    )
end

function normalize(correlator::VeryEarlyPromptLateCorrelator, integrated_samples)
    EarlyPromptLateCorrelator(
        get_very_early(correlator) / integrated_samples,
        get_early(correlator) / integrated_samples,
        get_prompt(correlator) / integrated_samples,
        get_late(correlator) / integrated_samples,
        get_very_late(correlator) / integrated_samples
    )
end

# TODO: correlate and dump for very early, very late
#=
Base.@propagate_inbounds @inline function correlate_iteration(
    ::Type{S},
    correlator::VeryEarlyPromptLateCorrelator,
    current_signal,
    early_late_sample_shift,
    carrier,
    prn,
    total_code_length,
    prompt_code_phase
) where S <: AbstractGNSSSystem
    early_code_phase = prompt_code_phase + code_phase_delta * early_late_sample_shift
    early_code_phase += (early_code_phase < 0) * total_code_length
    late_code_phase = prompt_code_phase - code_phase_delta * early_late_sample_shift
    late_code_phase -= (late_code_phase >= total_code_length) * total_code_length
    early_code = get_code_unsafe(S, early_code_phase, prn)
    prompt_code = get_code_unsafe(S, code_phase, prn)
    late_code = get_code_unsafe(S, late_code_phase, prn)
    early = get_early(correlator) + current_signal * carrier * early_code
    prompt = get_prompt(correlator) + current_signal * carrier * prompt_code
    late = get_late(correlator) + current_signal * carrier * late_code
    VeryEarlyPromptLateCorrelator(early, prompt, late)
end
=#


text:
module Tachyons

using WebIO

export tachyons_css, @class_str

# const tachyons_css = Scope(imports=["/pkg/Tachyons/tachyons.min.css"])

const path = joinpath(dirname(@__FILE__), "..", "assets", "tachyons.min.css")
function __init__()
    global tachyons_css = Scope(imports=[path])
end

macro class_str(s)
    :(function (elem)
          if haskey(props(elem), :className)
              elem(className=props(elem)[:className] * " " * $(esc(s)))
          else
              elem(className=$(esc(s)))
          end
      end)
end

end # module


text:
import HorizonSideRobots
using HorizonSideRobots
sides = [Nord, Ost, Sud, West]
nx = [0, 1, 0, -1]
ny = [1, 0, -1, 0]
used = [Int[0 for i in 1:50] for i in 1:50]

mutable struct MyRobot
    r::Robot
    x::Int
    y::Int
    max_x::Int
    min_x::Int
    spawn::Int
    max_y::Int
    min_y::Int
    MyRobot() = new(r, 25, 25, 25, 25, 25, 25, 25)
end

function clear()
    for i in 1:50
        for j in 1:50
            used[i][j] = 0
        end
    end
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function mark(R)
    if (abs(R.x + R.y - R.spawn - R.spawn) % 2 == 0)
        putmarker!(R.r)
    end
end

function up(R)
    R.min_x = min(R.x, R.min_x)
    R.min_y = min(R.y, R.min_y)
    R.max_x = max(R.x, R.max_x)
    R.max_y = max(R.y, R.max_y)
end

function dfs_find_bord(R)
    for i in 1:4
        x1 = R.x + nx[i]
        y1 = R.y + ny[i]
        if (used[x1][y1] == 0 && !isborder(R.r, sides[i]))
            used[x1][y1] = 1
            move!(R.r, sides[i])
            R.x += nx[i]
            R.y += ny[i]
            up(R)
            dfs_find_bord(R)
            R.x -= nx[i]
            R.y -= ny[i]
            move!(r, inverse(sides[i]))
        end
    end
end

function dfs_move_to(R, x, y)
    for i in 1:4
        x1 = R.x + nx[i]
        y1 = R.y + ny[i]
        if (used[x1][y1] == 0 && !isborder(r, sides[i]))
            used[x1][y1] = 1
            move!(r, sides[i])
            R.x += nx[i]
            R.y += ny[i]
            dfs_move_to(R, x, y)
            if (R.x == x && R.y == y)
                return
            end
            move!(r, inverse(sides[i]))
            R.x -= nx[i]
            R.y -= ny[i]
        end
    end
end

function try_move(R, side)
    if (!isborder(R.r, side))
        move!(R.r, side)
    end
end

function cnt_horizontal(R)
    x = R.min_x
    y = R.min_y
    ans = 0
    while (y != R.max_y)
        flag = false
        x = R.min_x
        clear()
        dfs_move_to(R, x, y)
        while (x != R.max_x)
            if (isborder(R.r, Nord))
                if (!flag)
                    flag = true
                    ans += 1
                end
            else
                flag = false
            end
            clear()
            dfs_move_to(R, x + 1, y)
            x += 1
        end
        y += 1
    end
    return ans
end

function main(r)
    R = MyRobot()
    clear()
    dfs_find_bord(R)
    println(cnt_horizontal(R))
end

r = Robot(animate = false, "mirea-progs/Task20/temp.sit")
main(r)
show(r)

#=
Посчитать число всех горизонтальных прямолинейных перегородок (вертикальных - нет)
=#

max_stars_repo_name: Sergej-Kiryakov/mirea-progs

max_issues_repo_name: Sergej-Kiryakov/mirea-progs

max_forks_repo_name: Sergej-Kiryakov/mirea-progs


text:
export cmaes, pure_cmaes

function cmaes(ctype::DataType, nin::Int64, nout::Int64, fitness::Function)
    randcgp = ctype(nin, nout)
    cgpfit(g::Array{Float64}) = -fitness(ctype(mod.(g, 1.0), nin, nout))
    fitmin, best_genes = pure_cmaes(
        cgpfit, randcgp.genes, 0.5*ones(length(randcgp.genes)),
        lambda=Config.cmaes_population, stopeval=Config.cmaes_stopeval)
    -fitmin, best_genes
end

# Below is a minimal implementation of the CMA-ES based on purecmaes.m
# The following are the original comments
#
# (mu/mu_w, lambda)-CMA-ES
# CMA-ES: Evolution Strategy with Covariance Matrix Adaptation for
# nonlinear function minimization. To be used under the terms of the
# GNU General Public License (http://www.gnu.org/copyleft/gpl.html).
# Copyright: Nikolaus Hansen, 2003-09.
# e-mail: hansen[at]lri.fr
#
# This code is an excerpt from cmaes.m and implements the key parts
# of the algorithm. It is intendend to be used for READING and
# UNDERSTANDING the basic flow and all details of the CMA-ES
# *algorithm*. Use the cmaes.m code to run serious simulations: it
# is longer, but offers restarts, far better termination options,
# and supposedly quite useful output.
#
# URL: http://www.lri.fr/~hansen/purecmaes.m
# References: See end of file. Last change: October, 21, 2010

function pure_cmaes(objFun::Function, pinit::Array{Float64}, sigma::Array{Float64};
                    lambda=0,stopeval=0,stopDeltaFitness=1e-12)

#   objFun(x) = sum( (x-linspace(0,100,length(x))).^2 )
#   objFun(x) = sum( 0.1*(x[1]-1).^2 + (x[2]-2).^2 )
#    N = 3                # number of objective variables/problem dimension
#    xmean = rand(N)      # objective variables initial point
#    sigma = 0.5          # coordinate wise standard deviation (step size)

    display = 1
    min_fit = Inf
    xbest = deepcopy(pinit)

    # Input checks
    N = length(pinit)

    if(size(sigma) != size(pinit))
        error("sigma and pinit have different sizes!")
    end

    initialValue = objFun(pinit)
    if(! (typeof(initialValue) <: AbstractFloat) )
        error("objFun should return a scalar Float")
    end

    xmean = pinit

    stopfitness = 1e-10  # stop if fitness < stopfitness (minimization)

    if(stopeval==0) stopeval = 1e3*N^2 end # stop after stopeval number of function evaluations

    #########
    # Strategy parameter setting: Selection
    if(lambda==0)
        lambda = convert(Integer,4+floor(3*log(N)))  # population size, offspring number
    end
    mu = lambda/2                   # number of parents/points for recombination
    weights = log.(mu+1/2).-log.(1:mu) # muXone array for weighted recombination
    mu = convert(Integer,floor(mu))
    weights = weights/sum(weights)     # normalize recombination weights array
    mueff=sum(weights)^2/sum(weights.^2) # variance-effectiveness of sum w_i x_i

    #########
    # Strategy parameter setting: Adaptation

    cc = 4.0/(N+4.0) # time constant for cumulation for C
    cs = (mueff+2.0) / (N+mueff+3.0)  # t-const for cumulation for sigma control

    ccov1 = 2.0 / ((N+1.3)^2+mueff)    # learning rate for rank-one update of C
    ccovmu = min(1-ccov1, 2 * (mueff-2+1/mueff) / ((N+2)^2+mueff) )   # and for rank-mu update

    damps = 1.0 + 2.0*max(0.0, sqrt((mueff-1)/(N+1))-1) + cs # damping for sigma, usually close to 1

    #########
    # Initialize dynamic (internal) strategy parameters and constants

    pc = zeros(N); ps = zeros(N)   # evolution paths for C and sigma

    diagD = sigma/maximum(sigma);      # diagonal matrix D defines the scaling
    diagC = diagD.^2;

    sigma = maximum(sigma)

    B = eye(N,N)                        # B defines the coordinate system

    BD = B.*repmat(diagD',N,1)          # B*D for speed up only
    C = diagm(diagC);                   # covariance matrix == BD*(BD)'

    chiN=N^0.5*(1-1/(4*N)+1/(21*N^2))  # expectation of  ||N(0,I)|| == norm(randn(N,1))

    #init a few things
    arx = zeros(N,lambda)
    arz = zeros(N,lambda)
    arfitness = zeros(lambda)
    arindex = zeros(lambda)

    previousFitness = 0

    iter = 0

    @printf("%i-%i CMA-ES\n",lambda,mu);

    # -------------------- Generation Loop --------------------------------
    counteval = 0  # the next 40 lines contain the 20 lines of interesting code
    while counteval < stopeval

        iter = iter+1

        # Generate and evaluate lambda offspring
        for k=1:lambda
            #      arx[:,k] = xmean + sigma * B * (D .* randn(N,1)) # m + sig * Normal(0,C)

            arz[:,k] = randn(N,1); # resample
            arx[:,k] = xmean + sigma * (BD * arz[:,k]);

            arfitness[k] = objFun( arx[:,k] ) # objective function call
            counteval = counteval+1
            if arfitness[k] < min_fit
                min_fit = arfitness[k]
                xbest = deepcopy(arx[:,k])
                Logging.info(@sprintf("CMAES: %d %0.4f %s", counteval, -min_fit,
                                      string(mod.(xbest, 1.0))))
            end
        end

        # Sort by fitness and compute weighted mean into xmean
        arindex = sortperm(arfitness)
        arfitness = arfitness[arindex]  # minimization

        # Calculate new xmean, this is selection and recombination
        xold = xmean; # for speed up of Eq. (2) and (3)
        xmean = arx[:,arindex[1:mu]]*weights
        zmean = arz[:,arindex[1:mu]]*weights # ==D^-1*B'*(xmean-xold)/sigma

        # Cumulation: Update evolution paths
        ps = (1-cs)*ps + sqrt(cs*(2-cs)*mueff) * (B*zmean)          # Eq. (4)
        hsig = norm(ps)/sqrt(1-(1-cs)^(2*iter))/chiN < 1.4 + 2/(N+1)
        pc = (1-cc)*pc + hsig*(sqrt(cc*(2-cc)*mueff)/sigma) * (xmean-xold)

        # Adapt covariance matrix C
        # artmp = (1/sigma) * (arx[:,arindex[1:mu]] - repmat(xold,1,mu))  # mu difference vectors

        if(ccov1+ccovmu>0)

          C =( (1-ccov1-ccovmu+(1-hsig)*ccov1*cc*(2-cc)) * C     # regard old matrix
              + ccov1 * pc*pc'                             # plus rank one update
              + ccovmu                                     # plus rank mu update
              * sigma^-2 * (arx[:,arindex[1:mu]]-repmat(xold,1,mu))
              * (repmat(weights,1,N) .* (arx[:,arindex[1:mu]]-repmat(xold,1,mu))')  )

        end

        # Adapt step size sigma
        sigma = sigma * exp((norm(ps)/chiN - 1)*cs/damps)  #Eq. (5)

        # Update B and D from C
        # if counteval - eigeneval > lambda/(ccov1+ccovmu)/N/10  # to achieve O(N^2)
        if( mod(iter, 1/(ccov1+ccovmu)/N/10) < 1 && ccov1+ccovmu>0 )

            C = triu(C) + triu(C,1)'    # enforce symmetry
            (tmp,B) = eig(C)            # eigen decomposition, B==normalized eigenvectors

            diagD = sqrt.(tmp); # D contains standard deviations now

            diagC = diag(C);

            BD = B.*repmat(diagD',N,1); # O(n^2)
        end

        #Stop conditions:
        # Break, if fitness is good enough or condition exceeds 1e14, better termination methods are advisable
        # if arfitness[1] <= stopfitness || maximum(diagD) > 1e7 * minimum(diagD)
        #   break
        # end

        #break if fitness doesn't change much
        # if(iter > 1)
        #     if( abs(arfitness[1]-previousFitness) < stopDeltaFitness)
        #         break
        #     end
        # end
        previousFitness = arfitness[1]

    end # while, end generation loop

    xmin = arx[:, arindex[1]];

    function corcov(C::AbstractMatrix)
        #should check if C is positive semidefinite

        sigma = sqrt(diag(C))
        return C ./ (sigma*sigma')
    end

    if( N<10 )
        print("\nCorrelation matrix:\n")
        println(corcov(C))
    end

    return min_fit, mod.(xbest, 1.0)

end

# ---------------------------------------------------------------
### REFERENCES
#
# Hansen, N. and S. Kern (2004). Evaluating the CMA Evolution
# Strategy on Multimodal Test Functions.  Eighth International
# Conference on Parallel Problem Solving from Nature PPSN VIII,
# Proceedings, pp. 282-291, Berlin: Springer.
# (http://www.bionik.tu-berlin.de/user/niko/ppsn2004hansenkern.pdf)
#
# Further references:
# Hansen, N. and A. Ostermeier (2001). Completely Derandomized
# Self-Adaptation in Evolution Strategies. Evolutionary Computation,
# 9(2), pp. 159-195.
# (http://www.bionik.tu-berlin.de/user/niko/cmaartic.pdf).
#
# Hansen, N., S.D. Mueller and P. Koumoutsakos (2003). Reducing the
# Time Complexity of the Derandomized Evolution Strategy with
# Covariance Matrix Adaptation (CMA-ES). Evolutionary Computation,
# 11(1).  (http://mitpress.mit.edu/journals/pdf/evco_11_1_1_0.pdf).


text:
# Estimates the basic Nerlove Cobb-Douglas model
using Econometrics, DelimitedFiles
function main()
data = readdlm("../Data/nerlove.data")
data = data[:,2:6]
data = log.(data)
n = size(data,1)
y = data[:,1]
x = data[:,2:end]
x = [ones(n,1) x]
names = ["constant", "output", "labor", "fuel", "capital"]
b, junk, junk, junk = ols(y, x, names=names, vc="ols")
return
end
main()

max_stars_repo_name: mcreel/Econometrics

max_issues_repo_name: mcreel/Econometrics

max_forks_repo_name: mcreel/Econometrics


text:
module GeoArrays

using ArchGDAL
using CoordinateTransformations
using StaticArrays

include("geoarray.jl")
include("geoutils.jl")
include("io.jl")
include("utils.jl")
include("interpolate.jl")
include("crs.jl")
include("operations.jl")
include("plot.jl")

export GeoArray

export coords
export centercoords
export indices

export bbox
export bboxes

export compose!

export interpolate!

export epsg!

export -,+,*,/

end


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

const _VT = NamedTuple{NAMES,NTuple{N,Int32}} where {NAMES,N}

throw_wrong_value(kind::AbstractString, val) = throw(ArgumentError("Invalid $kind $val"))

check_value(val, valid::_VT) = throw_wrong_value(Val(valid), val)
check_value(val::AbstractString, valid::_VT) = check_value(Symbol(lowercase(val)), valid)
check_value(val::Symbol, valid::_VT) = haskey(valid, val) ? val : throw_wrong_value(Val(valid), val)
check_value(val::Integer, valid::_VT) = val ∈ valid ? val : throw_wrong_value(Val(valid), val)

val_to_symbol(val, valid::_VT) = check_value(val, valid)
val_to_symbol(val::Integer, valid::_VT) = keys(valid)[indexin([check_value(val, valid)], [valid...])[1]]

val_to_int(val, valid::_VT) = valid[check_value(val, valid)]
val_to_int(val::Integer, valid::_VT) = check_value(val, valid)

const access_mode = (;
    readonly=Int32(1),
    create=Int32(2),
    overwrite=Int32(3),
    update=Int32(4),
    shared=Int32(5),
    write=Int32(6),
    direct_write=Int32(7)
)
throw_wrong_value(::Val{access_mode}, val) = throw_wrong_value("access mode", val)
check_mode(val) = check_value(val, access_mode)

# Chli and FameDatabase each depends on the other.
# We need AbstractFameDatabase only to make that work.
abstract type AbstractFameDatabase end;

"""
    struct Chli ⋯ end

The internal state of FAME.jl. There is no reason for users to play with this.
"""
struct Chli{AFD<:AbstractFameDatabase}
    lib::Union{Nothing,Ptr{Nothing}}
    opendb::Vector{AFD}
    workdb::Ref{AFD}
    Chli{AFD}(lib, od=AFD[], wd=Ref{AFD}()) where {AFD} = new(lib, od, wd)
end

"""
    global chli

The internal state of FAME.jl. There is no reason for users to play with this.
"""
global chli = Chli{AbstractFameDatabase}(nothing)
@inline getsym(sym::Symbol) = chli.lib !== nothing ? Libdl.dlsym(chli.lib, sym) : error("FAME not initialized.")

function push_db(db::AbstractFameDatabase)
    global chli
    push!(chli.opendb, db)
    return db
end

export FameDatabase
"""
    mutable struct FameDatabase … end

Fame database. See also [`opendb`](@ref), [`closedb!`](@ref), [`postdb`](@ref),
[`listdb`](@ref).
"""
mutable struct FameDatabase <: AbstractFameDatabase
    key::Int32
    name::String
    mode::Symbol
    function FameDatabase(key::Integer, name::AbstractString="", mode=:readonly)
        push_db(new(key, name, val_to_symbol(mode, access_mode)))
    end
end

Base.isopen(db::FameDatabase) = db.key >= 0
Base.isequal(l::FameDatabase, r::FameDatabase) = l.key == r.key

function Base.show(io::IO, x::FameDatabase)
    if x.key >= 0
        print(io, "$(x.name): open in $(x.mode) mode")
    else
        print(io, "$(x.name): closed")
    end
end

const fame_class = (;
    # refer to hli.h for the values
    series=Int32(1),
    scalar=Int32(2),
    formula=Int32(3),
    glname=Int32(5),
    glformula=Int32(6)
)
throw_wrong_value(::Val{fame_class}, val) = throw_wrong_value("class", val)
check_class(val) = check_value(val, fame_class)

#/* FAME BASIS Attribute Settings  */
const fame_basis = (;
    undefined=Int32(0),
    daily=Int32(1),
    business=Int32(2)
)
throw_wrong_value(::Val{fame_basis}, val) = throw_wrong_value("basis", val)
check_basis(val) = check_value(val, fame_basis)

#/* FAME OBSERVED Attribute Settings  */
const fame_observed = (;
    undefined=Int32(0),
    beginning=Int32(1),
    ending=Int32(2),
    averaged=Int32(3),
    summed=Int32(4),
    annualized=Int32(5),
    formula=Int32(6),
    high=Int32(7),
    low=Int32(8)
)
throw_wrong_value(::Val{fame_observed}, val) = throw_wrong_value("observed", val)
check_observed(val) = check_value(val, fame_observed)

#/* Relation */
const fame_relation = (;
    before=Int32(1),
    after=Int32(2),
    contains=Int32(3)
)
throw_wrong_value(::Val{fame_relation}, val) = throw_wrong_value("relation", val)
check_relation(val) = check_value(val, fame_relation)

const fame_freq = (;
    undefined=Int32(0),               #define HUNDFX    0	/* Undefined			*/
    daily=Int32(8),                   #define HDAILY    8	/* DAILY			*/
    business=Int32(9),                #define HBUSNS    9	/* BUSINESS			*/
    weekly_sunday=Int32(16),          #define HWKSUN   16	/* WEEKLY (SUNDAY)		*/
    weekly_monday=Int32(17),          #define HWKMON   17	/* WEEKLY (MONDAY)		*/
    weekly_tuesday=Int32(18),         #define HWKTUE   18	/* WEEKLY (TUESDAY)		*/
    weekly_wednesday=Int32(19),       #define HWKWED   19	/* WEEKLY (WEDNESDAY)		*/
    weekly_thursday=Int32(20),        #define HWKTHU   20	/* WEEKLY (THURSDAY)		*/
    weekly_friday=Int32(21),          #define HWKFRI   21	/* WEEKLY (FRIDAY)		*/
    weekly_saturday=Int32(22),        #define HWKSAT   22	/* WEEKLY (SATURDAY)		*/
    tenday=Int32(32),                 #define HTENDA   32	/* TENDAY			*/
    biweekly_asunday=Int32(64),       #define HWASUN   64	/* BIWEEKLY (ASUNDAY)		*/
    biweekly_amonday=Int32(65),       #define HWAMON   65	/* BIWEEKLY (AMONDAY)		*/
    biweekly_atuesday=Int32(66),      #define HWATUE   66	/* BIWEEKLY (ATUESDAY		*/
    biweekly_awednesday=Int32(67),    #define HWAWED   67	/* BIWEEKLY (AWEDNESDAY)	*/
    biweekly_athursday=Int32(68),     #define HWATHU   68	/* BIWEEKLY (ATHURSDAY)		*/
    biweekly_afriday=Int32(69),       #define HWAFRI   69	/* BIWEEKLY (AFRIDAY)		*/
    biweekly_asaturday=Int32(70),     #define HWASAT   70	/* BIWEEKLY (ASATURDAY)		*/
    biweekly_bsunday=Int32(71),       #define HWBSUN   71	/* BIWEEKLY (BSUNDAY)		*/
    biweekly_bmonday=Int32(72),       #define HWBMON   72	/* BIWEEKLY (BMONDAY)		*/
    biweekly_btuesday=Int32(73),      #define HWBTUE   73	/* BIWEEKLY (BTUESDAY)		*/
    biweekly_bwednesday=Int32(74),    #define HWBWED   74	/* BIWEEKLY (BWEDNESDAY)	*/
    biweekly_bthursday=Int32(75),     #define HWBTHU   75	/* BIWEEKLY (BTHURSDAY)		*/
    biweekly_bfriday=Int32(76),       #define HWBFRI   76	/* BIWEEKLY (BFRIDAY)		*/
    biweekly_bsaturday=Int32(77),     #define HWBSAT   77	/* BIWEEKLY (BSATURDAY)		*/
    twicemonthly=Int32(128),          #define HTWICM  128	/* TWICEMONTHLY			*/
    monthly=Int32(129),               #define HMONTH  129	/* MONTHLY			*/
    bimonthly_november=Int32(144),    #define HBMNOV  144	/* BIMONTHLY (NOVEMBER)		*/
    bimonthly_december=Int32(145),    #define HBIMON  145	/* BIMONTHLY (DECEMBER)		*/
    quarterly_october=Int32(160),     #define HQTOCT  160	/* QUARTERLY (OCTOBER)		*/
    quarterly_november=Int32(161),    #define HQTNOV  161	/* QUARTERLY (NOVEMBER)		*/
    quarterly_december=Int32(162),    #define HQTDEC  162	/* QUARTERLY (DECEMBER)		*/
    annual_january=Int32(192),        #define HANJAN  192	/* ANNUAL (JANUARY)		*/
    annual_february=Int32(193),       #define HANFEB  193	/* ANNUAL (FEBRUARY)		*/
    annual_march=Int32(194),          #define HANMAR  194	/* ANNUAL (MARCH)		*/
    annual_april=Int32(195),          #define HANAPR  195	/* ANNUAL (APRIL)		*/
    annual_may=Int32(196),            #define HANMAY  196	/* ANNUAL (MAY)			*/
    annual_june=Int32(197),           #define HANJUN  197	/* ANNUAL (JUNE)		*/
    annual_july=Int32(198),           #define HANJUL  198	/* ANNUAL (JULY)		*/
    annual_august=Int32(199),         #define HANAUG  199	/* ANNUAL (AUGUST)		*/
    annual_september=Int32(200),      #define HANSEP  200	/* ANNUAL (SEPTEMBER)		*/
    annual_october=Int32(201),        #define HANOCT  201	/* ANNUAL (OCTOBER)		*/
    annual_november=Int32(202),       #define HANNOV  202	/* ANNUAL (NOVEMBER)		*/
    annual_december=Int32(203),       #define HANDEC  203	/* ANNUAL (DECEMBER)		*/
    semiannual_july=Int32(204),       #define HSMJUL  204	/* SEMIANNUAL (JULY)		*/
    semiannual_august=Int32(205),     #define HSMAUG  205	/* SEMIANNUAL (AUGUST)		*/
    semiannual_september=Int32(206),  #define HSMSEP  206	/* SEMIANNUAL (SEPTEMBER)	*/
    semiannual_october=Int32(207),    #define HSMOCT  207	/* SEMIANNUAL (OCTOBER)		*/
    semiannual_november=Int32(208),   #define HSMNOV  208	/* SEMIANNUAL (NOVEMBER)	*/
    semiannual_december=Int32(209),   #define HSMDEC  209	/* SEMIANNUAL (DECEMBER)	*/
    ypp=Int32(224),                   #define HAYPP   224	/* YPP				*/
    ppy=Int32(225),                   #define HAPPY   225	/* PPY				*/
    secondly=Int32(226),              #define HSEC    226	/* SECONDLY			*/
    minutely=Int32(227),              #define HMIN    227	/* MINUTELY			*/
    hourly=Int32(228),                #define HHOUR   228	/* HOURLY			*/
    millisecondly=Int32(229),         #define HMSEC   229	/* MILLISECONDLY		*/
    case=Int32(232),                  #define HCASEX  232	/* CASE				*/
    weekly_pattern=Int32(233)        #define HWEEK_PATTERN 233 /* generic weekly pattern     */
)
throw_wrong_value(::Val{fame_freq}, val) = throw_wrong_value("frequency", val)
check_freq(val) = check_value(val, fame_freq)

iscase(x::Symbol) = x === :case
iscase(x::Integer) = x === fame_freq.case

const fame_type = (;
    undefined=Int32(0),   #define HUNDFT    0	/* Undefined	*/
    numeric=Int32(1),     #define HNUMRC    1	/* NUMERIC	*/
    namelist=Int32(2),    #define HNAMEL    2	/* NAMELIST	*/
    boolean=Int32(3),     #define HBOOLN    3	/* BOOLEAN	*/
    string=Int32(4),      #define HSTRNG    4	/* STRING	*/
    precision=Int32(5),   #define HPRECN    5	/* PRECISION	*/
    date=Int32(6)        #define HDATE     6	/* General DATE	*/
)
throw_wrong_value(::Val{fame_type}, val) = throw_wrong_value("type", val)
check_type(val) = check_value(val, fame_type)

"""
    FameIndex

Integer 64-bit type that represents the internal index FAME uses to access
elements of series.
"""
const FameIndex = Int64

"""
    FameDate

64-bit integer type that FAME uses to encode dates of different frequencies.
"""
const FameDate = FameIndex

## More types are defined in FameObjects. We need stuff from ChliLibrary.jl in
## order to define those, so we can't include them here.


text:
abstract type AbstractPolicy end

abstract type AbstractValuePolicy <: AbstractPolicy end

struct ϵGreedy <: AbstractValuePolicy
    ϵ::Float64
end

function sample(ap::ϵGreedy, values, rng)
    if rand(rng) > ap.ϵ
        return findmax(values)[2]
    else
        return rand(rng, 1:length(values))
    end
end

function get_prob(ap::ϵGreedy, values, action)
    if action == findmax(values)[2]
        return 1 - ap.ϵ + (ap.ϵ / length(values))
    else
        return ap.ϵ / length(values)
    end
end


text:
using Gadfly, DataFrames

i = indexin([Float64], [eltypes...])[1]
tconstr = tconstr[:,:,i]
teval = teval[:,:,i]

# Convert to dataframes, and for teval use points/second
dim = ones(Int,size(tconstr,1))*(1:size(tconstr,2))'
namelist = collect(take(cycle(testnames), length(dim)))
ordlist =  collect(take(cycle(testord), length(dim)))
dfconstr = DataFrame(t = vec(tconstr), name = namelist, ord = ordlist, dim = vec(dim))
dfeval   = DataFrame(rate = vec(10^6./teval), name = namelist, ord = ordlist, dim = vec(dim))

flagBspline = (dfconstr[:name] .== "IBSpline") | (dfconstr[:name] .== "Grid")
dfBconstr = dfconstr[flagBspline, :]
dfBeval   = dfeval[flagBspline, :]
dfGconstr = dfconstr[!flagBspline, :]
dfGeval   = dfeval[!flagBspline, :]

yscale = Scale.y_continuous

# Plot the B-spline methods
plotsize = (6inch, 8inch)
set_default_plot_size(plotsize...)
pBc = plot(dfBconstr, xgroup="dim", ygroup="ord", x="name", y="t", Geom.subplot_grid(Geom.point), yscale, Guide.ylabel("Time (s) by order"))
draw(PNG("constructionB.png", plotsize...), pBc)

pBr = plot(dfBeval, xgroup="dim", ygroup="ord", x="name", y="rate", Geom.subplot_grid(Geom.point), yscale, Guide.ylabel("Throughput (pts/s) by order"))
draw(PNG("rateB.png", plotsize...), pBr)

# Plot the gridded methods
plotsize = (10inch, 8inch)
set_default_plot_size(plotsize...)

pGc = plot(dfGconstr, xgroup="dim", ygroup="ord", x="name", y="t", Geom.subplot_grid(Geom.point), yscale, Guide.ylabel("Time (s) by order"))
draw(PNG("constructionG.png", plotsize...), pGc)

pGr = plot(dfGeval, xgroup="dim", ygroup="ord", x="name", y="rate", Geom.subplot_grid(Geom.point), Scale.y_sqrt, Guide.ylabel("Throughput (pts/s) by order"))
draw(PNG("rateG.png", plotsize...), pGr)


text:
using GeoPhyInv
using Statistics
#jl @init_parallel_stencil(2, false, Float32, 4)
#md using Plots; gr();
#md using ColorSchemes

# ### Medium
medium = Medium(:acou_homo2D, 5); # load a simple homogeneous acoustic medium from the gallery
update!(medium, [:vp, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)

# ### AGeom
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)

# ### SrcWav
tgrid = range(0.0, stop = 2.0, length = 2500); # generate a time grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # ricker wavelet
srcwav = SrcWav(tgrid, ageom, [:p]);
update!(srcwav, [:p], wav);

# ### Plotting the `vp` model
#md p1=heatmap(medium, :vp, seriescolor=cgrad(:roma))
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1)

# ## Acoustic
# Now we have all the required variables to create `SeisForwExpt` object and
# prepare the forward modeling.
# While creating, we switched the `snaps_flag` on, and instruct recording field at
# `tsnaps`.
# Once the `Expt` object is created, do the modeling without *additional* any
# memory allocations using `update!`

tsnaps=tgrid[1:div(length(tgrid),20):end] # store 20 snapshots
pa = SeisForwExpt(
    FdtdAcoustic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :p,
    tsnaps = tsnaps,
    tgrid = tgrid,
    rfields = [:p],
    verbose = true,
);

# ### Update `pa` to perform modelling; prints time and allocations
@time update!(pa);

# ### Extracting snaps
snaps1 = pa[:snaps, 1]; # extracting snaps of the first supersource
snaps2 = pa[:snaps, 2]; # second supersource

# ### Plotting pressure snapshots after acoustic simulation
#md function plot_snapshots(name)
#md     clip_perc = 90
#md     pmax = maximum([maximum(pa[:snaps][it]) for it = 1:length(tsnaps)])
#md     pmax -= pmax * 0.01 * clip_perc
#md
#md     anim = @animate for it = 1:length(tsnaps)
#md         plot(
#md             [
#md                 (
#md                     f = pa[:snaps, is][it];
#md                     heatmap(
#md                         f,
#md                         aspect_ratio = 1,
#md                         xlims = (1, size(f, 2)),
#md                         ylims = (1, size(f, 1)),
#md                         yflip = true,
#md                         seriescolor = cgrad(:seismic),
#md                         clims = (-pmax, pmax),
#md                         legend = false,
#md                     )
#md                 ) for is = 1:2
#md             ]...,
#md             title = string("snapshot at: ", round(tsnaps[it], digits = 5), " s"),
#md         )
#md     end
#md     return gif(anim, string(name,".gif"), fps = 1)
#md end
#md plot_snapshots("acoustic_snaps")

# ## Elastic
# Similarly, we can do the elastic wave-equation modeling
medium = Medium(:elastic_homo2D, 5); # load a simple homogeneous elastic medium from the gallery
update!(medium, [:vp, :vs, :rho], randn_perc = 5); # add some random noise to the medium
println(medium)
ageom = AGeom(medium.mgrid, :xwell, SSrcs(2)); # load a simple acquisition using `mgrid` of the medium
println(ageom)
# Add source term on `:vz` grid
srcwav = SrcWav(tgrid, ageom, [:vz]);
update!(srcwav, [:vz], wav);
# Perform modelling
pa = SeisForwExpt(
    FdtdElastic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    snaps_field = :vz,
    tsnaps = tsnaps,
    rfields = [:vz],
    tgrid = tgrid,
    verbose = true,
)
@time update!(pa);

snaps = pa[:snaps, 1];
snaps = pa[:snaps, 2];

# ### Plotting is vs model
#md p1=heatmap(medium, :vs, seriescolor=cgrad(:roma))
#md scatter!(ageom, SSrcs())
#md scatter!(ageom, Recs())
#md plot(p1)

# ### Plotting snapshots after elastic simulation
#md plot_snapshots("elastic_snaps") # P and S waves!


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

using Markdown
using InteractiveUtils

# ╔═╡ 4aebe8f8-502d-11eb-1925-694f35421317
md"""
# Recursion

Recursive function calls itself with a smaller instance of the same problem and ends with a base case. Three required ingredients:
1. Calls itself
2. with a smaller input parameter
3. end with a base case

Be careful with Infinite Recursion, which can occur when
- parameter don't decrease/increase towards the base case criteria
- base case cannot be reached

"""

# ╔═╡ 8c30ca4a-502d-11eb-2c0b-398047ab685f
function addrecursive(x)
	x >= 1 && error("x < 1")  # Defense against infinite recursion
    if x == 1  # if x is less than 1 this becomes an infinite recursion
        return 1
    else
        return x + addrecursive(x - 1)
    end
end

# ╔═╡ 0d1581d2-502e-11eb-2b2b-e3cc836fc2b0
# The Call Stack of Recursion

# ╔═╡ 130ca5e8-502e-11eb-3f82-5ff3930f7421
# Recurrence Relations

# ╔═╡ 1b0fe3ae-502e-11eb-2d34-c5a885ec3791
md"""
Applications
- Factorial
- Reverse a string
- Palindrome Checking
- List permutation
- String Permutations
- keypad Combinations
- deep reverse ?
- Tower of Hanoi
- Staircase ?

"""

# ╔═╡ Cell order:
# ╟─4aebe8f8-502d-11eb-1925-694f35421317
# ╠═8c30ca4a-502d-11eb-2c0b-398047ab685f
# ╠═0d1581d2-502e-11eb-2b2b-e3cc836fc2b0
# ╠═130ca5e8-502e-11eb-3f82-5ff3930f7421
# ╠═1b0fe3ae-502e-11eb-2d34-c5a885ec3791


text:
using ValueOrientedRiskManagementInsurance
using DataFrames
using Base.Test

### TESTS
println("start test ChainLadder_Test.jl...")
# We use the example from Mack's 1993-paper
# Cumulated claim amounts
cum_triangle =
  convert(Array{Real,2},
          readtable("test/ChainLadderMack_Input.csv", header = false))

triangle =cum2claims(cum_triangle)
@test_approx_eq(cum_triangle, claims2cum(triangle))

mack = Mack(cum_triangle, cum=true)

# Test: development factor
@test_approx_eq(
  [3.491, 1.747, 1.457, 1.174, 1.104, 1.086, 1.054, 1.077, 1.018],  round(mack.f,3))
# Test: future claims
@test_approx_eq(sum(mack.futureclaims), mack.tot_res)
# Test: standard error
@test_approx_eq(
  [0.80, 0.26, 0.19, 0.27, 0.29, 0.26, 0.22, 0.23, 0.29],
  round(√(mack.mse) ./ mack.res,2)[2:end])
# Test: total standard error
@test_approx_eq(0.13, round(√(mack.tot_mse)/mack.tot_res, 2))
println("...end test ChainLadder_Test.jl")


text:
using PyPlot

using mvIMPACT

DMR.Init()
DevCnt = DMR.getDeviceCount()
println("$(DevCnt) devices detected")

hDev = DMR.GetDeviceNr(0)

info = DMR.GetDeviceInfoEx(hDev, DMR.dmditDeviceInfoStructure)
println(info.serial)
println(info.family)
println(info.product)
println(info.firmwareVersion)

println("Opening device $hDev")
hDrv = DMR.OpenDevice(hDev)

hList = DMR.FindList(hDrv, DMR.dmltSetting)
println("hList of type: ", OBJ.as_str(OBJ.GetType(hList)))
println(OBJ.GetContentDesc(hList))
hGain = OBJ.FindProperty(hList, "Gain_dB")
println(OBJ.GetName(hGain), " of type ", OBJ.as_str(OBJ.GetType(hGain)),
        " has value ", OBJ.ReadProperty(hGain))
hExposure = OBJ.FindProperty(hList, "Expose_us")
println(OBJ.GetName(hExposure), " of type ", OBJ.as_str(OBJ.GetType(hExposure)),
        " has value ", OBJ.ReadProperty(hExposure))
hPixelFormat = OBJ.FindProperty(hList, "PixelFormat")
println(OBJ.GetName(hPixelFormat), " of type ", OBJ.as_str(OBJ.GetType(hPixelFormat)),
        " has value ", OBJ.ReadProperty(hPixelFormat))

hW = OBJ.FindProperty(hList, "W")
println(OBJ.GetName(hW), " of type ", OBJ.as_str(OBJ.GetType(hW)),
        " has value ", OBJ.ReadProperty(hW))
hH = OBJ.FindProperty(hList, "H")
println(OBJ.GetName(hH), " of type ", OBJ.as_str(OBJ.GetType(hH)),
        " has value ", OBJ.ReadProperty(hH))

println("Gain dict size: ", OBJ.GetDictSize(hGain))
println("Exposure dict size: ", OBJ.GetDictSize(hExposure))
println("PixelFormat dict size: ", OBJ.GetDictSize(hPixelFormat))

OBJ.WriteProperty(hGain, 0)
OBJ.WriteProperty(hExposure, 1000)
println(OBJ.GetIDictEntries(hPixelFormat))
OBJ.WriteProperty(hPixelFormat, "Mono10")

println("\t", OBJ.ReadProperty(hGain), "\n",
        "\t", OBJ.ReadProperty(hExposure), "\n",
        "\t", OBJ.ReadProperty(hPixelFormat))

l_imgs = []

@time for i in 1:10
    img = DMR.GetImage(hDrv)
    push!(l_imgs, img)
end

figure()
imshow(l_imgs[1])
colorbar()
figure()
imshow(l_imgs[end])
colorbar()

show()

println("Closing device $hDev with driver $hDrv")
DMR.CloseDevice(hDrv, hDev)
DMR.Close()


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2360bd70-1a15-11eb-3b28-75fdf85c0706
using NeXLDatabase

# ╔═╡ 290a7a40-1a15-11eb-338e-1d2c5e513ab3
using Gadfly, DataFrames, PlutoUI

# ╔═╡ bede52a0-1a27-11eb-1206-2bd005ea4d01
using SQLite

# ╔═╡ 30233060-1a15-11eb-2d1f-8b7537f66d45
db=openNeXLDatabase("C:\\Users\\nritchie\\Documents\\DrWatson\\K-Ratio Project\\data\\exp_pro\\kratio.db")

# ╔═╡ 7cfccdb0-1a15-11eb-37b7-4f1f1aff8f2f
begin
	labs = findall(db, DBLaboratory)
	md"""
	#### Laboratories
	Select a lab: $(@bind lab Select([ repr(i) => repr(lab) for (i, lab) in enumerate(labs) ]))
	"""
end

# ╔═╡ 3bded950-1a23-11eb-3d86-c1d125393e99
begin
	thelab = labs[parse(Int,lab)]
	members = findall(db, DBPerson, thelab)
	md"""
	#### Lab Members
	Select a lab member: $(@bind lmi Select([ repr(i) => person.name for (i, person) in enumerate(members) ]))
	"""
end

# ╔═╡ c5eb0000-1a15-11eb-0f46-99d1ca2092f6
begin
	theperson = members[parse(Int,lmi)]
	projects = findall(db, DBProject, theperson)
	md"""
	#### Projects
	Select a project: $(@bind pidx Select([ repr(i) => repr(proj) for (i, proj) in enumerate(projects) ]))
	"""
end

# ╔═╡ 2d6bb0b0-1a45-11eb-2233-c92381d65e13

# ╔═╡ c602c890-1a27-11eb-19fc-c79cbf0026c1
begin
	theproject = projects[parse(Int,pidx)]
	fspecs = findall(db, DBCampaign, project=theproject)
	asa(DataFrame, fspecs)
end

# ╔═╡ Cell order:
# ╟─2360bd70-1a15-11eb-3b28-75fdf85c0706
# ╟─290a7a40-1a15-11eb-338e-1d2c5e513ab3
# ╟─30233060-1a15-11eb-2d1f-8b7537f66d45
# ╟─7cfccdb0-1a15-11eb-37b7-4f1f1aff8f2f
# ╟─3bded950-1a23-11eb-3d86-c1d125393e99
# ╟─c5eb0000-1a15-11eb-0f46-99d1ca2092f6
# ╟─bede52a0-1a27-11eb-1206-2bd005ea4d01
# ╟─2d6bb0b0-1a45-11eb-2233-c92381d65e13
# ╟─c602c890-1a27-11eb-19fc-c79cbf0026c1


text:
include("diff.jl")
include("node.jl")


text:
# Compute graph statistics
export sample, edge_density, degree_distribution

## Graph statistics on concrete graphs

function sample(A::Array{Float64,2}, nsamples::Int)
    @assert size(A,1) == size(A,2)
    n = size(A,1)
    indices = (1:n)[int(ceil(n*rand(nsamples)))]
    return A[indices, indices]
end

function edge_density(A::AbstractArray)
    return mean(A)
end

function degree_distribution(A::AbstractArray)
    # sorted normalized degrees
    sort(vec(sum(A,1)/size(A,1)))
end

## Graph statistics on graphons

degree_distribution(g::Graphon, k::Int) = degree_distribution(sample(g, k))
edge_density(g::Graphon, k::Int) = edge_density(sample(g, k))


text:
#!/usr/bin/env julia
using Printf
using Distributed
using TimerOutputs
using CUDAdrv
using CUDAnative
using CuArrays

function ζ_2_sqr_gpu_range(N)
	MAX_THREADS_PER_BLOCK = 1024
    Nthr    = min(N, MAX_THREADS_PER_BLOCK)
    Nblocks = ceil.(Int, N / Nthr)

    psum_d = CuArrays.fill(0.0::Float64, Nthr)
	p  = zeros(Float64,Nthr)
	@printf("Nblocks = %i\nNthr = %i\n",Nblocks,Nthr)

    @cuda blocks=Nblocks threads=Nthr  innergpu_1!(psum_d)
    copyto!(p,psum_d)
    s = sum(p)
    s
end

function ζ_gen_gpu_range(N)
	MAX_THREADS_PER_BLOCK = 1024
    Nthr    = min(N, MAX_THREADS_PER_BLOCK)
    Nblocks = ceil.(Int, N / Nthr)

    psum_d = CuArrays.fill(0.0::Float64, Nthr)
	p  = zeros(Float64,Nthr)
	@printf("Nblocks = %i\nNthr = %i\n",Nblocks,Nthr)

    @cuda blocks=Nblocks threads=Nthr  innergpu_2!(psum_d,2)
    copyto!(p,psum_d)
    s = sum(p)
    s
end

function innergpu_1!(ps)
	ndx     = threadIdx().x
	i = (blockIdx().x-1) * blockDim().x + threadIdx().x
	r = 1.0/i
	@atomic ps[ndx] += r*r
	return nothing
end

function innergpu_2!(ps,x)
	ndx     = threadIdx().x
	i = (blockIdx().x-1) * blockDim().x + threadIdx().x
	r = 1.0/i
	@atomic ps[ndx] += CUDAnative.pow.(r,x)
	return nothing
end
to = TimerOutput()

N = 16000000000
st = 1
s =  zeros(Float64,32)
pd = zeros(Float64,32)
ndx = 1
#ENV["JULIA_DEBUG"] = "CUDAnative"

# first bit of code to handle the compilation stage, otherwise you
# see this in the time output

s[1] = ζ_2_sqr_gpu_range(N)  #
s[2] = ζ_gen_gpu_range(N)

for i=1:2
  @printf("s[%i] = %f\n",i,s[i])
end

@timeit to "ζ_2_sqr_gpu_range"     pd[ndx] = ζ_2_sqr_gpu_range(N)
ndx+=1

@timeit to "ζ_gen_gpu_range"     pd[ndx] = ζ_gen_gpu_range(N)

for i=1:2
  @printf("pd[%i] = %f\n",i,pd[i])
end

show(to)
@printf("\n");

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
using Colors: RGB, N0f8, Gray
using FastAI
using FastAI: ParamGroups, IndexGrouper, getgroup, DiscriminativeLRs, decay_optim
using FastAI: Image, Keypoints, Mask, testencoding, Label, OneHot, ProjectiveTransforms,
    encodedblock, decodedblock, encode, decode, mockblock
using FilePathsBase
using FastAI.Datasets
using DLPipelines
import DataAugmentation
import DataAugmentation: getbounds
using Flux
using Flux.Optimise: Optimiser, apply!
using StaticArrays
using Test
using TestSetExtensions
using DataFrames
using Tables
using CSV

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
include("testdata.jl")


text:
using ADSeismic
using ADCME
using PyPlot
using PyCall
using DelimitedFiles
using MAT
using Statistics
np = pyimport("numpy")
matplotlib.use("Agg")
use_gpu()

output_dir = "data/acoustic"
if !ispath(output_dir)
  mkpath(output_dir)
end
figure_dir = "figure/acoustic/rupture/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/acoustic/rupture/"
if !ispath(result_dir)
  mkpath(result_dir)
end
# if isfile(joinpath(result_dir, "loss.txt"))
#   rm(joinpath(result_dir, "loss.txt"))
# end

################### Generate synthetic data #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/layer-model-rupture.mat"; IT_DISPLAY=0)
src = load_acoustic_source("models/layer-model-rupture.mat")
rcv_sim = load_acoustic_receiver("models/layer-model-rupture.mat")
Rs_ = compute_forward_GPU(ap_sim, src, rcv_sim)

source_scale = maximum(src[1].srcv)
vmax = maximum(src[1].srcv)
# vmin = minimum(src[1].srcv)
vmin = -vmax

sess = Session(); init(sess)
Rs = run(sess, Rs_)

for i = 1:length(src)
    writedlm(joinpath(output_dir, "rupture-r$i.txt"), Rs[i])
end

u = run(sess, ap_sim(src[1]).u)
visualize_file(u, ap_sim.param, dirs=figure_dir)

################### Inversion using Automatic Differentiation #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/layer-model-rupture.mat"; IT_DISPLAY=0)
rcv_sim = load_acoustic_receiver("models/layer-model-rupture.mat")

ix = matread("models/layer-model-rupture.mat")["rupture"]["ix"][:]
iy = matread("models/layer-model-rupture.mat")["rupture"]["iy"][:]
src_ = [variable_source(ap_sim.param, ix, iy, source_scale=source_scale)]

Rs = Array{Array{Float64,2}}(undef, length(src_))
for i = 1:length(src_)
    Rs[i] = readdlm(joinpath(output_dir, "rupture-r$i.txt"))
end

sess = Session(); init(sess)

vars = [get("src")]
losses, gs = compute_loss_and_grads_GPU(ap_sim, src_, rcv_sim, Rs, vars)
g = Array{PyObject}(undef, length(gs[1]))
for i = 1:length(gs[1])
  g[i] = sum([pp[i] for pp in gs])
end

loss = sum(losses)
sess = Session(); init(sess)

function callback(v, iter, loss)
    if iter%10==0
      V = Array(reshape(v[1],  ap_sim.param.NSTEP+1, length(ix)))*source_scale
      close("all")
      t = collect(0:ap_sim.param.NSTEP)* ap_sim.param.DELTAT
      n = collect(1:length(ix)).-0.5
      n, t = np.meshgrid(n, t)
      pcolormesh(n, t, V, vmax=vmax, vmin=vmin, cmap="seismic")
      xlabel("Location index")
      ylabel("Time (s)")
      savefig(joinpath(figure_dir, "inv_$(lpad(iter,4,"0")).png"))
      writedlm(joinpath(result_dir, "inv_$(lpad(iter,5,"0")).txt"), V)
      open(joinpath(result_dir, "loss.txt"), "a") do io
        writedlm(io, loss)
      end
    else
      return
    end
end

callback(run(sess, vars), 0, run(sess, loss))
LBFGS!(sess, loss, vars=vars, callback=callback)
# LBFGS!(sess, loss, g, vars, callback=callback)