text:
const dev = CUDAdrv.device()
const ctx = CUDAdrv.CuContext(dev)

macro thread_local_index()
    return :(
        (threadIdx().z - 1) * blockDim().y * blockDim().x +
        (threadIdx().y - 1) * blockDim().x +
        threadIdx().x
    )
end

macro total_threads_per_block()
    return :(blockDim().z * blockDim().y * blockDim().x)
end

macro block_index()
    return :(
        blockIdx().x +
        (blockIdx().y - 1) * gridDim().x +
        (blockIdx().z - 1) * gridDim().x * gridDim().y
    )
end

macro total_blocks()
    return :(gridDim().z * gridDim().x * gridDim().y)
end

macro thread_global_index()
    return :(
        (@block_index() - 1) * (blockDim().x * blockDim().y * blockDim().z) +
        @thread_local_index()
    )
end

macro total_threads()
    return :(@total_blocks() * @total_threads_per_block())
end

function callkernel(dev, kernel, args)
    blocks, threads = getvalidconfig(dev, kernel, args)
    @cuda blocks = blocks threads = threads kernel(args...)

    return nothing
end

function getvalidconfig(dev, kernel, parallel_args)
    R = parallel_args[1]
    Rlength = length(R)
    Ssize = size(R)
    Slength = prod(Ssize)
    GC.@preserve parallel_args begin
        parallel_kargs = cudaconvert.(parallel_args)
        parallel_tt = Tuple{Core.Typeof.(parallel_kargs)...}
        parallel_kernel = cufunction(kernel, parallel_tt)

        # we are limited in how many threads we can launch...
        ## by the kernel
        kernel_threads = CUDAnative.maxthreads(parallel_kernel)
        ## by the device
        block_threads = (
            x=CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_X),
            y=CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_Y),
            total=CUDAdrv.attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK),
        )

        # figure out a legal launch configuration
        y_thr = min(nextpow(2, Rlength ÷ 512 + 1), 512, block_threads.y, kernel_threads)
        x_thr = min(
            512 ÷ y_thr,
            Slength,
            block_threads.x,
            ceil(Int, block_threads.total / y_thr),
            ceil(Int, kernel_threads / y_thr),
        )

        blk, thr = (Rlength - 1) ÷ y_thr + 1, (x_thr, y_thr, 1)
        blk = min(blk, ceil(Int, Rlength / prod(thr)))
    end
    return blk, thr
end


text:
# A quadratic quadrilateral in 2D space.
struct Quadrilateral8_2D <: Face_2D
    # The points are assumed to be ordered  in counter clockwise order as follows
    # p₁ = vertex A
    # p₂ = vertex B
    # p₃ = vertex C
    # p₄ = vertex D
    # p₅ = point on the quadratic segment from A to B
    # p₆ = point on the quadratic segment from B to C
    # p₇ = point on the quadratic segment from C to D
    # p₈ = point on the quadratic segment from D to A
    points::SVector{8, Point_2D}
end

# Constructors
# -------------------------------------------------------------------------------------------------
Quadrilateral8_2D(p₁::Point_2D, p₂::Point_2D, p₃::Point_2D, p₄::Point_2D,
                  p₅::Point_2D, p₆::Point_2D, p₇::Point_2D, p₈::Point_2D
                 ) = Quadrilateral8_2D(SVector(p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))

# Methods
# -------------------------------------------------------------------------------------------------
function (quad8::Quadrilateral8_2D)(r::Real, s::Real)
    # See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
    # Chapter 8, Advanced Data Representation, in the interpolation functions section
    ξ = 2*Float64(r) - 1; η = 2*Float64(s) - 1
    return (1 - ξ)*(1 - η)*(-ξ - η - 1)/4*quad8[1] +
           (1 + ξ)*(1 - η)*( ξ - η - 1)/4*quad8[2] +
           (1 + ξ)*(1 + η)*( ξ + η - 1)/4*quad8[3] +
           (1 - ξ)*(1 + η)*(-ξ + η - 1)/4*quad8[4] +
                      (1 - ξ^2)*(1 - η)/2*quad8[5] +
                      (1 - η^2)*(1 + ξ)/2*quad8[6] +
                      (1 - ξ^2)*(1 + η)/2*quad8[7] +
                      (1 - η^2)*(1 - ξ)/2*quad8[8]
end

# Interpolation with a point, instead of (r,s)
function (quad8::Quadrilateral8_2D)(p::Point_2D)
    r, s = p
    ξ = 2r - 1; η = 2s - 1
    return (1 - ξ)*(1 - η)*(-ξ - η - 1)/4*quad8[1] +
           (1 + ξ)*(1 - η)*( ξ - η - 1)/4*quad8[2] +
           (1 + ξ)*(1 + η)*( ξ + η - 1)/4*quad8[3] +
           (1 - ξ)*(1 + η)*(-ξ + η - 1)/4*quad8[4] +
                      (1 - ξ^2)*(1 - η)/2*quad8[5] +
                      (1 - η^2)*(1 + ξ)/2*quad8[6] +
                      (1 - ξ^2)*(1 + η)/2*quad8[7] +
                      (1 - η^2)*(1 - ξ)/2*quad8[8]
end

area(quad8::Quadrilateral8_2D) = area(quad8, Val(2))

function area(quad8::Quadrilateral8_2D, ::Val{N}) where {N}
    # Numerical integration required. Gauss-Legendre quadrature over a quadrilateral is used.
    # Let Q(r,s) be the interpolation function for quad8,
    #                             1  1
    # A = ∬ ||∂Q/∂r × ∂Q/∂s||dA = ∫  ∫ ||∂Q/∂r × ∂Q/∂s|| ds dr
    #      D                      0  0
    #
    #       N   N
    #   =   ∑   ∑  wᵢwⱼ||∂Q/∂r(rᵢ,sⱼ) × ∂Q/∂s(rᵢ,sⱼ)||
    #      i=1 j=1
    # N is the square root of the number of points used in the quadrature.
    # See tuning/Quadrilateral8_2D_area.jl for more info on how N = 3 was chosen.
    w, r = gauss_legendre_quadrature(Val(N))
    a = 0.0
    for i = 1:N, j = 1:N
        ∂Q_∂r, ∂Q_∂s = gradient(quad8, r[i], r[j])
        a += w[i]*w[j]*abs(∂Q_∂r × ∂Q_∂s)
    end
    return a
end

centroid(quad8::Quadrilateral8_2D) = centroid(quad8, Val(3))
function centroid(quad8::Quadrilateral8_2D, ::Val{N}) where {N}
    # Numerical integration required. Gauss-Legendre quadrature over a quadrilateral is used.
    # Let Q(r,s) be the interpolation function for quad8,
    #                             1  1
    # A = ∬ ||∂Q/∂r × ∂Q/∂s||dA = ∫  ∫ ||∂Q/∂r × ∂Q/∂s|| ds dr
    #      D                      0  0
    #
    #       N   N
    #   =   ∑   ∑  wᵢwⱼ||∂Q/∂r(rᵢ,sⱼ) × ∂Q/∂s(rᵢ,sⱼ)||
    #      i=1 j=1
    #
    # C_x = (∫∫ x dA)/A, C_y = (∫∫ y dA)/A
    #         D                  D
    w, r = gauss_legendre_quadrature(Val(N))
    A = 0.0
    C = Point_2D()
    for i = 1:N, j = 1:N
        ∂Q_∂r, ∂Q_∂s = gradient(quad8, r[i], r[j])
        weighted_val = w[i]*w[j]*abs(∂Q_∂r × ∂Q_∂s)
        A += weighted_val
        C += weighted_val * quad8(r[i], r[j])
    end
    return C/A
end

function gradient(quad8::Quadrilateral8_2D, r::Real, s::Real)
    # Chain rule
    # ∂Q   ∂Q ∂ξ     ∂Q      ∂Q   ∂Q ∂η     ∂Q
    # -- = -- -- = 2 -- ,    -- = -- -- = 2 --
    # ∂r   ∂ξ ∂r     ∂ξ      ∂s   ∂η ∂s     ∂η
    ξ = 2Float64(r) - 1; η = 2Float64(s) - 1
    ∂Q_∂ξ = (1 - η)*(2ξ + η)/4*quad8[1] +
            (1 - η)*(2ξ - η)/4*quad8[2] +
            (1 + η)*(2ξ + η)/4*quad8[3] +
            (1 + η)*(2ξ - η)/4*quad8[4] +
                    -ξ*(1 - η)*quad8[5] +
                   (1 - η^2)/2*quad8[6] +
                    -ξ*(1 + η)*quad8[7] +
                  -(1 - η^2)/2*quad8[8]

    ∂Q_∂η = (1 - ξ)*( ξ + 2η)/4*quad8[1] +
            (1 + ξ)*(-ξ + 2η)/4*quad8[2] +
            (1 + ξ)*( ξ + 2η)/4*quad8[3] +
            (1 - ξ)*(-ξ + 2η)/4*quad8[4] +
                   -(1 - ξ^2)/2*quad8[5] +
                     -η*(1 + ξ)*quad8[6] +
                    (1 - ξ^2)/2*quad8[7] +
                     -η*(1 - ξ)*quad8[8]

    return 2*∂Q_∂ξ, 2*∂Q_∂η
end

function in(p::Point_2D, quad8::Quadrilateral8_2D)
    # If the point is to the left of every edge
    #  4<-----3
    #  |      ^
    #  | p    |
    #  |      |
    #  |      |
    #  v----->2
    #  1
    return isleft(p, QuadraticSegment_2D(quad8[1], quad8[2], quad8[5])) &&
           isleft(p, QuadraticSegment_2D(quad8[3], quad8[4], quad8[7])) &&
           isleft(p, QuadraticSegment_2D(quad8[2], quad8[3], quad8[6])) &&
           isleft(p, QuadraticSegment_2D(quad8[4], quad8[1], quad8[8]))
end

# function in(p::Point_2D{F}, quad8::Quadrilateral8_2D{F}, N::Int64) where {F <: AbstractFloat}
#     # Determine if the point is in the triangle using the Newton-Raphson method
#     # N is the max number of iterations of the method.
#     p_rs = real_to_parametric(p, quad8, N)
#     ϵ = parametric_coordinate_ϵ
#     if (-ϵ ≤ p_rs[1] ≤ 1 + ϵ) &&
#        (-ϵ ≤ p_rs[2] ≤ 1 + ϵ)
#         return true
#     else
#         return false
#     end
# end

function intersect(l::LineSegment_2D, quad8::Quadrilateral8_2D)
    # Create the 3 quadratic segments that make up the triangle and intersect each one
    edges = SVector(QuadraticSegment_2D(quad8[1], quad8[2], quad8[5]),
                    QuadraticSegment_2D(quad8[2], quad8[3], quad8[6]),
                    QuadraticSegment_2D(quad8[3], quad8[4], quad8[7]),
                    QuadraticSegment_2D(quad8[4], quad8[1], quad8[8]))
    ipoints = MVector(Point_2D(), Point_2D(), Point_2D(),
                      Point_2D(), Point_2D(), Point_2D())
    n_ipoints = 0x00000000
    # We need to account for 6 points returned
    for k ∈ 1:4
        npoints, points = l ∩ edges[k]
        for i = 1:npoints
            n_ipoints += 0x00000001
            ipoints[n_ipoints] = points[i]
        end
    end
    return n_ipoints, SVector(ipoints)
end

function jacobian(quad8::Quadrilateral8_2D, r::Real, s::Real)
    # Return the 2 x 2 Jacobian matrix
    ∂Q_∂r, ∂Q_∂s = gradient(quad8, r, s)
    return SMatrix{2, 2}(∂Q_∂r.x, ∂Q_∂r.y,
                         ∂Q_∂s.x, ∂Q_∂s.y)
end

function real_to_parametric(p::Point_2D, quad8::Quadrilateral8_2D)
    return real_to_parametric(p, quad8, 30)
end

function real_to_parametric(p::Point_2D, quad8::Quadrilateral8_2D, N::Int64)
    # Convert from real coordinates to the triangle's local parametric coordinates using the
    # the Newton-Raphson method. N is the max number of iterations
    # If a conversion doesn't exist, the minimizer is returned.
    r = 0.5 # Initial guess at centroid
    s = 0.5
    for i = 1:N
        err = p - quad8(r, s)
        # Inversion is faster for 2 by 2 than \
        Δr, Δs = inv(jacobian(quad8, r, s)) * err
        r += Δr
        s += Δs
        if abs(Δr) + abs(Δs) < 1e-6
            break
        end
    end
    return Point_2D(r, s)
end

function triangulate(quad8::Quadrilateral8_2D, N::Int64)
    # N is the number of divisions of each edge
    triangles = Vector{Triangle_2D}(undef, 2*(N+1)*(N+1))
    if N === 0
        triangles[1] = Triangle_2D(quad8[1], quad8[2], quad8[3])
        triangles[2] = Triangle_2D(quad8[3], quad8[4], quad8[1])
    else
        for j = 0:N, i = 0:N
            triangles[2*(N+1)*j + 2i + 1] = Triangle_2D(quad8(    i/(N+1),     j/(N+1)),
                                                        quad8((i+1)/(N+1),     j/(N+1)),
                                                        quad8(    i/(N+1), (j+1)/(N+1)))
            triangles[2*(N+1)*j + 2i + 2] = Triangle_2D(quad8(    i/(N+1), (j+1)/(N+1)),
                                                        quad8((i+1)/(N+1),     j/(N+1)),
                                                        quad8((i+1)/(N+1), (j+1)/(N+1)))
        end
    end
    return triangles
end

# Plot
# -------------------------------------------------------------------------------------------------
if enable_visualization
    function convert_arguments(LS::Type{<:LineSegments}, quad8::Quadrilateral8_2D)
        q₁ = QuadraticSegment_2D(quad8[1], quad8[2], quad8[5])
        q₂ = QuadraticSegment_2D(quad8[2], quad8[3], quad8[6])
        q₃ = QuadraticSegment_2D(quad8[3], quad8[4], quad8[7])
        q₄ = QuadraticSegment_2D(quad8[4], quad8[1], quad8[8])
        qsegs = [q₁, q₂, q₃, q₄]
        return convert_arguments(LS, qsegs)
    end

    function convert_arguments(LS::Type{<:LineSegments}, Q::Vector{Quadrilateral8_2D})
        point_sets = [convert_arguments(LS, quad8) for quad8 ∈ Q]
        return convert_arguments(LS, reduce(vcat, [pset[1] for pset ∈ point_sets]))
    end

    function convert_arguments(M::Type{<:Mesh}, quad8::Quadrilateral8_2D)
        triangles = triangulate(quad8, 13)
        return convert_arguments(M, triangles)
    end

    function convert_arguments(M::Type{<:Mesh}, Q::Vector{Quadrilateral8_2D})
        triangles = reduce(vcat, triangulate.(Q, 13))
        return convert_arguments(M, triangles)
    end
