text:
"""
    mutable struct Structure_Wat  <: Str

This contains information for model constructing of a series of water models.

# Supported list:
 - structure_tip3p()
 - structure_spc()
 - structure_spce()
 - structure_tip4p_cut()
 - structure_tip4p_long()
 - structure_tip4p_ew()
 - structure_tip4p_fq()
 - structure_tip4p_2005()
 - structure_tip4p_ice()
 - structure_tip5p()
 - structure_tip4p_2005()
 - structure_tip5p_2018()
 - structure_tip7p()

# Parameters source:
- https://en.wikipedia.org/wiki/Water_model
- http://www1.lsbu.ac.uk/water/water_models.html
- https://lammps.sandia.gov/doc/Howto_tip3p.html
- https://lammps.sandia.gov/doc/Howto_tip4p.html
- https://lammps.sandia.gov/doc/Howto_spc.html
"""
mutable struct Structure_Wat <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_name
    atom_charge
    para_mass
    num_atoms
    num_atom_types
    bond_mode
    num_bonds
    num_bond_types
    angle_mode
    num_angles
    num_angle_types
    vec_type_id
    res_name
end

"""
    structure_tip3p()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip3p moedel.
"""
function structure_tip3p()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.25;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572;          # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    charge_vec = [0.4170, -0.8340]

    atom_type = [1, 2, 1]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994]
    atom_name = split("H O")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    atom_vec = [
        0          0          0
        1/ratio[1] 1/ratio[2] 0
        2/ratio[1] 0          0
    ]
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1
    res_name = "TIP3"

    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_spc()

This will generate a `Structure_Wat` type which contains all information needed to build a SPC moedel.
"""
function structure_spc()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 109.47;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 1.0;          # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    charge_vec = [0.41, -0.82]

    atom_type = [1, 2, 1]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994]
    atom_name = split("H O")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    atom_vec = [
        0          0          0
        1/ratio[1] 1/ratio[2] 0
        2/ratio[1] 0          0
    ]
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1
    res_name="SPC"

    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_spce()

This will generate a `Structure_Wat` type which contains all information needed to build a SPCE moedel.
"""
function structure_spce()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 109.47;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 1.0;          # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    charge_vec = [0.4238, -0.8476]

    atom_type = [1, 2, 1]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994]
    atom_name = split("H O")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    atom_vec = [
        0          0          0
        1/ratio[1] 1/ratio[2] 0
        2/ratio[1] 0          0
    ]
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1
    res_name = "SPCE"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_cut()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p moedel used in lammps with lj/cut/tip4p/cut. Detail info: https://lammps.sandia.gov/doc/Howto_tip4p.html

# Model info
- O mass = 15.9994
- H mass = 1.008
- O charge = -1.040
- H charge = 0.520
- ``r_0`` of OH bond = 0.9572
- ``θ`` of HOH angle = 104.52
- OM distance = 0.15
- LJ ``ε`` of O-O = 0.1550
- LJ `` σ`` of O-O = 3.1536
- LJ ``ε``, `` σ`` of OH, HH = 0.0
- Coulomb cutoff = 8.5
"""
function structure_tip4p_cut()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    charge_vec = [-1.040, 0.520]

    atom_type = [1, 2, 2]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [15.9994, 1.008]
    atom_name = split("O H")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        1/ratio[1] 1/ratio[2] 0
        0          0          0
        2/ratio[1] 0          0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [1, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [2, 1, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1
    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_long()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p moedel used in lammps with lj/cut/tip4p/long. Detail info: https://lammps.sandia.gov/doc/Howto_tip4p.html

# Model info
- O mass = 15.9994
- H mass = 1.008
- O charge = -1.0484
- H charge = 0.5242
- ``r_0`` of OH bond = 0.9572
- ``θ`` of HOH angle = 104.52
- OM distance = 0.1250
- LJ ``ε`` of O-O = 0.16275
- LJ `` σ`` of O-O = 3.16435
- LJ ``ε``, `` σ`` of OH, HH = 0.0
- Coulomb cutoff = 8.5
"""
function structure_tip4p_long()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    charge_vec = [-1.0484, 0.5242]

    atom_type = [1, 2, 2]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [15.9994, 1.008]
    atom_name = split("O H")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        1/ratio[1] 1/ratio[2] 0
        0          0          0
        2/ratio[1] 0          0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [1, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [2, 1, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1

    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_ew()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p_Ew moedel. Detail info: http://www1.lsbu.ac.uk/water/water_models.html
"""
function structure_tip4p_ew()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.125           # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    charge_vec = [0.52422, 0.0, -1.04844]

    atom_type = [1, 2, 1, 3]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994, 0]
    atom_name = split("H O Dummy")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        0          0          0
        1/ratio[1] 1/ratio[2] 0
        2/ratio[1] 0          0
        1/ratio[1] 1/cell_vec[2]*dummy_len 0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 3)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    bond_mode[3] = Bond(0, 1, [2, 4])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 2)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    angle_mode[2] = Angle(2, 1, [1, 2, 4])
    num_angles = length(angle_mode)
    num_angle_types = 2

    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_fq()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p_FQ moedel. Detail info: http://www1.lsbu.ac.uk/water/water_models.html
"""
function structure_tip4p_fq()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.15            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    charge_vec = [0.63, 0.0, -1.26]

    atom_type = [1, 2, 1, 3]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994, 0]
    atom_name = split("H O Dummy")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        0          0          0
        1/ratio[1] 1/ratio[2] 0
        2/ratio[1] 0          0
        1/ratio[1] 1/cell_vec[2]*dummy_len 0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 3)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    bond_mode[3] = Bond(0, 1, [2, 4])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 2)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    angle_mode[2] = Angle(2, 1, [1, 2, 4])
    num_angles = length(angle_mode)
    num_angle_types = 2

    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_2005()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p_2005 moedel used in lammps with lj/cut/tip4p/cut. Detail info: https://lammps.sandia.gov/doc/Howto_tip4p.html

# Model info
- O mass = 15.9994
- H mass = 1.008
- O charge = -1.1128
- H charge = 0.5564
- ``r_0`` of OH bond = 0.9572
- ``θ`` of HOH angle = 104.52
- OM distance = 0.1546
- LJ ``ε`` of O-O = 0.1852
- LJ `` σ`` of O-O = 3.1589
- LJ ``ε``, `` σ`` of OH, HH = 0.0
- Coulomb cutoff = 8.5
"""
function structure_tip4p_2005()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.15            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    charge_vec = [0.5564, -1.1128]

    atom_type = [1, 2, 2]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [15.9994, 1.008]
    atom_name = split("O H")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        1/ratio[1] 1/ratio[2] 0
        0          0          0
        2/ratio[1] 0          0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [1, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [2, 1, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1

    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip4p_ice()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip4p_Ice moedel used in lammps with lj/cut/tip4p/cut. Detail info: https://lammps.sandia.gov/doc/Howto_tip4p.html

# Model info
- O mass = 15.9994
- H mass = 1.008
- O charge = -1.1794
- H charge = 0.5897
- ``r_0`` of OH bond = 0.9572
- ``θ`` of HOH angle = 104.52
- OM distance = 0.1577
- LJ ``ε`` of O-O = 0.21084
- LJ `` σ`` of O-O = 3.1668
- LJ ``ε``, `` σ`` of OH, HH = 0.0
- Coulomb cutoff = 8.5"""
function structure_tip4p_ice()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.15            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    charge_vec = [0.5897, -1.1794]

    atom_type = [1, 2, 2]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [15.9994, 1.008]
    atom_name = split("O H")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        1/ratio[1] 1/ratio[2] 0
        0          0          0
        2/ratio[1] 0          0
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 2)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [1, 3])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 1)
    angle_mode[1] = Angle(1, 1, [2, 1, 3])
    num_angles = length(angle_mode)
    num_angle_types = 1

    res_name = "TIP4P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip5p()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip5p moedel. Detail info: http://www1.lsbu.ac.uk/water/water_models.html
"""
function structure_tip5p()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.70            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    dummy_ang = 109.47 / 180 * pi
    charge_vec = [0.2410, 0.0, -0.2410]

    atom_type = [1, 2, 1, 3, 3]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994, 0]
    atom_name = split("H O Dummy")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        0          0          0.5
        1/ratio[1] 1/ratio[2] 0.5
        2/ratio[1] 0          0.5
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 + 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 - 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 4)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    bond_mode[3] = Bond(0, 1, [2, 4])
    bond_mode[4] = Bond(0, 1, [2, 5])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 2)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    angle_mode[2] = Angle(2, 1, [2, 4, 5])
    num_angles = length(angle_mode)
    num_angle_types = 2

    res_name = "TIP5P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip5p_2018()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip5p_2018 moedel. Detail info: http://www1.lsbu.ac.uk/water/water_models.html
"""
function structure_tip5p_2018()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.70            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    dummy_ang = 109.47 / 180 * pi
    charge_vec = [0.394137, -0.641114, -0.07358]

    atom_type = [1, 2, 1, 3, 3]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994, 0]
    atom_name = split("H O Dummy")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        0          0          0.5
        1/ratio[1] 1/ratio[2] 0.5
        2/ratio[1] 0          0.5
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 + 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 - 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 4)
    bond_mode[1] = Bond(0, 1, [1, 2])
    bond_mode[2] = Bond(0, 1, [2, 3])
    bond_mode[3] = Bond(0, 1, [2, 4])
    bond_mode[4] = Bond(0, 1, [2, 5])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 2)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    angle_mode[2] = Angle(2, 1, [2, 4, 5])
    num_angles = length(angle_mode)
    num_angle_types = 2

    res_name = "TIP5P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end

"""
    structure_tip7p()

This will generate a `Structure_Wat` type which contains all information needed to build a Tip7p moedel. Detail info: http://www1.lsbu.ac.uk/water/water_models.html
"""
function structure_tip7p()
    # Parameters of water
    density = 1 / Const_cm2an^3;      # Unit: g/A^3
    angle = 104.52;
    angle = angle / 180 * pi;   # Convert int to rad
    ratio = [4, 5]              # Ratio to determing shape of box
    bond_len = 0.9572           # Unit: Anstrom. https://en.wikipedia.org/wiki/Water_model for more details
    dummy_len = 0.41            # Unit: Anstrom. http://www1.lsbu.ac.uk/water/water_models.html for more details
    dummy_ang = 109.47 / 180 * pi
    charge_vec = [0.58014, 0.11094, -0.17724, 0.45837]

    atom_type = [1, 2, 1, 3, 3, 4, 4]
    atom_charge = [charge_vec[atom_type[n]] for n=1:length(atom_type)]
    para_mass = [1.00784, 15.9994, 0, 0]
    atom_name = split("H O Dummy01 Dummy02")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)

    # Setting Atom Args
    cell_vec = [
        ratio[1] * bond_len * sin(angle/2)
        ratio[2] * bond_len * cos(angle/2)
        ( (para_mass[1]*2+para_mass[2]) * Const_gm2g) / (density*ratio[1]*bond_len*sin(angle/2) * ratio[2]*bond_len*cos(angle/2))
    ]
    atom_vec = [
        0          0          0.5
        1/ratio[1] 1/ratio[2] 0.5
        2/ratio[1] 0          0.5
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 + 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
        1/ratio[1] 1/ratio[2]+1/cell_vec[2]*dummy_len*cos(dummy_ang/2) 0.5 - 1/cell_vec[3]*dummy_len*sin(dummy_ang/2)
        1/ratio[1]/2 1/ratio[2]/2 0.5
        3/ratio[1]/2 1/ratio[2]/2 0.5
    ]
    cell_vec = diag(cell_vec)

    # Setting Bond Modes
    bond_mode = Vector{Bond}(undef, 6)
    bond_mode[1] = Bond(0, 1, [1, 6])
    bond_mode[2] = Bond(0, 1, [3, 7])
    bond_mode[3] = Bond(0, 1, [2, 4])
    bond_mode[4] = Bond(0, 1, [2, 5])
    bond_mode[5] = Bond(0, 1, [2, 6])
    bond_mode[6] = Bond(0, 1, [2, 7])
    num_bonds = length(bond_mode)
    num_bond_types = 1

    # Setting Angle Modes
    angle_mode = Vector{Angle}(undef, 2)
    angle_mode[1] = Angle(1, 1, [1, 2, 3])
    angle_mode[2] = Angle(2, 1, [2, 4, 5])
    num_angles = length(angle_mode)
    num_angle_types = 2

    res_name = "TIP7P"
    Structure_Wat(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms,
    num_atom_types, bond_mode, num_bonds, num_bond_types, angle_mode, num_angles, num_angle_types, [1, 2], res_name)
end


text:
using BNBKnapsack
using Test

@testset "BNBKnapsack.jl" begin
    # Execute the BNB
    #@test BNBKnapsack.execute(BNBKnapsack.Problem(4,[50, 40, 45, 55],[5, 4, 3, 4],13,1e-4, 1e-3)) == 0
end

n_objects = 4               # Also the number of variables.
values = [16, 22, 12, 8]
weights = [5, 7, 4, 3]
max_weight = 14

integer_precision = 1e-4    # Integers are represented as floating-point values: if they are close enough to an integer, they are considered integer.
gap_precision = 1e-3

problem = BNBKnapsack.Problem(n_objects,values,weights,max_weight,integer_precision,gap_precision)
BNBKnapsack.execute(problem)

ext: jl


text:
function plotrg!(
    ax::Axis,
    r::AbstractMatrix,
    colnames::AbstractVector;
    p::Union{AbstractMatrix, Nothing} = nothing,
    circle::Bool = true,
    diag::Bool = false,
    )

    n = size(r, 1)
    ps = Vector{Polygon}(undef, n^2)
    counter = 1
    for i in 1:n
        for j in 1:n
            ps[counter] = Polygon([
                Point2f(i - 1, 1 - j),
                Point2f(i - 1, -j),
                Point2f(i, -j),
                Point2f(i, 1 - j)
            ])
            counter += 1
        end
    end
    psd = Vector{Polygon}(undef, n)
    for i in 1:n
        psd[i] = Polygon([
            Point2f(i - 1, 1 - i),
            Point2f(i - 1, -i),
            Point2f(i, -i),
            Point2f(i, 1 - i)
        ])
    end
    if circle
        csl = Vector{Circle}(undef, binomial(n, 2))
        csu = Vector{Circle}(undef, binomial(n, 2))
        colorl = Vector{Float32}(undef, binomial(n, 2))
        coloru = Vector{Float32}(undef, binomial(n, 2))
        counter = 1
        for i in 1:n
            for j in (i + 1):n
                csl[counter] = Circle(Point2f(j - 0.5, 0.5 - i), r[i, j] * 0.45)
                csu[counter] = Circle(Point2f(i - 0.5, 0.5 - j), r[j, i] * 0.45)
                colorl[counter] = r[i, j]
                coloru[counter] = r[j, i]
                counter += 1
            end
        end
        if diag
            csd = Vector{Circle}(undef, n)
            colord = Vector{Float32}(undef, n)
            for i in 1:n
                csd[i] = Circle(Point2f(i - 0.5, 0.5 - i), r[i, i] * 0.45)
                colord[i] = r[i, i]
            end
        end
    else
        csl = Vector{Rect}(undef, binomial(n, 2))
        csu = Vector{Rect}(undef, binomial(n, 2))
        colorl = Vector{Float32}(undef, binomial(n, 2))
        coloru = Vector{Float32}(undef, binomial(n, 2))
        counter = 1
        for i in 1:n
            for j in (i + 1):n
                widthl = r[i, j] * 0.95
                widthu = r[j, i] * 0.95
                csl[counter] = Rect(j - 1 + (1 - widthl) / 2, 1 - i - (1 - widthl) / 2, widthl, -widthl)
                csu[counter] = Rect(i - 1 + (1 - widthu) / 2, 1 - j - (1 - widthu) / 2, widthu, -widthu)
                colorl[counter] = r[i, j]
                coloru[counter] = r[j, i]
                counter += 1
            end
        end
        if diag
            csd = Vector{Circle}(undef, n)
            colord = Vector{Float32}(undef, n)
            for i in 1:n
                widthd = r[i, i] * 0.95
                csd[i] = Rect(i - 1 + (1 - withd) / 2, 1 - i - (1 - withd) / 2, withd, -withd)
                colord[i] = r[i, i]
            end
        end
    end
    poly!(ax, ps, color = :white, strokewidth = 0.5, strokecolor = "#C3C3C3")
    # poly!(ax, psd, color = "#C3C3C3", strokewidth = 0.5, strokecolor = "#C3C3C3")
    poly!(ax, csl, color = colorl, colorrange = (-1, 1), colormap = :RdBu_10)
    poly!(ax, csu, color = coloru, colorrange = (-1, 1), colormap = :RdBu_10)
    if !isnothing(p)
        pvalues = Matrix{Float64}(undef, n^2, 3)
        counter = 1
        for i in 1:n
            for j in 1:n
                j == i ? pvalues[counter, 1] = 1 : pvalues[counter, 1] = p[j, i]
                pvalues[counter, 2] = i - 0.5
                pvalues[counter, 3] = -j + 0.5
                counter += 1
            end
        end
        ind = findall(x -> x < 0.05, pvalues[:, 1])
        scatter!(ax, pvalues[ind, 2], pvalues[ind, 3], marker = '✳', markersize = 7, color = :black)
    end
    if diag
        poly!(ax, csd, color = colord, colorrange = (-1, 1), colormap = :RdBu_10)
    end
    ax.xticks = (0.5:(n - 0.5), colnames)
    ax.yticks = (-0.5:-1:-(n - 0.5), colnames)
    ax.xticklabelsize = 6
    ax.yticklabelsize = 6
    ax.xticklabelpad = 0
    ax.yticklabelpad = 1
    ax.xaxisposition = :top
    ax.yaxisposition = :left
    hidedecorations!(ax, ticklabels = false)
    xlims!(ax, 0, n)
    ylims!(ax, -n, 0)
    ax.aspect = DataAspect()
end


text:
using SurrealNumbers: EmptySurrealSet, Singleton, WrappedSurrealSet
using SurrealNumbers: issetfinite, all_lt, all_lte

@testset "Conversion" begin
	@test SurrealSet() === ∅

	# Wraps single numbers in Singleton
	@test SurrealSet(1) === Singleton(1)

	# Wraps collections in WrappedSurrealSet
	s = [1, 2, 3]
	@test SurrealSet(s) === WrappedSurrealSet(s)

	# Passes existing instances through
	for s in [∅, Singleton(1), WrappedSurrealSet(1:3)]
		@test SurrealSet(s) === s
	end
end

@testset "EmptySurrealSet" begin
	@test Base.IteratorSize(∅) == Base.HasLength()
	@test length(∅) == 0
	@test isempty(∅)
	@test issetfinite(∅)
	@test eltype(∅) === Union{}

	@test iterate(∅) === nothing

	@test -∅ === ∅
end

@testset "Singleton" begin
	for value in [1, SurrealForm(0)]
		s = Singleton(value)

		@test Base.IteratorSize(s) == Base.HasLength()
		@test length(s) == 1
		@test issetfinite(s)
		@test eltype(s) == typeof(value)
		@test collect(s) == [value]
		@test only(s) === value
	end
end

@testset "WrappedSurrealSet" begin
	for s in []
		w = WrappedSurrealSet(s)

		@test Base.IteratorSize(w) == Base.IteratorSize(s)
		@test length(w) == length(s)
		@test eltype(w) == eltype(s)
		@test issetequal(w, s)
	end
end

@testset "Ordering" begin
	s1 = SurrealSet(1:3)
	s2 = SurrealSet(3:6)
	s3 = SurrealSet(4:6)
	s4 = SurrealSet(2:6)

	@test !all_lt(s1, s2)
	@test all_lte(s1, s2)
	@test !all_lt(s2, s1)
	@test !all_lte(s2, s1)

	@test all_lt(s1, s3)
	@test all_lte(s1, s3)
	@test !all_lt(s3, s1)
	@test !all_lte(s3, s1)

	@test !all_lt(s1, s4)
	@test !all_lte(s1, s4)
	@test !all_lt(s4, s1)
	@test !all_lte(s4, s1)
end


text:
export ActivateFuncs

module ActivateFuncs
import Base.sign
using QuaternionArgs
using Quaternions

export sign, quantization

# split type sign functions {{{
function sign{T<:Real}(z::Complex{T})
    return complex(sign(real(z)) ,sign(imag(z)))
end

function sign{T<:Real}(q::Quaternion{T})
    return Quaternion{T}(sign(q.q0),sign(q.q1),sign(q.q2),sign(q.q3))
end
# }}}

# quantization functions for multistate neurons {{{
function quantization{T<:Real}(z::Complex{T},K::Integer)
    halfphi = pi/K
    phi = 2.0*halfphi
    theta = arg(z) + halfphi
    q = floor(theta / phi)
    w = q*phi
    return complex(cos(w),sin(w))
end
function arg{T<:Real}(z::Complex{T})
    return atan2(imag(z),real(z))
end

function quantization_{T<:Real}(q::Quaternion{T},A::Integer,B::Integer,C::Integer)
    qarg = QuaternionArg(q)
    if randbool()
        phi = qarg.phi
        theta = qsign(qarg.theta, 0.5, B)
        psi = qsign(qarg.psi, 0.25, C)
    else
        phi = qsign(qarg.phi, 1.0, A)
        theta = qsign(qarg.theta, 0.5, B)
        psi = qarg.psi
    end
    return Quaternion(QuaternionArg(one(T),phi,theta,psi))
end

function quantization{T<:Real}(q::Quaternion{T},A::Integer,B::Integer,C::Integer)
    qarg = QuaternionArg(q)

    phi = qsign(qarg.phi, 1.0, A)
    theta = qsign(qarg.theta, 0.5, B)
    psi = qsign(qarg.psi, 0.25, C)

    return Quaternion(QuaternionArg(one(T),phi,theta,psi))
end

function qsign(phase::Real, coef::Real, K::Integer)
    dphase = 2*pi/K*coef
    phase0 = -pi*coef
    flg = false
    p = phase0

    # phase regulation
    phase = phase < phase0 ? phase0 : (phase >= -phase0 ? -phase0-0.5*dphase : phase)
    for i=1:K
        if p <= phase && phase < p+dphase
            flg = true
            break
        end
        p = phase0+dphase*i
    end
    @assert(flg)

    return p+0.5*dphase
end
# }}}

end


text:
import MCMCChains
using CmdStan

