text:
using JuliaChem.JCModules

using Base.Threads
using LinearAlgebra

#=
"""
	 index(a::Int64,b::Int64)
Summary
======
Triangular indexing determination.

Arguments
======
a = row index

b = column index
"""
=#
@inline function triangular_index(a::Int,b::Int)
  index = (a*(a-1)) >> 1 #bitwise divide by 2
  index += b
  return index
end

@inline function triangular_index(a::Int)
  return (a*(a-1)) >> 1
end

@inline function decompose(input::Int)
  #return ceil(Int,(-1.0+√(1+8*input))/2.0)
  return Base.fptosi(Int, Base.ceil_llvm((-1.0 +
    Base.Math.sqrt_llvm(float(1+8*input)))/2.0))
    #return ccall((:decompose, "/export/home/david/projects/Julia/JuliaChem.jl/src/eri/libjeri.so"),
  #  Int64, (Int64,), input)
end

function compute_nuc_grad(mol::Molecule)
  nuc_grad = zeros(Float64,(length(mol.atoms),3))
  for iatom in 1:length(mol.atoms), jatom in 1:length(mol.atoms)
    if iatom != jatom
      ix = mol.atoms[iatom].atom_center[1]
      jx = mol.atoms[jatom].atom_center[1]

      iy = mol.atoms[iatom].atom_center[2]
      jy = mol.atoms[jatom].atom_center[2]

      iz = mol.atoms[iatom].atom_center[3]
      jz = mol.atoms[jatom].atom_center[3]

      iatm = mol.atoms[iatom].atom_id
      jatm = mol.atoms[jatom].atom_id

      distance = √((jx-ix)^2 + (jy-iy)^2 + (jz-iz)^2)

      nuc_grad[jatom,1] -= iatm*jatm*(jx-ix)/(distance^3)
      nuc_grad[jatom,2] -= iatm*jatm*(jy-iy)/(distance^3)
      nuc_grad[jatom,3] -= iatm*jatm*(jz-iz)/(distance^3)
    end
  end

  return nuc_grad
end

function compute_overlap_grad(mol::Molecule,
  basis::Basis, W::Matrix{Float64}, jeri_oei_grad_engine)

  #== define initial variables ==#
  natoms = length(mol.atoms)
  WS_grad = zeros(Float64,(natoms,3))

  ncoord = length(WS_grad)
  shell_set = 2

  #== generate S derivative matrices ==#
  S_deriv = Vector{Matrix{Float64}}([ zeros(Float64,(basis.norb, basis.norb)) for i in 1:ncoord ])
  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    iatom = basis[ash].atom_id
    jatom = basis[bsh].atom_id

    S_block_JERI = zeros(Float64,(abas*bbas*ncoord))

    JERI.compute_overlap_grad_block(jeri_oei_grad_engine, S_block_JERI, ash, bsh,
      abas*bbas)

    axial_normalization_factor(S_block_JERI, basis[ash], basis[bsh], ncoord)

    #println("Julia Shells: $ash, $bsh")
    #println("Julia Atoms: $iatom, $jatom")

    shlset_idx = 1
    for ishlset in 1:shell_set
      atom = ishlset == 1 ? iatom : (ishlset == 2 ? jatom : ishlset-2)

      for icoord in 1:3
        idx = 1
        op = 3*(atom-1) + icoord

        #println("$ishlset, $icoord => $op, $shlset_idx")

        for ibas in 0:abas-1, jbas in 0:bbas-1
          iorb = apos + ibas
          jorb = bpos + jbas

          S_deriv[op][max(iorb,jorb),min(iorb,jorb)] += S_block_JERI[abas*bbas*(op-1) + idx]
          idx += 1
        end
        shlset_idx += 1
      end
    end
  end

  for ideriv in S_deriv
    for iorb in 1:basis.norb, jorb in 1:(iorb-1)
      ideriv[min(iorb,jorb),max(iorb,jorb)] = ideriv[max(iorb,jorb),min(iorb,jorb)]
    end
    #display(ideriv); println()
  end

  #== contract with energy-weighted density ==#
  for iatom in 1:natoms
    for icoord in 1:3
      deriv_idx = 3*(iatom-1)  + icoord
      for ibas in 1:basis.norb, jbas in 1:basis.norb
        scale = ibas == jbas ? 0.5 : 1.0

        WS_grad[iatom,icoord] -= scale * W[ibas,jbas] * S_deriv[deriv_idx][ibas, jbas]
      end
    end
  end
  return WS_grad
end

function compute_kinetic_grad(mol::Molecule,
  basis::Basis, P::Matrix{Float64}, jeri_oei_grad_engine)

  #== define initial variables ==#
  natoms = length(mol.atoms)
  PT_grad = zeros(Float64,(natoms,3))

  ncoord = length(PT_grad)
  shell_set = 2

  #== generate T derivative matrices ==#
  T_deriv = Vector{Matrix{Float64}}([ zeros(Float64,(basis.norb, basis.norb)) for i in 1:ncoord ])
  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    iatom = basis[ash].atom_id
    jatom = basis[bsh].atom_id

    T_block_JERI = zeros(Float64,(abas*bbas*ncoord))

    JERI.compute_kinetic_grad_block(jeri_oei_grad_engine, T_block_JERI, ash, bsh,
      abas*bbas)

    axial_normalization_factor(T_block_JERI, basis[ash], basis[bsh], ncoord)

    #println("Julia Shells: $ash, $bsh")
    #println("Julia Atoms: $iatom, $jatom")

    shlset_idx = 1
    for ishlset in 1:shell_set
      atom = ishlset == 1 ? iatom : (ishlset == 2 ? jatom : ishlset-2)

      for icoord in 1:3
        idx = 1
        op = 3*(atom-1) + icoord

        #println("$ishlset, $icoord => $op, $shlset_idx")

        for ibas in 0:abas-1, jbas in 0:bbas-1
          iorb = apos + ibas
          jorb = bpos + jbas

          T_deriv[op][max(iorb,jorb),min(iorb,jorb)] += T_block_JERI[abas*bbas*(op-1) + idx]
          idx += 1
        end
        shlset_idx += 1
      end
    end
  end

  for ideriv in T_deriv
    for iorb in 1:basis.norb, jorb in 1:(iorb-1)
      ideriv[min(iorb,jorb),max(iorb,jorb)] = ideriv[max(iorb,jorb),min(iorb,jorb)]
    end
    #display(ideriv); println()
  end

  #== contract with energy-weighted density ==#
  for iatom in 1:natoms
    for icoord in 1:3
      deriv_idx = 3*(iatom-1)  + icoord
      for ibas in 1:basis.norb, jbas in 1:basis.norb
        scale = ibas == jbas ? 0.5 : 1.0

        PT_grad[iatom,icoord] += scale * P[ibas,jbas] * T_deriv[deriv_idx][ibas, jbas]
      end
    end
  end
  return PT_grad
end

function compute_nuc_attr_grad(mol::Molecule,
  basis::Basis, P::Matrix{Float64}, jeri_oei_grad_engine)

  #== define initial variables ==#
  natoms = length(mol.atoms)
  PV_grad = zeros(Float64,(natoms,3))

  ncoord = length(PV_grad)
  shell_set = 2 + length(mol.atoms)

  #== generate T derivative matrices ==#
  V_deriv = Vector{Matrix{Float64}}([ zeros(Float64,(basis.norb, basis.norb)) for i in 1:ncoord ])
  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    iatom = basis[ash].atom_id
    jatom = basis[bsh].atom_id

    V_block_JERI = zeros(Float64,(abas*bbas*ncoord))

    JERI.compute_nuc_attr_grad_block(jeri_oei_grad_engine, V_block_JERI, ash, bsh,
      abas*bbas)

    axial_normalization_factor(V_block_JERI, basis[ash], basis[bsh], ncoord)

    #println("Julia Shells:", ash-1, ", ", bsh-1)
    #println("Julia Atoms:", iatom-1, ", ", jatom-1)

    shlset_idx = 1
    for ishlset in 1:shell_set
      atom = ishlset == 1 ? iatom : (ishlset == 2 ? jatom : ishlset-2)

      for icoord in 1:3
        idx = 1
        op = 3*(atom-1) + icoord
        #println(ishlset-1, ", ", icoord-1, " => ", op-1, ", ", shlset_idx-1)

        for ibas in 0:abas-1, jbas in 0:bbas-1
          iorb = apos + ibas
          jorb = bpos + jbas

          #println(abas*bbas*(op-1) + idx - 1)

          V_deriv[op][max(iorb,jorb),min(iorb,jorb)] += V_block_JERI[abas*bbas*(op-1) + idx]
          idx += 1
        end
        shlset_idx += 1
      end
    end
  end

  for ideriv in V_deriv
    for iorb in 1:basis.norb, jorb in 1:(iorb-1)
      ideriv[min(iorb,jorb),max(iorb,jorb)] = ideriv[max(iorb,jorb),min(iorb,jorb)]
    end
    #display(ideriv); println()
  end

  #== contract with energy-weighted density ==#
  for iatom in 1:natoms
    for icoord in 1:3
      deriv_idx = 3*(iatom-1)  + icoord

      for ibas in 1:basis.norb, jbas in 1:basis.norb
        scale = ibas == jbas ? 0.5 : 1.0

        PV_grad[iatom,icoord] += scale * P[ibas,jbas] * V_deriv[deriv_idx][ibas, jbas]
      end
    end
  end
  return PV_grad
end

function compute_schwarz_bounds(schwarz_bounds::Matrix{Float64},
  basis::Basis, nsh::Int64)

  max_am = max_ang_mom(basis)
  eri_quartet_batch = Vector{Float64}(undef,eri_quartet_batch_size(max_am))
  jeri_schwarz_engine = JERI.TEIEngine(basis.basis_cxx, basis.shpdata_cxx)

  for ash in 1:nsh, bsh in 1:ash
    fill!(eri_quartet_batch, 0.0)

    abas = basis[ash].nbas
    bbas = basis[bsh].nbas
    abshp = triangular_index(ash, bsh)

    JERI.compute_eri_block(jeri_schwarz_engine, eri_quartet_batch,
      ash, bsh, ash, bsh, abshp, abshp, abas*bbas, abas*bbas)

    schwarz_bounds[ash, bsh] = sqrt(maximum(abs.(eri_quartet_batch)) )
  end

  for ash in 1:nsh, bsh in 1:ash
    if ash != bsh
      schwarz_bounds[min(ash,bsh),max(ash,bsh)] =
        schwarz_bounds[max(ash,bsh),min(ash,bsh)]
    end
  end