end

ext: jl

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
module test_unsplashimages_unsplash_image

using Test
using UnsplashImages

yannic= unsplash_image("yannic-laderach-Dqx4XWuXu7w-unsplash.small.jpg")
@test size(yannic) == (800, 640)

yannic_mask = unsplash_image("yannic-laderach-Dqx4XWuXu7w-unsplash.small.u2net_human_seg.png")
@test size(yannic_mask) == (800, 640)

tyler = unsplash_image("tyler-nix-k5R4cdIAiuk-unsplash.small.jpg")
@test size(tyler) == (960, 640)

end # module test_unsplashimages_unsplash_image


text:
function julia_toolsdir()
    prefix = dirname(dirname(Base.julia_cmd().exec[1]))
    toolsdir = joinpath(prefix, "tools")
    isfile(joinpath(toolsdir, "opt")) && return toolsdir
    return nothing
end

Base.@kwdef mutable struct ShowCodeConfig
    toolsdir::Union{String,Nothing} = julia_toolsdir()
    opt = nothing
    llc = nothing
    dot = nothing
    pygmentize = nothing
end

getcmd(name::Symbol) = getcmd(ShowCode.CONFIG, name)
function getcmd(config::ShowCodeConfig, name::Symbol)
    cmd = getfield(config, name)
    cmd === nothing || return cmd
    toolsdir = config.toolsdir
    if toolsdir !== nothing
        candidate = joinpath(toolsdir, string(name))
        isfile(candidate) && return `$candidate`
    end
    try
        return getfield(LLVM_jll, name)()::Base.AbstractCmd
    catch
    end
    try
        return getfield(LLVM_jll, name)(identity)::Base.AbstractCmd
    catch
    end
    prog = string(string(name))
    if Sys.which(prog) === nothing && config === ShowCode.CONFIG
        msg = (
            "Program `$name` not found. If it exists outside the `\$PATH, " *
            "set `ShowConfig.CONFIG.$name`."
        )
        if name in (:opt, :llc)
            msg *= (
                " You can also set `ShowConfig.CONFIG.toolsdir` to `usr/tools` of your" *
                "Julia's build directory (if any)."
            )
        end
        @error "$msg"
    end
    return `$prog`
end

abstract type AbstractCode end

function Base.show(io::IO, ::MIME"text/plain", code::AbstractCode)
    T = nameof(typeof(code))
    print(io, T, "(...)")
    return
end

InteractiveUtils.edit(code::AbstractCode) = InteractiveUtils.edit(abspath(code))

struct Fields{T}
    object::T
end

Base.propertynames(::Fields{T}) where {T} = fieldnames(T)
Base.getproperty(fields::Fields, name::Symbol) = getfield(getfield(fields, :object), name)

function write_silently(cmd, input; stdout = devnull)
    errio = IOBuffer()
    proc = open(pipeline(cmd, stderr = errio, stdout = stdout), write = true)
    try
        write(proc, input)
    finally
        close(proc)
    end
    wait(proc)
    if proc.exitcode != 0
        cmd0 = setenv(cmd)
        error(
            "Command $cmd0 (cwd: $(cmd.dir)) failed with code $(proc.exitcode) and error:\n",
            String(take!(errio)),
        )
    end
    return proc
end

function finalize_module()
    doc = read(joinpath(dirname(@__DIR__), "README.md"), String)
    @eval ShowCode $Base.@doc $doc ShowCode
end


text:
module EMLgibbs
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...reactors.EMLtank_basic
	import EMLtypes.length
	include("gibbs/gibbs_vap.jl")
	include("gibbs/gibbs_liq.jl")
end


text:
"""Betina Robust Sampled based Formulation"""
struct SampleBased <: AbstractSampleBased
    sampled_returns::Array{Float64,2}
    number_of_assets::Int
    number_of_samples::Int
end

function SampleBased(;
    sampled_returns::Array{Float64,2}
)
    number_of_samples, number_of_assets = size(sampled_returns)

    return SampleBased(
        sampled_returns, number_of_assets, number_of_samples
    )
end

function predicted_portfolio_return!(::JuMP.Model, w, formulation::AbstractSampleBased;
    predicted_mean = mean(formulation.sampled_returns, dims=1)'[:,1]
)
    return sum(w'predicted_mean)
end

function portfolio_return!(model::JuMP.Model, w, formulation::AbstractSampleBased; kwargs...)
    return predicted_portfolio_return!(model, w, formulation; kwargs...)
end

function predicted_portfolio_variance!(::JuMP.Model, w, formulation::AbstractSampleBased;
    predicted_covariance = cov(formulation.sampled_returns, dims=1)
)
    return sum(w'predicted_covariance * w)
end

function portfolio_variance!(model::JuMP.Model, w, formulation::AbstractSampleBased; kwargs...)
    return predicted_portfolio_variance!(model, w, formulation)
end

function po_max_predicted_return_limit_return(formulation::AbstractPortfolioFormulation, minimal_return;
    rf::Real = 0.0, current_wealth::Real = 1.0,
    model::JuMP.Model = base_model(formulation.number_of_assets; current_wealth=current_wealth),
    w = model[:w],
    kwargs...
)
    # auxilary variables
    @variable(model, R)
    sum_invested = create_sum_invested_variable(model, w)

    # model
    @constraint(model, R == portfolio_return!(model, w, formulation) + rf * (current_wealth - sum_invested))
    @constraint(model, R >= minimal_return * current_wealth)

    # objective function
    @objective(model, Max, predicted_portfolio_return!(model, w, formulation; kwargs...) + rf * (current_wealth - sum_invested))

    return model
end


text:
# Copyright (c) 2022, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.

include("../src/sudoku.jl") # Defines Sudoku module

# Puzzle example from
# https://sourceforge.net/projects/winsudoku/ accessed Jan 30, 2022
# Example taken from BSD-licensed project
puzzle_values = [
    0 2 1 0 6 0 0 0 8;
    8 0 0 0 0 0 7 5 0;
    0 0 0 0 0 1 0 0 0;
    0 0 4 9 0 3 0 0 0;
    9 0 0 1 0 0 0 3 0;
    1 0 0 0 0 8 0 0 0;
    0 0 0 0 0 0 5 0 0;
    0 7 0 0 0 9 0 0 6;
    0 8 3 7 0 0 1 0 4]
puzzle = Sudoku.SolvablePuzzle(3)
Sudoku.assign_values!(puzzle,puzzle_values)
println("Unsolved puzzle")
print(Sudoku.as_text_grid(puzzle))
# Backtrack solve using
#   * up to 8-fold compound rules to simplify backtracking
#   * up to 57 recursive steps (24 knowns)
#   * solution cap of 2 to detect a unique solution
results = Sudoku.backtrack_solve(puzzle,8,81-24,2)
if length(results) == 0
    println("Failure: puzzle solution is incomplete")
elseif length(results) == 1
    print("\nSolved puzzle\n")
    print(Sudoku.as_text_grid(results[1]))
    if Sudoku.valid_puzzle(Sudoku.as_values(results[1]))
        println("Verified that solution satisfies Suduko rules")
    else
        println("Solution is not a valid Sudoku")
    end
    if Sudoku.satisfies(Sudoku.as_values(puzzle),Sudoku.as_values(results[1]))
        println("Verified that solution matches original puzzle")
    else
        println("Failure: solution does not match the original puzzle")
    end
else
    println("Failure: puzzle solution has multiple solutions")
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
using JLD2, Plots, Printf, LinearAlgebra
include(pwd() * "/analysis_scripts/" * "post_analysis.jl")
searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
mesoscale_dir = pwd()
checkpoints = searchdir(mesoscale_dir, "iteration")
file = jldopen( mesoscale_dir * "/" * checkpoints[1])
b = file["tracers"]["b"]["data"][2:end-1, 2:end-1, 2:end-1]
u = file["velocities"]["u"]["data"][2:end-1, 2:end-1, 2:end-1]
v = file["velocities"]["v"]["data"][2:end-1, 2:end-1, 2:end-1]
w = file["velocities"]["w"]["data"][2:end-1, 2:end-1, 2:end-1]
v = (v[:, 1:end-1, :] + v[:, 2:end, :]) .* 0.5 # put everything on cell centers
w = (w[:, :, 1:end-1] + w[:,:, 2:end]) .* 0.5  # put everything on cell centers
close(file)
n = length(b)
avg_b = sum(b) / n
tke = @. u^2 + v^2 + w^2
avg_tke = sum(tke) / n
println("the average buoyancy is " * string(avg_b))
println("the average tke is " * string(avg_tke))
##
const Nx = size(b)[1]
const Ny = size(b)[2]
const Nz = size(b)[3]
const Lz = 3000.0
const Lx = Ly = 10^6 * 1.0
##
@inline function zC(k)
    return - Lz + (k-0.5) / Nz * Lz
end
@inline function yC(j)
    return  (j-0.5) / Ny * Ly
end
@inline function xC(j)
    return  (j-0.5) / Nx * Lx
end

x = xC.( collect(1:Nx) )
y = yC.( collect(1:Ny) )
z = zC.( collect(1:Nz) )

x2 = (x[1:end-1] + x[2:end]) / 2
y2 = (y[1:end-1] + y[2:end]) / 2
z2 = (z[1:end-1] + z[2:end]) / 2
##
∂x = PartialDerivative((x, 1))
∂y = PartialDerivative((y, 2))
∂z = PartialDerivative((z, 3))
∇ = [∂x, ∂y, ∂z]
##
field_label = [(u, "u"), (v, "v"), (w, "w"), (b, "b"), ( w .* w, "ww")]
selection = 4
field = field_label[selection][1]
label = field_label[selection][2]
cmax = maximum(field)
cmin = minimum(field)
clims = (cmin, cmax)
p1 = contourf(y, z, field[ 1, :, :]',
    color = :thermometer, title = "Zonal Slice " * label,
    xlabel = "Meridional [m]", ylabel = "Depth [m]"
    , clims = clims, linewidth = 0)
##
field_label = [(u, "u"), (v, "v"), (w, "w"), (b, "b"), ( w .* w, "ww")]
selection = 4
field = field_label[selection][1]
label = field_label[selection][2]
cmax = maximum(field)
cmin = minimum(field)
clims = (cmin, cmax)
ϕ = sum(field, dims = 1)[1,:,:] ./ Nx
p1 = contourf(y, z, ϕ',
    color = :thermometer, title = "Zonal Average " * label,
    xlabel = "Meridional [m]", ylabel = "Depth [m]"
    , clims = clims, linewidth = 0)

##
field_label = [(u, "u"), (v, "v"), (w, "w"), (b, "b"), (w .* w, "ww"), (u .* u, "uu"), (v .* v, "vv"), (∂z( w .* w), "d(w .* w) / dz")]
selection = 1
field = sum(field_label[selection][1], dims = (1,2)) ./ (Nx * Ny)
label = field_label[selection][2]
cmax = maximum(field)
cmin = minimum(field)
clims = (cmin, cmax)
p1 = scatter(field[1,1,:],  z,
    color = :blue, xlabel = "Field", label = "Horizontally Averaged " * label,
    ylabel = "Depth [m]", xlims = (cmin, cmax)
    , legend = :bottomright)

##
# day_label = @sprintf("%.2f ", sim_day[i])
# surface values
field_label = [(u, "u"), (v, "v"), (w, "w"), (b, "b"), (u .* b, "ub")]
selection = 4
field = field_label[selection][1]
label = field_label[selection][2]
cmax = maximum(field)
cmin = minimum(field)
clims = (cmin, cmax)
p1 = contourf(x, y, field[ :, :, end]',
    color = :thermometer, title = "Surface " * label,
    xlabel = "Zonal [m]", ylabel = "Meridional [m]"
    , clims = clims, linewidth = 0)

## Ertel potential vorticity
U  = [u, v, w]
ω  = ∇ × U
∇b = [∂x(b), ∂y(b), ∂z(b)]
f  = -1e-4
β  = 1e-11
zero_ϕ = ω[1] .* 0
Ω = [zero_ϕ, zero_ϕ, zero_ϕ .+ -1e-4 .+ β .* reshape(y2, (1,191,1))]
@. ω +=  Ω
vec_pv = [ω[1] .* ∇b[1] ,  ω[2] .* ∇b[2] , ω[3] .* ∇b[3]]
pv = vec_pv[1] + vec_pv[2] + vec_pv[3]
magnitude = zeros(size(ω[1]))
@inbounds @simd for i in 1:length(ω[1])
    ω1  = sqrt(ω[1][i]^2  + ω[2][i]^2  + ω[3][i]^2)
    ∇b1 = sqrt(∇b[1][i]^2 + ∇b[2][i]^2 + ∇b[3][i]^2)
    magnitude[i] = ω1 * ∇b1
end
alignment = pv ./ magnitude

##
field = pv
label = "Alignment"
ϕ = field[ :, :, end]
cmax = maximum(ϕ)
cmin = minimum(ϕ)
clims = (cmin, cmax)
p1 = contourf(x2, y2, ϕ',
    color = :thermometer, title = "Potential Vorticity ",
    xlabel = "Zonal [m]", ylabel = "Meridional [m]"
    , clims = clims, linewidth = 0)

##
sf_escale = 32
sf = @. exp( -(y - Ly/2)^2 / (Ly^2 / sf_escale) ) - exp( -(Ly - Ly/2)^2 / (Ly^2 / sf_escale) )
plot(y, sf, ylims = (0,1))

nw_escale = 2*10^4
nw = @. exp( (x - Ly) / nw_escale )
plot!(y, nw, ylims = (0,1))
##

bd_escale = 1 / 200
bd = @. exp(- bd_escale * ( z + Lz))
plot(z, bd, ylims = (0,1))

##
# Checking relaxation profile
Δb = 10 * 2e-3
h = 1000
Nx = 192
Nz = 32
@inline function zC(k)
    return - Lz + (k-0.5) / Nz * Lz
end

@inline function yC(j)
    return  (j-0.5) / Nx * Ly
end

relaxation_profile(j) = Δb * (yC(j)/ Ly)
relaxation_profile_north_2(k) = Δb * ( exp(zC(k)/h) - exp(-Lz/h) ) / (1 - exp(-Lz/h))

tmp1 = relaxation_profile.( collect(1:192))
tmp2 = relaxation_profile_north_2.( collect(1:32))

zonal_average_b = sum(b, dims = (1, )) ./ Nx
tmp1 = reshape(tmp1, (1, 192, 1))
plot( b[1, : , end])
plot!(tmp1)
new_field = zonal_average_b .- tmp1
plot(y, new_field[1, :, end])
##
plot( b[2, end , :], legend = false)
for i in 1:10:192
    plot!(b[i, end, :])
end

plot!(tmp2)

##
# Check Geostrophy
∇pʰ = hydrostatic_pressure(b,x,y,z)
zero_ϕ = zeros(size(∇pʰ[1]))
f  = -1e-4
β  = 1e-11
Ω = [zero_ϕ, zero_ϕ, zero_ϕ .+ f .+ β .* reshape(y2, (1,length(y2),1))]
u_avg = (avg_other(u,1)[1:end-1,:,:] + avg_other(u,1)[2:end,:,:]) ./ 2
v_avg = (avg_other(v,1)[1:end-1,:,:] + avg_other(v,1)[2:end,:,:]) ./ 2
w_avg = (avg_other(w,1)[1:end-1,:,:] + avg_other(w,1)[2:end,:,:]) ./ 2
U = [u_avg, v_avg, w_avg]
coriolis_force = [-Ω[3] .* v_avg,  Ω[3] .* u_avg, zero_ϕ]

norm(coriolis_force[1] .+ ∇pʰ[1]) / norm(coriolis_force[1])
norm(coriolis_force[2] .+ ∇pʰ[2]) / norm(coriolis_force[2])

##
field_label = [(coriolis_force[1], "-Ω v"), (-∇pʰ[1], "-∂x(pʰ)"), (coriolis_force[2], "Ω u"), (-∇pʰ[2], "-∂y(pʰ)")]
selection = 2+2
field = field_label[selection][1]
label = field_label[selection][2]
cmax = maximum(field)
cmin = minimum(field)
clims = (cmin, cmax)
p1 = contourf(x2, y2, field[ :, :, end]',
    color = :thermometer, title = "100 [m] depth " * label,
    xlabel = "Zonal [m]", ylabel = "Meridional [m]"
    , clims = clims, linewidth = 0)

selection = 1+2
field = field_label[selection][1]
label = field_label[selection][2]
p2 = contourf(x2, y2, field[ :, :, end]',
    color = :thermometer, title = "100 [m] depth " * label,
    xlabel = "Zonal [m]", yaxis = false
    , clims = clims, linewidth = 0)
plot(p1,p2)
##
gr(size=(700,400))
p3 = plot(p1,p2)
savefig(p3, pwd() * "/geostrophic_1.pdf")
##
p = ∫dz(b, z)
h_p = sum(p, dims = (1,2)) ./ ( (Nx -1) * (Ny-1))
plot(h_p[:], z2)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
@testset "Test integrated function Nx = 1" begin

Nx = 1
Ne = 50

atol = 1e-6

X = randn(Nx, Ne)

Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
for b in Blist
    B = MultiBasis(b, Nx)

    idx = reshape([0; 1; 2; 3; 4], (5, 1))
    truncidx = idx[1:2:end,:]
    Nψ = 5

    coeff = randn(Nψ)

    f = ExpandedFunction(B, idx, coeff)
    R = IntegratedFunction(f)

    # Test evaluate
    ψt = zeros(Ne)

    for i=1:Ne
        x = X[:,i]
        ψt[i] = R.f(vcat(x[1:end-1], 0.0)) +
                quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(x[1:end-1],t))[end]), 0, x[end])[1]
    end

    ψ = evaluate(R, X)

    @test norm(ψ - ψt)<atol

    # Test grad_x
    dψ = grad_x(R, X)

    function Gradient(x)
        return ForwardDiff.gradient(y->begin
                               y[end] = 0.0
                               return  R.f(y)
            end, x) + quadgk(t->ForwardDiff.gradient(
                      z->R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1],t))[end]),
                x), 0.0, x[end])[1] + vcat(zeros(size(x,1)-1), R.g(ForwardDiff.gradient(y->R.f(y), x)[end]))
    end

    for i=1:Ne
        xi = X[:,i]
        dψt = Gradient(xi)
        @test norm(dψ[i,:] - dψt)<1e-8
    end

    # Test hess_x
    H = hess_x(R, X)
    Ht = zeros(Nx, Nx)
    function Hessian!(H, x)
        fill!(H, 0.0)

        H[1:Nx-1,1:Nx-1] .= ForwardDiff.hessian(y->begin
                            y[Nx] = 0.0
                            return R.f(y) end, x)[1:Nx-1,1:Nx-1]

        # Hessian of the integral term
        H[1:Nx-1,1:Nx-1] .+= (quadgk(t -> ForwardDiff.hessian(
        z-> R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1], t))[end]), x),
                0.0, x[end])[1])[1:Nx-1,1:Nx-1]

        # H[Nx,:] and H[:,Nx]
        H[:, Nx] .= ForwardDiff.gradient(z->R.g(ForwardDiff.gradient(y->R.f(y), z)[end]), x)
        H[Nx, :] .= H[:, Nx]
        return H
    end

    @inbounds for i=1:Ne
        @test norm(H[i,:,:] - Hessian!(Ht, X[:,i]))<atol
    end

    # Test grad_xd
    gdψ =  grad_xd(R, X)

    gdψt = zeros(Ne)

    for i=1:Ne
        gdψt[i] = R.g(ForwardDiff.gradient(R.f, X[:,i])[end])
    end

    @test norm(gdψ - gdψt) <atol

    # Test grad_coeff_grad_xd
    dψ_xd_dc = grad_coeff_grad_xd(R, X)
    dψ_xd_dct = zeros(Ne, Nψ)

    ∂kf(y) = ForwardDiff.gradient(R.f, y)[end]

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for i=1:Ne
            dψ_xd_dct[i,j] = ∂kfj(X[:,i])*grad_x(R.g, ∂kf(X[:,i]))
        end
    end

    @test norm(dψ_xd_dc - dψ_xd_dct)<atol

    # Test hess_coeff_grad_xd
    d2ψ_xd_dc = hess_coeff_grad_xd(R, X)
    d2ψ_xd_dct = zeros(Ne, Nψ, Nψ)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for l=1:Nψ
        fl = MultiFunction(R.f.B, f.idx[l,:])
        ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                d2ψ_xd_dct[i,j,l] = 0.0*grad_x(R.g, ∂kf(X[:,i])) +
                    hess_x(R.g, ∂kf(X[:,i]))*∂kfj(X[:,i])*∂kfl(X[:,i])
            end
        end
    end

    @test norm(d2ψ_xd_dc - d2ψ_xd_dct)<atol

    # Test integrate_xd
    intψt = zeros(Ne)

    for i=1:Ne
        xi = X[:,i]
        intψt[i] = quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(xi[1:end-1],t))[end]), 0, xi[end])[1]
    end

    intψ = integrate_xd(R, X)

    @test norm(intψt - intψ)<atol

    # Test grad_coeff_integrate_xd

    dcintψt = zeros(Ne, Nψ)
    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for i=1:Ne
            xi .= X[:,i]
            dcintψt[i,j] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
        end
    end

    dcintψ = grad_coeff_integrate_xd(R, X)
    @test norm(dcintψ - dcintψt)<atol

    # Test hess_coeff_integrate_xd
    d2cintψt = zeros(Ne, Nψ, Nψ)
    xi = zeros(Nx)
    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for l=1:Nψ

            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                xi .= X[:,i]
    #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                    d2cintψt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end
    end

    d2cintψ = hess_coeff_integrate_xd(R, X)

    @test norm(d2cintψt - d2cintψ)<atol

    # Test grad_coeff
    dcRt = zeros(Ne, Nψ)

    X0 = deepcopy(X)
    X0[end,:] .= zeros(Ne)
    # ∂_c f(x_{1:k-1},0)
    dcRt += evaluate_basis(R.f, X0)
    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for i=1:Ne
            xi .= X[:,i]
            dcRt[i,j] += quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]

        end
    end

    dcR = grad_coeff(R, X)

    @test norm(dcRt - dcR)<atol

    # Test hess_coeff
    d2cRt = zeros(Ne, Nψ, Nψ)

    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for l=1:Nψ

            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                xi .= X[:,i]
    #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                    d2cRt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end
    end

    d2cR = hess_coeff(R, X)

    @test norm(d2cRt - d2cR)<atol

    end
