text:
module Fmt

export @f_str

# types and utils
include("argument.jl")
include("spec.jl")
include("field.jl")
include("parser.jl")
include("compiler.jl")
include("interface.jl")
include("misc.jl")

# formatting functions
include("any.jl")
include("nothing.jl")
include("missing.jl")
include("char.jl")
include("string.jl")
include("int.jl")
include("bool.jl")
include("ptr.jl")
include("rational.jl")
include("float.jl")
include("bigfloat.jl")
include("irrational.jl")
include("complex.jl")

end


text:
using Test
include("../workflowtest.jl")

ext: jl


text:
"Map from PCRaster LDD value to a CartesianIndex"
const pcrdir = [
    CartesianIndex(1, -1),  # 1
    CartesianIndex(1, 0),  # 2
    CartesianIndex(1, 1),  # 3
    CartesianIndex(0, -1),  # 4
    CartesianIndex(0, 0),  # 5
    CartesianIndex(0, 1),  # 6
    CartesianIndex(-1, -1),  # 7
    CartesianIndex(-1, 0),  # 8
    CartesianIndex(-1, 1),  # 9
]

const mv = NaN

# timestep that the parameter units are defined in
const basetimestep = Second(Day(1))

# in case the input data is permuted the other way
permute_indices(inds) = [CartesianIndex(i[2], i[1]) for i in inds]

"Set at indices pit values (default = 5) in a gridded local drainage direction vector"
function set_pit_ldd(pits_2d, ldd, indices; pit = 5)
    pits = pits_2d[indices]
    index = filter(i -> isequal(pits[i], true), 1:length(indices))
    ldd[index] .= pit
    return ldd
end

"""
    active_indices(subcatch_2d, nodata)

Takes a 2D array of the subcatchments. And derive forward and reverse indices.

1: Get a list of `CartesianIndex{2}`` that are active, based on a nodata value.
These map from the 1D internal domain to the 2D external domain.

2: Make a reverse index, a `Matrix{Int}``, which maps from the 2D external domain to
the 1D internal domain, providing an Int which can be used as a linear index. Values of 0
represent inactive cells.
"""
function active_indices(subcatch_2d::AbstractMatrix, nodata)
    A = subcatch_2d
    all_inds = CartesianIndices(size(A))
    indices = filter(i -> !isequal(A[i], nodata), all_inds)

    reverse_indices = zeros(Int, size(A))
    for (i, I) in enumerate(indices)
        reverse_indices[I] = i
    end

    return indices, reverse_indices
end

function active_indices(network, key::Tuple)
    if :reservoir in key
        return network.reservoir.indices_outlet
    elseif :lake in key
        return network.lake.indices_outlet
    elseif :river in key
        return network.river.indices
    else
        return network.land.indices
    end
end

function lattometres(lat::Float64)
    m1 = 111132.92     # latitude calculation term 1
    m2 = -559.82       # latitude calculation term 2
    m3 = 1.175         # latitude calculation term 3
    m4 = -0.0023       # latitude calculation term 4
    p1 = 111412.84     # longitude calculation term 1
    p2 = -93.5         # longitude calculation term 2
    p3 = 0.118         # longitude calculation term 3

    # Calculate the length of a degree of latitude and longitude in meters
    latlen = m1 + (m2 * cosd(2.0 * lat)) + (m3 * cosd(4.0 * lat)) + (m4 * cosd(6.0 * lat))
    longlen = (p1 * cosd(lat)) + (p2 * cosd(3.0 * lat)) + (p3 * cosd(5.0 * lat))

    return longlen, latlen
end

"""
    set_states(instate_path, model, state_ncnames; <keyword arguments>)

Read states contained in `Dict` `state_ncnames` from NetCDF file located in `instate_path`,
and set states in `model` object. Active cells are selected with the corresponding network's
(`Vector{CartesianIndex}`) from the NetCDF file.

# Arguments
- `type = nothing`: type to convert data to after reading. By default no conversion is done.
"""
function set_states(
    instate_path,
    model,
    state_ncnames;
    type = nothing,
)
    @unpack network = model
    # states in NetCDF include dim time (one value) at index 3 or 4, 3 or 4 dims are allowed
    NCDataset(instate_path) do ds
        for (state, ncname) in state_ncnames
            sel = active_indices(network, state)
            dims = length(dimnames(ds[ncname]))
            # 4 dims, for example (x,y,layer,time) where dim layer is an SVector for soil layers
            if dims == 4
                A = permutedims(ds[ncname][sel, :, 1])
                # note that this array is allowed to have missings, since not every vertical
                # column is `maxlayers` layers deep
                A = replace!(A, missing => NaN)
                # Convert to desired type if needed
                if !isnothing(type)
                    if eltype(A) != type
                        A = map(type, A)
                    end
                end
                # set state in model object
                param(model, state) .= svectorscopy(A, Val{size(A)[1]}())
                # 3 dims (x,y,time)
            elseif dims == 3
                A = ds[ncname][sel, 1]
                A = nomissing(A)
                # Convert to desired type if needed
                if !isnothing(type)
                    if eltype(A) != type
                        A = map(type, A)
                    end
                end
                # set state in model object
                param(model, state) .= A
            else
                error("Number of state dims should be 3 or 4, number of dims = ", string(dims))
            end
        end
    end
end

"""
    ncread(nc, var; <keyword arguments>)

Read parameter `var` from NetCDF file `nc`. Supports various keyword arguments to get
selections of data in desired types, with or without missing values.

# Arguments
- `key=identity`: `key` is a function which maps `var` to the name of the parameter in the
        NetCDF file. The default `identity` keeps `var` as is.
- `sel=nothing`: a selection of indices, such as a `Vector{CartesianIndex}` of active cells,
        to return from the NetCDF. By default all cells are returned.
- `defaults=nothing`: a dictionary in which default values are looked up if `var` is not
        in `nc`. By default it gives an error in this case.
- `type=nothing`: type to convert data to after reading. By default no conversion is done.
- `allow_missing=false`: Missing values within `sel` is not allowed by default. Set to
        `true` to allow missing values.
- `fill=nothing`: Missing values are replaced by this fill value if `allow_missing` is `false`.
- `dimname` : name of third dimension of parameter `var`. By default no third dimension is expected.
"""
function ncread(
    nc,
    var;
    sel = nothing,
    defaults = nothing,
    type = nothing,
    allow_missing = false,
    fill = nothing,
    dimname = nothing,
)

    if isnothing(var)
        @assert !isnothing(defaults)
        if isnothing(dimname)
            return Base.fill(defaults, length(sel))
        else
            return Base.fill(defaults, (nc.dim[dimname], length(sel)))
        end
    end

    # Read the entire variable into memory, applying scale, offset and
    # set fill_values to missing.
    A = nc[var][:]

    # Take out only the active cells
    if !isnothing(sel)
        if isnothing(dimname)
            A = A[sel]
        else
            dim = findfirst(==(dimname), dimnames(nc[var]))
            if dim == 3
                A = permutedims(A[sel, :])
            elseif dim == 1
                A = A[:, sel]
            end
        end
    end

    if !allow_missing
        if isnothing(fill)
            # errors if missings are found
            A = nomissing(A)
        else
            # replaces missings with a fill value
            A = nomissing(A, fill)
        end
    end

    # Convert to desired type if needed
    if !isnothing(type)
        if eltype(A) != type
            A = convert(Array{type}, A)
        end
    end

    return A
end

"""
    set_layerthickness(d::Float64, sl::SVector)

Calculate actual soil thickness of layers based on a reference depth (e.g. soil depth or water table depth) `d`,
a SVector `sl` with cumulative soil depth starting at soil surface (0), and a SVector `tl` with actual thickness
per soil layer.
"""
function set_layerthickness(d::Float64, sl::SVector, tl::SVector)

    act_d = tl .* mv
    for i = 1:length(act_d)
        if d > sl[i+1]
            act_d = setindex(act_d, tl[i], i)
        elseif d - sl[i] > 0.0
            act_d = setindex(act_d, d - sl[i], i)
        end
    end

    nlayers = length(act_d) - sum(isnan.(act_d))
    return act_d, nlayers
end

"""
    detdrainlength(ldd, xl, yl)

Determines the drainaige length for a non square grid. Input `ldd` (drainage network), `xl` (length of cells in x direction),
`yl` (length of cells in y direction). Output is drainage length.
"""
function detdrainlength(ldd, xl, yl)
    # take into account non-square cells
    # if ldd is 8 or 2 use ylength
    # if ldd is 4 or 6 use xlength
    if ldd == 2 || ldd == 8
        yl
    elseif ldd == 4 || ldd == 6
        xl
    else
        hypot(xl, yl)
    end
end

"""
    detdrainwidth(ldd, xl, yl)

Determines the drainaige width for a non square grid. Input `ldd` (drainage network), `xl` (length of cells in x direction),
`yl` (length of cells in y direction). Output is drainage width.
"""
function detdrainwidth(ldd, xl, yl)
    # take into account non-square cells
    # if ldd is 8 or 2 use xlength
    # if ldd is 4 or 6 use ylength
    slantwidth = (xl + yl) * 0.5
    if ldd == 2 || ldd == 8
        xl
    elseif ldd == 4 || ldd == 6
        yl
    else
        slantwidth
    end
end

# 2.5x faster power method
"Faster method for exponentiation"
pow(x, y) = exp(y * log(x))

"Return the sum of the array `A` at indices `inds`"
function sum_at(A, inds)
    v = zero(eltype(A))
    for i in inds
        v += A[i]
    end
    return v
end

# https://juliaarrays.github.io/StaticArrays.jl/latest/pages/api/#Arrays-of-static-arrays-1
function svectorscopy(x::Matrix{T}, ::Val{N}) where {T,N}
    size(x, 1) == N || error("sizes mismatch")
    isbitstype(T) || error("use for bitstypes only")
    copy(reinterpret(SVector{N,T}, vec(x)))
end

"""
    fraction_runoff_toriver(graph, ldd, index_river, slope, n)

Determine ratio `frac` between `slope` river cell `index_river` and `slope` of each upstream neighbor (based on directed acyclic graph
`graph`).
"""
function fraction_runoff_toriver(graph, ldd, index_river, slope, n)
    frac = zeros(n)
    for i in index_river
        nbs = inneighbors(graph, i)
        for j in nbs
             if ldd[j] != ldd[i]
                 frac[j] = slope[j] / (slope[i] + slope[j])
             end
        end
    end
    return frac
end

"""
    equal_size_vectors(x)

Used in the structs of arrays to ensure all vectors are of equal length.

`equal_size_vectors(([1,2], [1,2,3]))` would throw an ArgumentError.
`equal_size_vectors(([4,5], [4,5]))` would pass.
`equal_size_vectors((1, [4,5], [4,5]))` would also pass, since `1` is not an AbstractVector.
"""
function equal_size_vectors(x)
    # all vectors in this struct should be the same size
    inds_vec = findall(arg -> isa(arg, AbstractVector), x)
    n = length(x[inds_vec[1]])
    x_vec = x[inds_vec]

    for arr in x_vec
        if length(arr) != n
            throw(ArgumentError("Not all vectors are of equal length"))
        end
    end
    return x
end

"""
    tosecond(x::Period)

Convert a Period into a Float64, which represents the number of seconds. Will fail if this
is not well defined, such as for Month.

# Examples
```julia-repl
julia> tosecond(Day(1))
86400.0
```
"""
tosecond(x::Hour) = Float64(Dates.value(Second(x)))
tosecond(x::Minute) = Float64(Dates.value(Second(x)))
tosecond(x::T) where {T<:DatePeriod} = Float64(Dates.value(Second(x)))
tosecond(x::T) where {T<:TimePeriod} = x / convert(T, Second(1))