const noncentered_schools_stan_model = """
    data {
        int<lower=0> J;
        real y[J];
        real<lower=0> sigma[J];
    }
    parameters {
        real mu;
        real<lower=0> tau;
        real eta[J];
    }
    transformed parameters {
        real theta[J];
        for (j in 1:J)
            theta[j] = mu + tau * eta[j];
    }
    model {
        mu ~ normal(0, 5);
        tau ~ cauchy(0, 5);
        eta ~ normal(0, 1);
        y ~ normal(theta, sigma);
    }
    generated quantities {
        vector[J] log_lik;
        vector[J] y_hat;
        for (j in 1:J) {
            log_lik[j] = normal_lpdf(y[j] | theta[j], sigma[j]);
            y_hat[j] = normal_rng(theta[j], sigma[j]);
        }
    }
"""

function makechains(names, ndraws, nchains; seed = 42, internal_names = [])
    rng = MersenneTwister(seed)
    nvars = length(names)
    vals = randn(rng, ndraws, nvars, nchains)
    chns = MCMCChains.Chains(vals, names, Dict(:internals => internal_names))
    return chns
end

function makechains(nvars::Int, args...; kwargs...)
    names = ["var$(i)" for i = 1:nvars]
    return makechains(names, args...; kwargs...)
end

dimsizes(ds) = ds._dims
convertindex(x::AbstractArray) = x
convertindex(o::PyObject) = o.array.values
vardict(ds) = Dict(k => convertindex(v._data) for (k, v) in ds._variables)
dimdict(ds) = Dict(k => v._dims for (k, v) in ds._variables)
attributes(ds) = ds.attrs

function cmdstan_noncentered_schools(data, draws, chains; proj_dir = pwd())
    model_name = "school8"
    stan_model = Stanmodel(
        name = model_name,
        model = noncentered_schools_stan_model,
        nchains = chains,
        num_warmup = draws,
        num_samples = draws,
    )
    rc, chns, cnames = stan(stan_model, data, proj_dir, summary = false)
    outfiles = []
    for i = 1:chains
        push!(outfiles, "$(proj_dir)/tmp/$(model_name)_samples_$(i).csv")
    end
    return (model = stan_model, files = outfiles, chains = chns)
end

function test_chains_data(chns, idata, group, names; coords = Dict(), dims = Dict())
    ndraws, nvars, nchains = size(chns)
    @test idata isa InferenceData
    @test group in propertynames(idata)
    ds = getproperty(idata, group)
    sizes = dimsizes(ds)
    @test length(sizes) == 2 + length(coords)
    vars = vardict(ds)
    for name in names
        @test name in keys(vars)
        dim = get(dims, name, [])
        s = (x -> length(get(coords, x, []))).(dim)
        @test size(vars[name]) == (nchains, ndraws, s...)
    end
    @test attributes(ds)["inference_library"] == "MCMCChains"
end

@testset "from_mcmcchains" begin
    @testset "posterior" begin
        nvars, nchains, ndraws = 2, 4, 20
        chns = makechains(nvars, ndraws, nchains)
        idata = from_mcmcchains(chns)
        test_chains_data(chns, idata, :posterior, names(chns))
    end

    @testset "prior" begin
        nvars, nchains, ndraws = 2, 4, 20
        chns = makechains(nvars, ndraws, nchains)
        idata = from_mcmcchains(prior = chns)
        test_chains_data(chns, idata, :prior, names(chns))
    end

    @testset "posterior + prior" begin
        nvars, nchains, ndraws = 2, 4, 20
        chns = makechains(nvars, ndraws, nchains)
        idata = from_mcmcchains(chns; prior = chns)
        test_chains_data(chns, idata, :posterior, names(chns))
        test_chains_data(chns, idata, :prior, names(chns))
    end

    @testset "coords/dim" begin
        names = ["a[1]", "a[2]", "b[1]", "b[2]"]
        coords = Dict("ai" => 1:2, "bi" => ["b1", "b2"])
        dims = Dict("a" => ["ai"], "b" => ["bi"])
        nchains, ndraws = 4, 20
        chns = makechains(names, ndraws, nchains)
        idata = from_mcmcchains(chns; coords = coords, dims = dims)
        test_chains_data(chns, idata, :posterior, ["a", "b"]; coords = coords, dims = dims)
        vardims = dimdict(idata.posterior)
        @test vardims["a"] == ("chain", "draw", "ai")
        @test vardims["b"] == ("chain", "draw", "bi")
    end

    @testset "multivariate" begin
        names = ["a[1][1]", "a.2.2", "a[2,1]", "a[1, 2]"]
        coords = Dict("ai" => 1:2, "aj" => ["aj1", "aj2"])
        dims = Dict("a" => ["ai", "aj"])
        nchains, ndraws = 4, 20
        chns = makechains(names, ndraws, nchains)
        idata = from_mcmcchains(chns; coords = coords, dims = dims)
        test_chains_data(chns, idata, :posterior, ["a"]; coords = coords, dims = dims)
        arr = vardict(idata.posterior)["a"]
        @test arr[:, :, 1, 1] == permutedims(chns.value[:, names[1], :], [2, 1])
        @test arr[:, :, 2, 2] == permutedims(chns.value[:, names[2], :], [2, 1])
        @test arr[:, :, 2, 1] == permutedims(chns.value[:, names[3], :], [2, 1])
        @test arr[:, :, 1, 2] == permutedims(chns.value[:, names[4], :], [2, 1])
    end

    @testset "complete" begin
        nchains, ndraws = 4, 20
        posterior = prior = ["a[1]", "a[2]"]
        posterior_predictive = prior_predictive = ["ahat[1]", "ahat[2]"]
        observed_data = ["xhat[1]", "xhat[2]"]
        constant_data = ["x[1]", "x[2]"]
        sample_stats = ["stat"]
        log_likelihood = "log_lik"
        post_names = [
            posterior
            posterior_predictive
            observed_data
            constant_data
            sample_stats
            log_likelihood
        ]
        prior_names = [prior; prior_predictive; sample_stats]
        chns = makechains(post_names, ndraws, nchains; internal_names = ["stat"])
        chns2 = makechains(prior_names, ndraws, nchains; internal_names = ["stat"])
        idata = from_mcmcchains(
            chns;
            posterior_predictive = "ahat",
            prior = chns2,
            prior_predictive = "ahat",
            observed_data = "xhat",
            constant_data = "x",
            log_likelihood = "log_lik",
        )
        for group in (
            :posterior,
            :prior,
            :posterior_predictive,
            :prior_predictive,
            :observed_data,
            :constant_data,
            :sample_stats,
            :sample_stats_prior,
        )
            @test group in propertynames(idata)
        end
        @test length(dimdict(idata.posterior)) == 4
        @test "a" ∈ keys(dimdict(idata.posterior))
        @test length(dimdict(idata.posterior_predictive)) == 4
        @test "ahat" ∈ keys(dimdict(idata.posterior_predictive))
        @test length(dimdict(idata.prior)) == 4
        @test "a" ∈ keys(dimdict(idata.prior))
        @test length(dimdict(idata.prior_predictive)) == 4
        @test "ahat" ∈ keys(dimdict(idata.prior_predictive))
        @test length(dimdict(idata.sample_stats)) == 4
        @test "stat" ∈ keys(dimdict(idata.sample_stats))
        @test "log_likelihood" ∈ keys(dimdict(idata.sample_stats))
        @test length(dimdict(idata.sample_stats_prior)) == 3
        @test "stat" ∈ keys(dimdict(idata.sample_stats_prior))
        @test length(dimdict(idata.observed_data)) == 4
        @test "xhat" ∈ keys(dimdict(idata.observed_data))
        @test length(dimdict(idata.constant_data)) == 4
        @test "x" ∈ keys(dimdict(idata.constant_data))
    end

    @testset "missing -> NaN" begin
        rng = MersenneTwister(42)
        nvars, nchains, ndraws = 2, 4, 20
        vals = randn(rng, ndraws, nvars, nchains)
        vals = Array{Union{Float64,Missing},3}(vals)
        vals[1, 1, 1] = missing
        names = ["var$(i)" for i = 1:nvars]
        chns = MCMCChains.Chains(vals, names)
        @test Missing <: eltype(chns.value)
        idata = from_mcmcchains(chns)
        vdict = vardict(idata.posterior)
        @test eltype(vdict["var1"]) <: Real
        @test isnan(vdict["var1"][1, 1, 1])
    end

    @testset "info -> attributes" begin
        nvars, nchains, ndraws = 2, 4, 20
        chns = makechains(nvars, ndraws, nchains)
        idata = from_mcmcchains(chns; library = "MyLib")
        attrs = idata.posterior.attrs
        @test attrs isa Dict
        @test attrs["inference_library"] == "MyLib"
        @test attrs["mcmcchains_summary"] isa Dict
        mcmcchains_summary = idata.posterior.attrs["mcmcchains_summary"]
        @test first(values(mcmcchains_summary)).__class__.__name__ == "DataFrame"
    end
end

@testset "convert_to_dataset(::MCMCChains.Chains)" begin
    nvars, nchains, ndraws = 2, 4, 20
    chns = makechains(nvars, ndraws, nchains)
    ds = ArviZ.convert_to_dataset(chns; library = "MyLib")
    @test ds isa ArviZ.Dataset
    attrs = ds.attrs
    @test "inference_library" ∈ keys(attrs)
    @test attrs["inference_library"] == "MyLib"
end

@testset "convert_to_inference_data(::MCMCChains.Chains)" begin
    nvars, nchains, ndraws = 2, 4, 20
    chns = makechains(nvars, ndraws, nchains)
    idata = convert_to_inference_data(chns; group = :posterior)
    @test idata isa InferenceData
    @test :posterior in propertynames(idata)
    idata = convert_to_inference_data(chns; group = :prior)
    @test idata isa InferenceData
    @test :prior in propertynames(idata)
end

if VERSION.minor > 0
    @testset "from_cmdstan" begin
        data = noncentered_schools_data()
        output = cmdstan_noncentered_schools(data, 500, 4)
        posterior_predictive = prior_predictive = ["y_hat"]
        log_likelihood = "log_lik"
        coords = Dict("school" => 1:8)
        dims = Dict(
            "theta" => ["school"],
            "y" => ["school"],
            "log_lik" => ["school"],
            "y_hat" => ["school"],
            "eta" => ["school"],
        )
        idata1 = from_cmdstan(
            output.chains;
            posterior_predictive = posterior_predictive,
            log_likelihood = log_likelihood,
            prior = output.chains,
            prior_predictive = prior_predictive,
            coords = coords,
            dims = dims,
        )
        idata2 = from_cmdstan(
            output.files;
            posterior_predictive = posterior_predictive,
            log_likelihood = log_likelihood,
            prior = output.files,
            prior_predictive = prior_predictive,
            coords = coords,
            dims = dims,
        )
        @testset "idata.$(group)" for group in Symbol.(idata2._groups)
            ds1 = getproperty(idata1, group)
            ds2 = getproperty(idata2, group)
            for f in (vardict, dimsizes)
                d1 = f(ds1)
                d2 = f(ds2)
                for k in keys(d1)
                    @test k in keys(d2)
                    @test d1[k] ≈ d2[k]
                end
            end
        end
    end
end


text:
"""
    abstract type Equation{m} end

Abstract supertype for equations that can be solved with this package.

# Type parameters
- `m`: number of spatial dimensions:
    - 1 for a non-radial one-dimensional equation (default);
    - 2 for a radial equation in polar or cylindrical coordinates;
    - 3 for a radial equation in spherical coordinates.
"""
abstract type Equation{m} end

"""
    isindomain(eq::Equation, val) -> Bool

`true` if `eq` is well defined for the solution value `val`; `false` otherwise.
"""
function isindomain end

"""
    DiffusionEquation(D; symbol=:θ) <: Equation{1}
    DiffusionEquation{m}(D; symbol=:θ) <: Equation{m}

Nonlinear diffusion equation.

    DiffusionEquation(model) <: Equation{1}
    DiffusionEquation{m}(model) <: Equation{m}

Nonlinear diffusion equation describing flow in a porous medium, with the diffusivity defined by a porous model.

# Arguments
- `D`: diffusivity function.
- `model::PorousModels.UnsaturatedFlowModel`: unsaturated flow model from which to obtain the diffusivity function.

# Keyword arguments
- `symbol::Symbol=:θ`: optional symbol used to represent the unknown function in the output.

# Type parameters
- `m::Int=1`: number of spatial dimensions:
    - 1 for non-radial one-dimensional diffusion (default);
    - 2 for radial diffusion in polar or cylindrical coordinates;
    - 3 for radial diffusion in spherical coordinates.

# Examples
```jldoctest; setup = :(using Fronts)
julia> D(θ) = θ^4
D (generic function with 1 method)

julia> eq = Fronts.DiffusionEquation(D)
∂θ/∂t = ∂(D(θ)*∂θ/∂r)/∂r

julia> eq = Fronts.DiffusionEquation{2}(D)
∂θ/∂t = 1/r*∂(r*D(θ)*∂θ/∂r)/∂r

julia> eq = Fronts.DiffusionEquation{3}(D, symbol=:c)
∂c/∂t = 1/r²*∂(r²*D(c)*∂c/∂r)/∂r
```

See also: [`PorousModels.UnsaturatedFlowModel`](@ref)
"""
struct DiffusionEquation{m,_TD} <: Equation{m}
    D::_TD
    symbol::Symbol

    function DiffusionEquation{m}(D; symbol::Symbol=:θ) where m
        @argcheck typeof(m) == Int TypeError(:m, Int, m)
        @argcheck m::Int in (1,2,3)
        new{m,typeof(D)}(D, symbol)
    end
end

DiffusionEquation(D; symbol::Symbol=:θ) = DiffusionEquation{1}(D, symbol=symbol)

function DiffusionEquation{m}(model::PorousModels.UnsaturatedFlowModel) where m
    function D(θ)
        PorousModels.Dθ(model, θ)
    end
    DiffusionEquation{m}(D)
end

DiffusionEquation(model::PorousModels.UnsaturatedFlowModel) = DiffusionEquation{1}(model)

function Base.show(io::IO, eq::DiffusionEquation{1})
    print(io, "∂", eq.symbol, "/∂t = ∂(", eq.D, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

function Base.show(io::IO, eq::DiffusionEquation{2})
    print(io, "∂", eq.symbol, "/∂t = 1/r*∂(r*", eq.D, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

function Base.show(io::IO, eq::DiffusionEquation{3})
    print(io, "∂", eq.symbol, "/∂t = 1/r²*∂(r²*", eq.D, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

flux(eq::DiffusionEquation, θ, r, t) = -eq.D(θ(r, t))*∂_∂r(θ, r, t)

function isindomain(eq::DiffusionEquation, θ)
    D = NaN
    dD_dθ = NaN
    try
        D, dD_dθ = value_and_derivative(eq.D, θ)
    catch e
        isa(e, DomainError) || rethrow()
    end

    return isfinite(D) && D>zero(D) && isfinite(dD_dθ)
end

"""
    RichardsEquation(; C, K, symbol=:h) <: Equation{1}
    RichardsEquation{m}(; C, K, symbol=:h) <: Equation{m}

Horizontal Richards equation, pressure-based formulation.

    RichardsEquation(model) <: Equation{1}
    RichardsEquation{m}(model) <: Equation{m}

Horizontal Richards equation, pressure-based formulation, with properties defined by a porous model.

# Arguments
- `pm::PorousModels.UnsaturatedFlowModel`: unsaturated flow model from which to obtain the relevant functions.

# Keyword arguments
- `C`: hydraulic capacity function, defined in terms of the unknown.
- `K`: hydraulic conductivity function, defined in terms of the unknown.
- `symbol::Symbol=:h`: optional symbol used to represent the unknown function in the output.

# Type parameters
- `m::Int=1`: number of spatial dimensions:
    - 1 for non-radial one-dimensional flow (default);
    - 2 for radial flow in polar or cylindrical coordinates;
    - 3 for radial flow in spherical coordinates.

See also: [`PorousModels.UnsaturatedFlowModel`](@ref)
"""
struct RichardsEquation{m,_TC,_TK} <: Equation{m}
    C::_TC
    K::_TK
    symbol::Symbol

    function RichardsEquation{m}(; C, K, symbol::Symbol=:h) where m
        @argcheck typeof(m) == Int TypeError(:m, Int, m)
        @argcheck m::Int in (1,2,3)
        new{m,typeof(C),typeof(K)}(C, K, symbol)
    end
end

RichardsEquation(; C, K, symbol::Symbol=:h) = RichardsEquation{1}(C=C, K=K, symbol=symbol)

function RichardsEquation{m}(model::PorousModels.UnsaturatedFlowModel) where m
    function C(h)
        PorousModels.Ch(model, h)
    end
    function K(h)
        PorousModels.Kh(model, h)
    end
    RichardsEquation{m}(; C=C, K=K)
end

RichardsEquation(model::PorousModels.UnsaturatedFlowModel) = RichardsEquation{1}(model)

function Base.show(io::IO, eq::RichardsEquation{1})
    print(io, eq.C, "*∂", eq.symbol, "/∂t = ∂(", eq.K, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

function Base.show(io::IO, eq::RichardsEquation{2})
    print(io, eq.C, "*∂", eq.symbol, "/∂t = 1/r*∂(r*", eq.K, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

function Base.show(io::IO, eq::RichardsEquation{3})
    print(io, eq.C, "*∂", eq.symbol, "/∂t = 1/r²*∂(r²*", eq.K, "(", eq.symbol, ")*∂", eq.symbol, "/∂r)/∂r")
end

flux(eq::RichardsEquation, h, r, t) = -eq.K(h(r, t))*∂_∂r(h, r, t)

function isindomain(eq::RichardsEquation, h)
    C = NaN
    K = NaN
    dK_dh = NaN
    try
        C = eq.C(h)
        K, dK_dh = value_and_derivative(eq.K, h)
    catch e
        isa(e, DomainError) || rethrow()
    end

    return isfinite(C) && isfinite(K) && K>zero(K) && isfinite(dK_dh)
end


text:
"""
    Baselet.Specialized

Functions that are specialized in Baselet.jl.
"""
baremodule Specialized
function getindex end
function filter end
function zip end
function enumerate end
function accumulate end
function cumsum end
function cumprod end
function flatten end
function foreach end
function any end
function all end
function findfirst end
function findlast end
function findall end
function findmax end
function findmin end
function argmax end
function argmin end
function maximum end
function minimum end
function extrema end
function in end
function unique end
function union end
function intersect end
function setdiff end
function symdiff end
function issubset end
function isdisjoint end
function sort end
end


text:
export make_atoms
export collect_atom_groups

function make_atoms(atom_xy)
    atom_points = map(Point2D, atom_xy[1,:], atom_xy[2,:])
    indexed_atoms = pairs(IndexLinear(), atom_points)
    return indexed_atoms
end

function collect_atom_groups(atom_xy::Array{Float64,2}; max_bondlength=12)
    kdtree = KDTree(atom_xy, leafsize = 10)
    natoms = size(atom_xy, 2)
    indices, bondlengths = knn(kdtree, atom_xy, min(natoms, 4), true)
    indices_inrange = inrange(kdtree, atom_xy, max_bondlength, true)
    atoms_knn = zip(indices, indices_inrange)
    map(atom_inrange, atoms_knn)
end

function collect_atom_groups(atom_xy::Array{Float64,1}; max_bondlength=12)
    return [[1]]
end
#
function atom_inrange(atom_knn)
    index, index_inrange = atom_knn
    atom_selected = sort!(delete_nothing(indexin(index_inrange, index)))
    index_new = index[atom_selected]
    return index_new
end

function delete_nothing(x)
    if typeof(x) <: Array{Union{Nothing, T}} where T <: Number
        convert(Array{Int64}, filter!(x -> !isnothing(x), x))
    else
        x
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

const _default_PT = Float32 # Default type for parameter values
const _default_float_WT = Float64 # Default type for float weights
const _default_int_WT = Int # Default type for int weights
const _default_LDT = Float64 # Default type for log-density values

const RandSampleable = Union{
    DistLikeDensity,
    MultivariateDistribution,
}

const AnyPosterior = Union{
    PosteriorDensity,
    DensitySampleVector,
    RandSampleable,
}

const AnyNSamples = Union{
    Integer,
    Tuple{Integer,Integer},
}

"""
    BAT.AbstractSamplingAlgorithm

Abstract type for BAT sampling algorithms. See [`bat_sample`](@ref).
"""
abstract type AbstractSamplingAlgorithm end

"""
    BAT.default_sampling_algorithm(posterior)

Get BAT's default sampling algorithm for `posterior`.
"""
function default_sampling_algorithm end

"""
    bat_sample(
        [rng::AbstractRNG],
        posterior::BAT.AnyPosterior,
        n::BAT.AnyNSamples,
        [algorithm::BAT.AbstractSamplingAlgorithm]
    )::DensitySampleVector

Draw `n` samples from `posterior`.

Returns a NamedTuple of the shape

```julia
(
    samples = X::DensitySampleVector,...
    stats = s::@test stats isa NamedTuple{(:mode,:mean,:cov,...)},
    ...
)
```

Result properties not listed here are algorithm-specific and are not part
of the stable BAT API.

`posterior` may be a

* [`BAT.AbstractPosteriorDensity`](@ref)

* [`BAT.DistLikeDensity`](@ref)

* [`BAT.DensitySampleVector`](@ref)

* `Distributions.MultivariateDistribution`

Depending on the type of posterior, `n` may be of type

* `Integer`: Number of samples

* `Tuple{Integer,Integer}`: Tuple of number of samples per sample source
   and number of sample sources (e.g. number of MCMC chains). The total number
   of samples is `product(n)`.

Depending on the type of `posterior`, the number of samples returned may be
somewhat larger or smaller than specified by `product(n)`.

Also depending on the `posterior` type, the samples may be independent or
correlated (e.g. when using MCMC).
"""
function bat_sample end
export bat_sample

@inline function bat_sample(
    posterior::AnyPosterior, n::AnyNSamples;
    kwargs...
)
    rng = bat_rng()
    bat_sample(rng, posterior, n; kwargs...)
end

@inline function bat_sample(
    posterior::AnyPosterior, n::AnyNSamples, algorithm::AbstractSamplingAlgorithm;
    kwargs...
)
    rng = bat_rng()
    bat_sample(rng, posterior, n, algorithm; kwargs...)
end

@inline function bat_sample(
    rng::AbstractRNG, posterior::AnyPosterior, n::AnyNSamples;
    kwargs...
)
    algorithm = default_sampling_algorithm(posterior)
    bat_sample(rng, posterior, n, algorithm; kwargs...)
end

"""
    RandSampling

Constructors:

    RandSampling()

Sample via `Random.rand`. Only supported for posteriors of type
`Distributions.MultivariateDistribution` and `BAT.DistLikeDensity`.
"""
struct RandSampling <: AbstractSamplingAlgorithm end
export RandSampling

default_sampling_algorithm(posterior::RandSampleable) = RandSampling()

function bat_sample(rng::AbstractRNG, posterior::RandSampleable, n::Integer, algorithm::RandSampling)
    npar = length(posterior)
    samples = DensitySampleVector{_default_PT,_default_LDT,_default_int_WT,Nothing,Nothing}(undef, n, npar)

    rand!(rng, sampler(posterior), flatview(samples.params))
    let logdensity = samples.logdensity, params = samples.params
        @uviews logdensity .= logpdf.(Ref(posterior), params)
    end
    samples.weight .= 1

    stats = bat_stats(samples)

    (samples = samples, stats = stats)
end

"""
    RandResampling <: AbstractSamplingAlgorithm

Constructors:

    RandResampling()

Resamples from a given set of samples.
"""
struct RandResampling <: AbstractSamplingAlgorithm end
export RandResampling

default_sampling_algorithm(posterior::DensitySampleVector) = RandResampling()

function bat_sample(rng::AbstractRNG, posterior::DensitySampleVector, n::Integer, algorithm::RandResampling)
    orig_idxs = eachindex(posterior)
    weights = FrequencyWeights(float(posterior.weight))
    resampled_idxs = sample(orig_idxs, weights, n, replace=true, ordered=false)

    samples = posterior[resampled_idxs]
    samples.weight .= 1

    stats = bat_stats(samples)

    (samples = samples, stats = stats)
end


text:
"""
Although the BSE exposes a special SOAP-based API for programmatic access to its database, I had too much trouble getting the desired results from it. Therefore the solution in this module is based around regexing of the
human readable HTML output.
"""
module BasisSetExchangeModule

export obtainBasisSetExchangeEntries, downloadBasisSetBasisSetExchange, computeBasisSetExchangeEntry, BasisSetExchangeEntry
using Requests
import Base.display

type BasisSetExchangeEntry
  url::UTF8String
  name::UTF8String
  availtype::UTF8String
  elts::UTF8String
  status::UTF8String
  hasEcp::UTF8String
  hasSpin::UTF8String
  lastModifiedDate::UTF8String
  contributionPI::UTF8String
  contributorName::UTF8String
  bsAbstract::UTF8String
end

function display(bseEntry::BasisSetExchangeEntry)
  @printf("%40s    |    %s\n", bseEntry.name, bseEntry.bsAbstract)
end

function display(arr::Array{BasisSetExchangeEntry,1})
  for entry in arr
    display(entry)
  end
end

function obtainBasisSetExchangeEntries()
  retval = Array{BasisSetExchangeEntry,1}()
  main = get("https://bse.pnl.gov/bse/portal/user/anon/js_peid/11535052407933/action/portlets.BasisSetAction/template/courier_content/panel/Main")
  for bsline in matchall(r"basisSets\[\d*\].*",Requests.text(main))
    regmatch = match(r"basisSet\(\"(?P<url>.*)\", \"(?P<name>.*)\", \"(?P<type>.*)\", \"\[(?P<elts>.*)\]\", \"(?P<status>.*)\", \"(?P<hasEcp>.*)\", \"(?P<hasSpin>.*)\", \"(?P<lastModifiedDate>.*)\", \"(?P<contributionPI>.*)\", \"(?P<contributorName>.*)\", \"(?P<bsAbstract>.*)\"\)",bsline)
    eltsList = replace(regmatch[:elts],",","")
    push!(retval,BasisSetExchangeEntry(regmatch[:url],regmatch[:name],regmatch[:type],eltsList,regmatch[:status],regmatch[:hasEcp],regmatch[:hasSpin],regmatch[:lastModifiedDate],regmatch[:contributionPI],regmatch[:contributorName],regmatch[:bsAbstract]))
  end
  return retval
end

function downloadBasisSetBasisSetExchange(entry::BasisSetExchangeEntry, filename::AbstractString, format::AbstractString="TX93")
  requestresponse = post("https://bse.pnl.gov:443/bse/portal/user/anon/js_peid/11535052407933/action/portlets.BasisSetAction/template/courier_content/panel/Main/eventSubmit_doDownload/true";
			data=Dict("bsurl" => entry.url, "bsname" => entry.name, "elts" => entry.elts, "format" => format, "minimize" => "false"), allow_redirects=false)
  sessioncookie = requestresponse.cookies["JSESSIONID"].value
  finalresponse = post("https://bse.pnl.gov/bse/portal/user/anon/panel/Main/template/courier_content/js_peid/11535052407933"; cookies = Dict("JSESSIONID" => sessioncookie))
  basSetDef = replace(Requests.text(finalresponse),r".*<pre style.*>\n*(.*)</pre>.*"s,s"\1")
  outfd = open(filename, "w")
  write(outfd, basSetDef)
  close(outfd)
end

function computeBasisSetExchangeEntry(name::AbstractString, arr::Array{BasisSetExchangeEntry,1})
  result = Void
  for entry in arr
    if ismatch(Regex(name,"i"),entry.name)
      if result == Void
        result = entry
      elseif isa(result,BasisSetExchangeEntry)
        result = [result,entry]
      else #isa(result,Array{BasisSetExchangeEntry})
        push!(result,entry)
      end
    end
  end
  if result==Void
    error("$name is not contained within the given BasisSetExchangeEntry list.")
  end
  return result
end

end # module


text:
module SGMCMC

    abstract SamplerState
    export AdamState, SGLDAdamState
    export SamplerState,HMCState,RelHMCState,RelSGHMCState, RelSGNHTState, SGHMCState, sample!
    export SGLDState, SGNHTState, mSGNHTState, SGDState, RelFrozenState, SantaState
    export SGLDAdagradState, AdagradState
    export SGLDRMSPropState

    add_inplace! = Base.LinAlg.axpy!
    #full gradient MCMC methods
    include("mcmc/hmc.jl")
    include("utils/ars.jl") #includes adaptive rejection sampling for momentum distribution.
    include("mcmc/relhmc.jl")

    #stochastic gradient MCMC methods
    include("sgmcmc/sgnht.jl")
    include("sgmcmc/msgnht.jl")
    include("sgmcmc/sghmc.jl")
    include("sgmcmc/relsghmc.jl")
    include("sgmcmc/relsgnht.jl")

    include("sgmcmc/sgld.jl")
    include("sgmcmc/sgldadam.jl")
    include("sgmcmc/sgldrmsprop.jl")
    include("sgmcmc/sgldadagrad.jl")

    include("sgd/sgd.jl")
    include("sgd/adam.jl")
    include("sgd/adagrad.jl")
    include("sgd/relfrozen.jl")

end


text:
using Reel
using Test
using Plots

let
    function render(t, dt)
        # t is the time into the sequence
        # dt is the time to advance for the next frame

        # any expression that results in an object which can be
        # rendered as png or jpg
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end

    film = roll(render, fps=30, duration=2)

    write("output.gif", film) # Write to a gif file
    write("file.webm", film) # Write to a webm video
    write("file.mp4", film)  # An mp4 formatted video

    @show typeof((render(t,0.1) for t in 0:0.1:2))
    roll((render(t,0.1) for t in 0:0.1:2), fps=30)
    write("output2.gif", film) # Write to a gif file
end

let
    film = roll(fps=30, duration=2) do t, dt
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end

    write("output.gif", film)
end

# Currently commented out because of problems even in Julia 0.6
# using Compose
# let
#     Compose.set_default_graphic_size(3inch, 3inch) # Square

#     # draw a regular n-gon
#     ngon(n) = compose(context(units=UnitBox(-1, -1, 2, 2)), fill("lightblue"),
#                        polygon([(cos(x), sin(x)) for x in π/2:2π/n:3.5π]))

#     roll(map(ngon, [3:10, 9:-1:3])), fps=5)

#     write("output.gif", film)
# end

let
    frames = Frames(MIME("image/png"), fps=2)

    for i=1:31
        t = i * 0.1
        push!(frames, plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6))
    end

    write("output.gif", frames)
    # write("output.webm", frames)
end

# issue 16  https://github.com/shashi/Reel.jl/issues/16
let
    ts = []
    film = roll(fps=1, duration=5) do t, dt
        push!(ts, t)
        plot([x -> sin(x+t*π), x -> cos(x+t*π)], 0, 6)
    end
    @test all(diff(ts).==1.0)
end


text:
#global

#adaptative
function niblack( image::Array{T,2}, width::Integer; k = 0.1 ) where {T<:Real}

    intA  = integralArray( image, type=Float64 );
    intA2 = integralArray( image, type=Float64, fun=(x)->(x^2));

    imgc  = zeros( Float32, size(image) );

    h, w = size(image);
    npix = width*width;

    for col in 1:w
        for row in 1:h

            TL = max( 1, row - width ), max( 1, col - width )
            BR = min( h, row + width ), min( w, col + width )

            npix2 = (length(TL[1]:BR[1])+1)*(length(TL[2]:BR[2])+1);

            sumI  = integralArea(  intA, TL, BR );
            sumI2 = integralArea( intA2, TL, BR );

            mean = sumI/npix2;
            m2   =  mean*mean;

            std  = sqrt( (sumI2 + npix2*m2 - 2*sumI*mean)/npix2 );

            th = mean + k*std;

            imgc[row,col] = ( image[row,col] <= th ) ? 1 : 0;
        end
    end

    return imgc;
end

function sauvola( image::Array{T,2}, width::Integer; R=1, k=0.3 ) where {T<:Real}

    intA  = integralArray( image, type=Float64 );
    intA2 = integralArray( image, type=Float64, fun=(x)->(x^2));

    imgc  = zeros( Float32, size(image) );

    h, w = size(image);
    npix = width*width;

    for col in 1:w
        for row in 1:h

            TL = max( 1, row - width ), max( 1, col - width )
            BR = min( h, row + width ), min( w, col + width )

            npix2 = (length(TL[1]:BR[1])+1)*(length(TL[2]:BR[2])+1);

            sumI  = integralArea(  intA, TL, BR );
            sumI2 = integralArea( intA2, TL, BR );

            mean = sumI/npix2;
            m2   =  mean*mean;

            std  = sqrt( (sumI2 + npix2*m2 - 2*sumI*mean)/npix2 );

            th = mean + k*( std/R - 1 )*(mean-1);

            imgc[row,col] = ( image[row,col] <= th ) ? 1 : 0;
            #imgc[row,col] = th
        end
    end

    return imgc;
end

heavensideStep( num ) = ( num <= 0 ) ? 0 : 1;

function adaptOtsu( image::Array{T,2}, width::Integer, step::Integer, bindata::AbstractRange; R=0.1 ) where {T<:Real}

    imgc = zeros( Float32, size(image) )
    img  = copy(image)
    h, w = size(image)

    len  = 2*width + 1;
    #step = 2*width + 1;
    centers_r = [ row - width for row in len:step:h ]
    centers_c = [ col - width for col in len:step:w ]

    thresholds = zeros( Float32, length(len:step:h), length(len:step:w) );

    nhImg,_,_ = normhistogram( imgc, bindata );
    imgOt = otsu( nhImg );

    # computing threshold for each grid element
    contc = 1;
    for col in len:step:w
        contr = 1;
        for row in len:step:h

            TL = max( 1, row - width ), max( 1, col - width )
            BR = min( h, row + width ), min( w, col + width )

            TL = ( row - len + 1, col - len + 1 )
            BR = ( row, col )

            subimg = image[ TL[1]:BR[1], TL[2]:BR[2] ]
            nh,_,_ = normhistogram( subimg, bindata )
            subOt  = otsu( nh )

            th = heavensideStep( R/abs(imgOt - subOt) - 1 )*150 + subOt

            thresholds[contr,contc] = th;
            contr += 1;
        end
        contc += 1;
    end

    ch = length( centers_r );
    cw = length( centers_c );

    # Interpolation
    dist = step + 1;

    # centers
    for idxc in 1:length(centers_c)
        col = centers_c[idxc]
        for idxr in 1:length(centers_r)
            row = centers_r[idxr]
            #imgc[row,col] = thresholds[idxr,idxc]
            imgc[row,col] = img[row,col] > thresholds[idxr,idxc]
        end
    end

    # vertically in line
    for idxc in 1:length(centers_c)
        col = centers_c[idxc];
        for idxr in 1:length(centers_r)-1
            U   = centers_r[ idxr ];
            B   = centers_r[idxr+1];
            thU = thresholds[ idxr , idxc ]
            thB = thresholds[idxr+1, idxc ]
            for row in U+1:B-1
                thI = ( B - row )/dist * thU + ( row - U )/dist * thB
                #imgc[ row, col ] = ( B - row )/dist * thU + ( row - U )/dist * thB
                imgc[ row, col ] = img[ row, col ] > thI
            end
        end
    end

    # laterally in line
    for idxc in 1:length(centers_c)-1
        L = centers_c[idxc]
        R = centers_c[idxc+1];
        for col in L+1:R-1
            fL = ( R - col )/dist;
            fR = ( col - L )/dist;
            for idxr in 1:length(centers_r)
                row = centers_r[idxr]
                thL = thresholds[idxr,  idxc ]
                thR = thresholds[idxr, idxc+1]
                #imgc[ row, col ] = fL * thL + fR * thR;
                imgc[ row, col ] = img[ row, col ] > fL * thL + fR * thR
            end
        end
    end

    # middle regions
    for idxc in 1:length(centers_c)-1
        L = centers_c[idxc]
        R = centers_c[idxc+1]
        for idxr in 1:length(centers_r)-1
            U = centers_r[idxr]
            B = centers_r[idxr+1]
            thUL = thresholds[  idxr ,  idxc  ]
            thUR = thresholds[  idxr , idxc+1 ]
            thBL = thresholds[ idxr+1,  idxc  ]
            thBR = thresholds[ idxr+1, idxc+1 ]

            for col in L+1:R-1
                fR  = ( col - L )/dist;
                fL  = ( R - col )/dist;

                for row in U+1:B-1
                    fU = ( B - row )/dist
                    fB = ( row - U )/dist

                    itpL  = fB * thBL + fU * thUL
                    itpR  = fB * thBR + fU * thUR
                    itpTH = fR * itpR + fL * itpL
                    # imgc[row,col] = itpTH
                    imgc[ row, col ] = img[ row, col ] > itpTH
                end
            end
        end
    end

    # first row
    fr = 1:centers_r[1]-1
    for col in 1:centers_c[1]
        #imgc[fr,col] .= thresholds[1,1]
        imgc[fr,col] .= img[fr,col] .> thresholds[1,1]
    end
    for col in centers_c[end]:w
        #imgc[fr,col] .= thresholds[1,end]
        imgc[fr,col] .= img[fr,col] .> thresholds[1,end]
    end
    for cidx in 1:length(centers_c)-1
        L   = centers_c[cidx]
        R   = centers_c[cidx+1]
        thL = thresholds[1,cidx]
        thR = thresholds[1,cidx+1]

        #imgc[fr,L] .= thL
        imgc[fr,L] .= img[fr,L] .> thL

        for col in L+1:R-1
            fR = ( col - L )/dist;
            fL = ( R - col )/dist;
            #imgc[fr,col] .= fR * thR + fL * thL
            imgc[fr,col] .= img[fr,col] .> ( fR * thR + fL * thL )
        end
    end

    # first col
    fc = 1:centers_c[1]-1
    for row in 1:centers_r[1]
        #imgc[row,fc] .= thresholds[1,1]
        imgc[row,fc] .= img[row,fc] .> thresholds[1,1]
    end
    for row in centers_r[end]:h
        #imgc[row,fc] .= thresholds[end,1]
        imgc[row,fc] .= img[row,fc] .> thresholds[end,1]
    end
    for ridx in 1:length(centers_r)-1
        U = centers_r[ridx]
        B = centers_r[ridx+1]

        thU = thresholds[ ridx , 1]
        thB = thresholds[ridx+1, 1]

        #imgc[U,fc] .= thU
        imgc[U,fc] .= img[U,fc] .> thU

        for row in U+1:B-1
            fB = ( row - U )/dist;
            fU = ( B - row )/dist;
            #imgc[row,fc] .= fU * thU + fB * thB
            imgc[row,fc] .= img[row,fc] .> fU * thU + fB * thB
        end
    end

    # last row
    lr = centers_c[end]:w
    for col in 1:centers_c[1]
        #imgc[lr,col] .= thresholds[end,1]
        imgc[lr,col] .= img[lr,col] .> thresholds[end,1]
    end
    for col in centers_c[end]:w
        #imgc[lr,col] .= thresholds[end,end]
        imgc[lr,col] .= img[lr,col] .> thresholds[end,end]
    end
    for cidx in 1:length(centers_c)-1
        L = centers_c[cidx]
        R = centers_c[cidx+1]

        thL = thresholds[end,cidx]
        thR = thresholds[end,cidx+1]

        #imgc[lr,L] .= thL
        imgc[lr,L] .= img[lr,L] .> thL

        for col in L+1:R-1
            fR = ( col - L )/dist;
            fL = ( R - col )/dist;
            #imgc[ lr,col] .= fR * thR + fL * thL
            imgc[ lr, col ] .= img[ lr, col ] .> ( fR * thR + fL * thL )
        end
    end

    # last col
    lc = centers_r[end]:h

    for row in 1:centers_r[1]
        #imgc[row,lc] .= thresholds[1,end]
        imgc[row,lc] .= img[row,lc] .> thresholds[1,end]
    end
    for row in centers_r[end]:h
        #imgc[row,lc] .= thresholds[end,end]
        imgc[row,lc] .= img[row,lc] .> thresholds[end,end]
    end
    for ridx in 1:length(centers_r)-1
        U = centers_r[ridx]
        B = centers_r[ridx+1]

        thU = thresholds[ ridx , end]
        thB = thresholds[ridx+1, end]

        #imgc[U,lc] .= thU
        imgc[U,lc] .= img[U,lc] .> thU

        for row in U+1:B-1
            fB = ( row - U )/dist;
            fU = ( B - row )/dist;
            #imgc[row,lc] .= fU * thU + fB * thB
            imgc[row,lc] .= img[row,lc] .> ( fU * thU + fB * thB )
        end
    end

    return imgc, thresholds
end


text:
struct LineIterator
    buf::Vector{UInt8}
    filestart::Int
    fileend::Int
    structtype::Union{Nothing, DataType}
    LineIterator(buf::Vector{UInt8},
        filestart::Int = 1,
        structtype::Union{Nothing, DataType} = nothing,
        fileend::Int = lastindex(buf)) =
        new(buf, filestart, fileend, structtype)
end

"""
    LineIterator(path::String; filestart = 1, structtype = nothing)

Create an iterator of a JSONLines file at `path`.

* Keyword Arguments:
    * `filestart=1`: Row at which to start the iterator
    * `structtype=nothing`: StructType passed to `JSON3.read` for each row
"""
LineIterator(path::String; filestart = 1, structtype = nothing) = LineIterator(Mmap.mmap(path), filestart, structtype)

## Iteration interface
function Base.iterate(lines::LineIterator, state::Int = lines.filestart; parse::Bool = true)
    state == lines.fileend && (return nothing)
    eol = _findnexteol(lines.buf, lines.fileend, state+1)
    ret = @inbounds(@view(lines.buf[state:eol]))
    return (parse ? parserow(ret, lines.structtype) : ret, eol)
end

function Base.eltype(lines::LineIterator)
    lines.structtype !== nothing && (return lines.structtype)
    return typeof(lines[1])
end

Base.IteratorSize(::Type{LineIterator}) = Base.SizeUnknown()

## Indexing interface
function _getindex(lines::LineIterator, i::Int)
    @warn "Indexin LineIterators is slow. Consider using LineIndex instead." maxlog=1
    linestart = skiprows(lines.buf, lines.fileend, i-1, lines.filestart) + (i>1)
    linestart >= lines.fileend && (throw(BoundsError(lines, i)))
    lineend = _findnexteol(lines.buf, lines.fileend, linestart)
    return linestart:lineend
end

function Base.getindex(lines::LineIterator, i::Int; parsed::Bool = true)
    index = _getindex(lines, i)
    ret = @inbounds(@view(lines.buf[index]))
    return parsed ? parserow(ret, lines.structtype) : ret
end

Base.firstindex(lines::LineIterator) = lines.filestart

function Base.getindex(lines::LineIterator, i::Int, col::Symbol)
    return lines[i][col]
end

function Base.getindex(lines::LineIterator, I::Vararg{Int, 2})
    row = lines[I[1]]
    col = collect(keys(row))[I[2]]
    return row[col]
end

## Tables.jl interface

Tables.istable(lines::LineIterator) = true
Tables.rowaccess(lines::LineIterator) = true
Tables.rows(lines::LineIterator) = lines
#Tables.materializer(lines::LineIterator) = Tables.rowtable


text:
function fit(
    df::DataFrame;
    max_time::Int = 10,
    fps = 0.022,
    loc_error = 0.03,
    p0::AbstractVector = [0.5, 0.5],
)
    @argcheck loc_error >= 0.0

    @. model(x, p) = 4 * p[1] * x^p[2] + 4 * loc_error^2
    result = curve_fit(
        model,
        df.delta_t[1:max_time] .* fps,
        df.msd[1:max_time],
        p0,
        lower = [0.0, 0.0],
        upper = [10.0, 2.0],
    )
    result
end

function estimate_α(df::DataFrame, δ1::Int, δ2::Int)
    @argcheck 1 <= δ2 <= δ1
    log(mean(df[df.Δt.==δ1, :msd]) / mean(df[df.Δt.==δ2, :msd])) / log(δ1 / δ2)
end

function fit_msd(
    df;
    max_time::Int64 = 10,
    loc_error::Float64 = 0.03,
    p0::AbstractVector = [0.5, 0.5],
)
    @argcheck loc_error >= 0.0
    @. model(x, p) = 4 * p[1] * x^p[2] + 4 * loc_error^2
    fit = curve_fit(
        model,
        df.delta_t[1:max_time] .* (1 / 45),
        df.msd[1:max_time],
        p0,
        lower = [0.0, 0.0],
        upper = [10.0, 2.0],
    )
    fit.param
end

function fit_msd(df::DataFrame, δ1::Int, δ2::Int)
    @argcheck 1 <= δ2 <= δ1
    log(mean(df[df.delta_t.==δ1, :msd]) / mean(df[df.delta_t.==δ2, :msd])) / log(δ1 / δ2)
end

fit_msd(df::DataFrame, δ::Int) = fit_msd(df::DataFrame, δ, 1)

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
using ProximalAlgorithms

x = Variable(10)
A = randn(5, 10)
y = Variable(7)
B = randn(5, 7)
b = randn(5)

println("\nTesting @minimize \n")
~x .= 0.
~y .= 0.
slv, = @minimize ls(A*x - B*y + b) st norm(x, 2) <= 1e4, norm(y, 1) <= 1.0 with ForwardBackward()
~x .= 0.
slv, = @minimize ls(A*x - b) st norm(x, 1) <= 1.0 with ForwardBackward()
~x .= 0.
slv, = @minimize ls(A*x - b) st norm(x, 1) <= 1.0
~x .= 0.
slv, = @minimize ls(A*x - b) + norm(x, 1) with ForwardBackward()
~x .= 0.
slv, = @minimize ls(A*x - b) + norm(x, 1)
~x .= 0.
slv, = @minimize ls(A*x - b)

#TODO many many more tests
x = Variable(5)
A = randn(10, 5)
b = randn(10)

println("\nTesting @minimize nonlinear \n")
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with ForwardBackward(tol = 1e-6)
xpg = copy(~x)
~x .= 0.
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with ZeroFPR(tol = 1e-6)
xz = copy(~x)
~x .= 0.
slv, = @minimize ls(sigmoid(A*x,10) - b)+norm(x,1) with PANOC(tol = 1e-6)
xp = copy(~x)
~x .= 0.

@test norm(xz-xpg) <= 1e-4
@test norm(xp-xpg) <= 1e-4

# test nonconvex Rosenbrock function with known minimum
solvers = [ZeroFPR(tol = 1e-6), PANOC(tol = 1e-6)]
for solver in solvers
    x = Variable(1)
    y = Variable(1)
    a,b = 2.0, 100.0

    cf = norm(x-a)^2+b*norm(pow(x,2)-y)^2
    @minimize cf+1e-10*norm(x,1)+1e-10*norm(y,1) with solver

    @test norm(~x-[a]) < 1e-4
    @test norm(~y-[a^2]) < 1e-4
end


text:
using Core.Compiler: NativeInterpreter

function process(@nospecialize(f), @nospecialize(TT=()); optimize=true)
    (interp, mi) = Cthulhu.mkinterp(NativeInterpreter(), f, TT)
    (ci, rt, infos, slottypes) = Cthulhu.lookup(interp, mi, optimize; allow_no_codeinf=true)
    if ci !== nothing
        ci = Cthulhu.preprocess_ci!(ci, mi, optimize, Cthulhu.CthulhuConfig(dead_code_elimination=true))
    end
    interp, ci, infos, mi, rt, slottypes
end

function find_callsites_by_ftt(@nospecialize(f), @nospecialize(TT=Tuple{}); optimize=true)
    interp, ci, infos, mi, _, slottypes = process(f, TT; optimize)
    ci === nothing && return Cthulhu.Callsite[]
    callsites = Cthulhu.find_callsites(interp, ci, infos, mi, slottypes, optimize)
    return callsites
end

macro find_callsites_by_ftt(ex0...)
    return InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :find_callsites_by_ftt, ex0)
end

function strip_base_linenums(s)
    s = replace(s, r"((boot|int|refvalue|refpointer)\.jl:)\d+" => s"\1")
    return s
end

function irshow_filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations)
    return "test_output/foo-$(optimize ? :opt : :unopt)-$(debuginfo)-$(iswarn ? :warn : :nowarn)\
        -$(hide_type_stable ? :hide_type_stable : :show_type_stable)-$(inline_cost ? :inline_cost : :nocost)\
        -$(type_annotations ? :types : :notypes)"
end

# to generate test cases for IRShow tests
function generate_test_cases(foo)
    outputs = Dict()
    tf = (true, false)
    for optimize in tf
        _, src, infos, mi, rt, slottypes = process(foo, (Int, Int); optimize);
        for (debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations) in Iterators.product(
            instances(Cthulhu.DInfo.DebugInfo), tf, tf, tf, tf,
        )
            !optimize && debuginfo === Cthulhu.DInfo.compact && continue
            !optimize && inline_cost && continue

            s = sprint(; context=:color=>true) do io
                Cthulhu.cthulhu_typed(io, debuginfo,
                                      src, rt, mi;
                                      iswarn, hide_type_stable, inline_cost, type_annotations)
            end
            s = strip_base_linenums(s)
            write(irshow_filename(optimize, debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations), s)

            in(s, values(outputs)) && (@show((optimize, debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations)); @show(findfirst(==(s), pairs(outputs))))
            @assert !in(s, values(outputs))
            outputs[(optimize, debuginfo, iswarn, hide_type_stable, inline_cost, type_annotations)] = s
        end
    end
end


text:
#Perform mulivariate regression on ilr transformed compositions

#Kristen R. Hunter-Cevera, Marine Biological Laboratory, Spring 2020

using JLD2, DataFrames, CSV, Dates, Query
#using StatsPlots, Colors
using LinearAlgebra, Clustering, GLM
using Gadfly, Colors, Compose
using RCall

using Distributions

include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

## Need to have oligo_df_plus loaded! This has zeros replaced with small amount:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"

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

#F and Chi-squared approximations:
function approx_F(Λ,p,vE,vH)
    if (p * vH) == 2
        t = 1
    else
        t = sqrt((p^2 * vH^2 - 4) / (p^2 + vH^2 -5))
    end
    w = vE + vH - 0.5 * (p + vH +1)
    df1 = p * vH
    df2 = w * t - 0.5 * (p * vH -2)

    F = ((1 - Λ.^(1/t)) / (Λ.^(1/t))) .* (df2 / df1)
    return F, df1, df2
end

function approx_Chi(Λ,p,vE,vH)
    df = p * vH
    chi = -(vE - 0.5 *(p - vH +1)) * log(Λ)
    return chi, df
end

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

#transform compositions to multivariate real numbers with ilr:
y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)]))
sbp0=[1 -1 1 1 1 1; 1 0 1 -1 -1 1; 0 0 0 1 -1 0; 1 0 1 0 0 -1; 1 0 -1 0 0 0]
lb0=["O1,O3,O4,O5,O6 | O2", "O1,O3,O6 | O4,O5", "O4 | O5", "O1,03 | O6", "O1 | O3"] #labels

ilr_y, bal=ilr(y,sbp0) #ilr transform
lb=lb0 #set labels

## Full model regression test - manually choose the variable combination:

vars=[:temp,:lightweek,:PO4,:NH4,:SiOH,:NO3]

#break up year into seasons:
wisp=findall(dayofyear.(oligo_df_plus[:,:date]) .< 166)
su=findall((dayofyear.(oligo_df_plus[:,:date]) .>= 166) .& (dayofyear.(oligo_df_plus[:,:date]) .< 258))
fa=findall(dayofyear.(oligo_df_plus[:,:date]) .>= 258)

X_wisp = convert(Array{Union{Float64,Missing}},oligo_df_plus[wisp,vars])
Y_wisp = ilr_y[wisp,:]
#jj=findall(ismissing.(X_wisp)) #-> this should be zero

X_su = convert(Array{Union{Float64,Missing}},oligo_df_plus[su,vars])
Y_su = ilr_y[su,:]
jj=findall(ismissing.(X_su))
qq=unique(map(x->x[1], jj)) #no nutrient values for one of the samples, exclude
ii=setdiff(collect(1:size(X_su,1)),qq)
X_su=X_su[ii,:]
Y_su=Y_su[ii,:]

X_fa = convert(Array{Union{Float64,Missing}},oligo_df_plus[fa,vars])
Y_fa = ilr_y[fa,:]
#jj=findall(ismissing.(X_fa)) #-> this also should be zero

## Step wise progression to test significance of added variables, one at a time:

#find the minimum number of variables that explains the most variability for each season:

#record in a dataframe
stepwise_df=DataFrame(season=String[], step=Int64[], base=String[], var=String[], Lambda=Float64[], F=Float64[], df1=Float64[], df2=Float64[], F_pval=Float64[], B_hat=Array{Union{Missing, Float64},2}[])

# first, find out which variable results in lowest lambda value:
seasons=["winter/spring", "summer", "fall"]

for j=1:3

    global season=seasons[j]
    #season = "winter/spring"
    if season == "winter/spring"
        X0=X_wisp
        Y=Y_wisp
    elseif season == "summer"
        X0=X_su
        Y=Y_su
    elseif season == "fall"
        X0=X_fa
        Y=Y_fa
    end

    intercept=ones(size(X0,1),1) #setup design matrix
    n=size(X0,1)
    y_bar=mean(Y, dims =1)

    @show season

    for q=1:6 # number of predictor variables to run through

        #@show q
        X=hcat(intercept,X0[:,q]) #build design matrix with each variable

        B_hat= inv(X' * X) * (X' * Y) #retrieve parameters

        H_full = (Y' * Y) .- (B_hat' * X' * Y) #calculate lambda
        H_null = Y' * Y .- (n .* (y_bar' * y_bar)) #E+H in book parlance
        Λf = det(H_full) ./ det(H_null)

        vH = length(q) #number of X vars, intercept excluded
        vE = n-(length(q))-1 #n - x vars -1
        p = size(Y,2) #num y vars

        f, df1 ,df2=approx_F(Λf,p,vE,vH)
        pvalF=1-cdf(FDist(df1,df2),f)

        push!(stepwise_df,[season,0,string(vars[q]), "-", Λf, f, df1, df2, pvalF, B_hat])

end

#once we have a baseline, find which var explains the most variance and add from there!
r=findall(stepwise_df[!,:season] .== season)
global qbase=argmin(stepwise_df[r,:Lambda])
global qtest=setdiff(1:6,qbase)
global stepp=1

