Solving Generalized Polyomino Puzzles Using the Ising Model

In the polyomino puzzle, the aim is to fill a finite space using several polyomino pieces with no overlaps or blanks. Because it is an NP-complete combinatorial optimization problem, various probabilistic and approximated approaches have been applied to find solutions. Several previous studies embedded the polyomino puzzle in a QUBO problem, where the original objective function and constraints are transformed into the Hamiltonian function of the simulated Ising model. A solution to the puzzle is obtained by searching for a ground state of Hamiltonian by simulating the dynamics of the multiple-spin system. However, previous methods could solve only tiny polyomino puzzles considering a few combinations because their Hamiltonian designs were not efficient. We propose an improved Hamiltonian design that introduces new constraints and guiding terms to weakly encourage favorable spins and pairs in the early stages of computation. The proposed model solves the pentomino puzzle represented by approximately 2000 spins with >90% probability. Additionally, we extended the method to a generalized problem where each polyomino piece could be used zero or more times and solved it with approximately 100% probability. The proposed method also appeared to be effective for the 3D polycube puzzle, which is similar to applications in fragment-based drug discovery.


Introduction
A polyomino, introduced by Golomb in 1954 [1], is a polygon formed by joining one or more squares edge to edge. Golomb studied the tiling problem using polyominoes, in which the goal of the problem is to cover an infinite or finite plane with replicas of a set of polyominoes [2][3][4]. The "polyomino puzzle" is a finite version of the tiling problem in which the goal is to cover a finite space using a number of replicas of a set of polyominoes with no overlaps or blanks. One of the most well-known polyomino puzzles is the pentomino puzzle, which is composed of twelve different pentomino pieces (formed by joining five squares) covering a 6 × 10 rectangular space, where each piece should be used only once.
The polyomino puzzle is known to be an NP-complete combinatorial optimization problem, and there is no polynomial-time exact algorithm to find its optimal solution. Therefore, various probabilistic and approximated algorithms have been proposed. Although heuristic tree search algorithms are conventionally used in such scenarios in computer science, they have difficulty finding optimal solutions with the desired probability as the scale of the problem increases.
In previous studies [5][6][7][8], the polyomino puzzle has been embedded in a quadratic unconstrained binary optimization (QUBO) problem and solved by simulating the dynamics of a multiple-spin system model inspired by the Ising model [9] in physics. The QUBO problem is a mathematical problem for which the objective is to find a vector of binary variables that minimizes the target quadratic polynomial function without any constraints on the variables. It is closely related to the Ising model with the quadratic Hamiltonian function, which has been thoroughly studied in statistical mechanics for over a century.
In the QUBO representation of the polyomino puzzle, each possible spatial placement of a specific polyomino is represented by a dedicated binary (or continuous) variable called a "spin". Further, any mutually enhancing or repressing relationship between a pair of candidate spatial placements is represented by a connection coefficient between the corresponding spins. The solution to the puzzle is then obtained by searching for a ground state of the Hamiltonian function of the multiple-spin system. This is achieved by updating the system status based on the specific dynamics equations, which guarantees the local minimization of the Hamiltonian in time.
For the QUBO representation of the polyomino puzzle, the design of the Hamiltonian function and the dynamics model have been studied over more than three decades, including by Akiyama et al. [5], Takefuji et al. [6], and Manabe et al. [7]. Recently, the QUBO approach has been gaining increased attention because of the successes of quantum annealing, which has been applied to find global minimum solutions effectively. Eagle et al. [8] solved polyomino puzzles using a quantum annealing simulator and actual hardware DW2000Q (D-Wave Systems Inc., Burnaby, BC, Canada).
In the QUBO approach, the standard pentomino puzzle requires approximately 2000 spins at a minimum to correctly represent all possible placements of pieces, considering their rotation and inversion. However, the models used in previous studies could only deal with lower complexity problems than the pentomino puzzle, such as covering a small rectangular space and limiting the rotation and inversion of pieces. This is because their Hamiltonian designs were not efficient enough and overlooked factors such as the generation of unacceptable small blanks (bubble) where no other polyominoes can be placed. The quality of a Hamiltonian design strongly influences the possibility of reaching a global minimum with the QUBO approach whether using a digital computer or a quantum annealer. Therefore, the development of an improved Hamiltonian design strategy for such an application is necessary.
In this study, we propose a Hamiltonian of the QUBO model that effectively represents the polyomino puzzle with consideration of the bubbles and the introduction of novel guiding terms that weakly enhance favorable spins and pairs of spins in the early stages of calculation. We used the Hopfield neural network model [10][11][12] and evaluated the performance of the proposed Hamiltonian formulation for the standard pentomino puzzle, which is difficult to solve with previous methods (Section 3).
Furthermore, we expanded our method to the generalized version of the polyomino puzzle, in which each polyomino piece can be used zero or more times (not exactly once), allowing the target space to have an arbitrary shape (Section 4), and also accepting mixed polyomino sizes: tromino, tetromino, pentomino, hexomino, etc. (Section 5). These generalization efforts are important to enhance the applicability of this technique and pave the way for future real-world applications such as rapid virtual screening in fragment-based drug discovery, which can be regarded as a complicated version of the generalized polyomino puzzle proposed in this study.

The Ising Model and the QUBO Problem
The Ising model-a simplified version of the Heisenberg model-is a mathematical model of ferromagnetism in statistical mechanics. In the Ising model, spins that can be in either of two states (+1 or −1) are typically located on the lattice points, and there are interactions between adjacent spins. The system has a quadratic Hamiltonian function and the spin status changes such that it minimizes the Hamiltonian spontaneously through the interaction between adjacent spins.
The Ising model is mathematically equivalent to the QUBO problem of finding a binary vector that minimizes the target quadratic polynomial function without any constraints on the variables. Lucas formulated various NP-hard combinatorial optimization problems, such as the maximum cut problem and the traveling salesperson problem, as Ising models [13]. He embedded the objective function and constraints of the combinatorial

