Merge dev

README.md

@@ -37,13 +37,13 @@ The library provides:
 
 First, create a conda environment (skip if you already have one):
 ```
-conda create --name opt python=3.11
+conda create --name opt python=3.12
 conda activate opt
 ```
 
 Then install the package:
 ```
-conda install "numpy<2.0"
+conda install numpy
 pip install sb-arch-opt
 ```
 

docs/index.md

@@ -31,13 +31,13 @@ The library provides:
 
 First, create a conda environment (skip if you already have one):
 ```
-conda create --name opt python=3.11
+conda create --name opt python=3.12
 conda activate opt
 ```
 
 Then install the package:
 ```
-conda install "numpy<2.0"
+conda install numpy
 pip install sb-arch-opt
 ```
 

sb_arch_opt/__init__.py

@@ -1 +1 @@
-__version__ = '1.6.1'
+__version__ = '1.6.3'

sb_arch_opt/algo/arch_sbo/infill.py

@@ -283,7 +283,7 @@ def _g(_):
             else:
                 x_optimized.append(x_ref_i)
 
-        return Population.new(X=np.row_stack(x_optimized))
+        return Population.new(X=np.vstack(x_optimized))
 
     @staticmethod
     def get_pareto_front(f: np.ndarray) -> np.ndarray:

sb_arch_opt/algo/egor_interface/algo.py

@@ -217,9 +217,9 @@ def _run_infills(self, n_infills: int):
             x, x_failed, y = self._get_xy(pop)
 
             # Update
-            self._x = np.row_stack([self._x, x])
-            self._y = np.row_stack([self._y, y])
-            self._x_failed = np.row_stack([self._x_failed, x_failed])
+            self._x = np.vstack([self._x, x])
+            self._y = np.vstack([self._y, y])
+            self._x_failed = np.vstack([self._x_failed, x_failed])
 
             # Store results
             if self._results_folder is not None:

sb_arch_opt/algo/segomoe_interface/algo.py

@@ -316,9 +316,9 @@ def _dummy_f_grouped(_):
         return sego.get_x(i=-1)
 
     def _tell_infill(self, x, x_failed, y):
-        self._x = np.row_stack([self._x, x]) if self._x is not None else x
-        self._y = np.row_stack([self._y, y]) if self._y is not None else y
-        self._x_failed = np.row_stack([self._x_failed, x_failed]) if self._x_failed is not None else x_failed
+        self._x = np.vstack([self._x, x]) if self._x is not None else x
+        self._y = np.vstack([self._y, y]) if self._y is not None else y
+        self._x_failed = np.vstack([self._x_failed, x_failed]) if self._x_failed is not None else x_failed
         self._save_results()
 
     def _get_sego(self, f_grouped):
@@ -465,10 +465,10 @@ def get_population(self, x: np.ndarray, y: np.ndarray, x_failed: np.ndarray = No
         f, g, h = self._split_y(y)
 
         if x_failed is not None and len(x_failed) > 0:
-            x = np.row_stack([x, x_failed])
-            f = np.row_stack([f, np.zeros((x_failed.shape[0], f.shape[1]))*np.inf])
-            g = np.row_stack([g, np.zeros((x_failed.shape[0], g.shape[1]))*np.inf])
-            h = np.row_stack([h, np.zeros((x_failed.shape[0], h.shape[1]))*np.inf])
+            x = np.vstack([x, x_failed])
+            f = np.vstack([f, np.zeros((x_failed.shape[0], f.shape[1]))*np.inf])
+            g = np.vstack([g, np.zeros((x_failed.shape[0], g.shape[1]))*np.inf])
+            h = np.vstack([h, np.zeros((x_failed.shape[0], h.shape[1]))*np.inf])
 
         kwargs = {'X': x, 'F': f, 'G': g, 'H': h}
         pop = Population.new(**kwargs)

sb_arch_opt/algo/tpe_interface/api.py

@@ -170,7 +170,7 @@ def __init__(self):
         super().__init__(sampling=None)
 
     def do(self, problem, n_samples, **kwargs):
-        x_init = np.row_stack([self.interface.ask_init() for _ in range(n_samples)])
+        x_init = np.vstack([self.interface.ask_init() for _ in range(n_samples)])
         return Population.new(X=x_init)
 
 

sb_arch_opt/design_space.py

@@ -647,8 +647,8 @@ def _get_all_discrete_x_by_trial_and_imputation(self):
                 x_repair = x_repair[is_not_repaired, :]
                 is_active = is_active[is_not_repaired, :]
 
-                x_discr = np.row_stack([x_discr, x_repair])
-                is_act_discr = np.row_stack([is_act_discr, is_active.astype(bool)])
+                x_discr = np.vstack([x_discr, x_repair])
+                is_act_discr = np.vstack([is_act_discr, is_active.astype(bool)])
 
             # Impute continuous values
             self.impute_x(x_discr, is_act_discr)

sb_arch_opt/pareto_front.py

@@ -140,9 +140,9 @@ def _calc_pareto_set_front(self, *_, pop_size=None, n_gen_min=10, n_repeat=4, n_
                         ps = res.X
                         pf = res.F
                     else:
-                        pf_merged = np.row_stack([pf, res.F])
+                        pf_merged = np.vstack([pf, res.F])
                         i_non_dom = NonDominatedSorting().do(pf_merged, only_non_dominated_front=True)
-                        ps = np.row_stack([ps, res.X])[i_non_dom, :]
+                        ps = np.vstack([ps, res.X])[i_non_dom, :]
                         pf = pf_merged[i_non_dom, :]
 
         # Reduce size of Pareto front to a predetermined amount to ease Pareto-front-related calculations

sb_arch_opt/problems/gnc.py

@@ -175,8 +175,8 @@ def _gen_all_discrete_x(self) -> Optional[Tuple[np.ndarray, np.ndarray]]:
             x_rows.append(x_combs)
             is_active_rows.append(is_act_combs)
 
-        x_all = np.row_stack(x_rows)
-        is_active_all = np.row_stack(is_active_rows)
+        x_all = np.vstack(x_rows)
+        is_active_all = np.vstack(is_active_rows)
         return x_all, is_active_all
 
     def _get_discrete_x_combs_type(self, x_base, j, n_objs):
@@ -234,7 +234,7 @@ def _iter_conns(n_src_, n_tgt_):
                     n_combinations += 1
 
             if return_conns:
-                n_combinations = np.row_stack(n_combinations)
+                n_combinations = np.vstack(n_combinations)
             n_comb_conn[n_src, n_tgt] = n_comb_conn[n_tgt, n_src] = n_combinations
 
             if return_conns:

sb_arch_opt/problems/rocket.py

@@ -51,6 +51,8 @@ class RocketArch(HierarchyProblemBase):
     _head_shapes = [HeadShape.CONE, HeadShape.ELLIPTICAL, HeadShape.SPHERE]
 
     _less_constrained = False
+    normalized_constraints = True
+    fail_if_cannot_reach_orbit = False
 
     def __init__(self):
         check_dependency()
@@ -85,7 +87,13 @@ def _arch_evaluate(self, x: np.ndarray, is_active_out: np.ndarray, f_out: np.nda
         lc = self._less_constrained
         rockets = self._get_rockets(x)
         for i, rocket in enumerate(rockets):
-            perf = RocketEvaluator.evaluate(rocket)
+
+            perf = RocketEvaluator.evaluate(
+                rocket,
+                normalize_constraints=self.normalized_constraints,
+                fail_if_cannot_reach_orbit=self.fail_if_cannot_reach_orbit,
+            )
+
             f_out[i, :] = (np.log10(perf.cost), -np.log10(max(1., perf.payload_mass)))
             if lc:
                 g_out[i, :] = (perf.delta_structural, perf.delta_payload, perf.delta_delta_v)
@@ -166,7 +174,7 @@ def _gen_all_discrete_x(self) -> Optional[Tuple[np.ndarray, np.ndarray]]:
             x_stage[:, 1:1+x_engines.shape[1]] = x_engines
             x_stages.append(x_stage)
 
-        x_stages = np.row_stack(x_stages)
+        x_stages = np.vstack(x_stages)
         x_all = np.repeat(x_stages, 3, axis=0)
         x_all[:, [11]] = np.tile(np.array([np.arange(3)]).T, (x_stages.shape[0], 1))  # Head shape
 
@@ -221,15 +229,15 @@ def __repr__(self):
     from pymoo.core.population import Population
     from sb_arch_opt.sampling import HierarchicalSampling
 
-    # problem = RocketArch()
+    problem = RocketArch()
     # x_pf = problem.pareto_set()
     # f_pf = problem.pareto_front()
     # problem = LCRocketArch()
     # problem = SOLCRocketArch(obj=RocketObj.OBJ_COST)
-    problem = SOLCRocketArch(obj=RocketObj.OBJ_PAYLOAD)
+    # problem = SOLCRocketArch(obj=RocketObj.OBJ_PAYLOAD)
     # problem = SOLCRocketArch(obj=RocketObj.OBJ_WEIGHTED)
 
-    problem.plot_pf()
+    # problem.plot_pf()
     # f_pf = problem.pareto_front()
     # f_so = f_pf[:, 0] + f_pf[:, 1]
 

sb_arch_opt/problems/rocket_eval.py

@@ -93,7 +93,7 @@ class RocketEvaluator:
     _rho_interp = None
 
     @classmethod
-    def evaluate(cls, rocket: Rocket) -> RocketPerformance:
+    def evaluate(cls, rocket: Rocket, normalize_constraints=False, fail_if_cannot_reach_orbit=False) -> RocketPerformance:
         check_dependency()
 
         # Calculate geometrical data
@@ -164,7 +164,7 @@ def evaluate(cls, rocket: Rocket) -> RocketPerformance:
         length_ratio = rocket.ellipse_l_ratio if rocket.head_shape == HeadShape.ELLIPTICAL else 0
         payload_mass, h_vector, v_vector, delta_delta_v = cls.calculate_trajectory(
             cone_angle, length_ratio, diameter, stage_thrusts, stage_structural_masses, stage_prop_masses, stage_mdots,
-            rocket.orbit_altitude, m_payload_fix=rocket.payload_mass)
+            rocket.orbit_altitude, m_payload_fix=rocket.payload_mass, normalize_constraint=normalize_constraints)
 
         # Cost estimation
         cost = cls.calculate_cost(stage_n_engines, stage_engine_masses, stage_solid_prop_masses, stage_h2_masses,
@@ -174,6 +174,17 @@ def evaluate(cls, rocket: Rocket) -> RocketPerformance:
         delta_structural = cls.calculate_max_q_constraint(rocket.max_q, h_vector, v_vector)
         delta_payload = cls.calculate_payload_constraint(volume_available, payload_mass, rocket.payload_density)
 
+        if normalize_constraints:
+            delta_structural /= rocket.max_q
+            delta_payload /= volume_available
+
+        if fail_if_cannot_reach_orbit and delta_delta_v < 0:
+            cost = math.nan
+            payload_mass = math.nan
+            delta_structural = math.nan
+            delta_payload = math.nan
+            delta_delta_v = math.nan
+
         return RocketPerformance(
             cost=cost, payload_mass=payload_mass,
             delta_structural=delta_structural, delta_payload=delta_payload, delta_delta_v=-delta_delta_v,
@@ -390,12 +401,13 @@ def modified_atmosphere(cls, x):
 
     @classmethod
     def calculate_trajectory(cls, cone_angle, length_ratio, diameter, T_stages, m_structural_stages, mp_stages,
-                             mdot_stages, h_orbit_target, m_payload_fix=None):
+                             mdot_stages, h_orbit_target, m_payload_fix=None, normalize_constraint=False):
         """Calculation of the launcher trajectory."""
         check_dependency()
 
         mu = 3.986004418e14
         r_earth = 6378e3
+        v_orbit = (mu / (r_earth + h_orbit_target)) ** 0.5
 
         # Drag coefficient calculation depending on head shape
         if cone_angle > 0:
@@ -483,6 +495,10 @@ def stage_state_eq(t_, y):    # Trajectory equation
                 h_first += h.tolist()
                 v_first += v.tolist()
 
+            if len(pos) == 0:
+                # No second stage available
+                return v_first[-1], h_first, v_first
+
             # Second stage
             h_0 = h[pos[0]]
             v_0 = v[pos[0]]
@@ -550,6 +566,10 @@ def stage_state_eq(t_, y):    # Trajectory equation
                 h_second += h.tolist()
                 v_second += v.tolist()
 
+            if len(pos2) == 0:
+                # No third stage available
+                return v_second[-1], h_first+h_second, v_first+v_second
+
             # Third stage
             v_0 = v[pos2[0]]
             t_0 = t[pos2[0]]
@@ -577,28 +597,35 @@ def stage_state_eq(t_, y):    # Trajectory equation
             return v_orbit_final, h_vector_, v_vector_
 
         def try_payload(m_payload):
-            v_orbit = (mu / (r_earth + h_orbit_target)) ** 0.5
             try:
                 v_final_, h_vector_, v_vector_ = simulate_trajectory(m_payload)
 
                 # Orbit minimum speed
                 v_target_diff = v_final_ - v_orbit
 
+                if normalize_constraint:
+                    v_target_diff /= v_orbit
+
                 return v_target_diff, v_final_, h_vector_, v_vector_
             except (IndexError, ValueError):
-                return -v_orbit, 0, [], []
+
+                v_diff = -v_orbit
+                if normalize_constraint:
+                    v_diff /= v_orbit
+
+                return v_diff, 0, [], []
 
         # Evaluate for fixed payload mass
         if m_payload_fix:
             v_tgt_diff, _, h_vector, v_vector = try_payload(m_payload_fix)
             if v_tgt_diff < 0:
-                return 0, [], []
+                return 0, [], [], v_tgt_diff
             return m_payload_fix, h_vector, v_vector, v_tgt_diff
 
         # Check if rocket could be feasible even without payload
         v_tgt_diff, _, _, _ = try_payload(0)
-        if v_tgt_diff <= 0:
-            return 0, [], [], v_tgt_diff
+        # if v_tgt_diff < 0:
+        #     return 0, [], [], v_tgt_diff
 
         m_payload, res = opt.newton(
             lambda mp_: try_payload(mp_)[0], 100, tol=1., maxiter=50, full_output=True, disp=False)

sb_arch_opt/sampling.py

@@ -339,8 +339,8 @@ def _choice(n_choose, n_from, replace=True):
             is_active.append(is_act_all[x_all_choose, :])
             i_x_sampled[x_all_choose] = True
 
-        x = np.row_stack(x)
-        is_active = np.row_stack(is_active)
+        x = np.vstack(x)
+        is_active = np.vstack(is_active)
 
         # Uniformly add discrete vectors if there are not enough (can happen if some groups are very small and there
         # are no continuous dimensions)
@@ -354,8 +354,8 @@ def _choice(n_choose, n_from, replace=True):
             else:
                 i_from_group = np.arange(x_available.shape[0])
 
-            x = np.row_stack([x, x_available[i_from_group, :]])
-            is_active = np.row_stack([is_active, is_act_available[i_from_group, :]])
+            x = np.vstack([x, x_available[i_from_group, :]])
+            is_active = np.vstack([is_active, is_act_available[i_from_group, :]])
 
         return x, is_active
 
@@ -489,7 +489,7 @@ def _choice(n_choose, n_from, replace=True):
                         i_opt_sampled = _choice(n_add, len(opt_values[i_dv]))
                         x_add[:, i_dv] = opt_values[i_dv][i_opt_sampled]
 
-            x = x_add if x is None else np.row_stack([x, x_add])
+            x = x_add if x is None else np.vstack([x, x_add])
 
             # Correct and remove duplicates
             x, is_active = self._correct(problem, repair, x)
@@ -509,8 +509,8 @@ def _choice(n_choose, n_from, replace=True):
         if x.shape[0] < n_samples and has_x_cont:
             n_add = n_samples-x.shape[0]
             i_select_dup = _choice(n_add, x.shape[0])
-            x = np.row_stack(x, x[i_select_dup, :])
-            is_active = np.row_stack(is_active, is_active[i_select_dup, :])
+            x = np.vstack(x, x[i_select_dup, :])
+            is_active = np.vstack(is_active, is_active[i_select_dup, :])
 
         return x, is_active