end

#=
"""
		get_oei_matrix(oei::Array{Float64,2})
Summary
======
Extract one-electron integrals from data file object. Kinetic energy integrals,
overlap integrals, and nuclear attraction integrals can all be extracted.

Arguments
======
oei = array of one-electron integrals to extract
"""
=#
function read_in_oei(oei::Vector{T}, nbf::Int) where T
	nbf2 = (nbf*(nbf+1)) >> 1

	oei_matrix = Matrix{Float64}(undef,(nbf,nbf))
	for ibf in 1:nbf2
    i = decompose(ibf)
    j = ibf - triangular_index(i)

		oei_matrix[i,j] = float(oei[ibf])
		oei_matrix[j,i] = oei_matrix[i,j]
	end

	return oei_matrix
end

function DIIS(F::Matrix{Float64}, e_array::Vector{Matrix{Float64}},
  F_array::Vector{Matrix{Float64}}, B_dim::Int64)

  B = Matrix{Float64}(undef,B_dim+1,B_dim+1)
  for i in 1:B_dim, j in 1:B_dim
    B[i,j] = LinearAlgebra.dot(e_array[i], e_array[j])

	  B[i,B_dim+1] = -1
	  B[B_dim+1,i] = -1
	  B[B_dim+1,B_dim+1] =  0
  end
  #DIIS_coeff::Vector{Float64} = [ fill(0.0,B_dim)..., -1.0 ]
  DIIS_coeff::Vector{Float64} = vcat(zeros(B_dim), [-1.0])

  #DIIS_coeff[:], B[:,:], ipiv = LinearAlgebra.LAPACK.gesv!(B, DIIS_coeff)
  DIIS_coeff[:], B[:,:], ipiv = LinearAlgebra.LAPACK.sysv!('U', B, DIIS_coeff)

  fill!(F, zero(Float64))
  for index in 1:B_dim
    F .+= DIIS_coeff[index] .* F_array[index]
  end
end

function axial_normalization_factor(oei, ash, bsh, ncoord)
  ama = ash.am
  amb = bsh.am

  na = ash.nbas
  nb = bsh.nbas

  for icoord in 1:ncoord
    ab = (icoord-1)*na*nb
    for asize::Int64 in 0:(na-1), bsize::Int64 in 0:(nb-1)
      ab += 1

      anorm = axial_norm_fact[asize+1,ama]
      bnorm = axial_norm_fact[bsize+1,amb]

      abnorm = anorm*bnorm
      oei[ab] *= abnorm
    end
  end
end

function compute_dipole_moment(mol::Molecule,
  basis::Basis, P::Matrix{Float64}, jeri_prop_engine)

  #== define initial variables ==#
  natoms = length(mol.atoms)
  dipole = zeros(Float64,(3,))

  ncoord = length(dipole)

  #== compute nuclear contribution to dipole ==#
  nuc_dipole = similar(dipole)
  for atom in mol
    nuc_dipole .+= atom.atom_id .* atom.atom_center
  end

  #== compute electronic contribution to dipole ==#
  elec_dipole = similar(dipole)

  elec_dipole_matrix = Vector{Matrix{Float64}}([ zeros(Float64,(basis.norb, basis.norb)) for i in 1:ncoord ])
  for ash in 1:length(basis), bsh in 1:ash
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    iatom = basis[ash].atom_id
    jatom = basis[bsh].atom_id

    elec_dipole_block_JERI = zeros(Float64,(abas*bbas*ncoord))

    JERI.compute_dipole_block(jeri_prop_engine, elec_dipole_block_JERI, ash, bsh,
      abas*bbas)

    axial_normalization_factor(elec_dipole_block_JERI, basis[ash], basis[bsh], ncoord)

    #println("Julia Shells: $ash, $bsh")
    #println("Julia Atoms: $iatom, $jatom")

    for icoord in 1:ncoord
      idx = 1

      #println("$ishlset, $icoord => $op, $shlset_idx")

      for ibas in 0:abas-1, jbas in 0:bbas-1
        iorb = apos + ibas
        jorb = bpos + jbas

        elec_dipole_matrix[icoord][max(iorb,jorb),min(iorb,jorb)] += elec_dipole_block_JERI[abas*bbas*(icoord-1) + idx]
        idx += 1
      end
    end
  end

  for imatrix in elec_dipole_matrix
    for iorb in 1:basis.norb, jorb in 1:(iorb-1)
      imatrix[min(iorb,jorb),max(iorb,jorb)] = imatrix[max(iorb,jorb),min(iorb,jorb)]
    end
    #display(ideriv); println()
  end

  #== contract with energy-weighted density ==#
  for icoord in 1:ncoord
    for ibas in 1:basis.norb, jbas in 1:basis.norb
      elec_dipole[icoord] += P[ibas,jbas] * elec_dipole_matrix[icoord][ibas, jbas]
    end
  end

  dipole .= 2.5417 .* (nuc_dipole .- elec_dipole)
  return dipole
end

function eri_quartet_batch_size(max_am)
  return am_to_nbas_cart(max_am)^4
end


text:
transform_coordinates(xyz, twist, section) = RotY(-twist) * xyz + section

function chop_sections(scaled_foils, twisties, leading_xyz, span_num, spacings, flip = false)
    # Reverse direction if left side
    if flip
        scaled_foils = reverse(scaled_foils, dims = 2)
        twisties     = (permutedims ∘ reverse)(twisties)
    end

    # Rotate and translate coordinates
    foil_coords  = reduce(hcat, transform_coordinates.(foil, Ref(twist), Ref(section)) for (foil, twist, section) in zip(eachcol(scaled_foils), twisties, leading_xyz))

    # Chop up spanwise sections
    permutedims(reduce(hcat, chop_coordinates(coords, span_num, spacings, flip) for coords in eachrow(foil_coords)))
end

"""
    surface_coordinates(wing :: HalfWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the coordinates of a `HalfWing` consisting of `Foil`s and relevant geometric quantities, given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels, with an option to flip the signs of the ``y``-coordinates.
"""
function surface_coordinates(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = ["cosine"], flip = false)
    leading_xyz  = leading_edge(wing, flip)
    scaled_foils = reduce(hcat, @. wing.chords * (extend_yz ∘ cosine_foil)(wing.foils, chord_num))
    chop_sections(scaled_foils, twists(wing), leading_xyz, span_num, spacings, flip)
end

"""
    camber_coordinates(wing :: HalfWing, n_s :: Integer, n_c :: Integer, flip = false)

Compute the coordinates of a `HalfWing` consisting of camber distributions of `Foil`s and relevant geometric quantities, given numbers of spanwise ``n_s`` and chordwise ``n_c`` panels, with an option to flip the signs of the ``y``-coordinates.
"""
function camber_coordinates(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = ["cosine"], flip = false)
    leading_xyz  = leading_edge(wing, flip)
    scaled_foils = reduce(hcat, @. wing.chords * (camber_coordinates ∘ camber_thickness)(wing.foils, chord_num))
    chop_sections(scaled_foils, twists(wing), leading_xyz, span_num, spacings, flip)
end

"""
    mesh_horseshoes(wing :: HalfWing, n_s :: Integer, n_c :: Integer; flip = false)

Mesh a `HalfWing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for lifting-line/vortex lattice analyses using horseshoe elements.
"""
mesh_horseshoes(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer, spacings = ["cosine"]; flip = false) = make_panels(coordinates(wing, span_num, chord_num; span_spacing = spacings, flip = flip))

"""
    mesh_wing(wing :: HalfWing, n_s :: Integer, n_c :: Integer; flip = false)

Mesh a `HalfWing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for 3D analyses using doublet-source elements or equivalent formulations. TODO: Tip meshing.
"""
mesh_wing(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; flip = false) = (make_panels ∘ surface_coordinates)(wing, span_num, chord_num, flip)

"""
    mesh_cambers(wing :: HalfWing, n_s :: Integer, n_c :: Integer; flip = false)

Mesh the camber distribution of a `HalfWing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for vortex lattice analyses.
"""
mesh_cambers(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer, spacings = ["cosine"]; flip = false) = make_panels(camber_coordinates(wing, span_num, chord_num; spacings = spacings, flip = flip))

chop_leading_edge(obj :: HalfWing, span_num; flip = false)  = chop_coordinates(leading_edge(obj, flip), span_num)
chop_trailing_edge(obj :: HalfWing, span_num; flip = false) = chop_coordinates(leading_edge(obj, flip), span_num)

chop_leading_edge(obj :: Wing, span_num :: Integer)  = let (lead, trail) = wing_bounds(obj); chop_spans(lead, trail, span_num)[1] end
chop_trailing_edge(obj :: Wing, span_num :: Integer) = let (lead, trail) = wing_bounds(obj); chop_spans(lead, trail, span_num)[2] end

"""
    mesh_horseshoes(wing :: Wing, n_s :: Integer, n_c :: Integer)

Mesh a `Wing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for lifting-line/vortex lattice analyses using horseshoe elements.
"""
function mesh_horseshoes(wing :: Wing, span_num, chord_num, spacings = ["cosine"])
    left_panels  = mesh_horseshoes(left(wing), reverse(span_num), chord_num, reverse(spacings), flip = true)
    right_panels = mesh_horseshoes(right(wing), span_num, chord_num, spacings)

    [ left_panels right_panels ]
end

"""
    mesh_wing(wing :: Wing, n_s :: Integer, n_c :: Integer)

Mesh a `Wing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for 3D analyses using doublet-source elements or equivalent formulations. TODO: Tip meshing.
"""
function mesh_wing(wing :: Wing, span_num, chord_num)
    left_panels  = mesh_wing(left(wing), reverse(span_num), chord_num, flip = true)
    right_panels = mesh_wing(right(wing), span_num, chord_num)

    [ left_panels right_panels ]
end

"""
    mesh_cambers(wing :: Wing, n_s :: Integer, n_c :: Integer)

Mesh the camber distribution of a `Wing` into panels of ``n_s`` spanwise divisions per section and ``n_c`` chordwise divisions meant for vortex lattice analyses.
"""
function mesh_cambers(wing :: Wing, span_num, chord_num, spacings = ["cosine"])
    left_panels  = mesh_cambers(left(wing), reverse(span_num), chord_num, reverse(spacings), flip = true)
    right_panels = mesh_cambers(right(wing), span_num, chord_num, spacings)

    [ left_panels right_panels ]
end

