text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

import Statistics

"""
    @unzip xs, ys, ... = us

will expand the assignment into the following code
    xs, ys, ... = map(x -> x[1], us), map(x -> x[2], us), ...
"""
macro unzip(args)
    args.head != :(=) && error("Expression needs to be of form `xs, ys, ... = us`")
    lhs, rhs = args.args
    items = isa(lhs, Symbol) ? [lhs] : lhs.args
    rhs_items = [:(map(x -> x[$i], $rhs)) for i in 1:length(items)]
    rhs_expand = Expr(:tuple, rhs_items...)
    esc(Expr(:(=), lhs, rhs_expand))
end
export @unzip

"""Contains the information defining a lenslet array."""
struct LensletSystem{T<:Real}
    eyebox_rectangle::Rectangle{T}
    subdivisions_of_eyebox::Tuple{Int64,Int64}
    lenses::Vector{ParaxialLens{T}}
    lattice_coordinates::Vector{Tuple{Int64,Int64}}
    lenslet_colors::Vector{String}
    projected_eyeboxes::Vector{SMatrix{3,4,T}}
    displayplanes::Vector{Plane{T,3}}
    lenslet_eye_boxes::Vector{SMatrix{3, 4, T, 12}}
    lenslet_eyebox_numbers::Vector{Int64}
    lenslet_eyebox_centers::Vector{SVector{3,T}}
    system_properties::Dict
    subdivided_eyebox_polys::Vector{SMatrix{3,4,T}}
    projected_eyebox_polygons::Vector{ConvexPolygon{N,T}} where{N}
end

# """ Computes the location of the optical center of the lens that will project the centroid of the display to the centroid of the eyebox. Normally the display centroid will be aligned with the geometric centroid of the lens, rather than the optical center of the lens."""
# function compute_optical_center(eyeboxcentroid,displaycenter,lensplane)
#     v = displaycenter - eyeboxcentroid
#     r = Ray(eyeboxcentroid,v)
#     intsct = point(closestintersection(surfaceintersection(lensplane,r),false)) #this seems complicated when you don't need CSG
#     @assert intsct !== nothing
#     return intsct
# end

""" Computes the location of the optical center of the lens that will project the centroid of the display to the centroid of the eyebox. Normally the display centroid will be aligned with the geometric centroid of the lens, rather than the optical center of the lens."""
function compute_optical_center(eyeboxcentroid,display_center,lens)
    lens_geometric_center = centroid(lens) #geometric and optical center coincide at this point
    v = eyeboxcentroid - lens_geometric_center
    r = Ray(display_center,v)
    #optical center may not lie inside the shape of the lens. surfaceintersection(lens,r) will return nothing in this case which will cause the setup_system code to crash. Instead intersect with the plane of the shape of the lens, which has infinite extent.
    surf = OpticSim.plane(lens)
    intsct = point(closestintersection(surfaceintersection(surf,r),false)) #this seems complicated when you don't need CSG
    @assert intsct !== nothing
    return intsct
end
export compute_optical_center

function localframe(a::ParaxialLens)
    if OpticSim.shape(a) isa ConvexPolygon
        return OpticSim.localframe(OpticSim.shape(a))
    else
        return nothing
    end
end

"""Paraxial lenses are created with a local transform of identity. This is confusing because the lens has a local transform and the shape in the lens has a local transform."""
function replace_optical_center(eyeboxcentroid,displaycenter,lens)
    world_optic_center = compute_optical_center(eyeboxcentroid,displaycenter,lens)
    world_to_local_lens = inv(localframe(lens))
    local_optic_center = world_to_local_lens * world_optic_center
    @assert isapprox(0.0, local_optic_center[3],atol = 1e-12) "z should be zero in the local frame but instead it is $(local_optic_center[3])"
    local_optic_center = SVector{2}((local_optic_center)[1:2]) #this is defined in the 2D lens plane so z = 0
    # need the untransformed convexpoly vertices
    return ParaxialLensConvexPoly(focallength(lens),shape(lens),local_optic_center)
end

"""returns display plane represented in world coordinates, and the center point of the display"""
function display_plane(lens)
    center_point = centroid(lens) + -OpticSim.normal(lens)* OpticSim.focallength(lens)
    pln = Plane(OpticSim.normal(lens), center_point, vishalfsizeu = .5, vishalfsizev = .5,interface = opaqueinterface())
    return pln,center_point
end
export display_plane

""" eyebox rectangle represented as a 3x4 SMatrix with x,y,z coordinates in rows 1,2,3 respectively. By default the eyebox z is 0.0"""
function eyeboxpolygon(xsize,ysize)::SMatrix{3,4}
    xext,yext = (xsize,ysize) ./ 2.0 #divide by 2 to get min and max extensions
    return SMatrix{3,4}([xext;yext;0.0 ;; xext;-yext;0.0 ;; -xext;-yext;0.0 ;; -xext;yext;0.0])
end

subpoints(pts::AbstractVector,subdivisions) = reshape(collect(range(extrema(pts)...,subdivisions)),1,subdivisions)
export subpoints

subdivide(poly::AbstractMatrix,xsubdivisions,ysubdivisions) = subdivide(SMatrix{3,4}(poly),xsubdivisions,ysubdivisions)

"""poly is a two dimensional rectangle represented as a 3x4 SMatrix with x values in row 1 and y values in row 2. Returns a vector of length xsubdivisions*ysubdivisions containing all the subdivided rectangles."""
function subdivide(poly::T, xsubdivisions, ysubdivisions)::Vector{T} where{T<:SMatrix{3,4}}
    result = Vector{T}(undef,0) #efficiency not an issue here so don't need to preallocate vector

    x = subpoints(poly[1,:],xsubdivisions+1) #range function makes one less subdivision that people will be expecting, e.g., collect(range(1,3,2)) returns [1,3]
    y = subpoints(poly[2,:],ysubdivisions+1)

    for i in 1:length(x) -1
        for j in 1:length(y) -1
            push!(result,T([x[i];y[j];0.0mm ;; x[i+1];y[j];0.0mm ;; x[i+1];y[j+1];0.0mm ;; x[i];y[j+1];0.0mm])) #this is klunky but the polygon vertices have units of mm so the 0.0 terms must as well.
        end
    end
    return result
end
export subdivide

function eyebox_number(tilecoords::NTuple{2,Int64},cluster::R,num_eyeboxes::Integer)::Int64 where {R<:AbstractLatticeCluster}
    _, tile_index = cluster_coordinates_from_tile_coordinates(cluster,tilecoords)
    eyeboxnum = mod(tile_index-1,num_eyeboxes) + 1
    return eyeboxnum
end

function eyebox_assignment(tilecoords::NTuple{2,Int64},cluster::R,eyeboxes) where{R<:AbstractLatticeCluster}
    #compute cluster
    eyeboxnum = eyebox_number(tilecoords,cluster,length(eyeboxes))
    return eyeboxes[eyeboxnum] #use linear index for Matrix.
end

"""given the eye box polygon and how it is to be subdivided computes the subdivided eyeboxes and centroids"""
function compute_lenslet_eyebox_data(eyeboxtransform,eyeboxpoly::SMatrix{3,4,T,12},subdivisions::Tuple{Int64,Int64})::Vector{SMatrix{3,4,Float64}} where{T} #can't figure out how to get the system to accept a type for an SMatrix of a Unitful quantity which is what eyeboxpoly is. Want to extract the number type from the Quantity type but this is not happening.
    subdivided_eyeboxpolys = subdivide(eyeboxpoly,subdivisions...)
    strippedpolys =  map(x-> ustrip.(mm,x),  subdivided_eyeboxpolys)
    subdivided_eyeboxpolys =[eyeboxtransform * x for x in strippedpolys] #these have units of mm which don't interact well with Transform.
    return subdivided_eyeboxpolys
end

function test_compute_lenslet_eyebox_data(eyeboxtransform,eyeboxpoly,subdivisions)
    subpolys = compute_lenslet_eyebox_data(eyeboxtransform,eyeboxpoly,subdivisions)
end

function project_eyebox_to_display_plane(eyeboxpoly::AbstractMatrix{T},lens,displayplane) where{T<:Real}
    rowdim,coldim = size(eyeboxpoly)
    rays = [Ray(SVector{rowdim}(point),opticalcenter(lens)-SVector{rowdim}(point)) for point in eachcol(eyeboxpoly)]

    points = collect([point(closestintersection(surfaceintersection(displayplane,ray),false)) for ray in rays])

    eyebox = SMatrix{rowdim,coldim}(reinterpret(Float64,points)...)
    threeDpts, toworld, _ = projectonbestfitplane(eyebox,[0.0,0.0,-1.0])
    twoDpts = reshape(threeDpts[1:2,:],2*size(threeDpts)[2])
    twoDpts = collect(reinterpret(SVector{2,T},twoDpts))
    polygon = ConvexPolygon(toworld,twoDpts,opaqueinterface())
    for (eyept,polypt) in zip(eachcol(eyebox),eachcol(vertices(polygon)))
        @assert isapprox(eyept,polypt)
    end

    return eyebox,polygon
end
export project_eyebox_to_display_plane

"""System parameters for a typical HMD."""
function systemparameters()
    return (eye_box = (10.0mm,9.0mm),fov = (90.0°,60.0°),eye_relief = 20.0mm, pupil_diameter = 4.0mm, display_sphere_radius = 40.0mm,min_fnumber = 1.6, pixel_pitch = .7μm) #these numbers give a reasonable system, albeit with 2000 lenslets. However displays are quite small (184x252)μm
end

"""System parameters which will yield fewer lenslets which will draw much faster"""
function smallsystemparameters()
    return (eye_box = (10.0mm,9.0mm),fov = (20.0°,20.0°),eye_relief = 20.0mm, pupil_diameter = 2mm, display_sphere_radius = 40.0mm,min_fnumber = 1.6, pixel_pitch = .7μm) #these numbers give a reasonable system, albeit with 2000 lenslets. However displays are quite small (184x252)μm
end
export smallsystemparameters

"""Coordinate frames for the eye/display system. The origin of this frame is at the geometric center of the eyeball. Positive Z axis is assumed to be the forward direction, the default direction of the eye's optical axis when looking directly ahead."""
function setup_coordinate_frames()
    eyeballframe = Transform()
    corneavertex = OpticSim.Data.cornea_to_eyecenter()
    eyeboxtransform = eyeballframe*OpticSim.translation(0.0,0.0,ustrip(mm,corneavertex))  #unfortunately can't use unitful values in transforms because the rotation and translation components would have different types which is not allowed in a Matrix.

    return (eyeball_frame = eyeballframe,eye_box_frame = eyeboxtransform)
end

"""repeat vec1 enough times to fill a vector of length(vec2)."""
function extend(vec1,vec2)
    numrepeats = Int64(ceil(length(vec1)/length(vec2)))
    remaining = mod(length(vec1),length(vec2))
    subdivs= similar(vec2,numrepeats*length(vec2)+remaining)
    empty!(subdivs)

    for i in 1:numrepeats-1
        append!(subdivs,vec2)
    end
    append!(subdivs,vec2[1:remaining])

    return subdivs
end

"""Create lenslet system that will cover the eyebox and fov. If the numbers in the fov tuple do not have units the system assumes they represent an angle in radians. Is is better to use unitful angles though, either rad or °:
```
setup_system(eye_box,(1rad,1.2rad),...)
setup_system(eye_box,(1°,1.2°),...)
```
If you do this
```
setup_system(eye_box,(1,1.2),...)
```
the system will assume the angle is in radians."""
function setup_system(eye_box,fov,eye_relief,pupil_diameter,display_sphere_radius,min_fnumber,pixel_pitch)
    #All coordinates are ultimately transformed into the eyeball_frame coordinate systems
    (eyeball_frame,eye_box_frame) = setup_coordinate_frames()

    @info "Input parameters\n"

    parameters = (eye_box = eye_box,fov = fov,eye_relief = eye_relief,pupil_diameter = pupil_diameter,display_sphere_radius = display_sphere_radius,min_fnumber = min_fnumber,pixel_pitch = pixel_pitch)

    for key in keys(parameters) @info "$key = $(parameters[key])" end
    println("\n\n")

    eyeboxz = (eye_box_frame*SVector(0.0,0.0,0.0))[3]
    eyebox_plane = Plane([0.0,0.0,1.0],[0.0,0.0,eyeboxz])
    #get system properties
    props = system_properties(eye_relief,eye_box,fov,pupil_diameter,.2,11,pixelpitch = pixel_pitch, minfnumber = min_fnumber)

    subdivisions = props[:subdivisions] #tuple representing how the eyebox can be subdivided given the cluster used for the lenslets

    clusterdata = props[:cluster_data]
    cluster = clusterdata[:cluster] #cluster that is repeated across the display to ensure continuous coverage of the eyebox and fov.
    focallength = ustrip(mm,props[:focal_length]) #strip units off because these don't work well with Transform

    #compute lenslets based on system properties. lattice_coordinates are the (i,j) integer lattice coordinates of the hexagonal lattice making up the display. These coordinates are used to properly assign color and subdivided eyebox to the lenslets.
    lenses,lattice_coordinates = spherelenslets(eyebox_plane,eye_relief,focallength,[0.0,0.0,1.0],display_sphere_radius,fov[1],fov[2],elementbasis(cluster))

    @info "$(length(lenses)) lenses required for this design"
    temp = display_plane.(lenses)
    displayplanes = [x[1] for x in temp]
    display_plane_centers = [x[2] for x in temp]

    lensletcolors = pointcolor.(lattice_coordinates,Ref(cluster))

    #compute subdivided eyebox polygons and assign to appropriate lenslets
    eyeboxpoly::SMatrix{3,4}  =  mm * (eye_box_frame * eyeboxpolygon(ustrip.(mm,eye_box)...)) #four corners of the eyebox frame which is assumed centered around the positive Z axis. Transformed to the eyeballframe. Have to switch back and forth between Unitful and unitless quantities because Transform doesn't work with Unitful values.

    @info "Subdividing eyebox polygon into $subdivisions sub boxes"
    subdivided_eyeboxpolys::Vector{SMatrix{3,4}} = compute_lenslet_eyebox_data(eye_box_frame,eyeboxpoly,subdivisions)
    test_compute_lenslet_eyebox_data(eye_box_frame,eyeboxpoly,subdivisions)

    #assign each lenslet one of the subdivided eyebox polygons
    lenslet_eye_boxes = eyebox_assignment.(lattice_coordinates,Ref(cluster),Ref(subdivided_eyeboxpolys))

    #assign each lenslet the index into subdivided_eyeboxpolys that corresponds to the subdivided eyebox polygon the lenslet is supposed to cover.
    lenslet_eyebox_numbers = eyebox_number.(lattice_coordinates,Ref(cluster),Ref(length(subdivided_eyeboxpolys)))

    #compute centroids of the subdivided eyeboxes on the eyebox plane
    lensleteyeboxcenters = [Statistics.mean(eachcol(x)) for x in lenslet_eye_boxes]

    offset_lenses = replace_optical_center.(lensleteyeboxcenters,display_plane_centers,lenses)  #make new lenses with optical centers of lenses to put display centers directly behind lenslets

    projected = project_eyebox_to_display_plane.(lenslet_eye_boxes,offset_lenses,displayplanes) #repeate subdivided_eyeboxpolys enough times to cover all lenses
    projected_eyeboxes = [x[1] for x in projected]
    projected_polygons = [x[2] for x in projected]
    @info "typeof polygon $(eltype(projected_polygons))"

    eyebox_rectangle = Rectangle(ustrip(mm,eye_box[1]/2),ustrip(mm,eye_box[2]/2),[0.0,0.0,1.0],[0.0,0.0,eyeboxz], interface = opaqueinterface())

    return LensletSystem{Float64}(eyebox_rectangle,subdivisions,offset_lenses,lattice_coordinates,lensletcolors,projected_eyeboxes,displayplanes,lenslet_eye_boxes,lenslet_eyebox_numbers,lensleteyeboxcenters,props,subdivided_eyeboxpolys,projected_polygons)

end
export setup_system


text:
struct EvalFunc{F} <: Function
  f::F
end
(f::EvalFunc)(args...) = f.f(args...)

NO_TSPAN_PROBS = Union{AbstractLinearProblem, AbstractNonlinearProblem,
                       AbstractQuadratureProblem,
                       AbstractSteadyStateProblem,AbstractJumpProblem}

has_kwargs(_prob::DEProblem) = has_kwargs(typeof(_prob))
Base.@pure has_kwargs(::Type{T}) where T = :kwargs ∈ fieldnames(T)

function init_call(_prob,args...;merge_callbacks = true,kwargs...)

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs,:callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
      values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}( (DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback ),) )
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if hasfield(typeof(_prob),:f) && hasfield(typeof(_prob.f),:f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__init,_prob,args...; kwargs...)#::T
  else
    __init(_prob,args...;kwargs...)#::T
  end
end

function init(prob::DEProblem,args...;kwargs...)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _prob = get_concrete_problem(prob,isadaptive(alg);kwargs...)
    init_call(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _prob = get_concrete_problem(prob,isadaptive(alg);kwargs...)
    init_call(_prob,args...;kwargs...)
  else
    _prob = get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);kwargs...)
    init_call(_prob,args...;kwargs...)
  end
end

function solve_call(_prob,args...;merge_callbacks = true, kwargs...)
  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs,:callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
      values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}( (DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback ),) )
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if hasfield(typeof(_prob),:f) && hasfield(typeof(_prob.f),:f) && typeof(_prob.f.f) <: EvalFunc
    Base.invokelatest(__solve,_prob,args...; kwargs...)#::T
  else
    __solve(_prob,args...;kwargs...)#::T
  end
end

# save_idxs and saveat are here due to https://github.com/FluxML/Zygote.jl/issues/664
function solve(prob::DEProblem,args...;sensealg=nothing,
               u0 = nothing, p = nothing, kwargs...)
  u0 = u0 !== nothing ? u0 : prob.u0
  p  = p  !== nothing ? p  : prob.p
  if sensealg === nothing && haskey(prob.kwargs,:sensealg)
    sensealg = prob.kwargs[:sensealg]
  end
  solve_up(prob,sensealg,u0,p,args...;kwargs...)
end

function solve_up(prob::DEProblem,sensealg,u0,p,args...;kwargs...)

  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    _alg = prepare_alg(alg,u0,p,prob)
    _prob = get_concrete_problem(prob,isadaptive(_alg);u0=u0,p=p,kwargs...)
    solve_call(_prob,_alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    _alg = prepare_alg(alg,u0,p,prob)
    _prob = get_concrete_problem(prob,isadaptive(_alg);u0=u0,p=p,kwargs...)
    solve_call(_prob,_alg,Base.tail(args)...;kwargs...)
  elseif isempty(args) # Default algorithm handling
    _prob = get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
    solve_call(_prob,args...;kwargs...)
  else
    _prob = get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
    solve_call(_prob,args...;kwargs...)
  end
end

function solve(prob::EnsembleProblem,args...;kwargs...)
  if isempty(args) || length(args) == 1 && typeof(args[1])<:EnsembleAlgorithm
    __solve(prob,nothing,args...;kwargs...)
  else
    __solve(prob,args...;kwargs...)
  end
end

function solve(prob::AbstractNoiseProblem,args...; kwargs...)
  __solve(prob,args...;kwargs...)
end

function solve(prob::AbstractJumpProblem,args...; kwargs...)
  __solve(prob,args...;kwargs...)
end

function get_concrete_problem(prob::AbstractJumpProblem, isadapt; kwargs...)
  prob
end

function get_concrete_problem(prob::SteadyStateProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, Inf, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0 = u0)
end

function get_concrete_problem(prob::NonlinearProblem, isadapt; kwargs...)
  u0 = get_concrete_u0(prob, isadapt, nothing, kwargs)
  u0 = promote_u0(u0, prob.p, nothing)
  remake(prob; u0 = u0)
end

function get_concrete_problem(prob::AbstractEnsembleProblem, isadapt; kwargs...)
  prob
end

function solve(prob::PDEProblem,alg::DiffEqBase.DEAlgorithm,args...;
                                          kwargs...)
    solve(prob.prob,alg,args...;kwargs...)
end

function init(prob::PDEProblem,alg::DiffEqBase.DEAlgorithm,args...;
                                          kwargs...)
    init(prob.prob,alg,args...;kwargs...)
end

function get_concrete_problem(prob, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  u0_promote = promote_u0(u0, p, tspan[1])
  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
                  prob.tspan == tspan && typeof(tspan) === typeof(tspan_promote) &&
                  p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f = f_promote, u0 = u0_promote, p=p, tspan = tspan_promote)
  end
end

function get_concrete_problem(prob::DAEProblem, isadapt; kwargs...)

  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
  du0 = get_concrete_du0(prob, isadapt, tspan[1], kwargs)

  u0_promote = promote_u0(u0, p, tspan[1])
  du0_promote = promote_u0(du0, p, tspan[1])

  f_promote = promote_f(prob.f, u0_promote)
  tspan_promote = promote_tspan(u0_promote, p, tspan, prob, kwargs)
  if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
     isconcretedu0(prob, tspan[1], kwargs) && typeof(du0_promote) === typeof(prob.du0) &&
                  prob.tspan == tspan && typeof(tspan) === typeof(tspan_promote) &&
                  p === prob.p && f_promote === prob.f
    return prob
  else
    return remake(prob; f = f_promote, du0 = du0_promote, u0 = u0_promote, p=p, tspan = tspan_promote)
  end
end

function get_concrete_problem(prob::DDEProblem, isadapt; kwargs...)
  p = get_concrete_p(prob, kwargs)
  tspan = get_concrete_tspan(prob, isadapt, kwargs, p)
  u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)

  if prob.constant_lags isa Function
    constant_lags = prob.constant_lags(p)
  else
    constant_lags = prob.constant_lags
  end

  u0 = promote_u0(u0, p, tspan[1])
  tspan = promote_tspan(u0, p, tspan, prob, kwargs)

  remake(prob; u0 = u0, tspan = tspan, p=p, constant_lags = constant_lags)
end

function promote_f(f,u0)
    # Ensure our jacobian will be of the same type as u0
    uElType = u0 === nothing ? Float64 : eltype(u0)
    if isdefined(f, :jac_prototype) && f.jac_prototype isa AbstractArray
        f = @set f.jac_prototype = similar(f.jac_prototype, uElType)
    end
    return f
end

promote_f(f::SplitFunction,u0) = typeof(f.cache) === typeof(u0) && isinplace(f) ? f : remake(f,cache=zero(u0))
prepare_alg(alg,u0,p,f) = alg

function get_concrete_tspan(prob, isadapt, kwargs, p)
  if prob.tspan isa Function
    tspan = prob.tspan(p)
  elseif haskey(kwargs, :tspan)
      tspan = kwargs[:tspan]
  elseif prob.tspan === (nothing, nothing)
    error("No tspan is set in the problem or chosen in the init/solve call")
  else
    tspan = prob.tspan
  end

  isadapt && eltype(tspan) <: Integer && (tspan = float.(tspan))

  tspan
end

function isconcreteu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.u0 !== nothing && !isdistribution(prob.u0)
end

function isconcretedu0(prob, t0, kwargs)
  !eval_u0(prob.u0) && prob.du0 !== nothing && !isdistribution(prob.du0)
end

function get_concrete_u0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif haskey(kwargs,:u0)
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  isadapt && eltype(u0) <: Integer && (u0 = float.(u0))

  handle_distribution_u0(u0)
end

function get_concrete_du0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.du0)
    du0 = prob.du0(prob.p, t0)
  elseif haskey(kwargs,:du0)
    du0 = kwargs[:du0]
  else
    du0 = prob.du0
  end

  isadapt && eltype(du0) <: Integer && (du0 = float.(du0))

  handle_distribution_u0(du0)
end

function get_concrete_p(prob, kwargs)
  if haskey(kwargs,:p)
    p = kwargs[:p]
  else
    p = prob.p
  end
end

handle_distribution_u0(_u0) = _u0
handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
isdistribution(_u0::Distributions.Sampleable) = true

eval_u0(u0::Function) = true
eval_u0(u0) = false

function __solve(prob::DEProblem,args...;default_set=false,second_time=false,kwargs...)
  if second_time
    error("""
      Default algorithm choices require DifferentialEquations.jl.
      Please specify an algorithm (e.g., `solve(prob, Tsit5())` for an ODE)
      or import DifferentialEquations directly.

      You can find the list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/
      and its associated pages.
      """
    )
  elseif length(args) > 0 && !(typeof(args[1]) <: Union{Nothing,DEAlgorithm})
    error("Inappropiate solve command. The arguments do not make sense. Likely, you gave an algorithm which does not actually exist (or does not `<:DiffEqBase.DEAlgorithm`)")
  else
    __solve(prob::DEProblem,nothing,args...;default_set=false,second_time=true,kwargs...)
  end
end

################### Differentiation

struct SensitivityADPassThrough <: SciMLBase.DEAlgorithm end

function ChainRulesCore.frule(::typeof(solve_up),prob,
                              sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
                              u0,p,args...;
                              kwargs...)
  _solve_forward(prob,sensealg,u0,p,args...;kwargs...)
end

function ChainRulesCore.rrule(::typeof(solve_up),prob::SciMLBase.DEProblem,
                              sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
                              u0,p,args...;
                              kwargs...)
  _solve_adjoint(prob,sensealg,u0,p,args...;kwargs...)
end

###
### Legacy Dispatches to be Non-Breaking
###

@deprecate concrete_solve(prob::SciMLBase.DEProblem,alg::Union{SciMLBase.DEAlgorithm,Nothing},
                        u0=prob.u0,p=prob.p,args...;kwargs...) solve(prob,alg,args...;u0=u0,p=p,kwargs...)

function _solve_adjoint(prob,sensealg,u0,p,args...;merge_callbacks = true, kwargs...)

  _prob = if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob,isadaptive(alg);u0=u0,p=p,kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob,isadaptive(alg);u0=u0,p=p,kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
  else
    get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs,:callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
      values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}( (DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback ),) )
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_adjoint(_prob,nothing,sensealg,u0,p;kwargs...)
  else
    _concrete_solve_adjoint(_prob,args[1],sensealg,u0,p,Base.tail(args)...;kwargs...)
  end
