text:
#=##############################################################################
# DESCRIPTION
    Methods for generation of surface geometries.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jul 2018
  * License   : MIT License
=###############################################################################

# LOFTING ######################################################################
"""
  `generate_loft(crosssections, bscale, b_pos, chords, twists, LE_x, LE_z)`

  Generates a lofted surface geometry. The name of the arguments are taken from
  its initial application to the lofting of an aircraft wing, so bear with it.

  **Arguments**
  * `bscale::Float64`         : Semi-span scale. It will scale the entire
                                geometry by this factor.
  * `b_pos::Array{Float64, 1}`: Normalized span positions y/bscale of the
                                following distributions.
  * `chords::Array{Float64,1}`: Chord c/bscale distribution.
  * `twists::Array{Float64,1}`: Twist (deg) distribution.
  * `LE_x::Array{Float64,1}`  : x-position (chordwise) x/bscale of leading edge
                                distribution.
  * `LE_z::Array{Float64,1}`  : z-position (dihedral-wise) z/bscale of leading
                                edge distribution.
  * `crosssections::Array{Tuple{T,Array{T,2}}, 1}`    : cross sections along the
                                span in the form [(y/bscale, crosssection)],
                                where `crosssection` is a matrix that contains
                                all points of the airfoil contour indexed by
                                row. In order to define the resulting normals
                                pointing out of the geometry, the points in
                                `crosssection` must start from the trailing
                                edge, go around the top side towards the leading
                                edge, and back to the trailing edge around the
                                bottom side. ALL CROSS SECTIONS MUST HAVE THE
                                SAME NUMBER OF POINTS.

  **Optional Arguments**
  * `tilt_z::Array{Float64,1}`: Tilting (deg) about the z-axis of every span
                                cross section. This is also a distribution.
  * `symmetric::Bool`         : Whether to consider the `crosssections` to be
                                symmetric about the semi-span. If true, only
                                positive y/bscale are neeeded.

  NOTE: The resulting geometry will be a structured surface mesh with the first
  index going in chordwise direction starting from the TE around the bottom
  surface and around back to the TE. The second index is the in the spanwise
  direction starting at left end (lowest position along y axis) to right end.
"""
function generate_loft(crosssections::Array{Tuple{T,Array{T,2}}, 1},
                        bscale::Real, b_pos::Array{T,1}, chords::Array{T,1},
                        twists::Array{T,1},
                        LE_x::Array{T,1}, LE_z::Array{T,1};
                        # MORE GEOMETRIC OPTIONS
                        tilt_z=nothing, tilt_y=nothing,
                        symmetric::Bool=false,
                        loop_dim::Int64=1,
                        # OUTPUT OPTIONS
                        save_path=nothing, paraview::Bool=true,
                        file_name::String="myloft"
                       ) where{T<:Real}

  sec_NDIVS = size(crosssections[1][2], 1)-1    # Cross-section cells
  b_NDIVS = size(b_pos, 1)-1                    # Span cells

  # ERROR CASES
  for (pos, sec) in crosssections
    if sec_NDIVS!=size(sec, 1)-1
      error("All airfoil sections must have the same number of points.")
    end
  end

  for (nam,val) in [("chords",chords), ("twists",twists), ("LE_x",LE_x),
                    ("LE_z",LE_z), ("tilt_z",tilt_z), ("tilt_y",tilt_y)]
    if (nam=="tilt_z" && tilt_z==nothing) || (nam=="tilt_y" && tilt_y==nothing)
      nothing
    else
      if size(b_pos,1)!=size(val,1)
          error("The size of `b_pos` doesn't match the size of `$nam`."*
            " Expected size $(size(b_pos,1)), got $(size(val,1)).")
      end
    end
  end

  # ----------------- PARAMETRIC GRID ------------------------------------------
  P_min = [0, 0, 0]               # Lower boundary arclength, span, dummy
  P_max = [1, 1, 0]               # Upper boundary arclength, span, dummy
  # loop_dim = 1                  # Loops the arclength dimension

  # Adds dummy division
  NDIVS = [sec_NDIVS, b_NDIVS, 0]

  # Generates parametric grid
  grid = Grid(P_min, P_max, NDIVS, loop_dim)

  # ----------------- SURFACE GRID ---------------------------------------------
  # Auxiliary function for weighting values across span
  function calc_vals(span, array)

    # Finds bounding airfoil position
    val_in, val_out = nothing, array[1]
    for val in array[2:end]
        val_in = val_out
        val_out = val
        if val[1]>=span; break; end
    end
    pos_in = val_in[1]
    val_in = val_in[2]
    pos_out = val_out[1]
    val_out = val_out[2]

    weight = (span-pos_in)/(pos_out-pos_in)

    return weight, val_in, val_out
  end

  # Space transformation function
  function my_space_transform(inds)
    span = b_pos[inds[2]]           # y/bscale span position
    chord = chords[inds[2]]         # c/bscale chord length
    twist = twists[inds[2]]         # deg twist
    le_x = LE_x[inds[2]]            # x/bscale LE position
    le_z = LE_z[inds[2]]            # z/bscale LE position

    inds1 = (sec_NDIVS+1) - (inds[1]-1) # Here it flips the airfoil contour

    # Merges airfoil contours at this span position
    weight, sec_in, sec_out = calc_vals(span*sign(span)^symmetric, crosssections)

    # Point over airfoil contour
    point = weight*sec_out[inds1, :]+(1-weight)*sec_in[inds1, :]
    point = vcat(point, 0)

    # Scales the airfoil contour by the normalized chord length
    point = chord*point

    # Applies twist to the airfoil point
    tlt_z = tilt_z!=nothing ?  tilt_z[inds[2]] : 0.0
    tlt_y = tilt_y!=nothing ?  tilt_y[inds[2]] : 0.0
    point = rotation_matrix(-twist, -tlt_z, -tlt_y)*point

    # Places the point relative to LE and scales by span scale
    point = [point[1]+le_x, span+point[3], point[2]+le_z]*bscale

    return point
  end

  # Transforms the quasi-two dimensional parametric grid into the wing surface
  transform2!(grid, my_space_transform)

  if save_path!=nothing
    # Outputs a vtk file
    save(grid, file_name; path=save_path)

    if paraview
      # Calls paraview
      run(`paraview --data=$(joinpath(save_path,file_name)).vtk`)
    end
  end

  return grid::Grid
end

"""
`generate_loft(crosssections, bscale, b_low, b_up, b_NDIVS, chords, twists,
LE_x, LE_z; optargs...)`

  Generates a lofted surface geometry. The name of the arguments are taken from
  its initial application to the lofting of an aircraft wing, so bear with it.

  **Arguments**
  * `bscale::Float64`         : Semi-span scale. It will scale the entire
                                geometry by this factor. All y/bscale values in
                                the following arguments must go from 0 to 1.
  * `b_low::Float64`          : Scaled lower bound of the span.
  * `b_up::Float64`           : Scaled upper bound of the span. To generate
                                a symmetric wing, give it b_low=-1, b_up=1 and .
                                symmetric=true; for a semi-span, give it
                                b_low=0, b_up=1. If generating a prop blade,
                                give it b_low=Rhub/Rtip, b_up=1.
  * `b_NDIVS`                 : Number of divisions (cells) along span. This
                                if either an Int or an array of sections in the
                                format of `multidiscretize()`.
  * `chords::Array{Float64,2}`: Chord distribution along the span in the form
                                [y/bscale c/bscale].
  * `twists::Array{Float64,2}`: Twist distribution along the span in the form
                                [y/bscale deg].
  * `LE_x::Array{Float64,2}`  : x-position (chordwise) of leading edge along the
                                span in the form [y/bscale x/bscale].
  * `LE_z::Array{Float64,2}`  : z-position (dihedral-wise) of leading edge along
                                the span in the form [y/bscale z/bscale].
  * `crosssections::Array{Tuple{T,Array{T,2}}, 1}`    : cross sections along the
                                span in the form [(y/bscale, crosssection)],
                                where `crosssection` is a matrix that contains
                                all points of the airfoil contour indexed by
                                row. In order to define the resulting normals
                                pointing out of the geometry, the points in
                                `crosssection` must start from the trailing
                                edge, go around the top side towards the leading
                                edge, and back to the trailing edge around the
                                bottom side. ALL CROSS SECTIONS MUST HAVE THE
                                SAME NUMBER OF POINTS.

  **Optional Arguments**
  * `tilt_z::Array{Float64,2}`: Tilting about the z-axis of every span cross
                                section in the form [(y/bscale, deg)].
  * `spl_k`, `spl_bc`, `spl_s`: Spline parameters with k the degree of the
                                spline, bc the boundary condition, and s the
                                smoothing or error of the spline.

  NOTE: The resulting geometry will be a structured surface mesh with the first
  index going in chordwise direction starting from the TE around the bottom
  surface and around back to the TE. The second index is the in the spanwise
  direction starting at left end (lowest position along y axis) to right end.
"""
function generate_loft(crosssections::Array{Tuple{T,Array{T,2}}, 1},
                        bscale::Real, b_low::Real, b_up::Real, b_NDIVS,
                        chords::Array{T,2}, twists::Array{T,2},
                        LE_x::Array{T,2}, LE_z::Array{T,2};
                        # MORE GEOMETRIC OPTIONS
                        tilt_z=nothing, symmetric::Bool=false,
                        # SPLINE OPTIONS
                        spl_k::Int64=5, spl_bc::String="extrapolate",
                        spl_s::Real=0.001, verify_spline::Bool=true,
                        optargs...
                       ) where{T<:Real}

 # Normalized span positions
 if typeof(b_NDIVS)==Int64
   b_pos = discretize(x->x, b_low, b_up, b_NDIVS, 1.0)
 elseif typeof(b_NDIVS)==multidisctype
   b_pos = multidiscretize(x->x, b_low, b_up, b_NDIVS)
 else
   error("Invalid b_NDIVS=$b_NDIVS."*
          " Expected type $(Int64) or $(multidiscrtype), got $(typeof(b_NDIVS))")
 end
 b_pos = T[val for val in b_pos]

  # ----------------- GEOMETRY SPLINES -----------------------------------------
  # Splines all distributions for a smooth geometry
  _spl_chord = Dierckx.Spline1D(chords[:, 1], chords[:, 2]; s=spl_s, bc=spl_bc,
                          k=size(chords,1)>=spl_k ? spl_k : size(chords,1)-1 )
  _spl_twist = Dierckx.Spline1D(twists[:, 1], twists[:, 2];s=spl_s, bc=spl_bc,
                          k=size(twists,1)>=spl_k ? spl_k : size(twists,1)-1)
  _spl_LE_x = Dierckx.Spline1D(LE_x[:, 1], LE_x[:, 2]; s=spl_s, bc=spl_bc,
                          k=size(LE_x,1)>=spl_k ? spl_k : size(LE_x,1)-1)
  _spl_LE_z = Dierckx.Spline1D(LE_z[:, 1], LE_z[:, 2]; s=spl_s, bc=spl_bc,
                          k=size(LE_z,1)>=spl_k ? spl_k : size(LE_z,1)-1)
  if tilt_z!=nothing
   _spl_tlt_z = Dierckx.Spline1D(tilt_z[:, 1], tilt_z[:, 2]; s=spl_s, bc=spl_bc,
                          k=size(tilt_z,1)>=spl_k ? spl_k : size(tilt_z,1)-1)
    new_tlt_z = [_spl_tlt_z(symmetric ? abs(b) : b) for b in b_pos]
  else
    new_tlt_z = nothing
  end

  new_chords = [_spl_chord(symmetric ? abs(b) : b) for b in b_pos]
  new_twists = [_spl_twist(symmetric ? abs(b) : b) for b in b_pos]
  new_LE_x = [_spl_LE_x(symmetric ? abs(b) : b) for b in b_pos]
  new_LE_z = [_spl_LE_z(symmetric ? abs(b) : b) for b in b_pos]

  # ----------------- SPLINE VERIFICATION --------------------------------------
  if verify_spline
    if isdefined(Main, :PyPlot)
      plot_loft_splines(b_pos, LE_x, LE_z, new_LE_x, new_LE_z, twists, new_twists,
        tilt_z, new_tlt_z, chords, new_chords)
    else
      @eval Main import PyPlot
      Base.invokelatest(plot_loft_splines, b_pos, LE_x, LE_z, new_LE_x, new_LE_z, twists, new_twists,
        tilt_z, new_tlt_z, chords, new_chords)
    end

  end

  # ----------------- LOFTING --------------------------------------------------
  return generate_loft(crosssections, bscale, b_pos, new_chords, new_twists,
                        new_LE_x, new_LE_z; tilt_z=new_tlt_z,
                        symmetric=symmetric, optargs...)

end