vlmesh_wing(wing :: Union{HalfWing, Wing}, span_num, chord_num, spacings = ["cosine"]) = mesh_horseshoes(wing, span_num, chord_num, spacings), mesh_cambers(wing, span_num, chord_num, spacings)

function paneller(wing :: Union{Wing, HalfWing}, span_num :: Union{Integer, Vector{<: Integer}}, chord_num :: Integer; rotation = one(RotMatrix{3, Float64}), translation = zeros(3), spacings = ["cosine"])
    horseshoes, cambers = vlmesh_wing(wing, ifelse(typeof(span_num) <: Integer, [span_num], span_num), chord_num, spacings)
    horseshoe_panels    = [ transform(panel, rotation, translation)               for panel in horseshoes ]
    panel_normals       = [ panel_normal(transform(panel, rotation, translation)) for panel in cambers 	  ]

    horseshoe_panels, panel_normals
end

panel_wing(comp :: Union{Wing, HalfWing}, span_panels :: Union{Integer, Vector{<: Integer}}, chord_panels :: Integer; position = zeros(3), angle = 0., axis = [1., 0., 0.], spacing = spanwise_spacing(comp)) = paneller(comp, span_panels, chord_panels, rotation = AngleAxis{Float64}(angle, axis...), translation = position, spacings = ifelse(typeof(spacing) <: String, [spacing], spacing))

number_of_spanwise_panels(wing :: HalfWing, span_num :: Integer) = ceil.(Int, span_num .* spans(wing) / span(wing))
number_of_spanwise_panels(wing :: Wing,     span_num :: Integer) = number_of_spanwise_panels(right(wing), ceil(Int, span_num / 2))

function number_of_spanwise_panels(wing :: HalfWing, span_num :: Vector{<: Integer})
    @assert (length ∘ spans)(wing) > 1 "Provide an integer number of spanwise panels for 1 wing section."
    span_num
end

function number_of_spanwise_panels(wing :: Wing, span_num :: Vector{<: Integer})
    @assert (length ∘ spans ∘ right)(wing) > 1 "Provide an integer number of spanwise panels for 1 wing section."
    ceil.(Int, span_num ./ 2)
end

spanwise_spacing(wing :: HalfWing) = [ "sine"; fill("cosine", (length ∘ spans)(wing)         - 1) ]
spanwise_spacing(wing :: Wing)     = [ "sine"; fill("cosine", (length ∘ spans ∘ right)(wing) - 1) ]

coordinates(wing :: Union{HalfWing, Wing}) = let (lead, trail) = wing_bounds(wing); permutedims([ lead trail ]) end

function coordinates(wing :: HalfWing, span_num, chord_num; span_spacing = ["cosine"], flip = false)
    lead, trail = wing_bounds(wing, flip)
    reduce(hcat, chop_wing(lead, trail, span_num, chord_num; span_spacing = span_spacing, flip = flip))
end

coordinates(wing :: Wing, span_num, chord_num; span_spacing = ["cosine"], rotation = one(RotMatrix{3, Float64}), translation = zeros(3)) = [ coordinates(left(wing), span_num, chord_num; span_spacing = span_spacing, flip = true)[:,1:end-1] coordinates(right(wing), span_num, chord_num; span_spacing = span_spacing) ]

ext: jl

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
module MatDeforI1RivlinADModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D, nstressstrain, nthermstrain
import FinEtoolsDeforLinear.MatDeforModule: AbstractMatDefor, stressvtot!, stressttov!, strainttov!, strainvdet, strainvtr
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, totlag2currsymm!
using LinearAlgebra: Transpose, Diagonal, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
using LinearAlgebra: eigen, eigvals, norm, cholesky, cross, dot, log, diagm, det
using ForwardDiff: gradient, hessian

"""
	MatDeforI1RivlinAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear

Type for triaxial two-term I1-based Rivlin hyperelastic material. Described in
K Bertoldi and others, Journal of Mechanics and Physics of Solids 56 (2008)
2642-2668.

Material parameters suggested for one specific type of rubber in that paper
was specified by parameters
```
c1 = 0.55 MPa
c2 = 0.3 MPa
K = 50 * mu = 55 MPa
```
Initial Young's modulus 3.25 MPa, which corresponds to nu = 0.4772727
"""
struct  MatDeforI1RivlinAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear
	mr::Type{MR} # model reduction type
	mass_density::FFlt # mass density
	c1::FFlt # material constant: c1 = mu/2 (half the shear modulus).
	c2::FFlt # material constant: neo-hookean results if c2 = 0.
	K::FFlt # bulk modulus
	tangentmoduli!::MTAN # Function to return the tangent moduli matrix
	update!::MUPD # Function to update the material state
end

################################################################################
# 3-D solid model
################################################################################

"""
    MatDeforI1RivlinAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)

Create triaxial two-term I1-based Rivlin hyperelastic material.
"""
function MatDeforI1RivlinAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, c1::FFlt, c2::FFlt, K::FFlt)
	_mI = diagm(0=>[1, 1, 1, 0.5, 0.5, 0.5]);
	_m1 = vec(FFlt[1 1 1 0 0 0]) ;
	_m1m1 = _m1*_m1';
	_I3 = [1.0 0 0; 0 1.0 0; 0 0 1.0]
	# neo-hookean:  mu/2*(trC-3) - mu*lJ + lambda/2*(lJ)^2;
	# two-term I1-based Rivlin:
	function strainenergy(Cv, mr, c1, c2, K)
		trC = strainvtr(mr, Cv)
		J = sqrt(strainvdet(mr, Cv))
		lJ = log(J)
		return c1*(trC-3) + c2*(trC-3)^2 - 2*c1*(lJ) + K/2*(J - 1)^2;
	end
	function tangentmoduli3d!(self::MatDeforI1RivlinAD, D::FFltMat, statev::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat, label::FInt)
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strainttov!(mr, Cv, C)
		Dtotal = 4 .* hessian(Cv -> strainenergy(Cv, mr, self.c1, self.c2, self.K), Cv);
		return totlag2currsymm!(D, Dtotal, Fn1)
	end
	function update3d!(self::MatDeforI1RivlinAD, statev::FFltVec, cauchy::FFltVec, output::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
		@assert length(cauchy) == nstressstrain(self.mr)
		# Note: the code below allocates. This may not be a big deal when the
		# code gets called only a few times, but in explicit dynamics
		# situations  this would be very expensive. The material type needs to
		# provide these temporary arrays as buffers.
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strainttov!(mr, Cv, C)
		Sv = 2 * gradient(Cv -> strainenergy(Cv, mr, self.c1, self.c2, self.K), Cv);
		S = fill(0.0, 3, 3)
		stressvtot!(mr, S, Sv)
		cauchyt = Fn1*(S/det(Fn1))*Fn1'; # Cauchy stress
		stressttov!(mr, cauchy, cauchyt)
		if quantity == :nothing
			#Nothing to be copied to the output array
		elseif quantity == :cauchy || quantity == :Cauchy
			(length(output) >= 6) || (output = zeros(6)) # make sure we can store it
			copyto!(output, cauchy);
		elseif quantity == :pressure || quantity == :Pressure
			output[1]  =  -sum(cauchy[1:3])/3.
		elseif quantity == :princCauchy || quantity == :princcauchy
			ep = eigen(sigma);
			(length(output) >= 3) || (output = zeros(3)) # make sure we can store it
			copyto!(output,  sort(ep.values, rev=true));
		elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
			s1=cauchy[1]; s2=cauchy[2]; s3=cauchy[3];
			s4=cauchy[4]; s5=cauchy[5]; s6=cauchy[6];
			(length(output) >= 1) || (output = zeros(1)) # make sure we can store it
			output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
		end
		return output
	end
	return MatDeforI1RivlinAD(mr, mass_density, c1, c2, K,
		tangentmoduli3d!, update3d!)
end

"""
    MatDeforI1RivlinAD(mr::Type{MR}, c1::FFlt, c2::FFlt, K::FFlt) where {MR}

Create two-term I1-based Rivlin hyperelastic material.

The mass density is the default value.
"""
function MatDeforI1RivlinAD(mr::Type{MR}, c1::FFlt, c2::FFlt, K::FFlt) where {MR}
	mass_density = 1.0
	return MatDeforI1RivlinAD(mr, mass_density, c1, c2, K)
end

"""
    MatDeforI1RivlinAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}

Create two-term I1-based Rivlin hyperelastic material equivalent to a simple neo-hookean material.

The mass density is the default value.
"""
function MatDeforI1RivlinAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}
	mass_density = 1.0
	c1 = E / 2 / (1 + nu) / 2 # half of the shear modulus
	c2 = 0.0
	K = E / 3 / (1 - 2*nu) # bulk modulus
	return MatDeforI1RivlinAD(mr, mass_density, c1, c2, K)
end

end


text:
module RocketFilterTypeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: filter_type()" begin

    println("Testing: operator filter_type()")

    run_proxyshowcheck("FilterType", filter_type(Int))

    @test_throws Exception from(1:5) |> filter_type(Float64)

    run_testset([
        (
            source      = from(1:5) |> filter_type(Int),
            values      = @ts([ 1, 2, 3, 4, 5, c ]),
            source_type = Int
        ),
        (
            source      = from(Real[ 1, 2.0, 3, 4.0, 5, 6.0 ]) |> filter_type(Int),
            values      = @ts([ 1, 3, 5, c ]),
            source_type = Int
        ),
        (
            source      = from(Union{Int, Float64}[ 1, 2.0, 3, 4.0, 5, 6.0 ]) |> filter_type(Float64),
            values      = @ts([ 2.0, 4.0, 6.0, c ]),
            source_type = Float64
        ),
        (
            source      = completed() |> filter_type(Int),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted("e") |> filter_type(Int),
            values      = @ts(e("e")),
            source_type = Int
        ),
        (
            source      = never() |> filter_type(Int),
            values      = @ts(),
            source_type = Int
        )
    ])

end

end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function Control(t,x,rate_vector,data_dictionary)
# ---------------------------------------------------------------------- #
# Control.jl was generated using the Kwatee code generation system.
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 03-31-2016 14:11:33
#
# Arguments:
# t  - current time
# x  - state vector
# rate_vector - vector of reaction rates
# data_dictionary  - Data dictionary instance (holds model parameters)
# ---------------------------------------------------------------------- #

# Set a default value for the allosteric control variables -
EPSILON = 1.0e-3;
number_of_reactions = length(rate_vector);
control_vector = ones(number_of_reactions);
control_parameter_array = data_dictionary["CONTROL_PARAMETER_ARRAY"];