while stepp <= 6
    global qtest
    global qbase
    global stepp

    for q in qtest

        # @show q
        # @show qbase
        X=hcat(intercept,X0[:,cat(qbase,q,dims=1)]) #full model with variable q as candidate variable
        B_hat= inv(X' * X) * (X' * Y) #retrieve parameters
        H_full = (Y' * Y) .- (B_hat' * X' * Y) #calculate lambda
        H_null = Y' * Y .- (n .* (y_bar' * y_bar)) #E+H in book parlance
        Λf = det(H_full) ./ det(H_null)

        Xr=hcat(intercept,X0[:,qbase]) #Reduced model
        B_hat_r = inv(Xr' * Xr) * (Xr' * Y)
        H_r = (Y' * Y) .- (B_hat_r' * Xr' * Y)
        Λr = det(H_r) ./ det(H_null) #overall test of signifance, 0.86-0.89

        Λ_test = Λf ./ Λr #Lamba to compare addition of extra variable

        # approximate F and χ statistics:
        vH = length(q) #additional varsglobal step=1
        vE = n-(length(qbase))-1 #n - x vars in reduced model -1
        p = size(Y,2) #num y vars

        f, df1, df2 = approx_F(Λ_test,p,vE,vH)
        #chi, df = approx_Chi(Λ_test,p,vE,vH)

        #p-values:
        pvalF=1-cdf(FDist(df1,df2),f)
        #pvalchi = 1-cdf(Chisq(df),chi)
        if length(qbase) == 1
            push!(stepwise_df,[season, stepp, string(vars[qbase]), string(vars[q]), Λ_test, f, df1, df2, pvalF, B_hat])
        else
            push!(stepwise_df,[season, stepp, join(string.(vars[qbase]),"-"), string(vars[q]), Λ_test,f, df1, df2, pvalF, B_hat])
        end

    end #end qtest loop

    r=findall((stepwise_df[!,:step] .== stepp) .& (stepwise_df[!,:season] .== season))

    if any(stepwise_df[r,:F_pval] .< 0.05) #alright- another variable does add significant explanation!

        qbase=vcat(qbase,qtest[argmin(stepwise_df[r,:Lambda])])
        qtest=setdiff(1:6,qbase)
        stepp=stepp .+ 1

    else
        break
    end

end

end #seaons
## Sweet!  Now, backtransform for plotting if needed:

#if needed (see below for saving)
@load  "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_ilr_regression_output2.jld2"

##Parameter values back to simplex:

#sumwinter/spring: temp-lightweek
rr=findall(stepwise_df[!,:season] .== "winter/spring")
final_step=maximum(stepwise_df[rr,:step])-1
rr2=findall((stepwise_df[!,:season] .== "winter/spring") .& (stepwise_df[!,:step] .== final_step))
best_model_wisp=rr2[argmin(stepwise_df[rr2,:Lambda])]

bhat_wisp_simplex=invilr(stepwise_df[best_model_wisp,:B_hat],bal)

#summer: temp-PO4-lightweek
rr=findall(stepwise_df[!,:season] .== "summer")
final_step=maximum(stepwise_df[rr,:step])-1
rr2=findall((stepwise_df[!,:season] .== "summer") .& (stepwise_df[!,:step] .== final_step))
best_model_su=rr2[argmin(stepwise_df[rr2,:Lambda])]

bhat_summer_simplex=invilr(stepwise_df[best_model_su,:B_hat],bal)

#fall: temp-lightweek
rr=findall(stepwise_df[!,:season] .== "fall")
final_step=maximum(stepwise_df[rr,:step])-1
rr2=findall((stepwise_df[!,:season] .== "fall") .& (stepwise_df[!,:step] .== final_step))
best_model_fa=rr2[argmin(stepwise_df[rr2,:Lambda])]

bhat_fall_simplex=invilr(stepwise_df[best_model_fa,:B_hat],bal) #bhat_wisp = Bhat

## Fitted data back to simplex:

Y_hat = Array{Union{Missing, Float64}}(missing, size(ilr_y))

#populate by seasons:
x =hcat(ones(size(wisp,1),1), X_wisp[:,[1,2]])
b = stepwise_df[best_model_wisp,:B_hat]
temp = b' * x' #B_hat is full model
Y_hat[wisp,:] .= temp'

x =hcat(ones(size(ii,1),1), X_su[:,[1,3,2]]) #initial inexing into summer as one sample does not have all the nutrients!
b = stepwise_df[best_model_su,:B_hat]
temp = b' * x' #B_hat is full model
Y_hat[su[ii],:] .= temp'

x =hcat(ones(size(fa,1),1), X_fa[:,[1,2]])
b = stepwise_df[best_model_fa,:B_hat]
temp = b' * x' #B_hat is full model
Y_hat[fa,:] .= temp'

#and back transform to the simplex:
yhat_simplex = invilr(Y_hat,bal)

#Only temperature as parameter:
Y_hat_temp = Array{Union{Missing, Float64}}(missing, size(ilr_y))

#populate by seasons:
x =hcat(ones(size(wisp,1),1), X_wisp[:,[1]])
b = stepwise_df[1,:B_hat] #temperature only model
temp = b' * x' #B_hat is full model
Y_hat_temp[wisp,:] .= temp'

x =hcat(ones(size(ii,1),1), X_su[:,1]) #initial inexing into summer as one sample does not have all the nutrients!
b = stepwise_df[16,:B_hat]
temp = b' * x' #B_hat is full model
Y_hat_temp[su[ii],:] .= temp'

x =hcat(ones(size(fa,1),1), X_fa[:,1])
b = stepwise_df[34,:B_hat]
temp = b' * x' #B_hat is full model
Y_hat_temp[fa,:] .= temp'

#and back transform to the simplex:
yhat_simplex = invilr(Y_hat,bal)
yhat_simplex_temp = invilr(Y_hat_temp,bal)
## and save!

@save "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_ilr_regression_output2.jld2" stepwise_df fa wisp su ii X_wisp X_fa X_su Y_hat yhat_simplex bhat_fall_simplex bhat_summer_simplex bhat_wisp_simplex

## if need to add to oligo_df_zero_plus dataframe:

# oligo_df_plus.ilr1 = ilr_y[:,1]
# oligo_df_plus.ilr2 = ilr_y[:,2]
# oligo_df_plus.ilr3 = ilr_y[:,3]
# oligo_df_plus.ilr4 = ilr_y[:,4]
# oligo_df_plus.ilr5 = ilr_y[:,5]

#oligo_df_plus.yearday = dayofyear.(oligo_df_plus[!,:date])
#CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/rel_abnd_analysis_scripts/matlab_scripts/olgio_df_plus_with_ilr.csv", oligo_df_plus, delim = ',')

b=convert(DataFrame,[X_su y[su[ii],:]]) #Y_fa
bnames=names(b)
for j=1:6
    rename!(b,bnames[j] => vars[j])
end

# rename!(b,:x7 => :ilr1)
# rename!(b,:x8 => :ilr2)
# rename!(b,:x9 => :ilr3)
# rename!(b,:x10 => :ilr4)
# rename!(b,:x11 => :ilr5)

rename!(b,:x7 => :O1)
rename!(b,:x8 => :O2)
rename!(b,:x9 => :O3)
rename!(b,:x10 => :O4)
rename!(b,:x11 => :O5)
rename!(b,:x12 => :O6)

CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/R_scripts/su_data.csv", b, delim = ',')

# a=DataFrame(yhat_1=yhat_simplex[:,1],
# yhat_2=yhat_simplex[:,2],
# yhat_3=yhat_simplex[:,3],
# yhat_4=yhat_simplex[:,4],
# yhat_5=yhat_simplex[:,5],
# yhat_6=yhat_simplex[:,6],
# yhat_temp1=yhat_simplex_temp[:,1],
# yhat_temp2=yhat_simplex_temp[:,2],
# yhat_temp3=yhat_simplex_temp[:,3],
# yhat_temp4=yhat_simplex_temp[:,4],
# yhat_temp5=yhat_simplex_temp[:,5],
# yhat_temp6=yhat_simplex_temp[:,6],
# yhat_ilr1=Y_hat[:,1],
# yhat_ilr2=Y_hat[:,2],
# yhat_ilr3=Y_hat[:,3],
# yhat_ilr4=Y_hat[:,4],
# yhat_ilr5=Y_hat[:,5])
#
# CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/rel_abnd_analysis_scripts/matlab_scripts/yhat_simplex.csv",a, delim = ',')

## If want to just models individually:

X0=X_fa, #X_wisp, X_su
Y=Y_fa # Y_su, Y_wisp
n=size(X0,1)
y_bar=mean(Y, dims =1)

intercept=ones(size(X0,1),1) #setup design matrix
X=hcat(intercept,X0[:,1]) #vars here

B_hat= inv(X' * X) * (X' * Y) #retrieve parameters

H_full = (Y' * Y) .- (B_hat' * X' * Y) #calculate lambda
H_null = Y' * Y .- (n .* (y_bar' * y_bar)) #E+H in book parlance
Λf = det(H_full) ./ det(H_null)

# Λ is distributed as q, p,n− p−1 as well as p,q,n−q−1, where:
# p is number of response variables
# q is the number of independent predictor variables
# n is the number of observations

# If Λf < Λ(p,q,n-q-1), then model is significant

#Can also estimate p-values with F or Chi-squared approximation:
## for full model:

vH = size(X,2)-1 #number of X vars, intercept excluded
vE = n-(size(X,2)-1)-1 #n - x vars -1
p = size(Y,2) #num y vars

f, df1 ,df2=approx_F(Λf,p,vE,vH)
pvalF=1-cdf(FDist(df1,df2),f)

#chi, df = approx_Chi(Λ_test,p,vE,vH)
#pvalchi = 1-cdf(Chisq(df),chi)

## parameter fits and projected model:

#If model is significant, see model estimates: populate once have run through each season:
seas_yhat=zeros(size(ilr_y))

Y_hat = B_hat' * X' #B_hat is full model
seas_yhat[ii[seas],:] = Y_hat'

# can also back transfrom fitted estimates to simplex:
# comp_yhat=invilr(seas_yhat,bal)
# comp_yhat_temp=invilr(seas_yhat_temp,bal)

#Can also back transform parameter fits into the simplex:

comp_bhat_wisp=invilr(bhat_wisp,bal) #bhat_wisp = Bhat

## correlations with phosphate:

using Statistics
using StatsBase

jj=findall(ismissing.(oligo_df_plus[!,:PO4]))
 #no nutrient values for one of the samples, exclude
ii=setdiff(collect(1:size(oligo_df_plus,1)),jj)

cor(oligo_df_plus[ii,:PO4],ilr_y[ii,5])
corspearman(convert(Array{Float64},oligo_df_plus[ii,:PO4]),ilr_y[ii,5])
corspearman(convert(Array{Float64},oligo_df_plus[ii,:PO4]),ilr_y[ii,4])

plot(x=oligo_df_plus[ii,:PO4],y=ilr_y[ii,4], Geom.point)
plot(x=oligo_df_plus[ii,:PO4],y=ilr_y[ii,5], Geom.point)

## plots for Wilk's Lambda...

#Sanity check just with interpolation....
#for:
#p=5 #number of response variables
#h=1 #number of additional variables compared to reduced - 1 in our case if want to test significance of added variable

#vE = n - q- 1 #number of samples, number of vars in full model -1

#testing for light signficance:
wispE = length(wisp) - 2 - 1
suE = length(ii)- 3 -1 #had two significant vars
faE = length(fa)- 2 -2

vE=[10:30; 40; 60; 80; 100]
cval = [0.215; 0.261; 0.303; 0.341; 0.376; 0.407; 0.436; 0.462; 0.486; 0.508; 0.529; 0.548; 0.565; 0.581; 0.596; 0.610; 0.623; 0.635; 0.647; 0.658;0.668; 0.744; 0.825; 0.867; 0.893]

l1 = layer(x=vE,y=cval,Geom.point)
l2 = layer(x=vE,y=cval,Geom.line)
l3= layer(x=[wispE], y=[stepwise_df[7,:Lambda]], Geom.point, style(default_color=RGB(0,0,0),point_size=2mm),shape=[Shape.square]) #,style(default_color=RGB(0,0,0), point_shapes=[square]))
p1=plot(l1,l2,l3)

l1 = layer(x=vE,y=cval,Geom.point)
l2 = layer(x=vE,y=cval,Geom.line)
l3= layer(x=[suE], y=[stepwise_df[27,:Lambda]], Geom.point, style(default_color=RGB(0,0,0),point_size=2mm),shape=[Shape.square]) #,style(default_color=RGB(0,0,0), point_shapes=[square]))
p2=plot(l1,l2,l3)

l1 = layer(x=vE,y=cval,Geom.point)
l2 = layer(x=vE,y=cval,Geom.line)
l3= layer(x=[faE], y=[stepwise_df[40,:Lambda]], Geom.point, style(default_color=RGB(0,0,0),point_size=2mm),shape=[Shape.square]) #,style(default_color=RGB(0,0,0), point_shapes=[square]))
p3=plot(l1,l2,l3)

Gadfly.set_default_plot_size(30cm,10cm)
gridstack([p1 p2 p3])


text:
module MakeFshared

export printout, calc, foo, moo

import Libdl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("shared library not installed properly; run Pkg.build(\"MakeFshared\"), restart Julia, and try again.")
end

include(depsjl_path)

function __init__()
    check_deps()
end

function printout()
    ccall((:display_, libextraf), Nothing, ())
end

function calc(x, y)
    ccall((:__eins_MOD_calc, libextraf), Int32,
        (Ptr{Int32}, Ptr{Int32},), x, y)
end

function foo(x, y)
    ccall((:__zwei_MOD_foo, libextraf), Nothing,
        (Ptr{Int32}, Ptr{Int32},), x, y)
end

function moo(x, y)
    ccall((:__zwei_MOD_moo, libextraf), Nothing,
        (Ptr{Int32}, Ptr{Int32},), x, y)
end

end # module