VMLS: Lighting problem

The data for this plot comes from Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares (Boyd and Vandenberghe, 2018).

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 480

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt

def get_pattern(power, weights):
    return np.array([power[i] * weights[i] for i in range(len(power))]).sum(axis=0).reshape(power.shape[1:3])

def rms_error(pattern):
    return np.linalg.norm(pattern - np.ones(m * m)) / m

m = 25
n = 10
lamps = {
    1: [ 4.13, 4.63, 4.0],
    2: [14.11, 3.72, 3.5],
    3: [22.60, 7.92, 6.0],
    4: [5.53, 12.71, 4.0],
    5: [12.21, 15.30, 4.0],
    6: [15.31, 11.21, 6.0],
    7: [21.30, 14.50, 5.5],
    8: [3.93, 21.69, 5.0],
    9: [13.12, 20.69, 5.0],
    10: [20.30, 20.79, 4.5]
    }

x = [loc[0] for loc in lamps.values()]
y = [loc[1] for loc in lamps.values()]

power = np.zeros((n, m, m))

for i in range(m):
    for j in range(m):
        for l in lamps:
            power[l - 1, i, j] = 1 / np.linalg.norm(np.array([i - 0.5, j + 0.5, 0]) - np.array(lamps[l]))**2

power = power / power.sum(axis=0).mean()

A = np.array([power[i].reshape(m * m) for i in range(n)]).T

app_ui = ui.page_fluid(
    ui.layout_sidebar(
        ui.sidebar(
            ui.input_radio_buttons(  
                'power',  
                'Lighting strategy',
                {'0s': 'Use 0 vector',
                '1s': 'Use 1s vector',
                'ls': 'Use least squares',
                '1': 'Turn light 1 on only',
                '110': 'Turn lights 1 and 10 on only',
                'corners': 'Turn corner lights on only',
                '1-opt': 'Turn light 1 on only (optimally)',
                '110-opt': 'Turn lights 1 and 10 on only (optimally)',
                'corners-opt': 'Turn corner lights (optimally)'},
                selected='0s'
            )  
        ),
        ui.layout_columns(
        ui.output_plot("mainplot"),
        ui.card(
            ui.output_plot("histplot", height="45%"),
            ui.output_plot("powerplot", height="45%"),
            ui.output_text_verbatim("RMS"), 
        ),
        ),
    ),
)


def server(input, output, session):

    @reactive.calc
    def pattern():
        if input.power() == '0s':
            p = np.zeros(n)
        elif input.power() == '1s':
            p = np.ones(n)
        elif input.power() == '1':
            p = np.zeros(n)
            p[0] = 1
        elif input.power() == '110':
            p = np.zeros(n)
            p[0] = 1
            p[9] = 1
        elif input.power() == 'corners':
            p = np.zeros(n)
            p[0] = 1
            p[2] = 1
            p[7] = 1
            p[9] = 1
        elif input.power() == 'ls':
            p = np.linalg.pinv(A) @ np.ones(m * m)
        elif input.power() == '1-opt':
            p = np.zeros(n)
            p[0] = np.linalg.pinv(A[:, [0]]) @ np.ones(m * m)
        elif input.power() == '110-opt':
            theta = np.linalg.pinv(A[:, [0, 9]]) @ np.ones(m * m)
            p = np.zeros(n)
            p[0] = theta[0]
            p[9] = theta[1]
        elif input.power() == 'corners-opt':
            theta = np.linalg.pinv(A[:, [0, 2, 7, 9]]) @ np.ones(m * m)
            p = np.zeros(n)
            p[0] = theta[0]
            p[2] = theta[1]
            p[7] = theta[2]
            p[9] = theta[3]
        return get_pattern(power, p), p

    # ========================================================================

    @render.plot
    def mainplot():
        fig, ax = plt.subplots()
        im = ax.imshow(pattern()[0], cmap='grey', vmin=0, vmax=3)
        plt.colorbar(im, ax=ax)
        ax.axis('off')
        plt.scatter(x, y, c='red')
        for i, loc in lamps.items():
            ax.annotate(f'{i}: ({loc[2]:.1f}m)', (loc[0] - 0.35, loc[1] - 0.45), c='red')
        return fig

    # ========================================================================

    @render.plot
    def histplot():
        fig, ax = plt.subplots()
        ax.hist(pattern()[0].reshape(m * m), bins=50)
        ax.set_xlim(0, 2)
        ax.set_yticks([])
        ax.set_title('Histogram of luminosity')
        return fig

    # ========================================================================

    @render.plot
    def powerplot():
        fig, ax = plt.subplots()
        ax.bar([i + 1 for i in range(n)], pattern()[1])
        ax.set_xticks([i + 1 for i in range(n)])
        ax.set_ylim(0, max(3, np.ceil(np.max(pattern()[1]))))
        ax.set_title('Lamp powers')
        return 

    # ========================================================================

    @render.text
    def RMS():
        return f'The RMS error is {rms_error(pattern()[0].reshape(m * m)) :.2f}.'

    # ========================================================================

app = App(app_ui, server)

References

Boyd, S. P., & Vandenberghe, L. (2018). Introduction to Applied Linear Algebra : Vectors, Matrices, and Least Squares (First edition.). Cambridge University Press.