Update based on review

examples/tutorial_3_LUPA.ipynb

@@ -1000,7 +1000,7 @@
     "    fp = 1/tp\n",
     "    spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
     "    efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
-    "    return wot.waves.long_crested_wave(efth,nrealizations=5)\n",
+    "    return wot.waves.long_crested_wave(efth,nrealizations)\n",
     "\n",
     "for case, data in wave_cases.items():\n",
     "    waves[case] = irregular_wave(data['Hs'], data['Tp'])"

examples/tutorial_4_Pioneer.ipynb

@@ -90,7 +90,7 @@
     "fp = 1/Tp\n",
     "spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, Hs)\n",
     "efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
-    "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations = nrealizations)"
+    "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations)"
    ]
   },
   {

tests/test_waves.py

@@ -190,7 +190,7 @@ def elevation(self, ndbc_omnidirectional, direction, nrealizations):
         """Complex sea state elevation amplitude [m] indexed by
         frequency and direction."""
         elev = wot.waves.long_crested_wave(
-            ndbc_omnidirectional.efth, direction, nrealizations)
+            ndbc_omnidirectional.efth, nrealizations, direction)
         return elev
 
     @pytest.fixture(scope="class")
@@ -244,7 +244,7 @@ def test_direction(self, elevation, direction):
         dir_out = elevation.wave_direction.values.item()
         assert np.isclose(dir_out, wot.degrees_to_radians(direction))
 
-    def test_realizations(self, elevation, direction):
+    def test_realizations(self, elevation):
         """Test that the number of realizations is correct."""
         realization_out = elevation.realization.values
         assert (realization_out == [0,1]).all()
@@ -259,7 +259,7 @@ def test_time_series(self, pm_spectrum, pm_f1, pm_nfreq):
         # create time-series
         direction = 0.0
         nrealizations = 1
-        wave = wot.waves.long_crested_wave(pm_spectrum, direction, nrealizations)
+        wave = wot.waves.long_crested_wave(pm_spectrum, nrealizations, direction)
         wave_ts = wot.fd_to_td(wave.sel(realization=0).values, pm_f1, pm_nfreq, False)
         # calculate the spectrum from the time-series
         t = wot.time(pm_f1, pm_nfreq)
@@ -289,12 +289,17 @@ def ndbc_spectrum(self,):
         files = [f'41013{i}2020.txt' for i in markers]
         spec = ws.read_ndbc_ascii([os.path.join(dir, file) for file in files])
         return spec.sel(time=time).interp(freq=freq)
+
+    @pytest.fixture(scope="class")
+    def nrealizations(self):
+        """Number of wave realizations."""
+        return 2
 
     @pytest.fixture(scope="class")
-    def elevation(self, ndbc_spectrum):
+    def elevation(self, ndbc_spectrum, nrealizations):
         """Complex sea state elevation amplitude [m] indexed by
         frequency and direction."""
