diff --git a/doc/common/img/limitation_beta_max.pdf b/doc/common/img/limitation_beta_max.pdf
new file mode 100644
index 00000000..49e43d3e
Binary files /dev/null and b/doc/common/img/limitation_beta_max.pdf differ
diff --git a/doc/common/img/limitation_beta_max.png b/doc/common/img/limitation_beta_max.png
new file mode 100644
index 00000000..4446b42d
Binary files /dev/null and b/doc/common/img/limitation_beta_max.png differ
diff --git a/doc/common/img/limitation_beta_min.pdf b/doc/common/img/limitation_beta_min.pdf
new file mode 100644
index 00000000..f130258d
Binary files /dev/null and b/doc/common/img/limitation_beta_min.pdf differ
diff --git a/doc/common/img/limitation_beta_min.png b/doc/common/img/limitation_beta_min.png
new file mode 100644
index 00000000..7360c798
Binary files /dev/null and b/doc/common/img/limitation_beta_min.png differ
diff --git a/doc/ja/source/input.rst b/doc/ja/source/input.rst
index d8606232..ea2bf59d 100644
--- a/doc/ja/source/input.rst
+++ b/doc/ja/source/input.rst
@@ -118,10 +118,11 @@ py2dmat は入力ファイルの形式に `TOML `_ を採
[``runner``] セクション
************************
``Algorithm`` と ``Solver`` を橋渡しする要素である ``Runner`` の設定を記述します。
-サブセクションとして ``mapping`` と ``log`` を持ちます。
+サブセクションとして ``mapping`` 、 ``limitation`` 、 ``log`` を持ちます。
+
[``mapping``] セクション
-************************
+************************************************
``Algorithm`` で探索している :math:`N` 次元のパラメータ :math:`x` から ``Solver`` で使う :math:`M` 次元のパラメータ :math:`y` への写像を定義します。
:math:`N \ne M` となる場合には、 ``solver`` セクションにも ``dimension`` パラメータを指定してください。
@@ -166,8 +167,96 @@ py2dmat は入力ファイルの形式に `TOML `_ を採
を表します。
+[``limitation``] セクション
+************************************************
+
+``Algorithm`` で探索している :math:`N` 次元のパラメータ :math:`x` に、制約条件を課すことが出来ます。
+``Algorithm`` ごとに定義する探索範囲（例：``exchange`` の ``min_list`` や ``max_list`` ） に加えて課すことが出来ます。
+現在は :math:`M` 行 :math:`N` 列の行列Aと :math:`M` 次元の縦ベクトルbから定義される :math:`Ax+b>0` の制約式が利用可能です。具体的に
+
+.. math::
+
+ A_{1,1} x_{1} + A_{1,2} x_{2} + &... + A_{1,N} x_{N} + b_{1} > 0\\
+ A_{2,1} x_{1} + A_{2,2} x_{2} + &... + A_{2,N} x_{N} + b_{2} > 0\\
+ &...\\
+ A_{M,1} x_{1} + A_{M,2} x_{2} + &... + A_{M,N} x_{N} + b_{M} > 0
+
+という制約をかけることが出来ます。
+ここで :math:`M` は制約式の個数（任意）となります。
+
+- ``co_a``
+
+ 形式: リストのリスト、あるいは文字列 (default: [])
+
+ 説明: 制約式の行列 :math:`A` を設定します。
+
+ 行数は制約式数 :math:`M` 列数は探索変数の数 :math:`N` である必要があります。
+
+ ``co_b`` を同時に定義する必要があります。
+
+- ``co_b``
+
+ 形式: リストのリスト、あるいは文字列 (default: [])
+
+ 説明: 制約式の縦ベクトル :math:`b` を設定します。
+
+ 次元数が制約式数 :math:`M` の縦ベクトルを設定する必要があります。
+
+ ``co_a`` を同時に定義する必要があります。
+
+行列の指定方法について、[``mapping``] セクションと同様で、例えば、 ::
+
+ A = [[1,1], [0,1]]
+
+と ::
+
+ A = """
+ 1 1
+ 0 1
+ """
+
+はともに
+
+.. math::
+
+ A = \left(
+ \begin{matrix}
+ 1 & 1 \\
+ 0 & 1
+ \end{matrix}
+ \right)
+
+を表します。また、 ::
+
+ co_b = [[0], [-1]]
+
+と ::
+
+ co_b = """0 -1"""
+
+と ::
+
+ co_b = """
+ 0
+ -1
+ """
+
+はともに
+
+.. math::
+
+ b = \left(
+ \begin{matrix}
+ 0 \\
+ -1
+ \end{matrix}
+ \right)
+
+を表します。
+``co_a`` と ``co_b`` のどちらも定義しない場合、制約式を課さずに探索します。
+
[``log``] セクション
-************************
+************************************************
solver 呼び出しのlogging に関する設定です。
diff --git a/doc/ja/source/tutorial/index.rst b/doc/ja/source/tutorial/index.rst
index 133470a5..8b3ebded 100644
--- a/doc/ja/source/tutorial/index.rst
+++ b/doc/ja/source/tutorial/index.rst
@@ -38,6 +38,7 @@ TRHEPDでは原子座標を与えた場合に、回折データがシミュレ
の5つのアルゴリズムが用意されています。
本チュートリアルでは、最初に順問題プログラム ``sim_trhepd_rheed`` の実行方法、
その後に ``minsearch`` , ``mapper_mpi``, ``bayes``, ``exchange``, ``pamc`` の実行方法について順に説明します。
+また、制約式を用いて探索範囲を制限出来る ``[runner.limitation]`` セクションを使用した実行方法も説明しています。
最後に、自分で順問題ソルバーを定義する簡単な例について説明します。
.. toctree::
@@ -48,6 +49,7 @@ TRHEPDでは原子座標を与えた場合に、回折データがシミュレ
mpi
bayes
exchange
+ limitation
pamc
solver_simple
diff --git a/doc/ja/source/tutorial/limitation.rst b/doc/ja/source/tutorial/limitation.rst
new file mode 100644
index 00000000..dd53f3be
--- /dev/null
+++ b/doc/ja/source/tutorial/limitation.rst
@@ -0,0 +1,192 @@
+制約式を適用したレプリカ交換モンテカルロ法による探索
+==========================================================================
+
+ここでは、 ``[runner.limitation]`` セクションに設定できる制約式機能のチュートリアルを示します。
+例として、レプリカ交換モンテカルロ法を、Himmelblauを対象に探索する計算に制約式を適用します。
+
+サンプルファイルの場所
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+サンプルファイルは ``sample/analytical/limitation`` にあります。
+フォルダには以下のファイルが格納されています。
+
+- ``ref.txt``
+
+ 計算が正しく実行されたか確認するためのファイル(本チュートリアルを行うことで得られる ``best_result.txt`` の回答)。
+
+- ``input.toml``
+
+ メインプログラムの入力ファイル。
+
+- ``do.sh``
+
+ 本チュートリアルを一括計算するために準備されたスクリプト
+
+以下、これらのファイルについて説明したあと、実際の計算結果を紹介します。
+
+入力ファイルの説明
+~~~~~~~~~~~~~~~~~~~
+
+メインプログラム用の入力ファイル ``input.toml`` について説明します。
+``input.toml`` の詳細については入力ファイルに記載されています。
+以下は、サンプルファイルにある ``input.toml`` の中身になります。
+
+.. code-block::
+
+ [base]
+ dimension = 2
+ output_dir = "output"
+
+ [algorithm]
+ name = "exchange"
+ seed = 12345
+
+ [algorithm.param]
+ max_list = [6.0, 6.0]
+ min_list = [-6.0, -6.0]
+ unit_list = [0.3, 0.3]
+
+ [algorithm.exchange]
+ Tmin = 1.0
+ Tmax = 100000.0
+ numsteps = 10000
+ numsteps_exchange = 100
+
+ [solver]
+ name = "analytical"
+ function_name = "himmelblau"
+
+ [runner]
+ [runner.limitation]
+ co_a = [[1, -1],[1, 1]]
+ co_b = [[0], [-1]]
+
+ここではこの入力ファイルを簡単に説明します。
+詳細は入力ファイルのレファレンスを参照してください。
+
+``[base]`` セクションはメインプログラム全体のパラメータです。
+``dimension`` は最適化したい変数の個数で、今の場合は2つの変数の最適化を行うので、``2`` を指定します。
+
+``[algorithm]`` セクションは用いる探索アルゴリズムを設定します。
+交換モンテカルロ法を用いる場合には、 ``name`` に ``"exchange"`` を指定します。
+``label_list`` は、``value_0x`` (x=1,2) を出力する際につけるラベル名のリストです。
+``seed`` は擬似乱数生成器に与える種です。
+
+``[algorithm.param]`` サブセクションは、最適化したいパラメータの範囲などを指定します。
+``min_list`` は最小値、 ``max_list`` は最大値を示します。
+
+``[algorithm.exchange]`` サブセクションは、交換モンテカルロ法のハイパーパラメータを指定します。
+
+- ``numstep`` はモンテカルロ更新の回数です。
+- ``numsteps_exchange`` で指定した回数のモンテカルロ更新の後に、温度交換を試みます。
+- ``Tmin``, ``Tmax`` はそれぞれ温度の下限・上限です。
+- ``Tlogspace`` が ``true`` の場合、温度を対数空間で等分割します。このオプションはデフォルト値が ``true`` であるため、今回の ``input.toml`` に指定は無いですが、 ``true`` になっています。
+
+``[solver]`` セクションではメインプログラムの内部で使用するソルバーを指定します。
+今回は ``analytical`` ソルバーを指定しています。 ``analytical`` ソルバーは ``function_name`` パラメータを用いて関数を設定します。
+今回はHimmelblau 関数を指定しています。
+``analytical`` ソルバーに関してはチュートリアル「順問題ソルバーの追加」を参照してください。
+
+``[runner]`` セクションは ``[runner.limitation]`` サブセクションを持ち、この中に制約式を設定します。
+現在、制約式は :math:`N` 次元のパラメータ :math:`x` 、 :math:`M` 行 :math:`N` 列の行列 :math:`A` 、
+:math:`M` 次元の縦ベクトル :math:`b` から定義される :math:`Ax+b>0` の制約式が利用可能です。
+パラメータとしては、以下の項目が設定可能です。
+
+- ``co_a`` は行列 :math:`A` を設定します。
+- ``co_b`` は縦ベクトル :math:`b` を設定します。
+
+パラメータの詳しい設定方法はマニュアル内「入力ファイル」項の「 [``limitation``] セクション」を参照してください。
+今回は
+
+.. math::
+
+ x_{1} − x_{2} > 0\\
+ x_{1} + x_{2} − 1 > 0
+
+の制約式を課して実行しています。
+
+計算実行
+~~~~~~~~~~~~
+
+最初にサンプルファイルが置いてあるフォルダへ移動します(以下、本ソフトウェアをダウンロードしたディレクトリ直下にいることを仮定します).
+
+.. code-block::
+
+ cd sample/analytical/limitation
+
+そのあとに、メインプログラムを実行します(計算時間は通常のPCで20秒程度で終わります)。
+
+.. code-block::
+
+ mpiexec -np 10 python3 ../../../src/py2dmat_main.py input.toml | tee log.txt
+
+ここではプロセス数10のMPI並列を用いた計算を行っています。
+(Open MPI を用いる場合で、使えるコア数よりも要求プロセス数の方が多い時には、
+``mpiexec`` コマンドに ``--oversubscribed`` オプションを追加してください。)
+実行すると、 ``output`` フォルダが生成され、その中に各ランクのフォルダが作成されます。
+更にその中には、各モンテカルロステップで評価したパラメータおよび目的関数の値を記した ``trial.txt`` ファイルと、
+実際に採択されたパラメータを記した ``result.txt`` ファイルが作成されます。
+ともに書式は同じで、最初の2列がステップ数とプロセス内のwalker 番号、次が温度、3列目が目的関数の値、4列目以降がパラメータです。
+以下は、 ``output/0/result.txt`` ファイルの冒頭部分です。
+
+.. code-block::
+
+ # step walker T fx x1 x2
+ 0 0 1.0 187.94429125133564 5.155393113805774 -2.203493345018569
+ 1 0 1.0 148.23606736778044 4.9995614992887525 -2.370212436322816
+ 2 0 1.0 148.23606736778044 4.9995614992887525 -2.370212436322816
+ 3 0 1.0 148.23606736778044 4.9995614992887525 -2.370212436322816
+
+最後に、 ``output/best_result.txt`` に、目的関数が最小となったパラメータとそれを得たランク、モンテカルロステップの情報が書き込まれます。
+
+.. code-block::
+
+ nprocs = 10
+ rank = 2
+ step = 4523
+ walker = 0
+ fx = 0.00010188398524402734
+ x1 = 3.584944906595298
+ x2 = -1.8506985826548874
+
+なお、一括計算するスクリプトとして ``do.sh`` を用意しています。
+``do.sh`` では ``best_result.txt`` と ``ref.txt`` の差分も比較しています。
+以下、説明は割愛しますが、その中身を掲載します。
+
+.. code-block::
+
+ #!/bin/bash
+ mpiexec -np 10 --oversubscribe python3 ../../../src/py2dmat_main.py input.toml
+
+ echo diff output/best_result.txt ref.txt
+ res=0
+ diff output/best_result.txt ref.txt || res=$?
+ if [ $res -eq 0 ]; then
+ echo TEST PASS
+ true
+ else
+ echo TEST FAILED: best_result.txt and ref.txt differ
+ false
+ fi
+
+計算結果の可視化
+~~~~~~~~~~~~~~~~~~~
+
+``result.txt`` を図示して、制約式を満たした座標のみを探索しているかを確認します。
+今回の場合は、以下のコマンドを打つことで2次元パラメータ空間の図が ``<実行日>_histogram`` フォルダ内に作成されます。
+生成されるヒストグラムは、burn-in期間として最初の1000ステップ分の探索を捨てたデータを使用しています。
+
+.. code-block::
+
+ python3 hist2d_limitation_sample.py -p 10 -i input.toml -b 0.1
+
+作成された図には2本の直線 :math:`x_{1} − x_{2} = 0, x_{1} + x_{2} − 1 = 0` と
+探索結果(事後確率分布のヒストグラム)を図示しています。
+図を見ると :math:`x_{1} − x_{2} > 0, x_{1} + x_{2} − 1 > 0` の範囲のみ探索をしていることが確認できます。
+以下に図の一部を掲載します。
+
+.. figure:: ../../../common/img/limitation_beta_min.*
+
+.. figure:: ../../../common/img/limitation_beta_max.*
+
+ サンプルされたパラメータと確率分布。横軸は ``value_01 (x1)`` , 縦軸は ``value_02 (x2)`` を表す。
diff --git a/sample/analytical/limitation/do.sh b/sample/analytical/limitation/do.sh
new file mode 100644
index 00000000..8de55d57
--- /dev/null
+++ b/sample/analytical/limitation/do.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+mpiexec -np 10 --oversubscribe python3 ../../../src/py2dmat_main.py input.toml
+
+echo diff output/best_result.txt ref.txt
+res=0
+diff output/best_result.txt ref.txt || res=$?
+if [ $res -eq 0 ]; then
+ echo TEST PASS
+ true
+else
+ echo TEST FAILED: best_result.txt and ref.txt differ
+ false
+fi
diff --git a/sample/analytical/limitation/hist2d_limitation_sample.py b/sample/analytical/limitation/hist2d_limitation_sample.py
new file mode 100644
index 00000000..0183c8a1
--- /dev/null
+++ b/sample/analytical/limitation/hist2d_limitation_sample.py
@@ -0,0 +1,210 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as clr
+import toml as tm
+import sys
+import datetime
+import os
+import math as m
+import collections as cl
+import argparse
+import decimal
+
+def read_toml():
+ #input.tomlを取得
+ dict_toml= tm.load(inputfile)
+
+ #変数の範囲の最大値・最小値とTがlogスケールかを取得
+ x1_min,x2_min= dict_toml['algorithm']['param']['min_list']
+ x1_max,x2_max= dict_toml['algorithm']['param']['max_list']
+ T_min= dict_toml['algorithm']['exchange']['Tmin']
+ T_max= dict_toml['algorithm']['exchange']['Tmax']
+ Tlog= True #dict_toml['algorithm']['exchange']['Tlogspace']
+ #Tlogはブーリアン型
+
+ #結果を出力したディレクトリ名を取得
+ output_dicname= dict_toml['base']['output_dir']
+
+ #探索回数を取得
+ numsteps= int(dict_toml['algorithm']['exchange']['numsteps'])
+ numsteps_exchange= int(dict_toml['algorithm']['exchange']['numsteps_exchange'])
+
+ name_for_simulate = ''
+ return x1_min, x2_min, x1_max, x2_max, T_max, T_min, Tlog, output_dicname, numsteps, numsteps_exchange, name_for_simulate
+
+
+def read_result():
+
+ #各ディレクトリから得た値を格納する変数
+ x1= []
+ x2= []
+ fx= []
+ T= []
+ step= range(1,numsteps+1,1)
+
+ #各探索プロセッサにおけるfxの最大値・最小値を格納する変数
+ fx_max_list= []
+ fx_min_list= []
+
+ #以下よりファイルの読み込み
+ i= 0
+ while i <= number_of_replica-1 :
+ print('{0}/{1} ファイル読み込み開始'.format(i+1,number_of_replica))
+ #リストの中にi番目のディレクトリの値を格納するリストを作成
+ x1.append([])
+ x2.append([])
+ fx.append([])
+ T.append([])
+
+ #ファイルを開き、改行区切りのリストにする。
+ f=open('{0}/{1}/result.txt'.format(output_dicname,i))
+ frs= f.read().splitlines()
+ f.close()
+
+ #繰り返し処理で値を抽出、float型へ変換
+ #jが1からスタートなのはresult.txtの1行目には結果の値が入っていないため
+ j= 1
+ while j< len(frs):
+ sp= frs[j].split()
+ x1[i].append(float(sp[-2]))
+ x2[i].append(float(sp[-1]))
+ fx[i].append(float(sp[-3]))
+ T[i].append(float(sp[-4]))
+ j+= 1
+
+ #i番目のresult.txtにおけるfxの最大値・最小値を格納
+ fx_max_list.append(max(fx[i]))
+ fx_min_list.append(min(fx[i]))
+ i+= 1
+
+ print('読み込み完了')
+
+ #読み込みができてるかの確認
+ #print(x1[0][0:5],'\n\n',x1[1][0:5],'\n\n',x1[2][0:5],'\n\n',x1[3][0:5])
+ return x1, x2, fx, T, step, max(fx_max_list), min(fx_min_list)
+
+
+def sort_T():
+ sorted_x1= []
+ sorted_x2= []
+ sorted_fx= []
+ replica_num= []
+ #original_x1= cp.deepcopy(x1)
+ #original_x2= cp.deepcopy(x2)
+ #original_fx= cp.deepcopy(fx)
+
+ Tvalue= make_Tvalue_list()
+
+ i= 0
+ while i < len(Tvalue):
+ sorted_x1.append([])
+ sorted_x2.append([])
+ sorted_fx.append([])
+ replica_num.append([])
+ i+= 1
+
+ i = int(numsteps*burn_in)
+ while i < numsteps:
+ j= 0
+ while j < number_of_replica:
+ k= 0
+ while k < len(Tvalue):
+ if Tvalue[k] == T[j][i]:
+ sorted_x1[k].append(x1[j][i])
+ sorted_x2[k].append(x2[j][i])
+ sorted_fx[k].append(fx[j][i])
+ replica_num[k].append(j)
+ k+= 1
+ j+= 1
+ i+= 1
+
+ sorted_T= True
+ return sorted_x1, sorted_x2, sorted_fx, replica_num, Tvalue, sorted_T
+
+
+def make_Tvalue_list():
+ Tvalue= []
+ i= 0
+ while i < number_of_replica :
+ Tvalue.extend(list(cl.Counter(T[i]).keys()))
+ i+= 1
+ Tvalue= list(cl.Counter(Tvalue).keys())
+ return Tvalue
+
+def soloplot():
+ #繰り返し処理でプロット・図の保存
+ print('グラフ描画開始')
+ min_l = []
+ max_l = []
+
+ l= 0
+ while l< number_of_replica:
+ fig = plt.figure(figsize= (8,8))
+ ax = fig.add_subplot()
+
+ num_of_sample = len(sorted_x1[l][:])
+ weight_l = np.ones(num_of_sample)/num_of_sample
+
+ hst2d = ax.hist2d(
+ sorted_x1[l][:], sorted_x2[l][:],
+ norm= clr.LogNorm(vmin=10**-4, vmax=10**-1),
+ range=[ [x1_min,x1_max] , [x2_min,x2_max] ],
+ cmap= 'Reds', weights=weight_l , bins=100)
+
+ sum_allspace = np.sum(hst2d[0])
+ max_in_one_space = hst2d[0].max()
+ min_in_one_space = hst2d[0][np.nonzero(hst2d[0])].min()
+ min_l.append(min_in_one_space)
+ max_l.append(max_in_one_space)
+
+ cb = fig.colorbar(hst2d[3])
+
+ temp_for_title = decimal.Decimal(1/Tvalue[l]).quantize(
+ decimal.Decimal('0.00001'), rounding=decimal.ROUND_HALF_UP
+ )
+ figtitlestr = f'beta = {temp_for_title}'
+
+ ax.set_title(figtitlestr)
+
+ line_x1 = np.arange(x1_min,x1_max+1,1)
+ line_x2_1 = line_x1
+ line_x2_2 = -line_x1 + 1
+ ax.plot(line_x1,line_x2_1,c="black",alpha=0.3,lw=0.5)
+ ax.plot(line_x1,line_x2_2,c="black",alpha=0.3,lw=0.5)
+ #'''
+ ax.set_xlabel("x1")
+ ax.set_ylabel("x2")
+ ax.set_xlim(x1_min,x1_max)
+ ax.set_ylim(x2_min,x2_max)
+ ax.set_aspect('equal')
+ #軸サイズ設定と図の保存
+ fig.savefig(f'{dirname}/{l}_beta_{temp_for_title}_burn_in_{burn_in}.png',dpi=300)
+ l+= 1
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-p','--mpiprocess',help='MPI process')
+ parser.add_argument('-i','--inputfile',help='input toml file')
+ parser.add_argument('-b','--burn_in',help='burn-in ratio')
+ args = parser.parse_args()
+
+ number_of_replica = int(args.mpiprocess)
+ inputfile = args.inputfile
+ burn_in = float(args.burn_in)
+
+ x1_min, x2_min, x1_max, x2_max, T_max, T_min, Tlog, output_dicname, numsteps, numsteps_exchange, name_for_simulate= read_toml()
+
+ #result.txtから読み取り
+ x1, x2, fx, T, step, fx_max, fx_min = read_result()
+
+ #グラフタイトル・ファイル名用に実行した時間の年日時間分を取得
+ time= datetime.date.today().strftime("%Y%m%d")
+
+ sorted_x1, sorted_x2, sorted_fx, replica_num, Tvalue, sorted_T= sort_T()
+
+ dirname= '{0}_histogram'.format(time)
+ if os.path.exists(dirname) == False:
+ os.makedirs(dirname)
+
+ soloplot()
diff --git a/sample/analytical/limitation/input.toml b/sample/analytical/limitation/input.toml
new file mode 100644
index 00000000..05e4163a
--- /dev/null
+++ b/sample/analytical/limitation/input.toml
@@ -0,0 +1,27 @@
+[base]
+dimension = 2
+output_dir = "output"
+
+[algorithm]
+name = "exchange"
+seed = 12345
+
+[algorithm.param]
+max_list = [6.0, 6.0]
+min_list = [-6.0, -6.0]
+unit_list = [0.3, 0.3]
+
+[algorithm.exchange]
+Tmin = 1.0
+Tmax = 100000.0
+numsteps = 10000
+numsteps_exchange = 100
+
+[solver]
+name = "analytical"
+function_name = "himmelblau"
+
+[runner]
+[runner.limitation]
+co_a = [[1, -1],[1, 1]]
+co_b = [[0], [-1]]
diff --git a/sample/analytical/limitation/ref.txt b/sample/analytical/limitation/ref.txt
new file mode 100644
index 00000000..63b22eae
--- /dev/null
+++ b/sample/analytical/limitation/ref.txt
@@ -0,0 +1,7 @@
+nprocs = 10
+rank = 2
+step = 4523
+walker = 0
+fx = 0.00010188398524402734
+x1 = 3.584944906595298
+x2 = -1.8506985826548874
diff --git a/src/py2dmat/_runner.py b/src/py2dmat/_runner.py
index 7e21e21f..f60818fb 100644
--- a/src/py2dmat/_runner.py
+++ b/src/py2dmat/_runner.py
@@ -184,47 +184,55 @@ def __init__(
self.ndim = info.base["dimension"]
if limitation is None:
info_limitation = info.runner.get("limitation",{})
- co_a: Optional[np.ndarray] = py2dmat.util.read_matrix.read_matrix(
+ co_a: np.ndarray = py2dmat.util.read_matrix.read_matrix(
info_limitation.get("co_a", [])
)
- co_b: Optional[np.ndarray] = py2dmat.util.read_matrix.read_matrix(
+ co_b: np.ndarray = py2dmat.util.read_matrix.read_matrix(
info_limitation.get("co_b", [])
)
- if co_a is not None:
- if co_a.size == 0:
- co_a = None
- else:
- if co_a.ndim != 2:
- raise InputError("co_a should be a matrix")
- if co_a.shape[1] != self.ndim:
- msg ='The number of columns in co_a should be equal to'
- msg+='the value of "dimension" in the [base] section'
- raise InputError(msg)
- n_row_co_a = co_a.shape[0]
- if co_b is not None:
- if co_b.size == 0:
- if co_a is None:
- co_b = None
- else: # co_a is not None:
- msg = "ERROR: co_a is defined but co_b is not."
- raise InputError(msg)
- elif co_b.ndim == 2:
- if co_a is not None:
- if co_b.shape[0] == 1 or co_b.shape[1] == 1:
- co_b = co_b.reshape(-1)
- else:
- raise InputError("co_b should be a vector")
- if co_b.size != n_row_co_a:
- msg ='The number of row in co_a should be equal to'
- msg+='the number of size in co_b'
- raise InputError(msg)
- else: # co_a is None:
- msg = "ERROR: co_b is defined but co_a is not."
+ if co_a.size == 0:
+ is_set_co_a = False
+ else:
+ is_set_co_a = True
+ if co_a.ndim != 2:
+ raise InputError("co_a should be a matrix")
+ if co_a.shape[1] != self.ndim:
+ msg ='The number of columns in co_a should be equal to'
+ msg+='the value of "dimension" in the [base] section'
+ raise InputError(msg)
+ n_row_co_a = co_a.shape[0]
+ if co_b.size == 0:
+ if not is_set_co_a :
+ is_set_co_b = False
+ else: # is_set_co_a is True
+ msg = "ERROR: co_a is defined but co_b is not."
+ raise InputError(msg)
+ elif co_b.ndim == 2:
+ if is_set_co_a:
+ if co_b.shape[0] == 1 or co_b.shape[1] == 1:
+ is_set_co_b = True
+ co_b = co_b.reshape(-1)
+ else:
+ raise InputError("co_b should be a vector")
+ if co_b.size != n_row_co_a:
+ msg ='The number of row in co_a should be equal to'
+ msg+='the number of size in co_b'
raise InputError(msg)
+ else: # not is_set_co_a:
+ msg = "ERROR: co_b is defined but co_a is not."
+ raise InputError(msg)
+ elif co_b.ndim > 2:
+ raise InputError("co_b should be a vector")
+
+ if is_set_co_a and is_set_co_b:
+ is_limitation = True
+ elif (not is_set_co_a) and (not is_set_co_b):
+ is_limitation = False
+ else:
+ msg = "ERROR: Both co_a and co_b must be defined."
+ raise InputError(msg)
- elif co_b.ndim > 2:
- raise InputError("co_b should be a vector")
- self.limitation = py2dmat.util.limitation.Inequality(co_a, co_b)
+ self.limitation = py2dmat.util.limitation.Inequality(co_a, co_b, is_limitation)
def prepare(self, proc_dir: Path):
self.logger.prepare(proc_dir)
diff --git a/src/py2dmat/util/limitation.py b/src/py2dmat/util/limitation.py
index 0c973f7f..5a1f03e4 100644
--- a/src/py2dmat/util/limitation.py
+++ b/src/py2dmat/util/limitation.py
@@ -6,23 +6,14 @@
# for type hints
from pathlib import Path
-from typing import List, Optional, TYPE_CHECKING, Dict, Tuple
+from typing import List, Optional, Dict
class LimitationBase(metaclass=ABCMeta):
@abstractmethod
- def __init__ (self, a: Optional[np.ndarray] = None,
- b: Optional[np.ndarray] = None):
- if (a is not None) or (b is not None):
- if a is None:
- msg = "ERROR: b is defined but a is not."
- raise RuntimeError(msg)
- if b is None:
- msg = "ERROR: a is defined but b is not."
- raise RuntimeError(msg)
- if (a is None) and (b is None):
- self.islimitary = False
- else:
- self.islimitary = True
+ def __init__ (self, a: np.ndarray, b: np.ndarray, is_limitary: bool):
+
+ self.islimitary = is_limitary
+ if self.islimitary:
self.a = a
self.b = b
self.n_formura = a.shape[0]
@@ -33,9 +24,8 @@ def judge(self, x: np.ndarray) -> bool:
pass
class Inequality(LimitationBase):
- def __init__ (self, a: Optional[np.ndarray] = None,
- b: Optional[np.ndarray] = None):
- super().__init__(a, b)
+ def __init__ (self, a: np.ndarray, b: np.ndarray, is_limitary: bool):
+ super().__init__(a, b, is_limitary)
def judge(self, x: np.ndarray) -> bool:
if self.islimitary :