"""
`generate_loft(crosssections, upper_rfl_NDIVS, lower_rfl_NDIVS, args...;
rflspl_k::Int64=5, rflspl_s::Real=0.001, verify_rflspline::Bool=true,
rfloptargs=[], optargs...)`

  This function also rediscretizes the cross sections as indicated by
  `upper_rfl_NDIVS`, and `lower_rfl_NDIVS` (upper and lower surface sections,
  respectively). Hence, the original cross sections need not to have the same
  number of points.
"""
function generate_loft(crosssections::Array{Tuple{T,Array{T,2}}, 1},
                            upper_rfl_NDIVS, lower_rfl_NDIVS,
                            bscale::Real, b_low::Real, b_up::Real, b_NDIVS,
                            chords::Array{T,2}, twists::Array{T,2},
                            LE_x::Array{T,2}, LE_z::Array{T,2};
                            # AIRFOIL SPLINE OPTIONS
                            rflspl_k::Int64=5, rflspl_s::Real=0.001,
                            verify_rflspline::Bool=true, rfloptargs=[],
                            optargs...
                           ) where{T<:Real}

  # Formats the NDIVS
  if typeof(upper_rfl_NDIVS)==Int64
    urfl_NDIVS = [(1.0, upper_rfl_NDIVS, 1.0, false)]
  elseif typeof(upper_rfl_NDIVS)==multidisctype
    urfl_NDIVS = upper_rfl_NDIVS
  else
    error("Invalid upper_rfl_NDIVS=$upper_rfl_NDIVS. Expected type $(Int64)"*
            " or $(multidiscrtype), got $(typeof(upper_rfl_NDIVS))")
  end
  if typeof(lower_rfl_NDIVS)==Int64
    lrfl_NDIVS = [(1.0, lower_rfl_NDIVS, 1.0, false)]
  elseif typeof(lower_rfl_NDIVS)==multidisctype
    lrfl_NDIVS = lower_rfl_NDIVS
  else
    error("Invalid lower_rfl_NDIVS=$lower_rfl_NDIVS. Expected type $(Int64)"*
            " or $(multidiscrtype), got $(typeof(lower_rfl_NDIVS))")
  end

  # Rediscretizes cross sections
  rediscr = []

  for (pos, M) in crosssections
    new_xy = rediscretize_airfoil(M[:, 1], M[:, 2], urfl_NDIVS, lrfl_NDIVS;
                                    spl_k=rflspl_k, spl_s=rflspl_s,
                                    verify_spline=verify_rflspline,
                                    rfloptargs...)
    push!(rediscr, [new_xy[j][i] for i in 1:size(new_xy[1],1), j in 1:2])
  end

  # Formates them back to the original format
  new_crosssections = [ (pos, rediscr[i])
                                  for (i,(pos, M)) in enumerate(crosssections)]

  return generate_loft(new_crosssections, bscale, b_low, b_up, b_NDIVS,
                          chords, twists, LE_x, LE_z; optargs...)
end

"""
`generate_loft(crosssections::Array{Tuple{T,String}, 1},
                        data_path::String, args...; header_len::Int64=1,
                        delim::String=" ", optargs...) where{T<:Real}`

Loft a geometry where the cross sections are read from the files indicated by
`crosssections` found in `data_path`.
"""
function generate_loft(crosssections::Array{Tuple{T,String}, 1},
                        data_path::String, args...; header_len::Int64=1,
                        delim::String=" ", optargs...) where{T<:Real}

 secs = [(pos, readcontour(f_name; path=data_path, header_len=header_len,
                              delim=delim, output="matrix"))
                              for (pos, f_name) in crosssections]

 return generate_loft(secs, args...; optargs...)
end

# REVOLTUION ###################################################################
"""
  `surface_revolution(profile, thetaNDIVS; loop_dim=0, axis_angle=0, low_a=0,
up_a=360, save_path=nothing, paraview=true, file_name="myrev")`

Receives a contour to revolve around an axis generating a surface of revolution.

**ARGUMENTS**
  * `profile::Array{Float64,2}`   : Contour to revolve. These are two-
                                    dimensional points in the YZ-plane that will
                                    get revolve around the Z-axis.
  * `thetaNDIVS::Int64`           : Number of angle-sections (cells) in the
                                    revolution.

**OPTIONAL ARGUMENTS**
  * `loop_dim::Int64=0`           : Whether to loop any dimension of the
                                    parametric grid.
  * `axis_angle::Float64=0`       : Tilting angle (deg) about the Z-axis to
                                    revolve the contour.
  * `low_a::Float64=0`            : Lower bound of angle (deg) of revolution.
  * `up_a::Float64=0`             : Upper bound of angle (deg) of revolution.
"""
function surface_revolution(profile::AbstractMatrix{T}, thetaNDIVS::Integer;
                              loop_dim::Integer=0, axis_angle::Number=0,
                              low_a::Number=0, up_a::Number=360,
                              # OUTPUT OPTIONS
                              save_path=nothing, paraview::Bool=true,
                              file_name::AbstractString="myrev"
                              ) where{T<:Real}
  #ERROR CASES
  if size(profile,2)!=2
    error("Invalid point dimensions in `profile`."*
          " Expected 2 dimensions, got $(size(profile,2))")
  elseif profile[1,:]==profile[end,:] && loop_dim==0
    @warn("Received a closed contour but parametric grid wasn't declared to"*
          " loop, resulting in overlaping start/end points. Give it"*
          " `loop_dim=1` to fix that.")
  end

  # First coordinate is a map to the points in the countor of revolution
  # Seciond coordinate is the angle of revolution
  # Third is a dummy
  NDIVS = [size(profile,1)-1, thetaNDIVS, 0]# Divisions in every coordinate
  P_min = [0, low_a, 0]                     # Lower boundaries
  P_max = [1, up_a, 0 ]                     # Upper boundaries
  # loop_dim = 1                            # Loops the countour of revolution

  # Parametric grid
  grid = Grid(P_min, P_max, NDIVS, loop_dim)

  # Rotates the profile according the revolution tilting
  if axis_angle!=0
    M = rotation_matrix(0,0,-axis_angle)
    M2D = M[2:3,2:3]
    M3D = collect(M)'
    points = collect(hcat([M2D*profile[i,:] for i in 1:size(profile,1)]...))'
  else
    M3D = I
    points = profile
  end

  # Space transformation function
  function my_space_transform(X, ind)
      p_ind = ind[1]                  # Point index
      angle = X[2]                    # (deg) angle

      # Places the contour in the YZ-plane
      point = [0, points[p_ind,1], points[p_ind,2]]

      # Rotates the contour around the Z-axis
      return M3D*axis_rotation(Float64[0,0,1], Float64(angle))*point
  end

  # Applies the space transformation to the parametric grid
  transform3!(grid, my_space_transform)

  if save_path!=nothing
    # Outputs a vtk file
    save(grid, file_name; path=save_path)

    if paraview
      # Calls paraview
      run(`paraview --data=$(joinpath(save_path,file_name)).vtk`)
    end
  end

  return grid::Grid
end

################################################################################


text:
# Perform Sanity Tests for Agents, Arm models and Experiments

@testset "Sanity Checks" begin

# ---------------------------------- Tests for Agents --------------------------------------
# Get all algorithms
list_algorithms     = subtypes( MAB.BanditAlgorithmBase )

basefn__get_arm_index   = which( get_arm_index,
                                    Tuple{MAB.Algorithms.BanditAlgorithmBase} )
basefn__update_reward   = which( update_reward!,
                                    Tuple{MAB.Algorithms.BanditAlgorithmBase,Real} )
basefn__reset           = which( reset!,
                                    Tuple{MAB.Algorithms.BanditAlgorithmBase} )
basefn__info_str        = which( info_str,
                                    Tuple{MAB.Algorithms.BanditAlgorithmBase,Bool} )

@testset "Algorithms" begin
    for a ∈ list_algorithms
        @testset "$(a)" begin
            # Test whether each of the algorithm implements all the required methods
            fn__get_arm_index   = which( get_arm_index, Tuple{a} )
            fn__update_reward   = which( update_reward!, Tuple{a,Real} )
            fn__reset           = which( reset!, Tuple{a} )
            fn__info_str        = which( info_str, Tuple{a,Bool} )
            @test fn__get_arm_index != basefn__get_arm_index
            # Implement atleast on of these methods
            @test basefn__update_reward != which( update_reward!, Tuple{a,Real} ) ||
                  basefn__update_reward != which( update_reward!, Tuple{a,Int} ) ||
                  basefn__update_reward != which( update_reward!, Tuple{a,Float64} )
            @test fn__reset != basefn__reset
            @test fn__info_str != basefn__info_str
        end
    end
end

# -------------------------------- Tests for Arm Models ------------------------------------
list_arms           = subtypes( MAB.Arms.BanditArmBase )

basefn__pull!   = which( Arms.pull!, Tuple{Arms.BanditArmBase} )
basefn__tick!   = which( Arms.tick!, Tuple{Arms.BanditArmBase} )
basefn__reset!  = which( Arms.reset!, Tuple{Arms.BanditArmBase} )

@testset "Arms" begin
    for arm ∈ list_arms
        @testset "$(arm)" begin
        fn__pull!   = which( Arms.pull!, Tuple{arm} )
        fn__tick!   = which( Arms.tick!, Tuple{arm} )
        fn__reset!  = which( Arms.reset!, Tuple{arm} )

        @test basefn__pull! != fn__pull!
        @test basefn__tick! != fn__tick!
        @test basefn__reset!    != fn__reset!
        end
    end
end

# ------------------------------- Tests for Experiments ------------------------------------
list_experiments    = subtypes( MAB.Experiments.BanditExpBase )
#TODO

end

ext: jl


text:
"""
    macro stack(into, structs...)

Stack the fields of multiple structs and create a new one.
The first argument is the name of the new struct followed by the
ones to be stacked. Parametric types are not supported and
the fieldnames needs to be unique.

Example:

    @stack Baz Foo Bar

Creates `Baz` with the concatenated fields of `Foo` and `Bar`

"""
macro stack(into, structs...)
    fields = []
    for _struct in structs
        names = fieldnames(getfield(__module__, _struct))
        types = [fieldtype(getfield(__module__, _struct), field) for field in names]
        for (n, t) in zip(names, types)
            push!(fields, :($n::$t))
        end
    end

    esc(
        quote
            struct $into
                $(fields...)
            end
        end
    )
end

"""
    unpack(x::CompressionHeader)

Return the following information:
- Name of compression algorithm
- Level of the compression
- compressedbytes and uncompressedbytes according to [uproot3](https://github.com/scikit-hep/uproot3/blob/54f5151fb7c686c3a161fbe44b9f299e482f346b/uproot3/source/compressed.py#L132)
"""
function unpack(x::CompressionHeader)
    algname = String(x.algo)
    # shift without casting to `Int` will give you 0x00 because we're shifting 0 bits into UInt8
    compressedbytes = x.c1 + (Int(x.c2) << 8) + (Int(x.c3) << 16)
    uncompressedbytes = x.u1 + (Int(x.u2) << 8) + (Int(x.u3) << 16)

    return algname, x.method, compressedbytes, uncompressedbytes
end

abstract type JaggType end
struct Nojagg      <:JaggType  end
struct Nooffsetjagg<:JaggType  end
struct Offsetjagg  <:JaggType  end
struct Offsetjaggjagg  <:JaggType  end

function JaggType(f, branch, leaf)
    # https://github.com/scikit-hep/uproot3/blob/54f5151fb7c686c3a161fbe44b9f299e482f346b/uproot3/interp/auto.py#L144
    (match(r"\[.*\]", leaf.fTitle) !== nothing) && return Nooffsetjagg
    leaf isa TLeafElement && leaf.fLenType==0 && return Offsetjagg
    !hasproperty(branch, :fClassName) && return Nojagg

    try
        streamer = streamerfor(f, branch.fClassName).streamer.fElements.elements[1]
        (streamer.fSTLtype == Const.kSTLvector) && return Offsetjagg
    catch
    end

    return Nojagg
end

"""
    parseTH(th::Dict{Symbol, Any})

Parse the output of [`TH`](@ref) into a tuple of `counts`, `edges`, and `sumw2`.
A `StatsBase.Histogram` can then be constructed with `Histogram(edges, counts)`.
TH1 and TH2 inputs are supported.
"""
function parseTH(th::Dict{Symbol, Any})
    xmin = th[:fXaxis_fXmin]
    xmax = th[:fXaxis_fXmax]
    xnbins = th[:fXaxis_fNbins]
    xbins = isempty(th[:fXaxis_fXbins]) ? range(xmin, xmax, length=xnbins+1) : th[:fXaxis_fXbins];
    counts = th[:fN]
    sumw2 = th[:fSumw2]
    if th[:fYaxis_fNbins] > 1
        ymin = th[:fYaxis_fXmin]
        ymax = th[:fYaxis_fXmax]
        ynbins = th[:fYaxis_fNbins]
        ybins = isempty(th[:fYaxis_fXbins]) ? range(ymin, ymax, length=ynbins+1) : th[:fYaxis_fXbins];
        counts = reshape(counts, (xnbins+2, ynbins+2))[2:end-1, 2:end-1]
        sumw2 = reshape(sumw2, (xnbins+2, ynbins+2))[2:end-1, 2:end-1]
        edges = (xbins, ybins)
    else
        counts = counts[2:end-1]
        sumw2 = sumw2[2:end-1]
        edges = (xbins,)
    end
    return counts, edges, sumw2
end

function samplefile(filename::AbstractString)
    return ROOTFile(normpath(joinpath(@__DIR__, "../test/samples", filename)))
end


text:
# src path
src_path="C:\\Users\\Koike\\Desktop\\git\\general_2d_euler\\src\\"

# main (変更しないこと)
src_read="read_grid.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="cal_time_step.jl"
include(src_path*src_read)
src_read="converge.jl"
include(src_path*src_read)
src_read="output.jl"
include(src_path*src_read)
src_read="setup.jl"
include(src_path*src_read)
src_read="boundary.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_main="main.jl"
include(src_path*src_main)


text:
for (fn) in (:floor, :ceil, :round)
  @eval begin
    function ($fn)(a::Float120)
        hi = ($fn)(a.hi)
        lo = 0.0
        if (hi == a.hi)
            lo = ($fn)(a.lo)
            hi,lo = eftAddGTE(hi,lo)
        end
        DD(hi,lo)
    end
  end
end

function (trunc)(a::Float120)
    ifelse(a.hi < zero(Float64), ceil(a), floor(a))
end

"""
stretch is the opposite of trunc()
it extends to the nearest integer away from zero
"""
function (stretch)(a::Float120)
    ifelse(a.hi < zero(Float64), floor(a), ceil(a))