end

function _solve_forward(prob,sensealg,u0,p,args...;merge_callbacks = true, kwargs...)

  _prob = if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    get_concrete_problem(prob,isadaptive(alg);u0=u0,p=p,kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    get_concrete_problem(prob,isadaptive(alg);u0=u0,p=p,kwargs...)
  elseif isempty(args) # Default algorithm handling
    get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
  else
    get_concrete_problem(prob,!(typeof(prob)<:DiscreteProblem);u0=u0,p=p,kwargs...)
  end

  if has_kwargs(_prob)
    if merge_callbacks && haskey(_prob.kwargs,:callback) && haskey(kwargs, :callback)
      kwargs_temp = NamedTuple{Base.diff_names(Base._nt_names(
      values(kwargs)), (:callback,))}(values(kwargs))
      callbacks = NamedTuple{(:callback,)}( (DiffEqBase.CallbackSet(_prob.kwargs[:callback], values(kwargs).callback ),) )
      kwargs = merge(kwargs_temp, callbacks)
    end
    kwargs = isempty(_prob.kwargs) ? kwargs : merge(values(_prob.kwargs), kwargs)
  end

  if isempty(args)
    _concrete_solve_forward(prob,nothing,sensealg,u0,p;kwargs...)
  else
    _concrete_solve_forward(prob,args[1],sensealg,u0,p,Base.tail(args)...;kwargs...)
  end
end

function _concrete_solve_adjoint(args...;kwargs...)
  error("No adjoint rules exist. Check that you added `using DiffEqSensitivity`")
end

function _concrete_solve_forward(args...;kwargs...)
  error("No sensitivity rules exist. Check that you added `using DiffEqSensitivity`")
end

ext: jl


text:
using Colors, Distributions, LaTeXStrings, PGFPlotsX

lprob(d, x) = cdf(d, x) - cdf(d, x-1/1000)
xs = map(x -> x/1000, 1:10000)

function seir_var_susp(r0, d, g, k, ep, days, conn)
    gamma_d = Gamma(k, 1/k)
    lprobs = map((x) -> lprob(gamma_d, x), xs[1:9999])
    push!(lprobs, 1 - sum(lprobs))
    st = [(1-ep) .* lprobs]
    et = [ep .* lprobs]
    it = [fill(0.0, 10000)]
    rt = [fill(0.0, 10000)]
    if conn
        b = r0/(1+1/k) * g
    else
        b = r0 * g
    end
    for t in 2:days
        if conn
            l = b * sum(xs .* it[t-1])
        else
            l = b * sum(it[t-1])
        end
        sred = xs .* l .* st[t-1]
        ered = d .* et[t-1]
        ired = g .* it[t-1]
        push!(st, st[t-1] .- sred)
        push!(et, et[t-1] .+ sred .- ered)
        push!(it, it[t-1] .+ ered .- ired)
        push!(rt, rt[t-1] .+ ired)
    end
    Dict(:st => st, :et => et, :it => it, :rt => rt,
        :r0 => r0, :d => d, :g => g, :k => k,
        :ep => ep, :conn => conn)
end

sedec(x) = replace("$x", "."=>"{,}")

function seirset(r0, d, g, kexps, ep, days)
    seirs_s = map(kexp -> seir_var_susp(r0, d, g, 2.0^kexp, ep, days, false), kexps)
    seirs_c = map(kexp -> seir_var_susp(r0, d, g, 2.0^kexp, ep, days, true), kexps)
    Dict(:seirs_s => seirs_s, :seirs_c => seirs_c, :kexps => kexps)
end

seirset_cends(seirset, c) = map(r -> r[end], map(s -> map(sum, s[c]), seirset))

function seirplot(seir)
    r0lab = latexstring(*("\\mathcal{R}_0=", sedec(seir[:r0])))
    klab = latexstring(*("k=", sedec(seir[:k])))
    if seir[:conn]
        connlab = "kontakt"
    else
        connlab = "mottaglighet"
    end
    days = 1:length(seir[:st])
    p = @pgf Axis({"yticklabel style"={"/pgf/number format/use comma"},
        legend_pos="outer north east",
        title="SEIR för $(r0lab), $(klab), varierande $(connlab)", xlabel=L"t",
        ylabel=L"p",
        xmajorgrids, ymajorgrids})
    cs = [("S", seir[:st]), ("E", seir[:et]), ("I", seir[:it]), ("R", seir[:rt])]
    plotcolors = distinguishable_colors(length(cs)+1, [RGB(1,1,1)])[2:end]
    for (i, c) in enumerate(cs)
        @pgf push!(p, PlotInc({no_markers, color = plotcolors[i]},
            Table([days, map(sum, c[2])])), LegendEntry(c[1]))
    end
    p
end

function seirset_rendsplot(seirset)
    seir = seirset[:seirs_s][1]
    r0lab = latexstring(*("\\mathcal{R}_0=", sedec(seir[:r0])))
    rendlab = latexstring("\\mathrm{R}($(length(seir[:st])))")
    p = @pgf Axis({"yticklabel style"={"/pgf/number format/use comma"},
        legend_pos="outer north east",
        title="SEIR för $(r0lab)", xlabel=L"\lg(k)",
        ylabel=rendlab,
        xmajorgrids, ymajorgrids})
    @pgf push!(p, PlotInc(Table([seirset[:kexps], seirset_cends(seirset[:seirs_s], :rt)])),
        LegendEntry("mottaglighet"))
    @pgf push!(p, PlotInc(Table([seirset[:kexps], seirset_cends(seirset[:seirs_c], :rt)])),
        LegendEntry("kontakt"))
    p
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# A Caesar cipher implementation in Julia

function caesar()
    println("Please enter the text: ")
    text = readline()
    println("Please enter key")
    key = parse(Int, readline())
    for char in text
        newchar = char + key
        if(newchar > 'z')
            newchar -= 26
        end
        if(newchar < 'a')
            newchar += 26
        end
        print(newchar)
    end
end
caesar()


text:
using Plots, FluxRC
using FastGaussQuadrature

"""
N polynomial degree -> Np solution points
"""
N = 2
Np = (N + 1) * (N + 2) ÷ 2
points, weights = tri_quadrature(N)
scatter(points[:, 1], points[:, 2], ratio = 1)

"""
Vandermonde matrix bridges Lagrange polynomials and modal polynomials
Vᵀℓ(r) = P(r)
"""
V = vandermonde_matrix(N, points[:, 1], points[:, 2])
Vr, Vs = ∂vandermonde_matrix(N, points[:, 1], points[:, 2]) # (r_i, ψ_j)

# Lagrange polynomials
∂l = zeros(Np, Np, 2)
for i = 1:Np
    ∂l[i, :, 1] .= V' \ Vr[i, :]
    ∂l[i, :, 2] .= V' \ Vs[i, :]
end

# test solution (derivatives at solution point)
u = [3.0, 2.0, 3.0, 2.0, 1.0, 2.0]
uhat = V \ u

uhat .* Vr[1, :] |> sum
u .* ∂l[1, :, 1] |> sum
uhat .* Vs[1, :] |> sum
u .* ∂l[1, :, 2] |> sum

# interface interpolation
function triface_quadrature(N)
    Δf = [1.0, √2, 1.0]

    pf = Array{Float64}(undef, 3, N + 1, 2)
    wf = Array{Float64}(undef, 3, N + 1)

    p0, w0 = gausslegendre(N + 1)

    pf[1, :, 1] .= p0
    pf[2, :, 1] .= p0
    pf[3, :, 1] .= -1.0
    pf[1, :, 2] .= -1.0
    pf[2, :, 2] .= -p0
    pf[3, :, 2] .= p0

    wf[1, :] .= w0 .* Δf[1]
    wf[2, :] .= w0 .* Δf[2]
    wf[3, :] .= w0 .* Δf[3]

    return pf, wf
end

pf, wf = triface_quadrature(N)

ψf = zeros(3, N + 1, Np)
for i = 1:3
    ψf[i, :, :] .= vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

∂ψf = zeros(3, N + 1, Np, 2)
for i = 1:3
    ∂ψf[i, :, :, 1], ∂ψf[i, :, :, 2] = ∂vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

lf = zeros(3, N + 1, Np)
for i = 1:3, j = 1:N+1
    lf[i, j, :] .= V' \ ψf[i, j, :]
end

u .* lf[2, 2, :] |> sum
uhat .* ψf[2, 2, :] |> sum

∂lf = zeros(3, N + 1, Np, 2)
for i = 1:3, j = 1:N+1, k = 1:2
    ∂lf[i, j, :, k] .= V' \ ∂ψf[i, j, :, k]
end

uhat .* ∂ψf[1, 2, :, 1] |> sum
u .* ∂lf[1, 2, :, 1] |> sum

# correction field
ϕ = zeros(3, N + 1)
σ = zeros(3, N + 1, Np)

for k = 1:Np
    for j = 1:N+1
        for i = 1:3
            σ[i, j, k] = wf[i, j] * ψf[i, j, k] * Δf[i]
        end
    end
end

for j = 1:N+1
    for i = 1:3
        ϕ[i, j] = sum(σ[i, j, :] .* ψf[i, j, :])
    end
end


text:
const CKSP = Ptr{Cvoid}
const CKSPType = Cstring

abstract type AbstractKSP{T, PetscLib} <: Factorization{T} end

Base.@kwdef mutable struct KSP{T, PetscLib} <: AbstractKSP{T, PetscLib}
    ptr::CKSP = C_NULL
    opts::Options{PetscLib}
    # Stuff to keep around so that they don't get gc'ed
    _A = nothing
    _P = nothing
    _dm = nothing
    # Function pointers
    ComputeRHS! = nothing
    ComputeOperators! = nothing
end

struct WrappedKSP{T, PetscLib} <: AbstractKSP{T, PetscLib}
    ptr::CKSP
end

scalartype(::KSP{T}) where {T} = T
Base.eltype(::KSP{T}) where {T} = T

LinearAlgebra.transpose(ksp) = LinearAlgebra.Transpose(ksp)
LinearAlgebra.adjoint(ksp) = LinearAlgebra.Adjoint(ksp)

"""
    KSPSetComputeRHS!(
        ksp::KSP{Number},
        ComputeRHS!,
        ctx = C_NULL,
    )

Set the right-hand side function `ComputeRHS!` for the `ksp` using the user
`ctx`. `ComputeRHS!` should be callable with three arguments of type
`(::KSP{Number}, ::Vec, ::Ptr)`;
see [PETSc manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetComputeRHS.html)
"""
function KSPSetComputeRHS! end

"""
    struct Fn_KSPComputeRHS{T} end

Type used to wrap `ComputeRHS!` functions in KSP
"""
struct Fn_KSPComputeRHS{T} end

"""
    KSPSetComputeOperators!(
        ksp::KSP{Number},
        ComputeOperators!,
        ctx = C_NULL,
    )

Set the linear operators function `ComputeOperators!` for the `ksp` using the
user `ctx`. `ComputeOperators!` should be callable with four arguments of type
`(::KSP{Number}, ::Mat, ::Mat, ::Ptr)`;
see [PETSc manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetComputeOperators.html)
"""
function KSPSetComputeOperators! end

"""
    struct Fn_KSPComputeOperators{T} end

Type used to wrap `ComputeOperators!` functions in KSP
"""
struct Fn_KSPComputeOperators{T} end

@for_libpetsc begin

    function KSP{$PetscScalar}(comm::MPI.Comm; kwargs...)
        @assert initialized($petsclib)
        opts = Options($petsclib, kwargs...)
        ksp = KSP{$PetscScalar, $PetscLib}(opts=opts)
        with(ksp.opts) do
          @chk ccall((:KSPCreate, $libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CKSP}), comm, ksp)
        end
        if comm == MPI.COMM_SELF
            finalizer(destroy, ksp)
        end
        return ksp
    end

    function destroy(ksp::KSP{$PetscScalar})
        finalized($petsclib) ||
        @chk ccall((:KSPDestroy, $libpetsc), PetscErrorCode, (Ptr{CKSP},), ksp)
        return nothing
    end

    function setoperators!(ksp::KSP{$PetscScalar}, A::AbstractMat{$PetscScalar}, P::AbstractMat{$PetscScalar})
        @chk ccall((:KSPSetOperators, $libpetsc), PetscErrorCode, (CKSP, CMat, CMat), ksp, A, P)
        ksp._A = A
        ksp._P = P
        return nothing
    end

    function (::Fn_KSPComputeRHS{$PetscScalar})(
        new_ksp_ptr::CKSP,
        cb::CVec,
        ksp_ptr::Ptr{Cvoid}
       )::$PetscInt
        ksp = unsafe_pointer_to_objref(ksp_ptr)
        new_ksp = WrappedKSP{$PetscScalar, $PetscLib}(new_ksp_ptr)
        b = Vec{$PetscScalar}(cb)
        ierr = ksp.ComputeRHS!(new_ksp, b)
        return $PetscInt(ierr)
    end

    function KSPSetComputeRHS!(
        ksp::KSP{$PetscScalar},
        ComputeRHS!
    )
        ptr_ksp = pointer_from_objref(ksp)
        fptr = @cfunction(Fn_KSPComputeRHS{$PetscScalar}(),
                          $PetscInt,
                          (CKSP, CVec, Ptr{Cvoid}))
        with(ksp.opts) do
            @chk ccall((:KSPSetComputeRHS, $libpetsc), PetscErrorCode,
                (CKSP, Ptr{Cvoid}, Ptr{Cvoid}),
                ksp, fptr, ptr_ksp)
        end
        ksp.ComputeRHS! = ComputeRHS!
        return ksp
    end

    function (::Fn_KSPComputeOperators{$PetscScalar})(
        new_ksp_ptr::CKSP,
        cA::CMat,
        cP::CMat,
        ksp_ptr::Ptr{Cvoid}
       )::$PetscInt
        ksp = unsafe_pointer_to_objref(ksp_ptr)
        new_ksp = WrappedKSP{$PetscScalar, $PetscLib}(new_ksp_ptr)
        A = Mat{$PetscScalar}(cA)
        P = Mat{$PetscScalar}(cP)
        ierr = ksp.ComputeOperators!(new_ksp, A, P)
        return $PetscInt(ierr)
    end

    function KSPSetComputeOperators!(
        ksp::KSP{$PetscScalar},
        ComputeOperators!
    )
        ptr_ksp = pointer_from_objref(ksp)
        fptr = @cfunction(Fn_KSPComputeOperators{$PetscScalar}(),
                          $PetscInt,
                          (CKSP, CMat, CMat, Ptr{Cvoid}))
        with(ksp.opts) do
            @chk ccall((:KSPSetComputeOperators, $libpetsc), PetscErrorCode,
                (CKSP, Ptr{Cvoid}, Ptr{Cvoid}),
                ksp, fptr, ptr_ksp)
        end
        ksp.ComputeOperators! = ComputeOperators!
        return ksp
    end

    function KSPSetDM!(ksp::KSP{$PetscScalar}, dm::AbstractDM{$PetscLib})
        with(ksp.opts) do
            @chk ccall((:KSPSetDM, $libpetsc), PetscErrorCode, (CKSP, CDM), ksp, dm)
        end
        ksp._dm = dm
        return nothing
    end

    function DMDA(ksp::AbstractKSP{$PetscScalar})
        t_dm = Ref{CDM}()
        @chk ccall(
            (:KSPGetDM, $libpetsc),
            PetscErrorCode,
            (CKSP, Ptr{CDM}),
            ksp,
            t_dm,
        )
        dm = DMDA{$PetscLib}(t_dm[])
        return dm
    end

    function settolerances!(ksp::KSP{$PetscScalar}; rtol=PETSC_DEFAULT, atol=PETSC_DEFAULT, divtol=PETSC_DEFAULT, max_it=PETSC_DEFAULT)
        @chk ccall((:KSPSetTolerances, $libpetsc), PetscErrorCode,
                    (CKSP, $PetscReal, $PetscReal, $PetscReal, $PetscInt),
                    ksp, rtol, atol, divtol, max_it)
        return nothing
    end

    function setfromoptions!(ksp::KSP{$PetscScalar})
        with(ksp.opts) do
            @chk ccall((:KSPSetFromOptions, $libpetsc), PetscErrorCode, (CKSP,), ksp)
        end
    end

    function gettype(ksp::KSP{$PetscScalar})
        t_r = Ref{CKSPType}()
        @chk ccall((:KSPGetType, $libpetsc), PetscErrorCode, (CKSP, Ptr{CKSPType}), ksp, t_r)
        return unsafe_string(t_r[])
    end

    function iters(ksp::KSP{$PetscScalar})
        r_its = Ref{$PetscInt}()
        @chk ccall((:KSPGetIterationNumber, $libpetsc), PetscErrorCode,
        (KSP, Ptr{$PetscInt}), ksp, r_its)
        return r_its[]
    end

    function view(ksp::KSP{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, getcomm(ksp)))
        @chk ccall((:KSPView, $libpetsc), PetscErrorCode,
                    (CKSP, CPetscViewer),
                ksp, viewer);
        return nothing
    end

    function resnorm(ksp::KSP{$PetscScalar})
        r_rnorm = Ref{$PetscReal}()
        @chk ccall((:KSPGetResidualNorm, $libpetsc), PetscErrorCode,
        (KSP, Ptr{$PetscReal}), ksp, r_rnorm)
        return r_rnorm[]
    end

    function solve!(x::AbstractVec{$PetscScalar}, ksp::KSP{$PetscScalar}, b::AbstractVec{$PetscScalar})
        with(ksp.opts) do
            @chk ccall((:KSPSolve, $libpetsc), PetscErrorCode,
            (CKSP, CVec, CVec), ksp, b, x)
        end
        return x
    end

    function solve!(ksp::KSP{$PetscScalar})
        with(ksp.opts) do
            @chk ccall((:KSPSolve, $libpetsc), PetscErrorCode,
            (CKSP, CVec, CVec), ksp, C_NULL, C_NULL)
        end
        return nothing
    end

    function solve!(x::AbstractVec{$PetscScalar}, tksp::Transpose{T,K}, b::AbstractVec{$PetscScalar}) where {T,K <: KSP{$PetscScalar}}
        ksp = parent(tksp)
        with(ksp.opts) do
            @chk ccall((:KSPSolveTranspose, $libpetsc), PetscErrorCode,
            (CKSP, CVec, CVec), ksp, b, x)
        end
        return x
    end

end

# no generic Adjoint solve defined, but for Real we can use Adjoint
solve!(x::AbstractVec{T}, aksp::Adjoint{T,K}, b::AbstractVec{T}) where {K <: KSP{T}} where {T<:Real} =
    solve!(x, transpose(parent(aksp)), b)

const KSPAT{T, LT} = Union{KSP{T, LT}, Transpose{T, KSP{T, LT}}, Adjoint{T, KSP{T, LT}}}

LinearAlgebra.ldiv!(x::AbstractVec{T}, ksp::KSPAT{T, LT}, b::AbstractVec{T}) where {T, LT} = solve!(x, ksp, b)
function LinearAlgebra.ldiv!(x::AbstractVector{T}, ksp::KSPAT{T, LT}, b::AbstractVector{T}) where {T, LT}
    parent(solve!(AbstractVec(x), ksp, AbstractVec(b)))
end
Base.:\(ksp::KSPAT{T, LT}, b::AbstractVector{T}) where {T, LT} = ldiv!(similar(b), ksp, b)

"""
    KSP(A, P; options...)

Construct a PETSc Krylov subspace solver.

Any PETSc options prefixed with `ksp_` and `pc_` can be passed as keywords.
"""
function KSP(A::AbstractMat{T}, P::AbstractMat{T}=A; kwargs...) where {T}
    ksp = KSP{T}(getcomm(A); kwargs...)
    setoperators!(ksp, A, P)
    setfromoptions!(ksp)
    return ksp
end

"""
    KSP(da::AbstractDM; options...)

Construct a PETSc Krylov subspace solver from the distributed mesh

Any PETSc options prefixed with `ksp_` and `pc_` can be passed as keywords.

see [PETSc manual](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetDM.html)
"""
function KSP(dm::AbstractDM{PetscLib}; kwargs...) where {PetscLib}
    T = scalartype(PetscLib)
    ksp = KSP{T}(getcomm(dm); kwargs...)
    KSPSetDM!(ksp, dm)
    setfromoptions!(ksp)
    return ksp
end

Base.show(io::IO, ksp::KSP) = _show(io, ksp)

"""
    iters(ksp::KSP)

Gets the current iteration number; if the `solve!` is complete, returns the number of iterations used.

# External Links
$(_doc_external("KSP/KSPGetIterationNumber"))
"""
iters

"""
    resnorm(ksp::KSP)

Gets the last (approximate preconditioned) residual norm that has been computed.

# External Links
$(_doc_external("KSP/KSPGetResidualNorm"))
"""
resnorm


text:
module ObservableFunctions

using PyCall

import ..@pyfunc

#qulacs = PyNULL()
observable = PyNULL()

const observable_functions = [
    :create_observable_from_openfermion_file,
    :create_observable_from_openfermion_text,
    :create_split_observable,
]

#=
for func in observable_functions
    @eval begin
        function $(func)(args...; kwargs...)
            qulacs.observable.$(func)(args...; kwargs...)
        end
        export $(func)
    end
end
=#

for func in observable_functions
    @eval begin
        @pyfunc observable $(func)
        export $(func)
    end
end

function __init__()
    #copy!(qulacs, pyimport("qulacs"))
    copy!(observable, pyimport("qulacs.observable"))
end

end # module ObservableFunctions


text:
abstract type AbstractBoolDomain <: AbstractDomain end

"""
    struct BoolDomain <: AbstractDomain

Boolean domain, uses a IntDomain in it. (true is 1 and false is 0)
"""

struct BoolDomain <: AbstractBoolDomain
    inner::IntDomain

    function BoolDomain(trailer::Trailer)
        return new(IntDomain(trailer, 2, -1))
    end
end

"""
    reset_domain!(dom::BoolDomain)

Used in `reset_model!`.
"""
reset_domain!(dom::BoolDomain) = reset_domain!(dom.inner)

function Base.show(io::IO, dom::BoolDomain)
    print(io, "[", join(dom, " "), "]")
end

function Base.show(io::IO, ::MIME"text/plain", dom::BoolDomain)
    print(io, typeof(dom), ": [", join(dom, " "), "]")
end

"""
    isempty(dom::BoolDomain)

Return `true` iff `dom` is an empty set. Done in constant time.
"""
Base.isempty(dom::BoolDomain) = Base.isempty(dom.inner)

"""
    length(dom::BoolDomain)

Return the size of `dom`. Done in constant time.
"""
Base.length(dom::SeaPearl.BoolDomain) = Base.length(dom.inner)

"""
    Base.in(value::Int, dom::BoolDomain)

Check if an integer is in the domain. Done in constant time.
"""
function Base.in(value::Bool, dom::BoolDomain)
    intValue = convert(Int, value)
    return Base.in(intValue, dom.inner)
end

"""
    remove!(dom::BoolDomain, value::Int)

Remove `value` from `dom`. Done in constant time.
"""
function remove!(dom::BoolDomain, value::Bool)
    if !(value in dom)
        return Bool[]
    end

    intValue = convert(Int, value)

    remove!(dom.inner, intValue)
    return [value]
end

"""
    removeAll!(dom::BoolDomain)

Remove every value from `dom`. Return the removed values. Done in constant time.
"""
removeAll!(dom::BoolDomain) = convert.(Bool, removeAll!(dom.inner))

"""
    assign!(dom::BoolDomain, value::Int)

Remove everything from the domain but `value`. Return the removed values. Return the pruned values.
Done in *constant* time.
"""
function assign!(dom::BoolDomain, value::Bool)
    @assert value in dom

    return convert.(Bool, assign!(dom.inner, convert(Int, value)))
end

"""
    Base.iterate(dom::BoolDomain, state=1)

Iterate over the domain in an efficient way. The order may not be consistent.
WARNING: Do **NOT** update the domain you are iterating on.
"""
function Base.iterate(dom::BoolDomain, state=1)
    returned = iterate(dom.inner, state)
    if isnothing(returned)
        return nothing
    end

    value, newState = returned
    return convert(Bool, value), newState
end


text:
@reset
@testset "CI" begin
    wf = @energy casci
    @test isapprox(wf.energy,-75.01775104333774,rtol=tol)
    @set ci_alg sparse
    wf = @energy casci
    @test isapprox(wf.energy,-75.0188434735,rtol=tol)
end


text:
abstract type AbstractObject end
abstract type AbstractCustomObject <: AbstractObject end


text:
"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::SequentialResampling)

Resample `x` according to a sequential resampling constraint.

This way of resampling introduces some serial dependence between the elements of `x` -
*beyond* what might already be present in the dataset. This is because imposing a
sequential constraint (e.g. `StrictlyIncreasing`) to the `i`-th value of the dataset
imposes constraints on what is possible to sample from the `i+1`th value.

## Example

```julia
# Some example data
N = 50
x_uncertain = [UncertainValue(Normal, x, rand(Uniform(0.1, 0.8))) for x in rand(N)]
y_uncertain = [UncertainValue(Normal, y, rand(Uniform(0.1, 0.8))) for y in rand(N)]
x = UncertainValueDataset(x_uncertain)
y = UncertainValueDataset(y_uncertain)

time_uncertain = [UncertainValue(Normal, i, 1) for i = 1:length(x)];
time_certain = [CertainValue(i) for i = 1:length(x)];
timeinds_x = UncertainIndexDataset(time_uncertain)
timeinds_y = UncertainIndexDataset(time_certain)

X = UncertainIndexValueDataset(timeinds_x, x)
Y = UncertainIndexValueDataset(timeinds_y, y);

# Resample
seq_resampling = SequentialResampling(StrictlyIncreasing())
resample(X, seq_resampling)
```
"""
function resample(x::AbstractUncertainIndexValueDataset,
        resampling::SequentialResampling{S}) where {S}
    resample(x, resampling.sequential_constraint)
end

"""
    resample(x::AbstractUncertainIndexValueDataset, resampling::SequentialInterpolatedResampling)

Resample `x` according to a sequential resampling constraint, then interpolate the draw(s) to
some specified grid.

This way of resampling introduces some serial dependence between the elements of `x` -
*beyond* what might already be present in the dataset. This is because imposing a
sequential constraint (e.g. `StrictlyIncreasing`) to the `i`-th value of the dataset
imposes constraints on what is possible to sample from the `i+1`th value.

## Example

```julia
# Some example data
N = 50
x_uncertain = [UncertainValue(Normal, x, rand(Uniform(0.1, 0.8))) for x in rand(N)]
y_uncertain = [UncertainValue(Normal, y, rand(Uniform(0.1, 0.8))) for y in rand(N)]
x = UncertainValueDataset(x_uncertain)
y = UncertainValueDataset(y_uncertain)

time_uncertain = [UncertainValue(Normal, i, 1) for i = 1:length(x)];
time_certain = [CertainValue(i) for i = 1:length(x)];
timeinds_x = UncertainIndexDataset(time_uncertain)
timeinds_y = UncertainIndexDataset(time_certain)

X = UncertainIndexValueDataset(timeinds_x, x)
Y = UncertainIndexValueDataset(timeinds_y, y);

# Resample
seqintp_resampling = SequentialInterpolatedResampling(StrictlyIncreasing(), RegularGrid(0:2:N))
resample(X, seqintp_resampling)
```
"""
function resample(x::AbstractUncertainIndexValueDataset,
        resampling::SequentialInterpolatedResampling{S, G}) where {S, G}
    resample(x, resampling.sequential_constraint, resampling.grid)
end


text:
# Options structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: May 2017
#

export Options, subsample

# Object for velocity/slowness models
mutable struct Options
    space_order::Integer
    free_surface::Bool
    limit_m::Bool
    buffer_size::Real
    save_data_to_disk::Bool
    file_path::String
    file_name::String
    sum_padding::Bool
    optimal_checkpointing::Bool
    num_checkpoints::Union{Integer, Nothing}
    checkpoints_maxmem::Union{Real, Nothing}
    frequencies::Array
    isic::Bool
    subsampling_factor::Integer
    dft_subsampling_factor::Integer
    return_array::Bool
    dt_comp::Union{Real, Nothing}
end

"""
    Options
        space_order::Integer
        free_surface::Bool
        limit_m::Bool
        buffer_size::Real
        save_rate::Real
        save_data_to_disk::Bool
        file_path::String
        file_name::String
        sum_padding::Bool
        optimal_checkpointing::Bool
        num_checkpoints::Integer
        checkpoints_maxmem::Real
	    frequencies::Array
	    isic::Bool
        subsampling_factor::Integer
	    dft_subsampling_factor::Integer
        return_array::Bool
        dt_comp::Real

Options structure for seismic modeling.

`space_order`: finite difference space order for wave equation (default is 8, needs to be multiple of 4)

`free_surface`: set to `true` to enable a free surface boundary condition.

`limit_m`: for 3D modeling, limit modeling domain to area with receivers (saves memory)

`buffer_size`: if `limit_m=true`, define buffer area on each side of modeling domain (in meters)

`save_data_to_disk`: if `true`, saves shot records as separate SEG-Y files

`file_path`: path to directory where data is saved

`file_name`: shot records will be saved as specified file name plus its source coordinates

`sum_padding`: when removing the padding area of the gradient, sum into boundary rows/columns for true adjoints

`optimal_checkpointing`: instead of saving the forward wavefield, recompute it using optimal checkpointing

`num_checkpoints`: number of checkpoints. If not supplied, is set to log(num_timesteps).

`checkpoints_maxmem`: maximum amount of memory that can be allocated for checkpoints (MB)

`frequencies`: calculate the FWI/LS-RTM gradient in the frequency domain for a given set of frequencies

`subsampling_factor`: compute forward wavefield on a time axis that is reduced by a given factor (default is 1)

`dft_subsampling_factor`: compute on-the-fly DFTs on a time axis that is reduced by a given factor (default is 1)

`isic`: use linearized inverse scattering imaging condition

`return_array`: return data from nonlinear/linear modeling as a plain Julia array.

`dt_comp`: overwrite automatically computed computational time step with this value.

Constructor
==========

All arguments are optional keyword arguments with the following default values:

    Options(;space_order=8, free_surface=false,
            limit_m=false, buffer_size=1e3,
            save_data_to_disk=false, file_path="",
            file_name="shot", sum_padding=false,
            optimal_checkpointing=false,
            num_checkpoints=nothing, checkpoints_maxmem=nothing,
            frequencies=[], isic=false,
            subsampling_factor=1, dft_subsampling_factor=1, return_array=false,
            dt_comp=nothing)

"""
Options(;space_order=8,
		 free_surface=false,
         limit_m=false,
		 buffer_size=1e3,
		 save_data_to_disk=false,
		 file_path="",
		 file_name="shot",
         sum_padding=false,
		 optimal_checkpointing=false,
		 num_checkpoints=nothing,
		 checkpoints_maxmem=nothing,
		 frequencies=[],
		 isic=false,
		 subsampling_factor=1,
		 dft_subsampling_factor=1,
         return_array=false,
         dt_comp=nothing) =
		 Options(space_order,
		 		 free_surface,
		         limit_m,
				 buffer_size,
				 save_data_to_disk,
				 file_path,
				 file_name,
				 sum_padding,
				 optimal_checkpointing,
				 num_checkpoints,
				 checkpoints_maxmem,
				 frequencies,
				 isic,
				 subsampling_factor,
				 dft_subsampling_factor,
                 return_array,
                 dt_comp)

function subsample(options::Options, srcnum)
    if isempty(options.frequencies)
        return options
    else
        opt_out = deepcopy(options)
        floc = options.frequencies[srcnum]
        typeof(floc) <: Array && (opt_out.frequencies = floc)
        return opt_out
    end
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
function reformulate(m, disj, bin_var, reformulation, param)
    vars = setdiff(all_variables(m), m[bin_var])
    @expression(m, original_model_variables, vars)
    for (i,constr) in enumerate(disj)
        if constr isa Vector || constr isa Tuple
            for (j,constr_j) in enumerate(constr)
                init_reformulation(m, constr_j, bin_var, reformulation, param, i, j)
            end
        elseif constr isa ConstraintRef || typeof(constr) <: Array || constr isa JuMP.Containers.DenseAxisArray
            init_reformulation(m, constr, bin_var, reformulation, param, i)
        end
    end
    if reformulation == :CHR
        add_disaggregated_constr(m, disj, vars)
    end
end

function init_reformulation(m, constr, bin_var, reformulation, param, i, j = missing)
    param = get_reform_param(param, i, j) #M or eps
    if constr isa ConstraintRef
        eval(:($reformulation($m, $constr, $bin_var, $i, $j, missing, $param)))
    elseif typeof(constr) <: Array
        for k in Iterators.product([1:s for s in size(constr)]...)
            eval(:($reformulation($m, $constr, $bin_var, $i, $j, $k, $param)))
        end
    elseif constr isa JuMP.Containers.DenseAxisArray
        for k in Iterators.product([s for s in constr.axes]...)
            eval(:($reformulation($m, $constr, $M, $bin_var, $i, $j, $k, $param)))
        end
    end
end

function get_reform_param(param, i, j)
    if param isa Number || ismissing(param)
        param = param
    elseif param isa Vector || param isa Tuple
        if param[i] isa Number
            param = param[i]
        elseif param[i] isa Vector || param[i] isa Tuple
            @assert !ismissing(j) "If constraint specific param values are provided, there must be more than one constraint in disjunct $i."
            @assert j <= length(param[i]) "If constraint specific param values are provided, a value must be provided for each constraint in disjunct $i."
            param = param[i][j]
        else
            error("Invalid param parameter provided for disjunct $i.")
        end
    else
        error("Invalid param parameter provided for disjunct $i.")
    end
end

function BMR(m, constr, bin_var, i, j, k, M)
    if ismissing(k)
        @assert is_valid(m,constr) "$constr is not a valid constraint in the model."
        ref = constr
    else
        @assert is_valid(m,constr[k...]) "$constr is not a valid constraint in the model."
        ref = constr[k...]
    end
    if ismissing(M)
        M = apply_interval_arithmetic(ref)
        @warn "No M value passed for $ref. M = $M was inferred from the variable bounds."
    end
    add_to_function_constant(ref, -M)
    bin_var_ref = variable_by_name(ref.model, string(bin_var[i]))
    set_normalized_coefficient(ref, bin_var_ref , M)
end

function apply_interval_arithmetic(ref)
    if ref isa NonlinearConstraintRef
        ref_str = string(ref)
        @assert length(findall(r"[<>]", ref_str)) <= 1 "$ref must be one of the following: GreaterThan, LessThan, or EqualTo."
        ref_func = replace(split(ref_str, r"[=<>]")[1], " " => "")
        ref_type = occursin(">", ref_str) ? :lower : :upper
        ref_rhs = 0 #Could be calculated with: parse(Float64,split(ref_str, " ")[end]). NOTE: @NLconstraint will always have a 0 RHS.
    elseif ref isa ConstraintRef
        ref_obj = constraint_object(ref)
        @assert ref_obj.set isa MOI.LessThan || ref_obj.set isa MOI.GreaterThan || ref_obj.set isa MOI.EqualTo "$ref must be one the following: GreaterThan, LessThan, or EqualTo."
        ref_func = string(ref_obj.func)
        ref_type = fieldnames(typeof(ref_obj.set))[1]
        ref_rhs = normalized_rhs(ref)
    end
    vars = all_variables(ref.model) #get all variable names
    for var in vars
        ub = has_upper_bound(var) ? upper_bound(var) : (is_binary(var) ? 1 : Inf)
        lb = has_lower_bound(var) ? lower_bound(var) : (is_binary(var) ? 0 : Inf)
        ref_func = replace(ref_func, "$var" => "($lb..$ub)")
    end
    func_bounds = eval(Meta.parse(ref_func))
    if ref_type == :lower
        M = func_bounds.lo - ref_rhs

    else
        M = func_bounds.hi - ref_rhs
    end
    if isinf(M)
        error("M parameter for $ref cannot be infered due to lack of variable bounds.")
    else
        return M
    end
end

function CHR(m, constr, bin_var, i, j, k, eps)
    if ismissing(k)
        if constr isa NonlinearConstraintRef #NOTE: CAN'T CHECK IF NL CONSTR IS VALID
            # @assert constr in keys(object_dictionary(m)) "$constr is not a named reference in the model."
        elseif constr isa ConstraintRef
            @assert is_valid(m,constr) "$constr is not a valid constraint in the model."
        end
        ref = constr
    else
        if constr isa NonlinearConstraintRef #NOTE: CAN'T CHECK IF NL CONSTR IS VALID
            # @assert constr in keys(object_dictionary(m)) "$constr is not a named reference in the model."
        elseif constr isa ConstraintRef
            @assert is_valid(m,constr[k...]) "$constr is not a valid constraint in the model."
        end
        ref = constr[k...]
    end
    bin_var_ref = variable_by_name(ref.model, string(bin_var[i]))
    for var in m[:original_model_variables]
        #get bounds for disaggregated variable
        @assert has_upper_bound(var) || is_binary(var) "Variable $var does not have an upper bound."
        @assert has_lower_bound(var) || is_binary(var) "Variable $var does not have a lower bound."
        UB = is_binary(var) ? 1 : upper_bound(var)
        LB = is_binary(var) ? 0 : lower_bound(var)
        #create disaggregated variable
        var_i = Symbol("$(var)_$i")
        if !(var_i in keys(object_dictionary(m)))
            eval(:(@variable($m, $LB <= $var_i <= $UB)))
            eval(:(@constraint($m, $LB * $bin_var_ref <= $var_i)))
            eval(:(@constraint($m, $var_i <= $UB * $bin_var_ref)))
            is_binary(var) && set_binary(m[var_i])
        end
    end
    #create convex hull constraint
    if ref isa NonlinearConstraintRef
        nl_perspective_function(ref, bin_var_ref, i, j, k, eps)
    elseif ref isa ConstraintRef
        lin_perspective_function(ref, bin_var_ref, i, j, k, eps)
    end
end

function lin_perspective_function(ref, bin_var_ref, i, j, k, eps)
    #check constraint type
    ref_obj = constraint_object(ref)
    @assert ref_obj.set isa MOI.LessThan || ref_obj.set isa MOI.GreaterThan || ref_obj.set isa MOI.EqualTo "$ref must be one the following: GreaterThan, LessThan, or EqualTo."
    for var in ref.model[:original_model_variables]
        #check var is present in the constraint
        coeff = normalized_coefficient(ref,var)
        iszero(coeff) && continue
        #modify constraint using convex hull
        rhs = normalized_rhs(ref) #get rhs
        var_i_ref = variable_by_name(ref.model, "$(var)_$i")
        set_normalized_rhs(ref,0) #set rhs to 0
        set_normalized_coefficient(ref, var, 0) #remove original variable
        set_normalized_coefficient(ref, var_i_ref, coeff) #add disaggregated variable
        set_normalized_coefficient(ref, bin_var_ref, -rhs) #add binary variable
    end
end

function nl_perspective_function(ref, bin_var_ref, i, j, k, eps)
    #extract info
    m = ref.model
    disj_name = replace("$(bin_var_ref)", "_binary" => "")
    vars = m[:original_model_variables]
    j = ismissing(j) ? "" : "_$j"
    k = ismissing(k) ? "" : "_$k"

    #check function has a single comparrison operator (<=, >=, ==)
    ref_str = string(ref)
    @assert length(findall(r"[<>]", ref_str)) <= 1 "$ref must be one of the following: GreaterThan, LessThan, or EqualTo."

    #create symbolic variables (using Symbolics.jl v0.1.25)
    sym_vars = []
    for var in vars
        var_sym = Symbol(var)
        push!(sym_vars, eval(:(Symbolics.@variables($var_sym)))[1])
    end
    ϵ = eps #epsilon parameter for perspective function (See Furman, Sawaya, Grossmann [2020] perspecive function)
    λ = Num(Symbolics.Sym{Float64}(Symbol(bin_var_ref)))
    FSG1 = Num(Symbolics.Sym{Float64}(gensym())) #this will become: [(1-ϵ)⋅λ + ϵ] (See Furman, Sawaya, Grossmann [2020] perspecive function)
    FSG2 = Num(Symbolics.Sym{Float64}(gensym())) #this will become: ϵ⋅(1-λ) (See Furman, Sawaya, Grossmann [2020] perspecive function)

    #convert ref_str into an Expr and extract comparrison operator (<=, >=, ==),
    #constraint function, and RHS
    ref_sym = Meta.parse(ref_str)
    ref_expr = ref_sym.args
    op = eval(ref_expr[1]) #comparrison operator
    rhs = ref_expr[3] #RHS of constraint
    gx = eval(ref_expr[2]) #convert the LHS of the constraint into a Symbolic function

    #use symbolic substitution to obtain the following expression:
    #[(1-ϵ)⋅λ + ϵ]⋅g(v/[(1-ϵ)⋅λ + ϵ]) - ϵ⋅g(0)⋅(1-λ) <= 0
    #first term
    g1 = FSG1*substitute(gx, Dict(var => var/FSG1 for var in sym_vars))
    #second term
    g2 = FSG2*substitute(gx, Dict(var => 0 for var in sym_vars))
    #create perspective function and simplify
    pers_func = simplify(g1 - g2, expand = true)
    #replace FSG expressions & simplify
    pers_func = substitute(pers_func, Dict(FSG1 => (1-ϵ)*λ+ϵ,
                                           FSG2 => ϵ*(1-λ)))
    pers_func = simplify(pers_func)

    #convert pers_func to Expr
    #use build_function from Symbolics.jl to convert pers_func into an Expr
    #where any operators (i.e. exp, *, <=) are replaced by the actual
    #operator (not the symbol).
    #This is done to later replace the symbolic variables with JuMP variables,
    #without messing with the math operators.
    pers_func_expr = Base.remove_linenums!(build_function(op(pers_func,rhs))).args[2].args[1]

    #replace symbolic variables by their JuMP variables
    replace_JuMPvars!(pers_func_expr, m)
    #replace the math operators by symbols
    replace_operators!(pers_func_expr)
    #add the constraint
    add_NL_constraint(m, pers_func_expr)

    #NOTE: the NLconstraint defined by `ref` needs to be deleted. However, this
    #   is not currently possible: https://github.com/jump-dev/JuMP.jl/issues/2355.
    #   As of today (5/12/21), JuMP is behind on its support for nonlinear systems.

    #NOTE: the new NLconstraint cannot be assigned a name (not an option in add_NL_constraint)
    # pers_func_name = Symbol("perspective_func_$(disj_name)$(j)$(k)")
end

function replace_JuMPvars!(expr, model)
    if expr isa Symbol
        return variable_by_name(model, string(expr))
    elseif expr isa Expr
        for i in eachindex(expr.args)
            expr.args[i] = replace_JuMPvars!(expr.args[i], model)
        end
    end
    expr
end

function replace_operators!(expr)
    if expr isa Expr
        for i in eachindex(expr.args)
            expr.args[i] = replace_operators!(expr.args[i])
        end
    elseif !isa(expr, Symbol) && !isa(expr, Number) && !isa(expr, VariableRef)
        return Symbol(expr)
    end
    expr
end

function add_disaggregated_constr(m, disj, vars)
    for var in vars
        d_vars = []
        for i in 1:length(disj)
            var_i = Symbol("$(var)_$i")
            var_i in keys(object_dictionary(m)) && push!(d_vars, m[var_i])
        end
        !isempty(d_vars) && eval(:(@constraint($m, $var == sum($d_vars))))
    end
end

################################################################################
###                             DEPRECATED                                   ###
################################################################################
function infer_BigM(ref)
    constr_set = constraint_object(ref).set
    set_fields = fieldnames(typeof(constr_set))
    @assert length(set_fields) == 1 "A reformulation cannot be done on constraint $ref because it is not one of the following GreaterThan, LessThan, or EqualTo."
    @assert :value in set_fields || :lower in set_fields || :upper in set_fields "$ref must be one the following: GreaterThan, LessThan, or EqualTo."
    vars = all_variables(ref.model) #get all variable names
    M = 0 #initialize M
    for var in vars
        coeff = normalized_coefficient(ref,var)
        iszero(coeff) && continue
        has_bounds = (has_lower_bound(var), has_upper_bound(var))
        if coeff > 0
            if :lower in set_fields && has_bounds[1]
                bound = lower_bound(var)
            elseif :upper in set_fields || :value in set_fields && has_bounds[2]
                bound = upper_bound(var)
            else
                error("M parameter cannot be infered due to lack of variable bounds for variable $var.")
            end
            M += coeff*bound
        elseif coeff < 0
            if :lower in set_fields && has_bounds[2]
                bound = upper_bound(var)
            elseif :upper in set_fields || :value in set_fields && has_bounds[1]
                bound = lower_bound(var)
            else
                error("M parameter cannot be infered due to lack of variable bounds for variable $var.")
            end
            M += coeff*bound
        end
    end
    if :lower in set_fields
        M -= constr_set.lower
    elseif :upper in set_fields
        M -= constr_set.upper
    elseif :value in set_fields
        M -= constr_set.value
    end
end


text:
import Base.convert

export BEDGraph

type BEDGraph <: Region
    contig :: String
    range :: Range1{Int}
    score :: Float64
end

register_file_extension!(BEDGraph, ".bg")
register_file_extension!(BEDGraph, ".bg.gz")
register_file_extension!(BEDGraph, ".bedGraph")
register_file_extension!(BEDGraph, ".bedGraph.gz")

function Base.read(io :: IO, ::Type{BEDGraph})
    while !eof(io)
        try
            return convert(BEDGraph, readline(io))
        catch
            continue
        end
    end
    throw(EOFError())
end

function Base.convert(::Type{BEDGraph}, s :: String)
    fields = split(strip(s), '\t')
    BEDGraph(fields[1], int(fields[2]):int(fields[3]), float(fields[4]))
end

strand (r :: BEDGraph) = STRAND_UNDEFINED
name (r :: BEDGraph) = ""


text:
using Turing, Test
using Turing: link, invlink, logpdf_with_trans
using ForwardDiff: derivative, jacobian
using LinearAlgebra: logabsdet

# logabsdet doesn't handle scalars.
_logabsdet(x::AbstractArray) = logabsdet(x)[1]
_logabsdet(x::Real) = log(abs(x))

# Generate a (vector / matrix of) random number(s).
_rand_real(::Real) = randn()
_rand_real(x) = randn(size(x))

# Standard tests for all distributions involving a single-sample.
function single_sample_tests(dist, jacobian)

    # Do the regular single-sample tests.
    single_sample_tests(dist)

    # Check that the implementation of the logpdf agrees with the AD version.
    x = rand(dist)
    logpdf_ad = logpdf(dist, x) - _logabsdet(jacobian(x->link(dist, x), x))
    @test logpdf_ad ≈ logpdf_with_trans(dist, x, true)
end

# Standard tests for all distributions involving a single-sample. Doesn't check that the
# logpdf implementation is consistent with the link function for technical reasons.
function single_sample_tests(dist)

    # Check that invlink is inverse of link.
    x = rand(dist)
    @test invlink(dist, link(dist, copy(x))) ≈ x atol=1e-9

    # Check that link is inverse of invlink. Hopefully this just holds given the above...
    y = link(dist, x)
    @test link(dist, invlink(dist, copy(y))) ≈ y atol=1e-9

    # This should probably be exact.
    @test logpdf(dist, x) == logpdf_with_trans(dist, x, false)

    # Check that invlink maps back to the apppropriate constrained domain.
    @test all(isfinite, logpdf.(Ref(dist), [invlink(dist, _rand_real(x)) for _ in 1:100]))

    # This is a quirk of the current implementation, of which it would be nice to be rid.
    @test typeof(x) == typeof(y)
end

# Standard tests for all distributions involving multiple samples. xs should be whatever
# the appropriate repeated version of x is for the distribution in question. ie. for
# univariate distributions, just a vector of identical values. For vector-valued
# distributions, a matrix whose columns are identical.
function multi_sample_tests(dist, x, xs, N)
    ys = link(dist, copy(xs))
    @test invlink(dist, link(dist, copy(xs))) ≈ xs atol=1e-9
    @test link(dist, invlink(dist, copy(ys))) ≈ ys atol=1e-9
    @test logpdf_with_trans(dist, xs, true) == fill(logpdf_with_trans(dist, x, true), N)
    @test logpdf_with_trans(dist, xs, false) == fill(logpdf_with_trans(dist, x, false), N)

    # This is a quirk of the current implementation, of which it would be nice to be rid.
    @test typeof(xs) == typeof(ys)
end

# Scalar tests
@testset "scalar" begin
let
    # Tests with scalar-valued distributions.
    uni_dists = [
        Arcsine(2, 4),
        Beta(2,2),
        BetaPrime(),
        Biweight(),
        Cauchy(),
        Chi(3),
        Chisq(2),
        Cosine(),
        Epanechnikov(),
        Erlang(),
        Exponential(),
        FDist(1, 1),
        Frechet(),
        Gamma(),
        InverseGamma(),
        InverseGaussian(),
        # Kolmogorov(),
        Laplace(),
        Levy(),
        Logistic(),
        LogNormal(1.0, 2.5),
        Normal(0.1, 2.5),
        Pareto(),
        Rayleigh(1.0),
        TDist(2),
        TruncatedNormal(0, 1, -Inf, 2),
    ]
    for dist in uni_dists

        single_sample_tests(dist, derivative)

        # specialised multi-sample tests.
        N = 10
        x = rand(dist)
        xs = fill(x, N)
        multi_sample_tests(dist, x, xs, N)
    end
end
end

# Tests with vector-valued distributions.
@testset "vector" begin
let
    vector_dists = [
        Dirichlet(2, 3),
        MvNormal(randn(10), exp.(randn(10))),
        MvLogNormal(MvNormal(randn(10), exp.(randn(10)))),
    ]
    for dist in vector_dists

        if dist isa Dirichlet
            single_sample_tests(dist)

            # This should fail at the minute. Not sure what the correct way to test this is.
            x = rand(dist)
            logpdf_turing = logpdf_with_trans(dist, x, true)
            J = jacobian(x->link(dist, x), x)
            @test_broken logpdf(dist, x) - _logabsdet(J) ≈ logpdf_turing
        else
            single_sample_tests(dist, jacobian)
        end

        # Multi-sample tests. Columns are observations due to Distributions.jl conventions.
        N = 10
        x = rand(dist)
        xs = repeat(x, 1, N)
        multi_sample_tests(dist, x, xs, N)
    end
end
end

# Tests with matrix-valued distributions.
@testset "matrix" begin
let
    matrix_dists = [
        Wishart(7, [1 0.5; 0.5 1]),
        InverseWishart(2, [1 0.5; 0.5 1]),
    ]
    for dist in matrix_dists

        single_sample_tests(dist)

        # This should fail at the minute. Not sure what the correct way to test this is.
        x = rand(dist)
        logpdf_turing = logpdf_with_trans(dist, x, true)
        J = jacobian(x->link(dist, x), x)
        @test_broken logpdf(dist, x) - _logabsdet(J) ≈ logpdf_turing

        # Multi-sample tests comprising vectors of matrices.
        N = 10
        x = rand(dist)
        xs = [x for _ in 1:N]
        multi_sample_tests(dist, x, xs, N)
    end
end
end

################################## Miscelaneous old tests ##################################

# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1000., -1000., -1000.]), true)
# NaN
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1000., -1000., -1000.], true, true)
# -1999.30685281944
#
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1., -2., -3.]), true)
# -3.006450206744678
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1., -2., -3.], true, true)
# -3.006450206744678
d  = Dirichlet([1., 1., 1.])
r  = [-1000., -1000., -1000.]
r2 = [-1., -2., -3.]

# test link
#link(d, r)

# test invlink
@test invlink(d, r) ≈ [0., 0., 1.] atol=1e-9

# test logpdf_with_trans
#@test logpdf_with_trans(d, invlink(d, r), true) -1999.30685281944 1e-9 ≈ # atol=NaN
@test logpdf_with_trans(d, invlink(d, r2), true) ≈ -3.760398892580863 atol=1e-9


text:
#!/bin/env julia

using Gadfly, DataFrames

Gadfly.set_default_plot_size(28cm, 16cm)

fname = ARGS[1]

df = readtable(fname * ".csv")

p = plot(df, x=:Days_since_ref, y=:Obj, ymin=:Obj_min, ymax=:Obj_max, Geom.line, Geom.ribbon);

draw(PDF("$(fname).pdf", 5inch, 3inch), p)
draw(PNG("$(fname).png", 5inch, 3inch), p)

#fo = open("$(fname).df", "w")
#println(fo, "$df")

writetable("$(fname).dat", df, separator = '|')

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
using Test, ParametricAirfoils, Ipopt, JuMP
include("testnaca.jl")


text:
Lar = LinearAlgebraicRepresentation

using IntervalTrees
using SparseArrays
using NearestNeighbors
using DataStructures

#-------------------------------------------------------------------------------
#   PAIRWISE DECOMPOSITION
#-------------------------------------------------------------------------------

"""
	facesplitting(model::CAGD.Model)::CAGD.Model

Performs the pairwise decomposition of all faces with the possibly intersecting.
"""

function facesplitting(model; atol = 10e-8)
    sp_idx = CAGD.spaceIndex(model, 2)
    de_models = [
        CAGD.face_decomposition(model, face_idx, sp_idx[face_idx], atol = atol)
        for face_idx = 1 : size(model, 2, 1)
    ]
    return CAGD.uniteMultipleModels(de_models)
end

function face_decomposition(
		model::CAGD.Model,
		face_idx::Int,
		sp_idx_face::Array{Int,1};
		atol = 1e-8
	)::CAGD.Model

	function face_intersection(													# Inserire il modello in depth
			V::Lar.Points,
			EV::Lar.ChainOp,
			face::Lar.Cell
		)::CAGD.Model

		# Get indices of `face` vertices in traversal order
	    Vidx = Lar.buildFV(EV, face)
	    retV = Lar.Points(undef, 2, 0)

	    visited_verts = []
		for i in 1:length(Vidx)
			# vector to vertex
			o = V[:, Vidx[i]]
			# next vertex (circular ordering)
			j = i < length(Vidx) ? i+1 : 1
			# vector from i-th to j-th
	        d = V[:, Vidx[j]] - o

			err = atol

			# if the point is coplanar, it should be added
			if (-err < d[3] < err)
				#retV = [retV V[1:2, Vidx[i]]]                                  Da aggiungere solo se la faccia è tutta parallela
				#retV = [retV V[1:2, Vidx[j]]]
			else

				# check whether the z-ratio between o and d is ∈ [-1, 0]
				# then there is an intersextion point `p = o + α * d`
	            alpha = -o[3] / d[3]

	            if -err <= alpha <= 1+err
	                p = o + alpha*d

	                if -err < alpha < err || 1-err < alpha < 1+err
	                    if !(Lar.vin(p, visited_verts))
	                        push!(visited_verts, p)
	                        retV = [retV reshape(p, 3, 1)[1:2, 1]]
	                    end
	                else
	                    retV = [retV reshape(p, 3, 1)[1:2, 1]]
	                end
	            end
	        end

	    end

	    vnum = size(retV, 2)

	    if vnum == 1
	        vnum = 0
	        retV = Lar.Points(undef, 2, 0)
	    end
	    enum = (÷)(vnum, 2)
	    retEV = spzeros(Int8, enum, vnum)

	    for i in 1:enum
	        retEV[i, 2*i-1:2*i] = [-1, 1]
	    end

	    retModel = CAGD.Model(retV)
		CAGD.addModelCells!(retModel, 1, retEV)
		return retModel
	end

	#== THE FUNCTION ==#

	Gface, Gfaceidx = CAGD.getModelCellGeometry(model, 2, face_idx, true)
	M = CAGD.build_projection_matrix(Gface)
	PG = (M * [model.G; ones(1, size(model, 0, 2))])[1:3, :]

	# Generate face Model (in 2D)
	Pmodel = CAGD.Model(PG[1:2, Gfaceidx])
	CAGD.addModelCells!(Pmodel, 1, model.T[1][CAGD.getModelLoCell(model, 2, face_idx), Gfaceidx])

	# Each other face adds new pieces
	for f in sp_idx_face
		# println("Intersecting $face_idx with $f")
		f_FE = model.T[2][f, :]
		f_EV = model.T[1][CAGD.getModelLoCell(model, 2, f), :]
		_, vs, _ = findnz(f_EV);
		if *((abs.(PG[3, vs]) .< atol)...)
			newVidx = sort(unique(vs));
			oldVidx = SparseArrays.spzeros(Int, size(PG, 2), 1)
			for i = 1 : length(newVidx)  oldVidx[newVidx[i]] = i  end
			newModel = CAGD.Model(PG[1:2, newVidx])
			I, J, Val = findnz(f_EV)
			f_EV = SparseArrays.sparse(I, map(x -> oldVidx[x], J), Val)#, max(I...), f_EV.n)
			CAGD.addModelCells!(newModel, 1, f_EV)
			#newModel = CAGD.getSubModel(model, 2, f)
		else
			newModel = face_intersection(PG, model.T[1], f_FE)
		end
		CAGD.uniteModels!(Pmodel, newModel)
	end

	Pmodel = CAGD.mergeModelVertices(Pmodel, signed_merge = false)

	model = CAGD.planar_arrangement(Pmodel, sparsevec(ones(Int8, length(Gfaceidx))))
	@assert !isnothing(model) "UNEXPECTED ERROR: a face should be mapped to itself"

	n = size(model, 0, 2)
	retModel = CAGD.Model((inv(M)*[model.G; zeros(1, n); ones(1, n)])[1:3, :])
	CAGD.addModelCells!(retModel, 1, model.T[1])
	CAGD.addModelCells!(retModel, 2, model.T[2])

	return retModel
end


text:
function plot(data::modwt)
	if data.series!=1
		error("Specify which series decomposition you want to plot.")
	else
		plot(data, 1)
	end
end

function plot(data::modwt, ser::Int)
	using Gadfly
	using DataFrames

	import Gadfly.plot

	df = convert(DataFrame, [reshape(data.W[:,:,ser], data.L, data.level) data.V[:,data.level,ser] data.original[:,ser] [1:data.L]])
	names!(df, [[symbol("W$i") for i in 1:data.level], symbol("V$levels"), symbol("Originalseries"), symbol("Time")])
	plotdata = stack(df, [1:(data.level+2)], [data.level+3])
	plot(plotdata, ygroup="variable", y="value", x="Time", Geom.subplot_grid(Geom.line))
end


text:
function lag0(x::AbstractVector, p::Int = 1)
    T = length(x);
    return [zeros(eltype(x), p); x[1:T - p]]
end

function lag0(x::AbstractArray, p::Int = 1)
    T, n = size(x);
    return [zeros(eltype(x), p, n); x[1:T - p,:]]
end

function eye(m::Int, n::Int = m)
    return Matrix{Float64}(I, m, n)
end

maxabs(x) = maximum(abs, x)

function vech(A::AbstractMatrix, k::Int = 0)
    # return A[tril!(trues(size(A)), 0)]
    return A[tril(trues(size(A)), k)]
end

function vech(A::AbstractArray{T,3}, k::Int = 0) where {T}
    B = vech(A[:,:,1], k)
    a = size(A, 3)

    if a > 1
        for i = 2:a
            append!(B, vech(A[:,:,i], k))
        end
    end

    return B
end

function vechRshp(A::AbstractArray{T,3}, k::Int = 0) where {T}
    a = size(A, 3)
    return reshape(vech(A, k), (:, a))
end

# this one is type unstable
function vec2ltri2(v::AbstractVector)
    n = length(v)
    s = round(Int, (sqrt(8n + 1) - 1) / 2)
    @assert s * (s + 1) / 2 == n "length of vector is not triangular"
    k = 0; [ i >= j ? (k += 1; v[k]) : zero(eltype(v)) for i = 1:s, j = 1:s ]
end

function vec2ltri(v::AbstractVector)
    n = length(v)
    s = round(Int, (sqrt(8n + 1) - 1) / 2)
    @assert s * (s + 1) / 2 == n "length of vector is not triangular"
    A = zeros(eltype(v), s, s)
    k = 0
    for j in 1:s, i in j:s
        k += 1
        A[i,j] = v[k]
    end
    return A
end

function vec2ltriW1(v::AbstractVector{T}, dm::Int = 1)::Array{T,2} where {T}
    n = length(v)
    s = round(Int, (sqrt(8n + 1) + 1) / 2)
    @assert s * (s - 1) / 2 == n "length of vector is not triangular below diagonal "
    A = Matrix{T}(I, s, s)
    k = 0
    if dm == 1
        for j in 1:s, i in j:s
            if i > j
                k += 1
                A[i,j] = v[k]
            end
        end
    elseif dm == 2
        for i in 1:s, j in 1:i
            if i > j
                k += 1
                A[i,j] = v[k]
            end
        end
    end
    return A
end

function vec2sym(v::AbstractVector)
    return Symmetric(vec2ltri(v), :L)
end

function vec2sym(v::AbstractMatrix, dm::Int = 1)
    if dm == 2
        v = v'
    end

    a, b = size(v)
    s = round(Int, (sqrt(8a + 1) - 1) / 2)

    A = similar(v, s, s, b)
    for i = 1:b
        A[:,:,i] = vec2sym(v[:,i])
    end
    return A
end

function companionf(beta::AbstractMatrix, n::Int, p::Int, c::Bool = false)
    return [beta[:, 1:end - c];eye(n * p - n) zeros(n * p - n, n)]
end

function companionf(beta::AbstractVector, n::Int, p::Int, c::Bool = false)
    @assert n == 1 throw(ArgumentError("provided beta is univariate"))
    if c || p > 1
        return [reshape(beta[1:end - c], 1, :); eye(p - 1) zeros(p - 1)]
    else
        return Matrix{Float64}(reshape(beta, 1, 1))
    end
end

function stabcheck(beta::AbstractVecOrMat, n::Int, p::Int, c::Bool = false)::Bool
    # beta is (Nx(N*L+c))
    FF = companionf(beta, n, p, c)
    return maximum(abs, eigvals(FF)) >= 1.
end

function stabcheckC(beta::AbstractVecOrMat, n::Int, p::Int, c::Bool = false)::Bool
    # beta is (Nx(N*L+c))
    # stability check which ignores complex numbers
    FF = companionf(beta, n, p, c)
    D = eigvals(FF)
    if isreal(D)
        S = maximum(abs, D) >= 1.
    else
        S = true
    end
    return S
end

function stabcheck(beta::Number, n::Int = 1, p::Int = 1, c::Bool = false)::Bool
    return abs(beta) >= one(eltype(beta))
end
stabcheckC(beta::Number, n::Int = 1, p::Int = 1, c::Bool = false) = stabcheck(beta, n, p, c)

function preparexy(data::AbstractVector, p::Int, c::Bool = false)
    T = length(data);
    Y = data[p + 1:end];
    X = similar(data, T, p + c);
    for i = 1:p
        X[:,(i - 1) + 1:i] = lag0(data, i);
    end
    if c
        X[:,end] = ones(eltype(data), T);
    end
    X = X[p + 1:end,:]
    return Y, X
end

function preparexy(data::AbstractMatrix, p::Int, c::Bool = false)
    T, n = size(data);
    Y = data[p + 1:end,:];
    X = similar(data, T, n * p + c);
    for i = 1:p
        X[:,n * (i - 1) + 1:n * i] = lag0(data, i);
    end
    if c
        X[:,end] = ones(eltype(data), T);
    end
    X = X[p + 1:end,:]
    return Y, X
end

function sumsqr(a::AbstractArray)
    return mapreduce(x -> x^2, +, a)
end

# function wish(h, n::Int)
#     A = cholesky!(h).L * randn(size(h, 1), n)
#     A = A * A'
#     return Symmetric(A)
# end

function wish(h::AbstractArray{S,2}, n::Int)   where {S <: AbstractFloat}
    A = randn(n, size(h, 1)) * cholesky(h).U
    return A'A
end

function wish(h::AbstractArray{S,1}, n::Int)   where {S <: AbstractFloat}
    A = randn(n) .* sqrt.(h)
    return dot(A', A)
end

function iwish(h::AbstractMatrix, n::Int)
    B = randn(n, size(h, 1))
    A = cholesky(h).L / qr!(B).R
    return A * A'
end

function iwish(h::AbstractVector, n::Int)
    B = randn(n, 1)
    A = sqrt.(h) / qr!(B).R
    return dot(A, A')
end

function genPSDmat(K::Int = rand(2:9))
    A = randn(K, K)
    return A'A
end

function genPSDmatStrict(K::Int = rand(2:9); no::Int = rand(1:K - 1))
    A = randn(K, K)
    A = A'A
    vals, vecs = eigen(A)
    vals[1:no] .= 0.
    return vecs * Diagonal(vals) * vecs'
end

function genPDmat(K::Int = rand(2:9))
    A = zeros(K, K)
    while isposdef(A) == false
        A .= randn(K, K)
        A .= A'A
    end
    return A
end

function cholPSD(A::AbstractMatrix{T})::Array{T,2}  where {T <: AbstractFloat}
    B = copy(Symmetric(A));
    if isposdef(B)

        # B = cholesky(B)
        # return convert(Array{Float64,2}, B.U)

        return cholesky!(B).U

    else
        # A = sqrt(A)
        # # D,U = eig(A);

        # A = eigen(A);

        # tol = length(A.values) * eps(maximum(A.values));
        # # D1 = A.values[A.values .> tol];
        # A = diagm(sqrt.(A.vectors)) * A.vectors[:, A.values .> tol]';

        # M = eigen(B)
        # D = M.values
        # U = M.vectors

        D, U = eigen(B)
        tol = length(D) * eps(maximum(D));
        # D1 = D[D .> tol];
        # # A = diagm(0 => sqrt.(D1)) * U[:, D .> tol]';
        # C = Diagonal(sqrt.(D1)) * U[:, D .> tol]';

        # return Diagonal(sqrt.(filter(z -> z > tol, D))) * U[:, D .> tol]'
        cord = findall(z -> z > tol, D)
        # return Diagonal(sqrt.(filter(z -> z > tol, D))) * U[:, findall(z -> z > tol, D)]'
        return Diagonal(sqrt.(D[cord])) * U[:, cord]'

    end
end

"the function regularises the matrix to be positive definite given that there are no eigenvalues of zeros"
function regMat2PD!(A::Matrix)
    D, U = eigen(A)
    A .= Symmetric(U * Diagonal(abs.(D)) * U')
    return nothing
end

function randnPSD(mu::AbstractMatrix{T}, sigma::AbstractMatrix{T})::Array{T}  where {T <: AbstractFloat}
    # input:
    # mu is row vector
    #
    # output:
    # draw is a row vector
    sigmaN = cholPSD(sigma)
    draw = mu + randn(1, size(sigmaN, 1)) * sigmaN
    return draw
end

function randnPSD(mu::AbstractVector{T}, sigma::AbstractMatrix{T})::Array{T}  where {T <: AbstractFloat}
    # input:
    # mu is row vector
    #
    # output:
    # draw is a row vector
    sigmaN = cholPSD(sigma)
    draw = mu' + randn(1, size(sigmaN, 1)) * sigmaN
    return draw
end

function mcSS(PI::Matrix)
    a = eigen(Matrix(PI'))
    b = a.vectors[:,a.values .≈ 1.]
    return vec(b) ./ sum(b)
end

function mcSS2(PI::Matrix, tol::Float64 = eps())
    l = size(PI, 1)
    q = ones(l) ./ l
    q1 = similar(q)
    a = 2.
    while a > tol
        q1 .= PI' * q
        a = maximum(abs, q1 .- q)
        q .= q1
    end
    return q
end

function quantileArr(v::AbstractArray{T,2}, p, dms::Int = 1) where {T}
    a, b = size(v)
    w = length(p)
    p = vec(p)
    if dms == 1
        A = similar(v, w, b)
        for i = 1:b
            A[:,i] = quantile(v[:,i], p)
        end
    elseif dms == 2
        A = similar(v, a, w)
        for i = 1:a
            A[i,:] = quantile(v[i,:], p)
        end
    end
    return A
end

function quantileArr(v::AbstractArray{T,3}, p, dms::Int = 1) where {T}
    a, b, c = size(v)
    w = length(p)
    p = vec(p)

    if dms == 1
        A = similar(v, w, b, c)
        for i in 1:b, j in 1:c
            A[:,i,j] = quantile(v[:,i,j], p)
        end
    elseif dms == 2
        A = similar(v, a, w, c)
        for i in 1:a, j in 1:c
            A[i,:,j] = quantile(v[i,:,j], p)
        end
    elseif dms == 3
        A = similar(v, a, b, w)
        for i in 1:a, j in 1:b
            A[i,j,:] = quantile(v[i,j,:], p)
        end
    end

    return A
end

ismyapprox(a, b, dgts::Int) = isapprox(a, b, atol = dgts)

function normpdf(x::Float64, mu::Float64 = 0.0, sigma::Float64 = 1.0)
    return exp(-0.5 * ((x - mu) ./ sigma).^2) ./ (sqrt(2. * pi) .* sigma)
end

function acf(y::AbstractVector, q::Int)
# function assumes that there are no missings

    z = y .- mean(y)
    T = length(y)

    ac0 = dot(z, z) ./ T
    ac = dot(z[1:end - q], z[q + 1:end]) ./ T
    ac /= ac0; # Normalize to correllation

    return ac
end

function acf(y::AbstractVector, q::AbstractRange = 1:10)
    # function assumes that there are no missings

    z = y .- mean(y)
    T = length(y)

    ac0 = dot(z, z) ./ T
    ac = Array{Float64}(undef, 0);
    for j = q
        push!(ac, dot(z[1:end - j], z[j + 1:end]) ./ T)
    end
    ac ./= ac0; # Normalize to correllation

    return ac
end

function acf(y::AbstractMatrix, q = 1:10)

    T, n = size(y)

    ac = Array{Float64}(undef, 0);

    for i = 1:n
        append!(ac, acf(y[:,i], q))
    end

    ac = reshape(ac, length(q), n)

    return ac
end

function decVCV(VCV::AbstractVector, getAinv::Bool = false)
    Omega = vec2sym(VCV)
    Omega = cholesky(Omega).L
    sigma = diag(Omega)
    A = Omega / Diagonal(sigma)
    if getAinv
        A = inv(A)
    end
    A = vech(A, -1)

    return A, sigma
end

function decVCV(VCV::Matrix{Float64}, getAinv::Bool = false)

    a, T = size(VCV)

    s = round(Int, (sqrt(8a + 1) - 1) / 2)
    s * (s + 1) / 2 == a || error("fromVCV: length of vector is not triangular")

    sigma = Array{Float64}(undef, s, T)
    A = Array{Float64}(undef, round(Int, s * (s - 1) / 2), T)

    for i = 1:T
        A[:,i], sigma[:,i] = decVCV(VCV[:,i], getAinv)
    end

    return A, sigma
end

function ols1(y::AbstractVector, x::AbstractMatrix)
    nobs, nvar = size(x)
    bhat = x \ y
    yhat = x * bhat
    res = y - yhat
    sig2hat = (res'res) / (nobs - nvar)
    sigbhat = sig2hat .* eye(nvar) / (x' * x)
    R2 = var(yhat) / var(y)

    return (nobs = nobs, nvar = nvar, bhat = bhat, yhat = yhat, res = res, sig2hat = sig2hat, sigbhat = sigbhat, R2 = R2)
end

function VARols1(data::AbstractVecOrMat{Float64}, cnst::Bool, P::Int)
    y, x = preparexy(data, P, cnst)
    nobs, nvar = size(x)
    bhat = x \ y
    yhat = x * bhat
    res = y - yhat
    sig2hat = (res'res) / (nobs - nvar)

    return (nobs = nobs, nvar = nvar, bhat = bhat, yhat = yhat, res = res, sig2hat = sig2hat)
end

function nanmean(x::AbstractArray{T}) where T <: AbstractFloat
    sum = zero(eltype(x))
    count = 0
    for i in x
        if !isnan(i)
            sum += i
            count += 1
        end
    end
    return sum / count
end

nanmean(x::AbstractArray, y::Int) = mapslices(nanmean, x; dims = y)

function inbetween(a::Matrix, b::AbstractArray{T,3})  where {T}
    return (view(b, :, :, 1) .<= a) .& (a .<= view(b, :, :, 2))
end

"The function multiplies matrices A and B and extracts the diagonal to a prealoccated vector v"
function getmultdiag!(v::Vector, A::Matrix, B::Matrix)
    # v .= dot.(eachrow(A), eachcol(B))
    tmp = A .* B'
    sum!(v, tmp)
    return nothing
end

"the function that transforms the value from domain (min,max) to range with minus/plus infinity"
function transf1To(a::Float64, minn::Float64, maxx::Float64)
   return -log((maxx - a) ./ (a - minn))
end

"the function that transforms the value back from minus/plus infinity  to range with domain (min,max)"
function transf1Back(a::Float64, minn::Float64, maxx::Float64)
    return minn + (maxx - minn) / (1. + exp(-a))
end

"Function to split samples."
function split_data(df, at = 0.70)
    r = size(df, 1)
    index = Int(round(r * at))
    train = df[1:index, :]
    test  = df[(index + 1):end, :]
    return train, test
end

    "A handy helper function to standardize the data"
function standardize(x;dims::Int = 1)
    return (x .- mean(x, dims = dims)) ./ std(x, dims = dims), x
end

"A handy helper function to unstandardize the data"
function unstandardize(x, orig;dims::Int = 1)
    return (x .+ mean(orig, dims = dims)) .* std(orig, dims = dims)
end

"a function that prints the loop progress"
function showprog(iter::Int, Nreps::Int; it_print::Int = floor(Nreps / iter))
    if mod(iter, it_print) == 0
        println("Iter: $iter / $Nreps ")
    end
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
using HOODESolver
using LinearAlgebra
using Plots
using Plots.PlotMeasures
function getmindif(tab::Vector{Vector{BigFloat}})
    nmmin = Inf
    ret = (0, 0)
    for i=1:size(tab,1)
        for j=(i+1):size(tab,1)
            nm = norm(tab[i]-tab[j], Inf)
            if nm < nmmin
                nmmin = nm
                ret = i, j
            end
        end
    end
    return ret
end
function fctmain(n_tau, prec)
    setprecision(prec)
    u0=BigFloat.([90, -44, 83, 13]//100)
    t_max = big"1.0"
    epsilon=big"0.0017"
    println("epsilon=$epsilon")
    nbmaxtest=12
    ordmax=17
    debord=3
    pasord=2
    y = ones(Float64, nbmaxtest, div(ordmax-debord,pasord)+1 )
    res_gen = Array{ Array{BigFloat,1}, 2}(undef, nbmaxtest, div(ordmax-debord,pasord)+1 )
    x=zeros(Float64,nbmaxtest)
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    prob = HOODEProblem(henon_heiles, u0, (big"0.0",t_max), missing, A, epsilon)
    tabsol = Array{Array{BigFloat,1},1}(undef,0)
    println("ordmax=$ordmax")
    indref = 1
    ind=1
    for order=debord:pasord:ordmax
        println("eps=$epsilon order=$order")
        nb = 100
        indc = 1
        labels=Array{String,2}(undef, 1, div(order-debord, pasord)+1)
        resnorm=0
        resnormprec=1
        sol =undef
        println("preparation ordre $order + 2")
        par_u0 = missing
        borindc = (order < ordmax) ? nbmaxtest : nbmaxtest+1
        while indc <= borindc
            if indc > nbmaxtest
                nb = 51*2^nbmaxtest
            end
            res = solve(prob, nb_tau=n_tau, order=order, nb_t=nb,par_u0=par_u0, dense=false)
            par_u0=res.par_u0
            sol = res[end]
            push!(tabsol, sol)
            if indc <= nbmaxtest
                res_gen[indc,ind] = sol
                x[indc] = 1.0/nb
            end
            (a, b) = getmindif(tabsol)
            if a != 0
                indref = a
                solref = tabsol[b]
                for i=1:ind
                    borne = (i <ind) ? size(res_gen,1) : min(indc, size(res_gen,1))
                    for j = 1:borne
                        nm = norm(res_gen[j,i] - solref, Inf)
                        if nm == 0
                            nm = norm(res_gen[j,i] - tabsol[a], Inf)
                        end
                        y[j,i] = (nm < 1) ? nm : NaN
                   end
                end
            end
            println("result=$y")
            println("log2(y)=$(log2.(y))")
            nb *= 2
            indc += 1
        end
        for i=debord:pasord:order
            labels[1,(i-debord)÷pasord+1] = " order = $i "
        end
        p=Plots.plot(
                        x,
                        view(y,:,1:ind),
                        xlabel="Δt",
                        xaxis=:log,
                        ylabel="error",
                        yaxis=:log,
                        legend=:bottomright,
                        label=labels,
                        marker=2
                    )
        pp=Plots.plot(
                        x,
                        view(y,:,1:ind),
                        xlabel="Δt",
                        xaxis=:log,
                        ylabel="error",
                        yaxis=:log,
                        legend=:bottomright,
                        label=labels,
                        marker=2,
                        bottom_margin=30px,
                        left_margin=60px
            )
        prec_v = precision(BigFloat)
        eps_v = convert(Float32,epsilon)
        Plots.savefig(p,"out/res4p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.pdf")
        Plots.savefig(p,"out/res4p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.png")
        Plots.savefig(pp,"out/res4pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.pdf")
        Plots.savefig(pp,"out/res4pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_henon_heiles.png")
        ind+= 1
    end
end

# testODESolver()

# for i=3:9
#     fctMain(2^i)
# end
fctmain(32,512)