text:
module CrossSectionalCharts

using Colors
using Base64

struct CrossSectionalChartSVG <: Main.SVG
    buf::IOBuffer
end

struct Axis
    index::Int # compornent index
    label::String
    range::AbstractRange
end
Base.first(a::Axis) = first(a.range)
Base.last(a::Axis) = last(a.range)

function crosssection(::Type{C};
                      x::Tuple=(1, "X", 0:1),
                      y::Tuple=(2, "Y", 0:1),
                      z::Tuple=(3, "Z", 0:1)) where C <: Color
    crosssection(C, Axis(x...), Axis(y...), Axis(z...))
end
function crosssection(::Type{C}, x::Axis, y::Axis, z::Axis) where C <: Color
    io = IOBuffer()
    id = String(nameof(C))

    write(io,
        """
        <svg xmlns="http://www.w3.org/2000/svg"
             xmlns:xlink="http://www.w3.org/1999/xlink"
             id="svg_$id" version="1.1"
             viewBox="0 0 40 30" width="27.09mm" height="20.32mm"
             stroke="none" style="display:inline; margin-left:2em; margin-bottom:1em">
        <defs>
        <filter id="filter_g_$id" filterUnits="userSpaceOnUse" x="0" y="-0.5" width="15.5" height="16">
          <feGaussianBlur stdDeviation="0.1" result="blur"/>
          <feMerge>
            <feMergeNode in="SourceGraphic"/>
            <feMergeNode in="blur"/>
          </feMerge>
        </filter>
        <style type="text/css"><![CDATA[
          #svg_$id g, #svg_$id image {
            transition: all 400ms ease 200ms;
          }
          #svg_$id path.b {
            stroke: currentColor;
            stroke-width: 0.5;
            opacity: 0;
            transition: all 200ms ease 0ms;
          }
          #svg_$id image {
            opacity: 0;
            filter:url(#filter_g_$id);
          }
          #svg_$id rect:active ~ image {
            opacity: 1;
            transition: all 200ms ease 0ms;
          }
          #svg_$id rect:hover ~ image {
            opacity: 1;
            transition: all 200ms ease 0ms;
          }
          #svg_$id rect:active ~ path.b {
            opacity: 0.8;
          }
          #svg_$id rect:hover ~ path.b {
            opacity: 0.8;
          }
          #svg_$id text {
            fill: currentColor;
            fill-opacity: 0.8;
            stroke: #aaa;
            stroke-width: 0.2;
            stroke-opacity: 0.4;
            font-size: 3px;
          }
          #svg_$id text.n {
            opacity:0;
          }
          #svg_$id:hover text.n {
            opacity:1;
          }
          #svg_$id:active text.n {
            opacity:1;
          }
        ]]></style>
        </defs>
        """)

    xs = [xv for xv in range(first(x), stop=last(x), length=16)]
    ys = [yv for yv in range(last(y), stop=first(y), length=16)]
    zs = [zv for zv in range(first(z), stop=last(z), length=11)]
    xmidf = (first(x) + last(x)) * 0.5
    ymidf = (first(y) + last(y)) * 0.5
    xmid = isinteger(xmidf) ? Int(xmidf) : xmidf
    ymid = isinteger(ymidf) ? Int(ymidf) : ymidf

    vec = [0.0, 0.0, 0.0]
    function col(xv, yv, zv)
        vec[x.index] = xv
        vec[y.index] = yv
        vec[z.index] = zv
        # TODO: Add more appropriate out-of-gamut color handling
        xyz = convert(XYZ, C(vec...))
        rgb = convert(RGB, XYZ(max(xyz.x,0), max(xyz.y,0), max(xyz.z,0)))
    end

    # add swatches of color bar and planes by layer
    for i = 1:11
        zi = isodd(i) ? 6 - i÷2 : 6 + i÷2 # zigzag order
        plane = [col(xs[xi], ys[yi], zs[zi]) for yi = 1:16, xi = 1:16]
        ccolor = col(xmid, ymid, zs[zi]) # center color (for color bar)
        barh = i == 1 ? 30 : 16.5 - 3*(i÷2)
        op = i == 1 ? "style=\"opacity:1;\"" : ""
        write(io,
            """
            <g>
              <rect fill="#$(hex(ccolor))" width="4" height="$(barh)" x="36" y="$(isodd(i) ? 30 - barh : 0)" />
              <image width="16" height="16" transform="scale(2) translate(-.5,0)" $op xlink:href="data:image/png;base64,""")
        b64enc = Base64EncodePipe(io)
        write_png(b64enc, plane)
        close(b64enc)
        write(io, "\" />\n")
        write(io,
            """
              <path d="m35,$(33-3zi) h 5" class="b"/>
            </g>
            """)
    end
    # add labels
    if first(x.range) * last(x.range) < 0
        write(io,
            """
            <path d="M0,15 h30 M15,0 v30" style="stroke:currentColor;stroke-width:0.125"/>
            <text x="29.5" y="14" style="text-anchor:end;">$(x.label)</text>
            <text x="29.5" y="18" style="text-anchor:end;" class="n">$(last(x))</text>
            <text x="16" y="3" style="text-anchor:start;">$(y.label)</text>
            <text x="14" y="3" style="text-anchor:end;" class="n">$(last(y))</text>
            <text x="14" y="18" style="text-anchor:end;" class="n">0</text>
            """)
    else
        write(io,
            """
            <text x="29.5" y="26" style="text-anchor:end;">$(x.label)</text>
            <text x="29.5" y="29" style="text-anchor:end;" class="n">$(last(x))</text>
            <text x="15" y="29" style="text-anchor:middle;" class="n">$xmid</text>
            <text x="6" y="3" style="text-anchor:start;">$(y.label)</text>
            <text x="0.5" y="3" style="text-anchor:start;" class="n">$(last(y))</text>
            <text x="0.5" y="16" style="text-anchor:start;" class="n">$ymid</text>
            <text x="0.5" y="29" style="text-anchor:start;" class="n">0</text>
            """)
    end
    write(io,
        """
        <text style="text-anchor:middle;" transform="translate(35,15) rotate(-90)">$(z.label)</text>
        <text x="36" y="3" style="text-anchor:end;" class="n">$(last(z))</text>
        <text x="36" y="29" style="text-anchor:end;" class="n">$(first(z))</text>
        <text x="2" y="26" style="fill:#fff;fill-opacity:1;text-anchor:start;">$id</text>
        <path d="m0,0 h40 v30 h-40 z" style="fill:none;stroke:none;" />
        </svg>""")
    CrossSectionalChartSVG(io)
end

function write_png(io::IO, cs::AbstractArray{T}) where T <: Color
    buf = IOBuffer() # to calculate chunk CRCs
    n = 16 # 16 x 16
    u8(x) = write(buf, UInt8(x & 0xFF))
    u16(x) = (u8((x & 0xFFFF)>>8); u8(x))
    u32(x) = (u16((x & 0xFFFFFFFF)>>16); u16(x))
    b(bstr) = write(buf, bstr)
    function palette(c::Color)
        rgb24 = convert(RGB24,c)
        u8(rgb24.color>>16); u8(rgb24.color>>8); u8(rgb24.color)
    end
    crct(x) = (for i = 1:8; x = x & 1==1 ? 0xEDB88320 ⊻ (x>>1) : x>>1 end; x)
    table = UInt32[crct(i) for i = 0x00:0xFF]
    function crc32()
        seekstart(buf)
        crc = 0xFFFFFFFF
        while !eof(buf)
            crc = (crc>>8) ⊻ table[(crc&0xFF) ⊻ read(buf, UInt8) + 1]
        end
        u32(crc ⊻ 0xFFFFFFFF)
    end
    flush() = write(io, take!(seekstart(buf)))

    # The following is a pre-encoded 256-indexed-color PNG with size of 16x16.
    # We only rewrite "pallets".
    b(b"\x89PNG\x0D\x0A\x1A\x0A")
    # Image header
    u32(13); flush(); b(b"IHDR"); u32(n); u32(n); u8(8); u8(3); u8(0); u8(0); u8(0); crc32()
    # Palette
    u32(n * n * 3); flush();
    b(b"PLTE")
    for y = 1:n, x = 1:n
        palette(cs[y,x])
    end
    crc32()
    # Image data
    u32(58); flush(); b(b"IDAT")
    b(b"\x78\xDA\x63\x64\x60\x44\x03\x02\xE8\x02\x0A\xE8\x02\x06\xE8\x02")
    b(b"\x0E\xE8\x02\x01\xE8\x02\x09\xE8\x02\x05\xE8\x02\x0D\xE8\x02\x13")
    b(b"\xD0\x05\x16\xA0\x0B\x6C\x40\x17\x38\x80\x2E\x70\x01\x5D\xE0\x01")
    b(b"\xBA\xC0\x07\x34\x3E\x00\x54\x4D\x08\x81"); crc32()
    # Image trailer
    u32(0); flush(); b(b"IEND"); crc32()
    flush()
end

"""
# Image data
using CodecZlib
raw = IOBuffer()
for y = 0:15
    write(raw, UInt8(1)) # filter: SUB
    write(raw, UInt8(y*16)) # line head
    write(raw, UInt8[1 for i=1:15]) # left + 1
end
flush(raw)
cd = ZlibCompressorStream(raw,level=9)
flush(cd)
seekstart(cd)
@show read(cd) # UInt8[0x78, 0xda, 0x63, 0x64, ...
"""

crosssection(::Type{HSV}) = crosssection(HSV, x=(2, "S", 0:1),
                                              y=(3, "V", 0:1),
                                              z=(1, "H", 0:360))
crosssection(::Type{HSL}) = crosssection(HSL, x=(2, "S", 0:1),
                                              y=(3, "L", 0:1),
                                              z=(1, "H", 0:360))
crosssection(::Type{HSI}) = crosssection(HSI, x=(2, "S", 0:1),
                                              y=(3, "I", 0:1),
                                              z=(1, "H", 0:360))

crosssection(::Type{Lab}) = crosssection(Lab, x=(2, "a*", -100:100),
                                              y=(3, "b*", -100:100),
                                              z=(1, "L*", 0:100))
crosssection(::Type{Luv}) = crosssection(Luv, x=(2, "u*", -100:100),
                                              y=(3, "v*", -100:100),
                                              z=(1, "L*", 0:100))
crosssection(::Type{LCHab}) = crosssection(LCHab, x=(2, "C*", 0:100),
                                                  y=(1, "L*", 0:100),
                                                  z=(3, "H", 0:360))
crosssection(::Type{LCHuv}) = crosssection(LCHuv, x=(2, "C*", 0:100),
                                                  y=(1, "L*", 0:100),
                                                  z=(3, "H", 0:360))

crosssection(::Type{YIQ}) = crosssection(YIQ, x=(2, "I", -1:1),
                                              y=(3, "Q", -1:1),
                                              z=(1, "Y", 0:1))
crosssection(::Type{YCbCr}) = crosssection(YCbCr, x=(2, "Cb", 0:256),
                                                  y=(3, "Cr", 0:256),
                                                  z=(1, "Y", 0:256))
end


text:
function check_folder_integrity(folder::String)
    folder_files = readdir(folder)
    alien_files = filter(!∈(KNOWN_SIMULATION_PATHS), folder_files)
    if isempty(alien_files)
        return true
    end
    if "data_store" ∉ folder_files
        error("The file path doesn't contain any data_store folder")
    end
    return false
end

function _fill_result_value_container(fields)
    return FieldResultsByTime(x => ResultsByTime() for x in fields)
end

struct SimulationResults <: PSIResults
    path::String
    params::SimulationStoreParams
    problem_results::Dict{String, SimulationProblemResults}
end

"""
Construct SimulationResults from a simulation output directory.

# Arguments
- `path::AbstractString`: Simulation output directory
- `execution::AbstractString`: Execution number. Default is the most recent.
- `ignore_status::Bool`: If true, return results even if the simulation failed.
"""
function SimulationResults(path::AbstractString, execution = nothing; ignore_status = false)
    # path will be either the execution_path or the directory containing all executions.
    contents = readdir(path)
    if "data_store" in contents
        execution_path = path
    else
        if execution === nothing
            executions = [parse(Int, f) for f in contents if occursin(r"^\d+$", f)]
            if isempty(executions)
                error("There are no simulation results in the path")
            end
            execution = maximum(executions)
        end
        execution_path = joinpath(path, string(execution))
        if !isdir(execution_path)
            error("Execution $execution not in the simulations results")
        end
    end

    status = deserialize_status(joinpath(execution_path, RESULTS_DIR))
    if status != RunStatus.SUCCESSFUL
        if ignore_status
            @warn "Simulation was not successful: $status. Results may not be valid."
        else
            error(
                "Simulation was not successful: status = $status. Set ignore_status = true to override.",
            )
        end
    end

    if !check_folder_integrity(execution_path)
        @warn "The results folder $(execution_path) is not consistent with the default folder structure. " *
              "This can lead to errors or unwanted results."
    end

    simulation_store_path = joinpath(execution_path, "data_store")
    check_file_integrity(simulation_store_path)

    return h5_store_open(simulation_store_path, "r") do store
        problem_results = Dict{String, SimulationProblemResults}()
        sim_params = get_params(store)
        for (name, problem_params) in sim_params.problems
            name = string(name)
            problem_result = SimulationProblemResults(
                store,
                name,
                problem_params,
                sim_params,
                execution_path,
            )
            problem_results[name] = problem_result
        end

        return SimulationResults(execution_path, sim_params, problem_results)
    end
end

