Density plot of chaos game using Python matplotlib?

Question

Essentially all I'm trying to do is produce as set of points via an IFS and use a color map to show the multiplicity of each point. In other words, if we assume a color map where high values are more yellow and lower ones are more red, then values repeatedly produced by the IFS will be more yellow.

I'm struggling to get correct results for this. Each thing I've tried has resulted in an image that looks interesting, but is clearly incorrect as it differs wildly from what you get from simply plotting the points without color mapping.

Below is the base code that I'm comfortable with, without the failed attempts at color mapping. What can I do to get a proper color map?

The basic strategy, I think, is to make a matrix 'mat' holding the point multiplicities and do something like plt.imshow(xs, ys, c=mat. cmap="..."). I've tried different approaches to this but keep coming up with incorrect results.

import numpy as np
import matplotlib.pyplot as plt
import random

def f(x, y, n):
    N = np.array([[x, y]])

    M = np.array([[1, 0], [0, 1]])

    b = np.array([[.5], [0]])
    b2 = np.array([[0], [.5]])

    if n == 0:
        return np.dot(M, N.T)

    elif n == 1:
       return np.dot(M, N.T) + b

    elif n == 2:
        return np.dot(M, N.T) + b2

    elif n == 3:
        return np.dot(M, N.T) - b

    elif n == 4:
        return np.dot(M, N.T) - b2


xs = [] # x coordinates
ys = [] # y coordinates
D = {}  # point multiplicities 

random.seed()

x = 1
y = 1

for i in range(0, 100000):
    n = random.randint(1, 4)

    V = f(x, y, n)

    x = V.item(0)
    y = V.item(1)

    xs.append(x)
    ys.append(y)

    xi = round(x, 3)
    yi = round(y, 3)

    if (xi, yi) in D:
        D[(xi, yi)] += 1
    else:   
        D[(xi, yi)] = 1

plt.xlabel('x')
plt.ylabel('y')

plt.scatter(xs,ys, s=.05)
plt.autoscale(True, True, True)
plt.show()

Answer 1

If I understand your problem, it sounds like you want to use a 2D histogram to get the density of points,

H, x, y = np.histogram2d(xs,ys,bins=100)
X, Y = np.meshgrid(x[:-1],y[:-1],indexing='ij')
plt.pcolormesh(X,Y,H,alpha=0.8, cmap = plt.cm.YlOrRd_r)
plt.colorbar()

Which gives,

This is a transparent colormesh plotted over the scatter plot. You could also colour your scatter plot by the value at point,

pc = some_fn_to_get_color_at_points(X, Y, H, xs, yx)
plt.scatter(xs,ys, s=.05, c=pc)