# Alias the species vector -
gene_Smad = x[1];
gene_Ecadherin = x[2];
gene_Bcatenin = x[3];
gene_PLCy = x[4];
gene_PKC = x[5];
gene_RAF = x[6];
gene_RAS = x[7];
gene_AP1_SP1 = x[8];
gene_tubulin = x[9];
gene_SNAIL_SLUG = x[10];
gene_Calcineurin = x[11];
gene_GSK3 = x[12];
mRNA_Smad = x[13];
mRNA_Ecadherin = x[14];
mRNA_Bcatenin = x[15];
mRNA_PLCy = x[16];
mRNA_PKC = x[17];
mRNA_RAF = x[18];
mRNA_RAS = x[19];
mRNA_AP1_SP1 = x[20];
mRNA_Calcineurin = x[21];
mRNA_SNAIL_SLUG = x[22];
mRNA_tubulin = x[23];
mRNA_GSK3 = x[24];
PLCy = x[25];
Active_PLCy = x[26];
Calcineurin = x[27];
Active_Calcineurin = x[28];
NFATc = x[29];
Active_NFATc = x[30];
PKC = x[31];
Active_PKC = x[32];
RAF = x[33];
RAF_P = x[34];
Smad = x[35];
pSmad = x[36];
RAS = x[37];
RAS_GTP = x[38];
MAPK = x[39];
pMAPK = x[40];
AP1_SP1 = x[41];
AP1_SP1_P = x[42];
Ecadherin = x[43];
Bcatenin = x[44];
Ecadherin_Bcatenin = x[45];
TCF4 = x[46];
Active_TCF4 = x[47];
LEF1 = x[48];
Active_LEF1 = x[49];
PPase = x[50];
YREG1 = x[51];
SNAIL_SLUG = x[52];
vimentin = x[53];
GSK3 = x[54];
tubulin = x[55];
RNAP = x[56];
RIBOSOME = x[57];
U0126 = x[58];
TGFB3_Ab = x[59];
ss_virus = x[60];
smad_virus = x[61];
lef1_virus = x[62];
gene_TCF4 = x[63];
gene_TGFB3 = x[64];
gene_YREG1 = x[65];
gene_LEF1 = x[66];
gene_vimentin = x[67];
gene_R = x[68];
gene_VEGFR = x[69];
gene_NFATc = x[70];
gene_MAPK = x[71];
mRNA_TCF4 = x[72];
mRNA_TGFB3 = x[73];
mRNA_YREG1 = x[74];
mRNA_LEF1 = x[75];
mRNA_vimentin = x[76];
mRNA_R = x[77];
mRNA_VEGFR = x[78];
mRNA_NFATc = x[79];
mRNA_MAPK = x[80];
R = x[81];
TGFB12 = x[82];
TGFB3 = x[83];
VEGF = x[84];
TGFB12_R = x[85];
TGFB3_R = x[86];
VEGFR = x[87];
VEGF_VEGFR = x[88];

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition
if (TGFB3_Ab<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[1,1]*(TGFB3_Ab)^control_parameter_array[1,2])/(1+control_parameter_array[1,1]*(TGFB3_Ab)^control_parameter_array[1,2]));
end

control_vector[3] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction
push!(transfer_function_vector,(control_parameter_array[2,1]*(Active_NFATc)^control_parameter_array[2,2])/(1+control_parameter_array[2,1]*(Active_NFATc)^control_parameter_array[2,2]));
# repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression
if (pSmad<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[3,1]*(pSmad)^control_parameter_array[3,2])/(1+control_parameter_array[3,1]*(pSmad)^control_parameter_array[3,2]));
end

# repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression
if (SNAIL_SLUG<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[4,1]*(SNAIL_SLUG)^control_parameter_array[4,2])/(1+control_parameter_array[4,1]*(SNAIL_SLUG)^control_parameter_array[4,2]));
end

# repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression
if (Active_LEF1<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[5,1]*(Active_LEF1)^control_parameter_array[5,2])/(1+control_parameter_array[5,1]*(Active_LEF1)^control_parameter_array[5,2]));
end

control_vector[135] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_TGFB3_by_Bcatenin target: induction_gene_TGFB3 actor: Bcatenin type: induction
push!(transfer_function_vector,(control_parameter_array[6,1]*(Bcatenin)^control_parameter_array[6,2])/(1+control_parameter_array[6,1]*(Bcatenin)^control_parameter_array[6,2]));
# activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction
push!(transfer_function_vector,(control_parameter_array[7,1]*(Active_TCF4)^control_parameter_array[7,2])/(1+control_parameter_array[7,1]*(Active_TCF4)^control_parameter_array[7,2]));
control_vector[138] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression
if (SNAIL_SLUG<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[8,1]*(SNAIL_SLUG)^control_parameter_array[8,2])/(1+control_parameter_array[8,1]*(SNAIL_SLUG)^control_parameter_array[8,2]));
end

control_vector[139] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction
push!(transfer_function_vector,(control_parameter_array[9,1]*(pSmad)^control_parameter_array[9,2])/(1+control_parameter_array[9,1]*(pSmad)^control_parameter_array[9,2]));
# repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression
if (YREG1<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[10,1]*(YREG1)^control_parameter_array[10,2])/(1+control_parameter_array[10,1]*(YREG1)^control_parameter_array[10,2]));
end
# DB constituatively inactive enforcement

if pSmad < 0.4
	mod = 0.0
else
	mod = 1.0
end
control_vector[140] = mean(transfer_function_vector*mod);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction
push!(transfer_function_vector,(control_parameter_array[11,1]*(Active_LEF1)^control_parameter_array[11,2])/(1+control_parameter_array[11,1]*(Active_LEF1)^control_parameter_array[11,2]));
# activation_vimentin_by_pSmad target: induction_gene_vimentin actor: pSmad type: induction
push!(transfer_function_vector,(control_parameter_array[12,1]*(pSmad)^control_parameter_array[12,2])/(1+control_parameter_array[12,1]*(pSmad)^control_parameter_array[12,2]));
# activation_vimentin_by_Bcatenin target: induction_gene_vimentin actor: Bcatenin type: induction
push!(transfer_function_vector,(control_parameter_array[13,1]*(Bcatenin)^control_parameter_array[13,2])/(1+control_parameter_array[13,1]*(Bcatenin)^control_parameter_array[13,2]));
# activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction
push!(transfer_function_vector,(control_parameter_array[14,1]*(AP1_SP1_P)^control_parameter_array[14,2])/(1+control_parameter_array[14,1]*(AP1_SP1_P)^control_parameter_array[14,2]));
control_vector[142] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction
push!(transfer_function_vector,(control_parameter_array[15,1]*(AP1_SP1_P)^control_parameter_array[15,2])/(1+control_parameter_array[15,1]*(AP1_SP1_P)^control_parameter_array[15,2]));
# activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction
push!(transfer_function_vector,(control_parameter_array[16,1]*(ss_virus)^control_parameter_array[16,2])/(1+control_parameter_array[16,1]*(ss_virus)^control_parameter_array[16,2]));
control_vector[152] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# Modify the rate_vector with the control variables -
rate_vector = rate_vector.*control_vector;

# Return the modified rate vector -
return rate_vector;
end


text:
# Autogenerated wrapper script for libgeotiff_jll for x86_64-unknown-freebsd
export applygeo, geotifcp, libgeotiff, listgeo, makegeo

using PROJ_jll
using Libtiff_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("libgeotiff")
JLLWrappers.@declare_library_product(libgeotiff, "libgeotiff.so.5")
JLLWrappers.@declare_executable_product(applygeo)
JLLWrappers.@declare_executable_product(geotifcp)
JLLWrappers.@declare_executable_product(listgeo)
JLLWrappers.@declare_executable_product(makegeo)
function __init__()
    JLLWrappers.@generate_init_header(PROJ_jll, Libtiff_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libgeotiff,
        "lib/libgeotiff.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        applygeo,
        "bin/applygeo",
    )

    JLLWrappers.@init_executable_product(
        geotifcp,
        "bin/geotifcp",
    )

    JLLWrappers.@init_executable_product(
        listgeo,
        "bin/listgeo",
    )

    JLLWrappers.@init_executable_product(
        makegeo,
        "bin/makegeo",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SingleFloats

export Single32

using Base: IEEEFloat
import Base: promote_rule, promote_type, convert, hash, show, string, repr, parse, tryparse,
    cmp, ==, !=, >=, >, <=, <, isless, isequal, iseven, isfinite, isnan, isodd, issubnormal,
    +, -, *, /, \, %, √, ÷, inv,
    prevfloat, nextfloat, eps, exponent, significand, zero, one, iszero, isone, isinteger, decompose,
    typemin, typemax, floatmin, floatmax, maxintfloat, widen, precision,
    ^, sqrt, cbrt, ∛, div, rem, fld, mod, cld, hypot, min, max, minmax, maximum, minimum, rad2deg, deg2rad, mod2pi, rem2pi,
    abs, abs2, signbit, sign, copysign, flipsign, floor, ceil, trunc, round, divrem, fldmod, clamp

import Base.Math: exp, exp2, exp10, expm1, log, log2, log10, log1p,
    sin, cos, tan, csc, sec, cot, asin, acos, atan, acsc, asec, acot,
    sinh, cosh, tanh, csch, sech, coth, asinh, acosh, atanh, acsch, asech, acoth

using LinearAlgebra
import LinearAlgebra: rank, cond, norm, opnorm, det, tr, inv, pinv,
                      transpose, adjoint, eigvals, eigvals!, svdvals, svdvals!, svd,
                      iszero, isone, isdiag, issymmetric, ishermitian,
                      isposdef, isposdef!, istril, istriu,
                      lu, lu!, qr, qr!,
                      schur, schur!, hessenberg, hessenberg!,
                      cholesky, cholesky!, factorize
using Random

include("type/Single32.jl")
include("math/arith.jl")
include("math/elementary.jl")
include("math/linearalgebra.jl")

end # SingleFloats


text:
import Random123
using KernelAbstractions.Extras
using FastClosures

"""
Return a tuple of 4 random UInt32
Uses the philox random CBRNG 4x32_10
"""
function randbits_4x32(key::UInt64, ctr1::UInt64, ctr2::UInt64)::NTuple{4, UInt32}
    r= Random123.Philox4x{UInt32, 10}(0,0,0,0,
                                    key%UInt32, UInt32(key>>>32),
                                    ctr1%UInt32, UInt32(ctr1>>>32),
                                    ctr2%UInt32, UInt32(ctr2>>>32), 0)
    Random123.random123_r(r)
end

"""
Return a tuple of 2 random UInt64
Uses the philox random CBRNG 4x32_10
"""
function randbits_2x64(key::UInt64, ctr1::UInt64, ctr2::UInt64)::NTuple{2, UInt64}
    x= randbits_4x32(key, ctr1, ctr2)
    (UInt64(x[1])|UInt64(x[2])<<32, UInt64(x[3])|UInt64(x[4])<<32)
end

"""
Return a tuple of 2 normal random Float64
Uses the philox random CBRNG 4x32_10 and box-muller transform
"""
function randn_2x64(key::UInt64, ctr1::UInt64, ctr2::UInt64)::NTuple{2, Float64}
    x= randbits_2x64(key, ctr1, ctr2)
    a= signed(x[1]) * 2.0^-63 # uniform [-1,1]
    b= x[2] * 2.0^-64 + 2.0^-65 # uniform (0,1]
    r= sqrt(-2*log(b))
    r .* sincos(a*pi)
end

"""
Return a tuple of 4 normal random Float32
Uses the philox random CBRNG 4x32_10 and box-muller transform
"""
function randn_4x32(key::UInt64, ctr1::UInt64, ctr2::UInt64)::NTuple{4, Float32}
    x= randbits_4x32(key, ctr1, ctr2)
    a= signed(x[1]) * 2.0f0^-31 # uniform [-1,1]
    b= x[2] * 2.0f0^-32 + 2.0f0^-33 # uniform (0,1]
    r= sqrt(-2*log(b))
    n1= r*sin(a*pi)
    n2= r*cos(a*pi)
    a= signed(x[3]) * 2.0f0^-31 # uniform [-1,1]
    b= x[4] * 2.0f0^-32 + 2.0f0^-33 # uniform (0,1]
    r= sqrt(-2*log(b))
    n3= r*sin(a*pi)
    n4= r*cos(a*pi)
    (n1,n2,n3,n4)
end

"""
Return a tuple of n normal random Float32, incrementing ctr2 for each 4.
Uses the philox random CBRNG 4x32_10 and box-muller transform
"""
function randn_32(key::UInt64, ctr1::UInt64, ctr2::UInt64, n::Val{N}) where N
    ntuple(@closure(i -> @inbounds(randn_4x32(key, ctr1, ctr2 + (i-1)÷4)[(i-1)%4+1])),n)
end


text:
function gen_window_names()
    s_set = {}
    for m in ["count", "sum", "mean", "median", "var", "std", "min", "max", "corr", "corr_pairwise", "cov", "skew", "kurt", "apply", "quantile"]
        for f in ["rolling", "expanding"]
            s = string(":", f, "_", m)
            push!(s_set, s)
        end
    end
    return join(s_set, ", ")
end


text:
using Iterators

function outer(f, vals...)
    shape = map(length, vals)
    out = Array{Any}(prod(shape)) # can I do better than Any
    for (idx, val) in enumerate(product(vals...))
        out[idx] = f(val...)
    end
    reshape(out, shape)
end

@show outer((x, y, z)->[x, y, z], [1, 2, 3], [4, 5, 6], [7, 8, 9])

# Outer[StringJoin, {"", "re", "un"}, {"cover", "draw", "wind"}, {"", "ing", "s"}]
@show outer(string, ["", "re", "un"], ["cover", "draw", "wind"], ["", "ing", "s"])

#TODO: so I have some column major issues due to how `product` iterates. I need to think
# about what this means for how I should input things

max_stars_repo_name: gabrielgellner/Scatchpad

max_issues_repo_name: gabrielgellner/Scatchpad

max_forks_repo_name: gabrielgellner/Scatchpad


text:
###############################################################################
# url file
function checkURL(features, metadata_path)
  url_path = joinpath(metadata_path, "url")
  features[:URL_EXISTS] = isfile(url_path)
  gitpath = readall(url_path)
  # Remove starting git:// or https://
  if gitpath[1:3] == "git"
    gitpath = gitpath[7:(end-5)]
  elseif gitpath[1:5] == "https"
    gitpath = gitpath[9:(end-5)]
  end
  features[:URL] = string("http://",gitpath)
end

###############################################################################
# DESCRIPTION.md file
function checkDesc(features, metadata_path)
  desc_path = joinpath(metadata_path, "DESCRIPTION.md")
  features[:DESC_EXISTS] = isfile(desc_path)
end

###############################################################################
# require files
function checkRequire(features, metadata_path)
  all_requires_ok = true
  features[:REQUIRES_FAILS] = ASCIIString[]
  features[:REQUIRES_PASSES] = ASCIIString[]

  versions_folder = joinpath(metadata_path, "versions")
  # What is the OS independent way to do this?
  version_list = split(readall(`ls $versions_folder`), "\n")
  for version in version_list[1:(end-1)]
    version_folder = joinpath(versions_folder, version)
    require_file = joinpath(version_folder, "requires")
    file_contents = ""
    try
      file_contents = readall(require_file)
    end
    passes = ismatch(r"julia", file_contents)
    if !passes
      all_requires_ok = false
      push!(features[:REQUIRES_FAILS], version)
    else
      push!(features[:REQUIRES_PASSES], version)
    end
  end
  features[:REQUIRES_OK] = all_requires_ok
  features[:VERSION] = version_list[end-1]
end


text:
using Statistics
using JLD2
using Printf
using GLMakie
using Oceananigans
using Oceananigans.Units

using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.Architectures: arch_array
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.TurbulenceClosures: HorizontallyCurvilinearAnisotropicDiffusivity, VerticallyImplicitTimeDiscretization
using CUDA: @allowscalar
using Oceananigans.Operators: Δzᵃᵃᶜ

include("cyclic_interpolate_utils.jl")

#####
##### Grid
#####

latitude = (-84.375, 84.375)
Δφ = latitude[2] - latitude[1]

# 2.8125 degree resolution
Nx = 128
Ny = 60
Nz = 18

output_prefix = "annual_cycle_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_temp"

arch = CPU()
reference_density = 1035

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

bathymetry_path = "bathy_128x60var4.bin"
east_west_stress_path = "off_TAUXvar1.bin"
north_south_stress_path = "off_TAUY.bin"
sea_surface_temperature_path="sst25_128x60x12.bin"

Nmonths = 12
bytes = sizeof(Float32) * Nx * Ny

bathymetry = reshape(bswap.(reinterpret(Float32, read(bathymetry_path, bytes))), (Nx, Ny))
τˣ = - reshape(bswap.(reinterpret(Float32, read(east_west_stress_path, Nmonths * bytes))), (Nx, Ny, Nmonths)) ./ reference_density
τʸ = - reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, Nmonths * bytes))), (Nx, Ny, Nmonths)) ./ reference_density
target_sea_surface_temperature = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, Nmonths * bytes))), (Nx, Ny, Nmonths))

bathymetry = arch_array(arch, bathymetry)
τˣ = arch_array(arch, τˣ)
τʸ = arch_array(arch, τʸ)
target_sea_surface_temperature = T★ = arch_array(arch, target_sea_surface_temperature)

H = 3600.0
# H = - minimum(bathymetry)

# Uncomment for a flat bottom:
# bathymetry = - H .* (bathymetry .< -10)

# A spherical domain
@show underlying_grid = LatitudeLongitudeGrid(arch,
                                              size = (Nx, Ny, Nz),
                                              longitude = (-180, 180),
                                              latitude = latitude,
                                              halo = (3, 3, 3),
                                              z = (-H, 0))

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

#####
##### Physics and model setup
#####

νh = 1e+5
νz = 1e+1
κh = 1e+3
κz = 1e-4

background_diffusivity =
    HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz,
                                                  time_discretization = VerticallyImplicitTimeDiscretization())

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0)

#####
##### Boundary conditions / constant-in-time surface forcing
#####

Δz_top = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.grid)
Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.grid)

@inline function surface_temperature_relaxation(i, j, grid, clock, fields, p)
    time = clock.time

    n₁ = current_time_index(time)
    n₂ = next_time_index(time)

    @inbounds begin
        T★₁ = p.T★[i, j, n₁]
        T★₂ = p.T★[i, j, n₂]
        T_surface = fields.T[i, j, grid.Nz]
    end

    T★ = cyclic_interpolate(T★₁, T★₂, time)

    return p.λ * (T_surface - T★)
end

T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation,
                                                discrete_form = true,
                                                parameters = (λ = Δz_top/3days, T★ = target_sea_surface_temperature))

@inline function wind_stress(i, j, grid, clock, fields, τ)
    time = clock.time
    n₁ = current_time_index(time)
    n₂ = next_time_index(time)

    @inbounds begin
        τ₁ = τ[i, j, n₁]
        τ₂ = τ[i, j, n₂]
    end

    return cyclic_interpolate(τ₁, τ₂, time)
end

u_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τˣ)
v_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τʸ)

@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

# Linear bottom drag:
μ = Δz_bottom / 10days

u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ)
v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ)

u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc)
v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc)
T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc)

# @inline function u_immersed_bottom_drag(i, j, k, grid, clock, fields, ν)

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = ExplicitFreeSurface(),
                                    #free_surface = ImplicitFreeSurface(maximum_iterations=10),
                                    #free_surface = ImplicitFreeSurface(),
                                    momentum_advection = VectorInvariant(),
                                    tracer_advection = WENO5(),
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs),
                                    buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=0.0)),
                                    tracers = (:T, :S),
                                    closure = (background_diffusivity, convective_adjustment))

#####
##### Initial condition:
#####

u, v, w = model.velocities
η = model.free_surface.η
T = model.tracers.T
T .= -1
S = model.tracers.S
S .= 30

#####
##### Simulation setup
#####

# Time-scale for gravity wave propagation across the smallest grid cell
g = model.free_surface.gravitational_acceleration
gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

minimum_Δx = abs(grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ[1:grid.Ny])) * deg2rad(grid.Δλ))
minimum_Δy = abs(grid.radius * deg2rad(grid.Δφ))
wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed

if model.free_surface isa ExplicitFreeSurface
    Δt = 60seconds
else
    Δt = 20minutes
end

simulation = Simulation(model, Δt = Δt, stop_time = 30year)

start_time = [time_ns()]

function progress(sim)
    wall_time = (time_ns() - start_time[1]) * 1e-9

    η = model.free_surface.η

    if model.free_surface isa ExplicitFreeSurface
        @info @sprintf("Time: % 12s, iteration: %d, max(|η|): %.2e m, wall time: %s",
                       prettytime(sim.model.clock.time),
                       sim.model.clock.iteration,
                       maximum(abs, η),
                       prettytime(wall_time))
    else
        @info @sprintf("Time: % 12s, iteration: %d, free surface iterations: %d, max(|η|): %.2e m, wall time: %s",
                       prettytime(sim.model.clock.time),
                       sim.model.clock.iteration,
                       sim.model.free_surface.implicit_step_solver.preconditioned_conjugate_gradient_solver.iteration,
                       maximum(abs, η),
                       prettytime(wall_time))
    end

    start_time[1] = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

u, v, w = model.velocities
T, S = model.tracers
η = model.free_surface.η

output_fields = (; u, v, T, S, η)
save_interval = 5days

simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, T, S, η),
                                                              schedule = TimeInterval(save_interval),
                                                              prefix = output_prefix * "_surface",
                                                              field_slicer = FieldSlicer(k=grid.Nz),
                                                              force = true)

simulation.output_writers[:bottom_fields] = JLD2OutputWriter(model, (; u, v, T, S),
                                                             schedule = TimeInterval(save_interval),
                                                             prefix = output_prefix * "_bottom",
                                                             field_slicer = FieldSlicer(k=1),
                                                             force = true)

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(1year),
                                                        prefix = output_prefix * "_checkpoint",
                                                        cleanup = true,
                                                        force = true)

# Let's goo!
@info "Running with Δt = $(prettytime(simulation.Δt))"

run!(simulation)

@info """

    Simulation took $(prettytime(simulation.run_time))
    Background diffusivity: $background_diffusivity
    Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale))
    Free surface: $(typeof(model.free_surface).name.wrapper)
    Time step: $(prettytime(Δt))
"""

#####
##### Visualize solution
#####

surface_file = jldopen(output_prefix * "_surface.jld2")
bottom_file = jldopen(output_prefix * "_bottom.jld2")

iterations = parse.(Int, keys(surface_file["timeseries/t"]))

iter = Node(0)

ηi(iter) = surface_file["timeseries/η/" * string(iter)][:, :, 1]
ui(iter) = surface_file["timeseries/u/" * string(iter)][:, :, 1]
vi(iter) = surface_file["timeseries/v/" * string(iter)][:, :, 1]
Ti(iter) = surface_file["timeseries/T/" * string(iter)][:, :, 1]
ti(iter) = string(surface_file["timeseries/t/" * string(iter)] / day)

ubi(iter) = bottom_file["timeseries/u/" * string(iter)][:, :, 1]
vbi(iter) = bottom_file["timeseries/v/" * string(iter)][:, :, 1]

η = @lift ηi($iter)
u = @lift ui($iter)
v = @lift vi($iter)
T = @lift Ti($iter)

ub = @lift ubi($iter)
vb = @lift vbi($iter)

max_η = 4
min_η = - max_η
max_u = 0.2
min_u = - max_u
max_T = 32
min_T = 0

fig = Figure(resolution = (1200, 900))

ax = Axis(fig[1, 1], title="Free surface displacement (m)")
hm = heatmap!(ax, η, colorrange=(min_η, max_η), colormap=:balance)
cb = Colorbar(fig[1, 2], hm)

ax = Axis(fig[2, 1], title="Sea surface temperature (ᵒC)")
hm = heatmap!(ax, T, colorrange=(min_T, max_T), colormap=:thermal)
cb = Colorbar(fig[2, 2], hm)

ax = Axis(fig[1, 3], title="East-west surface velocity (m s⁻¹)")
hm = heatmap!(ax, u, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[1, 4], hm)

ax = Axis(fig[2, 3], title="North-south surface velocity (m s⁻¹)")
hm = heatmap!(ax, v, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[2, 4], hm)

ax = Axis(fig[3, 1], title="East-west bottom velocity (m s⁻¹)")
hm = heatmap!(ax, ub, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[3, 2], hm)

ax = Axis(fig[3, 3], title="North-south bottom velocity (m s⁻¹)")
hm = heatmap!(ax, vb, colorrange=(min_u, max_u), colormap=:balance)
cb = Colorbar(fig[3, 4], hm)

title_str = @lift "Earth day = " * ti($iter)
ax_t = fig[0, :] = Label(fig, title_str)

GLMakie.record(fig, output_prefix * ".mp4", iterations, framerate=8) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

display(fig)

close(surface_file)


text:
using Metal
using BenchmarkTools

function kernel_fma(a, b, c, out)
    i = thread_position_in_grid_1d()
    a_val = a[i]
    b_val = b[i]
    c_val = c[i]

    for j in 1:99
        a_val = fma(a_val, b_val, c_val)
        b_val = fma(a_val, b_val, c_val)
        c_val = fma(a_val, b_val, c_val)
    end
    out[i] = fma(a_val, b_val, c_val)

    return
end
"Return calculated TFLOPS of Metal device"
function peakflops(_size=1024*1024*100)
    @info "Device: $(device().name)"
    a = round.(rand(Float32, _size) * 100)
    d_a = MtlArray(a)
    a = round.(rand(Float32, _size) * 100)
    d_b = MtlArray(a)
    a = round.(rand(Float32, _size) * 100)
    d_c = MtlArray(a)
    out = similar(a)
    d_out = MtlArray(out)

    _grid = cld(_size, 1024)

    # Warmup
    Metal.@sync @metal threads=1024 grid=_grid kernel_fma(d_a, d_b, d_c, d_out)

    bench = @benchmark Metal.@sync @metal threads=1024 grid=$_grid kernel_fma($d_a, $d_b, $d_c, $d_out)

    # Cleanup memory
    finalize(d_a)
    finalize(d_b)
    finalize(d_c)
    finalize(d_out)

    flopcount = (99*3+1) * 2 * _size
    secs = minimum(bench.times) * 1e-9
    flops = flopcount / secs

    @info "TFlops: $(round(flops/1e12; digits=2))"
    return
end

peakflops()


text:
struct WeightedUnionFind
    par::Vector{Int}
    rank::Vector{Int}
    diff_weight::Vector{Number}

    WeightedUnionFind(n) = new(collect(1:n), zeros(Int, n), zeros(Number, n))
end

function root!(wuf::WeightedUnionFind, x::Int)
    if wuf.par[x] == x
        return x
    else
        r = root!(wuf, wuf.par[x])
        wuf.diff_weight[x] += wuf.diff_weight[wuf.par[x]]
        return wuf.par[x] = r
    end
end

function weight(wuf::WeightedUnionFind, x)
    root!(wuf, x) # Calculate weight
    return wuf.diff_weight[x]
end

diff(wuf::WeightedUnionFind, x::Int, y::Int) = weight(wuf, y) - weight(wuf, x)
issame(wuf::WeightedUnionFind, x::Int, y::Int) = root!(wuf, x) == root!(wuf, y)

function union!(wuf, x::Int, y::Int, w::Number)
    w += weight(wuf, x)
    w -= weight(wuf, y)

    x = root!(wuf, x)
    y = root!(wuf, y)
    x == y && return false

    if wuf.rank[x] < wuf.rank[y]
        x, y = y, x
        w = -w
    elseif wuf.rank[x] == wuf.rank[y]
        wuf.rank[x] += 1
    end

    wuf.par[y] = x

    wuf.diff_weight[y] = w

    return true
end

function main()
    N, M = parse.(Int, split(readline()))

    wuf = WeightedUnionFind(N)
    ans = "Yes"
    for _ in 1:M
        L, R, D = parse.(Int, split(readline()))

        if issame(wuf, L, R)
            d = diff(wuf, L, R)
            if diff(wuf, L, R) ≠ D
                ans = "No"
                break
            end
        else
            union!(wuf, L, R, D)
        end
    end

    println(ans)
end

main()


text:
struct Clayton <: ArchimedeanCopula
    ϑ::Real

    function Clayton(ϑ::Real)
        # bivariate clayton can go as low as -1
        @assert -1 <= ϑ <= Inf

        if ϑ < -0.5
            @warn "Clayton returns a W copula for ϑ < -0.5"
            return W()
        elseif ϑ == 0.0
            @warn "Clayton returns an independence copula for ϑ == 0.0"
            return Independence()
        elseif ϑ == Inf
            @warn "Clayton returns an M copula for ϑ == Inf"
            return M()
        end

        return new(ϑ)
    end
end

function φ(x::Real, c::Clayton)
    return (1 + x)^(-1 / c.ϑ)
end

function φ⁻¹(x::Real, c::Clayton)
    return x^(-c.ϑ) - 1
end

function φ²(x::Real, c::Clayton)
    α = 1 / c.ϑ
    β = α * (1 + α)
    return β * (x + 1)^(-(2 + α))
end

function ∇φ⁻¹(x::Real, c::Clayton)
    return -c.ϑ * x^(-c.ϑ - 1)
end

function rosenblatt(M::AbstractMatrix, c::Clayton)
    u = 1 .+ cumsum(φ⁻¹.(M, c); dims=2)
    a = -1 / c.ϑ - 1

    return [M[:, 1] u[:, 2] .^ a ./ u[:, 1] .^ a]
end

function inverse_rosenblatt(U::AbstractMatrix, c::Clayton)
    u1 = U[:, 1]
    u2 = U[:, 2]

    return [u1 (u1 .^ (-c.ϑ) .* (u2 .^ (-c.ϑ / (c.ϑ + 1)) .- 1) .+ 1) .^ (-1 / c.ϑ)]
end

function τ(c::Clayton)
    return c.ϑ / (c.ϑ + 2)
end