"""
Construct SimulationResults from a simulation.
"""
SimulationResults(sim::Simulation; kwargs...) =
    SimulationResults(get_simulation_dir(sim); kwargs...)

Base.empty!(res::SimulationResults) = foreach(empty!, values(res.problem_results))
Base.isempty(res::SimulationResults) = all(isempty, values(res.problem_results))
Base.length(res::SimulationResults) = mapreduce(length, +, values(res.problem_results))
get_exports_folder(x::SimulationResults) = joinpath(x.path, "exports")

function get_problem_results(results::SimulationResults, problem)
    if !haskey(results.problem_results, problem)
        throw(IS.InvalidValue("$problem is not stored"))
    end

    return results.problem_results[problem]
end

"""
Return the problem names in the simulation.
"""
list_problems(results::SimulationResults) = collect(keys(results.problem_results))

"""
Export results to files in the results directory.

# Arguments
- `results::SimulationResults`: simulation results
- `exports`: SimulationResultsExport or anything that can be passed to its constructor.
  (such as Dict or path to JSON file)

An example JSON file demonstrating possible options is below. Note that `start_time`,
`end_time`, `path`, and `format` are optional.

```
{
  "problems": [
    {
      "name": "ED",
      "variables": [
        "P__ThermalStandard",
        "E__HydroEnergyReservoir"
      ],
      "parameters": [
        "all"
      ]
    },
    {
      "name": "UC",
      "variables": [
        "On__ThermalStandard"
      ],
      "parameters": [
        "all"
      ],
      "duals": [
        "all"
      ]
    }
  ],
  "start_time": "2020-01-01T04:00:00",
  "end_time": null,
  "path": null,
  "format": "csv"
}
```
"""
function export_results(results::SimulationResults, exports)
    simulation_store_path = joinpath(results.path, "data_store")
    h5_store_open(simulation_store_path, "r") do store
        export_results(results, exports, store)
    end
end

