text:
export tikz,
       tikz_output,
       tikz_output_and_reset

take(out::IO) =
  let contents = String(take!(out))
    print(out, contents)
    contents
 end

tikz_e(out::IO, arg) =
  begin
    print(out, arg)
    println(out, ";")
  end

export use_wireframe
use_wireframe = Parameter(false)

tikz_draw(out::IO, filled=false) =
  print(out, filled && ! use_wireframe() ? "\\fill " : "\\draw ")

tikz_number(out::IO, x::Real) =
  isinteger(x) ?
    print(out, x) :
    (abs(x) < 0.0001 ?
      print(out, 0) :
      print(out, round(x*10000.0)/10000.0))

tikz_cm(out::IO, x::Real) = begin
  tikz_number(out, x)
  print(out, "cm")
end

tikz_2d_coord(out::IO, c::Loc) =
  let c = in_world(c)
    print(out, "(")
    tikz_number(out, c.x)
    print(out, ",")
    tikz_number(out, c.y)
    print(out, ")")
  end

tikz_3d_coord(out::IO, c::Loc) =
  let c = in_world(c)
    print(out, "(")
    tikz_number(out, c.x)
    print(out, ",")
    tikz_number(out, c.y)
    print(out, ",")
    tikz_number(out, c.z)
    print(out, ")")
  end
#
tikz_coord(out::IO, c::Loc) =
  let c = in_world(c)
    print(out, "(")
    tikz_number(out, c.x)
    print(out, ",")
    tikz_number(out, c.y)
    if !iszero(c.z)
      print(out, ",")
      tikz_number(out, c.z)
    end
    print(out, ")")
  end

tikz_pgfpoint(out::IO, c::Loc) =
  let c = in_world(c)
    if ! iszero(c.z)
      error("Can't handle 3D coords")
    end
    print(out, "\\pgfpoint{")
    tikz_cm(out, c.x)
    print(out, "}{")
    tikz_cm(out, c.y)
    print(out, "}")
  end

