text:
module NPCFs
export NPCF, compute_npcf_simple, compute_npcf_pairwise, summarize, angular_indices

### Load packages
using Printf: @printf
using Distributed: @sync, @async, @spawnat, nworkers, workers, Future, myid
using Parameters: @with_kw
using Rotations: AngleAxis
using WignerSymbols: wigner3j, wigner6j
using SphericalHarmonics: computePlmcostheta
using LinearAlgebra: cross, norm
using LegendrePolynomials: Pl
using HypergeometricFunctions: _₂F₁
using ClassicalOrthogonalPolynomials: chebyshevt
using SpecialFunctions: gamma

### Define basic NPCF structure

"""
    NPCF(N, D, coords = "cartesian", periodic = "false", volume = 1e9, nbins = 10, r_min = 50, r_max = 200, lmax = 2, verb = false, _dr = -1, _boxsize = -1)

Create the NPCF object, containing all required attributes and output arrays.
Variables can be specified by their keywords, and take default values if not given.
The key variables are the NPCF order `N` and the dimension `D`.

# Arguments
- `N::Int64`: which NPCF to use (N = 2, 3, 4, 5 currently supported)
- `D::Int64`: number of dimensions

# Optional Arguments
- `coords::String = "cartesian"`: whether input coordinates are "cartesian" or "spherical"
- `periodic::Bool = false`: whether the particles use a periodic wrapping
- `volume::Float64 = 1e9`: volume of input (used for normalization, and periodic wrapping)

- `nbins::Int64 = 10`: number of radial bins
- `r_min::Float64 = 50.`: minimum bin
- `r_max::Float64 = 200`: # maximum bin

- `lmax::Int64 = 2`: maximum ell
- `verb::Bool = false`: whether to print throughout runtime

- `_dr::Float64 = -1:` bin-size, this will be reset on initialization
- `_boxsize::Float64 = -1:` periodic boxsize, this will be reset on initialization

# Examples
```julia-repl
julia> npcf = NPCF(N=4, D=3); # create a 4PCF object in 3D
...
julia> summarize(npcf); # summarize the NPCF inputs
...
julia> inputs = hcat(rand(100,3)*1000.,ones(100)); # generate 100 randomly distributed particles in [0, 1000], plus weights
...
julia> output_simple = compute_npcf_simple(npcf, inputs); # run the simple NPCF estimator
julia> output_pairwise = compute_npcf_pairwise(npcf, inputs); # run the efficient NPCF estimator
...
julia> l1, l2, l3  = angular_indices(npcf); # return a list of ell indices corresponding to the one-dimensional indexing used in the code.
```
"""
@with_kw struct NPCF
    # NPCF Structure
    N::Int64 # which NPCF to use (N = 2, 3, 4, 5 currently supported)
    D::Int64 # number of dimensions

    # Coordinates
    coords::String = "cartesian" # whether input coordinates are "cartesian" or "spherical"

    # Periodicity
    periodic::Bool = false # whether the particles use a periodic wrapping
    volume::Float64 = 1e9 # volume of input (used for normalization, and periodic wrapping)

    # Radial Binning
    nbins::Int64 = 10 # number of radial bins
    r_min::Float64 = 50. # minimum bin
    r_max::Float64 = 200. # maximum bin

    # Angular Binning
    lmax::Int64 = 2; # maximum ell

    # Verbosity
    verb::Bool = false # whether to print throughout runtime

    # Other parameters, set on initialization
    _dr::Float64 = -1; # bin-size, computed from binning parameters
    _boxsize::Float64 = -1; # boxsize, computed from volume

    # Check bounds and create object
    function NPCF(N,D,coords,periodic,volume,nbins,r_min,r_max,lmax,verb,_dr,_boxsize)

        # Check that derived parameters haven't been set
        @assert _boxsize==-1 "This will be reset on initialization!"
        @assert _dr==-1 "This will be reset on initialization!"

        # Set additional parameters
        _dr = (r_max-r_min)/nbins;
        if periodic
            _boxsize = volume^(1/D)
        end

        # Check input parameters
        @assert (D==4 || D==3 || D==2) "Only 2-, 3- and 4-dimensions implemented"
        @assert (N==5 || N==4 || N==3 || N==2) "Only 2PCF, 3PCF, 4PCF and 5PCF implemented"
        @assert r_max>r_min "Maximum radius must be greater than minimum radius"
        if coords=="spherical"
            @assert D==2 "Spherical coordinates only implemented in 2-dimensions"
            @assert periodic==false "Periodic wrapping not available for spherical coordinates"
            @assert (r_min>=-1 && r_max<=1) "Radial bins must be in range [-1,1]"
            if volume!=4pi
                println("Resetting volume to 4pi")
                volume=4pi
            end
        elseif coords=="cartesian"
            @assert (r_min>=0 && r_max>=0) "Radial bins must be in range [0,inf)"
            if periodic
                @assert (r_max<_boxsize/2) "r-max must be smaller than boxsize/2"
                @assert _boxsize>0 "Boxsize must be greater than zero"
            end
        end

        npcf=new(N,D,coords,periodic,volume,nbins,r_min,r_max,lmax,verb,_dr,_boxsize)
        if verb; summarize(npcf); end
        return npcf
    end
end

### Load utility functions

include("utils.jl")

### Key NPCF functions

"""
    compute_npcf_simple(inputs, npcf)

Compute the NPCF by summing over `N`-tuplets of particles given in the `inputs` dataset.
For `n` particles, the computational cost scales as `n^N`, and is thus not recommended for large `n`.
Work is distributed approximately equally over all ```nworkers()``` workers provided.

The function is initialized with the `NPCF` class, which specifies all relevant parameters.
The `inputs` array gives a list of particle positions in `D`-dimensions, and a `D+1`-th column specifying the particle weight.
The output is a vector array, where the first `N-1` indices specify the radial bins, and the `N`-th index gives the angular bin, collapsed into one dimension (see ```angular_indices``` for conversion to ells).
Note that we only fill up radial bins satisfying bin1<bin2<...<binN.
The actual computation is performed in the ```_compute_npcf_simple()``` function; this is just a wrapper function to facilitate multiprocessing.

# Examples
```julia-repl
julia> inputs = hcat(rand(100,3)*1000.,ones(100)); # 100 randomly distributed 3D particles in [0, 1000], plus weights
julia> npcf = NPCF(N=4, D=3); # create a 4PCF object in 3D
julia> output_simple = compute_npcf_simple(inputs, npcf); # run the simple NPCF estimator
```
"""
function compute_npcf_simple(inputs::Matrix{Float64}, npcf::NPCF)

    # First check inputs
    check_inputs(inputs, npcf)

    Npart = length(inputs[:,1])
    if npcf.verb; @printf("\nAnalyzing dataset with %d particles using %d worker(s) by summing over all N-tuplets of galaxies\n",Npart,nworkers()); end

    # Set up output array to hold NPCF
    output = create_npcf_array(npcf)

    # Work out how many workers are present and chunk memory
    tasks_per_worker = div(Npart,nworkers())+1

    # Now perform the chunking and run the algorithm across separate workers
    if npcf.verb;
        if nworkers()>1
            println("Starting distributed computation");
        else
            println("Starting computation");
        end
    end
    chunked_output = Dict{Int,Future}()
    @sync for (w_i,w) in enumerate(workers())
        imin = tasks_per_worker*(w_i-1)+1
        imax = min(imin+tasks_per_worker,Npart+1)-1
        @async chunked_output[w_i] = @spawnat w _compute_npcf_simple(imin, imax, inputs, npcf)
    end

    # Now combine the output together
    for w_i in 1:nworkers()
        output .+= chunked_output[w_i][]
    end

    # Apply normalization to the output array
    normalize_npcf!(output, npcf)

    if npcf.verb; @printf("Calculations complete!\n"); end

    return output
end;

"""
    compute_npcf_pairwise(inputs, npcf)

Compute the NPCF of the `inputs` dataset by summing over pairs of particles and utilizing hyperspherical harmonic decompositions.
For `n` particles, the computational cost scales as `n^2`, and thus can be extended to large `n`.
Work is distributed approximately equally over all ```nworkers()``` workers provided.

The function is initialized with the `NPCF` class, which specifies all relevant parameters.
The `inputs` array gives a list of particle positions in `D`-dimensions, and a `D+1`-th column specifying the particle weight.
The output is a vector array, where the first `N-1` indices specify the radial bins, and the `N`-th index gives the angular bin, collapsed into one dimension (see ```angular_indices``` for conversion to ells).
Note that we only fill up radial bins satisfying bin1<bin2<...<binN.
The actual computation is performed in the ```_compute_npcf_pairwise()``` function; this is just a wrapper function to facilitate multiprocessing.

# Examples
```julia-repl
julia> inputs = hcat(rand(100,3)*1000.,ones(100)); # 100 randomly distributed 3D particles in [0, 1000], plus weights
julia> npcf = NPCF(N=4, D=3); # create a 4PCF object in 3D
julia> output_pairwise = compute_npcf_pairwise(inputs, npcf); # run the pairwise NPCF estimator
```
"""
function compute_npcf_pairwise(inputs::Matrix{Float64}, npcf::NPCF)

    # First check inputs
    check_inputs(inputs, npcf)

    Npart = length(inputs[:,1])
    if npcf.verb; @printf("\nAnalyzing dataset with %d particles using %d worker(s) by summing over all N-tuplets of galaxies\n",Npart,nworkers()); end

    # Set up output array to hold NPCF
    output = create_npcf_array(npcf)

    # Compute array of weights
    if npcf.N>2
        if npcf.verb; println("Loading coupling weights"); end
        coupling_weights = create_weight_array(npcf)
    else
        # create an empty array for compatibility
        coupling_weights = []
    end

    # Work out how many workers are present and chunk memory
    tasks_per_worker = div(Npart,nworkers())+1

    # Now perform the chunking and run the algorithm across separate workers
    if npcf.verb;
        if nworkers()>1
            println("Starting distributed computation");
        else
            println("Starting computation");
        end
    end
    chunked_output = Dict{Int,Future}()
    @sync for (w_i,w) in enumerate(workers())
        imin = tasks_per_worker*(w_i-1)+1
        imax = min(imin+tasks_per_worker,Npart+1)-1
        @async chunked_output[w_i] = @spawnat w _compute_npcf_pairwise(imin, imax, inputs, coupling_weights, npcf)
    end

    # Now combine the output together
    for w_i in 1:nworkers()
        output .+= chunked_output[w_i][]
    end

    # Apply normalization to the output array
    normalize_npcf!(output, npcf)

    if npcf.verb; @printf("Calculations complete!\n"); end

    return output
end;

## Worker Functions