function export_results(results::SimulationResults, exports, store::SimulationStore)
    if !(exports isa SimulationResultsExport)
        exports = SimulationResultsExport(exports, results.params)
    end

    file_type = get_export_file_type(exports)

    for problem_results in values(results.problem_results)
        problem_exports = get_problem_exports(exports, problem_results.problem)
        path =
            exports.path === nothing ? problem_results.results_output_folder : exports.path
        for timestamp in get_existing_timestamps(problem_results)
            !should_export(exports, timestamp) && continue

            export_path = mkpath(joinpath(path, problem_results.problem, "variables"))
            for name in get_existing_variables(problem_results)
                if should_export_variable(problem_exports, name)
                    dfs = read_variable(
                        problem_results,
                        name;
                        initial_time = timestamp,
                        count = 1,
                        store = store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end

            export_path = mkpath(joinpath(path, problem_results.problem, "parameters"))
            for name in get_existing_parameters(problem_results)
                if should_export_parameter(problem_exports, name)
                    dfs = read_parameter(
                        problem_results,
                        name;
                        initial_time = timestamp,
                        count = 1,
                        store = store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end

            export_path = mkpath(joinpath(path, problem_results.problem, "duals"))
            for name in get_existing_duals(problem_results)
                if should_export_dual(problem_exports, name)
                    dfs = read_dual(
                        problem_results,
                        name;
                        initial_time = timestamp,
                        count = 1,
                        store = store,
                    )
                    export_result(file_type, export_path, name, timestamp, dfs[timestamp])
                end
            end
        end

        if problem_exports.optimizer_stats
            export_path = joinpath(path, problem_results.problem, "optimizer_stats.csv")
            df = read_optimizer_stats(problem_results, store = store)
            export_result(file_type, export_path, df)
        end
    end
end

function export_result(
    ::Type{CSV.File},
    path,
    name,
    timestamp::Dates.DateTime,
    df::DataFrames.DataFrame,
)
    filename = joinpath(path, string(name) * "_" * convert_for_path(timestamp) * ".csv")
    export_result(CSV.File, filename, df)
end

function export_result(::Type{CSV.File}, path, name, df::DataFrames.DataFrame)
    filename = joinpath(path, string(name) * ".csv")
    export_result(CSV.File, filename, df)
end

function export_result(
    ::Type{CSV.File},
    path,
    timestamp::Dates.DateTime,
    df::DataFrames.DataFrame,
)
    filename = joinpath(path, convert_for_path(timestamp) * ".csv")
    export_result(CSV.File, filename, df)
end

function export_result(::Type{CSV.File}, filename, df::DataFrames.DataFrame)
    open(filename, "w") do io
        CSV.write(io, df)
    end

    @debug "Exported $filename"
end

ext: jl


text:
# install PkgFix if not present
if !isdir(joinpath(Pkg.dir(), "PkgFix"))
  Pkg.clone("https://github.com/OptimalDesignLab/PkgFix.jl.git")
end
start_dir=pwd()
cd(Pkg.dir("PkgFix"))
run(`git checkout upgrade_0.6`)
cd(start_dir)

using PkgFix  # from now on, use PkgFix instead of Pkg for everything

pkgs = ["https://github.com/OptimalDesignLab/PumiInterface.jl.git" "v0.10";
#        "https://github.com/JuliaParallel/MPI.jl.git"  "v0.5.0";
        "https://github.com/OptimalDesignLab/SummationByParts.jl.git"  "jcwork";
#        "https://github.com/JuliaLang/ArrayViews.jl.git" "93e80390aeedb1dbcd90281b6dff7f760f430bc8";
#        "https://github.com/jipolanco/WriteVTK.jl.git"  "v0.6.1";
        "https://github.com/OptimalDesignLab/ODLCommonTools.jl.git" "master"]

pkg_dict = PkgFix.installed()

for i=1:size(pkgs, 1)
  pkg_name = PkgFix.getRepoName(pkgs[i, 1])
  if !haskey(pkg_dict, pkg_name)
    PkgFix.add(pkgs[i, 1], branch_ish=pkgs[i, 2])
  end
end

global const ARRAYVIEWS_URL = "https://github.com/JaredCrean2/ArrayViews.jl.git"
global const ARRAYVIEWS_VER = "work"

if !haskey(pkg_dict, "ArrayViews")
  PkgFix.add(ARRAYVIEWS_URL, branch_ish=ARRAYVIEWS_VER)
else
  PkgFix.checkout("ArrayViews", ARRAYVIEWS_VER)
end


text:
mutable struct Sphere{T} <: Intersectable{T}
    pos::SVector{3,T}
    radius::T
end

#=
\[d=-[\mathbf{l}\cdot(\mathbf{o}-\mathbf{c})]\pm\sqrt{[\mathbf{l}\cdot(\mathbf{o}-\mathbf{c})]^2-||\mathbf{o}-\mathbf{c}||^2+r^2}\]
=#

function extents(sphere::Sphere{T}) where T
    e = √(3one(T))
    [(pos[i]-e,pos[i]+e) for i=1:3]
end

function Base.intersect(ray::Ray{T}, sphere::Sphere{T}) where T
    oc = ray.pos - sphere.pos
    loc = dot(ray.dir, oc)
    loc^2 - dot(oc, oc) + sphere.radius^2 >= 0
end

function Base.intersect!(i::Intersection{T}, sphere::Sphere{T}) where T
    oc = i.ray.pos - sphere.pos
    loc = dot(i.ray.dir, oc)
    D = loc^2 - dot(oc, oc) + sphere.radius^2
    (D < 0) && return
    t = D != 0 ? -loc + sign(loc)*sqrt(D) : -loc
    (t > i.t) && return
    i.o = sphere
    i.t = t
end

normal(sphere::Sphere{T}, p::SVector{3,T}, ::Intersection{T}) where T =
    normalize(p-sphere.pos)

aabb(s::Sphere{T}) where T = AABB{T}(s.pos-s.radius,s.pos+s.radius)

export Sphere, extents, normal, aabb


text:
@testset "Resampling" begin
  c=Cube()
  d=subsetcube(c,variable="gross_primary_productivity",region="Germany",time=(Date(2001,1,1),Date(2001,1,15)))
  newlons=5.875:0.05:14.875
  newlats=54.875:-0.05:47.375
  shigh = spatialinterp(d,newlons,newlats)
  @test d[37,31,1]==shigh[181,151,1]
  @test d[37,30,1]==shigh[181,146,1]
end


text:
const LR = 0.0001
const dtau = 34 # sec
const ptau = 17 # sec

# is positive when Δt > 0
# note, this is defined on a ms scale
function stdp(Δt)
	# stdp function
    if Δt > 0
        return LR*exp(-Δt/ptau)
    elseif Δt < 0
        return -LR*exp(Δt/dtau)
    else # Δt=0, pre and post are simultaneous - no modification
        return 0
    end
end

# in the paper, t_pre = time when presynaptic cell activated
# t_post is when postsynaptic cell is activated. Therefore, Δt computed from
# spike times, not last spike arrival time of pre
function updateSTDP!(syn::SpikingSynapse, source::AbstractSpikingNeuron, target::AbstractSpikingNeuron)
    if source.didSpike || target.didSpike
        syn.w += stdp(target.lastSpike - source.lastSpike)
        syn.w = max(0,syn.w) # make sure weights dont go negative
    end
end


text:
    1+1

## --- Load external packages

    using StatGeochem
    using Plots; gr();
    using ProgressMeter: @showprogress

    if VERSION>=v"0.7"
        using Statistics, DelimitedFiles, SpecialFunctions
    end

## --- Load the OCTOPUS kml file
    # Unzip the kml file

    function load_octopus_kml(filename)

        # Check number and ordering of coordinate lists
        fid = open("data/octopus/crn_basins_global.kml")

        # Count the number of lines in the file
        i = 0
        for line in eachline(fid)
            i += 1
        end
        isfirstcoord = Array{Bool}(undef, i)

        # Parse the coords (basin polygon outlines) in each placemark
        i = 0
        ncoords = 0
        lastcoord = 0
        lastplacemark = 0
        seekstart(fid)
        for line in eachline(fid)
            i += 1
            if occursin("</ExtendedData>", line)
                lastplacemark = i
            elseif occursin("<coordinates>", line)
                ncoords += 1
                # If we've had a new placemark more recently than a cordinate, then
                # this is the first coord of the placemark
                if lastplacemark > lastcoord
                    isfirstcoord[ncoords] = true
                else
                    isfirstcoord[ncoords] = false
                end
                lastcoord = i
            end
        end
        isfirstcoord = isfirstcoord[1:ncoords]

        # Read entire string
        seekstart(fid)
        str = read(fid,String)

        # Close the file
        close(fid)

        nbasins = 0
        subbasins = Array{Any}(undef, sum(isfirstcoord))
        for i = 1:length(isfirstcoord)
            if isfirstcoord[i]
                nbasins += 1
                subbasins[nbasins] = i
            else
                subbasins[nbasins] = vcat(subbasins[nbasins],i)
            end
        end

        return (str, isfirstcoord, nbasins, subbasins)
    end

    run(`gunzip data/octopus/crn_basins_global.kml.gz`); # Decompress

    (str, isfirstcoord, nbasins, subbasins) = load_octopus_kml("data/octopus/crn_basins_global.kml")

    run(`gzip data/octopus/crn_basins_global.kml`); # Recompress

## --- Parse the file

    function parse_octopus_kml_variables(str)
        # Find a list of all the string variables
        lm = eachmatch(r"<SimpleField type=\"string\" name=\"(.*?)\"", str);
        stringvars = unique(map(x -> x[1], lm))

        # Find a list of all the numeric variables
        lm = eachmatch(r"<SimpleField type=\"double\" name=\"(.*?)\"", str);
        numvars = unique(map(x -> x[1], lm))

        # Make a dict to store all our parsed data
        data = Dict();

        # Parse all the string variables
        for i = 1:length(stringvars)
            lm = eachmatch(Regex("<SimpleData name=\"$(stringvars[i])\">(.*?)</SimpleData>"),str);
            data[stringvars[i]] = map(x -> x[1], lm);
        end

        # Parse all the numeric variables
        for i = 1:length(numvars)
            lm = eachmatch(Regex("<SimpleData name=\"$(numvars[i])\">(.*?)</SimpleData>"),str);
            data[numvars[i]] = map(x -> parse(Float64, x[1]), lm)
        end

        needsNaNs = ["alcorr", "be10ep", "beprod", "eal_err", "elev_std", "errbe_prod", "al26ep", "alprod", "be10ep_err", "beself", "eal_gcmyr", "erral_ams", "errbe_tot", "al26ep_err", "alself", "be10nc", "besnow", "eal_mmkyr", "erral_muon", "sizemax", "al26nc", "alsnow", "be10nc_err", "betopo", "ebe_err", "erral_prod", "sizemin", "al26nc_err", "altopo", "be10np", "betots", "ebe_gcmyr", "erral_tot", "slp_ave", "al26np", "altots", "be10np_err", "ebe_mmkyr", "errbe_ams", "slp_std", "al26np_err", "area", "becorr", "dbver", "elev_ave", "errbe_muon", ]

        for i = 1:length(needsNaNs)
            data[needsNaNs[i]][data[needsNaNs[i]] .< 0] .= NaN
        end

        return data
    end

    data = parse_octopus_kml_variables(str)

    # Parse the basin polygon outlines
    function parse_octopus_polygon_outlines(str,isfirstcoord)
        i = 0
        n = 0
        basin_polygon_lat = Array{Array}(undef, length(isfirstcoord))
        basin_polygon_lon = Array{Array}(undef, length(isfirstcoord))
        basin_polygon_n = Array{Int}(undef, length(isfirstcoord))
        for m in eachmatch(r"<coordinates>\n\t*(.*?)\n",str)
            i += 1 # Number of parsed values
            parsed = delim_string_function(x -> delim_string_parse(x, ',', Float64), m[1], ' ', Array{Float64,1})
            nparsed = length(parsed)

            basin_polygon_n[i] = nparsed;
            basin_polygon_lon[i] = parsed .|> x -> x[1]
            basin_polygon_lat[i] = parsed .|> x -> x[2]
        end
        return (basin_polygon_n, basin_polygon_lat, basin_polygon_lon)
    end

    (basin_polygon_n, basin_polygon_lat, basin_polygon_lon) = parse_octopus_polygon_outlines(str,isfirstcoord)

## --- Calculate precipitation for each basin

    using NetCDF
    prate = ncread("data/prate.sfc.mon.ltm.nc","prate")
    lat = ncread("data/prate.sfc.mon.ltm.nc","lat")
    lon = ncread("data/prate.sfc.mon.ltm.nc","lon")

    precip = mean(prate, dims=3)[:,:,1]
    imsc(collect(precip'), viridis)

    function find_precip(lat,lon)
        prate = ncread("prate.sfc.mon.ltm.nc","prate")
        precip = mean(prate, dims=3)[:,:,1]

        out = Array{Float64}(undef,size(lat))
        for i=1:length(lat)
            if isnan(lat[i]) || isnan(lon[i]) || lat[i] > 88.542 || lat[i] < -88.542
                out[i] = NaN
            else
                row = trunc(Int, mod(lon[i]*192/360 + 0.5, 192.0)) + 1
                col = trunc(Int, 46.99999 - lat[i]*94/177.084) + 1
                out[i] = precip[row, col]
            end
        end

        return out
    end

    data["precip"] = find_precip(data["y_wgs84"],data["x_wgs84"])

    # # Test find_precip to make sure it's not flipped
    # lat = repeat(-90:90,1,361)
    # lon = repeat((-180:180)',181,1)
    # imsc(find_precip(lat,lon),viridis)

## --- Recalculate slopes for each basin

    function get_basin_srtm15plus_aveslope(srtm::Dict,nbasins,subbasins,basin_polygon_lat,basin_polygon_lon)
        slope = srtm["slope"]
        x_lon_cntr = srtm["x_lon_cntr"]
        y_lat_cntr = srtm["y_lat_cntr"]

        basin_srtm15plus_aveslope = Array{Float64}(undef, nbasins)
        basin_N = Array{Int64}(undef, nbasins)
        for i = 1:nbasins
            rowsinbasin = Array{Int}(undef, 0)
            columnsinbasin = Array{Int}(undef, 0)
            for j = 1:length(subbasins[i])
                k = subbasins[i][j]
                subbasin_lon = basin_polygon_lon[k]
                subbasin_lat = basin_polygon_lat[k]
                (row, column) = find_grid_inpolygon(x_lon_cntr, y_lat_cntr, subbasin_lon, subbasin_lat)
                rowsinbasin = vcat(rowsinbasin,row)
                columnsinbasin = vcat(columnsinbasin,column)
            end

            pointsinbasin = unique(hcat(rowsinbasin,columnsinbasin),dims=1);
            basin_slopes = Array{UInt16}(undef, length(rowsinbasin))
            for j=1:size(pointsinbasin,1)
                basin_slopes[j] = slope[pointsinbasin[j,1],pointsinbasin[j,2]]
            end

            # Average all the slopes
            basin_srtm15plus_aveslope[i] = mean(basin_slopes)
            basin_N[i] = length(basin_slopes)

            print("Basin: $i, slope: $(round(basin_srtm15plus_aveslope[i],digits=3)), N: $(basin_N[i]) \n")
        end

        return basin_srtm15plus_aveslope
    end

    # using HDF5
    # srtm = h5read("data/srtm15plus_aveslope.h5","vars/")
    # basin_srtm15plus_aveslope = get_basin_srtm15plus_aveslope(srtm,nbasins,subbasins,basin_polygon_lat,basin_polygon_lon)

    using JLD
    # save("basin_srtm15plus_aveslope.jld","basin_srtm15plus_aveslope",basin_srtm15plus_aveslope)
    basin_srtm15plus_aveslope =  load("basin_srtm15plus_aveslope.jld")["basin_srtm15plus_aveslope"]

    # save("OctopusSlopeRecalc.jld","basin_srtm15plus_aveslope",basin_srtm15plus_aveslope,"basin_polygon_n",basin_polygon_n,"basin_polygon_lat",basin_polygon_lat,"basin_polygon_lon",basin_polygon_lon,"subbasins",subbasins)

## --- Plot raw Lat and Lon
    h = plot(data["x_wgs84"],data["y_wgs84"],seriestype=:scatter,label="basin locations")
    plot!(h, xlabel="Longitude", ylabel="Latitude")
    savefig(h, "Basin_Locations.pdf")
    display(h)

## --- Plot raw erosion rate vs basin area
    h = plot(data["area"],data["ebe_mmkyr"],seriestype=:scatter,label="OCTOPUS Be-10 data")
    plot!(h, data["area"],data["eal_mmkyr"],seriestype=:scatter,label="OCTOPUS Al-26 data")
    plot!(h, xlabel="Basin area (km^2)", ylabel="Erosion rate (mm/kyr)",yscale=:log10, xscale=:log10, legend=:topleft, fg_color_legend=:white)
    savefig(h, "Area_vs_erosion.pdf")
    display(h)

## --- Plot raw erosion rate vs latitude
    h = plot(data["y_wgs84"],data["ebe_mmkyr"],seriestype=:scatter,label="OCTOPUS Be-10 data")
    plot!(h,data["y_wgs84"],data["eal_mmkyr"],seriestype=:scatter,label="OCTOPUS Al-26 data")
    plot!(h, xlabel="Latitude", ylabel="Erosion rate (mm/kyr)", yscale=:log10,legend=:topleft, fg_color_legend=:white)
    savefig(h, "Latitude_vs_erosion.pdf")
    display(h)

## --- Plot raw erosion rate vs precipitation
    h = plot(data["precip"]*31557600,data["ebe_mmkyr"],seriestype=:scatter,label="OCTOPUS Be-10 data")
    plot!(h,data["precip"]*31557600,data["eal_mmkyr"],seriestype=:scatter,label="OCTOPUS Al-26 data")
    plot!(h, xlabel="Precipitation (kg/m^2/yr)", ylabel="Erosion rate (mm/kyr)", yscale=:log10, fg_color_legend=:white)
    savefig(h, "Precipitation_vs_erosion.pdf")
    display(h)

## --- Plot raw erosion rate vs elevation
    h = plot(data["elev_ave"],data["ebe_mmkyr"],seriestype=:scatter,label="OCTOPUS Be-10 data")
    plot!(h,data["elev_ave"],data["eal_mmkyr"],seriestype=:scatter,label="OCTOPUS Al-26 data")
    plot!(h, xlabel="Elevation (m)", ylabel="Erosion rate (mm/kyr)", yscale=:log10,legend=:topleft, fg_color_legend=:white)
    savefig(h, "Elevation_vs_erosion.pdf")
    display(h)

## --- Plot new slope vs erosion rate
    h = plot(basin_srtm15plus_aveslope,data["ebe_mmkyr"],seriestype=:scatter,label="OCTOPUS Be-10 data")
    plot!(h,basin_srtm15plus_aveslope,data["eal_mmkyr"],seriestype=:scatter,label="OCTOPUS Al-26 data")
    plot!(h, xlabel="SRTM15 Slope (m/km)", ylabel="Erosion rate (mm/kyr)", yscale=:log10,legend=:topleft, fg_color_legend=:white)
    savefig(h, "Slope_vs_Erosion.pdf")
    display(h)

    # log10(E) = slope/140 + 0.7
    # E = 10^(slope/140 + 0.7)

## --- Plot new slope vs precipitation
    h = plot(basin_srtm15plus_aveslope,data["precip"]*31557600,seriestype=:scatter,label="")
    plot!(h, xlabel="SRTM15 Slope (m/km)", ylabel="Precipitation (kg/m^2/yr)", yscale=:log10,legend=:topleft, fg_color_legend=:white)
    savefig(h, "Slope_vs_Precipitation.pdf")
    display(h)

## --- Fit raw erosion rate as a function of slope (m/km)

    using LsqFit: curve_fit

    function linear(x,p)
        y = p[1] .+ x * p[2]
    end
    p = [0.5, 1/100]
    t = .~isnan.(basin_srtm15plus_aveslope) .& .~isnan.(data["ebe_mmkyr"]) .& (basin_srtm15plus_aveslope .< 300)
    x = basin_srtm15plus_aveslope[t]
    y = log10.(data["ebe_mmkyr"][t])
    t = .~isnan.(basin_srtm15plus_aveslope) .& .~isnan.(data["eal_mmkyr"]) .& (basin_srtm15plus_aveslope .< 300)
    x = append!(x, basin_srtm15plus_aveslope[t])
    y = append!(y, log10.(data["eal_mmkyr"][t]))
    fobj = curve_fit(linear, x, y, p);

    # fobj.param[1]: 0.987237
    # fobj.param[2]: 0.00555159
    function Emmkyr(slp)
        return 10^(slp*0.00567517 + 0.971075)
    end

    h = plot(basin_srtm15plus_aveslope,data["ebe_mmkyr"], seriestype=:scatter, label="OCTOPUS Be-10 data")
    plot!(h, basin_srtm15plus_aveslope,data["eal_mmkyr"], seriestype=:scatter, label="OCTOPUS Al-26 data")
    plot!(h, xlabel="SRTM15 Slope (m/km)", ylabel="Erosion rate (mm/kyr)", yscale=:log10, legend=:topleft)
    plot!(h, 1:500, Emmkyr.(1:500), label = "E = 10^(slp/176.2 + 0.971)", fg_color_legend=:white)
    savefig(h, "Slope_vs_Erosion_Fitted.pdf")
    display(h)

## --- Plot Be-10 vs Al-26 erosion rate
    h = plot(data["ebe_mmkyr"],data["eal_mmkyr"],seriestype=:scatter,label="")
    plot!(h, xlabel="Be-10 erosion rate (t1/2=1.387E6) (mm/kyr) ", ylabel="Al-26 erosion rate (t1/2=7.17E5) (mm/kyr)",xscale=:log,yscale=:log)
    plot!(h,10.0.^(-1:4),10.0.^(-1:4),label="1:1",legend=:topleft,fg_color_legend=:white)

    function linslp1(x,p)
        y = p[1] .+ x
    end
    p = [0.0]
    t = .~isnan.(data["ebe_mmkyr"]) .& .~isnan.(data["eal_mmkyr"])
    x = log10.(data["ebe_mmkyr"][t])
    y = log10.(data["eal_mmkyr"][t])
    fobj = curve_fit(linslp1, x, y, p);

    plot!(h,10.0.^(-1:4), 10.0.^((-1:4) .+ fobj.param[1]), label="1:$(10.0^fobj.param[1])")
    savefig(h, "Be_vs_Al_rate.pdf")
    display(h)

## --- Plot residual erosion rate as a function of latitude
    resid = data["ebe_mmkyr"] - Emmkyr.(basin_srtm15plus_aveslope)
    h = plot(abs.(data["y_wgs84"][resid.<0]), abs.(resid[resid.<0]), seriestype=:scatter, label="negative residuals", color=:red)
    plot!(h, abs.(data["y_wgs84"][resid.>0]), resid[resid.>0], seriestype=:scatter, label="positive residuals", color=:blue)
    plot!(h, xlabel="Latitude", ylabel="Residual erosion rate (mm/kyr)", yscale=:log10, legend=:topleft, fg_color_legend=:white)
    savefig(h, "Latitude_vs_erosion_residual.pdf")
    display(h)

## --- Plot residual erosion rate as a function of precipitation
    resid = data["ebe_mmkyr"] - Emmkyr.(basin_srtm15plus_aveslope)
    h = plot(data["precip"][resid.<0]*31557600, abs.(resid[resid.<0]), seriestype=:scatter, label="negative residuals", marker=(0.3,:red))
    plot!(h, data["precip"][resid.>0]*31557600, resid[resid.>0], seriestype=:scatter, label="positive residuals", marker=(0.3,:blue))
    plot!(h, xlabel="Precipitation (kg/m^2/yr)", ylabel="Residual erosion rate (mm/kyr)", yscale=:log10, fg_color_legend=:white)
    savefig(h, "Precipitation_vs_erosion_residual.pdf")
    display(h)

## --- Plot residual erosion rate as a function of elevation
    resid = data["ebe_mmkyr"] - Emmkyr.(basin_srtm15plus_aveslope)
    h = plot(abs.(data["elev_ave"][resid.<0]), abs.(resid[resid.<0]), seriestype=:scatter, label="negative residuals",  marker=(0.3,:red))
    plot!(h, abs.(data["elev_ave"][resid.>0]), resid[resid.>0], seriestype=:scatter, label="positive residuals", marker=(0.3,:blue))
    # h = plot(data["elev_ave"],abs.(resid),seriestype=:scatter,label="OCTOPUS 10Be data",legend=:topleft)
    plot!(h, xlabel="Elevation", ylabel="Residual erosion rate (mm/kyr)",yscale=:log10,fg_color_legend=:white)
    savefig(h, "Elevation_vs_erosion_residual.pdf")
    display(h)

## --- Plot residual erosion ratio as a function of elevation
    resid = data["ebe_mmkyr"] ./ Emmkyr.(basin_srtm15plus_aveslope)
    h = plot(data["elev_ave"],abs.(resid),seriestype=:scatter,label="OCTOPUS 10Be data",legend=:topleft)
    plot!(h, xlabel="Elevation", ylabel="Residual erosion ratio",yscale=:log10)
    # savefig(h, "Latitude_vs_erosion_residual.pdf")
    display(h)

## ---
#
#     # nbins = 20;
#     # nsims = 10000;
#     # (c, m, el, eu) = bin_bsr_means(data["slp_ave"], log10.(data["ebe_mmkyr"]), 0, 800, nbins, data["slp_std"]/10, nsims)
#     # plot(c, m, yerror=(el,eu), seriestype=:scatter, label="",xlabel="Slope", ylabel="Log10 Erosion rate (mm/kyr)")
#
#     (c,m,sigma) = binmeans(data["slp_ave"], log10.(data["ebe_mmkyr"]), 0, 800, 20)
#     plot(c,m,yerror=2*sigma,seriestype="scatter")
#     plot!(xlabel="Slope", ylabel="Log10 Erosion rate (mm/kyr)")
#
# ## ---
#     (c,m,sigma) = binmeans(data["elev_ave"], data["ebe_mmkyr"], 0, 6000, 12)
#     plot(c,m,yerror=2*sigma,seriestype="scatter")
#     plot!(xlabel="Elevation (m)", ylabel="Erosion rate (mm/kyr)")
#
# ## --- Plot raw erosion rate vs elevation scaled to snowline
#     plot(data["elev_ave"]./data["Snowline"],data["ebe_mmkyr"],seriestype=:scatter)
#     plot!(xlabel="Elevation/ELA", ylabel="Erosion rate (mm/kyr)",yscale=:log10)
#
# ## ---
#     (c,m,sigma) = binmeans(data["elev_ave"]./data["Snowline"], data["ebe_mmkyr"], 0, 1.2, 12)
#     plot(c,m,yerror=2*sigma,seriestype="scatter")
#     plot!(xlabel="Elevation/ELA", ylabel="Erosion rate (mm/kyr)")
#
# ## --- Plot resampled elevation vs erosion rate
#     nbins = 20;
#     nsims = 10000;
#     (c, m, el, eu) = bin_bsr_means(data["elev_ave"], data["ebe_mmkyr"], 0, 6000, nbins, data["elev_std"], nsims)
#     plot(c, m, yerror=(el,eu), seriestype=:scatter, label="",xlabel="Elevation (m)", ylabel="Erosion rate (mm/kyr)")
#
# ## --- Plot resampled elevation / snowline vs erosion rate
#     nbins = 20;
#     nsims = 10000;
#     (c, m, el, eu) = bin_bsr_means(data["elev_ave"]./data["Snowline"], data["ebe_mmkyr"], 0, 1.2, nbins, data["elev_std"]./data["Snowline"], nsims)
#     plot(c, m, yerror=(el,eu), seriestype=:scatter, label="")
#     plot!(xlabel="Elevation/ELA", ylabel="Erosion rate (mm/kyr)")

## --- End of File


text:
# SPDX-License-Identifier: X11
# 2020-11-15
# Collinearity (300pt)

function main()
  N = parse(Int, readline())
  A = Vector{Tuple{Int, Int}}()
  for i ∈ 1:N
    x,y = parse.(Int, (split ∘ readline)())
    push!(A, (x,y))
  end

  for i ∈ 1:N - 2, j ∈ i + 1:N - 1, k ∈ j + 1:N
    if trianglearea(A[i], A[j], A[k]) == 0
      println("Yes")
      return
    end
  end

  println("No")
end

trianglearea((xa,ya)::Tuple{<:Real, <:Real},
             (xb,yb)::Tuple{<:Real, <:Real},
             (xc,yc)::Tuple{<:Real, <:Real}) =
  abs(xa * yb - xa * yc + xb * yc - xb * ya + xc * ya - xc * yb) / 2

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Documenter, Primes

makedocs(
    modules = [Primes],
    clean = false,
    format = :html,
    sitename = "Primes.jl",
    pages = Any[
        "Home" => "index.md",
        "Functions" => "api.md"
    ],
)

deploydocs(
    julia = "nightly",
    repo = "github.com/JuliaMath/Primes.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
"""
    replace_missing(a::AbstractRaster, newmissingval)
    replace_missing(a::AbstractRasterStack, newmissingval)

Replace missing values in the array or stack with a new missing value,
also updating the `missingval` field/s.

# Keywords

- `filename`: a filename to write to directly, useful for large files.
- `suffix`: a string or value to append to the filename.
    A tuple of `suffix` will be applied to stack layers. `keys(st)` are the default.

# Example

```jldoctest
using Rasters
A = Raster(WorldClim{Climate}, :prec; month=1) |> replace_missing
missingval(A)
# output
missing
```

"""
replace_missing(x; missingval=missing, kw...) = replace_missing(x, missingval; kw...)
function replace_missing(A::AbstractRaster{T}, missingval::MV;
    filename=nothing, suffix=nothing
) where {T,MV}
    MT = if ismissing(missingval)
        promote_type(T, Missing)
    else
        nonmissingtype(T)
    end
    old_missingval = Rasters.missingval(A)
    missingval = convert(MT, missingval)
    repmissing(x) = isequal(x, old_missingval) || ismissing(x) ? missingval : x
    # Disk-backed arrays need to be lazy, memory-backed don't.
    # But in both cases we make sure we return an array with the missingval
    # in the eltype, even if there are no missing values in the array.
    if !isnothing(filename)
        A1 = create(filename, MT, dims(A);
            parent=parent(A), suffix, missingval, name=name(A), metadata=metadata(A)
        )
        open(A1; write=true) do O
            O .= repmissing.(A)
        end
        return A1
    else
        # We need to force T of Union{T,Missing} for DiskArrays broadcasts
        if isdisk(A)
            data = repmissing.(parent(A))
            if missingval isa Missing
                data = MissingDiskArray(MT, data)
            end
        else
            data = similar(parent(A), MT)
            data .= repmissing.(parent(A))
        end
        return rebuild(A; data, missingval)
    end
end
function replace_missing(st::AbstractRasterStack, args...; suffix=keys(st), kw...)
    mapargs(st, suffix) do A, s
        replace_missing(A, args...; suffix=s, kw...)
    end
end
function replace_missing(s::AbstractRasterSeries, args...; kw...)
    map(x -> replace_missing(x, args...; kw...), s)
end

function replace_missing!(s::RasterSeriesOrStack, args...; kw...)
    map(x -> replace_missing!(x, args...; kw...), s)
end
function replace_missing!(A::AbstractRaster, missingval=missing)
    repmissing(x) = isequal(x, Rasters.missingval(A)) ? missingval : x
    A .= repmissing.(A)
    return rebuild(A; missingval)
end


text:
#
# Computes the moving average of the data
#

struct MovingAverage
  n :: Int64
  x :: Vector{Float64}
  R :: Float64
  residues :: Vector{Float64}
end

"""

`movingaverage(x,n)` or `movavg(x,n)`

Computes a moving average of `x[i]` in the range `i±(n-1)/2`. If `n` is even, we do `n ←  n + 1`.

# Example

```jldoctest
julia> x = rand(10);

julia> movingaverage(x,3)

 ------------------- Moving Average ----------

 Number of points averaged: 3 (± 1 points)

 Pearson correlation coefficient, R = 0.3532754137104625

 Averaged X: x = [0.5807828672543551, 0.40496733381946143...
 residues = [-0.22791917753944557, 0.4037347109743393...

 --------------------------------------------
```

"""
function movingaverage( X :: AbstractArray{<:Real}, n :: Int )
  if n == 0
    error(" Please set n with movavg(x,n=10), for example. ")
  end
  if ! isodd(n)
    n = n + 1
  end
  delta = round(Int64,(n-1)/2)
  nX = length(X)
  x = similar(X)
  for i in 1:nX
    jmin = max(i-delta,1)
    jmax = min(i+delta,nX)
    x[i] = 0.
    for j in jmin:jmax
      x[i] = x[i] + X[j]
    end
    x[i] = x[i] / (jmax-jmin+1)
  end
  residues = similar(X)
  for i in 1:nX
    residues[i] = X[i] - x[i]
  end
  R = Statistics.cor(x,X)
  return MovingAverage(n,x,R,residues)
end
movingaverage(X::AbstractArray{<:Real};n::Int=0) = movingaverage(X,n)
movavg = movingaverage

function Base.show( io :: IO, fit :: MovingAverage )
  println("")
  println(" ------------------- Moving Average ---------- ")
  println("")
  println(" Number of points averaged: $(fit.n) (± $(round(Int64,(fit.n-1)/2)) points)")
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println("")
  println(" Averaged X: x = [",fit.x[1],", ",fit.x[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" -------------------------------------------- ")
end

export movingaverage, movavg


text:
#= ==========================================================================
This solution is the translation from Python to Julia
of the wondeful solution by Søren Fuglede Jørgensen:

https://github.com/fuglede/adventofcode/blob/master/2020/day20/solutions.py

Thx Søren!
========================================================================== =#

using LinearAlgebra
using LightGraphs

Tile = BitArray{2}
Tiles = Dict{Int,Tile}

function readinput(str)
    tiles = Tiles()

    blocks = split(str, "\n\n")
    for b ∈ blocks
        bls = split(b, "\n")
        m = match(r"Tile (\d+):", bls[1])
        n = parse(Int, m[1])
        tiles[n] = Tile(hcat([[c == '#' for c ∈ l] for l ∈ bls[2:end]]...)')
    end
    tiles
end

struct D8
    tile::Tile
end

Base.iterate(S::D8) = (S.tile, (1, S.tile))

function Base.iterate(S::D8, state)

    if state[1] == 8
        nothing
    elseif isodd(state[1])
        (reverse(state[2], dims=2), (state[1] + 1, rotr90(state[2])))
    else
        (state[2], (state[1] + 1, state[2]))
    end
end

Base.length(S::D8) = 8

const N, E, S, W = 1, 2, 3, 4
const opposite = [S,W,N,E]

const edge = [[1, :], [:, 10], [10, :], [:, 1]]
# N, E, S, W edges of a tile

const offset = [(-1, 0), (0, 1), (1, 0), (0, -1)]
# N, E, S ,W offset of facing tiles

function Base.match(tile₁, tile₂)
    for tile ∈ D8(tile₂)
        for dir ∈ [N,E,S,W]
            if tile₁[edge[dir]...] == tile[edge[opposite[dir]]...]
                return offset[dir], tile
            end
        end
    end
    nothing
end

# First Half

puzzleinput = open("inputs/day20.txt") do file
    strip(read(file, String))
end

images = readinput(puzzleinput)

ids = collect(keys(images))
tiles = collect(values(images))

G = SimpleGraph(length(tiles))
combs = collect(Iterators.product(tiles, tiles))
for i in CartesianIndices(combs)
    if i[1] > i[2] && !isnothing(match(combs[i][1], combs[i][2]))
        add_edge!(G, i[1], i[2])
    end
end

prod(ids[degree(G) .== 2])
# Tiles at the edge of the image also have this border,
# but the outermost edges won't line up with any other tiles.
# It follows that each corner tile had only 2 tiles that match its edges.

# Second Half

puzzleinput = open("inputs/day20.txt") do file
    strip(read(file, String))
end

images = readinput(puzzleinput)

ids = collect(keys(images))
tiles = collect(values(images))
vrtcs = Dict([id => v for (v, id) ∈ enumerate(ids)])

G = SimpleGraph(length(tiles))
combs = collect(Iterators.product(tiles, tiles))
for i in CartesianIndices(combs)
    if i[1] > i[2] && !isnothing(match(combs[i][1], combs[i][2]))
        add_edge!(G, i[1], i[2])
    end
end

root, _ = first(images) # Anyone will be to start...
stack = [(root, (0, 0))]
locs = Dict((0, 0) => root)

while !isempty(stack)
    # Place the tiles on a grid, starting from root

    id₁, loc = pop!(stack)
    v₁ = vrtcs[id₁]
    for v₂ ∈ neighbors(G, v₁)
        id₂ = ids[v₂]
        if id₂ ∈ values(locs)
            continue
        end
        off, tile = match(images[id₁], images[id₂])
        images[id₂] = tile # Make sure to update the values of images as we apply a D8 transformation
        newloc = loc .+ off
        push!(stack, (id₂, newloc))
        locs[newloc] = id₂
    end
end

bigsz = Int(sqrt(length(images)))
bigimg = BitArray(undef, bigsz * 8, bigsz * 8)
# bigimag will contain the final image; its size is bigsz * 8 and not bigsz * 10,
# because the borders of each tile are not part of the actual image

offsetrc = (minimum([r for (r, _) ∈ keys(locs)]), minimum([c for (_, c) ∈ keys(locs)]))
# The root tile not necessarily is at top,left corner of the big image matrix

for i ∈ 0:bigsz - 1
    for j ∈ 0:bigsz - 1
        id = locs[(i, j) .+ offsetrc]

        tile = images[id][2:end - 1,2:end - 1]
        bigimg[1 + i * 8:(i + 1) * 8, 1 + j * 8:(j + 1) * 8] = tile
        # The borders of each tile are not part of the actual image

    end
end

# using Images
# using ImageView

# imshow(.!(bigimg))

monster = """
                  #
#    ##    ##    ###
 #  #  #  #  #  #   """
monster = hcat((collect.(split(monster, "\n")))...) |> rotr90
mask = BitArray(map(c -> c == '#', hcat(monster)))

# imshow(.!(mask))

function countmonsters(big)
    s = 0
    masksz = size(mask) .- 1
    max = size(big) .- masksz
    for i ∈ 1:max[1]
        for j ∈ 1:max[2]
            masked = big[i:i + masksz[1],j:j + masksz[2]]
            if all(masked .≥ mask)
                s += 1
            end
        end
    end
    s
end

count = maximum(map(countmonsters, [b for b ∈ D8(bigimg)]))
sum(bigimg) - count * sum(mask)


text:
using Gadfly

# wide-form plotting of heterogenous arrays of arrays

set_default_plot_size(4inch, 4inch)

xs = Any[Float32[0.547994, 0.819285, 0.567737, 0.557336, 0.27934, 0.777828, 0.8135, 0.00389743, 0.699683, 0.0958536],
         [0.951916, 0.999905, 0.251662, 0.986666, 0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512, 0.0769509],
         Float32[0.625995, 0.902331, 0.10869, 0.735902, 0.908017]]
plot(xs, x=Row.index, y=Col.value, color=Col.index, Scale.color_discrete, Geom.line)


text:
include("pif3d.jl");

Np=Int(1e5) #Number of particles

using BenchmarkTools
#Initial poisson solve
# @benchmark rho=accum_osde2(pm,zn[:,1],zn[:,2],zn[:,3],zn[:,4])
#
# b=@belapsed rho=accum_osde($pm,$zn[:,1],$zn[:,2],$zn[:,3],$zn[:,4])
#
# time(b)

x1n=rand(Float64,Np)
x2n=rand(Float64,Np)
x3n=rand(Float64,Np)
x4n=rand(Float64,Np)

pm=pif3d{Float64}(1,1,1,1.0)
rho=Array(Float64,1,1,1)

mem=Array(Float64,0) #bytes
walltime=Array(Float64,0) #nano seconds (10^-9 s)
dofs=Array(Float64,0)
idx=2
for idx=0:11
  Nx1=2^(idx)
  Nx2=2
  Nx3=0
  #Particle mesh coupling
  pm=pif3d{Float64}(Nx1,Nx2,Nx3,1.0)
  rho=Array(Float64, pm.N1,pm.N2,pm.N3)

  result=@benchmark rho=accum_osde(pm,x1n,x2n,x3n,x4n)

  # @benchmark rho=accum_osde($pm,$zn[:,1],$zn[:,2],$zn[:,3],$zn[:,4])
  push!(walltime,minimum(result).time)
  push!(mem, minimum(result).memory)
  push!(dofs, pm.N1*pm.N2*pm.N3)

  print(result,"\n")
  print("dofs:", dofs[end],", memory:", mem[end],"\n")

end

pm=pif3d{Float64}(1,0,0,1.0)

using PyPlot
plot(mem/1e3,walltime/1e6,"-o")
plot(64*[1., 1.], [walltime[1] ,walltime[end]]/1e6)
xlabel("memory [kB]")
ylabel("wall time [ms]")
grid()

# loglog(dofs,walltime/1e6,"-o")

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
using Knet, Base.Test

# Uncomment these if you want lots of messages:
# import Base.Test: default_handler, Success, Failure, Error
# default_handler(r::Success) = info("$(r.expr)")
# default_handler(r::Failure) = warn("FAIL: $(r.expr)")
# default_handler(r::Error)   = warn("$(r.err): $(r.expr)")

load_only = true
isapprox3(a,b,c)=all(map((x,y,z)->isapprox(x,y;rtol=z), a,b,c))

include(Knet.dir("examples","linreg.jl"))
@test LinReg.main("") == (0.0005497846637255735,32.77257400591496,0.11265426200775067)

include(Knet.dir("examples","mnist2d.jl"))
@test MNIST2D.main("--epochs 1") == (0.3204898f0,32.93997f0,4.614684f0)
@test MNIST2D.main("--epochs 1 --ysparse") == (0.3204898f0,32.93997f0,4.614684f0)
@test isapprox3(MNIST2D.main("--epochs 1 --xsparse"), (0.32048982f0,32.93997f0,4.6146836f0), (0.01,0.001,0.1))
@test isapprox3(MNIST2D.main("--epochs 1 --xsparse --ysparse"), (0.3204898f0,32.93997f0,4.614684f0), (0.01,0.001,0.1))

include(Knet.dir("examples","mnist4d.jl"))
@test isapprox3(MNIST4D.main("--epochs 1"), (0.20520441331629022,65.88371276855469,114.9615707397461), (0.02,0.001,0.2))

include(Knet.dir("examples","adding.jl"))
@test Adding.main("--epochs 1") == (0.22713459f0,3.3565507f0,5.3267756f0)
@test Adding.main("--epochs 1 --nettype lstm") == (0.24768005f0,3.601481f0,1.2290705f0)

seqdata = Knet.dir("data","seqdata.txt")

include(Knet.dir("examples","rnnlm.jl"))
@test RNNLM.main("--max_max_epoch 1 --dense $seqdata") == (30.65637197763184,110.81809997558594,29.36899185180664)

include(Knet.dir("examples","copyseq.jl"))
@test CopySeq.main("--epochs 1 --dense $seqdata") == (40.00286169097269,30.352935791015625,1.646486520767212)

### DEAD CODE:
# using Compat
# using Knet
# @date include("testdense.jl")
# @date include("testsparse.jl")
# @date include("testconvert.jl")
# @date include("testlinalg.jl")
# @date include("testcolops.jl")
# @date include("testmnist.jl")
# @date include("testlayers.jl")
# @date include("testperceptron.jl")
# @date include("testkperceptron.jl")
# # @date include("tutorial.jl")


text:
#=
multiply.jl
Multiplication of a LinearMapAX object with other things
2020-06-16 add 5-arg mul!
2020-06-27 Jeff Fessler
=#

export mul!

using LinearMaps
using LinearAlgebra: UniformScaling, I
import LinearAlgebra: mul! #, lmul!, rmul!
#import LinearAlgebra # AdjointAbsVec

# multiply with I or s*I (identity or scaled identity)
Base.:(*)(A::LinearMapAX, B::UniformScaling{Bool}) = B.λ ? A : (A * B.λ)
Base.:(*)(A::LinearMapAX, B::UniformScaling) = (B.λ == 1) ? A : (A * B.λ)
Base.:(*)(B::UniformScaling{Bool}, A::LinearMapAX) = B.λ ? A : (B.λ * A)
Base.:(*)(B::UniformScaling, A::LinearMapAX) = (B.λ == 1) ? A : (B.λ * A)

# multiply with scalars
Base.:(*)(s::Number, A::LinearMapAX) =
    LinearMapAA((s*I) * A._lmap ; prop=A._prop, idim=A._idim, odim=A._odim,
        operator = isa(A, LinearMapAO))
Base.:(*)(A::LinearMapAX, s::Number) =
    LinearMapAA(A._lmap * (s*I) ; prop=A._prop, idim=A._idim, odim=A._odim,
        operator = isa(A, LinearMapAO))

# multiply LMAX objects, if compatible
Base.:(*)(A::LinearMapAO{Ta,Do,D}, B::LinearMapAO{Tb,D,Di}) where {Ta,Tb,Do,D,Di} =
    lm_obj_mul(A, B)
Base.:(*)(A::LinearMapAO{T,Do,1}, B::LinearMapAM) where {T,Do} = lm_obj_mul(A, B)
Base.:(*)(A::LinearMapAM, B::LinearMapAO{T,1,Di}) where {T,Di} = lm_obj_mul(A, B)
Base.:(*)(A::LinearMapAM, B::LinearMapAM) = lm_obj_mul(A, B)

function lm_obj_mul(A::LinearMapAX, B::LinearMapAX)
    (A._idim != B._odim) && throw("dim mismatch")
    LinearMapAA(A._lmap * B._lmap ;
        prop = (prod=(A._prop,B._prop),),
        idim = B._idim, odim = A._odim,
        operator = isa(A, LinearMapAO) | isa(B, LinearMapAO),
    )
end

# multiply with a matrix
# subtle: AM * Matrix and Matrix * AM yield a new AM
# whereas AO * M and M * AO yield a matrix of numbers!
function AM_M(A::LinearMapAM, B::AbstractMatrix)
    (A._idim == (size(B,1),)) || throw("$(A._idim) * $(size(B,1)) mismatch")
    LinearMapAA(A._lmap * LinearMap(B) ; prop=A._prop, odim=A._odim)
end
function M_AM(A::AbstractMatrix, B::LinearMapAM)
    (B._odim == (size(A,2),)) || throw("$(B._odim) * $(size(A,1)) mismatch")
    LinearMapAA(LinearMap(A) * B._lmap ; prop=B._prop, idim=B._idim)
end
Base.:(*)(A::LinearMapAM, B::AbstractMatrix) = AM_M(A, B)
Base.:(*)(A::AbstractMatrix, B::LinearMapAM) = M_AM(A, B)

# LMAM case is easy!

# multiply with vectors (in-place)
# pass to lmam_mul! to handle composite maps (products) effectively
# 5-arg mul! requires julia 1.3 or later
# mul!(y, A, x, α, β) ≡ y .= A*(α*x) + β*y

mul!(y::AbstractVector, A::LinearMapAM, x::AbstractVector) =
    lm_mul!(y, A._lmap, x, 1, 0)
#    mul!(y, A._lmap, x)

mul!(y::AbstractVector, A::LinearMapAM, x::AbstractVector, α::Number, β::Number) =
    lm_mul!(y, A._lmap, x, α, β)
#    mul!(y, A._lmap, x, α, β)

# treat LinearMaps.CompositeMap as special case for in-place operations
function lm_mul!(y::AbstractVector, Lm::LinearMaps.CompositeMap,
    x::AbstractVector, α::Number, β::Number)
    LinearMaps.mul!(y, Lm, x, α, β) # todo: composite buffer
end

# 5-arg mul! for any other type
lm_mul!(y::AbstractVector, Lm::LinearMap,
    x::AbstractVector, α::Number, β::Number) =
    LinearMaps.mul!(y, Lm, x, α, β)

# with array
#=
these are unused because AM * array becomes a new AM
mul!(Y::AbstractArray, A::LinearMapAM, X::AbstractArray, α::Number, β::Number) =
    lmao_mul!(Y, A._lmap, X, α, β ; idim=A._idim, odim=A._odim)

mul!(Y::AbstractArray, A::LinearMapAM, X::AbstractArray) =
    LinearMapsAA.mul!(Y, A, X, 1, 0)
=#

# LMAO case

# 3-arg O*X
mul!(Y::AbstractArray, A::LinearMapAO, X::AbstractArray) =
    mul!(Y, A, X, 1, 0)

# 3-arg X*O
mul!(Y::AbstractArray, X::AbstractArray, A::LinearMapAO) =
    mul!(Y, X, A, 1, 0)

# 5-arg O*X
mul!(Y::AbstractArray, A::LinearMapAO, X::AbstractArray, α::Number, β::Number) =
    lmao_mul!(Y, A._lmap, X, α, β ; idim=A._idim, odim=A._odim)

# 5-arg X*O (todo: test complex case)
mul!(Y::AbstractArray, X::AbstractArray, A::LinearMapAO, α::Number, β::Number) =
    lmao_mul!(Y, A._lmap', X, α, β ; idim=A._odim, odim=A._idim) # note!

"""
     lmao_mul!(Y, A, X, α, β ; idim, odim)

Core routine for 5-arg multiply.
If `A._idim = (2,3,4)` and `A._odim = (5,6)` and
if input `X` has size `(2,3,4, 7,8)`
then output `Y` will have size `(5,6, 7,8)`
"""
function lmao_mul!(Y::AbstractArray, Lm::LinearMap, X::AbstractArray,
    α::Number, β::Number ;
    idim = (size(Lm,2),),
    odim = (size(Lm,1),),
)

    Di = length(idim)
    Do = length(odim)
    (Di > ndims(X) || (idim != size(X)[1:Di])) &&
        throw("idim=$(idim) vs size(RHS)=$(size(X))")
    (Do > ndims(Y) || (odim != size(Y)[1:Do])) &&
        throw("odim=$(odim) vs size(LHS)=$(size(Y))")
    size(X)[(Di+1):end] == size(Y)[(Do+1):end] ||
        throw("size(LHS)=$(size(Y)) vs size(RHS)=$(size(X))")

    x = reshape(X, prod(idim), :)
    y = reshape(Y, prod(odim), :)
    K = size(x,2)
    size(y,2) == K || throw("mismatch $(size(y,2)) K=$K")

    for k=1:K
        xk = selectdim(x,2,k)
        yk = selectdim(y,2,k)
        lm_mul!(yk, Lm, xk, α, β)
    end
    return Y
end

# multiply by array, with allocation

Base.:(*)(A::LinearMapAO, X::AbstractArray) = lmax_mul(A, X)
Base.:(*)(X::AbstractArray, A::LinearMapAO) = lmax_mul(A', X) # note!

#=
# this next line caused ambiguous method errors:
# Base.:(*)(A::LinearMapAM, X::AbstractArray) = lmax_mul(A, X)

# so i resort to this awful kludge:
Base.:(*)(A::LinearMapAM, X::AbstractArray{T,4}) where {T} = lmax_mul(A, X)

nah, too difficult, so revert to the AM*M returning an object, per above
=#

function lmax_mul(A::LinearMapAX{T}, X::AbstractArray) where {T}
    Di = length(A._idim)
    Do = length(A._odim)
    (Di > ndims(X) || (A._idim != size(X)[1:Di])) &&
         throw("idim=$(A._idim) vs size(RHS)=$(size(X))")
    extra = size(X)[(Di+1):end]
    Ty = promote_type(T, eltype(X))
    Y = Array{T}(undef, A._odim..., extra...) # allocate
    lmao_mul!(Y, A._lmap, X, 1, 0; idim=A._idim, odim=A._odim)
end

# multiply with vector

# O*v
Base.:(*)(A::LinearMapAO{T,Do,1}, v::AbstractVector) where {T,Do} =
    reshape(A._lmap * v, A._odim)
# u'*O (no, use general X*O above because unclear what this would mean)
#Base.:(*)(u::LinearAlgebra.AdjointAbsVec, A::LinearMapAO) =
#    reshape(A._lmap' * u', A._idim)

# A*v
Base.:(*)(A::LinearMapAM, v::AbstractVector) =
    A._lmap * v
# u'*A (nah, not worth it)
#Base.:(*)(u::LinearAlgebra.AdjointAbsVec, A::LinearMapAM) =
#    (A._lmap' * u')'

#= bad kludge
LMAAmanyFromOne = Union{
    LinearMapAA{T,2,1},
    LinearMapAA{T,3,1},
    LinearMapAA{T,4,1},
    } where {T}

Base.:(*)(A::LMAAmanyFromOne, v::AbstractVector{<:Number}) where {T} =
    reshape(A._lmap * v, A._odim)
Base.:(*)(A::LinearMapAA, v::AbstractVector{<:Number}) where {T} =
    reshape(A._lmap * v, A._odim)
=#

#= these are pointless; see multiplication with scalars above
lmul!(s::Number, A::LinearMapAA) = lmul!(s, A._lmap)
rmul!(A::LinearMapAA, s::Number) = rmul!(A._lmap, s)
=#


text:
# # Opinion spread

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../opinion.mp4" type="video/mp4">
# </video>
# ```

# This is a simple model of how an opinion spreads through a community.
# Each individual has a number of opinions as a list of integers.
# They can change their opinion by changing the numbers in the list.

# Agents can change their opinion at each step.
# They choose one of their neighbors randomly, and adopt one of the neighbor's opinion.
# They are more likely to adopt their neighbors opinion if the share more opinions with each other.

using Agents
using InteractiveDynamics # plotting agents
using CairoMakie # for static plotting
CairoMakie.activate!() # hide
using Random # hide

# ## Building the model
# ### 1. Model creation

mutable struct Citizen <: AbstractAgent
    id::Int
    pos::Dims{2}
    stabilized::Bool
    opinion::Array{Int,1}
    prev_opinion::Array{Int,1}
end

function create_model(; dims = (10, 10), nopinions = 3, levels_per_opinion = 4)
    space = GridSpace(dims)
    properties = Dict(:nopinions => nopinions)
    model = AgentBasedModel(
        Citizen,
        space,
        scheduler = Schedulers.randomly,
        properties = properties,
    )
    for pos in positions(model)
        add_agent!(
            pos,
            model,
            false,
            rand(model.rng, 1:levels_per_opinion, nopinions),
            rand(model.rng, 1:levels_per_opinion, nopinions),
        )
    end
    return model
end

# ### 2. Stepping functions

function adopt!(agent, model)
    neighbor = rand(model.rng, collect(nearby_ids(agent, model)))
    matches = model[neighbor].opinion .== agent.opinion
    nmatches = count(matches)
    if nmatches < model.nopinions && rand(model.rng) < nmatches / model.nopinions
        switchId = rand(model.rng, findall(x -> x == false, matches))
        agent.opinion[switchId] = model[neighbor].opinion[switchId]
    end
end

function update_prev_opinion!(agent, model)
    for i in 1:(model.nopinions)
        agent.prev_opinion[i] = agent.opinion[i]
    end
end

function is_stabilized!(agent, model)
    if agent.prev_opinion == agent.opinion
        agent.stabilized = true
    else
        agent.stabilized = false
    end
end

function agent_step!(agent, model)
    update_prev_opinion!(agent, model)
    adopt!(agent, model)
    is_stabilized!(agent, model)
end

# ## Running the model

# First, we create a stopping condition, which runs the model until all agents stabilize.

rununtil(model, s) = count(a -> a.stabilized, allagents(model)) == length(positions(model))

# Then we create our model, run it and collect some information

model = create_model(nopinions = 3, levels_per_opinion = 4)

agentdata, _ = run!(model, agent_step!, dummystep, rununtil, adata = [(:stabilized, count)])

# ## Plotting

# The plot shows the number of stable agents, that is, number of agents whose opinions
# don't change from one step to the next. Note that the number of stable agents can
# fluctuate before the final convergence.

f = Figure(resolution = (600, 400))
ax =
    f[1, 1] = Axis(
        f,
        xlabel = "Generation",
        ylabel = "# of stabilized agents",
        title = "Population Stability",
    )
lines!(ax, 1:size(agentdata, 1), agentdata.count_stabilized, linewidth = 2, color = :blue)
f

# ### Animation

# Here is an animation that shows change of agent opinions over time.
# The first three opinions of an agent determines its color in RGB.
Random.seed!(648) # hide
levels_per_opinion = 3
ac(agent) = CairoMakie.RGB((agent.opinion[1:3] ./ levels_per_opinion)...)
model = create_model(nopinions = 3, levels_per_opinion = levels_per_opinion)

abm_video(
    "opinion.mp4",
    model,
    agent_step!;
    ac = ac,
    am = '■',
    as = 20,
    framerate = 20,
    frames = 265,
    title = "Opinion Spread",
)
nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../opinion.mp4" type="video/mp4">
# </video>
# ```


text:
using MAT
using MinMaxFilter

# ######
#
# 2D
#
# ######

filen = matopen(joinpath(dirname(@__FILE__), "data", "2d_array.mat"))
A = read(filen, "Z")
A = convert(Array{FloatingPoint}, A)
close(filen)

minval, maxval = minmax_filter(A, 11)

matching = A[2:size(maxval)[1]+1, 2:size(maxval)[2]+1]
matching = matching .== maxval

peaks = maxval[matching]

for l = 1:length(peaks)
    x, y = ind2sub(size(A), find(A .== peaks[l]))
    println("x=$x, y=$y, val=$(A[x,y])")
end

# ######
#
# 3D
#
# #####

using Winston
using EEG

filen = matopen(joinpath(dirname(@__FILE__), "data", "3d_array.mat"))
A = read(filen, "image")
A = convert(Array{FloatingPoint}, A)
close(filen)

p = plot_dat(A)

minval, maxval = minmax_filter(A, 6)

matching = A[2:size(maxval)[1]+1, 2:size(maxval)[2]+1, 2:size(maxval)[3]+1]
matching = matching .== maxval

peaks = maxval[matching]
peaks = peaks[peaks .>= 0.1 * maximum(peaks)]

for l = 1:length(peaks)
    x, y, z = ind2sub(size(A), find(A .== peaks[l]))
    println("x=$x, y=$y, z=$z, val=$(A[x,y,z])")
    oplot(p, x, y, z)
end

file(p, "images/3D-dipole.png", width=1000, height=1000)


text:
module Ptrace
  using ..Gallium
  using ..Gallium: RemotePtr, X86_64
  using ..Gallium.Registers: RegisterValue
  using ObjFileBase
  import Base: start, next, done, iteratorsize

  include(joinpath(dirname(@__FILE__),"constants","ptrace.jl"))

  immutable Session
      pid::Cint
      mem_fd::Base.RawFD
  end

  immutable MapsIterator
      file::String
  end
  start(x::MapsIterator) = open(x.file)
  function next(x::MapsIterator, file)
      line = readline(file)
      parts = split(line, ' ', keep = false)
      range = UnitRange{UInt64}(parse.([UInt64],split(parts[1],'-'),16)...)
      # Our ranges are [,], kernel ones are [,)
      range = first(range):last(range)-1
      fname = length(parts) == 6 ? strip(parts[6]) : ""
      ((range, fname, line), file)
  end
  done(x::MapsIterator, file) = eof(file)
  iteratorsize(::Type{MapsIterator}) = Base.SizeUnknown()

  function waitpid(s::Session)
    status = Ref{Cint}(0)
    rpid = ccall(:waitpid, Cint, (Cint, Ptr{Cint}, Cint), s.pid, status, 0)
    @assert rpid == s.pid
  end

  function ptrace(req, pid, addr, data)
    Libc.errno(0) # Make sure to set errno to 0, so we can detect failure
    res = ccall(:ptrace, Clong, (Cint, Cint, Ptr{Void}, Ptr{Void}),
      req, pid, addr, data)
    Base.systemerror("ptrace", res == -1 && Libc.errno() != 0)
    (res % Culong)
  end

  function Gallium.read_exe(s::Session)
      readmeta(IOBuffer(open(Base.Mmap.mmap, "/proc/$pid/exe")))
  end

  function attach(pid)
    pid = convert(Cint, pid)
    ptrace(PTRACE_ATTACH, pid, C_NULL, C_NULL)
    sess = Session(pid, RawFD(
      ccall(:open, Cint, (Ptr{UInt8}, Cint), "/proc/$pid/mem", Base.JL_O_RDWR)))
    imageh = Gallium.read_exe(sess)
    modules = Gallium.GlibcDyldModules.load_library_map(sess, imageh)
    sess, modules
  end

  function Gallium.load{T}(s::Session, ptr::RemotePtr{T})
    r = Ref{T}()
    Base.systemerror("pread(/proc/mem)",
      ccall(:pread, Cssize_t, (Cint, Ptr{Void}, Csize_t, UInt64),
        s.mem_fd, r, sizeof(T), convert(UInt64, ptr)) == -1)
    r[]
  end

  function Gallium.load{T}(s::Session, ptr::RemotePtr{T}, count)
    r = Vector{T}(count)
    Base.systemerror("pread(/proc/mem)",
      ccall(:pread, Cssize_t, (Cint, Ptr{Void}, Csize_t, UInt64),
        s.mem_fd, r, sizeof(r), convert(UInt64, ptr)) == -1)
    r
  end

  function Gallium.mapped_file(s::Session, ptr::RemotePtr)
    first(filter(map->UInt64(ptr) ∈ map[1], MapsIterator("/proc/$(s.pid)/maps")))[2]
  end

  immutable iovec
      iov_base::Ptr{Void}
      iov_len::Csize_t
  end

  const NT_PRSTATUS = 1
  function Gallium.getregs(s::Session)
    regs = Array(UInt64, length(X86_64.kernel_order))
    iov = Ref{iovec}(iovec(pointer(regs), sizeof(regs)))
    ptrace(PTRACE_GETREGSET, s.pid, Ptr{Void}(NT_PRSTATUS), iov)
    @assert iov[].iov_len == sizeof(regs)
    rregs = X86_64.BasicRegs()
    for i in X86_64.basic_regs
      setfield!(rregs, i+1,
        RegisterValue{UInt64}(
          regs[X86_64.kernel_numbering[X86_64.dwarf_numbering[i]]]))
    end
    rregs
  end

  function Gallium.continue!(s::Session)
    ptrace(PTRACE_CONT, s.pid, C_NULL, C_NULL)
    waitpid(s)
  end

  function Gallium.single_step!(s::Session)
    ptrace(PTRACE_SINGLESTEP, s.pid, C_NULL, C_NULL)
    waitpid(s)
  end

end


text:
module StatPlots

using Reexport
@reexport using Plots
import Plots: _cycle
using StatsBase
using Distributions
using DataFrames

import KernelDensity
import Loess
@recipe f(k::KernelDensity.UnivariateKDE) = k.x, k.density
@recipe f(k::KernelDensity.BivariateKDE) = k.x, k.y, k.density

@shorthands cdensity

export groupapply
export get_groupederror

include("dataframes.jl")
include("corrplot.jl")
include("cornerplot.jl")
include("distributions.jl")
include("boxplot.jl")
include("violin.jl")
include("hist.jl")
include("marginalhist.jl")
include("bar.jl")
include("shadederror.jl")
include("groupederror.jl")

end # module


text:
wider_js_type(::Type{<:Integer}) = Float64  # Javascript thinks everything is a `double`
wider_js_type(::Type{Float64}) = Float64
wider_js_type(x) = x

function _setprop!(clip::AnimationClip, frame::Integer, prop::AbstractString, jstype::AbstractString, value)
    T = wider_js_type(typeof(value))
    track = get!(clip.tracks, prop) do
        AnimationTrack{T}(prop, jstype)
    end
    insert!(track, frame, value)
    return nothing
end

function getclip!(animation::Animation, path::Path)
    get!(animation.clips, path) do
        AnimationClip(fps=animation.default_framerate)
    end
end

js_quaternion(m::AbstractMatrix) = js_quaternion(RotMatrix(SMatrix{3, 3, eltype(m)}(m)))
js_quaternion(q::Quat) = [q.x, q.y, q.z, q.w]
js_quaternion(::UniformScaling) = js_quaternion(Quat(1., 0., 0., 0.))
js_quaternion(r::Rotation) = js_quaternion(Quat(r))
js_quaternion(tform::Transformation) = js_quaternion(transform_deriv(tform, SVector(0., 0, 0)))

function js_scaling(tform::AbstractAffineMap)
    m = transform_deriv(tform, SVector(0., 0, 0))
    SVector(norm(SVector(m[1, 1], m[2, 1], m[3, 1])),
            norm(SVector(m[1, 2], m[2, 2], m[3, 2])),
            norm(SVector(m[1, 3], m[2, 3], m[3, 3])))
end

js_position(t::Transformation) = convert(Vector, t(SVector(0., 0, 0)))

Cassette.@context AnimationCtx

function Cassette.overdub(ctx::AnimationCtx, ::typeof(settransform!), vis::Visualizer, tform::Transformation)
    animation, frame = ctx.metadata
    clip = getclip!(animation, vis.path)
    _setprop!(clip, frame, "scale", "vector3", js_scaling(tform))
    _setprop!(clip, frame, "position", "vector3", js_position(tform))
    _setprop!(clip, frame, "quaternion", "quaternion", js_quaternion(tform))
end

function Cassette.overdub(ctx::AnimationCtx, ::typeof(setprop!), vis::Visualizer, prop::AbstractString, value)
    animation, frame = ctx.metadata
    clip = getclip!(animation, vis.path)
    _setprop!(clip, frame, prop, "number", value)
end

function Cassette.overdub(ctx::AnimationCtx, ::typeof(setprop!), vis::Visualizer, prop::AbstractString, jstype::AbstractString, value)
    animation, frame = ctx.metadata
    clip = getclip!(animation, vis.path)
    _setprop!(clip, frame, prop, jstype, value)
end

function atframe(f, animation::Animation, frame::Integer)
    Cassette.overdub(AnimationCtx(metadata=(animation, frame)), f)
    return animation
end

function atframe(f::Function, anim::Animation, path::Path, frame::Integer)
    error("""
    atframe(f::Function, anim::Animation, path::Path, frame::Integer) is no longer supported.
    Please see the updated animation example notebook.
    """)
end

function atframe(f::Function, anim::Animation, vis::Visualizer, frame::Integer)
    Base.depwarn("""
    atframe(f::Function, anim::Animation, vis::Visualizer, frame::Integer) is deprecated.
    Please use atframe(g, anim, frame) instead, where g is similar to f but takes
    no arguments and should call methods on vis.
    See also the updated animation example notebook.
    """, :atframe)
    atframe(() -> f(vis), anim, frame)
end


text:
struct DDEProblem{uType,tType,lType,lType2,isinplace,P,F,H,C} <:
                          AbstractDDEProblem{uType,tType,lType,isinplace}
  f::F
  u0::uType
  h::H
  tspan::tType
  p::P
  constant_lags::lType
  dependent_lags::lType2
  callback::C
  neutral::Bool
  @add_kwonly function DDEProblem{iip}(f::AbstractDDEFunction{iip},
                                 u0,h,tspan,p=nothing;
                                 constant_lags=[],
                                 dependent_lags=[],
                                 neutral = f.mass_matrix == I ?
                                           false : det(f.mass_matrix)!=1,
                                 callback = nothing) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),
               typeof(constant_lags),typeof(dependent_lags),
               isinplace(f),typeof(p),
               typeof(f),typeof(h),typeof(callback)
               }(f,u0,h,_tspan,p,
                 constant_lags,
                 dependent_lags,callback,
                 neutral)
  end

  function DDEProblem{iip}(f,u0,h,tspan,p=nothing;kwargs...) where {iip}
    DDEProblem(convert(DDEFunction{iip},f),u0,h,tspan,p;kwargs...)
  end

end

function DDEProblem(f::AbstractDDEFunction,u0,h,tspan,p=nothing;kwargs...)
  DDEProblem{isinplace(f)}(f,u0,h,tspan,p;kwargs...)
end

function DDEProblem(f,u0,h,tspan,p=nothing;kwargs...)
  DDEProblem(convert(DDEFunction,f),u0,h,tspan,p;kwargs...)
end


text:
using JuliaVariables
using Test
using NameResolution
using MLStyle

rmlines(ex::Expr) = begin
    hd = ex.head
    tl = map(rmlines, filter(!islinenumbernode, ex.args))
    Expr(hd, tl...)
end
rmlines(@nospecialize(a)) = a
islinenumbernode(@nospecialize(x)) = x isa LineNumberNode

function first_scope_except_top(ex)
    inner = ex.args[2]
    find(ex) = @match ex begin
        Expr(:scoped, scope, inner) => (scope, inner)
        Expr(head, args...) =>
            for each in args
                res = find(each)
                res !== nothing && return res
            end
        _ => nothing
    end
    find(inner)
end

get_fn_scope(ex) =
    @match ex begin
        Expr(
            :scoped,
            _,
            Expr(:function, _,
                Expr(:scoped, scope, _)
            )
        ) => scope
        _ => error("malformed $ex")
    end

get_fn_body(ex) =
    @match ex begin
        Expr(
            :scoped,
            _,
            Expr(:function, _,
                Expr(:scoped, _, body)
            )
        ) => body
        _ => error("malformed $ex")
    end

@testset "JuliaVariables.jl" begin
    func = solve(:(function f(x)
        let y = x + 1
            y
        end
    end))
    println(func |> rmlines)

    func = solve(:(function f(x)
        y = x + 1
        z -> z + y
    end))

    println(func |> rmlines)

    func = solve(:(function f(x)
        y = x + 1
        let y = y + 1
            for z in 1:10
                x + y + z
            end
        end
    end))
    println(func |> rmlines)

    func = solve(
        macroexpand(@__MODULE__,:(
        @inline function f(x)
            y = x + 1
            let y = y + 1
                (x + y  + z for z in 1:10)
            end
        end
    )))
    println(func |> rmlines)
    # @test map(haskey(func.scope.freevars, _), [])

    func = solve_from_local(quote
        a -> a + z
    end)
    println(func |> rmlines)

    func = solve_from_local(quote
    q = a
    function z(x, k=1)
        (x=x, k=k, q=q)
    end
    end)
    println(func |> rmlines)
    # Write your own tests here.

    a = solve(:(function k(x::T) where T
                (1, T)
              end
    ))

    @test any(x -> x.name == :T, get_fn_scope(a).bounds)

    # a = solve(:(2 .^ [2, 3]))
    # @test eval(a) == [4, 8]

    a = solve(:(function z(x, k=1)
                   x + 20 + a + k
               end
    ))

    @test any(x -> x.name == :k, get_fn_scope(a).bounds)

     a = solve_from_local(quote
            x  = 1
            y  = 2
            .^ = 2
            function z()
                   x .^ y
            end
    end)
    scope = first_scope_except_top(a)[1]
    @test length(scope.freevars) == 2

    a = solve(:(function z(x, k=1)
                   (k=k, )
               end
    ))

    @test any(x -> x.name == :k, get_fn_scope(a).bounds)

    @test @when :(k=$_, ) = get_fn_body(a).args[2] begin
        true
    @otherwise
        false
    end

    func = solve(:(function (x)
              z = x + 1
              y -> begin
                  z += 1
                  z
              end
          end))

    println(func |> rmlines)
    @test any(x -> x.name == :z && x.is_mutable, get_fn_scope(func).bounds)

    func = solve(macroexpand(@__MODULE__, :(@inline function (x)
              z = x + 1
              @inline function (y)
                  z += 1
                  z
              end
          end)))
    println(func |> rmlines)
end

@testset "uninitialized bound variable" begin
        ast = quote
        let x
            x = 5
            function ()
                x = 2
            end
        end
        end |> solve
        scope = first_scope_except_top(ast)[1]
        @test isempty(scope.bound_inits)
    end

@testset "inplace binary" begin
    ast = quote
        local x
        function ()
            x ^= 5
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    var = scope.freevars[1]
    @test var.name == :x && var.is_mutable
end

@testset "type variables" begin
    ast = quote
        function (::T) where T <: Number
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    var = scope.bounds[1]
    @test var.name == :T

    ast = quote
        function (::T) where {T <: Number, A <: T}
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    # A is not used inside function scope
    @test Set([e.name for e in scope.bounds]) == Set([:T])

    ast = quote
        function (::A) where {T <: Number, T <: A <: GlobalType}
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test Set([e.name for e in scope.bounds]) == Set([:T, :A])
end

@testset "default argument is free variables" begin
    ast = quote
        z = 2
        function (x::T = z) where T <: Number
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    var = scope.freevars[1]
    @test var.name == :z
end

@testset "global mark" begin
    ast = quote
        x = 5
        function ()
            global x
            x = 2
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test isempty(scope.bounds) && isempty(scope.freevars)
end

@testset "while" begin
    ast = quote
        x = 5
        while cond
            x = 10
        end
    end |> solve
    scope = first_scope_except_top(ast)[1]
    x = scope.bounds[1]
    @test x.name === :x
    @test x.is_mutable === false
    @test x.is_global === false
end

@testset "let" begin
    ast = quote
        x = 5
        let x
            x = 3
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test isempty(scope.freevars)

    ast = quote
        x = 5
        let x = 3
            x
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test isempty(scope.freevars)

    ast = quote
        function ()
            x = 5
            let
                y = 1
            end
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test scope.bounds[1].name == :x
    @test all(scope.freevars) do var
        var.name !== :y
    end
    @test all(scope.bounds) do var
        var.name !== :y
    end
end

@testset "namedtuple" begin
    ast = quote
        function (c)
            (x = 3, )
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test [e.name for e in scope.bounds] == [:c]

    ast = quote
        function (c)
            (x = 3, y=2; z=3)
        end
    end |> solve_from_local
    scope = first_scope_except_top(ast)[1]
    @test [e.name for e in scope.bounds] == [:c]
end

@testset "broadcast" begin
    ast = quote
        x = 1
        .^ = ^
        y = 2
        function()
            x .^ y
        end
    end |> solve_from_local
    tp = first_scope_except_top(ast)
    scope = tp[1]
    @test length(scope.freevars) == 2
    @test Set([e.name for e in scope.freevars]) == Set([:x, :y])
    @test @match tp[2].args[2] begin
        :($_ .^ $_) => true
        _ => false
    end
end

@testset "broadcast fusion" begin
    ast = quote
        x = [1, 2]
        function()
            x .= 1
        end
    end |> solve_from_local
    tp = first_scope_except_top(ast)
    scope = tp[1]
    @test length(scope.freevars) == 1
    @test Set([e.name for e in scope.freevars]) == Set([:x])
    x = scope.freevars[1]
    @test x.is_mutable == false
end

@testset "simplification" begin

    ast = quote
        for i = I, j = J
            let x = 1, y
                c(x) = begin 3 end
                f(x) do
                    1
                end
            end
        end
    end |> rmlines |> simplify_ex

    expec = quote
        for i in I
            for j in J
                let x = 1
                    let y
                        function c(x)
                            3
                        end
                        f(function () 1 end, x)
                    end
                end
            end
        end
    end |> rmlines
    @test string(ast) == string(expec)
end

@testset "no effect" begin
    ast = macroexpand(@__MODULE__(), quote
        function ()
            @label x
            @goto x
        end
    end) |> solve_from_local
    tp = first_scope_except_top(ast)
    scope = tp[1]
    @test isempty(scope.bounds)
    @test isempty(scope.freevars)

    ast = macroexpand(@__MODULE__(), quote
        @inline function ()
        end
    end) |> solve_from_local
    tp = first_scope_except_top(ast)
    scope = tp[1]
    @test isempty(scope.bounds)
    @test isempty(scope.freevars)
end


text:
#=

These are fast algorithms for maniuplating Julia sparse matrices.
It is part of the Laplacians package.

Started by Dan Spielman.

=#

""" Permute the rows and columns of a symmetric matrix, a little faster than the default
routine in Julia.  It basically does a radix sort on rows than columns. It does require `a`
to be symmetrix.  But, the same ideas could probably speed up many sparse matrix routines.
`sympermute(a, perm)` is the same as `a[perm,perm]`.
"""
function symPermuteCSC(a::SparseMatrixCSC{Tv, Ti}, perm::Vector{Ti}) where {Tv,Ti}

    numnz = nnz(a)

    Ti1 = Ti
    Tv1 = Tv

    permi = Vector{Ti1}(undef, a.n)
    permi[perm] = collect(1:a.n)

    J = Vector{Ti1}(undef, numnz)

    deg = Vector{Ti1}(undef, a.n)

    ptr = 1
    for col = 1:a.n
        deg[permi[col]] = a.colptr[col+1] - a.colptr[col]
        for k = a.colptr[col] : (a.colptr[col+1]-1)
            J[ptr] = permi[col]
            ptr += 1
        end
    end

    I = permi[a.rowval];

    I1 = Vector{Ti1}(undef, numnz)
    J1 = Vector{Ti1}(undef, numnz)
    V1 = Vector{Tv1}(undef, numnz)

    cumdeg = cumsum(deg)
    col = [1;cumdeg .+ 1]

    cumdeg1 = copy(cumdeg)

    for i in numnz:-1:1
        ptr = cumdeg1[I[i]]
        cumdeg1[I[i]] -= 1
        I1[ptr] = I[i]
        J1[ptr] = J[i]
        V1[ptr] = a.nzval[i]
    end

    I2 = Vector{Ti1}(undef, numnz)
    V2 = Vector{Tv1}(undef, numnz)

    for i in numnz:-1:1
        ptr = cumdeg[J1[i]]
        cumdeg[J1[i]] -= 1
        I2[ptr] = I1[i]
        V2[ptr] = V1[i]
    end

    aperm = SparseMatrixCSC(a.n, a.n, col, I2, V2)

    return aperm
end

"""Compute the transpose of a matrix with symmetric nonzero structure"""
function symTransposeCSC(a::SparseMatrixCSC{Tv, Ti}) where {Tv,Ti}

    numnz = nnz(a)

    cumdeg = a.colptr[2:end] .- 1

    I = a.rowval;

    V1 = Vector{Tv}(undef,numnz)

    for i in numnz:-1:1
        ptr = cumdeg[I[i]]
        cumdeg[I[i]] -= 1
        V1[ptr] = a.nzval[i]
    end

    return SparseMatrixCSC(a.n, a.n, copy(a.colptr),I, V1)

end

"""Given a set of integers, `set` between 1 and n, return a sorted version of them"""
function sortSet(set::Vector{Ti},n::Ti) where Ti
    v = zeros(Bool,n)
    for i in set
        v[i] = true
    end
    out = Vector{Ti}(undef, length(set))
    ptr = 1
    for i in 1:n
        if v[i]
            out[ptr] = i
            ptr += 1
        end
    end
    out
end

"""For a the submatrix of a with the entries of ijv indexed by list.  The list must be sorted.
This is equivalent to

~~~
(ai,aj,av) = findnz(a)
sparse(ai[list],aj[list],av[list],a.m,a.n)
~~~
"""
function submatrixCSC(a::SparseMatrixCSC{Tv, Ti}, list::Vector{Ti}) where {Tv,Ti}

    list = sortSet(list,nnz(a))

    colptr = Vector{Ti}(undef, 1 .+ a.n)
    colptr[1] = 1

    numnz = length(list)

    rowval = Vector{Ti}(undef, numnz)
    nzval = Vector{Tv}(undef, numnz)

    ptr = 1
    col = 1

    for ind in list

        while a.colptr[col+1] <= ind
            col += 1
            colptr[col] = ptr
        end

        rowval[ptr] = a.rowval[ind]
        nzval[ptr] = a.nzval[ind]
        ptr += 1

    end

    colptr[(col+1):end] .= ptr

    return SparseMatrixCSC(a.m, a.n, colptr, rowval, nzval)
end

#=
  This is a version of sympermute that is not quite as fast.
  But, it could be the basis of a faster algorithm.
  =#
#=
function sympermuteB{Tv,Ti}(a::SparseMatrixCSC{Tv, Ti}, perm::Array{Ti,1})

    numnz = nnz(a)

    Ti1 = Ti
    Tv1 = Tv

    permi = Array(Ti1, a.n)
    permi[perm] = collect(1:a.n)

    J = Array(Ti1, numnz)
    I = Array(Ti1, numnz)
    deg = Array(Ti1, a.n)

    ptr = 1
    for col = 1:a.n
        deg[permi[col]] = a.colptr[col+1] - a.colptr[col]
        for k = a.colptr[col] : (a.colptr[col+1]-1)
            J[ptr] = permi[col]
            I[ptr] = permi[a.rowval[k]]
            ptr += 1
        end
    end

    I1 = Array(Ti1, numnz)
    J1 = Array(Ti1, numnz)
    V1 = Array(Tv1, numnz)

    cumdeg = cumsum(deg)
    for i in numnz:-1:1
        ptr = cumdeg[I[i]]
        I1[ptr] = I[i]
        J1[ptr] = J[i]
        V1[ptr] = a.nzval[i]
        cumdeg[I[i]] -= 1
    end

    I2 = Array(Ti1, numnz)
    V2 = Array(Tv1, numnz)
    cumdeg = cumsum(deg)
    for i in numnz:-1:1
        ptr = cumdeg[J1[i]]
        I2[ptr] = I1[i]
        V2[ptr] = V1[i]
        cumdeg[J1[i]] -= 1
    end

    cumdeg = cumsum(deg)
    col = [1;cumdeg+1]
    aperm = SparseMatrixCSC(a.n, a.n, col, I2, V2)

    return aperm
end

=#


text:
"""
    SimpleAxis(v)

Povides an `AbstractAxis` interface for any `AbstractUnitRange`, `v `. `v` will
be considered both the `values` and `keys` of the return instance.

## Examples

A `SimpleAxis` is useful for giving a standard set of indices the ability to use
the filtering syntax for indexing.
```jldoctest
julia> using AxisIndices, StaticRanges

julia> x = SimpleAxis(2:10)
SimpleAxis(2:10)

julia> x[2]
2

julia> x[==(2)]
2

julia> x[2] == x[==(2)]  # keys and values are same
true

julia> x[>(2)]
SimpleAxis(3:10)

julia> x[>(2)]
SimpleAxis(3:10)

julia> x[1]
ERROR: BoundsError: attempt to access SimpleAxis(2:10) at index [1]
[...]
```
"""
struct SimpleAxis{I,Inds<:AbstractUnitRange{I}} <: AbstractAxis{I,Inds}
    parent::Inds

    function SimpleAxis{I,Inds}(inds::AbstractUnitRange) where {I,Inds}
        if inds isa Inds
            return new{I,Inds}(inds)
        else
            return SimpleAxis{I,Inds}(Inds(inds))
        end
    end

    SimpleAxis{I}(inds::AbstractUnitRange{I}) where {I} = SimpleAxis{I,typeof(inds)}(inds)
    SimpleAxis{I}(inds::AbstractUnitRange) where {I} = SimpleAxis{I}(AbstractUnitRange{I}(inds))
    function SimpleAxis{I}(inds::SimpleAxis) where {I}
        if eltype(inds) <: I
            return inds
        else
            return SimpleAxis{I}(parent(inds))
        end
    end

    SimpleAxis(inds::SimpleAxis) = inds
    SimpleAxis(inds::AbstractUnitRange{I}) where {I} = SimpleAxis{I}(inds)
    SimpleAxis(inds::IdentityUnitRange) = SimpleAxis(inds.indices)
    SimpleAxis() = new{Int,OneToMRange{Int}}(OneToMRange(0))

    """
        SimpleAxis(start::Integer, stop::Integer) -> SimpleAxis{UnitRange{Integer}}

    Passes `start` and `stop` arguments to `UnitRange` to construct the values of `SimpleAxis`.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> SimpleAxis(1, 10)
    SimpleAxis(1:10)
    ```
    """
    SimpleAxis(start::Integer, stop::Integer) = SimpleAxis(start:stop)

    """
        SimpleAxis(stop::Integer) -> SimpleAxis{Base.OneTo{Integer}}

    Passes `stop` `Base.OneTo` to construct the values of `SimpleAxis`.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> SimpleAxis(10)
    SimpleAxis(Base.OneTo(10))
    ```
    """
    SimpleAxis(stop::Integer) = SimpleAxis(StaticInt(1):stop)
end

@inline Base.getproperty(axis::SimpleAxis, k::Symbol) = getproperty(parent(axis), k)

ArrayInterface.unsafe_reconstruct(axis::SimpleAxis, inds; kwargs...) = SimpleAxis(inds)

# FIXME this should be deleted once https://github.com/SciML/ArrayInterface.jl/issues/79 is resolved
@propagate_inbounds function Base.getindex(axis::SimpleAxis, arg::StepRange{I}) where {I<:Integer}
    @boundscheck checkbounds(axis, arg)
    return maybe_unsafe_reconstruct(axis, arg)
end

function print_axis(io::IO, axis::SimpleAxis)
    print(io, "SimpleAxis($(Int(first(axis)):Int(last(axis))))")
end


text:
using StandardizedPredictors
using Documenter

DocMeta.setdocmeta!(StandardizedPredictors, :DocTestSetup, :(using StandardizedPredictors);
                    recursive=true)

makedocs(; modules=[StandardizedPredictors],
         authors="Beacon Biosignals, Inc.",
         repo="https://github.com/beacon-biosignals/StandardizedPredictors.jl/blob/{commit}{path}#{line}",
         sitename="StandardizedPredictors.jl",
         format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://beacon-biosignals.github.io/StandardizedPredictors.jl",
                                assets=String[]),
         pages=["Home" => "index.md",
                "API" => "api.md"])

deploydocs(; repo="github.com/beacon-biosignals/StandardizedPredictors.jl",
           devbranch="main",
           push_preview=true)


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Analytic Solutions.
#Written by F.X. Giraldo on 10/2003
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,time,icase,DFloat)

    #Initialize
    qe=zeros(DFloat,Npoin,1)
    ue=zeros(DFloat,2,Npoin)

    #Constants
    xmin=minimum(coord[1,:])
    xmax=maximum(coord[1,:])
    ymin=minimum(coord[2,:])
    ymax=maximum(coord[2,:])
    xl=xmax-xmin
    yl=ymax-ymin
    xm=0.5*(xmax+xmin)
    ym=0.5*(ymax+ymin)
    xc=xmin + 0.25*xl
    yc=ymin + 0.5*yl
    σ=32.0

    #Initial Condition
    for I=1:Npoin
        x=coord[1,I]
        y=coord[2,I]

        if (icase == 1)  #Gaussian in CCW direction
            ue[1,I]=+(y-ym);
            ue[2,I]=-(x-xm);
            xx=x - xc*cos(time) - yc*sin(time)
            yy=y + xc*sin(time) - yc*cos(time)
            qe[I]=exp( - σ*(xx^2 + yy^2 ) )
        end
    end #I

    return (qe,ue)

end


text:
# developed in Julia 1.0.1
# -----------------------------
# prepare memories for variables (not parameters)
# NOTE: using Dict which can be quite easily converted to a DataFrame (by DataFrames package)
# NOTE: requires "env" defined
# ==============================================================================

# ---------------------------------------------------- Vector data 序列数据（年）
# NOTE: use Dict to conveniently convert it to a DataFrame
# Data (by year)
Dt = Dict(
    :Year => Array( range(env.START_YEAR,length=env.T) ),  # Years  年份索引（从真实年份开始）
    # ---------------- Economic States (stocks) 经济状态（存量）
    :K => zeros(env.T),  # Urban capital supply  城镇资本存量供应
    :L => zeros(env.T),  # Urban labour supply  城镇劳动力供应
    # ---------------- Economic States (flows & prices) 经济状态（流量&价格）
    :Y => zeros(env.T),  # Urban output (GDP)  城镇产出
    :r  => zeros(env.T),  # Net interest rate (investment returns)  净投资利率
    :w̄ => zeros(env.T),  # Urban average wage level  城镇平均工资水平
    :C  => zeros(env.T),  # Aggregated consumption  总消费
    :I  => zeros(env.T),  # Investment  总投资
    :G  => zeros(env.T),  # Government purchase  政府购买
    # ---------------- Social Security Benefits 社会保障给付
    :Λ => zeros(env.T),  # Urban average pension benefits  城镇平均养老金发放量
    # ---------------- Fiscal 财政
    :D => zeros(env.T),  # Government outstanding debt  政府未偿债务余额 (verbose, not used)
    :D2Y => zeros(env.T),  # Government outstanding debt to GDP  政府未偿债务余额/GDP (verbose, not used)
    :TRc => zeros(env.T),  # Urban total consumption tax revenues  城镇总消费税收入
    :TRw => zeros(env.T),  # Urban total income tax revenues  城镇总收入税收入
    :LI  => zeros(env.T),  # the gap/surplus of the urban pooling medical account (UE-BMI)  城镇医保统筹账户缺口
    # ---------------- Welfare 福利
    :U => zeros(env.T),  # aggregated cross-sectional social utility 当年城镇居民效用之和
    # ---------------- Abstract & Intermediate 抽象变量和中间变量
    :o => zeros(env.T),  # wage scaling coefficient 工资放缩系数
    :𝕡 => zeros(env.T),  # transfer amount from firm contirbution to retired households 企业医保缴纳对退休人群个人账户的转移支付量
)

# ----------------------------------------------------- Matrix data (t*s) 矩阵数据（年*岁）
# Data (by year * age)
Dst = Dict(
    :Year => Array( range(env.START_YEAR,length=env.T) ),  # Years  年份索引（从真实年份开始）
    :Age => Array( range(env.START_AGE,length=env.S) ),  # Ages  年龄索引（从真实年龄开始）
    # -------------------- Full ages: Stocks 全年龄存量
    :𝒜 => zeros(env.T,env.S),  # Urban agents Wealth = a + Φ 城镇居民个人财富
    :a => zeros(env.T,env.S),  # Urban agents Asset (not wealth) 城镇居民个人资产
    :Φ => zeros(env.T,env.S),  # Urban agents Individual medical accounts  城镇居民个人医保账户
    # -------------------- Full ages: Flows 全年龄流量
    :c => zeros(env.T,env.S),  # total consumption 消费
    :m => zeros(env.T,env.S),  # total medical expenses 总医疗支出
    :MA => zeros(env.T,env.S),  # outpatient expenses 门诊医疗支出
    :MB => zeros(env.T,env.S),  # inpatient expenses 住院医疗支出
    :ΦGaps => zeros(env.T,env.S),  # the transfer payments from agent's own a_{s} to the individual account of UEBMI Φ_{s}
    # -------------------- Working ages 工作时期
    :w => zeros(env.T,env.Sr),  # profiled age-specific wage level 经过profiling的工资水平
    :Lab => zeros(env.T,env.Sr),  # labor supply 劳动
)

#


text:
"""
    readDataJRmodel(::Val{false}, filename; x0 = NaN)
read partial observations of the Jansen and Rit model according
to observation scheme given in https://arxiv.org/abs/1903.01138
"""
function readDataJRmodel(::Val{false}, filename; x0 = NaN)
    df = CSV.read(filename)
    obs = ℝ{1}.(df.x1)
    obsTime = Float64.(df.time)
    if isnan(NaN)
        x0 = ℝ{1}(df.x1[1])
    end
    fpt = [NaN for _ in obsTime[2:end]]
    fptOrPartObs = PartObs()
    df, x0, obs, obsTime, fpt, fptOrPartObs
end