end


text:
module GLExtensions

using ModernGL

using MathManager
using GLDebug
using GLExtender

########################################################

GLExtender.openLib()
@glfunc glBindBuffersBase(target::GLenum, first::GLuint, count::GLsizei, buffers::Ptr{GLuint})::Cvoid
@glfunc glBindBuffersRange(target::GLenum, first::GLuint, count::GLsizei, buffers::Ptr{GLuint}, offsets::Ptr{GLintptr}, sizes::Ptr{GLintptr})::Cvoid
GLExtender.closeLib()

export glBindBuffersBase
export glBindBuffersRange

################################################################################

""" TODO """
function validateShader(shader)
  success = GLint[0]
  glGetShaderiv(shader, GL_COMPILE_STATUS, success)
  success[] == GL_TRUE
end
export validateShader

#loadShaderSource(shaderID::GLuint, source::String) = (shadercode=Vector{UInt8}(string(source,"\x00")); glShaderSource(shaderID, 1, Ptr{UInt8}[pointer(shadercode)], Ref{GLint}(length(shadercode))))
#glShaderSource(shader, 1, convert(Ptr{UInt8}, pointer([convert(Ptr{GLchar}, pointer(source))])), C_NULL)

""" TODO """
function compileShader(shader::GLuint, source::String)
  @GLCHECK glShaderSource(shader, 1, convert(Ptr{UInt8}, pointer([convert(Ptr{GLchar}, pointer(string(source,"\x00")))])), C_NULL)
  @GLCHECK glCompileShader(shader)
  validateShader(shader)
end
export compileShader

################################################################################

""" TODO """
glArrayFunc(ids::Array, func::Function) = begin len=length(ids); if len>0 func(len, ids) end; end
#glCreate(count::Number, arr::AbstractArray, func::Function, args...) = for i=1:count push!(arr, func(args...)) end
#glDelete(ids::Array, func::Function, args...) = for id in ids func(id, args...) end
export glArrayFunc

""" TODO """
glGenOne(glGenFn) = begin id=GLuint[0]; @GLCHECK glGenFn(1, id); id[] end
export glGenOne

""" TODO """
glDelOne(glDelFn, id) = begin ids=GLuint[id]; @GLCHECK glDelFn(1, ids) end
export glDelOne

""" TODO """
glGenBuffer() = glGenOne(glGenBuffers)
export glGenBuffer

""" TODO """
glDeleteBuffer(id) = glDelOne(glDeleteBuffers,id)
export glDeleteBuffer

""" TODO """
glGenFramebuffer() = glGenOne(glGenFramebuffers)
export glGenFramebuffer

""" TODO """
glDeleteFramebuffer(id) = glDelOne(glDeleteFramebuffers,id)
export glDeleteFramebuffer

""" TODO """
glGenRenderbuffer() = glGenOne(glGenRenderbuffers)
export glGenRenderbuffer

""" TODO """
glDeleteRenderbuffer(id) = glDelOne(glDeleteRenderbuffers,id)
export glDeleteRenderbuffer

""" TODO """
glGenVertexArray() = glGenOne(glGenVertexArrays)
export glGenVertexArray

""" TODO """
glDeleteVertexArray(id) = glDelOne(glDeleteVertexArrays,id)
export glDeleteVertexArray

""" TODO """
glGenTexture() = glGenOne(glGenTextures)
export glGenTexture

""" TODO """
glDeleteTexture(id) = glDelOne(glDeleteTextures, id)
export glDeleteTexture

""" TODO """
glGenSampler() = glGenOne(glGenSamplers)
export glGenSampler

""" TODO """
glDeleteSampler(id) = glDelOne(glDeleteSamplers, id)
export glDeleteSampler

################################################################################
export glGetValue

""" TODO """
glGetValue(f::Function, typ::DataType, t::Tuple) = begin value=typ[0]; @GLCHECK f(t...,value); value[] end

""" TODO """
function glGetValue(typ::DataType, name::GLenum, index::Union{Nothing, GLuint}=nothing)

  noindex = index == nothing
  t = noindex ? (name,) : (name,index)
	f = nothing

	if typ == GLboolean f = noindex ? glGetBooleanv : glGetBooleani_v
	elseif typ == GLint f = noindex ? glGetIntegerv : glGetIntegeri_v
	elseif typ == GLint64 f = noindex ? glGetInteger64v : glGetInteger64i_v
	elseif typ == GLfloat f = noindex ? glGetFloatv : glGetFloati_v
	elseif typ == GLdouble f = noindex ? glGetDoublev : glGetDoublei_v
  elseif typ == String return unsafe_string((noindex ? glGetString : glGetStringi)(t...))
	end

	if f != nothing return glGetValue(f, typ, t) end
	nothing
end

################################################################################

""" TODO """
glSetFragLocation(id::GLuint, name::String, colorNumber::GLuint) = glFragLocation(id, name, colorNumber)
export glSetFragLocation

""" TODO """
glFragLocation(id::GLuint, name::String, colorNumber::GLuint) = glBindFragDataLocation(id, colorNumber, name)
export glFragLocation

""" TODO """
glAttribLocation(id::GLuint, name::String) = glGetAttribLocation(id, name)
export glAttribLocation

""" TODO """
glUniformLocation(id::GLuint, name::String) = glGetUniformLocation(id, name)
export glUniformLocation

################################################################################

export glUniform

""" TODO """
glUniform(location::GLint, value::Bool) = glUniform1ui(location, GLuint(value ? 1 : 0))

""" TODO """
glUniform(location::GLint, value::UInt32) = glUniform1ui(location, GLuint(value))

""" TODO """
glUniform(location::GLint, value::UInt64) = glUniform1ui(location, GLuint(value))

""" TODO """
glUniform(location::GLint, value::Int32) = glUniform1i(location, GLint(value))

""" TODO """
glUniform(location::GLint, value::Int64) = glUniform1i(location, GLint(value))

""" TODO """
glUniform(location::GLint, value::Float32) = glUniform1f(location, GLfloat(value))

""" TODO """
glUniform(location::GLint, value::Float64) = glUniform1d(location, GLdouble(value))

""" TODO """
glUniform(location::GLint, value::Array{Float32,1}) = glUniform1fv(location, length(value), pointer(value))

""" TODO """
glUniform(location::GLint, value::Array{Float64,1}) = glUniform1dv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{2, Float32},1}) = glUniform2fv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{3, Float32},1}) = glUniform3fv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{4, Float32},1}) = glUniform4fv(location, length(value), pointer(value))

""" TODO """
glUniform(location::GLint, value::Vec2f) = glUniform2fv(location, 1, pointer(convert(Array, value)))

""" TODO """
glUniform(location::GLint, value::Vec3f) = glUniform3fv(location, 1, pointer(convert(Array, value)))

""" TODO """
glUniform(location::GLint, value::Vec4f) = glUniform4fv(location, 1, pointer(convert(Array, value)))

""" TODO """
glUniform(location::GLint, value::Mat2x2f, transpose=false) = glUniformMatrix2fv(location, 1, transpose, pointer(convert(Array, value)))

""" TODO """
glUniform(location::GLint, value::Mat3x3f, transpose=false) = glUniformMatrix3fv(location, 1, transpose, pointer(convert(Array, value)))

""" TODO """
glUniform(location::GLint, value::Mat4x4f, transpose=false) = glUniformMatrix4fv(location, 1, transpose, pointer(convert(Array, value)))

################################################################################

#glGetSamplerParameterIiv
#glGetPointerv
#glGetQueryObjectui64v
#glGetProgramiv
#glGetBufferPointerv
#glGetTexParameterIuiv
#glGetShaderiv
#glGetRenderbufferParameteriv
#glGetUniformuiv
#glGetProgramPipelineiv
#glGetVertexAttribfv
#glGetVertexAttribLdv
#glGetVertexAttribiv
#glGetSamplerParameteriv
#glGetActiveUniformBlockiv
#glGetVertexAttribdv
#glGetQueryObjectuiv
#glGetProgramResourceiv
#glGetUniformfv
#glGetUniformdv
#glGetProgramInterfaceiv
#glGetVertexAttribIuiv
#glGetFramebufferParameteriv
#glGetActiveSubroutineUniformiv
#glGetVertexAttribPointerv
#glGetBufferParameteriv
#glGetUniformiv
#glGetQueryObjecti64v
#glGetTexLevelParameteriv
#glGetBufferParameteri64v
#glGetQueryObjectiv
#glGetActiveUniformsiv
#glGetTexParameterfv
#glGetTexLevelParameterfv
#glGetVertexAttribIiv
#glGetFramebufferAttachmentParameteriv
#glGetActiveAtomicCounterBufferiv
#glGetSynciv
#glGetSamplerParameterfv
#glGetQueryiv
#glGetTexParameterIiv
#glGetUniformSubroutineuiv
#glGetQueryIndexediv
#glGetProgramStageiv
#glGetSamplerParameterIuiv
#glGetTexParameteriv
#glGetMultisamplefv

end #GLExtensions

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module USDAQuickStats

import HTTP: request

export
    set_api_key,
    get_counts,
    get_param_values,
    get_nass

const usda_url = "http://quickstats.nass.usda.gov"

include("set_api_key.jl")
include("nass.jl")
include("counts.jl")
include("param_values.jl")

end # module


text:
### Tests for each OPF objective ###

@testset "mincost objective" begin
    base_data = "data/case14/base"
    result = VariableFrequencyOPF.multifrequency_opf(base_data, "mincost", output_to_files=false)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 8081.5; atol = 1e0)

    variablef_data = "data/case14_twoarea/two_area/"
    result = VariableFrequencyOPF.multifrequency_opf(variablef_data, "mincost", output_to_files=false)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 7565.2; atol = 1e0)
    @test isapprox(result[1]["cost"][2], 553.6; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][1], 1.0; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][2], 60.0; atol = 1e0)

    base_data = "data/case14/base"
    result = VariableFrequencyOPF.multifrequency_opf(base_data, "mincost", output_to_files=true)
    result_output_dict = JSON.parsefile("results/case14/base/output_values.json")
    @test result_output_dict["status"] == "LOCALLY_SOLVED"
    @test isapprox(result_output_dict["cost"][1], 8081.5; atol = 1e0)
    try
        rm("results",force=true,recursive=true)
    catch IOError
        println("Unable to delete the results files generated by this test.")
    end
end

@testset "mincost objective_regularized" begin
    variablef_data = "data/case14_twoarea/two_area/"
    result = VariableFrequencyOPF.multifrequency_opf(variablef_data, "mincost", output_to_files=false, regularize_f=1e-4)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 7565.2; atol = 1e0)
    @test isapprox(result[1]["cost"][2], 553.6; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][1], 1.0; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][2], 60.0; atol = 1e0)

    variablef_data = "data/case14_twoarea/two_area/"
    result = VariableFrequencyOPF.multifrequency_opf(variablef_data, "mincost", output_to_files=false, regularize_f=1e-2)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 7565.2; atol = 1e0)
    @test isapprox(result[1]["cost"][2], 553.6; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][1], 1.0; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][2], 60.0; atol = 1e0)

    variablef_data = "data/case14_twoarea/two_area/"
    result = VariableFrequencyOPF.multifrequency_opf(variablef_data, "mincost", output_to_files=false, regularize_f=1e0)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 7565.2; atol = 1e0)
    @test isapprox(result[1]["cost"][2], 553.6; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][1], 1.0; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][2], 60.0; atol = 1e0)

    variablef_data = "data/case14_twoarea/two_area/"
    result = VariableFrequencyOPF.multifrequency_opf(variablef_data, "mincost", output_to_files=false, regularize_f=1e3)
    @test result[1]["status"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(result[1]["cost"][1], 7345.4; atol = 1e0)
    @test isapprox(result[1]["cost"][2], 825.7; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][1], 60.0; atol = 1e0)
    @test isapprox(result[1]["frequency (Hz)"][2], 60.0; atol = 1e0)

    try
        rm("results",force=true,recursive=true)
    catch IOError
        println("Unable to delete the results files generated by this test.")
    end
end

@testset "areagen objective" begin
    original_network = "data/nordic_fault/base/fault_4032_4044.m"
    # We choose to put the new network data in the same folder because
    # we also want to use some data that is already there. Each upgrade
    # gets its own folder.
    new_data_directory = "results/inputs/nordic_fault/"
    standard_frequency = 50.0
    lfac_branch_upgrades = [27,33]

    VariableFrequencyOPF.upgrade_branches(
        original_network,
        new_data_directory,
        standard_frequency,
        indices=lfac_branch_upgrades
    )

    # Minimize generation in areas 2 and 3 with plots including
    # the power flow between areas 1 and 2
    objective = "areagen"
    gen_areas = [2,3]
    area_transfer = [1,2]
    # We want to plot the results of all branch upgrades,
    # so we set `enum_branch` to true
    enum_branches = true

    solution =  VariableFrequencyOPF.run_series(
        new_data_directory,
        objective;
        gen_areas=gen_areas,
        area_transfer=area_transfer,
        enum_branches=enum_branches
    )
    results_dict = solution[1]
    n_subnets = solution[2]
    subnet_array = solution[3]
    idx_sorted = solution[4]
    series_output_folder = solution[5]
    plot_best_x = solution[6]

    subnet1 = findall(results_dict["subnet"]["br27"].==1)[1]
    subnet2 = 3-subnet1

    @test results_dict["status"]["br27"] == MathOptInterface.LOCALLY_SOLVED
    @test results_dict["status"]["br33"] == MathOptInterface.LOCALLY_SOLVED
    @test isapprox(results_dict["cost"]["br27"][subnet1], 4.40154e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br27"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet1], 50.0; atol = 1e0)
    # @test isapprox(results_dict["frequency (Hz)"]["br27"][subnet2], 20.92; atol = 1e0)
    @test isapprox(results_dict["cost"]["br33"][subnet1], 4.36446e6; rtol = 1e-2)
    @test isapprox(results_dict["cost"]["br33"][subnet2], 0.0; rtol = 1e-2)
    @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet1], 50.0; atol = 1e0)
    # @test isapprox(results_dict["frequency (Hz)"]["br33"][subnet2], 22.43; atol = 1e0)

    results_dict_allplots = [results_dict]

    plot_output_folder = "results/nordic_fault/"
    series_labels = ["LFAC upgrades"]

    # VariableFrequencyOPF.plot_results_dicts_bar(
    #     results_dict_allplots,
    #     n_subnets,
    #     subnet_array,
    #     idx_sorted,
    #     plot_output_folder,
    #     plot_best_x,
    #     series_labels,
    #     color_palette=:Dark2_8
    # )
    #
    # @test isfile("results/nordic_fault/plots/total/total generation.pdf")
    # @test isfile("results/nordic_fault/plots/subnet1/subnet 1 line loss.pdf")
    # @test isfile("results/nordic_fault/plots/subnet2/subnet 2 frequency (Hz).pdf")

    try
        rm("results",force=true,recursive=true)
    catch IOError
        println("Unable to delete the results files generated by this test.")
    end