text:
"""
    MatrixFDist(n1, n2, B)
```julia
n1::Real          degrees of freedom (greater than p - 1)
n2::Real          degrees of freedom (greater than p - 1)
B::AbstractPDMat  p x p scale
```
The [matrix *F*-distribution](https://projecteuclid.org/euclid.ba/1515747744)
(sometimes called the matrix beta type II distribution) generalizes the
*F*-Distribution to ``p\\times p`` real, positive definite matrices ``\\boldsymbol{\\Sigma}``.
If ``\\boldsymbol{\\Sigma}\\sim \\textrm{MF}_{p}(n_1/2,n_2/2,\\mathbf{B})``,
then its probability density function is

```math
f(\\boldsymbol{\\Sigma} ; n_1,n_2,\\mathbf{B}) =
\\frac{\\Gamma_p(\\frac{n_1+n_2}{2})}{\\Gamma_p(\\frac{n_1}{2})\\Gamma_p(\\frac{n_2}{2})}
|\\mathbf{B}|^{n_2/2}|\\boldsymbol{\\Sigma}|^{(n_1-p-1)/2}|\\mathbf{B}+\\boldsymbol{\\Sigma}|^{-(n_1+n_2)/2}.
```

If the joint distribution ``p(\\boldsymbol{\\Psi},\\boldsymbol{\\Sigma})=p(\\boldsymbol{\\Psi})p(\\boldsymbol{\\Sigma}|\\boldsymbol{\\Psi})``
is given by

```math
\\begin{aligned}
\\boldsymbol{\\Psi}&\\sim \\textrm{W}_p(n_1, \\mathbf{B})\\\\
\\boldsymbol{\\Sigma}|\\boldsymbol{\\Psi}&\\sim \\textrm{IW}_p(n_2, \\boldsymbol{\\Psi}),
\\end{aligned}
```

then the marginal distribution of ``\\boldsymbol{\\Sigma}`` is
``\\textrm{MF}_{p}(n_1/2,n_2/2,\\mathbf{B})``.
"""
struct MatrixFDist{T <: Real, TW <: Wishart} <: ContinuousMatrixDistribution
    W::TW
    n2::T
    logc0::T
end

#  -----------------------------------------------------------------------------
#  Constructors
#  -----------------------------------------------------------------------------

function MatrixFDist(n1::Real, n2::Real, B::AbstractPDMat)
    p = dim(B)
    n1 > p - 1 || throw(ArgumentError("first degrees of freedom must be larger than $(p - 1)"))
    n2 > p - 1 || throw(ArgumentError("second degrees of freedom must be larger than $(p - 1)"))
    logc0 = matrixfdist_logc0(n1, n2, B)
    T = Base.promote_eltype(n1, n2, logc0, B)
    prom_B = convert(AbstractArray{T}, B)
    W = Wishart(T(n1), prom_B)
    MatrixFDist{T, typeof(W)}(W, T(n2), T(logc0))
end

MatrixFDist(n1::Real, n2::Real, B::Union{AbstractMatrix, LinearAlgebra.Cholesky}) = MatrixFDist(n1, n2, PDMat(B))

#  -----------------------------------------------------------------------------
#  REPL display
#  -----------------------------------------------------------------------------

 show(io::IO, d::MatrixFDist) = show_multline(io, d, [(:n1, d.W.df), (:n2, d.n2), (:B, Matrix(d.W.S))])

#  -----------------------------------------------------------------------------
#  Conversion
#  -----------------------------------------------------------------------------

function convert(::Type{MatrixFDist{T}}, d::MatrixFDist) where T <: Real
    W = convert(Wishart{T}, d.W)
    MatrixFDist{T, typeof(W)}(W, T(d.n2), T(d.logc0))
end
Base.convert(::Type{MatrixFDist{T}}, d::MatrixFDist{T}) where {T<:Real} = d

function convert(::Type{MatrixFDist{T}}, W::Wishart, n2, logc0) where T <: Real
    WW = convert(Wishart{T}, W)
    MatrixFDist{T, typeof(WW)}(WW, T(n2), T(logc0))
end

#  -----------------------------------------------------------------------------
#  Properties
#  -----------------------------------------------------------------------------

dim(d::MatrixFDist) = dim(d.W)

size(d::MatrixFDist) = size(d.W)

rank(d::MatrixFDist) = dim(d)

insupport(d::MatrixFDist, Σ::AbstractMatrix) = isreal(Σ) && size(Σ) == size(d) && isposdef(Σ)

params(d::MatrixFDist) = (d.W.df, d.n2, d.W.S)

function mean(d::MatrixFDist)
    p = dim(d)
    n1, n2, B = params(d)
    n2 > p + 1 || throw(ArgumentError("mean only defined for df2 > dim + 1"))
    return (n1 / (n2 - p - 1)) * Matrix(B)
end

@inline partype(d::MatrixFDist{T}) where {T <: Real} = T

#  Konno (1988 JJSS) Corollary 2.4.i
function cov(d::MatrixFDist, i::Integer, j::Integer, k::Integer, l::Integer)
    p = dim(d)
    n1, n2, PDB = params(d)
    n2 > p + 3 || throw(ArgumentError("cov only defined for df2 > dim + 3"))
    n = n1 + n2
    B = Matrix(PDB)
    n1*(n - p - 1)*inv((n2 - p)*(n2 - p - 1)*(n2 - p - 3))*(2inv(n2 - p - 1)*B[i,j]*B[k,l] + B[j,l]*B[i,k] + B[i,l]*B[k,j])
end

function var(d::MatrixFDist, i::Integer, j::Integer)
    p = dim(d)
    n1, n2, PDB = params(d)
    n2 > p + 3 || throw(ArgumentError("var only defined for df2 > dim + 3"))
    n = n1 + n2
    B = Matrix(PDB)
    n1*(n - p - 1)*inv((n2 - p)*(n2 - p - 1)*(n2 - p - 3))*((2inv(n2 - p - 1) + 1)*B[i,j]^2 + B[j,j]*B[i,i])
end

#  -----------------------------------------------------------------------------
#  Evaluation
#  -----------------------------------------------------------------------------

function matrixfdist_logc0(n1::Real, n2::Real, B::AbstractPDMat)
    #  returns the natural log of the normalizing constant for the pdf
    p = dim(B)
    term1 = -logmvbeta(p, n1 / 2, n2 / 2)
    term2 = (n2 / 2) * logdet(B)
    return term1 + term2
end

function logkernel(d::MatrixFDist, Σ::AbstractMatrix)
    p = dim(d)
    n1, n2, B = params(d)
    ((n1 - p - 1) / 2) * logdet(Σ) - ((n1 + n2) / 2) * logdet(pdadd(Σ, B))
end

#  -----------------------------------------------------------------------------
#  Sampling
#  -----------------------------------------------------------------------------

function _rand!(rng::AbstractRNG, d::MatrixFDist, A::AbstractMatrix)
    Ψ = rand(rng, d.W)
    A .= rand(rng, InverseWishart(d.n2, Ψ) )
end

#  -----------------------------------------------------------------------------
#  Transformation
#  -----------------------------------------------------------------------------

inv(d::MatrixFDist) = ( (n1, n2, B) = params(d); MatrixFDist(n2, n1, inv(B)) )

#  -----------------------------------------------------------------------------
#  Test utils
#  -----------------------------------------------------------------------------

function _univariate(d::MatrixFDist)
    check_univariate(d)
    n1, n2, B = params(d)
    μ = zero(partype(d))
    σ = (n1 / n2) * Matrix(B)[1]
    return AffineDistribution(μ, σ, FDist(n1, n2))
end

