ternarydiagrams.jl

## ---
    using StatGeochem, Plots

    # Calculate a histogram of three variables in ternary space
    function ternaryhist(a,b,c; nbins::Int=10)
        # Normalize input data
        d = a + b + c
        a ./= d
        b ./= d
        c ./= d
        # abc = collect(zip(a./d,b./d,c./d))

        # Count them up
        n = 0 # triangle number
        edges = collect(0:1/nbins:1)
        N = Array{Int64}(undef,nbins*nbins)
        intriangle = BitArray(undef,length(d))
        # rightside-up triangles
        for i = 1:nbins
            for j = 1:i
                n += 1
                # > > >
                @inbounds intriangle .= (b .>= edges[1+nbins-i]) .& (a .>= edges[1+i-j]) .& (c .>= edges[j])
                N[n] = count(intriangle)
                # intriangle = 0
                # @inbounds for k = 1:length(a)
                #     if (b[k] >= edges[1+nbins-i]) && (a[k] >= edges[1+i-j]) && (c[k] >= edges[j])
                #         intriangle += 1
                #     end
                # end
                # N[n] = intriangle
                # N[n] = count(x -> (x[2] >= edges[1+nbins-i]) & (x[1] >= edges[1+i-j]) & (x[3] >= edges[j]), abc)
            end
        end

        # upside-down triangles
        for i = 1:(nbins-1)
            for j = 1:i
                n += 1
                # < < <
                @inbounds intriangle .= (b .< edges[1+nbins-i]) .& (a .< edges[2+i-j]) .& (c .< edges[1+j])
                N[n] = count(intriangle)
                # intriangle = 0
                # @inbounds for k = 1:length(a)
                #     if (b[k] < edges[1+nbins-i]) && (a[k] < edges[2+i-j]) && (c[k] < edges[1+j])
                #         intriangle += 1
                #     end
                # end
                # N[n] = intriangle
                # N[n] = count(x -> (x[2] < edges[1+nbins-i]) & (x[1] < edges[2+i-j]) & (x[3] < edges[1+j]), abc)
            end
        end

        return N
    end

    # Plot a histogram of three variables in ternary space
    function ternaryplot(N; cornerlabels=["","",""], nbins::Int=isqrt(length(N)), cmap=viridis, clims=(0,maximum(N)))

        # Scale the counts
        cind = ceil.(Int, (N .- clims[1]) * (length(cmap)-1) / clims[2]) .+ 1
        cind[cind .> length(cmap)] .= length(cmap)
        cind[cind .< 1] .= 1

        # Fill in null color by default
        h = plot(size=(600,519), xlims=(-1/2, 1/2), ylims=(0, sqrt(3)/2), xticks=[], yticks=[], framestyle=:none)
        plot!(h, [-1/2, 0, 1/2, -1/2,],[0, sqrt(3)/2, 0, 0,],label="",fill=true,fillcolor=cmap[1],linewidth=0)

        # Plot the histogram
        n = 0 # triangle number
        edges = collect(0:1/nbins:1)

        # rightside-up triangles
        for i = 1:nbins
            for j = 1:i
                n += 1
                # > > > If there's anything there
                if cind[n] > 1
                    y = 0.5*sqrt(3)*[edges[1+nbins-i], edges[2+nbins-i], edges[1+nbins-i], edges[1+nbins-i]]
                    x = -0.5*[edges[2+i-j], edges[1+i-j], edges[1+i-j], edges[2+i-j]] + 0.5*[edges[j], edges[j], edges[1+j], edges[j]]
                    plot!(h,x,y,label="",fill=true,color=cmap[cind[n]],linewidth=0)
                end
            end
        end

        # Upside-down triangles
        for i = 1:(nbins-1)
            for j = 1:i
                n += 1
                # < < < If there's anything there
                if cind[n] > 1
                    y = 0.5*sqrt(3)*[edges[1+nbins-i], edges[nbins-i], edges[1+nbins-i], edges[1+nbins-i]]
                    x = -0.5*[edges[2+i-j], edges[1+i-j], edges[1+i-j], edges[2+i-j]] + 0.5*[edges[j], edges[j], edges[1+j], edges[j]]
                    plot!(h,x,y,label="",fill=true,color=cmap[cind[n]],linewidth=0)
                end
            end
        end

        # Add lines on margins
        plot!(h, [-1/2, 0, 1/2, -1/2,], [0, sqrt(3)/2, 0, 0,], label="", color=:black)

        # Add corner labels
        plot!(h, [-0.515, 0, 0.515], [-0.015, sqrt(3)/2+0.015, -0.015], text=cornerlabels, seriestype=:scatter, label="")

        return h
    end

    # Plot three variables in ternary space
    function ternaryplot(a,b,c; cornerlabels=["","",""], markersize=0.5, color=:auto, label="")
        # Calculate sum of components
        d = a + b + c

        # x-axis of ABC diagram
        x = -0.5 * a./d + 0.5 * c./d
        # y-axis of ABC diagram
        y = 0.5 * sqrt(3) * b./d

        # Plot dataset
        h = plot(xlims=(-1/2, 1/2), ylims=(0, sqrt(3)/2), xticks=[], yticks=[], framestyle=:none, size=(600, 519))
        plot!(h,x,y,seriestype=scatter, label=label, markersize=markersize, color=color, mswidth=0)
        # Add lines on margins
        plot!(h, [-1/2, 0, 1/2, -1/2,],[0, sqrt(3)/2, 0, 0,],color=:black, label="")
        # Add corner labels
        plot!(h, [-0.515, 0, 0.515], [-0.015, sqrt(3)/2+0.015, -0.015], text=cornerlabels, seriestype=:scatter, label="")
        return h
    end

## --- Test it out

    npoints = 10^5
    a = abs.(randn(npoints))
    b = abs.(randn(npoints))
    c = abs.(randn(npoints))

    @time h = ternaryplot(a,b,c)
    savefig(h,"TernaryPlotTest.pdf")

    @time N = ternaryhist(a,b,c,nbins=100)
    colors = copy(viridis)
    colors[1] = eltype(viridis)(0,0,0)
    @time h = ternaryplot(N,cornerlabels=["Ca","Na","K"],cmap=colors)
    savefig(h,"TernaryHistTest.pdf")
    display(h)

## ---