text:
"""
An embedding in which the last point is guaranteed to lie within the convex
hull of the preceding points.
"""
mutable struct LinearlyInvariantEmbedding{D, T} <: AbstractEmbedding{D, T}
    points::AbstractArray{T, 2}
    embeddingdata::EmbeddingData{D, T}
end

export LinearlyInvariantEmbedding


text:
using Flux, Test, CUDA

@info "Testing GPU Support"
CUDA.allowscalar(false)

include("cuda.jl")
include("losses.jl")
include("layers.jl")

if CUDA.has_cudnn()
  @info "Testing Flux/CUDNN"
  include("cudnn.jl")
  include("curnn.jl")
else
  @warn "CUDNN unavailable, not testing GPU DNN support"
end


text:
using Dates
using FileIO
using MPI
using NCDatasets
using Printf
using Random
using StaticArrays
using Test

using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using Thermodynamics
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.Writers
using ClimateMachine.GenericCallbacks

using CLIMAParameters
using CLIMAParameters.Planet: grav, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

include("sin_init.jl")

function config_sin_test(FT, N, resolution, xmax, ymax, zmax)
    config = ClimateMachine.AtmosLESConfiguration(
        "Diagnostics SIN test",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_sin_test!,
    )

    return config
end

function config_diagnostics(driver_config)
    interval = "100steps"
    dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        interval,
        replace(driver_config.name, " " => "_"),
        writer = NetCDFWriter(),
    )
    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

function main()
    # Disable driver diagnostics as we're testing it here
    ClimateMachine.Settings.diagnostics = "never"

    FT = Float64

    # DG polynomial order
    N = 4

    # Domain resolution and size
    Δh = FT(50)
    Δv = FT(20)
    resolution = (Δh, Δh, Δv)

    xmax = FT(1500)
    ymax = FT(1500)
    zmax = FT(1500)

    t0 = FT(0)
    dt = FT(0.01)
    timeend = dt

    driver_config = config_sin_test(FT, N, resolution, xmax, ymax, zmax)

    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK54CarpenterKennedy,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        ode_dt = dt,
        init_on_cpu = true,
    )
    dgn_config = config_diagnostics(driver_config)

    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    Q = solver_config.Q
    solver = solver_config.solver

    outdir = mktempdir()
    currtime = ODESolvers.gettime(solver)
    starttime = replace(string(now()), ":" => ".")
    Diagnostics.init(mpicomm, param_set, dg, Q, starttime, outdir, false)
    GenericCallbacks.init!(
        dgn_config.groups[1],
        nothing,
        nothing,
        nothing,
        currtime,
    )

    ClimateMachine.invoke!(solver_config)

    # Check results
    mpirank = MPI.Comm_rank(mpicomm)
    if mpirank == 0
        dgngrp = dgn_config.groups[1]
        nm = @sprintf("%s_%s.nc", dgngrp.out_prefix, dgngrp.name)
        ds = NCDataset(joinpath(outdir, nm), "r")
        ds_u = ds["u"][:]
        ds_cov_w_u = ds["cov_w_u"][:]
        N = size(ds_u, 1)
        err = 0
        err1 = 0
        for i in 1:N
            u = ds_u[i]
            cov_w_u = ds_cov_w_u[i]
            err += (cov_w_u - 0.5)^2
            err1 += (u - 5)^2
        end
        close(ds)
        err = sqrt(err / N)
        @test err1 <= 1e-16
        @test err <= 2e-15
    end
end

main()


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

mutable struct IoK8sApiNetworkingV1IPBlock <: SwaggerModel
    cidr::Any # spec type: Union{ Nothing, String } # spec name: cidr
    except::Any # spec type: Union{ Nothing, Vector{String} } # spec name: except

    function IoK8sApiNetworkingV1IPBlock(;cidr=nothing, except=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1IPBlock, Symbol("cidr"), cidr)
        setfield!(o, Symbol("cidr"), cidr)
        validate_property(IoK8sApiNetworkingV1IPBlock, Symbol("except"), except)
        setfield!(o, Symbol("except"), except)
        o
    end
end # type IoK8sApiNetworkingV1IPBlock

const _property_map_IoK8sApiNetworkingV1IPBlock = Dict{Symbol,Symbol}(Symbol("cidr")=>Symbol("cidr"), Symbol("except")=>Symbol("except"))
const _property_types_IoK8sApiNetworkingV1IPBlock = Dict{Symbol,String}(Symbol("cidr")=>"String", Symbol("except")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1IPBlock }) = collect(keys(_property_map_IoK8sApiNetworkingV1IPBlock))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1IPBlock }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiNetworkingV1IPBlock[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1IPBlock }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1IPBlock[property_name]

function check_required(o::IoK8sApiNetworkingV1IPBlock)
    (getproperty(o, Symbol("cidr")) === nothing) && (return false)
    true
end

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


text:
#==============================================================================#
# AWSSDK/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2017 - All rights reserved
#==============================================================================#

using Base.Test
using AWSSDK

println("There are no tests defined AWSSDK.jl. See tests in AWSCore.jl")

#==============================================================================#
# End of file.
#==============================================================================#


text:
using GLFW
using ModernGL
using CSyntax
using Printf

# set up OpenGL context version
# it seems OpenGL 4.1 is the highest version supported by MacOS.
@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

@static if Sys.isapple()
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
else
    GLFW.DefaultWindowHints()
end

# _update_fps_counter
let previous = time()
    frame_count = 0
    global function updatefps(window::GLFW.Window)
        current = time()
        elapsed = current - previous
        if elapsed > 0.25
            previous = current
            fps = frame_count / elapsed
            GLFW.SetWindowTitle(window, @sprintf("opengl @ fps: %.2f", fps))
            frame_count = 0
        end
        frame_count += 1
    end
end

# set up GLFW error callbacks
@info "starting GLFW ..."
@info GLFW.GetVersionString()

# error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# set up GLFW key callbacks : press Esc to escape
function key_callback(window::GLFW.Window, key::GLFW.Key, scancode::Cint, action::GLFW.Action, mods::Cint)
    # GL_TRUE was original code, but julia true works just fine
    key == GLFW.KEY_ESCAPE && action == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# tell GLFW to run this function whenever the framebuffer size is changed
function framebuffer_size_callback(window::GLFW.Window, w::Cint, h::Cint)
	global fb_width = w
	global fb_height = h
	println("framebuffer size: ($w, $h)")
    return nothing
end

# tell GLFW to run this function whenever the window size is changed
function window_size_callback(window::GLFW.Window, w::Cint, h::Cint)
	global width = w
	global height = h
	println("window size: ($w, $h)")
    return nothing
end

width, height = fb_width, fb_height = 640, 480
window = GLFW.CreateWindow(width, height, "Extended Init.")
@assert window != C_NULL "could not open window with GLFW3."

GLFW.SetKeyCallback(window, key_callback)
GLFW.SetWindowSizeCallback(window, window_size_callback)
GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)
GLFW.MakeContextCurrent(window)
GLFW.WindowHint(GLFW.SAMPLES, 4)

# enable depth test
glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# get version info
renderer = unsafe_string(glGetString(GL_RENDERER))
version = unsafe_string(glGetString(GL_VERSION))
@info "Renderder: $renderer"
@info "OpenGL version supported: $version"

# hard-coded shaders
const vert_source = """
#version 410 core
in vec3 vp;
void main(void)
{
    gl_Position = vec4(vp, 1.0);
}"""
const frag_source = """
#version 410 core
out vec4 frag_colour;
void main(void)
{
    frag_colour = vec4(0.5, 0.0, 0.5, 1.0);
}"""

# compile shaders
vert_shader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vert_shader, 1, Ptr{GLchar}[pointer(vert_source)], C_NULL)
glCompileShader(vert_shader)
frag_shader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(frag_shader, 1, Ptr{GLchar}[pointer(frag_source)], C_NULL)
glCompileShader(frag_shader)

# create and link shader program
shader_prog = glCreateProgram()
glAttachShader(shader_prog, vert_shader)
glAttachShader(shader_prog, frag_shader)
glLinkProgram(shader_prog)

# vertex data
points = GLfloat[ 0.0,  0.5, 0.0,
                  0.5, -0.5, 0.0,
                 -0.5, -0.5, 0.0]

# create buffers located in the memory of graphic card
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)

# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, width, height)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
abstract type AbstractLayer end

abstract type ParamLayer <: AbstractLayer end


text:
# Testing version which uses Transit_Struct:
using PyPlot

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

npts = 1000
iu = 30
u_n = ones(iu) / iu
r = 0.1
b0 = zeros(npts)
for i=1:npts
    b0[i] = 3.0 * (i / float(npts) - 0.5)
end

lc_ana = zeros(npts); lc_str=zeros(npts)
trans = transit_init(r,0.0,u_n,false)
for i=1:length(b0)
  b=b0[i]
  lc_ana[i] = transit_poly(r,abs(b),u_n)
  trans.b = abs(b)
  lc_str[i] = transit_poly!(trans)
  if mod(i,10) == 0
    println("i: ",i,"b: ",b," lc_ana: ",lc_ana[i]," lc_str: ",lc_str[i])
  end
end

using PyPlot

plot(b0,lc_ana, linewidth=2, label="Normal")
plot(b0,lc_str, linewidth=1, label="Struct")

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# this file create julia functions from the python modules

empty!(type_map)  # for behaving nicely in system image

for (ee_expr, julia_type) in pre_type_map
    type_map[ee_expr()] = julia_type
end

# first loop over the ee types and set functions that take type
for key in collect(keys(type_map))
    # get list of methods
    methods = py"dir"(key)
    # loop over all methods and create julia methods
    for method in methods
        if ~startswith(method, "_")
            pyattr_set([type_map[key]], Symbol(method))
            if ~(method in ee_exports)
                push!(ee_exports, method)
            end
        end
    end
end

# next do two passes on the individual modules
# first pass is to get the methods for each module
# second pass is to map the methods to julia funcs

modules = [
    :Array,
    :Classifier,
    :Clusterer,
    :ConfusionMatrix,
    :Date,
    :DateRange,
    :Dictionary,
    :ErrorMargin,
    :Feature,
    :FeatureCollection,
    :Filter,
    :Geometry,
    :Image,
    :ImageCollection,
    :Join,
    :Kernel,
    :List,
    :Model,
    :Number,
    :PixelType,
    :Projection,
    :Reducer,
    :String,
    :Terrain,
    :data,
]

mod_methods = OrderedDict()

# loop over the modules
for mod in modules
    @eval begin
        submethods = collect(py"dir"(ee.$(string(mod))))
        # get the module functions and add to dict
        mod_methods[$(string(mod))] = submethods
    end
end

# loop over the modules again
for (k, v) in mod_methods
    # loop over the module methods
    for submethod in v
        # check if it is not private
        if ~startswith(submethod, "_")
            # create a julia function of the public methods
            m = Symbol(submethod)
            @eval begin
                function $m(args...; kwargs...)
                    method = ee.$k.$(string(submethod))
                    result = pycall(method, PyObject, args...; kwargs...)
                    ee_wrap(result)
                end
            end
            # add to list of methods for module export
            if ~(submethod in ee_exports)
                push!(ee_exports, submethod)
            end
        end
    end
end

# for some reason the ee.Algorithms module is a Dict....why???
# function to recusively search though the ee.Algorithms module and wrap all python functions
function wrap_eealgorithms(dict)
    for (k, v) in dict
        if typeof(v) == PyCall.PyObject
            m = Symbol(k)
            @eval begin
                function $m(args...; kwargs...)
                    method = $v
                    result = pycall(method, PyObject, args...; kwargs...)
                    ee_wrap(result)
                end
            end
            # add to list of methods for module export
            if ~(k in ee_exports)
                push!(ee_exports, k)
            end
        else
            wrap_eealgorithms(v)
        end
    end
end

# get the dict of ee algorithms
algorithms_dict = ee.Algorithms
# remove the GeometryConstructors functions that clobber with others
delete!(algorithms_dict, "GeometryConstructors")

# apply recusive wrapping
wrap_eealgorithms(algorithms_dict)


text:
@test max_smooth(-0.1, 100, 0.0001) ≈ 100
@test max_smooth(-0.1, -100, 0.0001) ≈ -0.1
@test max_smooth(-0.1, 100, 0) == max(-0.1,100)

@test min_smooth(-0.1, 100, 0.0001) ≈ -0.1
@test min_smooth(-0.1, -100, 0.0001) ≈ -100


text:
#### ExtremumFuncs

export Extremum, Min, Max
export n_hat, binary, ghost_vec, ghost_dual

abstract type Extremum end

"""
    Min

A type for dispatching on the minimum of
the domain along a particular direction
"""
struct Min <: Extremum end

"""
    Max

A type for dispatching on the maximum of
the domain along a particular direction
"""
struct Max <: Extremum end

"""
    n_hat(::Extremum)

The outward normal vector to the boundary
"""
n_hat(::Min) = -1
n_hat(::Max) = 1

"""
    binary(::Extremum)

Returns 0 for `Min` and 1 for `Max`
"""
binary(::Min) = 0
binary(::Max) = 1

"""
    ghost_vec(::Extremum)

A 3-element vector near the boundary with 1's
on ghost cells and 0's on interior cells
"""
ghost_vec(::Min) = (1, 0, 0)
ghost_vec(::Max) = (0, 0, 1)

"""
    ghost_dual(::Extremum)

A 2-element vector near the boundary with 1's
on ghost cells and 0's on interior cells
"""
ghost_dual(::Min) = (1, 0)
ghost_dual(::Max) = (0, 1)


text:
function verify_gcc(gcc)
    try
        return success(`$gcc --version`)
    catch
        return false
    end
end

# TODO: remove once Julia v0.7 is released
const julia_v07 = VERSION > v"0.7-"
if julia_v07
    const isunix = Sys.isunix
    const iswindows = Sys.iswindows
else
    const isunix = is_unix
    const iswindows = is_windows
end

if iswindows()
    using WinRPM
end

function build()
    gccpath = ""
    if isfile("deps.jl")
        include("deps.jl")
        if verify_gcc(gcc)
            info("gcc already installed and package already build")
            return
        else
            rm("deps.jl")
        end
    end

    if haskey(ENV, "CC")
        if !verify_gcc(`$(ENV["CC"]) -v`)
            error("Using compiler override from environment variable CC = $(ENV["CC"]), but unable to run `$(ENV["CC"]) -v`.")
        end
        gccpath = ENV["CC"]
        info("using $gccpath as a compiler from environment variable CC")
    end

    info("installing gcc")

    if verify_gcc("cc")
        gccpath = "cc"
        info("using cc as a compiler")
    elseif iswindows()
        gccpath = joinpath(
            WinRPM.installdir, "usr", "$(Sys.ARCH)-w64-mingw32",
            "sys-root", "mingw", "bin", "gcc.exe"
        )
        if !isfile(gccpath)
            WinRPM.install("gcc", yes = true)
        end
        if !isfile(gccpath)
            error("Couldn't install gcc via winrpm")
        end
        info("using gcc from WinRPM as a compiler")
    elseif isunix() && verify_gcc("gcc")
        gccpath = "gcc"
        info("using gcc as a compiler")
    end

    if isempty(gccpath)
        error("Please make sure to provide a working gcc in your path!")
    end
    open("deps.jl", "w") do io
        print(io, "const gcc = ")
        println(io, '"', escape_string(gccpath), '"')
    end
end

build()


text:
# Helper types to represent various combinations of vectors/matrices/scalings
const AbstractNumOrVector  = Union{Number, AbstractVector}
const AbstractNumOrMatrix  = Union{Number, AbstractVector, AbstractMatrix}
const AbstractNumOrUniform = Union{Number, AbstractVector, AbstractMatrix, UniformScaling}

const AbstractNumOrMatrixOrSymbol  = Union{Number, AbstractVector, AbstractMatrix, Symbol}
const AbstractNumOrUniformOrSymbol = Union{Number, AbstractVector, AbstractMatrix, UniformScaling, Symbol}

# Helper type to represent a Hermitian view of a block matrix
const HermitianBlockMatrix = Hermitian{T, M} where {T <: Number, M <: BlockMatrix{<:T}}

# Helper functions to convert various types to a 2d matrix
to_matrix(T, A::AbstractVector) = Matrix{T}(reshape(A, length(A), 1))
to_matrix(T, A::AbstractMatrix) = T.(A)  # Fallback
to_matrix(T, A::Number) = fill(T(A), 1, 1)
to_matrix(T, A::UniformScaling, s::Integer) = Matrix{T}(A*I, s, s)

function dftmatrix( n::Integer; unitary::Bool = false )

    M = fft( 1.0*I(n), 1 )

    # A scaling of 1/(√n) needs to be applied to make the DFT matrix unitary
    if( unitary )
        M = 1 / sqrt( n ) * M
    end

    return M
end


text:
const API_ROOT = "/api/v1"
const CLUSTERING_ROOT = "$(homedir())/datasets/participation/clusterings"
const DATASET_ROOT = "$(homedir())/datasets/participation/databases"
const DISTANCE_ROOT = "$(homedir())/datasets/participation/distances"
const HEADERS = [
    "Access-Control-Allow-Headers" => "*",
    "Access-Control-Allow-Methods" => "GET; POST; OPTIONS",
    "Access-Control-Allow-Origin" => "*",
    "Content-Type" => "application/json"
]

"""
    getdatasets(req::HTTP.Request)
"""
function getdatasets(::HTTP.Request)
    names = readdir(DATASET_ROOT)
    filter!(endswith(".sqlite"), names)
    return names
end

"""
    getcontributions(req::HTTP.Request)
"""
function getcontributions(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    database = SQLite.DB("$DATASET_ROOT/$datasetname")
    return DBInterface.execute(database, "SELECT * FROM contribution ORDER BY id") |> DataFrame
end

"""
    getcontribution(req::HTTP.Request)
"""
function getcontribution(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    contributionid = HTTP.URIs.splitpath(req.target)[6]
    database = SQLite.DB("$DATASET_ROOT/$datasetname")
    return DBInterface.execute(database, "SELECT * FROM contribution WHERE id = $contributionid") |> DataFrame
end

"""
    getclusterings(req::HTTP.Request)
"""
function getclusterings(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    names = readdir("$CLUSTERING_ROOT/$datasetname")
    filter!(endswith(".json"), names)
    return names
end

"""
    getclustering(req::HTTP.Request)
"""
function getclustering(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    clustering = HTTP.URIs.splitpath(req.target)[6]
    return read("$CLUSTERING_ROOT/$datasetname/$clustering", String)
end

"""
    getdistances(req::HTTP.Request)
"""
function getdistances(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    names = readdir("$DISTANCE_ROOT/$datasetname")
    filter!(endswith(".json"), names)
    return names
end

"""
    getdistance(req::HTTP.Request)
"""
function getdistance(req::HTTP.Request)
    datasetname = HTTP.URIs.splitpath(req.target)[4]
    distance = HTTP.URIs.splitpath(req.target)[6]
    return read("$DISTANCE_ROOT/$datasetname/$distance", String)
end

router = HTTP.Router()
HTTP.@register(router, "GET", "$API_ROOT/datasets", getdatasets)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/clusterings", getclusterings)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/clustering/*", getclustering)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/contributions", getcontributions)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/contribution/*", getcontribution)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/distances", getdistances)
HTTP.@register(router, "GET", "$API_ROOT/datasets/*/distance/*", getdistance)

"""
    handlejson(req::HTTP.Request)
"""
function handlejson(req::HTTP.Request)
    res = HTTP.handle(router, req)
    code = 200
    if isa(res, HTTP.Response) && HTTP.iserror(res)
        code = convert(Int, HTTP.status(res))
        body = JSON3.write(HTTP.body(res))
    elseif Tables.istable(res)
        body = arraytable(res)
    elseif isa(res, String)
        body = res
    else
        body = JSON3.write(res)
    end
    return HTTP.Response(code, HEADERS; body=body)
end

"""
    handlecors(req::HTTP.Request)
"""
function handlecors(req::HTTP.Request)
    if HTTP.hasheader(req, "OPTIONS")
        return HTTP.Response(200, HEADERS)
    end
    return handlejson(req)
end

"""
    startrest(port::Integer)
"""
function startrest(port::Integer=3010)
    @async HTTP.serve(handlecors, HTTP.Sockets.localhost, port)
    println("REST server is running at http://$(HTTP.Sockets.localhost):$port")
end


text:
module TestExamples2D

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "2D" begin

# Run basic tests
@testset "Examples 2D" begin
  # Linear advection
  include("test_examples_2d_advection.jl")

  # Hyperbolic diffusion
  include("test_examples_2d_hypdiff.jl")

  # Compressible Euler
  include("test_examples_2d_euler.jl")

  # Compressible Euler Multicomponent
  include("test_examples_2d_eulermulti.jl")

  # MHD
  include("test_examples_2d_mhd.jl")

  # Lattice-Boltzmann
  include("test_examples_2d_lbm.jl")
end

# Coverage test for all initial conditions
@testset "Tests for initial conditions" begin
  # Linear scalar advection
  @testset "elixir_advection_extended.jl with initial_condition_sin_sin" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [0.0001424424872539405],
      linf = [0.0007260692243253875],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_sin_sin)
  end

  @testset "elixir_advection_extended.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [3.2933000250376106e-16],
      linf = [6.661338147750939e-16],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @testset "elixir_advection_extended.jl with initial_condition_linear_x_y" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [2.478798286796091e-16],
      linf = [7.105427357601002e-15],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_linear_x_y,
      boundary_conditions = Trixi.boundary_condition_linear_x_y,
      periodicity=false)
  end

  @testset "elixir_advection_extended.jl with initial_condition_linear_x" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [1.475643203742897e-16],
      linf = [1.5543122344752192e-15],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_linear_x,
      boundary_conditions = Trixi.boundary_condition_linear_x,
      periodicity=false)
  end

  @testset "elixir_advection_extended.jl with initial_condition_linear_y" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
      l2   = [1.5465148503676022e-16],
      linf = [3.6637359812630166e-15],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_linear_y,
      boundary_conditions = Trixi.boundary_condition_linear_y,
      periodicity=false)
  end

  # Compressible Euler
  @testset "elixir_euler_vortex.jl one step with initial_condition_density_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.003489659044164644, 0.0034896590441646494, 0.0034896590441646502, 0.003489659044164646],
      linf = [0.04761180654650543, 0.04761180654650565, 0.047611806546505875, 0.04761180654650454],
      maxiters = 1,
      initial_condition = initial_condition_density_pulse)
  end

  @testset "elixir_euler_vortex.jl one step with initial_condition_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.00021747693728234874, 0.0022010142997830533, 0.0022010142997830485, 0.010855273768135729],
      linf = [0.005451116856088789, 0.03126448432601536, 0.03126448432601536, 0.14844305553724624],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_pressure_pulse)
  end

  @testset "elixir_euler_vortex.jl one step with initial_condition_density_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.003473649182284682, 0.005490887132955628, 0.005490887132955635, 0.015625074774949926],
      linf = [0.046582178207169145, 0.07332265196082899, 0.07332265196082921, 0.2107979471941368],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_density_pressure_pulse)
  end

  @testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [7.89034964747135e-17, 8.095575651413758e-17, 1.0847287658433571e-16, 1.2897732640029767e-15],
      linf = [2.220446049250313e-16, 3.191891195797325e-16, 4.163336342344337e-16, 3.552713678800501e-15],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @testset "elixir_euler_sedov_blast_wave.jl one step" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0021037031798961914, 0.010667428589443025, 0.01066742858944302, 0.10738893384136498],
      linf = [0.11854059147646778, 0.7407961272348982, 0.7407961272348981, 3.92623931433345],
      maxiters=1)
  end

  @testset "elixir_euler_sedov_blast_wave.jl one step with initial_condition_medium_sedov_blast_wave" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0021025532272874827, 0.010661548568022292, 0.010661548568022284, 0.10734939168392313],
      linf = [0.11848345578926645, 0.7404217490990809, 0.7404217490990809, 3.9247328712525973],
      maxiters=1, initial_condition=initial_condition_medium_sedov_blast_wave)
  end

  # GLM-MHD
  @testset "elixir_mhd_alfven_wave.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [7.144325530681224e-17, 2.123397983547417e-16, 5.061138912500049e-16, 3.6588423152083e-17, 8.449816179702522e-15, 3.9171737639099993e-16, 2.445565690318772e-16, 3.6588423152083e-17, 9.971153407737885e-17],
      linf = [2.220446049250313e-16, 8.465450562766819e-16, 1.8318679906315083e-15, 1.1102230246251565e-16, 1.4210854715202004e-14, 8.881784197001252e-16, 4.440892098500626e-16, 1.1102230246251565e-16, 4.779017148551244e-16],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @testset "elixir_mhd_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
      l2   = [1.2429643014570362, 1.7996363685163963, 1.6899889382208215, 0.0, 2.263014495582255, 0.2127948919559293, 0.23341167012961367, 0.0, 0.003341061230142962],
      linf = [10.401662016997985, 14.061042946011094, 15.556382737749237, 0.0, 16.714999099689578, 1.3250449624793987, 1.4148467737020096, 0.0, 0.0878998114279951],
      tspan = (0.0, 0.05))
  end

  @testset "elixir_mhd_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"),
      l2   = [0.17569823349539196, 3.853292514951796, 2.474036084054808, 0.0, 355.36545703316915, 2.3525087027248084, 1.394705056983077, 0.0, 0.029910308624236454],
      linf = [1.5831869462980066, 44.20237543674303, 12.866364462538552, 0.0, 2237.8614138686, 13.066894956593798, 8.984965484247244, 0.0, 0.5226498664960756],
      tspan = (0.0, 0.003))
  end

  # LBM
  @testset "elixir_lbm_couette.jl with initial_condition_couette_steady" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_lbm_couette.jl"),
      l2   = [9.321369073400123e-16, 1.6498793963435488e-6, 5.211495843124065e-16,
              1.6520893954826173e-6, 1.0406056181388841e-5, 8.801606429417205e-6,
              8.801710065560555e-6, 1.040614383799995e-5, 2.6135657178357052e-15],
      linf = [1.4432899320127035e-15, 2.1821189867266e-6, 8.881784197001252e-16,
              2.2481261510165496e-6, 1.0692966335143494e-5, 9.606391697600247e-6,
              9.62138334279633e-6, 1.0725969916147021e-5, 3.3861802251067274e-15],
      initial_condition=initial_condition_couette_steady,
      tspan = (0.0, 1.0))
  end

  @testset "elixir_lbm_lid_driven_cavity.jl with stationary walls" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_lbm_lid_driven_cavity.jl"),
      l2   = [1.7198203373689985e-16, 1.685644347036533e-16, 2.1604974801394525e-16,
              2.1527076266915764e-16, 4.2170298143732604e-17, 5.160156233016299e-17,
              6.167794865198169e-17, 5.24166554417795e-17, 6.694740573885739e-16],
      linf = [5.967448757360216e-16, 6.522560269672795e-16, 6.522560269672795e-16,
              6.245004513516506e-16, 2.1163626406917047e-16, 2.185751579730777e-16,
              2.185751579730777e-16, 2.393918396847994e-16, 1.887379141862766e-15],
      boundary_conditions=boundary_condition_wall_noslip,
      tspan = (0, 0.1))
  end