tikz_circle(out::IO, c::Loc, r::Real, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    tikz_coord(out, c)
    print(out, "circle(")
    tikz_cm(out, r)
    tikz_e(out, ")")
  end

tikz_point(out::IO, c::Loc) =
  tikz_circle(out, c, 0.01, true)

tikz_ellipse(out::IO, c::Loc, r0::Real, r1::Real, fi::Real, filled=false) =
  begin
    tikz_draw(out, filled)
    print(out, "[shift={")
    tikz_coord(out, c)
    print(out, "}]")
    print(out, "[rotate=")
    tikz_number(out, rad2deg(fi))
    print(out, "]")
    print(out, "(0,0)")
    print(out, "ellipse(")
    tikz_cm(out, r0)
    print(out, " and ")
    tikz_cm(out, r1)
    tikz_e(out, ")")
  end

tikz_arc(out::IO, c::Loc, r::Real, ai::Real, af::Real, filled=false) =
  begin
    tikz_draw(out, filled)
    if filled
      tikz_coord(out, c)
      print(out, "--")
    end
    tikz_coord(out, c+vpol(r, ai))
    print(out, "arc(")
    tikz_number(out, rad2deg(ai))
    print(out, ":")
    tikz_number(out, rad2deg(ai > af ? af+2*pi : af))
    print(out, ":")
    tikz_cm(out, r)
    print(out, ")")
    if filled
      tikz_e(out, "--cycle")    end
    println(out, ";")
  end

tikz_maybe_arc(out::IO, c::Loc, r::Real, ai::Real, da::Real, filled=false) =
  if iszero(r)
    tikz_point(out, c)
  elseif iszero(da)
    tikz_point(out, c + vpol(r, ai))
  elseif abs(da) >= 2*pi
    tikz_circle(out, c, r, filled)
  else
    let af = ai + da
      if af > ai
        tikz_arc(out, c, r, ai, af, filled)
      else
        tikz_arc(out, c, r, af, ai, filled)
      end
    end
  end

tikz_line(out::IO, pts::Locs, options::String="") =
  begin
    tikz_draw(out, false)
    print(out, options)
    tikz_coord(out, first(pts))
    for pt in Iterators.drop(pts, 1)
      print(out, "--")
      tikz_coord(out, pt)
    end
    println(out, ";")
  end

tikz_dimension(out::IO, p::Loc, q::Loc, text::AbstractString) =
  begin
    print(out, "\\dimline{")
    tikz_coord(out, p)
    print(out, "}{")
    tikz_coord(out, q)
    print(out, "}{;")
    print(out, text)
    println(out, "};")
  end

tikz_dim_line(out::IO, p::Loc, q::Loc, text::AbstractString, outside) =
  begin
    print(out, "\\draw[fill=black,latex-latex, very thin]")
    tikz_coord(out, p)
    print(out, "--")
    tikz_coord(out, q)
    println(out, "node[midway,auto=left]{$text};")
  end

tikz_closed_line(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    for pt in pts
      tikz_coord(out, pt)
      print(out, "--")
    end
    tikz_e(out, "cycle")
  end

tikz_closed_lines(out::IO, ptss, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    for pts in ptss
      for pt in pts
        tikz_coord(out, pt)
        print(out, "--")
      end
      print(out, "cycle ")
    end
    println(out, ";")
  end

tikz_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "plot [smooth,tension=1] coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

tikz_closed_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "plot [smooth cycle,tension=1] coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

# HACK we need to handle the starting and ending vectors
tikz_hobby_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    print(out, "\\begin{scope}")
    println(out, "[use Hobby shortcut, tension=0.1]")
    tikz_draw(out, false)
    tikz_coord(out, first(pts))
    for pt in Iterators.drop(pts, 1)
      print(out, "..")
      tikz_coord(out, pt)
    end
    println(out, ";")
    tikz_e(out, "\\end{scope}")
#=
    tikz_draw(out, filled)
    print(out, "[hobby, tension=0.1]")
    print(out, "plot coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")=#
  end

# HACK we need to handle the starting and ending vectors
tikz_hobby_closed_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "[closed hobby]")
    print(out, "plot coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

tikz_rectangle(out::IO, p::Loc, w::Real, h::Real, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    tikz_coord(out, p)
    print(out, "rectangle")
    tikz_coord(out, add_xy(p, w, h))
    println(out, ";")
  end

# Assuming default Arial font for AutoCAD
tikz_text(out::IO, txt, p::Loc, h::Real) =
  let (scale_x, scale_y) = (3.7*h, 3.7*h)
    tikz_draw(out)
    print(out, "[anchor=base west]")
    tikz_coord(out, p)
    print(out, "node[font=\\fontfamily{phv}\\selectfont,outer sep=0pt,inner sep=0pt")
    print(out, ",xscale=")
    tikz_number(out, scale_x)
    print(out, ",yscale=")
    tikz_number(out, scale_y)
    print(out, "]{")
    print(out, txt)
    tikz_e(out, "}")
  end

tikz_transform(out::IO, f::Function, c::Loc) =
  let m = c.cs.transform,
      t = in_world(c),
      a = m[1,1], b = m[2,1], c = m[1,2], d = m[2,2],
      tx = t.x, ty = t.y
    print(out, "\\begin{scope}")
    println(out, "[cm={$a, $b, $c, $d, ($tx, $ty)}]")
    f(out)
    tikz_e(out, "\\end{scope}")
  end

tikz_set_view(out::IO, view, options) =
  let v = view.camera - view.target,
      contents = String(take(out)),
      out = IOBuffer()
    print(out, "\\tdplotsetmaincoords{")
    tikz_number(out, rad2deg(sph_psi(v)))
    print(out, "}{")
    tikz_number(out, rad2deg(sph_phi(v))+90)
    println(out, "}")
    println(out, "\\begin{tikzpicture}[tdplot_main_coords$(use_wireframe() ? "" : ",fill=gray")$(options=="" ? "" : ",")$options]") #)opacity=0.2")]")
    print(out, contents)
    println(out, "\\end{tikzpicture}")
    String(take!(out))
  end

#=
tikz_set_view(out::IO, camera::Loc, target::Loc, lens::Real) =
  let v = camera - target,
      contents = String(take(out)),
      out = IOBuffer()
    println(out, raw"\begin{tikzpicture}")
    #println(out, "\\begin{axis}[view={$(rad2deg(sph_phi(v))+90)}{$(rad2deg(sph_psi(v)))},axis equal image,hide axis,colormap/blackwhite]")
    println(out, "\\begin{axis}[axis equal image,hide axis,colormap/blackwhite]")
    print(out, contents)
    println(out, raw"\end{axis}")
    println(out, raw"\end{tikzpicture}")
    String(take!(out))
  end
=#

tikz_set_view_top(out::IO, options) =
  let contents = String(take(out)),
      out = IOBuffer()
    println(out, "\\begin{tikzpicture}[$options]")
    print(out, contents)
    println(out, "\\end{tikzpicture}")
    String(take!(out))
  end

#=
  \begin{tikzpicture}
  \begin{axis}[view={135}{45},axis equal image,scale=4,hide axis
=#
#

abstract type TikZKey end
const TikZId = Nothing
const TikZIds = Vector{TikZId}
const TikZRef = GenericRef{TikZKey, TikZId}
const TikZRefs = Vector{TikZRef}
const TikZEmptyRef = EmptyRef{TikZKey, TikZId}
const TikZUniversalRef = UniversalRef{TikZKey, TikZId}
const TikZNativeRef = NativeRef{TikZKey, TikZId}
const TikZUnionRef = UnionRef{TikZKey, TikZId}
const TikZSubtractionRef = SubtractionRef{TikZKey, TikZId}
const TikZ = IOBufferBackend{TikZKey, TikZId}

KhepriBase.void_ref(b::TikZ) = TikZNativeRef(nothing)

const tikz = TikZ(view=View(xyz(10,10,10), xyz(0,0,0), 0, 0))

KhepriBase.backend_name(b::TikZ) = "TikZ"

tikz_output(options="") =
  let b = tikz
    b.view.lens == 0 ?
      tikz_set_view_top(connection(b), options) :
      tikz_set_view(connection(b), b.view, options)
  end

tikz_output_and_reset(options="") =
  let b = tikz,
      out = tikz_output(options)
    b_delete_all_refs(b)
    b.view.lens = 0 # Reset the lens
    out
  end

withTikZXForm(f, out, c) =
  if is_world_cs(c.cs)
    f(out, c)
  else
    tikz_transform(out,
      out -> f(out, u0(world_cs)),
      c)
   end

KhepriBase.b_delete_all_refs(b::TikZ) =
  truncate(connection(b), 0)

KhepriBase.b_point(b::TikZ, p) =
  tikz_pgfpoint(connection(b), p)

KhepriBase.b_line(b::TikZ, ps, mat) =
  tikz_line(connection(b), ps)

KhepriBase.b_polygon(b::TikZ, ps, mat) =
  tikz_closed_line(connection(b), ps)

KhepriBase.b_spline(b::TikZ, ps, v0, v1, interpolator, mat) =
  if (v0 == false) && (v1 == false)
    #tikz_hobby_spline(connection(b), ps, false)
    tikz_spline(connection(b), ps, false)
  elseif (v0 != false) && (v1 != false)
    TikZInterpSpline(connection(b), ps, v0, v1)
  else
    TikZInterpSpline(connection(b),
                     ps,
                     v0 == false ? ps[2] - ps[1] : v0,
                     v1 == false ? ps[end-1] - ps[end] : v1)
  end

KhepriBase.b_closed_spline(b::TikZ, ps, mat) =
  tikz_hobby_closed_spline(connection(b), ps)

KhepriBase.b_circle(b::TikZ, c, r, mat) =
  withTikZXForm(connection(b), c) do out, cc
    tikz_circle(out, cc, r)
  end

KhepriBase.b_arc(b::TikZ, c, r, α, Δα, mat) =
  withTikZXForm(connection(b), c) do out, cc
    tikz_maybe_arc(out, cc, r, α, Δα, false)
  end

KhepriBase.b_rectangle(b::TikZ, c, dx, dy, mat) =
  withTikZXForm(connection(b), c) do out, cc
    tikz_rectangle(out, cc, dx, dy)
  end

# KhepriBase.b_trig(b::TikZ, p1, p2, p3, mat) =
#   tikz_closed_line(connection(b), [p1, p2, p3], true)

#
# KhepriBase.b_trig(b::TikZ, p1, p2, p3, mat) =
#   let io = connection(b)
#     println(io, raw"\addplot3[patch,table/row sep=\\,patch table={")
#     println(io, "0 1 2 \\")
#     println(io, raw")}] table [row sep=\\] {")
#      x y z c\\
#      0 1 0 0\\
#      0 0 -1 0\\
#      -1 0 0 0\\
#      0 0 1 0\\
#      1 0 0 0\\
#     };
#
# KhepriBase.b_quad(b::TikZ, p1, p2, p3, p4, mat) =
#   tikz_closed_line(connection(b), [p1, p2, p3, p4], true)

KhepriBase.b_trig(b::TikZ, p1, p2, p3, mat) =
  let io = connection(b)
    print(io, raw"\addplot3[patch,shader=interp] coordinates {")
    tikz_3d_coord(io, p1)
    tikz_3d_coord(io, p2)
    tikz_3d_coord(io, p3)
    println(io, "};")
  end

KhepriBase.b_quad(b::TikZ, p1, p2, p3, p4, mat) =
  invoke(b_quad, Tuple{Backend, Any, Any, Any, Any, Any}, b, p1, p2, p3, p4, mat)
  # let io = connection(b)
  #   print(io, raw"\addplot3[patch,shader=interp] coordinates {")
  #   tikz_3d_coord(io, p1)
  #   tikz_3d_coord(io, p2)
  #   tikz_3d_coord(io, p3)
  #   tikz_3d_coord(io, p4)
  #   println(io, "};")
  # end

KhepriBase.b_surface_polygon(b::TikZ, ps, mat) =
  tikz_closed_line(connection(b), ps, true)
  #=

KhepriBase.b_surface_polygon_with_holes(b::TikZ, ps, qss, mat) =
  tikz_closed_lines(connection(b), [ps, qss...], true)

KhepriBase.b_surface_circle(b::TikZ, c, r, mat) =
  withTikZXForm(connection(b), c) do out, cc
    tikz_circle(out, cc, r, true)
  end

KhepriBase.b_surface_arc(b::TikZ, c, r, α, Δα, mat) =
  withTikZXForm(connection(b), c) do out, cc
    tikz_maybe_arc(out, cc, r, α, Δα, true)
  end
=#
# realize(b::TikZ, s::Ellipse) =
#   withTikZXForm(connection(b), s.center) do out, c
#     tikz_ellipse(out, c, s.radius_x, s.radius_y, 0, false)
#   end
#
# realize(b::TikZ, s::SurfaceEllipse) =
#   withTikZXForm(connection(b), s.center) do out, c
#     tikz_ellipse(out, c, s.radius_x, s.radius_y, 0, true)
#   end
#
# realize(b::TikZ, s::EllipticArc) =
#   error("Finish this")

#realize(b::TikZ, s::SurfaceElliptic_Arc) = TikZCircle(connection(b),

# KhepriBase.b_surface_rectangle(b::TikZ, c, dx, dy, mat) =
#   withTikZXForm(connection(b), c) do out, cc
#     tikz_rectangle(out, cc, dx, dy, true)
#   end

KhepriBase.b_text(b::TikZ, str, p, size, mat) =
  #invoke(b_text, Tuple{Backend, Any, Any, Any, Any}, b, str, p, size, mat)
  withTikZXForm(connection(b), p) do out, c
    tikz_text(out, str, c, size)
  end

KhepriBase.b_dim_line(b::TikZ, p, q, tv, str, size, outside, mat) =
  #invoke(b_dim_line, Tuple{Backend, Any, Any, Any, Any, Any, Any, Any}, b, p, q, tv, str, size, outside, mat)
  tikz_dim_line(connection(b), p, q, str, outside)

KhepriBase.b_ext_line(b::TikZ, p, q, mat) =
  tikz_line(connection(b), [p, q], "[very thin]")

# realize(b::TikZ, s::SurfaceGrid) =
#   invoke(realize, Tuple{Backend, SurfaceGrid}, b, s)
  # let n = size(s.points,1),
  #     m = size(s.points,2)
  #   for i in 1:n
  #     tikz_hobby_spline(connection(b), s.points[i,:], false)
  #   end
  #   for j in 1:m
  #     tikz_hobby_spline(connection(b), s.points[:,j], false)
  #   end
  # end

# For extra non-portable stuff
export add_tikz
add_tikz(str) =
  println(connection(tikz), str)
###
# Dimensioning

#backend_dimension(b::TikZ, pa, pb, sep, scale, style) =
#  tikz_dimension(connection(b), pa, pb, )

#=openpdf()
      pdfname = makepdf(latex)
      if Sys.iswindows()
          command = `cmd /K start \"\" $pdfname`
          run(command)
      else
          run(`open $pdfname`)
      end
      pdfname
  end
=#
const miktex_folder = Parameter("C:/Users/aml/AppData/Local/Programs/MiKTeX/miktex/bin/x64/")
miktex_cmd(cmd::AbstractString="pdflatex") = miktex_folder() * cmd

process_tikz(path) =
  let contents = tikz_output(),
      path = path_replace_suffix(path, ".tex"),
      pdfpath = path_replace_suffix(path, ".pdf")
    rm(pdfpath, force=true)
    open(path, "w") do out
      println(out, raw"\documentclass{standalone}")
      println(out, raw"\usepackage{tikz}")
      println(out, raw"\usetikzlibrary{patterns}")
      println(out, raw"\usetikzlibrary{calc,fadings,decorations.pathreplacing}")
      println(out, raw"\usetikzlibrary{shapes,fit}")
      println(out, raw"\usetikzlibrary{hobby}")
      #println(out, raw"\usepackage{pgfplots}")
      #println(out, raw"\pgfplotsset{compat=1.17}")
      println(out, raw"\usepackage{tikz-3dplot}")
      println(out, raw"\begin{document}")
      println(out, contents)
      println(out, raw"\end{document}")
    end
    # cd(dirname(path)) do
    #   output = read(`$(miktex_cmd()) -shell-escape -halt-on-error $(path)`, String)
    #   occursin("Error:", output) && println(output)
    # end
    cd(dirname(path)) do
      run(`$(miktex_cmd("texify")) --pdf --engine=luatex --run-viewer $(path)`, wait=true)
      #output = read(`$(miktex_cmd("texify")) --run-viewer $(path)`, String)
      #occursin("Error:", output) && println(output)
    end
  end

export visualize_tikz
visualize_tikz(name="Test") =
  with(render_kind_dir, "TikZ",
       render_ext, ".tex") do
    process_tikz(prepare_for_saving_file(render_pathname(name)))
    @info "Tex file: $(render_pathname(name))"
  end

KhepriBase.b_render_pathname(b::TikZ, name::String) =
  path_replace_suffix(render_pathname(name), ".pdf")

KhepriBase.b_render_view(b::TikZ, path) =
  process_tikz(path)


text:
module DocumenterMarkdown

using Documenter

const Markdown = Documenter.Writers.MarkdownWriter.Markdown
export Markdown

function __init__()
    if !isdefined(Documenter.Writers, :enable_backend)
        @warn """Incompatible Documenter version.

        Documenter.Writers is missing the enable_backend() function.
        """
        return
    end
    Documenter.Writers.enable_backend(:markdown)
    nothing
end

end # module

ext: jl


text:
push!(LOAD_PATH,"../src/")

ENV["GKSwstype"]="100"

using Documenter
using Plots
using Literate
using WaterWaves1D

# generate example

DOC_OUTPUT  = joinpath(@__DIR__, "src", "examples")
NB_OUTPUT   = joinpath(@__DIR__, "src", "notebooks")

# Literate.markdown(joinpath(@__DIR__, "..",  "examples","two_problems.jl"), DOC_OUTPUT, documenter=true)
# Literate.notebook(joinpath(@__DIR__, "..",  "examples","two_problems.jl"), NB_OUTPUT, execute=false)
# Literate.markdown(joinpath(@__DIR__, "..",  "examples","animation.jl"), DOC_OUTPUT, documenter=true)
# Literate.notebook(joinpath(@__DIR__, "..",  "examples","animation.jl"), NB_OUTPUT, execute=false)

makedocs(modules=[WaterWaves1D],
         doctest = false,
         authors = "Vincent Duchene",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         sitename = "WaterWaves1D.jl",
         pages = [ "Quickstart" => "quickstart.md",
                   "Documentation" => "index.md",
                   "Code basics"   => "basics.md",
                   "Contents"      => "contents.md"])

#                  "Animation"     => "examples/animation.md",
#                 "Example"       => "examples/two_problems.md",

deploydocs(
    repo   = "github.com/WaterWavesModels/WaterWaves1D.jl.git"
 )


text:
# Stubs - Can be used as references
struct Account <: QBObject end
struct ItemBasedExpenseLineDetail; end
struct Employee <: QBObject end
struct Vendor <: QBObject end
struct Customer <: QBObject end
struct Item <: QBObject end

struct Company <: QBObject
    Id::Maybe{Int}
end

from_json(::Type{ItemBasedExpenseLineDetail}, data) = ItemBasedExpenseLineDetail()

# Class
struct Class <: QBObject
    Id::Maybe{Int64}
    SyncToken::Maybe{Int64}
    Parent::Optional{QboRef{Class}}
    Name::Maybe{String}
    Active::Maybe{Bool}
end
@eval @make_default_methods $Class

# PaymentMethod
mutable struct PaymentMethod
    Id::Maybe{Int}
    SyncToken::Maybe{Int}
    Name::Maybe{String}
    Type::Maybe{String}
end
@eval @make_default_methods $PaymentMethod

### Purchases

struct AccountBasedExpenseLineDetail
    Account::QboRef{Account}
    Class::Optional{QboRef{Class}}
    BillableStatus::String
end
@eval @make_default_methods $AccountBasedExpenseLineDetail

const DetailType = Union{ItemBasedExpenseLineDetail,
                         AccountBasedExpenseLineDetail}

struct Line
    Id::Maybe{Int}
    Amount::Maybe{Decimal}
    Detail::DetailType
end

Line(Amount::Decimal, Detail::DetailType) = Line(missing, Amount, Detail)

@eval from_json(::Type{Line}, data) = @from_json $Line $(Dict(
    :Detail => quote
        if data["DetailType"] == "AccountBasedExpenseLineDetail"
            from_json(AccountBasedExpenseLineDetail,
                data["AccountBasedExpenseLineDetail"])
        elseif data["DetailType"] == "ItemBasedExpenseLineDetail"
            from_json(ItemBasedExpenseLineDetail,
                data["ItemBasedExpenseLineDetail"])
        end
    end
))

@eval to_json(data::Line) = @to_json $Line $(Dict(
    :Detail => quote
        if isa(data.Detail, AccountBasedExpenseLineDetail)
            ret["DetailType"] = "AccountBasedExpenseLineDetail"
            ret["AccountBasedExpenseLineDetail"] = to_json(data.Detail)
        elseif isa(data.Detail, ItemBasedExpenseLineDetail)
            ret["DetailType"] = "ItemBasedExpenseLineDetail"
            ret["ItemBasedExpenseLineDetail"] = to_json(data.Detail)
        end
    end
    ))

mutable struct Purchase <: QBObject
    Id::Maybe{Int}
    SyncToken::Maybe{Int}
    Account::Maybe{QboRef{Account}}
    PaymentMethod::Optional{QboRef{PaymentMethod}}
    PaymentType::Maybe{String}
    Entity::Maybe{QboRef{<:QBObject}}
    TxnDate::Maybe{Date}
    Lines::Vector{Line}
end

@eval @make_default_methods $Purchase $(Dict(
    :Lines => "Line"
))

# Type Mapping
function key_to_type(key::String)
    Dict{String, Type}(
        "Purchase" => Purchase,
        "Class" => Class,
        "Employee" => Employee,
        "PaymentMethod" => PaymentMethod
    )[key]
end


text:
include("util.jl")
include("graph.jl")
include("opinion.jl")

function main(dataset="congress")

    delim = DATASETS[dataset]
    DATASET = dataset
    FILEPATH =  "./data/raw/$(DATASET).edgelist"
    PLOT_PATH = "./plots/traces/$(DATASET).png"
    (W, U, n) = read_edgelist(FILEPATH, delim)
    normalize_edgelist!(W)
    normalize_edgelist!(U)
    y = generate_opinions(n)
    α, β, γ = generate_factors(n)
    M = (β .* W) - (γ .* U)
    #M = Matrix(M)
    #Q = inv(I - M)

    times  = 20

    trace = fill(WEIGHT(0.0), n, times)
    trace[:,1] = y

    println("Computing traces")
    for i = 2:times
        z = trace[:,i - 1]
        next_opinion = (α .* y) + (M * z)
        trace[:,i] = next_opinion
    end

    println("Ploting traces")
    p = plot(legend = false)
    for i = 1:n
        x = 1:times
        to_plot = trace[i, :]
        plot!(p, x, to_plot, legend = false)
    end

    println("Saving plots of traces")
    savefig(PLOT_PATH)

end

main("twitterreferendum")

max_stars_repo_name: InzamamRahaman/SignedDegroot

max_issues_repo_name: InzamamRahaman/SignedDegroot

max_forks_repo_name: InzamamRahaman/SignedDegroot


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################

using Coverage, Base.Test

# test our filename matching. These aren't exported functions but it's probably
# a good idea to have explicit tests for them, as they're used to match files
# that get deleted

@test Coverage.iscovfile("test.jl.cov")
@test Coverage.iscovfile("test.jl.2934.cov")
@test Coverage.iscovfile("/home/somebody/test.jl.2934.cov")
@test !Coverage.iscovfile("test.ji.2934.cov")
@test !Coverage.iscovfile("test.jl.2934.cove")
@test !Coverage.iscovfile("test.jicov")
@test !Coverage.iscovfile("test.c.cov")
@test Coverage.iscovfile("test.jl.cov", "test.jl")
@test !Coverage.iscovfile("test.jl.cov", "other.jl")
@test Coverage.iscovfile("test.jl.8392.cov", "test.jl")
@test Coverage.iscovfile("/somedir/test.jl.8392.cov", "/somedir/test.jl")
@test !Coverage.iscovfile("/otherdir/test.jl.cov", "/somedir/test.jl")

cd(Pkg.dir("Coverage")) do
    datadir = joinpath("test", "data")
    # Process a saved set of coverage data...
    r = process_file(joinpath(datadir,"Coverage.jl"))
    # ... and memory data
    analyze_malloc(datadir)
    lcov = IOBuffer()
    # we only have a single file, but we want to test on the Vector of file results
    LCOV.write(lcov, FileCoverage[r])
    open(joinpath(datadir, "expected.info")) do f
        @test takebuf_string(lcov) == readall(f)
    end

    # Test a file from scratch
    srcname = joinpath("test", "data","testparser.jl")
    covname = srcname*".cov"
    # clean out any previous coverage files. Don't use clean_folder because we
    # need to preserve the pre-baked coverage file Coverage.jl.cov
    clean_file(srcname)
    cmdstr = "include(\"$srcname\"); using Base.Test; @test f2(2) == 4"
    run(`julia --code-coverage=user -e $cmdstr`)
    r = process_file(srcname, datadir)
    # The next one is the correct one, but julia & JuliaParser don't insert a line number after the 1-line @doc -> test
    # See https://github.com/JuliaLang/julia/issues/9663 (when this is fixed, can uncomment the next line on julia 0.4)
    target = Union{Int64,Void}[nothing, nothing, nothing, nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing]
    #target = Union{Int64,Void}[nothing, nothing, nothing, nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing]
    @test r.coverage[1:length(target)] == target

    covtarget = (sum(x->x != nothing && x > 0, target), sum(x->x != nothing, target))
    @test get_summary(r) == covtarget
    @test get_summary(process_folder(datadir)) != covtarget

    #json_data = Codecov.build_json_data(Codecov.process_folder("data"))
    #@test typeof(json_data["coverage"]["data/Coverage.jl"]) == Array{Union{Int64,Void},1}

end


text:
@testset "SPIR-V" begin

include("definitions/spirv.jl")

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

@testset "IR" begin

@testset "kernel functions" begin
@testset "calling convention" begin
    kernel() = return

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{}; dump_module=true))
    @test !occursin("spir_kernel", ir)

    ir = sprint(io->spirv_code_llvm(io, kernel, Tuple{};
                                    dump_module=true, kernel=true))
    @test occursin("spir_kernel", ir)
end
end

end

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

end


text:
"""
    inpaint(A)

Inpaints `missing` values if `A` contains some.

Inspired by MATLAB's `inpaint_nans`'s (by John d'Errico: [link](https://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint_nans)).

# Example
```jldoctest
julia> A = convert(Array{Union{Float64, Missing},2}, (1:3)*(1:4)') ; A[1:2, 1:2] .= missing ; A
3×4 Array{Union{Missing, Float64},2}:
  missing   missing  3.0   4.0
  missing   missing  6.0   8.0
 3.0       6.0       9.0  12.0

julia> inpaint(A)
3×4 Array{Union{Missing, Float64},2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```
"""
inpaint(A::AbstractArray{Union{Missing, T}}; method=1, cycledims=Int64[]) where {T<:AbstractFloat} = inpaint(ismissing, A, method=method, cycledims=cycledims)

"""
    inpaint(A)

Inpaints `NaN` values if the elements of `A` are all float (i.e., `eltype(A) <: AbstractFloat`).

# Example

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, 1:2] .= NaN ; A
3×4 Array{Float64,2}:
 NaN    NaN    3.0   4.0
 NaN    NaN    6.0   8.0
   3.0    6.0  9.0  12.0

julia> inpaint(A)
3×4 Array{Float64,2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```
"""
inpaint(A::AbstractArray{T}; method=1, cycledims=Int64[]) where {T<:AbstractFloat} = inpaint(isnan, A, method=method, cycledims=cycledims)

"""
    inpaint(A, missing)

Inpaints `missing` values.
Should be the same as `inpaint(A)`.
"""
inpaint(A::AbstractArray, ::Missing; method=1, cycledims=Int64[]) = inpaint(ismissing, A, method=method, cycledims=cycledims)

"""
    inpaint(A, value_to_fill)

Inpaints `value_to_fill` values (can be `NaN`).
Useful in the case data was generated by filling missing values with an
arbitrary chosen one.

# Examples

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, 1:2] .= 999 ; A
3×4 Array{Float64,2}:
 999.0  999.0  3.0   4.0
 999.0  999.0  6.0   8.0
   3.0    6.0  9.0  12.0

julia> inpaint(A, 999)
3×4 Array{Float64,2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, 1:2] .= NaN ; A
3×4 Array{Float64,2}:
 NaN    NaN    3.0   4.0
 NaN    NaN    6.0   8.0
   3.0    6.0  9.0  12.0

julia> inpaint(A, NaN)
3×4 Array{Float64,2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```
"""
function inpaint(A::AbstractArray{T}, value_to_fill; method=1, cycledims=Int64[]) where T<:AbstractFloat
    return if isnan(value_to_fill)
        inpaint(isnan, A, method=method, cycledims=cycledims)
    else
        inpaint(x -> x == value_to_fill, A, method=method, cycledims=cycledims)
    end
end

"""
    inpaint(f, A)

Inpaints values of `A` for which `f(A) == true`.

# Examples

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, 1:2] .= 999 ; A
3×4 Array{Float64,2}:
 999.0  999.0  3.0   4.0
 999.0  999.0  6.0   8.0
   3.0    6.0  9.0  12.0

julia> inpaint(x -> x == 999, A)
3×4 Array{Float64,2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, [1, end]] .= NaN ; A
3×4 Array{Float64,2}:
 NaN    2.0  3.0  NaN
 NaN    4.0  6.0  NaN
   3.0  6.0  9.0   12.0

julia> inpaint(isnan, A)
3×4 Array{Float64,2}:
 1.0  2.0  3.0   4.0
 2.0  4.0  6.0   8.0
 3.0  6.0  9.0  12.0
```

    inpaint(f, A; method=1, cycledims=[])

With optional arguments, you can chose the inpainting method and if dimensions are cyclic.

```jldoctest
julia> A = float((1:3)*(1:4)') ; A[1:2, [1, end]] .= NaN ; A
3×4 Array{Float64,2}:
 NaN    2.0  3.0  NaN
 NaN    4.0  6.0  NaN
   3.0  6.0  9.0   12.0

julia> inpaint(A, cycledims=[2])
3×4 Array{Float64,2}:
 2.16475  2.0  3.0   2.83525
 3.76245  4.0  6.0   6.23755
 3.0      6.0  9.0  12.0
```
"""
function inpaint(f, A; method=1, cycledims=Int64[])
    return @match method begin
        1 => inpaint_method1(f, A, cycledims=cycledims) # Default
        _ => error("Method $method not available yet. Suggest it on Github!")
    end
end

"""
    list_neighbors(A, idx, neighbors)

Lists all the neighbors of the indices in `idx` in Array `A`.
Neighbors already in `idx` are accepted.
The argument `neighnors` contains a list of the neighbors about the
origin coordinate `(0, 0, ...)`.
In other words, it is a `Vector` of `CartesianIndex` such that
the direct neighbors of index `i` are given by `i + n for n in neighbors`.
Inspired by MATLAB's `inpaint_nans`'s `identify_neighbors` (by John d'Errico).
See https://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint_nans.
"""
function list_neighbors(R, idx, neighbors)
    out = [i + n for i in idx for n in neighbors if i + n ∈ R]
    out = sort(unique([out; idx]))
end
list_neighbors(R, idx::Vector{Int64}, neighbors) = list_neighbors(R, CartesianIndex.(idx), neighbors)

"""
    inpaint_method1(f, A::Array, cycledims=Int64[])

Inpaints values in `A` that `f` gives `true` on by solving a simple diffusion PDE.
Default method for `inpaint`.
The partial differential equation (PDE) is defined by the standard Laplacian, `Δ = ∇^2`.
Inspired by MATLAB's `inpaint_nans`'s method `0` for matrices (by John d'Errico).
See https://www.mathworks.com/matlabcentral/fileexchange/4551-inpaint_nans.
The discrete 1D stencil used for `Δ` looks like
```
┌───┐ ┌───┐ ┌───┐
│ 1 ├─┤-2 ├─┤ 1 │
└───┘ └───┘ └───┘
```
and is applied where it fits at the borders in all dimensions.
However, the user can supply a list of dimensions that should be considered cyclic.
In this case, the sentil will be used also at the borders and "jump" to the other side.
This is particularly useful for, e.g., world maps with longitudes spanning the entire globe.
"""
function inpaint_method1(f, A::Array; cycledims=Int64[])
    I_unknown = findall(@. f(A))
    I_known = findall(@. !f(A))

    # direct neighbors
    Ns = [CartesianIndex(ntuple(x -> x == d ? 1 : 0, ndims(A))) for d in 1:ndims(A)]
    neighbors = [Ns; -Ns]
    # Add neighbors on opposite border if dimensions are cyclic
    N′s = [CartesianIndex(ntuple(x -> x == d ? size(A,d) - 1 : 0, ndims(A))) for d in 1:ndims(A)]
    for d in cycledims # Add neighbors in cyclic case
        neighbors = push!(neighbors, N′s[d], -N′s[d])
    end

    # "working" list: the indices of all unkowns and their neighbors
    R = CartesianIndices(size(A))
    Iwork = list_neighbors(R, I_unknown, neighbors)
    iwork = LinearIndices(size(A))[Iwork]

    # Preallocate the Laplacian (not the fastest way but easier-to-read code)
    nA = length(A)
    Δ = sparse([], [], Vector{Float64}(), nA, nA)

    # Fill the Laplacian within the borders
    for d in 1:ndims(A)
        N = Ns[d]
        # R′ is the range of indices without the borders in dimension `d`
        R′ = CartesianIndices(size(A) .- 2 .* N.I)
        R′ = [r + N for r in R′]
        # Convert to linear indices to build the Laplacian
        u = vec(LinearIndices(size(A))[R′])
        n = LinearIndices(size(A))[first(R) + N] - LinearIndices(size(A))[first(R)]
        # Build the Laplacian (not the fastest way but easier-to-read code)
        Δ += sparse(u, u     , -2.0, nA, nA)
        Δ += sparse(u, u .- n, +1.0, nA, nA)
        Δ += sparse(u, u .+ n, +1.0, nA, nA)
    end

    # Add Laplacian along border if cyclic along dimension `d`
    for d in cycledims
        N = Ns[d]   # usual neighbor
        N′ = N′s[d] # neighbor on the opposite border
        R′₊ = [r for r in R if r + N ∉ R] # one border
        R′₋ = [r for r in R if r - N ∉ R] # the opposite border
        # Convert to linear indices to build the Laplacian
        n = LinearIndices(size(A))[first(R) + N] - LinearIndices(size(A))[first(R)]
        n′ = LinearIndices(size(A))[first(R) + N′] - LinearIndices(size(A))[first(R)]
        u₊ = LinearIndices(size(A))[R′₊]
        u₋ = LinearIndices(size(A))[R′₋]
        # Add to the Laplacian
        Δ += sparse(u₊, u₊      , -2.0, nA, nA)
        Δ += sparse(u₊, u₊ .- n , +1.0, nA, nA)
        Δ += sparse(u₊, u₊ .- n′, +1.0, nA, nA)
        Δ += sparse(u₋, u₋      , -2.0, nA, nA)
        Δ += sparse(u₋, u₋ .+ n , +1.0, nA, nA)
        Δ += sparse(u₋, u₋ .+ n′, +1.0, nA, nA)
    end

    # Use linear indices to access the sparse Laplacian
    i_unknown = LinearIndices(size(A))[I_unknown]
    i_known = LinearIndices(size(A))[I_known]

    # Place the knowns to right hand side
    rhs = -Δ[iwork, i_known] * A[i_known]

    # Solve for the unknowns
    B = zeros(size(A))
    B[i_known] .= A[i_known]
    B[i_unknown] .= Δ[iwork, i_unknown] \ rhs
    return B
end

export inpaint


text:
"""
cohorts as 3 matrices: year_b - initial year, year_e - last year,  x - abundance, y - catch, w - weight (obsolete)
"""
mutable struct  cohorts_t
     year_b
    year_e
    x
     y
    w

end

"""
detail cohort representation
tmax - number of seasons
jmax - number of ages
yb - physical year assciated with cohort (when age=1)
jb - youngest age
l - length
jt - trusted age
xv,yv,wv - vectors with abundunce, catch, weight
"""
mutable struct cohort_v
    tmax
    jmax
    yb
    jb
    l
    jt
     xv::Array{Float64}
     yv::Array{Float64}
     wv::Array{Float64}
 end

"""
detail  representation of given cohorts
tmax - number of seasons
jmax - number of ages
cs -  vector of cohort of mutable struct cohort_v
"""
mutable struct cohorts_v
       tmax
    jmax
    cs::Array{cohort_v}
 end

"""
detail  representation of given cohorts
tmax - number of seasons
jmax - number of ages
cs -  vector of cohort of mutable struct cohort_v
"""

mutable struct cohorts_v1
    year_b
    tmax
    jmax
    cs::Array{cohort_v}
 end

"""
detail cohort representation
no- the number
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
jbeg - initial age
len - length
jtrust - trusted age
xv,yv,wv - vectors with abundunce, catch, weight
"""
mutable struct cohort_v1
    no
    yb
    tmax
    jmax
    jbeg
    len
    jtrust
     xv::Array{Float64}
     yv::Array{Float64}
     wv::Array{Float64}
 end

"""
cohorts as 3 tables (=cohorts_m)
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
x,y,w- tables with abundunce, catch, weight
"""
mutable struct  cohorts_t1
     year_b
   tmax
    jmax
    x
     y
    w
end

"""cohorts as cs: cohort_v1?
year_b - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
"""
mutable struct cohorts_v2
    year_b
    tmax
    jmax
    cs
end

"""
cohorts related matrix
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
mtx: x,y,w,m,f, etc.
"""
mutable struct csmtx
    year_b
    tmax
    jmax
    mtx
end

"""
detail mtx representation of a cohort as sequence of vecs
no- the number
yb - physical year associated with cohort (when age=1)
jbeg - youngest age
len - length
jtrust - trusted age
vec
"""
mutable struct cvec
    no
    yb
    jbeg
    len
    jtrust
    vec::Array{Float64}

 end

"""
cohort as vecs
"""
mutable struct csvecs
    year_b
    tmax
    jmax
    vecs
end

"""
cohorts_t1->cohorts_v2

cohorts_t1:
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
x,y,w -  matrices  with abundunce, catch, weight

cohorts_v2:
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
cs - cohorts of mutable struct cohort_v1

cohort_v1:
no- the number
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
jbeg - initial age
len - length
jtrust - trusted age
xv,yv,wv- vectors with abundunce, catch, weight
tmax+jmax-1 -- number of related cohorts
"""
function cst2csv2(cs)
   x=cs.x
    y=cs.y
    w=cs.w
    tmax=cs.tmax
    jmax=cs.jmax
    year1=cs.year_b
    year0=year1-jmax+1
                #println(year0)
                #println(year1)

#    csv1=Array{cohort_v1}(tmax+jmax-1)
csv1=Array{cohort_v1}(undef,tmax+jmax-1)

    #csv1=Array{Any}(undef,tmax+jmax-1)

    for year in year0:year1-1
        i=year-year0+1
                #println(year)
                #println(i)

        csv1[i]=cohort_v1(i,year,tmax,jmax,
        jmax-i+1,i,jmax-i+1,
        map(j->x[j,jmax-i+j],1:i),
        map(j->y[j,jmax-i+j],1:i),
        map(j->w[j,jmax-i+j],1:i))
                #println(csv1[i])
    end

    year2=year1+tmax-jmax
                #println(year2)

    if tmax>=jmax
        for year in year1:year2
            i=year-year1+1
               #println(year)
               #println(i)
               #println(i-year0+year1)
            csv1[i-year0+year1]=cohort_v1(i-year0+year1,year,tmax,jmax,
            1,jmax,1,
            map(j->x[i+j-1,j],1:jmax),
            map(j->y[i+j-1,j],1:jmax),
            map(j->w[i+j-1,j],1:jmax))
               #println( csv1[i-year0+year1])
        end
    end

    year3=year2+jmax-1
                #println(year3)

    for year in year2+1:year3
            i=year-year1+1
                #println(year)
                #println(i)
                #println(i-year1+year2)
        csv1[i-year1+year2+1]=cohort_v1(i-year1+year2+1,year,tmax,jmax,
        1,tmax-i+1,1,
        map(j->x[i+j-1,j],1:tmax-i+1),
        map(j->y[i+j-1,j],1:tmax-i+1),
        map(j->w[i+j-1,j],1:tmax-i+1))
                #println( csv1[i-year0+year1])
    end

    cohorts_v2(year1,tmax,jmax,csv1)
end

"""
cohorts_v2->cohirts_t1

cohorts_t1:
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
x,y,w -  matrices  with abundunce, catch, weight

cohorts_v2:
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
cs - cohorts of mutable struct cohort_v1

cohort_v1:
no- the number
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
jbeg - initial age
len - length
jtrust - trusted age
xv,yv,wv- vectors with abundunce, catch, weight
"""
function csv2cst(csv)

    year1=csv.year_b
    tmax=csv.tmax
    jmax=csv.jmax
    cs=csv.cs

    x=fill(1.0,tmax,jmax)
    y=fill(1.0,tmax,jmax)
    w=fill(1.0,tmax,jmax)

    ncs=tmax+jmax-2
    year0=year1-jmax+1
    year2=year1+tmax-jmax
    year3=year2+jmax-1

    #println(year0)
    #println(year1)
    #println(year2)
    #println(year3)

   #1
    for year in year0:year1-1
        i=year-year0+1
        #println(year)
        #println(i)

    cv1=cs[i]
       no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
        jtrust=cv1.jtrust
        xv=cv1.xv
    yv=cv1.yv
    wv=cv1.wv

     #println(' ',i,' ',no,' ',yb,' ',jbeg,' ',len)

        for j in 1:i
            #println(j,' ',jmax-i+j)
            x[j,jmax-i+j]=xv[j]
            y[j,jmax-i+j]=yv[j]
            w[j,jmax-i+j]=wv[j]
        end
    end
  #2
     if tmax>=jmax
        for year in year1:year2
            i=year-year0+1
            i1=year-year1+1
                 #println(year)
                 #println(i)

            cv1=cs[i]
       no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
        jtrust=cv1.jtrust
        xv=cv1.xv
    yv=cv1.yv
    wv=cv1.wv

               #println(i,' ',no,' ',yb,' ',jbeg,' ',len)

            for j in 1:jmax
                #println(i+j-1,' ',j)
            x[i1+j-1,j]=xv[j]
            y[i1+j-1,j]=yv[j]
            w[i1+j-1,j]=wv[j]

        end
    end
        end

  #3
  for year in year2+1:year3
            i=year-year0+1
        i1=year-year1+1

                #println(year)
                #println(i)

        cv1=cs[i]
       no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
        jtrust=cv1.jtrust
        xv=cv1.xv
    yv=cv1.yv
    wv=cv1.wv

             #println(i,' ',no,' ',yb,' ',jbeg,' ',len)

        for j in 1:tmax-i1+1
             #println(i1+j-1,' ',j)

        x[i1+j-1,j]=xv[j]
            y[i1+j-1,j]=yv[j]
            w[i1+j-1,j]=wv[j]
        end
    end

    cohorts_t1(year1,tmax,jmax,x,y,w)
end

"""
cohorts related matrix (of mutable struct csmtx: x,y,w,m,f,z,g)->sequence of vectors (of mutable struct csvecs)

csrm:cohorts related matrix (of mutable struct csrm)
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
x,y,w, m,f,g -  cohorts related matrices

cvec: detail cohort representation as vec
no- the number in the sequence
yb - physical year associated with cohort (when age=1)
jbeg - initial age
len - length
jtrust - trusted age
vec

csvecs: cohorts related mtx as a sequence of cvec
no- the number
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
jbeg - initial age
len - length
jtrust - trusted age
xv,yv,wv- vectors with abundunce, catch, weight
tmax+jmax-1 -- number of related cohorts
"""
function csmtx2csvs(csrmtx)
   year1=csrmtx.year_b
    tmax=csrmtx.tmax
    jmax=csrmtx.jmax
    mtx=csrmtx.mtx

    year0=year1-jmax+1
      year2=year1+tmax-jmax
     year3=year2+jmax-1
    #println(year0)
    #println(year1)
    #println(year2)
    #println(year3)

    vecs=Array{cvec,1}(undef,tmax+jmax-1)

#    vecs=Array{Any}(undef,tmax+jmax-1)

    for year in year0:year1-1
        i=year-year0+1
        #println(year)
        #println(i)

        vecs[i]=cvec(i,year,
        jmax-i+1,i,argmax(map(j->mtx[j,jmax-i+j],1:i)),
        map(j->mtx[j,jmax-i+j],1:i))
        #println(vecs[i])
    end

    if tmax>=jmax
        for year in year1:year2
            i=year-year1+1
             #println(year)
             #println(i)
        #println(i-year0+year1)
            vecs[i-year0+year1]=cvec(i-year0+year1,year,
            1,jmax, argmax(map(j->mtx[i+j-1,j],1:jmax)),
            map(j->mtx[i+j-1,j],1:jmax))
            #println( vecs[i-year0+year1])
        end
    end

   #println(3)

    for year in year2+1:year3
            i=year-year1+1
            #println(year)
            #println(i)
            #println(i-year1+year2)
        vecs[i-year1+year2+1]=cvec(i-year1+year2+1,year,
        1,tmax-i+1,argmax(map(j->mtx[i+j-1,j],1:tmax-i+1)),
        map(j->mtx[i+j-1,j],1:tmax-i+1))
        #println( vecs[i-year1+year2+1])
    end

    csvecs(year1,tmax,jmax,vecs)
end

"""
sequence of vectors (of mutable struct csvecs)->cohorts related matrix (of mutable struct csrm) (x,y,w,m,f,z,g)

csrm:cohorts related matrix (of mutable struct csrm)
year_b - initial physical year
tmax - number of seasons
jmax - number of ages
x,y,w, m,f,g -  cohorts related matrices

cvec: detail cohort representation as vec
no- the number in the sequence
yb - physical year associated with cohort (when age=1)
jbeg - initial age
len - length
jtrust - trusted age
vec

csvecs: cohorts related mtx as a sequence of cvec
no- the number
yb - physical year associated with cohort (when age=1)
tmax - number of seasons
jmax - number of ages
jbeg - initial age
len - length
jtrust - trusted age
xv,yv,wv- vectors with abundunce, catch, weight
tmax+jmax-2 -- number of related cohorts
"""
function csvs2csmtx(csvs)

    year1=csvs.year_b
    tmax=csvs.tmax
    jmax=csvs.jmax
    vs=csvs.vecs

    mtx=fill(1.0,tmax,jmax)

    ncs=tmax+jmax-1
    year0=year1-jmax+1
    year2=year1+tmax-jmax
    year3=year2+jmax-1

    #println(year0)
    #println(year1)
    #println(year2)
    #println(year3)

   #1
    #println(1)
    for year in year0:year1-1
        i=year-year0+1
        #println(year)
        #println(i)

    cv1=vs[i]
    no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
    jtrust=cv1.jtrust
    xv=cv1.vec

     #println(' ',i,' ',no,' ',yb,' ',jbeg,' ',len)

        for j in 1:i
            #println(j,' ',jmax-i+j)
            mtx[j,jmax-i+j]=xv[j]

        end
    end
  #2
    #println(2)
     if tmax>=jmax
        for year in year1:year2
            i=year-year0+1
            i1=year-year1+1
            #println(year)
            #println(i)

    cv1=vs[i]
    no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
    jtrust=cv1.jtrust
    xv=cv1.vec

           #println(i,' ',no,' ',yb,' ',jbeg,' ',len)

            for j in 1:jmax
           # println(i+j-1,' ',j)
            mtx[i1+j-1,j]=xv[j]

        end
    end
        end

  #3
  #println(3)
    for year in year2+1:year3
            i=year-year0+1
        i1=year-year1+1

            #println(year)
            #println(i)

    cv1=vs[i]
    no=cv1.no
    yb=cv1.yb
    jbeg=cv1.jbeg
    len=cv1.len
    jtrust=cv1.jtrust
     xv=cv1.vec

            #println(i,' ',no,' ',yb,' ',jbeg,' ',len)

        for j in 1:tmax-i1+1
        # println(i1+j-1,' ',j)

        mtx[i1+j-1,j]=xv[j]

        end
    end

    csmtx(year1,tmax,jmax,mtx)
end

"""
generates xv1,yv1 of csvecs mutable struct with use of xv::csvecs,yv::csvecs,wv::csvecs,zv::csvecs,gv::csvecs
first makes a deep copy, then modifies vectors
"""
function gvecs(xv,yvs,wv,zv,gv)
    xv1=deepcopy(xv)
    yv1=deepcopy(yv)

    year1=xv.year_b
    tmax=xv.tmax
    jmax=xv.jmax
    xvecs=xv.vecs
    yvecs=yv.vecs
    wvecs=yv.vecs
    zvecs=zv.vecs
    gvecs=gv.vecs

    x1vecs=xv1.vecs
    y1vecs=yv1.vecs

    #println(xv1)

    # println( mvecs)
 #1
    for i in 2:jmax-1
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
        #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))
        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))
                y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]
        end
        end

        #println(x1vecs[i])

    end
   #2
   for i in jmax:tmax
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
         #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))
        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))
                 y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]
        end
        end

        #println(x1vecs[i])

    end
    #3
    for i in tmax+1:tmax+jmax-2
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
         #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))
        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))
                 y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]
        end
        end

        #println(x1vecs[i])
        #println(y1vecs[i])

    end

    return(xv1,yv1)

end

"""
generates xv1,yv1 of csvecs mutable struct with use of xv::csvecs,yv::csvecs,wv::csvecs,zv::csvecs,gv::csvecs
first makes a deep copy, then modifies vectors
"""
function gvecs1(xv,yv,wv,zv,gv)
    xv1=deepcopy(xv)
    yv1=deepcopy(yv)

    year1=xv.year_b
    tmax=xv.tmax
    jmax=xv.jmax
    xvecs=xv.vecs
    yvecs=yv.vecs
    wvecs=yv.vecs
    zvecs=zv.vecs
    gvecs=gv.vecs

    x1vecs=xv1.vecs
    y1vecs=yv1.vecs

    #println(xv1)

    # println( mvecs)
 #1
    for i in 2:jmax-1
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
        #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))

        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))

        end
        end

        #println(x1vecs[i])

    end

   #2
   for i in jmax:tmax
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
         #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))

        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))

        end
        end

        #println(x1vecs[i])

    end
    #3
    for i in tmax+1:tmax+jmax-2
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
        if jtr >1
        for i1 in 1:jtr-1
         #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(sum(zi[i1:jtr-1]))

        end
        end
        if jtr+1<len
        for i1 in jtr+1:len
          #println(i1)
            x1vecs[i].vec[i1]=xijtr*exp(-sum(zi[jtr+1:i1]))

        end
        end

        #println(x1vecs[i])
        #println(y1vecs[i])

    end

    ######

     for i in 1:jmax-1
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)

        for i1 in 1:len
        #println(i1)
           y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]

        end
                end

   #2
   for i in jmax:tmax
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

        #println(xijtr)
       for i1 in 1:len
        #println(i1)
           y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]

        end
                end

    #3
    for i in tmax+1:tmax+jmax-1
        xvi=xvecs[i]
        yvi=yvecs[i]
        wvi=wvecs[i]
        zvi=zvecs[i]
        gvi=gvecs[i]

        x1vi=x1vecs[i]
        y1vi=y1vecs[i]
      #  println(mvi)

         no=xvi.no
        yb=xvi.yb
        jbeg=xvi.jbeg
         len=xvi.len
        jtr=xvi.jtrust
        xi=xvi.vec

        yi=yvi.vec
        wi=wvi.vec
        zi=zvi.vec
        gi=gvi.vec

        x1i=x1vi.vec
        y1i=y1vi.vec

        #println(jbeg,' ',jtr)

        xijtr=xi[jtr]

        x1vecs[i].vec[jtr]=xijtr

       for i1 in 1:len
        #println(i1)
           y1vecs[i].vec[i1]= wi[i1]*x1vecs[i].vec[i1]*gi[i1]

        end
                end

    return(xv1,yv1)