end

# @testset "frequency sweep" begin
#     original_network = "data/nordic_fault/base/fault_4032_4044.m"
#     # We choose to put the new network data in the same folder because
#     # we also want to use some data that is already there. Each upgrade
#     # gets its own folder.
#     new_data_directory = "results/inputs/nordic_fault/"
#     standard_frequency = 50.0
#     lfac_branch_upgrades = [27]
#     new_data_file = new_data_directory*"/br27"
#
#     VariableFrequencyOPF.upgrade_branches(
#         original_network,
#         new_data_directory,
#         standard_frequency,
#         indices=lfac_branch_upgrades
#     )
#
#     # Minimize generation in areas 2 and 3 with plots including
#     # the power flow between areas 1 and 2
#     objective = "areagen"
#     fmin = 0:10:30
#     fmax = 0:10:30
#     results_dicts = Array{Dict}(undef, 1)
#     output_plots = Array{Any}(undef, 1)
#     (results_dicts[1], output_plots[1]) = VariableFrequencyOPF.frequency_ranges(
#         fmin,
#         fmax,
#         2,
#         new_data_file,
#         objective,
#         [("frequency (Hz)",2)],
#         [
#         ("nw",2,"branch",1,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",2,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",3,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",4,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",5,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",5,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",6,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",7,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",8,["pt","pf","qt","qf"])
#         ];
#         gen_areas=[2,3],
#         area_transfer=[1,2],
#         no_converter_loss=true
#     )
#
#     @test isfile("results/results/nordic_fault/br27/plots/vs subnet 2 frequency (Hz)/total/total generation in areas 2 and 3 (p.u.) vs subnet 2 frequency (Hz).pdf")
#     @test isfile("results/results/nordic_fault/br27/plots/vs subnet 2 frequency (Hz)/subnet1/subnet 1 line loss vs subnet 2 frequency (Hz).pdf")
#     @test isfile("results/results/nordic_fault/br27/plots/vs subnet 2 frequency (Hz)/subnet2/subnet 2 converter loss vs subnet 2 frequency (Hz).pdf")
#     @test isfile("results/nordic_fault/multiterminal_sweep/plots/vs subnet 2 frequency (Hz)/nw/2/branch/1/nw 2 branch 1 pf vs subnet 2 frequency (Hz).pdf")
#
#     # Now test the direct application of a data dictionary
#     output_results_folder = "time_0"
#     mn_data = VariableFrequencyOPF.read_sn_data(new_data_file, no_converter_loss=false)
#
#     folder = abspath(new_data_file)
#     folder_split = splitpath(folder)
#     toplevels = folder_split[1:end-3]
#     output_folder = joinpath(toplevels...,"results/$(folder_split[end-1])/$(folder_split[end])")
#     output_folder = joinpath(output_folder, output_results_folder)
#     if !isdir(output_folder)
#         mkpath(output_folder)
#     end
#
#     (results_dicts[1], output_plots[1]) = VariableFrequencyOPF.frequency_ranges(
#         fmin,
#         fmax,
#         2,
#         mn_data,
#         output_folder,
#         objective,
#         [("frequency (Hz)",2)],
#         [
#         ("nw",2,"branch",1,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",2,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",3,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",4,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",5,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",5,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",6,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",7,["pt","pf","qt","qf"]),
#         ("nw",2,"branch",8,["pt","pf","qt","qf"])
#         ];
#         gen_areas=[2,3],
#         area_transfer=[1,2]
#     )
#     @test isfile("results/results/nordic_fault/br27/time_0/plots/vs subnet 2 frequency (Hz)/total/total generation in areas 2 and 3 (p.u.) vs subnet 2 frequency (Hz).pdf")
#     @test isfile("results/results/nordic_fault/br27/time_0/plots/vs subnet 2 frequency (Hz)/subnet1/subnet 1 line loss vs subnet 2 frequency (Hz).pdf")
#     @test isfile("results/results/nordic_fault/br27/time_0/plots/vs subnet 2 frequency (Hz)/subnet2/subnet 2 converter loss vs subnet 2 frequency (Hz).pdf")
#
#     try
#         rm("results",force=true,recursive=true)
#     catch IOError
#         println("Unable to delete the results files generated by this test.")
#     end
# end


text:
using StaticArrays
n,m = 10,5
x = @SVector rand(n)
u = @SVector rand(m)
dt = 0.1
z = KnotPoint(x,u,dt)
RobotDynamics.set_state!(z, 2x)
RobotDynamics.set_control!(z, 2u)
@test state(z) ≈ 2x
@test control(z) ≈ 2u
RobotDynamics.set_z!(z, 3*[x;u])
@test state(z) ≈ 3x
@test control(z) ≈ 3u
@test z isa KnotPoint
@test z isa RobotDynamics.GeneralKnotPoint{T,N,M,SVector{NM,T}} where {T,N,M,NM}

# Terminal time step
z = KnotPoint(x,u,dt)
@test RobotDynamics.is_terminal(z) == false
@test RobotDynamics.get_z(z) ≈ [x;u]
z.dt = 0.0
@test RobotDynamics.is_terminal(z) == true
@test RobotDynamics.get_z(z) ≈ x

zterm = KnotPoint(x, m, 1.2)
@test control(zterm) ≈ zero(u)
@test state(zterm) ≈ x
@test RobotDynamics.is_terminal(zterm)

# GeneralKnotPoint
z_ = Vector([x;u])
z = RobotDynamics.GeneralKnotPoint(n,m,z_,dt)
@test z isa RobotDynamics.GeneralKnotPoint{Float64,n,m,Vector{Float64}}
@test RobotDynamics.get_z(z) isa Vector
@test state(z) isa SVector{n}
@test control(z) isa SVector{m}

# StaticKnotPoint
z = StaticKnotPoint(x, u, dt, 1.2)
@test_throws ErrorException RobotDynamics.set_state!(z, x)
@test_throws ErrorException RobotDynamics.set_control!(z, u)
@test_throws ErrorException RobotDynamics.set_z!(z, [x;u])

z0 = KnotPoint(2x,2u,dt,1.3)
z = StaticKnotPoint(z0, [x;u])
@test state(z) ≈ x
@test control(z) ≈ u
@test z.dt ≈ dt
@test z.t ≈ 1.3
z = StaticKnotPoint(z0)
@test state(z) ≈ 2x
@test control(z) ≈ 2u

z = z0 + z_
@test state(z) ≈ 3x
@test control(z) ≈ 3u
@test z isa StaticKnotPoint
@test z.dt ≈ z0.dt
@test z.t ≈ z0.t
z = z0 + 2*z0
@test state(z) ≈ 6x
@test control(z) ≈ 6u
@test z isa StaticKnotPoint
@test z.dt ≈ z0.dt
@test z.t ≈ z0.t

z *= 10
@test state(z) ≈ 60x
@test control(z) ≈ 60u


text:
# Scatter
# begin
#     brd = board("brd", xlim=[0,1], ylim=[0,1])
#     scatter(rand(5), rand(5), strokecolor="blue", withlabel=false) |> brd
#     brd
# end

# Scatter + fun
# begin
#     brd = board("brd", xlim=[0,1], ylim=[0,1])
#     slider("a", [[0.1,0.1],[0.6,0.1],[0,0.2,1]]) |> brd
#     @jsf foo() = val(a)
#     scatter(rand(5), foo, strokecolor="blue", withlabel=false) |> brd
#     brd
# end


text:
@testset "undirected_paths" begin
    G = SimpleGraph(Edge.([(1, 2), (2, 3), (2, 4), (3, 4)]))
    expected = Set([[1, 2, 4], [1, 2, 3, 4]])
    @test Set(Causality.undirected_paths(G, 1, 4)) == expected
end


