first draft even, odd etc period in grating

stimuli/components/__init__.py

@@ -1,6 +1,6 @@
 import itertools
 import warnings
-
+from copy import deepcopy
 import numpy as np
 
 from stimuli.utils import resolution
@@ -128,6 +128,7 @@ def mask_elements(
     # Mark elements with integer idx-value
     mask = np.zeros(base["shape"], dtype=int)
     for idx, edge in zip(reversed(range(len(edges))), reversed(edges)):
+        edge = np.round(edge, 10)
         mask[distances <= edge] = int(idx + 1)
 
     # Assemble output
@@ -139,6 +140,7 @@ def mask_elements(
         "shape": base["shape"],
         "visual_size": base["visual_size"],
         "ppd": base["ppd"],
+        "distances": distances,
     }
 
 
@@ -191,7 +193,7 @@ def resolve_grating_params(
         dictionary with all six resolution & size parameters resolved.
     """
 
-    if period not in ["ignore", "full", "half"]:
+    if period not in ["ignore", "even", "odd", "either"]:
         raise TypeError(f"period not understood: {period}")
 
     # Try to resolve resolution
@@ -204,16 +206,26 @@ def resolve_grating_params(
         length = length
         visual_angle = visual_angle
 
+    old_frequency = deepcopy(frequency)
+    old_n_phases = deepcopy(n_phases)
+    old_phase_width = deepcopy(phase_width)
+
     # Try to resolve number and width(s) of phases:
     # Logic here is that phase_width expresses "degrees per phase",
     # which we can convert to "phases_per_degree"
     # phase_width = degrees_per_phase = 1 / phases_per_degree = 1 / (2*frequency)
     if phase_width is not None:
+        # Make sure that phase width is realisable given ppd
+        phase_width = np.round(phase_width * ppd) / ppd
         phases_pd = 1 / phase_width
         if frequency is not None and phases_pd != 2 * frequency:
             raise ValueError(f"phase_width {phase_width} and frequency {frequency} don't match")
     elif frequency is not None:
+        # Make sure that frequency is realisable given ppd
         phases_pd = 2 * frequency
+        phase_width = np.round(1 / phases_pd * ppd) / ppd
+        phases_pd = 1 / phase_width
+        frequency = phases_pd / 2
     else:  # both are None:
         phases_pd = None
 
@@ -224,7 +236,7 @@ def resolve_grating_params(
     # pix = ppd * n_degrees
     # Thus we can resolve the number and spacing of phases also as a resolution
     try:
-        n_phases, min_angle, phases_pd = resolution.resolve_1D(
+        n_phases, min_angle, phases_pd_calc = resolution.resolve_1D(
             length=n_phases,
             visual_angle=visual_angle,
             ppd=phases_pd,
@@ -233,54 +245,84 @@ def resolve_grating_params(
     except Exception as e:
         raise Exception("Could not resolve grating frequency, phase_width, n_phases") from e
 
-    # Ensure full/half period?
-    if period != "ignore":
-        # Round n_phases
-        if period == "full":  # n_phases has to be even
-            n_phases = np.round(n_phases / 2) * 2
-        elif period == "half":  # n_phases can be odd
-            n_phases = np.round(n_phases)
-
-        # Check if n_phases fit in length
-        if length is not None and n_phases > 0 and length % n_phases:
-            raise resolution.ResolutionError(f"Cannot fit {n_phases} phases in {length} pix")
-
-        # Recalculate phases_pd
-        n_phases, min_angle, phases_pd = resolution.resolve_1D(
-            length=n_phases,
-            visual_angle=visual_angle,
-            ppd=None,
-            round=False,
-        )
-
-    # Convert to frequency
-    old_phase_width = phase_width
-    old_frequency = frequency
-    phase_width = 1 / phases_pd
+    phase_width = 1 / phases_pd_calc
+    phase_width = np.round(phase_width * ppd) / ppd
+    phases_pd = 1 / phase_width
     frequency = phases_pd / 2
-
-    if (old_phase_width is not None and phase_width != old_phase_width) or (
-        old_frequency is not None and frequency != old_frequency
-    ):
-        warnings.warn(
-            f"Adjusted frequency and phase width to ensure {period} period: {old_frequency} ->"
-            f" {frequency}, {old_phase_width} -> {phase_width}"
-        )
+
+    n_phases, min_angle, phases_pd_calc = resolution.resolve_1D(
+        length=None,
+        visual_angle=visual_angle,
+        ppd=phases_pd,
+        round=False,
+    )
 
     # Now resolve resolution
     visual_angle = min_angle if visual_angle is None else visual_angle
     length, visual_angle, ppd = resolution.resolve_1D(
         length=length, visual_angle=visual_angle, ppd=ppd
     )
 
+    # Ensure full/half period?
+    if period != "ignore":
+        # Round n_phases
+        if period == "even":  # n_phases has to be even
+            n_phases = np.round(n_phases / 2) * 2
+            phase_candidates = [n_phases, n_phases-2, n_phases+2, n_phases-4, n_phases+4, n_phases-6, n_phases+6]
+            phase_candidates = [n_phases]
+        elif period == "either":  # n_phases can be odd
+            n_phases = np.round(n_phases)
+            phase_candidates = [n_phases, n_phases-1, n_phases+1, n_phases-2, n_phases+2, n_phases-3, n_phases+3]
+        elif period == "odd":  # n_phases can be odd
+            n_phases = np.round(n_phases / 2) * 2 + 1
+            phase_candidates = [n_phases, n_phases-2, n_phases+2, n_phases-4, n_phases+4, n_phases-6, n_phases+6]
+
+        for pc in phase_candidates:
+            # Recalculate phases_pd
+            n_phases, min_angle, phases_pd = resolution.resolve_1D(
+                length=pc,
+                visual_angle=visual_angle,
+                ppd=None,
+                round=False,
+            )
+
+            # Convert to frequency
+            phase_width = 1 / phases_pd
+            frequency = phases_pd / 2
+
+            if not (phase_width*ppd) % 1:
+                break
+
+        if (phase_width*ppd) % 1:
+            raise resolution.ResolutionError(f"Cannot fit an {period} number of phases in {length} pix while "
+                                              "considering other parameters. Change parameters or increase ppd.")
+
+    if old_n_phases is not None:
+        if (n_phases < (old_n_phases-1)) or (n_phases > (old_n_phases+1)):
+            raise resolution.ResolutionError(f"Cannot fit {old_n_phases} phases in {length} pix."
+                                             " Change parameters or increase ppd. Best possible at "
+                                             f"given resolution: {n_phases} phases")
+
+        if n_phases != old_n_phases:
+            warnings.warn(
+                f"Adjusted n_phases={old_n_phases}->{n_phases} because original"
+                f"phase_width {old_phase_width} -> {phase_width} did not fit"
+            )
+
+    if (old_phase_width is not None and phase_width != old_phase_width):
+        warnings.warn(f"Adjusted phase width because of poor resolution: {old_phase_width} -> {phase_width}")
+
+    if (old_frequency is not None and frequency != old_frequency):
+        warnings.warn(f"Adjusted frequency because of poor resolution: {old_frequency} -> {frequency}")
+
     # Check that frequency does not exceed Nyquist limit:
     if frequency > (ppd / 2):
         raise ValueError(
             f"Grating frequency ({frequency}) should not exceed Nyquist limit {ppd/2} (ppd/2)"
         )
 
     # Accumulate edges of phases
-    edges = [*itertools.accumulate(itertools.repeat(phase_width, int(n_phases)))]
+    edges = [*itertools.accumulate(itertools.repeat(phase_width, int(np.ceil(n_phases))))]
     if "period" == "ignore":
         edges += [visual_angle]
 
@@ -290,7 +332,7 @@ def resolve_grating_params(
         "ppd": ppd,
         "frequency": frequency,
         "phase_width": phase_width,
-        "n_phases": int(n_phases),
+        "n_phases": n_phases,
         "edges": edges,
         "period": period,
     }

stimuli/components/grating.py

@@ -159,3 +159,48 @@ def square_wave(
         "n_bars": params["n_phases"],
         "period": params["period"],
     }
+
+
+if __name__ == "__main__":
+    from stimuli.utils.plotting import plot_stimuli
+
+    p1 = {
+          "visual_size": 5,
+          "ppd": 10,
+          "frequency": 2,
+          "period": "even",
+          }
+
+    p2 = {
+          "visual_size": 5,
+          "ppd": 20,
+          "n_bars": 20,
+          "period": "even",
+          }
+
+    p3 = {
+          "visual_size": 5,
+          "ppd": 10,
+          "frequency": 5,
+          "period": "even",
+          }
+
+    p4 = {
+          "visual_size": 5,
+          "ppd": 40,
+          "frequency": 1.2,
+          "period": "ignore",
+          }
+
+    sw1 = square_wave(**p1)
+    sw2 = square_wave(**p2)
+    sw3 = square_wave(**p3)
+    sw4 = square_wave(**p4)
+
+    stims = {
+        "sw1": sw1,
+        "sw2": sw2,
+        "sw3": sw3,
+        "sw4": sw4,
+        }
+    plot_stimuli(stims)