Hopfield Neural Network for Solving Combinatorial Optimization Problems
In the field of computer science, mathematical models inspired by the Ising model, such as the Hopfield neural network [10][11][12], have been studied as an exotic computing approach to solve QUBO problems. The Hopfield neural network, proposed by Hopfield in 1982 [10], is similar to the electron spin-glass model and can be used to solve QUBO problems by utilizing the intrinsic characteristics of its quadratic Hamiltonian minimization. In the original model, the output of a neuron, which is a node of the network, had a binary output value of zero or one. In the later model, the output has been extended to have a continuous value from zero to one [11]. Further, the spin status is computed as continuous values, and a binarization operation is performed after reading out the final status. The Hopfield neural network is able to solve combinatorial optimization problems by embedding the objective function and constraints into a QUBO formulation and searching for a ground state of Hamiltonian defined for the network [12]. The procedure merely provides a local minimum of Hamiltonian depending on the initial states and random noise, and thus many trials are required in general to find the ground state. A carefully designed Hamiltonian function that efficiently represents the characteristics of the target problem and provides a more convex energy landscape can significantly increase the probability of reaching the global minimum. An appropriate procedure for updating the system status can also significantly influence the convergence performance for finding the optimal solution in the Hopfield neural network and its derivative algorithms, such as the Gaussian machine [20]. In this study, we evaluated the performance of the proposed method by using a software simulator that implements a Hopfield neural network with continuous output values.
Considering the Hopfield neural network, the energy ground state is searched for by iteratively updating the neurons' output based on specific update rules. Regarding the model with continuous output values, each neuron first calculates its internal value U i using the following differential equation: where n is the number of neurons in the network, W ij is the weight between neurons i and j, V j is the output of another neuron j, and θ i is the bias (i.e., the input to neuron i).
In addition, τ is the time constant parameter, which determines how well the internal values of the neuron are conserved, considering the time variation dt. The Hopfield neural network is an interconnected network; however, it does not contain any self-connections (W ii = 0), and the weight values are symmetric (W ij = W ji ). Thereafter, the output value V i of the neuron i is determined by the following sigmoid function: where T m f is a parameter that determines the gradient of the sigmoid function (also known as the mean-field temperature). A higher T m f makes the sigmoid function closer to a linear function, whereas a lower T m f makes it closer to a step function. The energy of the network in the continuous value model is defined as follows: The discrete value model minimizes the energy E in Equation (4) by repeatedly updating the neurons, whereas the continuous value model minimizes the energy F, which also considers the entropy term S as in Equation (3) [11]. Here, it is possible to search for a solution to a combinatorial optimization problem using the following steps.

1.
Express the objective function to be minimized in a quadratic form for n binary variables.

2.
Compare the quadratic form to be minimized with the energy function E in Equation (4) and obtain the values of W ij and θ i , which make them coincide.

3.
Construct a Hopfield neural network with the obtained values of W ij and θ i .

4.
Select one neuron randomly and update it using Equations (1) and (2). 5. Perform Step 4 several times to minimize the free energy F in Equation (3). 6.
Obtain the output values of the neurons and binarize them based on some criteria. The solution is obtained if the output values of the neurons after the binarization operation satisfy the constraint conditions. In this study, the binarization is performed based on whether the output value is greater than or equal to 0.5.
Regarding this case, for the binarized result of the minimum solution to F to be consistent with the minimum solution to E, the value of T m f must be considerably small at the end. Moreover, there must be no significant difference between the landscape of the basin created by F and that of E. Akiyama et al. proposed a method known as the sharpening technique, in which T m f is large and the gradient of the sigmoid function is gradual in the initial states and becomes progressively steeper over time [5]. This method increases the probability of convergence to the global minimum by gradually transforming the space of F from a state where the ground state is easy to find to a state where there is almost no difference from the original landscape of E. We employ this technique, and the details of T m f scheduling are presented in Section 3.3.1. Figure 1 shows a pentomino puzzle, which is a well-known polyomino puzzle example. Pentominoes are formed by joining five squares, and each of the twelve unique pentominoes is used only once in the pentomino puzzle. Therefore, the total size of the board is 60 squares, and a typical example is a 6 × 10 board, as shown in Figure 1. In this study, pentominoes can rotate every 90 degrees and can be inverted (flipped front and back). Considering these conditions, the pentomino puzzle has 2339 unique solutions [21].

Improvement of the Hamiltonian Function
Applying the Hopfield neural network to the pentomino puzzle, we referred to previous studies to design the neurons; specifically, a neuron corresponds to a single pattern of placement of each pentomino on the board ( Figure 2) [5,22]. If the output of a neuron is one, it indicates that a pentomino is at the spatial corresponding position of the neuron. For instance, there are a total of 296 possible placement patterns on the board, considering the rotation and inversion of the specific pentomino piece shown in Figure 2. Similarly, considering the number of placement patterns for the other eleven types of pentominoes, we observe that a total of 1928 neurons are required to represent all the placement patterns of the 12 pentominoes. Considering the design of the Hamiltonian function, which includes the biases applied to neurons and the connection weights between neurons, we improve the formulation by adding new terms, compared to previous studies [5,22]. Here, we first describe the two conventional constraint terms. Subsequently, we describe four newly proposed terms.

Constraint on the Number of Times Each Pentomino Is Used
All pentominoes must be placed somewhere on the board only once. This constraint can be expressed by the following equation: where is the number of types of pentominoes and is the set of neurons that express all the placement patterns of a pentomino . The first term of Equation (6) becomes zero only when the summation of included in is one for each pentomino . An ideal state in which the output of only a single neuron is one and others are zero for each set of neurons satisfies the condition. The second term is to eliminate (the weight from a neuron back to itself) generated by the first term; therefore, = 0.

Improvement of the Hamiltonian Function
Applying the Hopfield neural network to the pentomino puzzle, we referred to previous studies to design the neurons; specifically, a neuron corresponds to a single pattern of placement of each pentomino on the board ( Figure 2) [5,22]. If the output of a neuron is one, it indicates that a pentomino is at the spatial corresponding position of the neuron. For instance, there are a total of 296 possible placement patterns on the board, considering the rotation and inversion of the specific pentomino piece shown in Figure 2. Similarly, considering the number of placement patterns for the other eleven types of pentominoes, we observe that a total of 1928 neurons are required to represent all the placement patterns of the 12 pentominoes.

Improvement of the Hamiltonian Function
Applying the Hopfield neural network to the pentomino puzzle, we referred to previous studies to design the neurons; specifically, a neuron corresponds to a single pattern of placement of each pentomino on the board ( Figure 2) [5,22]. If the output of a neuron is one, it indicates that a pentomino is at the spatial corresponding position of the neuron. For instance, there are a total of 296 possible placement patterns on the board, considering the rotation and inversion of the specific pentomino piece shown in Figure 2. Similarly, considering the number of placement patterns for the other eleven types of pentominoes, we observe that a total of 1928 neurons are required to represent all the placement patterns of the 12 pentominoes. Considering the design of the Hamiltonian function, which includes the biases applied to neurons and the connection weights between neurons, we improve the formulation by adding new terms, compared to previous studies [5,22]. Here, we first describe the two conventional constraint terms. Subsequently, we describe four newly proposed terms.

Constraint on the Number of Times Each Pentomino Is Used
All pentominoes must be placed somewhere on the board only once. This constraint can be expressed by the following equation: where is the number of types of pentominoes and is the set of neurons that express all the placement patterns of a pentomino . The first term of Equation (6) becomes zero only when the summation of included in is one for each pentomino . An ideal state in which the output of only a single neuron is one and others are zero for each set of neurons satisfies the condition. The second term is to eliminate (the weight from a neuron back to itself) generated by the first term; therefore, = 0. Considering the design of the Hamiltonian function, which includes the biases applied to neurons and the connection weights between neurons, we improve the formulation by adding new terms, compared to previous studies [5,22]. Here, we first describe the two conventional constraint terms. Subsequently, we describe four newly proposed terms.

Constraint on the Number of Times Each Pentomino Is Used
All pentominoes must be placed somewhere on the board only once. This constraint can be expressed by the following equation: where M is the number of types of pentominoes and P m is the set of neurons that express all the placement patterns of a pentomino m. The first term of Equation (6) becomes zero only when the summation of V i included in P m is one for each pentomino m. An ideal state in which the output of only a single neuron is one and others are zero for each set of neurons P m satisfies the condition. The second term is to eliminate W ii (the weight from a neuron back to itself) generated by the first term; therefore, W ii = 0. All pentominoes must be placed with no overlaps or blanks as a solution to the pentomino puzzle. This constraint can be expressed by the following equation: where N is the number of squares constituting the board, and L n is the set of neurons that express all the placement patterns of the pentominoes occupying a square n on the board. The first term of Equation (7) becomes zero only when the summation of V i , included in L n , is one for each square n on the board. An ideal state where the output of only a single neuron is one and the others are zero for each set of neurons L n satisfies this condition. If multiple neurons in a set of neurons corresponding to a specific square have an output value of one, it indicates that pentominoes are overlapping on that square. Furthermore, all neurons in the set have an output value of zero, indicating that the square is left blank. The second term in the equation is to eliminate W ii , which is the same as the second term of the H 1 described in Section 3.2.1.

Inhibition of Bubbles
In addition to the two constraints presented above, we added a new constraint to inhibit bubbles, which are small blanks where there is no room for another pentomino to be inserted, as demonstrated in the shaded areas in Figure 3. Bubbles can be generated by a single pentomino or by a combination of several pentominoes. The placement of a single pentomino that generates a bubble cannot be included in any solution; thus, it can be easily excluded beforehand. In contrast, it is desirable to restrict a pair of neurons that generate bubbles; therefore, the following term is added to the Hamiltonian: where Q is the total number of neurons, and the function f (i, j) represents whether the bubbles are generated or not. If all the pairs in the set of neurons with an output value of one do not generate bubbles, then this term is minimized to zero.
Entropy 2022, 24, 354 6 of 21 All pentominoes must be placed with no overlaps or blanks as a solution to the pentomino puzzle. This constraint can be expressed by the following equation: where is the number of squares constituting the board, and is the set of neurons that express all the placement patterns of the pentominoes occupying a square on the board. The first term of Equation (7) becomes zero only when the summation of , included in , is one for each square on the board. An ideal state where the output of only a single neuron is one and the others are zero for each set of neurons satisfies this condition. If multiple neurons in a set of neurons corresponding to a specific square have an output value of one, it indicates that pentominoes are overlapping on that square. Furthermore, all neurons in the set have an output value of zero, indicating that the square is left blank. The second term in the equation is to eliminate , which is the same as the second term of the described in Section 3.2.1.

Inhibition of Bubbles
In addition to the two constraints presented above, we added a new constraint to inhibit bubbles, which are small blanks where there is no room for another pentomino to be inserted, as demonstrated in the shaded areas in Figure 3. Bubbles can be generated by a single pentomino or by a combination of several pentominoes. The placement of a single pentomino that generates a bubble cannot be included in any solution; thus, it can be easily excluded beforehand. In contrast, it is desirable to restrict a pair of neurons that generate bubbles; therefore, the following term is added to the Hamiltonian: where is the total number of neurons, and the function ( , ) represents whether the bubbles are generated or not. If all the pairs in the set of neurons with an output value of one do not generate bubbles, then this term is minimized to zero.

Encouragement for Combinations of Pentominoes That Contact Each Other
Considering the described in Section 3.2.3, we restrict inappropriate combinations of neurons that generate bubbles. In contrast, it is also possible to improve the convergence performance to the optimal solution by encouraging appropriate combinations of neurons. One of the examples is the appropriate mutual contact of pentominoes. The more the edges of two pentominoes touch, the better the contact between them. Therefore, the following term is added to the Hamiltonian:

Encouragement for Combinations of Pentominoes That Contact Each Other
Considering the H 3 described in Section 3.2.3, we restrict inappropriate combinations of neurons that generate bubbles. In contrast, it is also possible to improve the convergence performance to the optimal solution by encouraging appropriate combinations of neurons. One of the examples is the appropriate mutual contact of pentominoes. The more the edges of two pentominoes touch, the better the contact between them. Therefore, the following term is added to the Hamiltonian: where µ ij is the number of edges of squares shared by the two pentominoes corresponding to the neurons i, j, and λ ij is the number of overlapping squares of the two pentominoes. The function g(i, j) evaluates the goodness of the combination for all pairs of neurons. If the placements of pentominoes corresponding to neurons i, j do not overlap, then the pair is encouraged in proportion to the number of edges they share with each other.

Encouragement for Pentominoes That Contact the Borders of the Board
Regarding the H 4 described in Section 3.2.4, we evaluate the goodness of the combinations of pairs of neurons. Similarly, we can also evaluate the goodness of a single neuron. To weakly encourage pentominoes that contact the borders of the board well, the following term is added to the Hamiltonian: where h i is the number of edges where the placement of a pentomino corresponding to neuron i contacts any borderline of the board.

Broad Restriction on All Pairs of Neurons
The encouragement terms described above will promote the finding of an optimal solution; however, this tends to make output values of several neurons fixed to one at the early stages in the simulation. It reduces the simulation runtime. Nonetheless, it increases the probability of being trapped in a local minimum because the Hopfield neural network is sensitive to the initial states of neurons initialized by uniform random numbers. To avoid this unintentional behavior, all pairs of neurons are equally suppressed. Thus, the following term is added to the Hamiltonian:

Hamiltonian for the Pentomino Puzzle
Combining Equations (6)-(11), the modified Hamiltonian for the pentomino puzzle is proposed as follows: Whereas H 1 to H 3 are constraint terms that describe the hard constraints on the problem, H 4 to H 6 are the guiding terms that work to guide the system to the global minimum. Guiding terms that encourage or restrict neurons or pairs of neurons have a common issue. If the coefficient of the term is too large, the ground state of the free energy function does not necessarily match the state in which the target puzzle is solved. Therefore, we introduce the parameter γ, which is a decaying weight coefficient for guiding terms. The parameter γ is one at the start of the simulation and decays linearly to zero at the end. Therefore, the effect of guiding terms is weakened over the course of the simulation. The coefficients A, B, C, D, E, and F on each term are positive real numbers. These coefficients determine the balance of each term; therefore, it is very important to select the appropriate values.
Comparing Equation (4) to Equation (12), the weights W ij between the neurons and biases θ i are obtained as follows: where a i is the size of the polyomino corresponding to neuron i (a i is five for the pentomino puzzle).