end

@testset "Displaying components 2D" begin
  @test_nowarn include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"))

  # test both short and long printing formats
  @test_nowarn show(mesh); println()
  @test_nowarn println(mesh)
  @test_nowarn display(mesh)

  @test_nowarn show(equations); println()
  @test_nowarn println(equations)
  @test_nowarn display(equations)

  @test_nowarn show(solver); println()
  @test_nowarn println(solver)
  @test_nowarn display(solver)

  @test_nowarn show(solver.basis); println()
  @test_nowarn println(solver.basis)
  @test_nowarn display(solver.basis)

  @test_nowarn show(solver.mortar); println()
  @test_nowarn println(solver.mortar)
  @test_nowarn display(solver.mortar)

  @test_nowarn show(semi); println()
  @test_nowarn println(semi)
  @test_nowarn display(semi)

  @test_nowarn show(summary_callback); println()
  @test_nowarn println(summary_callback)
  @test_nowarn display(summary_callback)

  @test_nowarn show(amr_controller); println()
  @test_nowarn println(amr_controller)
  @test_nowarn display(amr_controller)

  @test_nowarn show(amr_callback); println()
  @test_nowarn println(amr_callback)
  @test_nowarn display(amr_callback)

  @test_nowarn show(stepsize_callback); println()
  @test_nowarn println(stepsize_callback)
  @test_nowarn display(stepsize_callback)

  @test_nowarn show(save_solution); println()
  @test_nowarn println(save_solution)
  @test_nowarn display(save_solution)

  @test_nowarn show(analysis_callback); println()
  @test_nowarn println(analysis_callback)
  @test_nowarn display(analysis_callback)

  @test_nowarn show(alive_callback); println()
  @test_nowarn println(alive_callback)
  @test_nowarn display(alive_callback)

  @test_nowarn println(callbacks)
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # 2D

end #module


text:
@testset "nl_write" begin
    # Turn on debug mode so files persist
    old_debug = JAMSDWriter.CONFIG[:debug]
    JAMSDWriter.setdebug(true)

#    filename = "test"
#    filepath = joinpath(JAMSDWriter.solverdata_dir, filename)
#    JAMSDWriter.clean_solverdata()

#    context("all temp files deleted successfully") do
#        @fact length(readdir(JAMSDWriter.solverdata_dir)) --> 1
#    end

    m = Model(solver=JAMSDWriter.JAMSDSolver())
    @variable(m, x >= 0)
    @objective(m, Min, x)
    solve(m)

    # Reset debug mode and clean up
    JAMSDWriter.setdebug(old_debug)
    JAMSDWriter.clean_solverdata()
end


text:
#ДАНО: Робот - в произвольной клетке поля с внутренними перегородками(могут быть прямоугольники)
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы
include("library.jl")
function moveAroundAndPut!(r::Robot)
    for side in instances(HorizonSide)
        putmarker!(r)
        marksLine!(r, side)
    end
end
function perimetr_with_obstacles(r::Robot)
    sides = moveAndReturnDirections!(r)
    moveAroundAndPut!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
## iop_online
#  Shi Yu, Chaosheng Dong
#  shi_yu@vanguard.com
#  chaosheng@pitt.edu
#  All rights reserved 2017-03-17
# function
@everywhere using JuMP, Gurobi, Distributions, Clustering, Distances, CSV, NPZ
@everywhere function f(Q,c,r,x)
    1/2*r*vecdot(x, Q*x) - vecdot(c,x)
end

@everywhere function GenerateData(n_f,Sigma_f,c,Af,b,T,ra)
   x = SharedArray{Float64,2}((n_f,T));
   @parallel (+) for i = 1:T
       genData = Model( solver=GurobiSolver(OutputFlag=0, MIPGap = 1e-6,TimeLimit = 600 ) )
       @variable(genData, wf[1:n_f])
       @objective(genData, Min, f(Sigma_f,c,ra,wf))
       @constraint(genData, Af*wf .>= b )
       solve(genData)
       x[:,i] = getvalue(wf)
   end
   return x
end

@everywhere function GenerateData2(n,Q,c,A,b,r)
    genData = Model( solver=GurobiSolver( OutputFlag=0 ) )
    @variable(genData, y[1:n] )
    @objective(genData, Min, f(Q,c,r,y) )
    @constraint(genData, A*y .>= b)
    solve(genData)
    return getvalue(y)
end

@everywhere function qp_r_online_ver2(m,nf,M,eta,Sigma_f,c,Af,b,wf,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 600) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
   @variable(problem, x[1:nf])
#    @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2 + eta*vecdot(wf - x,wf - x) )
   @constraint(problem, Af*x .>= b )
   @constraint(problem, u .<= M*z )
   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, Sigma_f*x - ra_t*c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t), getvalue(x)
end

@everywhere function qp_r_online_ver3(m,nf,M,eta,Sigma_f,c,Af,b,wf,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 3000) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
   @variable(problem, x[1:nf])
#  @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   #@variable(problem, e<=0.0001)
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2 + eta*vecdot(wf - x,wf - x))
   @constraint(problem, Af*x .>= b)
   @constraint(problem, u .<= M*z )
   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, Sigma_f*x - ra_t*c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t), getvalue(x)
end

@everywhere function qp_r_online_ver4(m,nf,M,eta,Sigma_f,c,Af,b,x,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 300) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
#    @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2)
   @constraint(problem, Af*x .>= b )
   @constraint(problem, u .<= M*z )
#   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, ra_t*Sigma_f*x - c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t)
end

@everywhere function qp_r_online_ver5(m,lambda,Q,c,A,b,x,r)
    problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-3, TimeLimit = 300, InfUnbdInfo=1) )
    @variable(problem, -10 <= r_t <= 10)
    @variable(problem, u[1:m] >= 0)
    @variable(problem, t[1:m],Bin)
    @variable(problem, eta1)
    @objective(problem, Min, eta1/2 )
    @constraint(problem, norm([2*(r - r_t); eta1-1]) <= eta1 + 1 )
    #@constraint(problem, A*x .>= b )
    @constraint(problem, u.<= M*t )
    #@constraint(problem, A*x - b .<= M*(1-t) )
    @constraint(problem, r_t*Q*x - c - A'*u .== zeros(m) )
    solve(problem)
    return getvalue(r_t)
end

@everywhere function qp_r_online_ver2(m,nf,M,eta,Sigma_f,c,Af,b,wf,ra)
   problem = Model(solver=GurobiSolver(OutputFlag = 0, MIPGap = 1e-6, TimeLimit = 300) )
   @variable(problem, 0.000001<= ra_t<= 100) # risk aversion value at time t
   @variable(problem, x[1:nf])
#    @variable(problem, 0<=x[1:n]<=1 )
   @variable(problem, u[1:m] >= 0 )
   @variable(problem, z[1:m], Bin )
   @objective(problem, Min, vecdot(ra - ra_t,ra - ra_t)/2 + eta*vecdot(wf - x,wf - x) )
   @constraint(problem, Af*x .>= b )
   @constraint(problem, u .<= M*z )
   @constraint(problem, Af*x - b .<= M*(1-z))
   @constraint(problem, Sigma_f*x - ra_t*c - Af'*u .== zeros(nf) )
   solve(problem)
   return getvalue(ra_t), getvalue(x)
end

@everywhere function decompose(A_return, num_comp)
    # spectral decompositon of the covariance matrix
    cov_mat = cor(A_return);
    Q = cov_mat;
    Q[isnan(Q)] = 0;
    Q = Q + 0.0001*I;

    factor = eigfact(Q);
    eigvalues = factor[:values];
    eigvectors = factor[:vectors];

    large_eigs = eigvalues[length(eigvalues)-num_comp+1:length(eigvalues)];

    F = eigvectors[:,length(eigvalues)-num_comp+1:length(eigvalues)];

    Sigma_f = Diagonal(large_eigs);

   return F, Sigma_f
end

@everywhere function compute_c(A_return)
#    A_return: matrix of size (num of time steps, num of assets)

    num_asset = size(A_return)[2];
    c = SharedArray{Float64,1}(num_asset);
    for ind_asset = 1:num_asset
        A_vec = A_return[:, ind_asset];
        nonzeroInd = find(x->x!=0, A_vec);
        non_zeroA_vec = A_vec[nonzeroInd];
        c[ind_asset] = mean(non_zeroA_vec);
    end
   return c
end

@everywhere function compute_c_max(A_return)
#    A_return: matrix of size (num of time steps, num of assets)

    num_asset = size(A_return)[2];
    #print("num_asset ",num_asset)
    c = SharedArray{Float64,1}(num_asset);
    for ind_asset = 1:num_asset
        A_vec = A_return[:, ind_asset];
        nonzeroInd = find(x->x>=0, A_vec);
        non_zeroA_vec = A_vec[nonzeroInd];
        c[ind_asset] = median(non_zeroA_vec);
    end
   return c
end

using CSV, DataFrames, DataFramesMeta, Base.Dates
@everywhere using JuMP, Gurobi, Distributions, Clustering, Distances, Plots, JLD

fund_name="arlio";

A = zeros(10, 425);

icount = 1;

lambda =100;
M=100;

X_obs = CSV.read("/home/syu/Documents/Projects/ipo_kdd/ipo/data/drl/drl_100_assets_target1.1/arlio_X.txt",delim=",",datarow=1);
X_obs = Matrix(X_obs);
X_obs = X_obs[:,1:end]
Xobsmat = Matrix(X_obs);
Xobsmat = transpose(Xobsmat);

for sind=1:10:100

    universe_results = Vector{Vector{Float64}}();
    universe_results_ratio = Vector{Vector{Float64}}();

    all_results = Vector{Vector{Float64}}();
    all_results2 = Vector{Vector{Float64}}();
    Time = Float64[];

    for ind=1:1:425
        A_return = CSV.read("/home/syu/Documents/Projects/ipo_kdd/ipo/data/drl/drl_100_assets_target1.1/arlio_A_"*string(ind)*".txt",datarow=1);
        A_return = Matrix(A_return);
        println("A matrix size ", size(A_return));
        #sind=size(A_return)[1]-lookback;

        A_return = A_return[sind:size(A_return)[1], :];
        #A_return = A_return/12;

        #A_return = remove_extreme_negative_return(A_return);
        #A_return = A_return./100;
        #c = maximum(A_return,2);

        ind_results = Vector{Float64}();

        ra =1;
        #c = compute_c(A_return);
        c = compute_c_max(A_return);
        #println(c);

        c[isnan.(c)] .= 0;
        #print(c);

        n=10;

        Q = cov(A_return);
        Q = Q + 0.001*I;

        xx_obs = X_obs[:, 1:ind];
        est_return = Vector{Float64}();
        est_return2 = Vector{Float64}();

        #println(size(xx_obs)[2]);

        tic();

        if size(xx_obs)[2]<=10
            s=1
        else
            s=size(xx_obs)[2]-10
        end

        for t = s:size(xx_obs)[2]

            #println(ind,ra,t);

            y = xx_obs[:,t];
            y = y/sum(y);

            y[isnan.(y)] .= 0;

            A_complex=-[eye(n);ones(1,n);-ones(1,n);-eye(n)]

            b_complex= -[ones(n,1);1;-1;zeros(n,1)];

            (m,n) = size(A_complex);
            eta = lambda*t^(-1/2);
            #print(ra, y,c)

            try
                ra, x = qp_r_online_ver2(m,n,M,eta,Q,c,A_complex,b_complex,y,ra);
            catch e
                println(e);
                push!(est_return, 0);
                push!(est_return2, transpose(y)*Q*y./(transpose(c)*y));
                continue
            end

            push!(est_return, ra);
            push!(est_return2, transpose(y)*Q*y./(transpose(c)*y));
        end
        print(est_return);
        push!(all_results, est_return);
        push!(all_results2, est_return2);
        t = toc();
        Time = push!(Time, t);
    end

    temp_result = Vector{Float64}();
    temp_result2 = Vector{Float64}();

    for i=1:425
        push!(temp_result, all_results[i][end]);
        push!(temp_result2, all_results2[i][end]);
    end
    push!(universe_results, temp_result);
    push!(universe_results_ratio, temp_result2);

    for i=1:425
        A[icount, i]=universe_results[1][i];
    end

    icount = icount+1;
end

save("/home/syu/Documents/Projects/ipo_kdd/ipo/results/A_100assetnew_10_obs_100L_100M.jld","A",A)
npzwrite("/home/syu/Documents/Projects/ipo_kdd/ipo/results/A_100assetnew_10_obs_100L_100M.npy", A)

max_stars_repo_name: cyberyu/ipo

max_issues_repo_name: cyberyu/ipo

max_forks_repo_name: cyberyu/ipo


text:
# Section 2.3.2: The dynamic table can contain duplicate entries (i.e., entries
# with the same name and same value). Therefore, duplicate entries MUST NOT be
# treated as an error by a decoder.

table = HPack.new_dynamic_table()
HPack.add_header!(table, (":method", "GET"))
HPack.add_header!(table, (":method", "GET"))

# Section 2.3.3: Indices strictly greater than the sum of the lengths of both
# tables MUST be treated as a decoding error.

table = HPack.new_dynamic_table()
HPack.set_max_table_size!(table, 1)
@test_throws DecodeError HPack.get_header(table, 100)


text:
#=
        ScalarQuadraticFunction -in- LessThan
        ScalarQuadraticFunction -in- GreaterThan
        ScalarQuadraticFunction -in- EqualTo

    TODO:
        - constraint sets
        - constraint functions
        - deleting constraints
=#
constrdict(m::LinQuadOptimizer, ::QCI{LE})  = cmap(m).q_less_than
constrdict(m::LinQuadOptimizer, ::QCI{GE})  = cmap(m).q_greater_than
constrdict(m::LinQuadOptimizer, ::QCI{EQ})  = cmap(m).q_equal_to

function MOI.addconstraint!(m::LinQuadOptimizer, func::Quad, set::S) where S <: Union{LE, GE, EQ}
    addquadraticconstraint!(m, func, set)
    m.last_constraint_reference += 1
    ref = QCI{S}(m.last_constraint_reference)
    dict = constrdict(m, ref)
    push!(m.qconstraint_primal_solution, NaN)
    push!(m.qconstraint_dual_solution, NaN)
    # dict[ref] = get_number_quadratic_constraints(m)
    dict[ref] = length(m.qconstraint_primal_solution)
    return ref
end

function addquadraticconstraint!(m::LinQuadOptimizer, func::Quad, set::S) where S<: Union{LE, GE, EQ}
    addquadraticconstraint!(m, func, backend_type(m,set), _getrhs(set))
end

function addquadraticconstraint!(m::LinQuadOptimizer, f::Quad, sense::Cchar, rhs::Float64)
    if abs(f.constant) > 0
        warn("Constant in quadratic function. Moving into set")
    end
    ri, ci, vi = reduceduplicates(
        getcol.(m, f.quadratic_rowvariables),
        getcol.(m, f.quadratic_colvariables),
        f.quadratic_coefficients
    )
    add_quadratic_constraint!(m,
        getcol.(m, f.affine_variables),
        f.affine_coefficients,
        rhs - f.constant,
        sense,
        ri, ci, vi
    )
end

function reduceduplicates(rowi::Vector{T}, coli::Vector{T}, vals::Vector{S}) where T where S
    @assert length(rowi) == length(coli) == length(vals)
    d = Dict{Tuple{T, T},S}()
    for (r,c,v) in zip(rowi, coli, vals)
        if haskey(d, (r,c))
            d[(r,c)] += v
        else
            d[(r,c)] = v
        end
    end
    ri = Vector{T}(length(d))
    ci = Vector{T}(length(d))
    vi = Vector{S}(length(d))
    for (i, (key, val)) in enumerate(d)
        ri[i] = key[1]
        ci[i] = key[2]
        vi[i] = val
    end
    ri, ci, vi
end

MOI.candelete(m::LinQuadOptimizer, c::QCI{<: LinSets}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::QCI{<: LinSets})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    row = dict[c]
    delete_quadratic_constraints!(m, row, row)
    deleteat!(m.qconstraint_primal_solution, row)
    deleteat!(m.qconstraint_dual_solution, row)
    # shift all the other references
    shift_references_after_delete_quadratic!(m, row)
    delete!(dict, c)
end


text:
#=
# Annealed importance sampling

We can compute the partition function of the RBM (and hence the log-likelihood) with
annealed importance sampling (AIS).
=#

import MLDatasets
import Makie
import CairoMakie
import RestrictedBoltzmannMachines as RBMs
using Statistics: mean, std, middle
using ValueHistories: MVHistory
using RestrictedBoltzmannMachines: Binary, BinaryRBM, initialize!, pcd!,
    aise, raise, logmeanexp, logstdexp, sample_v_from_v

# Load MNIST (0 digit only).

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .== 0] .> 0.5)
nothing #hide

# Train an RBM

rbm = BinaryRBM(Float, (28,28), 128)
initialize!(rbm, train_x)
@time pcd!(rbm, train_x; epochs=100, batchsize=128)
nothing #hide

# Get some equilibrated samples from model
v = train_x[:, :, rand(1:size(train_x, 3), 1000)]
v = sample_v_from_v(rbm, v; steps=1000)
nothing #hide

# Estimate Z with AIS and reverse AIS.

nsamples=100
ndists = [10, 100, 1000, 10_000, 100_000]
R_ais = Vector{Float64}[]
R_rev = Vector{Float64}[]
init = initialize!(Binary(zero(rbm.visible.θ)), v)
nothing #hide

for nbetas in ndists
    push!(R_ais,
        @time aise(rbm; nbetas, nsamples, init)
    )
    push!(R_rev,
        @time raise(rbm; nbetas, init, v=v[:,:,rand(1:size(v, 3), nsamples)])
    )
end

# Plots

fig = Makie.Figure()
ax = Makie.Axis(
    fig[1,1], width=700, height=400, xscale=log10, xlabel="interpolating distributions", ylabel="log(Z)"
)
Makie.band!(
    ax, ndists,
    mean.(R_ais) - std.(R_ais),
    mean.(R_ais) + std.(R_ais);
    color=(:blue, 0.25)
)
Makie.band!(
    ax, ndists,
    mean.(R_rev) - std.(R_rev),
    mean.(R_rev) + std.(R_rev);
    color=(:black, 0.25)
)
Makie.lines!(ax, ndists, mean.(R_ais); color=:blue, label="AIS")
Makie.lines!(ax, ndists, mean.(R_rev); color=:black, label="reverse AIS")
Makie.lines!(ax, ndists, logmeanexp.(R_ais); color=:blue, linestyle=:dash)
Makie.lines!(ax, ndists, logmeanexp.(R_rev); color=:black, linestyle=:dash)
Makie.lines!(ax, ndists, -logmeanexp.(-R_rev); color=:orange, linestyle=:dash)
Makie.hlines!(ax, middle(mean(R_ais[end]), mean(R_rev[end])), linestyle=:dash, color=:red, label="limiting estimate")
Makie.xlims!(extrema(ndists)...)
Makie.axislegend(ax, position=:rb)
Makie.resize_to_layout!(fig)
fig


text:
function operator_pass(x, state)
    if typof(x) === CSTParser.BinaryOpCall
        if CSTParser.precedence(x.args[2]) in (8, 13, 14, 16) || x.args[2].fullspan == 0
            ensure_no_space_after(x.args[1], state, state.offset)
            ensure_no_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        else
            ensure_single_space_after(x.args[1], state, state.offset)
            ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        end
    elseif typof(x) === CSTParser.WhereOpCall
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
        n = length(x.args)
        offset = state.offset + x.args[1].fullspan + x.args[2].fullspan
        for i = 3:n
            a = x.args[i]
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif typof(x) === CSTParser.ColonOpCall
        offset = state.offset
        n = length(x.args)
        for (i, a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif typof(x) === CSTParser.ChainOpCall || typof(x) == CSTParser.Comparison
        offset = state.offset
        n = length(x.args)
        for (i, a) in enumerate(x.args)
            if i != n
                ensure_single_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function tuple_pass(x, state)
    if typof(x) === CSTParser.TupleH
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            i == n && continue
            if typof(a) === CSTParser.PUNCTUATION && kindof(a) === Tokens.COMMA && !(typof(x.args[i + 1]) === CSTParser.PUNCTUATION)
                ensure_single_space_after(a, state, offset)
            elseif !(typof(x.args[i + 1]) === CSTParser.Parameters)
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function curly_pass(x, state)
    if typof(x) === CSTParser.Curly
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function call_pass(x, state)
    if typof(x) === CSTParser.Call
        if issameline(state.offset, state.offset + x.span, state.lines)
            offset = state.offset + x.args[1].fullspan
            n = length(x)
            for (i, a) in enumerate(x)
                i == 1 && continue
                if typof(a) === CSTParser.PUNCTUATION && kindof(a) === Tokens.COMMA
                    ensure_single_space_after(a, state, offset)
                elseif i != n && !(typof(x.args[i + 1]) === CSTParser.Parameters)
                    ensure_no_space_after(a, state, offset)
                end
                offset += a.fullspan
            end
        else
            # space holder for splitting calls across lines
        end
    elseif typof(x) === CSTParser.Kw
        ensure_single_space_after(x.args[1], state, state.offset)
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
    end
end

function forloop_pass(x, state)
    if typof(x) === CSTParser.For
        offset = state.offset + x.args[1].fullspan
        for a in x.args[2]
            # convert iter = I into iter in I
            if typof(a) === CSTParser.BinaryOpCall && CSTParser.is_eq(a.args[2])
                offset += a.args[1].fullspan
                push!(state.edits, Edit(offset + 1:offset + 2, "in "))
                offset += a.args[2].fullspan
                offset += a.args[3].fullspan
            else
                offset += a.fullspan
            end
        end
    end
end

# TODO: move this to CSTParser?
function str_value(x)
    if typof(x) === CSTParser.PUNCTUATION
        kindof(x) == Tokens.LPAREN && return "("
        kindof(x) == Tokens.LBRACE && return "{"
        kindof(x) == Tokens.LSQUARE && return "["
        kindof(x) == Tokens.RPAREN && return ")"
        kindof(x) == Tokens.RBRACE && return "}"
        kindof(x) == Tokens.RSQUARE && return "]"
        kindof(x) == Tokens.COMMA && return ","
        kindof(x) == Tokens.SEMICOLON && return ";"
        kindof(x) == Tokens.AT_SIGN && return "@"
        kindof(x) == Tokens.DOT && return "."
        return ""
    elseif typof(x) === CSTParser.IDENTIFIER || typof(x) === CSTParser.LITERAL || typof(x) === CSTParser.OPERATOR || typof(x) === CSTParser.KEYWORD
        return CSTParser.str_value(x)
    else
        s = ""
        for a in x
            s *= str_value(a)
        end
        return s
    end
end

function doc_pass(x, state)
    return
    if typof(x) === CSTParser.MacroCall && typof(x.args[1]) === CSTParser.GlobalRefDoc
        # Align global docstring to:
        #
        # """
        # doc
        # """
        #
        # If the doc is single quoted i.e. "doc", they will be replaced with triple quotes.
        offset = state.offset + x.args[1].fullspan
        doc = x.args[2]

        val = str_value(doc)

        # s = escape_string(strip(val, ['\n']), "\$")
        s = strip(val, ['\n'])
        ds = string("\"\"\"\n", s, "\n", "\"\"\"\n")

        # Check if docstring needs to be edited
        if length(ds) != doc.fullspan || s != val
            # Remove previous docstring
            push!(state.edits, Edit(offset + 1:offset + doc.fullspan, ""))
            # Append newly formatted docstring
            push!(state.edits, Edit(offset, ds))
        end
    end
end

function kw_pass(x, state)
    if typof(x) === CSTParser.KEYWORD &&
        kindof(x) in (CSTParser.Tokens.ABSTRACT,
                      CSTParser.Tokens.BAREMODULE,
                      CSTParser.Tokens.CONST,
                      CSTParser.Tokens.DO,
                      CSTParser.Tokens.ELSEIF,
                      CSTParser.Tokens.EXPORT,
                      CSTParser.Tokens.FOR,
                      CSTParser.Tokens.FUNCTION,
                      CSTParser.Tokens.GLOBAL,
                      CSTParser.Tokens.IF,
                      CSTParser.Tokens.IMPORT,
                      CSTParser.Tokens.LOCAL,
                      CSTParser.Tokens.MACRO,
                      CSTParser.Tokens.MODULE,
                      CSTParser.Tokens.MUTABLE,
                      CSTParser.Tokens.OUTER,
                      CSTParser.Tokens.PRIMITIVE,
                      CSTParser.Tokens.STRUCT,
                      CSTParser.Tokens.TYPE,
                      CSTParser.Tokens.USING,
                      CSTParser.Tokens.WHILE)
        ensure_exactly_single_space_after(x, state, state.offset)
    end
end

function comments_pass(text, state)
    ts = tokenize(text)
    while !Tokenize.Lexers.eof(ts)
        t = Tokenize.Lexers.next_token(ts)
        if Tokens.kind(t) == Tokens.COMMENT
            val = Vector{UInt8}(t.val)
            if length(val) > 1 && val[2] == 0x3d # multiline
                if !(val[3] in (0x20, 0x09)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 2, " "))
                end
                if length(val) > 5 # ensure single space at
                    for i = length(val) - 2:-1:3
                        if val[i] in (0x20, 0x09, 0x0a, 0x0d)
                            continue
                        else
                            push!(state.edits, Edit(t.startbyte .+ (i + 1:length(val) - 2), " "))
                            break
                        end
                    end
                end

            elseif length(val) > 1
                if !(val[2] in (0x20, 0x09, 0x23)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 1, " "))
                end
            end
        end
    end
end

function lineends_pass(text, x, state)
    n = sizeof(text)
    io = IOBuffer(reverse(text))
    while !eof(io)
        c = read(io, Char)
        if c === '\n' && !eof(io)
            Base.peek(io) == 0x0d && read(io, Char) # crlf
            i1 = i2 = position(io)
            pc = Base.peek(io)
            while !eof(io) && pc in (0x20, 0x09)
                i2 = position(io)
                pc = read(io, UInt8)
            end
            if i1 != i2 && (y = get_expr(x, n - i1); y isa CSTParser.EXPR ?
                !(y.typ == CSTParser.LITERAL && y.kind in (CSTParser.Tokens.STRING, CSTParser.Tokens.TRIPLE_STRING, CSTParser.Tokens.CMD, CSTParser.Tokens.TRIPLE_CMD)) : true)
                push!(state.edits, Edit((n - i2) + 1:(n - i1), ""))
            end
        end
    end
end


text:
"""
All types in the package are part of the abstract type `SimpleSDMLayer`. A
`SimpleSDMLayer` has five core fields: `grid` is a matrix storing the cells, and
`left`, `right`, `bottom` and `top` are floating point numbers specifying the
bounding box.

It is assumed that the missing values will be represented as `nothing`, so
internally the matrix will have type `Union{T, Nothing}`.
"""
abstract type SimpleSDMLayer end

"""
A predictor is a `SimpleSDMLayer` that is immutable, and so does not have
methods for `setindex!`, etc. It is a safe way to store values that should not
be modified by the analysis. Note that if you are in a bind, the values of the
`grid` field are not immutable, but don't tell anyone we told you. The correct
way of handling predictors you need to modify would be to use `convert` methods.
"""
struct SimpleSDMPredictor{T} <: SimpleSDMLayer
    grid::Matrix{Union{Nothing,T}}
    left::AbstractFloat
    right::AbstractFloat
    bottom::AbstractFloat
    top::AbstractFloat
    function SimpleSDMPredictor(grid::Matrix{Union{Nothing,T}}, l::K, r::K, b::K, t::K) where {T, K<:AbstractFloat}
        return new{T}(grid, l, r, b, t)
    end
end

"""
A response is a `SimpleSDMLayer` that is mutable, and is the usual type to store
analysis outputs. You can transform a response into a predictor using `convert`.
"""
mutable struct SimpleSDMResponse{T} <: SimpleSDMLayer
    grid::Matrix{Union{Nothing,T}}
    left::AbstractFloat
    right::AbstractFloat
    bottom::AbstractFloat
    top::AbstractFloat
    function SimpleSDMResponse(grid::Matrix{Union{Nothing,T}}, l::K, r::K, b::K, t::K) where {T, K<:AbstractFloat}
        return new{T}(grid, l, r, b, t)
    end
end

# Begin code generation for the constructors

simplesdm_types = (:SimpleSDMResponse, :SimpleSDMPredictor)

for simplesdm_type in simplesdm_types
    eval(quote
        """
            $($simplesdm_type)(grid::Matrix{Union{Nothing,T}}) where {T}

        Returns a `$($simplesdm_type)` spanning the entire globe.
        """
        function $simplesdm_type(grid::Matrix{Union{Nothing,T}}) where {T}
            return $simplesdm_type(grid, -180.0, 180.0, -90.0, 90.0)
        end

        """
            $($simplesdm_type)(grid::Matrix{Union{Nothing,T}}) where {T}

        Returns a `$($simplesdm_type)` spanning the entire globe by converting to the
        correct type, *i.e.* with `Nothing` as an acceptable value.
        """
        function $simplesdm_type(grid::Matrix{T}) where {T}
            return $simplesdm_type(convert(Matrix{Union{Nothing,T}}, grid), -180.0, 180.0, -90.0, 90.0)
        end

        function $simplesdm_type(grid::Matrix{T}, l::K, r::K, b::K, t::K) where {T, K<:AbstractFloat}
            return $simplesdm_type(convert(Matrix{Union{Nothing,T}}, grid), l, r, b, t)
        end

        function $simplesdm_type(grid::Matrix{T}, L::K) where {T, K<:SimpleSDMLayer}
            return $simplesdm_type(convert(Matrix{Union{Nothing,T}}, grid), L.left, L.right, L.bottom, L.top)
        end
    end)
end


text:
# ---
# title: 1432. Max Difference You Can Get From Changing an Integer
# id: problem1432
# author: Indigo
# date: 2021-02-06
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/max-difference-you-can-get-from-changing-an-integer/description/>
# hidden: true
# ---
#
# You are given an integer `num`. You will apply the following steps exactly
# **two** times:
#
#   * Pick a digit `x (0 <= x <= 9)`.
#   * Pick another digit `y (0 <= y <= 9)`. The digit `y` can be equal to `x`.
#   * Replace all the occurrences of `x` in the decimal representation of `num` by `y`.
#   * The new integer **cannot** have any leading zeros, also the new integer **cannot** be 0.
#
# Let `a` and `b` be the results of applying the operations to `num` the first
# and second times, respectively.
#
# Return _the max difference_ between `a` and `b`.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 555
#     Output: 888
#     Explanation: The first time pick x = 5 and y = 9 and store the new integer in a.
#     The second time pick x = 5 and y = 1 and store the new integer in b.
#     We have now a = 999 and b = 111 and max difference = 888
#
#
# **Example 2:**
#
#
#
#     Input: num = 9
#     Output: 8
#     Explanation: The first time pick x = 9 and y = 9 and store the new integer in a.
#     The second time pick x = 9 and y = 1 and store the new integer in b.
#     We have now a = 9 and b = 1 and max difference = 8
#
#
# **Example 3:**
#
#
#
#     Input: num = 123456
#     Output: 820000
#
#
# **Example 4:**
#
#
#
#     Input: num = 10000
#     Output: 80000
#
#
# **Example 5:**
#
#
#
#     Input: num = 9288
#     Output: 8700
#
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 10^8`
#
#
## @lc code=start
using LeetCode

function max_score_1432(card_points::Vector{Int}, k::Int)
    len, tt = length(card_points), sum(card_points)
    i, j = 1, len - k
    minn = sum(card_points[i:j])
    tmp = minn
    while j < len
        tmp -= card_points[i]
        i += 1
        tmp += card_points[j += 1]
        minn = min(tmp, minn)
    end
    tt - minn
end
## @lc code=end


text:
using PCA
using Test

@testset "PCA.jl" begin

end


text:
function _create_group(h5::H5DataStore, name)
    if haskey(h5, name)
        g = h5[name]
    else
        g = create_group(h5, name)
    end
    return g
end


text:
using Plots

# Specify which experiment you want
alpha = .6
seedtype = 1
# Load the data for that experiment
using MAT
ppr_string = "../PPR_Test_Output/Seed_"*string(seedtype)*"_alpha_"*string(round(Int64,alpha*100))*"sets_stats.mat"
matPPR = matread(ppr_string)

# stats = [tol;numel(PR);tPR;cond;pPR;rPR;fPR];
ppr_stats = matPPR["stats"]

fs_string = "../FlowSeed_Test_Output/Seed_"*string(seedtype)*".mat"
matFS = matread(fs_string)

# [Ssize;timer;Sstats[4];Rcond;pr;re;f1]
fs_stats = matFS["stats"]

flowseed_f1 = fs_stats[7,1:2:end]
ppr_f1 = ppr_stats[7,1:2:end]

mat = matread("../Brain_Regions/BrainRegions_Test.mat")
Sizes = vec(round.(Int64,mat["RegionSizeTest"]))
Sizes = Sizes[1:2:end]

x_label = "Number of Nodes in Region"
y_label = "F1 Score"
numregions = length(ppr_f1)
l_place = :bottomleft
lstyle = :solid
lw = 1.5
s1 = 400
s2 = 300
ms = 3
xt = 1:3:numregions
sz = Sizes[xt]

plot(flowseed_f1, labels = "FlowSeed",grid = false,size = (s1,s2),
    xlabel = x_label,ylabel = y_label,legend = l_place,linewidth = lw,
    xticks = (xt,sz), xrotation=50,markershape = :circle,markersize = ms,
    linestyle = lstyle,color = :blue,ylim = [0,1]
)

plot!(ppr_f1, labels = "Push",grid = false,size = (s1,s2),
    xlabel = x_label,ylabel = y_label,legend = l_place,linewidth = lw,
    xticks = (xt,sz), xrotation=50,markershape = :cross,markersize = ms,
    linestyle = lstyle, color = :green,ylim = [0,1]
)

savefig("TestPlots_Seed_"*string(seedtype)*".pdf")

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
@testset "SplitHLAData" begin
    data = Escape.HLADataset("Test").data[1]
    @test Escape.SplitHLAData(data.name, data.records, data.hla_types, data.tree,
        data.stan_input, 1, 1, data, 1:14, 1) isa Escape.SplitHLAData

    records = Escape.records(joinpath(@__DIR__, "data", "test.fasta"))
    hla_types = rand(HLAType, 5)

    split = Escape.SplitHLAData("test", records, hla_types, data.tree, data.stan_input,
        1, 1, data, 1:14, 1)
    @test typeof(split) <: AbstractHLAData

    wrong_tree = phylogenetic_tree("(X:1.2,5:1.5,(1:0.5,(2:1.6,3:0.7):0.5):0.6);")
    wrong_tree_2 = phylogenetic_tree("(3:0.3,4:0.4)5:0.5);")

    @test_throws ErrorException Escape.SplitHLAData("test", records, hla_types, wrong_tree,
        data.stan_input, 1, 1, data, 1:14, 1)
    @test_throws DimensionMismatch Escape.SplitHLAData("test", records, hla_types,
        wrong_tree_2, data.stan_input, 1, 1, data, 1:14, 1)
end

@testset "split_hla_data(::HLAData, ::Int)" begin
    data = Escape.HLADataset("Test").data[1]
    r = Escape.replacements(data)

    @test Escape.split_hla_data(data, 2) isa Array{Escape.SplitHLAData}
    @test length(Escape.split_hla_data(data, 2)) == 2
    @test Escape.replacements(Escape.split_hla_data(data, 2)[1]) == r[1:7]
end


text:
"""
    eom_nonlinear3d(uᵢ, a₁, a₂, a₃,  b₁, b₂, b₃,
        c₁₂, c₂₃, c₁₃, σ₁, σ₂, σ₃) -> DiscreteDynamicalSystem

Equations of motion for a 3d nonlinear system with nonlinear couplings
``x_1 \\to x_2``, ``x_2 \\to x_3`` and ``x_1 \\to x_3``. Modified from [1].

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x_1(t+1) &= a_1 x_1 (1-x_1(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{1}(t) \\\\
x_2(t+1) &= a_1 x_2 (1-x_2(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{2}(t) + b x_1 x_2 \\\\
x_3(t+1) &= a_3 x_3 (1-x_3(t))^2  e^{-x_3(t)^2} + 0.4 \\xi_{3}(t) + c x_{2}(t) + d x_{1}(t)^2.
\\end{aligned}
```

## References

1. Gourévitch, B., Le Bouquin-Jeannès, R., & Faucon, G. (2006). Linear and nonlinear
    causality between signals: methods, examples and neurophysiological
    applications. Biological Cybernetics, 95(4), 349–369.
"""
function eom_nonlinear3d(x, p, n)
    x₁, x₂, x₃ = (x...,)
    a₁, a₂, a₃, b₁, b₂, b₃, c₁₂, c₂₃, c₁₃, σ₁, σ₂, σ₃ = (p...,)
    ξ₁ = rand(Normal(0, σ₁))
    ξ₂ = rand(Normal(0, σ₂))
    ξ₃ = rand(Normal(0, σ₃))

    dx₁ = a₁*x₁*(1-x₁)^2 * exp(-x₁^2) + b₁*ξ₁
    dx₂ = a₂*x₂*(1-x₂)^2 * exp(-x₂^2) + b₂*ξ₂ + c₁₂*x₁*x₂
    dx₃ = a₃*x₃*(1-x₃)^2 * exp(-x₃^2) + b₃*ξ₃ + c₂₃*x₂ + c₁₃*x₁^2

    return SVector{3}(dx₁, dx₂, dx₃)
end

function nonlinear3d(uᵢ, a₁, a₂, a₃,  b₁, b₂, b₃, c₁₂, c₂₃, c₁₃, σ₁, σ₂, σ₃)
    p = [a₁, a₂, a₃,  b₁, b₂, b₃, c₁₂, c₂₃, c₁₃, σ₁, σ₂, σ₃]
    s = DiscreteDynamicalSystem(eom_nonlinear3d, uᵢ, p)
    return s
end

"""
    nonlinear3d(;uᵢ = rand(3),
        σ₁ = 1.0, σ₂ = 1.0, σ₃ = 1.0,
        a₁ = 3.4, a₂ = 3.4, a₃ = 3.4,
        b₁ = 0.4, b₂ = 0.4, b₃ = 0.4,
        c₁₂ = 0.5, c₂₃ = 0.3, c₁₃ = 0.5) -> DiscreteDynamicalSystem