text:
function Annotate(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)
    ccall((:ImPlot_AnnotateStr, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)
    ccall((:ImPlot_AnnotateClampedStr, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function Annotate(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)
    ccall((:ImPlot_AnnotateVec4, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)
    ccall((:ImPlot_AnnotateClampedVec4, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end


text:
# Load docstring from markdown files:
for (name, path) in [
    :SplittablesBase => joinpath(dirname(@__DIR__), "README.md"),
    :halve => joinpath(@__DIR__, "halve.md"),
    :amount => joinpath(@__DIR__, "amount.md"),
]
    try
        include_dependency(path)
        str = read(path, String)
        @eval @doc $str $name
    catch err
        @error "Failed to import docstring for $name" exception = (err, catch_backtrace())
    end
end

"""
    _unzip(xs::Tuple)

# Examples
```jldoctest; setup = :(using SplittablesBase.Implementations: _unzip)
julia> _unzip(((1, 2, 3), (4, 5, 6)))
((1, 4), (2, 5), (3, 6))
```
"""
_unzip(xs::Tuple{Vararg{NTuple{N,Any}}}) where {N} = ntuple(i -> map(x -> x[i], xs), N)

tail(xs) = _tail(xs...)
_tail(_, args...) = args

"""
    arguments(xs)

An "inverse" of a factory function.  It returns a tuple s.t.

```
args == arguments(f(args...))
```

where `f` is `zip`, `Base.Generator`, `Iterators.filter`, etc.
"""
function arguments end

if isdefined(Iterators, :Zip1)  # VERSION < v"1.1-"
    arguments(xs::Iterators.Zip1) = (xs.a,)
    arguments(xs::Iterators.Zip2) = (xs.a, xs.b)
    arguments(xs::Iterators.Zip) = (xs.a, arguments(xs.z)...)
else
    arguments(xs::Iterators.Zip) = xs.is
end
if isdefined(Iterators, :AbstractZipIterator)  # VERSION < v"1.1-"
    const _Zip = Iterators.AbstractZipIterator
else
    const _Zip = Iterators.Zip
end

arguments(xs::Base.Generator) = (xs.f, xs.iter)
arguments(xs::Iterators.Filter) = (xs.flt, xs.itr)
arguments(xs::Iterators.Flatten) = (xs.it,)
arguments(xs::Iterators.PartitionIterator) = (xs.c, xs.n)
arguments(xs::Base.SkipMissing) = (xs.x,)

safelength(xs) =
    Base.IteratorSize(xs) isa Union{Base.HasLength,Base.HasShape} ? length(xs) : nothing

amount(xs) = length(xs)
amount(xs::AbstractString) = lastindex(xs) - firstindex(xs) + 1
amount(xs::Base.EachStringIndex) = amount(xs.s)

struct WithOffset{Iter,Offset}
    iter::Iter
    offset::Offset
end

withoffset(iter, offset) = WithOffset(iter, offset)
withoffset(iter::WithOffset, offset) = WithOffset(iter.iter, iter.offset + offset)

_shift(wo, ::Nothing) = nothing
_shift(wo, (i, s)) = (convert(eltype(wo), wo.offset + i), s)
Base.iterate(wo::WithOffset) = _shift(wo, iterate(wo.iter))
Base.iterate(wo::WithOffset, s) = _shift(wo, iterate(wo.iter, s))

Base.IteratorEltype(::Type{<:WithOffset{Iter}}) where {Iter} =
    Base.IteratorEltype(Iter)
Base.IteratorSize(::Type{<:WithOffset{Iter}}) where {Iter} =
    Base.IteratorSize(Iter)
Base.eltype(::Type{WithOffset{Iter,Offset}}) where {Iter,Offset} =
    promote_type(eltype(Iter), Offset)
Base.length(wo::WithOffset) = length(wo.iter)
Base.size(wo::WithOffset) = size(wo.iter)
amount(wo::WithOffset) = amount(wo.offset)

function halve(xs::AbstractVector)
    mid = length(xs) ÷ 2
    left = @view xs[firstindex(xs):firstindex(xs)-1+mid]
    right = @view xs[firstindex(xs)+mid:end]
    return (left, right)
end

function halve(xs::AbstractArray)
    i = something(findlast(x -> x > 1, size(xs)), 1)
    mid = size(xs, i) ÷ 2
    leftranges = ntuple(ndims(xs)) do j
        if i == j
            firstindex(xs, j):firstindex(xs, j) + mid - 1
        else
            firstindex(xs, j):lastindex(xs, j)
        end
    end
    rightranges = ntuple(ndims(xs)) do j
        if i == j
            firstindex(xs, j) + mid:lastindex(xs, j)
        else
            firstindex(xs, j):lastindex(xs, j)
        end
    end
    return (view(xs, leftranges...), view(xs, rightranges...))
end

halve(xs::AbstractArray{<:Any,0}) = ((), xs)

# A view wrapper returned from `halve` on `Dict`, `keys`, `values`, and `Set`.
struct DictView{D,F}
    dict::D
    firstslot::Int
    lastslot::Int
    f::F  # identity, first, or last
end

DictView(xs::DictView, i::Int, j::Int) = DictView(xs.dict, i, j, xs.f)

Base.IteratorEltype(::Type{<:DictView{D}}) where {D} = Base.IteratorEltype(D)
Base.IteratorSize(::Type{<:DictView}) = Base.SizeUnknown()

Base.eltype(::Type{<:DictView{D,typeof(identity)}}) where {D} = eltype(D)
Base.eltype(::Type{<:DictView{D,typeof(first)}}) where {D} = keytype(D)
Base.eltype(::Type{<:DictView{D,typeof(last)}}) where {D} = valtype(D)

const DictWrapper = Union{Set,KeySet,ValueIterator}
const DictLike = Union{Dict,DictView,DictWrapper}

DictView(xs::AbstractDict, i::Int, j::Int) = DictView(xs, i, j, identity)
DictView(xs::Union{KeySet,Set}, i::Int, j::Int) = DictView(xs.dict, i, j, first)
DictView(xs::ValueIterator, i::Int, j::Int) = DictView(xs.dict, i, j, last)

# Note: this relies on the implementation detail of `iterate(::Dict)`.
@inline function Base.iterate(xs::DictView, i = xs.firstslot)
    i <= xs.lastslot || return nothing
    y = iterate(xs.dict, i)
    y === nothing && return nothing
    x, j = y
    # If `j` is `xs.lastslot + 1` or smaller, it means the current element is
    # within the range of this `DictView`:
    j <= xs.lastslot + 1 && return xs.f(x), j
    # Otherwise, we need to stop:
    return nothing
end

function Base.length(xs::DictView)
    n = 0
    for _ in xs
        n += 1
    end
    return n
end

firstslot(xs::Dict) = xs.idxfloor
lastslot(xs::Dict) = lastindex(xs.slots)

firstslot(xs::DictView) = xs.firstslot
lastslot(xs::DictView) = xs.lastslot

firstslot(xs::DictWrapper) = firstslot(xs.dict)
lastslot(xs::DictWrapper) = lastslot(xs.dict)

amount(xs::DictLike) = lastslot(xs) - firstslot(xs) + 1

function halve(xs::DictLike)
    i1 = firstslot(xs)
    i3 = lastslot(xs)
    i2 = (i3 - i1 + 1) ÷ 2 + i1
    left = DictView(xs, i1, i2 - 1)
    right = DictView(xs, i2, i3)
    return (left, right)
end

function halve(xs::Iterators.Pairs)
    left, right = halve(keys(xs))
    return Iterators.Pairs(values(xs), left), Iterators.Pairs(values(xs), right)
end

function halve(xs::AbstractString)
    offset = firstindex(xs) - 1
    mid = thisind(xs, (lastindex(xs) - offset) ÷ 2 + offset)
    left = SubString(xs, firstindex(xs):mid)
    right = SubString(xs, nextind(xs, mid):lastindex(xs))
    return (left, right)
end

function halve(xs::Base.EachStringIndex)
    offset = firstindex(xs.s) - 1
    mid = thisind(xs.s, (lastindex(xs.s) - offset) ÷ 2 + offset)
    mid2 = nextind(xs.s, mid)
    left = eachindex(SubString(xs.s, firstindex(xs.s):mid))
    right = withoffset(eachindex(SubString(xs.s, mid2:lastindex(xs.s))), mid2 - 1)
    return (left, right)
end

function halve(xs::WithOffset{<:Base.EachStringIndex})
    left, right = halve(xs.iter)
    return withoffset(left, xs.offset), withoffset(right, xs.offset)
end

@generated function halve(xs::NTuple{N,Any}) where {N}
    m = N ÷ 2
    quote
        (($([:(xs[$i]) for i in 1:m]...),), ($([:(xs[$i]) for i in m+1:N]...),))
    end
end

@inline function halve(xs::NamedTuple{names}) where {names}
    lnames, rnames = halve(names)
    return NamedTuple{lnames}(xs), NamedTuple{rnames}(xs)
end

amount(xs::Base.Generator) = amount(arguments(xs)[2])

function halve(xs::Base.Generator)
    f, coll = arguments(xs)
    left, right = halve(coll)
    return Base.Generator(f, left), Base.Generator(f, right)
end

amount(xs::Iterators.Filter) = amount(arguments(xs)[2])

function halve(xs::Iterators.Filter)
    f, coll = arguments(xs)
    left, right = halve(coll)
    return Iterators.filter(f, left), Iterators.filter(f, right)
end

amount(xs::Iterators.Flatten) = amount(arguments(xs)[1])

function halve(xs::Iterators.Flatten)
    coll, = arguments(xs)
    left, right = halve(coll)
    return Iterators.flatten(left), Iterators.flatten(right)
    # return _flatten(left), _flatten(right)
end

# Is this better?
# _flatten(xs) = safelength(xs) == 1 ? first(xs) : Iterators.flatten(xs)

shape(xs::Union{AbstractArray,Broadcast.Broadcasted}) = size(xs)
shape(xs::Base.Generator) = size(arguments(xs)[2])
shape(xs::Union{NamedTuple,Tuple}) = (length(xs),)

# Make sure that `halve` on the collections give consistent result.
checkshape(xs::_Zip) = length(arguments(xs)) == 1 || shape(xs)
function shape(xs::_Zip)
    args = arguments(xs)
    length(args) == 0 && throw(ArgumentError("empty `zip` does not have `size`."))
    sizes = map(shape, args)
    if !all(==(sizes[1]), tail(sizes))
        throw(ArgumentError("`halve(zip(...))` requires collections with identical `size`"))
    end
    return sizes[1]
end

function halve(xs::_Zip)
    checkshape(xs)
    lefts, rights = _unzip(map(halve, arguments(xs)))
    return zip(lefts...), zip(rights...)
end

shape(xs::Iterators.ProductIterator) = map(length, xs.iterators)

function halve(product::Iterators.ProductIterator)
    i = findlast(x -> length(x) > 1, product.iterators)
    if i === nothing
        # It doesn't matter which "dimension" is used.
        left, right = halve(product.iterators[1])
    else
        left, right = halve(product.iterators[i])
    end
    return (@set(product.iterators[i] = left), @set(product.iterators[i] = right))
end

function shape(xs::Iterators.PartitionIterator)
    coll, = arguments(xs)
    shape(coll)  # make sure it has a shape
    return (length(xs),)
end

function halve(xs::Iterators.PartitionIterator)
    coll, n = arguments(xs)
    m = n * cld(div(length(coll), n), 2)
    offset = firstindex(coll) - 1
    return (
        Iterators.partition(view(coll, offset .+ (1:m)), n),
        Iterators.partition(view(coll, offset .+ (m+1:length(coll))), n),
    )
end

shape(xs::Iterators.Enumerate) = shape(xs.itr)

function halve(xs::Iterators.Enumerate)
    left, right = halve(xs.itr)
    return enumerate(left), zip(reshape(length(left)+1:length(xs), size(right)), right)
end

function halve(xs::Iterators.Reverse)
    left, right = halve(xs.itr)
    return Iterators.reverse(right), Iterators.reverse(left)
end

amount(xs::Base.SkipMissing) = amount(arguments(xs)[1])

function halve(xs::Base.SkipMissing)
    coll, = arguments(xs)
    left, right = halve(coll)
    return skipmissing(left), skipmissing(right)
end


text:
########################## Stochastic Approximations ###########################
# f should be dependent on data x and parameters θ, f(θ, x)
# possibly extend to f(θ, x, y) for regression
# Idea: stochastic and online could be subsumed if we have an infinite vector
# whose getindex method always returns random elements (maybe that is not a good definition)
# computes stochastic approximation to f by randomly samling from x with replacement
# TODO: add option to sample without replacement
function stochastic(f, x::AbstractVector)
    return function sf(θ, batchsize::Int = 1)
        subset = sort!(rand(1:length(x), batchsize)) # without replacement # with replacement: sample(1:length(x), batchsize, replace = false)
        f(θ, x[subset])
    end
end

# stochastic approximation of f(θ, x) by going through x sequentially (in online fashion)
# combine this with PeriodicVector in MyLinearAlgebra to get infinite stream of data
function online(f, x::AbstractVector)
    return function sf(θ, t::Int, batchsize::Int = 1)
        subset = t*batchsize:(t+1)*batchsize
        f(θ, x[subset])
    end
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

@testset "ordinal" begin
    @test ordinal.([3, 32, 391, 2412, 1000000]) ==
        ["3rd", "32nd", "391st", "2412th", "1000000th"]
end

# Test rad2sec
@testset "rad2sec" begin
    @test @inferred(rad2sec(1)) ≈ 206264.80624709636
    @test @inferred(rad2sec(pi)) ≈ 648000.0
    @test @inferred(sec2rad(rad2sec(12.34))) ≈ 12.34
end

# Test sec2rad
@testset "sec2rad" begin
    @test @inferred(sec2rad(3600*30)) ≈ pi/6
    @test @inferred(1/sec2rad(1)) ≈ 206264.80624709636
    @test @inferred(rad2sec(sec2rad(56.78))) ≈ 56.78
end


text:
"""
Authors : Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020
"""

# Ne pas enlever
include("./main.jl")

println("Chargement de JuMP")
using JuMP
println("Chargement de CPLEX")
using CPLEX
println("Cplex chargé")

@doc """
Modifie la solution sol de l'instance inst en appliquant CPLEX sur un
programme linéaire modélisant le problème (G_ULS).
""" ->
function run(inst, sol)
    T = inst.T
    M = inst.M
    R = inst.R
    Emax = inst.Emax
    f = inst.f
    e = inst.e
    h = inst.h
    p = inst.p
    d = inst.d

	m = Model(with_optimizer(CPLEX.Optimizer))
	MOI.set(m, MOI.RawParameter("CPXPARAM_MIP_Display"),0)

	@variable(m, 0 <= x[1:T,1:M])
	@variable(m, 0 <= s[1:T])
	@variable(m, y[1:T,1:M], Bin)

	## OBJECTIF
	##
	@objective(m, Min, sum(sum(p[t][j]*x[t,j]+f[t][j]*y[t,j] for j in 1:M)+h[t]*s[t] for t in 1:T))

	## CONTRA1NTE 1
	##
	## Contrainte de conservation
	@constraint(m, constraint1[t in 2:T], sum( x[t,j] for j in 1:M)-s[t]+s[t-1] == d[t])

	## CONTRA1NTE 2
	##
	## Contrainte de conservation
	@constraint(m, constraint2, sum(x[1,j]-s[1] for j in 1:M) == d[1])

	## CONTRA1NTE 3
	##
	## Contrainte de demande
	@constraint(m, constraint3[t in 1:T, j in 1:M], sum(d[t1] for t1 in t:T)*y[t,j]>=x[t,j])

	## CONTRA1NTE 4
	##
	## Contrainte écologique
	@constraint(m, constraint4[t in R:T], sum(sum( (e[t1][j]-Emax[t1])*x[t1,j] for j in 1:M) for t1 in t-R+1:t) <=0)

	optimize!(m)

	return (m,x,s,y)
end

@doc """
Cette fonction est appelée après la fonction `run` et permet de faire de l'affichage et des traitement sur la sortie de la fonction `run` ; sans pour autant affecter son temps de calcul.

Le paramètre cpu time est le temps de calcul de `run`. Les valeurs de `inst` et `sol` sont les mêmes qu’à la sortie de la fonction run. Enfin, `others` est ce qui est renvoyé par la fonction `run`. Vous pouvez ainsi effectuer des tests et afficher des résultats sans affecter le temps de calcul.
""" ->
function post_process(cpu_time::Float64, inst, sol, others)

	# Run a renvoyé le modèle et ses variables, qui ont été mis dans others.
	m,x,s,y = others

	solx = value.(x)
	emission_moyenne = 0
	prod_moyenne = 0
	for i in 1:inst.T
	    for j in 1:inst.M
	    	emission_moyenne+=solx[i,j]*inst.e[i][j]
	    	prod_moyenne+=solx[i,j]
	    end
	end
	emission_moyenne = emission_moyenne/prod_moyenne

	println("OBJECTIF : $(objective_value(m))")
	println("emission moyenne : $(emission_moyenne)")
	# println("VALEURS de x : $(value.(x))")
	# println("VALEURS de s : $(value.(s))")
	# println("VALEURS de y : $(value.(y))")

	# println("Temps de calcul : $cpu_time.")
end

# Ne pas enlever
if length(ARGS) > 0
	input_file = ARGS[1]

	main(input_file)
end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
# ------------------------------------------------------------------------------------------
# Models
# ------------------------------------------------------------------------------------------
struct PatMat{S<:Surrogate, T<:Real} <: AbstractPatMat{S}
    τ::T
    C::T
    l1::S
    l2::S
    state::State{T}

    function PatMat(
        τ::Real;
        C::Real = 1,
        ϑ1::Real = 1,
        ϑ2::Real = 1,
        surrogate::Type{<:Surrogate} = Hinge,
        T = Float32,
    )
        l1 = surrogate(T(ϑ1))
        l2 = surrogate(T(ϑ2))
        return new{typeof(l1), T}(T(τ), T(C), l1, l2, State(T))
    end
end

struct PatMatNP{S<:Surrogate, T<:Real} <: AbstractPatMat{S}
    τ::T
    C::T
    l1::S
    l2::S
    state::State{T}

    function PatMatNP(
        τ::Real;
        C::Real = 1,
        ϑ1::Real = 1,
        ϑ2::Real = 1,
        surrogate::Type{<:Surrogate} = Hinge,
        T = Float32,
    )
        l1 = surrogate(T(ϑ1))
        l2 = surrogate(T(ϑ2))
        return new{typeof(l1), T}(T(τ), T(C), l1, l2, State(T))
    end
end

function parameters(m::AbstractPatMat{S}) where {S<:Surrogate}
    return (τ = m.τ, C = m.C, surrogate = S.name.name, ϑ1 = m.l1.ϑ, ϑ2 = m.l2.ϑ)
end

function threshold(model::AbstractPatMat, K::KernelMatrix)
    s_neg = .- model.state.s[inds_β(K)]
    f(t) = sum(value.(model.l2, s_neg .- t)) - K.nβ*model.τ

    return find_root(f, (-Inf, Inf))
end

# ------------------------------------------------------------------------------------------
# Auxilliary functions
# ------------------------------------------------------------------------------------------
function initialization!(model::AbstractPatMat{S}, K::KernelMatrix) where {S <: Surrogate}
    T = eltype(K)
    α, β = rand(T, K.nα), rand(T, K.nβ)
    if S <: Hinge
        δ = maximum(β)/model.l1.ϑ
        α, β, δ = projection(α, β, δ, model.l1.ϑ*model.C, model.l2.ϑ)
    elseif S <: Quadratic
        δ = sqrt(max(sum(abs2, β)/(4*K.nβ*model.τ*model.l2.ϑ^2), 0))
        α, β, δ = projection(α, β, δ)
    else
        throw(ArgumentError("$(S.name.name) surrogate function is not supported"))
    end
    αβ = vcat(α, β)

    model.state.s = K * αβ
    model.state.αβ = αβ
    model.state.δ = δ
    model.state.βsort = sort(αβ[inds_β(K)], rev = true)
    return
end

function permutation(::PatMat, y::BitVector)
    perm_α = findall(y)
    perm_β = 1:length(y)

    return length(perm_α), length(perm_β), vcat(perm_α, perm_β)
end

function permutation(::PatMatNP, y::BitVector)
    perm_α = findall(y)
    perm_β = findall(.~y)

    return length(perm_α), length(perm_β), vcat(perm_α, perm_β)
end

function extract_scores(model::PatMat, K::KernelMatrix)
    s = copy(model.state.s)
    return .- s[inds_β(K)]
end

function extract_scores(model::PatMatNP, K::KernelMatrix)
    s = copy(model.state.s)
    s[inds_β(K)] .*= -1
    return s[invperm(K.perm)]
end

function extract_state(model::AbstractPatMat, K::KernelMatrix)
    αβ = copy(model.state.αβ)
    return  (α = αβ[inds_α(K)], β = αβ[inds_β(K)], δ = model.state.δ)
end

function objective(model::AbstractPatMat{S}, K::KernelMatrix) where {S<:Surrogate}
    @unpack s, αβ, δ = model.state
    τ = model.τ
    C, ϑ1, ϑ2 = model.C, model.l1.ϑ, model.l2.ϑ

    α = αβ[inds_α(K)]
    β = αβ[inds_β(K)]
    n = K.nβ
    sα = s[inds_α(K)]
    w_norm = s'*αβ/2
    t = threshold(model, K)

    # objectives
    L_primal = w_norm + C*sum(value.(model.l1, t .- sα))
    if S <: Hinge
        L_dual = - w_norm + sum(α)/ϑ1 + sum(β)/ϑ2 - δ*n*τ
    else
        L_dual = - w_norm + sum(α)/ϑ1 - sum(abs2, α)/(4*C*ϑ1^2) + sum(β)/ϑ2 - sum(abs2, β)/(4*δ*ϑ2^2) - δ*n*τ
    end
    return L_primal, L_dual, L_primal - L_dual
end

# ------------------------------------------------------------------------------------------
# Update rules - Hinge loss
# ------------------------------------------------------------------------------------------
function rule_αα(model::AbstractPatMat{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, δ = model.state
    C, ϑ1 = model.C, model.l1.ϑ

    num = s[k] - s[l]
    den = K[k, k] - 2*K[k, l] + K[l, l]
    lb = max(- αβ[k], αβ[l] - ϑ1*C)
    ub = min(ϑ1*C - αβ[k], αβ[l])

    return UpdateRule(model, K, num, den, lb, ub, k, l, δ)
end

function rule_αβ(model::AbstractPatMat{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ = model.state
    C, ϑ1, ϑ2 = model.C, model.l1.ϑ, model.l2.ϑ
    βmax = find_βmax(model, αβ[l])

    num = s[k] + s[l] - 1/ϑ1 - 1/ϑ2
    den = K[k, k] + 2*K[k, l] + K[l, l]
    lb = max(- αβ[k], - αβ[l])
    ub = ϑ1*C - αβ[k]

    # solution 1
    Δ1 = compute_Δ(num, den, lb, ub)
    r1 = if αβ[l] + Δ1 <= βmax
        UpdateRule(model, K, num, den, lb, ub, k, l, βmax/ϑ2)
    else
        UpdateRule(k, l; T = eltype(K))
    end

    # solution 2
    c = K.nβ*model.τ/ϑ2
    Δ2 = compute_Δ(num + c, den, lb, ub)
    r2 = if αβ[l] + Δ2 >= βmax
        UpdateRule(model, K, num + c, den, lb, ub, k, l, (αβ[l] + Δ2)/ϑ2)
    else
        UpdateRule(k, l; T = eltype(K))
    end
    return select(r1, r2)
end

function rule_ββ(model::AbstractPatMat{<:Hinge}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ = model.state
    ϑ2 = model.l2.ϑ
    βmax = find_βmax(model, αβ[l])

    num = s[k] + s[l]
    den = K[k, k] - 2*K[k, l] + K[l, l]
    lb = - αβ[k]
    ub = αβ[l]

    # solution 1
    Δ1 = compute_Δ(num, den, lb, ub)
    r1 = if max(αβ[k] + Δ1, αβ[l] - Δ1) <= βmax
        UpdateRule(model, K, num, den, lb, ub, k, l, βmax/ϑ2)
    else
        UpdateRule(k, l; T = eltype(K))
    end

    # solution 2
    c = K.nβ*model.τ/ϑ2
    Δ2 = compute_Δ(num + c, den, lb, ub)
    r2 = if αβ[k] + Δ2 >= max(βmax, αβ[l] - Δ2)
        UpdateRule(model, K, num + c, den, lb, ub, k, l, (αβ[k] + Δ2)/ϑ2)
    else
        UpdateRule(k, l; T = eltype(K))
    end

    # solution 3
    Δ3 = compute_Δ(num - c, den, lb, ub)
    r3 = if αβ[l] - Δ3 >= max(βmax, αβ[k] + Δ2)
        UpdateRule(model, K, num - c, den, lb, ub, k, l, (αβ[l] - Δ3)/ϑ2)
    else
        UpdateRule(k, l; T = eltype(K))
    end
    return select(select(r1, r2), r3)
end

# ------------------------------------------------------------------------------------------
# Update rules - Quadratic loss
# ------------------------------------------------------------------------------------------
function rule_αα(model::AbstractPatMat{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, δ = model.state
    C, ϑ1 = model.C, model.l1.ϑ

    num = s[k] - s[l] + (αβ[k] - αβ[l])/(2*C*ϑ1^2)
    den = K[k, k] - 2*K[k, l] + K[l, l] + 1/(C*ϑ1^2)
    lb = - αβ[k]
    ub = αβ[l]

    return UpdateRule(model, K, num, den, lb, ub, k, l, δ)
end

function rule_αβ(model::AbstractPatMat{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, δ = model.state
    C, ϑ1, ϑ2 = model.C, model.l1.ϑ, model.l2.ϑ
    n, τ = K.nβ, model.τ

    num = s[k] + s[l] - 1/ϑ1 + αβ[k]/(2*C*ϑ1^2) - 1/ϑ2 + αβ[l]/(2*δ*ϑ2^2)
    den = K[k, k] + 2*K[k, l] + K[l, l] + 1/(2*C*ϑ2^2) + 1/(2*δ*ϑ2^2)
    lb = max(- αβ[k], - αβ[l])
    ub = eltype(s)(Inf)

    Δ = compute_Δ(num, den, lb, ub)
    δnew = sqrt(max(δ^2 + (Δ^2 + 2*Δ*αβ[l])/(4*ϑ2^2*n*τ), 0))
    # num -= (1/δ - 1/δnew)*αβ[l]/(2*ϑ2^2)
    # den -= (1/δ - 1/δnew)/(2*ϑ2^2)

    return UpdateRule(model, K, num, den, lb, ub, k, l, δnew)
end

function rule_ββ(model::AbstractPatMat{<:Quadratic}, K::KernelMatrix, k::Int, l::Int)
    @unpack s, αβ, δ = model.state
    ϑ2 = model.l2.ϑ
    n, τ = K.nβ, model.τ

    num = s[k] - s[l] + (αβ[k] - αβ[l])/(2*δ*ϑ2^2)
    den = K[k, k] - 2*K[k, l] + K[l, l] + 1/(δ*ϑ2^2)
    lb = - αβ[k]
    ub = αβ[l]

    Δ = compute_Δ(num, den, lb, ub)
    δnew = sqrt(max(δ^2 + (Δ^2 + Δ*(αβ[k] - αβ[l]))/(2*ϑ2^2*n*τ), 0))
    # num -= (1/δ - 1/δnew)*(αβ[k] - αβ[l])/(2*ϑ2^2)
    # den -= (1/δ - 1/δnew)/(ϑ2^2)

    return UpdateRule(model, K, num, den, lb, ub, k, l, δnew)
end


text:
# Defines the function that computes the transit:
include("../../../src/transit_poly_struct.jl")
using PyPlot
using PyCall

# System properties
npts = 10000
nu = 5
u_n = 0.5*ones(nu) / nu
r = 0.1
b = zeros(npts)
for i=1:npts
    b[i] = 3.0 * ((i-0.5) / float(npts) - 0.5)
end

# Julia flux
flux_julia = zeros(npts)
#flux_big = zeros(npts)
trans = transit_init(r, 0.0, u_n, true)
dfdx = zeros(nu+2,npts)
#trans_big = transit_init(big(r), big(0.0), big.(u_n), false)
for i=1:npts
    trans.b = abs(b[i])
    flux_julia[i] = transit_poly!(trans)
    dfdx[1:2,i] = trans.dfdrb
    dfdx[3:nu+2,i] = trans.dfdu
    if b[i] < 0.0
      dfdx[2,i] *= -1.0
    end
#    trans_big.b = big(abs(b[i]))
#    flux_big[i] = convert(Float64,transit_poly!(trans_big))
#    flux_julia[i] = transit_poly(r,abs(b[i]),u_n)
end

# Starry flux
@pyimport starry
map = starry.Map(nu)
for i=1:nu
    map[:__setitem__](i, u_n[i])
end
flux_starry, grad = map[:flux](xo=b, ro=r, gradient = true)

# Plot
fig,axes = subplots(2,2)
ax = axes[1]
ax[:plot](b, flux_julia, linewidth=2, label="julia")
ax[:plot](b, flux_starry, linewidth=1, label="starry")
ax[:set_xlabel]("Impact parameter")
ax[:set_ylabel]("Flux")
ax[:legend]()
ax[:axis]([-1.5,1.5,0.985, 1.01])

ax = axes[2]
#ax[:semilogy](b, abs.(flux_julia-flux_big), linewidth=2, label="julia-big")
ax[:semilogy](b, abs.(flux_julia-flux_starry), linewidth=2, label="julia-starry")
#ax[:plot](b, abs.(flux_starry-flux_big), linewidth=2, label="starry-big")
ax[:set_xlabel]("abs(Difference)")
ax[:legend]()
#ax[:axis]([-1.5,1.5,1e-10,1e-5])

ax = axes[3]
ax[:plot](b, grad["ro"],label=L"$df/dr$ (starry)")
ax[:plot](b, dfdx[1,:],label=L"$df/dr$ (limbdark)")
ax[:plot](b, grad["xo"],label=L"$df/db$ (starry)")
ax[:plot](b, dfdx[2,:],label=L"$df/db$ (limbdark)")
for iu=1:nu
  ax[:plot](b, grad["u"][iu,:],label=L"$df/du$ (starry)")
  ax[:plot](b, dfdx[2+iu,:],label=L"$df/du$ (limbdark)")
end
ax[:set_xlabel]("Derivative")
ax[:set_ylabel]("Flux")
#ax[:legend](fontsize=6,loc="lower right")

ax = axes[4]
ax[:semilogy](b, abs.(grad["ro"]-dfdx[1,:]),label=L"$df/dr$ diff")
ax[:plot](b, abs.(grad["xo"]-dfdx[2,:]),label=L"$df/db$ diff")
for iu=1:nu
  ax[:plot](b, abs.(grad["u"][iu,:]-dfdx[2+iu,:]),label=L"$df/du$ diff")
end
ax[:set_xlabel]("Derivative differences")
ax[:set_ylabel]("Flux")
ax[:legend](fontsize=4,loc="lower right")
#ax[:axis]([-1.5,1.5,1e-10,1e-5])

savefig("compare_to_starry_" * string(nu) * ".pdf", bbox_inches="tight")


text:
@testset "Parallel.BFS" begin

    g5 = SimpleDiGraph(4)
    add_edge!(g5, 1, 2); add_edge!(g5, 2, 3); add_edge!(g5, 1, 3); add_edge!(g5, 3, 4)

    for g in testdigraphs(g5)
        @test @inferred(Parallel.gdistances(g, 1))  == gdistances(g, 1)
        @test @inferred(Parallel.gdistances(g, [1, 3]))  == gdistances(g, [1, 3])
    end

    g6 = smallgraph(:house)

    for g in testgraphs(g6)
        @test @inferred(Parallel.gdistances(g, 2))  == gdistances(g, 2)
        @test @inferred(Parallel.gdistances(g, [1, 2])) == gdistances(g, [1, 2])
    end

end


text:
function DiffEqBase.__solve(
    prob::DiffEqBase.AbstractODEProblem{uType,tuptType,isinplace},
    alg::AlgType,
    timeseries=[],ts=[],ks=[];
    saveat = Float64[],
    verbose=true,save_everystep = isempty(saveat),
    save_on = true,
    save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[1] in saveat,
    timeseries_errors=true,dense_errors=false,
    callback=nothing, alias_u0=false, kwargs...) where
    {uType,tuptType,isinplace,AlgType<:ODEInterfaceAlgorithm}

    tType = eltype(tuptType)

    isstiff = alg isa ODEInterfaceImplicitAlgorithm
    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction) && isstiff
            if DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
    if max_len_cb isa VectorContinuousCallback
      callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltype,uBottomEltype)
    else
      callback_cache = nothing
    end

    tspan = prob.tspan

    o = KW(kwargs)

    u0 = prob.u0

    if typeof(u0) <: Number
        u = [u0]
    else
        if alias_u0
            u = u0
        else
            u = deepcopy(u0)
        end
    end

    tdir = sign(tspan[2]-tspan[1])

    saveat_internal =
      saveat_disc_handling(saveat,tdir,tspan,tType)

    sizeu = size(u)

    o[:RHS_CALLMODE] = ODEInterface.RHS_CALL_INSITU

    if save_everystep
        _timeseries = Vector{uType}(undef,0)
        ts = Vector{tType}(undef,0)
    else
        _timeseries = [copy(u0)]
        ts = [tspan[1]]
    end

    uprev = similar(u)

    sol = DiffEqBase.build_solution(prob,  alg, ts, _timeseries,
                         timeseries_errors = timeseries_errors,
                         calculate_error = false,
                         destats = DiffEqBase.DEStats(0),
                         retcode = :Default)

    opts = DEOptions(saveat_internal,save_on,save_everystep,callbacks_internal)
    if !isinplace && typeof(u)<:AbstractArray
        f! = (t,u,du) -> (du[:] = vec(prob.f(reshape(u,sizeu),integrator.p,t)); nothing)
    elseif !(typeof(u)<:Vector{Float64})
        f! = (t,u,du) -> (prob.f(reshape(du,sizeu),reshape(u,sizeu),integrator.p,t);
                          du=vec(du); nothing)
    else
        f! = (t,u,du) -> prob.f(du,u,integrator.p,t)
    end

    integrator = ODEInterfaceIntegrator(prob.f,u,uprev,tspan[1],tspan[1],prob.p,opts,
                                        false,tdir,sizeu,sol,
                                        (t)->[t],0,1,callback_cache,alg,0.)
    initialize_callbacks!(integrator)

    outputfcn = OutputFunction(integrator)
    o[:OUTPUTFCN] = outputfcn
    if !(typeof(callbacks_internal.continuous_callbacks)<:Tuple{}) || !isempty(saveat)
        if typeof(alg) <: Union{ddeabm,ddebdf}
            @warn("saveat and continuous callbacks ignored for ddeabm and ddebdf")
            o[:OUTPUTMODE] = ODEInterface.OUTPUTFCN_WODENSE
        else
            o[:OUTPUTMODE] = ODEInterface.OUTPUTFCN_DENSE
        end
    else
        o[:OUTPUTMODE] = ODEInterface.OUTPUTFCN_WODENSE
    end

    dict = buildOptions(o,
                        ODEINTERFACE_OPTION_LIST,
                        ODEINTERFACE_ALIASES,
                        ODEINTERFACE_ALIASES_REVERSED)
    if prob.f.mass_matrix != I
        if typeof(prob.f.mass_matrix) <: Matrix && isstiff
            dict[:MASSMATRIX] = prob.f.mass_matrix
        elseif !isstiff
            error("This solver does not support mass matrices")
        else
            error("This solver must use full or banded mass matrices.")
        end
    end
    if DiffEqBase.has_jac(prob.f)
        dict[:JACOBIMATRIX] = (t,u,J) -> prob.f.jac(J,u,prob.p,t)
    end

    if isstiff && alg.jac_lower !== nothing
        dict[:JACOBIBANDSSTRUCT] = (alg.jac_lower,alg.jac_upper)
    end

    # Convert to the strings
    opts = ODEInterface.OptionsODE([Pair(ODEINTERFACE_STRINGS[k],v) for (k,v) in dict]...)

    if typeof(alg) <: dopri5
        tend, uend, retcode, stats =
            ODEInterface.dopri5(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: dop853
        tend, uend, retcode, stats =
            ODEInterface.dop853(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: odex
        tend, uend, retcode, stats =
            ODEInterface.odex(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: seulex
        tend, uend, retcode, stats =
            ODEInterface.seulex(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: radau
        tend, uend, retcode, stats =
            ODEInterface.radau(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: radau5
        tend, uend, retcode, stats =
            ODEInterface.radau5(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: rodas
        tend, uend, retcode, stats =
            ODEInterface.rodas(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: ddeabm
        tend, uend, retcode, stats =
            ODEInterface.ddeabm(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    elseif typeof(alg) <: ddebdf
        tend, uend, retcode, stats =
            ODEInterface.ddebdf(f!, tspan[1], tspan[2], vec(integrator.u), opts)
    end

    if !save_everystep
        push!(ts,tend)
        save_value!(_timeseries,uend,uType,sizeu)
    end

    if retcode < 0
        if retcode == -1
            verbose && @warn("Input is not consistent.")
            return_retcode = :Failure
        elseif retcode == -2
            verbose && @warn("Interrupted. Larger maxiters is needed.")
            return_retcode = :MaxIters
        elseif retcode == -3
            verbose && @warn("Step size went too small.")
            return_retcode = :DtLessThanMin
        elseif retcode == -4
            verbose && @warn("Interrupted. Problem is probably stiff.")
            return_retcode = :Unstable
        end
    else
        return_retcode = :Success
    end

    if DiffEqBase.has_analytic(prob.f)
        DiffEqBase.calculate_solution_errors!(integrator.sol;
        timeseries_errors=timeseries_errors,
        dense_errors=dense_errors)
    end

    destats = sol.destats
    destats.nf = stats["no_rhs_calls"]
    if haskey(stats,"no_steps_rejected")
        destats.nreject = stats["no_steps_rejected"]
        destats.naccept = stats["no_steps_accepted"]
    end
    if haskey(stats,"no_jac_calls")
        destats.njacs = stats["no_jac_calls"]
    end
    if haskey(stats,"no_lu_decomp")
        destats.nw = stats["no_lu_decomp"]
    end
    DiffEqBase.solution_new_retcode(sol,return_retcode)
end

function save_value!(_timeseries,u,::Type{T},sizeu) where T<:Number
    push!(_timeseries,first(u))
end

function save_value!(_timeseries,u,::Type{T},sizeu) where T<:Vector
    push!(_timeseries,u)
end

function save_value!(_timeseries,u,::Type{T},sizeu) where T<:Array
    push!(_timeseries,reshape(u,sizeu))
end

function buildOptions(o, optionlist, aliases, aliases_reversed)
    dict1 = Dict{Symbol,Any}([Pair(k,o[k]) for k in (keys(o) ∩ optionlist)])
    dict2 = Dict([Pair(aliases_reversed[k],o[k]) for k in (keys(o) ∩ values(aliases))])
    merge(dict1,dict2)
end

function saveat_disc_handling(saveat,tdir,tspan,tType)

  if typeof(saveat) <: Number
    if (tspan[1]:saveat:tspan[end])[end] == tspan[end]
      saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:tspan[end]))
    else
      saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat)))
    end
  else
    saveat_vec = vec(collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x<tdir*tspan[end],saveat)))
  end

  if tdir>0
    saveat_internal = BinaryMinHeap(saveat_vec)
  else
    saveat_internal = BinaryMaxHeap(saveat_vec)
  end

  saveat_internal
end

const ODEINTERFACE_OPTION_LIST =
    Set([:RTOL, :ATOL, :OUTPUTFCN, :OUTPUTMODE, :MAXSTEPS, :STEST, :EPS, :RHO, :SSMINSEL,
         :SSMAXSEL, :SSBETA, :MAXSS, :INITIALSS, :MAXEXCOLUMN, :STEPSIZESEQUENCE,
         :MAXSTABCHECKS, :MAXSTABCHECKLINE, :DENSEOUTPUTWOEE, :INTERPOLDEGRE,
         :SSREDUCTION, :SSSELECTPAR1, :SSSELECTPAR2, :ORDERDECFRAC, :ORDERINCFRAC,
         :OPT_RHO, :OPT_RHO2, :RHSAUTONOMOUS, :M1, :M2, :LAMBDADENSE, :TRANSJTOH,
         :STEPSIZESEQUENCE, :JACRECOMPFACTOR, :MASSMATRIX, :JACOBIMATRIX, :JACOBIBANDSSTRUCT,
         :WORKFORRHS, :WORKFORJAC, :WORKFORDEC, :WORKFORSOL,
         :MAXNEWTONITER, :NEWTONSTARTZERO, :NEWTONSTOPCRIT, :DIMFIND1VAR,
         :MAXSTAGES, :MINSTAGES, :INITSTAGES, :STEPSIZESTRATEGY,
         :FREEZESSLEFT, :FREEZESSRIGHT, :ORDERDECFACTOR,
         :ORDERINCFACTOR, :ORDERDECCSTEPFAC1, :ORDERDECSTEPFAC2, :RHS_CALLMODE
         ])

const ODEINTERFACE_ALIASES =
    Dict{Symbol,Symbol}(:RTOL=>:reltol,
                        :ATOL=>:abstol,
                        :MAXSTEPS=> :maxiters,
                        :MAXSS=>:dtmax,
                        :INITIALSS=>:dt,
                        #:SSMINSEL=>:qmin,
                        :SSBETA=>:beta2,
                        :SSMAXSEL=>:qmax)

const ODEINTERFACE_ALIASES_REVERSED =
    Dict{Symbol,Symbol}([(v,k) for (k,v) in ODEINTERFACE_ALIASES])

const ODEINTERFACE_STRINGS = Dict{Symbol,String}(
  :LOGIO            => "logio",
  :LOGLEVEL         => "loglevel",
  :RHS_CALLMODE     => "RightHandSideCallMode",

  :RTOL             => "RelTol",
  :ATOL             => "AbsTol",
  :MAXSTEPS         => "MaxNumberOfSteps",
  :EPS              => "eps",

  :OUTPUTFCN        => "OutputFcn",
  :OUTPUTMODE       => "OutputFcnMode",

  :STEST            => "StiffTestAfterStep",
  :RHO              => "rho",
  :SSMINSEL         => "StepSizeMinSelection",
  :SSMAXSEL         => "StepSizeMaxSelection",
  :SSBETA           => "StepSizeBeta",
  :MAXSS            => "MaxStep",
  :INITIALSS        => "InitialStep",

  :MAXEXCOLUMN      => "MaxExtrapolationColumn",
  :MAXSTABCHECKS    => "MaxNumberOfStabilityChecks",
  :MAXSTABCHECKLINE => "MaxLineForStabilityCheck",
  :INTERPOLDEGREE   => "DegreeOfInterpolation",
  :ORDERDECFRAC     => "OrderDecreaseFraction",
  :ORDERINCFRAC     => "OrderIncreaseFraction",
  :STEPSIZESEQUENCE => "StepSizeSequence",
  :SSREDUCTION      => "StepSizeReduction",
  :SSSELECTPAR1     => "StepSizeSelectionParam1",
  :SSSELECTPAR2     => "StepSizeSelectionParam2",
  :RHO2             => "rho2",
  :DENSEOUTPUTWOEE  => "DeactivateErrorEstInDenseOutput",

  :TRANSJTOH        => "TransfromJACtoHess",
  :MAXNEWTONITER    => "MaxNewtonIterations",
  :NEWTONSTARTZERO  => "StartNewtonWithZeros",
  :DIMOFIND1VAR     => "DimensionOfIndex1Vars",
  :DIMOFIND2VAR     => "DimensionOfIndex2Vars",
  :DIMOFIND3VAR     => "DimensionOfIndex3Vars",
  :STEPSIZESTRATEGY => "StepSizeStrategy",
  :M1               => "M1",
  :M2               => "M2",
  :JACRECOMPFACTOR  => "RecomputeJACFactor",
  :NEWTONSTOPCRIT   => "NewtonStopCriterion",
  :FREEZESSLEFT     => "FreezeStepSizeLeftBound",
  :FREEZESSRIGHT    => "FreezeStepSizeRightBound",
  :MASSMATRIX       => "MassMatrix",
  :JACOBIMATRIX     => "JacobiMatrix",
  :JACOBIBANDSTRUCT => "JacobiBandStructure",

  :MAXSTAGES        => "MaximalNumberOfStages",
  :MINSTAGES        => "MinimalNumberOfStages",
  :INITSTAGES       => "InitialNumberOfStages",
  :ORDERINCFACTOR   => "OrderIncreaseFactor",
  :ORDERDECFACTOR   => "OrderDecreaseFactor",
  :ORDERDECSTEPFAC1 => "OrderDecreaseStepFactor1",
  :ORDERDECSTEPFAC2 => "OrderDecreaseStepFactor2",

  :RHSAUTONOMOUS    => "AutonomousRHS",
  :LAMBDADENSE      => "LambdaForDenseOutput",
  :WORKFORRHS       => "WorkForRightHandSide",
  :WORKFORJAC       => "WorkForJacobimatrix",
  :WORKFORDEC       => "WorkForLuDecomposition",
  :WORKFORSOL       => "WorkForSubstitution",

  :BVPCLASS         => "BoundaryValueProblemClass",
  :SOLMETHOD        => "SolutionMethod",
  :IVPOPT           => "OptionsForIVPsolver")

struct OutputFunction{T} <: Function
    integrator::T
end

function (f::OutputFunction)(reason::ODEInterface.OUTPUTFCN_CALL_REASON,
      tprev::Float64, t::Float64, u::Vector{Float64},
      eval_sol_fcn, extra_data::Dict)

  if reason == ODEInterface.OUTPUTFCN_CALL_STEP

      integrator = f.integrator

      integrator.uprev .= integrator.u
      if eltype(integrator.sol.u) <: Vector
          integrator.u .= u
      else
          integrator.u .= reshape(u,integrator.sizeu)
      end
      integrator.t = t
      integrator.tprev = tprev
      integrator.eval_sol_fcn = eval_sol_fcn

      handle_callbacks!(integrator,eval_sol_fcn)

      if integrator.u_modified

          if eltype(integrator.sol.u) <: Vector
              u .= integrator.u
          else
              tmp = reshape(u,integrator.sizeu)
              tmp .= integrator.u
          end

          return ODEInterface.OUTPUTFCN_RET_CONTINUE_XCHANGED
      else
          return ODEInterface.OUTPUTFCN_RET_CONTINUE
      end
      # TODO: ODEInterface.OUTPUTFCN_RET_STOP for terminate!

  end

  ODEInterface.OUTPUTFCN_RET_CONTINUE
end


text:
# src path
src_path = "C:\\Users\\Koike\\Desktop\\git\\test\\test_2d_general_NS\\src\\"

# main (変更しないこと)
src_read="read_grid.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="output.jl"
include(src_path*src_read)

src_read="advection_term.jl"
include(src_path*src_read)
src_read="boundary.jl"
include(src_path*src_read)
src_read="converge.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_read="rhs.jl"
include(src_path*src_read)
src_read="value_setup.jl"
include(src_path*src_read)
src_read="viscos.jl"
include(src_path*src_read)

src_read="setup.jl"
include(src_path*src_read)
src_read="cal_time_step.jl"
include(src_path*src_read)
src_main="main.jl"
include(src_path*src_main)

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
basis_matrix(L::ZLat) = L.basis_matrix

ambient_space(L::ZLat) = L.space

base_ring(L::ZLat) = FlintZZ

################################################################################
#
#  Creation
#
################################################################################

function Zlattice(B::fmpq_mat; gram = identity_matrix(FlintQQ, ncols(B)))
  V = quadratic_space(FlintQQ, gram)
  return lattice(V, B)
end

function Zlattice(B::fmpz_mat; gram = identity_matrix(FlintQQ, ncols(B)))
  V = quadratic_space(FlintQQ, gram)
  return lattice(V, B)
end

function Zlattice(;gram)
  n = nrows(gram)
  return lattice(quadratic_space(FlintQQ, gram), identity_matrix(FlintQQ, n))
end

function lattice(V::QuadSpace{FlintRationalField, fmpq_mat}, B::MatElem)
  Gc = change_base_ring(FlintQQ, B)
  if typeof(Gc) !== fmpq_mat
    throw(error("Cannot convert entries of the matrix to the rationals"))
  end
  return ZLat(V, Gc)
end

@doc Markdown.doc"""
    orthogonal_sum(L1::ZLat, L2::ZLat)

Return the orthogonal direct sum of the lattice L1 and L2.
"""
function orthogonal_sum(L1::ZLat, L2::ZLat)
  V1 = ambient_space(L1)
  V2 = ambient_space(L2)
  B1 = basis_matrix(L1)
  B2 = basis_matrix(L2)
  B = diagonal_matrix(B1, B2)
  V = orthogonal_sum(V1, V2)
  return lattice(V, B)
end

function gram_matrix(L::ZLat)
  b = basis_matrix(L)
  return b * gram_matrix(ambient_space(L)) * transpose(b)
end

gram_matrix_of_rational_span(L::ZLat) = gram_matrix(L)

function rational_span(L::ZLat)
  if isdefined(L, :rational_span)
    return L.rational_span
  else
    G = gram_matrix(L)
    V = quadratic_space(FlintQQ, G)
    L.rational_span = V
    return V
  end
end

################################################################################
#
#  String I/O
#
################################################################################

# if ambient_representation = true, they are given with respect to the ambient space
function Base.show(io::IO, L::ZLat)
  print(io, "Quadratic lattice of rank ", rank(L),
            " and degree ", degree(L), " over the rationals")
end

function assert_has_automorphisms(L::ZLat; redo::Bool = false,
                                           try_small::Bool = true)

  if !redo && isdefined(L, :automorphism_group_generators)
    return nothing
  end

  V = ambient_space(L)
  GL = gram_matrix(L)
  d = denominator(GL)
  res = fmpz_mat[change_base_ring(FlintZZ, d * GL)]
  # So the first one is either positive definite or negative definite
  # Make it positive definite. This does not change the automorphisms.
  if res[1][1, 1] < 0
    res[1] = -res[1]
  end
  Glll, T = lll_gram_with_transform(res[1])
  Ttr = transpose(T)
  res_orig = copy(res)
  res[1] = Glll

  bm = basis_matrix(L)

  # Make the Gram matrix small

  C = ZLatAutoCtx(res)
  if try_small
    fl, Csmall = try_init_small(C)
    if fl
      auto(Csmall)
      _gens, order = _get_generators(Csmall)
      gens = fmpz_mat[matrix(ZZ, g) for g in _gens]
    else
      init(C)
      auto(C)
      gens, order = _get_generators(C)
    end
  else
    init(C)
    auto(C)
    gens, order = _get_generators(C)
  end

  # Now translate back
  Tinv = inv(T)
  for i in 1:length(gens)
    gens[i] = Tinv * gens[i] * T
  end

  # Now gens are with respect to the basis of L
  @hassert all(gens[i] * GL * transpose(gens[i]) == GL for i in 1:length(gens))

  L.automorphism_group_generators = gens
  L.automorphism_group_order = order

  return nothing
end

# natural action = action on ambient_space

@doc Markdown.doc"""
    automorphism_group_generators(L::ZLat; ambient_representation::Bool = true)

Returns generators of the automorphism group of $L$. By default, the
automorphisms are acting on the coordinate vectors of lattice elements.
If `ambient_representation = true`, the automorphisms act on elements in the
ambient space of `L`.
"""
function automorphism_group_generators(L::ZLat; check::Bool = true,
                                                ambient_representation::Bool = true)

  @req isdefinite(L) "The lattice must be definite"
  assert_has_automorphisms(L)

  gens = L.automorphism_group_generators
  if !ambient_representation
    return fmpq_mat[ change_base_ring(FlintQQ, g) for g in gens]
  else
    # Now translate to get the automorphisms with respect to basis_matrix(L)
    bm = basis_matrix(L)
    bminv = inv(bm)
    gens = L.automorphism_group_generators
    V = ambient_space(L)
    if rank(L) == rank(V)
      res = fmpq_mat[bminv * change_base_ring(FlintQQ, g) * bm for g in gens]
    else
      # Extend trivially to the orthogonal complement of the rational span
      !isregular(V) &&
        throw(error(
          """Can compute ambient representation only if ambient space is
             regular"""))
      C = orthogonal_complement(V, basis_matrix(L))
      C = vcat(basis_matrix(L), C)
      Cinv = inv(C)
      D = identity_matrix(FlintQQ, rank(V) - rank(L))
      res = fmpq_mat[Cinv * diagonal_matrix(change_base_ring(FlintQQ, g), D) * C for g in gens]
    end
    if check
      for g in res
        @assert g * gram_matrix(V) * g' == gram_matrix(V)
      end
    end
    return res
  end
end

function automorphism_group_order(L::ZLat)
  assert_has_automorphisms(L)
  return L.automorphism_group_order
end

@doc Markdown.doc"""
    isisometric(L::ZLat, M::ZLat; ambient_representation::Bool = true
                                  check::Bool = true)
        -> (Bool, MatElem)

Tests if $L$ and $M$ are isometric. If this is the case, the second return value
is an isometry $T$ from $L$ to $M$.

By default, that isometry is represented with respect to the bases of the
ambient spaces, that is, $T V_M T^t = V_L$ where $V_L$ and $V_M$ are the gram
matrices of the ambient spaces of $L$ and $M$ respectively. If
`ambient_representation = true`, then the isometry is represented with respect
to the bases of $L$ and $M$, that is, $T G_M T^t = G_L$ where $G_M$ and $G_L$ are
the gram matrices of $L$ and $M$ respectively.
"""
function isisometric(L::ZLat, M::ZLat; ambient_representation::Bool = true,
                                       check::Bool = true)
  GL = gram_matrix(L)
  dL = denominator(GL)
  GLint = change_base_ring(FlintZZ, dL * GL)
  cL = content(GLint)
  GLint = divexact(GLint, cL)

  GM = gram_matrix(M)
  dM = denominator(GM)
  GMint = change_base_ring(FlintZZ, dM * GM)
  cM = content(GMint)
  GMint = divexact(GMint, cM)

  # GLint, GMint are integral, primitive scalings of GL and GM
  # If they are isometric, then the scalars must be identitcal.
  if dL//cL != dM//cM
    return false, zero_matrix(FlintQQ, 0, 0)
  end

  # Now compute LLL reduces gram matrices

  GLlll, TL = lll_gram_with_transform(GLint)
  @assert TL * change_base_ring(FlintZZ, GL) * TL' * dL == GLlll * cL
  GMlll, TM = lll_gram_with_transform(GMint)
  @assert TM * change_base_ring(FlintZZ, GM) * TM' * dM == GMlll * cM

  # Setup for Plesken--Souvignier

  G1 = fmpz_mat[GLlll]
  G2 = fmpz_mat[GMlll]

  fl, CLsmall, CMsmall = _try_iso_setup_small(G1, G2)
  if fl
    b, _T = isometry(CLsmall, CMsmall)
    T = matrix(FlintZZ, _T)
  else
    CL, CM = _iso_setup(fmpz_mat[GLlll], fmpz_mat[GMlll])
    b, T = isometry(CL, CM)
  end

  if b
    T = change_base_ring(FlintQQ, inv(TL)*T*TM)
    if !ambient_representation
      if check
        @assert T * gram_matrix(M) * T' == gram_matrix(L)
      end
      return true, T
    else
      V = ambient_space(L)
      W = ambient_space(L)
      if rank(L) == rank(V)
        T = inv(basis_matrix(L)) * T * basis_matrix(M)
      else
        (!isregular(V) || !isregular(W)) &&
          throw(error(
            """Can compute ambient representation only if ambient space is
               regular"""))
          (rank(V) != rank(W)) &&
          throw(error(
            """Can compute ambient representation only if ambient spaces
            have the same dimension."""))

        CV = orthogonal_complement(V, basis_matrix(L))
        CV = vcat(basis_matrix(L), CV)
        CW = orthogonal_complement(V, basis_matrix(M))
        CW = vcat(basis_matrix(M), CW)
        D = identity_matrix(FlintQQ, rank(V) - rank(L))
        T = inv(CV) * diagonal_matrix(T, D) * CW
      end
      if check
        @assert T * gram_matrix(ambient_space(M))  * T' ==
                  gram_matrix(ambient_space(L))
      end
      return true, T
    end
  else
    return false, zero_matrix(FlintQQ, 0, 0)
  end
end

################################################################################
#
#  Root lattice
#
################################################################################

@doc Markdown.doc"""
    root_lattice(R::Symbol, n::Int)

Determine the root lattice of type `R` with parameter `n`. At the moment only
type `:A` is supported.
"""
function root_lattice(R::Symbol, n::Int)
  if R === :A
    return Zlattice(gram = _root_lattice_A(n))
  else
    error("Type (:$R) must be :A")
  end
end

function _root_lattice_A(n::Int)
  n < 0 && error("Parameter ($n) for root lattice of type :A must be positive")
  z = zero_matrix(FlintQQ, n, n)
  for i in 1:n
    z[i, i] = 2
    if i > 1
      z[i, i - 1] = -1
    end
    if i < n
      z[i, i + 1] = -1
    end
  end
  return z
end

################################################################################
#
#  Dual lattice
#
################################################################################

function dual(L::ZLat)
  G = gram_matrix(ambient_space(L))
  Gi = inv(G)
  new_bmat = transpose(inv(basis_matrix(L)) * Gi)
  return lattice(ambient_space(L), new_bmat)
end


text:
function get防区编号和FL(sc,t0)
   global    CalChn1_2Obj
    meters    =  CalChn1_2Obj.CfgInfoObj.meterperSeg
   segIDChn2 =  CalChn1_2Obj.CfgInfoObj.chn2SegBegIdx
   chnid = 1
   fl = 0
   segid = floor(Int, sc/10)
   if segid >= segIDChn2
        chnid =2
        fl = (segid - segIDChn2 )  * meters
    else
        chnid = 1
        fl = (segid )  * meters
    end
    return ( chnid,fl , UInt2DTStr( t0 ))
# @info "----SegParaAlarmProc.CalChn1_2Obj.CfgInfoObj.meterperSeg-- ???:",SegParaAlarmProc.CalChn1_2Obj.CfgInfoObj.meterperSeg
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool