Simulate spectra (take 2)

[thomased]

Fix https://github.com/rmaia/pavo/issues/151. Moving discussion here from #185 because I started from scratch after completely forgetting that you'd already had a go at it @Bisaloo, then tried to merge them & screwed everything up, so just started again here (I'll learn git one day...).

Stuff to to discuss (partly from #185):

• [X] Should it just be a function in pavo, or would there value in thinking of a little package dedicated to simulating 'natural' spectra at this early stage? Could also split off and expand sensmodel() from pavo, and explore possibility of functionality for radiance and/or irradiance spectra (trickier...). So it could be a package for simulating 'biological' reflectance, absorbance/absorbtance, radiance (maybe), or irradiance spectra? Alternatively, could add a few new related functions to pavo. Like sim.reflect(), sim.absorb(), sim.irrad() or something. A: Keep in pavo for immediacy & convenience benefits (& lack of strong reason not to at this stage), though reserve the right to roll into its own package for fun later on.
• [X] Argument names? A: Use snake_case for argument names. Snake has always been my preference anyway - easy. 711f77c
• [x] Function name? simulate_spec() or sim_spec() more consistent with current norms? I'd reach for the former since it's not that much longer, but is more clear. Should it be simulate_reflect() in case of future simulate_irrad and the like though? Or I guess you could use this to simulate irradiance spectra of a sort already - maybe it's more a case of fleshing out the examples to show as much. In which case simulate_spec is sensible. Added an irradiance example to the docs - I'm leaning toward keeping simulate_spec as the name.
• [X] Currently it's just simple curve-generation. Should it also implement physical models? Thin-film, multilayer structures etc. A: Agree it's too big a can of worms to include physical models etc. Let's stick to simple sims for now, which will cover a good set of use cases for users as-is.
• [X] Offer the option to simulate several spectra with different characteristics at once? I'm inclined not to, and just have people lapply() it or whatever. A: No, leave it to users, but added an example. ee0f7a9f7a55dd229f2cfbabdcc435147555620b
• [X] Offer the option to add noise? A: No. I didn't actually think of the potential for fakery, but fair point. Easy enough for users to add if they want anyway - not worth an additional argument at this stage in lieu of clear use case.
• [X] Return 'ideal' (flat) reflectance when no options specified? A: Yes, for now (still open to alternative/error). Returns an ideal reflector with the max as defined by ylim. Added example to support. 5b43030f0624fa4523548d03d6134241b02aeb82
• [X] Multiple peaks/inflections are currently combined additively, then normalised to ylim() range. This is fine, and does what you'd expect when the changes in reflectance are separated (like a bell curve in the UV, and sigmoid in the red, in the example below). But when changes are close together, the resulting curve can be hard to predict, and the 'shape' of the underlying constituent curves isn't preserved. Various alternatives, but then the y range is often blown out. I think it's ultimately a trade-off, but need to think about & understand it better. A: Agree I think current behaviour is what you'd most likely expect. May revise on further testing/playing, but happy for now.
• [x] Include ylim()? A: I'm currently inclined to keep it (still open to removing though). It'd make simulations of luminance/brightness variation simpler. i.e. users could lapply() through a set of ylim()s to simulate certain kinds of limited-view iridescence, without having to manually renormalise along the way. The direction of any sigmoidal curve is also now controlled by this, as is the maximum of the (default) 'ideal' spectrum, so I'm increasingly inclined to keep it. Added an example - 2fc8a4966af84299cf4c4107bcf9376c6222916c
• [x] More informative spectrum names (like sensmodel() does)? A: Added unique naming convention, which includes information on specified inflection points and/or peaks. 3cfca35c0166006aed168818d7c08c20f9bc478c.
• [X] Add ability to generate (biologically silly?) high-to-low sigmoidal curves? Perhaps controlled by ylim(), such that ylim = c(0, 100), wl_inflect = 550 would generate a low-to-high curve, whereas ylim = c(100, 0), wl_inflect = 550 would generate a high-to-low? A: Can't see any harm? Done here - 494c4f2e0aca0ac68f1c3b7a1bde54c8035b320b
• [x] Add further flexibility for simulating non-reflectance spectra (irradiance, radiance etc.)? What other parameters might be useful to control?
  • A: Added skew_gauss parameter so Gaussian curves are drawn from skew-normal distribution. Useful for simulating convincing irrad/rad spectra - https://github.com/rmaia/pavo/pull/247/commits/7ba9c127dfb7f35502317f00ddd9de9c5646dab1
• [X] Add tests A: Roughed in - 95758311683d13a03300a4daa16f8357d2a65f13
• [X] Add news - https://github.com/rmaia/pavo/pull/247/commits/5088504c3285e2f3461ded430f2495766b787e46
• [ ] Add to handbook (will do once finalised & merged)

Thoughts welcome!

Working example:

# Single sigmoidal spectrum, with an inflection at 500 nm.
spec <- as.rspec(simulate_spec2(wl.inflect = 500))

# Single Gaussian spectrum, with a peak at 400 nm
spec2 <- as.rspec(simulate_spec2(wl.peak = 400))

# Combination of both Gaussian (with peak at 340 nm) and sigmoidal (with inflection
# at 560 nm)
spec3 <- as.rspec(simulate_spec2(wl.peak = 340, wl.inflect = 560))

# Double-Gaussian peaks of differing widths
spec4 <- as.rspec(simulate_spec2(wl.peak = c(340, 560), width.gauss = c(12, 40)))

# Complex spectrim with multi-Gaussian and single sigmoidal peak
spec5 <- as.rspec(simulate_spec2(wl.peak = c(340, 430), width.gauss = c(20, 60), wl.inflect = 575))

# Merge and plot
spectra <- merge(spec, spec2)
spectra <- merge(spectra, spec3)
spectra <- merge(spectra, spec4)
spectra <- merge(spectra, spec5)
explorespec(spectra, ylim = c(0, 100))

[thomased]

Awesome, many thanks for those edits & your thoughts. I'm with you on all points — I'll keep adjusting things along those lines and we can see how we like it.

[thomased]

Latest examples, from the help docs:

# Single ideal 'grey' reflectance spectrum, with 50% reflectance across 300 - 700 nm.
reflect0 <- simulate_spec(ylim = c(0, 50))

# Single sigmoidal spectrum, with a low-to-high inflection at 550 nm.
reflect1 <- simulate_spec(wl_inflect = 550)

# Single Gaussian spectrum, with a peak at 400 nm
reflect2 <- simulate_spec(wl_peak = 400)

# Combination of both Gaussian (with peak at 340 nm) and sigmoidal (with inflection
# at 560 nm)
reflect3 <- simulate_spec(wl_inflect = 560, wl_peak = 340)

# Double-Gaussian peaks of differing widths
reflect4 <- simulate_spec(wl_peak = c(340, 560), width_gauss = c(12, 40))

# Complex spectrum with single sigmoidal peak and multi-Gaussian peaks
reflect5 <- simulate_spec(wl_inflect = 575, wl_peak = c(340, 430), width_gauss = c(20, 60))

# Merge and plot
spectra <- Reduce(merge, list(reflect0, reflect1, reflect2, reflect3, reflect4, reflect5))
plot(spectra, lty = 1:5)

---

# Simulate a set of Gaussian reflectance curves with peaks varying between 400-600nm
# in increments of 10, then combine into a single rspec object, and plot the result
peaks <- seq(400, 600, 10)  # Peak locations
reflectances <- lapply(seq_along(peaks), function(x) simulate_spec(wl_peak = peaks[x]))  # Simulate
reflectances <- Reduce(merge, reflectances)  # Combine
plot(reflectances)  # Plot

---

# Simulate a set of Gaussian reflectance curves with a single peak at 500 nm, but
# with maximum reflectance varying from 30 to 90% in 10% increments, then combine
# into a single rspec object, and plot the result
ymax <- lapply(seq(30, 90, 10), function(x) c(0, x))  # Varying reflectance maxima
reflectances2 <- lapply(ymax, function(x) simulate_spec(wl_peak = 500, ylim = x))  # Simulate
reflectances2 <- Reduce(merge, reflectances2)  # Combine
plot(reflectances2)  # Plot

---

# To simulate non-reflectance spectra (like irradiances or radiances), it's often useful
# to explore more 'extreme' parameters. Here's a simple example which simulates
# natural daylight, as represented by the D65 standard daylight spectrum.
D65_real <- procspec(sensdata(illum = 'D65'), opt = 'smooth')  # Official D65 daylight spectrum
D65_sim <- simulate_spec(wl_peak = 400,
                         width_gauss = 1300,
                         skew_gauss = 10,
                         ylim = c(0, 1))  # Simulated D65
cor.test(D65_real$D65, D65_sim$spec_p400)  # >0.99 correlation
plot(merge(D65_real, D65_sim), lty = 1:2, ylab = 'Irradiance (%)')  # Merge and plot the two spectra

[thomased]

Absolutely no rush at all @Bisaloo, but if you want to cast another eye over it sometime it'd be appreciated. I think I'm basically happy with where it's at - just about ready to merge. Can always do some further fine-tuning before the next release of course. And happy to hear alternative views on any of the above points - I've just jotted down my thinking.