Merge branch 'limitation' into TRHEPD_manybeam_with_limitation

src/py2dmat/algorithm/_algorithm.py

@@ -170,22 +170,28 @@ def _read_param(
         if initial_list.ndim == 1:
             initial_list = initial_list.reshape(1, -1)
         if initial_list.size == 0:
-            # Repeat until an "initial_list" is generated 
-            # that satisfies the constraint formula.
+            initial_list = min_list + (max_list - min_list) * self.rng.rand(
+                num_walkers, self.dimension
+            )
+            # Repeat until an "initial_list" is generated
+            # that satisfies the constraint expression.
+            # If "co_a" and "co_b" are not set in [runner.limitation], 
+            # all(isOK_judge) = true and do not repeat.
             loop_count = 0
+            isOK_judge = np.full(num_walkers, False)
             while True:
-                judge_result = []
-                initial_list = min_list + (max_list - min_list) * self.rng.rand(
-                    num_walkers, self.dimension
-                )
-                for walker_index in range(num_walkers):
-                    judge = self.runner.limitation.judge(
-                                    initial_list[walker_index,:])
-                    judge_result.append(judge)
-                if all(judge_result):
+                for index in np.where(~isOK_judge)[0]:
+                    isOK_judge[index] = self.runner.limitation.judge(
+                                            initial_list[index,:]
+                                            )
+                if np.all(isOK_judge):
                     break
                 else:
-                    loop_count += 1 
+                    initial_list[~isOK_judge] = (
+                            min_list + (max_list - min_list) * self.rng.rand(
+                                np.count_nonzero(~isOK_judge), self.dimension
+                                    )   )
+                    loop_count += 1   
         if initial_list.shape[0] != num_walkers:
             raise exception.InputError(
                 f"ERROR: initial_list.shape[0] != num_walkers ({initial_list.shape[0]} != {num_walkers})"

src/py2dmat/algorithm/min_search.py

@@ -91,6 +91,14 @@ def _f_calc(x_list: np.ndarray, extra_data: bool = False) -> float:
                         )
                     )
                     out_of_range = True
+
+            if not self.runner.limitation.judge(x_list):
+                msg ="Warning: "
+                msg+="Variables do not satisfy the constraint formula.\n"
+                for index in range(dimension):
+                    msg+="{} = {}\n".format(label_list[index],x_list[index])
+                print(msg,end="")
+                out_of_range = True
 
             for index in range(dimension):
                 x_list[index] /= unit_list[index]

src/py2dmat/algorithm/montecarlo.py

@@ -302,7 +302,17 @@ def local_update(
         x_old = copy.copy(self.x)
         if self.iscontinuous:
             self.x = self.propose(x_old)
-            in_range = ((self.xmin <= self.x) & (self.x <= self.xmax)).all(axis=1)
+            #judgement of "in_range"
+            in_range_xmin = self.xmin <= self.x
+            in_range_xmax = self.x <= self.xmax
+            in_range_limitation = np.full(self.nwalkers, False)
+            for index_walker in range(self.nwalkers):
+                in_range_limitation[index_walker] = self.runner.limitation.judge(
+                                                        self.x[index_walker]
+                                                            )
+
+            in_range = (in_range_xmin & in_range_xmax).all(axis=1) \
+                       &in_range_limitation 
         else:
             i_old = copy.copy(self.inode)
             self.inode = self.propose(self.inode)
@@ -313,7 +323,7 @@ def local_update(
         fx_old = self.fx.copy()
         self._evaluate(in_range)
         self._write_result(file_trial, extra_info_to_write=extra_info_to_write)
-
+        
         fdiff = self.fx - fx_old
 
         # Ignore an overflow warning in np.exp(x) for x >~ 710
@@ -322,9 +332,9 @@ def local_update(
         # old_setting = np.seterr(over="ignore")
         old_setting = np.seterr(all="ignore")
         probs = np.exp(-beta * fdiff)
-        # probs[np.isnan(probs)] = 0.0
+        #probs[np.isnan(probs)] = 0.0
         np.seterr(**old_setting)
-
+        
         if not self.iscontinuous:
             probs *= self.ncandidates[i_old] / self.ncandidates[self.inode]
         tocheck = in_range & (probs < 1.0)