A 3d nonlinear system with nonlinear couplings ``x_1 \\to x_2``,
``x_2 \\to x_3`` and ``x_1 \\to x_3``. Modified from [1].

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x_1(t+1) &= a_1 x_1 (1-x_1(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{1}(t) \\\\
x_2(t+1) &= a_1 x_2 (1-x_2(t))^2  e^{-x_2(t)^2} + 0.4 \\xi_{2}(t) + b x_1 x_2 \\\\
x_3(t+1) &= a_3 x_3 (1-x_3(t))^2  e^{-x_3(t)^2} + 0.4 \\xi_{3}(t) + c x_{2}(t) + d x_{1}(t)^2.
\\end{aligned}
```

## References

1. Gourévitch, B., Le Bouquin-Jeannès, R., & Faucon, G. (2006). Linear and nonlinear
    causality between signals: methods, examples and neurophysiological
    applications. Biological Cybernetics, 95(4), 349–369.
"""
nonlinear3d(;uᵢ = rand(3),
        σ₁ = 1.0, σ₂ = 1.0, σ₃ = 1.0,
        a₁ = 3.4, a₂ = 3.4, a₃ = 3.4,
        b₁ = 0.4, b₂ = 0.4, b₃ = 0.4,
        c₁₂ = 0.5, c₂₃ = 0.3, c₁₃ = 0.5) =
    nonlinear3d(uᵢ, a₁, a₂, a₃,  b₁, b₂, b₃, c₁₂, c₂₃, c₁₃, σ₁, σ₂, σ₃)


text:
module HDFTables

using HDF5
using DataFrames

struct Table{T,N}
    columns::NTuple{N,Symbol}
    rows::Vector{T}
end
export Table

function table_type(obj::HDF5Dataset)
    """ Builds a Tuple type for a compound dataset """
    objtype = HDF5.datatype(obj)
    class_id = HDF5.h5t_get_class(objtype.id)
    if class_id == HDF5.H5T_COMPOUND
        n = Int(HDF5.h5t_get_nmembers(objtype.id))
        field_names = Symbol[]
        field_types = DataType[]
        for i in 0:n-1
            field_tid = HDF5.h5t_get_member_type(objtype.id, i)
            field_type = HDF5.hdf5_to_julia_eltype(HDF5Datatype(field_tid))
            field_name = HDF5.h5t_get_member_name(objtype.id, i)
            push!(field_names, Symbol(field_name))
            push!(field_types, field_type)
        end
        T = Tuple{field_types...}
        return T, tuple(field_names...)
    else
        error("Dataset is not a table: $obj")
    end
end

function read_table(obj::HDF5Dataset)
    # Build a Julia Tuple type for the table elements:
    T, field_names = table_type(obj)
    @assert isbits(T)

    # Build a HDF5 datatype that reflects the memory layout of the Julia type:
    memtype_id = HDF5.h5t_create(HDF5.H5T_COMPOUND, sizeof(T))
    try
        for (i, ftype) in enumerate(T.types)
            HDF5.h5t_insert(memtype_id, string(field_names[i]), fieldoffset(T,i), HDF5.datatype(ftype))
        end

        # Allocate the destination array:
        space = HDF5.dataspace(obj);
        dims, maxdims = HDF5.get_dims(space)
        data = Array{T}(dims...)

        # Read the data using HDF5's built-in data layout conversion:
        status = HDF5.readarray(obj, memtype_id, data)
        status < 0 && error("Couldn't read table: $status")

        return Table(field_names, data)
    finally
        HDF5.h5t_close(memtype_id)
    end
end
export read_table

column(table::Table, i::Int) = collect(row[i] for row in table.rows)

function Base.convert(::Type{DataFrame}, t::Table)
    cols = [column(t, i) for i in eachindex(t.columns)]
    DataFrame(cols, Symbol[t.columns...])
end

end # module


text:
using Colors_2021
using Test

@testset "Colors_2021.jl" begin
    # Write your tests here.
end


text:
# Module for constructing lattices
# Things to add:
#   - support for unit cell
#   - support for generic lattices (so far only have square)

module Lattices

using LinearAlgebra
using LightGraphs
using SimpleWeightedGraphs
import LightGraphs.cartesian_product

abstract type Lattice end

export Lattice

struct SquareLattice <: Lattice
    dims::Tuple
    graph::LightGraphs.SimpleGraph{Int64}
    lattice_vectors::Array{Int64}
end

struct MixedBCSquareLattice <: Lattice
    dims::Tuple
    graph::SimpleWeightedGraphs.SimpleWeightedGraph
end

function get_MixedBoundarySquareLattice(dims::Tuple)::MixedBCSquareLattice
    """ This is an extension of the function Grid from LightGraphs to
    make a weighted graph with antiperiodic boundary conditions in one direction
    and periodic boundary conditions in the other, for the purposes
    of obtaining a nondegenerate fermi surface.
    https://journals.jps.jp/doi/pdf/10.1143/JPSJ.57.2482 """

   # make anti periodic cycle in x direction
   g = SimpleWeightedGraph(CycleGraph(dims[1]))
   add_edge!(g,1,dims[1],-1)

   for d in dims[2:end]
       g = cartesian_product(SimpleWeightedGraph(CycleGraph(d)), g)
   end
   return MixedBCSquareLattice(dims,g)
end

function get_SquareLattice(dims;pbc=true)::SquareLattice
    graph = LightGraphs.Grid([d for d in dims],periodic=pbc)
    # lattice vectors for a d dimensional square lattice are just
    # v_i = \hat e_i
    lattice_vectors = get_SquareLatticeVectors(dims)
    return SquareLattice(dims,graph,lattice_vectors)
end

function get_SquareLatticeVectors(dims::Tuple)::Array{Int64}
    """ returns array of lattice vectors
    v[i,:,:] is the i'th lattice vector"""
    ndims = length(dims)
    m0 = Matrix{Int}(I,ndims,ndims)
    return m0
end

struct TriangularLattice <: Lattice
    dims::Tuple
    graph::LightGraphs.SimpleGraph{Int64}
end

function get_Tightbinding_Wavefunctions(lattice::SquareLattice)::Eigen
    ham = -Matrix(adjacency_matrix(lattice.graph))
    F = eigen(ham) # F.values are evals, F.vectors[:,k] is k'th evec.
    return F
end

function get_Tightbinding_Wavefunctions(lattice::MixedBCSquareLattice)::Eigen
    ham = -Matrix(weights(lattice.graph))
    F = eigen(ham) # F.values are evals, F.vectors[:,k] is k'th evec.
    return F
end

function get_2DCSL_Wavefunctions(lattice::SquareLattice,Δ)::Eigen
    lg = lattice.graph
    dims = lattice.dims
    nsites = prod(dims)
    ham = zeros(ComplexF64,(nsites,nsites))
    Lx = dims[1]
    Ly = dims[2]

    for y = 1:dims[2]
        # loop over y layers. odd layers have hopping -t in the x dir
        # even layers have hopping +t in the x dir
        for x = 1:dims[1]
            id0 = sub2ind(dims,x,y)

            if (y % 2) == 1
                idp = sub2ind(dims,x,wrap(y+1,Ly))
                ham[idp,id0] = -1

                idp = sub2ind(dims,wrap(x+1,Lx),y)
                ham[idp,id0] = +1

                idp = sub2ind(dims,wrap(x+1,Lx),wrap(y+1,Ly))
                ham[idp,id0] = +1im*Δ

                idp = sub2ind(dims,wrap(x+1,Lx),wrap(y-1,Ly))
                ham[idp,id0] = -1im*Δ
            else
                idp = sub2ind(dims,x,wrap(y+1,Ly))
                ham[idp,id0] = -1

                idp = sub2ind(dims,wrap(x+1,Lx),y)
                ham[idp,id0] = -1

                idp = sub2ind(dims,wrap(x+1,Lx),wrap(y+1,Ly))
                ham[idp,id0] = -1im*Δ

                idp = sub2ind(dims,wrap(x+1,Lx),wrap(y-1,Ly))
                ham[idp,id0] = +1im*Δ
            end
        end
    end

    ham += ham'
    # print(ham)
    F = eigen(ham)
    return F
end

function wrap(x,a)
    """ wraps x to between 1 and a """
    return fldmod(x-1,a)[2]+1
end

"""transition to julia 1.0"""
function ind2sub(dims, ind)
    return Tuple(CartesianIndices(dims)[ind])
end

"""transition to julia 1.0"""
function sub2ind(dims, inds...)
    return (LinearIndices(dims))[inds...]
end

function cartesian_product(g::G, h::G) where G <: AbstractSimpleWeightedGraph
    """ Function overload for cartesian product for weighted graphs. For
    some reason, not correctly loading override functions from package..."""
    z = G(nv(g) * nv(h))
    id(i, j) = (i - 1) * nv(h) + j
    for e in edges(g)
        i1, i2 = Tuple(e)
        for j = 1:nv(h)
            add_edge!(z, id(i1, j), id(i2, j), weight(e))
        end
    end

    for e in edges(h)
        j1, j2 = Tuple(e)
        for i in vertices(g)
            add_edge!(z, id(i, j1), id(i, j2), weight(e))
        end
    end
    return z
end

end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("datasets", "iris"),
     layer(x="SepalLength", y="SepalWidth", Geom.point),
     layer([sin, cos], 0, 25))


text:
# You finally decode the Elves' message. HI, the message says. You continue
# searching for the sleigh keys.
#
# Ahead of you is what appears to be a large ocean trench. Could the keys have
# fallen into it? You'd better send a probe to investigate.
#
# The probe launcher on your submarine can fire the probe with any integer
# velocity in the x (forward) and y (upward, or downward if negative)
# directions. For example, an initial x,y velocity like 0,10 would fire the
# probe straight up, while an initial velocity like 10,-1 would fire the probe
# forward at a slight downward angle.
#
# The probe's x,y position starts at 0,0. Then, it will follow some trajectory
# by moving in steps. On each step, these changes occur in the following order:
#
# - The probe's x position increases by its x velocity.
# - The probe's y position increases by its y velocity.
# - Due to drag, the probe's x velocity changes by 1 toward the value 0; that
#   is, it decreases by 1 if it is greater than 0, increases by 1 if it is less
#   than 0, or does not change if it is already 0.
# - Due to gravity, the probe's y velocity decreases by 1.
#
# For the probe to successfully make it into the trench, the probe must be on
# some trajectory that causes it to be within a target area after any step. The
# submarine computer has already calculated this target area (your puzzle
# input). For example:
#
# target area: x=20..30, y=-10..-5
#
# This target area means that you need to find initial x,y velocity values such
# that after any step, the probe's x position is at least 20 and at most 30,
# and the probe's y position is at least -10 and at most -5.
#
# Given this target area, one initial velocity that causes the probe to be
# within the target area after any step is 7,2:
#
# .............#....#............
# .......#..............#........
# ...............................
# S........................#.....
# ...............................
# ...............................
# ...........................#...
# ...............................
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTT#TT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
#
# In this diagram, S is the probe's initial position, 0,0. The x coordinate
# increases to the right, and the y coordinate increases upward. In the bottom
# right, positions that are within the target area are shown as T. After each
# step (until the target area is reached), the position of the probe is marked
# with #. (The bottom-right # is both a position the probe reaches and a
# position in the target area.)
#
# Another initial velocity that causes the probe to be within the target area
# after any step is 6,3:
#
# ...............#..#............
# ...........#........#..........
# ...............................
# ......#..............#.........
# ...............................
# ...............................
# S....................#.........
# ...............................
# ...............................
# ...............................
# .....................#.........
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................T#TTTTTTTTT
# ....................TTTTTTTTTTT
#
# Another one is 9,0:
#
# S........#.....................
# .................#.............
# ...............................
# ........................#......
# ...............................
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTT#
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
# ....................TTTTTTTTTTT
#
# One initial velocity that doesn't cause the probe to be within the target
# area after any step is 17,-4:
#
# S..............................................................
# ...............................................................
# ...............................................................
# ...............................................................
# .................#.............................................
# ....................TTTTTTTTTTT................................
# ....................TTTTTTTTTTT................................
# ....................TTTTTTTTTTT................................
# ....................TTTTTTTTTTT................................
# ....................TTTTTTTTTTT..#.............................
# ....................TTTTTTTTTTT................................
# ...............................................................
# ...............................................................
# ...............................................................
# ...............................................................
# ................................................#..............
# ...............................................................
# ...............................................................
# ...............................................................
# ...............................................................
# ...............................................................
# ...............................................................
# ..............................................................#
#
# The probe appears to pass through the target area, but is never within it
# after any step. Instead, it continues down and to the right - only the first
# few steps are shown.
#
# If you're going to fire a highly scientific probe out of a super cool probe
# launcher, you might as well do it with style. How high can you make the probe
# go while still reaching the target area?
#
# In the above example, using an initial velocity of 6,9 is the best you can
# do, causing the probe to reach a maximum y position of 45. (Any higher
# initial y velocity causes the probe to overshoot the target area entirely.)
#
# Find the initial velocity that causes the probe to reach the highest y
# position and still eventually be within the target area after any step. What
# is the highest y position it reaches on this trajectory?

# read input
(_, _, targety1, _) = parse.(
  Int,
  match(
    r"x=(-?\d+)..(-?\d+), y=(-?\d+)..(-?\d+)",
    readline("input.txt")
  )
)

# If our probe starts at a y velocity of Y, it will reach its maximum height
# and then fall back to 0 with a velocity of -Y because, at each step, the
# velocity will be Y, Y-1, Y-2, ..., 0 (at its maximum height), -1, -2, -3,
# ..., -Y (at a y position of 0), and so on. Assuming the minimum target y is
# negative, the velocity that produces the maximum height will be
# abs(min(targety))-1. Then we can just find the maximum height as the sum of
# (Y, Y-1, Y-2, ..., 0).
inity = abs(targety1) - 1
maxheight = inity * (inity + 1) / 2
println("Maximum Height: ", maxheight)


text:
using OrdinaryDiffEq, Test

jac_called = Ref(false)
tgrad_called = Ref(false)

function Lotka(du,u,p,t)
  du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
  du[2] = -3 * u[2] + 1 * u[1] * u[2]
  nothing
end

function Lotka_jac(J,u,p,t)
  jac_called.x = true
  J[1,1] = 1.0 - u[2]
  J[1,2] = -u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

function Lotka_tgrad(grad,u,p,t)
  tgrad_called.x = true
  grad[1] = 1 * 0
  grad[2] = 1 * 0
  nothing
end

Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad)
prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))

good_sol = solve(prob,Rosenbrock23())

@test jac_called[]
@test tgrad_called[]

prob2 = ODEProblem(Lotka,ones(2),(0.0,10.0))

sol = solve(prob2,Rosenbrock23(autodiff=true))
@test ≈(good_sol[end],sol[end],rtol=1e-2)

sol = solve(prob2,Rosenbrock23(autodiff=true,chunk_size=1))
@test ≈(good_sol[end],sol[end],rtol=1e-2)

sol = solve(prob2,Rosenbrock23(autodiff=false))
@test ≈(good_sol[end],sol[end],rtol=1e-2)


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

function glavnaia(r)
    ygolok(r, West, Sud)
    n = 0
    flag = 1
    while !isborder(r, Nord)
        while !isborder(r, Ost)
            move!(r, Ost)
            if(flag == 0 && isborder(r, Nord))
                flag = 1
                n += 1
            elseif(flag == 1 && !isborder(r, Nord))
                flag = 0
            end
        end
        move!(r, Nord)
        kray(r,West)
    end
print(n)
end

function ygolok(r,side1,side2)
    while !isborder(r,side1) || !isborder(r,side2)
        if !isborder(r,side1)
            move!(r,side1)
        else
            move!(r,side2)
        end
    end
end
function kray(r,side) #функция перемещает робота к краю
    while isborder(r,side) == false
        move!(r,side)
    end
end

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs