<jupyter_start><jupyter_text>Stimulation devices===================To stimulate one class of retinal photoreceptor without changing the activation in others&mdash;which is to perform the method of silent substitution&mdash;requires a multiprimary stimulation device with at least as many primaries as there are photoreceptors in the retina. This generally means that [5 primaries are needed](./01_background.rst), although 4 primaries may suffice when working in the photopic range as rod photoreceptors become saturated and incapable of signalling above 300 cd/m$^2$ ([<NAME>, 1954](https://doi.org/10.1080/713818657); [Adelson, 1982](https://doi.org/10.1016/0042-6989(82)90143-2); but see [Shapiro, 2002](https://iovs.arvojournals.org/article.aspx?articleid=2200107); [Kremers et al., 2009](https://doi.org/10.1007/s10633-008-9159-0)). The primaries should be independantly addressable, additive, and ideally stable over time with a linear gamma function. Peak wavelength and bandwidth are key considerations that will ultimately define the gamut and available contrast ([Evéquoz et al., 2021](https://doi.org/10.1364/josaa.420373)), and the light source will also need to be integrated into an optical setup for stimulus delivery&mdash;usually either a Ganzfeld (e.g., [Martin et al., 2021](https://doi.org/10.3758/s13428-021-01759-3)) or Maxwellian (e.g., [Cao et al., 2015](https://doi.org/10.1167/15.1.27)).[Conus and Geiser (2020)](https://doi.org/10.3390/photonics7040121) reviewed stimulation devices from a range of silent substitution studies and found that in most cases the device had 4 or 5 primaries and was built from scratch using LEDs, optical bench components, and microprocessors, such as Arduino, for pulse width modulation control of intensity. Only a few devices were commercially bought. Whatever the device and setup, it will be necessary to have a forward calibration model that is representative of what an observer sees, and which can be used to predict spectral output for any combination of settings. Typically this involves sampling each of the primaries at a range of itentsities with an external spectrometer and using interpolation methods.  `pysilsub.device.StimulationDevice`-----------------------------------*PySilSub* has a generic `StimulationDevice` class that will turn spectral calibration data for any multiprimary system into a predictive model. To get started, you will need a calibration file containing spectral measurements. Each row should be a spectral measurement for a device primary at a given setting, and there should be column labels 'Primary' and 'Setting' to indicate this. Note that it will be important to include a dark measurement where LEDs are turned off (i.e., Setting = 0) to account for ambient light in an experimental setup. If you have a [JETI](https://www.jeti.com/) or [OceanOptics](https://www.oceaninsight.com/) spectrometer, *[PyPlr](https://github.com/PyPlr/cvd_pupillometry)* has some Python interfaces that may help you obtain the spectral measurements.The following demonstration uses calibration data from a linear 5-primary device with 8-bit resolution depth, which is to say that the intensity of each LED is set with a number between 0 and 255. The full range of intensities was sampled in steps of 8, so there are 33 measurements for each LED. The first row of the calibration file contains the wavelengths of the spectral measurements as well as columns to identify the primary and setting.<jupyter_code>import pandas as pd

data = pd.read_csv(
    '../../data/BCGAR_5_Primary_8_bit_linear.csv', 
    index_col=['Primary', 'Setting']
)
data.columns = data.columns.astype('int64')
data.columns.name = 'Wavelength'
data<jupyter_output><empty_output><jupyter_text>After reading the CSV file into a `DataFrame` using pandas, we can pass the data to `StimulationDevice` along with some further information to gain access to many useful features.<jupyter_code>from pysilsub.device import StimulationDevice

device = StimulationDevice(
    resolutions=[255]*5,
    colors=['blue', 'cyan', 'green', 'orange', 'red'],
    spds=data,
    spd_binwidth=1,
    name='BCGAR (8-bit, linear)'
)<jupyter_output><empty_output><jupyter_text>For example, plot the spectral power distributions of all measurements.<jupyter_code>spd_fig = device.plot_spds()<jupyter_output><empty_output><jupyter_text>Or the gamut of the device, showing the boundaries of the colour space that can be achieved on the [CIE *xy* chromaticity horseshoe](https://en.wikipedia.org/wiki/CIE_1931_color_space).<jupyter_code>gamut_fig = device.plot_gamut()<jupyter_output><empty_output><jupyter_text>More useful, we can use `device.predict_primary_spd(...)` to predict the spectral output of the stimulation device for a chosen primary at a given setting. This method simply interpolates between the measured spectra that were provided in the calibration data. It is a workhorse method that gets called repeatedly.<jupyter_code>primary_spd = device.predict_primary_spd(
    primary=2, 
    setting=.5,  # Could also say 128
    name='Primary 2 (half power)'
)
print(primary_spd)
primary_spd.plot(ylabel='W/m$^2$/nm', c=device.colors[2]);<jupyter_output>Wavelength
380    0.001115
381    0.001191
382    0.001198
383    0.001231
384    0.001003
         ...   
776    0.000710
777    0.000693
778    0.000562
779    0.000733
780    0.000826
Name: Primary 2 (half power), Length: 401, dtype: float64<jupyter_text>We can also use `device.predict_multiprimary_spd(...)` to predict the spectral output for a combination of primaries and settings. This method calls the previous method for each primary/setting pair and sums the spectral power distributions.<jupyter_code>half_max = [.5, .5, .5, .5, .5]  # Same as [128, 128, 128, 128, 128]
device_spd = device.predict_multiprimary_spd(half_max, 'Half max')
device_spd.plot(legend=True, ylabel='W/m$^2$/nm', color='k')
print(f'Predicted output for device settings: {half_max}')<jupyter_output>Predicted output for device settings: [0.5, 0.5, 0.5, 0.5, 0.5]<jupyter_text>`device.predict_multiprimary_aopic(...)` will return the photoreceptor-specific irradiances for a spectrum produced by any combination of settings. This method calls the previous method and then weights the spectrum by the photoreceptor spectral sensitivities, before summing up the values for each photoreceptor.<jupyter_code>device_ao = device.predict_multiprimary_aopic(half_max, 'Half max')
ao_colors = list(device.aopic_colors.values())
device_ao.plot(kind='bar', color=ao_colors, ylabel='W/m$^2$')
device_ao<jupyter_output><empty_output><jupyter_text>We can also search for a spectrum based on *xy* chromaticity and luminance.<jupyter_code>target_xy=[0.31271, 0.32902]
target_luminance=600.
bg = device.find_settings_xyY(target_xy, target_luminance)<jupyter_output>basinhopping step 0: f 1.11019e-11
basinhopping step 1: f 5.7635e-08 trial_f 5.7635e-08 accepted 1  lowest_f 1.11019e-11
Requested LMS: [ 0.87374046  0.78105394  0.44556496]
Solution LMS: [ 0.87373816  0.78105248  0.44556688]
{'fraction_non_alphanumeric': 0.0954745353950915, 'fraction_numerical': 0.06142715278762945, 'mean_word_length': 5.157205240174672, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 12, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '364498', 'license': 'MIT', 'n_tokens_mistral': 2404, 'n_tokens_neox': 2038, 'n_words': 783, 'path': 'docs/source/03a_stimulation_devices.ipynb', 'repo_name': 'PySilentSubstitution/silent-sub', 'size': '7062'}