end

@testset "Test integrated function Nx = 2" begin

Nx = 2
Ne = 50

atol = 1e-6

X = randn(Nx, Ne)

Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
for b in Blist
    B = MultiBasis(b, Nx)

    idx = [0 0; 0 1; 1 0; 1 1; 1 2]
    truncidx = idx[1:2:end,:]
    Nψ = 5

    coeff = randn(Nψ)

    f = ExpandedFunction(B, idx, coeff)
    R = IntegratedFunction(f)

    # Test evaluate
    ψt = zeros(Ne)

    for i=1:Ne
        x = X[:,i]
        ψt[i] = R.f(vcat(x[1:end-1], 0.0)) +
                quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(x[1:end-1],t))[end]), 0, x[end])[1]
    end

    ψ = evaluate(R, X)

    @test norm(ψ - ψt)<atol

    # Test grad_x
    dψ = grad_x(R, X)

    function Gradient(x)
        return ForwardDiff.gradient(y->begin
                               y[end] = 0.0
                               return  R.f(y)
            end, x) + quadgk(t->ForwardDiff.gradient(
                      z->R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1],t))[end]),
                x), 0.0, x[end])[1] + vcat(zeros(size(x,1)-1), R.g(ForwardDiff.gradient(y->R.f(y), x)[end]))
    end

    for i=1:Ne
        xi = X[:,i]
        dψt = Gradient(xi)
        @test norm(dψ[i,:] - dψt)<1e-8
    end

    # Test hess_x
    H = hess_x(R, X)
    Ht = zeros(Nx, Nx)
    function Hessian!(H, x)
        fill!(H, 0.0)

        H[1:Nx-1,1:Nx-1] .= ForwardDiff.hessian(y->begin
                            y[Nx] = 0.0
                            return R.f(y) end, x)[1:Nx-1,1:Nx-1]

        # Hessian of the integral term
        H[1:Nx-1,1:Nx-1] .+= (quadgk(t -> ForwardDiff.hessian(
        z-> R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1], t))[end]), x),
                0.0, x[end])[1])[1:Nx-1,1:Nx-1]

        # H[Nx,:] and H[:,Nx]
        H[:, Nx] .= ForwardDiff.gradient(z->R.g(ForwardDiff.gradient(y->R.f(y), z)[end]), x)
        H[Nx, :] .= H[:, Nx]
        return H
    end

    @inbounds for i=1:Ne
        @test norm(H[i,:,:] - Hessian!(Ht, X[:,i]))<atol
    end

    # Test grad_xd
    gdψ =  grad_xd(R, X)

    gdψt = zeros(Ne)

    for i=1:Ne
        gdψt[i] = R.g(ForwardDiff.gradient(R.f, X[:,i])[end])
    end

    @test norm(gdψ - gdψt) <atol

    # Test grad_coeff_grad_xd
    dψ_xd_dc = grad_coeff_grad_xd(R, X)
    dψ_xd_dct = zeros(Ne, Nψ)

    ∂kf(y) = ForwardDiff.gradient(R.f, y)[end]

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for i=1:Ne
            dψ_xd_dct[i,j] = ∂kfj(X[:,i])*grad_x(R.g, ∂kf(X[:,i]))
        end
    end

    @test norm(dψ_xd_dc - dψ_xd_dct)<atol

    # Test hess_coeff_grad_xd
    d2ψ_xd_dc = hess_coeff_grad_xd(R, X)
    d2ψ_xd_dct = zeros(Ne, Nψ, Nψ)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for l=1:Nψ
        fl = MultiFunction(R.f.B, f.idx[l,:])
        ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                d2ψ_xd_dct[i,j,l] = 0.0*grad_x(R.g, ∂kf(X[:,i])) +
                    hess_x(R.g, ∂kf(X[:,i]))*∂kfj(X[:,i])*∂kfl(X[:,i])
            end
        end
    end

    @test norm(d2ψ_xd_dc - d2ψ_xd_dct)<atol

    # Test integrate_xd
    intψt = zeros(Ne)

    for i=1:Ne
        xi = X[:,i]
        intψt[i] = quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(xi[1:end-1],t))[end]), 0, xi[end])[1]
    end

    intψ = integrate_xd(R, X)

    @test norm(intψt - intψ)<atol

    # Test grad_coeff_integrate_xd

    dcintψt = zeros(Ne, Nψ)
    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
        for i=1:Ne
            xi .= X[:,i]
            dcintψt[i,j] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
        end
    end

    dcintψ = grad_coeff_integrate_xd(R, X)
    @test norm(dcintψ - dcintψt)<atol

    # Test hess_coeff_integrate_xd
    d2cintψt = zeros(Ne, Nψ, Nψ)
    xi = zeros(Nx)
    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for l=1:Nψ

            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                xi .= X[:,i]
    #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                    d2cintψt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end
    end

    d2cintψ = hess_coeff_integrate_xd(R, X)

    @test norm(d2cintψt - d2cintψ)<atol

    # Test grad_coeff
    dcRt = zeros(Ne, Nψ)

    X0 = deepcopy(X)
    X0[end,:] .= zeros(Ne)
    # ∂_c f(x_{1:k-1},0)
    dcRt += evaluate_basis(R.f, X0)
    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for i=1:Ne
            xi .= X[:,i]
            dcRt[i,j] += quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]

        end
    end

    dcR = grad_coeff(R, X)

    @test norm(dcRt - dcR)<atol

    # Test hess_coeff
    d2cRt = zeros(Ne, Nψ, Nψ)

    xi = zeros(Nx)

    for j=1:Nψ
        fj = MultiFunction(R.f.B, f.idx[j,:])
        ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

        for l=1:Nψ

            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

            for i=1:Ne
                xi .= X[:,i]
    #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                    d2cRt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end
    end

    d2cR = hess_coeff(R, X)

    @test norm(d2cRt - d2cR)<atol

    end
