Simplex heatmap

bryce · July 1, 2022, 7:53pm

I’m trying to generate a heatmap to visualize the output of a function evaluated on a grid of points over the unit 2-simplex. Here’s an example of the sort of data I’m working with:

using Plots, Combinatorics, StatsBase

function simplex_grid(num_outcomes, points_per_dim)
    num_points = multinomial(num_outcomes-1, points_per_dim)
    points = Array{Float64,2}(undef, num_outcomes, num_points)
    for (p,comb) in enumerate(with_replacement_combinations(1:num_outcomes, points_per_dim))
        distr = counts(comb, 1:num_outcomes) ./ points_per_dim
        points[:,p] = distr
    end
    return points
end

function bary2cart(points)
    x = (2 .* points[3,:] .+ points[1,:]) ./ dropdims(sum(points, dims=1), dims=1) ./ 2
    y = √3 / 2 .* points[1,:] ./ dropdims(sum(points, dims=1), dims=1);
    return (x,y)
end

points_per_dim = 30
num_outcomes = 3
grid = simplex_grid(num_outcomes, points_per_dim)
(x,y) = bary2cart(grid)
z = rand(length(x))

The best solutions I’ve come up with so far are 1) a scatter plot with hexagonal markers:

plot(showaxis=false, grid=false, ticks=false, xlims=(-.08,1.08), ylims=(-.08,.√.75+.08))
scatter!(x,y, zcolor=z, markerstrokewidth=0, markershape=:hexagon, markersize=8, label="")
plot!([0,.5,1,0],[0,√.75,0,0], color=:black, markeralpha=0, label="")

annotate!(.5, √.75+.07, text("A", 12))
annotate!(-.05, -.05, text("B", 12))
annotate!(1.05,-.05, text("C", 12))

hexagons

or 2) a surface with the camera directly above:

plot(showaxis=false, grid=false, ticks=false, xlims=(-.08,1.08), ylims=(-.08,.√.75+.08))
surface!(x,y,z, zcolor=z, camera=(0,90))
plot!([0,.5,1,0],[0,√.75,0,0], color=:black, markeralpha=0, label="")

surface

But both of these have significant drawbacks. The hexagons tend to have gaps or overlaps, and don’t allow any interpolation. The surface doesn’t allow annotations, and has weirdly jagged edges. Any suggestions for how to improve on (or replace) either of these would be appreciated!

Related but less important: I also tried plotting the barycentric coordinates directly, but couldn’t figure out how to get the view aligned orthogonal to the plane of the simplex. Is there anywhere I can find a basic explanation of how the “camera” keyword works?

juliohm · July 1, 2022, 8:09pm

Are you after a ternary plot?

There are open issues in plotting packages to support these maps. I know that TernaryPlots.jl had something working at some point:

Also, Zulip has a discussion with various more advanced tricks:

https://julialang.zulipchat.com/#narrow/stream/225542-helpdesk/topic/Can.20a.20recipe.20be.20used.20for.20ternary.20plots.3F

joa-quim · July 1, 2022, 10:04pm

Examples of ternary plots Ternary plots with GMT.jl

bryce · July 6, 2022, 5:22pm

I had looked into TernaryPlots.jl, but heatmaps are listed under “work in progress”.

juliohm · July 6, 2022, 5:29pm

What about the thread in Zulip? Maybe there is a heatmap there? It is been a while I read it…

bryce · July 6, 2022, 5:42pm

GMT seems powerful, but its Julia interface seems very poorly documented. I’ve managed to get the following figure, but I can’t figure out how to annotate it. I’d like to add lots of labels, starting with something like the A/B/C annotations from the example above. Any suggestions?

using GMT
grid = simplex_grid(num_outcomes, points_per_dim)
m = Matrix{Float64}(undef,496,4)
m[:,1:3] .= grid'
m[:,4] = generate_data(grid)
C = makecpt(T=(0,maximum(m[:,4])*1.001))
ternary(m, show=true, marker=:p, image=true, cmap=C, markersize=.01, region=(0,1,0,1,0,1), frame=(grid=0.1, ticks=0.1, annot=0.5))

If possible, I’d also like to add annotations using LatexStrings, such as:

using LatexStrings
annotate!(0.5, -.05, text(L"\left< 0, \frac{1}{2}, \frac{1}{2}\right>", 12))

but since that doesn’t seem to work right in any of the plotting libraries I wouldn’t be surprised if it were just totally infeasible.

rafael.guerra · July 6, 2022, 8:35pm

Fyi, it seems to work well with Plots.jl’s pgfplotsx() backend:

using LaTeXStrings, Plots; pgfplotsx()
heatmap(rand(10,10));
annotate!(5.5, 5.5, text(L"\left< 0, \frac{1}{2}, \frac{1}{2}\right>", 16))

joa-quim · July 6, 2022, 9:46pm

I know that the documentation is lagging behind but no need to exaggerate also. Doesn’t the link that I posted above show numerous examples on how to add … A/B/C labeling and other examples? (posting it again)

https://www.generic-mapping-tools.org/GMTjl_doc/examples/ternary/ternary_examples/

All that respect frames is described in the manual … and we can’t say it’s a short description.

Regarding LaTex, that’s a different issue. Yes, it should be possible. GMT does it but depends on a full LaTex installation, which something that doesn’t really pleases me. But so far it seems the only option. Put the latex string in a label="@[\nabla^4 \psi - \Delta \sigma_{xx}^2@[ (MPa)" like in the GMT example I linked to. Never tried but no reason it should not work … as long as the GMT conditions are met.

And

help?> ternary

gives also some info

bryce · July 7, 2022, 3:53pm

I appreciate your help, and I don’t want to get into an argument over tone on the internet, but please understand that it’s incredibly frustrating to be told to RTFM when you’ve already been trying to understand that manual for many hours. There’s certainly tons of information in the general GMT manual, but it’s quite challenging to parse if you don’t already know the ins and outs of GMT. Specifically, it’s very hard to figure out how things will translate to the Julia interface and it’s very hard to track down the right search terms for features one doesn’t already know.

Here’s my current status:

ternary(m, show=true, image=true, region=(0,1,0,1,0,1), frame=(grid=0, ticks=0.5, annot=0), vertex_labels="3/1/2", labels=(["0, .5, .5"],[".5, 0, .5"],[".5, .5, 0"]))

There are lots of things I’d like to add, but the most pressing is to make it clear that the labels like [0, .5, .5] refer to the point on the edge of the simplex. I could imagine several ways of doing this, but haven’t managed to get any of them working. One option might be various GMT config options. It seems like this line should make the ticks larger (and if it worked, there’d be a number of other interesting config options I’d like to explore), but it doesn’t appear to change anything:

gmtset(MAP_TICK_LENGTH=2)

Or if I could overlay a scatter plot on top of the heatmap, I could do something like this to emphasize those points:

ternary([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0], frame=(grid=0, ticks=0, annot=0), noclip=true, marker=:p, markersize=.2, show=true)

But what seems to me like the obvious way to do this in Julia (calling ternary!) doesn’t work, and my searching hasn’t turned up any alternatives.

Many of the other things I’d like to do could be managed if I had a way to add arbitrary text to an arbitrary location in the figure, but again, I can’t find the right GMT command to do that. Any tips would be wonderful; thanks again!

joa-quim · July 7, 2022, 4:36pm

Sorry for your doc-frustation. I share it because I try to keep adding layers of it and seem to be failing.

The gmtset(MAP_TICK_LENGTH=2,) should work, but note the trailing comma needed because it’s a named tuple, but looks it’s not. I will have to look at that. But you can obtain a similar effect (not exactly equal because the gmtset form is supposed to be permanent until the figure is generated) with the par=( put-your-settings-here ) alternative.

Regarding plotting a scatter plot on an image. Isn’t what that example does? Or you want a different image?

ternary("@ternary.txt", marker=:p, image=true, clockwise=true,
        frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt",
        clabel="Sand", suffix=" %"), par=(MAP_TICK_LENGTH=2,), show=true)

joa-quim · July 7, 2022, 4:47pm

I can obtain the same figure as above by splitting it into two commands. This means you can get an image with one dataset and do a scatter plot with another.

ternary("@ternary.txt", image=true, clockwise=true,
        frame=(annot=:auto, grid=:a, ticks=:a, alabel="Clay", blabel="Silt", 
        clabel="Sand", suffix=" %"), par=(MAP_TICK_LENGTH=2,))

ternary!("@ternary.txt", marker=:p, show=true)

bryce · July 7, 2022, 8:27pm

I had show=true in both commands , and the result was an uninformative error message.

joa-quim · July 7, 2022, 11:16pm

Nothing I can do to improve that error message. The ! form tells to append to a figure. If that figure was not initialized with a call to a method without the ! form the resulting postscript file is invalid and the error you saw is issued by ghostscript, not GMT.

joa-quim · July 8, 2022, 4:12pm

@bryce, ternary man page. Feel free to propose improvements to the parts that made you suffer more.

bryce · July 8, 2022, 11:47pm

Thanks, I’ll take a look and see if I have any suggestions.

On a related note, when I try to add another series to the plot with ternary!, I’m sometimes getting the following error:

psbasemap [ERROR]: Bad interval in -B option (x-component, a-info): (null) gave interval = 0
psbasemap [ERROR]: Bad interval in -B option (y-component, a-info): (null) gave interval = 0
psbasemap [ERROR]: Option -B parsing failure. Correct syntax:
	-B[p|s][x|y|z]<intervals>[+a<angle>|n|p][+l|L<label>][+p<prefix>][+u<unit>] -B[<axes>][+b][+g<fill>][+o<lon>/<lat>][+s<subtitle>][+t<title>][+w[<pen>]][+x<fill>][+y<fill>][+z<fill>] OR
	-B[p|s][x|y|z][a|f|g]<tick>[m][l|p] -B[p|s][x|y|z][+l<label>][+p<prefix>][+u<unit>] -B[<axes>][+b][+g<fill>][+o<lon>/<lat>][+s<subtitle>][+t<title>][+w[<pen>]][+x<fill>][+y<fill>][+z<fill>]
psbasemap [ERROR]: Offending option -B(null)
psternary [ERROR]: Unable to plot A axis

Here’s a minimal example that produces it:

ternary([.1 .3 .6], marker=:a)
ternary!([.2 .4 .4], marker=:c, show=true)

Any idea what’s going on here or any tips on how to avoid it?

joa-quim · July 8, 2022, 11:51pm

This?
https://www.generic-mapping-tools.org/GMTjl_doc/examples/ternary/ternary_examples/#diagram_with_lines_and_text