end

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

"""
look how to evaluate closeness of data and model!!!!
as L1+L2
no prints out
"""
function Lp(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0

    for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+abs(log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t])^p
           #println(t," ",j," ",L1t)
        end
        L1=L1+L1t

    end

   #println(L1)

     local   L2=0.0
    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+abs(log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t])))^p

           #println(t," ",j," ",L2t)
        end;
        L2=L2+L2t;

    end
    #println(L2)
  # println(L1+L2)

    if L1+L2==Inf 10000.0 else L1+L2 end
end

"""
L as L1+L2 with dispersions for the normal distribution, p=2
no prints out
"""
function Lmine2(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0
    local  L11=0.0
    local mean1

    for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+
                    log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t]
           #println(t," ",j," ",L1t)
        end
      L1=L1+L1t
     end

        mean1=L1/((tmax-1)*(jmax-1))

     for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+
                    (log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t]-mean1)^2
           #println(t," ",j," ",L1t)
        end
        L11=L11+L1t
      end

   #println(L1)

     local   L2=0.0
        local  L22=0.0
        local mean2

    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+
                log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t]))

           #println(t," ",j," ",L2t)
        end;
        L2=L2+L2t;

    end

     mean2=L2/(tmax*jmax)

    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+
                    (log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t]))-mean2)^2

           #println(t," ",j," ",L2t)
        end;
        L22=L22+L2t;

    end

    #println(L2)
  # println(L1+L2)
    local result
    result=lambda *L11/((tmax-1)*(jmax-1)-1)+(1. - lambda)*L22/(tmax*jmax)

    if result==Inf 10000.0 else result end
end

"""
L as max
no prints out
"""
function Lmine2max(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0
    local  L11=0.0
    local mean1

    for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+
                    log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t]
           #println(t," ",j," ",L1t)
        end
      L1=L1+L1t
     end

        mean1=L1/((tmax-1)*(jmax-1))

     for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+
                    (log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t]-mean1)^2
           #println(t," ",j," ",L1t)
        end
        L11=L11+L1t
      end

   #println(L1)

     local   L2=0.0
        local  L22=0.0
        local mean2

    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+
                log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t]))

           #println(t," ",j," ",L2t)
        end;
        L2=L2+L2t;

    end

     mean2=L2/(tmax*jmax)

    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+
                    (log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t]))-mean2)^2

           #println(t," ",j," ",L2t)
        end;
        L22=L22+L2t;

    end

    #println(L2)
  # println(L1+L2)
    local result
    global lambda

    result=max(lambda*L11/((tmax-1)*(jmax-1)-1),(1. - lambda)*L22/(tmax*jmax))

    if result==Inf 10000.0 else result end
end

"""
look how to evaluate closeness of data and model!!!!
as L1+L2
"""
function Lp12(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0

    for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+abs(log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t])^p
           #println(t," ",j," ",L1t)
        end
        L1=L1+L1t

    end

   println(L1)

     local   L2=0.0
    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+abs(log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t])))^p

           #println(t," ",j," ",L2t)
        end;
        L2=L2+L2t;

    end
    println(L2)
   println(L1+L2)

    if L1+L2==Inf 10000.0 else L1+L2 end
end

function Lp1(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0

    for t=1:tmax-1
        L1t=0.0
        for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1t=L1t+abs(log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t])^p
           #println(t," ",j," ",L1t)
        end
        L1=L1+L1t

    end

   println(L1)

     local   L2=0.0
    for t=1:tmax
        L2t=0.0
        for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2t=L2t+abs(log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t])))^p

           #println(t," ",j," ",L2t)
        end;
        L2=L2+L2t;

    end
    println(L2)
   println(L1+L2)

    if L1+L2==Inf 10000.0 else L1+L2 end
end

function Lp2(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)

    function g(hh,mm,ss,ff)
   ss*ff/(hh*mm+ss*ff)*(1-exp(-hh*mm-ss*ff))
    end

    local  x=cst.x

      local  y=cst.y
     local   w=cst.w

     (tmax,jmax)=size(x)

   #println(tmax," ",jmax)
   # println(a1m,a2m,b1m,b2m,cm,d1m,d2m,as,bs)
   #println(htj)
   #println(ft)
     local  L1=0.0

    for t=1:tmax-1
                for j=1:jmax-1
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            L1=L1+abs(log(x[t,j]/x[t+1,j+1])-htj[t,j]*mj-sj*ft[t])^p
           #println(t," ",j," ",L1t)
        end

    end

   println(L1)

     local   L2=0.0
    for t=1:tmax
              for j= 1:jmax
            sj= sigma(j,as,bs,jmax)
            mj=ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)

            L2=L2+abs(log(y[t,j])-log(kw*w[t,j]*kx*x[t,j]*g(htj[t,j],mj,sj,ft[t])))^p

           #println(t," ",j," ",L2t)
        end;

    end
    println(L2)
   println(L1+L2)

    if L1+L2==Inf 10000.0 else L1+L2 end
end

