text:
function _save_shapes!(shapes::Vector{<:Mesh}, dir_name)
    for i in eachindex(shapes)
        fname = lpad(i,4,"0")
        fname = "t" * fname * ".ply"
        shape = shapes[i]
        save_mesh(joinpath(dir_name, fname), shape)
    end
    0
end

function main(data_folder::String, track_id::Int, t_start::Int, t_end::Int)
    t0 = convert(Float64, t_start)
    t1 = convert(Float64, t_end)
    df = create_group_df(data_folder)
    cell = get_cell(df, track_id; t0 = t0, t1 = t1)
    coordinate = cell.coordinates
    length_step = 0.01
    add_moving_frame_data!(cell, step = length_step)
    if t0 == t1 == -1.0
        selected_df = filter(row -> row[:trackid] == track_id, df)
    else
        @assert t1 > t0 >= 1.0
        selected_df =
            filter(row -> (row.trackid == track_id) & (t0 <= row.t <= t1), df)
    end

    shape_folder = joinpath(data_folder, "hires")
    add_shape_data!(cell, selected_df, data_dir = shape_folder)
    str_t0::String = string(t0)
    str_t1::String = string(t1)
    str_id::String = string(track_id)

    dir_name = "movement_analysis" * "_" * str_id * "_" * str_t0 * "-" * str_t1
    dir_name = joinpath(data_folder, dir_name)
    if !isdir(dir_name)
        mkdir(dir_name)
    end
    movement_df = DataFrame(
        x = cell.smooth_coordinates[:, 1],
        y = cell.smooth_coordinates[:, 2],
        z = cell.smooth_coordinates[:, 3],
        time = cell.smooth_time,
        T = matrix_to_array_of_array(cell.T),
        N = matrix_to_array_of_array(cell.N),
        B = matrix_to_array_of_array(cell.B),
        curvature = cell.curvature,
        torsion = cell.torsion,
        speed = cell.speed,
    )
    fname = joinpath(dir_name, "result.csv")
    CSV.write(fname, movement_df)

    ori_dir_name =
        "original_shapes" * "_" * str_id * "_" * str_t0 * "-" * str_t1
    ori_dir_name = joinpath(data_folder, ori_dir_name)
    if !isdir(ori_dir_name)
        mkdir(ori_dir_name)
    end
    reori_dir_name =
        "reoriented_shapes" * "_" * str_id * "_" * str_t0 * "-" * str_t1
    reori_dir_name = joinpath(data_folder, reori_dir_name)
    if !isdir(reori_dir_name)
        mkdir(reori_dir_name)
    end
    _save_shapes!(cell.original_shapes, ori_dir_name)
    _save_shapes!(cell.reoriented_shapes, reori_dir_name)
    return 0
end


text:
"""
    ABINITPseudopotentialInfo

Information about a pseudopotential used in an ABINIT calculation.
"""
struct ABINITPseudopotentialInfo
    title::String           # some info about the pseudopotential
    znuclpsp::Float64       # atomic number
    zionpsp::Float64        # charge
    pspso::Int32            # spin-orbit coupling considerations?
    pspdat::Int32           # revision date
    pspcod::Int32           # might be weird if 4 or 5, otherwise ignore
    pspxc::Int32            # XC functional used with pseudopotential
    lmn_size::Int32         # Not sure what this is
    md5_pseudos::String      # Really, not as a UInt?
end

function Base.show(io::IO, ::MIME"text/plain", pspinfo::ABINITPseudopotentialInfo)
    println(io, typeof(pspinfo), ":")
    println(io, "  Title:  \"", pspinfo.title)
    println(io, "  Atom:    ", pspinfo.znuclpsp, " (", ELEMENTS[Int(pspinfo.znuclpsp)], ")")
    println(io, "  Charge:  ", pspinfo.zionpsp)
end

# Making this thing mutable is probably the best idea here
# That way we can initialize the struct and then update all the data
"""
    ABINITHeader

Header information from an abinit FORTRAN binary output file.

This data structure is mutable.
"""
Base.@kwdef mutable struct ABINITHeader
    # abinit version
    codvsn::VersionNumber = v"0.0.0"
    # Header format version
    headform::Int32 = 0
    # Unknown
    fform::Int32 = 0
    # Product of nkpt and nband
    bantot::Int32 = 0
    # Calculation date
    date::Int32 = 0
    # Grid interpolation scheme used for exchange-correlation energy
    intxc::Int32 = 0
    # Exchange-correlation functional
    ixc::Int32 = 0
    # Number of atoms
    natom::Int32 = 0
    # Size of the FFT grid
    ngfft::SVector{3,Int32} = zeros(SVector{3,Int32})
    # Number of k-points
    nkpt::Int32 = 0
    # Number of spin density components
    nspden::Int32 = 1
    # Number of spinor components
    nspinor::Int32 = 1
    # Number of spin polarizations
    nsppol::Int32 = 1
    # Number of symmetry operations
    nsym::Int32 = 1
    # Number of pseudopotentials used
    npsp::Int32 = 1
    # number of types of atoms
    ntypat::Int32 = 0
    # Method for determing band occupation
    occopt::Int32 = 0
    # Probably set to 1 if perturbative DFT is used
    pertcase::Int32 = 0
    # 1 if PAW is used
    usepaw::Int32 = 0
    # Energy cutoff (in Hartree)
    ecut::Float64 = 0
    # Energy cutoff for the second grid (used in PAW)
    ecutdg::Float64 = 0
    # Smearing of the energy cutoff (used in optimizations)
    ecutsm::Float64 = 0
    # ecut*dilatmx^2 (probably only set or different from ecut in optimizations)
    ecut_eff::Float64 = 0
    # q-points (needed for phonon calculations only)
    qptn::SVector{3,Float64} = zeros(SVector{3,Int32})
    # Translation vectors for the primitive cell
    # TODO: what exactly are the units for this thing?
    rprimd::SMatrix{3,3,Float64,9} = zeros(SMatrix{3,3,Float64})
    # Bias voltage used in electron microscopy
    stmbias::Float64 = 0
    # Physical temperature of electrons (not sure of the units)
    tphysel::Float64 = 0
    # Smearing temperature for band occupation
    tsmear::Float64 = 0
    # Use of a wavelet basis set (we probably never use this)
    usewvl::Int32 = 0
    # Unknown for now
    nshiftk_orig::Int32 = 0                             # New for abinit 8.10
    # Number of k-point grid shifts
    nshiftk::Int32 = 0                                  # New for abinit 8.10
    # Maximum number of bands
    mband::Int32 = 0                                    # New for abinit 8.10
    # Set to 1 if the wavefunction should be stored at that k-point
    istwfk::Vector{Int32} = Int32[]    # size nkpt
    # Number of bands at each k-point
    nband::Vector{Int32} = Int32[]    # size nkpt
    # Number of planewave arrays at each k-point
    npwarr::Vector{Int32} = Int32[]    # size nkpt
    # Whether spin-orbit coupling is considered for each pseudopotential
    so_psp::Vector{Int32} = Int32[]    # size npsp
    # Antiferromagnetic symmetry components (-1 for those that reverse magnetization)
    symafm::Vector{Int32} = Int32[]    # size nsym
    # Symmetry operation matrices
    symrel::Vector{SMatrix{3,3,Int32,9}} = SMatrix{3,3,Int32}[]        # size nsym
    # Atom types (each integer corresponds to znucltypat)
    typat::Vector{Int32} = Int32[]    # size natom
    # List of k-points
    kpt::Vector{SVector{3,Float64}} = SVector{3,Float64}[]
    # Occupancies for each band in each k-point
    # TODO: changes this to a matrix or a vector of vectors
    occ::Vector{Float64} = Float64[]    # size bantot (nband*nkpt)
    # "Translation non-symmorphic vectors"
    tnons::Vector{SVector{3,Float64}} = SVector{3,Float64}[]    # size nsym
    # Atomic numbers for each nucleus
    znucltypat::Vector{Float64} = Float64[]    # size ntypat
    # k-point weights
    wtk::Vector{Float64} = Float64[]    # size nkpt
    # Setup method for k-point list
    kptopt::Int32 = 1
    # 2 if complex rho(i,j) occupancies are used (only for PAW calculations)
    pawcpxocc::Int32 = 1
    # Number of electrons in the cell
    nelect::Int32 = 0
    # Charge of the full cell
    cellcharge::Int32 = 0
    # 1 if wavelet solver is used/non-standard boundary conditions
    icoulomb::Int32 = 0
    # Lattice describing the k-point grid (ngkpt is along diagonal)
    kptrlatt::SMatrix{3,3,Int32,9} = zeros(SMatrix{3,3,Int32})
    # Unknown
    kptrlatt_orig::SMatrix{3,3,Int32,9} = zeros(SMatrix{3,3,Int32})
    # Unknown
    shiftk_orig::SVector{3,Float64} = zeros(SVector{3,Float64})
    # Shift of k-point grid off origin
    shiftk::SVector{3,Float64} = zeros(SVector{3,Float64})
    # Pseudopotential info
    pspinfo::Vector{ABINITPseudopotentialInfo} = ABINITPseudopotentialInfo[]    # size npsp
    # Residual density (a metric for calculation quality)
    residm::Float64 = 0
    # Positions of atoms (types given by typat matrix)
    xred::Vector{SVector{3,Float64}} = SVector{3,Float64}[]        # size natom
    # Total energy as determined by the calculation
    etotal::Float64 = 0
    # Fermi energy as determined by the calculation
    fermie::Float64 = 0
    # Masses of each type of atom (in amu)
    amu::Vector{Float64} = Float64[]
    # 2 if complex density information is used in PAW calculations
    cplex::Int32 = 1
    # rho(i,j) matrices for each atom
    pawdata::Matrix{Matrix{Float64}} = Matrix{Matrix{Float64}}(undef, 0, 0)    # size nspden*natom
end

# Index notation for this thing, just in case that's easier
Base.getindex(h::ABINITHeader, name::Symbol) = getfield(h, name)
Base.setindex!(x, h::ABINITHeader, name::Symbol) = setfield!(x, h, name)

"""
    symrel_to_sg(symrel::AbstractVector{<:AbstractMatrix{<:Integer}}) -> Int

Converts a list of ABINIT symmetry operations to the corresponding space group number and setting.
"""
function symrel_to_sg(symrel::AbstractVector{<:AbstractMatrix{<:Integer}})
    # any parity inverting operations?
    has_mir = !any(x -> x < 0, det.(symrel))
    # inversion center?
    has_inv = diagm([-1,-1,-1]) in symrel
    # All twofold axes
    twofold = [diagm([-1, -1, -1])]
    return 1 # for now
end

symrel_to_sg(h::ABINITHeader) = symrel_to_sg(h.symrel)

function Crystal{3}(h::ABINITHeader; convert=:P)
    # Lattice vectors converted to angstroms
    latt = BasisVectors{3}(BOHR2ANG*h.rprimd)
    atomlist = AtomList{3}(
        latt,
        AtomPosition{3}.(
            Int.(h.znucltypat[h.typat]),
            h.xred
        )
    )
    return Crystal{3}(
        # ABINIT seems to use bohr by default
        RealLattice{3}(BOHR2ANG*h.rprimd, prim=true, ctr=convert),
        symrel_to_sg(h),        # TODO: get space group from header?
        [0, 0, 0],              # TODO: get origin setting from header?
        atomlist,
        atomlist,
    )
end

# This function gets indices for the rhoij entries in a PAW calculation
function triang_index(n)
    row = floor(Int,(1 + sqrt(8*n - 7))/2)
    col = row + n - Int(row*(row + 1)/2)
    return (row, col)
end

"""
    get_abinit_version(io::IO) -> NamedTuple{}

Gets the ABINIT version information from a calculation output header.
"""
function get_abinit_version(io::IO)
    # Get the size of the first data entry
    sz = read(io, Int32)
    codvsn = VersionNumber(String(read(io, sz - 8)))
    headform = read(io, Int32)
    fform = read(io, Int32)
    # Skip to the next field - assume the next thing to happen is one of the header reading methods
    skip(io, 4)
    return (codvsn=codvsn, headform=headform, fform=fform)
end

"""
    read_abinit_header_57(io::IO) -> ABINITHeader

Reads in an abinit header from the outputs of calculations made by versions up to 7.10. These files
will contain a `headform` value of 57.
"""
function read_abinit_header_57(io::IO)
    # Variable names in this function correspond to names given here:
    # https://docs.abinit.org/guide/abinit/#header
    # Generate the new header
    h = ABINITHeader()
    # Skip the header for the first entry
    skip(io, 4)
    # Get a bunch of data now
    (
        h.bantot,     # total number of bands used in calculation (nband*nkpt*nsppol)
        h.date,       # date of calculation
        h.intxc,      # grid interpolation for XC
        h.ixc,        # XC functional index (usually 1 for Teter93 LDA)
        h.natom,      # number of atoms in unit cell
        ngfft1, ngfft2, ngfft3,    # FFT grid dimensions
        h.nkpt,       # number of k-points
        h.nspden,     # 1 for no spinors, 2 for scalar spin, 4 for vector spin
        h.nspinor,    # 1 for scalar wavefunction, 2 for spinor wavefunction
        h.nsppol,     # 1 for unpolarized, 2 for polarized
        h.nsym,       # number of symmetry operations
        h.npsp,       # number of pseudopotentials read
        h.ntypat,     # number of types of atoms
        h.occopt,     # occupation setting - usually 3
        h.pertcase,   # something about perturbations - probably not important
        h.usepaw      # 1 if the projector augmented wave method is used
    ) = tuple((read(io, Int32) for n = 1:18)...)
    # Reinterpret the FFT grid
    h.ngfft = SVector{3,Int}(ngfft1, ngfft2, ngfft3)
    # Get more data, but this time it's Float64s
    (
        h.ecut,       # Energy cutoff
        h.ecutdg,     # Energy cutoff for the second grid (used in PAW)
        h.ecutsm,     # Smearing of ecut for unit cell changes in optimization
        h.ecut_eff,   # ecut*dilatmx^2
    ) = tuple((read(io, Float64) for n in 1:4)...)
    # q-points (needed for phonon calculations only)
    h.qptn = SVector{3,Float64}((read(io, Float64) for n in 1:3)...)
    # primitive cell vectors in terms of acell
    # TODO: can we assume that the units are in Å for rprimd vectors?
    h.rprimd = SMatrix{3,3,Float64}((read(io, Float64) for n in 1:9)...)
    (
        h.stmbias,    # STM bias voltage: probably not used
        h.tphysel,    # Physical electron temperature
        h.tsmear      # Band smearing factor
    ) = tuple((read(io, Float64) for n in 1:3)...)
    h.usewvl = read(io, Int32)
    # next entry
    skip(io, 8)
    # storage of wavefunction by k-points
    h.istwfk = [read(io, Int32) for n in 1:h.nkpt]
    # number of bands (these elements should probably all be the same)
    h.nband = [read(io, Int32) for n in 1:(h.nkpt*h.nsppol)]
    # number of planewave arrays (or size ?)
    h.npwarr = [read(io, Int32) for n in 1:h.nkpt]
    # spin-orbit treatment of pseudopotentials (1 uses the psp file info)
    h.so_psp = [read(io, Int32) for n in 1:h.npsp]
    # symmetries for ferromagnetic materials (-1 if magnetization is altered)
    h.symafm = [read(io, Int32) for n in 1:h.nsym]
    # matrices for symmetry operations
    h.symrel = [SMatrix{3,3,Int32}(read(io, Int32) for n in 1:9) for n in 1:h.nsym]
    # types of atoms
    h.typat = [read(io, Int32) for n in 1:h.natom]
    # k-point coordinates
    h.kpt = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.nkpt]
    # Band occupancy
    h.occ = [read(io, Float64) for n in 1:h.bantot]
    # "translation non-symmorphic vectors"
    h.tnons = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.nsym]
    # Atomic numbers for each typat entry
    h.znucltypat = [read(io, Float64) for n in 1:h.ntypat]
    # k-point weights
    h.wtk = [read(io, Float64) for n in 1:h.nkpt]
    # Get all the pseudopotential info
    title = Vector{String}(undef, h.npsp)         # probably can be ignored
    znuclpsp = Vector{Float64}(undef, h.npsp)     # atomic number
    zionpsp = Vector{Float64}(undef, h.npsp)      # charge
    pspso = Vector{Int32}(undef, h.npsp)          # spin-orbit coupling considerations?
    pspdat = Vector{Int32}(undef, h.npsp)         # revision date
    pspcod = Vector{Int32}(undef, h.npsp)         # might be weird if 4 or 5, otherwise ignore
    pspxc = Vector{Int32}(undef, h.npsp)          # XC functional used with pseudopotential
                                                # probably matches ixc
    lmn_size = Vector{Int32}(undef, h.npsp)       # size of spherical harmonic components? Idk
    for m in 1:h.npsp
        skip(io, 8)
        title[m] = join(read(io, Char) for n in 1:132)
        (
            znuclpsp[m],
            zionpsp[m]
        ) = tuple((read(io, Float64) for n in 1:2)...)
        (
            pspso[m],
            pspdat[m],
            pspcod[m],
            pspxc[m],
            lmn_size[m]
        ) = tuple((read(io, Int32) for n in 1:5)...)
    end
    h.pspinfo = [
        ABINITPseudopotentialInfo(
            title[n],
            znuclpsp[n],
            zionpsp[n],
            pspso[n],
            pspdat[n],
            pspcod[n],
            pspxc[n],
            lmn_size[n],
            ""
        ) for n in 1:h.npsp
    ]
    # One more field
    skip(io, 8)
    # Residual density (gives the calculation quality)
    h.residm = read(io, Float64)
    # Atomic coordinates
    h.xred = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.natom]
    # Total energy and Fermi energy
    (h.etotal, h.fermie) = tuple((read(io, Float64) for n in 1:2)...)
    # Vector containing PAW data
    h.pawdata = [zeros(Float64, 8, 8) for s in 1:h.nspden, a in 1:h.natom]
    # Fields for a PAW calculation
    if h.usepaw == 1
        # First field
        skip(io, 8)
        # Number of entries in the rho_ij matrix (which is sparse and symmetric)
        nrhoijsel = [read(io, Int32) for s in 1:h.nspden, a in 1:h.natom]
        # Complex components of the density
        h.cplex = read(io, Int32)
        nsp_new = read(io, Int32)
        # Second field
        skip(io, 8)
        # Non-zero values of the rho_ij matrix
        rhoijselect = [
            [read(io, Int32) for n in 1:nrhoijsel[s,a]] for s in 1:h.nspden, a in 1:h.natom
        ]
        # Convert to matrix indices
        rhoijinds = [triang_index.(v) for v in rhoijselect]
        # Get the matrix elementwise
        # TODO: assuming that rhoijp consists of 8×8 matrices, but does it always?
        # rhoij = [zeros(Float64, 8, 8) for s in 1:nspden, a in 1:natom]
        for s in 1:h.nspden, a in 1:h.natom
            for ind in rhoijinds[s,a]
                x = read(io, Float64)
                h.pawdata[s,a][ind...] = x
                h.pawdata[s,a][reverse(ind)...] = x
            end
        end
    end
    # Skip the last record entry
    skip(io, 4)
    # Return output variables
    # TODO: perhaps try to calculate nelect assuming a neutral crystal?
    return h
end

"""
    read_abinit_header_80(io::IO) -> ABINITHeader

Reads in an abinit header from the outputs of calculations made by versions up to 8.10. These files
will contain a `headform` value of 80.
"""
function read_abinit_header_80(io::IO)
    # Variable names in this function correspond to names given here:
    # https://docs.abinit.org/guide/abinit/#header
    # Generate the new header
    h = ABINITHeader()
    # Records start and end with the same Float32 containing the entry length
    skip(io, 4)                                           # seek to 26
    # Get a bunch of data now
    (
        h.bantot,     # total number of bands used in calculation (nband*nkpt*nsppol)
        h.date,       # date of calculation
        h.intxc,      # grid interpolation for XC
        h.ixc,        # XC functional index (usually 1 for Teter93 LDA)
        h.natom,      # number of atoms in unit cell
        ngfft1, ngfft2, ngfft3,    # FFT grid dimensions
        h.nkpt,       # number of k-points
        h.nspden,     # 1 for no spinors, 2 for scalar spin, 4 for vector spin
        h.nspinor,    # 1 for scalar wavefunction, 2 for spinor wavefunction
        h.nsppol,     # 1 for unpolarized, 2 for polarized
        h.nsym,       # number of symmetry operations
        h.npsp,       # number of pseudopotentials read
        h.ntypat,     # number of types of atoms
        h.occopt,     # occupation setting - usually 3
        h.pertcase,   # something about perturbations - probably not important
        h.usepaw      # 1 if the projector augmented wave method is used
    ) = tuple((read(io, Int32) for n = 1:18)...)
    # Reinterpret the FFT grid
    h.ngfft = SVector{3,Int}(ngfft1, ngfft2, ngfft3)
    # Get more data, but this time it's Float64s
    (
        h.ecut,       # Energy cutoff
        h.ecutdg,     # Energy cutoff for the second grid (used in PAW)
        h.ecutsm,     # Smearing of ecut for unit cell changes in optimization
        h.ecut_eff,   # ecut*dilatmx^2
    ) = tuple((read(io, Float64) for n in 1:4)...)
    # q-points (needed for phonon calculations only)
    h.qptn = SVector{3,Float64}((read(io, Float64) for n in 1:3)...)
    # primitive cell vectors in terms of acell
    # TODO: can we assume that the units are in Å for rprimd vectors?
    h.rprimd = SMatrix{3,3,Float64}((read(io, Float64) for n in 1:9)...)
    (
        h.stmbias,    # STM bias voltage: probably not used
        h.tphysel,    # Physical electron temperature
        h.tsmear      # Band smearing factor
    ) = tuple((read(io, Float64) for n in 1:3)...)
    (
        h.usewvl,
        h.nshiftk_orig,
        h.nshiftk,
        h.mband
    ) = tuple((read(io, Int32) for n in 1:4)...)
    # next entry
    skip(io, 8)
    # storage of wavefunction by k-points
    h.istwfk = [read(io, Int32) for n in 1:h.nkpt]
    # number of bands (these elements should probably all be the same)
    h.nband = [read(io, Int32) for n in 1:(h.nkpt*h.nsppol)]
    # number of planewave arrays (or size ?)
    h.npwarr = [read(io, Int32) for n in 1:h.nkpt]
    # spin-orbit treatment of pseudopotentials (1 uses the psp file info)
    h.so_psp = [read(io, Int32) for n in 1:h.npsp]
    # symmetries for ferromagnetic materials (-1 if magnetization is altered)
    h.symafm = [read(io, Int32) for n in 1:h.nsym]
    # matrices for symmetry operations
    h.symrel = [SMatrix{3,3,Int32}(read(io, Int32) for n in 1:9) for n in 1:h.nsym]
    # types of atoms
    h.typat = [read(io, Int32) for n in 1:h.natom]
    # k-point coordinates
    h.kpt = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.nkpt]
    # Band occupancy
    h.occ = [read(io, Float64) for n in 1:h.bantot]
    # "translation non-symmorphic vectors"
    h.tnons = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.nsym]
    # Atomic numbers for each typat entry
    h.znucltypat = [read(io, Float64) for n in 1:h.ntypat]
    # k-point weights
    h.wtk = [read(io, Float64) for n in 1:h.nkpt]
    # This is new to the version 80 header
    @debug string("position before new v80 data: ", position(io))
    # One more field
    skip(io, 8)
    # Residual density (gives the calculation quality)
    h.residm = read(io, Float64)
    # Atomic coordinates
    h.xred = [SVector{3,Float64}(read(io, Float64) for n in 1:3) for n in 1:h.natom]
    # Total energy and Fermi energy
    (h.etotal, h.fermie) = tuple((read(io, Float64) for n in 1:2)...)
    # Masses of atoms used
    h.amu = [read(io, Float64) for n in 1:h.ntypat]
    @debug string("position after amu: ", position(io))
    # This data is changed from version 57
    skip(io, 8)
    h.kptopt = read(io, Int32)
    h.pawcpxocc = read(io, Int32)
    h.nelect = read(io, Float64)
    h.cellcharge = read(io, Float64)
    h.icoulomb = read(io, Int32)
    # A lot of useful k-point info!
    @debug string("position before kptrlatt: ", position(io))
    h.kptrlatt = SMatrix{3,3,Int32}((read(io, Int32) for n in 1:9)...)
    h.kptrlatt_orig = SMatrix{3,3,Int32}((read(io, Int32) for n in 1:9)...)
    h.shiftk_orig = SVector{3,Float64}((read(io, Float64) for n in 1:3)...)
    h.shiftk = SVector{3,Float64}((read(io, Float64) for n in 1:3)...)
    # Get all the pseudopotential info
    title = Vector{String}(undef, h.npsp)           # probably can be ignored
    znuclpsp = Vector{Float64}(undef, h.npsp)       # atomic number
    zionpsp = Vector{Float64}(undef, h.npsp)        # charge
    pspso = Vector{Int32}(undef, h.npsp)            # spin-orbit coupling considerations?
    pspdat = Vector{Int32}(undef, h.npsp)           # revision date
    pspcod = Vector{Int32}(undef, h.npsp)           # might be weird if 4 or 5, otherwise ignore
    pspxc = Vector{Int32}(undef, h.npsp)            # XC functional used with pseudopotential
                                                    # probably matches ixc
    lmn_size = Vector{Int32}(undef, h.npsp)         # size of spherical harmonic components? Idk
    md5_pseudos = Vector{String}(undef, h.npsp)     # MD5 hash of pseudopotential used
    @debug string("before getting psp data: file pointer is at ", position(io))
    for m in 1:h.npsp
        skip(io, 8)
        @debug string(
            "Character in first field (should be H, probably):\n",
            repr(MIME("text/plain"), peek(io, Char))
        )
        title[m] = join(read(io, Char) for n in 1:132)
        (
            znuclpsp[m],
            zionpsp[m]
        ) = tuple((read(io, Float64) for n in 1:2)...)
        (
            pspso[m],
            pspdat[m],
            pspcod[m],
            pspxc[m],
            lmn_size[m]
        ) = tuple((read(io, Int32) for n in 1:5)...)
        md5_pseudos[m] = String(read(io, 32))
    end
    @debug string("after getting psp data: file pointer is at ", position(io))
    h.pspinfo = [
        ABINITPseudopotentialInfo(
            title[n],
            znuclpsp[n],
            zionpsp[n],
            pspso[n],
            pspdat[n],
            pspcod[n],
            pspxc[n],
            lmn_size[n],
            md5_pseudos[n]
        ) for n in 1:h.npsp
    ]
    # Vector containing PAW data
    h.pawdata = [zeros(Float64, 8, 8) for s in 1:h.nspden, a in 1:h.natom]
    # Fields for a PAW calculation
    if h.usepaw == 1
        # First field
        skip(io, 8)
        # Number of entries in the rho_ij matrix (which is sparse and symmetric)
        nrhoijsel = [read(io, Int32) for s in 1:h.nspden, a in 1:h.natom]
        # Complex components of the density
        h.cplex = read(io, Int32)
        nsp_new = read(io, Int32)
        # Second field
        skip(io, 8)
        # Non-zero values of the rho_ij matrix
        rhoijselect = [
            [read(io, Int32) for n in 1:nrhoijsel[s,a]] for s in 1:h.nspden, a in 1:h.natom
        ]
        # Convert to matrix indices
        rhoijinds = [triang_index.(v) for v in rhoijselect]
        # Get the matrix elementwise
        # TODO: assuming that rhoijp consists of 8×8 matrices, but does it always?
        # rhoij = [zeros(Float64, 8, 8) for s in 1:nspden, a in 1:natom]
        for s in 1:h.nspden, a in 1:h.natom
            for ind in rhoijinds[s,a]
                x = read(io, Float64)
                h.pawdata[s,a][ind...] = x
                h.pawdata[s,a][reverse(ind)...] = x
            end
        end
    end
    # Skip the last record entry
    skip(io, 4)
    # Return output variables
    return h
end

"""
    read_abinit_header(io::IO)

Reads the header of an ABINIT output file, automatically determining the format of the header
from the first few digits.
"""
function read_abinit_header(io::IO)
    # Select which function to use based on the headform value
    fdict = Dict{Int, Function}(
        57 => read_abinit_header_57,
        80 => read_abinit_header_80,
    )
    # Get the info stored in the header
    (codvsn, headform, fform) = get_abinit_version(io)
    # @debug string("abinit version ", codvsn, " (header version ", headform, ")")
    # Read the rest of the header
    # Select the appropriate function by version number
    # header::ABINITHeader = fdict[headform](io)::ABINITHeader
    header = if headform == 57
        read_abinit_header_57(io)
    elseif headform == 80
        read_abinit_header_80(io)
    else
        error("Unsupported header format.")
    end
    # Add in the rest of the data and return
    header.codvsn = codvsn
    header.headform = headform
    header.fform = fform
    return header
end

"""
    function read_abinit_datagrids(T, io, nspden, ngfft) -> Vector{Matrix{T}}

Reads the datagrid portion of an abinit density output, following the header.

Electron density values are given in electrons/Bohr³ in the abinit output files. These are
converted to electrons/Å³ in the resulting grids.

Depending on the value of `cplex` from the header, the resulting matrix may either contain
`Float64` or `Complex{Float64}` data.
"""
function read_abinit_datagrids(
    ::Type{T},
    io::IO,
    nspden::Integer,
    ngfft::AbstractVector{<:Integer};
    conversion::Real = 1
) where T
    # Density matrices
    rho = [Array{T,3}(undef, ngfft[1], ngfft[2], ngfft[3]) for s in 1:nspden]
    # Assuming that every grid is stored in its own
    for s in 1:nspden
        @debug "Before assertion: file pointer is at " * string(position(io))
        sz1 = sizeof(T) * prod(ngfft)
        sz2 = read(io, Int32)
        @assert sz1 == sz2 string(
            "The reported FORTRAN formatted data entry doesn't match the grid size.\n",
            "Expected to read a datagrid of type ", T,
            " and size ", join(string.(ngfft), "×"),
            " ( total size ", sz1 ,")\n",
            "Reported size from file: ", sz2
        )
        @debug "After assertion: file pointer is at " * string(position(io))
        for a in 1:ngfft[1], b in 1:ngfft[2], c in 1:ngfft[3]
        # Convert to electrons/Ang^3
            rho[s][a,b,c] = read(io, T) * conversion
        end
        @debug "Read a datagrid: file pointer is at " * string(position(io))
        skip(io, 4)
    end
    return rho
end

function read_abinit_density(io::IO)
    # Get the header from the file
    header = read_abinit_header(io)
    # Get the type of the datagrid (set by cplex in the header)
    T = (Float64, Complex{Float64})[header.cplex]
    rho = read_abinit_datagrids(T, io, header.nspden, header.ngfft, conversion=BOHR2ANG^3)
    # Add each dataset to the dictionary
    data = Dict{String, RealSpaceDataGrid{3,T}}()
    # Convert the basis
    basis = BasisVectors(BOHR2ANG * header.rprimd)
    # Fill the dictionary
    data["density_total"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    if header.nspden == 2
        data["density_spinup"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    elseif header.nspden == 4
        data["density_spinup_x"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["density_spinup_y"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["density_spinup_z"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    end
    return CrystalWithDatasets{3,String,RealSpaceDataGrid{3,T}}(
        Crystal{3}(header),
        data
    )
end

read_abinit_density(filename::AbstractString) = open(read_abinit_density, filename)

function read_abinit_potential(io::IO)
    # Get the header from the file
    header = read_abinit_header(io)
    # Get the type of the datagrid (set by cplex in the header)
    T = (Float64, Complex{Float64})[header.cplex]
    # No conversion will occur here: assume units of Hartree
    rho = read_abinit_datagrids(T, io, header.nspden, header.ngfft)
    data = Dict{String, RealSpaceDataGrid{3,T}}()
    basis = BasisVectors{3}(BOHR2ANG * header.rprimd)
    if header.nspden == 1
        data["potential_total"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    elseif header.nspden == 2
        data["potential_spinup"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["potential_spindown"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    elseif header.nspden == 4
        data["potential_up_up"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["potential_down_down"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["potential_up_down_real"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
        data["potential_up_down_imag"] = RealSpaceDataGrid(basis, [0, 0, 0], rho[1])
    end
    return CrystalWithDatasets{3,String,RealSpaceDataGrid{3,T}}(
        Crystal{3}(header),
        data
    )
end

read_abinit_potential(filename::AbstractString) = open(read_abinit_potential, filename)

"""
    readHGH(io::IO)

Reads information for an abinit-formatted HGH pseudopotential. Pseudopotential information may be
downloaded from here: https://www.abinit.org/sites/default/files/PrevAtomicData/psp-links/hgh.html
"""
function readHGH(io::IO)
    # Function to generate the off-diagonal h and k coefficients
    function gen_h_k(m::AbstractMatrix, lmax::Integer)
        h = cat([diagm(m[l,:]) for l in 1:(lmax+1)]..., dims=3)
        for l = 0:lmax
            # TODO: there must be a better way to do this
            h[1,2,l+1] = -1/2*sqrt((2*l+3)/(2*l+5)) * h[2,2,l+1]
            h[1,3,l+1] = 1/2*sqrt((2*l+3)*(2*l+5)/((2*l+7)*(2*l+9))) * h[3,3,l+1]
            h[2,3,l+1] = -sqrt((2*l+5)^2/((2*l+7)*(2*l+9))) * h[3,3,l+1]
            h[2,1,l+1] = h[1,2,l+1]
            h[3,1,l+1] = h[1,3,l+1]
            h[3,2,l+1] = h[2,3,l+1]
        end
        return h
    end
    # Skip the first line (just has a line of text)
    ln = (readline(io); readline(io))
    # Get zatom and zion
    (zatom, zion) = parse.(Int8, split(ln)[1:2])
    # Get the exchange-correlation functional used and maximum l
    (pspcod, pspxc, lmax) = parse.(Int8, split(readline(io))[1:3])
    # pspcod should *always* equal 3 for valid HGH pseudopotentials
    if pspcod != 3
        error("pspcod was not 3")
    end
    (rloc, c1, c2, c3, c4) = parse.(Float64, split(readline(io))[1:5])
    c = SVector{4,Float64}(c1, c2, c3, c4)
    r = zeros(Float64, lmax + 1)
    # Generate matrices to store the diagonal h and k components for each l
    hdiag = zeros(Float64, lmax + 1, 3)
    kdiag = zeros(Float64, lmax + 1, 3)
    # Get the s components (which have no corresponding spin-orbit parameters)
    data = parse.(Float64, split(readline(io))[1:4])
    r[1] = data[1]
    # Assign h11s, h22s, h33s
    for n in 1:3
        hdiag[1,n] = data[n+1]
    end
    # Assign the diagonal h and k components
    for l in 1:lmax
        data = parse.(Float64, split(readline(io))[1:4])
        r[l+1] = data[1]
        for n in 1:3
            hdiag[l+1,n] = data[n+1]
        end
        data = parse.(Float64, split(readline(io))[1:3])
        for n in 1:3
            kdiag[l+1,n] = data[n]
        end
    end
    # Calculate the off-diagonal components
    h = gen_h_k(hdiag, lmax)
    k = gen_h_k(kdiag, lmax)
    return HGHPseudopotential(zatom, zion, XCFunctional(pspxc), lmax, rloc, c, r, h, k)
end

readHGH(filename::AbstractString) = open(readHGH, filename)

ext: jl

max_stars_repo_name: brainandforce/Xtal

max_issues_repo_name: brainandforce/Xtal

max_forks_repo_name: brainandforce/Xtal


text:
"""
SPIR-V extension which may have been promoted to a core version or be enabled implicitly by enabled Vulkan extensions.
"""
struct SpecExtensionSPIRV
    "Name of the SPIR-V extension."
    name::String
    "Core version of the Vulkan API in which the extension was promoted, if promoted."
    promoted_to::Optional{VersionNumber}
    "Vulkan extensions that implicitly enable the SPIR-V extension."
    enabling_extensions::Vector{String}
end

"""
Condition that a feature needs to satisfy to be considered enabled.
"""
struct FeatureCondition
    "Name of the feature structure relevant to the condition."
    type::Symbol
    "Member of the structure which must be set to true to enable the feature."
    member::Symbol
    "Core version corresponding to the structure, if any."
    core_version::Optional{VersionNumber}
    "Extension required for the corresponding structure, if any."
    extension::Optional{String}
end

"""
Device property that enables a SPIR-V capability when supported.
"""
struct PropertyCondition
    "Name of the property structure relevant to the condition."
    type::Symbol
    "Member of the property structure to be tested."
    member::Symbol
    "Required core version of the Vulkan API, if any."
    core_version::Optional{VersionNumber}
    "Required extension, if any."
    extension::Optional{String}
    "Whether the property to test is a boolean. If not, then it will be a bit from a bitmask."
    is_bool::Bool
    "Name of the bit enum that must be included in the property, if the property is not a boolean."
    bit::Optional{Symbol}
end

"""
SPIR-V capability with information regarding various requirements to consider it enabled.
"""
struct SpecCapabilitySPIRV
    "Name of the SPIR-V capability."
    name::Symbol
    "Core version of the Vulkan API in which the SPIR-V capability was promoted, if promoted."
    promoted_to::Optional{VersionNumber}
    "Vulkan extensions that implicitly enable the SPIR-V capability."
    enabling_extensions::Vector{String}
    "Vulkan features that implicitly enable the SPIR-V capability."
    enabling_features::Vector{FeatureCondition}
    "Vulkan properties that implicitly enable the SPIR-V capability."
    enabling_properties::Vector{PropertyCondition}
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Variants
using Base.Test

# Constructors

x1 = Variant{Int,Int}(1, Val{true})
x2 = Variant{Int,Int}(2, true)
x3 = Variant{Int,Int}(3, Val{false})
x4 = Variant{Int,Int}(4, false)
@test_throws VariantException Variant{Int,Int}(5)
x6 = Variant{Int,Char}(6)
x7 = Variant{Int,Char}(7, Val{true})
x8 = Variant{Int,Char}(8, true)
@test_throws MethodError Variant{Int,Char}(9, Val{false})
@test_throws TypeError Variant{Int,Char}(10, false)
xa = Variant{Int,Char}('a')
@test_throws MethodError Variant{Int,Char}('b', Val{true})
@test_throws TypeError Variant{Int,Char}('c', true)
xd = Variant{Int,Char}('d', Val{false})
xe = Variant{Int,Char}('e', false)

@test lefttype(Variant{Int,Int}) === Int
@test righttype(Variant{Int,Int}) === Int
@test lefttype(Variant{Int,Char}) === Int
@test righttype(Variant{Int,Char}) === Char

# Basic properties

@test isleft(x1) && isleft(x2)
@test isright(x3) && isright(x4)
@test isleft(x6) && isleft(x7) && isleft(x8)
@test isright(xa) && isright(xd) && isright(xe)

@test getleft(x1) == 1
@test_throws VariantException getright(x1)
@test getleft(x1, 42) == 1
@test getright(x1, 42) == 42
@test isequal(getleft_null(x1), Nullable(1))
@test isequal(getright_null(x1), Nullable{Int}())

@test_throws VariantException getleft(x3)
@test getright(x3) == 3
@test getleft(x3, 42) == 42
@test getright(x3, 42) == 3
@test isequal(getleft_null(x3), Nullable{Int}())
@test isequal(getright_null(x3), Nullable(3))

@test getleft(x6) == 6
@test_throws VariantException getright(x6)
@test getleft(x6, 42) == 6
@test getright(x6, 42) == 42
@test isequal(getleft_null(x6), Nullable(6))
@test isequal(getright_null(x6), Nullable{Int}())

@test_throws VariantException getleft(xa)
@test getright(xa) == 'a'
@test getleft(xa, 'z') == 'z'
@test getright(xa, 'z') == 'a'
@test isequal(getleft_null(xa), Nullable{Char}())
@test isequal(getright_null(xa), Nullable('a'))

@test isequal(x1, Variant{Int,Int}(1, true))
@test !isequal(x1, Variant{Int,Int}(1, false))
@test !isequal(x3, Variant{Int,Int}(3, true))
@test isequal(x3, Variant{Int,Int}(3, false))

@test isequal(x6, Variant{Int,Char}(6))
@test !isequal(x6, Variant{Int,Char}('a'))
@test !isequal(xa, Variant{Int,Char}(6))
@test isequal(xa, Variant{Int,Char}('a'))

@test hash(x1) == hash(Variant{Int,Int}(1, true))
@test hash(x1) != hash(Variant{Int,Int}(1, false))
@test hash(x3) != hash(Variant{Int,Int}(3, true))
@test hash(x3) == hash(Variant{Int,Int}(3, false))

@test hash(x6) == hash(Variant{Int,Char}(6))
@test hash(x6) != hash(Variant{Int,Char}('a'))
@test hash(xa) != hash(Variant{Int,Char}(6))
@test hash(xa) == hash(Variant{Int,Char}('a'))

@test !isless(x1, Variant{Int,Int}(1, true))
@test isless(x1, Variant{Int,Int}(1, false))
@test !isless(x3, Variant{Int,Int}(3, true))
@test !isless(x3, Variant{Int,Int}(3, false))

@test !isless(x6, Variant{Int,Char}(6))
@test isless(x6, Variant{Int,Char}('a'))
@test !isless(xa, Variant{Int,Char}(6))
@test !isless(xa, Variant{Int,Char}('a'))

int = string(Int)   # either "Int32" or "Int64"
@test string(x1) == "Variant{$int,$int}(left:1)"
@test string(x3) == "Variant{$int,$int}(right:3)"
@test string(x6) == "Variant{$int,Char}(left:6)"
@test string(xa) == "Variant{$int,Char}(right:'a')"

# Non-trivial nesting

typealias V1 Variant{Int, Char}
typealias V2 Variant{Float64, AbstractString}
typealias V3 Variant{V1, V2}
v1a = V1(1)
v1b = V1('b')
v2a = V2(3.0)
v2b = V2("d")
v3aa = V3(v1a)
v3ab = V3(v1b)
v3ba = V3(v2a)
v3bb = V3(v2b)
@test getleft(getleft(v3aa)) == 1
@test getright(getleft(v3ab)) == 'b'
@test getleft(getright(v3ba)) == 3.0
@test getright(getright(v3bb)) == "d"

# Higher order functions

@test variant(Variant{Int,Int}(4,true)) == 4
@test variant(Variant{Int,Int}(4,false)) == 4

@test variant(+, -, Variant{Int,Int}(4,true)) == 4
@test variant(+, -, Variant{Int,Int}(4,false)) == -4

v1arr = V1[v1a, v1b, v1b]
v2arr = V2[v2b, v2a]
@test lefts(v1arr) == [1]
@test rights(v1arr) == ['b', 'b']
@test lefts(v2arr) == [3.0]
@test rights(v2arr) == ["d"]
@test partition_variants(v1arr) == (lefts(v1arr), rights(v1arr))
@test partition_variants(v2arr) == (lefts(v2arr), rights(v2arr))

typealias V4 Variant{Vector{V1}, Vector{V2}}
v4a = V4(v1arr)
v4b = V4(v2arr)
@test variant(length, sum, v4a) == 3
@test variant(max, length, v4b) == 2


text:
function sum_by_col(x)
    s = zero(eltype(x))
    for j in 1:size(x, 2)
        for i in 1:size(x, 1)
            s += x[i, j]
        end
    end
    s
end

function sum_by_row(x)
    s = zero(eltype(x))
    for i in 1:size(x, 1)
        for j in 1:size(x, 2)
            s += x[i, j]
        end
    end
    s
end

max_stars_repo_path: Chapter04/01. Traversing matrices efficiently/sums.jl


text:
using Gadfly

set_default_plot_size(7inch, 3inch)

p1 = plot(x=[0], y=[0], Scale.x_discrete, Geom.histogram2d)
p2 = plot(x=[0], y=[0], Scale.y_discrete, Geom.histogram2d)

hstack(p1, p2)


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 33012b19-0f0b-4943-aa99-e4a65a1de0f0
md"Hello"

# ╔═╡ e9d8b87c-f192-4ddc-8068-f25d132ad75f
	md"""
Oka
"""

# ╔═╡ f98090e8-23b4-44c1-99e3-da43571814bb

# ╔═╡ 43e5dfe4-90b2-4f7b-b467-3107b78520f4
x = 1 + 1

# ╔═╡ 036fdfbc-94c9-42e2-be5e-54e8c958ce65

# ╔═╡ Cell order:
# ╠═33012b19-0f0b-4943-aa99-e4a65a1de0f0
# ╠═e9d8b87c-f192-4ddc-8068-f25d132ad75f
# ╠═f98090e8-23b4-44c1-99e3-da43571814bb
# ╠═43e5dfe4-90b2-4f7b-b467-3107b78520f4
# ╠═036fdfbc-94c9-42e2-be5e-54e8c958ce65


text:
using Test

using BEAST
using CompScienceMeshes

V = [
    point(0,0,0),
    point(1,0,0),
    point(0,1,0)]
F = [index(1,2,3)]

m1 = Mesh(V,F)
m2 = CompScienceMeshes.rotate(m1, 0.5pi*point(1,0,0))
m3 = CompScienceMeshes.rotate(m1, 1.0pi*point(1,0,0))

m = weld(m1,m2,m3)

@test numcells(skeleton(m,0)) == 5
@test numcells(m) == 3

Nd = BEAST.nedelec(m)
@test numfunctions(Nd) == 7

function interior(mesh::Mesh, edges=skeleton(mesh,1))
    @assert dimension(mesh) == 2
    @assert vertices(mesh) === vertices(edges)

    C = connectivity(edges, mesh)
    @assert size(C) == (numcells(mesh), numcells(edges))

    nn = vec(sum(abs.(C), dims=1))
    T = CompScienceMeshes.celltype(edges)
    interior_edges = Vector{T}()
    for (i,edge) in pairs(cells(edges))
        nn[i] > 1 && push!(interior_edges, edge)
    end
    Mesh(vertices(mesh), interior_edges)
end

interior_edges = interior(m)
@test numcells(interior_edges) == 1
Nd = BEAST.nedelec(m, interior_edges)
@test numfunctions(Nd) == 1


text:
using Base.Test

using DataFrames
import FITSIO

import Celeste: AccuracyBenchmark
import Celeste: DeterministicVI
import Celeste: Model

@testset "flux <-> mags conversion" begin
    # reference values based on `nmgy2lups()` from
    # https://github.com/esheldon/sdsspy/blob/683d6e0f16a998240a129942f80ad3ce6e7d5dfe/sdsspy/util.py
    @test isapprox(AccuracyBenchmark.flux_to_mag(15.0, 1), 19.559677, atol=1e-5)
    @test isapprox(AccuracyBenchmark.flux_to_mag(15.0, 3),19.559702, atol=1e-5)
    @test isapprox(AccuracyBenchmark.mag_to_flux(19.559677, 1), 15.0, atol=1e-5)
    @test isapprox(AccuracyBenchmark.mag_to_flux(19.559702, 3), 15.0, atol=1e-5)
end

@testset "color calculations" begin
    @test isapprox(AccuracyBenchmark.color_from_fluxes(15.0, 20.0), log(20 / 15))
    @test isna(AccuracyBenchmark.color_from_fluxes(15.0, 0.0))
    fluxes = AccuracyBenchmark.fluxes_from_colors(10.0, [-1.0, 0.0, 1.0, 2.0])
    @test isapprox(fluxes[1], exp(1.0) * 10.0)
    @test isapprox(fluxes[2], 10.0)
    @test isapprox(fluxes[3], 10.0)
    @test isapprox(fluxes[4], exp(1.0) * 10.0)
    @test isapprox(fluxes[5], exp(3.0) * 10.0)
end

@testset "angle canonicalization" begin
    @test AccuracyBenchmark.canonical_angle(95.0) == 95.0
    @test AccuracyBenchmark.canonical_angle(195.0) == 15.0
    @test AccuracyBenchmark.canonical_angle(-20.0) == 160.0

    @test isapprox(AccuracyBenchmark.degrees_to_diff(20., -30.), 50.)
    @test isapprox(AccuracyBenchmark.degrees_to_diff(-10., 190.), 20.)
end

@testset "variational params -> data frame -> catalog entry conversion" begin
    variational_params = DeterministicVI.generic_init_source([1.0, 2.0])
    variational_params[Model.ids.gal_ab] = 0.5
    variational_params[Model.ids.gal_scale] = 10.0
    variational_params[Model.ids.gal_angle] = -pi / 4
    variational_params[Model.ids.flux_loc[2]] = log(20.0)
    variational_params[Model.ids.is_star[1]] = 0.01
    variational_params[Model.ids.is_star[2]] = 0.99

    data = AccuracyBenchmark.variational_parameters_to_data_frame_row(
        variational_params)
    # we'll just check a few particularly troublesome parameters :)
    @test isapprox(data[1, :half_light_radius_px], 10 * sqrt(0.5))
    @test isapprox(data[1, :angle_deg], 135.0)
    @test isapprox(data[1, :reference_band_flux_nmgy], 20.0)

    data[1, :is_star] = false
    catalog_entry = AccuracyBenchmark.make_catalog_entry(first(eachrow(data)))
    @test isapprox(catalog_entry.gal_scale, 10.0)
    @test isapprox(catalog_entry.gal_angle, 3 * pi / 4)
    @test isapprox(catalog_entry.gal_fluxes[3], 20.0)
end

@testset "PSF serialize/deserialize" begin
    psf = PsfComponent[
        PsfComponent(0.3, StaticArrays.@SVector([1., 2.]), StaticArrays.@SMatrix([3. 1.; 1. 4.])),
        PsfComponent(0.7, StaticArrays.@SVector([-2., -3.]), StaticArrays.@SMatrix([2. 0.; 0. 1.]))
    ]
    header = FITSIO.FITSHeader(String[], [], String[])

    AccuracyBenchmark.serialize_psf_to_header(psf, header)
    new_psf = AccuracyBenchmark.make_psf_from_header(header)

    @test new_psf == psf
end

const FITS_HEADER_STRING = "SIMPLE  =                    T /                                                BITPIX  =                  -32 / 32 bit floating point                          NAXIS   =                    2                                                  NAXIS1  =                 2048                                                  NAXIS2  =                 1489                                                  EXTEND  =                    T /Extensions may be present                       BZERO   =              0.00000 /Set by MRD_SCALE                                BSCALE  =              1.00000 /Set by MRD_SCALE                                TAI     =        4537928038.58 / 1st row - Number of seconds since Nov 17 1858  RA      =            359.76830 / 1st row - Right ascension of telescope boresighDEC     =            0.000000  / 1st row - Declination of telescope boresight (dSPA     =              90.000  / 1st row - Camera col position angle wrt north (IPA     =             102.229  / 1st row - Instrument rotator position angle (deIPARATE =              0.0000  / 1st row - Instrument rotator angular velocity (AZ      =            14.483496 / 1st row - Azimuth  (encoder) of tele (0=N?) (deALT     =            56.247810 / 1st row - Altitude (encoder) of tele        (deFOCUS   =           -436.80000 / 1st row - Focus piston (microns?)              DATE-OBS= '2002-09-05'         / 1st row - TAI date                             TAIHMS  = '07:33:58.57'        / 1st row - TAI time (HH:MM:SS.SS) (TAI-UT = apprCOMMENT  TAI,RA,DEC,SPA,IPA,IPARATE,AZ,ALT,FOCUS at reading of col 0, row 0     ORIGIN  = 'SDSS    '                                                            "

@testset "parse FITS header" begin
    header = AccuracyBenchmark.parse_fits_header_from_string(FITS_HEADER_STRING)
    @test header["SIMPLE"] == true
    @test header["NAXIS1"] == 2048
    @test isapprox(header["RA"], 359.7683)
    @test strip(header["ORIGIN"]) == "SDSS"
end

@testset "filtering rows for scoring errors" begin
    function make_data()
        (
            DataFrame(
                is_saturated=false,
                half_light_radius_px=10.0,
                de_vaucouleurs_mixture_weight=0.99,
                minor_major_axis_ratio=0.8,
            ),
            DataFrame(
                is_saturated=false,
                minor_major_axis_ratio=0.5,
                angle_deg=10.0,
                declination_deg=0.0,
            )
        )
    end

    truth, error = make_data()

    function check_row(column_name=:minor_major_axis_ratio)
        AccuracyBenchmark.is_good_row(first(eachrow(truth)), first(eachrow(error)), column_name)
    end

    @test check_row()

    truth, error = make_data()
    error[1, :minor_major_axis_ratio] = NaN
    @test !check_row()

    truth, error = make_data()
    error[1, :minor_major_axis_ratio] = NA
    @test !check_row()

    truth, error = make_data()
    truth[1, :half_light_radius_px] = NA
    @test check_row()

    truth, error = make_data()
    truth[1, :half_light_radius_px] = 25.0
    @test !check_row()

    truth, error = make_data()
    truth[1, :is_saturated] = true
    @test !check_row()

    truth, error = make_data()
    error[1, :is_saturated] = true
    @test !check_row()

    truth, error = make_data()
    truth[1, :de_vaucouleurs_mixture_weight] = 0.5
    @test !check_row()
    @test check_row(:declination_deg)

    truth, error = make_data()
    @test !check_row(:angle_deg)
    truth[1, :minor_major_axis_ratio] = 0.2
    @test check_row(:angle_deg)
end

@testset "match catalogs" begin
    ra = [0.0, 1.0, 2.0, 3.0]
    dec = [50.0, 51.0, 52.0, 53.0]
    truth = DataFrame(Any[ra, dec], [:right_ascension_deg, :declination_deg])

    off = 0.2 / 3600.0
    ra1 = [0.0 + off, 1.0 - off, 5.0, 3.0 + off, 4.0]
    dec1 = [50.0 - off, 51.0 + off, 60.0, 53.0 - off, 50.0]
    pred1 = DataFrame(Any[ra1, dec1], [:right_ascension_deg, :declination_deg])

    ra2 = [3.0 - off, 0.0 - off, 7.0, 5.0, 4.0]
    dec2 = [53.0 + off, 50.0 + off, 51.0, 60.0, 50.0]
    pred2 = DataFrame(Any[ra2, dec2], [:right_ascension_deg, :declination_deg])

    truth_matched, predictions_matched =
        AccuracyBenchmark.match_catalogs(truth, [pred1, pred2])

    @test nrow(truth_matched) == 2
    @test length(predictions_matched) == 2
    for pred in predictions_matched
        @test nrow(pred) == 2
    end
    @test truth_matched[:right_ascension_deg] == [0.0, 3.0]
    @test truth_matched[:declination_deg] == [50.0, 53.0]
    @test predictions_matched[1][:right_ascension_deg] == [0.0 + off, 3.0 + off]
    @test predictions_matched[1][:declination_deg] == [50.0 - off, 53.0 - off]
    @test predictions_matched[2][:right_ascension_deg] == [0.0 - off, 3.0 - off]
    @test predictions_matched[2][:declination_deg] == [50.0 + off, 53.0 + off]
end


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

mutable struct FrontendIPConfigurationPropertiesFormat <: SwaggerModel
    inboundNatRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: inboundNatRules
    inboundNatPools::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: inboundNatPools
    outboundRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: outboundRules
    loadBalancingRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: loadBalancingRules
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, IPAllocationMethod } # spec name: privateIPAllocationMethod
    privateIPAddressVersion::Any # spec type: Union{ Nothing, IPVersion } # spec name: privateIPAddressVersion
    subnet::Any # spec type: Union{ Nothing, Subnet2 } # spec name: subnet
    publicIPAddress::Any # spec type: Union{ Nothing, PublicIPAddress } # spec name: publicIPAddress
    publicIPPrefix::Any # spec type: Union{ Nothing, SubResource } # spec name: publicIPPrefix
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function FrontendIPConfigurationPropertiesFormat(;inboundNatRules=nothing, inboundNatPools=nothing, outboundRules=nothing, loadBalancingRules=nothing, privateIPAddress=nothing, privateIPAllocationMethod=nothing, privateIPAddressVersion=nothing, subnet=nothing, publicIPAddress=nothing, publicIPPrefix=nothing, provisioningState=nothing)
        o = new()
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("inboundNatRules"), inboundNatRules)
        setfield!(o, Symbol("inboundNatRules"), inboundNatRules)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("inboundNatPools"), inboundNatPools)
        setfield!(o, Symbol("inboundNatPools"), inboundNatPools)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("outboundRules"), outboundRules)
        setfield!(o, Symbol("outboundRules"), outboundRules)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("loadBalancingRules"), loadBalancingRules)
        setfield!(o, Symbol("loadBalancingRules"), loadBalancingRules)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        setfield!(o, Symbol("privateIPAddressVersion"), privateIPAddressVersion)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("publicIPPrefix"), publicIPPrefix)
        setfield!(o, Symbol("publicIPPrefix"), publicIPPrefix)
        validate_property(FrontendIPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type FrontendIPConfigurationPropertiesFormat

const _property_map_FrontendIPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("inboundNatRules")=>Symbol("inboundNatRules"), Symbol("inboundNatPools")=>Symbol("inboundNatPools"), Symbol("outboundRules")=>Symbol("outboundRules"), Symbol("loadBalancingRules")=>Symbol("loadBalancingRules"), Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("privateIPAddressVersion")=>Symbol("privateIPAddressVersion"), Symbol("subnet")=>Symbol("subnet"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("publicIPPrefix")=>Symbol("publicIPPrefix"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_FrontendIPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("inboundNatRules")=>"Vector{SubResource}", Symbol("inboundNatPools")=>"Vector{SubResource}", Symbol("outboundRules")=>"Vector{SubResource}", Symbol("loadBalancingRules")=>"Vector{SubResource}", Symbol("privateIPAddress")=>"String", Symbol("privateIPAllocationMethod")=>"IPAllocationMethod", Symbol("privateIPAddressVersion")=>"IPVersion", Symbol("subnet")=>"Subnet2", Symbol("publicIPAddress")=>"PublicIPAddress", Symbol("publicIPPrefix")=>"SubResource", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ FrontendIPConfigurationPropertiesFormat }) = collect(keys(_property_map_FrontendIPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ FrontendIPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FrontendIPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ FrontendIPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_FrontendIPConfigurationPropertiesFormat[property_name]

function check_required(o::FrontendIPConfigurationPropertiesFormat)
    true
end

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


text:
module SemianalyticalWrapper

using JavaCall

include("DsstWrapper/DsstWrapper.jl")

end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module Interface

using Random
using ExtractMacro
using ..Common

# this is silly but it's required for correct cross-linking in docstrings, apparently
using ...RRRMC

export Config, AbstractGraph, SimpleGraph, DiscrGraph, SingleGraph, DoubleGraph,
       spinflip!, energy, delta_energy, energy_type, neighbors, getN, allΔE, inner_graph,
       delta_energy_residual, update_cache!, update_cache_residual!,
       cenergy, distances

import Base: length, copy, copy!, ==

@static VERSION < v"1.1" && using Future

struct Config
    N::Int
    s::BitVector
    function Config(N::Integer; init::Bool = true)
        s = BitArray(undef, N)
        init && rand!(s)
        return new(N, s)
    end
end

@doc """
    Config(N::Integer)

The object storing the configuration for an Ising model. Although
the spin values are \$σ_i ∈ {-1,1}\$, internally they are stored in
a `BitArray`, in the type field `s`, so that to obtain the real
value one needs to perform the transformation \$σ_i = 2s_i - 1\$.
""" -> Config(N::Integer)

length(C::Config) = C.N

#spinflip!(C::Config, move::Int) = (C.s[move] $= 1)
spinflip!(C::Config, move::Int) = unsafe_bitflip!(C.s, move)

function copy!(dst::Config, src::Config)
    length(dst) == length(src) || throw(ArgumentError("Incompatible Config lengths, dst=$(length(dst)) src=$(length(src))"))
    M = @static VERSION < v"1.1" ? Future : Base
    M.copy!(dst.s, src.s)
    return dst
end

copy(C::Config) = copy!(Config(length(C), init=false), C)

(==)(C1::Config, C2::Config) = length(C1) == length(C2) && C1.s == C2.s

"""
    AbstractGraph{ET<:Real}

An abstract type representing an Ising spin model. The `ET` parameter
is the type returned by the [`energy`](@ref) and [`delta_energy`](@ref)
functions.

See also [`SimpleGraph`](@ref), [`DiscrGraph`](@ref), [`SingleGraph`](@ref)
and [`DoubleGraph`](@ref).
"""
abstract type AbstractGraph{ET<:Real} end

"""
    update_cache!(X::AbstractGraph, C::Config, move::Int)

A function which is called every time a spin is flipped. This may happen:

1. when a move is accepted, in [`standardMC`](@ref), [`rrrMC`](@ref), [`bklMC`](@ref) and
   [`wtmMC`](@ref).
2. when a move is attempted to evaluate the effect on the neighbors, in [`rrrMC`](@ref).

`move` is the spin index. By default, this function does nothing, but it may be overloaded
by particular graph types.

When `X` is a [`DoubleGraph`](@ref), there is a default implementation which first calls
`update_cache!` on [`inner_graph`](@ref)`(X)`, then
calls [`update_cache_residual!`](@ref) on `X`.

*Note:* this function is always invoked *after* the flip has been performed, unlike in [`delta_energy`](@ref)
and [`delta_energy_residual`](@ref).
"""
update_cache!(X::AbstractGraph, C::Config, move::Int) = nothing

function spinflip!(X::AbstractGraph, C::Config, move::Int)
    spinflip!(C, move)
    update_cache!(X, C, move)
end

"""
    energy(X::AbstractGraph, C::Config)

Returns the energy of graph `X` in the configuration `C`. This is always invoked at the
beginning of [`standardMC`](@ref), [`rrrMC`](@ref), [`bklMC`](@ref) and [`wtmMC`](@ref).
Subsequently, [`delta_energy`](@ref) is used instead.

All graphs must implement this function.

It *must* also be used to initialize/reset the cache for a given graph, if any (see [`update_cache!`](@ref)).
"""
energy(::AbstractGraph, ::Config) = error("not implemented")

"""
    energy_type(X::AbstractGraph{ET}) where ET

Returns the type `ET`, that is the return type of the [`energy`](@ref) and [`delta_energy`](@ref) functions.
"""
energy_type(::AbstractGraph{ET}) where ET = ET

"""
    delta_energy(X::AbstractGraph, C::Config, move::Int)

Returns the energy difference that would be associated to flipping the spin `move`.

A default fallback implementation based on `energy` is provided, to be used for debugging,
but having an efficient implementation for each graph is critical for performance.

*Note:* when `X` is a [`DiscrGraph`](@ref), the absolute value of the result must be contained in the
tuple returned by [`allΔE`](@ref) – no approximations are allowed, and missing values will cause crashes
(unless Julia is run with the `--check-bounds=yes` option, in which case they will cause errors).

*Note:* this function is always invoked *before* performing the flip, unlike in [`update_cache!`](@ref)
and [`update_cache_residual!`](@ref).
"""
function delta_energy(X::AbstractGraph, C::Config, move::Int)
    @extract C : s
    oldn = energy(X, C)
    s[move] ⊻= 1
    newn = energy(X, C)
    s[move] ⊻= 1
    Δn0 = newn - oldn
    return Δn0
end

"""
    getN(X::AbstractGraph)

Returns the number of spins for a graph. The default implementation just returns `X.N`.
"""
getN(X::AbstractGraph) = X.N

"""
    neighbors(X::AbstractGraph, i::Int)

Returns an iterable with all the neighbors of spin `i`. This is required
by [`rrrMC`](@ref), [`bklMC`](@ref) and [`wtmMC`](@ref) since those methods need to
evaluate the effect of flipping a spin on its neighbors' local fields. It is
not required by [`standardMC`](@ref).

For performance reasons, it is best if the returned value is stack-allocated
rather than heap-allocated, e.g. it is better to return an `NTuple` than a `Vector`.
"""
neighbors(::AbstractGraph, i::Int) = error("not implemented")

"""
    SimpleGraph{ET} <: AbstractGraph{ET}

An abstract type representing a generic graph.

The `ET` parameter is the type returned by [`energy`](@ref) and [`delta_energy`](@ref).
"""
abstract type SimpleGraph{ET} <: AbstractGraph{ET} end

"""
    DiscrGraph{ET} <: AbstractGraph{ET}

An abstract type representing a graph in which the [`delta_energy`](@ref) values
produced when flipping a spin belong to a finite discrete set, and thus can be
sampled more efficiently with [`rrrMC`](@ref) or [`bklMC`](@ref) if the set is
small.

The `ET` parameter is the type returned by [`energy`](@ref) and [`delta_energy`](@ref).

See also [`allΔE`](@ref).
"""
abstract type DiscrGraph{ET} <: AbstractGraph{ET} end

"""
    allΔE{P<:DiscrGraph}(::Type{P})
    allΔE(X::DiscrGraph)

Returns a tuple of all possible *non-negative* values that can be returned
by [`delta_energy`](@ref). This must be implemented by all `DiscrGraph`
objects in order to use [`rrrMC`](@ref) or [`bklMC`](@ref).

The second form above, called directly on the `DiscrGraph` object, takes
precedence; however, for performance reasons, it is best if the result can
be computed from the type of the graph alone (possibly using a generated
function), i.e. it is best to implement the first definition above — the
second one then falls back automatically to the first.
"""
allΔE(::Type{P}) where {P<:DiscrGraph} = error("not implemented")
allΔE(X::DiscrGraph) = allΔE(typeof(X))

"""
    SingleGraph{ET}

A type alias representing either a [`SimpleGraph`](@ref) or a
[`DiscrGraph{ET}`](@ref). See also [`DoubleGraph`](@ref).
"""
const SingleGraph{ET} = Union{SimpleGraph{ET},DiscrGraph{ET}}

"""
    DoubleGraph{GT<:SingleGraph,ET} <: AbstractGraph{ET}

An abstract type representing a graph in which the energy is the sum of two
contributions, one of which is encoded in a graph of type `GT` (see
[`SingleGraph`](@ref)). This allows [`rrrMC`](@ref) to sample values
more efficiently.

The `ET` parameter is the type returned by the [`energy`](@ref) and
[`delta_energy`](@ref) functions. Note that it can be different from
the energy type of the internal `GT` object (e.g., one can have
a `DoubleGraph{DiscrGraph{Int},Float64}` object).

*Note:* When you declare a type as subtype of this, `GT` should *not* be the
concrete type of the inner graph, but either `SimpleGraph{T}` or
`DiscrGraph{T}` for some `T`.

See also [`inner_graph`](@ref), [`delta_energy_residual`](@ref) and
[`update_cache_residual!`](@ref).
"""
abstract type DoubleGraph{GT,ET} <: AbstractGraph{ET} end

"""
    inner_graph(X::DoubleGraph)

Returns the internal graph used by the given [`DoubleGraph`](@ref).
The default implementation simply returns `X.X0`.
"""
inner_graph(X::DoubleGraph) = X.X0
inner_graph(X::SingleGraph) = X

"""
    delta_energy_residual(X::DoubleGraph, C::Config, move::Int)

Returns the residual part of the energy difference produced if the spin `move` would
be flipped, excluding the contribution from the internal [`SimpleGraph`](@ref)
(see [`inner_graph`](@ref)).

See also [`delta_energy`](@ref). There is a default fallback implementation, but
it should be overloaded for efficiency.
"""
delta_energy_residual(X::DoubleGraph, C::Config, move::Int) =
    delta_energy(X, C, move) - delta_energy(inner_graph(X), C, move)

"""
    update_cache_residual!(X::DoubleGraph, C::Config, move::Int)

Called internally by the default [`update_cache!`](@ref) when the argument is a [`DoubleGraph`](@ref).
Can be useful to overload this if the residual part of the graph has an indipendent cache.

By default, it does nothing.
"""
update_cache_residual!(X::DoubleGraph, C::Config, move::Int) = nothing

function update_cache!(X::DoubleGraph, C::Config, move::Int)
    update_cache!(inner_graph(X), C, move)
    update_cache_residual!(X, C, move)
end

allΔE(X::DoubleGraph{ET,GT}) where {ET,GT<:DiscrGraph} = allΔE(typeof(inner_graph(X)))

"""
    cenergy(X::DoubleGraph)

Returns the individual energy (as defined by the original model)
of the reference configuration in a [`GraphLocalEntropy`](@ref RRRMC.GraphLocalEntropy) graph or a
[`GraphTopologicalLocalEntropy`](@ref RRRMC.GraphTopologicalLocalEntropy) graph.
"""
cenergy(::AbstractGraph) = error("not implemented")

"""
    distances(X::DoubleGraph)

Returns the matrix of the Hamming distances between replicas, in graphs with
replicas: [`GraphQuant`](@ref RRRMC.GraphQuant), [`GraphLocalEntropy`](@ref RRRMC.GraphLocalEntropy),
[`GraphTopologicalLocalEntropy`](@ref RRRMC.GraphTopologicalLocalEntropy) or [`GraphRobustEnsemble`](@ref RRRMC.GraphRobustEnsemble).
"""
distances(::AbstractGraph) = error("not implemented")

end # module


text:
using SeisDvv
using Test

@testset "SeisDvv.jl" begin
    # Write your own tests here.
end


text:
__precompile__()
module Neo4Julia

# Import/Export
importall PackageStream

# Includes
include("connections.jl")
include("session_api.jl")

end # module


text:
module ZippedArrays

export ZippedArray

# Alias for a tuple of arrays.
const ArrayTuple{L,N} = NTuple{L,AbstractArray{<:Any,N}}

"""
    ZippedArray(A,B,C,...) -> Z

yields a zipped array instance `Z` based on arrays `A`, `B`, `C`, etc.  such
that the syntax `Z[i]` yields a tuple of values `(A[i],B[i],C[i],...)` while
the syntax `Z[i] = (a,b,c,...)` is equivalent to `(A[i],B[i],C[i],...) =
(a,b,c,...)`.

Any number of arrays can be zipped together, they must however have the same
indices (as given by calling the `axes` method).

""" ZippedArray

# Parameters:
#     T = array element type;
#     N = number of dimensions;
#     L = number of zipped arrays
#     I = indexing style (true for linear, false otherwise)
#     S = type of L-tuple of zipped arrays
struct ZippedArray{T,N,L,I,S<:ArrayTuple{L,N}} <: AbstractArray{T,N}
    args::S
end

# Alias for zipped arrays with fast linear indexing.
const FastZippedArray{T,N,L,S} = ZippedArray{T,N,L,true,S}

ZippedArray(args::AbstractArray{<:Any,N}...) where {N} =
    ZippedArray{Tuple{map(eltype,args)...}, N, length(args),
                get_index_style(args...) === IndexLinear(),
                typeof(args)}(args)

ZippedArray() = error("at least one array argument must be provided")

Base.length(A::ZippedArray) = length(A.args[1])

Base.size(A::ZippedArray) = size(A.args[1])
Base.size(A::ZippedArray, i) = size(A.args[1], i)

Base.axes1(A::ZippedArray) = Base.axes1(A.args[1])
Base.axes(A::ZippedArray) = axes(A.args[1])
Base.axes(A::ZippedArray, i::Integer) = axes(A.args[1], i)

Base.IndexStyle(::Type{<:FastZippedArray}) = IndexLinear()
Base.IndexStyle(::Type{<:ZippedArray}) = IndexCartesian()

@generated function Base.getindex(A::FastZippedArray{T,N,L},
                                  i::Int) where {T,N,L}
    lhs = Expr(:tuple, ntuple(j -> :(A.args[$j][i]), Val(L))...)
    quote
        $(Expr(:meta, :inline))
        @boundscheck checkbounds(A, i)
        @inbounds r = $lhs
        r
    end
end

@generated function Base.setindex!(A::FastZippedArray{T,N,L},
                                   val, i::Int) where {T,N,L}
    lhs = Expr(:tuple, ntuple(j -> :(A.args[$j][i]), Val(L))...)
    quote
        $(Expr(:meta, :inline))
        @boundscheck checkbounds(A, i)
        @inbounds $lhs = val
        A
    end
end

@inline function Base.checkbounds(::Type{Bool},
                                  A::FastZippedArray{T,N,L},
                                  i::Int) where {T,N,L}
    checkbounds(Bool, A.args[1], i)
end

@generated function Base.getindex(A::ZippedArray{T,N,L},
                                  i::Vararg{Int,N}) where {T,N,L}
    lhs = Expr(:tuple, ntuple(j -> :(A.args[$j][i...]), Val(L))...)
    quote
        $(Expr(:meta, :inline))
        @boundscheck checkbounds(A, i...)
        @inbounds r = $lhs
        r
    end
end

@generated function Base.setindex!(A::ZippedArray{T,N,L},
                                   val, i::Vararg{Int,N}) where {T,N,L}
    lhs = Expr(:tuple, ntuple(j -> :(A.args[$j][i...]), Val(L))...)
    quote
        $(Expr(:meta, :inline))
        @boundscheck checkbounds(A, i...)
        @inbounds $lhs = val
        A
    end
end

@inline function Base.checkbounds(::Type{Bool},
                                  A::ZippedArray{T,N,L},
                                  i::Vararg{Int,N}) where {T,N,L}
    checkbounds(Bool, A.args[1], i...)
end

"""
    all_match(val, f, args...) -> bool

yields whether `f(arg) == val` for all `arg ∈ args`.

"""
all_match(val, f::Function) = true
all_match(val, f::Function, A) = f(A) == val
@inline all_match(val, f::Function, A, B...) =
    all_match(val, f, A) && all_match(val, f::Function, B...)

"""
    get_index_style(A...) -> IndexLinear() or IndexCartesian()

yields the most efficient indexing style for accessing arrays `A...` with the
same index.  An exception is thrown if arguments do not have the same axes.

"""
get_index_style(A::AbstractArray) = IndexStyle(A)

@inline get_index_style(A::AbstractArray, B::AbstractArray) =
    get_index_style(IndexStyle(A, B), A, B)
@inline get_index_style(A::AbstractArray, B::AbstractArray...) =
    get_index_style(IndexStyle(A, B...), A, B...)

@inline function get_index_style(I::Union{IndexLinear,IndexCartesian},
                                 A::AbstractArray,
                                 B::AbstractArray...)
    all_match(axes(A), axes, B...) || throw_indices_mismatch(A, B...)
    return I
end

@noinline throw_indices_mismatch(A::AbstractArray...) =
    throw_indices_mismatch(Base.has_offset_axes(A...), A...)

@noinline function throw_indices_mismatch(show_axes::Bool, A::AbstractArray...)
    io = IOBuffer()
    write(io, "all arguments of `ZippedArray` must have the same ",
          (show_axes ? "axes" : "shape"), ", got ")
    m = length(A)
    for i in 1:m
        write(io, (i == 1 ? "(" : i == m ? " and (" : ", )"))
        if show_axes
            inds = axes(A[i])
            n = length(inds)
            for j in 1:n
                j > 1 && write(io, ",")
                print(io, first(inds[j]), ":", last(inds[j]))
            end
        else
            dims = size(A[i])
            n = length(dims)
            for j in 1:n
                j > 1 && write(io, ",")
                print(io, dims[j])
            end
        end
        write(io, (n == 1 ? ",)" : ")"))
    end
    throw(DimensionMismatch(String(take!(io))))
end

end # module


text:
#
# file: tensorgrid.jl
#
# Creation of 1D,2D,3D tensor grids
#

export tensor_grid

# This meshgrid functions are extracted from Julia (License is MIT) and modified
# by Dirk Peschka: Unfortunately this seemed not to be a callable part of Julia
# 0.4.0 base. Extracted only meshgrid functionality from the ndgrid.jl file available at
#
# https://github.com/JuliaLang/julia/blob/master/examples/ndgrid.jl
#
# and modified swapping x <-> y in the meshgrid output, that
# is probably a (strange) artifact from MATLAB.

"""
**Summary:**     `meshgrid{T}(vx, vy) :: Array{T,2}`

**Description:** `meshgrid{T}(vx,vy)` returns as type `T` two arrays xx,yy, so
that x[i,j]=vx[i], y[i,j]=vy[j] for all i,j.
"""
function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
    m, n = length(vy), length(vx)
    vx = reshape(vx, n, 1)
    vy = reshape(vy, 1, m)
    (repmat(vx, 1, m), repmat(vy, n, 1))
end

"""
**Summary:**     `meshgrid{T}(vx, vy, vz) :: Array{T,3}`

**Description:** `meshgrid{T}(vx,vy)` returns as type `T` two arrays x,y,z so
that x[i,j,k]=vx[i], y[i,j,k]=vy[j], z[i,j,k]=vz[k] for all i,j,k.
"""
function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T},
                     vz::AbstractVector{T})
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, n, 1, 1)
    vy = reshape(vy, 1, m, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end

"""
**Summary:** `isincreasing(x) :: Bool`

**Description:** The `function isincreasing(x)` returns true, if the elements
of the iterable object `x` are in strict increasing order, otherwise it returns
false.

Examples:

    * isincreasing([1 2 3]) => true
    * isincreasing([1 2 2]) => false
"""
function isincreasing(x)
  cnt = start(x)
  done(x,cnt) && return true
  prev, cnt = next(x, cnt)
  while !done(x, cnt)
    this, cnt = next(x, cnt)
    !isless(prev,this) && return false
    prev = this
  end
  return true
end

"""
**Summary:** `function tensor_grid(x::AbstractArray,pattr=1,eattr=1)`

**Description:** Creates a 1D tensor grid out of vertex positions in `x`.

Keyword arguments:

  * `pattr::Int` number of vertex attributes
  * `eattr::Int` number of element attributes
Examples:

    * g = tensor_grid(linspace(0,1,5))
    * g = tensor_grid([0. 0.1 0.3 0.6 1.0])
"""
function tensor_grid(xx::AbstractArray,pattr=1,eattr=1)

  # check if x strictly increasing
  x = collect(float(xx))
  if !isincreasing(x)
    error("Requires strictly increasing input!")
  end

  # create indices for topology
  n = length(x)
  m = pattr  # point attributes as a parameter (with default=1)?
  k = eattr  # element attributes (see above)
  x    = reshape(x,n,1)
  v    = [collect(1:n) collect(1:n) zeros(Int64,n,m)]
  # mark boundaries by default? or rather not
  v[1,3] = 1
  v[n,3] = 1
  elem = [collect(1:n-1) collect(1:(n-1)) collect(2:n) zeros(Int64,n-1,k)]

  # create geometry and topology using standard constructors
  vertices = VertexTopology(v,1,m)
  geometry = Geometry{1}(x,vertices)
  topology = IntervalTopology(elem,2,k)

  # create Grid type
  Grid{1,IntervalTopology}(geometry,topology)

end

"""
 Creates 2D tensor grid
"""
function tensor_grid(x::AbstractArray, y::AbstractArray)

  # check if x,y strictly increasing
  if !isincreasing(x)|!isincreasing(y)
    error("Requires strictly increasing input!")
  end

end

"""
Creates 3D tensor grid
"""
function tensor_grid(x::AbstractArray, y::AbstractArray, z::AbstractArray)
  # check whether arrays are monotone
  # check for singular subintervals
  # create arrays for coordinates and cells and do nested loops to fill them
  # create array for boundary facets and fill it
  if !isincreasing(x)|!isincreasing(y)|!isincreasing(z)
    error("Requires strictly increasing input!")
  end

end

#=
Assigns region markers to certain cells of the grid
=#
#function mark_cells(g::AbstractGrid)
#end

#=
Assigns region markers to certain facets of the grid
=#
#function mark_facets(g::AbstractGrid)
#end

#=
Remove parts of the grid
=#
#function extract_subgrid(g::AbstractGrid)
#end

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
# Issue #14
stack = @make_stack map(x->2x, 1:10)
state = dummy_state(stack)
execute_command(state, state.stack[1], Val{:finish}(), "finish")
@test isempty(state.stack)
@test state.overall_result == 2 .* [1:10...]

# Issue #12
function complicated_keyword_stuff(args...; kw...)
    args[1] == args[1]
    (args..., kw...)
end
stack = @make_stack complicated_keyword_stuff(1)
state = dummy_state(stack)
execute_command(state, state.stack[1], Val{:n}(), "n")
execute_command(state, state.stack[1], Val{:finish}(), "finish")
@test isempty(state.stack)

@test runframe(ASTInterpreter2.enter_call(complicated_keyword_stuff, 1, 2)) ==
      runframe(@make_stack(complicated_keyword_stuff(1, 2))[1])
@test runframe(ASTInterpreter2.enter_call(complicated_keyword_stuff, 1, 2; x=7, y=33)) ==
      runframe(@make_stack(complicated_keyword_stuff(1, 2; x=7, y=33))[1])

# Issue #22
f22() = string(:(a+b))
@test step_through(enter_call_expr(:($f22()))) == "a + b"
f22() = string(QuoteNode(:a))
@test step_through(enter_call_expr(:($f22()))) == ":a"


text:
#' @description
#'
#' An AbstractDataArray is an Array whose entries can take on
#' values of type `T` or the value `NA`.
abstract AbstractDataArray{T, N} <: AbstractArray{T, N}

#' @description
#'
#' An AbstractDataVector is an AbstractDataArray of order 1.
typealias AbstractDataVector{T} AbstractDataArray{T, 1}

#' @description
#'
#' An AbstractDataMatrix is an AbstractDataArray of order 2.
typealias AbstractDataMatrix{T} AbstractDataArray{T, 2}

#' @description
#' Determine the type of the elements of an AbstractDataArray.
#'
#' @param ada::AbstractDataArray{T} The AbstractDataArray whose
#'        element type is desired.
#'
#' @returns T::DataType The type of the non-NA elements of `ada`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' T = eltype(dv)
Base.eltype{T, N}(d::AbstractDataArray{T, N}) = T

# Generic iteration over AbstractDataArray's

# TODO: Document
Base.start(x::AbstractDataArray) = 1

# TODO: Document
Base.next(x::AbstractDataArray, state::Integer) = (x[state], state + 1)

# TODO: Document
Base.done(x::AbstractDataArray, state::Integer) = state > length(x)

#' @description
#'
#' Determine if any of the entries of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are any of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' anyna(a)
anyna(a::AbstractArray) = false # -> Bool

#' @description
#'
#' Determine if any of the entries of an DataArray are `NA`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are any of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' anyna(da)
anyna(d::AbstractDataArray) = any(isna, d) # -> Bool

#' @description
#'
#' Determine if all of the entries of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are all of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' allna(a)
allna(a::AbstractArray) = false # -> Bool

#' @description
#'
#' Determine if all of the entries of an DataArray are `NA`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are all of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' allna(da)
allna(d::AbstractDataArray) = all(isna, d) # -> Bool

#' @description
#'
#' Determine if the values of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose missingness will
#'        be assessed.
#'
#' @returns na::BitArray{N} Elementwise Boolean whether entry is missing.
#'
#' @examples
#'
#' a = [1, 2, 3]
#' isna(a)
isna(a::AbstractArray) = falses(size(a)) # -> BitArray

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' AbstractArray is `NA`.
#'
#' @param a::AbstractArray The AbstractArray whose missingness will
#'        be assessed.
#' @param i::Integer The index of the element to be checked for `NA`.
#'
#' @returns na::Bool Is the element `NA` or not?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' isna(a, 1)
isna(a::AbstractArray, i::Real) = false # -> Bool

# TODO: Remove this method? It seems like a trap to me.
#' @description
#'
#' Insert a set of values, `vals`, into the positions described
#' by `inds`. The inserted values, `vals`, are echoed back as the return
#' value.
#'
#' NB: The indices in `inds` must be exhaustive.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param vals::AbstractVector A set of values to be inserted into the
#'        specified entries of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to
#'        be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[false, true, true]] = [5.0, 5.0]
function Base.setindex!(da::AbstractDataArray,
                        vals::AbstractVector,
                        inds::AbstractVector{Bool})
	if size(inds) != size(da)
		throw(ArgumentError("Boolean indices must be exhaustive"))
	end
    return setindex!(da, vals, find(inds))
end

#' @description
#'
#' Insert a set of values, `vals`, into the positions described
#' by `inds`. The inserted values, `vals`, are echoed back as the return
#' value.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param vals::AbstractVector A set of values to be inserted into the
#'        specified entries of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to
#'        be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[1, 3]] = [5.0, 5.]
function Base.setindex!(da::AbstractDataArray,
                        vals::AbstractVector,
                        inds::AbstractVector)
	if length(vals) != length(inds)
		throw(ArgumentError("Values and indices must have the same length"))
	end
    for (val, ind) in zip(vals, inds)
        da[ind] = val
    end
    return da
end

#' @description
#'
#' Insert a value, `val`, of type `Any` into the positions described
#' by `inds`. The inserted value, `val`, is echoed back as the return
#' value.
#'
#' NB: The indices in `inds` must be exhaustive.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::Any A value to be inserted into the specified entries of `da`.
#' @param inds::AbstractVector{Bool} A Boolean vector of indices of `da` to
#'        be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[false, true]] = 5.0
function Base.setindex!{T}(da::AbstractDataArray{T},
                           val::Any,
                           inds::AbstractVector{Bool})
    return setindex!(da, convert(T, val), find(inds))
end

#' @description
#'
#' Insert a value, `val`, of type `Any` into the positions described
#' by `inds`. The inserted value, `val`, is echoed back as the return
#' value. Note that `val` be coerced to the element type of `da`, which
#' may fail.
#'
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::Any A value to be inserted into the specified entries of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[1:3] = 5.0
function Base.setindex!{T}(da::AbstractDataArray{T},
                           val::Any,
                           inds::AbstractVector)
    val = convert(T, val)
    for ind in inds
        da[ind] = val
    end
    return da
end

# TODO: Include BitArray indexing here.
#' @description
#'
#' Insert a value, `val`, of type `Any` into the positions described
#' by `inds`. The inserted value, `val`, is echoed back as the return
#' value.
#'
#' NB: The vector of indices, `inds`, must exhausitively state for every
#' element of `da` whether it is being modified or not.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::Any A value to be inserted into the specified entries of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[true, false, false]] = 5.0
function Base.setindex!(da::AbstractDataArray,
                        val::Any,
                        inds::AbstractVector{Bool})
    # TODO: Make inds a BooleanIndex
    if size(inds) != size(da)
        throw(ArgumentError("Boolean indices must be exhausitive"))
    end
    return setindex!(da, val, find(inds))
end

#' @description
#'
#' Insert a value, `val`, of type `Any` into the positions described
#' by `inds`. The inserted value, `val`, is echoed back as the return
#' value.
#'
#' @param da::DataArray{T, N} A DataArray whose entries will be modified.
#' @param val::Any A value to be inserted into the specified entries of `da`.
#' @param inds::AbstractVector A vector of indices of `da` to be modified.
#'
#' @returns da::DataArray{T, N} The modified DataArray.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' da[[1, 2]] = 5.0
function Base.setindex!{T}(da::AbstractDataArray{T},
                           val::Any,
                           inds::AbstractVector)
    val = convert(T, val)
    for ind in inds
        da[ind] = val
    end
    return da
end


text:
using SafeTestsets

@safetestset "Complex transformation math tests" begin include("transform_tests.jl") end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
using OrdinaryDiffEq: ODEProblem, Rodas5, solve!, DiscreteCallback, CallbackSet, _initialize_dae!, u_modified!
using DiffEqCallbacks: GeneralDomain
using DiffEqBase: solve
using Setfield

import PowerDynamics: simulate

"""
```Julia
modified_copy(component, key, value)
```

provides a copy of an immutable type and modifies a given field

# Arguments
- `component`: type instance
- `key`: field name that should be modified
- `value`: new value for field `key`

Note that `value` needs to have the appropriate field type.
"""
function modified_copy(component, key::Symbol, value)
    kwargs = Dict{Symbol,Any}()
    # collect values
    for fn in component |> typeof |> fieldnames
        kwargs[fn] = getfield(component, fn)
    end
    # modify field
    @assert value isa fieldtype(typeof(component), key)
    kwargs[key] = value
    return eval(typeof(component))(; kwargs...)
end

using NLsolve: nlsolve, converged
using LinearAlgebra: pinv

struct RootRhs_ic
    rhs
    mpm
end
function (rr::RootRhs_ic)(x)
    dx = similar(x)
    rr.rhs(dx, x, nothing, 0.0)
    rr.mpm * dx .- dx
end

function RootRhs_ic(of::ODEFunction)
    mm = of.mass_matrix
    @assert mm != nothing
    mpm = pinv(mm) * mm
    RootRhs_ic(of.f, mpm)
end

function find_valid_initial_condition(pg, ic_guess)
    rr = RootRhs_ic(rhs(pg))
    nl_res = nlsolve(
        rr,
        ic_guess,
        autodiff = :forward,
        show_trace = false,
        method=:trust_region,
    )
    if converged(nl_res) == true
        return nl_res.zero #State(pg, nl_res.zero)
    else
        println("Failed to find initial conditions on the constraint manifold!")
        println("Try running nlsolve with other options.")
    end
end

"""
```Julia
    NodeShortCircuitAlt(;node_number,Y,tspan_fault)
```
# Keyword Arguments
- `node_number`: number  of the node
- `Y`: admittance of the short circuit
- `tspan_fault`: short circuit timespan
"""
Base.@kwdef struct NodeShortCircuitAlt
    node_number
    Y
    tspan_fault
    shunt_symbol = :Y_n
end

"Error to be thrown if something goes wrong during short circuit"
struct NodeShortCircuitAltError <: PowerDynamicsError
    msg::String
end

function (nsc::NodeShortCircuitAlt)(powergrid)
    nodes = copy(powergrid.nodes)
    sc_node = powergrid.nodes[nsc.node_number]

    # add shunt to the node component at the fault location
    if !(hasproperty(sc_node, nsc.shunt_symbol))
        throw(NodeShortCircuitAltError("Node number: $(nsc.node_number) must have a shunt field called $(nsc.shunt_symbol)."))
    end
    lens = Setfield.compose(Setfield.PropertyLens{nsc.shunt_symbol}())
    faulted_component = Setfield.set(sc_node, lens, nsc.Y)

    nodes[nsc.node_number] = faulted_component
    ControlledPowerGrid(
        powergrid.graph,
        nodes,
        powergrid.lines,
        powergrid.controller,
    )
end

"""
```Julia
simulate(nsc::NodeShortCircuitAlt, powergrid, x1, timespan)
```
Simulates a [`NodeShortCircuitAlt`](@ref)
"""
function simulate(nsc::NodeShortCircuitAlt, powergrid, x1, timespan)
    @assert first(timespan) <= nsc.tspan_fault[1] "fault cannot begin in the past"
    @assert nsc.tspan_fault[2] <= last(timespan) "fault cannot end in the future"

    nsc_powergrid = nsc(powergrid)

    problem = ODEProblem{true}(rhs(powergrid), x1, timespan)

    init_alg =  BrownFullBasicInit()

    function errorState(integrator)
        sol1 = integrator.sol
        println("searching fault state")
        #x2 = find_valid_initial_condition(nsc_powergrid, sol1[end]) # Jump the state to be valid for the new system.
        integrator.f = rhs(nsc_powergrid)
        DiffEqBase.initialize_dae!(integrator, init_alg)
        println("found fault state")
        println(State(powergrid, integrator.u)[:, :v])
        #integrator.u = x2
        #u_modified!(integrator,true)
    end

    function regularState(integrator)
        sol2 = integrator.sol
        println("searching post-fault state")
        #x3 = find_valid_initial_condition(powergrid, sol2[end]) # Jump the state to be valid for the new system.
        integrator.f = rhs(powergrid)
        DiffEqBase.initialize_dae!(integrator, init_alg)
        println("found post-fault state")
        println(State(powergrid, integrator.u)[:, :v])
        #integrator.u = x3
        #u_modified!(integrator,true)
    end

    t1 = nsc.tspan_fault[1]
    t2 = nsc.tspan_fault[2]
    dt = 1e-7

    # discrete callback requires tstop
    cb1 = DiscreteCallback(
        ((u, t, integrator) -> t in nsc.tspan_fault[1]),
        errorState,
    )
    cb2 = DiscreteCallback(
        ((u, t, integrator) -> t in nsc.tspan_fault[2]),
        regularState,
    )

    # cb1 = ContinuousCallback(
    #     ((u, t, integrator) -> t - nsc.tspan_fault[1]),
    #     errorState,
    # )
    # cb2 = ContinuousCallback(
    #     ((u, t, integrator) -> t - nsc.tspan_fault[2]),
    #     regularState,
    # )

    # Vlim = 0.1 #pu
    # v_idx = Vector{Int}()
    # offset = 0.
    # for node in cpg.nodes
    #     push!(v_idx, offset + 1)
    #     push!(v_idx, offset + 2)
    #     offset += dimension(node)
    # end
    # function is_valid_voltage(u, resid)
    #     #v = State(powergrid, u)[:, :v]
    #     #resid .= sum(map(x-> x<Vlim ? Vlim - x : 0., v))
    #     for i in 1:length(u)
    #         if i ∈ v_idx
    #             resid[i] =  abs(u[i]) < Vlim ? abs(u[i]) - Vlim : 0.
    #         else
    #             resid[i] = 0.
    #         end
    #     end
    # end
    #
    # cb3 = GeneralDomain(is_valid_voltage; autonomous=true)

    sol = solve(
        problem,
        Rodas5(autodiff = true),
        #force_dtmin = true,
        initializealg=BrownFullBasicInit(),
        callback = CallbackSet(cb1, cb2),
        tstops = [t1, t2],
        #tstops = [t1-dt, t1+dt, t2-dt, t2+dt],
    )
    return PowerGridSolution(sol, powergrid)
end

simulate(nsc::NodeShortCircuitAlt, op::State, timespan) =
    simulate(nsc, op.grid, op.vec, timespan)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
@testset "bitboard.jl" begin
    @printf("> bitboard\n")

    b = Bb.setBoard()
    @test b.black.P == Bb.MASK_RANKS[7]
    @test b.white.P == Bb.MASK_RANKS[2]
    @test b.white.K == p2u["e1"]
    @test b.black.K == p2u["e8"]
    @test b.white.Q == p2u["d1"]
    @test b.black.Q == p2u["d8"]
    @test b.white.R == p2u["a1"] | p2u["h1"]
    @test b.black.R == p2u["a8"] | p2u["h8"]
    @test b.white.B == p2u["c1"] | p2u["f1"]
    @test b.black.B == p2u["c8"] | p2u["f8"]
    @test b.white.N == p2u["b1"] | p2u["g1"]
    @test b.black.N == p2u["b8"] | p2u["g8"]

    [@test Bb.INT2UINT[i] & b.taken == Bb.EMPTY for i = 17:48]
    [@test Bb.INT2UINT[i] & ~b.taken == Bb.EMPTY for i = 1:16]
    [@test Bb.INT2UINT[i] & ~b.taken == Bb.EMPTY for i = 49:64]

    @test b.active
    @test b.castling == UInt8(15)
    @test b.enpassant == Bb.EMPTY
    @test b.halfmove == 0
    @test b.fullmove == 1
end


text:
#=
Given a string, find the longest palindromic contiguous substring. If there are more than one with the maximum length, return any one.

For example, the longest palindromic substring of "aabcdcb" is "bcdcb". The longest palindromic substring of "bananas" is "anana".
=#

using Test
include("Solutions/problem46_findlongestpalindrome.jl")

@test longestpalindromicsubstring("aabcdcb") == "bcdcb"
@test longestpalindromicsubstring("bananas") == "anana"

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
@recipe(RecipePlot) do scene
    th = merge(
        default_theme(scene),
        Attributes(palette = Palette(rwong))
    )
    th.color = automatic
    return th
end

function plot!(p::T) where T <: RecipePlot

    # What happens here is that I want to lift on every available observable,
    # so they need to be splatted.  This also means that nested attributes
    # will not be lifted on, but that's an acceptable tradeoff.
    #
    # After lifting on everything,

    # Node(1) is a dummy observable for dispatch.
    lift(Node(1), p.attributes, p.converted, p.converted..., values(p.attributes)...) do _, attrs, args, __lifted...

        !isempty(p.plots) && empty!(p.plots)

        RecipesPipeline.recipe_pipeline!(
            p,
            Dict{Symbol, Any}(keys(attrs) .=> to_value.(values(attrs))),
            to_value.(args)
        )

    end

    return nothing

end