end

@testset "Test integrated function Nx = 3" begin

    Nx = 3
    Ne = 50

    X = randn(Nx, Ne)

    atol = 1e-6

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        # idx = [0 0; 0 1; 1 0; 1 1; 1 2]
        idx = [0 0 0 ;0  0 1; 0 1 0; 1 0 0;1 1 0; 0 1 1; 1 0 1; 1 1 1; 1 2 0; 2 1 0; 2 2 2]

        truncidx = idx[1:2:end,:]
        Nψ = 11

        coeff = randn(Nψ)

        f = ExpandedFunction(B, idx, coeff)
        R = IntegratedFunction(f)

        # Test grad_x
        dψ = grad_x(R, X)

        function Gradient(x)
            return ForwardDiff.gradient(y->begin
                                   y[end] = 0.0
                                   return  R.f(y)
                end, x) + quadgk(t->ForwardDiff.gradient(
                          z->R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1],t))[end]),
                    x), 0.0, x[end])[1] + vcat(zeros(size(x,1)-1), R.g(ForwardDiff.gradient(y->R.f(y), x)[end]))
        end

        for i=1:Ne
            xi = X[:,i]
            dψt = Gradient(xi)
            @test norm(dψ[i,:] - dψt)<atol
        end

        # Test grad_xd
        gdψ =  grad_xd(R, X)

        gdψt = zeros(Ne)

        for i=1:Ne
            gdψt[i] = R.g(ForwardDiff.gradient(R.f, X[:,i])[end])
        end

        @test norm(gdψ - gdψt) <atol

        # Test grad_coeff_grad_xd
        dψ_xd_dc = grad_coeff_grad_xd(R, X)
        dψ_xd_dct = zeros(Ne, Nψ)

        ∂kf(y) = ForwardDiff.gradient(R.f, y)[end]

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                dψ_xd_dct[i,j] = ∂kfj(X[:,i])*grad_x(R.g, ∂kf(X[:,i]))
            end
        end

        @test norm(dψ_xd_dc - dψ_xd_dct)<atol

        # Test hess_coeff_grad_xd
        d2ψ_xd_dc = hess_coeff_grad_xd(R, X)
        d2ψ_xd_dct = zeros(Ne, Nψ, Nψ)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for l=1:Nψ
            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

                for i=1:Ne
                    d2ψ_xd_dct[i,j,l] = 0.0*grad_x(R.g, ∂kf(X[:,i])) +
                        hess_x(R.g, ∂kf(X[:,i]))*∂kfj(X[:,i])*∂kfl(X[:,i])
                end
            end
        end

        @test norm(d2ψ_xd_dc - d2ψ_xd_dct)<atol

        # Test integrate_xd
        intψt = zeros(Ne)

        for i=1:Ne
            xi = X[:,i]
            intψt[i] = quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(xi[1:end-1],t))[end]), 0, xi[end])[1]
        end

        intψ = integrate_xd(R, X)

        @test norm(intψt - intψ)<atol

        # Test grad_coeff_integrate_xd

        dcintψt = zeros(Ne, Nψ)
        xi = zeros(Nx)
        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                xi .= X[:,i]
                dcintψt[i,j] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end

        dcintψ = grad_coeff_integrate_xd(R, X)
        @test norm(dcintψ - dcintψt)<atol

        # Test hess_coeff_integrate_xd
        d2cintψt = zeros(Ne, Nψ, Nψ)
        xi = zeros(Nx)
        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

            for l=1:Nψ

                fl = MultiFunction(R.f.B, f.idx[l,:])
                ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

                for i=1:Ne
                    xi .= X[:,i]
        #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                        d2cintψt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                end
            end
        end

        d2cintψ = hess_coeff_integrate_xd(R, X)

        @test norm(d2cintψt - d2cintψ)<atol

        # Test grad_coeff
        dcRt = zeros(Ne, Nψ)

        X0 = deepcopy(X)
        X0[end,:] .= zeros(Ne)
        # ∂_c f(x_{1:k-1},0)
        dcRt += evaluate_basis(R.f, X0)
        xi = zeros(Nx)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                xi .= X[:,i]
                dcRt[i,j] += quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]

            end
        end

        dcR = grad_coeff(R, X)

        @test norm(dcRt - dcR)<atol

        # Test hess_coeff
        d2cRt = zeros(Ne, Nψ, Nψ)

        xi = zeros(Nx)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

            for l=1:Nψ

                fl = MultiFunction(R.f.B, f.idx[l,:])
                ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

                for i=1:Ne
                    xi .= X[:,i]
        #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                    d2cRt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                end
            end
        end

        d2cR = hess_coeff(R, X)

        @test norm(d2cRt - d2cR)<atol
    end
end

@testset "Test integrated function Nx = 4" begin

    Nx = 4
    Ne = 50

    X = randn(Nx, Ne)

    atol = 5e-6

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        # idx = [0 0; 0 1; 1 0; 1 1; 1 2]
        idx = [0 0 0 0; 0 2 0 1; 0 2 3 0; 2 0 2 1; 0 0 1 2; 1 2 0 2;3 2 2 2]

        truncidx = idx[1:2:end,:]
        Nψ = 7

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)
        R = IntegratedFunction(f)

        # Test grad_x
        dψ = grad_x(R, X)

        function Gradient(x)
            return ForwardDiff.gradient(y->begin
                                   y[end] = 0.0
                                   return  R.f(y)
                end, x) + quadgk(t->ForwardDiff.gradient(
                          z->R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1],t))[end]),
                    x), 0.0, x[end])[1] + vcat(zeros(size(x,1)-1), R.g(ForwardDiff.gradient(y->R.f(y), x)[end]))
        end

        for i=1:Ne
            xi = X[:,i]
            dψt = Gradient(xi)
            @test norm(dψ[i,:] - dψt)<atol
        end

        # Test hess_x
        H = hess_x(R, X)
        Ht = zeros(Nx, Nx)
        function Hessian!(H, x)
            fill!(H, 0.0)

            H[1:Nx-1,1:Nx-1] .= ForwardDiff.hessian(y->begin
                                y[Nx] = 0.0
                                return R.f(y) end, x)[1:Nx-1,1:Nx-1]

            # Hessian of the integral term
            H[1:Nx-1,1:Nx-1] .+= (quadgk(t -> ForwardDiff.hessian(
            z-> R.g(ForwardDiff.gradient(y->R.f(y), vcat(z[1:end-1], t))[end]), x),
                    0.0, x[end])[1])[1:Nx-1,1:Nx-1]

            # H[Nx,:] and H[:,Nx]
            H[:, Nx] .= ForwardDiff.gradient(z->R.g(ForwardDiff.gradient(y->R.f(y), z)[end]), x)
            H[Nx, :] .= H[:, Nx]
            return H
        end

        @inbounds for i=1:Ne
            @test norm(H[i,:,:] - Hessian!(Ht, X[:,i]))<atol
        end

        # Test grad_xd
        gdψ =  grad_xd(R, X)

        gdψt = zeros(Ne)

        for i=1:Ne
            gdψt[i] = R.g(ForwardDiff.gradient(R.f, X[:,i])[end])
        end

        @test norm(gdψ - gdψt) <atol

        # Test grad_coeff_grad_xd
        dψ_xd_dc = grad_coeff_grad_xd(R, X)
        dψ_xd_dct = zeros(Ne, Nψ)

        ∂kf(y) = ForwardDiff.gradient(R.f, y)[end]

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                dψ_xd_dct[i,j] = ∂kfj(X[:,i])*grad_x(R.g, ∂kf(X[:,i]))
            end
        end

        @test norm(dψ_xd_dc - dψ_xd_dct)<atol

        # Test hess_coeff_grad_xd
        d2ψ_xd_dc = hess_coeff_grad_xd(R, X)
        d2ψ_xd_dct = zeros(Ne, Nψ, Nψ)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for l=1:Nψ
            fl = MultiFunction(R.f.B, f.idx[l,:])
            ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

                for i=1:Ne
                    d2ψ_xd_dct[i,j,l] = 0.0*grad_x(R.g, ∂kf(X[:,i])) +
                        hess_x(R.g, ∂kf(X[:,i]))*∂kfj(X[:,i])*∂kfl(X[:,i])
                end
            end
        end

        @test norm(d2ψ_xd_dc - d2ψ_xd_dct)<atol

        # Test integrate_xd
        intψt = zeros(Ne)

        for i=1:Ne
            xi = X[:,i]
            intψt[i] = quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(xi[1:end-1],t))[end]), 0, xi[end])[1]
        end

        intψ = integrate_xd(R, X)

        @test norm(intψt - intψ)<atol

        # Test grad_coeff_integrate_xd

        dcintψt = zeros(Ne, Nψ)
        xi = zeros(Nx)
        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                xi .= X[:,i]
                dcintψt[i,j] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
            end
        end

        dcintψ = grad_coeff_integrate_xd(R, X)
        @test norm(dcintψ - dcintψt)<atol

        # Test hess_coeff_integrate_xd
        d2cintψt = zeros(Ne, Nψ, Nψ)
        xi = zeros(Nx)
        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for l=1:Nψ
                fl = MultiFunction(R.f.B, f.idx[l,:])
                ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]
                for i=1:Ne
                    xi .= X[:,i]
                    d2cintψt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                end
            end
        end

        d2cintψ = hess_coeff_integrate_xd(R, X)

        @test norm(d2cintψt - d2cintψ)<atol

        # Test grad_coeff
        dcRt = zeros(Ne, Nψ)

        X0 = deepcopy(X)
        X0[end,:] .= zeros(Ne)
        # ∂_c f(x_{1:k-1},0)
        dcRt += evaluate_basis(R.f, X0)
        xi = zeros(Nx)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]
            for i=1:Ne
                xi .= X[:,i]
                dcRt[i,j] += quadgk(t->∂kfj(vcat(xi[1:end-1],t))*grad_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]

            end
        end

        dcR = grad_coeff(R, X)

        @test norm(dcRt - dcR)<atol

        # Test hess_coeff
        d2cRt = zeros(Ne, Nψ, Nψ)

        xi = zeros(Nx)

        for j=1:Nψ
            fj = MultiFunction(R.f.B, f.idx[j,:])
            ∂kfj(y) = ForwardDiff.gradient(fj, y)[end]

            for l=1:Nψ

                fl = MultiFunction(R.f.B, f.idx[l,:])
                ∂kfl(y) = ForwardDiff.gradient(fl, y)[end]

                for i=1:Ne
                    xi .= X[:,i]
        #                 @show quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                        d2cRt[i,j,l] = quadgk(t->∂kfj(vcat(xi[1:end-1],t))*∂kfl(vcat(xi[1:end-1],t))*hess_x(R.g, ∂kf(vcat(xi[1:end-1],t))),0,xi[end])[1]
                end
            end
        end

        d2cR = hess_coeff(R, X)

        @test norm(d2cRt - d2cR)<atol
    end
end


text:
### Define the target distribution and its gradient
using Distributions: logpdf, MvNormal
using DiffResults: GradientResult, value, gradient
using ForwardDiff: gradient!

D = 10
target = MvNormal(zeros(D), ones(D))
ℓπ(θ) = logpdf(target, θ)

function ∂ℓπ∂θ(θ)
    res = GradientResult(θ)
    gradient!(res, ℓπ, θ)
    return (value(res), gradient(res))
end

### Build up a HMC sampler to draw samples
using AdvancedHMC

# Sampling parameter settings
n_samples, n_adapts = 10_000, 2_000

# Draw a random starting points
θ_init = randn(D)

# Define metric space, Hamiltonian, sampling method and adaptor
metric = DiagEuclideanMetric(D)
h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
int = Leapfrog(find_good_eps(h, θ_init))
prop = NUTS{MultinomialTS,GeneralisedNoUTurn}(int)
adaptor = StanHMCAdaptor(n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, int))

# Draw samples via simulating Hamiltonian dynamics
# - `samples` will store the samples
# - `stats` will store statistics for each sample
samples, stats = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; progress=true)


text:
using JLBoost, RDatasets
iris = dataset("datasets", "iris")

iris[!, :is_setosa] = iris[!, :Species] .== "setosa"
target = :is_setosa

features = setdiff(names(iris), ["Species", "is_setosa"])

# fit one tree
# ?jlboost for more details
xgtreemodel = jlboost(iris, target)

typeof(trees(xgtreemodel))

typeof(xgtreemodel.loss)

typeof(xgtreemodel.target)

xgtreemodel2 = jlboost(iris, target; nrounds = 2, max_depth = 2)

xgtreemodel3 = jlboost(iris, target; nrounds = 1, max_leaves = 4, max_depth = 0)

iris.pred1 = predict(xgtreemodel, iris)
iris.pred2 = predict(xgtreemodel2, iris)
iris.pred1_plus_2 = predict(vcat(xgtreemodel, xgtreemodel2), iris)

AUC(-iris.pred1, iris.is_setosa)

gini(-iris.pred1, iris.is_setosa)

feature_importance(xgtreemodel, iris)

nrow(df) # returns the number of rows
ncol(df)
view(df, rows, cols)

using DataFrames
using JLBoost
df = DataFrame(x = rand(100) * 100)

df[!, :y] = 2*df.x .+ rand(100)

target = :y
features = [:x]
warm_start = fill(0.0, nrow(df))

using LossFunctions: L2DistLoss
loss = L2DistLoss()
jlboost(df, target, features, warm_start, loss; max_depth=2) # default max_depth = 6

JLBoost.save(xgtreemodel, "model.jlb")
JLBoost.save(trees(xgtreemodel), "model_tree.jlb")

JLBoost.load("model.jlb")
JLBoost.load("model_tree.jlb")

using JLBoost, RDatasets, JDF
iris = dataset("datasets", "iris")

iris[!, :is_setosa] = iris[!, :Species] .== "setosa"
target = :is_setosa

features = setdiff(names(iris), [:Species, :is_setosa])

savejdf("iris.jdf", iris)
irisdisk = JDFFile("iris.jdf")

# fit using on disk JDF format
xgtree1 = jlboost(irisdisk, target, features)
xgtree2 = jlboost(iris, target, features; nrounds = 2, max_depth = 2)

# predict using on disk JDF format
iris.pred1 = predict(xgtree1, irisdisk)
iris.pred2 = predict(xgtree2, irisdisk)

