text:
"""
    Mattsson2017(version::Symbol)

Coefficients of the upwind SBP operators given in
- Mattsson (2017)
  Diagonal-norm upwind SBP operators.
  Journal of Computational Physics 335, pp. 283-310.

You can choose between the different versions `:central`, `:plus`, and `:minus`.
"""
struct Mattsson2017 <: SourceOfCoefficients
  kind::Symbol

  function Mattsson2017(kind::Symbol)
    if (kind !== :plus) && (kind !== :minus) && (kind !== :central)
      throw(ArgumentError("The only choices are :plus, :minus, and :central, not :$kind."))
    end
    new(kind)
  end
end

function Base.show(io::IO, source::Mattsson2017)
  if get(io, :compact, false)
    summary(io, source)
  else
    print(io,
        "Mattsson (2017) \n",
        "  Diagonal-norm upwind SBP operators. \n",
        "  Journal of Computational Physics 335, pp. 283-310. \n",
        "  (upwind coefficients ", source.kind, ")")
  end
end

function first_derivative_coefficients(source::Mattsson2017, order::Int,
                                       T=Float64, mode=FastMode())
    if order == 2
        left_boundary_plus = (
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3),
                                                    T(5),
                                                    T(-2), )),
            DerivativeCoefficientRow{T,1,4}(SVector(T(-1//5),
                                                    T(-1),
                                                    T(8//5),
                                                    T(-2//5), )),
        )
        right_boundary_plus = (
            DerivativeCoefficientRow{T,1,2}(SVector(T(1),
                                                    T(-1), )),
            DerivativeCoefficientRow{T,1,2}(SVector(T(1),
                                                    T(-1), )),
        )
        upper_coef_plus = SVector( T(2),
                                   T(-1//2), )
        central_coef_plus = T(-3//2)
        lower_coef_plus = SVector{0,T}()
        left_weights = SVector( T(1//4),
                                T(5//4), )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        left_boundary_minus = (
            DerivativeCoefficientRow{T,1,2}(SVector(T(-1),
                                                    T(1), )),
            DerivativeCoefficientRow{T,1,2}(SVector(T(-1),
                                                    T(1), )),
        )
        right_boundary_minus = (
            DerivativeCoefficientRow{T,1,3}(SVector(T(3),
                                                    T(-5),
                                                    T(2), )),
            DerivativeCoefficientRow{T,1,4}(SVector(T(1//5),
                                                    T(1),
                                                    T(-8//5),
                                                    T(2//5), )),
        )
        upper_coef_minus     = .- lower_coef_plus
        central_coef_minus   = .- central_coef_plus
        lower_coef_minus     = .- upper_coef_plus

        left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
        right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
        upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
        central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
        lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

        if source.kind === :plus
          left_boundary  = left_boundary_plus
          right_boundary = right_boundary_plus
          upper_coef     = upper_coef_plus
          central_coef   = central_coef_plus
          lower_coef     = lower_coef_plus
        elseif source.kind === :minus
          left_boundary  = left_boundary_minus
          right_boundary = right_boundary_minus
          upper_coef     = upper_coef_minus
          central_coef   = central_coef_minus
          lower_coef     = lower_coef_minus
        elseif source.kind === :central
          left_boundary  = left_boundary_central
          right_boundary = right_boundary_central
          upper_coef     = upper_coef_central
          central_coef   = central_coef_central
          lower_coef     = lower_coef_central
        end
        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 3
      left_boundary_plus = (
          DerivativeCoefficientRow{T,1,3}(SVector(T(-7//5),
                                                  T(9//5),
                                                  T(-2//5), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(-5//13),
                                                  T(-5//13),
                                                  T(12//13),
                                                  T(-2//13), )),
      )
      right_boundary_plus = (
          DerivativeCoefficientRow{T,1,2}(SVector(T(1),
                                                  T(-1), )),
          DerivativeCoefficientRow{T,1,3}(SVector(T(9//13),
                                                  T(-5//13),
                                                  T(-4//13), )),
      )
      upper_coef_plus = SVector( T(1),
                                 T(-1//6), )
      central_coef_plus = T(-1//2)
      lower_coef_plus = SVector( T(-1//3), )
      left_weights = SVector( T(5//12),
                              T(13//12), )
      right_weights = left_weights
      left_boundary_derivatives = Tuple{}()
      right_boundary_derivatives = left_boundary_derivatives

      left_boundary_minus = (
          DerivativeCoefficientRow{T,1,2}(SVector(T(-1),
                                                  T(1), )),
          DerivativeCoefficientRow{T,1,3}(SVector(T(-9//13),
                                                  T(5//13),
                                                  T(4//13), )),
      )
      right_boundary_minus = (
          DerivativeCoefficientRow{T,1,3}(SVector(T(7//5),
                                                  T(-9//5),
                                                  T(2//5), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(5//13),
                                                  T(5//13),
                                                  T(-12//13),
                                                  T(2//13), )),
      )
      upper_coef_minus     = .- lower_coef_plus
      central_coef_minus   = .- central_coef_plus
      lower_coef_minus     = .- upper_coef_plus

      left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
      right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
      upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
      central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
      lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

      if source.kind === :plus
        left_boundary  = left_boundary_plus
        right_boundary = right_boundary_plus
        upper_coef     = upper_coef_plus
        central_coef   = central_coef_plus
        lower_coef     = lower_coef_plus
      elseif source.kind === :minus
        left_boundary  = left_boundary_minus
        right_boundary = right_boundary_minus
        upper_coef     = upper_coef_minus
        central_coef   = central_coef_minus
        lower_coef     = lower_coef_minus
      elseif source.kind === :central
        left_boundary  = left_boundary_central
        right_boundary = right_boundary_central
        upper_coef     = upper_coef_central
        central_coef   = central_coef_central
        lower_coef     = lower_coef_central
      end
      DerivativeCoefficients(left_boundary, right_boundary,
                              left_boundary_derivatives, right_boundary_derivatives,
                              lower_coef, central_coef, upper_coef,
                              left_weights, right_weights, mode, 1, order, source)
    elseif order == 4
      left_boundary_plus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(-75//49),
                                                  T(205//98),
                                                  T(-29//49),
                                                  T(3//98), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(-169//366),
                                                  T(-11//61),
                                                  T(99//122),
                                                  T(-43//183),
                                                  T(4//61), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(11//123),
                                                  T(-39//82),
                                                  T(-29//41),
                                                  T(389//246),
                                                  T(-24//41),
                                                  T(4//41), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(9//298),
                                                  T(-11//149),
                                                  T(-65//298),
                                                  T(-117//149),
                                                  T(216//149),
                                                  T(-72//149),
                                                  T(12//149), )),
      )
      right_boundary_plus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(69//49),
                                                  T(-169//98),
                                                  T(11//49),
                                                  T(9//98), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(205//366),
                                                  T(-11//61),
                                                  T(-39//122),
                                                  T(-11//183), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(-29//123),
                                                  T(99//82),
                                                  T(-29//41),
                                                  T(-65//246), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(3//298),
                                                  T(-43//149),
                                                  T(389//298),
                                                  T(-117//149),
                                                  T(-36//149), )),
      )
      upper_coef_plus = SVector( T(3//2),
                                 T(-1//2),
                                 T(1//12), )
      central_coef_plus = T(-5//6)
      lower_coef_plus = SVector( T(-1//4), )
      left_weights = SVector( T(49//144),
                              T(61//48),
                              T(41//48),
                              T(149//144), )
      right_weights = left_weights
      left_boundary_derivatives = Tuple{}()
      right_boundary_derivatives = left_boundary_derivatives

      left_boundary_minus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(-69//49),
                                                  T(169//98),
                                                  T(-11//49),
                                                  T(-9//98), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(-205//366),
                                                  T(11//61),
                                                  T(39//122),
                                                  T(11//183), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(29//123),
                                                  T(-99//82),
                                                  T(29//41),
                                                  T(65//246), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(-3//298),
                                                  T(43//149),
                                                  T(-389//298),
                                                  T(117//149),
                                                  T(36//149), )),
      )
      right_boundary_minus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(75//49),
                                                  T(-205//98),
                                                  T(29//49),
                                                  T(-3//98), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(169//366),
                                                  T(11//61),
                                                  T(-99//122),
                                                  T(43//183),
                                                  T(-4//61), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-11//123),
                                                  T(39//82),
                                                  T(29//41),
                                                  T(-389//246),
                                                  T(24//41),
                                                  T(-4//41), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-9//298),
                                                  T(11//149),
                                                  T(65//298),
                                                  T(117//149),
                                                  T(-216//149),
                                                  T(72//149),
                                                  T(-12//149), )),
      )
      upper_coef_minus     = .- lower_coef_plus
      central_coef_minus   = .- central_coef_plus
      lower_coef_minus     = .- upper_coef_plus

      left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
      right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
      upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
      central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
      lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

      if source.kind === :plus
        left_boundary  = left_boundary_plus
        right_boundary = right_boundary_plus
        upper_coef     = upper_coef_plus
        central_coef   = central_coef_plus
        lower_coef     = lower_coef_plus
      elseif source.kind === :minus
        left_boundary  = left_boundary_minus
        right_boundary = right_boundary_minus
        upper_coef     = upper_coef_minus
        central_coef   = central_coef_minus
        lower_coef     = lower_coef_minus
      elseif source.kind === :central
        left_boundary  = left_boundary_central
        right_boundary = right_boundary_central
        upper_coef     = upper_coef_central
        central_coef   = central_coef_central
        lower_coef     = lower_coef_central
      end
      DerivativeCoefficients(left_boundary, right_boundary,
                              left_boundary_derivatives, right_boundary_derivatives,
                              lower_coef, central_coef, upper_coef,
                              left_weights, right_weights, mode, 1, order, source)
    elseif order == 5
      left_boundary_plus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(-366//251),
                                                  T(941//502),
                                                  T(-94//251),
                                                  T(-21//502), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(-869//1794),
                                                  T(-22//299),
                                                  T(375//598),
                                                  T(-86//897),
                                                  T(8//299), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(58//633),
                                                  T(-255//422),
                                                  T(-58//211),
                                                  T(1309//1266),
                                                  T(-60//211),
                                                  T(8//211), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(45//1478),
                                                  T(-22//739),
                                                  T(-661//1478),
                                                  T(-234//739),
                                                  T(720//739),
                                                  T(-180//739),
                                                  T(24//739), )),
      )
      right_boundary_plus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(354//251),
                                                  T(-869//502),
                                                  T(58//251),
                                                  T(45//502), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(941//1794),
                                                  T(-22//299),
                                                  T(-255//598),
                                                  T(-22//897), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(-94//633),
                                                  T(375//422),
                                                  T(-58//211),
                                                  T(-661//1266),
                                                  T(12//211), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-21//1478),
                                                  T(-86//739),
                                                  T(1309//1478),
                                                  T(-234//739),
                                                  T(-360//739),
                                                  T(36//739), )),
      )
      upper_coef_plus = SVector( T(1),
                                 T(-1//4),
                                 T(1//30), )
      central_coef_plus = T(-1//3)
      lower_coef_plus = SVector( T(-1//2),
                                 T(1//20), )
      left_weights = SVector( T(251//720),
                              T(299//240),
                              T(211//240),
                              T(739//720), )
      right_weights = left_weights
      left_boundary_derivatives = Tuple{}()
      right_boundary_derivatives = left_boundary_derivatives

      left_boundary_minus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(-354//251),
                                                  T(869//502),
                                                  T(-58//251),
                                                  T(-45//502), )),
          DerivativeCoefficientRow{T,1,4}(SVector(T(-941//1794),
                                                  T(22//299),
                                                  T(255//598),
                                                  T(22//897), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(94//633),
                                                  T(-375//422),
                                                  T(58//211),
                                                  T(661//1266),
                                                  T(-12//211), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(21//1478),
                                                  T(86//739),
                                                  T(-1309//1478),
                                                  T(234//739),
                                                  T(360//739),
                                                  T(-36//739), )),
      )
      right_boundary_minus = (
          DerivativeCoefficientRow{T,1,4}(SVector(T(366//251),
                                                  T(-941//502),
                                                  T(94//251),
                                                  T(21//502), )),
          DerivativeCoefficientRow{T,1,5}(SVector(T(869//1794),
                                                  T(22//299),
                                                  T(-375//598),
                                                  T(86//897),
                                                  T(-8//299), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-58//633),
                                                  T(255//422),
                                                  T(58//211),
                                                  T(-1309//1266),
                                                  T(60//211),
                                                  T(-8//211), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-45//1478),
                                                  T(22//739),
                                                  T(661//1478),
                                                  T(234//739),
                                                  T(-720//739),
                                                  T(180//739),
                                                  T(-24//739), )),
      )
      upper_coef_minus     = .- lower_coef_plus
      central_coef_minus   = .- central_coef_plus
      lower_coef_minus     = .- upper_coef_plus

      left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
      right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
      upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
      central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
      lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

      if source.kind === :plus
        left_boundary  = left_boundary_plus
        right_boundary = right_boundary_plus
        upper_coef     = upper_coef_plus
        central_coef   = central_coef_plus
        lower_coef     = lower_coef_plus
      elseif source.kind === :minus
        left_boundary  = left_boundary_minus
        right_boundary = right_boundary_minus
        upper_coef     = upper_coef_minus
        central_coef   = central_coef_minus
        lower_coef     = lower_coef_minus
      elseif source.kind === :central
        left_boundary  = left_boundary_central
        right_boundary = right_boundary_central
        upper_coef     = upper_coef_central
        central_coef   = central_coef_central
        lower_coef     = lower_coef_central
      end
      DerivativeCoefficients(left_boundary, right_boundary,
                              left_boundary_derivatives, right_boundary_derivatives,
                              lower_coef, central_coef, upper_coef,
                              left_weights, right_weights, mode, 1, order, source)
    elseif order == 6
      left_boundary_plus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(-58148100//36496453),
                                                  T(1146190567//547446795),
                                                  T(-14369571//52137790),
                                                  T(-55265831//182482265),
                                                  T(26269819//1094893590),
                                                  T(9858004//182482265), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-1116490567//2422752675),
                                                  T(-954612//32303369),
                                                  T(190538869//484550535),
                                                  T(102705469//969101070),
                                                  T(4964892//161516845),
                                                  T(-191689861//4845505350), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(9869571//102452500),
                                                  T(-135385429//215150250),
                                                  T(-2198412//7171675),
                                                  T(45137333//35858375),
                                                  T(-253641811//430300500),
                                                  T(70665929//358583750),
                                                  T(-72//2675), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(66965831//723199750),
                                                  T(-208765789//867839700),
                                                  T(-17623253//72319975),
                                                  T(-657684//2066285),
                                                  T(410905829//433919850),
                                                  T(-477953317//1446399500),
                                                  T(576//5395),
                                                  T(-72//5395), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(-49219819//3153258150),
                                                  T(3519588//105108605),
                                                  T(26422771//630651630),
                                                  T(-141938309//315325815),
                                                  T(-12476988//21021721),
                                                  T(2217185207//1576629075),
                                                  T(-4320//7841),
                                                  T(1152//7841),
                                                  T(-144//7841), )),
          DerivativeCoefficientRow{T,1,10}(SVector(T(-9498004//587634985),
                                                  T(142906261//3525809910),
                                                  T(-3137129//587634985),
                                                  T(-29884283//1175269970),
                                                  T(-630168407//1762904955),
                                                  T(-9609300//16789571),
                                                  T(57600//43837),
                                                  T(-21600//43837),
                                                  T(5760//43837),
                                                  T(-720//43837), )),
      )
      right_boundary_plus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(57671100//36496453),
                                                  T(-1116490567//547446795),
                                                  T(9869571//52137790),
                                                  T(66965831//182482265),
                                                  T(-49219819//1094893590),
                                                  T(-9498004//182482265), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(1146190567//2422752675),
                                                  T(-954612//32303369),
                                                  T(-135385429//484550535),
                                                  T(-208765789//969101070),
                                                  T(3519588//161516845),
                                                  T(142906261//4845505350), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-14369571//102452500),
                                                  T(190538869//215150250),
                                                  T(-2198412//7171675),
                                                  T(-17623253//35858375),
                                                  T(26422771//430300500),
                                                  T(-3137129//358583750), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-55265831//723199750),
                                                  T(102705469//867839700),
                                                  T(45137333//72319975),
                                                  T(-657684//2066285),
                                                  T(-141938309//433919850),
                                                  T(-2298791//111261500), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(26269819//3153258150),
                                                  T(4964892//105108605),
                                                  T(-253641811//630651630),
                                                  T(410905829//315325815),
                                                  T(-12476988//21021721),
                                                  T(-630168407//1576629075),
                                                  T(288//7841), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(9858004//587634985),
                                                  T(-191689861//3525809910),
                                                  T(70665929//587634985),
                                                  T(-477953317//1175269970),
                                                  T(2217185207//1762904955),
                                                  T(-9609300//16789571),
                                                  T(-17280//43837),
                                                  T(1440//43837), )),
      )
      upper_coef_plus = SVector( T(4//3),
                                 T(-1//2),
                                 T(2//15),
                                 T(-1//60), )
      central_coef_plus = T(-7//12)
      lower_coef_plus = SVector( T(-2//5),
                                 T(1//30), )
      left_weights = SVector( T(13613//43200),
                              T(12049//8640),
                              T(535//864),
                              T(1079//864),
                              T(7841//8640),
                              T(43837//43200), )
      right_weights = left_weights
      left_boundary_derivatives = Tuple{}()
      right_boundary_derivatives = left_boundary_derivatives

      left_boundary_minus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(-57671100//36496453),
                                                  T(1116490567//547446795),
                                                  T(-9869571//52137790),
                                                  T(-66965831//182482265),
                                                  T(49219819//1094893590),
                                                  T(9498004//182482265), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-1146190567//2422752675),
                                                  T(954612//32303369),
                                                  T(135385429//484550535),
                                                  T(208765789//969101070),
                                                  T(-3519588//161516845),
                                                  T(-142906261//4845505350), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(14369571//102452500),
                                                  T(-190538869//215150250),
                                                  T(2198412//7171675),
                                                  T(17623253//35858375),
                                                  T(-26422771//430300500),
                                                  T(3137129//358583750), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(55265831//723199750),
                                                  T(-102705469//867839700),
                                                  T(-45137333//72319975),
                                                  T(657684//2066285),
                                                  T(141938309//433919850),
                                                  T(2298791//111261500), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-26269819//3153258150),
                                                  T(-4964892//105108605),
                                                  T(253641811//630651630),
                                                  T(-410905829//315325815),
                                                  T(12476988//21021721),
                                                  T(630168407//1576629075),
                                                  T(-288//7841), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(-9858004//587634985),
                                                  T(191689861//3525809910),
                                                  T(-70665929//587634985),
                                                  T(477953317//1175269970),
                                                  T(-2217185207//1762904955),
                                                  T(9609300//16789571),
                                                  T(17280//43837),
                                                  T(-1440//43837), )),
      )
      right_boundary_minus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(58148100//36496453),
                                                  T(-1146190567//547446795),
                                                  T(14369571//52137790),
                                                  T(55265831//182482265),
                                                  T(-26269819//1094893590),
                                                  T(-9858004//182482265), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(1116490567//2422752675),
                                                  T(954612//32303369),
                                                  T(-190538869//484550535),
                                                  T(-102705469//969101070),
                                                  T(-4964892//161516845),
                                                  T(191689861//4845505350), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-9869571//102452500),
                                                  T(135385429//215150250),
                                                  T(2198412//7171675),
                                                  T(-45137333//35858375),
                                                  T(253641811//430300500),
                                                  T(-70665929//358583750),
                                                  T(72//2675), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(-66965831//723199750),
                                                  T(208765789//867839700),
                                                  T(17623253//72319975),
                                                  T(657684//2066285),
                                                  T(-410905829//433919850),
                                                  T(477953317//1446399500),
                                                  T(-576//5395),
                                                  T(72//5395), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(49219819//3153258150),
                                                  T(-3519588//105108605),
                                                  T(-26422771//630651630),
                                                  T(141938309//315325815),
                                                  T(12476988//21021721),
                                                  T(-2217185207//1576629075),
                                                  T(4320//7841),
                                                  T(-1152//7841),
                                                  T(144//7841), )),
          DerivativeCoefficientRow{T,1,10}(SVector(T(9498004//587634985),
                                                  T(-142906261//3525809910),
                                                  T(3137129//587634985),
                                                  T(29884283//1175269970),
                                                  T(630168407//1762904955),
                                                  T(9609300//16789571),
                                                  T(-57600//43837),
                                                  T(21600//43837),
                                                  T(-5760//43837),
                                                  T(720//43837), )),
      )
      upper_coef_minus     = .- lower_coef_plus
      central_coef_minus   = .- central_coef_plus
      lower_coef_minus     = .- upper_coef_plus

      left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
      right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
      upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
      central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
      lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

      if source.kind === :plus
        left_boundary  = left_boundary_plus
        right_boundary = right_boundary_plus
        upper_coef     = upper_coef_plus
        central_coef   = central_coef_plus
        lower_coef     = lower_coef_plus
      elseif source.kind === :minus
        left_boundary  = left_boundary_minus
        right_boundary = right_boundary_minus
        upper_coef     = upper_coef_minus
        central_coef   = central_coef_minus
        lower_coef     = lower_coef_minus
      elseif source.kind === :central
        left_boundary  = left_boundary_central
        right_boundary = right_boundary_central
        upper_coef     = upper_coef_central
        central_coef   = central_coef_central
        lower_coef     = lower_coef_central
      end
      DerivativeCoefficients(left_boundary, right_boundary,
                              left_boundary_derivatives, right_boundary_derivatives,
                              lower_coef, central_coef, upper_coef,
                              left_weights, right_weights, mode, 1, order, source)
    elseif order == 7
      left_boundary_plus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(-81216540//51172247),
                                                  T(1587945773//767583705),
                                                  T(-17337249//73103210),
                                                  T(-84398989//255861235),
                                                  T(48781961//1535167410),
                                                  T(13716476//255861235), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-1570125773//3386062785),
                                                  T(-2863836//225737519),
                                                  T(240029831//677212557),
                                                  T(202934303//1354425114),
                                                  T(1418484//225737519),
                                                  T(-231357719//6772125570), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(14637249//144536540),
                                                  T(-206937767//303526734),
                                                  T(-6595236//50587789),
                                                  T(49602727//50587789),
                                                  T(-218919665//607053468),
                                                  T(51815011//505877890),
                                                  T(-216//18869), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(91418989//1008619010),
                                                  T(-266570495//1210342812),
                                                  T(-33094279//100861901),
                                                  T(-1973052//14408843),
                                                  T(440626231//605171406),
                                                  T(-365711063//2017238020),
                                                  T(2016//37621),
                                                  T(-216//37621), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(-62551961//4426143330),
                                                  T(9588//385217),
                                                  T(82588241//885228666),
                                                  T(-279245719//442614333),
                                                  T(-37430964//147538111),
                                                  T(2312302333//2213071665),
                                                  T(-18144//55031),
                                                  T(4032//55031),
                                                  T(-432//55031), )),
          DerivativeCoefficientRow{T,1,10}(SVector(T(-13500476//822302915),
                                                  T(202087559//4933817490),
                                                  T(-11297731//822302915),
                                                  T(61008503//1644605830),
                                                  T(-1360092253//2466908745),
                                                  T(-5765580//23494369),
                                                  T(60480//61343),
                                                  T(-18144//61343),
                                                  T(4032//61343),
                                                  T(-432//61343), )),
      )
      right_boundary_plus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(80930340//51172247),
                                                  T(-1570125773//767583705),
                                                  T(14637249//73103210),
                                                  T(91418989//255861235),
                                                  T(-62551961//1535167410),
                                                  T(-13500476//255861235), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(1587945773//3386062785),
                                                  T(-2863836//225737519),
                                                  T(-206937767//677212557),
                                                  T(-266570495//1354425114),
                                                  T(9588//589393),
                                                  T(202087559//6772125570), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-17337249//144536540),
                                                  T(240029831//303526734),
                                                  T(-6595236//50587789),
                                                  T(-33094279//50587789),
                                                  T(82588241//607053468),
                                                  T(-11297731//505877890), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-84398989//1008619010),
                                                  T(202934303//1210342812),
                                                  T(49602727//100861901),
                                                  T(-1973052//14408843),
                                                  T(-279245719//605171406),
                                                  T(61008503//2017238020),
                                                  T(-288//37621), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(48781961//4426143330),
                                                  T(1418484//147538111),
                                                  T(-218919665//885228666),
                                                  T(440626231//442614333),
                                                  T(-37430964//147538111),
                                                  T(-1360092253//2213071665),
                                                  T(6048//55031),
                                                  T(-576//55031), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(13716476//822302915),
                                                  T(-231357719//4933817490),
                                                  T(51815011//822302915),
                                                  T(-365711063//1644605830),
                                                  T(2312302333//2466908745),
                                                  T(-5765580//23494369),
                                                  T(-36288//61343),
                                                  T(6048//61343),
                                                  T(-576//61343), )),
      )
      upper_coef_plus = SVector( T(1),
                                 T(-3//10),
                                 T(1//15),
                                 T(-1//140), )
      central_coef_plus = T(-1//4)
      lower_coef_plus = SVector( T(-3//5),
                                 T(1//10),
                                 T(-1//105), )
      left_weights = SVector( T(19087//60480),
                              T(84199//60480),
                              T(18869//30240),
                              T(37621//30240),
                              T(55031//60480),
                              T(61343//60480), )
      right_weights = left_weights
      left_boundary_derivatives = Tuple{}()
      right_boundary_derivatives = left_boundary_derivatives

      left_boundary_minus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(-80930340//51172247),
                                                  T(1570125773//767583705),
                                                  T(-14637249//73103210),
                                                  T(-91418989//255861235),
                                                  T(62551961//1535167410),
                                                  T(13500476//255861235), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(-1587945773//3386062785),
                                                  T(2863836//225737519),
                                                  T(206937767//677212557),
                                                  T(266570495//1354425114),
                                                  T(-9588//589393),
                                                  T(-202087559//6772125570), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(17337249//144536540),
                                                  T(-240029831//303526734),
                                                  T(6595236//50587789),
                                                  T(33094279//50587789),
                                                  T(-82588241//607053468),
                                                  T(11297731//505877890), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(84398989//1008619010),
                                                  T(-202934303//1210342812),
                                                  T(-49602727//100861901),
                                                  T(1973052//14408843),
                                                  T(279245719//605171406),
                                                  T(-61008503//2017238020),
                                                  T(288//37621), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(-48781961//4426143330),
                                                  T(-1418484//147538111),
                                                  T(218919665//885228666),
                                                  T(-440626231//442614333),
                                                  T(37430964//147538111),
                                                  T(1360092253//2213071665),
                                                  T(-6048//55031),
                                                  T(576//55031), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(-13716476//822302915),
                                                  T(231357719//4933817490),
                                                  T(-51815011//822302915),
                                                  T(365711063//1644605830),
                                                  T(-2312302333//2466908745),
                                                  T(5765580//23494369),
                                                  T(36288//61343),
                                                  T(-6048//61343),
                                                  T(576//61343), )),
      )
      right_boundary_minus = (
          DerivativeCoefficientRow{T,1,6}(SVector(T(81216540//51172247),
                                                  T(-1587945773//767583705),
                                                  T(17337249//73103210),
                                                  T(84398989//255861235),
                                                  T(-48781961//1535167410),
                                                  T(-13716476//255861235), )),
          DerivativeCoefficientRow{T,1,6}(SVector(T(1570125773//3386062785),
                                                  T(2863836//225737519),
                                                  T(-240029831//677212557),
                                                  T(-202934303//1354425114),
                                                  T(-1418484//225737519),
                                                  T(231357719//6772125570), )),
          DerivativeCoefficientRow{T,1,7}(SVector(T(-14637249//144536540),
                                                  T(206937767//303526734),
                                                  T(6595236//50587789),
                                                  T(-49602727//50587789),
                                                  T(218919665//607053468),
                                                  T(-51815011//505877890),
                                                  T(216//18869), )),
          DerivativeCoefficientRow{T,1,8}(SVector(T(-91418989//1008619010),
                                                  T(266570495//1210342812),
                                                  T(33094279//100861901),
                                                  T(1973052//14408843),
                                                  T(-440626231//605171406),
                                                  T(365711063//2017238020),
                                                  T(-2016//37621),
                                                  T(216//37621), )),
          DerivativeCoefficientRow{T,1,9}(SVector(T(62551961//4426143330),
                                                  T(-9588//385217),
                                                  T(-82588241//885228666),
                                                  T(279245719//442614333),
                                                  T(37430964//147538111),
                                                  T(-2312302333//2213071665),
                                                  T(18144//55031),
                                                  T(-4032//55031),
                                                  T(432//55031), )),
          DerivativeCoefficientRow{T,1,10}(SVector(T(13500476//822302915),
                                                  T(-202087559//4933817490),
                                                  T(11297731//822302915),
                                                  T(-61008503//1644605830),
                                                  T(1360092253//2466908745),
                                                  T(5765580//23494369),
                                                  T(-60480//61343),
                                                  T(18144//61343),
                                                  T(-4032//61343),
                                                  T(432//61343), )),
      )
      upper_coef_minus     = .- lower_coef_plus
      central_coef_minus   = .- central_coef_plus
      lower_coef_minus     = .- upper_coef_plus

      left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
      right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
      upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
      central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
      lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

      if source.kind === :plus
        left_boundary  = left_boundary_plus
        right_boundary = right_boundary_plus
        upper_coef     = upper_coef_plus
        central_coef   = central_coef_plus
        lower_coef     = lower_coef_plus
      elseif source.kind === :minus
        left_boundary  = left_boundary_minus
        right_boundary = right_boundary_minus
        upper_coef     = upper_coef_minus
        central_coef   = central_coef_minus
        lower_coef     = lower_coef_minus
      elseif source.kind === :central
        left_boundary  = left_boundary_central
        right_boundary = right_boundary_central
        upper_coef     = upper_coef_central
        central_coef   = central_coef_central
        lower_coef     = lower_coef_central
      end
      DerivativeCoefficients(left_boundary, right_boundary,
                              left_boundary_derivatives, right_boundary_derivatives,
                              lower_coef, central_coef, upper_coef,
                              left_weights, right_weights, mode, 1, order, source)
    # elseif order == 6
    #   left_boundary_plus = (
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,9}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,10}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #   )
    #   right_boundary_plus = (
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #   )
    #   upper_coef_plus = SVector( T(),
    #                              T(),
    #                              T(), )
    #   central_coef_plus = T()
    #   lower_coef_plus = SVector( T(),
    #                              T(),
    #                              T(), )
    #   left_weights = SVector( T(),
    #                           T(),
    #                           T(),
    #                           T(),
    #                           T(),
    #                           T(), )
    #   right_weights = left_weights
    #   left_boundary_derivatives = Tuple{}()
    #   right_boundary_derivatives = left_boundary_derivatives

    #   left_boundary_minus = (
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #   )
    #   right_boundary_minus = (
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,9}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #       DerivativeCoefficientRow{T,1,10}(SVector(T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(),
    #                                               T(), )),
    #   )
    #   upper_coef_minus     = .- lower_coef_plus
    #   central_coef_minus   = .- central_coef_plus
    #   lower_coef_minus     = .- upper_coef_plus

    #   left_boundary_central  = (left_boundary_plus  .+ left_boundary_minus)  ./ 2
    #   right_boundary_central = (right_boundary_plus .+ right_boundary_minus) ./ 2
    #   upper_coef_central     = widening_plus(upper_coef_plus, upper_coef_minus) / 2
    #   central_coef_central   = (central_coef_plus   + central_coef_minus) / 2
    #   lower_coef_central     = widening_plus(lower_coef_plus, lower_coef_minus) / 2

    #   if source.kind === :plus
    #     left_boundary  = left_boundary_plus
    #     right_boundary = right_boundary_plus
    #     upper_coef     = upper_coef_plus
    #     central_coef   = central_coef_plus
    #     lower_coef     = lower_coef_plus
    #   elseif source.kind === :minus
    #     left_boundary  = left_boundary_minus
    #     right_boundary = right_boundary_minus
    #     upper_coef     = upper_coef_minus
    #     central_coef   = central_coef_minus
    #     lower_coef     = lower_coef_minus
    #   elseif source.kind === :central
    #     left_boundary  = left_boundary_central
    #     right_boundary = right_boundary_central
    #     upper_coef     = upper_coef_central
    #     central_coef   = central_coef_central
    #     lower_coef     = lower_coef_central
    #   end
    #   DerivativeCoefficients(left_boundary, right_boundary,
    #                           left_boundary_derivatives, right_boundary_derivatives,
    #                           lower_coef, central_coef, upper_coef,
    #                           left_weights, right_weights, mode, 1, order, source)
  else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

# function second_derivative_coefficients(source::Mattsson2017, order::Int, T=Float64, mode=FastMode())
#     if order == 2
#         left_boundary = (
#             # d1
#             DerivativeCoefficientRow{T,1,3}(SVector(T(1),
#                                                     T(-2),
#                                                     T(1) )),
#         )
#         right_boundary = left_boundary
#         upper_coef = SVector(T(1))
#         central_coef = T(-2)
#         lower_coef = upper_coef
#         left_weights = SVector(T(1//2))
#         right_weights = left_weights
#         left_boundary_derivatives = (
#             DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
#                                                     T(2),
#                                                     T(-1//2) )),
#         )
#         right_boundary_derivatives = (-left_boundary_derivatives[1], )

#         DerivativeCoefficients(left_boundary, right_boundary,
#                                 left_boundary_derivatives, right_boundary_derivatives,
#                                 lower_coef, central_coef, upper_coef,
#                                 left_weights, right_weights, mode, 2, order, source)
#     else
#         throw(ArgumentError("Order $order not implemented/derived."))
#     end
# end


text:
module SparseND

import ...BEAST
import FillArrays

struct MatrixOfConvolutions{T} <: AbstractArray{Vector{T},2}
    # banded::AbstractArray{T,3}
    convop::BEAST.ConvOp{T}
end

function Base.eltype(x::MatrixOfConvolutions{T}) where {T}
    Vector{T}
end
Base.size(x::MatrixOfConvolutions) = size(x.convop)[1:2]
function Base.getindex(x::MatrixOfConvolutions, m, n)
    return x.convop[m,n,:]
end

BEAST.convolve!(y, Z::MatrixOfConvolutions, x, csx, i, k_start, k_stop) = BEAST.convolve!(y, Z.convop, x, csx, i, k_start, k_stop)
BEAST.convolve!(y, Z::FillArrays.Zeros, x, csx, i, k_start, k_stop) = nothing
BEAST.convolve!(y, Z::FillArrays.Fill, x, csx, i, k_start, k_stop) = nothing

struct SpaceTimeData{T} <: AbstractArray{Vector{T},1}
    data::Array{T,2}
end

Base.eltype(x::SpaceTimeData{T}) where {T} = Vector{T}
Base.size(x::SpaceTimeData) = (size(x.data)[1],)
Base.getindex(x::SpaceTimeData, i::Int) = x.data[i,:]

end # module

ext: jl


text:
abstract type Field{T} end

# ******************************************************************************
# T
# ******************************************************************************
struct FieldT{T, TFFT} <: Field{T}
    Eu :: T
    Iu :: T
    w0 :: T
    E :: AbstractVector{Complex{T}}
    FFT :: TFFT
end

function Field(grid::GridT{T}, w0, E; Eu=1, Iu=1) where T
    (; arch) = grid

    E = Array{Complex{T}}(E)
    E = adapt_array(arch, E)

    FFT = plan_fft!(E)
    rsig2asig!(E, FFT)   # real signal -> analytic signal
    return FieldT{T, typeof(FFT)}(Eu, Iu, w0, E, FFT)
end

# ******************************************************************************
# RT
# ******************************************************************************
struct FieldRT{T, TDHT, TFFT} <: Field{T}
    Eu :: T
    Iu :: T
    w0 :: T
    E :: AbstractMatrix{Complex{T}}
    DHT :: TDHT
    FFT :: TFFT
end

function Field(grid::GridRT{T}, w0, E; Eu=1, Iu=1) where T
    (; arch, Nr, Nt, rmax) = grid

    E = Array{Complex{T}}(E)
    E = adapt_array(arch, E)

    Nwr = iseven(Nt) ? div(Nt, 2) : div(Nt+1, 2)
    region = (Nr, Nwr)
    DHT = plan_dht(rmax, E; region=region, save=false)

    FFT = plan_fft!(E, [2])
    rsig2asig!(E, FFT)   # real signal -> analytic signal

    return FieldRT{T, typeof(DHT), typeof(FFT)}(Eu, Iu, w0, E, DHT, FFT)
end

# ******************************************************************************
# Tools
# ******************************************************************************
function field2intensity(E, n0)
    return real(n0) * EPS0 * C0 * abs2(E) / 2
end

function intensity2field(I, n0)
    return sqrt(I / (real(n0) * EPS0 * C0 / 2))
end

function space2frequency!(E, DHT)
    DHT * E
    return nothing
end

function frequency2space!(E, DHT)
    DHT \ E
    return nothing
end

function time2frequency!(E, FFT)
    FFT \ E   # time -> frequency [exp(-i*w*t)]
    return nothing
end

function frequency2time!(E, FFT)
    FFT * E  # frequency -> time [exp(-i*w*t)]
    return nothing
end


text:
# functions related to studentized range distribution

# R implementations
using .RFunctions:
    srdistcdf,
    srdistccdf,
    srdistlogcdf,
    srdistlogccdf,
    srdistinvcdf,
    srdistinvccdf,
    srdistinvlogcdf,
    srdistinvlogccdf


text:
module Day2

abstract type Instruction end

struct Forward <: Instruction
    x::Int64
end

struct Down <: Instruction
    x::Int64
end

struct Up <: Instruction
    x::Int64
end

function instruction(s::String)
    (cmd, num) = split(s, " ")
    num = parse(Int64, num)
    if cmd == "forward"
        return Forward(num)
    end
    if cmd == "up"
        return Up(num)
    end
    if cmd == "down"
        return Down(num)
    end
    error("Bad instruction")
end

struct Position
    horizontal::Int64
    depth::Int64
end

Position() = Position(0, 0)

Base.:+(up::Up, pos::Position) = Position(pos.horizontal, pos.depth - up.x)
Base.:+(down::Down, pos::Position) = Position(pos.horizontal, pos.depth + down.x)
Base.:+(fwd::Forward, pos::Position) = Position(pos.horizontal + fwd.x, pos.depth)

read(f::String) = map(instruction, readlines(f))

run_part1(data_dir::String) = joinpath(data_dir, "day2.txt") |>
    read |>
    arr -> foldr(+, reverse(arr); init=Position()) |>
    pos -> pos.horizontal * pos.depth

struct PositionWithAim
    pos::Position
    aim::Int64
end

PositionWithAim() = PositionWithAim(Position(), 0)

Base.:+(up::Up, p::PositionWithAim) = PositionWithAim(p.pos, p.aim - up.x)
Base.:+(down::Down, p::PositionWithAim) = PositionWithAim(p.pos, p.aim + down.x)
Base.:+(fwd::Forward, p::PositionWithAim) = PositionWithAim(Down(p.aim * fwd.x) + (fwd + p.pos), p.aim)

run_part2(data_dir::String) = joinpath(data_dir, "day2.txt") |>
    read |>
    arr -> foldr(+, reverse(arr); init=PositionWithAim()) |>
    pos -> pos.pos.horizontal * pos.pos.depth

end


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

mutable struct IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec <: SwaggerModel
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec(;namespace=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec

const _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec = Dict{Symbol,Symbol}(Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec = Dict{Symbol,String}(Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1SelfSubjectRulesReviewSpec)
    true
end

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


text:
using DispatchError
using Base.Test

# this fails with error "boo!"
DispatchError.Foo(rand(3), rand(3), 1)


text:
if ENV["HOME"] == "/mnt/juliabox"
    Pkg.dir(path...)=joinpath("/home/jrun/.julia/v0.6",path...)
else
    for p in ("JLD","Knet")
        Pkg.installed(p) == nothing && Pkg.add(p)
    end
end

using JLD,Knet

info("Loading Shakespeare data")
include(Knet.dir("data","gutenberg.jl"))
trn,tst,chars = shakespeare()
shake_text = String(chars[vcat(trn,tst)])

info("Loading Shakespeare model")
isfile("shakespeare.jld") || download("http://people.csail.mit.edu/deniz/models/nlp-demos/shakespeare.jld","shakespeare.jld")
shake_model = load("shakespeare.jld","model")

info("Reading Julia files")
base = joinpath(Base.JULIA_HOME, Base.DATAROOTDIR, "julia", "base")
julia_text = ""
for (root,dirs,files) in walkdir(base)
    for f in files
        f[end-2:end] == ".jl" || continue
        julia_text *= readstring(joinpath(root,f))
    end
    # println((root,length(files),all(f->contains(f,".jl"),files)))
end

info("Loading Julia model")
isfile("juliacharlm.jld") || download("http://people.csail.mit.edu/deniz/models/nlp-demos/juliacharlm.jld","juliacharlm.jld")
julia_model = load("juliacharlm.jld","model")

# Given the current character, predict the next character
function predict(ws,xs,hx,cx;pdrop=0)
    r,wr,wx,wy,by = ws
    x = wx[:,xs]                                    # xs=(B,T) x=(X,B,T)
    x = dropout(x,pdrop)
    y,hy,cy = rnnforw(r,wr,x,hx,cx,hy=true,cy=true) # y=(H,B,T) hy=cy=(H,B,L)
    y = dropout(y,pdrop)
    y2 = reshape(y,size(y,1),size(y,2)*size(y,3))   # y2=(H,B*T)
    return wy*y2.+by, hy, cy
end

# Sample from trained model
function generate(model,n)
    function sample(y)
        p,r=Array(exp.(y-logsumexp(y))),rand()
        for j=1:length(p); (r -= p[j]) < 0 && return j; end
    end
    h,c = nothing,nothing
    chars = model[end]
    x = findfirst(chars,'\n')
    for i=1:n
        y,h,c = predict(model,[x],h,c)
        x = sample(y)
        print(chars[x])
    end
    println()
end

nothing

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
hexDic = Dict{Char,String}(
  '0' => "0000",
  '1' => "0001",
  '2' => "0010",
  '3' => "0011",
  '4' => "0100",
  '5' => "0101",
  '6' => "0110",
  '7' => "0111",
  '8' => "1000",
  '9' => "1001",
  'A' => "1010",
  'B' => "1011",
  'C' => "1100",
  'D' => "1101",
  'E' => "1110",
  'F' => "1111",
)

hex2bin(s::String) = string(map(x -> hexDic[x], collect(s))...)

# Converts from binary to decimal (binary given as array of bools).
function toInt(v :: Array{Bool,1}) :: Int
  num = 0
  for i in 1:length(v)
    num += (v[i] ? 1 : 0) * 2^(length(v)-i)
  end
  return num
end

toInt(s::String) :: Int = toInt(map(c -> (c == '0') ? false : true, collect(s)))

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

abstract type Packet end

struct NumberPacket <: Packet
  version :: Int
  type    :: Int
  value   :: Int
end

struct OperatorPacket <: Packet
  version  :: Int
  type     :: Int
  nPackets :: Int
  packets  :: Vector{Packet}
end

function parsePacket(s::String) :: Tuple{Packet,Int}
  println(" > Parsing Packet")
  (toInt(s[4:6]) == 4) ? (return parseNumber(s)) : (return parseOperator(s))
end

function parseNumber(s::String) :: Tuple{NumberPacket,Int}
  println("   > Parsing Number")
  version :: Int = toInt(s[1:3])
  type    :: Int = toInt(s[4:6])

  lastPack :: Bool = false
  pos :: Int = 7
  packs :: Vector{String} = []
  while !lastPack
    lastPack = (s[pos] == '0')
    push!(packs,s[pos+1:pos+4])
    pos += 5
  end
  value :: Int = toInt(string(packs...))
  return NumberPacket(version,type,value), pos
end

function parseOperator(s::String) :: Tuple{OperatorPacket,Int}
  println("   > Parsing Operator")
  version :: Int = toInt(s[1:3])
  type    :: Int = toInt(s[4:6])

  pos :: Int = 7
  packets :: Vector{Packet} = []
  if s[7] == '0'
    lenPackets :: Int = toInt(s[8:22])
    pos = 23
    while pos-23 < lenPackets
      packet, step = parsePacket(s[pos:end])
      pos += step-1
      push!(packets,packet)
    end
  else
    numPackets :: Int = toInt(s[8:18])
    pos = 19
    for i in 1:numPackets
      println(pos, " - ", s[pos:end])
      packet, step = parsePacket(s[pos:end])
      pos += step-1
      push!(packets,packet)
    end
  end

  return OperatorPacket(version, type, length(packets), packets), pos
end

versionSum(p::NumberPacket) :: Int = p.version

versionSum(p::OperatorPacket) :: Int = p.version + sum(map(x -> versionSum(x), p.packets))

compute(p::NumberPacket) :: Int = p.value

function compute(p::OperatorPacket) :: Int
  if p.type == 0
    return sum(map(x->compute(x), p.packets))
  elseif p.type == 1
    return prod(map(x->compute(x), p.packets))
  elseif p.type == 2
    return minimum(map(x->compute(x), p.packets))
  elseif p.type == 3
    return maximum(map(x->compute(x), p.packets))
  elseif p.type == 5
    return (compute(p.packets[1]) > compute(p.packets[2])) ? (return 1) : (return 0)
  elseif p.type == 6
    return (compute(p.packets[1]) < compute(p.packets[2])) ? (return 1) : (return 0)
  elseif p.type == 7
    return (compute(p.packets[1]) == compute(p.packets[2])) ? (return 1) : (return 0)
  else
    @error "Operator type not recognised!"
  end
end

function readInput()
  hex = readline("inputs/day16.txt")
  bin = hex2bin(hex)
  data, pos = parsePacket(bin)
  return data
end

p = readInput()

println("Problem 1:")
res = versionSum(p)
println(res)

prinln("Problem 2:")
res = compute(p)
println(res)


text:
# Histogram of Violation

using Plots

function plotDynConstVio(dynVio, constVio, useLog10 = true)

    if useLog10
        dV = log10.(dynVio)
        cV = log10.(constVio)
        header = "Log10 of "
    else
        dV = dynVio
        cV = constVio
        header = ""
    end

    plt = scatter(dV, cV)
    xlabel!(header * "Dynamics Constraint Violation")
    ylabel!(header * "Inequality Constraint Violation")
    title!("Clustering of Constraint Violation")

    display(plt)

    return plt
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
export BehaviorCloningPolicy

mutable struct BehaviorCloningPolicy{A} <: AbstractPolicy
    approximator::A
    explorer::AbstractExplorer
    sampler::BatchSampler{(:state, :action)}
    min_reservoir_history::Int
end

"""
    BehaviorCloningPolicy(;kw...)

# Keyword Arguments

- `approximator`: calculate the logits of possible actions directly
- `explorer=GreedyExplorer()`
- `batch_size::Int = 32`
- `min_reservoir_history::Int = 100`, number of transitions that should be experienced before updating the `approximator`.
- `rng = Random.GLOBAL_RNG`
"""
function BehaviorCloningPolicy(;
    approximator::A,
    explorer::AbstractExplorer = GreedyExplorer(),
    batch_size::Int = 32,
    min_reservoir_history::Int = 100,
    rng = Random.GLOBAL_RNG,
) where {A}
    sampler = BatchSampler{(:state, :action)}(batch_size; rng = rng)
    BehaviorCloningPolicy(approximator, explorer, sampler, min_reservoir_history)
end

function (p::BehaviorCloningPolicy)(env::AbstractEnv)
    s = state(env)
    s_batch = Flux.unsqueeze(s, ndims(s) + 1)
    s_batch = send_to_device(device(p.approximator), s_batch)
    logits = p.approximator(s_batch) |> vec |> send_to_host # drop dimension
    typeof(ActionStyle(env)) == MinimalActionSet ? p.explorer(logits) :
    p.explorer(logits, legal_action_space_mask(env))
end

function RLBase.update!(p::BehaviorCloningPolicy, batch::NamedTuple{(:state, :action)})
    s = send_to_device(device(p.approximator), batch.state)
    a = send_to_device(device(p.approximator), batch.action)
    m = p.approximator
    gs = gradient(params(m)) do
        ŷ = m(s)
        y = Flux.OneHotMatrix(a, size(ŷ, 1))
        logitcrossentropy(ŷ, y)
    end
    update!(m, gs)
end

function RLBase.update!(p::BehaviorCloningPolicy, t::Any)
    (length(t) <= p.min_reservoir_history || length(t) <= p.sampler.batch_size) && return

    _, batch = p.sampler(t)
    update!(p, batch)
end

function RLBase.prob(p::BehaviorCloningPolicy, env::AbstractEnv)
    s = state(env)
    m = p.approximator
    s_batch = send_to_device(device(m), Flux.unsqueeze(s, ndims(s) + 1))
    values = m(s_batch) |> vec |> send_to_host
    typeof(ActionStyle(env)) == MinimalActionSet ? prob(p.explorer, values) :
    prob(p.explorer, values, legal_action_space_mask(env))
end

function RLBase.prob(p::BehaviorCloningPolicy, env::AbstractEnv, action)
    A = action_space(env)
    P = prob(p, env)
    @assert length(A) == length(P)
    if A isa Base.OneTo
        P[action]
    else
        for (a, p) in zip(A, P)
            if a == action
                return p
            end
        end
        @error "action[$action] is not found in action space[$(action_space(env))]"
    end
end


text:
@testset "65.valid-number.jl" begin
    true_data = [
        "2",
        "0089",
        "-0.1",
        "+3.14",
        "4.",
        "-.9",
        "2e10",
        "-90E3",
        "3e+7",
        "+6e-1",
        "53.5e93",
        "-123.456e789",
    ]
    false_data = [
        "abc", "1a", "1e", "e3", "99e2.5", "--6", "-+3", "95a54e53", "3.5e+3.5e+3.5"
    ]
    @test all(is_number_65(s) for s in true_data)
    @test !any(is_number_65(s) for s in false_data)
end


text:
using Symbolics: get_variables
"""
    inputs(sys)

Return all variables that mare marked as inputs. See also [`unbound_inputs`](@ref)
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref)
"""
inputs(sys) = [filter(isinput, states(sys)); filter(isinput, parameters(sys))]

"""
    outputs(sys)

Return all variables that mare marked as outputs. See also [`unbound_outputs`](@ref)
See also [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
function outputs(sys)
    o = observed(sys)
    rhss = [eq.rhs for eq in o]
    lhss = [eq.lhs for eq in o]
    unique([filter(isoutput, states(sys))
            filter(isoutput, parameters(sys))
            filter(x -> x isa Term && isoutput(x), rhss) # observed can return equations with complicated expressions, we are only looking for single Terms
            filter(x -> x isa Term && isoutput(x), lhss)])
end

"""
    bound_inputs(sys)

Return inputs that are bound within the system, i.e., internal inputs
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
bound_inputs(sys) = filter(x -> is_bound(sys, x), inputs(sys))

"""
    unbound_inputs(sys)

Return inputs that are not bound within the system, i.e., external inputs
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
unbound_inputs(sys) = filter(x -> !is_bound(sys, x), inputs(sys))

"""
    bound_outputs(sys)

Return outputs that are bound within the system, i.e., internal outputs
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
bound_outputs(sys) = filter(x -> is_bound(sys, x), outputs(sys))

"""
    unbound_outputs(sys)

Return outputs that are not bound within the system, i.e., external outputs
See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
unbound_outputs(sys) = filter(x -> !is_bound(sys, x), outputs(sys))

"""
    is_bound(sys, u)

Determine whether or not input/output variable `u` is "bound" within the system, i.e., if it's to be considered internal to `sys`.
A variable/signal is considered bound if it appears in an equation together with variables from other subsystems.
The typical usecase for this function is to determine whether the input to an IO component is connected to another component,
or if it remains an external input that the user has to supply before simulating the system.

See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
"""
function is_bound(sys, u, stack = [])
    #=
    For observed quantities, we check if a variable is connected to something that is bound to something further out.
    In the following scenario
    julia> observed(syss)
        2-element Vector{Equation}:
        sys₊y(tv) ~ sys₊x(tv)
        y(tv) ~ sys₊x(tv)
    sys₊y(t) is bound to the outer y(t) through the variable sys₊x(t) and should thus return is_bound(sys₊y(t)) = true.
    When asking is_bound(sys₊y(t)), we know that we are looking through observed equations and can thus ask
    if var is bound, if it is, then sys₊y(t) is also bound. This can lead to an infinite recursion, so we maintain a stack of variables we have previously asked about to be able to break cycles
    =#
    u ∈ Set(stack) && return false # Cycle detected
    eqs = equations(sys)
    eqs = filter(eq -> has_var(eq, u), eqs) # Only look at equations that contain u
    # isout = isoutput(u)
    for eq in eqs
        vars = [get_variables(eq.rhs); get_variables(eq.lhs)]
        for var in vars
            var === u && continue
            if !same_or_inner_namespace(u, var)
                return true
            end
        end
    end
    # Look through observed equations as well
    oeqs = observed(sys)
    oeqs = filter(eq -> has_var(eq, u), oeqs) # Only look at equations that contain u
    for eq in oeqs
        vars = [get_variables(eq.rhs); get_variables(eq.lhs)]
        for var in vars
            var === u && continue
            if !same_or_inner_namespace(u, var)
                return true
            end
            if is_bound(sys, var, [stack; u]) && !inner_namespace(u, var) # The variable we are comparing to can not come from an inner namespace, binding only counts outwards
                return true
            end
        end
    end
    false
end

"""
    same_or_inner_namespace(u, var)

Determine whether or not `var` is in the same namespace as `u`, or a namespace internal to the namespace of `u`.
Example: `sys.u ~ sys.inner.u` will bind `sys.inner.u`, but `sys.u` remains an unbound, external signal. The namepsaced signal `sys.inner.u` lives in a namspace internal to `sys`.
"""
function same_or_inner_namespace(u, var)
    nu = get_namespace(u)
    nv = get_namespace(var)
    nu == nv ||           # namespaces are the same
        startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
        occursin('₊', var) && !occursin('₊', u) # or u is top level but var is internal
end

function inner_namespace(u, var)
    nu = get_namespace(u)
    nv = get_namespace(var)
    nu == nv && return false
    startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
        occursin('₊', var) && !occursin('₊', u) # or u is top level but var is internal
end

"""
    get_namespace(x)

Return the namespace of a variable as a string. If the variable is not namespaced, the string is empty.
"""
function get_namespace(x)
    sname = string(x)
    parts = split(sname, '₊')
    if length(parts) == 1
        return ""
    end
    join(parts[1:(end - 1)], '₊')
end

"""
    has_var(eq, x)

Determine whether or not an equation or expression contains variable `x`.
"""
function has_var(eq::Equation, x)
    has_var(eq.rhs, x) || has_var(eq.lhs, x)
end

has_var(ex, x) = x ∈ Set(get_variables(ex))

# Build control function

"""
    (f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); implicit_dae = false, ddvs = if implicit_dae

For a system `sys` that has unbound inputs (as determined by [`unbound_inputs`](@ref)), generate a function with additional input argument `in`
```
f_oop : (u,in,p,t)      -> rhs
f_ip  : (uout,u,in,p,t) -> nothing
```
The return values also include the remaining states and parameters, in the order they appear as arguments to `f`.

# Example
```
using ModelingToolkit: generate_control_function, varmap_to_vars, defaults
f, dvs, ps = generate_control_function(sys, expression=Val{false}, simplify=true)
p = varmap_to_vars(defaults(sys), ps)
x = varmap_to_vars(defaults(sys), dvs)
t = 0
f[1](x, inputs, p, t)
```
"""
function generate_control_function(sys::AbstractODESystem;
                                   implicit_dae = false,
                                   has_difference = false,
                                   simplify = true,
                                   kwargs...)
    ctrls = unbound_inputs(sys)
    if isempty(ctrls)
        error("No unbound inputs were found in system.")
    end

    # One can either connect unbound inputs to new parameters and allow structural_simplify, but then the unbound inputs appear as states :( .
    # One can also just remove them from the states and parameters for the purposes of code generation, but then structural_simplify fails :(
    # To have the best of both worlds, all unbound inputs must be converted to `@parameters` in which case structural_simplify handles them correctly :)
    sys = toparam(sys, ctrls)

    if simplify
        sys = structural_simplify(sys)
    end

    dvs = states(sys)
    ps = parameters(sys)

    dvs = setdiff(dvs, ctrls)
    ps = setdiff(ps, ctrls)
    inputs = map(x -> time_varying_as_func(value(x), sys), ctrls)

    eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
    check_operator_variables(eqs, Differential)
    # substitute x(t) by just x
    rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
           [eq.rhs for eq in eqs]

    # TODO: add an optional check on the ordering of observed equations
    u = map(x -> time_varying_as_func(value(x), sys), dvs)
    p = map(x -> time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)

    # pre = has_difference ? (ex -> ex) : get_postprocess_fbody(sys)

    args = (u, inputs, p, t)
    if implicit_dae
        ddvs = map(Differential(get_iv(sys)), dvs)
        args = (ddvs, args...)
    end
    pre, sol_states = get_substitutions_and_solved_states(sys)
    f = build_function(rhss, args...; postprocess_fbody = pre, states = sol_states,
                       kwargs...)
    f, dvs, ps
end

"""
    toparam(sys, ctrls::AbstractVector)

Transform all instances of `@varibales` in `ctrls` appearing as states and in equations of `sys` with similarly named `@parameters`. This allows [`structural_simplify`](@ref)(sys) in the presence unbound inputs.
"""
function toparam(sys, ctrls::AbstractVector)
    eqs = equations(sys)
    subs = Dict(ctrls .=> toparam.(ctrls))
    eqs = map(eqs) do eq
        substitute(eq.lhs, subs) ~ substitute(eq.rhs, subs)
    end
    ODESystem(eqs, name = nameof(sys))
end

function inputs_to_parameters!(state::TransformationState)
    @unpack structure, fullvars, sys = state
    @unpack var_to_diff, graph, solvable_graph = structure
    @assert solvable_graph === nothing

    inputs = BitSet()
    var_reidx = zeros(Int, length(fullvars))
    ninputs = 0
    nvar = 0
    new_parameters = []
    input_to_parameters = Dict()
    new_fullvars = []
    for (i, v) in enumerate(fullvars)
        if isinput(v) && !is_bound(sys, v)
            if var_to_diff[i] !== nothing
                error("Input $(fullvars[i]) is differentiated!")
            end
            push!(inputs, i)
            ninputs += 1
            var_reidx[i] = -1
            p = toparam(v)
            push!(new_parameters, p)
            input_to_parameters[v] = p
        else
            nvar += 1
            var_reidx[i] = nvar
            push!(new_fullvars, v)
        end
    end
    ninputs == 0 && return state

    nvars = ndsts(graph) - ninputs
    new_graph = BipartiteGraph(nsrcs(graph), nvars, Val(false))

    for ie in 1:nsrcs(graph)
        for iv in 𝑠neighbors(graph, ie)
            iv = var_reidx[iv]
            iv > 0 || continue
            add_edge!(new_graph, ie, iv)
        end
    end

    new_var_to_diff = DiffGraph(nvars, true)
    for (i, v) in enumerate(var_to_diff)
        new_i = var_reidx[i]
        (new_i < 1 || v === nothing) && continue
        new_v = var_reidx[v]
        @assert new_v > 0
        new_var_to_diff[new_i] = new_v
    end
    @set! structure.var_to_diff = new_var_to_diff
    @set! structure.graph = new_graph

    @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
    @set! sys.states = setdiff(states(sys), keys(input_to_parameters))
    @set! sys.ps = [parameters(sys); new_parameters]

    @set! state.sys = sys
    @set! state.fullvars = new_fullvars
    @set! state.structure = structure
end


text:
using Revise, PyPlot, PyCall
using Kwant
kwant = Kwant; pyplot = PyPlot

######### SPIN ORBIT ############

include("spin_orbit.jl")

syst = make_system()

# Check that the system looks as intended.
kwant.plot(syst)

# Finalize the system.
system = syst.finalized()

# We should see non-monotonic conductance steps.
plot_conductance(system, [0.01 * i - 0.3 for i in range(0,length=100)])
# CHANGE from oringal: energies not specified as keyword

######### QUANTUM WELL ############
include("quantum_well.jl")

syst = make_system()

# Check that the system looks as intended.
kwant.plot(syst)

# Finalize the system.
system = syst.finalized()

# We should see conductance steps.
plot_conductance(system, 0.2, [0.01 * i for i in 0:100-1]) ## note 1:100 not range(100)

################# AHARONOV-BOHM RING ################
include("ab_ring.jl")

syst = make_system()

# Check that the system looks as intended.
kwant.plot(syst)

# Finalize the system.
system = syst.finalized()

# We should see a conductance that is periodic with the flux quantum
plot_conductance(system, energy=0.15, fluxes=[0.01*i*3*2pi for i in range(0,length=100-1)])


text:
module IUP_IM

export
	# From libim.jl
	imFileOpen,
	imFileOpenAs,
	imFileNew,
	imFileClose,
	imFileHandle,
	imFileGetInfo,
	imFileSetInfo,
	imFileSetAttribute,
	imFileGetAttribute,
	imFileGetAttributeList,
	imFileGetPalette,
	imFileSetPalette,
	imFileReadImageInfo,
	imFileWriteImageInfo,
	imFileReadImageData,
	imFileWriteImageData,
	imFormatRegisterInternal,
	imFormatRemoveAll,
	imFormatList,
	imFormatInfo,
	imFormatInfoExtra,
	imFormatCompressions,
	imFormatCanWriteImage,
	imEncodeColor,
	imDecodeColor,
	imFileFormat,
	imImageInfo,
	imLoadRGB,
	imSaveRGB,
	imLoadMap,
	imSaveMap,
	imRGB2Map,
	imMap2RGB,
	imRGB2Gray,
	imMap2Gray,
	imResize,
	imStretch,
	imRegisterCallback,
	#
	# From im_image.jl		-- A submodule perhaps?
	#
	imImageCreate,
	imImageInit,
	imImageCreateBased,
	imImageDestroy,
	imImageAddAlpha,
	imImageSetAlpha,
	imImageRemoveAlpha,
	imImageReshape,
	imImageCopy,
	imImageCopyData,
	imImageCopyAttributes,
	imImageMergeAttributes,
	imImageCopyPlane,
	imImageDuplicate,
	imImageClone,
	imImageSetAttribute,
	imImageGetAttribute,
	imImageGetAttributeList,
	imImageClear,
	imImageIsBitmap,
	imImageSetPalette,
	imImageMatchSize,
	imImageMatchColor,
	imImageMatchDataType,
	imImageMatchColorSpace,
	imImageMatch,
	imImageSetMap,
	imImageSetBinary,
	imImageSetGray,
	imImageMakeBinary,
	imImageMakeGray,
	imFileLoadImage,
	imFileLoadImageFrame,
	imFileLoadBitmap,
	imFileLoadImageRegion,
	imFileLoadBitmapFrame,
	imFileSaveImage,
	imFileImageLoad,
	imFileImageLoadBitmap,
	imFileImageLoadRegion,
	imFileImageSave,
	#
	# From libim_h.jl
	#
	IM_INTERRUPTED,
	IM_ALL,
	IM_COUNTER_CB,
	IM_RESOLUTION_CB,
	IM_GIF_TRANSPARENT_COLOR_CB,
	IM_TIF_IMAGE_DESCRIPTION_CB,
	IM_BYTE,
	IM_SHORT,
	IM_USHORT,
	IM_INT,
	IM_FLOAT,
	IM_CFLOAT,
	IM_RGB,
	IM_MAP,
	IM_GRAY,
	IM_BINARY,
	IM_CMYK,
	IM_YCBCR,
	IM_LAB,
	IM_LUV,
	IM_XYZ,
	IM_ALPHA,
	IM_PACKED,
	IM_TOPDOWN,
	IM_ERR_NONE,
	IM_ERR_OPEN,
	IM_ERR_ACCESS,
	IM_ERR_FORMAT,
	IM_ERR_DATA,
	IM_ERR_COMPRESS,
	IM_ERR_MEM,
	IM_ERR_COUNTER,
	IM_ERR_READ,
	IM_ERR_WRITE,
	IM_ERR_TYPE,
	IM_ERR_COMP,
	IM_BMP,
	IM_PCX,
	IM_GIF,
	IM_TIF,
	IM_RAS,
	IM_SGI,
	IM_JPG,
	IM_LED,
	IM_TGA,
	IM_NONE,
	IM_DEFAULT,
	IM_COMPRESSED,
	IM_RES_NONE,
	IM_RES_DPI,
	IM_RES_DPC,
	# From iupim.jl
	IupLoadImage, IupSaveImage, IupGetNativeHandleImage, IupGetImageNativeHandle, IupImageFromImImage,
	# From im_image_h.jl
	imImage	# A composite type

include("libim_h.jl")
include("libim.jl")
include("im_image_h.jl")
include("im_image.jl")
include("iupim.jl")

end  # module