"""
generates x,y with M,F,W,kx,kw,x1j,xt1
"""
function gm1(tmax,jmax,M,F,W,kx,kw,x1j,xt1)
##0
    x=fill(0.0,tmax,jmax);

    for j in [1:jmax;] x[1,j]=x1j[j] end

    for t in 1:tmax x[t,1]=xt1[t] end

 ##1
        for t in [1:tmax-jmax+1;]
           for i in [0:jmax-2;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end

 ##2
    for t in [tmax-jmax+2:tmax-1;]
           for i in [0:tmax-t-1;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end
   ##3

    for j in [2:jmax;]
           for i in [0:jmax-j-1;]
           #println(i+1+1,"! ",j+i+1,"<-", i+1," ",j+i)
            x[i+1+1,j+i+1]=x[i+1,j+i]*exp(-M[i+1,j+i]-F[i+1,j+i])
           end
        end

 ##0
     y=fill(0.0,tmax,jmax);

##1
   for t in 1:tmax
        for j in 1:jmax
            #println(t+i," ",i+1)
            Z=M[t,j]+F[t,j]
            y[t,j]=
            (kw*W[t,j])*(kx*x[t,j])*F[t,j]/Z*(1- exp(-Z))
        end
    end

    return x,y
end

"""
returns ovsyannikov's weight for j
"""
function ww_osp(t,j)
   wv=[0.012 0.052 0.11 0.179 0.253 0.323 0.483 0.677 0.908 1.116 1.297 1.541 1.758 1.971 1.997 2.358 2.475];
    return wv[j]
end

"""
returns total mortality z as sum of natural m and fishing f mortality matrices
"""
function gz(m::Array{Float64,2},f::Array{Float64,2},tmax,jmax)
    z=Array{Float64}(undef,tmax,jmax)
    z=m+f
    return(z)
end

"""
generates g matrix with m and f, tmax,jmax
"""
function gg(m::Array{Float64,2},f::Array{Float64,2},tmax,jmax)
    z=Array{Float64}(undef,tmax,jmax)
    g=Array{Float64}(undef,tmax,jmax)
    z=gz(m,f,tmax,jmax)
    for t in 1:tmax
        for j in 1:jmax
            g[t,j]=f[t,j]/z[t,j]*(1-exp(-z[t,j]))
        end
    end
    return(g)
end

"""
evaluates criteria
"""
function Lpv(cst,a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,as,bs,ft,kx,kw,p)
    local  x=cst.x

      local  y=cst.y
     local   w=cst.w
   local year1= cst.year_b

     (tmax,jmax)=size(x)

     M0=gmtj(a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,tmax,jmax);
    F0=gftj(as,bs,ft,tmax,jmax);
    Z0=M0+F0;
    G0=gg(M0,F0,tmax,jmax);

    xv=csmtx2csvs(csmtx(year1,tmax,jmax,x));
yv=csmtx2csvs(csmtx(year1,tmax,jmax,y));
    mv=csmtx2csvs(csmtx(year1,tmax,jmax,M0));
fv=csmtx2csvs(csmtx(year1,tmax,jmax,F0));
zv=csmtx2csvs(csmtx(year1,tmax,jmax,Z0));
gv=csmtx2csvs(csmtx(year1,tmax,jmax,G0));
wv=csmtx2csvs(csmtx(year1,tmax,jmax,w));
 xv0,yv0=gvecs(xv,yv,wv,zv,gv);

   x11=csvs2csmtx(xv0);
y11=csvs2csmtx(yv0);

    L=0;
for i in 1:tmax+jmax-1
    lenxv0=xv0.vecs[i].len;
   vecxv0=xv0.vecs[i].vec;
    vecyv0=yv0.vecs[i].vec;
    vecxv=xv.vecs[i].vec;
    vecyv=yv.vecs[i].vec;
           for n in 1:lenxv0
        L=L+abs(log(vecxv0[n])-log(vecxv[n]))^p+abs(log(vecyv0[n])-log(vecyv[n]))^p
    end

end

    if L==Inf 10000.0 else L end
end

"""
optimizes Lpv
"""
function optimizeLp_v(x,tmax,jmax,a,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

   global count = 0
opt = Opt(ALGS[a], 4*tmax+9)
#l_opt= Opt(ALGS[a1], 4*tmax+9)
    lower_bounds!(opt, fill(lb,4*tmax+9))
    #lower_bounds!(l_opt, fill(lb1,4*tmax+9))
    xtol_rel!(opt,xtol)
    #xtol_rel!(l_opt,xtol1)

min_objective!(opt, buragofunc_v)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

#local_optimizer!(opt,l_opt)

is=fill(0.4,4*tmax+9)
    is[1:3*tmax]=fill(is1,3*tmax)#h
    is[3*tmax+1:4*tmax]=fill(is2,tmax)#f
    is[4*tmax+1:4*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
#population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
returns Lpv
"""
function buragofunc_v(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1.,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lpv(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,1)

    #println("f_$count($x)=$y")

end

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
using QuasiGeostrophy, LinearAlgebra, Plots, BenchmarkTools
# Polynomial Order
n = 16

# Domain [a b]
Ω = [-1.0 1.0]
b = Ω[end]
a = Ω[1]
L = b - a

D, cx = cheb(n)

x =  (b-a) .* (cx .+ 1) ./ 2 .+ a

∂ˣ = 2 / L .* D

# Equation Dy = y, y[x=-1] = 1
A = ∂ˣ - I

# Boundary Conditions
A[end, :] .= 0.0
A[end, end] = 1.0

# Solve
rhs = 0.0 * x
rhs[end] = 1.0
luA = lu(A)
y = A \ rhs

# Exact solution
solution = exp.( x .+ 1)

plot(x, solution, label = "analytic", xlabel = "x", ylabel = "solution")
scatter!(x, y, label = "numerical")
plot!(legend = :topleft)

##
# performance check
n = m = 32
A = lu(randn(m*n,m*n))
b = ones(m*n)
x = copy(b)
@btime x = A \ b
@btime ldiv!(x, A, b)
@btime ldiv!(A, b)

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
module TestDegreeFunctions
  using Test
  using EcologicalNetworks

  # Generate some data

  N = UnipartiteProbabilisticNetwork([0.22 0.4; 0.3 0.1], [:a, :b])

  d_in = vec([0.52 0.5])
  d_out = vec([0.62 0.4])
  d_tot = vec([1.14 0.9])

  Q = BipartiteQuantitativeNetwork(rand(Float64, (3, 5)))
  D = convert(BinaryNetwork, Q)
  @test degree(Q) == degree(D)
  @test degree(Q; dims=1) == degree(D; dims=1)
  @test degree(Q; dims=2) == degree(D; dims=2)

  Dov = degree_var(N; dims=1)
  @test Dov[:b] ≈ 0.3

  Dv = degree_var(N; dims=2)
  @test Dv[:b] ≈ 0.33

  @test_throws ArgumentError degree(N; dims=3)
  @test_throws ArgumentError degree(N; dims=0)

  # Variance

  N = UnipartiteProbabilisticNetwork([0.22 0.4; 0.3 0.1])

  Din = degree(N; dims=2)
  Dout = degree(N; dims=1)
  Dtot = degree(N)

  Dov = degree_var(N; dims=1)
  @test Dov["s2"] ≈ 0.3 atol=0.001
  @test Dov["s1"] ≈ 0.41159 atol=0.001

  Div = degree_var(N; dims=2)
  @test Div["s1"] ≈ 0.3816 atol=0.01

  Dv = degree_var(N)
  @test Dv["s2"] ≈ 0.63 atol=0.001

end


text:
# This file was generated, do not modify it. # hide
combine(gdf, names(iris, Not(:Species)) .=> std)


text:
#=##############################################################################
# DESCRIPTION
    Tools for aerodynamic calculations.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

"""
    Calculates the aerodynamic force at each element of a VLM system using it's
current `Gamma` solution, the velocity induced by the particle field `pfield`,
and the kinematic velocity calculated from `prev_vlm_system`. It saves the force
as the field `vlm_system.sol["Ftot"]`
"""
function calc_aerodynamicforce(vlm_system::Union{vlm.Wing, vlm.WingSystem},
                                prev_vlm_system, pfield, Vinf, dt, rho; t=0.0,
                                vpm_solver="ExaFMM", per_unit_span=false,
                                Vout=nothing, lenout=nothing,
                                lencrit=-1)

    m = vlm.get_m(vlm_system)    # Number of horseshoes

    # Nodes of every horseshoe
    Ap = _get_Xs(vlm_system, "Ap")
    A = _get_Xs(vlm_system, "A")
    B = _get_Xs(vlm_system, "B")
    Bp = _get_Xs(vlm_system, "Bp")

    # Midpoints of bound vortices
    ApA = (Ap .+ A)/2
    AB = (A .+ B)/2
    BBp = (B .+ Bp)/2

    # Evaluate VPM on each midpoint
    Xs = vcat(ApA, AB, BBp)
    Vvpm = Vvpm_on_Xs(pfield, Xs; dt=dt)

    # Evaluate VLM on each midpoint
    Vvlm = vlm.Vind.(Ref(vlm_system), Xs; t=t, ign_col=true, ign_infvortex=true)

    # Evaluate Vinf on each midpoint
    Vinfs = Vinf.(Xs, t)

    # Evaluate kinematic velocity on each node
    Vtran = _Vkinematic(vlm_system, prev_vlm_system, dt; t=t,
                                                targetX=["Ap", "A", "B", "Bp"])

    Ftot = [zeros(3) for i in 1:m]

    for i in 1:m                 # Iterate over horseshoes
        for j in 1:3             # Iterate over bound vortices (1==ApA, 2==AB, 3==BBp)

            # Bound vortex' length vector
            if j==1
                l = (A[i]-Ap[i])
            elseif j==2
                l = (B[i]-A[i])
            else
                l = (Bp[i]-B[i])
            end

            # Bound vortex' midpoint
            # X = Xs[i + m*(j-1)]

            # Kinematic velocity at midpoint
            Vtranmid = (Vtran[i + m*(j-1)] + Vtran[i + m*j])/2

            # Effective velocity at midpoint
            V = Vvpm[i + m*(j-1)] + Vvlm[i + m*(j-1)] + Vinfs[i + m*(j-1)] + Vtranmid
            if Vout != nothing
                push!(Vout, [Vvpm[i + m*(j-1)], Vvlm[i + m*(j-1)], Vinfs[i + m*(j-1)], Vtranmid])
            end

            # Circulation
            Gamma = vlm_system.sol["Gamma"][i]

            # Normalization factor
            if per_unit_span
                # len = norm(cross(V/norm(V),l))
                unitV = V/sqrt(V[1]*V[1]+V[2]*V[2]+V[3]*V[3])
                projl = zeros(3)
                projl[1] = unitV[2]*l[3] - unitV[3]*l[2]
                projl[2] = unitV[3]*l[1] - unitV[1]*l[3]
                projl[3] = unitV[1]*l[2] - unitV[2]*l[1]
                len = sqrt(projl[1]*projl[1]+projl[2]*projl[2]+projl[3]*projl[3])
            else
                len = 1
            end
            if lenout != nothing
                push!(lenout, per_unit_span==false || len>lencrit ? len : -10*lencrit)
            end

            # Kutta–Joukowski theorem: F = rho * V × vecGamma
            if per_unit_span==false || len > lencrit # NOTE: Use lencrit to avoid dividing by zero
                Ftot[i][1] += rho * Gamma * (V[2]*l[3] - V[3]*l[2]) / len
                Ftot[i][2] += rho * Gamma * (V[3]*l[1] - V[1]*l[3]) / len
                Ftot[i][3] += rho * Gamma * (V[1]*l[2] - V[2]*l[1]) / len
            end

        end
    end

    vlm._addsolution(vlm_system, per_unit_span ? "ftot" : "Ftot", Ftot; t=t)

    return Ftot
end

"""
Receives an array of vectors `Fs` and decomposes them into components in the
`ihat`, `jhat`, `khat` unitary orthogonal directions (`khat` is deduced).
"""
function decompose(Fs, ihat, jhat)

    if abs(ihat[1]*ihat[1]+ihat[2]*ihat[2]+ihat[3]*ihat[3] - 1) >= 1e-12
        error("Got non-unitary vector ihat=$ihat")
    elseif abs(jhat[1]*jhat[1]+jhat[2]*jhat[2]+jhat[3]*jhat[3] - 1) >= 1e-12
        error("Got non-unitary vector jhat=$jhat")
    end

    khat = zeros(3)
    khat[1] = ihat[2]*jhat[3] - ihat[3]*jhat[2]
    khat[2] = ihat[3]*jhat[1] - ihat[1]*jhat[3]
    khat[3] = ihat[1]*jhat[2] - ihat[2]*jhat[1]

    Fis = [zeros(3) for i in 1:length(Fs)]
    Fjs = [zeros(3) for i in 1:length(Fs)]
    Fks = [zeros(3) for i in 1:length(Fs)]

    for (n, F) in enumerate(Fs)
        Fis[n] .= (F[1]*ihat[1]+F[2]*ihat[2]+F[3]*ihat[3]) * ihat
        Fjs[n] .= (F[1]*jhat[1]+F[2]*jhat[2]+F[3]*jhat[3]) * jhat
        Fks[n] .= (F[1]*khat[1]+F[2]*khat[2]+F[3]*khat[3]) * khat
    end

    return Fis, Fjs, Fks
end

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
# Equations:
# http://www.virtualtrebuchet.com/#documentation_EquationsOfMotion

function simulate(t)
   time = 0.0
   while !isa(t.stage, Val{:End})
      s = simulate_(t, time)
      time = s.t[end]
      derive!(t, s)
      transition(t, s[end])
   end
   return t.sol
end

simulate_(t::TrebuchetState, time) = simulate_(t, time, t.stage)

function terminate(it, tr)
    tr.terminated = tr.stage
    terminate!(it)
end

function simulate_(t::TrebuchetState, time, ::Val{:Ground})
    a = t.a
    aw = t.aw
    u0 = oftype.(t.m.w, [a.aq, a.wq, a.sq, aw.aw, aw.ww, aw.sw])
    ti = oftype.(t.m.w, (time, time + 1.0))
    prob = ODEProblem(stage1!, u0, ti, t)

    string_tension = (u, time, it, θ) -> begin
        m = t.m.p
        time == it.tprev && return 0
        acc = (sling_velocity(t, u) - sling_velocity(t, it.uprev))/(time - it.tprev)
        Tn = -acc.x*m/sin(θ)
    end

    ccb = ContinuousCallback(
        (u, time, it) -> begin
            mg = t.m.p * t.c.Grav
            θ = u[1] + u[3] - π
            Tn = string_tension(u, time, it, θ)
            t.Tn = Tn
            mg - Tn*cos(θ)
        end,
        (it) -> terminate(it, t),
        abstol=1e-2, reltol=1e-2)

    dcb = DiscreteCallback(
        (u, time, it) -> t.Tn < zero(typeof(t.Tn)), # loose string case
        (i) -> terminate!(i))

    solve(prob, Euler(), dt = 1/t.rate, callback=CallbackSet(ccb, dcb))
end

function simulate_(t::TrebuchetState, time, ::Val{:Hang})
    a, aw = t.a, t.aw
    r = t.c.r
    u0 = oftype.(t.m.w, [a.aq, a.wq, a.sq, aw.aw, aw.ww, aw.sw])
    ti = oftype.(t.m.w, (time, time + 1.0))
    prob = ODEProblem(stage2!, u0, ti, t)

    cb = ContinuousCallback(
        (u, time, it) -> sin(projectile_angle(t, u) - r),
        (it) -> terminate(it, t))

    solve(prob, Euler(), dt = 1/t.rate, callback=cb)
end

function simulate_(t::TrebuchetState, time, ::Val{:Released})
    a = t.l.a
    u0 = oftype.(t.m.w, [t.p.x, t.p.y, t.v.x, t.v.y])
    ti = oftype.(t.m.w, (time, 5.0 + time))
    prob = ODEProblem(stage3!, u0, ti, t)

    cb = ContinuousCallback(
        (u, time, it) -> u[2] + a,
        (it) -> terminate(it, t))

    solve(prob, Euler(), dt = 1/t.rate, callback=cb)
end

function stage1!(du, u, p::TrebuchetState, t)
    l, a, m, i = p.l, p.a, p.m, p.i
    LAl, LAs, LW, LS, h = l.b, l.c, l.d, l.e, l.a
    mA, mW, mP = m.a, m.w, m.p
    IA3, IW3 = i.iw, i.ia
    LAcg = (LAl - LAs)/2
    Grav = p.c.Grav
    SIN = sin
    COS = cos

    Aq = u[1]
    Wq = u[2]
    Sq = u[3]
    Aw = u[4]
    Ww = u[5]
    Sw = u[6]

    M11 = -mP*LAl^2*(-1+2*SIN(Aq)*COS(Sq)/SIN(Aq+Sq)) + IA3 + IW3 + mA*LAcg^2 + mP*LAl^2*SIN(Aq)^2/SIN(Aq+Sq)^2 + mW*(LAs^2+LW^2+2*LAs*LW*COS(Wq))
    M12 = IW3 + LW*mW*(LW+LAs*COS(Wq))
    M21 = IW3 + LW*mW*(LW+LAs*COS(Wq))
    M22 = IW3 + mW*LW^2
    r1 = Grav*LAcg*mA*SIN(Aq) + LAl*LS*mP*(SIN(Sq)*(Aw+Sw)^2+COS(Sq)*(COS(Aq+Sq)*Sw*(Sw+2*Aw)/SIN(Aq+Sq)+(COS(Aq+Sq)/SIN(Aq+Sq)+LAl*COS(Aq)/(LS*SIN(Aq+Sq)))*Aw^2)) + LAl*mP*SIN(Aq)*(LAl*SIN(Sq)*Aw^2-LS*(COS(Aq+Sq)*Sw*(Sw+2*Aw)/SIN(Aq+Sq)+(COS(Aq+Sq)/SIN(Aq+Sq)+LAl*COS(Aq)/(LS*SIN(Aq+Sq)))*Aw^2))/SIN(Aq+Sq) - Grav*mW*(LAs*SIN(Aq)+LW*SIN(Aq+Wq)) - LAs*LW*mW*SIN(Wq)*(Aw^2-(Aw+Ww)^2)
    r2 = -LW*mW*(Grav*SIN(Aq+Wq)+LAs*SIN(Wq)*Aw^2)

    dAq = Aw
    dWq = Ww
    dSq = Sw
    dAw = (r1*M22-r2*M12)/(M11*M22-M12*M21)
    dWw = -(r1*M21-r2*M11)/(M11*M22-M12*M21)
    dSw = -COS(Aq+Sq)*dSq*(dSq+2*dAq)/SIN(Aq+Sq) - (COS(Aq+Sq)/SIN(Aq+Sq)+LAl*COS(Aq)/(LS*SIN(Aq+Sq)))*dAq^2 - (LAl*SIN(Aq)+LS*SIN(Aq+Sq))*dAw/(LS*SIN(Aq+Sq))

    du[1] = dAq
    du[2] = dWq
    du[3] = dSq
    du[4] = dAw
    du[5] = dWw
    du[6] = dSw
end

function stage2!(du, u, p::TrebuchetState, t)
    l, a, m, i = p.l, p.a, p.m, p.i
    LAl, LAs, LW, LS, h = l.b, l.c, l.d, l.e, l.a
    mA, mW, mP = m.a, m.w, m.p
    IA3, IW3 = i.iw, i.ia

    LAcg = (LAl - LAs)/2
    Grav = p.c.Grav
    COS = cos
    SIN = sin

    Aq = u[1]
    Wq = u[2]
    Sq = u[3]
    Aw = u[4]
    Ww = u[5]
    Sw = u[6]

    M11 = IA3 + IW3 + mA*LAcg^2 + mP*(LAl^2+LS^2+2*LAl*LS*COS(Sq)) + mW*(LAs^2+LW^2+2*LAs*LW*COS(Wq))
    M12 = IW3 + LW*mW*(LW+LAs*COS(Wq))
    M13 = LS*mP*(LS+LAl*COS(Sq))
    M21 = IW3 + LW*mW*(LW+LAs*COS(Wq))
    M22 = IW3 + mW*LW^2
    M31 = LS*mP*(LS+LAl*COS(Sq))
    M33 = mP*LS^2

    r1 = Grav * LAcg * mA * SIN(Aq) + Grav * mP * ( LAl * SIN(Aq) + LS * SIN(Aq+Sq) ) - Grav * mW * ( LAs * SIN(Aq) + LW* SIN(Aq+Wq)) - LAl * LS * mP * SIN(Sq) * ( Aw^2 - ( Aw + Sw )^2 ) - LAs * LW * mW * SIN(Wq) * ( Aw^2 - (Aw+Ww)^2)
    r2 = -LW * mW * (Grav * SIN(Aq+Wq) + LAs * SIN(Wq) * Aw^2)
    r3 = LS * mP * ( Grav * SIN(Aq + Sq) - LAl * SIN(Sq) * Aw^2 )

    dAq = Aw
    dWq = Ww
    dSq = Sw
    dAw = -( r1 * M22 * M33 - r2 * M12 * M33 - r3 * M13 * M22 ) / ( M13 * M22 * M31 - M33 * ( M11 * M22 - M12 * M21))
    dWw = ( r1 * M21 * M33 - r2 * (M11 * M33 - M13 * M31 ) - r3 * M13 * M21) / ( M13 * M22 * M31 - M33 * ( M11 * M22 - M12 * M21))
    dSw = ( r1 * M22 * M31 - r2 * M12 * M31 - r3 * (M11 * M22 - M12 * M21)) / (M13 * M22 * M31 - M33 * ( M11 * M22 - M12 * M21))

    du[1] = dAq
    du[2] = dWq
    du[3] = dSq
    du[4] = dAw
    du[5] = dWw
    du[6] = dSw
end

function stage3!(du, u, pr::TrebuchetState, t)
    c = pr.c
    Grav = c.Grav
    ρ = c.ρ
    WS = c.w
    Aeff = π*(pr.l.z)^2
    mP = pr.m.p
    Cd = c.Cd

    Pvx = u[3]
    Pvy = u[4]

    dPx = Pvx
    dPy = Pvy
    dPvx = -(ρ*Cd*Aeff*(Pvx-WS)*sqrt(Pvy^2+(WS-Pvx)^2))/(2*mP)
    dPvy = -Grav - (ρ*Cd*Aeff*Pvy*sqrt(Pvy^2+(WS-Pvx)^2))/(2*mP)

    du[1] = dPx
    du[2] = dPy
    du[3] = dPvx
    du[4] = dPvy
end


text:
try
	cd(Pkg.dir("EGR"))
catch
	cd("/Users/stepa/Google Drive/Research/EGRproject/EGR.jl")
end

using EGR

include("TestM1P1.jl"     )
include("TestRedisReadWrite.jl")

Oracles = Tuple{Function,Int64,Int64,DataType,ASCIIString,Outputter,Bool,ASCIIString,Function}[]
include("LoadTestToyData.jl")
include("LoadTestToySparseData.jl")
include("LoadTestAgaricusData.jl")

include("TestVerify.jl")
include("TestEGRlikeSG.jl")
include("TestFindBest.jl")
include("TestRedisRuns.jl")
include("TestSSVRG.jl")
include("TestSAG.jl")
include("TestSAGA.jl")
include("TestSAGAinit.jl")
include("TestSAGinit.jl")
include("TestUniversal.jl")


text:
# "Escaping the swirling vortex"
#
# Nonautonomous ODE from pg 121 of Hale's ODE book, credited to
# Markus, L. and H. Yamabe "Global stability criteria for diferential systems",
# Osaka Math. J. 12 (1960), 305-317.
#
# A nonautonomous linear ODE of the form `du/dt = A(t)u`, with all the eigenvalues
# of `A(t)` being complex with negative real part.
#
# However, origin is an unstable fixed point!
#
# The vector field visualization was based on
# https://gist.github.com/Wikunia/4ca58b4c08da90978c0dcdee5b0c908a

using Javis
using Colors
using LinearAlgebra

"""
    ground(args...)

Define the background and the initial "foreground" color.
"""
function ground(args...)
    background("black")
    sethue("white")
    return nothing
end

"""
dudt(x, y, t)

Return the right hand side of the linear equation `u'=A(t)u`,
where u=[x,y]. In this case, the equation takes the from
    x' = (-1 + 1.5cos(t)^2)x + (1 - 1.5cos(t)sin(t))y
    y' = (-1 - 1.5sin(t)cos(t))x + (-1 + 1.5sin(t)^2)y
"""
function dudt(x, y, t)
    ct = cos(t)
    st = sin(t)
    dxdt = (-1.0 + 1.5 * ct^2) * x + (1.0 - 1.5 * ct * st) * y
    dydt = (-1.0 - 1.5 * st * ct) * x + (-1.0 + 1.5 * st^2) * y
    return dxdt, dydt
end

"""
    jacobian!(mat, t)

Return (and modify in place) the linear part `A(t)` of the equation `u'=A(t)u`
(see [`dudt`](@ref)). This is also seen as the Jacobian of `u -> f(t,u) = A(t)u`.
"""
function jacobian!(mat, t)
    ct = cos(t)
    st = sin(t)
    mat[1, 1] = -1.0 + 1.5 * ct^2
    mat[1, 2] = 1.0 - 1.5 * ct * st
    mat[2, 1] = -1.0 - 1.5 * st * ct
    mat[2, 2] = -1.0 + 1.5 * st^2
    return mat
end

"""
    solution!(u,t)

Return a particular solution of the equation `u'=A(t)u`, given by
`u = exp(t/2)  * [-cos(t), sin(t)] / 20.0` corresponding to the initial condition
`u(0) = [-0.05, 0.0]`
"""
function solution!(u, t)
    et2 = exp(t / 2) / 20.0
    u[1] = -et2 * cos(t)
    u[2] = et2 * sin(t)
    return u
end

"""
    space_to_frame_coordinates(u, frm_scl)

Transform a point in xy coordinates to the frame coordinates XY.
They are related according to `X = frm_scl * x` and `Y = - frm_scl * y`,
for a given scale factor `frm_scl`.
"""
space_to_frame_coordinates(u, frm_scl) = frm_scl * Point(u[1], -u[2])
space_to_frame_coordinates(x, y, frm_scl) = space_to_frame_coordinates([x, y], frm_scl)

"""
    ellipse_level_set!(u, t, mat, σ₁, σ₂, frm_scl)

Draw a "level-set"-like ellipse associated with the solution at time `t`.

For a given time `t`, the solutions `w=w(s)` of the autonomous equation
`dw(s)/ds = A(t)w(s)` spiral towards the origin, decreasing along the level
sets of an elliptic paraboloid. This function computes and draws the ellipse
level set containing the solution `u(t)`. Note, however, that `u=u(t)` is
the solution of the nonautonomous equation `du(t)/dt = A(t)u(t)` given
by [`solution!`](@ref), which spirals away from the origin.

The ellipse is computed by finding the complex eigenvector `v=v(t)` of the
complex eigenvalue `α + iβ` of `A(t)` and considering the real and invertible
matrix `P(t) = [Re(v(t)) Im(v(t))]` which transforms `A(t)` into its real
Jordan form `P(t)^-1 * A(t) P(t) = [α β, -β α]`. In this particular case, the
eigenvalues `α ± iβ` are time-independent. The ellipse is obtained from the
SVD decomposition `P = UΣVᵗ`, where `Σ = diag(σ₁, σ₂)` is diagonal with
nonnegative elements `σ₁ ≥ σ₂ ≥ 0`, and `U` and `V` are orthogonal matrices.
In general, all these terms may change in time, but in the particular
example we are considering, only `U=U(t)` and `V=V(t)` are time-dependent,
with `σ₁ > σ₂ > 0` being constant. The elements `σ₁` and `σ₂` are, precisely,
the semi-major and semi-minor axes of the ellipse associated with the level
set of the elliptic paraboloid at level one. It determines the eccentricity
of the ellipse. The orientation of the ellipse is given by the rotation
determined by the orthogonal matrix `U=U(t)`.

At each time `t`, we find and draw the level set that contains the solution `u(t)`
as follows. We rotate `u(t)` via `Uᵗ = inv(U)`, compute the semi-major
and semi-minor axis of the ellipse that has the same level as `Uᵗu(t)`, compute
the foci of this ellipse, rotate the foci back to the right orientation by applying
the orthogonal matrix `U`, and then find the ellipse, in frame coordinates, that
defined by the foci, in frame coordinates, and passing by the solution `u(t)`, in
frame coordinates.
"""
function ellipse_level_set!(u, t, mat, σ₁, σ₂, frm_scl)
    sethue("white")
    jacobian!(mat, t)
    v = eigen(mat).vectors[:, 1]
    U = svd([real(v) imag(v)]).U
    x, y = transpose(U) * solution!(u, t)
    γ = sqrt(x^2 / σ₁^2 + y^2 / σ₂^2) # level
    σ̃₁ = γ * σ₁ # semi-major
    σ̃₂ = γ * σ₂ # semi-minor
    fd = sqrt(abs(σ̃₁^2 - σ̃₂^2)) # focal distance to the origin/center
    f = U * [fd, 0.0] # one focus point; the other is -f
    poly(
        ellipse(
            space_to_frame_coordinates(f, frm_scl), # focus in frame coordinates
            -space_to_frame_coordinates(f, frm_scl), # the other focus
            space_to_frame_coordinates(u, frm_scl), # ellipse/level containing u(t)
        ),
        :stroke,
        close = true,
    )
    return nothing
end

"""
    ball(p=O, size=4, color="red")

Draw a filled circle of size `size` and color `color` around the point `p`
"""
function ball(p = O, size = 4, color = "red")
    sethue(color)
    circle(p, size, :fill)
    return p
end

"""
    path!(points, pos, color)

Update and draw the trail of the solution.

Add the point `pos` to the vector `points` containing the trail of the trajectory
and draw the trail.
"""
function path!(points, pos, color)
    sethue(color)
    push!(points, pos) # add pos to points
    prettypoly(points) # draw dotted path with circles
    return nothing
end

"""
    ode_arrow(x, y, t, frm_scl, color_min, color_max)

Create one arrow at the base point `(x,y)` and arrow vector given by `dudt(x, y, t)`.
The arrow is properly scaled by the factor `frm_scl` and the orientation takes into
account that the vertical coordinate increases downwards. Depending on the magnitude
of the vector field at that point, the color ranges from `color_min` to `color_max`.
"""
function ode_arrow(x, y, t, frm_scl, color_min, color_max)
    @JShape begin
        translate(space_to_frame_coordinates(x, y, frm_scl))
        dx, dy = dudt(x, y, t)
        delta_t = 0.1
        ξ = dx * delta_t
        η = dy * delta_t
        l = sqrt(ξ^2 + η^2)
        len_colors = 100
        color_range = range(color(color_min), stop = color(color_max), length = len_colors)
        idx = clamp(floor(Int, l * len_colors), 1, len_colors)
        sethue(color_range[idx])
        if ξ != 0 || η != 0
            arrow(
                O,
                space_to_frame_coordinates(ξ, η, frm_scl);
                linewidth = 2,
                arrowheadlength = 8,
            )
        end
    end x = x y = y t = t frm_scl = frm_scl
end

"""
    horizontal_grid_line(video, c, frm_scl)

Create an horizontal grid line at `y=c`
"""
function horizontal_grid_line(video, r, frm_scl)
    if r == 0
        setline(1.5)
    end
    sethue("gray")
    line(Point(-video.width / 2, r * frm_scl), Point(video.width / 2, r * frm_scl), :stroke)
    return nothing
end

"""
    vertical_grid_line(video, c, frm_scl)

Create a vertical grid line at `x=c`.
"""
function vertical_grid_line(video, c, frm_scl)
    if c == 0
        setline(1.5)
    end
    sethue("gray")
    line(
        Point(c * frm_scl, -video.height / 2),
        Point(c * frm_scl, video.height / 2),
        :stroke,
    )
    return nothing
end

"""
    animate()

Create the animation. 😃
"""
function animate()
    # initial setup
    vid = Video(500, 500)
    nframes = 200
    xgrid = -2.0:0.25:2.0 # grid for bars and arrow base points
    ygrid = -2.0:0.25:2.0 # grid for bars and arrow base points
    t₀, tf = 0.0, 5π / 2 # time interval for the simulation
    u = [0.0, 0.0] # dummy initizaliation of the solution vector
    solution!(u, t₀) # initial vector
    mat = zeros(2, 2) # dummy initialization of the jacobian
    jacobian!(mat, t₀) # jacobian at time t₀
    v = eigen(mat).vectors[:, 1] # first eigenvector
    svd_P = svd([real(v) imag(v)]) # P = [real(v) imag(v)] => inv(P) * mat * P = [α β; -β α]
    σ₁, σ₂ = svd_P.S # the semi-major and semi-minor axes are constant in this problem
    t = t₀ # set t to initial time
    frm_scl = 120 # frame scale
    solution_trail = Point[] # initialize history vector for the solution path

    # objects
    ## set background for all frames
    Background(1:nframes, ground)
    ## grids
    [
        Object(1:nframes, (args...) -> horizontal_grid_line(args[1], x, frm_scl)) for
        x in xgrid
    ]
    [Object(1:nframes, (args...) -> vertical_grid_line(args[1], y, frm_scl)) for y in ygrid]
    ## vector field
    arrows = [
        Object(1:nframes, ode_arrow(x, y, t, frm_scl, "yellow", "red")) for x in xgrid,
        y in ygrid
    ]
    ## solution at time t
    sol_pos = Object(
        1:nframes,
        (args...; t) -> ball(
            space_to_frame_coordinates(solution!(u, t), frm_scl),
            6,
            "lightskyblue",
        ),
    )
    ## solution trail
    Object(1:nframes, (args...) -> path!(solution_trail, pos(sol_pos), "lightskyblue"))
    ## ellipse
    obj_ellipse =
        Object(1:nframes, (args...; t) -> ellipse_level_set!(u, t, mat, σ₁, σ₂, frm_scl))

    # animation
    ## update solution, vector field and ellipse
    ## rmk: animation of the solution trail is a byproduct of updating solution
    act!(sol_pos, Action(1:nframes, change(:t, t₀ => tf)))
    act!(arrows, Action(1:nframes, change(:t, t₀ => tf)))
    act!(obj_ellipse, Action(1:nframes, change(:t, t₀ => tf)))

    # rendering
    ## render and generate gif
    render(
        vid;
        framerate = 15,
        pathname = joinpath(@__DIR__, "gifs/escaping_the_swirling_vortex.gif"),
    )

    return nothing
end

animate()


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb02fd(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    ITMAX = 10
    KMAX = 20
    NMAX = 20
    LDR = NMAX*KMAX
    LDT = KMAX
    LDWORK = ( NMAX + 1 )*KMAX
    S = Array{BlasInt,1}(undef, ITMAX)
    tdim = NMAX*KMAX
    #T = Array{Float64,2}(undef, LDT,NMAX*KMAX)
    T = fill(zero(Float64), LDT,NMAX*KMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    #R = Array{Float64,2}(undef, LDR,NMAX*KMAX)
    R = fill(zero(Float64),LDR,NMAX*KMAX)
    V = Array{Float64,1}(undef, NMAX*KMAX)
    W = Array{Float64,1}(undef, NMAX*KMAX)
    Z = Array{Float64,1}(undef, NMAX*KMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    K = parse(BlasInt, vs[2])
    IT = parse(BlasInt, vs[3])
    TYPET = 'R'
    M = N*K
    local NNRM
    if ( N<=0 || N>NMAX )
        @error "Illegal N=$N"
    elseif ( K<=0 || K>KMAX )
        @error "Illegal K=$K"
    elseif ( IT<=0 || IT>ITMAX )
        @error "Illegal IT=$IT"
    end

    vs = String[]
    _isz = IT
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    S[1:_isz] .= parsex.(BlasInt, vs)

    vs = String[]
    _isz,_jsz = (K,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       T[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
    P = 0
    POS = 1
    for SCIT in 1:IT # do 90
# interp call 1

        #Rtmp = R[POS:LDR,POS:tdim]
        INFO = SLICOT.mb02fd!(TYPET, K, N, P, S[SCIT], view(T,1:LDT,POS:tdim),
                              view(R,POS:LDR,POS:tdim), LDWORK)
        #R[POS:LDR,POS:tdim] .= Rtmp

        @test INFO == 0
        if ( INFO!=0 )
            return
        end # if
        S1 = S[SCIT] + P
        if ( S1==0 )
            LEN = N*K
            # interp call 2

            #FOREIGN.dlaset!( 'All', LEN, 1, ONE, ONE, V, 1 )
            V[1:LEN] .= ONE

            # interp call 3
            for PIT in 1:5 # do 30
                for i in 1:N
                    #FOREIGN.dgemv!( 'NoTranspose', K, LEN-(I-1)*K, ONE, T, LDT, V((I-1)*K+1), 1, ZERO, W((I-1)*K+1), 1 )
                    #W[(i-1)*K+1:i*K] .= T[1:K,1:LEN-(i-1)*K] * V[(i-1)*K+1:LEN]
                    _nc = LEN-(i-1)*K
                    BLAS.gemv!('N',ONE,T[1:K,1:_nc],V[(i-1)*K+1:LEN],ZERO,view(W,(i-1)*K+1:i*K))
                end

                # interp call 4

                for i in 1:N-1
                    #FOREIGN.dgemv!( 'Transpose', K, (N-I)*K, ONE, T(1,K+1), LDT, V((I-1)*K+1), 1, ONE, W(I*K+1), 1 )
                    #W[i*K+1:N*K] .= T[1:K,K+1:(N-i+1)*K]' * V[(i-1)*K+1:i*K]
                    BLAS.gemv!('T',ONE,T[1:K,K+1:(N-i+1)*K],V[(i-1)*K+1:i*K],ONE,view(W,i*K+1:N*K))
                end
                # interp call 5

                #FOREIGN.dcopy!( LEN, W, 1, V, 1 )
                V[1:LEN] .= W[1:LEN]
            end # for PIT
            #NNRM = DNRM2( LEN, V, 1 )
            NNRM = norm(V[1:LEN])
            # interp call 6

            #FOREIGN.dscal!( LEN, ONE/NNRM, V, 1 )
            V[1:LEN] .*= ONE / NNRM

        else
            LEN = ( N - S1 )*K
            # interp call 7

            #FOREIGN.dlaset!( 'All', LEN, 1, ONE, ONE, V, 1 )
            V[1:LEN] .= ONE

            for PIT in 1:5 # do 80
                POSR = ( S1 - 1 )*K + 1
                for i in 1:N-S1 # do 40
                    # interp call 8

                    #FOREIGN.dgemv!( 'NoTranspose', K, LEN-(I-1)*K, ONE, T(1,POSR+K), LDT, V((I-1)*K+1), 1, ZERO, W((I-1)*K+1), 1 )
                    _nc = LEN-(i-1)*K
                    #W[(i-1)*K+1:i*K] .= T[1:K,POSR+K:POSR+K+_nc-1] * V[(i-1)*K+1:LEN]
                    BLAS.gemv!('N',ONE,T[1:K,POSR+K:POSR+K+_nc-1],V[(i-1)*K+1:LEN],ZERO,
                               view(W,(i-1)*K+1:i*K))
                end

                for i in 1:N-S1 # do 50
                    # interp call 9

                    #FOREIGN.dtrmv!( 'U', 'N', 'N', K, R[POSR,POSR], LDR, V[(I-1)*K+1], 1 )
                    #V[(i-1)*K+1:i*K] .= UpperTriangular(R[POSR:POSR+K-1,POSR:POSR+K-1]) * V[(i-1)*K+1:i*K]
                    BLAS.trmv!('U','N','N',R[POSR:POSR+K-1,POSR:POSR+K-1],view(V,(i-1)*K+1:i*K))

                    # interp call 10

                    #FOREIGN.dgemv!( 'N', K, LEN-I*K, ONE, R[POSR,POSR+K], LDR, V[I*K+1], 1, ONE, V[(I-1)*K+1], 1 )
                    #V[(i-1)*K+1:i*K] .= R[POSR:POSR+K-1,POSR+K:POSR+K+LEN-i*K-1] * V[i*K+1:LEN]
                    BLAS.gemv!('N',ONE,R[POSR:POSR+K-1,POSR+K:POSR+K+LEN-i*K-1],
                               V[i*K+1:LEN],ONE,view(V,(i-1)*K+1:i*K))
                end
                # interp call 11

                #FOREIGN.dlaset!( 'All', LEN, 1, ZERO, ZERO, Z, 1 )
                Z[1:LEN] .= ZERO

                for i in 1:N-S1 # do 60
                    # interp call 12

                    #FOREIGN.dgemv!( 'T', K, LEN-I*K, ONE, R[POSR,POSR+K], LDR, V[(I-1)*K+1], 1, ONE, Z[I*K+1], 1 )
                    #Z[i*K+1:LEN] .+= R[POSR:POSR+K-1,POSR:POSR+LEN-i*K-1]' * V[(i-1)*K+1:i*K]
                    BLAS.gemv!('T',ONE,R[POSR:POSR+K-1,POSR:POSR+LEN-i*K-1],V[(i-1)*K+1:i*K],
                               ONE,view(Z,i*K+1:LEN))

                    # interp call 13

                    #FOREIGN.dtrmv!( 'U', 'T', 'N', K, R[POSR,POSR], LDR, V[(I-1)*K+1], 1 )
                    #V[(i-1)*K+1:i*K] .= UpperTriangular(R[POSR:POSR+K-1,POSR:POSR+K-1])' * V[(i-1)*K+1:i*K]
                    BLAS.trmv!('U','T','N',R[POSR:POSR+K-1,POSR:POSR+K-1],
                               view(V,(i-1)*K+1:i*K))
                    # interp call 14

                    #FOREIGN.daxpy!( K, ONE, V((I-1)*K+1), 1, Z((I-1)*K+1), 1 )
                    Z[(i-1)*K+1:i*K] .+= V[(i-1)*K+1:i*K]
                end

                # interp call 15

                # FOREIGN.dlaset!( 'All', LEN, 1, ZERO, ZERO, V, 1 )
                V[1:LEN] .= ZERO

                for i in 1:N-S1 # do 70
                    # interp call 16

                    #FOREIGN.dgemv!( 'T', K, LEN-(I-1)*K, ONE, T(1,POSR+K), LDT, W((I-1)*K+1), 1, ONE, V((I-1)*K+1), 1 )
                    #V[(i-1)*K+1:LEN] .+= T[1:K,POSR+K:POSR+K+LEN-(i-1)*K-1]' * W[(i-1)*K+1:i*K]
                    BLAS.gemv!('T',ONE,T[1:K,POSR+K:POSR+K+LEN-(i-1)*K-1],W[(i-1)*K+1:i*K],
                               ONE,view(V,(i-1)*K+1:LEN))
                end
                # interp call 17

                #FOREIGN.daxpy!( LEN, -ONE, Z, 1, V, 1 )
                #V[1:LEN] .-= Z[1:LEN]
                BLAS.axpy!(-ONE,Z[1:LEN],view(V,1:LEN))

                # NNRM = DNRM2( LEN, V, 1 )
                NNRM = norm(V[1:LEN])
                # interp call 18

                # FOREIGN.dscal!( LEN, -ONE/NNRM, V, 1 )
                V[1:LEN] .*= (-ONE / NNRM)
            end # for / do 80
            POS = ( S1 - 1 )*K + 1
            P = S1
        end # if
        println(io, "rows: $(P*K); norm of Schur complement: $NNRM")

    end # for / do 90
    # interp output 1
    println(io, "R:")
    _nc = M
    _nr = P*K
    show(io, "text/plain", R[1:_nr,1:_nc])
    println(io)

end # run_mb02fd()


text:
abstract type StaticInjection <: Device end
supports_services(d::Device) = true
get_dynamic_injector(d::StaticInjection) = nothing


text:
using CRlibm
using FlexFloat

using Polynomials
using Quaternions
using DualNumbers