Experimental Parameters
Considering the simulation, we randomly selected one of the n neurons and calculated its output value V i based on the update rule described in Section 2.2. In other words, we chose the asynchronous update of neurons, which is only feasible in software simulations and has a better effect for faster convergence. In this study, we define one time step as performing this update process n times, which is the number of neurons. Because the selection of a neuron was performed randomly, there were neurons that were updated multiple times and some that were never updated during a single step. Moreover, ∆t = 1 and τ = 1 are used in Equation (1), which indicates that the influence of the previous value of U i is completely forgotten in each update. In addition, the output value V i of each neuron is initialized by a uniform random number in the range [0, 2/n] such that the sum of the outputs of the n neurons approaches one. This is because if each output value is too large, such as 0.5, it will take longer for the network to stabilize at the earlier stages of the simulation.
We employed the sharpening technique that decays T m f , which determines the gradient of the sigmoid function over the course of the simulation. The scheduling of T m f is as follows: m f is the initial value of T m f , which is fixed at one in this experiment. Moreover, t is the number of steps in the simulation, and t max is the total number of steps to be specified. The computer experiments presented in this section were performed by varying t max from 500 to 5000.
The coefficients of the proposed Hamiltonian were determined via Bayesian optimization to obtain the values shown in Table 1. The details are described in Section 6.1.  Figure 4 shows the experimental results obtained by changing the total number of simulation steps, t max . Considering each setting of t max , we ran the simulation 1000 times with different random seeds for generating the initial output value of the neurons. Figure 4a shows the distribution of the Hamiltonian energy obtained in the final state of the simulation. As t max increases, the system discovers smaller energy states, indicating better solutions. Figure 4b shows the probability of finding the optimal solution, where the probability is 90.5% at t max = 5000, indicating that the proposed Hamiltonian has a better performance in searching for the solution.
The CPU time for a simulation with t max = 5000 using a single CPU core (Intel Xeon E5-2680 2.4 GHz, Intel Corporation, Santa Clara, CA, USA) was approximately 30 s. The CPU time for a simulation with = 5000 using a single CPU core (Intel Xeon E5-2680 2.4 GHz) was approximately 30 s.   Figure 5 shows the intermediate states of conversion to an optimal solution. Considering any trial, the pentominoes in contact with the borders were placed first, mainly because of the term. Additionally, the term, which restricts the overlap of pentominoes, suppresses the placement to the center of the board where pentominoes are easily overlapped. Thereafter, the placement of pentominoes that have been in contact with the already placed pentominoes can be selected. This is because the term encourages good pairs of placements and when complex spaces are produced by the placement of pentominoes, only specific pentominoes can cover the space properly. Thus, because the proposed method is not formulated in such a manner as to increase the diversity of solutions, it tends to converge to the same optimal solution in multiple trials. When the number of steps = 5000, the optimal solution was found in 905 out of 1000 trials, with only 4 unique solutions.   Figure 5 shows the intermediate states of conversion to an optimal solution. Considering any trial, the pentominoes in contact with the borders were placed first, mainly because of the H 5 term. Additionally, the H 2 term, which restricts the overlap of pentominoes, suppresses the placement to the center of the board where pentominoes are easily overlapped. Thereafter, the placement of pentominoes that have been in contact with the already placed pentominoes can be selected. This is because the H 4 term encourages good pairs of placements and when complex spaces are produced by the placement of pentominoes, only specific pentominoes can cover the space properly. Thus, because the proposed method is not formulated in such a manner as to increase the diversity of solutions, it tends to converge to the same optimal solution in multiple trials. When the number of steps t max = 5000, the optimal solution was found in 905 out of 1000 trials, with only 4 unique solutions.   Figure 5 shows the intermediate states of conversion to an optimal solution. Considering any trial, the pentominoes in contact with the borders were placed first, mainly because of the term. Additionally, the term, which restricts the overlap of pentominoes, suppresses the placement to the center of the board where pentominoes are easily overlapped. Thereafter, the placement of pentominoes that have been in contact with the already placed pentominoes can be selected. This is because the term encourages good pairs of placements and when complex spaces are produced by the placement of pentominoes, only specific pentominoes can cover the space properly. Thus, because the proposed method is not formulated in such a manner as to increase the diversity of solutions, it tends to converge to the same optimal solution in multiple trials. When the number of steps = 5000, the optimal solution was found in 905 out of 1000 trials, with only 4 unique solutions.

Target Problem
The pentomino puzzle presented above can be considered a special case of the generalized problem. In this section, we extend the puzzle in two aspects as a first step:

1.
Each pentomino can be used zero or more times; 2.
An arbitrarily shaped board is the subject of the problem.
Hereinafter, the generalized model allows each polyomino to be used zero or more times. Additionally, we consider boards of any shape (all the squares making up the board need to be interconnected) as a crucial step for approaching real-world applications. Figure 6 shows an example of the generalized polyomino puzzle problem.

Target Problem
The pentomino puzzle presented above can be considered a special case of the generalized problem. In this section, we extend the puzzle in two aspects as a first step: 1. Each pentomino can be used zero or more times; 2. An arbitrarily shaped board is the subject of the problem.
Hereinafter, the generalized model allows each polyomino to be used zero or more times. Additionally, we consider boards of any shape (all the squares making up the board need to be interconnected) as a crucial step for approaching real-world applications. Figure 6 shows an example of the generalized polyomino puzzle problem.

Hamiltonian for the Generalized Polyomino Puzzle
For the generalized polyomino puzzle described in Section 4.1, a modification of the Hamiltonian in Equation (12), which is for pentomino puzzles, needs to be considered. Because none of the terms in Equation (12) depend on the shape of the target board, the same Hamiltonian can be used for the arbitrary shape of the board. In contrast, because the term restricts the number of times each pentomino can be used, this must be removed to place each pentomino zero or more times. Therefore, the Hamiltonian for the generalized polyomino puzzle is introduced as follows: The weights between the neurons and biases are obtained as follows:

Generation of Arbitrarily Shaped Boards
We generated several boards with complex shapes for the generalized polyomino puzzle other than the 6 × 10 board used in the pentomino puzzle. Each board was randomly generated by joining 60 squares without any hole on the board. Figure 7 shows the 101 different boards in an ascending order based on their perimeters, representing their complexity. This consists of the original 6 × 10 board and 100 randomly generated boards. We selected the 6 × 10 board with a perimeter of 32 as the original board (Original), and the board with the longest perimeter of 94 was considered as the most complex board (Complex). Furthermore, the board with a perimeter of 64, which is closest to the average of those two perimeters, was chosen as the board with the average complexity (Middle). The total number of neurons for the Original, Middle, and Complex boards was 1928, 857, and 129, respectively. These three boards are shown in Figure 7.

Hamiltonian for the Generalized Polyomino Puzzle
For the generalized polyomino puzzle described in Section 4.1, a modification of the Hamiltonian in Equation (12), which is for pentomino puzzles, needs to be considered. Because none of the terms in Equation (12) depend on the shape of the target board, the same Hamiltonian can be used for the arbitrary shape of the board. In contrast, because the H 1 term restricts the number of times each pentomino can be used, this H 1 must be removed to place each pentomino zero or more times. Therefore, the Hamiltonian for the generalized polyomino puzzle is introduced as follows: The weights W ij between the neurons and biases θ i are obtained as follows:

Generation of Arbitrarily Shaped Boards
We generated several boards with complex shapes for the generalized polyomino puzzle other than the 6 × 10 board used in the pentomino puzzle. Each board was randomly generated by joining 60 squares without any hole on the board. Figure 7 shows the 101 different boards in an ascending order based on their perimeters, representing their complexity. This consists of the original 6 × 10 board and 100 randomly generated boards. We selected the 6 × 10 board with a perimeter of 32 as the original board (Original), and the board with the longest perimeter of 94 was considered as the most complex board (Complex). Furthermore, the board with a perimeter of 64, which is closest to the average of those two perimeters, was chosen as the board with the average complexity (Middle). The total number of neurons for the Original, Middle, and Complex boards was 1928, 857, and 129, respectively. These three boards are shown in Figure 7.

Experimental Parameters
The experiments were performed with the total number of steps changing from 50 to 500. The steps were set shorter than in the experiment described in Section 3.3.1 because the number of simulation steps required to obtain an optimal solution was smaller. We used Bayesian optimization to determine the coefficients of the proposed Hamiltonian function, shown in Table 2. For the other parameters, we used the same settings as in Section 3.3.1.  Figure 8 shows the probability of finding an optimal solution for various . Similar to Section 3.3.2, considering each , we ran the simulation 1000 times with different random seeds for generating the initial output value of the neurons. More than 99% of the trials reached an optimal solution at = 200 for each board, indicating that the proposed Hamiltonian had a sufficient performance for these problems.
Moreover, the probability of the Original board was substantially higher than the standard pentomino puzzle shown in Section 3.3.2, even with fewer . The extension of the number of times each pentomino is used can be considered as the drastic relaxation

Experimental Parameters
The experiments were performed with the total number of steps t max changing from 50 to 500. The steps were set shorter than in the experiment described in Section 3.3.1 because the number of simulation steps required to obtain an optimal solution was smaller. We used Bayesian optimization to determine the coefficients of the proposed Hamiltonian function, shown in Table 2. For the other parameters, we used the same settings as in Section 3.3.1.  Figure 8 shows the probability of finding an optimal solution for various t max . Similar to Section 3.3.2, considering each t max , we ran the simulation 1000 times with different random seeds for generating the initial output value of the neurons. More than 99% of the trials reached an optimal solution at t max = 200 for each board, indicating that the proposed Hamiltonian had a sufficient performance for these problems.
Moreover, the probability of the Original board was substantially higher than the standard pentomino puzzle shown in Section 3.3.2, even with fewer t max . The extension of the number of times each pentomino is used can be considered as the drastic relaxation of the problem, resulting in easier achievement. In addition, the overall probability of finding an optimal solution in the Middle and Complex boards was higher than that in the Original. This is because of their complex shapes, where the total number of neurons is lower because there are areas with fewer candidates suitable to cover the space. Moreover, the output values of certain neurons are fixed to one in the early stages of the simulation, narrowing the search space.
Entropy 2022, 24, 354 12 of 21 of the problem, resulting in easier achievement. In addition, the overall probability of finding an optimal solution in the Middle and Complex boards was higher than that in the Original. This is because of their complex shapes, where the total number of neurons is lower because there are areas with fewer candidates suitable to cover the space. Moreover, the output values of certain neurons are fixed to one in the early stages of the simulation, narrowing the search space.

Set of Polyominoes
Popular polyomino puzzles use polyominoes of the same size. However, some problems deal with various sizes of polyominoes, such as the example of a commercialized  of the problem, resulting in easier achievement. In addition, the overall probability of finding an optimal solution in the Middle and Complex boards was higher than that in the Original. This is because of their complex shapes, where the total number of neurons is lower because there are areas with fewer candidates suitable to cover the space. Moreover, the output values of certain neurons are fixed to one in the early stages of the simulation, narrowing the search space.

Set of Polyominoes
Popular polyomino puzzles use polyominoes of the same size. However, some problems deal with various sizes of polyominoes, such as the example of a commercialized

Set of Polyominoes
Popular polyomino puzzles use polyominoes of the same size. However, some problems deal with various sizes of polyominoes, such as the example of a commercialized polyomino puzzle, where an 8 × 8 board is covered with twelve pentominoes and one tetromino. Figure 10 is an example of the problem discussed here. In addition to the extensions described in Section 4, we further extend the puzzle to deal with polyominoes of various sizes simultaneously as a subsequent step. However, covering a board with small polyominoes is generally easier than with large polyominoes. Regarding some industrial applications, such as fragment-based drug discovery, it is desirable to cover the space by combining a small number of larger and functional components, rather than simply filling it up with several tiny components. Therefore, we consider the problem of minimizing the total number of polyominoes used by enhancing the use of larger polyominoes.
Entropy 2022, 24, 354 13 of 21 polyomino puzzle, where an 8 × 8 board is covered with twelve pentominoes and one tetromino. Figure 10 is an example of the problem discussed here. In addition to the extensions described in Section 4, we further extend the puzzle to deal with polyominoes of various sizes simultaneously as a subsequent step. However, covering a board with small polyominoes is generally easier than with large polyominoes. Regarding some industrial applications, such as fragment-based drug discovery, it is desirable to cover the space by combining a small number of larger and functional components, rather than simply filling it up with several tiny components. Therefore, we consider the problem of minimizing the total number of polyominoes used by enhancing the use of larger polyominoes.  Table 3 shows the number of types of polyominoes for each size. Polyominoes that have the same shape when rotated or inverted are treated as the same type of polyomino. In addition to the twelve types of pentominoes discussed in the previous sections, we introduced new polyominoes of various sizes. However, if we include tiny polyominoes such as the monomino (size 1) and the domino (size 2) in the candidate set, there are trivial solutions, such as covering the board by simply repeating them. Therefore, only the tromino (size 3) and the tetromino (size 4) were included as small-size polyominoes. We also included the hexomino (size 6) as larger polyominoes. We have not included heptomino (size 7) and those larger than that because they have a large number of shapes and the total number of neurons required will be too large. The total number of types of polyominoes used was 54. When all the 54 polyominoes were used, the total number of neurons needed for the Original, Middle, and Complex boards was 8469, 4141, and 693, respectively.   Table 3 shows the number of types of polyominoes for each size. Polyominoes that have the same shape when rotated or inverted are treated as the same type of polyomino. In addition to the twelve types of pentominoes discussed in the previous sections, we introduced new polyominoes of various sizes. However, if we include tiny polyominoes such as the monomino (size 1) and the domino (size 2) in the candidate set, there are trivial solutions, such as covering the board by simply repeating them. Therefore, only the tromino (size 3) and the tetromino (size 4) were included as small-size polyominoes. We also included the hexomino (size 6) as larger polyominoes. We have not included heptomino (size 7) and those larger than that because they have a large number of shapes and the total number of neurons required will be too large. The total number of types of polyominoes used was 54. When all the 54 polyominoes were used, the total number of neurons needed for the Original, Middle, and Complex boards was 8469, 4141, and 693, respectively.

Experimental Parameters
Considering the generalized polyomino puzzle with various sizes of polyominoes, we used the same Hamiltonian as in Section 4 because the H 6 term could reduce the total number of polyominoes to be used. A detailed discussion is presented in Section 5.2.3. The experiments were performed by varying the total number of steps t max from 50 to 500. We used Bayesian optimization to determine the coefficients of the proposed Hamiltonian, which are shown in Table 4. For the other parameters, we used the same settings as in Section 3.3.1.  Figure 11 shows the probability of finding an optimal solution and the breakdown of the total number of polyominoes used in each solution. Considering each t max , we ran the simulation 1000 times with different random seeds to generate the initial output value of the neurons. There was a little difference in the probability of finding an optimal solution among the three boards, and the probabilities were reached at approximately 100% at t max = 200 or larger for all the boards. The probability of finding a solution does not differ significantly compared to the results using only the pentominoes shown in Figure 8. This is because although the total number of required neurons was increasing, the total number of possible solutions was also increasing simultaneously. In addition, considering the breakdown of the total number of polyominoes used in the obtained solutions, most of the solutions in the Original and Middle boards covered the 60 squares using only 10 hexominoes (polyominoes of size 6), minimizing the total number of polyominoes. There was no solution that covered the Complex board with ten polyominoes; therefore, the solutions with eleven polyominoes were obtained.

Experimental Parameters
Considering the generalized polyomino puzzle with various sizes of polyominoes, we used the same Hamiltonian as in Section 4 because the term could reduce the total number of polyominoes to be used. A detailed discussion is presented in Section 5.2.3. The experiments were performed by varying the total number of steps from 50 to 500. We used Bayesian optimization to determine the coefficients of the proposed Hamiltonian, which are shown in Table 4. For the other parameters, we used the same settings as in Section 3.3.1.  Figure 11 shows the probability of finding an optimal solution and the breakdown of the total number of polyominoes used in each solution. Considering each , we ran the simulation 1000 times with different random seeds to generate the initial output value of the neurons. There was a little difference in the probability of finding an optimal solution among the three boards, and the probabilities were reached at approximately 100% at = 200 or larger for all the boards. The probability of finding a solution does not differ significantly compared to the results using only the pentominoes shown in Figure  8. This is because although the total number of required neurons was increasing, the total number of possible solutions was also increasing simultaneously. In addition, considering the breakdown of the total number of polyominoes used in the obtained solutions, most of the solutions in the Original and Middle boards covered the 60 squares using only 10 hexominoes (polyominoes of size 6), minimizing the total number of polyominoes. There was no solution that covered the Complex board with ten polyominoes; therefore, the solutions with eleven polyominoes were obtained. Figure 11. Probability of finding an optimal solution and the breakdown of the number of placed polyominoes of the optimal solution for various . Figure 11. Probability of finding an optimal solution and the breakdown of the number of placed polyominoes of the optimal solution for various t max .

Size of the Placed Polyominoes
In the previous section, we showed that the total number of polyominoes used to cover the board was minimized by using the same Hamiltonian as in Section 4.4.1, which was not an intuitive result because the effect of the H 6 term suppresses all the pairs of neurons universally. The H 6 term penalty increases when the number of neurons to output one is increasing, indicating the use of more polyominoes. Therefore, the use of fewer and larger polyominoes was encouraged to minimize the Hamiltonian energy, and the total number of polyominoes was minimized. Figure 12 shows the breakdown of the size of the placed polyominoes when the coefficient F for the H 6 term varies from the optimum value obtained shown in Table 4. At F = 0, a large proportion of trominoes and tetrominoes are observed. However, as coefficient F increases, the proportion of such small-sized polyominoes decreases, and at the setting of F = 1, hexominoes are mostly used. The same is true for F > 1. For the Complex board, a certain percentage of tetrominoes and pentominoes were used, even at F = 1 or larger, because the board contained areas that were not coverable only with hexominoes.
cover the board was minimized by using the same Hamiltonian as in Section 4.4.1, which was not an intuitive result because the effect of the term suppresses all the pairs of neurons universally. The term penalty increases when the number of neurons to output one is increasing, indicating the use of more polyominoes. Therefore, the use of fewer and larger polyominoes was encouraged to minimize the Hamiltonian energy, and the total number of polyominoes was minimized. Figure 12 shows the breakdown of the size of the placed polyominoes when the coefficient for the term varies from the optimum value obtained shown in Table 4. At = 0, a large proportion of trominoes and tetrominoes are observed. However, as coefficient increases, the proportion of such small-sized polyominoes decreases, and at the setting of = 1, hexominoes are mostly used. The same is true for > 1. For the Complex board, a certain percentage of tetrominoes and pentominoes were used, even at = 1 or larger, because the board contained areas that were not coverable only with hexominoes. Figure 12. Breakdown of the size of the placed polyominoes for various . It includes the results of all the trials, not only those in which an optimal solution is found. Moreover, of all the trials was fixed at 500, which was sufficiently large, considering the probability of finding an optimal solution.

Performance with the Reduced Types of Polyominoes
When all the 54 types of polyominoes were available, an optimal solution was found with a high probability and a small number of simulation steps. Therefore, to verify the performance when the number of available polyomino types becomes smaller, we performed experiments by randomly removing the polyominoes from the set of polyominoes of the basic 54 types. Figure 13 shows the probability of finding an optimal solution and the breakdown of the number of placed polyominoes for the number of types of polyominoes. Furthermore, is fixed at 500, which was considered large enough (for reference, the results for various are shown in Figure S1 in the Supplementary Material). Because the shape of the board became more complicated from Original to Complex, the probability of finding an optimal solution decreased more rapidly, when gradually removing the polyominoes from the set. For the Middle and Complex boards, it was almost impossible to obtain any optimal solution with randomly selected ten types of polyominoes from the original set. This is because a complex local part of a board tends to have few candidate polyominoes to cover it, and all such candidates may have been removed. In contrast, considering the Original board, it was possible to cover the board with fewer types of polyominoes. Even when we randomly selected only ten types of polyominoes, the optimal solution was found with a probability of approximately 40%. Figure 12. Breakdown of the size of the placed polyominoes for various F. It includes the results of all the trials, not only those in which an optimal solution is found. Moreover, t max of all the trials was fixed at 500, which was sufficiently large, considering the probability of finding an optimal solution.

Performance with the Reduced Types of Polyominoes
When all the 54 types of polyominoes were available, an optimal solution was found with a high probability and a small number of simulation steps. Therefore, to verify the performance when the number of available polyomino types becomes smaller, we performed experiments by randomly removing the polyominoes from the set of polyominoes of the basic 54 types. Figure 13 shows the probability of finding an optimal solution and the breakdown of the number of placed polyominoes for the number of types of polyominoes. Furthermore, t max is fixed at 500, which was considered large enough (for reference, the results for various t max are shown in Figure S1 in the Supplementary Material). Because the shape of the board became more complicated from Original to Complex, the probability of finding an optimal solution decreased more rapidly, when gradually removing the polyominoes from the set. For the Middle and Complex boards, it was almost impossible to obtain any optimal solution with randomly selected ten types of polyominoes from the original set. This is because a complex local part of a board tends to have few candidate polyominoes to cover it, and all such candidates may have been removed. In contrast, considering the Original board, it was possible to cover the board with fewer types of polyominoes. Even when we randomly selected only ten types of polyominoes, the optimal solution was found with a probability of approximately 40%. 6. Discussion

Optimization of the Coefficients of the Hamiltonian
For the optimization of the coefficients of the Hamiltonian function, we used Optuna [23], a Bayesian optimization library in Python. Bayesian optimization makes it possible to perform the parameter search efficiently, which is extremely time-consuming when grid search is used. Based on the preliminary experiments, the search range of each parameter was set as shown in Table 5, and the number of iterations in the Bayesian optimization was set to 1000 to find the optimal parameters. In each experiment, the obtained parameters were rounded to three digits. Among the searched parameters, a positive correlation was found between coefficient of the term, which restricts the overlap of polyominoes, and coefficient of the term, which suppresses all the pairs of neurons. This is to maintain a balance between the term (which provides a positive bias value for each neuron) and the term (which provides a penalty for each pair of neurons equally). No significant correlation was observed among the other parameters.

Proposal of a New Term That Minimizes the Total Number of Polyominoes Placed
We consider the generalized polyomino puzzle with various sizes of polyominoes in Section 5, and the total number of placed polyominoes is minimized by the effect of the term, which restricts all the pairs of neurons. However, the term is a guiding term 6. Discussion

Optimization of the Coefficients of the Hamiltonian
For the optimization of the coefficients of the Hamiltonian function, we used Optuna [23], a Bayesian optimization library in Python. Bayesian optimization makes it possible to perform the parameter search efficiently, which is extremely time-consuming when grid search is used. Based on the preliminary experiments, the search range of each parameter was set as shown in Table 5, and the number of iterations in the Bayesian optimization was set to 1000 to find the optimal parameters. In each experiment, the obtained parameters were rounded to three digits. Among the searched parameters, a positive correlation was found between coefficient B of the H 2 term, which restricts the overlap of polyominoes, and coefficient F of the H 6 term, which suppresses all the pairs of neurons. This is to maintain a balance between the H 2 term (which provides a positive bias value for each neuron) and the H 6 term (which provides a penalty for each pair of neurons equally). No significant correlation was observed among the other parameters.

Proposal of a New Term That Minimizes the Total Number of Polyominoes Placed
We consider the generalized polyomino puzzle with various sizes of polyominoes in Section 5, and the total number of placed polyominoes is minimized by the effect of the H 6 term, which restricts all the pairs of neurons. However, the H 6 term is a guiding term whose effect is ignored at the end of the simulation by parameter γ. Therefore, considering the current design of the Hamiltonian, the final value of the Hamiltonian is the same, regardless of the size of the polyominoes that cover the board. To explicitly differentiate the final Hamiltonian value depending on the size of the placed polyominoes, we designed an additional term that is inversely proportional to the size of the polyominoes as follows: Adding Equation (14) to Equation (13), the Hamiltonian is modified as follows: The weights W ij between the neurons and biases θ i are obtained as follows: We performed the same experiment as in Section 5.2.2 using the modified Hamiltonian. For the coefficients of the modified Hamiltonian, we used Bayesian optimization to determine the values shown in Table 6.  Figure 14 shows the results of the experiment. Although the probability of finding an optimal solution at t max = 50 increases slightly, no significant improvement is observed overall. Moreover, there is no change in the total number of placed polyominoes, and the introduction of the H 7 term may assist in finding an optimal solution when the problem size increases.
Entropy 2022, 24, 354 17 of 21 whose effect is ignored at the end of the simulation by parameter . Therefore, considering the current design of the Hamiltonian, the final value of the Hamiltonian is the same, regardless of the size of the polyominoes that cover the board. To explicitly differentiate the final Hamiltonian value depending on the size of the placed polyominoes, we designed an additional term that is inversely proportional to the size of the polyominoes as follows: Adding Equation (14) to Equation (13), the Hamiltonian is modified as follows: The weights between the neurons and biases are obtained as follows: We performed the same experiment as in Section 5.2.2 using the modified Hamiltonian. For the coefficients of the modified Hamiltonian, we used Bayesian optimization to determine the values shown in Table 6.  Figure 14 shows the results of the experiment. Although the probability of finding an optimal solution at = 50 increases slightly, no significant improvement is observed overall. Moreover, there is no change in the total number of placed polyominoes, and the introduction of the term may assist in finding an optimal solution when the problem size increases.

Possibility of Application for Drug Discovery
Recently, fragment-based drug discovery (FBDD) has attracted much attention in the field of drug discovery. The method evaluates a drug candidate compound by the combination of fragments of the compound [15,24,25]. The number of unique fragments included is rather saturating, even when the size of the compound database grows very rapidly [24]. Considering this context, several methods have been proposed for solving the problems appearing in the FBDD as combinatorial optimization problems. For example, eHiTS performs a fragment-based docking simulation by formulating the problem as a weighted maximum clique-finding problem [25]. In addition, Banchi et al. solved a similar weighted maximum clique problem by emulating the Gaussian boson sampler (GBS), which is a type of optical quantum computer [15].
The problem of filling the pocket of a target protein with fragments is similar to the generalized polyomino puzzle of covering a board with polyominoes. Therefore, we consider that the proposed Hamiltonian can be applied to drug discovery in the future. The fragments have varying sizes, and each fragment can be used zero or more times. Thus, our extension of the generalized polyomino puzzle illustrated in Section 5 is similar to the real-world problem. We needed a better resolution for representing the spherical shape of a protein pocket; therefore, we adopted the same Hamiltonian function with polyhex tiles ( Figure 15).

Possibility of Application for Drug Discovery
Recently, fragment-based drug discovery (FBDD) has attracted much attention in the field of drug discovery. The method evaluates a drug candidate compound by the combination of fragments of the compound [15,24,25]. The number of unique fragments included is rather saturating, even when the size of the compound database grows very rapidly [24]. Considering this context, several methods have been proposed for solving the problems appearing in the FBDD as combinatorial optimization problems. For example, eHiTS performs a fragment-based docking simulation by formulating the problem as a weighted maximum clique-finding problem [25]. In addition, Banchi et al. solved a similar weighted maximum clique problem by emulating the Gaussian boson sampler (GBS), which is a type of optical quantum computer [15].
The problem of filling the pocket of a target protein with fragments is similar to the generalized polyomino puzzle of covering a board with polyominoes. Therefore, we consider that the proposed Hamiltonian can be applied to drug discovery in the future. The fragments have varying sizes, and each fragment can be used zero or more times. Thus, our extension of the generalized polyomino puzzle illustrated in Section 5 is similar to the real-world problem. We needed a better resolution for representing the spherical shape of a protein pocket; therefore, we adopted the same Hamiltonian function with polyhex tiles ( Figure 15). However, there are still huge differences between them. For example, regarding the FBDD, the space to be filled is 3D instead of 2D, and the pieces filling the space will take more complex shapes. Figure 16 shows examples of the Hamiltonian proposed in Section 4.2, applied to a problem with more complex pieces in 3D by using a cube instead of a square. We performed a series of computer experiments on these problems; however, the details are not presented in this paper. However, there are still huge differences between them. For example, regarding the FBDD, the space to be filled is 3D instead of 2D, and the pieces filling the space will take more complex shapes. Figure 16 shows examples of the Hamiltonian proposed in Section 4.2, applied to a problem with more complex pieces in 3D by using a cube instead of a square. We performed a series of computer experiments on these problems; however, the details are not presented in this paper.
One further difference is the necessity of satisfying constraints. For drug discovery, filling the 3D space, called the binding pocket, in the drug target protein is not the main aim, but rather designing more stable molecules in the binding pocket. This means that having some space remaining is acceptable. Furthermore, as the binding pocket has some flexibility, there is tolerance for overlapping pieces. This means that states that do not strictly satisfy the constraints may still be valuable. Considering this characteristic, some space on the board and overlapping of polycubes are allowed. One further difference is the necessity of satisfying constraints. For drug discovery, filling the 3D space, called the binding pocket, in the drug target protein is not the main aim, but rather designing more stable molecules in the binding pocket. This means that having some space remaining is acceptable. Furthermore, as the binding pocket has some flexibility, there is tolerance for overlapping pieces. This means that states that do not strictly satisfy the constraints may still be valuable. Considering this characteristic, some space on the board and overlapping of polycubes are allowed.
More importantly, we need to optimize the space-filling and the physicochemical interactions such as electrostatic interaction, hydrogen bonding, and hydrophobic interaction between the atoms. By defining a reward in each cube of the board (for each atom type), it is also possible to design the space-filling problem to maximize the sum of the rewards defined in each cube with soft constraints. When applying the Hamiltonian proposed in this study to this problem, we can introduce a reward for each cube in the term as follows: However, there are more difficult technical issues that need to be solved. First, the variety of possible fragments corresponding to the types of polycubes is enormous. For example, Yanagisawa et al. [24] demonstrated that the total number of types of fragments reached 263,319 for compounds registered in the ZINC12 database [26], a collection of commercially available chemicals. Therefore, instead of dealing with all these chemically observable fragments, it is necessary to select a set of core and important fragments that can be used to navigate the early stages of drug design. It is also necessary to consider the fine rotation of the fragments in the 3D space, and it may require a huge number of neurons or spins. Therefore, it is important to narrow down the number of candidates in advance, depending on the local information of the protein pocket space. Furthermore, although the polyomino puzzles only consider the contact between the polyominoes, the fragments on molecules have chemical bonds. In principle, these issues can be addressed by modifying the Ising model treated in this study; however, the number of spins will be too large unless it can be reduced using domain knowledge.

Conclusions
In this study, we considered a standard polyomino puzzle in which the board is covered with a given number of polyominoes, designed an improved Hamiltonian function, More importantly, we need to optimize the space-filling and the physicochemical interactions such as electrostatic interaction, hydrogen bonding, and hydrophobic interaction between the atoms. By defining a reward in each cube of the board (for each atom type), it is also possible to design the space-filling problem to maximize the sum of the rewards defined in each cube with soft constraints. When applying the Hamiltonian proposed in this study to this problem, we can introduce a reward R n for each cube n in the H 2 term as follows: However, there are more difficult technical issues that need to be solved. First, the variety of possible fragments corresponding to the types of polycubes is enormous. For example, Yanagisawa et al. [24] demonstrated that the total number of types of fragments reached 263,319 for compounds registered in the ZINC12 database [26], a collection of commercially available chemicals. Therefore, instead of dealing with all these chemically observable fragments, it is necessary to select a set of core and important fragments that can be used to navigate the early stages of drug design. It is also necessary to consider the fine rotation of the fragments in the 3D space, and it may require a huge number of neurons or spins. Therefore, it is important to narrow down the number of candidates in advance, depending on the local information of the protein pocket space. Furthermore, although the polyomino puzzles only consider the contact between the polyominoes, the fragments on molecules have chemical bonds. In principle, these issues can be addressed by modifying the Ising model treated in this study; however, the number of spins will be too large unless it can be reduced using domain knowledge.

Conclusions
In this study, we considered a standard polyomino puzzle in which the board is covered with a given number of polyominoes, designed an improved Hamiltonian function, and proposed a method to obtain solutions efficiently using the Ising model. Moreover, we considered the generalized polyomino puzzle with relaxed restrictions on the number of times the polyominoes could be used, the various polyomino sizes, and the flexibility of the board shape. To verify the performance of the proposed method in finding the optimal solution, we conducted evaluation experiments on a software simulator of the Hopfield neural network. Using our formulation, we were able to solve the standard pentomino puzzle-with the problem represented by 1928 spins-in more than 90% of the cases, whereas the models proposed in previous studies had difficulty solving them. Furthermore, several generalized polyomino puzzles were solved approximately 100% of the time with a small number of simulation steps. In addition, we consider the possibility of future applications to real problems such as fragment-based drug discovery. When dealing with combinatorial optimization problems inherent in real-world problems, the scale of the model to be solved can be enormous if we apply a simple modeling approach that treats all candidates equally. Therefore, it is necessary to keep the size of the problem compact by considering the specific domain knowledge and local properties of the problem. A theoretical study aimed at improving the design method of the Ising model and efficiently controlling the dynamics of the system to handle larger-scale optimization problems should also be considered.