"""
    _compute_npcf_simple(imin,imax,inputs,npcf)

Compute the NPCF by summing over `N`-tuplets of particles given in the `inputs` dataset using primary particles `imin` to `imax`.
This will run the code for a single worker; see ```compute_npcf_simple()``` for the wrapper and documentation.
"""
function _compute_npcf_simple(imin::Int64, imax::Int64, inputs::Matrix{Float64}, npcf::NPCF)

    Npart = length(inputs[:,1])
    if npcf.verb;
        if nworkers()>1
            println("Running on thread id $(myid()) for particles $imin to $imax");
        else
            println("Iterating over primary particles");
        end
    end

    # Set up output array for this worker
    output = create_npcf_array(npcf)

    # Define positions and weights
    positions = inputs[:,1:npcf.D]
    weights = inputs[:,npcf.D+1]

    # Iterate over first particle in required range
    for i in imin:imax

        p1 = positions[i,:]
        w1 = weights[i]

        if npcf.coords=="spherical"
            # Reorient points such that p1 is at the origin
            shift_positions = reorient_points(positions,p1)
        else
            # Translate positions such that p1 is at the origin
            shift_positions = positions
            for ii in 1:npcf.D
                shift_positions[:,ii] .-= p1[ii]
            end
            if npcf.periodic
                periodic_wrapping!(shift_positions,npcf)
            end
        end

        # Iterate over all secondary particles
        for j in 1:Npart
            p2 = shift_positions[j,:]
            w12 = w1*weights[j]

            # Find separation bin
            sep12 = sep(p2,npcf)
            bin12 = return_bin(sep12,npcf)
            if bin12 == -1; continue; end

            # Accumulate 2PCF
            if npcf.N==2
                output[bin12] += w12
                continue
            end

            # Define unit vector
            if npcf.coords=="cartesian"; p2 /= sep12; end

            # Iterate over third particle
            for k in 1:Npart
                p3 = shift_positions[k,:]
                w123 = w12*weights[k]

                if j==k; continue; end

                # Find separation bin
                sep13 = sep(p3,npcf)
                bin13 = return_bin(sep13,npcf)
                if bin13==-1 || bin13<=bin12; continue; end

                # Define unit vector
                if npcf.coords=="cartesian"; p3 /= sep13; end

                if npcf.N==3
                    # Compute NPCF contribution

                    for l in 0:npcf.lmax
                        # Find basis vector
                        basis = basis_3pcf(l,p2,p3,npcf)

                        # Add to output array
                        output[bin12,bin13,l+1] += w123*conj(basis)
                    end
                else
                    # Iterate over fourth particle
                    for l in 1:Npart
                        p4 = shift_positions[l,:]
                        w1234 = w123*weights[l]

                        if k==l; continue; end

                        # Find separation bin
                        sep14 = sep(p4,npcf)
                        bin14 = return_bin(sep14,npcf)
                        if bin14==-1 || bin14<=bin13; continue; end

                        # Define unit vector
                        if npcf.coords=="cartesian"; p4/=sep14; end

                        if npcf.N==4
                            # Compute NPCF contribution
                            if npcf.D==2
                                l_index = 1
                                for l1_1 in -npcf.lmax:npcf.lmax
                                    for l2_1 in -npcf.lmax:npcf.lmax
                                        l3_1 = -l1_1-l2_1
                                        if abs(l3_1)>npcf.lmax || l3_1>0; continue; end

                                        # Find basis vector
                                        basis = basis_4pcf(l1_1,l2_1,l3_1,p2,p3,p4,npcf)

                                        # Add to output array
                                        output[bin12,bin13,bin14,l_index] += w1234*conj(basis)
                                        l_index += 1
                                    end
                                end
                            elseif npcf.D>=3
                                l_index  = 1
                                for l1_2 in 0:npcf.lmax
                                    for l2_2 in 0:npcf.lmax
                                        for l3_2 in abs(l1_2-l2_2):min(l1_2+l2_2,npcf.lmax)

                                            # remove odd-parity NPCF contributions
                                            if (-1)^(l1_2+l2_2+l3_2)==-1; continue; end

                                            # Find basis vector
                                            basis = basis_4pcf(l1_2,l2_2,l3_2,p2,p3,p4,npcf)

                                            # Add to output array
                                            output[bin12,bin13,bin14,l_index] += w1234*conj(basis)
                                            l_index += 1
                                        end
                                    end
                                end
                            end
                        else
                            # Iterate over fifth particle
                            for m in 1:Npart
                                p5 = shift_positions[m,:]
                                w12345 = w1234*weights[m]

                                if l==m; continue; end

                                # Find separation bin
                                sep15 = sep(p5,npcf)
                                bin15 = return_bin(sep15,npcf)

                                if bin15==-1 || bin15<=bin14; continue; end

                                # Define unit vector
                                if npcf.coords=="cartesian"; p5/=sep15; end

                                if npcf.N==5
                                    # Compute NPCF contribution
                                    if npcf.D==2
                                        l_index = 1
                                        for l1_1 in -npcf.lmax:npcf.lmax
                                            for l2_1 in -npcf.lmax:npcf.lmax
                                                for l3_1 in -npcf.lmax:npcf.lmax
                                                    l4_1 = -l1_1-l2_1-l3_1
                                                    if abs(l4_1)>npcf.lmax || l4_1>0; continue; end

                                                    # Find basis vector
                                                    basis = basis_5pcf(l1_1,l2_1,l1_1+l2_1,l3_1,l4_1,p2,p3,p4,p5,npcf)

                                                    # Add to output array
                                                    output[bin12,bin13,bin14,bin15,l_index] += w12345*conj(basis)
                                                    l_index += 1
                                                end
                                            end
                                        end
                                    elseif npcf.D>=3
                                        l_index  = 1
                                        for l1_2 in 0:npcf.lmax
                                            for l2_2 in 0:npcf.lmax
                                                for l12_2 in abs(l1_2-l2_2):l1_2+l2_2
                                                    for l3_2 in 0:npcf.lmax
                                                        for l4_2 in abs(l12_2-l3_2):min(l12_2+l3_2,npcf.lmax)

                                                            # remove odd-parity NPCF contributions
                                                            if (-1)^(l1_2+l2_2+l3_2+l4_2)==-1; continue; end

                                                            # Find basis vector
                                                            basis = basis_5pcf(l1_2,l2_2,l12_2,l3_2,l4_2,p2,p3,p4,p5,npcf)

                                                            # Add to output array
                                                            output[bin12,bin13,bin14,bin15,l_index] += w12345*conj(basis)
                                                            l_index += 1
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return output
end;

"""
    _compute_npcf_pairwise(imin, imax, inputs, coupling_weights, npcf)

Compute the NPCF of the `inputs` dataset using primary particles `imin` to `imax` by summing over pairs of particles and utilizing hyperspherical harmonic decompositions.
This will run the code for a single worker; see ```compute_npcf_pairwise()``` for the wrapper and documentation.
"""
function _compute_npcf_pairwise(imin::Int64, imax::Int64, inputs::Matrix{Float64}, coupling_weights, npcf::NPCF)

    Npart = length(inputs[:,1])
    if npcf.verb;
        if nworkers()>1
            println("Running on thread id $(myid()) for particles $imin to $imax");
        else
            println("Iterating over primary particles");
        end
    end

    # Set up output array for this worker
    output = create_npcf_array(npcf)

    # Define positions and weights
    positions = inputs[:,1:npcf.D]
    weights = inputs[:,npcf.D+1]

    # Define a_l array
    nl = length(hyperspherical_harmonic(zeros(npcf.D),npcf))
    al = zeros(ComplexF64,npcf.nbins,nl)

    # Iterate over first particle in required range
    for i in imin:imax
        p1 = positions[i,:]
        w1 = weights[i]

        if npcf.coords=="spherical"
            # Rotate particles such that p1 is at the origin
            shift_positions = reorient_points(inputs,p1)
        else
            # Translate particles such that p1 is at the origin
            shift_positions = positions
            for ii in 1:npcf.D
                shift_positions[:,ii] .-= p1[ii]
            end
            if npcf.periodic
                periodic_wrapping!(shift_positions, npcf)
            end
        end

        # Empty array
        al .= 0

        # Iterate over all secondary particles and compute al array
        for j in 1:Npart
            p2 = shift_positions[j,:]
            w2 = weights[j]

            # Find separation bin
            sep12 = sep(p2,npcf)
            bin12 = return_bin(sep12,npcf)
            if bin12 == -1; continue; end

            # Accumulate 2PCF
            if npcf.N==2
                output[bin12] += w1*w2
                continue
            end

            # Define unit vector
            if npcf.coords=="cartesian"; p2/=sep12; end

            # Define spherical harmonics (all at once)
            yl = hyperspherical_harmonic(p2,npcf)

            # Add these to the al array
            al[bin12,:] += w2*conj(yl[:])
        end

        # Now fill up NPCF array
        if npcf.N==2
            # no further computations needed!
            continue
        elseif npcf.N==3
            for b1 in 1:npcf.nbins
                for b2 in b1+1:npcf.nbins
                    # iterate over primary ells
                    for l1 in 0:npcf.lmax

                        # Compute the NPCF contribution in this bin and set of primary ells
                        contribution = accumulate_npcf_contribution([b1,b2],l1,al,coupling_weights,npcf)

                        # Add to output array
                        if contribution!=0; output[b1,b2,l1+1] += w1*contribution; end

                    end
                end
            end
        elseif npcf.N==4
            # Iterate over radial bins
            for b1 in 1:npcf.nbins
                for b2 in b1+1:npcf.nbins
                    for b3 in b2+1:npcf.nbins
                        if npcf.D==2
                            # Iterate over primary ells
                            l_index = 1
                            for l1_1 in -npcf.lmax:npcf.lmax
                                for l2_1 in -npcf.lmax:npcf.lmax
                                    if abs(l1_1+l2_1)>npcf.lmax || l1_1+l2_1<0; continue; end

                                    # Compute the NPCF contribution in this bin and set of primary ells
                                    contribution = accumulate_npcf_contribution([b1,b2,b3],[l1_1,l2_1],al,coupling_weights,npcf)

                                    # Add to output array
                                    if contribution!=0; output[b1,b2,b3,l_index] += w1*contribution; end
                                    l_index += 1

                                end
                            end
                        else
                            # Iterate over primary ells for D>2
                            l_index = 1
                            for l1 in 0:npcf.lmax
                                for l2 in 0:npcf.lmax
                                    for l3 in abs(l1-l2):min(l1+l2,npcf.lmax)
                                        # skip parity-odd contributions
                                        if (-1)^(l1+l2+l3)==-1; continue; end

                                        # Compute the NPCF contribution in this bin and set of primary ells
                                        contribution = accumulate_npcf_contribution([b1,b2,b3],[l1,l2,l3],al,coupling_weights,npcf)

                                        # Add to output array
                                        if contribution!=0; output[b1,b2,b3,l_index] += w1*contribution; end
                                        l_index += 1
                                    end
                                end
                            end
                        end
                    end
                end
            end
        elseif npcf.N==5
            # Iterate over radial bins
            for b1 in 1:npcf.nbins
                for b2 in b1+1:npcf.nbins
                    for b3 in b2+1:npcf.nbins
                        for b4 in b3+1:npcf.nbins
                            if npcf.D==2
                                # Iterate over primary ells
                                l_index = 1
                                for l1_1 in -npcf.lmax:npcf.lmax
                                    for l2_1 in -npcf.lmax:npcf.lmax
                                        for l3_1 in -npcf.lmax:npcf.lmax

                                            if abs(l1_1+l2_1+l3_1)>npcf.lmax || l1_1+l2_1+l3_1<0; continue; end

                                            # Compute the NPCF contribution in this bin and set of primary ells
                                            contribution = accumulate_npcf_contribution([b1,b2,b3,b4],[l1_1,l2_1,l3_1],al,coupling_weights,npcf)

                                            # Add to output array
                                            if contribution!=0; output[b1,b2,b3,b4,l_index] += w1*contribution; end
                                            l_index += 1
                                        end
                                    end
                                end
                            else
                                # Iterate over primary ells for D>2
                                l_index = 1
                                for l1 in 0:npcf.lmax
                                    for l2 in 0:npcf.lmax
                                        for l12 in abs(l1-l2):(l1+l2)
                                            for l3 in 0:npcf.lmax
                                                for l4 in abs(l12-l3):min(l12+l3,npcf.lmax)
                                                    # skip parity-odd contributions
                                                    if (-1)^(l1+l2+l3+l4)==-1; continue; end

                                                    # Compute the NPCF contribution in this bin and set of primary ells
                                                    contribution = accumulate_npcf_contribution([b1,b2,b3,b4],[l1,l2,l12,l3,l4],al,coupling_weights,npcf)

                                                    # Add to output array
                                                    if contribution!=0; output[b1,b2,b3,b4,l_index] += w1*contribution; end
                                                    l_index += 1
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        else
            error("Not yet implemented!")
        end
    end

    return output
