extfieldsofstudy:
Mathematics
Computer Science

text:
Guidesort: Simpler Optimal Deterministic Sorting for the Parallel Disk Model

A new algorithm, Guidesort, for sorting in the uniprocessor variant of the parallel disk model (PDM) of Vitter and Shriver is presented. The algorithm is deterministic and executes a number of (parallel) I/O operations that comes within a constant factor $C$ of the optimum. The algorithm and its analysis are simpler than those proposed in previous work, and the achievable constant factor $C$ of essentially 3 appears to be smaller than for all other known deterministic algorithms, at least for plausible parameter values.

Introduction
Sorting is an important problem. In the seventies Knuth quoted an estimate that over 25% of computers' running time is spent on sorting [7]. Frequently the data sets to be sorted are so large that they do not fit in internal memory and must be held on external storage, often one or more magnetic disks or similar devices. In this setting of massive data, sorting acquires even more importance because a high number of algorithms that use external storage efficiently do so by reducing other problems to sorting, so that the time to sort almost can be viewed as playing the role that linear time has in RAM computation.
Accesses to magnetic disks are much slower than CPU operations and accesses to internal memory. A natural and increasingly popular way to sort faster is to use many disks in parallel. Magnetic disks have high latencies, so that efficiency dictates that an access to a magnetic disk must be used to transfer not one data item, but a whole block of many data items. The parallel disk model of Vitter and Shriver [14] tries to capture these characteristics of disk systems and has been used for much of the extensive research on sorting with several disks.

Model and Problem Statement
An instance of the uniprocessor variant of the parallel disk model or PDM of Vitter and Shriver [14] is specified via three positive integers, M, B and D with 2B ≤ M. It features a machine comprising an internal memory of M cells and D disks, each with an infinite number of cells. Every disk is linearly ordered and partitioned into block frames of B consecutive cells, each of which can accommodate a block of B data items. In slight deviation from the original definition of the PDM, we will assume that the internal memory is also linearly ordered and partitioned into block frames of B consecutive cells and-therefore-that M/B is an integer. An I/O operation or, for short, an I/O can copy the blocks stored in D pairwise distinct block frames in the internal memory to block frames on D pairwise distinct disks (an output operation) or vice versa (an input operation). Thus the disks are assumed to be synchronized. The description in [14] does not indicate explicitly whether it is also possible to transfer fewer than D blocks between pairwise distinct block frames in the internal memory and block frames on as many pairwise distinct disks; here we will assume this to be the case. Arbitrary computation (internal computation) can take place on data stored in the internal memory, whereas items stored on disks can participate in computation only after being input to the internal memory. An algorithm is judged primarily by the number of I/O operations that it executes, computation in the internal memory usually being considered free.
So that operations in the internal memory and I/O operations can specify their arguments, we assume that the disks, the cells in the internal memory and the block frames both in the internal memory and on each disk are numbered consecutively, starting at 0 (say). A number of conventions make the PDM convenient to argue about, but less precise. First, the meaning of a "cell" depends on the problem under consideration. In the context of sorting, as relevant here, a cell is the amount of memory needed to store one of the items to be sorted or a comparable object (e.g., for sorting it is generally assumed that one can ensure at no cost that the keys of the input items are pairwise distinct by appending to each its position in the input). Second, the only space accounted for is that taken up by "data items", not that needed to realize the control structure of an algorithm under execution. E.g., algorithms like those discussed in the following may want to manipulate such data structures as a recursion stack and various vectors indexed by disk numbers. Space for such bookkeeping information is assumed implicitly to be available whenever needed.
It is customary to restrict the parallel disk model by imposing an additional condition that says, informally, that M is sufficiently large relative to B and D. The condition varies from description to description, however, and seems to reflect the requirements of particular algorithms more than any fundamental deliberation concerning the model. E.g., [14] requires that D ≤ ⌊M/B⌋, [3] that M ≥ 2DB, and [11] that DB ≤ ⌊(M − M β )/2⌋ for some fixed β < 1. Since no more than M/B disks can take part in a (parallel) I/O operation, the condition D ≤ M/B seems somewhat more canonical than the others, and we will adopt it here (alternatively, D can be interpreted as the minimum of M/B and a true number of disks). It cannot be a priori excluded, however, that more than M/B disks can be put to good use in a PDM algorithm (cf. the RAMBO model of Fredman and Saks [5]).
A sequence of items is said to be stored in the striped format if its blocks are distributed over the disks in a round-robin fashion. More precisely, a sequence (x 0 , . . . , x N −1 ) is stored in the striped format if there are (known) nonnegative integers a 0 , . . . , a D−1 such that the following holds for i = 0, . . . , N − 1: If ⌊i/B⌋ = qD + r for integers q and r with q ≥ 0 and 0 ≤ r < D, then x i is stored on the disk numbered r in the block frame numbered a r + q. Every sequence for which nothing else is stated explicitly in the following is assumed to be stored in the striped format.

The Challenge
In the remainder of the paper, consider the problem of sorting a sequence of N > M items and take n = ⌈N/B⌉ (the number of block frames occupied by the input) and m = M/B (the number of block frames in the internal memory). Without loss of generality we will assume that the N items to be sorted have pairwise distinct keys. The last of the n blocks of input may contain a segment of data beyond the N items to be sorted. Such a "foreign" segment should be treated as a sorted sequence of dummy items larger than all real items, so that the segment will not be modified by the sorting.
In the sequential case, i.e., for D = 1, we can sort within A PDM machine with D disks can simulate one with a single disk but a larger block size of B ′ = BD without slowdown by operating the disks in lock-step, i.e., groups of corresponding block frames, one from each disk, are formed once and for all, and every I/O operation inputs from or outputs to the block frames in a group. Applying this to the sequential merging algorithm discussed above and assuming that m ′ = m/D is an integer and at least 2 yields a sorting algorithm that uses at most 2⌈n/D⌉⌈log m ′ (n/D)⌉ I/Os. Call this algorithm sorting by naive m ′ -way striping. If D is sufficiently small relative to m, there is no significant difference between logarithms to base m and logarithms to base m ′ . In particular, if D = O(m 1−ǫ ) for some fixed ǫ > 0, the bound of 2⌈n/D⌉⌈log m ′ (n/D)⌉ I/Os is within a constant factor of the second lower bound. As D approaches m, however, the problem becomes increasingly difficult. An attempt to parallelize the sequential merging algorithm simply by inputting the D next blocks in parallel meets with the difficulty that the D next blocks may not be known and, even if they are, may not be stored on D distinct disks. This could be called the problem of read contention.

The New Result
We present a new algorithm, Guidesort, for sorting in the PDM. The algorithm is deterministic, simple and easy to analyze and to implement. If the factor of 1 + F ′ (M, B, D, N) of the previous subsection is ignored, the number of I/Os executed by Guidesort is (C/D)Sort M,B (N), where C is approximately 3 for typical values of M, B and D-for brevity, the constant factor of Guidesort is 3. As D becomes large relative to m and B, C grows to a maximum of around 9.
Guidesort works by computing a guide that can be used to redistribute blocks to disks in such a way that the read-contention problem disappears. More details are provided in Section 2.

Previous Work
As befits a fundamental problem of great practical importance, a high number of algorithms for sorting in the PDM has been proposed. Many were qualified as simple. We are not aware, however, of any previously published algorithm that can be proved efficient using simple arguments. Some of the algorithms work via repeated merging [1,3,4,6,8,11,12] and hence bottom-up (this also applies to Guidesort), while others use the approximately inverse process of repeated splitting at chosen partitioning elements [6,10,13,14] and therefore operate in a top-down fashion.
Most published algorithms resort to randomization to cope with the readcontention problem. Some of them [3,6,13] are optimal and achieve constant factors close to 1 if D is sufficiently small relative to m (this is also when sorting by naive striping comes into its own), but as D approaches m their performance degrades.
Explicit constant factors were not indicated for any of the optimal deterministic algorithms published to date. The scheme of Aggarwal and Plaxton [1], based on Sharesort and bottom-up, but with some elements of top-down, is applicable to a variety of computational models. Its constant factor seems difficult to determine, but the generality of the approach lets one expect it to be quite large. Balance Sort by Nodine and Vitter [10] is a deterministic top-down algorithm that depends on subroutines for complicated tasks such as load balancing, matching and derandomization. Again, the constant factor appears to be large.
Greed Sort, also due to Nodine and Vitter [11], is perhaps closest in spirit to the Guidesort algorithm presented here. It is deterministic, based on repeated merging, and sufficiently simple that estimating its constant factor seems feasible. In order to merge, Greed Sort first carries out an approximate R-way merge, for a certain R, that brings each item to a position within some distance L of its rank in the sorted combined sequence, and then finishes by using Leighton's Columnsort [9] to sort locally within overlapping segments of 2L items each.
As indicated in [11], the number of I/Os executed by the approximate merge of Greed Sort is between 3 and 5 times the number of blocks involved. Subsequently each block participates in two applications of Columnsort. Columnsort, in turn, consists of 8 steps, each of which reads and writes every block at least once. Thus each recursive level needs at least 3 + 2 · 8 · 2 = 35 I/Os per block. In addition, R is chosen approximately as √ m, so that sorting based on R-way merging has about twice as many recursive levels as sorting by m-way merging. Since each recursive level in sequential sorting reads and writes each block just once (2 I/Os), this calculation indicates a constant factor for Greed Sort of at least 35. In fact, a more detailed study of [11] reveals the estimate of 35 to be optimistic, especially if D is not much smaller than m. The work described here was borne out of the author's desire to have an optimal sorting algorithm for the PDM simple enough to serve as the basis of a homework problem. Whereas this seemed out of the question for all previously published algorithms, the plan was carried out successfully for Guidesort.

Guidesort
Guidesort still sorts recursively, with the base case given by internal sorting of at most M items and each nonterminal call executing a multiway merge of recursively sorted sequences. Only the merge is done in a novel way.

The Main Idea
Call each sorted input sequence of a merge a run and consider a multiway merge of runs partitioned into blocks. Assume first that the blocks are input one by one. As long as there is enough internal memory to store one block from each run and to buffer the output, each input block must be read only once.
Define the leader of a block of items to be its smallest item. It is natural to read the next block from a run as soon as the last item in the previous block from that run has been consumed, i.e., moved to the output buffer. If we knew the leader x of the next block in the run, however, we could postpone reading that block until just before x is to be consumed. This shows that if the blocks of all runs are input in the order given by the sorted order of their leaders, it is still the case that each input block must be read only once. We can discover the appropriate canonical order by forming the sample of each run as the sorted sequence of its leaders and merging the samples, which creates the canonical sequence. Since the samples are far smaller than the full runs, the cost of merging them is usually negligible.
The computation of the canonical sequence was also considered by Hutchinson, Sanders and Vitter [6], who used the term "trigger" for what we call a leader.
Suppose now that we manage to input the blocks in the canonical order in batches of D ≤ D consecutive blocks each, where a batch is input in a single I/O operation. If we provide an additional buffer of D block frames to hold the latest input batch, it will still be the case that each input block must be read only once. In order for it to be possible to consume the input in batches, the blocks in each batch must be stored on different disks. In addition, to make it possible to produce each run D blocks at a time, we must ensure that if the run is partitioned into subsequences of D consecutive blocks each, the blocks within each subsequence are stored on different disks. Viewed more abstractly, we are faced with the problem of coloring the vertices of a graph G defined by the canonical sequence with at most D colors in such a way that no two adjacent vertices receive the same color. The vertices of G correspond to the leaders or the blocks in all runs, each edge joins two blocks that must be stored on distinct disks, and the D colors correspond to the D disks. The degree of G is bounded by 2(D − 1). If we choose D = D, coloring G may be difficult or impossible. With D ≤ ⌈D/2⌉, however, the obvious greedy algorithm can color the canonical sequence in a single pass. Call the resulting sequence of colored leaders the guide. While keeping the full guide, we also transfer the colors of leaders found to the original samples. This can be done with a recursive splitting that reverses the steps of the merge that created the canonical sequence.
To carry out the top-level merge, we first process each run, redistributing its blocks to the disks specified in its colored sample. This can be accomplished in a single pass over the run synchronized with a pass over its colored sample. Subsequently the top-level merge of the runs can proceed under the control of the guide, D blocks at a time, and this process can also create the corresponding sample for the merge at the next higher recursive level at essentially no additional cost.
To summarize, Guidesort comprises the following steps, illustrated in Fig. 1  While the back-of-the-envelope estimate essentially leads to the correct result, a more precise argument must account for the cost of Steps 1-3, which is not always negligible, and must also work without a number of assumptions that were made implicitly above. This is the topic of the next two subsections.

A Realistic Simple Special Case
Define the block length of a stored sequence to be the number of block frames that it occupies (fully or in part). Let us be more specific about the recursive sorting algorithm, which is parameterized by an integer r with 2 ≤ r ≤ m. When asked to sort a sequence of block length p, the algorithm first computes k = min{⌈p/m⌉, r}. If k = 1, the sorting takes place in the internal memory without recursive calls. Otherwise the input sequence is split into k subsequences, each of block length ⌊p/k⌋ or ⌈p/k⌉, and the k subsequences are sorted recursively and subsequently merged. The main point is that although the algorithm "nominally" uses r-way merging, a call that issues terminal calls issues only as many as necessary. Without increasing the depth of the recursion, this ensures that each terminal call deals with at least m/2 blocks and therefore that the number of terminal calls and the number of calls altogether are O(n/m). This is useful for dealing with rounding issues. E.g., if each merge executed as part of the overall sorting is associated with a quantity of the form ⌈t⌉, where t is some expression, we can upper-bound the sum of ⌈t⌉ over all merges by summing t over all merges and adding O(n/m).
Our analysis frequently sums the number of I/Os of a particular kind over all merges carried out as part of the overall sorting. We shall use the term "accumulated" to denote this situation, i.e., "accumulated" means "summed over all merges". Let z be the accumulated total block length of the runs of all merges. Exactly as in the sequential sorting algorithm, z ≤ n⌈log r n⌉. Similarly, let y be the accumulated block length of all samples produced by the algorithm. Since each leader in a sample "represents" B items, an application of the counting method set out at the end of the previous paragraph shows that y = z/B + O(n/m).
Step 1 merges at most r samples. If the block length of some sample is < m, and therefore the block lengths of all samples are ≤ m, the samples are first partitioned into bundles in such a way that the total block length of the runs in each bundle, except possibly the last bundle, is at least m/2 and at most m. The blocks in each bundle are then input, D at a time, sorted internally, and output as one sorted sequence, again D blocks at a time, at an accumulated cost of 2y/D + O(n/m) I/Os. Assume that the number of sorted runs (original samples or sorted bundles) at this point is r ′ . The r ′ runs are merged in a binary merge tree of height ⌈log 2 r ′ ⌉ in which each binary merge is carried out with D 1 disks that operate in lock-step to simulate a single disk with a block size of B ′ = D 1 B.
Observe that n/m ≤ n/D = O((z/D)(1/log m n)) (a convexity argument shows that we even have n/D ≤ (z/D)(1/log m n)). When r ′ runs are merged in a binary merge tree, each run, except possibly the last one, contains at least M/2 leaders that "represent" at least MB/2 input items. With respect to a block size of B ′ , the accumulated block length of all runs input to nontrivial binary merge trees-those with r ′ ≥ 2 leaves-is therefore O(z/(D 1 B)) = O(z/(DB)), and the merges can be carried out with O(z⌈log 2 r⌉/(DB)) = O((z ln m)/(DB)) I/Os. Altogether, the number of I/Os executed by Step 1 is O((z ln m)/(DB) + n/m). This bound also covers Step 3, which reverses the merging to transfer the information attached to leaders in Step 2 from the guide to the samples. Steps 1 and 3 need 2D 1 block frames of internal memory, which are available since D 1 ≤ m/2.
Step 2 colors the canonical sequence. Say that a color is used recently in a sequence of k colored objects if it is the color of one of the min{D − 1, k} last objects in the sequence. The leaders are processed in the order in which they occur in the canonical sequence, and each leader x is given a color in {0, . . . , D − 1} that is used recently neither in the complete sorted sequence of colored leaders nor in its (not necessarily contiguous) subsequence of leaders drawn from the same sample as x. To accomplish its task, the algorithm maintains at most r + 1 history sequences of the chronologically ordered recently used colors overall and within each sample. To process a leader, the algorithm identifies and assigns an appropriate color to the leader and updates two history sequences accordingly, all of which is straightforward. In addition, the leaders of each color are numbered consecutively in the order in which they are colored and the number of each leader, called its index, is attached to the leader. We will assume that the history sequences with their at most (r + 1)(D − 1) color values can be stored in q = ⌈(r + 1)(D − 1)/B⌉ block frames. Using an additional input/output buffer of D 2 block frames, we can execute Step 2 with an accumulated number of I/Os of 2y/D 2 + O(n/m) = 2z/(D 2 B) + O(n/m). We must ensure that q + D 2 ≤ m.
Step 4 redistributes blocks to new disks one run at a time. When processing a run of ℓ blocks, the algorithm inputs the ℓ blocks and, interleaved, a sample of ℓ colored leaders and, for i = 1, . . . , ℓ, stores the ith block on the disk corresponding to the color of the ith leader and in a block frame whose number is the index of the ith leader plus an offset chosen for that disk. Thus the blocks on each disk are stored compactly, but they are not necessarily written in the order of increasing frame numbers. With a "primary" input buffer of D 4 block frames for the input blocks, a "secondary" input buffer of D L block frames for the leaders and an output buffer of D block frames for the redistributed blocks, the accumulated number of I/Os spent in Step 4 is z/D 4 + z/D + y/D L + O(n/m) = z(1/D 4 + 1/D + 1/(D L B)) + O(n/m). We must ensure that D 4 + D + D L ≤ m.
Step 5 actually merges the runs with the aid of one block frame of input buffer for each run, a total of at most r block frames, an input buffer of D block frames for the latest batch (D − 1 block frames actually suffice), two buffers of D L block frames each for the guide, which is input, and the sample for the next recursive level, which is output, and finally a buffer of D 5 block frames for the primary output, the final outcome of the merge. The accumulated number of I/Os is z/D + z/D 5 + 2y/D L + O(n/m) = z(1/D + 1/D 5 + 2/(D L B)) + O(n/m), and we must ensure that r + D + D 5 + 2D L ≤ m.
Collecting the contributions identified above and not forgetting the O(n/D) I/Os consumed by terminal calls of the recursive sorting algorithm, we arrive at a total number of I/Os for the complete sorting of

The General Case: An Order-of-Magnitude Bound
In order to deal with situations in which B is smaller than D or ln m is not small relative to B, we generalize the algorithm by introducing an additional parameter, s, which must be a positive divisor of D. The algorithm described so far corresponds to the special case s = 1.
Whereas until now a leader was the smallest item within a block, we redefine it to be the smallest item within a segment of s consecutive blocks in some run. Thus each run is partitioned into segments of sB items each, except that the last segment may be smaller, and each segment contributes only a single leader to the sample of the run. As a consequence, the accumulated block length of all samples now is y = z/(sB) + O(n/m). We stipulate that the s blocks that form a segment are colored using consecutive colors that begin at a multiple of s. An intuitive view of this is that the task now is to color segments with one of 2D/s − 1 colors so as to avoid the at most D/s − 1 most recent colors both overall and within the same run. This has the beneficial effect of reducing the state information that must be kept by the greedy coloring in Step 2 from (r + 1)(D −1) colors to at most (r + 1)(D − 1)/s colors, which can be stored in q = ⌈(r + 1)(D − 1)/(sB)⌉ blocks.
On the other hand, in order for the leaders to fulfil their function, Step 5 must be changed to input runs whole segments at a time, which means that an input buffer of s block frames rather than of a single block frame must be provided for each run. The requirement for Step 5 accordingly becomes rs + D + D 5 + 2D L ≤ m for the revised sorting algorithm, while the total number of I/Os generalizes from (1) to For m ≤ 3 and therefore D ≤ 3, sequential sorting obviously works within  N)). If we want to prove a more precise bound, we must choose the parameters more carefully.

The General Case: Good Constant Factors
In order to obtain the best result, we change the algorithm slightly: In Step 4, instead of having separate primary input and output buffers of D 4 and D block frames, respectively, we use a single buffer of D 4 block frames for both input and output. This is trivial, but requires us to ensure that D | D 4 (i.e., D divides D 4 ). Besides r, s, D, D 2 , D 4 , D 5 and D L being positive integers, the conditions that the parameters must satisfy for the final algorithm are the following: rs Our task at this point is to (approximately) minimize the I/O bound of (2) subject to the constraints (3). As an aid in dealing with the divisibility requirements of the last constraint, we first prove a simple technical lemma that, informally, says that, given two positive integers a and b, we can make the smaller divide the larger without changing their sorted order by lowering each by less than half and b by less than its square root. Let N = {1, 2, . . .}.
There is a function f : N 2 → N 2 that can be evaluated in constant time and has the following properties : For all a, b ∈ N, if (a ′ , b ′ ) = f (a, b), then From now on assume that 2 ≤ a ≤ b − 2. Since then b ≥ 4 and √ b ≤ b/2, we need no longer verify the condition b ′ > b/2 explicitly. Consider two cases: If a ≤ √ b + 1, take a ′ = a and b ′ = ⌊b/a⌋a ≤ b and observe that if a < √ b + 1, , take a ′ = ⌊b/⌈b/a⌉⌋ ≤ a and observe that if a ≥ b/2, then a ′ = ⌊b/2⌋ ≥ ⌊(a + 2)/2⌋ > a/2, whereas if a < b/2, then we can successively conclude that a ≥ 4, b ≥ 9 and a ≥ 5 and then that a ′ ≥ ⌊b/(b/a + 1)⌋ = ⌊ab/(b + a)⌋ ≥ ⌊2a/3⌋ > a/2. Also take and note that b ′ ≤ b is obvious from the third expression in the previous line.
The stage is set for our main result: Proof. Because of the term 1/ √ D in the expression for F ′ , we can assume without loss of generality that m is larger than a certain constant m 0 , which we choose as m 0 = 2 12 . As justified earlier, we will also assume that D ≥ √ m. What remains is to prove that r ≥ 2 or, equivalently, that r 2 ≥ 4 and r 5 ≥ 2. First, and hence, by the second property of Lemma 2.1, s ≥ (1/2) 7D/(16B). Now (m − 1)sB/(D − 1) ≥ 2sB ≥ 7D/16 ≥ 5, where the last inequality follows from D ≥ √ m ≥ 2 6 , and this shows that r 2 ≥ 4.
In order to demonstrate that r 5 ≥ 2, we prove that m − D − D 5 − 2D L ≥ 2s.
We next bound the number of I/Os executed by the algorithm, essentially by estimating the values of its parameters. First By the third property of Lemma 2.1, it follows that As observed earlier, (m − 1)sB/(D − 1) ≥ r 2 +1. But r 2 +1 ≥ 2r+1 ≥ (3/2)(r+1) and therefore (r + 1)(D − 1)/(sB) ≤ (2/3)(m−1). As a consequence, D 2 = Ω(D). Because D 4 is just 2D, Noting that for x ≥ 1, we find that We are almost ready to sum the terms in (2). Observe first that since D 4 , D 5 and D are all larger than ǫD for some fixed ǫ > 0, our bound of the form D 4 = D(1 − O(· · · )) yields a bound of the form 1/D 4 = (1/D)(1 + O(· · · )), and analogously for D 5 and D. We already noted that sB = Ω( √ DB). Hence We also have z ln m sDB Now the total number of I/Os can be seen to be By the second property of Lemma 2. where F ′ is as in the theorem. To show that r 2 ≥ 4, the proof of Theorem 2.2 argued that 2sB ≥ 5, a relation that is a triviality for B ≥ 16. In the case of r 5 ≥ 2, the relation to show was m > 16s − 10. But for 8 ≤ m ≤ 16 this is clear since s = 1, while for m ≥ 16 it follows from s ≤ D/B ≤ (1/4) √ m.
It is easy to see that the asymptotic assertions of Theorem 2.2 continue to hold.

Conclusion
Guidesort is simple and based on a natural idea. If a reasonably programmed algorithm for sorting through multiway merging is available, the new algorithm can be grafted onto it by replacing the old subroutine for multiway merging by one that executes Steps 1-5.
Step 1 is standard merging. Because the sequences to be merged are smaller by a factor of B than the actual runs, in many practical situations the merging can be done any which way. The algorithm proposed here that can be viewed as the final part of sorting by naive 2-way striping is also very easy to implement. Step 3 is an even simpler approximately inverse operation.
Step 5 is standard merging, except that the tests to identify the disks from which to input the next blocks have been executed beforehand.
Step 2 implements a straightforward greedy algorithm, and Step 4 is a trivial redistribution of blocks according to a precomputed pattern.
While most or all previous algorithms are optimal with respect to internal computation, i.e., execute O(N ln N) steps of internal computation to sort N items, it is not obvious how to carry out the greedy coloring of our Step 2 in less than Θ(D/s) time per leader (or Θ(D/(sw) + 1) time per leader if w bits can be manipulated together in constant time), for which reason we can bound the amount of internal computation of Guidesort only by O(N ln N +(nD/s) log m n) = O(N(ln N + D/B 3 log m n)). For realistic parameter values, however, this is still O (N ln N).
It may seem deplorable that the function F ′ should include an "error term" of O(1/log m n), which can be quite large even for realistic values of the parameters n and m. However, recall that the term represents just O(n/D) I/Os, i.e., the cost of a single pass over the input at full parallel speed. Approximating the total depth of the leaves in the recursion tree of a bottom-up or top-down algorithm similar to those published to date by ⌈log m n⌉ times the number of leaves incurs an error of the same magnitude. An "error" of O(n/m), which is comparable if D ≈ m, seems even harder to avoid: The number of nodes in the recursion tree is Θ(n/m), so if Ω(1) I/Os are "lost" on average at every node, e.g., because of a last block that is only partially occupied by items, the accumulated waste amounts to Ω(n/m) I/Os.


extfieldsofstudy:
Physics
Medicine
Computer Science
Mathematics

text:
A Java Application to Characterize Biomolecules and Nanomaterials in Electrolyte Aqueous Solutions

The electrostatic, entropic and surface interactions between a macroion (nanoparticle or biomolecule), surrounding ions and water molecules play a fundamental role in the behavior and function of colloidal systems. However, the molecular mechanisms governing these phenomena are still poorly understood. One of the major limitations in procuring this understanding is the lack of appropriate computational tools. Additionally, only experts in the field with an extensive background in programming, who are trained in statistical mechanics, and have access to supercomputers are able to study these systems. To overcome these limitations, in this article, we present a free, multiplatform, portable Java software, which provides experts and non-experts in the field an easy and efficient way to obtain an accurate molecular characterization of electrical and structural properties of aqueous electrolyte mixture solutions around both cylindricaland spherical-like macroions under multiple conditions. The Java software does not require outstanding skills, and comes with detailed user-guide documentation. The application is based on the so-called Classical Density Functional Theory Solver (CSDFTS), which was successfully applied to a variety of rod-like biopolymers, rigidlike globular proteins, nanoparticles, and nano-rods. CSDFTS implements several electrolyte and macroion models, uses different level of approximation and takes advantage of high performance Fortran90 routines and optimized libraries. These features enable the software to run on single processor computers at low-to-moderate computational cost depending on the computer performance, the grid resolution, and the characterization of the macroion and the electrolyte solution, among other factors. As a unique feature, the software comes with a graphical user interface (GUI) that allows users to take advantage of the visually guided setup of the required input data to properly characterize the system and configure the solver. Several examples on nanomaterials and biomolecules are provided to illustrate the use of the GUI and the solver performance. aEmail: marcelo.marucho@utsa.edu, website: http://physics.utsa.edu/Faculty/Marucho/Marucho.html 1 ar X iv :1 80 1. 07 20 1v 1 [ ph ys ic s. co m pph ] 1 9 Ja n 20 18