# AUC
AUC(-predict(xgtree1, irisdisk), irisdisk[!, :is_setosa])

# gini
gini(-predict(xgtree1, irisdisk), irisdisk[!, :is_setosa])

# clean up
rm("iris.jdf", force=true, recursive=true)


text:
# ---
# title: 115. Distinct Subsequences
# id: problem115
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/distinct-subsequences/description/>
# hidden: true
# ---
#
# Given two strings `s` and `t`, return _the  number of distinct subsequences of
# `s` which equals `t`_.
#
# A string's **subsequence** is a new string formed from the original string by
# deleting some (can be none) of the characters without disturbing the relative
# positions of the remaining characters. (i.e., `"ACE"` is a subsequence of
# `"ABCDE"` while `"AEC"` is not).
#
# It's guaranteed the answer fits on a 32-bit signed integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "rabbbit", t = "rabbit"
#     Output: 3
#     Explanation:
#     As shown below, there are 3 ways you can generate "rabbit" from S.
#     **_rabb_** b ** _it_**
#     **_ra_** b ** _bbit_**
#     **_rab_** b ** _bit_**
#
#
# **Example 2:**
#
#
#
#     Input: s = "babgbag", t = "bag"
#     Output: 5
#     Explanation:
#     As shown below, there are 5 ways you can generate "bag" from S.
#     **_ba_** b _ **g**_ bag
#     **_ba_** bgba ** _g_**
#     _**b**_ abgb ** _ag_**
#     ba _ **b**_ gb _ **ag**_
#     babg ** _bag_**
#
#
#
# **Constraints:**
#
#   * `0 <= s.length, t.length <= 1000`
#   * `s` and `t` consist of English letters.
#
#
## @lc code=start
using LeetCode

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


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 88629870-3490-11eb-2215-d765d59a440c
using PlutoUI

# ╔═╡ 1be7c0c0-d4ba-11ea-1e43-ef66c35b0830
using Plots

# ╔═╡ 159f28e2-d4db-11ea-13e6-a97e2e7f7c60
#using DifferentialEquations
using OrdinaryDiffEq

# ╔═╡ 6cbb32f0-d665-11ea-3f24-41cbe1910b80
using LinearAlgebra

# ╔═╡ 6ba02e80-d4cf-11ea-19cb-33028d3dc67d
md"# _Primer tema_
+ ¿Qué es una ecuación diferencial y como se resuelve?
    * Euler
    * Métodos mejores
+ Sistema de segundo orden con dos polos
+ Explosión en tiempo finito
+ Linealización
+ Ejemplos
    * Generador
+ Ciclo límite
    * Van der Pool
+ Sistemas discretos
    * Ecuación logística, el camino hacia el Caos"

# ╔═╡ 7db3fdb0-d4d2-11ea-21fc-597c55dda572
md"# ¿Qué es una ecuación diferencial?

En la asignatura nos interesan sistemas como este:

$$\frac{dx}{dt}=f(x,t)$$

$$x(0)=x_0$$

Donde x es un vector y la derivada respecto al tiempo se suele escribir como $$\dot x$$

Las ecuaciones diferenciales son útiles cuando es fácil expresar cómo cambia una magnitud en el tiempo, pero no tanto cual es el estado $x(t)$ directamente, que es lo que ocurre en casi todos los sistemas físicos.
"

# ╔═╡ 76079a70-d4cf-11ea-2645-8f43c78eae52
md" ## Euler
Si la función $$f$$ tiene buen comportamiento (lo veremos más adelante cuando repasemos la existencia y unicidad) hay una solución de la ecuación que es continua y derivable y usando taylor:

$$x(\Delta t)=x(0) + f(x(0),0)\Delta t + o(\Delta t^2)$$

Si se toman pasos pequeños y repetimos el procedimiento:

$$x(2\Delta t)=x(\Delta t) + f(x(\Delta t),\Delta t) \Delta t + o(\Delta t^2)$$
...

y llamando $$t_n$$ a $$n\Delta t$$ y $$x_n$$ a $$x(t_n)$$

$$x_{n+1} \approx x_n + f(x_n,t_n)\Delta t + o(\Delta t^2)$$

### Ejemplo de implementación sencilla
Para el caso escalar

"

# ╔═╡ 66e9ffd0-d4d7-11ea-09c0-03ebaacab879
function euler(f,x0,tf,Δt)
	t=0:Δt:tf
	x=zeros(size(t))
	x[1]=x0
	for n=1:(length(t)-1)
		x[n+1]=x[n]+f(x[n],t[n])*Δt
	end
	return (x,t)
end

# ╔═╡ 0b94f1c0-7113-11eb-13da-c5c5b95ebbf5
md"Vamos a ver si funciona resolviendo
$\dot x=-x+sin(t)$ "

# ╔═╡ e885f5b0-d4d4-11ea-06d6-9ddf9619961f
begin
	f(x,t)=-x+sin(t)

	f(x,p,t)=f(x,t) #¡Despacho múltiple! (lo necesitamos para luego)
end

# ╔═╡ b99fc9b0-d4d9-11ea-0e18-cfd5af90188f
tf=10

# ╔═╡ 3be6b020-d4d8-11ea-3356-23afd148ebc9
@bind N Slider(1:100)

# ╔═╡ 220f6f12-d4d9-11ea-159e-c1e1551fd732
Δt=tf/N;

# ╔═╡ 23f5ece0-d4da-11ea-1330-c5d077bd3a89
md"N=$(N) $$\to \Delta t=$$ $(Δt)"

# ╔═╡ d8b5ed40-d4d7-11ea-27d4-bf594518802b
x,t =euler(f,0,10,Δt);

# ╔═╡ 2b7df180-d4d8-11ea-118f-75e8709e36c3
plot(t,x,label="x")

# ╔═╡ 455e0c30-d4d7-11ea-22ce-0d04b4e8d860
md"## Hay métodos mejores:

+ Ajustan los pasos automáticamente

+ Tienen más precisión para los mismos pasos

+ Están integrados con otras funciones como *plots*, permite interpolar los resultados automáticamente, etc...

Vamos a usar la librería OrdinaryDiffEq que es parte de DifferencialEquation.jl
"

# ╔═╡ c90980d0-d637-11ea-112d-c336b2080d87
begin
	x0 = 0
	tspan = (0.0,10.0)
	#ojo se necesita f(x,p,t) donde p son los parámetros (no tenemos)
	#de ahí la definición de antes
	prob = ODEProblem(f,x0,tspan);
	sol = solve(prob,Tsit5(), reltol=1e-8, abstol=1e-8);
	sol.u #su u es nuestro x
end

# ╔═╡ 1f44a88e-7114-11eb-0ac2-d12cd3fb40ef
sol.t

# ╔═╡ 4a7f5f72-0e0c-11eb-2b7c-4d7214c09c4b
sol(2.751) #¡se puede usar como una función continua de t!

# ╔═╡ 51059990-0e0c-11eb-2d55-89ac739619c7

# ╔═╡ a92de300-d63c-11ea-2446-5b099592d734
plot(sol, xaxis="t",yaxis="x(t)",label="x",legend=false)

# ╔═╡ 09980bd0-d63d-11ea-0752-1d033ee00cc9
md"## Sistema de segundo orden con dos polos

En un sistema lineal

$\dot x=Ax$

La solución es $x=x_0e^{At}$ y su carácter depende de los autovalores, es decir las soluciones de $|A- \lambda \mathbb{1} |=0$

La estabilidad depende de que la parte real de los $\lambda$ sea negativa (los autovalores estén en el semiplano izquierdo del plano complejo)

Nos interesa el sistema lineal segundo orden porque es simple y contiene todos los casos de interés:

$\dot x_1=x_2$
$\dot x_2=p_1x_1+p_2x_2$
"

# ╔═╡ e6afc70e-d63d-11ea-0d36-738b13f0de30
function f_lineal(x,p,t)
	dx=[0.0 ; 0.0]
	dx[1]=x[2];
	dx[2]=p[1]*x[1]+p[2]*x[2]
	return dx
end

# ╔═╡ 2ff44950-7116-11eb-3ad9-c3e145cdae4e
md"
Polos Reales? $(@bind son_reales CheckBox())
$(@bind slider1 Slider(-1:0.05:1))
$(@bind slider2 Slider(-1:0.05:1))
"

# ╔═╡ a0573fe0-d6a2-11ea-057e-5140c323c005
begin
	if son_reales==true
		# el polinomio es (s-real1)*(s-real2)
		# s^2+(-real1-real2)s + real1*real2 --> k2= real1 + real2 k1=real1*real2
		k2= slider1 + slider2
		k1=-slider1*slider2
	else # el polinomio es (s-real-imag*im)*(s-real+imag*im)
		 # es decir s^2 - 2*real*s + real^2 + imag^2 --> k2=-2*real k1=-real^2-imag^2
		k2=2*slider1
		k1=-slider1^2 - slider2^2
	end
	A=[0 1; k1 k2]
	autovalores,_=eigen(A)  #para esto necesito LinearAlgebra
	md"""
	k1= $(k1), k2= $(k2)

	autovalores= $(autovalores)"""
end

# ╔═╡ 5a802550-d665-11ea-3f9f-b103ce80822f
scatter(real(autovalores),imag(autovalores),
	grid=true,legend=false, framestyle = :origin,
	xaxis=[-1,1], yaxis=[-1,1] )

# ╔═╡ 1790c760-0e0d-11eb-1c47-f5e902e158fb
md" # Muy bien pero ¡Enséñame el código!
... vosotros lo habéis querido ☺	"

# ╔═╡ f5ca89e0-0e0c-11eb-32a1-fb22511f42ae
function flechas!(figura,f,p,rangos;N=10)
	#modifica una figura añadiendo las flechas por eso tiene una exclamación !
	# N es un parámetro opcional, si no se lo paso el valor es 10
	# los otros son posicionales
	xmin,xmax,ymin,ymax=rangos #desestructurar un iterable (explicar)
	longitud=max(abs(xmax-xmin),abs(ymax-ymin))/N

    #hago una maya de x e ys y pongo las flechas en cada punto
	x1s=xmin:((xmax-xmin)/N):xmax
	x2s=ymin:((ymax-ymin)/N):ymax

    #barrido para las flechas
	X1=[];
	X2=[];
	U=[];
	V=[];
	for x1 in x1s #barrido para calcular las flechas
		for x2 in x2s
			push!(X1,x1); #Push añade un componente al final de un vector
			push!(X2,x2);
			F=f([x1,x2],p,0)
			X=0.1*F/(0.1+sqrt(F[1]^2+F[2]^2)) #Saturo el tamaño de la flecha
			push!(U,X[1]);
			push!(V,X[2]);
		end
	end
	#Pinta las flechas
	quiver!(figura,X1,X2,quiver=(U,V), arrow=arrow(0.1*longitud, 0.1*longitud))
	return figura
end

# ╔═╡ 165464e2-0e0f-11eb-3851-3701a6c78502
#= Es un euler modificado que integra hacia delante y hacia atrás
    congtinúa integrando hasta que se sale o completa un número de pasos.
	   -hacia delante se detectan bien los PE y ciclos límite estables
	   -Hacia atrás permite ver también los inestables
	Ajusta dt para que el paso tenga un tamaño que sea aproximadamente constante=#

function trayectoria(x0,f,p,rangos;N=10,pasosMax=100)
	xmin,xmax,ymin,ymax=rangos
	long=max(abs(xmax-xmin),abs(ymax-ymin))/N
	x,y=x0
	pasos=0

	X1=[x]
	Y1=[y]
	X2=[x] #parte de la trayectoria hacia atrás en el tiempo
	Y2=[y]

	while (xmin<x<xmax) & (ymin<y<ymax) & (pasos<pasosMax)
		derivada=f([x,y],p,0)
		dt=0.05*long/(0.1+sqrt(derivada[1]^2+derivada[2]^2))
		Xn=[x;y]+derivada*dt
		x,y=Xn
		pasos=pasos+1
		push!(X1,x)
		push!(Y1,y)
	end
	pasos=0
	x,y=x0
	while (xmin<x<xmax) & (ymin<y<ymax) & (pasos<pasosMax)
		derivada=f([x,y],p,0)
		dt=-0.05*long/(0.1+sqrt(derivada[1]^2+derivada[2]^2))
		Xn=[x;y]+derivada*dt
		x,y=Xn
		pasos=pasos+1
		push!(X2,x)
		push!(Y2,y)
	end
	X=append!(reverse(X2),X1) #junto ambas partes
	Y=append!(reverse(Y2),Y1)
	return (X,Y)
end

# ╔═╡ b51d4d90-d7f4-11ea-07aa-adcb4f418c99
md"## Clasificación sistema de segundo orden"

# ╔═╡ 0b2775be-d7fb-11ea-0b82-0d2b4d544d31
md"### Estables autovalores en el semiplano izquierdo"

# ╔═╡ 4100c3f0-d7f5-11ea-10c3-593316528d50
md"Nodo **estable** (polos reales en -1 y -0.5)"

# ╔═╡ f8cc9600-d7f8-11ea-138d-5f40babd5dce
md"Foco **estable** (polos reales en $$-1\pm i$$)"

# ╔═╡ 34bbafa0-d7fb-11ea-2ed3-8b300f45fbac
md"### En el limite de la estabilidad (autovalores imaginarios puros)"

# ╔═╡ 345ec290-d7fb-11ea-1f16-d919242fefc3
md"Centro, polos en $$\pm i$$"

# ╔═╡ 8c6fa2e0-dd75-11ea-163b-c70c276c3f2a
md"### Caso *degenerado*, autovalor nulo"

# ╔═╡ ad48c140-dd75-11ea-03fe-75a0b19f91bd
md"""No hay un punto de equilibrio sino una "*línea de equilibrio*"
autovalores -1 y 0"""

# ╔═╡ 33f181d0-d7fb-11ea-0a51-d3a11d1613df
md"### Inestables, algún autovalor en el semiplano derecho"

# ╔═╡ fa4a8370-d7f8-11ea-12b0-fde507f287be
md"Nodo **inestable** (polos reales en +1 y +0.5)"

# ╔═╡ 08c0cb70-d7fa-11ea-0efd-ab3b7f76a8b8
md"Foco **inestable** (polos reales en $$+1\pm i$$)"

# ╔═╡ b6db1b70-d7fa-11ea-2bea-43295b3179a6
md"punto de silla (**Inestable**) polos en $$\pm1$$"

# ╔═╡ fd823200-d7fc-11ea-1af8-c91386a02ab5
md"""## Escape en tiempo finito
Consideremos la ecuación aparentemente *inocente*:

$\dot x=x^2$

$x_0=1$

Es separable y como se vio en teoría la solución es obviamente (basta con derivar y sustituir)

$x(t)=\frac{1}{1-t}$

Que en t=1 se hace infinita, veámoslo

"""

# ╔═╡ 2647d140-d7fd-11ea-201a-5744621fbe12
escape(x,p,t)=x^2

# ╔═╡ 7c13c6b0-d7fd-11ea-0cc5-55f5ab6106fe
begin
	prob2 = ODEProblem(escape,1,(0,2));
	sol2 = solve(prob2,Tsit5(), reltol=1e-5, abstol=1e-5);
	fig=plot(sol2.t,min.(sol2.u,10000),xlim=(0,1))
	#saturo a 10000 para evitar que "pete"
	ylims!(fig,(0.0,10))
end

# ╔═╡ 82023190-dd77-11ea-12d1-cfa1a2496348
md"# Linealización"

# ╔═╡ 16d7ee1e-dd7a-11ea-13a7-197edb480651
md""" Ejemplo del generador de corriente

$$\dot x_1 = x_2$$

$$\dot x_2 = \frac{P_m - P_e sin(x_1)}{H}$$  """

# ╔═╡ ab5bd370-dd77-11ea-034f-6153487b9658
function generador(x,p,t)
	x1,x2=x #así es más bonito, es la forma RECOMENDADA para la legibilidad
	H,Pm,Pe=p #Pm es u de la teoría
    dx1=x2
    dx2=(1/H)*(Pm-Pe*sin(x1))
	#el return es opcional (si no se pone devuelve la última línea
	#pero me gusta ser explícito
	return [dx1; dx2]
end

# ╔═╡ 1668b190-dd7a-11ea-3c80-c5f082ec8647
md"""## Puntos de Equilibrio
Los puntos *interesantes* es donde las flechas se **desvanecen**, es decir donde la derivada se anula:

$$0= x_2$$

$$0 = \frac{P_m - P_e sin(x_1)}{H} \to sin(x_1)=\frac{P_m}{P_e}$$

Las soluciones existen si ${|P_m| \le |P_e|}$

soluciones tipo 1: $$(arcsin\left(\frac{P_m}{P_e}\right)+2n\pi,0)$$

soluciones tipo 2: $$(\pi-arcsin\left(\frac{P_m}{P_e}\right)+2n\pi,0)$$

"""

# ╔═╡ 351f5332-dd7c-11ea-2fc4-c9dc71684f23
@bind sliderMotor Slider(-1.1:0.05:1.1)

# ╔═╡ 5fd238f0-dd7b-11ea-2a00-a587ee59c261
begin
	plot(-2*pi:0.01:2*pi,sin,legend=false)
	plot!([-2*pi,2*pi],[sliderMotor,sliderMotor ])
	title!("Pm/Pe=$sliderMotor")
end

# ╔═╡ db58def0-dd78-11ea-335b-4140d22d9a78
md""" En el punto de equilibrio $f(x^*)=0$ así que usando taylor:

definiendo $\delta x=x-x^*$

$\delta \dot x= \dot x=f(x^*+\delta x) \approx f(x^*)+\nabla f(x^*)\delta x=A\delta x$

donde A es la matriz jacobiana

$$A=\begin{pmatrix}
\frac{df_1}{dx_1}  & ... & \frac{df1}{dx_n} \\
\vdots             & \ddots   &  \vdots \\
\frac{df_n}{dx_1}  & ... &\frac{df_n}{dx_n}
\end{pmatrix}$$

que en el caso del motor es:

$$A=\begin{pmatrix}
0                     &  1 \\
\frac{P_ecos(x^*)}{H}  &  0
\end{pmatrix}$$

"""

# ╔═╡ 33af7e82-d7f5-11ea-12bb-5982575f35d3
md"""## Ciclo límite
Ecuación de Van der Pol

$\dot x_1=x_2$

$\dot x_2=\frac{-kx_1+c(1-x_1^2)x_2}{m}$

"""

# ╔═╡ 3d70984e-d7f5-11ea-3243-9d3932ee155b
function van_der_pol(x,p,t)
    m,c,k=p
	x1,x2=x
    dx1=x2
    dx2=(1/m)*(-k*x1+c*(1-x1^2)*x2)
	return [dx1;dx2]
end

# ╔═╡ 713de9b0-3490-11eb-041b-c7083c829ba7
md"""# El camino hacia el Caos
Presentado por el Dr. Chaos ☺

Ecuación logística

$y_{n+1}=ry_n(1-y_n)$
$r \in [0, 4]$
$y \in [0, 1]$

"""

# ╔═╡ 98165ef0-3490-11eb-0e57-259fcf7361c2
logistica(y,r)=r*y*(1-y)

# ╔═╡ 10625a80-7122-11eb-3601-273ea4390a23
md"Lo bueno de las ecuaciones en diferencias es que es trivial resolverlas, basta con iterar"

# ╔═╡ ea63b4f0-3490-11eb-26a8-6f5e1f2d3c05
function trayectoria(y_inicial,r,pasos)
	y=zeros(pasos)
	y[1]=y_inicial
	for i=1:(pasos-1)
		y[i+1]=logistica(y[i],r)
	end
	return y
end

# ╔═╡ 915aa560-0e0e-11eb-18d4-e3a899ee2c1c
function trayectorias!(figura,f,p,rangos;N=10)
	#modifica una figura añadiendo las trayectorias por eso tiene una exclamación !
	#barrido para las soluciones
	xmin,xmax,ymin,ymax=rangos
	N=round(N/2);
	x1s=xmin:((xmax-xmin)/N):xmax
	x2s=ymin:((ymax-ymin)/N):ymax
	T=10;
	for x1 in x1s
		for x2 in x2s
			X,Y=trayectoria([x1;x2],f,p,rangos,N=N,pasosMax=100)
			figura=plot!(figura,X,Y,legend=false)
		end
	end
	return figura
end

# ╔═╡ 7115307e-d660-11ea-3f2f-4d4963a60e39
function fases(f,p;rangos=[-1,1,-1,1],N=10)
	xmin,xmax,ymin,ymax=rangos
	x1s=xmin:((xmax-xmin)/N):xmax
	x2s=ymin:((ymax-ymin)/N):ymax

	p1=plot()

	p1=flechas!(p1,f,p,rangos,N=N) #argumento opcional
	p1=trayectorias!(p1,f,p,rangos,N=N)

	xlims!(p1,(xmin, xmax))
	ylims!(p1,(ymin, ymax))
	return p1
end

# ╔═╡ 90208f82-d4ba-11ea-11e6-7f51e6ffc5d5
fases(f_lineal,[k1,k2])

# ╔═╡ d846eed0-d7f8-11ea-0596-d56c65bbbb5e
fases(f_lineal,[-0.5,-1.5])

# ╔═╡ ee673070-d7f9-11ea-3cb9-330579db1f09
fases(f_lineal,[-2,-2])

# ╔═╡ 86cc4ac0-d7fb-11ea-2101-4b12fafd41f7
fases(f_lineal,[-1,0])

# ╔═╡ 8bc272f0-dd75-11ea-3f24-1766ca4f5a7c
fases(f_lineal,[0,-1])

# ╔═╡ dcb4f750-d7f8-11ea-0027-47491be6b131
fases(f_lineal,[-0.5,1.5])

# ╔═╡ 7dacbd90-d7fa-11ea-1063-e5741bf3f618
fases(f_lineal,[-2,2])

# ╔═╡ df4f1840-d7fa-11ea-1e89-3d6d103458d1
fases(f_lineal,[1,0])

# ╔═╡ f102ad80-dd78-11ea-08b5-2d95e6f1a678
begin
	H=1
	Pe=1
	Pm=Pe*sliderMotor
	dibujo=fases(generador,[H,Pm,Pe],rangos=[-2pi,2pi,-2,2],N=15)
	if abs(Pm/Pe)<=1
		x1a=asin(Pm/Pe)
	    x1b=pi-asin(Pm/Pe)
		scatter!([x1a,x1b],[0,0]) #puntos de equilibrio
	end
    dibujo
end

# ╔═╡ 46dd7a80-de3e-11ea-0e9d-abe7496efef7
begin
	if abs(Pm/Pe)<=1
	   A1=[0 1;
		  Pe*cos(x1a) 0]
	   autov1,_1=eigen(A1)
	   fases((x,p,t)->A1*x,[]) #otra forma e hacerlo con una función anónima
	   title!("""El primer PE es inestable (silla),
		λ1=$(autov1[1]) λ2=$(autov1[2])""")
	end
end

# ╔═╡ 44e6bce2-800b-11eb-3651-e3ba5df8617e
begin
	if abs(Pm/Pe)<=1
	   A2=[0 1;
		  Pe*cos(x1b) 0]
	   autov2,_2=eigen(A2)
	   fases((x,p,t)->A2*x,[]) #otra forma e hacerlo con una función anónima
       title!("""El primer PE es un centro (en la aproximación lineal),
	λ1=$(autov2[1]) λ2=$(autov2[2])""")
	end
end

# ╔═╡ afa04610-d7f4-11ea-1cd0-e3a91d91c5c9
fases(van_der_pol,[1.0,1.0,1.0],rangos=[-3,3,-3,3])

# ╔═╡ b8326ede-3490-11eb-1e0b-530752ed836e
@bind r s=Slider(0:0.05:4,default=2)

# ╔═╡ a07f3bc0-3490-11eb-11b3-6bedb92502cb
let  #introduce un "scope" nuevo no como begin-end lo que permite usar las mismas variables sin tener que inventar nombres nuevos
	x=0:0.01:1
	fig=plot(x,logistica.(x,r),xlim=(0,1),ylim=(0,1),
		label=:none, title="r=$r",ratio=1)
end

# ╔═╡ d43266e0-3490-11eb-0d9a-59ef7a663b50
@bind ver CheckBox(default=false)

# ╔═╡ 165343ee-3491-11eb-10a7-7b4d6fef56f4
function puntos(R)
	px=[]
	py=[]
	for r=0:0.02:R
		yt=trayectoria(0.5,r,100)
		k=length(yt)
		y_est=yt[(k÷2):k]
		append!(px,r*ones(length(y_est))) #append es como push pero añade un vector
		append!(py,y_est)
	end
	return (px,py)
end

# ╔═╡ e042eeee-3490-11eb-1406-058ee6753c1f
let
	y=trayectoria(0.5,r,50)
	p1=plot(y,title="r=$r",marker=:dot,markersize=:2,label=:none,ylimits=(0,1))
	k=length(y)
		p2=scatter()
	if ver
		px,py=puntos(4)
		p2=scatter!(p2,px,py,legend=:none,color=:grey, alpha=0.1,
		            linecolor=:grey, linealpha=0.1,marker=1)
	end
	y_estacionarios=y[(k÷2):k] #Nota ver cómo se escribe esto
	p2=scatter!(p2,r*ones(length(y_estacionarios)),y_estacionarios,
		        xlim=(0,4),ylim=(0,1),legend=:none,color=:red,marker=2)
	l = @layout [a ; b]
	plot(p1,p2, layout=l, size=(650,500))
end

# ╔═╡ 315a35f0-3491-11eb-0ce1-636a82ef3dd9
md"""# ¿Cómo se llega al caos?

Haciendo lo mismo una y otra vez...

Primer paso, periodo 1:"""

# ╔═╡ dd340bd0-3491-11eb-31ab-5b1bf6286cf0
@bind Niter Slider(1:50,default=1,show_value=true)

# ╔═╡ 3765e8e2-3491-11eb-00a5-274f3caf6a15
@bind y0 Slider(0:0.05:1,default=0.5,show_value=true)

# ╔═╡ de04c952-3491-11eb-276b-b530d61da953
@bind r2 Slider(0:0.05:4,default=2,show_value=true)

# ╔═╡ 39b45040-3492-11eb-1e48-998b1e4b1911
function iterar(f,N)
	x=0:0.001:1
	fig2=plot(x,f.(x),xlim=(0,1),ylim=(0,1),label=:none, title="r=$r2")
	plot!(fig2,[0,1],[0,1],color=:red,label=:none, ratio=1)

	x=y0 #punto en el plano de la figura
	y=0
	for i=1:N
		#levanto y hasta f(x)
		xn=x
		yn=f(x)
		plot!(fig2,[x,xn],[y,yn],color=:green,label=:none)
		#desplazo x hasta y manteniendo yn
		x=yn
		y=yn
		plot!(fig2,[xn,x],[yn,y],color=:green,label=:none)
	end
	fig2
end

# ╔═╡ 2e058370-729c-11eb-0822-793d3db17e74
begin
    f1(x)=logistica(x,r2)
	f2(x)=f1(f1(x))
	f4(x)=f2(f2(x))
	f8(x)=f4(f4(x))
	f16(x)=f8(f8(x))
	iterar(f1,Niter) # prueba con f1 f2 ...
end

# ╔═╡ 4c657840-3492-11eb-093d-59a4b3854b7a
md"""En un sistema discreto la estabilidad depende de la derivada, es parecido a los sistemas continuos.
Un punto de equilibrio $x^*$ es aquel en el que $x_{n+1}=x_n$, es decir

$x_{n+1}=f(x_n)=x_n=x^*$

Es decir que es un punto fijo de $f$, es decir $f(x^*)=x^*$

Si estamos cerca de un punto de equilibrio $x_n=x^*+\delta x_n$ luego haciendo taylor:

$f(x^*+\delta x_n)=f(x^*) + \nabla f(x^*)\delta x_n + o(\delta x_n^2)$

Pero entonces:

$\delta x_{n+1}=x_{n+1}-x^*=f(x^*)-x^* + \nabla f(x^*)\delta x_n + o(\delta x_n^2) \approx \nabla f(x^*) \delta x_n=A \delta x_n$

Cuya solución es $\delta x_n=A^n\delta x_0$ que es estable si los autovalores de $A$ están en el círculo unitario. Es parecido a los sistemas continuos cambiando el eje imaginario por el **círculo unidad**.

En el caso de una variable se reduce a que $|f'(x^*)|<1$ como se intuye gráficamente

"""