-        return wot.waves.irregular_wave(ndbc_spectrum.efth)
+        return wot.waves.irregular_wave(ndbc_spectrum.efth, nrealizations)
 
     def test_coordinates(self, elevation):
         """Test that the elevation dataArray has the correct
@@ -304,15 +309,20 @@ def test_coordinates(self, elevation):
         for icoord in coordinates:
             assert icoord in elevation.coords, f'missing coordinate: {icoord}'
 
-    def test_shape(self, ndbc_spectrum, elevation):
+    def test_shape(self, ndbc_spectrum, elevation, nrealizations):
         """Test that the elevation dataArray has the correct shape."""
         nfreq = len(ndbc_spectrum.freq)
         ndir = len(ndbc_spectrum.dir)
-        assert np.squeeze(elevation.values).shape == (nfreq, ndir)
+        assert np.squeeze(elevation.values).shape == (nfreq, ndir, nrealizations)
 
     def test_type(self, elevation):
         """Test that the elevation dataArray has the correct type."""
         assert np.iscomplexobj(elevation)
+
+    def test_realizations(self, elevation):
+        """Test that the number of realizations is correct."""
+        realization_out = elevation.realization.values
+        assert (realization_out == [0,1]).all()
 
 
 class TestRandomPhase:
@@ -323,7 +333,8 @@ def shape(self,):
         """Shape of the phase matrix, randomized each time the test is
         run.
         """
-        return (np.random.randint(10, 100), np.random.randint(10, 100))
+        return (np.random.randint(10, 100), np.random.randint(10, 100),
+                np.random.randint(10, 100))
 
     @pytest.fixture(scope="class")
     def phase_mat(self, shape):

wecopttool/core.py

@@ -683,11 +683,13 @@ def solve(self,
                                 obj_fun=pto.average_power,
                                 nstate_opt=2*nfreq+1)
 
-        To get the post-processed results for the
-        :py:class:`wecopttool.pto.PTO`, you may call
+        To get the post-processed results for the :py:class:`wecopttool.WEC`
+        and :py:class:`wecopttool.pto.PTO` for a single realization, you
+        may call
 
-        >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt[0])
-        >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[0])
+        >>> realization = 0 # realization index
+        >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt[realization])
+        >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[realization])
 
         See Also
         --------
@@ -712,41 +714,13 @@ def solve(self,
 
         # composite scaling vector
         scale = np.concatenate([scale_x_wec, scale_x_opt])
-
-        # objective function
-        sign = -1.0 if maximize else 1.0
-
-        def obj_fun_scaled(x):
-            x_wec, x_opt = self.decompose_state(x/scale)
-            return obj_fun(self, x_wec, x_opt, wave)*scale_obj*sign
-
-        # constraints
-        constraints = self.constraints.copy()
-
-        for i, icons in enumerate(self.constraints):
-            icons_new = {"type": icons["type"]}
-
-            def make_new_fun(icons):
-                def new_fun(x):
-                    x_wec, x_opt = self.decompose_state(x/scale)
-                    return icons["fun"](self, x_wec, x_opt, wave)
-                return new_fun
-
-            icons_new["fun"] = make_new_fun(icons)
-            if use_grad:
-                icons_new['jac'] = jacobian(icons_new['fun'])
-            constraints[i] = icons_new
-
-        # system dynamics through equality constraint, ma - Σf = 0
-        def scaled_resid_fun(x):
-            x_s = x/scale
-            x_wec, x_opt = self.decompose_state(x_s)
-            return self.residual(x_wec, x_opt, wave)
-
-        eq_cons = {'type': 'eq', 'fun': scaled_resid_fun}
-        if use_grad:
-            eq_cons['jac'] = jacobian(scaled_resid_fun)
-        constraints.append(eq_cons)
+
+        # decision variable initial guess
+        if x_wec_0 is None:
+            x_wec_0 = np.random.randn(self.nstate_wec)
+        if x_opt_0 is None:
+            x_opt_0 = np.random.randn(nstate_opt)
+        x0 = np.concatenate([x_wec_0, x_opt_0])*scale 
 
         # bounds
         if (bounds_wec is None) and (bounds_opt is None):
@@ -771,32 +745,60 @@ def scaled_resid_fun(x):
             bounds = Bounds(lb=np.hstack([le.lb for le in bounds_list])*scale,
                             ub=np.hstack([le.ub for le in bounds_list])*scale)
 
-        # callback
-        if callback is None:
-            def callback_scipy(x):
-                x_wec, x_opt = self.decompose_state(x)
-                max_x_opt = np.nan if np.size(x_opt)==0 else np.max(np.abs(x_opt))
-                _log.info("Scaled [max(x_wec), max(x_opt), obj_fun(x)]: "
-                          + f"[{np.max(np.abs(x_wec)):.2e}, "
-                          + f"{max_x_opt:.2e}, "
-                          + f"{obj_fun_scaled(x):.2e}]")
-        else:
-            def callback_scipy(x):
-                x_s = x/scale
-                x_wec, x_opt = self.decompose_state(x_s)
-                return callback(self, x_wec, x_opt, wave)
-
         for realization, wave in waves.groupby('realization'):
 
             _log.info("Solving pseudo-spectral control problem "
                       + f"for realization number {realization}.") 
+
+            # objective function
+            sign = -1.0 if maximize else 1.0
+
+            def obj_fun_scaled(x):
+                x_wec, x_opt = self.decompose_state(x/scale)
+                return obj_fun(self, x_wec, x_opt, wave)*scale_obj*sign
+
+            # constraints
+            constraints = self.constraints.copy()
+
+            for i, icons in enumerate(self.constraints):
+                icons_new = {"type": icons["type"]}
+
+                def make_new_fun(icons):
+                    def new_fun(x):
+                        x_wec, x_opt = self.decompose_state(x/scale)
+                        return icons["fun"](self, x_wec, x_opt, wave)
+                    return new_fun
 
-            # decision variable initial guess
-            if x_wec_0 is None:
-                x_wec_0 = np.random.randn(self.nstate_wec)
-            if x_opt_0 is None:
-                x_opt_0 = np.random.randn(nstate_opt)
-            x0 = np.concatenate([x_wec_0, x_opt_0])*scale 
+                icons_new["fun"] = make_new_fun(icons)
+                if use_grad:
+                    icons_new['jac'] = jacobian(icons_new['fun'])
+                constraints[i] = icons_new
+
+            # system dynamics through equality constraint, ma - Σf = 0
+            def scaled_resid_fun(x):
+                x_s = x/scale
+                x_wec, x_opt = self.decompose_state(x_s)
+                return self.residual(x_wec, x_opt, wave)
+
+            eq_cons = {'type': 'eq', 'fun': scaled_resid_fun}
+            if use_grad:
+                eq_cons['jac'] = jacobian(scaled_resid_fun)
+            constraints.append(eq_cons)
+
+            # callback
+            if callback is None:
+                def callback_scipy(x):
+                    x_wec, x_opt = self.decompose_state(x)
+                    max_x_opt = np.nan if np.size(x_opt)==0 else np.max(np.abs(x_opt))
+                    _log.info("Scaled [max(x_wec), max(x_opt), obj_fun(x)]: "
+                              + f"[{np.max(np.abs(x_wec)):.2e}, "
+                              + f"{max_x_opt:.2e}, "
+                              + f"{obj_fun_scaled(x):.2e}]")
+            else:
+                def callback_scipy(x):
+                    x_s = x/scale
+                    x_wec, x_opt = self.decompose_state(x_s)
+                    return callback(self, x_wec, x_opt, wave)
 
             # optimization problem
             optim_options['disp'] = optim_options.get('disp', True)

wecopttool/waves.py

@@ -95,7 +95,7 @@ def elevation_fd(
     freq = frequency(f1, nfreq, False)
     omega = freq*2*np.pi
 
-    dims = ('omega', 'wave_direction','realization')
+    dims = ('omega', 'wave_direction', 'realization')
     omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'}
     freq_attr = {'long_name': 'Frequency', 'units': 'Hz'}
     dir_attr = {'long_name': 'Wave direction', 'units': 'rad'}
@@ -107,11 +107,17 @@ def elevation_fd(
 
     if amplitudes is None:
         amplitudes = np.zeros([nfreq, ndirections, nrealizations])
+    else:
+        if amplitudes.shape == (nfreq, ndirections):
+            amplitudes = np.expand_dims(amplitudes,axis=2)
+        assert amplitudes.shape == (nfreq, ndirections, 1) or \
+                amplitudes.shape == (nfreq, ndirections, nrealizations)
 
     if phases is None:
-        phases = random_phase([nfreq, ndirections, nrealizations],seed)
+        phases = random_phase([nfreq, ndirections, nrealizations], seed)
     else:
         phases = degrees_to_radians(phases, False)
+        assert phases.shape == (nfreq, ndirections, nrealizations)
 
     camplitude = amplitudes * np.exp(1j*phases)
 
@@ -190,8 +196,8 @@ def regular_wave(
 
 def long_crested_wave(
     efth: DataArray,
+    nrealizations: int,
     direction: Optional[float] = 0.0,
-    nrealizations: Optional[float] = 1,
     seed: Optional[float] = None,
 ) -> DataArray:
     """Create a complex frequency-domain wave elevation from an