end;

### Other functions

"""
    summarize(npcf)

Print a summary of the NPCF hyperparameters
"""
function summarize(npcf::NPCF)
    @printf("Computing the %dPCF in %d dimensions\n\n",npcf.N,npcf.D)
    @printf("Coordinates = %s\n",npcf.coords)
    @printf("Periodic = %s\n",npcf.periodic)
    if npcf.periodic
        @printf("Periodic boxsize: %f\n",npcf._boxsize)
    end
    @printf("Radial bins: [%.2f, %.2f, %.2f]\n",npcf.r_min,npcf.r_max,npcf._dr)
    @printf("l-max: %d\n",npcf.lmax)
end

"""
    angular_indices(npcf)

Return a list of primary angular indices corresponding to the collapsed indices used in the output of the ```compute_npcf_simple``` and ```compute_npcf_pairwise``` functions.
For `N = 3`, we return a single list of `l1`, for `N = 4`, we return a list of `l1, l2, l3`, and for `N = 5`, we return a list of `l1, l2, l12, l3, l4`.
"""
function angular_indices(npcf::NPCF)
    if npcf.N==2
        return error("Isotropic 2PCF contains no angular indices!")
    elseif npcf.N==3
        l1s = Int64[]

        # Fill up array of ells
        l_index = 1
        for l1 in 0:npcf.lmax
            push!(l1s,l1)
        end

        return l1s

    elseif npcf.N==4
        l1s = Int64[]
        l2s = Int64[]
        l3s = Int64[]

        if npcf.D==2

            # Fill up array of ells
            for l1 in -npcf.lmax:npcf.lmax
                for l2 in -npcf.lmax:npcf.lmax
                    l3 = -l1-l2
                    if abs(l3)>npcf.lmax || l3>0; continue; end

                    push!(l1s,l1)
                    push!(l2s,l2)
                    push!(l3s,l3)
                end
            end

        elseif npcf.D>=3

            # Fill up array of ells
            for l1 in 0:npcf.lmax
                for l2 in 0:npcf.lmax
                    for l3 in abs(l1-l2):min(l1+l2,npcf.lmax)
                        # only even-parity NPCF contributions
                        if (-1)^(l1+l2+l3)==-1; continue; end

                        push!(l1s,l1)
                        push!(l2s,l2)
                        push!(l3s,l3)
                    end
                end
            end
        end

        return l1s, l2s, l3s

    elseif npcf.N==5

        if npcf.D==2
            l1s = Int64[]
            l2s = Int64[]
            l3s = Int64[]
            l4s = Int64[]

            # Fill up array of ells
            for l1 in -npcf.lmax:npcf.lmax
                for l2 in -npcf.lmax:npcf.lmax
                    for l3 in -npcf.lmax:npcf.lmax
                        l4 = -l1-l2-l3
                        if abs(l4)>npcf.lmax || l4>0; continue; end

                        push!(l1s,l1)
                        push!(l2s,l2)
                        push!(l3s,l3)
                        push!(l4s,l4)
                    end
                end
            end
            return l1s,l2s,l3s,l4s

        elseif npcf.D>=3
            l1s = Int64[]
            l2s = Int64[]
            l12s = Int64[]
            l3s = Int64[]
            l4s = Int64[]

            # Fill up array of ells
            for l1 in 0:npcf.lmax
                for l2 in 0:npcf.lmax
                    for l12 in abs(l1-l2):l1+l2
                        for l3 in 0:npcf.lmax
                            for l4 in abs(l12-l3):min(l12+l3,npcf.lmax)
                                # only even-parity NPCF contributions
                                if (-1)^(l1+l2+l3+l4)==-1; continue; end

                                push!(l1s,l1)
                                push!(l2s,l2)
                                push!(l12s,l12)
                                push!(l3s,l3)
                                push!(l4s,l4)
                            end
                        end
                    end
                end
            end
            return l1s,l2s,l12s,l3s,l4s
        end
    end
end

end


text:
Faker.language("es-MX")
Faker.seed(1234)

@test Faker.random_digit() == "7"
@test Faker.random_digit_not_null() == "4"
@test Faker.random_digit_or_empty() == " "
@test Faker.random_number(digits= 1) == "6"
@test Faker.random_int(min=0, max=9999) == "4994"
@test Faker.randomize_nb_elements(number=10, le=false, ge=false) == "14.9"
@test Faker.random_letter() == "d"
@test Faker.random_element(("a", "b", "b")) == "b"
@test Faker.lexify("????") == "VTSa"
@test Faker.numerify("####") == "7491"
@test Faker.bothify("## ??") == "26 Eu"

Faker.seed(1234)

@test Faker.state_abbr() == "HID"
@test Faker.latitude() == 2.505859375
@test Faker.street_name() == "Paseo Coquimatlán"
@test Faker.address() == "Escalinata Xochitepec 52 Querétaro, PUE 12657"
@test Faker.street_address() == "Mercado Silvia Espinoza S/N 7ª Planta"
@test Faker.postcode() == "77309"
@test Faker.longitude() == 17.015625
@test Faker.country() == "Kazajistán"
@test Faker.geo_coordinate(center="None", radius=0.001) == Float16(-16.98)
@test Faker.secondary_address() == "Depto. 894"
@test Faker.city_prefix() == "Sur"
@test Faker.city_suffix() == "los bajos"
@test Faker.building_number() == "8"
@test Faker.city() == "Asientos"
@test Faker.state() == "Quintana Roo"

Faker.seed(1234)

@test Faker.company()                                 == "Osorio, Peña y Andrade Asociados S. R.L. de C.V."
@test Faker.company_suffix()                          == "S.A.B."
@test Faker.catch_phrase()                            == "adaptador defectos cero actualizable "
@test Faker.bs()                                      == "optimiza interfaces magnéticas "

Faker.seed(1234)

@test Faker.months() ==8
@test Faker.am_pm() =="PM"
 Faker.iso8601() == "2076-06-20T03:19:24"
 Faker.date_time() =="2041-1-3 08:20:25"
 Faker.month_name() =="March"
 Faker.day_of_week() =="Saturday"
 Faker.day_of_month() == 13
 Faker.date_time_between("-30y", "now") =="1989-6-7 21:38:44"
 Faker.time("H:M:S") =="18:46:36"
 Faker.unix_time() ==454113618933
 Faker.timezone() =="Atlantic/Azores"
 Faker.century() =="V"
 Faker.date("Y-m-d") =="2023-2-19"
 Faker.year() ==2035
 Faker.date_time_this_century(before_now=true, after_now=false) =="2007-1-3 00:00:00"
 Faker.date_time_this_decade(before_now=true, after_now=false) =="2018-7-23 00:00:00"
 Faker.date_time_this_year(before_now=true, after_now=false) =="2018-5-17 00:00:00"
 Faker.date_time_this_month(before_now=true, after_now=false) =="2018-9-1 00:00:00"
Faker.seed(1234)

@test Faker.email()                                   == "kPeña@agyaasrdc.org"
@test Faker.free_email()                              == "rGranados@prodigy.net.mx"
@test Faker.safe_email()                              == "Contreras.Julio César@example.info"
@test Faker.user_name()                               == "odel Ángel"
@test Faker.domain_name()                             == "ddlssrdc.info"
@test Faker.domain_word()                             == "lcs"
@test Faker.domain_suffix()                           ==  "me"
@test Faker.ipv4()                                    ==  "242.225.18.190"
@test Faker.ipv6()                                    == "2f04:35d1:dc5c:57a8:8850:40a5:35e2:c8a7"
@test Faker.uri()                                     == "http://www.áos.mx/"
@test Faker.url()                                     == "http://www.mygs.com.mx/"
@test Faker.image_url()                               == "http://www.lorempixum.com/116/188"
@test Faker.uri_extension()                           == ".html"
@test Faker.uri_path()                                == "blog/app/blog"
@test Faker.uri_page()                                == "privacy"
@test Faker.mac_address()                             == "81:32:31:b1:2a:06"
@test Faker.company_email()                           == "Alfonso.Beltrán@mysa.info"

Faker.seed(1234)

@test Faker.text(number_chars=200) ==  "Laudantium molestias eos rerum. Cupiditate voluptatem exercitationem enim.. Qui qui suscipit aspernatur porro. Non. Officiis consequatur harum architecto quibusdam. Explicabo voluptas cumque distinctio maiores voluptatibus.."
@test Faker.sentence(number_words=6, variable_nb_words=true) == "Distinctio numquam velit voluptatem incidunt."
@test Faker.word() == "eaque"
@test Faker.paragraphs() == ["Numquam unde et dolorem. Rem aut harum perspiciatis. Veritatis..", "Consequuntur consequatur recusandae quis. Nobis dolorem ipsum non ipsam. Quia vel et minus sed. Molestiae velit iure..", "Culpa expedita. Quo voluptatibus voluptatem ea placeat.."]
@test Faker.words(number_words=3) == ["qui", "eos", "sit"]
@test Faker.paragraph(number=3) == "Earum exercitationem architecto deserunt dicta alias. Qui eos sint iusto aut.."
@test Faker.sentences(number_sentences=3) == ["In et voluptatem.", "Nostrum sapiente molestiae.", "Saepe ut in non quasi.", "Ab totam perspiciatis."]

Faker.seed(1234)

@test Faker.phone_number()  == "81-2935-2382"
@test Faker.cell_phone() ==  "044 (223) 714 5444"

Faker.seed(1234)

@test Faker.password() == "D3Wt#5(2&5&b#7'R%9BsA5Xr%*Uu%kRy99htT8Fn"
@test Faker.locale() == "cn_NGA"
@test Faker.null_boolean() == false
@test Faker.country_code() == "SLV"
@test Faker.language_code() == "es"
@test Faker.boolean() == true

Faker.seed(1234)

@test Faker.color_name()                              == "LightYellow"
@test Faker.safe_color_name()                         == "orchid"
@test Faker.hex_color()                               == "#FFFAFA"
@test Faker.rgb_color_list()                          == "(215,221,84)"
@test Faker.rgb_color()                               == "131,80,74"
@test Faker.rgb_css_color()                           == "rgb(68,107,109)"

Faker.seed(1234)

@test Faker.credit_card_security_code()               == "736"
Faker.credit_card_full()                        =="JCB 16 digit Beatriz Acosta 3096935238270743 CVC 571"
Faker.credit_card_number()                      == "3096544451199132"
@test Faker.credit_card_provider()                    == "laser"
@test Faker.generate_number("52",16)                  == "5276029874944399"

Faker.seed(1234)

@test Faker.mime_type()                               == "image/x-freehand"
@test Faker.file_extension()                          == "rmvb"
@test Faker.file_name()                               == "wTPd.pgn"

Faker.seed(1234)

@test Faker.first_name()                              == "Lorena"
@test Faker.last_name()                               == "Peña"
@test Faker.prefixe()                                 == "Dr."
@test Faker.simple_profile()                          == Dict("name"=>"Elizabeth","mail"=>"Omar23@hotmail.com","username"=>"Marcos29","birthdate"=>"2022-2-9","sex"=>"M","address"=>"Avenida Luis Alberto de los Santos 3 5ª Planta Cerritos, SON 34045")
@test Faker.job()                                     == "Administrador de alojamiento"
a= Faker.profile()
b= Dict(
 "username" => "Villegas.María Elena",
  "sex" => "M",
  "address" => "Huerta México 63143 Interior 080 Choix, COL 62719",
  "name" => "Erika",
  "job" => "Agrimensor, minero",
  "mail" => "Rentería.Adriana@prodigy.net.mx",
  "birthdate" => "2027-12-6",
  "current_location" => (-0.9951171875, 8.0078125),
  "blood_group" => "AB-",
  "residence" => "Barranco Querétaro 71590 Planta baja Nazareno Etla, TLA 69478",
  "company" => "Morán y Trujillo S.R.L",
  "website" => "http://www.neh.mx/"
  )