# ╔═╡ 55e2ab40-3492-11eb-2953-e9d8bd17735c
md"En cada paso el cambio de $r$ necesario para producir una bifurcación es más pequeño.

De hecho la bifurcación de orden 4 es como una versión pequeña de la de orden 2 como hemos visto, si llamamos $r_n$ al $r$ que produce una bifurcación tenemos que

$\frac{r_{n+2}-r_{n+1}}{r_{n+1}-r_{n}} \to 4.6992...$ que es la constante de Feigenbaum y que es **la misma** para muchos otros procesos esto se conoce como **Universalidad**.

Como el paso sigue aproximadamente una serie geométrica la suma de infinitas bifurcaciones se produce con un cambio finito de $r$ ya que:

$r_\infty-r_0 \approx \Delta r_0 +  \frac {\Delta r_0}{4.7} +\frac {\Delta r_0}{4.7^2}+\frac {\Delta r_0}{4.7^3}...=\Delta r_0 \frac {1}{1-1/4.7}\approx 1.27\Delta r_0$

"

# ╔═╡ Cell order:
# ╟─6ba02e80-d4cf-11ea-19cb-33028d3dc67d
# ╟─7db3fdb0-d4d2-11ea-21fc-597c55dda572
# ╟─76079a70-d4cf-11ea-2645-8f43c78eae52
# ╠═66e9ffd0-d4d7-11ea-09c0-03ebaacab879
# ╟─0b94f1c0-7113-11eb-13da-c5c5b95ebbf5
# ╠═e885f5b0-d4d4-11ea-06d6-9ddf9619961f
# ╠═b99fc9b0-d4d9-11ea-0e18-cfd5af90188f
# ╠═88629870-3490-11eb-2215-d765d59a440c
# ╠═3be6b020-d4d8-11ea-3356-23afd148ebc9
# ╠═220f6f12-d4d9-11ea-159e-c1e1551fd732
# ╟─23f5ece0-d4da-11ea-1330-c5d077bd3a89
# ╠═d8b5ed40-d4d7-11ea-27d4-bf594518802b
# ╠═1be7c0c0-d4ba-11ea-1e43-ef66c35b0830
# ╠═2b7df180-d4d8-11ea-118f-75e8709e36c3
# ╟─455e0c30-d4d7-11ea-22ce-0d04b4e8d860
# ╠═159f28e2-d4db-11ea-13e6-a97e2e7f7c60
# ╠═c90980d0-d637-11ea-112d-c336b2080d87
# ╠═1f44a88e-7114-11eb-0ac2-d12cd3fb40ef
# ╠═4a7f5f72-0e0c-11eb-2b7c-4d7214c09c4b
# ╠═51059990-0e0c-11eb-2d55-89ac739619c7
# ╠═a92de300-d63c-11ea-2446-5b099592d734
# ╟─09980bd0-d63d-11ea-0752-1d033ee00cc9
# ╠═e6afc70e-d63d-11ea-0d36-738b13f0de30
# ╠═6cbb32f0-d665-11ea-3f24-41cbe1910b80
# ╟─5a802550-d665-11ea-3f9f-b103ce80822f
# ╟─2ff44950-7116-11eb-3ad9-c3e145cdae4e
# ╟─a0573fe0-d6a2-11ea-057e-5140c323c005
# ╠═90208f82-d4ba-11ea-11e6-7f51e6ffc5d5
# ╟─1790c760-0e0d-11eb-1c47-f5e902e158fb
# ╠═f5ca89e0-0e0c-11eb-32a1-fb22511f42ae
# ╠═165464e2-0e0f-11eb-3851-3701a6c78502
# ╠═915aa560-0e0e-11eb-18d4-e3a899ee2c1c
# ╠═7115307e-d660-11ea-3f2f-4d4963a60e39
# ╟─b51d4d90-d7f4-11ea-07aa-adcb4f418c99
# ╟─0b2775be-d7fb-11ea-0b82-0d2b4d544d31
# ╟─4100c3f0-d7f5-11ea-10c3-593316528d50
# ╠═d846eed0-d7f8-11ea-0596-d56c65bbbb5e
# ╟─f8cc9600-d7f8-11ea-138d-5f40babd5dce
# ╠═ee673070-d7f9-11ea-3cb9-330579db1f09
# ╟─34bbafa0-d7fb-11ea-2ed3-8b300f45fbac
# ╟─345ec290-d7fb-11ea-1f16-d919242fefc3
# ╠═86cc4ac0-d7fb-11ea-2101-4b12fafd41f7
# ╟─8c6fa2e0-dd75-11ea-163b-c70c276c3f2a
# ╟─ad48c140-dd75-11ea-03fe-75a0b19f91bd
# ╠═8bc272f0-dd75-11ea-3f24-1766ca4f5a7c
# ╟─33f181d0-d7fb-11ea-0a51-d3a11d1613df
# ╟─fa4a8370-d7f8-11ea-12b0-fde507f287be
# ╠═dcb4f750-d7f8-11ea-0027-47491be6b131
# ╟─08c0cb70-d7fa-11ea-0efd-ab3b7f76a8b8
# ╠═7dacbd90-d7fa-11ea-1063-e5741bf3f618
# ╟─b6db1b70-d7fa-11ea-2bea-43295b3179a6
# ╠═df4f1840-d7fa-11ea-1e89-3d6d103458d1
# ╟─fd823200-d7fc-11ea-1af8-c91386a02ab5
# ╠═2647d140-d7fd-11ea-201a-5744621fbe12
# ╠═7c13c6b0-d7fd-11ea-0cc5-55f5ab6106fe
# ╟─82023190-dd77-11ea-12d1-cfa1a2496348
# ╟─16d7ee1e-dd7a-11ea-13a7-197edb480651
# ╠═ab5bd370-dd77-11ea-034f-6153487b9658
# ╟─1668b190-dd7a-11ea-3c80-c5f082ec8647
# ╟─351f5332-dd7c-11ea-2fc4-c9dc71684f23
# ╟─5fd238f0-dd7b-11ea-2a00-a587ee59c261
# ╠═f102ad80-dd78-11ea-08b5-2d95e6f1a678
# ╟─db58def0-dd78-11ea-335b-4140d22d9a78
# ╟─46dd7a80-de3e-11ea-0e9d-abe7496efef7
# ╟─44e6bce2-800b-11eb-3651-e3ba5df8617e
# ╟─33af7e82-d7f5-11ea-12bb-5982575f35d3
# ╠═3d70984e-d7f5-11ea-3243-9d3932ee155b
# ╠═afa04610-d7f4-11ea-1cd0-e3a91d91c5c9
# ╟─713de9b0-3490-11eb-041b-c7083c829ba7
# ╠═98165ef0-3490-11eb-0e57-259fcf7361c2
# ╟─10625a80-7122-11eb-3601-273ea4390a23
# ╠═ea63b4f0-3490-11eb-26a8-6f5e1f2d3c05
# ╠═a07f3bc0-3490-11eb-11b3-6bedb92502cb
# ╟─b8326ede-3490-11eb-1e0b-530752ed836e
# ╟─d43266e0-3490-11eb-0d9a-59ef7a663b50
# ╠═e042eeee-3490-11eb-1406-058ee6753c1f
# ╠═165343ee-3491-11eb-10a7-7b4d6fef56f4
# ╟─315a35f0-3491-11eb-0ce1-636a82ef3dd9
# ╠═dd340bd0-3491-11eb-31ab-5b1bf6286cf0
# ╠═3765e8e2-3491-11eb-00a5-274f3caf6a15
# ╠═de04c952-3491-11eb-276b-b530d61da953
# ╠═2e058370-729c-11eb-0822-793d3db17e74
# ╠═39b45040-3492-11eb-1e48-998b1e4b1911
# ╟─4c657840-3492-11eb-093d-59a4b3854b7a
# ╟─55e2ab40-3492-11eb-2953-e9d8bd17735c

max_stars_repo_path: Tema 1.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 1.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 1.jl

max_forks_repo_name: Dictino/CNL


text:
include("dkk17.jl")

function rcov_quantum_filter(reps, eps, k, α=4, τ=0.1; limit1=2, limit2=1.5)
    d, n = size(reps)
    reps_pca, U = pca(reps, k)
    if k == 1
        reps_estimated_white = reps_pca
        Σ′ = ones(1, 1)
    else
        selected = cov_estimation_iterate(reps_pca, eps/n, τ, nothing, limit=round(Int, limit1*eps))
        reps_pca_selected = reps_pca[:, selected]
        Σ = cov(reps_pca_selected', corrected=false)
        reps_estimated_white = Σ^(-1/2)*reps_pca
        Σ′ = cov(reps_estimated_white')
    end
    M = k > 1 ? exp(α*(Σ′- I)/(opnorm(Σ′) - 1)) : ones(1, 1)
    M /= tr(M)
    estimated_poison_ind = k_lowest_ind(
        -[x'M*x for x in eachcol(reps_estimated_white)],
        round(Int, limit2*eps)
    )
    return .! estimated_poison_ind
end

function rcov_auto_quantum_filter(reps, eps, α=4, τ=0.1; limit1=2, limit2=1.5)
    reps_pca, U = pca(reps, 100)
    best_opnorm, best_selected, best_k = -Inf, nothing, nothing
    for k in round.(Int, range(1, sqrt(100), length=10) .^ 2)
        selected = rcov_quantum_filter(reps, eps, k, α, τ; limit1=limit1, limit2=limit2)
        Σ = cov(reps_pca[:, selected]')
        Σ′ = cov((Σ^(-1/2)*reps_pca)')
        M = exp(α*(Σ′- I)/(opnorm(Σ′) - 1))
        M /= tr(M)
        op = tr(Σ′ * M) / tr(M)
        poison_removed = sum((.! selected)[end-eps+1:end])
        @show k, op, poison_removed
        if op > best_opnorm
            best_opnorm, best_selected, best_k = op, selected, k
        end
    end
    @show best_k, best_opnorm
    return best_selected
end

max_stars_repo_name: TDteach/spectre-defense

max_issues_repo_name: TDteach/spectre-defense

max_forks_repo_name: TDteach/spectre-defense


text:
@testset "1332.remove-palindromic-subsequences.jl" begin
    @test remove_palindromic_subsequences("ababa") == 1
    @test remove_palindromic_subsequences("abb") == 2
    @test remove_palindromic_subsequences("baabb") == 2
end


text:
module mbox1
using MeshSteward: initbox, updatebox!, inflatebox!, inbox, boxesoverlap
using Test
function test()
	box = Float64[]
	box = initbox([0.0])
	@test isapprox(box, [0.0, 0.0])

	box = updatebox!(box, [-1.0])
	@test isapprox(box, [-1.0, 0.0])

	box = Float64[]
	box = initbox([1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test isapprox(box, [0.99, 1.01])

	box = initbox([0.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test inbox(box, [-0.003]) == true

	box = initbox([2.0])
	updatebox!(box, [1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test inbox(box, [rand()+1.0]) == true

	box = Float64[]
	box = initbox([0.0; 1.0])
	@test isapprox(box, [0.0, 0.0, 1.0, 1.0])

	box = updatebox!(box, [1.0; -1.0])
	@test isapprox(box, [0.0, 1.0, -1.0, 1.0])

	box = Float64[]
	box = initbox([0.0; 1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test isapprox(box, [-0.01, 0.01, 0.99, 1.01])

	box = initbox([0.0; 1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test inbox(box, [-0.003; 0.995]) == true

	box = initbox([0.0; 0.0])
	updatebox!(box, [1.0; 1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test inbox(box, [-0.003; 0.995]) == true
	@test inbox(box, [rand(); rand()]) == true
	@test inbox(box, [0.5; 0.5]) == true
	@test inbox(box, [0.0-inflate; 0.0-inflate]) == true
	@test inbox(box, [0.0-inflate; 0.0-2*inflate]) == false
	@test inbox(box, [0.0+inflate; 1.0+inflate]) == true

	box = initbox([0.0; 0.0; 0.0].-1.0)
	updatebox!(box, [1.0; 1.0; 1.0])
	inflate = 0.01
	box = inflatebox!(box, inflate)
	@test inbox(box, [-0.003; 0.995; 1.0]) == true
	@test inbox(box, [rand(); rand(); rand()]) == true
	@test inbox(box, [-rand(); rand(); rand()]) == true
	@test inbox(box, [rand(); -rand(); rand()]) == true
	@test inbox(box, [rand(); rand(); -rand()]) == true
	@test inbox(box, [0.0-inflate; 0.0-inflate; 1.0]) == true
	@test inbox(box, [0.0-inflate; 0.0-inflate; -1.0]) == true

	return true
end
end
using .mbox1
mbox1.test()

module mbox2
using MeshSteward: initbox, updatebox!, inflatebox!, inbox, boxesoverlap, boundingbox
using LinearAlgebra
using Test
function test()
	a = [0.431345 0.611088 0.913161;
	    0.913581 0.459229 0.82186;
	    0.999429 0.965389 0.571139;
	    0.390146 0.711732 0.302495;
	    0.873037 0.248077 0.51149;
	    0.999928 0.832524 0.93455]
	    b1 = boundingbox(a)
	    @test norm(b1 - [0.390146, 0.999928, 0.248077, 0.965389, 0.302495, 0.93455]) < 1.0e-4
	    b2 = updatebox!(b1, a)
	    @test norm(b1 - b2) < 1.0e-4
	    b2 = initbox(a)
	    @test norm(b1 - b2) < 1.0e-4
	    c = [-1.0, 3.0, -0.5]
	    b3 = updatebox!(b1, c)
	    # # println("$(b3)")
	    @test norm(b3 - [-1.0, 0.999928, 0.248077, 3.0, -0.5, 0.93455]) < 1.0e-4
	    x = [0.25 1.1 -0.3]
	    @test inbox(b3, x)
	    @test inbox(b3, c)
	    @test inbox(b3, a[2, :])
	    b4 = boundingbox(-a)
	    # # println("$(b3)")
	    # # println("$(b4)")
	    # # println("$(boxesoverlap(b3, b4))")
	    @test !boxesoverlap(b3, b4)
	    b5 = updatebox!(b3, [0.0 -0.4 0.0])
	    # # println("$(b5)")
	    # # println("$(boxesoverlap(b5, b4))")
	    @test boxesoverlap(b5, b4)
end
end
using .mbox2
mbox2.test()

module mboxt3
using MeshSteward: intersectboxes
using LinearAlgebra
using Test
function test()
	b1 = [-1.0, 2.0, 0.5, 2.5]
	b2 = [-0.7, 1.2, -0.5, 3.5]
	ib = intersectboxes(b1, b2)
	@test isapprox(ib, [-0.7, 1.2, 0.5, 2.5])
	b1 = [-1.0, 2.0, 0.5, 2.5, 77.1, 1000.0]
	b2 = [-0.7, 1.2, -0.5, 3.5, 1077.1, 10000.0]
	ib = intersectboxes(b1, b2)
	@test isapprox(ib, Float64[])
end
end
using .mboxt3
mboxt3.test()


text:
function setup_polishing()
    options = Dict(:verbose => false,
                   :polish => true,
                   :eps_abs => 1e-03,
                   :eps_rel => 1e-03,
                   :verbose => false,
                   :max_iter => 5000)
    return options
end

tol = 1e-3

@testset "polishing" begin

    @testset "polishing_problem" begin
        P = sparse(Diagonal([11.; 0.]))
        q = [3.; 4.]
        A = sparse([-1. 0.; 0. -1.; -1. -3; 2. 5.; 3. 4.])
        u = [0.; 0.; -15.; 100.; 80]
        l = -Inf * ones(length(u))
        (m, n) = size(A)

        # Solve problem
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        x_test = [9.90341e-11; 5.0]
        y_test = [1.66667; 0.0; 1.33333; 1.20431e-14; 1.49741e-14]
        obj_test = 20.

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, y_test, atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status_polish == 1

    end

    @testset "polishing_unconstrained" begin

        srand(1)

        n = 10
        m = n
        P = sparse(Diagonal(randn(n)) + 1. * sparse(I, n, n))
        q = randn(n)
        A = sparse(I, n, n)
        l = -100 * ones(m)
        u = 100 * ones(m)

        # Solve problem
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        # Explicit solution
        invP = inv(Array(P))
        x_test = - invP * q
        y_test = zeros(m)
        obj_test = - .5 * q' * invP * q

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, zeros(m), atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status_polish == 1
    end

    @testset "polish_random" begin
        srand(1)

        n = 30
        m = 50
        Pt = sprandn(n, n, 0.5)
        P = Pt * Pt'
        q = randn(n)
        A = sprandn(m, n, 0.5)
        l = -3 .+ randn(m)
        u = 3 .+ randn(m)

        # Solve problem
        options = setup_polishing()
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        results = OSQP.solve!(model)

        # # Solve with Gurobi
        # using Gurobi
        # env = Gurobi.Env()
        # setparam!(env, "OutputFlag", 1)
        # model = Gurobi.Model(env, "model")
        # add_cvars!(model, q); update_model!(model)
        # add_constrs!(model, A, repmat(['<'], m), u); update_model!(model)
        # add_constrs!(model, A, repmat(['>'], m), l); update_model!(model)
        # add_qpterms!(model, P); update_model!(model)
        # optimize(model)
        # x_test = get_solution(model)
        # obj_test= get_objval(model)
        # y_test = -Gurobi.get_dblattrarray(model, "Pi", 1, 2 * m)
        # y_test = y_test[m+1:end] + y_test[1:m]
                #
        # println(x_test)
        # println(y_test)
        # println(obj_test)

        x_test = [0.369834, -0.0209277, 0.068939, -0.604151, 0.60773, -0.715965, -0.128837, -0.593642, 0.928932, 0.678794, 0.194334, 0.228124, 0.188779, 0.607246, 0.479752, 0.0878951, 0.275489, 0.761593, 0.0600479, -0.950663, -0.191354, 1.12286, -0.593261, -0.0932879, 0.0286826, -0.0257784, 0.358456, -0.0693514, -0.00358498, -0.0270959]
        y_test = [4.88779e-14, 0.441478, 5.89704e-14, 1.76439e-13, 8.10529e-15, 4.9186e-14, 5.00637e-14, 2.00889e-14, 5.39402e-13, -1.9796e-14, 1.52041e-13, 0.149979, 8.97844e-12, 5.93941e-14, -4.11719e-14, -2.38358e-13, 0.226114, -1.8933e-14, 1.54839e-13, 2.99422e-14, 7.16687e-15, -5.32536e-14, 0.116129, -2.53245e-14, 1.21586e-14, -1.66761e-13, 3.04905e-14, 7.35909e-14, 3.84878e-12, -3.09855e-15, 4.30889e-9, 4.09921e-14, -2.4863e-14, -2.46841e-14, 4.37156e-15, 0.219029, -5.14155e-14, 8.31705e-14, 4.51252e-15, -7.90941e-14, 0.206435, -0.196264, 1.06762e-14, 0.338783, 7.58781e-10, -0.105739, 3.30014e-13, -4.05467e-14, 3.3212e-14, -0.161755]

        obj_test = -6.021607803961248

        @test isapprox(results.x, x_test, atol=tol)
        @test isapprox(results.y, y_test, atol=tol)
        @test isapprox(results.info.obj_val, obj_test, atol=tol)
        @test results.info.status_polish == 1

    end
end


text:
module GoogleMaps

using ArgCheck
using HTTP
using JSON
using Dates

include("geocoding.jl")
export geocode

include("timezone.jl")
export timezone

end


text:
"""
A Sweeps objects holds information
about the various parameters controlling
a density matrix renormalization group (DMRG)
or similar matrix product state (MPS) calculation.

For a Sweeps object `sw` the available
parameters are:
* `nsweep(sw)` -- the number of sweeps to do
* `maxdim(sw,n)` -- maximum MPS bond dimension for sweep n
* `mindim(sw,n)` -- minimum MPS bond dimension for sweep n
* `cutoff(sw,n)` -- truncation error cutoff for sweep n
* `noise(sw,n)` -- noise term coefficient for sweep n
"""
mutable struct Sweeps
  nsweep::Int
  maxdim::Vector{Int}
  cutoff::Vector{Float64}
  mindim::Vector{Int}
  noise::Vector{Float64}

  function Sweeps(nsw::Int)
    return new(nsw, fill(1, nsw), zeros(nsw), fill(1, nsw), zeros(nsw))
  end
end

"""
    Sweeps(d::AbstractMatrix)

    Sweeps(nsweep::Int, d::AbstractMatrix)

Make a sweeps object from a matrix of input values.
The first row should be strings that define which
variables are being set ("maxdim", "cutoff", "mindim",
and "noise").

If the number of sweeps are not specified, they
are determined from the size of the input matrix.

# Examples
```julia
julia> Sweeps([
         "maxdim" "mindim" "cutoff" "noise"
          50       10       1e-12    1E-7
          100      20       1e-12    1E-8
          200      20       1e-12    1E-10
          400      20       1e-12    0
          800      20       1e-12    1E-11
          800      20       1e-12    0
         ])
Sweeps
1 cutoff=1.0E-12, maxdim=50, mindim=10, noise=1.0E-07
2 cutoff=1.0E-12, maxdim=100, mindim=20, noise=1.0E-08
3 cutoff=1.0E-12, maxdim=200, mindim=20, noise=1.0E-10
4 cutoff=1.0E-12, maxdim=400, mindim=20, noise=0.0E+00
5 cutoff=1.0E-12, maxdim=800, mindim=20, noise=1.0E-11
6 cutoff=1.0E-12, maxdim=800, mindim=20, noise=0.0E+00
```
"""
function Sweeps(nsw::Int, d::AbstractMatrix)
  sw = Sweeps(nsw)
  vars = d[1, :]
  for (n, var) in enumerate(vars)
    inputs = d[2:end, n]
    if var == "maxdim"
      maxdim!(sw, inputs...)
    elseif var == "cutoff"
      cutoff!(sw, inputs...)
    elseif var == "mindim"
      mindim!(sw, inputs...)
    elseif var == "noise"
      noise!(sw, float.(inputs)...)
    else
      error("Sweeps object does not have the field $var")
    end
  end
  return sw
end

Sweeps(d::AbstractMatrix) = Sweeps(size(d, 1) - 1, d)

"""
    nsweep(sw::Sweeps)
    length(sw::Sweeps)

Obtain the number of sweeps parameterized
by this sweeps object.
"""
nsweep(sw::Sweeps)::Int = sw.nsweep

Base.length(sw::Sweeps)::Int = sw.nsweep

"""
    maxdim(sw::Sweeps,n::Int)

Maximum MPS bond dimension allowed by the
Sweeps object `sw` during sweep `n`
"""
maxdim(sw::Sweeps, n::Int)::Int = sw.maxdim[n]

"""
    mindim(sw::Sweeps,n::Int)

Minimum MPS bond dimension allowed by the
Sweeps object `sw` during sweep `n`
"""
mindim(sw::Sweeps, n::Int)::Int = sw.mindim[n]

"""
    cutoff(sw::Sweeps,n::Int)

Truncation error cutoff setting of the
Sweeps object `sw` during sweep `n`
"""
cutoff(sw::Sweeps, n::Int)::Float64 = sw.cutoff[n]

"""
    noise(sw::Sweeps,n::Int)

Noise term coefficient setting of the
Sweeps object `sw` during sweep `n`
"""
noise(sw::Sweeps, n::Int)::Float64 = sw.noise[n]

get_maxdims(sw::Sweeps) = sw.maxdim
get_mindims(sw::Sweeps) = sw.mindim
get_cutoffs(sw::Sweeps) = sw.cutoff
get_noises(sw::Sweeps) = sw.noise

"""
    maxdim!(sw::Sweeps,maxdims::Int...)

Set the maximum MPS bond dimension for each
sweep by providing up to `nsweep(sw)` values.
If fewer values are provided, the last value
is repeated for the remaining sweeps.
"""
function setmaxdim!(sw::Sweeps, maxdims::Int...)::Nothing
  mdims = collect(maxdims)
  for i in 1:nsweep(sw)
    sw.maxdim[i] = get(mdims, i, maxdims[end])
  end
end
maxdim!(sw::Sweeps, maxdims::Int...) = setmaxdim!(sw, maxdims...)

"""
    mindim!(sw::Sweeps,maxdims::Int...)

Set the minimum MPS bond dimension for each
sweep by providing up to `nsweep(sw)` values.
If fewer values are provided, the last value
is repeated for the remaining sweeps.
"""
function setmindim!(sw::Sweeps, mindims::Int...)::Nothing
  mdims = collect(mindims)
  for i in 1:nsweep(sw)
    sw.mindim[i] = get(mdims, i, mindims[end])
  end
end
mindim!(sw::Sweeps, mindims::Int...) = setmindim!(sw, mindims...)

"""
    cutoff!(sw::Sweeps,maxdims::Int...)

Set the MPS truncation error used for each
sweep by providing up to `nsweep(sw)` values.
If fewer values are provided, the last value
is repeated for the remaining sweeps.
"""
function setcutoff!(sw::Sweeps, cutoffs::Real...)::Nothing
  cuts = collect(cutoffs)
  for i in 1:nsweep(sw)
    sw.cutoff[i] = get(cuts, i, cutoffs[end])
  end
end
cutoff!(sw::Sweeps, cutoffs::Real...) = setcutoff!(sw, cutoffs...)

"""
    noise!(sw::Sweeps,maxdims::Int...)

Set the noise-term coefficient used for each
sweep by providing up to `nsweep(sw)` values.
If fewer values are provided, the last value
is repeated for the remaining sweeps.
"""
function setnoise!(sw::Sweeps, noises::Real...)::Nothing
  nvals = collect(noises)
  for i in 1:nsweep(sw)
    sw.noise[i] = get(nvals, i, noises[end])
  end
end
noise!(sw::Sweeps, noises::Real...) = setnoise!(sw, noises...)

function Base.show(io::IO, sw::Sweeps)
  println(io, "Sweeps")
  for n in 1:nsweep(sw)
    @printf(
      io,
      "%d cutoff=%.1E, maxdim=%d, mindim=%d, noise=%.1E\n",
      n,
      cutoff(sw, n),
      maxdim(sw, n),
      mindim(sw, n),
      noise(sw, n)
    )
  end
end

struct SweepNext
  N::Int
  ncenter::Int
end

"""
    sweepnext(N::Int; ncenter::Int=2)

Returns an iterable object that evaluates
to tuples of the form `(b,ha)` where `b`
is the bond number and `ha` is the half-sweep
number. Takes an optional named argument
`ncenter` for use with an n-site MPS or DMRG
algorithm, with a default of 2-site.
"""
function sweepnext(N::Int; ncenter::Int=2)::SweepNext
  if ncenter < 0
    error("ncenter must be non-negative")
  end
  return SweepNext(N, ncenter)
end

function Base.iterate(sn::SweepNext, state=(0, 1))
  b, ha = state
  if ha == 1
    inc = 1
    bstop = sn.N - sn.ncenter + 2
  else
    inc = -1
    bstop = 0
  end
  new_b = b + inc
  new_ha = ha
  done = false
  if new_b == bstop
    new_b -= inc
    new_ha += 1
    if ha == 2
      return nothing
    end
  end
  return ((new_b, new_ha), (new_b, new_ha))
end


text:
export iauXys06a
"""
For a given TT date, compute the X,Y coordinates of the Celestial
Intermediate Pole and the CIO locator s, using the IAU 2006
precession and IAU 2000A nutation models.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double  TT as a 2-part Julian Date (Note 1)

Returned:
   x,y          double  Celestial Intermediate Pole (Note 2)
   s            double  the CIO locator s (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The Celestial Intermediate Pole coordinates are the x,y components
   of the unit vector in the Geocentric Celestial Reference System.

   3. The CIO locator s (in radians) positions the Celestial
   Intermediate Origin on the equator of the CIP.

   4. Series-based solutions for generating X and Y are also available:
   see Capitaine & Wallace (2006) and iauXy06.

Called:
   iauPnm06a    classical NPB matrix, IAU 2006/2000A
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauS06       the CIO locator s, given X,Y, IAU 2006

References:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauXys06a(double date1, double date2,
#                double *x, double *y, double *s)

function iauXys06a(date1::Real, date2::Real)
   ref_x = Ref{Float64}(0.0)
   ref_y = Ref{Float64}(0.0)
   ref_s = Ref{Float64}(0.0)

   ccall((:iauXys06a, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
         convert(Float64, date1), convert(Float64, date2),
         ref_x, ref_y, ref_s)

   return ref_x[], ref_y[], ref_s[]
end