text:
module TestStructure

using Base.Test

using Bio.Structure
using TestFunctions.get_bio_fmt_specimens
using Bio.Structure:
    fixlists!,
    parsestrict,
    parselenient,
    parsevalue,
    spacestring

get_bio_fmt_specimens()

# Access a PDB file in BioFmtSpecimens
function pdbfilepath(filename::AbstractString)
    return joinpath(dirname(@__FILE__), "..", "BioFmtSpecimens", "PDB", filename)
end

@testset "Model" begin
    # Test constructors and indexing
    struc = ProteinStructure("Test structure")
    struc[1] = Model(1, struc)
    mod = struc[1]
    @test isa(mod, Model)
    struc[3] = Model(3, struc)
    struc[1]['A'] = Chain('A', mod)
    ch = struc[1]['A']
    @test isa(ch, Chain)
    struc['B'] = Chain('B', mod)
    struc['A'][10] = Residue("ALA", 10, ' ', false, ch)
    res = struc['A'][10]
    @test isa(res, Residue)
    struc['A']["H_20A"] = DisorderedResidue(Dict(
        " VA"=> Residue(" VA", 20, 'A', true, ch),
        "ILE"=> Residue("ILE", 20, 'A', true, ch)
    ), " VA")
    dis_res = struc['A']["H_20A"]
    @test isa(dis_res, DisorderedResidue)
    struc['A'][10][" CA "] = Atom(
        100, " CA ", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res)
    at = struc['A'][10]["CA"]
    @test isa(at, Atom)
    struc['A'][10][" CB "] = DisorderedAtom(Dict(
        'A'=> Atom(200, " CB ", 'A', [10.0, 20.0, 30.0], 0.6, 20.0, " C", "1+", res),
        'B'=> Atom(201, " CB ", 'B', [11.0, 21.0, 31.0], 0.4, 30.0, " C", "1+", res)
    ), 'A')
    dis_at = struc['A'][10]["CB"]
    @test isa(dis_at, DisorderedAtom)
    struc['A']["H_20A"][" CG "] = Atom(
        300, " CG ", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", defaultresidue(dis_res))
    disorderedres(dis_res, "ILE")[" O  "] = Atom(
        400, " O  ", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " O", "  ", disorderedres(dis_res, "ILE"))
    fixlists!(struc)

    # Test alternate constructors
    ProteinStructure("struc", Dict(1=> Model()))
    ProteinStructure()

    Model(1, Dict('A'=> Chain('A')), ProteinStructure())
    Model(1, ProteinStructure())
    Model(1)
    Model()

    Chain('A', ["1"], Dict("1"=> Residue("ALA", 1, ' ', false, Chain('A'))), Model())
    Chain('A', Model())
    Chain('A')

    Residue("ALA", 1, ' ', false, ["CA"], Dict("CA"=>
        Atom(1, "CA", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, "  ", "  ",
        Residue("ALA", 1, ' ', false, Chain('A')))), Chain('A'))
    Residue("ALA", 1, ' ', false, Chain('A'))

    # Test show
    show(DevNull, at)
    show(DevNull, dis_at)
    show(DevNull, res)
    show(DevNull, dis_res)
    show(DevNull, ch)
    show(DevNull, mod)
    show(DevNull, struc)

    # Test getters/setters
    @test serial(at) == 100
    @test serial(dis_at) == 200

    @test atomname(at) == "CA"
    @test atomname(dis_at) == "CB"
    @test atomname(at, strip=false) == " CA "
    @test atomname(dis_at, strip=false) == " CB "

    @test altlocid(at) == ' '
    @test altlocid(dis_at) == 'A'

    @test x(at) == 1.0
    @test x(dis_at) == 10.0

    x!(at, 5.0)
    @test coords(at) == [5.0, 2.0, 3.0]
    # Only the coordinates on the default atom are changed
    x!(dis_at, 12.0)
    @test coords(dis_at) == [12.0, 20.0, 30.0]
    @test coords(dis_at['B']) == [11.0, 21.0, 31.0]

    @test y(at) == 2.0
    @test y(dis_at) == 20.0

    y!(at, 6.0)
    @test coords(at) == [5.0, 6.0, 3.0]
    y!(dis_at, 22.0)
    @test coords(dis_at) == [12.0, 22.0, 30.0]
    @test coords(dis_at['B']) == [11.0, 21.0, 31.0]

    @test z(at) == 3.0
    @test z(dis_at) == 30.0

    z!(at, 7.0)
    @test coords(at) == [5.0, 6.0, 7.0]
    z!(dis_at, 32.0)
    @test coords(dis_at) == [12.0, 22.0, 32.0]
    @test coords(dis_at['B']) == [11.0, 21.0, 31.0]

    coords!(at, [40.0, 50.0, 60.0])
    @test coords(at) == [40.0, 50.0, 60.0]
    coords!(dis_at, [40.0, 50.0, 60.0])
    @test coords(dis_at) == [40.0, 50.0, 60.0]
    @test coords(dis_at['B']) == [11.0, 21.0, 31.0]
    @test_throws ArgumentError coords!(at, [40.0, 50.0, 60.0, 70.0])
    x!(dis_at['A'], 100.0)
    @test coords(dis_at) == [100.0, 50.0, 60.0]
    @test coords(dis_at['B']) == [11.0, 21.0, 31.0]
    x!(dis_at['B'], 110.0)
    @test coords(dis_at) == [100.0, 50.0, 60.0]
    @test coords(dis_at['B']) == [110.0, 21.0, 31.0]

    @test occupancy(at) == 1.0
    @test occupancy(dis_at) == 0.6

    @test tempfactor(at) == 10.0
    @test tempfactor(dis_at) == 20.0

    @test element(at) == "C"
    @test element(dis_at) == "C"
    @test element(at, strip=false) == " C"
    @test element(dis_at, strip=false) == " C"

    @test charge(at) == ""
    @test charge(dis_at) == "1+"
    @test charge(at, strip=false) == "  "
    @test charge(dis_at, strip=false) == "1+"

    @test residue(at) == res
    @test residue(dis_at) == res
    @test residue(res) == res
    @test residue(dis_res) == dis_res

    @test !ishetero(res)
    @test !ishetero(at)
    @test !ishetero(dis_at)
    @test ishetero(dis_res)

    @test !isdisorderedatom(at)
    @test isdisorderedatom(dis_at)

    @test defaultaltlocid(dis_at) == 'A'

    dis_at_mod = DisorderedAtom(dis_at, 'B')
    @test defaultaltlocid(dis_at_mod) == 'B'
    @test serial(dis_at_mod) == 201
    @test_throws ArgumentError DisorderedAtom(dis_at, 'C')

    @test isa(defaultatom(dis_at), Atom)
    @test serial(defaultatom(dis_at)) == 200

    @test altlocids(at) == [' ']
    @test altlocids(dis_at) == ['A', 'B']

    @test atomid(at) == ("10:A", "ALA", "CA")
    @test atomid(dis_at) == ("10:A", "ALA", "CB")

    @test resname(res) == "ALA"
    @test resname(at) == "ALA"
    @test resname(dis_at) == "ALA"
    @test resname(dis_res) == "VA"
    @test resname(dis_res, strip=false) == " VA"

    @test resnumber(res) == 10
    @test resnumber(at) == 10
    @test resnumber(dis_at) == 10
    @test resnumber(dis_res) == 20

    @test inscode(res) == ' '
    @test inscode(at) == ' '
    @test inscode(dis_at) == ' '
    @test inscode(dis_res) == 'A'

    @test resid(at) == "10"
    @test resid(dis_at) == "10"
    @test resid(at, full=true) == "10:A"
    @test resid(dis_at, full=true) == "10:A"
    @test resid(res) == "10"
    @test resid(res, full=true) == "10:A"
    @test resid(dis_res) == "H_20A"
    @test resid(dis_res, full=true) == "H_20A:A"

    @test atomnames(res) == ["CA", "CB"]
    @test atomnames(dis_res) == ["CG"]
    @test atomnames(res, strip=false) == [" CA ", " CB "]
    @test atomnames(dis_res, strip=false) == [" CG "]

    @test isa(atoms(res), Dict{String, AbstractAtom})
    @test length(atoms(res)) == 2
    @test serial(atoms(res)[" CA "]) == 100
    @test isa(atoms(dis_res), Dict{String, AbstractAtom})
    @test length(atoms(dis_res)) == 1
    @test serial(atoms(dis_res)[" CG "]) == 300

    @test !isdisorderedres(res)
    @test isdisorderedres(dis_res)

    @test isa(disorderedres(dis_res, "ILE"), AbstractResidue)
    @test resname(disorderedres(dis_res, "ILE")) == "ILE"

    @test defaultresname(dis_res) == " VA"

    @test isa(defaultresidue(dis_res), Residue)
    @test resname(defaultresidue(dis_res)) == "VA"

    @test resnames(res) == ["ALA"]
    @test resnames(dis_res) == [" VA", "ILE"]

    dis_res_mod = DisorderedResidue(dis_res, "ILE")
    @test defaultresname(dis_res_mod) == "ILE"
    @test atomnames(dis_res_mod) == ["O"]
    @test_throws ArgumentError DisorderedResidue(dis_res, "SER")

    @test chain(at) == ch
    @test chain(dis_at) == ch
    @test chain(res) == ch
    @test chain(dis_res) == ch
    @test chain(ch) == ch

    @test chainid(at) == 'A'
    @test chainid(dis_at) == 'A'
    @test chainid(res) == 'A'
    @test chainid(dis_res) == 'A'
    @test chainid(ch) == 'A'

    @test resids(ch) == ["10", "H_20A"]

    @test isa(residues(ch), Dict{String, AbstractResidue})
    @test length(residues(ch)) == 2
    @test serial(residues(ch)["10"]["CA"]) == 100

    @test model(at) == mod
    @test model(dis_at) == mod
    @test model(res) == mod
    @test model(dis_res) == mod
    @test model(ch) == mod
    @test model(mod) == mod

    @test modelnumber(at) == 1
    @test modelnumber(dis_at) == 1
    @test modelnumber(res) == 1
    @test modelnumber(dis_res) == 1
    @test modelnumber(ch) == 1
    @test modelnumber(mod) == 1

    @test chainids(mod) == ['A', 'B']
    @test chainids(struc) == ['A', 'B']

    @test isa(chains(mod), Dict{Char, Chain})
    @test length(chains(mod)) == 2
    @test resname(chains(mod)['A']["H_20A"]) == "VA"
    @test length(chains(struc)) == 2

    @test structure(at) == struc
    @test structure(dis_at) == struc
    @test structure(res) == struc
    @test structure(dis_res) == struc
    @test structure(ch) == struc
    @test structure(mod) == struc
    @test structure(struc) == struc

    @test structurename(at) == "Test structure"
    @test structurename(dis_at) == "Test structure"
    @test structurename(res) == "Test structure"
    @test structurename(dis_res) == "Test structure"
    @test structurename(ch) == "Test structure"
    @test structurename(mod) == "Test structure"
    @test structurename(struc) == "Test structure"

    @test modelnumbers(struc) == [1, 3]

    @test isa(models(struc), Dict{Int, Model})
    @test length(models(struc)) == 2
    @test resids(models(struc)[1]['A']) == ["10", "H_20A"]

    @test isa(defaultmodel(struc), Model)
    @test modelnumber(defaultmodel(struc)) == 1

    # Test iteration over elements
    at_col = collect(at)
    @test isa(at_col, Vector{Atom})
    @test length(at_col) == 1
    @test serial(at_col[1]) == 100

    dis_at_col = collect(dis_at)
    @test isa(dis_at_col, Vector{Atom})
    @test length(dis_at_col) == 2
    @test serial(dis_at_col[2]) == 201

    res_col = collect(res)
    @test isa(res_col, Vector{AbstractAtom})
    @test length(res_col) == 2
    @test serial(res_col[2]) == 200

    dis_res_col = collect(dis_res)
    @test isa(dis_res_col, Vector{AbstractAtom})
    @test length(dis_res_col) == 1
    @test serial(dis_res_col[1]) == 300

    ch_col = collect(ch)
    @test isa(ch_col, Vector{AbstractResidue})
    @test length(ch_col) == 2
    @test resname(ch_col[2]) == "VA"

    mod_col = collect(mod)
    @test isa(mod_col, Vector{Chain})
    @test length(mod_col) == 2
    @test chainid(mod_col[2]) == 'B'

    struc_col = collect(struc)
    @test isa(struc_col, Vector{Model})
    @test length(struc_col) == 2
    @test modelnumber(struc_col[2]) == 3

    # Test element indices
    @test isa(dis_at['A'], Atom)
    @test serial(dis_at['A']) == 200
    @test serial(dis_at['B']) == 201
    @test_throws KeyError dis_at['C']
    @test_throws MethodError dis_at["A"]

    @test isa(res["CA"], AbstractAtom)
    @test serial(res["CA"]) == 100
    @test serial(res["CB"]) == 200
    @test_throws KeyError res["N"]

    @test isa(dis_res["CG"], AbstractAtom)
    @test serial(dis_res["CG"]) == 300
    @test_throws KeyError dis_res["N"]

    @test isa(ch["10"], AbstractResidue)
    @test isa(ch[10], AbstractResidue)
    @test_throws KeyError ch["H_10"]
    @test serial(ch["10"]["CA"]) == 100
    @test serial(ch[10]["CA"]) == 100
    @test resname(ch["H_20A"]) == "VA"
    @test_throws KeyError ch["11"]
    @test_throws KeyError ch[11]

    @test isa(mod['A'], Chain)
    @test resname(mod['A']["H_20A"]) == "VA"
    @test_throws KeyError mod['C']
    @test_throws MethodError mod["A"]

    @test isa(struc[1], Model)
    @test isa(struc[3], Model)
    @test isa(struc['A'], Chain)
    @test_throws KeyError struc[2]
    @test_throws KeyError struc['C']
    @test_throws MethodError struc["A"]
    @test ishetero(struc[1]['A']["H_20A"])

    # Test selector functions
    ch_a = Chain('A')
    ch_a["10"] = Residue("ALA", 10, ' ', false, ch_a)
    res_a = ch_a["10"]
    ch_a["H_11"] = Residue("MG", 11, ' ', true, ch_a)
    res_b = ch_a["H_11"]
    ch_a["H_100"] = Residue("HOH", 100, ' ', true, ch_a)
    res_c = ch_a["H_100"]
    ch_a["10"]["CA"] = Atom(
        100, " CA ", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res_a)
    at_a = ch_a["10"]["CA"]
    ch_a["H_11"]["MG"] = Atom(
        110, "MG", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res_b)
    at_b = ch_a["H_11"]["MG"]

    @test standardselector(at_a)
    @test !standardselector(at_b)
    @test standardselector(res_a)
    @test !standardselector(res_b)
    @test !heteroselector(at_a)
    @test heteroselector(at_b)
    @test !heteroselector(res_a)
    @test heteroselector(res_b)
    @test atomnameselector(at_a, Set(["CA", "N", "C"]))
    @test atomnameselector(at_a, ["CA", "N", "C"])
    @test !atomnameselector(at_b, Set(["CA", "N", "C"]))
    @test !atomnameselector(at_a, ["CA", "N", "C"], strip=false)
    @test calphaselector(at_a)
    @test !calphaselector(at_b)
    @test !calphaselector(Atom(
        100, "CA", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res_b))
    @test !cbetaselector(at_a)
    @test cbetaselector(Atom(
        100, "CB", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " H", "  ", res_a))
    @test backboneselector(at_a)
    @test !backboneselector(at_b)
    @test !backboneselector(Atom(
        100, "CA", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res_b))
    @test heavyatomselector(at_a)
    @test !heavyatomselector(at_b)
    @test !heavyatomselector(Atom(
        100, "H1", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " H", "  ", res_a))
    @test resnameselector(at_a, Set(["ALA"]))
    @test resnameselector(at_a, ["ALA"])
    @test !resnameselector(at_b, Set(["ALA"]))
    @test resnameselector(res_a, Set(["ALA"]))
    @test !resnameselector(res_b, Set(["ALA"]))
    @test !waterselector(res_a)
    @test waterselector(res_c)
    @test !waterselector(at_a)
    @test notwaterselector(res_a)
    @test !notwaterselector(res_c)
    @test notwaterselector(at_a)
    @test !disorderselector(at_a)
    @test disorderselector(dis_at)
    @test !disorderselector(res_a)
    @test disorderselector(dis_res)
    @test hydrogenselector(Atom(
        100, "H", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " H", "  ", res_a))
    @test !hydrogenselector(Atom(
        100, "H", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res_a))
    @test hydrogenselector(Atom(
        100, "H1", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, "  ", "  ", res_a))
    @test hydrogenselector(Atom(
        100, "1H", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, "  ", "  ", res_a))
    @test !hydrogenselector(Atom(
        100, "NH1", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, "  ", "  ", res_a))

    # Further tests for structural element ordering
    # Order when looping over a DisorderedAtom is the atom serial
    res = Residue("ALA", 1, ' ', false, Chain('A'))
    dis_at_ord = DisorderedAtom(Dict(
        'A' => Atom(102, "CA", 'A', [1.0, 2.0, 3.0], 0.3, 10.0, "C", "", res),
        'B' => Atom(101, "CA", 'B', [1.0, 2.0, 3.0], 0.4, 10.0, "C", "", res),
        'C' => Atom(100, "CA", 'C', [1.0, 2.0, 3.0], 0.3, 10.0, "C", "", res),
    ), 'B')
    @test altlocids(dis_at_ord) == ['C', 'B', 'A']

    # Order when sorting an atom list is the atom serial
    at_list_ord = AbstractAtom[
        DisorderedAtom(Dict(
            'A' => Atom(100, "CA", 'A', [1.0, 2.0, 3.0], 0.4, 10.0, "C", "", res),
            'B' => Atom(104, "CA", 'B', [1.0, 2.0, 3.0], 0.6, 10.0, "C", "", res)
        ), 'B'),
        Atom(102, "CB", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, "C", "", res),
        Atom(103, "CG", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, "C", "", res)
    ]
    @test map(atomname, sort(at_list_ord)) == ["CB", "CG", "CA"]
    sort!(at_list_ord)
    @test map(atomname, at_list_ord) == ["CB", "CG", "CA"]

    # Order when sorting a residue list is chain ID, then stdres/hetres,
    # then residue number, then insertion code
    res_ord = AbstractResidue[
        Residue("ALA", 201, 'A', false, Chain('A')),
        Residue("ALA", 203, ' ', false, Chain('A')),
        Residue("ALA", 200, ' ', true,  Chain('A')),
        Residue("ALA", 201, 'B', false, Chain('A')),
        Residue("ALA", 202, ' ', false, Chain('A')),
        Residue("ALA", 300, ' ', false, Chain('B')),
        Residue("ALA", 201, ' ', true,  Chain('A')),
        Residue("ALA", 201, ' ', false, Chain('A')),
        Residue("ALA", 201, 'A', true,  Chain('A')),
        Residue("ALA", 100, ' ', false, Chain('B')),
        Residue("ALA", 203, ' ', true,  Chain('A')),
        Residue("ALA", 200, ' ', false, Chain('A')),
    ]

    # Test sequentialresidues
    @test sequentialresidues(res_ord[3], res_ord[7])
    @test !sequentialresidues(res_ord[7], res_ord[3])
    @test sequentialresidues(res_ord[7], res_ord[9])

    @test map(resid, res_ord) == [
        "201A", "203", "H_200", "201B", "202", "300", "H_201", "201", "H_201A",
        "100", "H_203", "200"]
    @test map(res -> resid(res, full=true), res_ord) == [
        "201A:A", "203:A", "H_200:A", "201B:A", "202:A", "300:B", "H_201:A",
        "201:A", "H_201A:A", "100:B", "H_203:A", "200:A"]
    @test map(res -> resid(res, full=true), sort(res_ord)) == [
        "200:A", "201:A", "201A:A", "201B:A", "202:A", "203:A", "H_200:A",
        "H_201:A", "H_201A:A", "H_203:A", "100:B", "300:B"]
    sort!(res_ord)
    @test map(res -> resid(res, full=true), res_ord) == [
        "200:A", "201:A", "201A:A", "201B:A", "202:A", "203:A", "H_200:A",
        "H_201:A", "H_201A:A", "H_203:A", "100:B", "300:B"]

    # Order of listing residue names in a DisorderedResidue is default then alphabetical
    dis_res_ord = DisorderedResidue(Dict(
        "THR" => Residue("THR", 201, ' ', false, Chain('A')),
        "ALA" => Residue("ALA", 201, ' ', false, Chain('A')),
        "ILE" => Residue("ILE", 201, ' ', false, Chain('A')),
        "SER" => Residue("SER", 201, ' ', false, Chain('A')),
        "VAL" => Residue("VAL", 201, ' ', false, Chain('A'))
    ), "SER")
    @test defaultresname(dis_res_ord) == "SER"
    @test resnames(dis_res_ord) == ["SER", "ALA", "ILE", "THR", "VAL"]

    # Order when sorting chain IDs is character ordering with the empty chain ID at the end
    mod_ord = Model(1, Dict(
        'A' => Chain('A'),
        ' ' => Chain(' '),
        '1' => Chain('1'),
        'a' => Chain('a'),
        'X' => Chain('X'),
    ), ProteinStructure())
    @test chainids(mod_ord) == ['1', 'A', 'X', 'a', ' ']

    # Test sequence extraction
    struc = read(pdbfilepath("1AKE.pdb"), PDB)
    seq = AminoAcidSequence(struc['B'])
    @test seq == AminoAcidSequence(
        "MRIILLGAPGAGKGTQAQFIMEKYGIPQISTGDMLRAAVKSGSELGKQAKDIMDAGKLVTDELVIALVKERIAQEDCRNG" *
        "FLLDGFPRTIPQADAMKEAGINVDYVLEFDVPDELIVDRIVGRRVHAPSGRVYHVKFNPPKVEGKDDVTGEELTTRKDDQ" *
        "EETVRKRLVEYHQMTAPLIGYYSKEAEAGNTKYAKVDGTKPVAEVRADLEKILGX-------------------------" *
        "--------------------------------------------------------------------------------" *
        "--------------------------------------------------------------------------------" *
        "--------------------------------------------------------------------------------" *
        "---------------------------------------------------XX---------------------------" *
        "----------------------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" *
        "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX" *
        "XXXXXXXXXXXXXXX"
    )
    seq = AminoAcidSequence(struc['B'], standardselector)
    @test seq == AminoAcidSequence(
        "MRIILLGAPGAGKGTQAQFIMEKYGIPQISTGDMLRAAVKSGSELGKQAKDIMDAGKLVTDELVIALVKERIAQEDCRNG" *
        "FLLDGFPRTIPQADAMKEAGINVDYVLEFDVPDELIVDRIVGRRVHAPSGRVYHVKFNPPKVEGKDDVTGEELTTRKDDQ" *
        "EETVRKRLVEYHQMTAPLIGYYSKEAEAGNTKYAKVDGTKPVAEVRADLEKILG"
    )
    seq = AminoAcidSequence(AbstractResidue[
        Residue("VAL", 20, 'B', true, Chain('B')),
        Residue("ALA", 10, 'A', false, Chain('A')),
    ])
    @test seq == AminoAcidSequence("VA")
end

@testset "Parsing" begin
    # Test parsevalue
    line = "ATOM     40  CB  LEU A   5      22.088  45.547  29.675  1.00 22.23           C  "
    @test parsevalue(line, 7, 11, Int) == 40
    @test parsevalue(line, 31, 38, Float64) == 22.088
    @test parsevalue(line, 13, 16, String) == " CB "
    @test parsevalue(line, 22, 22, Char) == 'A'
    @test_throws ErrorException parsevalue(line, 1, 4, Int)
    @test_throws ErrorException parsevalue(line, 1, 4, Bool)
    @test_throws ErrorException parsevalue(line, 79, 100, Int)

    # Test parsestrict
    line =   "ATOM    591  C   GLY A  80      29.876  54.131  35.806  1.00 40.97           C  "
    line_a = "ATOM    591  C   GLY A  80              54.131  35.806  1.00 40.97           C  "
    line_b = "ATOM    591  C   GLY A  80 "
    @test parsestrict(line, 7, 11, Int, "could not read atom serial number", 10) == 591
    @test parsestrict(line, 13, 16, String, "could not read atom name", 20) == " C  "
    @test_throws PDBParseError parsestrict(line_a, 31, 38, Float64, "could not read x coordinate", 10)
    @test_throws PDBParseError parsestrict(line_b, 31, 38, Float64, "could not read x coordinate", 10)
    @test_throws PDBParseError parsestrict(line, 7, 11, Bool, "could not read atom serial number", 10)

    # Test parselenient
    line =   "ATOM     40  CB  LEU A   5      22.088  45.547  29.675  1.00 22.23           C  "
    line_a = "ATOM     40  CB  LEU A   5      22.088  45.547  29.675  1.00 22.23              "
    line_b = "ATOM     40  CB  LEU A   5      22.088  45.547  29.675  1.00 22.23  "
    line_c = "ATOM     40  CB  LEU A   5      22.088  45.547  29.675       22.23           C  "
    @test parselenient(line, 77, 78, String, "  ") == " C"
    @test parselenient(line_a, 77, 78, String, "  ") == "  "
    @test parselenient(line_b, 77, 78, String, "  ") == "  "
    @test parselenient(line_b, 77, 78, String, " C") == " C"
    @test parselenient(line_c, 55, 60, Float64, 1.0) == 1.0
    @test parselenient(line, 77, 78, Bool, " N") == " N"

    # Test AtomRecord constructor
    line_a = "ATOM    669  CA  ILE A  90      31.743  33.110  31.221  1.00 25.76           C  "
    line_b = "HETATM 3474  O  B XX A 334A      8.802  62.000   8.672  1.00 39.15           O1-"
    line_c = "ATOM    669  CA  ILE A  90      xxxxxx  33.110  31.221  1.00 25.76           C  "
    line_d = "ATOM    669  CA  ILE A  90      31.743   "
    line_e = "REMARK   1 REFERENCE 1                                                          "
    at_rec = AtomRecord(line_a, 10)
    show(DevNull, at_rec)
    @test !at_rec.het_atom
    @test at_rec.serial == 669
    @test at_rec.atom_name == " CA "
    @test at_rec.alt_loc_id == ' '
    @test at_rec.res_name == "ILE"
    @test at_rec.chain_id == 'A'
    @test at_rec.res_number == 90
    @test at_rec.ins_code == ' '
    @test at_rec.coords == [31.743, 33.110, 31.221]
    @test at_rec.occupancy == 1.00
    @test at_rec.temp_factor == 25.76
    @test at_rec.element == " C"
    @test at_rec.charge == "  "
    at_rec = AtomRecord(line_b)
    @test at_rec.het_atom
    @test at_rec.serial == 3474
    @test at_rec.atom_name == " O  "
    @test at_rec.alt_loc_id == 'B'
    @test at_rec.res_name == " XX"
    @test at_rec.chain_id == 'A'
    @test at_rec.res_number == 334
    @test at_rec.ins_code == 'A'
    @test at_rec.coords == [8.802, 62.0, 8.672]
    @test at_rec.occupancy == 1.00
    @test at_rec.temp_factor == 39.15
    @test at_rec.element == " O"
    @test at_rec.charge == "1-"
    @test_throws PDBParseError AtomRecord(line_c)
    @test_throws PDBParseError AtomRecord(line_d)

    # Test parsing 1AKE (multiple chains, disordered atoms)
    struc = read(pdbfilepath("1AKE.pdb"), PDB)
    @test structurename(struc) == "1AKE.pdb"
    @test countmodels(struc) == 1
    @test modelnumbers(struc) == [1]
    @test countchains(struc) == 2
    @test countchains(struc[1]) == 2
    @test chainids(struc) == ['A', 'B']
    @test chainids(struc[1]) == ['A', 'B']
    @test resname(struc['A'][10]) == "GLY"
    @test !isdisorderedres(struc['A'][10])
    @test serial(struc['A'][200]["NZ"]) == 1555
    @test !isdisorderedatom(struc['A'][200]["NZ"])
    @test isdisorderedatom(struc['A'][167]["CD"])
    @test altlocids(struc['A'][167]["CD"]) == ['A', 'B']
    @test x(struc['A'][167]["CD"]) == 24.502
    @test x(struc['A'][167]["CD"]['A']) == 24.502
    @test x(struc['A'][167]["CD"]['B']) == 24.69
    mods = collect(struc)
    @test modelnumber(mods[1]) == 1
    chs = collect(mods[1])
    @test map(chainid, chs) == ['A', 'B']
    res = collect(chs[1])
    @test length(res) == 456
    @test resid(res[20]) == "20"
    ats = collect(res[20])
    @test length(ats) == 8
    @test map(atomname, ats) == ["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"]

    # Test choosedefaultaltlocid
    res = Residue("ALA", 1, ' ', false, Chain('A'))
    at_a = Atom(100, "CA", 'A', [1.0, 2.0, 3.0], 0.4, 10.0, "C", "", res)
    at_b = Atom(101, "CA", 'B', [1.0, 2.0, 3.0], 0.6, 10.0, "C", "", res)
    @test choosedefaultaltlocid(at_a, at_b) == 'B'
    @test choosedefaultaltlocid(at_b, at_a) == 'B'
    at_a = Atom(100, "CA", 'A', [1.0, 2.0, 3.0], 0.5, 10.0, "C", "", res)
    at_b = Atom(101, "CA", 'B', [1.0, 2.0, 3.0], 0.5, 10.0, "C", "", res)
    @test choosedefaultaltlocid(at_a, at_b) == 'A'
    @test choosedefaultaltlocid(at_b, at_a) == 'A'

    # Test applyselectors
    ats = collectatoms(struc)
    # Not providing any selector functions just returns the input list
    ats_min = applyselectors(ats)
    @test length(ats_min) == length(ats)
    @test map(serial, ats_min) == map(serial, ats)
    applyselectors!(ats_min)
    @test length(ats_min) == length(ats)
    @test map(serial, ats_min) == map(serial, ats)
    ats_min = applyselectors(ats, standardselector)
    @test length(ats_min) == 3312
    @test serial(ats_min[2000]) == 2006
    applyselectors!(ats, standardselector)
    @test length(ats) == 3312
    @test serial(ats[2000]) == 2006
    ats = collectatoms(struc)
    ats_min = applyselectors(ats, standardselector, disorderselector)
    @test length(ats_min) == 5
    @test serial(ats_min[4]) == 1294
    applyselectors!(ats, standardselector, disorderselector)
    @test length(ats) == 5
    @test serial(ats[4]) == 1294

    res = collectresidues(struc)
    res_min = applyselectors(res)
    @test length(res_min) == length(res)
    @test map(res -> resid(res, full=true), res_min) == map(res -> resid(res, full=true), res)
    applyselectors!(res_min)
    @test length(res_min) == length(res)
    @test map(res -> resid(res, full=true), res_min) == map(res -> resid(res, full=true), res)
    res_min = applyselectors(res, waterselector)
    @test length(res_min) == 378
    @test resid(res_min[300], full=true) == "H_657:B"
    applyselectors!(res, waterselector)
    @test length(res) == 378
    @test resid(res[300], full=true) == "H_657:B"
    res = collectresidues(struc)
    # Test anonymous selector function
    res_min = applyselectors(res, standardselector, res -> chainid(res) == 'A')
    @test length(res_min) == 214
    @test resid(res_min[200], full=true) == "200:A"
    applyselectors!(res, standardselector, res -> chainid(res) == 'A')
    @test length(res) == 214
    @test resid(res[200], full=true) == "200:A"

    # Test parsing options
    struc = read(pdbfilepath("1AKE.pdb"), PDB, structure_name="New name")
    @test structurename(struc) == "New name"
    @test countatoms(struc) == 3804

    struc = read(pdbfilepath("1AKE.pdb"), PDB, read_het_atoms=false)
    @test countatoms(struc) == 3312
    @test serial(collectatoms(struc)[2000]) == 2006
    @test sum(map(ishetero, collectatoms(struc))) == 0

    struc = read(pdbfilepath("1AKE.pdb"), PDB, read_std_atoms=false)
    @test countatoms(struc) == 492
    @test serial(collectatoms(struc)[400]) == 3726
    @test sum(map(ishetero, collectatoms(struc))) == 492

    struc = read(pdbfilepath("1AKE.pdb"), PDB, read_het_atoms=false, read_std_atoms=false)
    @test countatoms(struc) == 0
    @test countresidues(struc) == 0
    @test countchains(struc) == 0
    @test countmodels(struc) == 0

    struc = read(pdbfilepath("1AKE.pdb"), PDB, remove_disorder=true)
    @test countatoms(struc) == 3804
    @test sum(map(isdisorderedatom, collectatoms(struc))) == 0
    @test tempfactor(struc['A'][167]["NE"]) == 23.32

    # Test parsing from stream
    open(pdbfilepath("1AKE.pdb"), "r") do file
        struc = read(file, PDB)
        @test countatoms(struc) == 3804
        @test countresidues(struc) == 808
    end

    # Test parsing 1EN2 (disordered residue)
    struc = read(pdbfilepath("1EN2.pdb"), PDB)
    @test modelnumbers(struc) == [1]
    @test chainids(struc[1]) == ['A']
    @test serial(struc['A'][48]["CA"]) == 394
    @test isdisorderedres(struc['A'][10])
    @test defaultresname(struc['A'][10]) == "SER"
    @test resname(disorderedres(struc['A'][10], "SER")) == "SER"
    @test resname(disorderedres(struc['A'][10], "GLY")) == "GLY"
    # Atoms in a disordered residue are not necessarily disordered
    @test !isdisorderedatom(struc['A'][10]["CA"])
    @test altlocid(struc['A'][10]["CA"]) == 'A'
    @test isdisorderedres(struc['A'][16])
    @test defaultresname(struc['A'][16]) == "ARG"
    @test resname(disorderedres(struc['A'][16], "ARG")) == "ARG"
    @test resname(disorderedres(struc['A'][16], "TRP")) == "TRP"
    @test !isdisorderedatom(disorderedres(struc['A'][16], "TRP")["CA"])
    @test altlocid(disorderedres(struc['A'][16], "TRP")["CA"]) == 'C'
    @test isdisorderedatom(struc['A'][16]["CA"])
    @test altlocids(struc['A'][16]["CA"]) == ['A', 'B']
    @test defaultaltlocid(struc['A'][16]["CA"]) == 'A'
    @test occupancy(struc['A'][16]["CA"]) == 0.22
    ats = collectatoms(DisorderedResidue[struc['A'][10], struc['A'][16]])
    @test length(ats) == 17
    @test isa(ats, Vector{AbstractAtom})
    @test serial(ats[10]) == 113
    @test isa(ats[10], DisorderedAtom)
    @test countatoms(DisorderedResidue[struc['A'][10], struc['A'][16]]) == 17
    res = collectresidues(DisorderedResidue[struc['A'][16], struc['A'][10]])
    @test length(res) == 2
    @test isa(res, Vector{DisorderedResidue})
    @test resnumber(res[1]) == 10
    @test countresidues(DisorderedResidue[struc['A'][16], struc['A'][10]]) == 2

    # Test parsing 1SSU (multiple models)
    struc = read(pdbfilepath("1SSU.pdb"), PDB)
    # Test countmodels
    @test countmodels(struc) == 20
    @test modelnumbers(struc) == collect(1:20)
    @test countchains(struc) == 1
    @test countchains(struc[5]) == 1
    @test chainids(struc) == ['A']
    @test chainids(struc[5]) == ['A']
    @test serial(struc[10]['A'][40]["HB2"]) == 574
    @test countatoms(struc) == 756
    @test map(countatoms, [struc[i] for i in modelnumbers(struc)]) == 756 * ones(Int, 20)
    @test countatoms(struc, hydrogenselector) == 357
    ats = collectatoms(Model[struc[5], struc[10]])
    @test length(ats) == 1512
    @test z(ats[20]) == -16.522
    @test z(ats[1000]) == -0.394
    @test countatoms(Model[struc[5], struc[10]]) == 1512
    res = collectresidues(Model[struc[5], struc[10]])
    @test length(res) == 102
    @test y(res[10]["O"]) == -6.421
    @test y(res[100]["O"]) == -15.66
    @test countresidues(Model[struc[5], struc[10]]) == 102
    chs = collectchains(Model[struc[5], struc[10]])
    @test length(chs) == 2
    @test map(chainid, chs) == ['A', 'A']
    @test z(chs[2][5]["CA"]) == -5.667
    @test countchains(Model[struc[5], struc[10]]) == 2
    mods = collectmodels(Model[struc[10], struc[5]])
    @test length(mods) == 2
    @test map(modelnumber, mods) == [5, 10]
    @test z(mods[2]['A'][5]["CA"]) == -5.667
    @test countmodels(Model[struc[10], struc[5]]) == 2

    # Test collectatoms
    struc = read(pdbfilepath("1AKE.pdb"), PDB)
    ats = collectatoms(struc)
    @test length(ats) == 3804
    @test isa(ats, Vector{AbstractAtom})
    @test isa(ats[70], Atom)
    @test isa(ats[1290], DisorderedAtom)
    @test serial(ats[1660]) == 1666
    ats = collectatoms(struc, heteroselector)
    @test length(ats) == 492
    @test serial(ats[80]) == 3406
    ats = collectatoms(struc, disorderselector)
    @test length(ats) == 12
    @test serial(ats[10]) == 3338
    ats = collectatoms(struc, standardselector, disorderselector)
    @test length(ats) == 5
    @test serial(ats[4]) == 1294
    ats = collectatoms(struc[1])
    @test length(ats) == 3804
    @test serial(ats[1660]) == 1666
    ats = collectatoms(struc['A'])
    @test length(ats) == 1954
    @test serial(ats[240]) == 240
    ats = collectatoms(struc['A'][50])
    @test length(ats) == 9
    @test serial(ats[4]) == 358
    ats = collectatoms(struc['A'][50]["CA"])
    @test length(ats) == 1
    @test isa(ats, Vector{AbstractAtom})
    @test serial(ats[1]) == 356
    ats = collectatoms(struc['A'][167]["CZ"])
    @test length(ats) == 1
    @test isa(ats, Vector{AbstractAtom})
    @test isa(ats[1], DisorderedAtom)
    @test serial(ats[1]) == 1292
    ats = collectatoms(Chain[struc['B'], struc['A']])
    @test length(ats) == 3804
    @test serial(ats[5]) == 5
    ats = collectatoms(Residue[struc['A'][51], struc['A'][50]])
    @test length(ats) == 17
    @test serial(ats[5]) == 359
    ats = collectatoms(AbstractResidue[struc['A'][51], struc['A'][50]])
    @test length(ats) == 17
    @test serial(ats[5]) == 359
    ats = collectatoms(Atom[struc['A'][51]["CA"], struc['A'][50]["CA"]])
    @test length(ats) == 2
    @test isa(ats, Vector{Atom})
    @test serial(ats[2]) == 365
    ats = collectatoms(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]])
    @test length(ats) == 2
    @test isa(ats, Vector{DisorderedAtom})
    @test serial(ats[2]) == 1292
    ats = collectatoms(AbstractAtom[struc['A'][50]["CA"], struc['A'][167]["CZ"]])
    @test length(ats) == 2
    @test isa(ats, Vector{AbstractAtom})
    @test serial(ats[2]) == 1292

    # Test countatoms
    @test countatoms(struc) == 3804
    @test countatoms(struc[1]) == 3804
    @test countatoms(struc['A']) == 1954
    @test countatoms(struc['A'][50]) == 9
    @test countatoms(struc['A'][50]["CA"]) == 1
    @test countatoms([struc['A'], struc['B']]) == 3804
    @test countatoms(AbstractResidue[struc['A'][50], struc['A'][51]]) == 17
    @test countatoms(Residue[struc['A'][50], struc['A'][51]]) == 17
    @test countatoms(collectatoms(struc['A'])) == 1954
    @test countatoms(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]]) == 2
    @test countatoms(Atom[struc['A'][51]["CA"], struc['A'][50]["CA"]]) == 2

    @test countatoms(struc['A'], standardselector) == 1656
    @test countatoms(struc['A'], heteroselector) == 298
    @test countatoms(struc['A'], standardselector, disorderselector) == 5

    @test countatoms(ProteinStructure()) == 0
    @test countatoms(Model()) == 0
    @test countatoms(Chain('X')) == 0
    @test countatoms(Residue("ALA", 100, ' ', false, Chain('A'))) == 0

    # Test collectresidues
    res = collectresidues(struc)
    @test length(res) == 808
    @test isa(res, Vector{AbstractResidue})
    @test isa(res[50], Residue)
    @test resnumber(res[220]) == 305
    res = collectresidues(struc, heteroselector)
    @test length(res) == 380
    @test resnumber(res[370]) == 725
    res = collectresidues(struc, standardselector, res -> chainid(res) == 'A')
    @test length(res) == 214
    @test resnumber(res[200]) == 200
    res = collectresidues(struc[1])
    @test length(res) == 808
    @test resnumber(res[220]) == 305
    res = collectresidues(struc['A'])
    @test length(res) == 456
    @test resnumber(res[220]) == 305
    res = collectresidues(struc['A'][50])
    @test length(res) == 1
    @test isa(res, Vector{AbstractResidue})
    @test resnumber(res[1]) == 50
    res = collectresidues(struc['A'][50]["CA"])
    @test length(res) == 1
    @test isa(res, Vector{AbstractResidue})
    @test isa(res[1], Residue)
    @test resnumber(res[1]) == 50
    res = collectresidues(struc['A'][167]["CZ"])
    @test length(res) == 1
    @test isa(res, Vector{AbstractResidue})
    @test isa(res[1], Residue)
    @test resnumber(res[1]) == 167
    res = collectresidues(Chain[struc['B'], struc['A']])
    @test length(res) == 808
    @test resid(res[5], full=true) == "5:A"
    res = collectresidues(Residue[struc['A'][51], struc['A'][50]])
    @test length(res) == 2
    @test isa(res, Vector{Residue})
    @test resnumber(res[1]) == 50
    res = collectresidues(AbstractResidue[struc['A'][51], struc['A'][50]])
    @test length(res) == 2
    @test isa(res, Vector{AbstractResidue})
    @test resnumber(res[1]) == 50
    res = collectresidues(Atom[struc['A'][51]["CA"], struc['A'][50]["CA"]])
    @test length(res) == 2
    @test isa(res, Vector{AbstractResidue})
    @test resnumber(res[2]) == 51
    res = collectresidues(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]])
    @test length(res) == 1
    @test isa(res, Vector{AbstractResidue})
    @test atomnames(res[1]) == ["N", "CA", "C", "O", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2"]
    res = collectresidues(AbstractAtom[struc['A'][50]["CA"], struc['A'][167]["CZ"]])
    @test length(res) == 2
    @test isa(res, Vector{AbstractResidue})
    @test atomnames(res[1]) == ["N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"]

    # Test countresidues
    @test countresidues(struc) == 808
    @test countresidues(struc[1]) == 808
    @test countresidues(struc['A']) == 456
    @test countresidues(struc['A'][50]) == 1
    @test countresidues(struc['A'][50]["CA"]) == 1
    @test countresidues([struc['A'], struc['B']]) == 808
    @test countresidues(AbstractResidue[struc['A'][50], struc['A'][51]]) == 2
    @test countresidues(Residue[struc['A'][50], struc['A'][51]]) == 2
    @test countresidues(collectatoms(struc['A'])) == 456
    @test countresidues(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]]) == 1
    @test countresidues(Atom[struc['A'][51]["CA"], struc['A'][50]["CA"]]) == 2

    @test countresidues(struc['A'], standardselector) == 214
    @test countresidues(struc['A'], heteroselector) == 242
    @test countresidues(struc, standardselector, res -> chainid(res) == 'A') == 214

    @test countresidues(ProteinStructure()) == 0
    @test countresidues(Model()) == 0
    @test countresidues(Chain('X')) == 0
    @test countresidues(Residue("ALA", 100, ' ', false, Chain('A'))) == 1

    # Test collectchains
    chs = collectchains(struc)
    @test length(chs) == 2
    @test isa(chs, Vector{Chain})
    @test chainid(chs[2]) == 'B'
    chs = collectchains(struc, ch -> chainid(ch) == 'B')
    @test length(chs) == 1
    @test chainid(chs[1]) == 'B'
    chs = collectchains(struc[1])
    @test length(chs) == 2
    @test chainid(chs[2]) == 'B'
    chs = collectchains(struc['A'])
    @test length(chs) == 1
    @test chainid(chs[1]) == 'A'
    chs = collectchains(struc['A'][50])
    @test length(chs) == 1
    @test chainid(chs[1]) == 'A'
    chs = collectchains(struc['A'][50]["CA"])
    @test length(chs) == 1
    @test chainid(chs[1]) == 'A'
    chs = collectchains(Chain[struc['B'], struc['A']])
    @test length(chs) == 2
    @test chainid(chs[2]) == 'B'
    chs = collectchains(Residue[struc['A'][51], struc['B'][50]])
    @test length(chs) == 2
    @test chainid(chs[2]) == 'B'
    chs = collectchains(AbstractResidue[struc['A'][51], struc['B'][50]])
    @test length(chs) == 2
    @test chainid(chs[2]) == 'B'
    chs = collectchains(Atom[struc['B'][51]["CA"], struc['A'][50]["CA"]])
    @test length(chs) == 2
    @test chainid(chs[2]) == 'B'
    chs = collectchains(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]])
    @test length(chs) == 1
    @test chainid(chs[1]) == 'A'
    chs = collectchains(AbstractAtom[struc['A'][50]["CA"], struc['A'][167]["CZ"]])
    @test length(chs) == 1
    @test chainid(chs[1]) == 'A'

    # Test countchains
    @test countchains(struc) == 2
    @test countchains(struc[1]) == 2
    @test countchains(struc['A']) == 1
    @test countchains(struc['A'][50]) == 1
    @test countchains(struc['A'][50]["CA"]) == 1
    @test countchains([struc['A'], struc['B']]) == 2
    @test countchains(AbstractResidue[struc['A'][50], struc['B'][51]]) == 2
    @test countchains(Residue[struc['A'][50], struc['B'][51]]) == 2
    @test countchains(collectatoms(struc['A'])) == 1
    @test countchains(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]]) == 1
    @test countchains(Atom[struc['A'][51]["CA"], struc['B'][50]["CA"]]) == 2

    @test countchains(struc, ch -> chainid(ch) == 'B') == 1

    @test countchains(ProteinStructure()) == 0
    @test countchains(Model()) == 0
    @test countchains(Chain('X')) == 1
    @test countchains(Residue("ALA", 100, ' ', false, Chain('A'))) == 1

    # Test collectmodels
    struc_1SSU = read(pdbfilepath("1SSU.pdb"), PDB)
    mods = collectmodels(struc_1SSU)
    @test length(mods) == 20
    @test isa(mods, Vector{Model})
    @test map(modelnumber, mods) == collect(1:20)
    mods = collectmodels(struc_1SSU, mod -> modelnumber(mod) < 4)
    @test length(mods) == 3
    @test modelnumber(mods[2]) == 2
    mods = collectmodels(struc_1SSU[10])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 10
    mods = collectmodels(struc_1SSU['A'])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(struc_1SSU[7]['A'][50])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 7
    mods = collectmodels(struc_1SSU['A'][50]["CA"])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(Chain[struc['B'], struc['A']])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(Residue[struc['A'][51], struc['B'][50]])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(AbstractResidue[struc_1SSU[12]['A'][51], struc_1SSU['A'][50]])
    @test length(mods) == 2
    @test modelnumber(mods[2]) == 12
    mods = collectmodels(Atom[struc_1SSU['A'][51]["CA"], struc_1SSU['A'][50]["CA"]])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1
    mods = collectmodels(AbstractAtom[struc['A'][50]["CA"], struc['A'][167]["CZ"]])
    @test length(mods) == 1
    @test modelnumber(mods[1]) == 1

    # Test countmodels
    @test countmodels(struc_1SSU) == 20
    @test countmodels(struc_1SSU[10]) == 1
    @test countmodels(struc_1SSU['A']) == 1
    @test countmodels(struc_1SSU['A'][50]) == 1
    @test countmodels(struc_1SSU[10]['A'][50]["CA"]) == 1
    @test countmodels([struc['A'], struc['B']]) == 1
    @test countmodels(AbstractResidue[struc_1SSU[1]['A'][50], struc_1SSU[2]['A'][51]]) == 2
    @test countmodels(Residue[struc['A'][50], struc['B'][51]]) == 1
    @test countmodels(collectatoms(struc_1SSU)) == 1
    @test countmodels(collectatoms(collect(struc_1SSU))) == 20
    @test countmodels(DisorderedAtom[struc['A'][167]["CZ"], struc['A'][167]["CD"]]) == 1
    @test countmodels(Atom[struc_1SSU[7]['A'][51]["CA"], struc_1SSU['A'][50]["CA"]]) == 2

    @test countmodels(struc_1SSU, mod -> modelnumber(mod) < 4) == 3

    @test countmodels(ProteinStructure()) == 0
    @test countmodels(Model()) == 1
    @test countmodels(Chain('X')) == 1
    @test countmodels(Residue("ALA", 100, ' ', false, Chain('A'))) == 1

    # Test parser error handling
    error = PDBParseError("message", 10, "line")
    showerror(DevNull, error)
    # Missing coordinate (blank string)
    @test_throws PDBParseError read(pdbfilepath("1AKE_err_a.pdb"), PDB)
    # Missing chain ID (line ends early)
    @test_throws PDBParseError read(pdbfilepath("1AKE_err_b.pdb"), PDB)
    # Bad MODEL record
    @test_throws PDBParseError read(pdbfilepath("1SSU_err.pdb"), PDB)
    # Duplicate atom names in same residue
    @test_throws ErrorException read(pdbfilepath("1AKE_err_c.pdb"), PDB)
    # Non-existent file
    @test_throws SystemError read(pdbfilepath("non_existent_file.pdb"), PDB)
end

@testset "Writing" begin
    # Test spacestring
    @test spacestring(1.5, 5) == "  1.5"
    @test spacestring("A", 3) == "  A"
    @test spacestring('A', 3) == "  A"
    @test_throws ArgumentError spacestring(1.456789, 5)
    @test_throws ArgumentError spacestring("ABCDEF", 3)

    # Test spaceatomname
    res = Residue("ALA", 1, ' ', false, Chain('A'))
    @test spaceatomname(Atom(1, " CA ", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " C", "  ", res)) == " CA "
    @test spaceatomname(Atom(1, "N",    ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " N", "  ", res)) == " N  "
    @test spaceatomname(Atom(1, "N",    ' ', [0.0, 0.0, 0.0], 1.0, 0.0, "  ", "  ", res)) == " N  "
    @test spaceatomname(Atom(1, "CA",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " C", "  ", res)) == " CA "
    @test spaceatomname(Atom(1, "NE1",  ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " N", "  ", res)) == " NE1"
    @test spaceatomname(Atom(1, "2HD1", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res)) == "2HD1"
    @test spaceatomname(Atom(1, "HH11", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res)) == "HH11"
    @test spaceatomname(Atom(1, "1H",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res)) == "1H  "
    @test spaceatomname(Atom(1, "MG",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, "MG", "  ", res)) == "MG  "
    @test spaceatomname(Atom(1, "MG",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, "  ", "  ", res)) == " MG "
    @test_throws ArgumentError spaceatomname(Atom(1, "11H",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res))
    @test_throws ArgumentError spaceatomname(Atom(1, "11H11", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res))
    @test_throws ArgumentError spaceatomname(Atom(1, "1MG",   ' ', [0.0, 0.0, 0.0], 1.0, 0.0, "MG", "  ", res))

    # Test pdbline
    ch_a = Chain('A')
    ch_a["1"] = Residue("ALA", 1, ' ', false, ch_a)
    ch_a["1"][" N  "] = Atom(10, " N  ", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " N", "  ", ch_a["1"])
    line = pdbline(ch_a["1"][" N  "])
    @test line == "ATOM     10  N   ALA A   1         0.0     0.0     0.0   1.0   0.0           N  "
    ch_b = Chain('B')
    ch_b["H_20"] = Residue("X", 20, ' ', true, ch_b)
    ch_b["H_20"]["C"] = Atom(101, "C", 'A', [10.5, 20.12345, -5.1227], 0.50, 50.126, "C", "1+", ch_b["H_20"])
    line = pdbline(ch_b["H_20"]["C"])
    @test line == "HETATM  101  C  A  X B  20        10.5  20.123  -5.123   0.5 50.13           C1+"
    ch_b["H_20"]["11H11"] = Atom(1, "11H11", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", ch_b["H_20"])
    @test_throws ArgumentError pdbline(ch_b["H_20"]["11H11"])

    line_a = "ATOM    669  CA  ILE A  90      31.743  33.110  31.221  1.00 25.76           C  "
    line_b = "HETATM 3474  O  B XX A 334A      8.802  62.000   8.672  1.00 39.15           O1-"
    at_rec = AtomRecord(line_a)
    @test pdbline(at_rec) == "ATOM    669  CA  ILE A  90      31.743   33.11  31.221   1.0 25.76           C  "
    at_rec = AtomRecord(line_b)
    @test pdbline(at_rec) == "HETATM 3474  O  B XX A 334A      8.802    62.0   8.672   1.0 39.15           O1-"

    # Test writepdb
    # Counts lines in a file
    function countlines(filename::AbstractString)
        counter = 0
        open(filename, "r") do file
            for line in eachline(file)
                counter += 1
            end
        end
        return counter
    end

    struc = read(pdbfilepath("1SSU.pdb"), PDB)
    # All writing is done to one temporary file which is removed at the end
    temp_filename = tempname()
    writepdb(temp_filename, struc)
    @test countlines(temp_filename) == 15160
    struc_written = read(temp_filename, PDB)
    @test isa(struc_written, ProteinStructure)
    @test modelnumbers(struc_written) == collect(1:20)
    @test countatoms(struc_written) == 756
    @test z(struc_written[4]['A']["30"]["OG"]) == -2.177
    @test atomnames(struc_written[15]['A']["39"]) == [
        "N", "CA", "C", "O", "CB", "SG", "H", "HA", "HB2", "HB3"]

    # Test writing to stream
    open(temp_filename, "w") do file
        writepdb(file, struc)
    end
    @test countlines(temp_filename) == 15160
    struc_written = read(temp_filename, PDB)
    @test modelnumbers(struc_written) == collect(1:20)
    @test countatoms(struc_written) == 756
    @test z(struc_written[4]['A']["30"]["OG"]) == -2.177
    @test atomnames(struc_written[15]['A']["39"]) == [
        "N", "CA", "C", "O", "CB", "SG", "H", "HA", "HB2", "HB3"]

    # Test selectors
    struc = read(pdbfilepath("1AKE.pdb"), PDB)
    writepdb(temp_filename, struc, heteroselector)
    @test countlines(temp_filename) == 499
    struc_written = read(temp_filename, PDB)
    @test modelnumbers(struc_written) == [1]
    @test countatoms(struc_written) == 492
    @test chainids(struc_written) == ['A', 'B']
    @test tempfactor(struc_written['B']["H_705"]["O"]) == 64.17
    writepdb(temp_filename, struc, standardselector, disorderselector)
    @test countlines(temp_filename) == 10
    struc_written = read(temp_filename, PDB)
    @test countatoms(struc_written) == 5
    @test sum(map(isdisorderedatom, collectatoms(struc_written))) == 5
    @test defaultaltlocid(struc_written['A'][167]["NH1"]) == 'A'

    # Test writing different element types
    writepdb(temp_filename, struc[1])
    @test countlines(temp_filename) == 3816
    struc_written = read(temp_filename, PDB)
    @test modelnumbers(struc_written) == [1]
    @test countatoms(struc_written) == 3804
    writepdb(temp_filename, struc['A'])
    @test countlines(temp_filename) == 1966
    struc_written = read(temp_filename, PDB)
    @test chainids(struc_written) == ['A']
    writepdb(temp_filename, struc['A'][50])
    @test countlines(temp_filename) == 9
    struc_written = read(temp_filename, PDB)
    @test chainids(struc_written) == ['A']
    @test countresidues(struc_written) == 1
    @test atomnames(struc_written['A'][50]) == [
        "N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"]
    writepdb(temp_filename, struc['A'][50]["CA"])
    @test countlines(temp_filename) == 1
    struc_written = read(temp_filename, PDB)
    @test countatoms(struc_written) == 1
    @test !isdisorderedatom(collectatoms(struc_written)[1])
    @test serial(collectatoms(struc_written)[1]) == 356
    writepdb(temp_filename, struc['A'][167]["CZ"])
    @test countlines(temp_filename) == 2
    struc_written = read(temp_filename, PDB)
    @test countatoms(struc_written) == 1
    @test isdisorderedatom(collectatoms(struc_written)[1])
    @test length(collectatoms(struc_written)[1]) == 2
    @test tempfactor(collectatoms(struc_written)[1]) == 16.77
    writepdb(temp_filename, Chain[struc['A'], struc['B']])
    @test countlines(temp_filename) == 3816
    struc_written = read(temp_filename, PDB)
    @test chainids(struc_written) == ['A', 'B']
    @test countatoms(struc_written['A']) == 1954
    @test countatoms(struc_written['B']) == 1850
    @test altlocids(struc_written['A']["H_215"]["O1G"]) == ['A', 'B']
    writepdb(temp_filename, AbstractResidue[struc['A'][51], struc['A'][50]])
    @test countlines(temp_filename) == 17
    struc_written = read(temp_filename, PDB)
    @test countresidues(struc_written) == 2
    @test map(resid, collectresidues(struc_written)) == ["50", "51"]
    @test countatoms(struc_written) == 17
    writepdb(temp_filename, AbstractAtom[struc['A'][51]["CA"], struc['A'][50]["CA"]])
    @test countlines(temp_filename) == 2
    struc_written = read(temp_filename, PDB)
    @test countatoms(struc_written) == 2
    @test !ishetero(struc_written['A'][51]["CA"])

    # Test multiple model writing
    struc = read(pdbfilepath("1SSU.pdb"), PDB)
    writepdb(temp_filename, Model[struc[10], struc[5]])
    @test countlines(temp_filename) == 1516
    struc_written = read(temp_filename, PDB)
    @test modelnumbers(struc_written) == [5, 10]
    @test modelnumber(defaultmodel(struc_written)) == 5
    @test countatoms(struc_written[5]) == 756
    @test countatoms(struc_written[10]) == 756
    @test_throws KeyError struc_written[1]

    # Test disordered residue writing
    struc = read(pdbfilepath("1EN2.pdb"), PDB)
    writepdb(temp_filename, struc)
    @test countlines(temp_filename) == 819
    struc_written = read(temp_filename, PDB)
    @test countatoms(struc_written) == 754
    @test isa(struc_written['A'][15], Residue)
    @test isa(struc_written['A'][16], DisorderedResidue)
    @test defaultresname(struc_written['A'][16]) == "ARG"
    @test isa(struc_written['A'][16]["N"], DisorderedAtom)
    @test defaultaltlocid(struc_written['A'][16]["N"]) == 'A'
    @test isa(disorderedres(struc_written['A'][16], "TRP")["N"], Atom)
    @test countatoms(struc_written['A'][16]) == 11
    @test countatoms(disorderedres(struc_written['A'][16], "TRP")) == 14
    writepdb(temp_filename, AbstractResidue[struc['A'][16], struc['A'][10]])
    @test countlines(temp_filename) == 46
    struc_written = read(temp_filename, PDB)
    @test countresidues(struc_written) == 2
    @test isa(struc_written['A'][10], DisorderedResidue)
    @test isa(struc_written['A'][16], DisorderedResidue)
    @test defaultresname(struc_written['A'][10]) == "SER"
    @test isa(disorderedres(struc_written['A'][10], "GLY")["O"], Atom)
    @test altlocid(disorderedres(struc_written['A'][10], "GLY")["O"]) == 'B'
    @test countatoms(struc_written['A'][10]) == 6
    @test countatoms(struc_written['A'][16]) == 11

    @test_throws ArgumentError writepdb(temp_filename, Atom(
        1, "11H11", ' ', [0.0, 0.0, 0.0], 1.0, 0.0, " H", "  ", res))

    # Delete temporary file
    rm(temp_filename)
end

@testset "Spatial" begin
    # Test coordarray
    res = Residue("ALA", 1, ' ', false, Chain('A'))
    at = Atom(100, "CA", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res)
    cs = coordarray(at)
    @test size(cs) == (3,1)
    @test cs[1] == 1.0
    @test cs[2] == 2.0
    @test cs[3] == 3.0

    struc_1AKE = read(pdbfilepath("1AKE.pdb"), PDB)
    cs = coordarray(struc_1AKE)
    @test size(cs) == (3,3804)
    @test cs[1,3787] == 20.135
    @test cs[2,3787] == -10.789
    @test cs[3,3787] == -1.732
    cs = coordarray(struc_1AKE['A'], calphaselector)
    @test size(cs) == (3,214)
    @test cs[1,10] == 17.487
    @test cs[2,10] == 42.426
    @test cs[3,10] == 19.756
    @test coordarray(cs) == cs

    # Test rmsd
    cs_one = [
        0.0 0.0
        1.0 0.0
        1.0 3.0
    ]
    cs_two = [
        0.0 2.0
        2.0 0.0
        1.0 3.0
    ]
    # This line gives a @simd warning when --inline=no
    @test isapprox(rmsd(cs_one, cs_two), sqrt(5/2))
    cs_one = [
        0.0 0.0 1.0
        1.0 0.0 2.0
        1.0 3.0 3.0
    ]
    cs_two = [
        0.0 2.0
        2.0 0.0
        1.0 3.0
    ]
    @test_throws ArgumentError rmsd(cs_one, cs_two)

    struc_1SSU = read(pdbfilepath("1SSU.pdb"), PDB)
    @test isapprox(rmsd(struc_1SSU[1], struc_1SSU[2], calphaselector), 4.1821925809691889)
    @test isapprox(rmsd(struc_1SSU[5], struc_1SSU[6], backboneselector), 5.2878196391279939)
    @test_throws ArgumentError rmsd(struc_1SSU[1]['A'][8], struc_1SSU[1]['A'][9])

    # Test displacements
    cs_one = [
        0.0 0.0
        1.0 0.0
        1.0 3.0
    ]
    cs_two = [
        0.0 2.0
        2.0 0.0
        1.0 4.0
    ]
    # This line gives a @simd warning when --inline=no
    @test isapprox(displacements(cs_one, cs_two), [1.0, sqrt(5)])
    cs_one = [
        0.0 0.0 1.0
        1.0 0.0 2.0
        1.0 3.0 3.0
    ]
    cs_two = [
        0.0 2.0
        2.0 0.0
        1.0 4.0
    ]
    @test_throws ArgumentError displacements(cs_one, cs_two)

    disps = displacements(struc_1SSU[5], struc_1SSU[10])
    @test isa(disps, Vector{Float64})
    @test length(disps) == 756
    @test isapprox(disps[20], sqrt(1.984766))
    disps = displacements(struc_1SSU[5], struc_1SSU[10], calphaselector)
    @test length(disps) == 51
    @test isapprox(disps[20], sqrt(0.032822))

    # Test sqdistance and distance
    at_a = Atom(100, "CA", ' ', [1.0, 2.0, 3.0], 1.0, 10.0, " C", "  ", res)
    at_b = Atom(110, "CA", ' ', [0.0, -1.0, 3.0], 1.0, 10.0, " C", "  ", res)
    @test sqdistance(at_a, at_b) == 10.0
    @test isapprox(distance(at_a, at_b), sqrt(10))

    @test isapprox(distance(struc_1AKE['A'], struc_1AKE['B']), sqrt(6.852947))
    @test isapprox(distance(struc_1AKE['A'], struc_1AKE['B'][50]), sqrt(530.645746))
    @test isapprox(distance(struc_1AKE['A'], struc_1AKE['B'][50]["CA"]), sqrt(574.699125))
    @test isapprox(distance(struc_1AKE['A'], struc_1AKE['B'], backboneselector), sqrt(17.350083))
    @test isapprox(distance(struc_1AKE['A'], struc_1AKE['B'], standardselector), sqrt(11.252973))
    @test isapprox(distance(struc_1AKE['A'][50]["CA"], struc_1AKE['B'][50]["CA"]), sqrt(2607.154834))

    # Test bondangle
    at_a = Atom(100, "CA", ' ', [1.0, 0.0, 1.0], 1.0, 10.0, " C", "  ", res)
    at_b = Atom(100, "CA", ' ', [0.0, 0.0, 0.0], 1.0, 10.0, " C", "  ", res)
    at_c = Atom(100, "CA", ' ', [3.0, 2.0, 1.0], 1.0, 10.0, " C", "  ", res)
    @test isapprox(bondangle(at_a, at_b, at_c), 0.713724378944765)
    vec_a = [2.0, 0.0, 0.0]
    vec_b = [2.0, 1.0, 1.0]
    @test isapprox(bondangle(vec_a, vec_b), 0.615479708670387)

    # Test dihedral functions
    at_a = Atom(100, "CA", ' ', [-1.0, -1.0, 0.0], 1.0, 10.0, " C", "  ", res)
    at_b = Atom(100, "CA", ' ', [0.0, 0.0, 0.0], 1.0, 10.0, " C", "  ", res)
    at_c = Atom(100, "CA", ' ', [1.0, 0.0, 0.0], 1.0, 10.0, " C", "  ", res)
    at_d = Atom(100, "CA", ' ', [2.0, 1.0, -1.0], 1.0, 10.0, " C", "  ", res)
    @test isapprox(dihedralangle(at_a, at_b, at_c, at_d), 2.356194490192345)
    vec_a = [1.0, 1.0, 0.0]
    vec_b = [1.0, 0.0, 0.0]
    vec_c = [1.0, -1.0, 1.0]
    @test isapprox(dihedralangle(vec_a, vec_b, vec_c), -0.785398163397448)

    @test isapprox(omegaangle(struc_1AKE['A'][20], struc_1AKE['A'][19]), -3.091913621551854, atol=1e-5)
    @test isapprox(phiangle(struc_1AKE['A'][7], struc_1AKE['A'][6]), 2.851151641716221, atol=1e-5)
    @test isapprox(psiangle(struc_1AKE['A'][8], struc_1AKE['A'][9]), 2.838265381719911, atol=1e-5)
    @test_throws ArgumentError omegaangle(struc_1AKE['A'][20], Residue("ALA", 19, ' ', false, Chain('A')))
    @test_throws ArgumentError phiangle(struc_1AKE['A'][7], Residue("ALA", 6, ' ', false, Chain('A')))
    @test_throws ArgumentError psiangle(struc_1AKE['A'][8], Residue("ALA", 9, ' ', false, Chain('A')))

    phis, psis = ramachandranangles(struc_1AKE['A'])
    @test size(phis) == (456,)
    @test size(psis) == (456,)
    @test isapprox(phis[5], -1.764512005880236, atol=1e-5)
    @test isapprox(psis[10], 0.4425094841355222, atol=1e-5)
    @test isnan(phis[1])
    @test isnan(psis[214])
    @test sum(map(x -> Int(isnan(x)), phis)) == 243
    @test sum(map(x -> Int(isnan(x)), psis)) == 243
    @test_throws ArgumentError ramachandranangles(struc_1AKE['A'][10]["CA"])

    # Test contactmap
    cas = collectatoms(struc_1AKE, calphaselector)[1:10]
    @test isa(contactmap(cas, 10), BitArray{2})
    @test contactmap(cas, 10) == [
        true  true  true  false false false false false false false
        true  true  true  true  false false false false false false
        true  true  true  true  true  true  false false false false
        false true  true  true  true  true  false false false false
        false false true  true  true  true  true  false false false
        false false true  true  true  true  true  true  true  false
        false false false false true  true  true  true  true  true
        false false false false false true  true  true  true  true
        false false false false false true  true  true  true  true
        false false false false false false true  true  true  true
    ]
    @test contactmap(struc_1AKE[1], 1.0) == [
        true  false
        false true
    ]
    cmap = contactmap(struc_1AKE['A'], 5.0)
    @test size(cmap) == (456, 456)
    @test cmap[196, 110]
    @test !cmap[15, 89]

    @test contactmap(struc_1AKE['A'][10], struc_1AKE['A'][11], 4.0) == [
        true  false false false false
        true  true  false false false
        true  true  true  false true
        true  true  true  false false
    ]
end

end # TestStructure


text:
img = open("lena.pgm");
magic = chomp(readline(img));
params = chomp(readline(img));

pm = split(params);
wd = int(pm[1]);
ht = int(pm[2]);

data = uint8(readbytes(img,wd*ht));
data = reshape(data,wd,ht);
close(img);

dout = copy(data);
for i = 2:wd-1
  for j = 2:ht-1
    temp = data[i-1:i+1, j-1:j+1];
    x = sum(Gx.*temp)
    y = sum(Gy.*temp)
    p = int(sqrt(x*x + y*y))
    dout[i,j] = (p < 256) ? uint8(p) : 0xff
  end
end
out = open("lenaX.pgm","w");
println(out,magic);
println(out,params);
write(out,dout);
close(out);

ext: jl

max_stars_repo_path: Module 3/Chapter06/lena.jl

max_issues_repo_path: Module 3/Chapter06/lena.jl

max_forks_repo_path: Module 3/Chapter06/lena.jl


text:
#!/usr/bin/env julia6
using ReadGlobal

function makenorm!(txx::AbstractArray,txy,txz,tyy,tyz,tzz,out)
    @inbounds begin
        Threads.@threads for i in linearindices(out)
            p = txx[i]^2 + tyy[i]^2 + tzz[i]^2 + 2 *(txy[i]^2 + txz[i]^2 + tyz[i]^2)
            out[i] = sqrt(2*p)
        end
    end
end

function makenorm!(tx::AbstractArray,ty,tz,out)
    @inbounds begin
        Threads.@threads for i in linearindices(out)
            p = tx[i]^2 + ty[i]^2 + tz[i]^2
            out[i] = sqrt(p)
        end
    end
end

function makenorm(tx::String,ty,tz,out)
    doinchunks(makenorm!,input=(tx,ty,tz),output=(out,))
end

function makenorm(txx::String,txy,txz,tyy,tyz,tzz,out)
    doinchunks(makenorm!,input=(txx,txy,txz,tyy,tyz,tzz),output=(out,))
end

makenorm(ARGS...)


text:
"""
The `Lexer` sub-module is concerned with transcribing a given
stream of characters into a sequence of domain specific lexical
units called "token".

Basic methodology:

1. Wrap a plain `IO` object into a [`Lexer.CharStream`](@ref).
2. Call [`Lexer.next_token`](@ref) to collect another
   [`Lexer.Token`](@ref) from the character stream.
3. Goto 2. unless end of file is reached.

For convenience the above process is simplified by providing the
type [`Lexer.TokenStream`](@ref), which supports `eof`, `read`
and [`Lexer.peek`](@ref).
"""
module Lexer

using Base.UTF8proc.isspace

"""
    LexicalError(msg)

The string or stream passed to [`Lexer.next_token`](@ref) was not
a valid sequence of characters according to the smart game
format.
"""
struct LexicalError <: Exception
    msg::String
end

"""
    Token(name::Char, [value::String])

A SGF specific lexical token. It can be either for the following:

- `Token('\\0')`: Empty token to denote trailing whitespaces.
- `Token(';')`: Separator for nodes.
- `Token('(')` and `Token(')')`: Delimiter for game trees.
- `Token('[')` and `Token(']')`: Delimiter for property values.
- `Token('I', "AB1")`: Identifier for properties. In general these
  are made up of one or more uppercase letters. However, with the
  exception of the first position, digits are also allowed to
  occur. This is done in order to supported older FF versions.
- `Token('S', "abc 23(\\)")`: Any property value between `'['`
  and `']'`. This includes moves, numbers, simple text, and text.
"""
struct Token
    name::Char
    value::String
end

Token(name::Char) = Token(name, "")

Base.hash(t::Token, h::UInt) = hash(t.name, hash(t.value, hash(:Token, h)))
Base.:(==)(a::Token, b::Token) = isequal(a.name, b.name) && isequal(a.value, b.value)

"""
    CharStream(io::IO)

Stateful decorator around `io` to keep track of some context
information, as well as allow the use of [`peek`](@ref) (i.e.
looking at the next character without consuming it).
"""
mutable struct CharStream{I<:IO}
    io::I
    inproperty::Bool
    next::Char
end

CharStream(io::IO) = CharStream(io, false, '\0')

Base.eof(cs::CharStream) = eof(cs.io) && (cs.next == '\0')

"""
    peek(cs::CharStream, ::Type{Char}) -> Char

Return the next `Char` in `cs` without consuming it, which means
that the next time `peek` or `read` is called, the same `Char`
will be returned.
"""
function peek(cs::CharStream, ::Type{Char})
    if cs.next == '\0'
        cs.next = read(cs.io, Char)::Char
    end
    cs.next
end

function Base.read(cs::CharStream, ::Type{Char})
    if cs.next == '\0'
        read(cs.io, Char)::Char
    else
        c = cs.next
        cs.next = '\0'
        c
    end
end

"""
    next_token(cs::CharStream) -> Token

Reads and returns the next [`Token`](@ref) from the given
character stream `cs`. If no more token are available, then a
`EOFError` will be thrown.

Note that the lexer should support FF[1]-FF[4] versions. In case
any unambiguously illegal character sequence is encountered, the
function will throw a [`LexicalError`](@ref).
"""
function next_token(cs::CharStream)
    eof(cs) && throw(EOFError())
    c = read(cs, Char)
    if cs.inproperty && c === ']'
        cs.inproperty = false
        Token(']')
    elseif cs.inproperty
        buf = UInt8[c]
        while true
            c = peek(cs, Char)
            if c === ']'
                break
            elseif c === '\\'
                read(cs, Char) # consume \\
                c = read(cs, Char)
                # newlines are removed after '\\'
                if c === '\r'
                    # consume additional \n if it exists
                    peek(cs, Char) == '\n' && read(cs, Char)
                elseif c === '\n'
                    # consume additional \r if it exists
                    peek(cs, Char) == '\r' && read(cs, Char)
                elseif isspace(c)
                    # all other whitespaces are converted to " "
                    push!(buf, ' ')
                else
                    # everything else is stored verbatim
                    # this includes ']', ')', and ':'
                    push!(buf, c)
                end
            elseif c === '\r'
                read(cs, Char) # consume \r
                # we represent newlines always as a single \n
                push!(buf, '\n')
                # consume additional \n if it exists
                peek(cs, Char) == '\n' && read(cs, Char)
            elseif c === '\n'
                read(cs, Char) # consume \n
                push!(buf, '\n')
                # consume additional \r if it exists
                peek(cs, Char) == '\r' && read(cs, Char)
            elseif isspace(c)
                read(cs, Char) # consume " "
                # since newlines are already processed at this point,
                # every other whitespace can be replaced with ' '
                push!(buf, ' ')
            else
                read(cs, Char) # consume character
                push!(buf, c)
            end
        end
        Token('S', String(buf))
    elseif c === '['
        cs.inproperty = true
        Token('[')
    elseif c === '('
        Token('(')
    elseif c === ')'
        Token(')')
    elseif c === ';'
        Token(';')
    elseif c in 'A':'Z'
        buf = UInt8[c]
        while true
            c = read(cs, Char)
            if c in 'A':'Z'
                push!(buf, c)
            elseif c in '0':'9' # compatibility with FF[1-3]
                push!(buf, c)
            else
                cs.next = c
                break
            end
        end
        Token('I', String(buf))
    elseif isspace(c)
        eof(cs) ? Token('\0') : next_token(cs)
    else
        throw(LexicalError("Invalid character: \"$c\""))
    end
end

"""
    TokenStream(cs::CharStream)

Stateful decorator around `cs` to allow the use of [`peek`](@ref)
(i.e. looking at the next [`Token`](@ref) without consuming it).

It uses the function [`next_token`](@ref) to create a new
[`Token`](@ref) from the current position of `cs` onwards.
"""
mutable struct TokenStream{I<:IO}
    stream::CharStream{I}
    next::Token
end

TokenStream(cs::CharStream) = TokenStream(cs, Token('\0'))
TokenStream(io::IO) = TokenStream(CharStream(io))

Base.eof(ts::TokenStream) = eof(ts.stream) && (ts.next == Token('\0'))

"""
    peek(ts::TokenStream, ::Type{Token}) -> Token

Return the next [`Token`](@ref) in `ts` without consuming it,
which means that the next time `peek` or `read` is called, the
same [`Token`](@ref) will be returned.
"""
function peek(ts::TokenStream, ::Type{Token})
    if ts.next == Token('\0')
        ts.next = next_token(ts.stream)
    end
    ts.next
end

function Base.read(ts::TokenStream, ::Type{Token})
    if ts.next == Token('\0')
        next_token(ts.stream)
    else
        t = ts.next
        ts.next = Token('\0')
        t
    end
end

end


text:
function TriangleFilter(nr, n12, uu, vv)

    pp = zeros(Float32,n12+nr-1)
    qq = zeros(Float32,n12+nr+nr-2)
    tt = zeros(Float32,n12)
    for i = 1 : n12
	qq[i] = uu[i]
    end
    rectangle_filter(nr,n12,qq,pp)
    rectangle_filter(nr,nr+n12-1,pp,qq)
    for i = 1 : n12
	tt[i] = qq[i+nr-1]
    end
    for i = 1 : nr-1
	tt[i] = tt[i] + qq[nr-i]
    end
    for i = 1 : nr-1
	tt[n12-i+1] = tt[n12-i+1] + qq[n12+(nr-1)+i]
    end
    for i = 1 : n12
	vv[i] = tt[i]
    end

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# This file illustrates how to create inferrable tree-iteration methods in circumstances
# where the children are not naturally indexable.
# See "binarytree_easy.jl" for a simpler approach.

if !isdefined(@__MODULE__, :BinaryNode)
    include("binarytree_core.jl")
end

## Enhancement of the "native" binary tree
# You might have the methods below even if you weren't trying to support AbstractTrees.

# Implement iteration over the immediate children of a node
function Base.iterate(node::BinaryNode)
    isdefined(node, :left) && return (node.left, false)
    isdefined(node, :right) && return (node.right, true)
    return nothing
end
function Base.iterate(node::BinaryNode, state::Bool)
    state && return nothing
    isdefined(node, :right) && return (node.right, true)
    return nothing
end
Base.IteratorSize(::Type{BinaryNode{T}}) where T = Base.SizeUnknown()
Base.eltype(::Type{BinaryNode{T}}) where T = BinaryNode{T}

## Things we need to define to leverage the native iterator over children
## for the purposes of AbstractTrees.
# Set the traits of this kind of tree
Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()
AbstractTrees.parentlinks(::Type{BinaryNode{T}}) where T = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{BinaryNode{T}}) where T = AbstractTrees.StoredSiblings()
# Use the native iteration for the children
AbstractTrees.children(node::BinaryNode) = node

Base.parent(root::BinaryNode, node::BinaryNode) = isdefined(node, :parent) ? node.parent : nothing

function AbstractTrees.nextsibling(tree::BinaryNode, child::BinaryNode)
    isdefined(child, :parent) || return nothing
    p = child.parent
    if isdefined(p, :right)
        child === p.right && return nothing
        return p.right
    end
    return nothing
end

# We also need `pairs` to return something sensible.
# If you don't like integer keys, you could do, e.g.,
#   Base.pairs(node::BinaryNode) = BinaryNodePairs(node)
# and have its iteration return, e.g., `:left=>node.left` and `:right=>node.right` when defined.
# But the following is easy:
Base.pairs(node::BinaryNode) = enumerate(node)

AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)

root = BinaryNode(0)
l = leftchild(1, root)
r = rightchild(2, root)
lr = rightchild(3, l)

print_tree(root)
collect(PostOrderDFS(root))
@static if isdefined(@__MODULE__, :Test)
    @testset "binarytree_infer.jl" begin
        @test @inferred(map(x->x.data, PostOrderDFS(root))) == [3, 1, 2, 0]
        @test @inferred(map(x->x.data, PreOrderDFS(root))) == [0, 1, 3, 2]
        @test @inferred(map(x->x.data, Leaves(root))) == [3, 2]
    end
end


text:
# implements a multi-agent policy assignment strategy
module PolicyAssignment

using DiscreteMDPs
using JuMP
using Gurobi

export AssignmentPolicy
export action
export merge

import DiscreteMDPs.Policy
import DiscreteMDPs.action

typealias SingleAgentPolicy Vector{Int64}
typealias SingleAgentValues Vector{Float64}

# this is an online policy
# TODO: only works for a single sniper and multiple target/monitor startegy
type AssignmentPolicy <: Policy
    nAgents::Int64
    policies::Vector{SingleAgentPolicy}
    values::Vector{SingleAgentValues}
    mdp::DiscreteMDP
end

function action(policy::AssignmentPolicy, sIdx::Int64)
    mdp      = policy.mdp
    nAgents  = policy.nAgents
    nActions = mdp.nFactoredActions # need this to get the final action
    p        = policy.policies
    v        = policy.values
    ds       = mdp.map.xSize

    # c[1,:] the values for all policies for a given agent
    c = zeros(nAgents, nAgents)
    saIdxs = zeros(Int64, nAgents) # single agent state indices
    saActions = zeros(Int64, nAgents) # single agent actions
    aSub = zeros(Int64, nAgents) + nActions # for combinig single agent actions
    dimSub = zeros(Int64, nAgents) + ds # only works for single sniper problems

    # get positions from sIdx
    i2x!(mdp, sIdx)
    x = mdp.tempSub1
    sp = [x[1], x[2]]
    # combine positions of sniper and each monitor to get single agent positions
    ai = 3 # agent counter assumes single sniper
    for i = 1:nAgents # assumes single sniper
        mp = [x[ai], x[ai+1]]
        cp = [sp, mp] # combined position
        # convert positions to idx
        sa = sub2ind(dimSub, cp) # state of single agent
        saIdxs[i] = sa
        for j = 1:nAgents
            # fill the value matrix
            c[i,j] = v[j][sa]
        end
        ai += 2
    end

    y = assignment(policy, c)
    for i = 1:nAgents
        pp = indmax(y[i,:]) # policy to be assigned to an agent
        sa = saIdxs[i]
        saActions[i] = p[pp][sa]
    end

    return sub2ind(aSub, saActions)
end

# return individual actions for each agent
function assignment(policy::AssignmentPolicy, c::Matrix)
    # c is value-policy matrix
    n = policy.nAgents
    m = Model(solver=GurobiSolver())
    # define the assignment matrix
    @defVar(m, 0 <= y[1:n,1:n] <= 1, Int)
    # the objective is double sum over the c and y values
    @setObjective(m, Max, sum{c[i,j]*y[i,j], i = 1:n, j = 1:n})
    # column and row constraints
    @addConstraints m begin
        col[i=1:n], sum(y[i,:]) == 1
        row[j=1:n], sum(y[:,j]) == 1
    end
    status = solve(m)
    # assign a policy to each agent based on y
    return getValue(y)
end

function merge{T}(a::Vector{T}...)
    l = length(a)
    n = [a[i] for i = 1:l]
    return n
end

end # module


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("Utils.jl")
include("CovFuncs.jl")
using LinearAlgebra
using RecursiveArrayTools
using Distributions
using Distances
using Random
using Match
using Profile
using JLD
using MAT
using BlockArrays
Random.seed!(123)

NTrain = 10000
NTest = 1
N = NTest + NTrain
d =  2

xTrain = rand(d, NTrain)
xTestClose = rand(d, NTest)
xTestFar = 1.5 * ones(d,NTest) .+ 0.1 * rand(d, NTest)
xTest = xTestFar
x = hcat(xTrain, xTest)

l = 0.5
σ = 1e-16 * Diagonal(2.0 * ones(N))
maternOrder = "12"

covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)

# Make sure to avoid NANs
@match maternOrder begin
  "12" => (covfunc! = covfunc12!)
  "32" => (covfunc! = covfunc32!)
  "52" => (covfunc! = covfunc52!)
end

ρ = 10.0
λ = 1.5

# =====================================
# Screw that, try different dataset:
# =====================================
h = 1 / NTrain
xTrain = h : h : 1
xTrain = vcat(sin.(2 * π * xTrain)',
              2 * cos.(2 * π * xTrain)')
xTest = 0.1 * randn(2,NTest)
x = hcat(xTrain, xTest)

# =====================================
# Prediction variables last
# =====================================

# reordering points and forming skeletons
colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
xOrd = x[:, P]
skeletons = construct_skeletons(colptr, rowval, distances, λ)

# Creating
Θ = (pairwise(Euclidean(), xOrd, dims=2))
covfunc!(Θ)
Θ = Θ

L = assembleL(skeletons, xOrd, covfunc!)

ΘPredLast = inv(Matrix(L * L'))

@show KL(Θ, ΘPredLast)
@show KL(ΘPredLast, Θ)

# =====================================
# Computing posterior mean for prediction variables first
# =====================================

colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
xTrainOrd = xTrain[:, P]
xOrd = hcat(xTrainOrd, xTest)

L = SparseMatrixCSC(NTrain, NTrain, colptr, rowval, ones(size(rowval, 1)))
L = vcat(L, ones(NTest, NTrain))
L = tril(hcat(L, ones(N, NTest)))

colptr = L.colptr
rowval = L.rowval
distances = vcat(distances, fill(Inf, NTest))

skeletons = construct_skeletons(colptr, rowval, distances, λ)

# Creating
Θ = (pairwise(Euclidean(), xOrd, dims=2))
covfunc!(Θ)
Θ = Θ

L = assembleL(skeletons, xOrd, covfunc!)

ΘPredFirst = inv(Matrix(L * L'))

@show KL(Θ, ΘPredFirst)
@show KL(ΘPredFirst, Θ)

# =====================================
# Don't include prediction variables at all
# =====================================

# reordering points and forming skeletons
colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
xOrdTrain = xTrain[:, P]
xOrd = hcat(xOrdTrain, xTest)
skeletons = construct_skeletons(colptr, rowval, distances, λ)
L = assembleL(skeletons, xOrdTrain, covfunc!, σ)

# Creating
Θ = (pairwise(Euclidean(), xOrd, dims=2))
covfunc!(Θ)
Θ = Θ

TMat = pairwise(Euclidean(), xOrdTrain, xTest, dims=2)
covfunc!(TMat)
CMat = pairwise(Euclidean(), xTest, dims=2)
covfunc!(CMat)

ΘNoPred = inv(Matrix(L * L'))
ΘNoPred = vcat(ΘNoPred, TMat')
ΘNoPred = hcat(ΘNoPred, vcat(TMat, CMat))

@show KL(Θ, ΘNoPred)
@show KL(ΘNoPred, Θ)


text:
module V3

using LinearAlgebra
import Base.:(==)
import Base.+
import Base.-
import Base.*
import Base./
import Base.isapprox
import Base.iterate

struct Vec3
    comp::Vector{Float64}
end

x(a::Vec3) = a.comp[1]
y(a::Vec3) = a.comp[2]
z(a::Vec3) = a.comp[3]
Base.iterate(a::Vec3; kwargs...) = Base.iterate(a.comp; kwargs...)
==(a::Vec3, b::Vec3) = a.comp == b.comp
+(a::Vec3, b::Vec3) = Vec3(a.comp + b.comp)
-(a::Vec3, b::Vec3) = Vec3(a.comp - b.comp)
-(a::Vec3) = Vec3(-a.comp)
*(a::Vec3, b::Number) = Vec3(a.comp * b)
/(a::Vec3, b::Number) = Vec3(a.comp / b)
LinearAlgebra.dot(a::Vec3, b::Vec3) = a.comp ⋅ b.comp
LinearAlgebra.cross(a::Vec3, b::Vec3) = Vec3(a.comp × b.comp)
LinearAlgebra.norm(a::Vec3) = norm(a.comp)
isapprox(a::Vec3, b::Vec3; kwargs...) = isapprox(a.comp, b.comp; kwargs...)
LinearAlgebra.normalize(a::Vec3) = Vec3(normalize(a.comp))
at_len(a::Vec3, b::Number) = normalize(a) * b

end # module V3

max_stars_repo_name: hdtx/raytracer-comparison

max_issues_repo_name: hdtx/raytracer-comparison

max_forks_repo_name: hdtx/raytracer-comparison


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const GDAL_VERSION_MAJOR = 1
const GDAL_VERSION_MINOR = 11
const GDAL_VERSION_REV = 2
const GDAL_VERSION_BUILD = 0

# Skipping MacroDefinition: GDAL_COMPUTE_VERSION ( maj , min , rev ) ( ( maj ) * 1000000 + ( min ) * 10000 + ( rev ) * 100 )
# Skipping MacroDefinition: GDAL_VERSION_NUM ( GDAL_COMPUTE_VERSION ( GDAL_VERSION_MAJOR , GDAL_VERSION_MINOR , GDAL_VERSION_REV ) + GDAL_VERSION_BUILD )

const GDAL_RELEASE_DATE = 20150210
const GDAL_RELEASE_NAME = "1.11.2"
const GDALMD_AREA_OR_POINT = "AREA_OR_POINT"
const GDALMD_AOP_AREA = "Area"
const GDALMD_AOP_POINT = "Point"
const CPLE_WrongFormat = 200
const GDAL_DMD_LONGNAME = "DMD_LONGNAME"
const GDAL_DMD_HELPTOPIC = "DMD_HELPTOPIC"
const GDAL_DMD_MIMETYPE = "DMD_MIMETYPE"
const GDAL_DMD_EXTENSION = "DMD_EXTENSION"
const GDAL_DMD_CREATIONOPTIONLIST = "DMD_CREATIONOPTIONLIST"
const GDAL_DMD_CREATIONDATATYPES = "DMD_CREATIONDATATYPES"
const GDAL_DMD_SUBDATASETS = "DMD_SUBDATASETS"
const GDAL_DCAP_CREATE = "DCAP_CREATE"
const GDAL_DCAP_CREATECOPY = "DCAP_CREATECOPY"
const GDAL_DCAP_VIRTUALIO = "DCAP_VIRTUALIO"

# Skipping MacroDefinition: SRCVAL ( papoSource , eSrcType , ii ) ( eSrcType == GDT_Byte ? ( ( GByte * ) papoSource ) [ ii ] : ( eSrcType == GDT_Float32 ? ( ( float * ) papoSource ) [ ii ] : ( eSrcType == GDT_Float64 ? ( ( double * ) papoSource ) [ ii ] : ( eSrcType == GDT_Int32 ? ( ( GInt32 * ) papoSource ) [ ii ] : ( eSrcType == GDT_UInt16 ? ( ( GUInt16 * ) papoSource ) [ ii ] : ( eSrcType == GDT_Int16 ? ( ( GInt16 * ) papoSource ) [ ii ] : ( eSrcType == GDT_UInt32 ? ( ( GUInt32 * ) papoSource ) [ ii ] : ( eSrcType == GDT_CInt16 ? ( ( GInt16 * ) papoSource ) [ ii * 2 ] : ( eSrcType == GDT_CInt32 ? ( ( GInt32 * ) papoSource ) [ ii * 2 ] : ( eSrcType == GDT_CFloat32 ? ( ( float * ) papoSource ) [ ii * 2 ] : ( eSrcType == GDT_CFloat64 ? ( ( double * ) papoSource ) [ ii * 2 ] : 0 ) ) ) ) ) ) ) ) ) ) )

const GMF_ALL_VALID = 0x01
const GMF_PER_DATASET = 0x02
const GMF_ALPHA = 0x04
const GMF_NODATA = 0x08

# Skipping MacroDefinition: GDAL_CHECK_VERSION ( pszCallingComponentName ) GDALCheckVersion ( GDAL_VERSION_MAJOR , GDAL_VERSION_MINOR , pszCallingComponentName )

# begin enum GDALDataType
typealias GDALDataType UInt32
const GDT_Unknown = (UInt32)(0)
const GDT_Byte = (UInt32)(1)
const GDT_UInt16 = (UInt32)(2)
const GDT_Int16 = (UInt32)(3)
const GDT_UInt32 = (UInt32)(4)
const GDT_Int32 = (UInt32)(5)
const GDT_Float32 = (UInt32)(6)
const GDT_Float64 = (UInt32)(7)
const GDT_CInt16 = (UInt32)(8)
const GDT_CInt32 = (UInt32)(9)
const GDT_CFloat32 = (UInt32)(10)
const GDT_CFloat64 = (UInt32)(11)
const GDT_TypeCount = (UInt32)(12)
# end enum GDALDataType

# begin enum GDALAsyncStatusType
typealias GDALAsyncStatusType UInt32
const GARIO_PENDING = (UInt32)(0)
const GARIO_UPDATE = (UInt32)(1)
const GARIO_ERROR = (UInt32)(2)
const GARIO_COMPLETE = (UInt32)(3)
const GARIO_TypeCount = (UInt32)(4)
# end enum GDALAsyncStatusType

# begin enum GDALAccess
typealias GDALAccess UInt32
const GA_ReadOnly = (UInt32)(0)
const GA_Update = (UInt32)(1)
# end enum GDALAccess

# begin enum GDALRWFlag
typealias GDALRWFlag UInt32
const GF_Read = (UInt32)(0)
const GF_Write = (UInt32)(1)
# end enum GDALRWFlag

# begin enum GDALColorInterp
typealias GDALColorInterp UInt32
const GCI_Undefined = (UInt32)(0)
const GCI_GrayIndex = (UInt32)(1)
const GCI_PaletteIndex = (UInt32)(2)
const GCI_RedBand = (UInt32)(3)
const GCI_GreenBand = (UInt32)(4)
const GCI_BlueBand = (UInt32)(5)
const GCI_AlphaBand = (UInt32)(6)
const GCI_HueBand = (UInt32)(7)
const GCI_SaturationBand = (UInt32)(8)
const GCI_LightnessBand = (UInt32)(9)
const GCI_CyanBand = (UInt32)(10)
const GCI_MagentaBand = (UInt32)(11)
const GCI_YellowBand = (UInt32)(12)
const GCI_BlackBand = (UInt32)(13)
const GCI_YCbCr_YBand = (UInt32)(14)
const GCI_YCbCr_CbBand = (UInt32)(15)
const GCI_YCbCr_CrBand = (UInt32)(16)
const GCI_Max = (UInt32)(16)
# end enum GDALColorInterp

# begin enum GDALPaletteInterp
typealias GDALPaletteInterp UInt32
const GPI_Gray = (UInt32)(0)
const GPI_RGB = (UInt32)(1)
const GPI_CMYK = (UInt32)(2)
const GPI_HLS = (UInt32)(3)
# end enum GDALPaletteInterp

typealias GDALMajorObjectH Ptr{Void}
typealias GDALDatasetH Ptr{Void}
typealias GDALRasterBandH Ptr{Void}
typealias GDALDriverH Ptr{Void}
typealias GDALProjDefH Ptr{Void}
typealias GDALColorTableH Ptr{Void}
typealias GDALRasterAttributeTableH Ptr{Void}
typealias GDALAsyncReaderH Ptr{Void}

type GDALOptionDefinition
    pszOptionName::Ptr{UInt8}
    pszValueType::Ptr{UInt8}
    pszDescription::Ptr{UInt8}
    papszOptions::Ptr{Ptr{UInt8}}
end

type GDAL_GCP
    pszId::Ptr{UInt8}
    pszInfo::Ptr{UInt8}
    dfGCPPixel::Cdouble
    dfGCPLine::Cdouble
    dfGCPX::Cdouble
    dfGCPY::Cdouble
    dfGCPZ::Cdouble
end

typealias GDALDerivedPixelFunc Ptr{Void}

immutable Array_20_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
    d7::Cdouble
    d8::Cdouble
    d9::Cdouble
    d10::Cdouble
    d11::Cdouble
    d12::Cdouble
    d13::Cdouble
    d14::Cdouble
    d15::Cdouble
    d16::Cdouble
    d17::Cdouble
    d18::Cdouble
    d19::Cdouble
    d20::Cdouble
end

zero(::Type{Array_20_Cdouble}) = begin  # /Users/fgans/.julia/v0.4/Clang/src/wrap_c.jl, line 264:
        Array_20_Cdouble(fill(zero(Cdouble),20)...)
    end

type GDALRPCInfo
    dfLINE_OFF::Cdouble
    dfSAMP_OFF::Cdouble
    dfLAT_OFF::Cdouble
    dfLONG_OFF::Cdouble
    dfHEIGHT_OFF::Cdouble
    dfLINE_SCALE::Cdouble
    dfSAMP_SCALE::Cdouble
    dfLAT_SCALE::Cdouble
    dfLONG_SCALE::Cdouble
    dfHEIGHT_SCALE::Cdouble
    adfLINE_NUM_COEFF::Array_20_Cdouble
    adfLINE_DEN_COEFF::Array_20_Cdouble
    adfSAMP_NUM_COEFF::Array_20_Cdouble
    adfSAMP_DEN_COEFF::Array_20_Cdouble
    dfMIN_LONG::Cdouble
    dfMIN_LAT::Cdouble
    dfMAX_LONG::Cdouble
    dfMAX_LAT::Cdouble
end

type GDALColorEntry
    c1::Int16
    c2::Int16
    c3::Int16
    c4::Int16
end

# begin enum GDALRATFieldType
typealias GDALRATFieldType UInt32
const GFT_Integer = (UInt32)(0)
const GFT_Real = (UInt32)(1)
const GFT_String = (UInt32)(2)
# end enum GDALRATFieldType

# begin enum GDALRATFieldUsage
typealias GDALRATFieldUsage UInt32
const GFU_Generic = (UInt32)(0)
const GFU_PixelCount = (UInt32)(1)
const GFU_Name = (UInt32)(2)
const GFU_Min = (UInt32)(3)
const GFU_Max = (UInt32)(4)
const GFU_MinMax = (UInt32)(5)
const GFU_Red = (UInt32)(6)
const GFU_Green = (UInt32)(7)
const GFU_Blue = (UInt32)(8)
const GFU_Alpha = (UInt32)(9)
const GFU_RedMin = (UInt32)(10)
const GFU_GreenMin = (UInt32)(11)
const GFU_BlueMin = (UInt32)(12)
const GFU_AlphaMin = (UInt32)(13)
const GFU_RedMax = (UInt32)(14)
const GFU_GreenMax = (UInt32)(15)
const GFU_BlueMax = (UInt32)(16)
const GFU_AlphaMax = (UInt32)(17)
const GFU_MaxCount = (UInt32)(18)
# end enum GDALRATFieldUsage

# begin enum GDALTileOrganization
typealias GDALTileOrganization UInt32
const GTO_TIP = (UInt32)(0)
const GTO_BIT = (UInt32)(1)
const GTO_BSQ = (UInt32)(2)
# end enum GDALTileOrganization

typealias GDALTransformerFunc Ptr{Void}

immutable Array_4_UInt8
    d1::UInt8
    d2::UInt8
    d3::UInt8
    d4::UInt8
end

zero(::Type{Array_4_UInt8}) = begin  # /Users/fgans/.julia/v0.4/Clang/src/wrap_c.jl, line 264:
        Array_4_UInt8(fill(zero(UInt8),4)...)
    end

type GDALTransformerInfo
    szSignature::Array_4_UInt8
    pszClassName::Ptr{UInt8}
    pfnTransform::GDALTransformerFunc
    pfnCleanup::Ptr{Void}
    pfnSerialize::Ptr{Void}
end

typealias GDALContourWriter Ptr{Void}
typealias GDALContourGeneratorH Ptr{Void}

immutable Array_6_Cdouble
    d1::Cdouble
    d2::Cdouble
    d3::Cdouble
    d4::Cdouble
    d5::Cdouble
    d6::Cdouble
end

zero(::Type{Array_6_Cdouble}) = begin  # /Users/fgans/.julia/v0.4/Clang/src/wrap_c.jl, line 264:
        Array_6_Cdouble(fill(zero(Cdouble),6)...)
    end

type OGRContourWriterInfo
    hLayer::Ptr{Void}
    adfGeoTransform::Array_6_Cdouble
    nElevField::Cint
    nIDField::Cint
    nNextID::Cint
end

# begin enum GDALGridAlgorithm
typealias GDALGridAlgorithm UInt32
const GGA_InverseDistanceToAPower = (UInt32)(1)
const GGA_MovingAverage = (UInt32)(2)
const GGA_NearestNeighbor = (UInt32)(3)
const GGA_MetricMinimum = (UInt32)(4)
const GGA_MetricMaximum = (UInt32)(5)
const GGA_MetricRange = (UInt32)(6)
const GGA_MetricCount = (UInt32)(7)
const GGA_MetricAverageDistance = (UInt32)(8)
const GGA_MetricAverageDistancePts = (UInt32)(9)
# end enum GDALGridAlgorithm

type GDALGridInverseDistanceToAPowerOptions
    dfPower::Cdouble
    dfSmoothing::Cdouble
    dfAnisotropyRatio::Cdouble
    dfAnisotropyAngle::Cdouble
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMaxPoints::UInt32
    nMinPoints::UInt32
    dfNoDataValue::Cdouble
end

type GDALGridMovingAverageOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMinPoints::UInt32
    dfNoDataValue::Cdouble
end

type GDALGridNearestNeighborOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    dfNoDataValue::Cdouble
end

type GDALGridDataMetricsOptions
    dfRadius1::Cdouble
    dfRadius2::Cdouble
    dfAngle::Cdouble
    nMinPoints::UInt32
    dfNoDataValue::Cdouble
end

const GMO_VALID = 0x0001
const GMO_IGNORE_UNIMPLEMENTED = 0x0002
const GMO_SUPPORT_MD = 0x0004
const GMO_SUPPORT_MDMD = 0x0008
const GMO_MD_DIRTY = 0x0010
const GMO_PAM_CLASS = 0x0020

# Skipping MacroDefinition: ARE_REAL_EQUAL ( dfVal1 , dfVal2 ) ( dfVal1 == dfVal2 || fabs ( dfVal1 - dfVal2 ) < 1e-10 || ( dfVal2 != 0 && fabs ( 1 - dfVal1 / dfVal2 ) < 1e-10 ) )
# Skipping MacroDefinition: DIV_ROUND_UP ( a , b ) ( ( ( a ) % ( b ) ) == 0 ? ( ( a ) / ( b ) ) : ( ( ( a ) / ( b ) ) + 1 ) )

const GDALSTAT_APPROX_NUMSAMPLES = 2500
const GCIF_GEOTRANSFORM = 0x01
const GCIF_PROJECTION = 0x02
const GCIF_METADATA = 0x04
const GCIF_GCPS = 0x08
const GCIF_NODATA = 0x00001000
const GCIF_CATEGORYNAMES = 0x00002000
const GCIF_MINMAX = 0x00004000
const GCIF_SCALEOFFSET = 0x00008000
const GCIF_UNITTYPE = 0x00010000
const GCIF_COLORTABLE = 0x00020000
const GCIF_COLORINTERP = 0x00020000
const GCIF_BAND_METADATA = 0x00040000
const GCIF_RAT = 0x00080000
const GCIF_MASK = 0x00100000
const GCIF_BAND_DESCRIPTION = 0x00200000
const GCIF_ONLY_IF_MISSING = 0x10000000
const GCIF_PROCESS_BANDS = 0x20000000
const GCIF_PAM_DEFAULT = (((((((((((((((GCIF_GEOTRANSFORM | GCIF_PROJECTION) | GCIF_METADATA) | GCIF_GCPS) | GCIF_NODATA) | GCIF_CATEGORYNAMES) | GCIF_MINMAX) | GCIF_SCALEOFFSET) | GCIF_UNITTYPE) | GCIF_COLORTABLE) | GCIF_COLORINTERP) | GCIF_BAND_METADATA) | GCIF_RAT) | GCIF_MASK) | GCIF_ONLY_IF_MISSING) | GCIF_PROCESS_BANDS) | GCIF_BAND_DESCRIPTION
const GPF_DIRTY = 0x01
const GPF_TRIED_READ_FAILED = 0x02
const GPF_DISABLED = 0x04
const GPF_AUXMODE = 0x08
const GPF_NOSAVE = 0x10

type GDALRasterBandPamInfo
    poParentDS::Ptr{Cint}
    bNoDataValueSet::Cint
    dfNoDataValue::Cdouble
    poColorTable::Ptr{Cint}
    eColorInterp::GDALColorInterp
    pszUnitType::Ptr{UInt8}
    papszCategoryNames::Ptr{Ptr{UInt8}}
    dfOffset::Cdouble
    dfScale::Cdouble
    bHaveMinMax::Cint
    dfMin::Cdouble
    dfMax::Cdouble
    bHaveStats::Cint
    dfMean::Cdouble
    dfStdDev::Cdouble
    psSavedHistograms::Ptr{Void} # Ptr{CPLXMLNode}
    poDefaultRAT::Ptr{Cint}
end

const MAX_ULPS = 10

# begin enum GDALBurnValueSrc
typealias GDALBurnValueSrc UInt32
const GBV_UserBurnValue = (UInt32)(0)
const GBV_Z = (UInt32)(1)
const GBV_M = (UInt32)(2)
# end enum GDALBurnValueSrc

# begin enum GDALRasterMergeAlg
typealias GDALRasterMergeAlg UInt32
const GRMA_Replace = (UInt32)(0)
const GRMA_Add = (UInt32)(1)
# end enum GDALRasterMergeAlg

type GDALRasterizeInfo
    pabyChunkBuf::Ptr{Cuchar}
    nXSize::Cint
    nYSize::Cint
    nBands::Cint
    eType::GDALDataType
    padfBurnValue::Ptr{Cdouble}
    eBurnValueSource::GDALBurnValueSrc
    eMergeAlg::GDALRasterMergeAlg
end

typealias llScanlineFunc Ptr{Void}
typealias llPointFunc Ptr{Void}
typealias GDALTransformDeserializeFunc Ptr{Void}

type GDALProxyPoolDatasetHS
end

typealias GDALProxyPoolDatasetH Ptr{GDALProxyPoolDatasetHS}

const RAT_MAX_ELEM_FOR_CLONE = 1000000

typealias VRTImageReadFunc Ptr{Void}
typealias VRTDriverH Ptr{Void}
typealias VRTSourceH Ptr{Void}
typealias VRTSimpleSourceH Ptr{Void}
typealias VRTAveragedSourceH Ptr{Void}
typealias VRTComplexSourceH Ptr{Void}
typealias VRTFilteredSourceH Ptr{Void}
typealias VRTKernelFilteredSourceH Ptr{Void}
typealias VRTAverageFilteredSourceH Ptr{Void}
typealias VRTFuncSourceH Ptr{Void}
typealias VRTDatasetH Ptr{Void}
typealias VRTWarpedDatasetH Ptr{Void}
typealias VRTRasterBandH Ptr{Void}
typealias VRTSourcedRasterBandH Ptr{Void}
typealias VRTWarpedRasterBandH Ptr{Void}
typealias VRTDerivedRasterBandH Ptr{Void}
typealias VRTRawRasterBandH Ptr{Void}

const EXIFOFFSETTAG = 0x8769
const INTEROPERABILITYOFFSET = 0xa005
const GPSOFFSETTAG = 0x8825
const MAXSTRINGLENGTH = 65535

type gpsname
    tag::Cint
    name::Ptr{UInt8}
end

type tagname
    tag::Cint
    name::Ptr{UInt8}
end

type intr_tag
    tag::Cint
    name::Ptr{UInt8}
end

# begin enum TIFFDataType
typealias TIFFDataType UInt32
const TIFF_NOTYPE = (UInt32)(0)
const TIFF_BYTE = (UInt32)(1)
const TIFF_ASCII = (UInt32)(2)
const TIFF_SHORT = (UInt32)(3)
const TIFF_LONG = (UInt32)(4)
const TIFF_RATIONAL = (UInt32)(5)
const TIFF_SBYTE = (UInt32)(6)
const TIFF_UNDEFINED = (UInt32)(7)
const TIFF_SSHORT = (UInt32)(8)
const TIFF_SLONG = (UInt32)(9)
const TIFF_SRATIONAL = (UInt32)(10)
const TIFF_FLOAT = (UInt32)(11)
const TIFF_DOUBLE = (UInt32)(12)
const TIFF_IFD = (UInt32)(13)
# end enum TIFFDataType

type TIFFDirEntry
    tdir_tag::Cint
    tdir_type::Cint
    tdir_count::Cint
    tdir_offset::Cint
end

typealias GDALGridFunction Ptr{Void}

# Skipping MacroDefinition: GCC_CPUID ( level , a , b , c , d ) __asm__ ( "xchgl %%ebx, %1\n" "cpuid\n" "xchgl %%ebx, %1" : "=a" ( a ) , "=r" ( b ) , "=c" ( c ) , "=d" ( d ) : "0" ( level ) )

type GDALGridXYArrays
    padfX::Ptr{Cdouble}
    padfY::Ptr{Cdouble}
end

type GDALGridPoint
    psXYArrays::Ptr{GDALGridXYArrays}
    i::Cint
end

type GDALGridExtraParameters
    hQuadTree::Ptr{Void} # Ptr{CPLQuadTree}
    dfInitialSearchRadius::Cdouble
    pafX::Ptr{Cfloat}
    pafY::Ptr{Cfloat}
    pafZ::Ptr{Cfloat}
end

type _GDALJP2GeoTIFFBox
end

typealias GDALJP2GeoTIFFBox Void

# begin enum GDALResampleAlg
typealias GDALResampleAlg UInt32
const GRA_NearestNeighbour = (UInt32)(0)
const GRA_Bilinear = (UInt32)(1)
const GRA_Cubic = (UInt32)(2)
const GRA_CubicSpline = (UInt32)(3)
const GRA_Lanczos = (UInt32)(4)
const GRA_Average = (UInt32)(5)
const GRA_Mode = (UInt32)(6)
# end enum GDALResampleAlg

typealias GDALMaskFunc Ptr{Void}
typealias GDALProgressFunc Ptr{Void}
typealias GDALTransformerFunc Ptr{Void}

type GDALWarpOptions
    papszWarpOptions::Ptr{Ptr{UInt8}}
    dfWarpMemoryLimit::Cdouble
    eResampleAlg::GDALResampleAlg
    eWorkingDataType::GDALDataType
    hSrcDS::GDALDatasetH
    hDstDS::GDALDatasetH
    nBandCount::Cint
    panSrcBands::Ptr{Cint}
    panDstBands::Ptr{Cint}
    nSrcAlphaBand::Cint
    nDstAlphaBand::Cint
    padfSrcNoDataReal::Ptr{Cdouble}
    padfSrcNoDataImag::Ptr{Cdouble}
    padfDstNoDataReal::Ptr{Cdouble}
    padfDstNoDataImag::Ptr{Cdouble}
    pfnProgress::GDALProgressFunc
    pProgressArg::Ptr{Void}
    pfnTransformer::GDALTransformerFunc
    pTransformerArg::Ptr{Void}
    papfnSrcPerBandValidityMaskFunc::Ptr{GDALMaskFunc}
    papSrcPerBandValidityMaskFuncArg::Ptr{Ptr{Void}}
    pfnSrcValidityMaskFunc::GDALMaskFunc
    pSrcValidityMaskFuncArg::Ptr{Void}
    pfnSrcDensityMaskFunc::GDALMaskFunc
    pSrcDensityMaskFuncArg::Ptr{Void}
    pfnDstDensityMaskFunc::GDALMaskFunc
    pDstDensityMaskFuncArg::Ptr{Void}
    pfnDstValidityMaskFunc::GDALMaskFunc
    pDstValidityMaskFuncArg::Ptr{Void}
    pfnPreWarpChunkProcessor::Ptr{Void}
    pPreWarpProcessorArg::Ptr{Void}
    pfnPostWarpChunkProcessor::Ptr{Void}
    pPostWarpProcessorArg::Ptr{Void}
    hCutline::Ptr{Void}
    dfCutlineBlendDist::Cdouble
end

typealias GDALWarpOperationH Ptr{Void}

# CPLErr Data Type: added by hand from cplerror.h
typealias CPLErr Cint
const CE_None = Cint(0)
const CE_Debug = Cint(1)
const CE_Warning = Cint(2)
const CE_Failure = Cint(3)
const CE_Fatal = Cint(4)