
using Statistics

function format(x)
    if abs(x) < 999
        return @sprintf("%12.7f", x)
    else
        return @sprintf("%12.5e", x)
    end
end

import Base.write
"""

```
write(R::ComplexMixtures.Result, filename::String, solute::Selection, solvent::Selection)
```

Function to write the final results to output files as simple tables
that are human-readable and easy to analyze with other software

If the solute and solvent selections are provides, pass on the atom
names.

"""
write(R::Result, filename::String, solute::Selection, solvent::Selection) =
    write(R, filename, solute_names = solute.names, solvent_names = solvent.names)


"""

```
write(R::ComplexMixtures.Result, filename::String; 
      solute_names::Vector{String} = ["nothing"], 
      solvent_names::Vector{String} = ["nothing"])
```

Optional passing of atom names.

"""
function write(
    R::Result,
    filename::String;
    solute_names::Vector{String} = ["nothing"],
    solvent_names::Vector{String} = ["nothing"],
)

    # Names of output files containing atomic contibutions
    atom_contrib_solvent =
        FileOperations.remove_extension(filename) *
        "-ATOM_CONTRIB_SOLVENT." *
        FileOperations.file_extension(filename)
    atom_contrib_solute =
        FileOperations.remove_extension(filename) *
        "-ATOM_CONTRIB_SOLUTE." *
        FileOperations.file_extension(filename)


    #
    # GMD computed with minimum distance
    #

    # Conversion factor for volumes (as KB integrals), from A^3 to cm^3/mol
    mole = 6.022140857e23
    convert = mole / 1.e24

    output = open(filename, "w")
    println(output, @sprintf("#"))
    println(output, @sprintf("# Output of ComplexMixtures - MDDF"))
    println(output, @sprintf("# Trajectory files and weights:"))
    for i = 1:length(R.files)
        println(output, "#  $(R.files[i]) - w = $(R.weights[i])")
    end
    println(output, @sprintf("#"))
    println(
        output,
        @sprintf(
            "# Density of solvent in simulation box (sites/A^3): %15.8f",
            R.density.solvent
        )
    )
    println(
        output,
        @sprintf(
            "# Density of solvent in bulk (estimated) (sites/A^3): %15.8f",
            R.density.solvent_bulk
        )
    )
    println(
        output,
        @sprintf(
            "# Molar volume of solvent in simulation (cc/mol): %15.8f",
            convert / R.density.solvent
        )
    )
    println(
        output,
        @sprintf(
            "# Molar volume of solvent in bulk (estimated) (cc/mol): %15.8f",
            convert / R.density.solvent_bulk
        )
    )
    println(output, @sprintf("#"))
    println(output, @sprintf("# Number of atoms solute: %i9", size(R.solute_atom, 2)))
    println(
        output,
        @sprintf("# Number of atoms of the solvent: %i9", size(R.solvent_atom, 2))
    )
    println(output, @sprintf("#"))

    if R.options.usecutoff
        ibulk = setbin(R.options.dbulk, R.options.binstep)
    else
        ibulk = round(Int, R.nbins - 1 / R.options.binstep)
    end
    bulkerror = Statistics.mean(R.mddf[ibulk:R.nbins])
    sdbulkerror = Statistics.std(R.mddf[ibulk:R.nbins])
    println(output, "#")
    println(
        output,
        @sprintf(
            "# Average and standard deviation of bulk-gmd: %12.5f +/- %12.5f",
            bulkerror,
            sdbulkerror
        )
    )
    println(output, "#")

    # Output table

    println(output, "# COLUMNS CORRESPOND TO: ")
    println(output, "#       1  Minimum distance to solute (dmin)")
    println(
        output,
        "#       2  GMD distribution (md count normalized by md count of random-solute distribution)",
    )
    println(
        output,
        "#       3  Kirwood-Buff integral (cc/mol) computed [(1/bulkdensity)*(col(6)-col(7))].",
    )
    println(output, "#       4  Minimum distance site count for each dmin.")
    println(
        output,
        "#       5  Minimum distance site count for each dmin for random solute distribution.",
    )
    println(
        output,
        "#       6  Cumulative number of molecules within dmin in the simulation.",
    )
    println(
        output,
        "#       7  Cumulative number of molecules within dmin for random solute distribution.",
    )
    println(output, "#       8  Volume of the shell of distance dmin and width binstep.")
    println(output, "#")
    println(
        output,
        "#   1-DISTANCE         2-GMD      3-KB INT    4-MD COUNT  5-COUNT RAND      6-SUM MD    7-SUM RAND   8-SHELL VOL",
    )

    for i = 1:R.nbins
        line = "  " * format(R.d[i])                                 #  1-DISTANCE
        line = line * "  " * format(R.mddf[i])                         #  2-GMD
        line = line * "  " * format(R.kb[i])                           #  3-KB INT
        line = line * "  " * format(R.md_count[i])                     #  4-MD COUNT
        line = line * "  " * format(R.md_count_random[i])              #  5-COUNT RAND
        line = line * "  " * format(R.sum_md_count[i])                 #  6-SUM MD
        line = line * "  " * format(R.sum_md_count_random[i])          #  7-SUM RAND
        line = line * "  " * format(R.volume.shell[i])                 #  8-SHELL VOL
        println(output, line)
    end
    close(output)

    # Writting gmd per atom contributions for the solvent

    output = open(atom_contrib_solvent, "w")
    println(output, "#")
    println(output, "# Solvent atomic contributions to total MDDF. ")
    println(output, "#")
    println(output, "# Trajectory files: ", R.files)
    println(output, "#")
    println(output, "# Atoms: ")
    for i = 1:size(R.solvent_atom, 2)
        if solvent_names[1] == "nothing"
            println(output, @sprintf("# %9i", i))
        else
            println(output, @sprintf("# %9i %5s", i, solvent_names[i]))
        end
    end
    println(output, "#")
    string = "#   DISTANCE     GMD TOTAL"
    for i = 1:size(R.solvent_atom, 2)
        string = string * @sprintf("  %12i", i)
    end
    println(output, string)
    for i = 1:R.nbins
        string = format(R.d[i])
        string = string * "  " * format(R.mddf[i])
        for j = 1:size(R.solvent_atom, 2)
            string = string * "  " * format(R.solvent_atom[i, j])
        end
        println(output, string)
    end
    close(output)

    # Writting gmd per atom contributions for the solute

    output = open(atom_contrib_solute, "w")
    println(output, "#")
    println(output, "# Solute atomic contributions to total MDDF. ")
    println(output, "#")
    println(output, "# Trajectory files: ", R.files)
    println(output, "#")
    println(output, "# Atoms:")
    for i = 1:size(R.solute_atom, 2)
        if solute_names[1] == "nothing"
            println(output, @sprintf("# %9i", i))
        else
            println(output, @sprintf("# %9i %5s", i, solute_names[i]))
        end
    end
    println(output, "#")
    string = "#   DISTANCE      GMD TOTAL"
    for i = 1:size(R.solute_atom, 2)
        string = string * @sprintf("  %12i", i)
    end
    println(output, string)
    for i = 1:R.nbins
        string = format(R.d[i]) * "  " * format(R.mddf[i])
        for j = 1:size(R.solute_atom, 2)
            string = string * "  " * format(R.solute_atom[i, j])
        end
        println(output, string)
    end
    close(output)

    # Write final messages with names of output files and their content

    println()
    println(" OUTPUT FILES: ")
    println()
    println(" Wrote solvent atomic GMD contributions to file: ", atom_contrib_solvent)
    println(" Wrote solute atomic GMD contributions to file: ", atom_contrib_solute)
    println()
    println(" Wrote main output file: ", filename)
    println()

    return nothing
end
