An Efficient Enumeration of Flat-Foldings : Study on Random Single Vertex Origami
Chihiro Nakajima


Gr-IoU: Ground-Intersection over Union for Robust Multi-Object Tracking with 3D Geometric Constraints
    Keisuke Toida10009-0006-4873-3651 Naoki Kato20009-0004-3815-0829 Osamu Segawa20009-0000-2469-6098 Takeshi Nakamura20009-0001-4991-3383 Kazuhiro Hotta10000-0002-5675-8713

    September 9, 2024
==============================================================================================================================================================================



§ ABSTRACT 
This paper deals with themes such as approximate counting/evaluation of the total number of flat-foldings for random origami diagrams, evaluation of the values averaged over various instances, obtaining forcing sets for general origami diagrams, and evaluation of average computational complexity. An approach to the above problems using a physical model and an efficient size reduction method for them is proposed.
Using a statistical mechanics model and a numerical method of approximate enumeration based on it, we give the result of approximate enumeration of the total number of flat-foldings of single-vertex origami diagram with random width of angles gathering anound the central vertex, and obtain its size dependence for an asymptotic prediction towards the limit of infinite size.

In addition, an outlook with respect to the chained determination of local stacking orders of facets caused by the constraint that prohibits the penetration of them is also provided from the viewpoint of organizing the terms included in the physical model.
A method to efficiently solve the problem of the determination or enumeration of flat-foldings is discussed based on the above perspectives. 
This is thought to be closely related to forcing sets.





§ INTRODUCTION





Statistical mechanics deals with configurations given by combinations of variables that take on a small number of states (for example, binary variables)<cit.>.
It introduces realization probabilities to configurations and discusses the properties of moments and large deviation functions, and is closely related to mathematics through probability theory and combinatorics.
In mathematical and computer science research on origami, a statistical mechanics perspective has been introduced to some problems (though it is sporadically, not based on a unified perspective)<cit.>.




§ PHYSICAL MODEL FOR FOLDING OF ORIGAMI DIAGRAM




In this study, the folding diagrams are generated according to probabilistically generating algorithm explained below.
With the total number of facets which is descried as n is fixed, large numbers of instances are generated by randomly giving widths of angles around the center, and discuss the statistical properties of the resulting set of instances.
We try to obtain asymptotic predictions on the behavior in the limit of an infinite number of facets n →∞ from the sequences of those with finite number of facets. 

[][ht]
  
    < g r a p h i c s >

  (a)Example of origami diagram. Each edge in the figure represents a crease. In this figure there are no overlaps of facets. (b)Corresponding pre-folded diagram, which describes the overlaps of facets when the figure (a) is folded along the creases. Each vertiex indicated by the same mark is the same as that in Fig. (a). (c)Schematic picture of introduction of the Ising variable to a local layer-ordering.
  


Attempts to introduce physical models into origami have been made in the past.
It includes one that deals with geometric constraints on the bending of flat structures in three-dimensional space<cit.>,
 one aimed at the physical properties of polymers and membranes<cit.>,
 one that extracts the phase transition phenomenon seen in self-folding origami<cit.>,
 and one that uses it in the context of mathematical research<cit.>.

In this paper we consider the model consisting of binary variables s_i,j∈{-1,+1} with the condition subjected with the form of products of two or four variables<cit.> to traet the stacking of facets and determinate flat foldability accurately.
The variables s_i,j represent the vertical relationship in stacking of two facets i and j in the origami diagram.
Each realization of s_i,j∈{-1,+1} represent the global stacking order of n pieces of facets.

The conditions to prohibit the situation that a crease is penetrated or intruded by other facets are imposed with the following terms in the energy function,

    E^(i)_i,j;k
         =1/2(1-J_(ik)(kj)s_i,ks_k,j)                
    
      E^(c)_i,j,k
         =1/4(1-L_(ij)(jk)s_i,js_j,k-L_(jk)(ki)s_j,ks_k,i-L_(ki)(ij)s_k,is_i,j), 
    
      E^(q)_i,j ; k,l
         =1/2(1-K_ijkls_i,ks_i,ls_j,ks_j,l),

where J_(ik)(kj)=-τ_ikτ_kj,
K_ijkl=τ_ikτ_ilτ_jkτ_jl,
L_(ij)(jk)=-τ_ijτ_jk,
and τ_ij is the sign of the difference between two indices of facets, namely τ_ij=sign(j-i).
Note that for the variable s_i,j, pay attention to the order in which subscripts are written, is assigned only when i<j and vice versa.
The term (<ref>), which prohibits a facet k from intruding between two facets i and j which are connected by a crease, is assigned for a geometry where a crease with (i,j) has an overlap with a facet k in the pre-folded diagram.
The term (<ref>), which prohibits cyclic stacking among the three facets i, j, k, is assigned for a geometry where three facets i,j and k simultaenously share an area.
The term (<ref>), which prohibits the unrealizable overlap between the four facets i, j, k, l that form two creases in (sometimes partially) coincident position, is assigned for a geometry where the two creases each consisting two facets i and j, k and l with respect are in coincident position.
The terms with the form (<ref>)-(<ref>) are respectively assigned to each geometries of overlap of a crease and a facet, simultaneously shared area, and coincident creases in the pre-folded diagram.
Hence, the total energy function is described as followings,

    H({s})=∑_(i,j ; k)E^(i)_i,j;k + ∑_(i,j,k)E^(c)_i,j,k + ∑_(i,j;k,l)E^(q)_i,j;k,l,

where the each summation is taken over all corresponding geometries in the pre-folded diagram.
Thus the problem of finding the flat folding is casted as an optimization problem,
that of finding the conbination of {s_i,j} with H=0.



§ SINGLE-VERTEX ORIGAMI DIAGRAM


Preparing a flat-foldable diagram in general is considered to be equivalent to a combinatorial optimization problem itself.
In fact, when an origami diagrams are generated at random angles, the frequency of the flat-foldable diagrams is extremely small, so that it is almost unrealistic to discuss the foldability problem.

Instead of the generation of general diagrams, generating ones of single-vertex structure such as exhibited in Fig.<ref>(a) is performed in this research.
The properties of single-vertex diagrams have been actively researched in the context of determining foldability under a fixed crease pattern (mauntain-valley assignment), enumerating foldable crease patterns, and enumerating flat-foldings (global stacking-order of facets)<cit.>.
It is known that the foldability of a single vertex diagram can be determined using the Kawasaki's theorem<cit.>.
This themrem makes to generate foldable diagrams with random angles around the center possible by maintaining the condition that the alternating sum of those is 0.
Regarding the enumeration of the total number of flat foldings at least for the case of random angles, even in single vertex diagrams there is still no known conclusion, as far as the narrow knowledge of the auther.



 §.§ Numerical procedure for generating each instance




For an origami diagram whose total number of facets is n, the total number of pairs of facets is n(n-1)/2, so the upper limit of N, the number of pairs which the local layering order should be considered, is also this value. 
However, depending on the instance, there are facets that do not have an overlapping area in the folded diagram. 
No variables are assigned to these pairs, as there is no need to consider direct hierarchical relationships. 
As a result, the value of N for the origami diagram of n facets roughly takes a value close to the upper limit, but takes various values depending on the details of the overlap.
Thus the number of variables s_i,j for the corresponding optimization problem is given by N.

The instances are generated with two method, the examples of diagrams are shown in Fig.<ref>.
In the first method, the width of each angle around the center is ramdomly given from continuous real numbers in the range 0 to 1 by uniform distribution.
The angles are given so that the alternating sum is 0, and eventually normalized so that the sum is 2π.
This can be expected to correspond to the infinitesimal limit of the width compared to the case where discrete unit width w is introduced to the angles, which will be described later.
[][ht]

    < g r a p h i c s >

Examples of of single vertex origami diagrams,
(a) An origami diagram with n=24 and randomly generated width of angles around a center vertex,
(b) A Corresponding pre-folded diagram for the diagram (a),
(c) An origami diagram with n=24 and w=24.
(d) A Corresponding pre-folded diagram for the diagram (b),







 §.§ Definition of quantities




  §.§.§ Entropy as Logarithm of Total Number of Flat-Foldings

In statistical mechanics, a quantity called a partition function Z(β) works as a generating function to obtain the expected value of various physical quantities, energy, each variable, or the sum of them at a fixed temperature.
Its definition is the summation of the quantity exp(-β H) over combinations of variables included in the energy function H, namely,
Z(β)=∑_s_i,j=± 1∑_s_i',j'=± 1⋯exp(-β H({s})),
where the summation symbol over many variables, ∑_s_i,j=± 1∑_s_i',j'=± 1⋯, means that it is taken over all combinations of the values of the variables of the system.
The H discussed in this paper obviously has a minimum value 0, at least for the origami diagrams that satisfy the Kawasaki theorem. Therefore, the value of Z(β) in the limit of β→∞, which is called zero-temperature limit in physics, is the sum of the value 1 by the total number of flat-foldings. 
The logarithm of the partition function, F(β)=βlog Z(β), is usually called the free energy.
In particular for the energy function with the minimum value 0, the zero temperature limit lim_β→∞ F(β) is also called zero-temperature entropy or ground-state entropy.
The value of the ground-state entropy is numerically obtained by integrating the expected value of energy at a finite value of β using the following formula.
This paper describes the results of an approximate evaluation of this absolute zero entropy using a numerical calculation or probabilistic sampling method based on the Markov chain Monte Carlo method, especially the replica exchange Monte Carlo method<cit.>.

In this paper, in order to discuss the average behavior for various folding diagrams, the zero-temperature entropy lim_β→∞ F(β) on each instance is further averaged over several instances.
Therefore we represent each instance of origami diagram with the symbol Δ and write the zero temperature entropy for a diagram Δ as S_tot^(Δ). 
In addition, these quantities averaged over all generated Δ is represented as [S_tot]=(∑_Δ S_tot^(Δ)) / (∑_Δ1).



  §.§.§ Number of facets sandwitched


If a certain facet k is sandwiched between two other facets i and j, the product s_i,ks_j,k takes the value τ_ikτ_kj=-J_(ik)(kj).
Hence, the number of facets which are sandwitched between i and j is obtained by calculating the following quantity,

    n^(sdw)_i,j=∑_k1-J_(ik)(kj)s_i,ks_j,k/2,

where the running suffix k is taken for every k for which a variable s is given between both of the facets i and j those which form a crease.
However, no variables are introduced for facet pairs that do not have a direct vertical relationship due to their overlapping positions. 
In other words, a facet l that does not have a vertical relationship with either of the two facets i or j is not included in the sum ∑_k in Eq. (<ref>).
Therefore, in the scope of this study, facets that are sandwiched indirectly or hidden are not counted in the number.

After computing n^(sdw)_i,j for all creases (i,j) included in the diagram,
the maximum and minimum of n^(sdw)_i,j included in each flat folding are written as n_max, n_min.
In addition, for each value of n_max and n_min, we approximately enumerated (the logarithm of) the number of flat foldings that have such values as the maximum and minimum of n^(sdw)_i,js, S(n_max) and S(n_min) resepctively, in the same way as S_tot.
Furthermore, these values are averaged over instances as well as S_tot. Those are written as [S(n_max)] and [S(n_min)].
In the below section <ref>, they are written [S(ν_max)] and [S(ν_min)] with ν_max=n_max/n and ν_min=n_min/n.

The number of facets sandwiched between a crease is also called the crease width. 
Awareness of this issue was introduced in the paper <cit.> and a research context called folding complexity has been developed.
It is known that the problem of finding a flat-folding that minimizes the maximum value of crease widths n_max in stamp-folding problem, in the single vertex origami-diagram with uniform facets and open boundary, is NP-hard.

In introducing [ S(n_max) ] or [ S(n_min) ], it must be mentioned that the maximum and minimum values of n_max or n_min are different for each instance.
There is a gap between the average value taken over only instances where non-zero contribution of S(n_max) is confirmed and the frequency itself of obtaining instances with foldings whose maximum crease width is n_max.
If an instance does not have a fold whose maximum crease width is n_max, its number is 0 and its logarithm is -∞, which is not suitable for averaging.
To deal with this problem with, the following two quantity is computed,

    [S(n_max)]=lim_α→ 0[ log{α+exp(S(n_max))}],
    [S(n_max)]=∑_Δ∈ D^(n_max)S^(Δ)(n_max)  /  ∑_Δ∈ D^(n_max)1,

where 
D^(n_max) is the set of instances that have foldings whose maximum crease width is n_max.



  §.§.§ Overlap between Two Configurations

As an indicator of how the solution space corresponding to flat folding is embedded in the variable space expressed by the spin variables, we conducted two independent numerical experiments at sufficiently low temperatures and computed the overlap of the spin variables. 
The configurations sampled from two independent numerical simulations at the same temperature are denoted by spin variables s_i,j and s'_i,j, respectively. 
We calculate the following quantities for these pair of simulations,

    Q=∑_i,j1-s_i,js'_i,j/2,

where the summation in Eq. (<ref>) is taken over all variables of the system composed of {s_i,j}.
Let the total number of variables be N and q=Q/N. If q=1, the two configurations are completely the same, and if 0, they are completely reversed. 
Also, when q ≃ 0.5, it means that about half of all variables are reversed, and the two configurations can be said to be almost uncorrelated or “unrelated”.

A normalized histogram h_β,n(Q) is obtained from the sequence of Qs sampled by a pair of numerical simulations with a certain value of (inverse) temperature β.
Using h_β,n(Q), the probability density distribution P_β(q) is derived as,

    P_β(q)=∑_Q=0^N h_β,n(Q) χ_q(Q/N)/∑_Q=0^N h_β,n(Q),

where χ_q(Q/N) is an indicator function that returens the value 1 when Q/N=q and returns 0 when Q/N ≠ q.
This quantity provides informations about how far apart the sampled configurations in a pair of simulations are from each other in variable space at a given temperature.
At the limit of β→∞, this asymptotes to the structure of the embedding of the solution in the configuration space.
If there is not one or only a small number of configurations that completely satisfy the constraints, then P_β(q) is sufficiently low and q=1, and the physical model in question in this paper is considering that it is symmetric about inversion, there is a sharp peak near q=0 and q=1. 
Furthermore, if there is another satisfied configurations in a region where the values of some variables are different, a peak similarly occurs in the region of q values corresponding to the difference.
The average over instances are taken as well as the entropies, namely [P_β(q)]=∑_Δ{P_β^(Δ)(q)× 1/(∑_Δ1)}.



§ RESULTS






 §.§ Average behavior of total numbers of foldings

The logarithm of the total number of flat-foldings in a origami diagram consisting of facets of uniform length is proportional to nlog n. 
Whereas, in diagrams with random facets generated by the procedure explained in Section <ref>, the entropy value which corresponds to the total number of folds exhibits the behavior shown in Fig..
In Fig. the values of entropy is proportional to n in average, which means that the dependence of the total number of flat-foldings on n is exponential.
This dependence is qualitatively different from the case with uniform facets.
In the exponential function, whose form is exp(γ nlog n), the value of γ can be readed as 0.098 from the figure, even though the figure is the result of a numerical experiment and is an averaged value over various instances.

[][ht]
  
    < g r a p h i c s >

  
The values of total entropy  averaged over instances [S_tot] is plotted against the inverse of the number of facets, 1/n, for n=18 (square), 24 (circle) and 36 (triangle).
Each point is overlaid with an error bar of the average value itself over instances and the variance over instances.
The error bars are given by the bootstrap method and drawn as solid lines, however they are smaller than the size of the mark. 
The variances are written with dashed lines.

  





  §.§.§ Distribution of the minimal and maximal numbers of sandwitched facets


Next, the results regarding the entropy of flat-foldings on the maximum and minimum numbers of facets sandwiched between each crease, respectively noted as n_max and n_min, are described.

In Figure <ref> the horizontal axis describes n_max or n_min divided by the total number of facets of the system n,
which are respectively represented as ν_max=n_max/n and ν_min=n_min/n,
 and the vertical axis is the logarithm of the total number of the flat foldings such that the value of n_max or n_min takes the value on the horizontal axis, [ S(ν_max) ] and [ S(ν_max) ] with respect devided by n log n.
A curve is shown for n=18,24 and 36.
For reference, the same quantities for a single vertex diagram with uniform facet angle are shown in the inset.


The same as [S_tot], the curves of [S(ν_max)] or [S(ν_min)] for various sizes appears to collapse into a master curve when divided by n log n.
The master curve is found in both cases with uniform and random facets. 
However, those two exhibit qualitatively different behaviors.
In the case of uniform facets, the curve of S(n_max) increases monotonically as the value of n_max increases. 
This suggests that various combinations of stacking orders are realized even among the facets sandwiched between the creases with the maximum width.
On the other hand, in the case of random facets, the value of S(ν_max) shows a behavior that decreases slightly in the region where the value of n_max is large.
At n_max=n-2, it roughly takes the maximum value of S(ν_max) ≃ 0.3 × nlog n. 
Whereas, in the case of random angle width [S(ν_max)] takes its maximum value, approximately [S(ν_max)]≃ 0.1 × nlog n, at ν_max≃ 0.3 and decreases once the value of ν_max exceeds that argument as shown in Fig.<ref>.
Also regarding [S(ν_max)], this scaling behavior, including the value of the coefficient, of the maximum value is the same as shown in Fig.<ref>, although the decrease of the value is rather significant in this case bacause the fraction becomes smaller as n_max goes to its large value.

Meanwhile, S(ν_min) or [S(ν_min)] have its contribution only on ν_min=0, both for origami diagrams with random and uniform angle width.
This is thought to be due to the fact that in the single-vertex diagram the array of facets is closed, that is, all facets have creases at both side of them and no have an open end.
In the case with exact stamp folding problem, there obviously exists a folding in which the minimum value of the crease width is a finite positive value.

Note here that there is nothing unnatural about the curve of [S(ν_max)] having a finite number of folds assigned to the position ν_max=0.
For example, if a random folding diagram contains two facets with fairly large angles and they are not connected by a crease, it is possible to create a folding where one of the two is located at the top in the vertical relationship of the overall stacking and the other is at the bottom.
In such a folding, between the top and bottom facets, two stacks may be formed at a considerable distance from each other.
Due to the circumstances described in Section <ref>, in such stacking there may be the case in any creases and two facets that compose them no other facets sandwiched between.



  §.§.§ Averaged Shape of Overlap Distribution


The results of this study are shown in Figure <ref>.

The contribution of the distribution is widely dispersed, with the ratio of the degree of overlap to the total number of variables q ranging from about 0.2 to 0.8. 
Comparing to the shape of the overlap distribution of the origami diagram with a uniform angular width is a Gaussian,
 that on diagrams with random angular widths clearly has the different shape. 
Instead, the contribution of [P(q)] is large even in regions where the absolute value of q is far from 0. 
Therefore, it can be seen that the set of flat foldings is composed of stacks of facets whose vertical relationships are very different from each other.

As a consequence of the theory of replica symmetry breaking in statistical mechanics, the three peak structure in the asymptotic shape of [P(q)] (at q=0.5 and symmetrical positions centered on it) implies that the average complexity class of the seach problem is NP-complete<cit.>.
When the size dependence shown in Fig. <ref> is focused based on this viewpoint, the shape found in 0.2<q<0.8 in the range up to n = 36 somewhat seems to lead the three peak structure.
Meanwhile,the width of the error bar is also large. 
Thus, it yet require the carefull observation to conclude that the average-case computational complexity of searching flat foldings of random single-vertex origami diagrams is NP-complete.



 §.§ Reduction of System Size for Computational Efficiency

In the physical model given by Eq. (<ref>), each term with the form (<ref>) in particular gives a relationship that should be satisfied in the ground state between variables respectively involved in each term.
Applying this relationship to a set of variables in (P) and rearranging them brings an efficiency for numerical approximations to the total number of flat folds and overlapping distributions we have seen in the above by mapping the ground state to a system described by fewer variables.
In addition, this procedure is expected to aggregate variables whose values are determined uniquely each other into a single cluster and to reveal the combinatorial structure that is the essential cause of the problem.

Here, the result of the reduction of variables by the contraction via the equations (<ref>) for the origami diagram instances ranging from n=24 to 68.
For the origami diagram obtained using the method in Sec.<ref>, the variables are reduced as demonstrated in Sec. <ref>.
This reduced number of variables is represented as Ñ and the ratio R=Ñ/N is discusssed for each instance.
Note that there are many cases that the reduction result as the collection of several independent component each consisting of variables {s,C}. 
The two variables s_i,j or C_l which belong to each different component are not involved in a same product term.
For such cases let Ñ be the number of variables contained in the largest component in the collection.
The density plot of the distribution of the ratio D(R) for s = 24,48,68, from 465, 518, 319 instances respectively, is shown in Fig.<ref>. 


In Fig.<ref>, for small n, for example n=24, the distribution has a unimodal peak roughly around 0.05<R<0.15.
However, in the case of a larger system size, another peak appears on the side with a smaller value of R.
As the system size further increases, the distribution becomes bimodal, with the height of two peaks swapping, and the peak with the smaller R value becoming dominant.

A similar change in the shape of D(R) also occurs in the case with sets of instances generated by the second method, the method with the minimum unit is introduced to the angular width.
As shown in Fig.<ref>, it is observed when the value of s is fixed and the values of w are increased.
From this observation, it is thought that when the number of terms with a product of four-variables like Eq.(<ref>) is large, D(R) has a unimodal shape with its peak in the region of relatively large values of R,
and as the number decreases the shape changes to one with a dominant peak on the small-R side via the bimodal shape.

However, if the angular width takes continuous random values for a small n, the number of terms like Eq.(<ref>) in the corresponding cost function still remains almost 0, which does not mean that there are a large number of terms.
For this point further consideration on the relationship to the number of the four-variables product is currently required.

Within the range of system sizes tested, the system size Ñ can be reduced to roughly 1/20 of the original system size.
Compared to this, in the diagram with a uniform angles width, we have to treat the system with n(n-1)/2 variables because there is no room for reducing the variables using the contraction procedure.
This suggests that the average computational amount for a set of random instances can behave quite differently compared to that for a uniform angular value, which is considered to be the corresponding to the worst-case computational complexity.
It is also interesting from a computational complexity perspective.



  §.§.§ Trying to obtain system-size dependence of each separated peak

In Fig. <ref>, it is found that the peak on the small R side in D(R) gradually becomes more prominent as the larger n. 
Although being aware that this may not be an accurate evaluation, an attempt to evaluate the asymptotics of this peak is made as shown in Fig.7. 
The dispersion σ_R of the small-R side peak, meaning the width of the peak, is computed using only instances whose R are less than or equal to a cutoff.


    < g r a p h i c s >


The dispersion σ_R for the lower side peak in D(R) with various thresholding values.



While the cutoff is located between the two peaks, the slope of the variance change is relatively small, and when the cutoff is located on the small R peak, the slope of the change is large. 
Using this behavior, we estimated the value at which the slope of the curve representing dependence begins to increase (again) as seen in the caption of Fig.7. 
The threshold value of R is estimated as Rt = 0.030, 0.037, 0.045 for n = 36, 48, 68, respectively. 
These values exhibit that it is also not possible to exclude the possibility that [Ñ] has a depencence of sub-linear order on [N], where [Ñ] and [N] is averaged over instances because the system size N is different among each instances.



§ SUMMARY


In this study based on the physical model formulation, approximate enumeration of the total number of foldings and its decomposition into those with fixed maximum and minimum values of the number of facets between each crease, i.e., the crease width.
The value of the logarithm of the total numbers of the configuration, computed by the replica exchange Monte Carlo method, is 0.098 × nlog n and rather smaller than that regarding the origami diagrams with uniform angle width.
Furthermore, the behaviors of the decomposed entropy with the maximum number of facets, n_max sandwithed by each crease of origami diagrams of respectively random or uniform angle width behave qualitatively different from  each other.
This comparison implys that the constraints on mutual penetration of facets sandwiched within the crease with the maximum width brings a limit of the diversity of combinations of their stacking orders.

In addition, we approach the average computational complexity of the problem of enumerating the number of flat folding of the single vertex origami diagram from two perspectives: observation of the behavior of overlap distribution functions from the viewpoint of replica symmetry breaking in spin glass theory, and reduction of system size based on chain relationships that fix the values of variables.
Given the range of system sizes observed in this study, it has still not been possible to draw conclusions regarding the results from both and their consistency with each other. However, it is expected that research in this direction will still have the potential to be a future topic.





§ APPENDIX







 §.§ Contraction of spin variables


Here, the method of variable reduction is intruduced based on the example of the origami diagram shown in Fig.<ref>.
First, from the diagram of this figure, the following cost function is obtained by the modeling described in Section <ref>.
In a term whose form is the same as Eq.(<ref>), the combination of the involved spin variables that brings the value 0 of the term is uniquely determined except for total inversion of the both two.
Therefore, we can translate the constraint from the term like Eq.(<ref>) into an allowed relationship between the spin variables.
By using such translation, several spin variables that appear in the energy function can be collectively re-expressed as a single variable. 
And the energy function can be re-given as a combination of fewer binary variables.


    < g r a p h i c s >


An example of diagram of single vertex origami with n=6. α, β, γ, δ at the bottom of the figure are the angle values around the center point of facets.
 two colored creases in (a) overlap when the diagram is flat-folded as in (b).



    ∑_(ijkl)E^(q)_ijkl   =1+s_24s_25s_34s_35/2,
    ∑_(ij,k)E^(i)_ij,k   =1-s_13s_23/2+1-s_14s_24/2+1-s_12s_13/2+1-s_26s_36/2
           +1-s_14s_15/2+1-s_46s_56/2+1-s_15s_16/2+1-s_25s_26/2
           +1-s_35s_36/2+1-s_45s_46/2,
    ∑_(ijk)E^(c)_ijk   =1+s_12s_23-s_23s_13-s_13s_12/4+⋯+1+s_35s_56-s_56s_36-s_36s_35/4,

where the detail of Eq.(<ref>) is written as Eq.(<ref>)-(<ref>), described later.
By applying the rewriting process based on the product of two spin variables, which is mentioned at the beginning of this section, from the 2nd to 11th term of the above cost function (<ref>) each leads the relationship between two variables.
The relationships are eventually summarized as follows,

    C_1=s_12=s_13=s_23,
    
    C_2=s_14=s_24=s_15=s_16,
    
    C_3=s_26=s_36=s_25=s_35,
    
    C_4=s_45=s_46=s_56.

With the above introduced cluster variables, the terms of Eq.(<ref>) is rewritten as

    1+s_12s_25-s_25s_15-s_15s_12/4   =1+C_1C_3-C_3C_2-C_2C_1/4,
    1+s_12s_26-s_26s_16-s_16s_12/4   =1+C_1C_3-C_3C_2-C_2C_1/4,
    1+s_34s_45-s_45s_35-s_35s_34/4   =1+s_34C_4-C_4C_3-C_3s_34/4,
    1+s_34s_46-s_46s_36-s_36s_34/4   =1+s_34C_4-C_4C_3-C_3s_34/4,
    1+s_12s_23-s_23s_13-s_13s_12/4   =1+C_1C_1-C_1C_1-C_1C_1/4=0
    
    
    1+s_12s_24-s_24s_14-s_14s_12/4   =1+C_1C_2-C_2C_2-C_2C_1/4=0 
    
    
    1+s_13s_34-s_34s_14-s_14s_13/4   =1+C_1s_34-s_34C_2-C_2C_1/4,
    
    
    1+s_13s_35-s_35s_15-s_15s_13/4   =1+C_1C_3-C_3C_2-C_2C_1/4
    1+s_13s_36-s_36s_16-s_16s_13/4   =1+C_1C_3-C_3C_2-C_2C_1/4,
    
    
    1+s_14s_45-s_45s_15-s_15s_14/4   =1+C_2C_4-C_4C_2-C_2C_2/4=0,


    1+s_14s_46-s_46s_16-s_16s_14/4   =1+C_2C_4-C_4C_2-C_2C_2/4=0,
    
    
    1+s_23s_35-s_35s_25-s_25s_23/4   =1+C_1C_3-C_3C_2-C_2C_1/4, 
    1+s_23s_36-s_36s_26-s_26s_23/4   =1+C_1C_3-C_3C_2-C_2C_1/4,
    1+s_24s_45-s_45s_25-s_25s_24/4   =1+C_2C_4-C_4C_3-C_3C_2/4,
    
    
    1+s_24s_46-s_46s_26-s_26s_24/4   =1+C_2C_4-C_4C_3-C_3C_2/4,
    1+s_15s_56-s_56s_16-s_16s_15/4   =1+C_2C_4-C_4C_2-C_2C_2/4=0, 
    
    
    1+s_45s_56-s_56s_46-s_46s_45/4   =1+C_4C_4-C_4C_4-C_4C_4/4=0,
    1+s_25s_56-s_56s_26-s_26s_25/4   =1+C_3C_4-C_4C_3-C_3C_3/4=0,
    1+s_35s_56-s_56s_36-s_36s_35/4   =1+C_3C_4-C_4C_3-C_3C_3/4=0.


In particular, we focus on the rewriting of the first term, Eq.(<ref>), which is

    1+s_24s_25s_24s_35/2=1+C_2C_3s_34C_3/2=1+C_2s_34/2,

whose form composes the product of two spin variables similar to Eq.(<ref>) again.
The reproduction of the term whose form is similar to Eq.(<ref>) induces the recursive application of the rewriting process.
Here, the following relationship,

    s_34=-C_2,

is lead from Eq.(<ref>).
In addition, further rewritings and relationships are induced from Eq.(<ref>) as follows,

    1-C_1C_2+C_2C_2-C_2C_1/4   =1-C_1C_2/2,
    
    1-C_2C_4-C_4C_3+C_3C_2/4+
    1+C_2C_4-C_4C_3-C_3C_2/4   =1-C_4C_3/2,

which results

    C_1=C_2,
    
    C_3=C_4,

where Eqs(<ref>) and (<ref>) are brought from Eqs.(<ref>) and (<ref>)+(<ref>)+(<ref>)+(<ref>), respectively.

The method of variable reduction explained here utilizes the relationship between variables, meaning the vertical relationship of facets, that must be satisfied in the ground state (state of zero energy) of the physical model. 
The concept that gives a cascading chain of decisions on the vertical relationship of facets is also studied in mathematics and information science, and is called a forcing set.
Algorithms for finding forcing sets for one-dimensional origami<cit.> and Miura-ori diagrams<cit.> have been proposed, however it appears that no algorithms exist for general origami at present.
Based on a series of studies on this subject, the original definition of a forcing set is considered to be the subset of minimal assignments of the crease patterns that determine the assignments for all creases included in the origami diagram. 
It is known that the problem of giving an entire forcing set is can be an NP-complete problem.
However, a chain of decisions for stacking relationships can be given for local relationships among facets. 
Therefore, it is an interesting application to make the determination of the possibility or enumeration of flat-folding more efficient by partially utilizing the subsets of chained decisions contained in forcing sets.
The method demonstrated in this section is expected to provide insight into the derivation of forced sets for general origami.




 §.§ Numerical method for approximately estimating the number of folding


The numerical simulations performed for this study are based on the replica-exchange Monte Carlo method and the multiple histogram reweighting technique.

Each variable s_i,j is updated using the Metropolis rule. That is, the accepting probability for local updating,  p_ud, is given by

    p_ud=min{1,exp(-β (E'-E))},

where β is the physical inverse temperature, E and E' is the value of Hamiltonian (<ref>) for configuration respectively before and after the update.
The exchange of replicas with indices l and l+1 is performed using the Metropolis rule with excahge probability p_exch, which is given as

    p_exch=min{1,exp((β^(l)-β^(l+1))(H^(l+1)-H^(l)))}.


The histograms sampled with each replica, h_l(E,Q), are integrated to estimate the number of states W(E,Q) with the multiple histogram reweighting method.
W(E,Q) is obtained via the equation

    W(E,Q)=∑_lh_l(E,Q)/∑_l(ω_l(E,Q)/z_l∑_E=0^n∑_Q=0^nh_l(E,Q)),

where ω_l(E,Q)=exp(-β_l E), and z_l=∑_E=0^n∑_Q=0^nW(E,Q)ω_l(E,Q).

Eq. (<ref>) is originally derived in the paper 
<cit.>.


osmebibstyle