I. INTRODUCTION
A rich and complex, yet not fully understood, electrical double layer (EDL) formation arises when a nanoparticle or biomolecule is immersed in a liquid solution [1][2][3]. In such conditions, the structural and thermodynamic properties of the liquid surrounding the macroion are shown to be very different from those corresponding to the bulk phase, which in turn, have a high impact on the behaviors and functions of the macroion. For instance, physicochemical properties of nanomaterials such as particle size (PS), shape, Zeta potential (ZP), and surface charge density (SCD) have significant effects on their stability, circulation in blood and absorption into cell membranes [4][5][6][7]. Alteration of the pH and electrolyte conditions modify the SCD and ZP, which affect the binding of nanomaterials to biological tissue and target them to specific sites inside the cell. Whereas, in biophysics [8][9][10], the formation of high order structures, including bundles and networks, involve interactions between highly charged linker proteins and biopolymers that are often modulated by the biological environment. The environment has been shown to have a profound impact on biological functions. For instance, ZP and SCD in microfilaments, nucleic acids, and globular proteins play a fundamental role in the stabilization of these systems due to charge residue groups on their surface. Extremely positive or negative ZP values cause larger repulsive forces. This repulsion between macroions with similar electric charge prevents aggregation of the macroions and thus ensures easy redispersion. Adversely, altering the ZP to the point at which it exhibits zero net charge may decrease stabilization forces and promote agglomeration.
Beyond the substantial progress done in the characterization of these colloidal systems, the complex interplay between the physicochemical properties of macroions and biological environment still remains elusive due to the lack of appropriate methodologies. Conventional computational tools and approaches are limited by their approximations and computational cost. Certainly, the current understanding of EDL properties of macroions is based mainly on mean-field theories like the non linear Poisson-Boltzmann (NLPB) formalism and its modifications [10], which consider electrostatic potential interactions only. Furthermore, the poor resolution and sensitivity of current experimental techniques limit our ability to obtain information on the molecular mechanisms governing these phenomena. Thus, they may be inappropriate to explore a large number of colloidal systems. Additionally, scientific software developed to characterize solvation and electrical double layers of nanomaterials and biomolecules might be helpful in understanding these phenomena. However, they usually require specialized training and expertise in computational biology, expensive commercial licenses, and access to clusters or supercomputers, which are often an obstacle for many researchers, experimentalists, even students lacking these requirements. Thus, it is essential to develop not only more accurate and efficient approaches, but also readily accessible, friendly implementation of a software package that provides experts and non-experts a visualized guide (graphical user interface) to perform these calculations without limitations [11].
In this article, we present a free, multi-platform and portable Java software which provides both experts and non-experts in the field an easy and efficient way to get an accurate molecular characterization of the EDL properties for biomolecules and nanomaterials at infinite dilution. The application is based on the so-called Classical Solvation Density Function Theory (CSDFT) and its modifications, which have been shown to be particularly useful in studying multiple environmental scenarios for a variety of rod-like [12,13] and spherical [14][15][16] rigid-like macroions without computational restriction. CSDFT extends the capabilities of NLPB formalism, eliminating the extremely high computational demands of full atomistic simulation calculations without losing important structural features of complex EDL properties (see Figure 1). It considers not only the electric, but also the entropic and many-body interactions. This feature has been particularly useful for identifying and characterizing dominant interactions and molecular mechanisms governing the behavior and function of In the next sections, we explain how to install, run and use the software. We also include examples to illustrate the solver performance and applicability. In Appendix A, we provide a brief summary of the approach and computational scheme implemented by the software to solve the CSDFT numerically.

II. OVERVIEW
The software allows the user to use a variety of models and approaches. It provides four macroion models: nanoparticle, nanorod, globular protein, and rod-like biopolymer. For each of these models, there are two approaches to characterize the macroion size and the surface charge density: experimental and protonation/deprotonation models for nanomaterials applications, and experimental and molecular structure models for biophysics applications. Additionally, the software comes with two electrolyte approaches: NLPB and CSDFT. These approaches include 3 ion size types (Crystal, Effective and Hydrated), and two solvent models (implicit and explicit). As a unique feature, the Java application comes with a graphical user interface (GUI) which eliminates the arduous and error-prone manual entry of data, and substantially reduces the time spent on the setup of the information required to characterize the macroion and the electrolyte solution. Additionally, each GUI screen provides helpful information about how to fill out the input data by simply holding the mouse pointer over the corresponding text or blank box. The GUI also provides default values for key input parameters and preselects some relevant algorithms to speed up the setup of the input data. However, they may be easily changed at any time. Moreover, the GUI tests all the input data before running the CSDFTS to avoid the incorrect use of the software and prevent meaningless results. At the end of the calculations, the GUI generates two-dimension plots of the selected output files to provide graphical representations of the structural and electrostatic properties of the macroion EDL. Finally, all the output data files are properly saved and organized in the designated folder for post-analysis purposes.

System Requirements
The software requires at least 500 MB of disk space to install. Additional disc space is required to save the output files generated by each simulation, which may be as large as 1Gb or more. The available RAM memory required to run the software depends on many factors. Based on the examples presented in this article, most simulations require only a few Gbs whereas highly demanding calculations may demand up to 12 Gbs.
The software is portable, e.g. it does not require external libraries, leave its files or settings on the host computer, or modify the existing system and its configuration. The software comes with Java 1.8 which is configured by the installer to run the GUI. The Fortran90 codes are distributed as binary files. For biophysics applications, the software requires the molecular structure visualization java tool "Jmol" (http://jmol.sourceforge.net/) [17] and either the java based protein volume "provol" (gmlab.bio.rpi.edu/) [18] or the Unix based "v3" (http://3vee.molmovdb.org/volumeCalc.php) [19] applications. These open sources are compiled and distributed along with the software. To use the latter application, the user's computer needs to have the GNU "gawk" application installed (https://www.gnu.org/software /gawk/). The software also includes "pdb2pqr" and "propka" binary files with the necessary libraries to perform molecular structure calculations (http://apbs-pdb2pqr.readthedocs.io /en/latest/pdb2pqr/index.html) [20]. Additionally, the GUI requires a variety of the user's computer applications, which are a part of the basic operating systems. In particular, the GUI uses "Activity Monitor", "procexp.exe", and "gnome-system-monitor" (or "top") applications for Mac, Windows and Linux users, respectively, to monitor the user's computer performance. It also uses "Terminal.app", "cmd.exe" and "gnome-terminal" (or "xterm") for Mac, Windows, and Linux users, respectively, to display CSDFT calculations.

Installation
The software is distributed as a single zip file containing the self-installer application. The user does not need "Root" or "Admin" permissions to install the software on the user's home directory. The user may unzip the file on the desktop directory, open the resulting folder, then double-click the self-installer application to run the wizard setup and follow a number of module steps (See Figure 2). The installer creates the folder "CSDFTS" in the chosen directory to save the software source files, including the GUI executable "CSDFTSsoftware". If selected, the installer also creates an icon on the desktop for easy access to run the GUI. The software was successfully installed and tested on a variety of operating systems (Windows 7 and 10, Sierra, Centos 7, Ubuntu 14, Fedora 24, Debian 8) and provides an uninstaller application for easy removal. During the first simulation the GUI creates a folder named "CSDFTS_Workspace" in "My Documents" and "$HOME" directory for Windows and Unix operating systems, respectively. Within this folder, the GUI creates one sub-folder per each simulation to save and organize all the corresponding output data files. The name of the sub-folders follows the convention: Analysis_<<macroion-electrolyte>>_<<Date>>__<<Time>>, where <<macroion-elec-trolyte>> is the macroion and electrolyte model selected by the user, <<Date>> is the current date, and <<Time>> is the real time when the model was selected (see Figure 3).
The GUI generates the following output files: "DensityProfile.dat (normalized ions and water density profile distributions"); "ElectrostaticPotential.dat" (mean electrostatic potential); "IntegrateCharge.dat" (the integrated charge), and "IonContributions.dat"(normalized electrostatic potential energy, particle crowding entropy energy, and ion-ion electrostatic correlation energy to the ionic potential of mean force). Explicit expressions for these properties are provided in Appendix A. The size of these files may vary depending on the number of ions species, number of grid points, and electrolyte model. These files are formatted in a multi-column arrangement where the row number corresponds to the number of grid points, the first column contains the discretized distance in unit of Å and the remaining column(s) provide the numerical solution(s) (see Table I). Other relevant data including SCD, ZP, and PS are saved in the log file "out.txt", whereas information on the input data and solver configuration are included in the files "inp.in" and "input.file" for spherical and cylindrical macroions, respectively. If the molecular structure option is selected for biophysical applications, the GUI generates the following additional files: "input.pdb" (copy of the molecular structure uploaded by the user); "input.pqr" (pdb2pqr output data); "input.propka" (propka output data) ; "out.pdb" (molecular structure information with heterogem atoms removed); "out.xyzr" (atomic positions in xyzr format); "vol.txt" (3v and provol output data including information on the macroion surface and volume as well as the number of residues and atoms); "lengthandcharge.txt" (total macroion charge and dimensions), and "radiusandcharge.txt" (parameters characterizing the macroion radius and charge density).

A. GUI Description
In this section we describe the screen sequence generated by the GUI. The first, second and third screens correspond to "project window", "model window", and "results visualization window", respectively. Each screen provides information to help the user fill out the input data by moving the mouse pointer over the corresponding text or blank box (see Figure 4).

Units Description
DensityProfile.dat The first column is for distances. The second column corresponds to the first ion species selected in ion model section, and so on.
If the CSDFT is selected, the last column corresponds to water ElectrostaticPotential.dat KT/e The first column is for distances. The order of the columns in each set of contributions is the following: first column corresponds to the first selected ion species, and so on. The Project window shown in Figure 5 is the main screen which provides user access to the Menu, Information and Model Sections.
The menu section, located at the top left corner of the window, contains the File, Tools and Help menus. The File Menu contains the Results Visualization and Exit options (see Figure 6(a)). The Results Visualization option may be used after running the simulation. This option allows users to select output file(s) and visualize the solution(s) in two dimensional plot(s).
The Tools Menu contains the Browse PDB files, Crystal Ion Library, Effective Ion Library and Hydrated Ion Library options (see Figure 6(b)). The tool Browse PDB files is used for biophysical applications only. It allows users to open a web browser and download protein molecular structure(s) from the protein data base web page (http://www.rcsb.org/). The user may use this tool at anytime. The Ion Library tools provide the user access to the tables pretabulated with specific information on ion species, valences and diameters that are required to characterize the electrolyte aqueous solution. These tools can be used to redefine existing ion species and define new ion species as well. The last menu in Figure 6(c) is the Help Menu which provides access to the user's guide, software license and contact information.
The Solute Geometry Model section is located at the central-bottom part of the window. To select the macroion (solute), the user has four options: Nanoparticles and nanorods for nanomaterials applications, or globular proteins and rod-biopolymers for biophysical applications.
The Electrolyte Aqueous Solution Model is located right below the Solute geometry model. CSDFTS offers two electrolyte aqueous solution theories. NLPB uses an implicit solvent model and considers electric interactions only, whereas, CSDFT uses an explicit solvent model (SPM) and considers not only the electric, but also the entropic and ion-ion correlation interactions. More information on solvent and ion models is provided in Appendix A. The selection of the EDL model mainly depends on the computational resources, the electrolyte conditions, physicochemical properties of the solute, and the accuracy required for the numerical solution (illustrative examples are provided below). As a rule of thumb, NLPB is more efficient, but less accurate than CSDFT. It is worth mentioning that CSDFT is a unique feature of the software, whereas NLPB theory is included for testing and comparison purposes, mostly. Indeed, there are other efficient programs based on implicit solvent models such as APBS [21] and MPBEC [22] who provide the solution for both symmetric and asymmetric macroion shapes.

Model Window: system and solver configuration
Once the Solute Geometry Model and the Electrolyte Aqueous Solution Model are selected in the project window, the GUI automatically opens a second screen, the model window. In this screen, most of the text based user interaction occurs (see Figure 7). It contains the following modules: Electrolyte Model, Solute Model, and Numerical Scheme and Runtime Options.
The Ion Model section shown in Figure 8 allows users to characterize the electrolyte solution by providing the bulk concentration [ρ 0 i ] and the valence z i for each ion species i. In NLPB theory, ions are represented as point-like particles, and consequently, there is no need to define the ion sizes. Whereas, in CSDFT each ionic species i is represented by  charged hard spheres of diameter d i . In this case, CSDFTS offers a pretabulated crystal, effective and hydrated ionic diameter types, which are estimated using different experimental techniques [23]. To characterize the electrolyte solution, the user has to select the first ion species and its properties, then click the button "Add", and subsequently repeat the procedure to add more ion species. The "Delete" and "Update" buttons allow users to remove and change the properties of a previously selected ion. Once the user selects all the ion species (e.g. valence and bulk concentration), the electroneutrality condition (e.g.    For nanomaterials it is required to provide the effective radius which is usually obtained from experiments. For biomolecules there are two ways to characterize the solute size. The user may either provide this information from experiments (Experimental data model ) or it can be estimated from the molecular structure (Molecular structure model ).
A novel feature of the software is the option to calculate the macroion surface charge density σ arising from the protonation/deprotonation reactions of the dissociable functional groups at the solid/liquid interface. The user can also click the Illustration buttons to visualize the different options and models (see Figure 11). For nanomaterials, this charging / discharging mechanism has traditionally been described using surface complexation models (SCMs), where titration of surface groups is described by the pH level and an ensemble of mass balance equilibria with associated equilibrium constants (see Appendix A). This approach is included in the protonation/deprotonation model where the user has to provide the equilibrium constants pKa and pKb for protonation and deprotonation of the active functional group, respectively, as well as the number density of total functional groups on the surface N, and the pH level (see Figures 10 a)). The resulting SCD depends on the ZP and vise-versa. A different scenario is presented for globular proteins and rod-biopolymers since the titration is usually described by a set of empirical rules relating the protein structure to the equilibrium constant values of ionizable groups (residues). This mechanism is included in the Molecular structure model option where the user has to upload the molecular structure (pdb file), select a force field, and define the pH level as shown in Figure 10 b). This uncharged molecular structure in pdb format is used by the software to automatically run the pdb2pqr / propka application (http://sourceforge.net/projects/pdb2pqr/). This application assigns atomic charges and sizes, adds missing hydrogens, optimizes the hydrogen bonding network, and renormalizes atomic charges of the residues exposed to the surface due to pH effects (protonation/deprotonation process). Afterwards, the GUI uses the molecular structure in pdb format to automatically run either the provol (default) or 3v application to estimate the total macroion volume by rolling a probe particle of radius 1.4Å and using a grid mesh of 0.5Å. Finally, the GUI uses the information obtained from the total protein charge and volume to estimate the effective solute radius as well as the uniform surface charge density and Bjerrum length for spherical and cylindrical shapes, respectively. For each of these approaches, the GUI opens an external terminal window to display these calculations (see Figure 12). It also opens a RAM memory window to monitor the user's computer performance (see Figure 13). The GUI closes these windows once the calculations are over and enables the execution buttons "Generate the input file and run CSDFT manually" and "run CSDFT automatically" to run CSDFT.
Another important feature for biophysics applications is the molecular visualization tool. After all the aforementioned calculations are over, the user can visualize the uploaded Figure 11: Illustrations molecular structure which provides a detailed molecular characterization including the amino acid sequence and the number and type of residues exposed to the electrolyte (see Figure 14). Alternatively, the user may select the experimental data option to provide the values of these parameters from experiments. In this case, there is no waiting time to run CSDFTS.
The Thermodynamic Conditions option is used to define the electrolyte solution temperature. The default value displayed in Figure 15 corresponds to the room temperature.
The Numerical Scheme and Runtime Option section shown in Figure 17 allows the user to configure key parameters that play a fundamental role in the solver performance. They control the accuracy and computational cost. The tolerance number represents the numerical error required by the user to numerically obtain the normalized density profile solutions. A minimum of six digits of precision is highly recommended (default value). The radial grid resolution represents the regular separation distance h between two consecutive points in the domain discretized of the radial distance to solve CSDFT numerically. The value recommended for this parameter is 0.05Å which has been shown to work for most applications on Figure 12: CSDFTS calculations computers without RAM memory restrictions. The domain of the solution ranges from the macroion surface (e.g. the radius R) to the cutoff L. The latter is determined automatically by CSDFTS to provide the correct (long range) asymptotic behavior of the mean electrostatic potential, and consequently, satisfy the electroneutrality condition of the system. The number of grid points is calculated as follows N = (L − R)/h + 1. The Mixing Parameter is a number between 0 and 1 which helps the solver to reach stability and convergence in the solution iteratively. The default value for this parameter is 0.35. It is required for the protonation/deprotonation model only.
We note that the high resolution and tolerance generate high computational cost, allocated RAM memory, solver stability and accuracy, but slow the rate of convergence on the calculations. The user may use lower resolution to speed up the calculations and reduce the allocated RAM memory at the risk of loosing convergence in the CSDFTS iterations. In that case, the solver will stop the simulation after 10000 iterations. The user may also kill the process using standard approaches (e.g. "Ctrl+c", etc). More details on the numerical solver scheme are provided in Appendix A.
After all changes on preselected parameters and options are completed, all boxes are filled out and precalculations are over, the user may either press the button Generate the input file and run CSDFT manually or run CSDFT automatically. The GUI tests all the input data, and sends the user warning messages if unusual/nonphysical values are assigned to the required parameters or missing information is detected (see Figure 16).
We highly recommend the user to use the option run CSDFT automatically. The other option is mainly for those users that have computer skills and prefer to run the simulations using command line. In this case, the GUI generates the input file "inp.in" for spherical solutes and "inputfile.in" for cylindrical solutes. To run CSDFTS manually in Linux operating systems, for instance, the user should open a terminal, change the current directory to the corresponding output sub-folder directory and type "~/CSDFT/Linux/sphere_linux.sh Figure 13: User's computer performance on Linux operating systems inp.in" for spherical solutes or "~/CSDFT/Linux/cylinder_linux.sh inputfile.in" for cylindrical solutes. A similar procedure should be performed for other operating systems.

Results Visualization
If the "run CSDFT automatically" option is selected, the results visualization window appears once the simulation is over. This screen allows the user to select the output file(s) and visualize the numerical solution(s) as a function of the radial separation distance in two dimension plot(s) (see Figures 18). To change the plot properties the user must right click on the desired plot (see Figure 18 c)). The user can also visualize the data on spreadsheets as shown in Figure 18 b).

B. Quick start
The user is recommended to follows these simple steps: • On the first screen, select the macroion type and subsequently the electrolyte model.
• On the second screen, select the ions and provide the corresponding concentrations.
Make sure to fulfill the electroneutrality condition. Next, define the solvent model.  The dielectric permittivity value may be changed if needed. If CSDFT was selected, the solvent concentration and molecule size may be also changed. Subsequently, configure the solver. The resolution and tolerance may be changed depending on the user's computer performance and available RAM memory. The default values may be used if there is no RAM memory restrictions. Otherwise, decrease the resolution accordingly. For nanomaterials applications and protonation / deprotonation models, the mixing parameter may be changed if needed. For biophysics applications and the molecular structure option, the protein volume application may be changed depending on the user's operating system. The default application works for all the operating systems whereas "3v" only works for Unix operating systems and requires the GNU gawk application installed on the user's computer. As a rule of thumb, the latter runs faster and requires less allocated RAM memory. Next, the temperature may be changed if needed. The last step is the selection of the solute model. For nanomaterials applications, select either the experimental or protonation / deprotonation option and provide the corresponding information. Now the software is ready to run the CSDFTS. For bio-  physics applications, select either the experimental or the molecular structure option and provide the corresponding information. In the latter case, select the force field and set up the pH before uploading the molecular structure. Wait until some calculations on the molecular structure are over. Afterwards, the molecular structure my be visualized. Finally, select either the "Generate the input file and run CSDFT manually" or "run CSDFT automatically".
• If the "run CSDFT automatically" was selected, an additional screen is available to visualize solutions in 2D plots.
• If there is any issue in the previous steps, close all screens and start over again. Double check that all the input data and selected options are appropriate for the required simulation and change those parameters to improve the CSDFTS convergence and the user's computer performance. In this section we analyze several macroions under different electrolyte conditions (see Table II). We ran these simulations on several operating systems and CPU processors. The results shown in the present article corresponds to those obtained on a single Intel Xeon 5680 and CentOS 7 operating system.
The   To illustrate numerical solutions, we show in Figure 19 the structural properties of the EDL predicted by CSDFT and NLPB for a myoglogin protein, whereas in Figure 20, we present the electric properties of the EDL corresponding to a actin filament . Certainly, this visualization is useful for analyzing the behavior of the electrostatic potential interactions and the ion-water distributions at short and long distances from the macroion. It is also relevant for analyzing the effective range, intensity, and nature of the driving force governing the potential of mean force for each ion species. An extensive analysis and discussion on structural and electrical properties of EDLs can be found elsewhere [10,12,[14][15][16]24].
Overall, the software offers diversity in applicability and analysis on many colloidal systems. For instance, running multiple simulations at different pH levels may provide a molecular understanding on the impact of pH on the ZP, SCD, and EDL properties of macroions, and consequently, their stability and aggregation. Furthermore, simulating different electrolyte solutions may be useful to understand the effects of biological fluids on biomolecules and their functions. It is also possible to perform multiple simulations for different macroion sizes, charges and shapes to study their influence on membrane absorption.

V. FUTURE DIRECTIONS
Users will be notified to update the software when a new release is available. Additional features to be included in coming releases are: 1) modeling more complex charging mech- anism on the surface of nanomaterials including more than one active functional group and the release /adsorption of ions; 2) modeling solvent polarization effects [24]. Future work will be focused on modeling: 1) inhomogenous surface charge density and solutes with irregular shapes; 2) nanomaterials having surface functionalization; 3) ionic currents along biopolymers [13].

ACKNOWLEDGMENTS
This work was supported by NIH Grant 1SC2GM112578-03. The author warmly thanks technicians Esteban Valderrama and Samrita Neogi for their participation in programming and testing the software. In this approach, we consider a rigid charged macroion of effective radius R and uniform surface charge density σ surrounded by an electrolyte solution comprised of m ionic species. We use the solvent particle model to characterize the electrolyte. Each ionic species i is represented by bulk Molar concentration [ρ 0 i ], a charged hard sphere of diameter d i , and total charge q i = ez i , where e is the electron charge and z i is the corresponding ionic valence.
Additionally, the solvent molecules are represented as a neutral ion species whereas the solvent electrostatics is considered implicitly by using the continuum dielectric environment with a dielectric constant . The macroion-liquid interaction induces inhomogeneous ion profiles [ρ i (r)] which are calculated using CSDFT as follows [12,14]: the integrated charge (see Figure 1). Expression (2) is the formal solution of the PB equation for an homogeneous anisotropic dielectric media


Optimal Release Time Estimation of Software System using Box-Cox Transformation and Neural Network

The determination of the software release time for a new software product is the most critical issue for designing and controlling software development processes. This paper presents an innovative technique to predict the optimal software release time using a neural network. In our approach, a three-layer perceptron neural network with multiple outputs is used, where the underlying software fault count data are transformed into the Gaussian data by means of the well-known Box-Cox power transformation. Then the prediction of the optimal software release time, which minimizes the expected software cost, is carried out using the neural network. Numerical examples with four actual software fault count data sets are presented, where we compare our approach with conventional Non-Homogeneous Poisson Process (NHPP) -based Software Reliability Growth Models (SRGMs). KeywordsSoftware cost model, Optimal software release time, Software reliability, Artificial neural network, Data transformation, Long-term prediction, Fault count data, Empirical validation.

Introduction
Software testing is an important feature in software development processes and plays an imperative role in defining the quality of software. Moreover, due to huge competitions in the market, software can be seldom survived during the long lifetime due to the competitions and always faces the version up. At the same time, newly added feature in software increases the complexity, which leads to the software faults if not evaluated appropriately. On the other hand, the main concern in software process management by practitioners is to determine when to stop software testing and release the software to the market or users. The determination of the optimal timing to stop software testing is called the optimal software release problem, which is directly related to software testing costs. The total testing cost is reduced successfully if the volume of software test work is less. The debugging cost is more significant in the operational phase than that in the testing phase after releasing software. On the contrary, the higher the software reliability desired, the longer testing period, which further accumulates higher testing cost. Hence, it is important to find as appropriate software release time taking account of the expected software cost. Many authors have discussed in the past on when to stop software testing and to release it for use. The recent trend shows that the quantitative software reliability assessment processes and the efficient software testing methodologies are becoming a firm prerequisite for software developers. Towards achieving the prerequisite, a great number of SRGMs have been developed in the literature (see Goel and Okumoto,1979;Yamada et al., 1983;Ohba, 1984;Littlewood, 1984;Goel, 1985;Abdel-Ghaly et al., 1986;Lyu, 1996;Achcar et al., 1998;Cai, 1998;Gokhale and Trivedi, 1998;Pham, 2000;Ohishi et al., 2009;. Many researchers have also studied the optimal software release policies as well. Okumoto and Goel (1980) derived the optimal software released time for their NHPP-based software reliability model (SRGM) such that the software reliability attains a certain requirement level. Pham and Zhang (1999) proposed a software cost model with warranty and risk costs. The influential involvement of Leung (1992), Dalal and McIntosh (1994), Zhang and Pham (1998), Pham and Zhang (1999), Pham (2003) has been seen in a number of optimization problems under different modeling assumptions. They focused on the completely known stochastic point process such as NHPP-based SRGM to estimate the optimal software testing time. However, it is obvious that there is a need for a better SRGM to fit every software fault count data, addressing the fact that the formulation of the optimal software release problems strongly depends on the accurate statistical estimation of software fault count process. Kaneishi and Dohi (2013), Xiao and Dohi (2013), and Saito and Dohi (2015) introduced new non-parametric assessment methods for software fault count process in testing. However, the disadvantages of these methods rises on the failure to predict the long-term behavior of software fault count process; therefore, these methods are not applicable to the actual software release problems. Recently, a non-parametric method was employed to predict approximately the optimal software release time by minimizing the upper or lower bound of the expected software cost by Saito et al. (2016).
In this paper, a neural network approach with fault count data is used to estimate the optimal software release time, which reduces the relevant software cost. The well-known data transformation technique called the Box-Cox power transform (see Box and Cox, 1964), is applied to transform the data with an arbitrary probability law into the Gaussian data, and to transform the underlying software fault prediction problem into a nonlinear Gaussian regression problem with the Multilayer Perceptron (MLP). An MLP is a feed forward artificial neural network model where all of the input data are mapped onto a set of appropriate outputs. The MLP exploits a supervised learning technique called back propagation. The back propagation is a common learning technique for training Artificial Neural Networks (ANNs) and is used in conjunction with an optimization method such as gradient descent. The algorithm repeats a two-phase cycle; even propagation and weight update. The problems of ANN exist in many applications including software fault prediction, so that there is no efficient way to determine the best neural network architecture. The number of input neurons and output neurons may be determined from physical requirements in many cases. On the other hand, the number of hidden layers and hidden neurons significantly influence the prediction performance in the MLP feed forward neural networks. Nevertheless, the MLP is widely applicable to the statistical identification, function approximation, and time series forecasting.
The idea on data transformation is rather simple, fit useful to handle the non-Gaussian data. Tukey (1957) introduced a family of power transformations where an arbitrary transformation parameter is involved in the power transformations. Box and Cox (1964) proposed a method to calculate the maximum likelihood as well as the Bayesian methods for estimation of the parameter, and derived the asymptotic distribution of the likelihood ratio to test some hypotheses about the parameter. The main contribution by Box and Cox (1964) are two-folds: the first one is to calculate the profile likelihood function and to acquire an approximate confidence interval from the asymptotic property of the maximum likelihood estimator; the second one is to ensure that the probability model is fully identifiable in the Bayesian approach. Dashora et.al. (2015) used the Box-Cox power transformation for data driven prediction models with single output neuron in MLP and analyzed the intermittent stream flow for Narmada river basin in the context of neural computation. It compared against the MLP with seasonal autoregressive integrated moving average and Thomas-Fiering models. Sennaroglu and Senvar (2015) estimated the process capability indices in industry and compared the Box and Cox power transformation with weighted variance methods in the Weibull distribution model in terms of the process variation on the product specifications. In these ways, the potential applicability of the Box and Cox power transform is acknowledged in several research fields. Begum and Dohi (2016a) proposed a one-stage look-ahead prediction method with an MLP where the software fault count data are used. The method is adequate to predict the cumulative number of software faults in practice. Dohi et al. (1999) considered an optimal software release problem by means of the MLP-based software fault prediction, and was able to estimate the optimal software release time. They proposed an interesting graphical method to find directly an optimal software testing time, which minimizes the relevant software cost. However, the model dealt with only software fault-detection time data, and failed to make the multi-stage look-ahead prediction. An idea to use Multiple-Input Multiple-Output (MIMO) came from Park et al. (2014) who proposed a refined ANN approach to predict the long-term behavior of a software fault count with the grouped data. Dohi (2016b, 2016c) imposed a plausible assumption that the underlying fault count process obeys the Poisson law with an unknown mean value function, and proposed to utilize three data transform methods from the Poisson count data to the Gaussian data; Bartlett transform (1936), Anscombe transform (1948) and Fisz transform (1955). They also proposed the idea to apply an MIMO type of MLP to an optimal software release problem with the grouped data. Recently, Begum and Dohi (2017) proposed another MLP-based software fault prediction method with Box-Cox power transformation.

NHPP-Based Software Reliability Modeling
Assume that software system test starts at time = 0. Let ( ) denote the cumulative number of software faults detected by time t. Then the counting process { ( ), ≥ 0} is said to be a nonhomogeneous Poisson process (NHPP) with intensity function ∅( ; ), if the following conditions hold: where (ℎ) is the higher term of infinitesimal time h, and ∅( ; ) is the intensity function of an NHPP which denotes the instantaneous fault detection rate per each fault. In the above definition, is the model parameter (vector) included in the intensity function. Then, the probability mass function (p.m.f.) of the NHPP is given by The identification problem of the NHPP is reduced to a statistical estimation problem of the unknown model parameter , if the mean value function ( ; ) or the intensity function ∅( ; ) is known. The parametric statistical estimation methods can then be applied when the parametric form of the mean value function or the intensity function is given. The maximum likelihood method and/or the EM (Expectation-Maximization) algorithm are commonly used to estimate the mean value function.  summarized eleven parametric NHPP-based SRGMs (see Table 1) and developed a parameter estimation tool, SRATS. The best SRGM with smallest AIC (Akaike Information Criterion) is automatically selected in SRATS, where the best SRGM fits the past observation data on software fault counts among the eleven models.
Suppose that realizations of ( ), ( = 1,2, … , ), are observed up to the observation point . We estimate the model parameter , by means of the maximum likelihood method. Then, the log likelihood function for the fault count data ( , ) ( = 1,2, … , ) is given by where (0; ) = 0, = 0 and = for simplification. By maximizing Equation (3) with respect to the model parameter , we can predict the future value of the intensity function or the mean value function at an arbitrary time ( = 1,2, … ), where l denotes the prediction length. In parametric modeling, the prediction at time is easily done by substituting estimated model parameter into the time evolution ( ; ), where the unconditional and conditional mean value functions at an arbitrary future time are given by If the parametric form of the mean value function or the intensity function is unknown, the identification problem of an NHPP becomes much more difficult. A limited number of nonparametric approaches have been developed by Kaneishi and Dohi (2013), Saito and Dohi (2015), etc. However, those approaches can deal with the fault-detection time data, but will not work for long-term prediction with the grouped fault count data. In addition, the wavelet-based method in Xiao and Dohi (2013) can handle the grouped data, but is unsuccessful to make the long-term prediction. To overcome the above problem we use a fundamental MLP for long-term software fault prediction.

MLP Architecture
An Artificial Neuron Network (ANN) is a computational model inspired by the structure and functions of biological neural networks. An ANN has several advantages but one of the most significant ones is that it can be trained from past observation. The simplest ANN has three layers that are interconnected. The first layer consists of input neurons. Those neurons send data to the second layer i.e. hidden layer, which sends the outputs to the third layer. Subsequently, the hidden neurons have no communication with the external world, so that the output layer of neurons sends the final output to the external world. The problem is how to get an appropriate number of hidden neurons in the neural computation. Here, we study an MIMO type of MLP with only one hidden layer. Similar to Section 2, suppose that n software fault count data ( , ) ( = 1,2, … , ) are observed at the observation point (= ). Our concern is about the future prediction of the cumulative number of software faults at time ( = 1,2, … ).

First Phase: Data Transformation
The simplest MLP with only one output neuron is considered as a nonlinear regression model, where the explanatory variables are randomized by the Gaussian white noise. In particular, the output data in the MLP are indirectly assumed to be realizations of a nonlinear Gaussian model. By contrast, the fault count data are integer values. Hence, the original data shall be transformed to the Gaussian data in advance. Xiao and Dohi (2013) used a pre-data processing which is common in the wavelet shrinkage estimation, where the underlying data follows the Poisson data.
In the same way, we apply the Box-Cox power transformation technique proposed by Box and Cox (1964), from an arbitrary random data to the Gaussian data. As mentioned in the above, they developed a method to find an appropriate exponent to transform data into a "normal shape". Table 2 presents the Box-Cox power transformation and its inverse transform formula where denotes the cumulative number of software faults detected at (1,2, … , )-th testing day. Then, we have the transformed data by means of the Box-Cox power transformation. The transformation parameter indicates the power to which all data should be raised, where the parameter has to be adjusted in the Box-Cox power transformation. Generally, the transformation parameter should be optimal that's why before time series prediction pre-experiments is essential. Let ( = 1,2, … , ) and ( = 1,2, … ) be the input and output for the MIMO type of MLP, respectively. Then, the prediction of the cumulative number of software faults are given by the inversion of the data transform. Fig 1 depicts the architecture of back propagation type MIMO, where is the number of software fault count data and is the prediction length. We suppose that there is only one hidden layer with (= 1,2, … ) hidden neurons in our MIMO type of MLP.  Vol. 3, No. 2, 177-194, 2018ISSN: 2455

Second Phase: Training the Neural Network
In Fig 1, denotes the number of input neurons in the input layer, of the hidden neurons and l indicates output neurons so we assume that all the connection weights ( weights from input to hidden layer, kl weights from hidden to output layer) are first given by the uniformly distributed pseudo-random varieties. In our MIMO type of MLP, output neuron can be calculated ( , … , ) from the previous transformed input ( , … , )if these weights are completely known. However, in principle of the common BP algorithm, it is impossible to train all the weights including ( + ) unknown patterns, as a result, it is required to improve the common BP algorithm for long-term prediction scheme.

(i) Long-Term Prediction Scheme
In Fig 2, we show the structure of our prediction scheme. For the long-term prediction, we assume that > without any loss of generality. In order to predict the cumulative number of software faults for testing days from the observation point , the prediction has to be made at the point . This indicates that only ( − ) + = weights can be estimated with the training data experienced for the period ( , ] and that the remaining ( + ) − weights are not trained at time . We call these ( + ) − weights the non-estimable weights in this paper. In addition, the prediction improbability also rises more when the prediction length is longer, and the number of non-estimable weights becomes greater. In this scheme, the Box-Cox transformed data ( , … , ) with given are used for the input in the MIMO, and the remaining data ( , … , ) are used for the teaching signals in the training phase.

(ii) Common BP Algorithm
Back Propagation (BP) is a common method for training a neural network. The BP algorithm is the well-known gradient descent method to update the connection weights, to minimize the squared error between the network output values and the teaching signals. There are two special inputs; bias units which always have the unit values. These inputs are used to evaluate the bias to the hidden neurons and output neurons, respectively. For the value coming out an input neuron, ( = 1,2, … , − ), let ∈ [−1, +1] be the connection weights from -th input neuron to -th hidden neuron, where and denote the bias weights for -th hidden neuron and -th output neuron, respectively, for the training phase with = 0,1, … , − , = 0,1, … , and = − + 1, − + 2, … , . Each hidden neuron calculates the weighted sum of the input neuron, ℎ , in the following equation: Since there is no universal method to determine the number of hidden neurons, we change in the pre-experiments and choose an appropriate value. After calculating ℎ for each hidden neuron, we apply a sigmoid function (ℎ ) = 1 exp (−ℎ ) ⁄ as a threshold function in the MIMO. Since ℎ are summative and weighted inputs from respective hidden neurons, the -th output ( = − + 1, − + 2, … , ) in the output layer is given by Because are also summative and weighted inputs from respective hidden neurons in the output layer, the weight is connected from j-th hidden neuron to s-th output neuron. The output value of the network in the training phase, , is calculated by ( ) = 1 exp (− ) ⁄ . In the BP algorithm, the error is propagated from an output layer to a successive hidden layer by updating the weights, where the error function is defined by with the prediction value and the teaching signal observed for the period ( , ]. Next we overview the BP algorithm. It updates the weight parameters so as to minimize SSE between the network output values ( = − + 1, − + 2, … , ) and the teaching signals , where each connection weight is adjusted using the gradient descent algorithm according to the contribution to SSE in Equation (8). The momentum, , and the learning rate, , are controlled to adjust the weights and the convergence speed in the BP algorithm, respectively. Since these are the most important tuning parameters in the BP algorithm, we carefully examine these parameters in pre-experiments. In this paper we set =0.25~0.90 and = 0.001~0.500. Then, the connection weights are updated in the following: where ℎ and are the output gradient of j-th hidden neuron and the output gradient in the output layer, and are defined by respectively. Also, the updated bias weights for hidden and output neurons are respectively given by The above procedure is repeated until the desired output is achieved.

Last Phase: Long-Term Prediction
Once the weights are estimated with the training data experienced for the period ( , ], we need to obtain the remaining ( + ) − non-estimable weights for prediction through the BP algorithm. Unfortunately, since these cannot be trained with the information at time , we need to give these values by the uniform pseudo random variates in the range [−1,1]. By giving the random connection weights, the output as the prediction of the cumulative number of software faults, ( , … , ), are calculated by replacing Equations (6) and (7) by respectively, for = 1,2, … , , = 1, … , and = 1,2, … , . Note that the resulting output is based on one sample by generating a set of uniform pseudo random variates. In order to obtain the prediction of the expected cumulative number of software faults, we generate m sets of random variaties and take the arithmetic mean of m predictions of ( , … , ), where = 1,000 is confirmed to be enough in our preliminary experiments. In other words, the prediction in the MIMO type of MLP is reduced to a combination of the BP learning and a Monte Carlo simulation on the connection weights.

Optimal Software Release Decision
Assume that the system test of a software product starts at = 0 and terminates at = . Let be the software lifetime or the upper limit of the software warranty period, where the time length ( , ] denotes the operational period of software after the release. When the software fault count data = { , … , } which are the cumulative number of detected faults at time ( = 1,2, … , ) are observed at time (0 < ≤ ), the model parameter is estimated with these data in parametric NHPP-based SRGMs. Define the following cost components: . (17) Note that Λ ; is independent of and that Λ ; = from Eq.(5). When = , l corresponds to the remaining testing length. Hence, our optimization problem is essentially reduced to When SRGMs are assumed, the mean value function Λ ; can be estimated from the underlying data. On the other hand, in the context of neural computation, the output of MLP in Eq. (16) is regarded as an estimate of Λ ; . Hence, output sequence( , , … , ) denotes Λ ; , Λ ; , … , Λ ; in the MIMO type MLP. Dohi et al. (1999) gave the essentially similar but somewhat different dual problem to Eq. (18). Consequently, the optimal release time ̂ is equivalent to the point with the maximum vertical distance between the predicted curve Λ ; and a straight line ( − ) ⁄ in the two-dimensional plane. Then, it can be graphically checked that there exists a finite optimal software release time. In fact, when ̂ is equal to the observation point , then it is optimal not to continue testing the software product any more. In other words, the optimization problem considered here is a prediction problem of the function Λ ; .

Numerical Experiments
We give numerical examples to predict the optimal software release time based on our neural network approach, where four data sets, DS1∼DS4, are used for analysis (see Lyu (1996)); which consist of the software fault count (grouped) data. In these data sets, the length of software testing and the total number of detected software faults are given by (62, 133), (41, 351), (46,266) and (109,535) respectively. To find out the desired output via the BP algorithm, we need much computation cost to calculate the gradient descent. The initial guess of weights, , , and , are given by the uniform random variates ranged in [−1, +1], The number of total iterations in the BP algorithm run is 1,000 and the convergence criterion on the minimum error is 0.001 which is same as our previous paper by Begum and Dohi (2016a). In our experiments, it is shown that the search range of the transformation parameter λ should be [−3, +2]. We define the error criterion in the following to evaluate the prediction performance of the optimal software release time via MIMO type of MLP or SRGM, where ̂ is a posteriori optimal testing time and is calculated by finding the point (= 0,1,2, … , ) which maximizes Λ ; − ( − ) ⁄ , and q satisfies = . This error criterion is used for the predictive performance when all the data , , … , , , , … , are given. In Table 3, we present posteriori optimal release times for DS1∼DS4, when = 4, = 1 = 1 and varies in the range of 2∼10.
Meanwhile the cost parameter is not sensitive to the posteriori optimal release time in Table 3, we focus on only two cases; = 2 and = 10 hereafter. In Tables 4 and 5 we give the prediction results on the optimal software release time for DS1 with = 2 and = 10, respectively, where the prediction of optimal release time and its associated prediction error are International Journal of Mathematical, Engineering andManagement Sciences Vol. 3, No. 2, 177-194, 2018 ISSN: 2455-7749 187 calculated at each observation point (50%∼90% points of the whole data). In these tables, the bold number implies the best prediction model by MIMO type of MLP in comparison with the best SRGM among eleven models in Table 1. In the MIMO type of MLP, we compare Box-Cox power transformation results with the non-transformed case (Normal) and the best SRGM. In the column of SRGM, we denote the best SRGMs in terms of prediction performance (in the sense of minimum average relative error; AE) and estimation performance (in the sense of minimum Akaike Information Criterion; AIC), denoted by P and E, respectively, where the capability of the prediction model is measured by the Average Error (AE); where is called the relative error for the future time t = n + s and is given by So we regard the prediction model with smaller AE as a better prediction model. In Table 4 and 5, it is seen that our approaches could provide smaller prediction error than the common SRGM in almost all observation points (50%∼90%).  Table 4. Prediction results of optimal software release time with DS1 with = 2 Table 5. Prediction results of optimal software release time with DS1 with = 10 International Journal of Mathematical, Engineering andManagement Sciences Vol. 3, No. 2, 177-194, 2018 ISSN: 2455-7749 189 In contrast, at 80% observation time the non-transformed case (Normal) showed small errors (see Table 5). Even in these cases, it should be noted that the best SRGM with the minimum AIC is not always equivalent to the best SRGM with the minimum AE. This fact tells that one cannot know exactly the best SRGM in terms of judgment on when to stop software testing.
Tables 6 and 7 give the optimal software release time for = 2 and = 10 with DS2. In the latter phase of software testing, i.e., 80%∼ 90% observation points, SRGMs, such as txvmin, txvmax and lxvmin, could offer less error than our MIMO type of MLP (see Table 6). Similarly Tables 8∼11 provide the optimal software release time for = 2 and = 10 with DS3 and DS4 respectively. In most of the cases our MIMO type of MLP could give the best timing of software release for the all testing phase. Throughout our numerical experiments, we can conclude that when the cost parameter increases, the resulting software release time also increases. Table 6. Prediction results of optimal software release time with DS2 with = 2 International Journal of Mathematical, Engineering andManagement Sciences Vol. 3, No. 2, 177-194, 2018 ISSN: 2455-7749 190 Table 7. Prediction results of optimal software release time with DS2 with = 10 Table 8. Prediction results of optimal software release time with DS3 with = 2 International Journal of Mathematical, Engineering andManagement Sciences Vol. 3, No. 2, 177-194, 2018 ISSN: 2455-7749 191 Table 9. Prediction results of optimal software release time with DS3 with = 10 Table 10. Prediction results of optimal software release time with DS4 with = 2 International Journal of Mathematical, Engineering andManagement Sciences Vol. 3, No. 2, 177-194, 2018 ISSN: 2455-7749 192 Table 11. Prediction results of optimal software release time with DS4 with = 10

Conclusion
Almost all software is complex in nature, which has a huge testing scope. All defects in the software are possible to find out but testing will be endless and testing cycles will continue until a decision is made when to stop. Now it becomes even more complicated to come to a decision to stop testing which is being minimized the relevant expected cost. We have given illustrative examples with four real software fault data to derive estimates of the optimal software release time and performed a sensitivity analysis for each the observation point. In the numerical examples are that our MIMO type of MLP could give the accurate optimal software release time in the all testing phases. Additionally, the method can also help project managers to decide when to stop the software testing for market release at the earlier time. Future work will include to develop a tool to support the optimal software testing decision by automating the BP learning and the size determination of a hidden layer.