@test a["username"]== b["username"]
@test a["sex"]== b["sex"]
@test a["address"]== b["address"]
@test a["name"]== b["name"]
@test a["job"]== b["job"]
@test a["mail"]== b["mail"]
@test a["birthdate"]== b["birthdate"]
@test a["current_location"]== b["current_location"]
@test a["blood_group"]== b["blood_group"]
@test a["residence"]== b["residence"]
@test a["company"]== b["company"]
@test a["website"]== b["website"]

@test Faker.profile("username")            ==  Dict{Any,Any}("username"=>"Delgado.Gerardo")
@test Faker.profile("username","blood_group")   == Dict{Any,Any}("blood_group"=>"AB+","username"=>"zRamos")
@test Faker.simple_profile("birthdate","address","sex")    ==  Dict{Any,Any}("birthdate"=>"2087-3-15","sex"=>"M","address"=>"Camino Luis Alberto Zárate S/N Int. 9 Bacalar, BCN 72662")
@test Faker.simple_profile("name")    ==  Dict{Any,Any}("name"=>"David")

Faker.seed(1234)

@test Faker.mac_processor()                           == "Intel"
@test Faker.linux_processor()                         == "x86_64"
@test Faker.user_agent()                              == "opera"
@test Faker.windows_platform_token()                  == "Windows CE"
@test Faker.linux_platform_token()                    == "X11; Linux i686"
@test Faker.mac_platform_token()                      ==  "Macintosh; U; Intel Mac OS X 10_5_3"
@test Faker.chrome()                                  =="Mozilla/5.0 (Windows 98; Win 9x 4.90) AppleWebKit/5331 (KHTML, like Gecko) Chrome/32.0.824.0 Safari/5331"
#@test Faker.firefox()                                 =="Mozilla/5.0 (Windows NT 4.0; sl-SI; rv:1.9.0.20) Gecko/2015-4-22 21:10:48 Firefox/6.0"
#@test Faker.safari()      == "Mozilla/5.0 (iPod; U; CPU iPhone OS 3_1 like Mac OS X; it-IT) AppleWebKit/534.28.6 (KHTML, like Gecko) Version/3.0.5 Mobile/8B115 Safari/6534.28.6"
#@test Faker.opera()                             ==  "Opera/9.26.(Windows NT 6.1; sl-SI) Presto/2.9.160 Version/10.00)"
#@test Faker.internet_explorer()                 == "Mozilla/5.0 (compatible; MSIE 9.0; Windows 95; Trident/3.1)"

ext: jl


text:
@enum CSVType singledata pairdata assocdata groupdata

const NumberOrString = Union{Union{T1,Missing},Union{T2,Missing}} where {T1 <: AbstractString, T2 <: Number}

"""
    getfileextension(filepath)

A quick helper to get the file extension of any given path (without the dot).

# Examples
```julia-repl
julia> getfileextension("~/Desktop/text.txt")
"txt"
```
"""
function getfileextension(filepath::AbstractString)
    path,ext = splitext(filepath)
    return String(chop(ext,head=1,tail=0))
end

"""
    getpaths(location; relativetodatabase=false)

Returns database paths that is optionally relative to Clapeyron.jl directory.
If path is a file, then return an Array containing a single path to that file.
If path is a directory, then return an Array containing paths to all csv files in that directory.

# Examples
```julia-repl
julia> getpaths("SAFT/PCSAFT"; relativetodatabase=true)
3-element Array{String,1}:
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_assoc.csv"
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_like.csv"
 "/home/user/.julia/packages/Clapeyron.jl/xxxxx/database/SAFT/PCSAFT/data_PCSAFT_unlike.csv"

```
"""
function getpaths(location::AbstractString; relativetodatabase::Bool=false)
    # We do not use realpath here directly because we want to make the .csv suffix optional.
    filepath = relativetodatabase ? normpath(dirname(pathof(Clapeyron)), "..", "database", location) : location
    isfile(filepath) && return [realpath(filepath)]
    isfile(filepath * ".csv") && return [realpath(filepath * ".csv")]
    if !isdir(filepath)
        relativetodatabase ? error("The path ", location, " does not exist in the Clapeyron database.") :
            error("The path ", location, " does not exist.")
    end
    files = joinpath.(filepath, readdir(filepath))
    return realpath.(files[isfile.(files) .& (getfileextension.(files) .== "csv")])
end

function flattenfilepaths(locations,userlocations)
    res =
                vcat(
                    reduce(vcat,getpaths.(locations; relativetodatabase=true),init = String[]),
                    reduce(vcat,getpaths.(userlocations),init = String[]),
                    String[]
                    )
    return res
end

#function getparams()

"""
    params, sites = getparams(components,locations;kwargs...)

returns a `Dict{String,ClapeyronParam}` containing all the parameters found for the list of components
in the available CSVs. `locations` are the locations relative to `Clapeyron` database. the available keywords are the ones used in [`ParamOptions`](@ref)

if `return_sites` is set to false, `getparams` will only return the found params.
"""
function getparams(components,
                    locations::Array{String,1}=String[];
                    userlocations::Vector{String}=String[],
                    asymmetricparams::Vector{String}=String[],
                    ignore_missing_singleparams::Vector{String}=String[],
                    ignore_headers::Vector{String} =  ["dipprnumber", "smiles"],
                    verbose::Bool=false,
                    species_columnreference::String="species",
                    source_columnreference::String="source",
                    site_columnreference::String="site",
                    group_columnreference::String="groups",
                    normalisecomponents::Bool=true,
                    return_sites::Bool = true,
                    component_delimiter::String = "~|~"
                    )

    options = ParamOptions(;userlocations,
                            asymmetricparams,
                            ignore_missing_singleparams,
                            ignore_headers,
                            verbose,
                            species_columnreference,
                            source_columnreference,
                            site_columnreference,
                            group_columnreference,
                            normalisecomponents,
                            return_sites,
                            component_delimiter)

    # locations is a list of paths relative to the Clapeyron database directory.
    # userlocations is a list of paths input by the user.
    # If parameters exist in multiple files, Clapeyron gives priority to files in later paths.
    # asymmetricparams is a list of parameters for which matrix reflection is disabled.
    # ignore_missingsingleparams gives users the option to disable component existence check in single params.
    return getparams(components,locations,options)

end
function getparams(components::Vector{String},locations::Vector{String},options)
    filepaths = flattenfilepaths(locations,options.userlocations)
    allcomponentsites = findsitesincsvs(components, filepaths,options)
    allparams, paramsourcecsvs, paramsources = createparamarrays(components, filepaths, allcomponentsites,options)
    result = packageparams(allparams, components, allcomponentsites, paramsourcecsvs, paramsources,options)
    if !options.return_sites
        return result
    end
    if any(x isa AssocParam for x in values(result))
        sites = buildsites(result,components,allcomponentsites,options)
        return result,sites
    else
        return result
    end
end

function buildsites(result,components,allcomponentsites,options)
    n_sites_columns = options.n_sites_columns
    v = String[]
    for sitei in allcomponentsites
        append!(v,sitei)
    end
    unique!(v)
    iszero(length(v)) && return SiteParam(components)
    if !any(haskey(result,n_sites_columns[vi]) for vi in v)
        @error """No columns containing number of sites were found. Supposing zero sites.
        If your model doesn't use sites, but some input CSV folders contain Association params. consider using return_sites=false.
        If there are columns containing number of sites, then those aren't recognized. Consider passing a Dict with the mappings between the sites and the name of the column containing the number of said sites in with the n_sites_columns keyword"
        """
        assoc_csv = Set(String[])
        for x in values(result)
            if x isa AssocParam
                for csvx in x.sourcecsvs
                push!(assoc_csv,csvx)
                end
            end
        end
        assoc_csv = collect(assoc_csv)
        @error "Parsed Association CSV were:"
        for csv in assoc_csv
            println(csv)
        end
        return SiteParam(components)
    end

    n_sites_dict = Dict{String,SingleParam{Int}}(vi => result[n_sites_columns[vi]] for vi in v)
    return SiteParam(n_sites_dict,allcomponentsites)
end

function getparams(groups::GroupParam, locations::Vector{String}=String[],options::ParamOptions=DefaultOptions)
    return getparams(groups.flattenedgroups, locations,options)
end

function getparams(components::String, locations::Vector{String}=String[],options::ParamOptions=DefaultOptions)
    return getparams([components],locations,options)
end

function packageparams(allparams::Dict,
    components::Vector{String},
    allcomponentsites::Array{Array{String,1},1},
    paramsourcecsvs::Dict{String,Set{String}},
    paramsources::Dict{String,Set{String}},
    options::ParamOptions = DefaultOptions)

    asymmetricparams = options.asymmetricparams
    ignore_missingsingleparams = options.ignore_missing_singleparams
    # Package params into their respective Structs.
    output = Dict{String,ClapeyronParam}()
    for (param, value) ∈ allparams
        output[param] = pkgparam(param,value,components,allcomponentsites,paramsourcecsvs,paramsources,options)
    end
    return output
end
#SingleParam
function pkgparam(param::String,
    value::Vector{<:NumberOrString},
    components::Vector{String},
    allcomponentsites::Array{Array{String,1},1},
    paramsourcecsvs::Dict{String,Set{String}},
    paramsources::Dict{String,Set{String}},
    options::ParamOptions = DefaultOptions)
    newvalue, ismissingvalues = defaultmissing(value)
    if param ∉ options.ignore_missing_singleparams && any(ismissingvalues)
        error("Missing values exist in single parameter ", param, ": ", value, ".")
    end
    return SingleParam(param, components, newvalue, ismissingvalues, collect(paramsourcecsvs[param]), collect(paramsources[param]))
end
#PairParam
function pkgparam(param::String,
    value::Matrix{<:NumberOrString},
    components::Vector{String},
    allcomponentsites::Array{Array{String,1},1},
    paramsourcecsvs::Dict{String,Set{String}},
    paramsources::Dict{String,Set{String}},
    options::ParamOptions = DefaultOptions)

    param ∉ options.asymmetricparams && mirrormatrix!(value)
    newvalue, ismissingvalues = defaultmissing(value)
    if (param ∉ options.ignore_missing_singleparams
        && !all([ismissingvalues[x,x] for x ∈ 1:size(ismissingvalues,1)])
        && any([ismissingvalues[x,x] for x ∈ 1:size(ismissingvalues,1)]))
        error("Partial missing values exist in diagonal of pair parameter ", param, ": ", [value[x,x] for x ∈ 1:size(ismissingvalues,1)], ".")
    end
    diagvalues = view(newvalue,diagind(newvalue))
    return PairParam(param,components, newvalue,diagvalues,ismissingvalues, collect(paramsourcecsvs[param]), collect(paramsources[param]))
end
#AssocParam
function pkgparam(param::String,
    value::Matrix{<:Matrix{<:NumberOrString}},
    components::Vector{String},
    allcomponentsites::Array{Array{String,1},1},
    paramsourcecsvs::Dict{String,Set{String}},
    paramsources::Dict{String,Set{String}},
    options::ParamOptions = DefaultOptions)

    param ∉ options.asymmetricparams && mirrormatrix!(value)
    newvalue_ismissingvalues = defaultmissing.(value)
    newvalue = first.(newvalue_ismissingvalues)
    return AssocParam(param,components, newvalue , allcomponentsites, collect(paramsourcecsvs[param]), collect(paramsources[param]))
end

function param_type(t1,t2)
    t_promoted = promote_type(t1,t2)
    local res
    if t_promoted !== Any
        res =  Union{Missing,t_promoted}
    else
        res = Union{Missing,t1,t2}
    end
    return res
end

#try not to promote to Int64. this will error if a user pass integer values different from 0 and 1
function param_type(::Type{Bool},::Type{Int})
    return Union{Bool,Missing}
end

function createparamarrays(components::Array{String,1}, filepaths::Array{String,1}, allcomponentsites::Array{Array{String,1},1}, options::ParamOptions = DefaultOptions)
    # Returns Dict with all parameters in their respective arrays.
    verbose = options.verbose
    check_clashingheaders(filepaths,options)
    allparams = Dict{String,Any}()
    paramsourcecsvs = Dict{String,Set{String}}()
    paramsources = Dict{String,Set{String}}()
    # Read the filepaths in reverse in order to ensure that unused sources do not get added.

    for filepath ∈ reverse(filepaths)
        csvtype = readcsvtype(filepath)
        if csvtype == groupdata
            verbose && @info("Skipping groupdata csv $filepath")
            continue
        end
        headerparams = readheaderparams(filepath,options)
        verbose && @info("Searching for $(Symbol(csvtype)) headers $headerparams for components $components at $filepath ...")
        foundparams, paramtypes, sources = findparamsincsv(components, filepath,options)
        foundparams = swapdictorder(foundparams)
        for headerparam ∈ headerparams
            if !haskey(allparams, headerparam)
                if ismissing(paramtypes[headerparam])
                    allparams[headerparam] = createemptyparamsarray(missing, csvtype, components, allcomponentsites)
                else
                    allparams[headerparam] = createemptyparamsarray(paramtypes[headerparam], csvtype, components, allcomponentsites)
                end
            end
            if !haskey(paramsourcecsvs, headerparam)
                paramsourcecsvs[headerparam] = Set{String}()
                paramsources[headerparam] = Set{String}()
            end
            if csvtype == singledata
                isempty(foundparams) && continue
                for (component, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(allparams[headerparam]))
                    tx = param_type(paramtypes[headerparam],currenttype)
                    allparams[headerparam] = convert(Array{tx}, allparams[headerparam])
                    idx = findfirst(isequal(component), components)
                    if allparams[headerparam] isa Matrix
                        !ismissing(allparams[headerparam][idx,idx]) && continue
                        allparams[headerparam][idx,idx] = value
                    else
                        !ismissing(allparams[headerparam][idx]) && continue
                        allparams[headerparam][idx] = value
                    end
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[component]) && push!(paramsources[headerparam], sources[component])
                end
            end
            if csvtype == pairdata
                if allparams[headerparam] isa Vector
                    allparams[headerparam] = singletopair(allparams[headerparam], missing)
                end
                isempty(foundparams) && continue
                for (componentpair, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(allparams[headerparam]))
                    tx = param_type(paramtypes[headerparam],currenttype)
                    allparams[headerparam] = convert(Matrix{tx}, allparams[headerparam])
                    idx1 = findfirst(isequal(componentpair[1]), components)
                    idx2 = findfirst(isequal(componentpair[2]), components)
                    !ismissing(allparams[headerparam][idx1,idx2]) && continue
                    allparams[headerparam][idx1,idx2] = value
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[componentpair]) && push!(paramsources[headerparam], sources[componentpair])
                end
            end
            if csvtype == assocdata
                isempty(foundparams) && continue
                for (assocpair, value) ∈ foundparams[headerparam]
                    currenttype = nonmissingtype(eltype(first(allparams[headerparam])))
                    tx = param_type(paramtypes[headerparam],currenttype)
                    allparams[headerparam] = convert(Array{Array{tx,2},2}, allparams[headerparam])
                    idx1 = findfirst(isequal(assocpair[1][1]), components)
                    idx2 = findfirst(isequal(assocpair[1][2]), components)
                    idx21 = findfirst(isequal(assocpair[2][1]), allcomponentsites[idx1])
                    idx22 = findfirst(isequal(assocpair[2][2]), allcomponentsites[idx2])
                    !ismissing(allparams[headerparam][idx1,idx2][idx21,idx22]) && continue
                    allparams[headerparam][idx1,idx2][idx21,idx22] = value
                    push!(paramsourcecsvs[headerparam], filepath)
                    !ismissing(sources[assocpair]) && push!(paramsources[headerparam], sources[assocpair])
                end
            end
        end
    end
    return allparams, paramsourcecsvs, paramsources
end

function defaultmissing(array::Array{<:Number},defaultvalue = zero(eltype(array)))
    return deepcopy(array),Array(ismissing.(array))
end

function defaultmissing(array::Array{<:AbstractString},defaultvalue = "")
    return string.(array),Array(ismissing.(array))
end

function defaultmissing(array::Array{String},defaultvalue = "")
    return deepcopy(array),Array(ismissing.(array))
end

function defaultmissing(array::Array{Union{Missing, T}},defaultvalue="") where T <:AbstractString
    return string.(coalesce.(array,defaultvalue)),Array(ismissing.(array))
end

function defaultmissing(array::Array{Union{Missing, T}},defaultvalue=zero(eltype(array))) where T<:Number
    return coalesce.(array,defaultvalue),Array(ismissing.(array))
end

function defaultmissing(array::Array{Union{Missing, Bool}},defaultvalue=zero(eltype(array)))
    return [coalesce(i,defaultvalue) for i in array],Array(ismissing.(array))
end

#if an array with only missings is passed, the Resulting ClapeyronParam will be
#of the type that this function returns
function defaultmissing(array::Array{Missing},defaultvalue=0.0)
    return coalesce.(array,defaultvalue),Array(ismissing.(array))
end

function defaultmissing(array::Array{Union{Missing,T1,T2}},defaultvalue="") where {T1 <:Number,T2<:AbstractString}
    return string.(coalesce.(array,defaultvalue)),Array(ismissing.(array))
end

function defaultmissing(array::Array{Any},defaultvalue="")
    return string.(coalesce.(array,defaultvalue)),Array(ismissing.(array))
end
function defaultmissing(array::Array{T},defaultvalue::T2) where T<:Union{T2,Missing} where T2
    coalesce.(array,Ref(defaultvalue)),Array(ismissing.(array))
end
function defaultmissing(array)
    throw("Unsupported array element type  $(typeof(array))")
end

function swapdictorder(dict)
    # Swap the first two levels in a nested dictionary.
    isempty(dict) && return dict
    output = Dict()
    outerkeys = keys(dict)
    try
        innerkeys = keys(first(values(dict)))
        for innerkey ∈ innerkeys, outerkey ∈ outerkeys
            if !haskey(output, innerkey)
                output[innerkey] = Dict{Any,Any}(outerkey => dict[outerkey][innerkey])
            end
            push!(output[innerkey], outerkey => dict[outerkey][innerkey])
        end
    catch e
        error("The format of the nested dictionary is not correct.")
    end
    return output
end

function col_indices(csvtype,headernames,options=DefaultOptions)
    columnreference = options.species_columnreference
    normalised_columnreference = normalisestring(columnreference)

    idx_species = 0
    idx_groups = 0
    idx_species1 = 0
    idx_species2 = 0
    idx_sites1 = 0
    idx_sites2 = 0

    if csvtype === singledata || csvtype === groupdata
        lookupcolumnindex = findfirst(isequal(normalised_columnreference), headernames)
        isnothing(lookupcolumnindex) && error("Header ", normalised_columnreference, " not found.")
        idx_species = lookupcolumnindex
        if csvtype === groupdata
            groupcolumnreference= options.group_columnreference
            normalised_groupcolumnreference = normalisestring(groupcolumnreference)
            lookupgroupcolumnindex = findfirst(isequal(normalised_groupcolumnreference), headernames)
            isnothing(lookupgroupcolumnindex) && error("Header ", normalised_groupcolumnreference, " not found.")
            idx_groups = lookupgroupcolumnindex
        end

    elseif csvtype === pairdata || csvtype == assocdata
        normalised_columnreference1 = normalised_columnreference * '1'
        normalised_columnreference2 = normalised_columnreference * '2'
        lookupcolumnindex1 = findfirst(isequal(normalised_columnreference1), headernames)
        lookupcolumnindex2 = findfirst(isequal(normalised_columnreference2), headernames)
        isnothing(lookupcolumnindex1) && error("Header ", normalised_columnreference1, " not found.")
        isnothing(lookupcolumnindex2) && error("Header ", normalised_columnreference2, " not found.")
        idx_species1 = lookupcolumnindex1
        idx_species2 = lookupcolumnindex2
        if csvtype == assocdata
            sitecolumnreference = options.site_columnreference
            normalised_sitecolumnreference = normalisestring(sitecolumnreference)
            normalised_sitecolumnreference1 = normalised_sitecolumnreference * '1'
            normalised_sitecolumnreference2 = normalised_sitecolumnreference * '2'
            lookupsitecolumnindex1 = findfirst(isequal(normalised_sitecolumnreference1), headernames)
            lookupsitecolumnindex2 = findfirst(isequal(normalised_sitecolumnreference2), headernames)
            isnothing(lookupsitecolumnindex1) && error("Header ", normalised_sitecolumnreference1, " not found.")
            isnothing(lookupsitecolumnindex2) && error("Header ", normalised_sitecolumnreference2, " not found.")
            idx_sites1 = lookupsitecolumnindex1
            idx_sites2 = lookupsitecolumnindex2
        end
    end

    _single = (idx_species,idx_groups)
    _pair = (idx_species1,idx_species2)
    _assoc = (idx_sites1,idx_sites2)
    return (_single,_pair,_assoc)
end

function findparamsincsv(components::Array{String,1},
    filepath::AbstractString,
    options::ParamOptions = DefaultOptions)

    headerparams = readheaderparams(filepath,options)
    sourcecolumnreference = options.source_columnreference
    verbose = options.verbose
    normalisecomponents = options.normalisecomponents
    component_delimiter = options.component_delimiter

    csvtype = readcsvtype(filepath)
    df = CSV.File(filepath; header=3, pool=0,silencewarnings=true)
    csvheaders = String.(Tables.columnnames(df))
    normalised_components = normalisestring.(components,normalisecomponents)
    normalised_csvheaders = normalisestring.(csvheaders)
    normalised_headerparams = normalisestring.(headerparams)
    normalised_headerparams ⊈ normalised_csvheaders && error("Headers ", setdiff(normalised_headerparams, normalised_csvheaders), " not present in csv header.")

    foundvalues = Dict()
    paramtypes = Dict(headerparams .=> [Tables.columntype(df, Symbol(x)) for x ∈ headerparams])
    sources = Dict{Any,Union{String,Missing}}()

    normalised_sourcecolumnreference = normalisestring(sourcecolumnreference)
    getsources = false
    if normalised_sourcecolumnreference ∈ normalised_csvheaders
        getsources = true
        sourcecolumn = findfirst(isequal(normalised_sourcecolumnreference), normalised_csvheaders)
    end

    single_idx,pair_idx,assoc_idx = col_indices(csvtype,normalised_csvheaders,options)
    lookupcolumnindex,_ = single_idx
    lookupcolumnindex1,lookupcolumnindex2 = pair_idx
    lookupsitecolumnindex1,lookupsitecolumnindex2 = assoc_idx
    headerparams_indices = [findfirst(isequal(i),normalised_csvheaders) for i in normalised_headerparams]
    if csvtype == singledata
        for row ∈ Tables.rows(df)
            component_split = split(row[lookupcolumnindex], component_delimiter, keepempty=false)
            for component ∈ component_split
                foundcomponentidx = findfirst(isequal(normalisestring(component,normalisecomponents)), normalised_components)
                isnothing(foundcomponentidx) && continue

                component = components[foundcomponentidx]
                foundvalues[component] = Dict{String,Any}()
                for (headerparam,idx) ∈ zip(headerparams,headerparams_indices)
                    foundvalues[component][headerparam] = row[idx]
                end
                verbose && @info("""Found component: $(component)
                with values: $(foundvalues[component])
                """)
                if getsources
                    source = row[sourcecolumn]
                    sources[component] = source
                else
                    sources[component] = missing
                end
            end
        end
    elseif csvtype == pairdata
        for row ∈ Tables.rows(df)
            component_split1 = split(row[lookupcolumnindex1], component_delimiter, keepempty=false)
            for component1 ∈ component_split1
                foundcomponentidx1 = findfirst(isequal(normalisestring(component1,normalisecomponents)), normalised_components)
                isnothing(foundcomponentidx1) && continue
                component_split2 = split(row[lookupcolumnindex2], component_delimiter, keepempty=false)
                for component2 ∈ component_split2
                    foundcomponentidx1 = findfirst(isequal(normalisestring(component1,normalisecomponents)), normalised_components)
                    foundcomponentidx2 = findfirst(isequal(normalisestring(component2,normalisecomponents)), normalised_components)
                    isnothing(foundcomponentidx2) && continue
                    componentpair = (components[foundcomponentidx1], components[foundcomponentidx2])
                    foundvalues[componentpair] = Dict{String,Any}()
                    for (headerparam,idx) ∈ zip(headerparams,headerparams_indices)
                        foundvalues[componentpair][headerparam] = row[idx]
                    end
                    verbose && @info("""Found component pair: ($(component1),$(component2))
                    with values: $(foundvalues[componentpair])
                    """)
                    if getsources
                        source = row[sourcecolumn]
                        sources[componentpair] = source
                    else
                        sources[componentpair] = missing
                    end
                end
            end
        end
    elseif csvtype == assocdata
        for row ∈ Tables.rows(df)
            component_split1 = split(row[lookupcolumnindex1], component_delimiter, keepempty=false)
            for component1 ∈ component_split1
                foundcomponentidx1 = findfirst(isequal(normalisestring(component1,normalisecomponents)), normalised_components)
                isnothing(foundcomponentidx1) && continue
                component_split2 = split(row[lookupcolumnindex2], component_delimiter, keepempty=false)
                for component2 ∈ component_split2
                    foundcomponentidx2 = findfirst(isequal(normalisestring(component2,normalisecomponents)), normalised_components)
                    isnothing(foundcomponentidx2) && continue
                    site1 = row[lookupsitecolumnindex1]
                    site2 = row[lookupsitecolumnindex2]
                    assocpair = ((components[foundcomponentidx1], components[foundcomponentidx2]), (site1, site2))
                    foundvalues[assocpair] = Dict{String,Any}()
                    for (headerparam,idx) ∈ zip(headerparams,headerparams_indices)
                        foundvalues[assocpair][headerparam] = row[idx]
                    end
                    verbose && @info("""Found assoc pair: ($(component1),$(component2)) , ($(site1),$(site2))
                    with values: $(foundvalues[assocpair])
                    """)
                    if getsources
                        source = row[sourcecolumn]
                        sources[assocpair] = source
                    else
                        sources[assocpair] = missing
                    end
                end
            end
        end
    else
        error("File is of type ", String(csvtype), " and cannot be read with this function.")
    end
    return foundvalues, paramtypes, sources
end

function normalisestring(str::AbstractString, isactivated::Bool=true; tofilter::Regex=r"[ \-\_]", changecase::Bool=true)
    !isactivated && return str
    normalisedstring = replace(str, tofilter => "")
    changecase && (normalisedstring = lowercase(normalisedstring))
    return normalisedstring
end

const readcsvtype_keywords  = ["like", "single", "unlike", "pair", "assoc", "group", "groups"]

function readcsvtype(filepath)
    # Searches for type from second line of CSV.
    keywords = readcsvtype_keywords
    words = split(lowercase(rstrip(getline(String(filepath), 2), ',')), ' ')
    foundkeywords = intersect(words, keywords)
    isempty(foundkeywords) && error("Unable to determine type of database", filepath, ". Check that keyword is present on Line 2.")
    length(foundkeywords) > 1 && error("Multiple keywords found in database ", filepath, ": ", foundkeywords)
    _readcsvtype(only(foundkeywords))
end

function _readcsvtype(key)
    key == "single" && return singledata
    key == "like" && return singledata
    key == "pair" && return pairdata
    key == "unlike" && return pairdata
    key == "assoc" && return assocdata
    key == "group" && return groupdata
    key == "groups" && return groupdata
    error("Unable to determine database type of $key")
end

function getline(filepath::AbstractString, selectedline::Int)
    # Simple function to return text from filepath at selectedline.
    open(filepath) do file
        linecount = 1
        for line ∈ eachline(file)
            linecount == selectedline && return line
            linecount += 1
        end
        error("Selected line number exceeds number of lines in file")
    end
end

function readheaderparams(filepath::AbstractString, options::ParamOptions = DefaultOptions,headerline::Int = 3)
    # Returns array of filtered header strings at line 3.
    ignorelist = deepcopy(options.ignore_headers)
    push!(ignorelist,options.species_columnreference)
    push!(ignorelist,options.source_columnreference)
    push!(ignorelist,options.site_columnreference)
    headers = split(getline(filepath, headerline), ',',keepempty = false)
    while last(headers) == ""
        pop!(headers)
    end
    return String.(filter(x -> normalisestring(x; tofilter=r"[ \-\_\d]") ∉ ignorelist, headers))
end

function check_clashingheaders(filepaths::Vector{String},
                                options::ParamOptions = DefaultOptions)
    # Raises an error if the header of any assoc parameter clashes with a non-assoc parameter.
    headerparams = String[]
    headerparams_assoc = String[]
    for filepath in filepaths
        csvtype = readcsvtype(filepath)
        if csvtype == singledata || csvtype == pairdata
            append!(headerparams, readheaderparams(filepath,options))
        elseif csvtype == assocdata
            append!(headerparams_assoc, readheaderparams(filepath,options))
        end
    end
    clashingheaders = intersect(headerparams, headerparams_assoc)
    !isempty(clashingheaders) && error("Headers ", clashingheaders, " appear in both loaded assoc and non-assoc files.")
end

function findsitesincsvs(components::Array{String,1},
                        filepaths::Array{String,1},
                        options::ParamOptions = DefaultOptions)

    verbose = options.verbose
    normalisecomponents = options.normalisecomponents
    # Look for all relevant sites in the database.
    # Note that this might not necessarily include all sites associated with a component.
    normalised_components = normalisestring.(components,normalisecomponents)
    sites = Dict(components .=> [Set{String}() for _ ∈ 1:length(components)])
    for filepath ∈ filepaths
        csvtype = readcsvtype(filepath)
        csvtype != assocdata && continue
        df = CSV.File(filepath; header=3,silencewarnings = !verbose)
        csvheaders = String.(Tables.columnnames(df))
        normalised_csvheaders = normalisestring.(csvheaders)
        _,pair_idx,assoc_idx = col_indices(csvtype,normalised_csvheaders,options)
        lookupcolumnindex1,lookupcolumnindex2 = pair_idx
        lookupsitecolumnindex1,lookupsitecolumnindex2 = assoc_idx
        for row ∈ Tables.rows(df)
            component1 = row[lookupcolumnindex1]
            component2 = row[lookupcolumnindex2]
            foundcomponentidx1 = findfirst(isequal(normalisestring(component1,normalisecomponents)), normalised_components)
            foundcomponentidx2 = findfirst(isequal(normalisestring(component2,normalisecomponents)), normalised_components)
            (isnothing(foundcomponentidx1) || isnothing(foundcomponentidx2)) && continue
            push!(sites[components[foundcomponentidx1]], row[lookupsitecolumnindex1])
            push!(sites[components[foundcomponentidx2]], row[lookupsitecolumnindex2])
        end
    end
output = Array{Array{String,1}}(undef, 0)
for component ∈ components
    push!(output, collect(sites[component]))
end
verbose && @info("Found sites for $components are $(output).")
return output
end

function findgroupsincsv(components::Vector{String},
                        filepath::String,
                        options::ParamOptions = DefaultOptions)

    verbose = options.verbose
    normalisecomponents = options.normalisecomponents
    csvtype = readcsvtype(filepath)
    csvtype != groupdata && return Dict{String,String}()
    normalised_components = normalisestring.(components,normalisecomponents)
    df = CSV.File(filepath; header=3,silencewarnings = !verbose)
    columns = Tables.columns(df)
    csvheaders = String.(Tables.columnnames(df))
    normalised_csvheaders = normalisestring.(csvheaders)
    single_idx,_,_ = col_indices(csvtype,normalised_csvheaders,options)
    lookupcolumnindex,lookupgroupcolumnindex = single_idx
    species_column = Tables.getcolumn(columns,lookupcolumnindex)
    groups_column = Tables.getcolumn(columns,lookupgroupcolumnindex)
    norm_species_column = normalisestring.(species_column,normalisecomponents)
    idx = findall(in(normalised_components),norm_species_column)
    found_comps = @view species_column[idx]
    found_groups = @view groups_column[idx]
    if verbose
        for i in 1:length(found_comps)
            @info("""Found component: $(found_comps[i])
            with groups: $(found_groups[i])
            """)
        end
    end
    foundgroups = Dict(comp_i => group_i for (comp_i,group_i) in zip(found_comps,found_groups))
    return foundgroups
end

function createemptyparamsarray(datatype::Type, csvtype::CSVType, components::Array{String,1})
    allcomponentsites = [String[] for _ in 1:length(components)]
    return createemptyparamsarray(datatype, csvtype, components, allcomponentsites)
end

function createemptyparamsarray(datatype::Type, csvtype::CSVType, components::Array{String,1}, allcomponentsites::Array{Array{String,1},1})
    # Creates a missing array of the appropriate size.
    componentslength = length(components)
    csvtype == singledata && return (Array{Union{Missing,datatype}}(undef, componentslength) .= missing)
    csvtype == pairdata && return (Array{Union{Missing,datatype}}(undef, componentslength, componentslength) .= missing)
    if csvtype == assocdata
        output = Array{Array{Union{Missing,datatype},2},2}(undef, componentslength, componentslength)
        for i ∈ 1:componentslength, j ∈ 1:componentslength
            output[i,j] = (Array{Union{Missing,datatype}}(undef, length(allcomponentsites[i]), length(allcomponentsites[j])) .= missing)
        end
        return output
    end
end

function createemptyparamsarray(csvtype::CSVType, components::Array{String,1})
    allcomponentsites = [String[] for _ in 1:length(components)]
    return createemptyparamsarray(missing, csvtype, components, allcomponentsites)
end

function createemptyparamsarray(csvtype::CSVType, components::Array{String,1}, allcomponentsites::Array{Array{String,1},1})
    return createemptyparamsarray(missing, csvtype, components, allcomponentsites)
end

_zero(t::Number) = zero(t)
_zero(::String) = ""
_zero(::Missing) = missing
function _zero(x::Type{T})  where T<:Number
    return zero(T)
end
_zero(::Type{String}) = ""
_zero(::Type{Missing}) = missing
function _zero(::Type{T}) where T <:Union{T1,Missing} where T1
    return missing
end

_iszero(t::Number) = iszero(t)
_iszero(::Missing) = true
_iszero(t::AbstractString) = isempty(t)
"""
    singletopair(params::Vector,outputmissing=zero(T))

Generates a square matrix, filled with "zeros" (considering the "zero" of a string, a empty string).
The generated matrix will have the values of `params` in the diagonal.
If missing is passed, the matrix will be filled with `missing`
"""
function singletopair(params::Vector{T1},::T2 =_zero(T1)) where {T1,T2}
    len = length(params)
    T = Union{T1,T2}
    output = Matrix{T}(undef,len,len)
    fill!(output,_zero(T))
    @inbounds  for i in 1:len
        output[i,i] = params[i]
    end
    return output
end

function mirrormatrix!(matrix::Matrix{T}) where T
    # Mirrors a square matrix.
    matrixsize = size(matrix)
    matrixsize[1] != matrixsize[2] && error("Matrix is not square.")
    for i ∈ 2:matrixsize[1], j ∈ 1:i-1
        lowervalue = matrix[i,j]
        uppervalue = matrix[j,i]
        if !ismissing(lowervalue) && !ismissing(uppervalue) && lowervalue != uppervalue
            error("Dissimilar non-zero entries exist across diagonal:", lowervalue, ", ", uppervalue)
        end
        !ismissing(lowervalue) && (matrix[j,i] = lowervalue)
        !ismissing(uppervalue) && (matrix[i,j] = uppervalue)
    end
    return matrix
end

function mirrormatrix!(matrix::Array{Array{T,2},2}) where T
    # Mirrors a square matrix.
    matrixsize = size(matrix)
    matrixsize[1] != matrixsize[2] && error("Matrix is not square.")
    [mirrormatrix!(matrix[i,i]) for i ∈ 1:matrixsize[1]]
    for i ∈ 2:matrixsize[1], j ∈ 1:i-1
        matrixsize2 = size(matrix[i,j])
        for a ∈ 1:matrixsize2[1], b ∈ 1:matrixsize2[2]
            lowervalue = matrix[i,j][a,b]
            uppervalue = matrix[j,i][b,a]
            if !ismissing(lowervalue) && !ismissing(uppervalue) && lowervalue != uppervalue
                error("Dissimilar non-zero entries exist across diagonal:", lowervalue, ", ", uppervalue)
            end
            !ismissing(lowervalue) && (matrix[j,i][b,a] = lowervalue)
            !ismissing(uppervalue) && (matrix[i,j][a,b] = uppervalue)
        end
    end
    return matrix
end

function GroupParam(gccomponents::Vector,
    grouplocations::Vector{String}=String[];
    usergrouplocations::Vector{String}=String[],
    verbose::Bool=false)
    options = ParamOptions(;usergrouplocations,verbose)
    return GroupParam(gccomponents,grouplocations,options)
end

function GroupParam(gccomponents,
    grouplocations::Array{String,1}=String[],
    options::ParamOptions = DefaultOptions)
    # The format for gccomponents is an arary of either the species name (if it
    # available in the Clapeyron database, or a tuple consisting of the species
    # name, followed by a list of group => multiplicity pairs.  For example:
    # gccomponents = ["ethane",
    #                ("hexane", ["CH3" => 2, "CH2" => 4]),
    #                ("octane", ["CH3" => 2, "CH2" => 6])]
    usergrouplocations = options.usergrouplocations
    verbose = options.verbose
    BuildSpeciesType = Union{Tuple{String, Array{Pair{String, Int64},1}}, String, Tuple{String}}
    any(.!(isa.(gccomponents, BuildSpeciesType))) && error("The format of the components is incorrect.")
    filepaths = flattenfilepaths(grouplocations,usergrouplocations)
    componentstolookup = String[]
    #componentstolookup = filter(x-> x isa Union{String,Tuple{String}},gccomponents)
    append!(componentstolookup, [x for x ∈ gccomponents[isa.(gccomponents, String)]])
    append!(componentstolookup, [first(x) for x ∈ gccomponents[isa.(gccomponents, Tuple{String})]])
    allfoundcomponentgroups = Dict{String,String}()
    groupsourcecsvs = String[]
    for filepath in filepaths
        csvtype = readcsvtype(filepath)
        if csvtype != groupdata
            verbose && @info("Skipping $csvtype csv at $filepath")
            continue
        end
        verbose && @info("Searching for groups for components $componentstolookup at $filepath ...")
        merge!(allfoundcomponentgroups,  findgroupsincsv(componentstolookup, filepath, options))
        append!(groupsourcecsvs, [filepath])
    end
    gccomponents_parsed = PARSED_GROUP_VECTOR_TYPE(undef,length(gccomponents))
    for (i,gccomponent) ∈ pairs(gccomponents)
        if gccomponent isa Tuple{String, Array{Pair{String, Int64},1}}
            gccomponents_parsed[i] = gccomponent
        elseif gccomponent isa String
            !haskey(allfoundcomponentgroups, gccomponent) && error("Predefined component ", gccomponent, " not found in any group input csvs.")
            groupsandngroups = eval(Meta.parse(allfoundcomponentgroups[gccomponent]))
            gccomponents_parsed[i] = (gccomponent,groupsandngroups)
        elseif gccomponent isa Tuple{String}
            !haskey(allfoundcomponentgroups, first(gccomponent)) && error("Predefined component ", gccomponent, " not found in any group input csvs.")
            groupsandngroups = eval(Meta.parse(allfoundcomponentgroups[first(gccomponent)]))
            gccomponents_parsed[i] = (first(gccomponent),groupsandngroups)
        end
    end
    return GroupParam(gccomponents_parsed,groupsourcecsvs)
end

export getparams


text:
abstract type System end

namefor(s::S) where {S<:System} = namefor(S)
namefor(S) = nameof(S)
typefor(s::S) where {S<:System} = S
Base.names(s::S) where {S<:System} = names(S)
Base.names(S::Type{<:System}) = (n = split(String(Symbol(S)), "."); [Symbol(join(n[i:end], ".")) for i in 1:length(n)])

Base.length(s::System) = length(fieldnamesunique(s))
Base.iterate(s::System) = iterate(s, 1)
Base.iterate(s::System, i) = begin
    F = fieldnamesunique(s)
    l = length(F)
    l == 0 ? nothing : (s[F[i]], l == i ? nothing : i+1)
end
Base.iterate(s::System, ::Nothing) = nothing

Base.broadcastable(s::System) = Ref(s)

Base.getindex(s::System, i) = getproperty(s, i)
Base.getindex(s::System, ::Nothing) = s

Base.getproperty(s::System, n::AbstractString) = begin
    isempty(n) && return s
    reduce((a, b) -> begin
        m = match(r"([^\[\]]+)(?:\[(.+)\])?", b)
        n, i = m[1], m[2]
        v = getfield(a, Symbol(n))
        if isnothing(i)
            v
        else
            #HACK: support symbol as-is (i.e. "a[:i]" vs. "a[i]")
            k = if startswith(i, ":")
                i
            elseif isprivatename(i)
                canonicalname(i, n)
            else
                try
                    #HACK: support indexing of non-Variable (i.e. "a[1]" for Vector{Layer})
                    parse(Int, i)
                catch
                    #HACK: support accessing index (i.e. "a[i]")
                    getfield(a, Symbol(i)) |> value
                end
            end
            v[k]
        end
    end, [s, split(n, ".")...])
end

Base.hasproperty(s::System, n::AbstractString) = begin
    try
        getproperty(s, n)
    catch
        return false
    end
    true
end

#HACK: calculate variable body with external arguments for debugging purpose
value(s::System, k::Symbol; kw...) = begin
    d = dependency(s)
    v = d.M[k]
    @assert v.state in (:Preserve, :Track)
    emit(a) = let p = extractfuncargpair(a), k = p[1]; :($k = $(kw[k])) end
    args = emit.(v.args)
    body = v.body
    eval(:(let $(args...); $body end))
end

Base.show(io::IO, s::System) = print(io, "<$(namefor(s))>")
Base.show(io::IO, ::MIME"text/plain", s::System) = look(io, s; header=false, doc=false)

#TODO: see if we can move it to util/look.jl
include("look.jl")

labelstring(v; maxlength=nothing) = begin
    l = repr(v; context=IOContext(devnull, :compact => true, :limit => true))
    n = length(l)
    i = findfirst('\n', l)
    i = isnothing(i) ? n : i-1
    x = isnothing(maxlength) ? n : maxlength
    i = min(i, x)
    i < n ? l[1:i] * "…" : l
end

export System


text:
module DGMethods

using MPI
using ..MPIStateArrays
using ..Mesh.Grids
using ..Mesh.Topologies
using StaticArrays
using ..VariableTemplates
using DocStringExtensions
using KernelAbstractions
using KernelAbstractions.Extras: @unroll

export BalanceLaw,
    DGModel, init_ode_state, restart_ode_state, restart_auxiliary_state

include("balance_law.jl")
include("NumericalFluxes.jl")
include("DGModel.jl")
include("DGModel_kernels.jl")

end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# return the indices of sample inputs in `obj::MwlsNaiveObject` within `dist` from `pt`
function inrange_data(obj::MwlsNaiveObject, pt::Point, dist::Real = obj.EPS)
  res::Vector{Int} = []

  @inbounds for p in 1:size(obj.inputs, 2)
    if norm(obj.inputs[:, p] - pt) < dist + dist / 1e8
      push!(res, p)
    end
  end

  return res
end

# return the indices of sample inputs in `obj::MwlsNaiveObject` within `dist` from `pt`
function inrange_data(obj::MwlsKdObject, pt::Point, dist::Real = obj.EPS)
  return inrange(obj.tree, pt, dist)
end

# return the indices of sample inputs in `obj::MwlsNaiveObject` within `dist` from `pt`
function inrange_data(obj::MwlsCllObject, pt::Point, dist::Real = obj.EPS)
  return cll_inrange(obj.cll, pt, dist)
end

getInrangeData(a...;b...) = warn("`getInrangeData` is deprecated, use` inrange_data` instead.")

"""
    mwls_coefficients(obj::MwlsObject, pt::Point, dist::Real)

Calculates the coefficients of the linear combination of polynomials used for approximation.
This is done for each dimension of output data.

!!! note
    If the matrix in the system of linear equations used to find the coefficients is singular, then zero coefficients are returned!
"""
function mwls_coefficients(obj::MwlsObject, pt::Point, dist::Real = obj.EPS)
  m = size(obj.b, 1)
  outputDim = size(obj.outputs, 1)
  firstTerm = zeros(m, m)
  secondTerm = zeros(outputDim, m)

  data = inrange_data(obj, pt, dist)

  for p in data
    curPt = obj.inputs[:, p]
    w = obj.weightfunc(norm(pt - curPt), obj.EPS)
    @inbounds for j in 1:m
      @inbounds for i in 1:m
        firstTerm[i, j] += w * obj.matrix[i, j](obj.vars => curPt)
      end
      secondTerm[:, j] += w * obj.b[j](obj.vars => curPt) * obj.outputs[:, p]
    end
  end

  # probably singular matrix
  # returns zeroes
  if cond(firstTerm, 1) >= 1e14 || cond(firstTerm, 2) >= 1e14
    return zeros(m, outputDim)
  end

  result = firstTerm \ transpose(secondTerm)
  return result
end

calcMwlsCoefficients(a...;b...) = warn("`calcMwlsCoefficients` is deprecated, use `mwls_coefficients` instead.")

# does the actual approximation
"""
    approximate(obj::MwlsObject, pt::Point)
    approximate(obj::MwlsObject, pt::Point; dist::Real = obj.EPS)

This calculates the approximated value at `pt` for each dimension of output data.
The approximated value is returned.
"""
function approximate(obj::MwlsObject, pt::Point, dist::Real = obj.EPS)
  cs = mwls_coefficients(obj, pt, dist)
  poly = [polynomial(cs[:, i], obj.b) for i in 1:size(cs, 2)]
  res = [p(obj.vars => pt) for p in poly]
  return length(res) == 1 ? res[1] : res
end

function (obj::MwlsNaiveObject)(pt::Point, dist::Real = obj.EPS)
  approximate(obj, pt, dist)
end

function (obj::MwlsNaiveObject)(pt::Real, dist::Real = obj.EPS)
  return obj([pt], dist)
end

function (obj::MwlsKdObject)(pt::Point, dist::Real = obj.EPS)
  approximate(obj, pt, dist)
end

function (obj::MwlsKdObject)(pt::Real, dist::Real = obj.EPS)
  return obj([pt, 0], dist)
end

function (obj::MwlsCllObject)(pt::Point, dist::Real = obj.EPS)
  approximate(obj, pt, dist)
end

# TODO: huh
function (obj::MwlsCllObject)(pt::Real, dist::Real = obj.EPS)
  return obj([pt, 0], dist)
end

"""
    mwls_diff_polys(obj::MwlsObject, inPt::Point, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}

Polynomials created by [mwls_coefficients](@ref) are differentiated according to dirs.
For an example if ``f`` is the polynomial used for approximation, then ``dirs = (1,2)`` returns ``\\frac{\\partial^3}{\\partial x_1 \\partial x_2^2}f``.
"""
function mwls_diff_polys(obj::MwlsObject, inPt::Point, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
  cs = mwls_coefficients(obj, inPt, dist)
  poly = [polynomial(cs[:, i], obj.b) for i in 1:size(cs, 2)]
  @inbounds for j in 1:length(poly)
    @inbounds for i in 1:length(dirs)
      poly[j] = differentiate(poly[j], obj.vars[i], dirs[i])
    end
  end
  return poly
end

calcDiffMwlsPolys(a...;b...) = warn("`calcDiffMwlsPolys` is deprecated, use `mwls_diff_polys` instead")

"""
    mwls_diff(obj::MwlsNaiveObject, inPt::Real, dirs::Integer; dist::Real = obj.EPS)
"""
function mwls_diff(obj::MwlsNaiveObject, inPt::Real, dirs::Integer; dist::Real = obj.EPS)
  mwls_diff(obj, [inPt], ntuple(x -> dirs, 1); dist = dist)
end

"""
    mwls_diff(obj::MwlsNaiveObject, inPt::Point, dirs::Integer; dist::Real = obj.EPS)
"""
function mwls_diff(obj::MwlsNaiveObject, inPt::Point, dirs::Integer; dist::Real = obj.EPS)
  mwls_diff(obj, inPt, ntuple(x -> dirs, 1); dist = dist)
end

"""
    mwls_diff(obj::MwlsNaiveObject, inPt::Real, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
"""
function mwls_diff(obj::MwlsNaiveObject, inPt::Real, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
  mwls_diff(obj, [inPt], dirs, dist)
end

"""
    mwls_diff(obj::MwlsObject, inPt::Real, dirs::Integer; dist::Real = obj.EPS)
"""
function mwls_diff(obj::MwlsObject, inPt::Real, dirs::Integer; dist::Real = obj.EPS)
  mwls_diff(obj, [inPt, 0], ntuple(x -> dirs, 1); dist = dist)
end

"""
    mwls_diff(obj::MwlsObject, inPt::Real, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
"""
function mwls_diff(obj::MwlsObject, inPt::Real, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
  mwls_diff(obj, [inPt, 0], ntuple(x -> dirs, 1); dist = dist)
end

"""
    mwls_diff(obj::MwlsObject, inPt::Point, dirs::Integer; dist::Real = obj.EPS)
"""
function mwls_diff(obj::MwlsObject, inPt::Point, dirs::Integer; dist::Real = obj.EPS)
  mwls_diff(obj, inPt, ntuple(x -> dirs, 1); dist = dist)
end

"""
    mwls_diff(obj::MwlsObject, inPt::Point, dirs::NTuple{N, Integer}; dist::Real = obj.EPS)

Calculates the approximated derivative at `inPt`, where `x[i]` is differentiated `dirs[i]` times.
"""
function mwls_diff(obj::MwlsObject, inPt::Point, dirs::NTuple{N, Integer}; dist::Real = obj.EPS) where {N}
  N != length(obj.vars) && error("Mismatch between tuple size and amount of variables")
  pl = mwls_diff_polys(obj, inPt, dirs; dist = dist)
  res = [p(obj.vars => inPt) for p in pl]
  return length(res) == 1 ? res[1] : res
end

mwlsDiff(a...;b...) = warn("`mwlsDiff` is deprecated, use `mwls_diff` instead")


text:
module Test_Dynamic_Pendulum_xAxis

using  Modia3D
import Modia3D.ModiaMath

groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=5*Ly, m=1.0) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(frame0; r=[0.0, 0.0, -Lx/2])
   frame2 = Modia3D.Object3D(frame0; r=[0.0, 0.0, Lx/2])
   cyl    = Modia3D.Cylinder(Ly/2,Lx; material=vmat2)
   cyl1   = Modia3D.Object3D(frame1, cyl; R=ModiaMath.rot2(-pi/2), visualizeFrame=false)
   cyl2   = Modia3D.Object3D(frame2, cyl; R=ModiaMath.rot2(-pi/2), visualizeFrame=false)
end

@assembly DoublePendulum(;Lx = 1.0, m=1.0) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))
   bar1  = Bar(Lx=Lx, m=m)
   bar2  = Bar(Lx=Lx, m=m)
   rev1  = Modia3D.Revolute(world, bar1.frame1; axis = 1)
   rev2  = Modia3D.Revolute(bar1.frame2, bar2.frame1; axis = 1)
end

gravField = Modia3D.UniformGravityField(n=[0,-1,0])
doublePendulum = DoublePendulum(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.3))
model = Modia3D.SimulationModel( doublePendulum)
result = ModiaMath.simulate!(model; stopTime=5.0, tolerance=1e-6,interval=0.001, log=false)

ModiaMath.plot(result, [("rev1.phi", "rev2.phi"),
                        ("rev1.w"  , "rev2.w"),
                        ("rev1.a"  , "rev2.a")])

println("... success of Test_Dynamic_Pendulum_xAxis.jl!")
end


text:
using Test
using WorldBankData

# checks the iso 3-letter country codes for updates
all_iso3c = ["ABW", "AFE", "AFG", "AFR", "AFW", "AGO", "ALB", "AND", "ARB", "ARE", "ARG", "ARM", "ASM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEA", "BEC", "BEL", "BEN", "BFA", "BGD", "BGR", "BHI", "BHR", "BHS", "BIH", "BLA", "BLR", "BLZ", "BMN", "BMU", "BOL", "BRA", "BRB", "BRN", "BSS", "BTN", "BWA", "CAA", "CAF", "CAN", "CEA", "CEB", "CEU", "CHE", "CHI", "CHL", "CHN", "CIV", "CLA", "CME", "CMR", "COD", "COG", "COL", "COM", "CPV", "CRI", "CSA", "CSS", "CUB", "CUW", "CYM", "CYP", "CZE", "DEA", "DEC", "DEU", "DFS", "DJI", "DLA", "DMA", "DMN", "DNF", "DNK", "DNS", "DOM", "DSA", "DSF", "DSS", "DZA", "EAP", "EAR", "EAS", "ECA", "ECS", "ECU", "EGY", "EMU", "ERI", "ESP", "EST", "ETH", "EUU", "FCS", "FIN", "FJI", "FRA", "FRO", "FSM", "FXS", "GAB", "GBR", "GEO", "GHA", "GIB", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GRL", "GTM", "GUM", "GUY", "HIC", "HKG", "HND", "HPC", "HRV", "HTI", "HUN", "IBB", "IBD", "IBT", "IDA", "IDB", "IDN", "IDX", "IMN", "IND", "INX", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAC", "LAO", "LBN", "LBR", "LBY", "LCA", "LCN", "LDC", "LIC", "LIE", "LKA", "LMC", "LMY", "LSO", "LTE", "LTU", "LUX", "LVA", "MAC", "MAF", "MAR", "MCO", "MDA", "MDE", "MDG", "MDV", "MEA", "MEX", "MHL", "MIC", "MKD", "MLI", "MLT", "MMR", "MNA", "MNE", "MNG", "MNP", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAC", "NAF", "NAM", "NCL", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRS", "NRU", "NXS", "NZL", "OED", "OMN", "OSS", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRE", "PRI", "PRK", "PRT", "PRY", "PSE", "PSS", "PST", "PYF", "QAT", "ROU", "RRS", "RUS", "RWA", "SAS", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SMR", "SOM", "SRB", "SSA", "SSD", "SSF", "SST", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SXM", "SXZ", "SYC", "SYR", "TCA", "TCD", "TEA", "TEC", "TGO", "THA", "TJK", "TKM", "TLA", "TLS", "TMN", "TON", "TSA", "TSS", "TTO", "TUN", "TUR", "TUV", "TWN", "TZA", "UGA", "UKR", "UMC", "URY", "USA", "UZB", "VCT", "VEN", "VGB", "VIR", "VNM", "VUT", "WLD", "WSM", "XKX", "XZN", "YEM", "ZAF", "ZMB", "ZWE"]

@test WorldBankData.download_countries(verbose=true)[!, :iso3c] == all_iso3c


text:
importall Bukdu

type WelcomeController <: ApplicationController
end

layout(::Layout, body) = body
index(::WelcomeController) = render(Text/Layout, "hello world")
show(::WelcomeController) = pi

Router() do
    get("/", WelcomeController, index)
    get("/pi", WelcomeController, show)
end

Endpoint() do
    plug(Plug.Logger, level=:error)
    plug(Router)
end

using Base.Test
import HttpCommon: Request, Response

req = Request()

req.method = "GET"
req.resource = "/"
res = Bukdu.handler(req, Response())
@test 200 == res.status
@test "hello world" == String(res.data)

req.method = "GET"
req.resource = "/pi"
res = Bukdu.handler(req, Response())
@test 200 == res.status
@test "π = 3.1415926535897..." == String(res.data)

req.method = "POST"
req.resource = "/"
res = Bukdu.handler(req, Response())
@test 404 == res.status

req.method = "GET"
req.resource = "/test"
res = Bukdu.handler(req, Response())
@test 404 == res.status

import Requests: statuscode, text

Bukdu.start([8082, 8083])

resp1 = Requests.get("http://localhost:8082/")
resp2 = Requests.get("http://localhost:8083/")
@test 200 == statuscode(resp1)
@test 200 == statuscode(resp2)
@test "hello world" == text(resp1)
@test "hello world" == text(resp2)

req.method = "GET"
req.resource = "/"

(server,task) = first(Bukdu.Farm.servers)
@test :runnable == task.state
@test "hello world" == text(server.http.handle(req, Response()))

sleep(0.1)
Bukdu.stop()

logs = []

before(::Request, ::Response) = push!(logs, :br)
before(::WelcomeController) = push!(logs, :bc)
after(::WelcomeController) = push!(logs, :ac)
after(::Request, ::Response) = push!(logs, :ar)

conn = (Router)(get, "/")
@test [:bc, :ac] == logs
empty!(logs)

before(render, Text/Layout) do text
    push!(logs, :bvl)
end
before(render, Text) do text
    push!(logs, :bv)
end
after(render, Text) do text
    push!(logs, :av)
end
after(render, Text/Layout) do text
    push!(logs, :avl)
end

Bukdu.start(8082)
resp1 = Requests.get("http://localhost:8082/")
@test_throws Base.UVError Requests.get("http://localhost:8083/")
@test 200 == statuscode(resp1)
@test "hello world" == text(resp1)
@test [:br,:bc,:bvl,:bv,:av,:avl,:ac,:ar] == logs
sleep(0.1)
Bukdu.stop()