@@ -210,11 +216,11 @@ def long_crested_wave(
     efth
         Omnidirection wave spectrum in units of m^2/Hz, in the format
         used by :py:class:`wavespectra.SpecArray`.
-    direction
-        Direction (in degrees) of the long-crested wave.
     nrealizations
         Number of wave phase realizations to be created for the 
         long-crested wave.
+    direction
+        Direction (in degrees) of the long-crested wave.
     seed
         Seed for random number generator. Used for reproducibility.
         Generally should not be used except for testing.
@@ -225,7 +231,6 @@ def long_crested_wave(
     values = efth.values
     values[values<0] = np.nan
     amplitudes = np.sqrt(2 * values * df)
-    amplitudes = np.expand_dims(amplitudes,axis=2)
 
     attr = {
         'Wave type': 'Long-crested irregular',
@@ -236,7 +241,7 @@ def long_crested_wave(
 
 
 def irregular_wave(efth: DataArray,
-                   nrealizations: Optional[float] = 1,
+                   nrealizations: int,
                    seed: Optional[float] = None,) -> DataArray:
     """Create a complex frequency-domain wave elevation from a spectrum.
 
@@ -270,7 +275,6 @@ def irregular_wave(efth: DataArray,
     values = efth.values
     values[values<0] = np.nan
     amplitudes = np.sqrt(2 * values * df * dd)
-    amplitudes = np.expand_dims(amplitudes,axis=2)
 
     attr = {'Wave type': 'Irregular'}