function _rand_params(::Type{MatrixFDist}, elty, n::Int, p::Int)
    n == p || throw(ArgumentError("dims must be equal for MatrixFDist"))
    n1 = elty( n + 1 + abs(10randn()) )
    n2 = elty( n + 3 + abs(10randn()) )
    B = (X = 2rand(elty, n, n) .- 1; X * X')
    return n1, n2, B
end


text:
import Base.vec
import Base.show
import Base.isapprox
import Base.length
import Base.transpose

export VectorSpace
export I
export dot
export norm
export outer
export normalize

"""
    Represents an abstract vector space.
"""
abstract type VectorSpace end

## Generic Functions ##

vec(vs::VectorSpace) = Base.vec(vs.a)
show(io::IO, vs::VectorSpace) = Base.show(io, vec(vs))
length(vs::VectorSpace) = Base.length(vec(vs))
transpose(vs::VectorSpace) = Base.transpose(vec(vs))
I(vs::VectorSpace) = reshape([Float64(i == j) for i in 1:length(vs) for j in 1:length(vs)],
                             length(vs),
                             length(vs))
# dot product with a given metric
dot(vs1::VectorSpace, vs2::VectorSpace, M::Array{Float64}) = transpose(vs1) * M * vec(vs2)
dot(vs1::VectorSpace, vs2::VectorSpace) = dot(vs1, vs2, I(vs1))
norm(vs::VectorSpace) = sqrt(dot(vs, vs))
outer(vs1::VectorSpace, vs2::VectorSpace) = vec(vs1) * transpose(vs2)
normalize(vs::VectorSpace) = vs * (1 / norm(vs))

Base.isapprox(vs1::VectorSpace, vs2::VectorSpace; atol::Float64 = TOLERANCE) = begin
    isapprox(vec(vs1), vec(vs2), atol = atol)
end
isapprox(a₁::Array{<:VectorSpace}, a₂::Array{<:VectorSpace}; atol::Real = TOLERANCE) = begin
    all(map((i, j) -> isapprox(i, j, atol = atol), zip(a₁, a₂)))
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# import Pkg; Pkg.add("GR")

println("using GR")
@time using GR  # takes 0.2 sec!

println("randn")
@time randn(10000)

println("histogram")
@time histogram(randn(10000))  # takes 3 sec!

# without these last lines, run this script with
#   julia -i test_gr.jl
#   quit with CTRL-D

println("done. press <ENTER> to quit")
readline()

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
#EasyPlotGrace base types & core functions
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
const NCOLORS_GRACEDFLT = 15 #Number of colors expected to be defined by Grace.
#NOTE: Assuming color 0 is white/background color (not to be used)

#==Base types
===============================================================================#
#Manages additional colors (leaves default ones intact):
#NOTE: 0 appears to be bkgnd color & 1: default frame color
mutable struct ColorMgr
	plt::GracePlot.Plot #Needed to send commands
	colormap::Dict{String, Int} #Maps #RRGGBB hex color string to registered color number
	nextcolor::Int
end
ColorMgr(p::GracePlot.Plot) = ColorMgr(p, Dict(), NCOLORS_GRACEDFLT+1)

mutable struct WfrmBuilder <: EasyPlot.AbstractWfrmBuilder
	ref::GracePlot.GraphRef #Axes reference
	theme::EasyPlot.Theme
	colormgr::ColorMgr
	fold::Optional{EasyPlot.FoldedAxis}
end

#==Registering new colors
===============================================================================#
const HEX_CODES = UInt8[
	'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F'
]
function int2hexcolorstr(v::UInt32)
	result = Array{UInt8}(undef, 7)
	result[1] = '#'
	for i in length(result):-1:2
		result[i] = HEX_CODES[(v & 0xF)+1]
		v >>= 4
	end
	return String(result)
end

function buildcmd_mapcolor(id::Int, colorval::UInt32, idstr::String)
	r = (colorval>>16) & 0xFF
	g = (colorval>>8) & 0xFF
	b = colorval & 0xFF
	result = "map color $id to ($r,$g,$b), \"$idstr\""
end

#Automatically adds color if does not exist:
function _getcoloridx(mgr::ColorMgr, v::Colorant)
	v = convert(RGB24, v)
	v = v.color
	colorid = int2hexcolorstr(v)
	idx = get(mgr.colormap, colorid, missing)
	if !ismissing(idx)
		return idx
	end

	idx = mgr.nextcolor
	cmd = buildcmd_mapcolor(idx, v, colorid)
	GracePlot.sendcmd(mgr.plt, cmd)
	mgr.colormap[colorid] = idx
	mgr.nextcolor += 1
	return idx
end

#==Mapping functions
===============================================================================#
mapcolor(mgr::ColorMgr, v::Colorant) = _getcoloridx(mgr, v)
mapfacecolor(mgr::ColorMgr, v) = mapcolor(mgr, v) #In case we want to diverge

#Linewidth:
maplinewidth(w) = w
maplinewidth(::Nothing) = maplinewidth(1) #default

#Glyph size:
mapglyphsize(sz) = sz/2
mapglyphsize(::Nothing) = mapglyphsize(1) #default

#TODO: Support ColorScheme:
function _graceline(attr::EasyPlot.WfrmAttributes, mgr::ColorMgr)
	return line(style=attr.linestyle,
	            width=maplinewidth(attr.linewidth),
	            color=mapcolor(mgr, attr.linecolor),
	           )
end
function _graceglyph(attr::EasyPlot.WfrmAttributes, mgr::ColorMgr)
	nofill = (colorant"transparent" == attr.glyphfillcolor)
	glyphfillcolor = nofill ? colorant"white" : attr.glyphfillcolor
	return glyph(shape=attr.glyphshape,
	             size=mapglyphsize(attr.glyphsize),
	             linewidth=maplinewidth(attr.linewidth),
	             color=mapcolor(mgr, attr.glyphlinecolor),
	             fillcolor=mapfacecolor(mgr, glyphfillcolor),
	             fillpattern=(nofill ? 0 : 1)
	            )
end

#==AbstractWfrmBuilder implementation
===============================================================================#
EasyPlot.needsfold(b::WfrmBuilder) = b.fold

#Called by EasyPlot, for each individual DataF1 ∈ DataMD.
function EasyPlot.addwfrm(b::WfrmBuilder, d::DataF1, id::String,
	la::EasyPlot.LineAttributes, ga::EasyPlot.GlyphAttributes, strip::Int)
	attr = EasyPlot.WfrmAttributes(b.theme, la, ga) #Apply theme to attributes
	_line = _graceline(attr, b.colormgr)
	_glyph = _graceglyph(attr, b.colormgr)
	add(b.ref, d.x, d.y, _line, _glyph, id=id)
end

#Last line


text:
module HighLevelInterface

using ..SolverInterface
using ..MathProgBase

function warn_no_inf(T)
    if !(isinf(typemin(T)) && isinf(typemax(T)))
        @warn("Element type $T does not have an infinite value. Note that this may artifically introduce ranged (two-sided) constraints. To avoid this, consider casting the problem data to Float64.")
    end
end

include("linprog.jl")
include("mixintprog.jl")
include("quadprog.jl")

end


text:
module StrongStabilityPreservingRungeKuttaMethod
export StrongStabilityPreservingRungeKutta
export SSPRK33ShuOsher, SSPRK34SpiteriRuuth

using GPUifyLoops
include("StrongStabilityPreservingRungeKuttaMethod_kernels.jl")

using ..ODESolvers
ODEs = ODESolvers
using ..SpaceMethods
using ..MPIStateArrays: device, realview

"""
    StrongStabilityPreservingRungeKutta(f, RKA, RKB, RKC, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a strong-stability-preserving Runge--Kutta scheme
based on the provided `RKA`, `RKB` and `RKC` coefficient arrays.

The available concrete implementations are:

  - [`SSPRK33ShuOsher`](@ref)
  - [`SSPRK34SpiteriRuuth`](@ref)
"""
mutable struct StrongStabilityPreservingRungeKutta{T, RT, AT, Nstages} <: ODEs.AbstractODESolver
  "time step"
  dt::RT
  "time"
  t::RT
  "rhs function"
  rhs!
  "Storage for RHS during the `StrongStabilityPreservingRungeKutta` update"
  Rstage::AT
  "Storage for the stage state during the `StrongStabilityPreservingRungeKutta` update"
  Qstage::AT
  "RK coefficient vector A (rhs scaling)"
  RKA::Array{RT,2}
  "RK coefficient vector B (rhs add in scaling)"
  RKB::Array{RT,1}
  "RK coefficient vector C (time scaling)"
  RKC::Array{RT,1}

  function StrongStabilityPreservingRungeKutta(rhs!, RKA, RKB, RKC,
                                               Q::AT; dt=0, t0=0) where {AT<:AbstractArray}
    T = eltype(Q)
    RT = real(T)
    new{T, RT, AT, length(RKB)}(RT(dt), RT(t0), rhs!, similar(Q), similar(Q), RKA, RKB, RKC)
  end
end

function StrongStabilityPreservingRungeKutta(spacedisc::AbstractSpaceMethod, RKA, RKB, RKC,
                                             Q::AT; dt=0, t0=0) where {AT<:AbstractArray}
  rhs! = (x...; increment) -> SpaceMethods.odefun!(spacedisc, x..., increment = increment)
  StrongStabilityPreservingRungeKutta(rhs!, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

ODEs.updatedt!(ssp::StrongStabilityPreservingRungeKutta, dt) = (ssp.dt = dt)
ODEs.updatetime!(lsrk::StrongStabilityPreservingRungeKutta, time) = (lsrk.t = time)

"""
    ODESolvers.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, p,
                       timeend::Real, adjustfinalstep::Bool)

Use the strong stability preserving Runge--Kutta method `ssp` to step `Q`
forward in time from the current time, to the time `timeend`. If
`adjustfinalstep == true` then `dt` is adjusted so that the step does not take
the solution beyond the `timeend`.
"""
function ODEs.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, p,
                      timeend::Real, adjustfinalstep::Bool)
  time, dt = ssp.t, ssp.dt
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
  end
  @assert dt > 0

  ODEs.dostep!(Q, ssp, p, time, dt)

  if dt == ssp.dt
    ssp.t += dt
  else
    ssp.t = timeend
  end
end

"""
    ODESolvers.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, p,
                       time::Real, dt::Real, [slow_δ, slow_rv_dQ, slow_scaling])

Use the strong stability preserving Runge--Kutta method `ssp` to step `Q`
forward in time from the current time `time` to final time `time + dt`.

If the optional parameter `slow_δ !== nothing` then `slow_rv_dQ * slow_δ` is
added as an additional ODE right-hand side source. If the optional parameter
`slow_scaling !== nothing` then after the final stage update the scaling
`slow_rv_dQ *= slow_scaling` is performed.
"""
function ODEs.dostep!(Q, ssp::StrongStabilityPreservingRungeKutta, p,
                      time::Real, dt::Real, slow_δ = nothing,
                      slow_rv_dQ = nothing, in_slow_scaling = nothing)

  RKA, RKB, RKC = ssp.RKA, ssp.RKB, ssp.RKC
  rhs! = ssp.rhs!
  Rstage, Qstage = ssp.Rstage, ssp.Qstage

  rv_Q = realview(Q)
  rv_Rstage = realview(Rstage)
  rv_Qstage = realview(Qstage)

  threads = 256
  blocks = div(length(rv_Q) + threads - 1, threads)

  rv_Qstage .= rv_Q
  for s = 1:length(RKB)
    rhs!(Rstage, Qstage, p, time + RKC[s] * dt, increment = false)

    slow_scaling = nothing
    if s == length(RKB)
      slow_scaling = in_slow_scaling
    end
    @launch(device(Q), threads = threads, blocks = blocks,
            update!(rv_Rstage, rv_Q, rv_Qstage, RKA[s,1], RKA[s,2], RKB[s], dt,
                    slow_δ, slow_rv_dQ, slow_scaling))
  end
  rv_Q .= rv_Qstage
end

"""
    SSPRK33ShuOsher(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the third-order, 3-stage, strong-stability-preserving, Runge--Kutta scheme
of Shu and Osher (1988)

### References
    @article{shu1988efficient,
      title={Efficient implementation of essentially non-oscillatory shock-capturing schemes},
      author={Shu, Chi-Wang and Osher, Stanley},
      journal={Journal of computational physics},
      volume={77},
      number={2},
      pages={439--471},
      year={1988},
      publisher={Elsevier}
    }
"""
function SSPRK33ShuOsher(F, Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  T = eltype(Q)
  RT = real(T)
  RKA = [ RT(1) RT(0); RT(3//4) RT(1//4); RT(1//3) RT(2//3) ]
  RKB = [ RT(1), RT(1//4), RT(2//3) ]
  RKC = [ RT(0), RT(1), RT(1//2) ]
  StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

"""
    SSPRK34SpiteriRuuth(f, Q; dt, t0 = 0)

This function returns a [`StrongStabilityPreservingRungeKutta`](@ref) time stepping object
for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the third-order, 4-stage, strong-stability-preserving, Runge--Kutta scheme
of Spiteri and Ruuth (1988)

### References
    @article{spiteri2002new,
      title={A new class of optimal high-order strong-stability-preserving time discretization methods},
      author={Spiteri, Raymond J and Ruuth, Steven J},
      journal={SIAM Journal on Numerical Analysis},
      volume={40},
      number={2},
      pages={469--491},
      year={2002},
      publisher={SIAM}
    }
"""
function SSPRK34SpiteriRuuth(F, Q::AT; dt=0, t0=0) where {AT <: AbstractArray}
  T = eltype(Q)
  RT = real(T)
  RKA = [ RT(1) RT(0); RT(0) RT(1); RT(2//3) RT(1//3); RT(0) RT(1) ]
  RKB = [ RT(1//2); RT(1//2); RT(1//6); RT(1//2) ]
  RKC = [ RT(0); RT(1//2); RT(1); RT(1//2) ]
  StrongStabilityPreservingRungeKutta(F, RKA, RKB, RKC, Q; dt=dt, t0=t0)
end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA