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/en/source/input.rst b/doc/en/source/input.rst index 3edb49c6..f28d61e9 100644 --- a/doc/en/source/input.rst +++ b/doc/en/source/input.rst @@ -24,6 +24,10 @@ The input file consists of the following six sections. - Define the mapping from a parameter searched by ``Algorithm`` . +- ``limitation`` + + - Define the limitation (constration) of parameter searched by ``Algorithm`` . + - ``log`` - Specify parameters related to logging of solver calls. @@ -157,6 +161,52 @@ mean \end{matrix} \right). + +[``limitation``] section +************************* + +This section defines the limitation (constraint) in an :math:`N` dimensional parameter searched by ``Algorithm``, :math:`x`, in addition of ``min_list`` and ``max_list``. + +In the current version, a linear inequation with the form :math:`Ax+b>0` is available. + +- ``co_a`` + + Format: List of list of float, or a string (default: ``[]``) + + Description: :math:`N \times M` matrix :math:`A`. An empty list ``[]`` is a shorthand of an identity matrix. + If you want to set it by a string, arrange the elements of the matrix separated with spaces and newlines (see the example). + + +- ``co_b`` + + Format: List of float, or a string (default: ``[]``) + + Description: :math:`M` dimensional vector :math:`b`. An empty list ``[]`` is a shorthand of a zero vector. + If you want to set it by a string, arrange the elements of the vector separated with spaces. + +For example, both :: + + A = [[1,1], [0,1]] + +and :: + + A = """ + 1 1 + 0 1 + """ + +mean + +.. math:: + + A = \left( + \begin{matrix} + 1 & 1 \\ + 0 & 1 + \end{matrix} + \right). + + [``log``] section ************************ diff --git a/doc/en/source/solver/sim-trhepd-rheed.rst b/doc/en/source/solver/sim-trhepd-rheed.rst index 4e2cb2f7..da3fc4d6 100644 --- a/doc/en/source/solver/sim-trhepd-rheed.rst +++ b/doc/en/source/solver/sim-trhepd-rheed.rst @@ -19,10 +19,20 @@ The ``surf.exe`` is called from ``py2dmat``. Input parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Input parameters can be specified in subcsections ``config``, ``post``, ``param``, ``reference`` in ``solver`` section. +Input parameters can be specified in subsections ``config``, ``post``, ``param``, ``reference`` in ``solver`` section. -[``config``] section +[``solver``] section ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +- ``generate_rocking_curve`` + + Format: boolean (default: false) + + Description: Whether to generate ``RockingCurve_calculated.txt``. + If ``true``, ``RockingCurve_calculated.txt`` will be generated in the working directory ``Log%%%_###``. + Note that if ``remove_work_dir`` (in [``post``] subsection) is ``true``, ``Log%%%_###`` will be removed. + +[``solver.config``] subsection +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``surface_exec_file`` @@ -60,30 +70,94 @@ Input parameters can be specified in subcsections ``config``, ``post``, ``param` Description: One of the parameters that specifies the range of output files to be read, calculated by the solver. This parameter specifies the last line to be read. -- ``row_number`` +- ``calculated_info_line`` - Format: integer (default: 8) + Format: integer (default: 2) - Description: One of the parameters that specifies the range of output files to be read, calculated by the solver. This parameter specifies the column to be read. + Description: One of the parameters that specifies the range of output files to be read, calculated by the solver. + This parameter specifies the line to be read, which contains the number of glancing angles (second column) and the number of beams (third column). + +- ``cal_number`` -[``post``] section -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + Format: List of integers -- ``normalization`` + Description: Columns of dataset to be read. Multiple columns can be specified (many-beam condition). - Format: string ("TOTAL" or "MAX", default: "TOTAL") - Description: This parameter specifies whether the experimental and computational data vectors are normalized by the sum of the whole values or by the maximum value. +[``solver.post``] subsection +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This subsection is used to the postprocess -- to specify how to calculate the objective function, that is, the deviation between the experimental and computational data, and to draw the rocking curve. - ``Rfactor_type`` Format: string ("A" or "B", default: "A") - Description: This parameter specifies how to calculate the R-factor. - The experimental and computational data vectors are denoted as :math:`u = (u_{1}, u_{2},...,u_{m})`, - :math:`v = (v_{1}, v_{2},...,v_{m})`, respectively. - When "A" type is chosen, R-factor is defined as :math:`R = (\sum_i^m (u_{i}-v_{i})^{2})^{1/2}`. - When "B" type is chosen, R-factor is defined as :math:`R = (\sum_i^m (u_{i}-v_{i})^{2})^{1/2}/( \sum_i^m u_{i}^2 + \sum_i^m v_{i}^2)`. + Description: This parameter specifies how to calculate the R-factor to be minimized. + Let :math:`n` be the number of dataset, :math:`m` be the number of glancing angles, and :math:`v^{(n)} = (v^{(n)}_{1}, v^{(n)}_{2}, \dots ,v^{(n)}_{m})` be the calculated data. + With the weights of the beams, :math:`w^{(j)}`, R-factors is defined as follows: + + - "A" type: + + .. math:: + + R = \sqrt{ \sum_{j}^{n} w^{(j)} \sum_{i}^{m} \left(u^{(j)}_{i}-v^{(j)}_{i}\right)^{2} } + + - "A2" type: + + .. math:: + + R^{2} = \sum_{j}^{n} w^{(j)} \sum_{i}^{m} \left(u^{(j)}_{i}-v^{(j)}_{i}\right)^{2} + + - "B" type: + + .. math:: + + R = \frac{\sum_{i}^{m} \left(u^{(1)}_{i}-v^{(1)}_{i}\right)^{2}}{\sum_{i}^{m} \left(u^{(1)}_{i}\right)^{2} + \sum_{i}^{m} (v^{(1)}_{i})^2} + + - "B" type is available only for the single dataset (:math:`n=1`). + + +- ``normalization`` + + Format: string ("TOTAL" or "MANY_BEAM") + + Description: This parameter specifies how to normalize the experimental and computational data vectors. + + - "TOTAL" + + - To normalize the data as the summation is 1. + - The number of dataset should be one (the number of ``cal_number`` should be one). + + - "MANY_BEAM" + + - To normalize with weights as specified by ``weight_type``. + + **NOTE: "MAX" is no longer available** + +- ``weight_type`` + + Format: string or ``None``. "calc" or "manual" (default: ``None``) + + Description: The weights of the datasets for the "MANY_BEAM" normalization. + + - "calc" + + .. math:: + + w^{(n)} = \left(\frac{\sum_i^m v^{(n)}_{i}}{\sum_j^n \sum_i^m v^{(j)}_i} \right)^2 + + - "manual" + + :math:`w^{(n)}` is specified by ``spot_weight``. + +- ``spot_weight`` + + Format: list of float (mandatory when ``weight_type`` is "manual") + + Description: The weights of the beams in the calculation of R-factor. + The weights are automatically normalized as the sum be 1. + For example, [3,2,1] means :math:`w^{(1)}=1/2, w^{(2)}=1/3, w^{(3)}=1/6`. - ``omega`` @@ -95,9 +169,9 @@ Input parameters can be specified in subcsections ``config``, ``post``, ``param` Format: boolean (default: false) - Description: Whether to remove working directories after reading R-factor or not + Description: Whether to remove working directories ``Log%%%_###`` after reading R-factor or not -[``param``] section +[``solver.param``] subsection ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``string_list`` @@ -106,14 +180,9 @@ Input parameters can be specified in subcsections ``config``, ``post``, ``param` Description: List of placeholders to be used in the reference template file to create the input file for the solver. These strings will be replaced with the values of the parameters being searched for. -- ``degree_max`` - - Format: float (default: 6.0) - - Description: Maximum angle (in degrees) -[``reference``] section -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +[``solver.reference``] subsection +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``path`` @@ -121,18 +190,24 @@ Input parameters can be specified in subcsections ``config``, ``post``, ``param` Description: Path to the experimental data file. -- ``first`` +- ``reference_first_line`` Format: integer (default: 1) Description: One of the parameters that specify the range of experimental data files to be read. This parameter specifies the first line of the experimental file to be read. -- ``last`` +- ``reference_last_line`` Format: integer (default: 56) Description: One of the parameters that specify the range of experimental data files to be read. This parameter specifies the last line of the experimental file to be read. +- ``exp_number`` + + Format: List of integers + + Description: Columns of dataset to be read. Multiple columns can be specified (many-beam condition). + Reference file ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -169,21 +244,22 @@ In this case, ``value_01``, ``value_02``, and ``value_03`` are the parameters to Target file ^^^^^^^^^^^^^^ This file (``experiment.txt``) contains the data to be targeted. -The first column contains the angle, and the second column contains the calculated value of the reflection intensity multiplied by the weight. +The first column contains the angle, and the second and following columns contain the calculated value of the reflection intensity multiplied by the weight. An example of the file is shown below. .. code-block:: - 0.100000 0.002374995 - 0.200000 0.003614789 - 0.300000 0.005023215 - 0.400000 0.006504978 - 0.500000 0.007990674 - 0.600000 0.009441623 - 0.700000 0.010839445 - 0.800000 0.012174578 - 0.900000 0.013439485 - 1.000000 0.014625579 + 3.00000e-01 8.17149e-03 1.03057e-05 8.88164e-15 ... + 4.00000e-01 1.13871e-02 4.01611e-05 2.23952e-13 ... + 5.00000e-01 1.44044e-02 1.29668e-04 4.53633e-12 ... + 6.00000e-01 1.68659e-02 3.49471e-04 7.38656e-11 ... + 7.00000e-01 1.85375e-02 7.93037e-04 9.67719e-10 ... + 8.00000e-01 1.93113e-02 1.52987e-03 1.02117e-08 ... + 9.00000e-01 1.92590e-02 2.53448e-03 8.69136e-08 ... + 1.00000e+00 1.86780e-02 3.64176e-03 5.97661e-07 ... + 1.10000e+00 1.80255e-02 4.57932e-03 3.32760e-06 ... + 1.20000e+00 1.77339e-02 5.07634e-03 1.50410e-05 ... + 1.30000e+00 1.80264e-02 4.99008e-03 5.53791e-05 ... ... @@ -210,20 +286,49 @@ An example is shown below. output-filename : surf-bulkP.s -``RockingCurve.txt`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +``RockingCurve_calculated.txt`` +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This file is located in the ``Log%%%%%_#####`` folder. -The first line is the header, and the second and subsequent lines are the angle, convoluted calculated/experimental values, normalized calculated/experimental values, and raw calculated values in that order. -An example is shown below. +At the beginning of the file, the lines beginning with ``#`` are headers. +The header contains the values of the input variables, the objective function value ``f(x)``, the parameters ``Rfactor_type``, ``normalization``, ``weight_type``, ``cal_number``, ``spot_weight``, and what is marked in the data portion columns (e.g. ``# #0 glancing_angle``). + +The header is followed by the data. +The first column shows the glancing angle, and the second and subsequent columns show the intensity of each data column. +You can see which data columns are marked in the header. For example, .. code-block:: - #degree convolution_I_calculated I_experiment convolution_I_calculated(normalized) I_experiment(normalized) I_calculated - 0.1 0.0023816127859192407 0.002374995 0.004354402952499057 0.005364578226620574 0.001722 - 0.2 0.003626530149456865 0.003614789 0.006630537795012198 0.008164993342397588 0.003397 - 0.3 0.00504226607469267 0.005023215 0.009218987407498791 0.011346310125551366 0.005026 - 0.4 0.006533558304296079 0.006504978 0.011945579793136154 0.01469327865677437 0.006607 - 0.5 0.00803056955158873 0.007990674 0.014682628499657693 0.018049130948243314 0.008139 - 0.6 0.009493271317558538 0.009441623 0.017356947736613827 0.021326497600946535 0.00962 - 0.7 0.010899633015118851 0.010839445 0.019928258053867838 0.024483862338931763 0.01105 - ... + # #0 glancing_angle + # #1 cal_number=1 + # #2 cal_number=2 + # #3 cal_number=4 + +shows that the first column is the glancing angle, and the second, third, and fourth columns are the calculated data corresponding to the first, second, and fourth columns of the calculated data file, respectively. + +Intencities in each column are normalized so that the sum of the intensity is 1. +To calculate the objective function value (R-factor and R-factor squared), the data columns are weighted by ``spot_weight`` and normalized by ``normalization``. + +.. code-block:: + + #value_01 = 0.00000 value_02 = 0.00000 + #Rfactor_type = A + #normalization = MANY_BEAM + #weight_type = manual + #fx(x) = 0.03686180462340505 + #cal_number = [1, 2, 4, 6, 8] + #spot_weight = [0.933 0.026 0.036 0.003 0.002] + #NOTICE : Intensities are NOT multiplied by spot_weight. + #The intensity I_(spot) for each spot is normalized as in the following equation. + #sum( I_(spot) ) = 1 + # + # #0 glancing_angle + # #1 cal_number=1 + # #2 cal_number=2 + # #3 cal_number=4 + # #4 cal_number=6 + # #5 cal_number=8 + 0.30000 1.278160358686800e-02 1.378767858296659e-04 8.396046839668212e-14 1.342648818357391e-30 6.697979700048016e-53 + 0.40000 1.778953628930054e-02 5.281839702773564e-04 2.108814173486245e-12 2.467220122612354e-28 7.252675318478533e-50 + 0.50000 2.247181148723425e-02 1.671115124520428e-03 4.250758278908295e-11 3.632860054842994e-26 6.291667506376419e-47 + ... + diff --git a/doc/en/source/tutorial/bayes.rst b/doc/en/source/tutorial/bayes.rst index bcf80073..bc5c1370 100644 --- a/doc/en/source/tutorial/bayes.rst +++ b/doc/en/source/tutorial/bayes.rst @@ -7,7 +7,7 @@ This tutorial subscribes how to estimate atomic positions from the experimental Sample files ~~~~~~~~~~~~~~~~~~~~~~~~ -Sample files are available from ``sample/sim-trhepd-rheed/bayes`` . +Sample files are available from ``sample/sim-trhepd-rheed/single_beam/bayes`` . This directory includes the following files: - ``bulk.txt`` @@ -123,14 +123,14 @@ First, move to the folder where the sample file is located (hereinafter, it is a .. code-block:: - cd sample/sim-trhepd-rheed/bayes + cd sample/sim-trhepd-rheed/single_beam/bayes Copy ``bulk.exe`` and ``surf.exe`` as the tutorial for the direct problem. .. code-block:: - cp ../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . - cp ../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . Execute ``bulk.exe`` to generate ``bulkP.b`` . @@ -142,7 +142,7 @@ Then, run the main program (it takes a few secondes) .. code-block:: - python3 ../../../src/py2dmat_main.py input.toml | tee log.txt + python3 ../../../../src/py2dmat_main.py input.toml | tee log.txt This makes a directory with the name of ``0`` . The following standard output will be shown: @@ -216,7 +216,7 @@ I will omit the explanation below, but I will post the contents. ./bulk.exe - time python3 ../../../src/py2dmat_main.py input.toml + time python3 ../../../../src/py2dmat_main.py input.toml echo diff BayesData.txt ref_BayesData.txt res=0 diff --git a/doc/en/source/tutorial/exchange.rst b/doc/en/source/tutorial/exchange.rst index 6d1b13af..e6dce52c 100644 --- a/doc/en/source/tutorial/exchange.rst +++ b/doc/en/source/tutorial/exchange.rst @@ -6,7 +6,7 @@ This tutorial subscribes how to estimate atomic positions from the experimental Sample files ~~~~~~~~~~~~~~~~~~ -Sample files are available from ``sample/sim-trhepd-rheed/exchange`` . +Sample files are available from ``sample/sim-trhepd-rheed/single_beam/exchange`` . This directory includes the following files: - ``bulk.txt`` @@ -126,14 +126,14 @@ First, move to the folder where the sample file is located (hereinafter, it is a .. code-block:: - cd sample/sim-trhepd-rheed/exchange + cd sample/sim-trhepd-rheed/single_beam/exchange Copy ``bulk.exe`` and ``surf.exe`` as the tutorial for the direct problem. .. code-block:: - cp ../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . - cp ../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . Execute ``bulk.exe`` to generate ``bulkP.b`` . @@ -145,7 +145,7 @@ Then, run the main program (it takes a few secondes) .. code-block:: - mpiexec -np 4 python3 ../../../src/py2dmat_main.py input.toml | tee log.txt + mpiexec -np 4 python3 ../../../../src/py2dmat_main.py input.toml | tee log.txt Here, the calculation is performed using MPI parallel with 4 processes. @@ -187,7 +187,7 @@ I will omit the explanation below, but I will post the contents. ./bulk.exe - time mpiexec --oversubscribe -np 4 python3 ../../../src/py2dmat_main.py input.toml + time mpiexec --oversubscribe -np 4 python3 ../../../../src/py2dmat_main.py input.toml echo diff best_result.txt ref.txt res=0 @@ -207,7 +207,7 @@ The ``result.txt`` in each rank folder records the data sampled by each replica, .. code-block:: - python3 ../../../script/separateT.py + python3 ../../../../script/separateT.py The data reorganized for each temperature point is written to ``result_T%.txt`` (``%`` is the index of the temperature point). diff --git a/doc/en/source/tutorial/index.rst b/doc/en/source/tutorial/index.rst index 6b33ffa1..b22f28c5 100644 --- a/doc/en/source/tutorial/index.rst +++ b/doc/en/source/tutorial/index.rst @@ -44,5 +44,6 @@ In this tutorial, we will first introduce how to run the sequential problem prog mpi bayes exchange + limitation pamc solver_simple diff --git a/doc/en/source/tutorial/limitation.rst b/doc/en/source/tutorial/limitation.rst new file mode 100644 index 00000000..baf41136 --- /dev/null +++ b/doc/en/source/tutorial/limitation.rst @@ -0,0 +1,175 @@ +Replica Exchange Monte Carlo search with limitation +========================================================================== + +This tutorial describes the constraint expression function that can be set in the ``[runner.limitation]`` section. +As an example, the replica exchange Monte Carlo method is applied to the calculation of Himmelblau with constraints. + +Sample files location +~~~~~~~~~~~~~~~~~~~~~~~~ + +Sample files are available in the ``sample/analytical/limitation`` directory. +This directory contains the following files. + +- ``ref.txt`` + + File to check if the calculation is executed correctly (answer to obtain by performing this tutorial). + +- ``input.toml`` + + Input file for the main program. + +- ``do.sh`` + + Script prepared to calculate this tutorial at once. + +In the following, we will explain these files, and then introduce the actual calculation results. + +Input files +~~~~~~~~~~~~~~~~~~~ + +The following ``input.toml`` is an input file for the main program. + +.. 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]`` section is the parameter of the main program. +``dimension`` is the number of variables to be optimized, and in this case, it is 2. + +``[algorithm]`` section is the section to set the search algorithm. +``name`` is the name of the search algorithm. In this case, specify ``"exchange"`` for the replica exchange Monte Carlo method. +``seed`` is the seed given to the pseudo-random number generator. + +``[algorithm.param]`` sub-section specifies the range of parameters to be optimized. +``min_list`` indicates the minimum value, and ``max_list`` indicates the maximum value. + +``[algorithm.exchange]`` sub-section specifies the hyperparameters of the replica exchange Monte Carlo method. + +- ``numstep`` is the number of Monte Carlo updates. +- ``numsteps_exchange`` specifies the number of times to attempt temperature exchange. +- ``Tmin`` and ``Tmax`` are the lower and upper limits of the temperature, respectively. +- If ``Tlogspace`` is ``true``, the temperature is divided equally in log space. This option is not specified in this ``input.toml`` because the default value is ``true``. + +``[solver]`` section specifies the solver used internally in the main program. +In this case, the ``analytical`` solver is specified. +The ``analytical`` solver sets the function using the ``function_name`` parameter, and in this case, the Himmelblau function is specified. + +``[runner]`` section has a sub-section ``[runner.limitation]``, and in this section, the constraint expression is set. +In the current version, the constraint expression is defined as :math:`Ax+b>0` where :math:`x` is :math:`N` dimensional input parameter, :math:`A` is a :math:`M \times N` matrix, and :math:`b` is a :math:`M` dimensional vector. +:math:`A` and :math:`b` are set by ``co_a`` and ``co_b``, respectively. +For details, see the ``[limitation]`` section in the input file in the manual. + +In this case, the following constraint is imposed: + +.. math:: + + x_{1} − x_{2} > 0\\ + x_{1} + x_{2} − 1 > 0 + + +Calculation +~~~~~~~~~~~~~~~~ + +First, move to the folder where the sample file is located (assuming that you are directly under the directory where you downloaded this software). + +.. code-block:: + + cd sample/analytical/limitation + +Then, execute the main program as follows (the calculation time will end in about 20 seconds on a normal PC). + +.. code-block:: + + mpiexec -np 10 python3 ../../../src/py2dmat_main.py input.toml | tee log.txt + +In this case, a calculation with 10 MPI parallel processes is performed. +(When using OpenMPI, if the number of processes to be used is greater than the number of available cores, add the ``--oversubscribed`` option to the ``mpiexec`` command.) +After executed, the ``output`` folder is generated, and in it, a subfolder for each rank is created. +Each subfolder contains the results of the calculation. +``trial.txt`` file, which contains the parameters and objective function values evaluated at each Monte Carlo step, and ``result.txt`` file, which contains the parameters actually adopted, are created. +Both files have the same format, with the first two columns being the step number and the walker number within the process, the next being the temperature, the third being the value of the objective function, and the fourth and subsequent being the parameters. +The following is the beginning of the ``output/0/result.txt`` file: + +.. 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 + +Finally, the best parameter and the rank and Monte Carlo step at which the objective function is minimized are written to ``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`` is available as a script to calculate all at once. +Additionally, in ``do.sh``, the difference between ``best_result.txt`` and ``ref.txt`` is also compared. + +.. 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 + +Visualization of the calculation result +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +By visualizing the ``result.txt`` file, we can confirm that the search is only for coordinates that satisfy the constraint expression. +``hist2d_limitation_sample.py`` is prepared to visualize the 2D parameter space. +This generates a histogram of the posterior probability distribution in the ``_histogram`` folder. +The histogram is generated using the data obtained by discarding the first 1000 steps of the search as a burn-in period. + +.. code-block:: + + python3 hist2d_limitation_sample.py -p 10 -i input.toml -b 0.1 + +The figure shows the posterior probability distribution and the two lines :math:`x_{1} − x_{2} = 0, x_{1} + x_{2} − 1 = 0`, +and it is confirmed that the search is only for the range where :math:`x_{1} − x_{2} > 0, x_{1} + x_{2} − 1 > 0`. + +.. figure:: ../../../common/img/limitation_beta_min.* + +.. figure:: ../../../common/img/limitation_beta_max.* diff --git a/doc/en/source/tutorial/minsearch.rst b/doc/en/source/tutorial/minsearch.rst index 3a668f21..3d20a633 100644 --- a/doc/en/source/tutorial/minsearch.rst +++ b/doc/en/source/tutorial/minsearch.rst @@ -22,7 +22,7 @@ In the main program, the Nelder-Mead method (using `scipy.optimize.fmin `_ を採 [``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` 列の行列:math:`A` と :math:`M` 次元の縦ベクトル:math:`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/solver/sim-trhepd-rheed.rst b/doc/ja/source/solver/sim-trhepd-rheed.rst index c769c489..1c9b2a7b 100644 --- a/doc/ja/source/solver/sim-trhepd-rheed.rst +++ b/doc/ja/source/solver/sim-trhepd-rheed.rst @@ -20,11 +20,20 @@ 入力パラメータ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -``solver`` セクション中のサブセクション +``solver`` セクションとセクション中のサブセクション ``config``, ``post``, ``param``, ``reference`` を利用します。 -[``config``] セクション +[``solver``] セクション ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +- ``generate_rocking_curve`` + + 形式: 真偽値 (default: false) + + 説明: ``RockingCurve_calculated.txt`` の生成をするかどうか。 ``RockingCurve_calculated.txt`` は作業ディレクトリ ``Log%%%_###`` の中に保存されます。このオプションが"true"であっても、 ``remove_work_dir`` （[``post``] セクション）が"true"であれば ``Log%%%_###`` は削除されてしまうため、注意してください。 + + +[``config``] セクション（サブセクション） +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``surface_exec_file`` @@ -58,35 +67,94 @@ - ``calculated_last_line`` - 形式: 整数型 (default: 60) + 形式: 整数型 (default: 5 + [``reference``] セクション（サブセクション）にて指定した実験データファイルの視写角数 - 1) 説明: ソルバーにより計算された出力ファイルを読み込む範囲を指定するパラメータ。読み込む最後の行を指定。 -- ``row_number`` +- ``calculated_info_line`` - 形式: 整数型 (default: 8) + 形式: 整数型 (default: 2) - 説明: ソルバーにより計算された出力ファイルを読み込む範囲を指定するパラメータ。読み込む列を指定。 + 説明: ソルバーにより計算された出力ファイルの行を指定するパラメータ。出力ファイル内に記されているglancing angle数やbeam数が記してある行を指定。 + +- ``cal_number`` -[``post``] セクション -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + 形式: 整数型のリスト (設定必須) -- ``normalization`` + 説明: ソルバーにより計算された出力ファイルを読み込む範囲を指定するパラメータ。読み込むデータ列を指定。複数指定が可能。設定した値は、後述の ``exp_number`` （[``reference``] セクション）の値とリストの長さを一致させる必要があります。 - 形式: string型。"TOTAL"または"MAX"のいずれかをとります。 (default: "TOTAL") +[``post``] セクション（サブセクション） +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - 説明: 実験・計算のデータベクトルの規格化の方法。全体の和で規格化するか最大値で規格化するかを指定します。 +実験と計算の「ズレ」の大きさを示し、最小化すべき関数である「目的関数」の定義や、ロッキングカーブを描くためのコンボリューション半値幅などを指定するセクションです。 - ``Rfactor_type`` - 形式: string型。"A"または"B"のいずれかをとります。 (default: "A") + 形式: string型。"A"、"A2"または"B"のいずれかをとります。 (default: "A") + + 説明: 目的関数値の計算方法の指定。 + 視写角のインデックスを :math:`i = 1,2,\dots,m`, ビームのインデックスを :math:`j = 1,2,\dots,n` として、 + 実験・計算のデータをそれぞれ :math:`u^{(j)}_i` と :math:`v^{(j)}_i` で表し、 + 各ビームの重みを :math:`w^{(j)}` と書くとき + + - "A"タイプ: + + .. math:: + + R = \sqrt{ \sum_{j}^{n} w^{(j)} \sum_{i}^{m} \left(u^{(j)}_{i}-v^{(j)}_{i}\right)^{2} } + + - "A2"タイプ: + + .. math:: - 説明: Rファクターの計算方法の指定。 - 実験・計算のデータベクトルを、それぞれ、:math:`u = (u_{1}, u_{2},...,u_{m})`, - :math:`v = (v_{1}, v_{2},...,v_{m})` と書いたとき、 - "A"タイプを指定すると、Rファクター値は :math:`R = (\sum_i^m (u_{i}-v_{i})^{2})^{1/2}` で計算され、 - "B"タイプを指定すると、Rファクター値は :math:`R = (\sum_i^m (u_{i}-v_{i})^{2})^{1/2}/( \sum_i^m u_{i}^2 + \sum_i^m v_{i}^2)` で計算されます。 + R^{2} = \sum_{j}^{n} w^{(j)} \sum_{i}^{m} \left(u^{(j)}_{i}-v^{(j)}_{i}\right)^{2} + + - "B"タイプ: + + .. math:: + + R = \frac{\sum_{i}^{m} \left(u^{(1)}_{i}-v^{(1)}_{i}\right)^{2}}{\sum_{i}^{m} \left(u^{(1)}_{i}\right)^{2} + \sum_{i}^{m} (v^{(1)}_{i})^2} + + - "B"タイプはデータ列が1つの実験・計算データを用いた実行( :math:`n=1` )のみ対応しています。 + +- ``normalization`` + + 形式: string型。"TOTAL"または"MANY_BEAM"のいずれかをとります。 (設定必須) + + 説明: 実験・計算のデータベクトルの規格化の方法。 + + - "TOTAL" + + - 全体の和で規格化をします。 + - 計算に使用するデータ列が1つのみ適用できる規格化方法であり、 ``cal_number`` （[``config``] セクション）および ``exp_number`` （[``reference``] セクション）で設定したリストの長さが1である必要が有ります。 + + - "MANY_BEAM" + + - "MANY_BEAM"はデータ列が2つ以上であるときに利用できる規格化方法です。後述の ``weight_type`` によって規格化方法が変わります。 + + **なお、 normalization="MAX" は廃止となりました。** + +- ``weight_type`` + + 形式: string型または ``None`` 。"calc"または"manual"のいずれかを設定する必要があります。 (default: ``None`` 、 ``normalization = "MANY_BEAM"`` としたとき設定必須) + + 説明: 目的関数値を計算するときの、ビームごとの重み :math:`w^{(j)}` の計算方法を指定します。 + "calc"とした場合、データ列ごとの重み :math:`w^{(n)}` は次の式で与えられます。 + + .. math:: + + w^{(j)} = \left(\frac{\sum_{i=1}^m v^{(j)}_{i}}{\sum_{k=1}^n \sum_{i=1}^m v^{(j)}_i} \right)^2 + + "manual"とした場合、オプション ``spot_weight`` を用いることで、ユーザーが重みを指定可能です。 + +- ``spot_weight`` + + 形式: float型のリスト。 (default: [] 、 ``weight_type = "manual"`` としたとき設定必須) + + 説明: 目的関数値を計算するときの、データ列ごとの重みを設定します。総和が1になるように自動的に規格化されます。 + 例えば、[3,2,1]を指定すると、 :math:`w^{(1)}=1/2, w^{(2)}=1/3, w^{(3)}=1/6` となります。 + - ``omega`` 形式: 実数型 (default: 0.5) @@ -97,11 +165,11 @@ 形式: 真偽値 (default: false) - 説明: R-factor を読み取った後に作業ディレクトリ ``Log%%%_###`` を削除するかどうか + 説明: R-factor を読み取った後に作業ディレクトリ ``Log%%%_###`` を削除するかどうか。なお、 ``generate_rocking_curve`` ([``solver``] セクション) が"true"であっても、本オプションが"true"ならば ``Log%%%_###`` を削除します。 -[``param``] セクション -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +[``param``] セクション（サブセクション） +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``string_list`` @@ -109,14 +177,8 @@ 説明: ソルバーの入力ファイルを作成するための参照用テンプレートファイルで利用するプレースホルダーのリスト。これらの文字列が探索中のパラメータの値に置換されます。 -- ``degree_max`` - - 形式: 実数型 (default: 6.0) - - 説明: 最大角度（度単位）の指定 - -[``reference``] セクション -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +[``reference``] セクション（サブセクション） +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - ``path`` @@ -124,18 +186,23 @@ 説明: 実験データファイルへのパス。 -- ``first`` +- ``reference_first_line`` 形式: 整数型 (default: 1) 説明: 実験データファイルを読み込む範囲を指定するパラメータ。実験ファイルを読み込む最初の行を指定。 -- ``last`` +- ``reference_last_line`` - 形式: 整数型 (default: 56) + 形式: 整数型 (default: 実験データファイルの最後の行の行数) 説明: 実験データファイルを読み込む範囲を指定するパラメータ。実験ファイルを読み込む最後の行を指定。 +- ``exp_number`` + + 形式: 整数型のリスト (default: []、設定必須) + + 説明: 実験データファイルを読み込む範囲を指定するパラメータ。読み込むデータ列を指定。複数指定が可能。設定した値は、前述の ``cal_number`` （[``config``] セクション）の値とリストの長さを一致させる必要があります。 ソルバー用補助ファイル ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -151,45 +218,45 @@ .. code-block:: 2 ,NELMS, -------- Ge(001)-c4x2 - 32,1.0,0.1 ,Ge Z,da1,sap + 32,1.2,0.15 ,Ge Z,da1,sap 0.6,0.6,0.6 ,BH(I),BK(I),BZ(I) - 32,1.0,0.1 ,Ge Z,da1,sap + 32,1.2,0.15 ,Ge Z,da1,sap 0.4,0.4,0.4 ,BH(I),BK(I),BZ(I) - 9,4,0,0,2, 2.0,-0.5,0.5 ,NSGS,msa,msb,nsa,nsb,dthick,DXS,DYS + 9,4,0,0,2,1.7,-0.5,0.5 ,NSGS,msa,msb,nsa,nsb,dthick,DXS,DYS 8 ,NATM - 1, 1.0, 1.34502591 1 value_01 ,IELM(I),ocr(I),X(I),Y(I),Z(I) - 1, 1.0, 0.752457792 1 value_02 - 2, 1.0, 1.480003343 1.465005851 value_03 - 2, 1.0, 2 1.497500418 2.281675 - 2, 1.0, 1 1.5 1.991675 - 2, 1.0, 0 1 0.847225 - 2, 1.0, 2 1 0.807225 - 2, 1.0, 1.009998328 1 0.597225 + 1, 1.0, value_01, 1.00000, 5.231000 ,IELM(I),ocr(I),X(I),Y(I),Z(I + 1, 1.0, value_02, 1.00000, 4.371000 + 2, 1.0, 1.50000, 1.50000, 3.596000 + 2, 1.0, 2.00000, 1.49751, 2.100000 + 2, 1.0, 1.00000, 1.50000, 2.000000 + 2, 1.0, 0.00000, 1.00000, 0.849425 + 2, 1.0, 2.00000, 1.00000, 0.809425 + 2, 1.0, 1.00997, 1.00000, 0.599425 1,1 ,(WDOM,I=1,NDOM) - -この場合、 ``value_01``, ``value_02``, ``value_03`` が動かすパラメータとなります。 +この場合、 ``value_01``, ``value_02`` が動かすパラメータとなります。 ターゲット参照ファイル ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ターゲットにするデータが格納されたファイル ``experiment.txt`` を指定します。 -第一列に角度、第二列に反射強度に重みをかけて計算した値が入ってます。 -以下、ファイルの例を指定します。 +第一列に角度、第二列以降に反射強度にコンボリューションを計算した値が入ってます。 +以下、ファイルの例を示します。 .. code-block:: - 0.100000 0.002374995 - 0.200000 0.003614789 - 0.300000 0.005023215 - 0.400000 0.006504978 - 0.500000 0.007990674 - 0.600000 0.009441623 - 0.700000 0.010839445 - 0.800000 0.012174578 - 0.900000 0.013439485 - 1.000000 0.014625579 + 3.00000e-01 8.17149e-03 1.03057e-05 8.88164e-15 ... + 4.00000e-01 1.13871e-02 4.01611e-05 2.23952e-13 ... + 5.00000e-01 1.44044e-02 1.29668e-04 4.53633e-12 ... + 6.00000e-01 1.68659e-02 3.49471e-04 7.38656e-11 ... + 7.00000e-01 1.85375e-02 7.93037e-04 9.67719e-10 ... + 8.00000e-01 1.93113e-02 1.52987e-03 1.02117e-08 ... + 9.00000e-01 1.92590e-02 2.53448e-03 8.69136e-08 ... + 1.00000e+00 1.86780e-02 3.64176e-03 5.97661e-07 ... + 1.10000e+00 1.80255e-02 4.57932e-03 3.32760e-06 ... + 1.20000e+00 1.77339e-02 5.07634e-03 1.50410e-05 ... + 1.30000e+00 1.80264e-02 4.99008e-03 5.53791e-05 ... ... @@ -219,21 +286,50 @@ output-filename : surf-bulkP.s -``RockingCurve.txt`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +``RockingCurve_calculated.txt`` +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +``generate_rocking_curve`` ([``solver``] セクション) が"true"の場合のみ ``Log%%%%%_#####`` フォルダに出力されます。 + +ファイル冒頭、 ``#`` で始まる行はヘッダーです。 +ヘッダーには探索変数の値、目的関数値 ``f(x)`` オプションで指定した ``Rfactor_type``, ``normalization`, ``weight_type``, ``cal_number``, オプションで指定またはプログラムが計算したデータ列ごとの重み ``spot_weight``, データ部分のどの列に何が記されているか(例： ``# #0 glanceing_angle`` など)が記されています。 + +``#`` が付いていない部分はデータ表記部分になります。1列目は視写角、2列目以降はデータ列ごとに強度が記しています。どのデータ列が記されているかはヘッダーの表記で確認できます。例えば + +.. code-block:: + + # #0 glancing_angle + # #1 cal_number=1 + # #2 cal_number=2 + # #3 cal_number=4 + +との記載があれば、1列目は視写角、2列目は計算データファイルの1列目に相当する反射強度、3列目は計算データファイルの2列目に相当する反射強度、4列目は計算データファイルの4列目に相当する反射強度が記されていることがわかります。 + +また、各列の反射強度は各列の総和が1になるように規格化されています。目的関数値（R-factor及びR-factorの二乗）を算出する際は、データ列ごとの重み ``spot_weight`` を加味して計算されています。 -``Log%%%%%_#####`` フォルダに出力されます。 -1行目にヘッダ、2行目以降は角度、コンボリューションされた計算値・実験値、規格化された計算値・実験値と、生の計算値が順に出力されます。 以下、出力例です。 .. code-block:: - #degree convolution_I_calculated I_experiment convolution_I_calculated(normalized) I_experiment(normalized) I_calculated - 0.1 0.0023816127859192407 0.002374995 0.004354402952499057 0.005364578226620574 0.001722 - 0.2 0.003626530149456865 0.003614789 0.006630537795012198 0.008164993342397588 0.003397 - 0.3 0.00504226607469267 0.005023215 0.009218987407498791 0.011346310125551366 0.005026 - 0.4 0.006533558304296079 0.006504978 0.011945579793136154 0.01469327865677437 0.006607 - 0.5 0.00803056955158873 0.007990674 0.014682628499657693 0.018049130948243314 0.008139 - 0.6 0.009493271317558538 0.009441623 0.017356947736613827 0.021326497600946535 0.00962 - 0.7 0.010899633015118851 0.010839445 0.019928258053867838 0.024483862338931763 0.01105 - ... + #value_01 = 0.00000 value_02 = 0.00000 + #Rfactor_type = A + #normalization = MANY_BEAM + #weight_type = manual + #fx(x) = 0.03686180462340505 + #cal_number = [1, 2, 4, 6, 8] + #spot_weight = [0.933 0.026 0.036 0.003 0.002] + #NOTICE : Intensities are NOT multiplied by spot_weight. + #The intensity I_(spot) for each spot is normalized as in the following equation. + #sum( I_(spot) ) = 1 + # + # #0 glancing_angle + # #1 cal_number=1 + # #2 cal_number=2 + # #3 cal_number=4 + # #4 cal_number=6 + # #5 cal_number=8 + 0.30000 1.278160358686800e-02 1.378767858296659e-04 8.396046839668212e-14 1.342648818357391e-30 6.697979700048016e-53 + 0.40000 1.778953628930054e-02 5.281839702773564e-04 2.108814173486245e-12 2.467220122612354e-28 7.252675318478533e-50 + 0.50000 2.247181148723425e-02 1.671115124520428e-03 4.250758278908295e-11 3.632860054842994e-26 6.291667506376419e-47 + ... + diff --git a/doc/ja/source/tutorial/bayes.rst b/doc/ja/source/tutorial/bayes.rst index 001061a5..3df2625d 100644 --- a/doc/ja/source/tutorial/bayes.rst +++ b/doc/ja/source/tutorial/bayes.rst @@ -8,7 +8,7 @@ サンプルファイルの場所 ~~~~~~~~~~~~~~~~~~~~~~~~ -サンプルファイルは ``sample/sim-trhepd-rheed/bayes`` にあります。 +サンプルファイルは ``sample/sim-trhepd-rheed/single_beam/bayes`` にあります。 フォルダには以下のファイルが格納されています。 - ``bulk.txt`` @@ -166,8 +166,8 @@ .. code-block:: - cp ../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . - cp ../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . 最初に ``bulk.exe`` を実行し、 ``bulkP.b`` を作成します。 @@ -179,7 +179,7 @@ .. code-block:: - python3 ../../../src/py2dmat_main.py input.toml | tee log.txt + python3 ../../../../src/py2dmat_main.py input.toml | tee log.txt 実行すると、各ランクのフォルダが作成されます。 以下の様な標準出力がされます。 @@ -255,7 +255,7 @@ ./bulk.exe - time python3 ../../../src/py2dmat_main.py input.toml + time python3 ../../../../src/py2dmat_main.py input.toml echo diff BayesData.txt ref_BayesData.txt res=0 diff --git a/doc/ja/source/tutorial/exchange.rst b/doc/ja/source/tutorial/exchange.rst index 8c50f52d..8222b00c 100644 --- a/doc/ja/source/tutorial/exchange.rst +++ b/doc/ja/source/tutorial/exchange.rst @@ -7,7 +7,7 @@ サンプルファイルの場所 ~~~~~~~~~~~~~~~~~~~~~~~~ -サンプルファイルは ``sample/sim-threpd-rheed/exchange`` にあります。 +サンプルファイルは ``sample/sim-threpd-rheed/single_beam/exchange`` にあります。 フォルダには以下のファイルが格納されています。 - ``bulk.txt`` @@ -116,14 +116,14 @@ .. code-block:: - cd sample/sim-trhepd-rheed/exchange + cd sample/sim-trhepd-rheed/single_beam/exchange 順問題の時と同様に、 ``bulk.exe`` と ``surf.exe`` をコピーします。 .. code-block:: - cp ../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . - cp ../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/bulk.exe . + cp ../../../../../sim-trhepd-rheed/src/TRHEPD/surf.exe . 最初に ``bulk.exe`` を実行し、``bulkP.b`` を作成します。 @@ -135,7 +135,7 @@ .. code-block:: - mpiexec -np 4 python3 ../../../src/py2dmat_main.py input.toml | tee log.txt + mpiexec -np 4 python3 ../../../../src/py2dmat_main.py input.toml | tee log.txt ここではプロセス数4のMPI並列を用いた計算を行っています。 (Open MPI を用いる場合で、使えるコア数よりも要求プロセス数の方が多い時には、 @@ -177,7 +177,7 @@ ./bulk.exe - time mpiexec --oversubscribe -np 4 python3 ../../../src/py2dmat_main.py input.toml + time mpiexec --oversubscribe -np 4 python3 ../../../../src/py2dmat_main.py input.toml echo diff best_result.txt ref.txt res=0 @@ -198,7 +198,7 @@ Py2DMat の実装では同一レプリカが様々な温度のサンプルを保 .. code-block:: - python3 ../../../script/separateT.py + python3 ../../../../script/separateT.py ``result_T%.txt`` に各温度点ごとにまとめなおされたデータが書き込まれます(``%`` は温度点のindex)。 1列目がステップ、2列めがランク、3列目が目的関数の値、4列目以降がパラメータです。 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..8ad1a8f4 --- /dev/null +++ b/doc/ja/source/tutorial/limitation.rst @@ -0,0 +1,191 @@ +制約式を適用したレプリカ交換モンテカルロ法による探索 +========================================================================== + +ここでは、 ``[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"`` を指定します。 +``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/doc/ja/source/tutorial/minsearch.rst b/doc/ja/source/tutorial/minsearch.rst index cdf0fcca..192ab9bb 100644 --- a/doc/ja/source/tutorial/minsearch.rst +++ b/doc/ja/source/tutorial/minsearch.rst @@ -24,7 +24,7 @@ Nelder-Mead法 (`scipy.optimize.fmin 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) + + self.limitation = py2dmat.util.limitation.Inequality(co_a, co_b, is_limitation) def prepare(self, proc_dir: Path): self.logger.prepare(proc_dir) @@ -185,14 +240,17 @@ def prepare(self, proc_dir: Path): def submit( self, message: py2dmat.Message, nprocs: int = 1, nthreads: int = 1 ) -> float: - x = self.mapping(message.x) - message_indeed = py2dmat.Message(x, message.step, message.set) - self.solver.prepare(message_indeed) - cwd = os.getcwd() - os.chdir(self.solver.work_dir) - self.solver.run(nprocs, nthreads) - os.chdir(cwd) - result = self.solver.get_results() + if self.limitation.judge(message.x): + x = self.mapping(message.x) + message_indeed = py2dmat.Message(x, message.step, message.set) + self.solver.prepare(message_indeed) + cwd = os.getcwd() + os.chdir(self.solver.work_dir) + self.solver.run(nprocs, nthreads) + os.chdir(cwd) + result = self.solver.get_results() + else: + result = np.inf self.logger.count(message, result) return result diff --git a/src/py2dmat/algorithm/_algorithm.py b/src/py2dmat/algorithm/_algorithm.py index cf503b7c..ab2b95dd 100644 --- a/src/py2dmat/algorithm/_algorithm.py +++ b/src/py2dmat/algorithm/_algorithm.py @@ -26,6 +26,7 @@ import numpy as np import py2dmat +import py2dmat.util.limitation from py2dmat import exception, mpi # for type hints @@ -172,14 +173,41 @@ def _read_param( 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: + 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: + 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})" ) if initial_list.shape[1] != self.dimension: raise exception.InputError( - f"ERROR: initial_list.shape[1] != dimension ({initial_list.shape[1]} != {self.dimension})" - ) + f"ERROR: initial_list.shape[1] != dimension ({initial_list.shape[1]} != {self.dimension})" ) + judge_result = [] + for walker_index in range(num_walkers): + judge = self.runner.limitation.judge( + initial_list[walker_index,:]) + judge_result.append(judge) + if not all(judge_result): + raise exception.InputError( + "ERROR: initial_list does not satisfy the constraint formula." + ) return initial_list, min_list, max_list, unit_list def _meshgrid( diff --git a/src/py2dmat/algorithm/min_search.py b/src/py2dmat/algorithm/min_search.py index 5b6f2c07..46495472 100644 --- a/src/py2dmat/algorithm/min_search.py +++ b/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] diff --git a/src/py2dmat/algorithm/montecarlo.py b/src/py2dmat/algorithm/montecarlo.py index 1e4cae1c..fd3cb3a2 100644 --- a/src/py2dmat/algorithm/montecarlo.py +++ b/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) diff --git a/src/py2dmat/solver/lib_make_convolution.py b/src/py2dmat/solver/lib_make_convolution.py new file mode 100644 index 00000000..a4cc12c4 --- /dev/null +++ b/src/py2dmat/solver/lib_make_convolution.py @@ -0,0 +1,50 @@ +import numpy as np +import copy + + +def calc(RC_data_org, number_of_beams, number_of_glancing_angles, omega, verbose_mode): + + sigma = 0.5 * omega / (np.sqrt(2.0 * np.log(2.0))) + + def g(x): + g = (1.0 / (sigma * np.sqrt(2.0 * np.pi))) * np.exp( + -0.5 * x**2.0 / sigma**2.0 + ) + + return g + + RC_data_cnv = np.zeros((number_of_glancing_angles, number_of_beams + 1)) + # copy glancing angle + RC_data_cnv[:, 0] = copy.deepcopy(RC_data_org[:, 0]) + + if verbose_mode: + print("RC_data_org =\n", RC_data_org) + print("RC_data_cnv =\n", RC_data_cnv) + + for beam_index in range(number_of_beams): + for index in range(number_of_glancing_angles): + integral = 0.0 + angle_interval = 0.0 + for index2 in range(number_of_glancing_angles): + if index2 < number_of_glancing_angles - 1: + angle_interval = RC_data_org[index2 + 1, 0] - RC_data_org[index2, 0] + integral += ( + RC_data_org[index2, beam_index + 1] + * g(RC_data_org[index, 0] - RC_data_org[index2, 0]) + * angle_interval + ) + if verbose_mode: + print( + "beam_index, index, index2, g(RC_data_org[index,0] - RC_data_org[index2,0]) =", + beam_index, + index, + index2, + g(RC_data_org[index, 0] - RC_data_org[index2, 0]), + ) + RC_data_cnv[index, beam_index + 1] = integral + + if verbose_mode: + print("RC_data_cnv =\n", RC_data_cnv) + + return RC_data_cnv + diff --git a/src/py2dmat/solver/sim_trhepd_rheed.py b/src/py2dmat/solver/sim_trhepd_rheed.py index 35cc92b1..6effbf05 100644 --- a/src/py2dmat/solver/sim_trhepd_rheed.py +++ b/src/py2dmat/solver/sim_trhepd_rheed.py @@ -1,81 +1,182 @@ -# 2DMAT -- Data-analysis software of quantum beam diffraction experiments for 2D material structure -# Copyright (C) 2020- The University of Tokyo -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see http://www.gnu.org/licenses/. - -from typing import List import itertools import os import os.path import shutil -import subprocess -from pathlib import Path +import time +import ctypes import numpy as np import py2dmat -from py2dmat import exception +from py2dmat import exception, mpi +from . import lib_make_convolution + +# for type hints +from pathlib import Path +from typing import List, Dict, Optional, TYPE_CHECKING + +if TYPE_CHECKING: + from mpi4py import MPI class Solver(py2dmat.solver.SolverBase): + mpicomm: Optional["MPI.Comm"] + mpisize: int + mpirank: int + run_scheme: str + isLogmode: bool + detail_timer: Dict path_to_solver: Path def __init__(self, info: py2dmat.Info): super().__init__(info) + self.mpicomm = mpi.comm() + self.mpisize = mpi.size() + self.mpirank = mpi.rank() + + self._name = "sim_trhepd_rheed_mb_connect" + + self.run_scheme = info.solver.get("run_scheme", "subprocess") + scheme_list = ["subprocess", "connect_so"] + + if self.run_scheme not in scheme_list: + raise exception.InputError("ERROR : solver.run_scheme should be 'subprocess' or 'connect_so'.") - self._name = "sim_trhepd_rheed" - p2solver = info.solver["config"].get("surface_exec_file", "surf.exe") - if os.path.dirname(p2solver) != "": - # ignore ENV[PATH] - self.path_to_solver = self.root_dir / Path(p2solver).expanduser() + if self.run_scheme == "connect_so": + self.load_so() + + elif self.run_scheme == "subprocess": + # path to surf.exe + p2solver = info.solver["config"].get("surface_exec_file", "surf.exe") + if os.path.dirname(p2solver) != "": + # ignore ENV[PATH] + self.path_to_solver = self.root_dir / Path(p2solver).expanduser() + else: + for P in itertools.chain( + [self.root_dir], os.environ["PATH"].split(":") + ): + self.path_to_solver = Path(P) / p2solver + if os.access(self.path_to_solver, mode=os.X_OK): + break + if not os.access(self.path_to_solver, mode=os.X_OK): + raise exception.InputError(f"ERROR: solver ({p2solver}) is not found") + + self.isLogmode = False + self.set_detail_timer() + + self.input = Solver.Input(info, self.isLogmode, self.detail_timer) + self.output = Solver.Output(info, self.isLogmode, self.detail_timer) + + def set_detail_timer(self) -> None: + # TODO: Operate log_mode with toml file. Generate txt of detail_timer. + if self.isLogmode: + self.detail_timer = {} + self.detail_timer["prepare_Log-directory"] = 0 + self.detail_timer["make_surf_input"] = 0 + self.detail_timer["launch_STR"] = 0 + self.detail_timer["load_STR_result"] = 0 + self.detail_timer["convolution"] = 0 + self.detail_timer["normalize_calc_I"] = 0 + self.detail_timer["calculate_R-factor"] = 0 + self.detail_timer["make_RockingCurve.txt"] = 0 + self.detail_timer["delete_Log-directory"] = 0 else: - for P in itertools.chain([self.root_dir], os.environ["PATH"].split(":")): - self.path_to_solver = Path(P) / p2solver - if os.access(self.path_to_solver, mode=os.X_OK): - break - if not os.access(self.path_to_solver, mode=os.X_OK): - raise exception.InputError(f"ERROR: solver ({p2solver}) is not found") - info_config = info.solver.get("config", {}) - - self.input = Solver.Input(info) - self.output = Solver.Output(info) - self.result = None + self.detail_timer = {} + + def default_run_scheme(self) -> str: + """ + Return + ------- + str + run_scheme. + """ + return self.run_scheme + + def command(self) -> List[str]: + """Command to invoke solver""" + return [str(self.path_to_solver)] def prepare(self, message: py2dmat.Message) -> None: fitted_x_list, subdir = self.input.prepare(message) self.work_dir = self.proc_dir / Path(subdir) + self.output.prepare(fitted_x_list) - self.result = None + + def get_results(self) -> float: + return self.output.get_results(self.work_dir) + + def load_so(self) -> None: + self.lib = np.ctypeslib.load_library("surf.so", os.path.dirname(__file__)) + self.lib.surf_so.argtypes = ( + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_int), + np.ctypeslib.ndpointer(), + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_int), + np.ctypeslib.ndpointer(), + np.ctypeslib.ndpointer(), + ) + self.lib.surf_so.restype = ctypes.c_void_p + + def launch_so(self) -> None: + + n_template_file = len(self.input.template_file) + m_template_file = self.input.surf_template_width_for_fortran + n_bulk_file = len(self.input.bulk_file) + m_bulk_file = self.input.bulk_out_width_for_fortran + + # NOTE: The "20480" is related to the following directive in surf_so.f90. + # character(c_char), intent(inout) :: surf_out(20480) + emp_str = " " * 20480 + self.output.surf_output = np.array([emp_str.encode()]) + self.lib.surf_so( + ctypes.byref(ctypes.c_int(n_template_file)), + ctypes.byref(ctypes.c_int(m_template_file)), + self.input.template_file, + ctypes.byref(ctypes.c_int(n_bulk_file)), + ctypes.byref(ctypes.c_int(m_bulk_file)), + self.input.bulk_file, + self.output.surf_output, + ) + self.output.surf_output = self.output.surf_output[0].decode().splitlines() def run(self, nprocs: int = 1, nthreads: int = 1) -> None: - self._run_by_subprocess([str(self.path_to_solver)]) + if self.isLogmode: + time_sta = time.perf_counter() - def get_results(self) -> float: - if self.result is None: - self.result = self.output.get_results(self.work_dir) - return self.result + if self.run_scheme == "connect_so": + self.launch_so() + elif self.run_scheme == "subprocess": + self._run_by_subprocess([str(self.path_to_solver)]) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["launch_STR"] += time_end - time_sta class Input(object): + mpicomm: Optional["MPI.Comm"] + mpisize: int + mpirank: int + root_dir: Path output_dir: Path dimension: int + run_scheme: str + generate_rocking_curve: bool string_list: List[str] bulk_output_file: Path surface_input_file: Path surface_template_file: Path + template_file_origin: List[str] + + def __init__(self, info, isLogmode, detail_timer): + self.mpicomm = mpi.comm() + self.mpisize = mpi.size() + self.mpirank = mpi.rank() + + self.isLogmode = isLogmode + self.detail_timer = detail_timer - def __init__(self, info): self.root_dir = info.base["root_dir"] self.output_dir = info.base["output_dir"] @@ -84,7 +185,17 @@ def __init__(self, info): else: self.dimension = info.base["dimension"] + # read info info_s = info.solver + self.run_scheme = info_s["run_scheme"] + self.generate_rocking_curve = info_s.get("generate_rocking_curve", False) + + # NOTE: + # surf_template_width_for_fortran: Number of strings per line of template.txt data for surf.so. + # bulk_out_width_for_fortran: Number of strings per line of bulkP.txt data for surf.so. + if self.run_scheme == "connect_so": + self.surf_template_width_for_fortran = 128 + self.bulk_out_width_for_fortran = 1024 info_param = info_s.get("param", {}) v = info_param.setdefault("string_list", ["value_01", "value_02"]) @@ -107,17 +218,58 @@ def __init__(self, info): f"ERROR: surface_template_file ({self.surface_template_file}) does not exist" ) - self._check_template() + if self.mpirank == 0: + self._check_template() + temp_origin = self.load_surface_template_file(filename) + else: + temp_origin = None + self.template_file_origin = self.mpicomm.bcast(temp_origin, root=0) + + if self.run_scheme == "connect_so": + filename = info_config.get("bulk_output_file", "bulkP.txt") + filename = Path(filename).expanduser().resolve() + self.bulk_output_file = self.root_dir / filename + if not self.bulk_output_file.exists(): + raise exception.InputError( + f"ERROR: bulk_output_file ({self.bulk_output_file}) does not exist" + ) + + if self.mpirank == 0: + bulk_f = self.load_bulk_output_file(filename) + else: + bulk_f = None + self.bulk_file = self.mpicomm.bcast(bulk_f, root=0) - filename = info_config.get("bulk_output_file", "bulkP.b") - filename = Path(filename).expanduser().resolve() - self.bulk_output_file = self.root_dir / filename - if not self.bulk_output_file.exists(): - raise exception.InputError( - f"ERROR: bulk_output_file ({self.bulk_output_file}) does not exist" - ) + else: + filename = info_config.get("bulk_output_file", "bulkP.b") + filename = Path(filename).expanduser().resolve() + self.bulk_output_file = self.root_dir / filename + if not self.bulk_output_file.exists(): + raise exception.InputError( + f"ERROR: bulk_output_file ({self.bulk_output_file}) does not exist" + ) + + def load_surface_template_file(self, filename): + template_file = [] + with open(self.surface_template_file) as f: + for line in f: + template_file.append(line) + return template_file + + def load_bulk_output_file(self, filename): + bulk_file = [] + with open(self.bulk_output_file) as f: + for line in f: + line = line.replace("\t", " ").replace("\n", " ") + line = line.encode().ljust(self.bulk_out_width_for_fortran) + bulk_file.append(line) + bulk_f = np.array(bulk_file) + return bulk_f def prepare(self, message: py2dmat.Message): + if self.isLogmode: + time_sta = time.perf_counter() + x_list = message.x step = message.step iset = message.set @@ -132,25 +284,77 @@ def prepare(self, message: py2dmat.Message): fitted_value += format(x_list[index], ".8f") fitted_value = fitted_value[: len(string_list[index])] fitted_x_list.append(fitted_value) - for index in range(dimension): - print(string_list[index], "=", fitted_x_list[index]) - folder_name = self._pre_bulk(step, bulk_output_file, iset) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["make_surf_input"] += time_end - time_sta + + if self.isLogmode: + time_sta = time.perf_counter() + + folder_name = "." + if self.generate_rocking_curve: + folder_name = self._pre_bulk(step, bulk_output_file, iset) + else: + if self.run_scheme == "connect_so": + folder_name = "." + + elif self.run_scheme == "subprocess": + # make workdir and copy bulk output file + folder_name = self._pre_bulk(step, bulk_output_file, iset) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["prepare_Log-directory"] += time_end - time_sta + + if self.isLogmode: + time_sta = time.perf_counter() + self._replace(fitted_x_list, folder_name) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["make_surf_input"] += time_end - time_sta + return fitted_x_list, folder_name + def _pre_bulk(self, Log_number, bulk_output_file, iset): + folder_name = "Log{:08d}_{:08d}".format(Log_number, iset) + os.makedirs(folder_name, exist_ok=True) + if self.run_scheme == "connect_so": + pass + else: # self.run_scheme == "subprocess": + shutil.copy( + bulk_output_file, os.path.join(folder_name, bulk_output_file.name) + ) + return folder_name + def _replace(self, fitted_x_list, folder_name): - with open(self.surface_template_file, "r") as file_input, open( - os.path.join(folder_name, self.surface_input_file), "w" - ) as file_output: - for line in file_input: - for index in range(self.dimension): - if line.find(self.string_list[index]) != -1: - line = line.replace( - self.string_list[index], - fitted_x_list[index], - ) + template_file = [] + if self.run_scheme == "subprocess": + file_output = open( + os.path.join(folder_name, self.surface_input_file), "w" + ) + for line in self.template_file_origin: + for index in range(self.dimension): + if line.find(self.string_list[index]) != -1: + line = line.replace( + self.string_list[index], fitted_x_list[index] + ) + + if self.run_scheme == "connect_so": + line = line.replace("\t", " ").replace("\n", " ") + line = line.encode().ljust(self.surf_template_width_for_fortran) + template_file.append(line) + + elif self.run_scheme == "subprocess": file_output.write(line) + if self.run_scheme == "connect_so": + self.template_file = np.array(template_file) + elif self.run_scheme == "subprocess": + file_output.close() + def _check_template(self) -> None: found = [False] * self.dimension with open(self.surface_template_file, "r") as file_input: @@ -166,81 +370,117 @@ def _check_template(self) -> None: msg += label raise exception.InputError(msg) - def _pre_bulk(self, Log_number, bulk_output_file, iset): - folder_name = "Log{:08d}_{:08d}".format(Log_number, iset) - os.makedirs(folder_name, exist_ok=True) - shutil.copy( - bulk_output_file, os.path.join(folder_name, bulk_output_file.name) - ) - return folder_name - class Output(object): """ Output manager. """ + mpicomm: Optional["MPI.Comm"] + mpisize: int + mpirank: int + + run_scheme: str + generate_rocking_curve: bool dimension: int - string_list: List[str] - surface_output_file: str normalization: str + weight_type: Optional[str] + spot_weight: List Rfactor_type: str + omega: float + remove_work_dir: bool + string_list: List[str] + reference_first_line: int + reference_last_line: Optional[int] + reference_path: str + exp_number: List + I_reference_normalized_l: np.ndarray + surface_output_file: str calculated_first_line: int calculated_last_line: int - row_number: int - degree_max: float - degree_list: List[float] + calculated_info_line: int + cal_number: List + + def __init__(self, info, isLogmode, detail_timer): + self.mpicomm = mpi.comm() + self.mpisize = mpi.size() + self.mpirank = mpi.rank() - reference: List[float] - reference_norm: float - reference_normalized: List[float] - degree_list: List[float] + self.isLogmode = isLogmode + self.detail_timer = detail_timer - def __init__(self, info): if "dimension" in info.solver: self.dimension = info.solver["dimension"] else: self.dimension = info.base["dimension"] info_s = info.solver + self.run_scheme = info_s["run_scheme"] - # solver.config - info_config = info_s.get("config", {}) - self.surface_output_file = info_config.get( - "surface_output_file", "surf-bulkP.s" - ) - - v = info_config.get("calculated_first_line", 5) - if not (isinstance(v, int) and v >= 0): - raise exception.InputError( - "ERROR: calculated_first_line should be non-negative integer" - ) - self.calculated_first_line = v - - v = info_config.get("calculated_last_line", 60) - if not (isinstance(v, int) and v >= 0): - raise exception.InputError( - "ERROR: calculated_last_line should be non-negative integer" - ) - self.calculated_last_line = v - - v = info_config.get("row_number", 8) - if not (isinstance(v, int) and v >= 0): - raise exception.InputError( - "ERROR: row_number should be non-negative integer" - ) - self.row_number = v + # If self.run_scheme == "connect_so", + # the contents of surface_output_file are retailned in self.surf_output. + self.surf_output = np.array([]) + self.generate_rocking_curve = info_s.get("generate_rocking_curve", False) # solver.post info_post = info_s.get("post", {}) - - v = info_post.get("normalization", "TOTAL") - if v not in ["TOTAL", "MAX"]: - raise exception.InputError("ERROR: normalization must be TOTAL or MAX") + v = info_post.get("normalization", "") + available_normalization = ["TOTAL", "MANY_BEAM"] + if v == "MAX": + raise exception.InputError( + 'ERROR: solver.post.normalization == "MAX" is no longer available' + ) + if v not in available_normalization: + msg = "ERROR: solver.post.normalization must be " + msg += "MANY_BEAM or TOTAL" + raise exception.InputError(msg) self.normalization = v + v = info_post.get("weight_type", None) + available_weight_type = ["calc", "manual"] + if self.normalization == "MANY_BEAM": + if v is None: + msg = 'ERROR: If solver.post.normalization = "MANY_BEAM", ' + msg += '"weight_type" must be set in [solver.post].' + raise exception.InputError(msg) + elif v not in available_weight_type: + msg = "ERROR: solver.post.weight_type must be " + msg += "calc or manual" + raise exception.InputError(msg) + else: + if v is not None: + if self.mpirank == 0: + msg = 'NOTICE: If solver.post.normalization = "MANY_BEAM" is not set, ' + msg += '"solver.post.weight_type" is NOT considered in the calculation.' + print(msg) + self.weight_type = None + self.weight_type = v + + v = info_post.get("spot_weight", []) + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": + if v == []: + msg = 'ERROR: With solver.post.normalization="MANY_BEAM" and ' + msg += 'solver.post.weight_type=="manual", the "spot_weight" option ' + msg += "must be set in [solver.post]." + raise exception.InputError(msg) + self.spot_weight = np.array(v) / sum(v) + else: + if v != []: + if self.mpirank == 0: + msg = 'NOTICE: With the specified "solver.post.normalization" option, ' + msg += 'the "spot_weight" you specify in the toml file is ignored, ' + msg += "since the spot_weight is automatically calculated within the program." + print(msg) + if self.normalization == "TOTAL": + self.spot_weight = np.array([1]) + v = info_post.get("Rfactor_type", "A") - if v not in ["A", "B"]: - raise exception.InputError("ERROR: Rfactor_type must be A or B") + if v not in ["A", "B", "A2"]: + raise exception.InputError("ERROR: solver.post.Rfactor_type must be A, A2 or B") + if self.normalization == "MANY_BEAM": + if (v != "A") and (v != "A2"): + msg = 'With solver.post.normalization="MANY_BEAM", ' + msg += 'only solver.post.Rfactor_type="A" or "A2" is valid.' + raise exception.InputError(msg) self.Rfactor_type = v v = info_post.get("omega", 0.5) @@ -259,49 +499,202 @@ def __init__(self, info): ) self.string_list = v - v = info_param.get("degree_max", 6.0) - self.degree_max = v - # solver.reference info_ref = info_s.get("reference", {}) - reference_path = info_ref.get("path", "experiment.txt") - - v = info_ref.setdefault("first", 1) + v = info_ref.setdefault("reference_first_line", 1) if not (isinstance(v, int) and v >= 0): raise exception.InputError( "ERROR: reference_first_line should be non-negative integer" ) firstline = v - v = info_ref.setdefault("last", 56) - if not (isinstance(v, int) and v >= firstline): + # None is dummy value + # If "reference_last_line" is not specified in the toml file, + # the last line of the reference file is used for the R-factor calculation. + v = info_ref.setdefault("reference_last_line", None) + if v is None: + reference_are_read_to_final_line = True + else: + reference_are_read_to_final_line = False + if not (isinstance(v, int) and v >= firstline): + raise exception.InputError( + "ERROR: reference_last_line < reference_first_line" + ) + lastline = v + + reference_path = info_ref.get("path", "experiment.txt") + data_experiment = self.read_experiment( + reference_path, firstline, lastline, reference_are_read_to_final_line + ) + self.angle_number_experiment = data_experiment.shape[0] + self.beam_number_exp_raw = data_experiment.shape[1] + + v = info_ref.get("exp_number", []) + + if len(v) == 0: + raise exception.InputError("ERROR: You have to set the 'solver.reference.exp_number'.") + + if not isinstance(v, list): + raise exception.InputError("ERROR: 'solver.reference.exp_number' must be a list type.") + + if max(v) > self.beam_number_exp_raw: + raise exception.InputError("ERROR: The 'solver.reference.exp_number' setting is wrong.") + + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": + if len(v) != len(self.spot_weight): + raise exception.InputError( + "ERROR:len('solver.reference.exp_number') and len('solver.post.spot_weight') do not match." + ) + + if self.normalization == "TOTAL" and len(v) != 1: + msg = 'When solver.post.normalization=="TOTAL" is specified, ' + msg += "only one beam data can be specified with " + msg += '"solver.reference.exp_number" option.' + raise exception.InputError(msg) + + self.exp_number = v + + # Normalization of reference data + self.beam_number_experiment = len(self.exp_number) + for loop_index in range(self.beam_number_experiment): + exp_index = self.exp_number[loop_index] + I_reference = data_experiment[:, exp_index] + if self.normalization == "TOTAL": + I_reference_norm = np.sum(I_reference) + I_reference_normalized = I_reference / I_reference_norm + I_reference_norm_l = np.array([I_reference_norm]) + self.I_reference_normalized_l = np.array([I_reference_normalized]) + elif self.normalization == "MANY_BEAM" and self.weight_type == "calc": + I_reference_norm = np.sum(I_reference) + I_reference_normalized = I_reference / I_reference_norm + if loop_index == 0: # first loop + I_reference_norm_l = np.array([I_reference_norm]) + self.I_reference_normalized_l = np.array( + [I_reference_normalized] + ) + else: # N-th loop + I_reference_norm_l = np.block( + [I_reference_norm_l, I_reference_norm] + ) + self.I_reference_normalized_l = np.block( + [[self.I_reference_normalized_l], [I_reference_normalized]] + ) + elif self.normalization == "MANY_BEAM" and self.weight_type == "manual": + I_reference_norm = np.sum(I_reference) + I_reference_normalized = I_reference / I_reference_norm + if loop_index == 0: # first loop + I_reference_norm_l = np.array([I_reference_norm]) + self.I_reference_normalized_l = np.array( + [I_reference_normalized] + ) + else: # N-th loop + I_reference_norm_l = np.block( + [I_reference_norm_l, I_reference_norm] + ) + self.I_reference_normalized_l = np.block( + [[self.I_reference_normalized_l], [I_reference_normalized]] + ) + else: + msg = "ERROR: solver.post.normalization must be " + msg += "MANY_BEAM or TOTAL" + raise exception.InputError(msg) + + # solver.config + info_config = info_s.get("config", {}) + self.surface_output_file = info_config.get( + "surface_output_file", "surf-bulkP.s" + ) + + v = info_config.get("calculated_first_line", 5) + if not (isinstance(v, int) and v >= 0): raise exception.InputError( - "ERROR: reference_last_line < reference_first_line" + "ERROR: solver.config.calculated_first_line should be non-negative integer" ) - lastline = v + self.calculated_first_line = v + + v = info_config.get( + "calculated_last_line", + self.calculated_first_line + self.angle_number_experiment - 1, + ) + if not (isinstance(v, int) and v >= 0): + raise exception.InputError( + "ERROR: solver.config.calculated_last_line should be non-negative integer" + ) + self.calculated_last_line = v + + # Number of lines in the computation file + # used for R-factor calculations. + self.calculated_nlines = ( + self.calculated_last_line - self.calculated_first_line + 1 + ) + + if self.angle_number_experiment != self.calculated_nlines: + raise exception.InputError( + "ERROR: The number of glancing angles in the calculation data does not match the number of glancing angles in the experimental data." + ) + + # Variable indicating whether the warning + # "self.calculated_nlines does not match + # the number of glancing angles in the calculated file" + # is displayed. + self.isWarning_calcnline = False + + v = info_config.get("calculated_info_line", 2) + if not (isinstance(v, int) and v >= 0): + raise exception.InputError( + "ERROR: solver.config.calculated_info_line should be non-negative integer" + ) + self.calculated_info_line = v + + v = info_config.get("cal_number", []) + if len(v) == 0: + raise exception.InputError("ERROR: You have to set the 'solver.config.cal_number'.") - # Read experiment-data - nline = lastline - firstline + 1 - self.degree_list = [] - self.reference = [] - with open(reference_path, "r") as fp: - for _ in range(firstline - 1): - fp.readline() - for _ in range(nline): - line = fp.readline() - words = line.split() - self.degree_list.append(float(words[0])) - self.reference.append(float(words[1])) - - self.reference_norm = 0.0 - if self.normalization == "TOTAL": - self.reference_norm = sum(self.reference) - else: # self.normalization == "MAX": - self.reference_norm = max(self.reference) - - self.reference_normalized = [ - I_exp / self.reference_norm for I_exp in self.reference - ] + if not isinstance(v, list): + raise exception.InputError("ERROR: 'solver.config.cal_number' must be a list type.") + + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": + if len(self.spot_weight) != len(v): + raise exception.InputError( + "len('solver.config.cal_number') and len('solver.post.spot_weight') do not match." + ) + if self.normalization == "TOTAL" and len(v) != 1: + msg = 'When solver.post.normalization=="TOTAL" is specified, ' + msg += "only one beam data can be specified with " + msg += '"solver.config.cal_number" option.' + raise exception.InputError(msg) + self.cal_number = v + + def read_experiment(self, ref_path, first, last, read_to_final_line): + # Read experiment data + if self.mpirank == 0: + file_input = open(ref_path, "r") + Elines = file_input.readlines() + + firstline = first + if read_to_final_line: + lastline = len(Elines) + else: + lastline = last + + n_exp_row = lastline - firstline + 1 + + # get value from data + for row_index in range(n_exp_row): + line_index = firstline + row_index - 1 + line = Elines[line_index] + data = line.split() + # first loop + if row_index == 0: + n_exp_column = len(data) + data_e = np.zeros((n_exp_row, n_exp_column)) # init value + + for column_index in range(n_exp_column): + data_e[row_index, column_index] = float(data[column_index]) + else: + data_e = None + data_exp = self.mpicomm.bcast(data_e, root=0) + return data_exp def prepare(self, fitted_x_list): self.fitted_x_list = fitted_x_list @@ -309,166 +702,314 @@ def prepare(self, fitted_x_list): def get_results(self, work_dir) -> float: """ Get Rfactor obtained by the solver program. - Returns ------- """ - - cwd = os.getcwd() # Calculate Rfactor and Output numerical results + cwd = os.getcwd() os.chdir(work_dir) Rfactor = self._post(self.fitted_x_list) os.chdir(cwd) - if self.remove_work_dir: - def rmtree_error_handler(function, path, excinfo): - print(f"WARNING: Failed to remove a working directory, {path}") + # delete Log-directory + if self.isLogmode: + time_sta = time.perf_counter() + + if self.run_scheme == "subprocess": + if self.remove_work_dir: + + def rmtree_error_handler(function, path, excinfo): + print(f"WARNING: Failed to remove a working directory, {path}") + + shutil.rmtree(work_dir, onerror=rmtree_error_handler) - shutil.rmtree(work_dir, onerror=rmtree_error_handler) + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["delete_Log-directory"] += time_end - time_sta return Rfactor def _post(self, fitted_x_list): - degree_list = self.degree_list - I_experiment_norm = self.reference_norm - I_experiment_list = self.reference + I_experiment_normalized_l = self.I_reference_normalized_l ( - convolution_I_calculated_list_normalized, - I_calculated_norm, - I_calculated_list, - convolution_I_calculated_list, + glancing_angle, + conv_number_of_g_angle, + conv_I_calculated_norm_l, + conv_I_calculated_normalized_l, ) = self._calc_I_from_file() - Rfactor = self._calc_Rfactor(convolution_I_calculated_list_normalized) - - print("R-factor =", Rfactor) + if self.isLogmode: + time_sta = time.perf_counter() + Rfactor = self._calc_Rfactor( + conv_number_of_g_angle, + conv_I_calculated_normalized_l, + I_experiment_normalized_l, + ) + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["calculate_R-factor"] += time_end - time_sta + # generate RockingCurve_calculated.txt dimension = self.dimension string_list = self.string_list - - with open("RockingCurve.txt", "w") as file_RC: - I_experiment_list_normalized = self.reference_normalized - # Write headers - file_RC.write("#") - for index in range(dimension): + cal_number = self.cal_number + spot_weight = self.spot_weight + weight_type = self.weight_type + Rfactor_type = self.Rfactor_type + normalization = self.normalization + if self.generate_rocking_curve: + if self.isLogmode: + time_sta = time.perf_counter() + with open("RockingCurve_calculated.txt", "w") as file_RC: + # Write headers + file_RC.write("#") + for index in range(dimension): + file_RC.write( + "{} = {} ".format(string_list[index], fitted_x_list[index]) + ) + file_RC.write("\n") + file_RC.write("#Rfactor_type = {}\n".format(Rfactor_type)) + file_RC.write("#normalization = {}\n".format(normalization)) + if weight_type is not None: + file_RC.write("#weight_type = {}\n".format(weight_type)) + file_RC.write("#fx(x) = {}\n".format(Rfactor)) + file_RC.write("#cal_number = {}\n".format(cal_number)) + file_RC.write("#spot_weight = {}\n".format(spot_weight)) file_RC.write( - "{} = {} ".format(string_list[index], fitted_x_list[index]) + "#NOTICE : Intensities are NOT multiplied by spot_weight.\n" ) - file_RC.write("\n") - file_RC.write("#R-factor = {}\n".format(Rfactor)) - if self.normalization == "TOTAL": - file_RC.write("#I_calculated_total={}\n".format(I_calculated_norm)) - file_RC.write("#I_experiment_total={}\n".format(I_experiment_norm)) - else: # self.normalization == "MAX" - file_RC.write("#I_calculated_max={}\n".format(I_calculated_norm)) - file_RC.write("#I_experiment_max={}\n".format(I_experiment_norm)) - file_RC.write("#") - for xname in ( - "degree", - "convolution_I_calculated", - "I_experiment", - "convolution_I_calculated(normalized)", - "I_experiment(normalized)", - "I_calculated", - ): - file_RC.write(xname) - file_RC.write(" ") - file_RC.write("\n") - - # Write rocking curve - for index in range(len(degree_list)): file_RC.write( - "{} {} {} {} {} {}\n".format( - degree_list[index], - convolution_I_calculated_list[index], - I_experiment_list[index], - convolution_I_calculated_list_normalized[index], - I_experiment_list_normalized[index], - I_calculated_list[index], - ) + "#The intensity I_(spot) for each spot is normalized as in the following equation.\n" + ) + file_RC.write("#sum( I_(spot) ) = 1\n") + file_RC.write("#\n") + + label_column = ["glancing_angle"] + fmt_rc = "%.5f" + for i in range(len(cal_number)): + label_column.append(f"cal_number={self.cal_number[i]}") + fmt_rc += " %.15e" + + for i in range(len(label_column)): + file_RC.write(f"# #{i} {label_column[i]}") + file_RC.write("\n") + g_angle_for_rc = np.array([glancing_angle]) + np.savetxt( + file_RC, + np.block([g_angle_for_rc.T, conv_I_calculated_normalized_l.T]), + fmt=fmt_rc, ) - return Rfactor - def _g(self, x): - g = (0.939437 / self.omega) * np.exp( - -2.77259 * (x ** 2.0 / self.omega ** 2.0) - ) - return g + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["make_RockingCurve.txt"] += time_end - time_sta + return Rfactor def _calc_I_from_file(self): + if self.isLogmode: + time_sta = time.perf_counter() + surface_output_file = self.surface_output_file calculated_first_line = self.calculated_first_line calculated_last_line = self.calculated_last_line - row_number = self.row_number - degree_max = self.degree_max - degree_list = self.degree_list - - nlines = calculated_last_line - calculated_first_line + 1 - # TODO: handling error - assert 0 < nlines - - # TODO: nlines == len(degree_list) ? - # assert nlines == len(degree_list) - - I_calculated_list = [] - with open(surface_output_file, "r") as file_result: - for _ in range(calculated_first_line - 1): - file_result.readline() - for _ in range(nlines): - line = file_result.readline() - words = line.replace(",", "").split() - I_calculated_list.append(float(words[row_number - 1])) - # TODO: degree_calc == degree_exp should be checked for every line? - degree_last = round(float(words[0]), 1) - if degree_last != degree_max: - print( - "WARNING : degree in lastline = {}, but {} expected".format( - degree_last, degree_max - ) + calculated_info_line = self.calculated_info_line + calculated_nlines = self.calculated_nlines + omega = self.omega + + cal_number = self.cal_number + + assert 0 < calculated_nlines + + if self.run_scheme == "connect_so": + Clines = self.surf_output + + elif self.run_scheme == "subprocess": + file_input = open(self.surface_output_file, "r") + Clines = file_input.readlines() + file_input.close() + + # Reads the number of glancing angles and beams + line = Clines[calculated_info_line - 1] + line = line.replace(",", "") + data = line.split() + + calc_number_of_g_angles_org = int(data[1]) + calc_number_of_beams_org = int(data[2]) + + if calc_number_of_g_angles_org != calculated_nlines: + if self.mpirank == 0 and not self.isWarning_calcnline: + msg = "WARNING:\n" + msg += "The number of lines in the calculated file " + msg += "that you have set up does not match " + msg += "the number of glancing angles in the calculated file.\n" + msg += "The number of lines (nline) in the calculated file " + msg += "used for the R-factor calculation " + msg += "is determined by the following values\n" + msg += ( + 'nline = "calculated_last_line" - "calculated_first_line" + 1.' ) + print(msg) + self.isWarning_calcnline = True + calc_number_of_g_angles = calculated_nlines + + # Define the array for the original calculated data. + RC_data_org = np.zeros( + (calc_number_of_g_angles, calc_number_of_beams_org + 1) + ) + for g_angle_index in range(calc_number_of_g_angles): + line_index = (calculated_first_line - 1) + g_angle_index + line = Clines[line_index] + line = line.replace(",", "") + data = line.split() + RC_data_org[g_angle_index, 0] = float(data[0]) + for beam_index in range(calc_number_of_beams_org): + RC_data_org[g_angle_index, beam_index + 1] = data[beam_index + 1] + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["load_STR_result"] += time_end - time_sta + + if self.isLogmode: + time_sta = time.perf_counter() + verbose_mode = False + # convolution + data_convolution = lib_make_convolution.calc( + RC_data_org, + calc_number_of_beams_org, + calc_number_of_g_angles, + omega, + verbose_mode, + ) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["convolution"] += time_end - time_sta - ##### convolution ##### - convolution_I_calculated_list = [] - for index in range(len(I_calculated_list)): - integral = 0.0 - degree_org = degree_list[index] - for index2 in range(len(I_calculated_list)): - integral += ( - I_calculated_list[index2] - * self._g(degree_org - degree_list[index2]) - * 0.1 + if self.isLogmode: + time_sta = time.perf_counter() + + angle_number_convolution = data_convolution.shape[0] + if self.angle_number_experiment != angle_number_convolution: + raise exception.InputError( + "ERROR: The number of glancing angles in the calculation data does not match the number of glancing angles in the experimental data." + ) + glancing_angle = data_convolution[:, 0] + + beam_number_reference = len(cal_number) + + # Normalization of calculated data. + for loop_index in range(beam_number_reference): + cal_index = cal_number[loop_index] + conv_I_calculated = data_convolution[:, cal_index] + if self.normalization == "TOTAL": + conv_I_calculated_norm = np.sum(conv_I_calculated) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm ) - convolution_I_calculated_list.append(integral) - - if self.normalization == "TOTAL": - I_calculated_norm = sum(convolution_I_calculated_list) - else: # self.normalization == "MAX" - I_calculated_norm = max(convolution_I_calculated_list) - convolution_I_calculated_list_normalized = [ - c / I_calculated_norm for c in convolution_I_calculated_list - ] + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) + conv_I_calculated_normalized_l = np.array( + [conv_I_calculated_normalized] + ) + elif self.normalization == "MANY_BEAM" and self.weight_type == "calc": + conv_I_calculated_norm = np.sum(conv_I_calculated) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm + ) + if loop_index == 0: # first loop + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) + conv_I_calculated_normalized_l = np.array( + [conv_I_calculated_normalized] + ) + else: # N-th loop + conv_I_calculated_norm_l = np.block( + [conv_I_calculated_norm_l, conv_I_calculated_norm] + ) + conv_I_calculated_normalized_l = np.block( + [ + [conv_I_calculated_normalized_l], + [conv_I_calculated_normalized], + ] + ) + if loop_index == beam_number_reference - 1: # first loop + # calculate spot_weight + self.spot_weight = ( + conv_I_calculated_norm_l / sum(conv_I_calculated_norm_l) + ) ** 2 + elif self.normalization == "MANY_BEAM" and self.weight_type == "manual": + conv_I_calculated_norm = np.sum(conv_I_calculated) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm + ) + if loop_index == 0: # first loop + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) + conv_I_calculated_normalized_l = np.array( + [conv_I_calculated_normalized] + ) + else: # N-th loop + conv_I_calculated_norm_l = np.block( + [conv_I_calculated_norm_l, conv_I_calculated_norm] + ) + conv_I_calculated_normalized_l = np.block( + [ + [conv_I_calculated_normalized_l], + [conv_I_calculated_normalized], + ] + ) + else: + msg = "ERROR: solver.post.normalization must be " + msg += "MANY_BEAM or TOTAL" + raise exception.InputError(msg) + + if self.isLogmode: + time_end = time.perf_counter() + self.detail_timer["normalize_calc_I"] += time_end - time_sta + return ( - convolution_I_calculated_list_normalized, - I_calculated_norm, - I_calculated_list, - convolution_I_calculated_list, + glancing_angle, + angle_number_convolution, + conv_I_calculated_norm_l, + conv_I_calculated_normalized_l, ) - def _calc_Rfactor(self, calc_result): - I_experiment_list_normalized = self.reference_normalized + def _calc_Rfactor(self, n_g_angle, calc_result, exp_result): + spot_weight = self.spot_weight + n_spot = len(spot_weight) if self.Rfactor_type == "A": R = 0.0 - for I_exp, I_calc in zip(I_experiment_list_normalized, calc_result): - R += (I_exp - I_calc) ** 2 + for spot_index in range(n_spot): + R_spot = 0.0 + for angle_index in range(n_g_angle): + R_spot += ( + exp_result[spot_index, angle_index] + - calc_result[spot_index, angle_index] + ) ** 2 + R_spot = spot_weight[spot_index] * R_spot + R += R_spot R = np.sqrt(R) + elif self.Rfactor_type == "A2": + R = 0.0 + for spot_index in range(n_spot): + R_spot = 0.0 + for angle_index in range(n_g_angle): + R_spot += ( + exp_result[spot_index, angle_index] + - calc_result[spot_index, angle_index] + ) ** 2 + R_spot = spot_weight[spot_index] * R_spot + R += R_spot else: # self.Rfactor_type == "B" - # Pendry's R-factor + all_exp_result = [] + all_calc_result = [] + for spot_index in range(n_spot): + for angle_index in range(n_g_angle): + all_exp_result.append(exp_result[spot_index, angle_index]) + all_calc_result.append(calc_result[spot_index, angle_index]) y1 = 0.0 y2 = 0.0 y3 = 0.0 - for I_exp, I_calc in zip(I_experiment_list_normalized, calc_result): + for I_exp, I_calc in zip(all_exp_result, all_calc_result): y1 += (I_exp - I_calc) ** 2 - y2 += I_exp ** 2 - y3 += I_calc ** 2 + y2 += I_exp**2 + y3 += I_calc**2 R = y1 / (y2 + y3) return R diff --git a/src/py2dmat/util/limitation.py b/src/py2dmat/util/limitation.py new file mode 100644 index 00000000..5a1f03e4 --- /dev/null +++ b/src/py2dmat/util/limitation.py @@ -0,0 +1,42 @@ +from abc import ABCMeta, abstractmethod + +import numpy as np + +import py2dmat + +# for type hints +from pathlib import Path +from typing import List, Optional, Dict + +class LimitationBase(metaclass=ABCMeta): + @abstractmethod + 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] + self.ndim = a.shape[1] + + @abstractmethod + def judge(self, x: np.ndarray) -> bool: + pass + +class Inequality(LimitationBase): + 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 : + judge_formula = [] + for formula_index in range(self.n_formura): + value_formula = 0 + for dim_index in range(self.ndim): + value_formula += self.a[formula_index, dim_index]*x[dim_index] + value_formula += self.b[formula_index] + judge_formula.append(value_formula>0) + judge_result = all(judge_formula) + else: + judge_result = True + return judge_result