extfieldsofstudy:
Physics
Computer Science

text:
Paraiso : An Automated Tuning Framework for Explicit Solvers of Partial Differential Equations

We propose Paraiso, a domain specific language embedded in functional programming language Haskell, for automated tuning of explicit solvers of partial differential equations (PDEs) on GPUs as well as multicore CPUs. In Paraiso, one can describe PDE solving algorithms succinctly using tensor equations notation. Hydrodynamic properties, interpolation methods and other building blocks are described in abstract, modular, re-usable and combinable forms, which lets us generate versatile solvers from little set of Paraiso source codes. We demonstrate Paraiso by implementing a compressive hydrodynamics solver. A single source code less than 500 lines can be used to generate solvers of arbitrary dimensions, for both multicore CPUs and GPUs. We demonstrate both manual annotation based tuning and evolutionary computing based automated tuning of the program.

Introduction
Today, computer architectures are becoming more and more complex, making it more and more difficult to predict the performance of a program before running it. Parallel architectures, GPUs being one example and distributed memory machines another, forces us to program in different ways and our programs tend to become longer. Moreover, to optimize them we are asked to write many different implementations of such programs, and the coding task becomes time-consuming. However, if we describe the problem we want to solve in a domain-specific language (DSL) from which a lot of possible implementations are generated and benchmarked, we can automate the tuning processes.
Examples of such automated-tuning DSL approaches are found in fast Fourier transformation library FFTW [1], linear algebra library ATLAS [2], digital signal processing library SPIRAL [3]. In these works the authors regarded automated tuning not just as a tool to avoid manual tuning, but as a necessary tool to have "portable performance" -the practical way of optimizing domain-specific codes for various complicated architecture we have today and in the future. Now, our domain is explicit solvers of partial differential equations (PDEs). The main portions of such solvers consist of stencil computations, and a few global reduction operations needed to calculate Courant-Friedrichs-Lewy (CFL) conditions. Stencil computations are computations that update mesh structures, and the next state of a mesh depends only on the states of its neighbor meshes. Due to its inherent and coherent parallelism, stencil computation DSLs have been actively studied. For example, [4] have demonstrated ×1.5 − ×5.6 speedup on various hardwares via hardware and memoryhierarchy aware automated tuning. [5] have demonstrated generating multi-GPU codes that weakly scales up to 256 GPUs.
Still, reports are limited to simple equations such as diffusion equations and Jacobi solvers of Poisson equations, and to implement solvers of more complicated equations such as compressive hydrodynamics, magneto-hydrodynamics or general relativity, and their higher-order versions, we are forced to manually decompose the solving algorithms to imperative instructions that are often tens and hundred thousands of lines.
It is a problem common to all the languages that is designed to be compatible with C or Fortran, such as OpenMP, CUDA, PGI Fortran and OpenACC -that it is hard, if not impossible at all, to avoid repeating yourself in these languages. Despite all the efforts made so far to extend these languages, the programs in those languages fundamentally lack the ability to manipulate programs themselves and other abstract concepts. Every time a new generation of language appear, we are forced to make painful choices of porting a huge amounts of legacy codes to the new language. It is also a pain that various numerical techniques are mixed in one code, and can hardly be reused. If we want to make computers automatically combine such numerical techniques, compose a variety of implementations, and search for the better ones, the abstraction power is necessary. A complementary approach is needed here, that works together with the parallel languages.
Our contribution, Paraiso (PARallel Automated Integration Scheme Organizer ), enables us describe PDE-solving algorithms succinctly using algebraic concepts such as tensors. Paraiso also enables us to describe various manipulations on algorithms such as introducing higher-order interpolations, in composable and reusable manner. That is, once we define a certain interpolation method in our language as a transformation of basic solver to a higher-order one, we can reuse the transformation for any equations. In this way Paraiso reduces the cost of rewriting when we search for better discretizations or interpolations; enables us to port the collections of algorithms to new parallel languages, without changing the Paraiso source codes but by updating the common code generator; allows yet another layer of automated tuning at translating Paraiso codes to other parallel languages.
In writing Paraiso, we wanted to define arithmetic operations between tensors and code generator fragments. We want our tensors to be polymorphic, but we don't want to allow addition between tensors of different dimensions. We wanted to make generalized functional applications. We wanted to traverse over data structures. We need to manage various contexts, like context of code generation and the context of serialization from a program to a genome.
Programming language Haskell supports all of these. Moreover, they are supported not as built-in features that users cannot change; they are libraries that Haskell users can freely combine or create their own. This flexibility of Haskell is based on its strong, static, higher-order type system with type classes and type inference. Thus Haskell essentially allows us to develop our own type-systems within it, which gives it a unique advantage as a platform for developing embedded DSLs. Many parallel and distributed programming languages has been implemented using Haskell [6]. Nepal [7] and Data Parallel Haskell [8] are implementations of NESL, a language for operating nested arrays in Haskell. Accelerate [9] and Nikola [10] are languages to manipulate arrays on GPUs written in Haskell. Finally, Liszt [11] is a DSL for solving mesh-based PDEs based on functional programming language Scala.
Our contributions are the following: • A domain-specific language (DSL) embedded in Haskell, with which one can describe explicit solvers of partial differential equations (PDEs) in a succinct and organized manner, using tensor notations.
• A code-generation mechanism that takes the DSL and generates OpenMP and CUDA programs.
• A compressible Euler equations solver implemented in the DSL, and tuning experiments using the solver. Written in tensor notations, the solver can be applied to problems of arbitrary dimension without changing the source code at all but a single type declaration that sets the dimension of the solver.
• An annotation mechanism with which one can give hints to code generators, which makes it drastically easy to search for better implementations manually. By adding just one line of hint to the solver causes overall refactoring, adding 4 subroutines to the generated code, making it consume 1.3 times more memory but 6.42 times faster.
• An automated benchmarking and tuning mechanism based on parallel simulated annealing and genetic algorithms, which generates a lot of different implementations of the PDE solvers and search for faster implementations. It speeds up the unannotated solver by factor of 5.26, and the annotated solver further by factor of 1.78.
Just for a comparison, the Paraiso framework is about 5'000 lines of code in Haskell, and the compressible Euler equations solver implemented in Paraiso is less than 500 lines. From that, we have generated more than 500'000 instances of the solver, each being 3'000 -10'000 lines of code in CUDA. The automatically tuned codes are faster than the manually tuned codes reported by other groups. Paraiso is ready to optimize solvers of other equations, and all the solvers written in Paraiso can migrate to new parallel languages and new hardwares once the Paraiso code generator supports them.
Our work, Paraiso is the first consistent system that combine those previously studied techniques of symbolic computations, DSLs, GPU computations and automated tuning. We demonstrate the utility of such a system in the domain of explicit solvers of PDEs.
This paper is organized as follows. In section 2, we describe the overall design of Paraiso as well as its components. In section 3, we describe our automated tuning mechanism, which is a combination of genetic algorithm and simulated annealing that is designed to utilize the varying number of available nodes in a shared computer system. In section 4, we introduce the compressible hydrodynamics solver we choose as the tuning target, and describe the manual and automated tuning experiments. In section 5, we analyze the automated tuning history. In section 6 are concluding remarks and discussions.

Overall Design
Paraiso is to tackle the ambitious problem of generating fast and massively-parallel codes from human-friendly notations of algorithms. To divide the problem into major components, which may be conquered by different set of people, we set the overall design of Paraiso, as illustrated in Fig. 1.
The users of Paraiso must manually invent a discretized algorithm for the partial differential equations they want to solve. The solving algorithms are described in Paraiso using Builder Monads. Builder Monads generate programs for Orthotope Machine (OM), a virtual parallel machine designed to denote parallel computations on multidimensional arrays. Then the back end OM compiler generates the native machine source codes such as C++ and CUDA. Finally, native compilers translate them to executable codes.
When someone translates a simulation algorithm in his mind, to different languages such as Fortran, C and CUDA, he starts from building mathematical notations in his mind, and then gradually decomposes it to machine level. There should be the last common level that is independent of the detail of the target languages or the target hardwares. The level consist of very primitive operations on arrays. Orthotope Machine (OM) is designed to capture this level: The instruction set of the OM is kept as compact as possible, while maintaining all the parallelisms found in the solving algorithms. In forthcoming technical report we plan to provide the formal definition of the Orthotope Machine.
In order to generate the OM data-flow graphs from tensor expressions, and to translate them as native programs and further to apply manual and automated tuning over them, Paraiso introduces a number of abstract concepts centered around the OM. These components are quite orthogonal to each other, and some may even be useful outside the context of Paraiso. The data-flow graph for the OM. c.f. §2.2, Fig. 2, Fig. 3.
type Annotation in Language.Paraiso.Annotation The collection of annotations that are added to each OM data-flow graph node. c.f. §2.6, §3.1.
data Plan in Language.Paraiso.Generator.Plan The fixed detail of the code to be generated such as amount of memories and what to do in each subroutine. c.f. Fig. 5. To see how a plan is fixed c.f. Fig. 7.
data Program in Language.Paraiso.Generator.Claris The subset of C++ and CUDA syntax tree which is sufficient in generating codes in scope of Paraiso. c.f. Fig. 5.
newtype Genome in Language.Paraiso.Tuning.Genetic The set of annotations that belongs to an individual encoded as a string of letters, with which one can mutate, cross and triangulate. The evolution algorithm is in §3.2. interface, we can combine various stencil computation applications with state-of-theart techniques developed so far, such as cache-friendly data structures [4], overlapping communication with computation [12,5] or heterogeneous utilization of CPU/GPU [13]. On the other hand, various concepts of numerical simulations has been decomposed to elemental calculations before the OM level, so the problem of how to build a parallel computation from components such as new spatial interpolations, time marching methods and approximate Riemann solvers, can be addressed and developed separately from the detail of the hardware. To achieve such orthogonality is one of the aims of the Paraiso project.

Outline of The Orthotope Machine(OM)
The Orthotope Machine (OM) is a virtual machine much like vector computers. Each register of OM is multidimensional array of infinite size. Arithmetic operations of OM work in parallel on each mesh, or loads from neighbor cells. We have no intention of building a real hardware: OM is a thought object to capture parallel algorithms to data-flow graphs without losing parallelism.
The instruction set of OM resembles those of historical parallel machines such as PAX computer [14], and is subset of partitioned global address space (PGAS) languages such as XcalableMP [15].
Each instance of OM have a specific dimension (e.g. a two-dimensional OM of size N 0 × N 1). The variables of OM are either arrays of that dimension, or a scalar value. All of the arrays must have a common size. The actual numbers (N 0, N 1) are fixed at native code generation phase. We say that arrays are variable with Local Realm , while the scalars are variable with Global Realm.
OM has a set of Static variables which denotes the current state of the simulation. Each Static variable has a string name. OM has a set of Kernels -they are subroutines for updating the Static variables. Inside each Kernel, you can generate Temporal values in static single assignment (SSA) manner.
To summarize, the lifetime of an OM variable is either Static or Temporal, and the Realm of an OM variable is either Local or Global. Static variables survive multiple Kernel calls, while Temporal variables are limited to one Kernel. Local variables are arrays. Global variables are scalar values; in other words, they are arrays whose elements are globally the same.
OM cannot handle array of structures; Local variables may contain only simple objects such as Bool or Double, but not composite ones such as complex numbers or vectors. Nevertheless we can easily use such composite concepts in Paraiso programs at Builder Monad level and translate them to OM instructions, for example by using the applicative programming style [16].
An OM Kernel is a directed bipartite graph consisting of NInst nodes and NValue nodes, as illustrated in Fig. 2. Each node has arity :: a -> (Int, Int), number of incoming and outgoing edges. NValue nodes have arity of (1,n), where n is the number of NInst nodes that uses the value and can be an arbitrary integer. Arities of NInst nodes are inherited from the instructions they carry.
OM has nine instructions: Here, we limit ourselves to C pseudo-code of OM instruction semantics. The formal definition of the OM will be in forthcoming technical reports. Note that we do not translate the OM instructions one by one to C codes listed here. Instead, the instructions are merged into sub-graphs and then translated to C loops (or CUDA kernels) with much larger bodies, to make efficient use of computation resource and memory bandwidth (c.f. §3.1).
Imm arity (0,  Arith perform various arithmetic operations. The arity of this NInst node is inherited from its operator. The realms of the inputs and outputs must match, and can be either of Global or Local. If Local, array elements at matching index are operated in parallel (i.e. zipWith).
For example,

An Example of Orthotope Machine Data-Flow Graph
Here we show the use case of the Orthotope Machine operators introduced in §2.2 within a simple PDE solver. We solve the following linear wave equation: where t is time, x is one-dimensional space coordinate and c is the signal speed. By introducing the time derivative g = ∂f /∂t, Eq. (1) is rewritten as following system of first-order PDEs: and satisfies the following conservation law: We will write a Paraiso code that simulate the discrete system of Eqs. (2,3) and tests the conservation law. Let f n [i] denote the discrete field where n and i are the discrete time and space coordinate, respectively. We choose the following 2nd-order, Lax-Wendorf and Leap Frog scheme to solve Eqs. (2,3) : and measure the following discrete form of the conserved quantity: Table 2 shows the Paraiso implementation for the algorithm Eq. (6,7). The corresponding OM data-flow graph is visualized in Fig. 3. See how the data-flows  Table 2. For brevity, NValue nodes are represented as little circles and only NInst nodes are shown; the Arith instruction nodes are replaced by arithmetic operators inside them, and the instruction arguments are not shown but for Load, Store and Reduce. proceed :: Builder Vec1 Int Annotation () proceed = do c <-bind $ imm 3.43 n <-bind $ loadSize TLocal (0::Double) $ Axis 0 f0 <-bind $ load TLocal (0::Double) $ fieldF g0 <-bind $ load TLocal (0::Double) $ fieldG dx <-bind $ 2 * pi / n dt <-bind $ dx / c f1 <-bind $ f0 + dt * g0 fR <-bind $ shift (Vec:~-1) f1 fL <-bind $ shift (Vec:~1) f1 g1 <-bind $ g0 + dt * c^2 / dx^2 * (fL + fR -2 * f1) store fieldF f1 store fieldG g1 dfdx <-bind $ (fR -fL) / (2*dx) store energy $ reduce Reduce.Sum $ 0.5 * (c^2 * dfdx^2 + ((g0+g1)/2)^2) * dx from Load node to Store node. Here, we have c = 3.43 and 0 ≤ x < 2π. We have performed numerical simulations using this program and the following initial condition: g(x) = cos(3x), with resolution varying from N = 8 to N = 3072, and have confirmed that the fluctuation of discrete conserved quantity, Eq. (8) was smaller than of the order of 10 −13 .

Typelevel Tensor
We introduce typelevel-tensor library, to abstract over the dimensions; the use of tensor notations allow us to describe the algorithms for different dimensions in a single source code. The benefit of having information on tensor dimensions and ranks at the type-level is that we can detect erroneous operations like adding two tensors of different dimensions at compile time. Also we have much less need for explicitly mentioning the tensor dimensions in our programs, because the dimensions will be type-inferred. Our approach is similar to that in [17], where two constructors Z and :. are used to inductively define multi-dimensional tuples. Their n-dimensional vectors are actually n-tuples, so it is possible that the set of n elements consist of different types. In contrast, we needed that all the n elements of a vector are of the same type, because we wanted to operate on them, for example by applying a same function to all of them.
So instead of [17]'s approach: infixl 3 :. data Z = Z data tail :. head = tail :. head We have these: Here, Vec is the type constructor for 0-dimensional vector, and :~is a type-level function that takes the type constructor for n-dimensional vector as an argument and returns the type constructor for n + 1-dimensional vector. We define type synonyms for successively higher vector: Vec0, Vec1, Vec2, and Vec3 are all types of kind * -> *, meaning that it takes one type (the element type) and returns another (the vector type). For example, the three-dimensional double precision vector type is Vec3 Double. Higher-rank tensors are defined as nested vectors; for example Vec3 (Vec3 Int) is a 3 × 3 matrix of integers.
Since we know that all the elements of our tensor are of the same type a, we can make our tensors instances of Traversable type class [16,18]: The benefit of making our tensors instances of Traversable is that we can traverse on them: traverse :: Here, suppose t is our tensor type-constructor and f is some context -for example, a code generation context. a and b are elements of our tensor. Then the type of the traverse function means that if we have code generators for the computation of one element (a -> f b), and we have a t a, a tensor whose elements are of type a, then we can deduce the code generator for computation of the entire tensor f (t b).

Builder Monad
Builder monad is a State monad whose state is the half-built data-flow graph of the Orthotope Machine. To represent the data-flow graph, we use Functional Graph Library (FGL) [19,20]. Authors learned from the Q monad [21] how to encapsulate the construction process. The graph carried by the State monad has the following type: Graph takes three type arguments. vector :: *->* is a type constructor that denotes the dimension of the OM. gauge is the type for the indices of the arrays, which is usually an Int. anot is the type the nodes of the graph are annotated with. Such annotations are used to analyze and optimize the data-flow graph. The three types vector, gauge, and anot are passed to the Nodes of the graph. Nodes are either NValue or NInst. Two types vector and gauge are further passed to the instruction type Inst, because instructions such as shift requires the information on the array dimension and indices. Every graph nodes are also annotated by type anot.
On the other hand, the edges contain none of the three types. They are just unordered edges EUnord or edges ordered by an integer EOrd Int. For example, we can exchange two edges going into addition instruction, but cannot exchange those into subtraction.
We define various mathematical operations between Builder Monad in a consistent manner (c.f. Fig. 4). For any operator ⊕, where Value i is the value computed by Program i which is generated by Builder i. For example, a helper function that takes an operator symbol op, two builders builder1 and builder2, and create a binary operator for builder, is as follows: mkOp2 :: (TRealm r, Typeable c) => A.Operator --^The operator symbol -> (Builder v g a (Value r c)) --^Input 1 -> (Builder v g a (Value r c)) --^Input 2 -> (Builder v g a (Value r c)) --^Output mkOp2 op builder1 builder2 = do We first extract the graph nodes from left hand side and right hand side builders in Builder context, add an NInst node that contains op symbol, add an NValue node after that, and return the node index.
In Haskell, defining mathematical operators between a data type is done by declaring the data type as an instance of the type class that manages the operator. In this sense type classes in Haskell are the parallels of algebraic structures such as group, ring and field, that manages addition, multiplication, and division, respectively. We use a Haskell package numeric-prelude that provides such algebraic structures [22].
For example, an Additive instance declaration of Builder is as follows: instance (TRealm r, Typeable c, Additive.C c) => Additive.C (Builder v g a (Value r c)) where zero = return $ FromImm unitTRealm Additive.zero These type class instances, together with typelevel-tensor library, allows us to write tensor equations used in Paraiso application programs. Moreover, such equations are good for arbitrary instances of the type class. For example, here are the definition of momentum and momentum flux in our Euler equations solver: These functions can be used to directly calculate momentum vector and momentum flux tensor whose components are of type Double. The very same functions are used to generate the solvers for CPUs and GPUs. At that time, their components are inferred to be Builder types. In addition to that, these functions can handle tensors of arbitrary dimensions.

Backend
The back end converts the data-flow graph of Kernels native codes, and create a C++ class corresponding to an OM. The code generation processes of Paraiso is illustrated in Fig. 5. First, various analysis and optimizations are applied. Analysis and optimizations are functions that takes an OM and returns an OM, so we can combine them in arbitrary ways. Though some analysis are mandatory for code generation.
Paraiso has an omnibus interface for analysis and optimization using dynamic programming library Data.Dynamic in Haskell [23] : Here, analyzers as well as human beings can add annotations of arbitrary type a to the graph. On the other hand, optimizers can read out annotations of what type they recognize and perform transformations on the data-flow graph. The set of annotations can be serialized to, and deserialized from genomes, which are binary strings used in automated tuning phase.
Examples of annotations are the choices of whether to store a value on memory and reuse or not to store and recompute it as is needed (c.f. Fig. 6); the boundary analysis result used for automatically adding ghost cells; dependency analysis; and labels used for dead code eliminations.
Once the analysis and optimizations are done, an OM is translated to a code generation Plan. Here, decisions are made for how many memory are used, what portion of computation goes into a same subroutine, and so on. The nodes in data-flow graph are greedily merged as long as they have no dependence and can be calculated in the (a) Some nodes need to be Manifest.
(b) For each of other nodes, the genome specifies whether it is Manifest or Delayed.
(c) The sets of nodes that can be calculated in the same loop are greedily merged into write groups.
(d) The Kernel is divided to multiple subKernels, each of which representing one write group. Each subKernel is in turn translated to a function in C or a global function in CUDA. same loop (Fig. 7). The Plan is further translated to CLARIS (C++-Like Abstract Representation of Intermediate Syntax). CLARIS is subset of C++ and CUDA syntax which is sufficient in generating codes in scope of Paraiso.
Finally, CLARIS is translated to native C++ or CUDA codes.

Tuning Targets
The objects we want to optimize are the implementations of a partial differential equations solving algorithm. We call them individuals, adopting the genetic algorithm terminology. Each individual has a genome that encodes how it have chosen to implement the algorithm. The choices are (1) how much computation speed and memory bandwidth to use, (2) where to synchronize the computation, and (3) the CUDA kernel execution configuration. The fitness of the genome is the benchmark score of the generated code, measured in cups (the number of fluid cells updated per second) which we want to maximize. A simple example program Fig. 6 indicates that we have implementation choices for intermediate variables: whether to store its entire contents on the memory (Manifest) or not to store them and recompute them as they're needed (Delayed). The terms Manifest and Delayed are inherited from REPA [17]. If we increase the number of Manifest nodes, we consume less arithmetic units but more memory and its bandwidth; decreasing the number of Manifest nodes have the opposite effect. There are two extreme configurations, one is making as many nodes Manifest as possible and the other is making as few nodes Manifest as possible. In most cases both of them result in poor performance and moderate configurations are faster. Paraiso generate codes for many possible combinations of Manifest / Delayed choices (Fig. 7) and searches for such configurations by automated tuning.
Another tuning done by Paraiso is the choice of synchronization points. In CUDA, inserting syncthreads(), especially before load/store instructions cause the next instruction to coalesce and increase the speed of the program. Inserting too much synchronization, on the other hand, is a waste of time. Again, Paraiso searches for better configurations by automated tuning.
The last tuning done by Paraiso is to find the optimal CUDA kernel execution configuration, i.e. how many CUDA threads and thread blocks to be launched simultaneously. These tuning items summed, the genome size is approximately 6000 bits for our hydrodynamics solver, which means that there are 2 6000 possible implementations. Brute-force searching for the fastest implementation from this space is inutile and in the first place impossible. Instead, the goal of Paraiso is to stochastically solve such a global optimization problem in this genome space, where the function to be optimized is the benchmark score of the PDE solver generated from the genome. proceed = do x <-bind $ loadLD density y <-bind $ x * x z <-bind $ y + y store density z proceed = do x <-bind $ Anot.add Sync.Pre <?> loadLD density y <-bind $ Anot.add Alloc.Manifest <?> x * x z <-bind $ y + y store density z The generated .cu program file, abbreviated:  Table 3. An example of changes in the genome and the implemented algorithm caused by annotations. Table 3 shows a simple Paraiso kernel and how its genome and implementation is altered by adding annotations. This kernel named proceed performs the following calculations: It first loads from an array named density, performs a multiplication, then an addition, and then stores the result to density again. The header file on the right side of the Table 3, compared to the left one, allocates an additional device vector and declares an additional subkernel as results of a Manifest annotation. The right .cu file, compared to the left one, differs in two points: one is that it performs the multiplication y = x * x and the addition z = y + y in separate CUDA kernels and stores the intermediate result to the additional device vector, which is yet another result of the Manifest annotation. Another difference is that it calls syncthreads() just before the load instruction, which is the result of the synchronization annotation.

Parallel Asynchronous Genetic Simulated Annealing
Simulated annealing has been widely used to solve global optimization problems. However, a standard simulated annealing method is not suitable for parallelization, and if the annealing schedule (how we cool down the temperature as function of time) is too quick, it tends to fall into local minima. Replica-exchange Monte Carlo method [24] solves these drawbacks by introducing multiple replicas of heat baths with different temperatures, and by allowing replica exchange between adjacent heat baths. Thus, replicas can be computed in parallel, and the annealing schedule is spontaneously managed by replica migration.
We further extend this method to fit into benchmark based tuning on shared cluster computer systems. First, in shared systems it is hard to maintain a fixed number of replicas because the available amount of nodes changes with time. Second, it is not efficient to perform replica-exchange synchronously, since the wall clock time required to calculate the fitness function varies with replicas. For these two reasons, we modify replica-exchange Monte Carlo method to a master/worker model (c.f. Fig. 8), where a master asynchronously launches varying number of workers in parallel.
The master holds the genomes and scores of all the past individuals in a database (DB). The role of the master is to draw individuals from the DB, to create new individuals using their genomes, and to launch the workers when the computer resource is available. The role of the workers, on the other hand, is to generate codes from the given genome, to take the benchmark, and to write the results into the DB.
In this way we can launch any number of workers asynchronously as long as the DB is not a bottleneck. In addition to that, we eliminate a common weak point of simulated annealing and genetic algorithms -that individuals of older generation are overwritten and are inaccessible.
For each individual I, the generated code is benchmarked 30 times. We record the mean µ(I) and the sample standard deviation σ(I) of the score. It is important to record the deviation. When we benchmark each individual only once, some individuals receive over-evaluated score by chance, infest the system and tend to stall the evolution.
At the end of each benchmark, the individual is briefly tested if the state of the simulation has developed substantially from the initial condition and there is no NaN (not a number). The test is chosen because not evolving at all and generating NaN are two dominant modes of failure, and the individual test time needs to be kept smaller than benchmark itself. If an individual fails a test, its score is 0. At the end of a tuning experiment, the champion individual is extensively tested if it can reproduce various analytic solutions, which takes hours.
A draw is the operation to randomly choose an individual from the DB. Each draw has a temperature T . In a draw of temperature T , the probability the individual I is where I is the individual with the largest µ(I).
Low-temperature draws strictly prefer high-score individuals while hightemperature draws does not care the score too much. And the difference small compared to σ(I ) + σ(I) is gracefully ignored.
Every time master creates a new individual, it chooses the draw temperature T randomly, so that the probability density of log T is uniformly distributed between log(σ(I )) and log(µ(I )).
To summarize, our method is a master/worker variant of replica-exchange Monte Carlo method [24], using genetic algorithms as neighbor generators. Thus, our method • can utilize parallel and dynamically varying computer resource.
• can find global maximum without hand-adjusted annealing schedule.
• can combine independently-found improvements.
An extensive review on the applications of the evolutionary computation in astronomy and astrophysics is found in [25].

Three Methods of Birth for Generating New Individuals
We use three different methods of birth to generate new individuals (c.f. Fig. 9).  The first method is mutation. We draw one individual from the database (DB), take its genome, overwrite it randomly and create a new genome. The second method is crossover. We draw two individuals from the DB, split their genomes at several random points and exchange the segments. The third method is triangulation, or three-parent crossover.
Multi-parent crossovers are crossovers that takes more than two parents and create one or multiple children. The proposed multi-parent crossover methods include taking bit-wise majority vote [26,27] and exchanging segments between multiple parents [28]. Multi-parent crossovers for real-number coded genomes are also studied [29]. The novelty of our three-parent crossover (triangulation) is to take the scores of the parents into consideration. Triangulation is designed to efficiently combine independent improvement found in sub-problems.
FFTW [1] uses the divide and conquer approach to its target problems. It first recursively decomposes large FFT problem into smaller ones and solves the optimization problem from smaller part. In contrast, Paraiso deals with monolithic data-flow graphs. It is not obvious how to decompose them into subgraphs -decompositions are actually the target of optimization. Therefore, we optimize subgraphs in vivo; we optimize subgraphs keeping them embedded into the entire graph.
Aiming to merge two distinct optimized subgraphs, we draw three individuals from the DB. Then we sort them in ascending order of the score and name them Base, Secondary and Primary. We perform bit-wise operation as in Table 4 to create their child. For each bit, if at least one of Primary or Secondary has been changed from Base, then we adopt the change.
When creating a new individual, the master chooses mutation, crossover, or triangulation at equal probability of 1/3. Then the master draws the needed number of parents from the DB. If two of the parents are the same the master retries the draws. Otherwise, the master creates a child with the chosen method of birth. Then if the genome of the newly created individual is already in the DB, the master discards the individual and retries the mutation until it gets a genome not found in the DB. With each retry, draw temperature T is multiplied by 1.2, so that the master eventually succeed in creating a new genome.

The Target Program
We implemented a 2nd order compressible hydrodynamics solver [30] in Paraiso, and then optimize it. Here we detail numerical scheme. Let where ρ , v a , m a , p and E are the density, the velocity, the momentum, the pressure and the total energy, respectively, The primitive and conserved variables are related by m a = ρv a , E = E i (p) + 1 2 ρv 2 where E i is the internal energy of the gas. By assuming the adiabatic equation of state for perfect gas, E i (p) = (γ − 1) −1 p where γ is the ratio of specific heats.
We numerically solve the Euler equations ∂ t U = ∂ a F a , or more specifically: The numerical method we used is as follows, based on a function Flux a (U0) which calculates the fluxes across the cell surfaces using a Riemann solver. Flux a (U0) is defined as follows. First, V0 are the primitive variables calculated from U0: The interpolated primitive variables VL and VR are where we use piecewise linear interpolation with minbee flux limiter [31,30]: Then, fluxes across the boundaries are defined using the HLLC Riemann solver [32].
This is the flux defined by Flux a (U0). Then, linear addition of flux to conserved variable AddFlux(∆t, F a , U) is defined by where ∆r a is the mesh-size along the a-axis. Using these notations, we construct the second-order time marching as follows, where ∆t is the time step determined by the CFL-condition: U2 is the set of conserved variables for the next generation. Notice how we rely on human mind flexibility when we convey the algorithms in the forms like above. In languages like Fortran or C, we are often forced to decompose those expression into elemental expressions and the source codes tend to become longer. In Haskell, we can express these ideas in well-defined machine readable forms while keeping the compactness and flexibility as in above.
We declare both the primitive variables Eq. (13) and conserved variables Eq. (14) as instances of the type class Hydrable, the set of variables large enough for calculating any other hydrodynamic variables. We used the general form by default, and bind either the primitive or conserved variables when we need the specific form. By doing so requisite minimum number of conversion code between primitive and conserved variables are generated.
When we defined the minbee interpolation for a triplet of real numbers Eq. (21) and then applied it to primitive variables Eq. (20), we implicitly extended a function on real numbers to function on set of real numbers. In Haskell, we can define how any such generalized function application over the set of hydrodynamic variables should behave, just by making it an instance of Applicative type class.  Table 5.
The codes annotated by hand and their performances. First four columns are the ID of the individuals, the CUDA kernel execution configuration, and the two choices of annotation -whether to Manifest or Delay (1) the result of spatial interpolation (2) the result of the Riemann solver. Next three columns are the properties of the generated codes as results of these choices made. They are namely the size of the code (number of lines), the number of subKernels in the code, the memory consumption in proportion to fluid cell number N . The last two columns are the speed of the generated codes for single precision (SP) and double precision (DP). The speed of the codes are measured by Mcups (10 6 cell update per second).
When we calculate the space-derivatives of the fluxes in Eq. (23), the component of the flux we access (index a in F a ) and the direction in which we differentiate (indices a in [i − 1 2 e a ] and [i + 1 2 e a ]) should match. Although the flux F and the spatial array indices i have very different types, Haskell's type inference guarantees that they are both tensors of the same dimensions, and we can access both of them by the common tensor index a. And we can sum over a, because tensors are Foldable and their components are Additive. Tensors are Traversable as mentioned before, and any type constructors that are Traversable are also Foldable.

Annotating By Hand
Before we start the automated tuning experiment, we generate several individuals by hand. The initial individual Izanagi is the one with the least number of Manifest nodes as possible. Izanami is the same individual with the CUDA execution configuration suggested by the CUDA occupancy calculator. Based on it, we add several Manifest annotation and create new individuals (c.f. Table 5). Manual annotations are not blindsearch process; each annotation has clear motivation such as "let us store the result of the Riemann solvers because it is computationally heavy", "let us store the result of interpolations because it moves a lot of data," etc. For example, the sample code in Table 6 shows the implementation of the piecewise linear interpolation with the minbee The interpolation code for Izanami: interpolateSingle order x0 x1 x2 x3 | order == 1 = do return (x1, x2) | order == 2 = do d01 <-bind $ x1-x0 d12 <-bind $ x2-x1 d23 <-bind $ x3-x2 let absmaller a b = select ((a*b) 'le' 0) 0 $ select (abs a 'lt' abs b) a b d1 <-bind $ absmaller d01 d12 d2 <-bind $ absmaller d12 d23 l <-bind $ x1 + d1/2 r <-bind $ x2 -d2/2 return (l,r) For Iwatsuchibiko, the last line is modified as follows: return (Anot.add Alloc.Manifest <?> l, Anot.add Alloc.Manifest <?> r) Table 6. An annotation made in Paraiso source code that causes the results of the interpolations to be stored on memory and reused. flux limiter, Eq. (21) in Paraiso, and how to annotate the return values of the limiter as Manifest.
We use Izanami as the base line of the benchmark, and use Izanagi as the initial individual of some experiment to see if the automated tuning can find out the optimal CUDA execution configuration by itself.
The codes are benchmarked on TSUBAME 2.0 cluster at Tokyo Institute of Technology. Each individual was benchmarked on a TSUBAME node with two Intel Xeon X5670 CPU(2.93GHz, 6 Cores × HT = 12 processors) and three M2050 GPU(1.15GHz 14MP x 32Cores/MP=448Cores).
The abstraction power of Builder Monad lets us change the code drastically with little modification. For example, Paraiso source code of Shinatsuhiko and Izanami differs by just one line of annotation. This introduces the 16 bits of difference in their genome, causing Shinatsuhiko generate a code that has 1/3 lines, which contains four more subroutines, consume 1.31 times more memory and is 6.42 times faster.

Automated Tuning
Next, we performed several automated tuning experiments c.f. Table 7.
To distinguish the contribution of the three methods to generate new individuals, we also performed automated tuning experiments with either crossover or triangulation turned off. Table 8 shows the initial possibility of the master node attempting each RunID prec.
initial score wct best ID/total high score GA-1 DP  Table 7.
The statistics of GB-* experiment series. The columns are RunID, precision, the score of initial individual, the ID of the best individual and the number of individuals generated, the high score (in Mcups). Experiments started with Shinatsuhiko. The nine experiments were performed in parallel, and took 92252 minutes of wall-clock time. method of birth. Note that the actual frequencies of crossover and triangulation are smaller than these value because the master may default to mutation. Table 9 shows the result of experiments. The resolutions were 1024 2 for GA-1 and GA-S1, problems, 100 3 for GA-3D problems and 512 2 for other GA-* series. The resolutions were 512 2 in GB-* series. The automated tuning system can generate and benchmark approximately 10 000 individuals per day. 20 − 100 workers were running at the same time. It takes a few days to tune up Izanami to speed comparable to Shinatsuhiko, or speed up Shinatsuhiko by another factor of 2. The best speed obtained was 35.3Mcups for double precision, and 53.7Mcups for single precision. Our automated tuning experiments on 3D solvers mark 42.4Mcups SP. These are competitive performances to hand-tuned codes for single GPUs; e.g. Schive et. al. [33] reports 30Mcups per C2050 card (single precision, note that their code is 3D. Asuncióna et.al. [34] reports 6.8Mcups per GTX580 card (single precision, 2D). Fig. 10 shows the result of convergence tests for individuals GA-1.20756, GA-4.33991, GA-D.59841, GA-DE.41250, GA-S1.33958, Izanami and Shinatsuhiko. The resolution N is varied from 16 2 to 1024 2 . The codes are tested for entropy wave propagation, sound wave propagation and Sod's shock-tube problem. The detail of test initial conditions are as follows. For all tests, numerically solved domains are (0 < x < 1, 0 < y < 1), out of which the analytic solutions are continuously substituted as boundary conditions. The system was numerically developed until t = 1 for entropy wave and sound wave problems, and t = 0.125 for shock tube problem, and then the numerical solutions were compared to analytic solutions using the norm of L 1 error vector.
The initial conditions are, for entropy wave problem: for sound wave problem: where the amplitude A = 10 −5 , and for Sod's shock-tube problem: x > 0.5 While the "coaxial" tests used the above initial conditions as they are, the "oblique" tests used the initial conditions rotated about (x, y) = (0.5, 0.5) for 60 • .

Overview of The Simulated Evolution and Analyses
The functions of score against evolution progress exhibit self-similar structure of repeated cliff and plateau. For example, see the evolution history of experiment GA-D illustrated in Fig. 11 and its two enlarged views Fig. 12 and Fig. 13, with the family tree of the best individual superimposed. The evolution of the high scores in the experiments are shown in Fig. 14.
To design more efficient auto-tuning strategies, we investigate what have happened and which part of the evolution contributed to create better individuals in our automated tuning experiments. In section 5.2, we measure the contributions from three different tuning items. In section 5.3, we classify the individuals by their aspects such as their method of birth, their distance to the champion in the family tree, and their fitness relative to their parents. In section 5.4 we study how these classes contribute to the evolution by performing correlation analyses among these classes. In section 5.5, we study how the method of birth of parents affect their children. In section 5.6 we summarize and conclude the analyses. The ID and the scores (in cups) of all the individuals generated in automated tuning experiment GA-D. The individual with the highest score as well as its ancestors are marked by crosses, and each of them is connected to its parents with lines.   The average of GB-370-* was always higher than that of GB-307-*, and the average of GB-333-* was always higher than that of GB-370-* for the 89% of the time.

Contributions of The Three Genome Parts
We assign the symbols to genome parts with different functions as follows: (C): the CUDA kernel execution configuration, (M): which data to store on the memory (to make them Manifest), and (S): when to synchronize the computation. We first investigate how these three components contributed to the score by component-wise artificial crossover between the initial individual and the best scoring individual (Table  10). In the case of GA-S1.33958, introducing improvement only in C, M, S part increase the score by 14%, 30%, and 0%, respectively. In the case of GA-4.33991, the increase are 0%, 85%, 0%. Both cases exhibit synergy effect. Introducing several modifications simultaneously have more effect than the sum of the separate effects, they multiply. So addition in log-space explains 98% and 88% of the progress for GA-S1.33958 and GA-4.33991, respectively, but the score of the final individuals are still slightly higher than predicted.
The score of the individuals created by artificial crossover between the initial individual I 0 and the best scoring individual I . The second to fourth columns indicate which component was taken from which individual. Columns C,M,S correspond to CUDA kernel execution configuration, Manifest/Delay choice, synchronization timing, respectively. For each individual I the fifth column shows µ(I) ± σ(I), the sixth column shows µ(I) µ(I ) − µ(I 0 ) ± σ(I) µ(I ) − µ(I 0 ) , and the seventh column shows log µ(I) − log µ(I 0 ) log µ(I ) − log µ(I 0 ) ± σ(I) (log µ(I ) − log µ(I 0 ))µ(I) .
100%, and 14%, respectively. Removing C,M,S part from GA-S1.33958 decreases the score by 71%, 87%, and −3%, respectively. Again we see the synergy effect, except that GA-S1.33958 were faster if S part were removed. We conclude that the Manifest/Delayed trade-off plays the central part in improving the score. Fixing the Manifest/Delayed nodes, determines the decomposition of the data-flow graph. Then tuning the synchronization timing and CUDA kernel execution configuration help further improve the score.

Classifications of Individuals
To measure how each individual I contributed in generating one of the best individuals, we define contribution distance d(I), and classify the individuals to those who contributed to the evolution and those who didn't. To begin with, let P(I) be the set of individual I's parents. (we use p for probability.) For individual I born by mutation, crossover, and triangulation, the size of the parent set n(P(I)) is 1,2, and 3, respectively.
We define d(I) as follows: • d(I ) = 0 where I is the individual whose µ(I) was the largest in the history.
• d(I) = 0 if one of I's children I c satisfies d(I c ) = 0.
We say that I is [ ] are children significantly faster than any of their parents; [ ] are children whose scores are comparable to their fastest parents; [ ] are children significantly slower than any of their parents; and [≤] are children significantly slower than their best parent but not significantly slower than their slowest parents. Note that children of rank [≤] are never born by mutations since there is only one parent. Table A8 -A14 shows the classifications for experiments.

Chi-squared Tests of Correlations between Classes
To study which set of individuals are contributing to produce best species, we perform Pearson's chi-squared test on pairs of predicates on individuals; see Table A15 and A16. We use 10 −3 significance, or X 2 > 10.83 unless otherwise mentioned. Our observations are as follows: • (row 1.) Significant negative correlations are detected between being born of mutation (n(P(I)) = 1) and being d(I) = 0 for all experiments.
• (row 2.) Significant positive correlations are detected between being born of crossover (n(P(I)) = 2) and being d(I) = 0 for all experiments but GA-S1.
• (row 3.) Significant positive correlations are detected between being born of triangulation (n(P(I)) = 3) and being d(I) = 0 for all experiments.

Correlation along the Family Line
Next, we analyze the correlation along the family line by interpreting the family tree as Markov processes (c.f. Table 11). For each individual I such that d(I) = 0, we trace back its family tree in all possible ways for n steps and obtain an n-letter word. For example, "312" means that I is born of crossover, and one of its parent was born of mutation from an individual who was born of triangulation. With such bag of words obtained, we investigate how the last letter of a word is correlated with letters in front of them. The second columns show X 2 for null hypothesis "The bag of words is a 0th order Markov chain," i.e. in the two-letter words the second letter is not correlated to the first. Formally written, the null hypothesis is: where a, b, c, ... denote characters and ab, abc, ... denote words. The third columns show X 2 for null hypothesis "The bag of words is a 0th order Markov chain," i.e. the probability of the occurrences of the letters depend only on their immediate predecessors and the distribution of three-letter words can be determined  Table 11. Chi-squared test of the family tree being lower-order Markov processes. The each column of the table shows the X 2 statistics of the null hypothesis the family tree being a n-th order Markov process and having no longer correlation.
from the distribution of two-letter words. Formally written, the null hypothesis is: Table 11 shows the result of the Markov chain analysis. Note that, the degree of freedom k for 0th-order and 1st-order Markov chains are 4 and 16, respectively, and the chi-squared values for 10 −3 significance are X 2 > 18.47 and X 2 > 37.70, respectively. From the table, we can deny the 0th-order Markov chain model for all of the experiments, and deny the 1st-order Markov chain model for 8 out of the 10 experiments. We also studied whether an occurence of a character is correlated to its prefix word. "3" was significantly likely to be followed by another "3" in all of the experiments, while "2" was significantly unlikely to be followed by another "2" in 6 out of 10 experiments. On the other hand, three consecutive occurences of the same letter was not significant at least in 8 out of 10 experiments.

Summary of The Analyses
Mutation is not efficient in directly producing individuals with good scores. Nevertheless mutations are indispensable part of the evolution because it is the only way that introduces new genomes to the genome pools and that have chance of finding new improvements.
Triangulations and crossover are two different methods to combine independentlyfound improvements, and are efficient in directly producing good individuals. Triangulations are relatively more efficient in producing d(I) = 0 individuals. However, crossover is relatively good at making large jumps that belongs to [ ] . Triangulations are working in another way; the statistics indicate that one triangulation are likely to be followed by another triangulation in a family tree, and the sequence of Triangulations work to accumulate minor improvements that belongs to [ ] to make a larger one. Fig. 14 -although the number of samples are too few to make a statistically significant statement -suggests that having both mode of combiner is better than having only one.
The main source of performance improvement was the change in memory layout and subroutine re-organization caused by making different Manifest/Delayed trade-offs. Tuning the synchronization timing and CUDA kernel execution configuration provided additional improvements.
The case of GA-S1.33958 is an evidence of the current genetic algorithm not finding the optimal individual. Random mutations tend to increase the genome entropy even under the selection pressure. Adding the rule-based mutations such as "remove all synchronization," which are improbable under random mutations, may contribute to further optimization.

Conclusion and Discussion
We have designed and implemented Paraiso, a domain specific language for describing explicit solvers of partial differential equations, embedded in Haskell. Using the typelevel-tensor and algebraic type classes, we can explain our algorithms with simplicity of mathematical equations. Then the algorithm is translated to OpenMP code or CUDA code. The generated code can be optimized both by applying annotations by hand, and by automated tuning based on a genetic algorithm. Although we present just one example here -a second order Euler equations solver -the front end of Paraiso with typelevel-tensor and Builder Monad readily accept various other algorithms. The code generation and automated tuning can revolutionize the way we invent, implement and optimize various partial differential equation solvers.
The website of Paraiso is http://www.paraiso-lang.org/wiki/, where you can find codes, documentations, slides and videos.
Making more efficient searches for fast individuals is an important future extension of Paraiso. One way is to use machine learning, for example to make suggestions for genomes to introduce, or to predict the scores of the individuals before benchmarking them. Another approach is careful selection of the tuning items, importance of them given either by hand or by inference from the syntactic structure of Paraiso source code. The computing time invested in automated tuning is only justified if the optimized program is used repeatedly, or if optimized single-GPU program is used as a "core" of a multi-GPU program.
Making the automated tuning results more flexible and re-usable is another future challenge. For example, the size of the array e.g. (N 0, N 1) in §2.2 for the twodimensional OM, is fixed before native code generation under the current design of Paraiso. To change the size N 0 or N 1, the users need to re-run the code generation process. One solution is to improve the representation of the constant values in Paraiso, so that the users can set the values that are constant within a run but may vary across different runs. Possible examples of such constants include the OM size and the parameters of the algorithm.
On the other hand, the change in the dimension of the OM e.g. to three; (N 0, N 1, N 2), or the change in the algorithm described by the Builder Monad will change the number of nodes in the OM data flow graph and the bits in the genome. Therefore, the automatically tuned genomes before the change become useless. The genome must be redesigned to perform more re-usable, general automated tunings.
Generating codes for distributed computers / distributed and accelerated computers is another important future extension of Paraiso. To this end, instead of generating MPI codes by ourselves, we plan to generate parallel languages such as X10 [35], Chapel [36], and XcalableMP [15], or domain specific languages for stencil computations such as Physis [5].
Finally, the potential users and developers of Paraiso are obstructed by the programming language Haskell, which is a very computer-science language and is far from being mainstream in the fields of simulation science and high-performance computing. This is a disadvantage of Paraiso compared to other embedded DSL approaches based on more popular syntax such as Python (e.g. [37].) We make two remarks on this.
On the one hand, we show in this paper that abstract concepts developed by computer scientists and found in Haskell, is useful and was almost necessary in implementing a DSL that consistently handles the translation from mathematical notations to evolutionary computation. We hope that such powerful DSLs appear in many different fields, with help of the findings in computer science. On the other hand, we want many other researchers to become developers and users of Paraiso without requiring much knowledge in computer science.
For promotion of such development and use of Paraiso, we plan to provide interfaces to the code generation and the tuning stages in Paraiso, in terms of strings or in lightweight markup languages, along with their specifications. Those will allow many other programming languages can access those stages, such as OM data-flow graphs, their annotations, and the genomes. Adding a simulation scientist friendly scripting language layer on top of Builder Monad is another future work, so that people not familiar with Haskell can access Paraiso. We also plan to publish a series of tutorials and example programs in the Paraiso website. Such collaborations are fundamental to the realization of the future works listed in this section and computer-assisted programming in many fields of simulation science.  Table A1. Contribution distance analysis for experiment GA-1 : The numbers of individuals with specific d(I) that is born of mutation, crossover, and triangulation. The last column shows the total numbers of individuals with specific d(I). The decimal numbers inside the parentheses are the ratio of the number of individuals with specific d(I) of specific birth, divided by total number of individuals of specific birth. Note that the grand total is slightly smaller than that in Table 7; this is because some worker nodes failed to write into DB.   Table A8. Children relative fitness for Experiment GA-1: The individuals classified based upon how they were born and their score with respect to their parents. For each class, the first row shows the number of individuals of that class, and the third row shows the number of individuals of that class with d(I) = 0. The second and fourth rows are first and third row divided by the total number of individuals born in that way. The fifth row is the contribution ratio, the third row divided by the first row (or fourth divided by second).


Construction of Community Web Directories based on Web usage Data

This paper support the concept of a community Web directory, as a Web directory that is constructed according to the needs and interests of particular user communities. Furthermore, it presents the complete method for the construction of such directories by using web usage data. User community models take the form of thematic hierarchies and are constructed by employing clustering approach. We applied our methodology to the ODP directory and also to an artificial Web directory, which was generated by clustering Web pages that appear in the access log of an Internet Service Provider. For the discovery of the community models, we introduced a new criterion that combines a priori thematic informativeness of the Web directory categories with the level of interest observed in the usage data. In this context, we introduced and evaluated new clustering method. We have tested the methodology using access log files which are collected from the proxy servers of an Internet Service Provider and provided results that indicates the usability of the community Web directories. The proposed clustering methodology is evaluated both on a specialized artificial and a community Web directory, indicating its value to the user of the web.

INTRODUCTION
The objective of the paper is building the Web directories according to the preferences of user communities. A Web directory, such as Microsoft (www.microsoft.com) and the Open Directory Project (ODP) (dmoz.org), users can find any topic that they are interested by searching the websites, starting with large categories and gradually moving down and choosing the category related to their interests. Searching for a particular topic ,user has to check deep inside the directory. Hence, the size and the complexity of the Web directory has increased and information overload problem will occur. To overcome the difficulties of web directories, we focus on the construction of community web directories. The users of a community can use the community directory as a starting point for navigating the Web, based on the topics that they are interested in, instead of accessing vast Web directories. The construction of community web directories with usage data collected from the Internet Service Provider raises a number of research issues, which are addressed and solved in this paper. User communities are formed using data collected from the log files of Internet Service Provider. The goal is constructing community web directories based on web usage data. The usage data that form the basis for the construction of community Web directories are collected by the proxy servers of Internet service provider.

Previous work
A Web directory, such as Microsoft (www.microsoft.com) and the Open Directory Project (ODP) (dmoz.org), users can find any topic that they are interested ,by searching the websites starting with large categories and gradually moving down and choosing the category related to their interests. However, the information for the particular topic ,user has to check deep inside the directory. Hence, the size and the complexity of the Web directory has increased and information overload problem will occur, i.e., it is often difficult to navigate to the information of interest to a particular user. On the other hand, Web Personalization, i.e., the task of making Web-based information systems changed to the needs and interests of individual users, or groups of users, to tackle information overload. However, in achieving personalization, we need to create the accurate and operational user models, but creating and maintaining these models are difficult for various reasons.

Present Work
In the concept of a community Web directory, a Web directory is created according to the needs and interests of particular user communities. Furthermore, it presents the complete method for the creation of such directories by using web usage data. User community models take the form of thematic hierarchies and are constructed by using clustering approach. We applied our methodology to the ODP directory, and also to an artificial Web directory, which was generated by clustering Web pages that appear in the log file of a proxy server. For the discovery of the community models, we introduced a new criterion that combines a priori thematic informativeness of the Web directory categories with the level of interest observed in the usage data. We have tested the methodology using access logs from the proxy servers of an Internet Service Provider(ISP) and provided results that indicates the usability of the community Web directories. The results presented here provide an initial measure of the benefits that we can obtain by constructing Web directories to the needs and interests of user communities. However, we have only approximated the gain of the end user and have not bothered about the cost of "losses" that could be encountered in the case of personalized web directory. This issue requires the evaluation of community Web directories according to the user preferences.

ISP Log Format
The Common Log Format is a standardized text file format used by web servers when generating access log files. Because the format is standardized, the files may be analyzed by a number of analysis programs. Each line in a file stored in the Common Log Format has the following syntax: host ident authuser date request status bytes  "GET /apache_pb.gif HTTP/1.0" is the request line from the client. The method GET, /apache_pb.gif the resource requested, and HTTP/1.0 the HTTP protocol.
200 is the HTTP status code returned to the client. 2xx indicates successful response, 3xx indicates redirection, 4xx indicates client error, and 5xx indicates server error.
2326 is the size of the object returned to the client, which is measured in bytes.

Community Web Directory Discovery Algorithm
In this section, we describe the new clustering methodology that we propose for the construction of the community Web directories. This method is used for the selection of categories of the Web directory .
Step2: Extract information from the log file.
Step3: Extract websites from log file.
Step4: Extract directories of the websites.
Step5: Cluster directories based on training set into various clusters. If cluster is not found then cluster is "unspecified"

Evaluation of Community Web Directories
The methodology introduced in this paper for the construction of community Web directories has been tested for analyzing the data collected from the proxy servers. Evaluation of community web directories involves the following steps.

Experimental Setup
Step1: Click on load file.
Step3: Click Get Site.

Conclusions and Future Work
We conclude that the concept of community Web directory specializes to the needs and interests of particular user communities. Furthermore, it presents the complete method for the creation of such directories by using web usage data. User community models are constructed by employing clustering approach. We applied our methodology to the ODP directory, as well as to an artificial Web directory, which was generated by clustering Web pages that appear in the access log files of a Web proxy. For the discovery of the community models, we introduced a new criterion that combines a priori thematic informativeness of the Web directory categories with the level of interest observed in the usage data. We have tested the methodology using access logs from the proxy servers of an Internet Service Provider(ISP) and provided results that indicates the usability of the community Web directories. The proposed methodology addresses the issues of existing method by reducing the dimensionality of the problem, through the classification of individual Web pages into the different categories of the web directory. The results presented here provide an initial measure of the benefits that we can obtain by creating Web directories to the needs and interests of user communities. However, we have only approximated the gain of the end user and have not bothered about the cost of "losses" that could be encountered in the case of personalized web directory. This issue requires the evaluation of community Web directories according to user preferences. Various components of the methodology could be replaced by a number of alternatives. Most importantly, more sophisticated methods for extracting the categories from usage data, in addition to the use of an existing Web directory, would make the mapping of pages to domains and then to categories more accurate and complete.


RELcat: a Relation Registry for ISOcat data categories

The ISOcat Data Category Registry contains basically a flat and easily extensible list of data category specifications. To foster reuse and standardization only very shallow relationships among data categories are stored in the registry. However, to assist crosswalks, possibly based on personal views, between various (application) domains and to overcome possible proliferation of data categories more types of ontological relationships need to be specified. RELcat is a first prototype of a Relation Registry, which allows storing arbitrary relationships. These relationships can reflect the personal view of one linguist or a larger community. The basis of the registry is a relation type taxonomy that can easily be extended. This allows on one hand to load existing sets of relations specified in, for example, an OWL (2) ontology or SKOS taxonomy. And on the other hand allows algorithms that query the registry to traverse the stored semantic network to remain ignorant of the original source vocabulary. This paper describes first experiences with RELcat and explains some initial design decisions.

Introduction
Since early 2009 the ISOcat Data Category Registry (DCR) is operational 1 . A DCR is an ISO 12620:2009 compliant registry for elaborate specifications of data categories (ISO 12620, 2009). Data categories are elementary descriptors in a linguistic resource or annotation scheme. The predecessor of ISOcat SYNTAX, which was based on a draft of the revised ISO 12620 DCR data model allowed some relationships to be included in the specification (Ide and Romary, 2004): 1) value of relationships between simple and closed data categories, and 2) broader generic concept relationships between data categories. The final DCR data model restricted the broader generic concept relationship to be between simple data categories only (Kemps-Snijders, Windhouwer et al., 2008). The reasons for these were: 1. It would already be hard to agree on specifications of standardized data categories, let alone of putting them into a single ideal linguistic ontology, 2. Data categories would be put in ontological relationships based on very domain or user specific views thus hampering reuse in different contexts, and 3. Allowing multiple ontologies to coexist in the registry could lead to endless ontological clutter. These are valid reasons and highlight the goal of the DCR to come to a standardized and reusable set of data categories. However, it was already clear that in the end ontological relationships, especially based on equivalence, would be important to provide crosswalks between various (application) domains and/or between various registries. Now that ISOcat gets nearer to a stable release with relatively complete functionality 2 a 1 See http://www.isocat.org/ 2 Active usage of ISOcat has of course revealed shortcomings which will have to be dealt with in subsequent versions, but version 1.0 will be functionally companion registry is under development whose aim is to store all kinds of relationships and manage the ontological clutter. This Relation Registry is called RELcat (Schuurman and Windhouwer, 2011). This paper will report on the first experiments with RELcat and explain some initial design choices. But the next section will first stress the need for this new registry based on experience gained with the use of ISOcat.

Proliferation of Data Categories
In an ideal world a linguistic concept would be represented by a single data category, whose specification can then be standardized and used by the linguistic community as a whole. However, this idealistic view is already compromised by the DCR data model itself. In this data model data categories can be of various types:  Complex data categories which have a conceptual domain,  Simple data categories which are values in such a conceptual domain, and  Container data categories 3 which group other container or complex data categories. Which data category type is appropriate is directly related to how the data category is used in a linguistic resource. For example, the linguistic concept noun phrase can be realized as either a container data category or a simple data category. This is illustrated in Figure 1 4 where in the feature structure noun phrase is a value of the category attribute, i.e., a simple data category which is part of the conceptual domain of a complex data category, and in the parse tree NP (a common abbreviation for noun phrase) is an inner node, i.e., a container data category.
complete with regard to the original aims. 3 The container data category type isn't part of the data model in ISO 12620:2009, but its addition to ISOcat has been sanctioned by ISO TC 37. 4 Source of the feature structure: http://en.wikipedia.org/wiki/Feature_structure, and source of the parse tree: http://en.wikipedia.org/wiki/Bottom-up_parsing The fact that these two data categories which both realize the concept noun phrase are closely related is valuable information which can be exploited by, for example, semantic search algorithms to find closely related resources. The Relation Registry, RELcat, allows storing and typing these kinds of relationships. Another source of proliferation is the fact that existing sets are bulk loaded into ISOcat, e.g., the NKJP (Patejuk and Przepiórkowski, 2010) or STTS tagsets or the concepts of the GOLD ontology (Farrar and Langendoen, 2010). Ideally these existing sets would all reuse the same concepts and indicate their own encoding using Data Element Names, e.g., for the data category /adverb/ (DC-1232) the following Data Element Names could exist: 1. Name adverbial, source: GOLD 2. Name adv, source: NKJP 3. Name ADV, source STTS However, there is currently no data category that is shared this way. Either due to fundamental differences between data categories, which seem equivalent just based on their names but are deemed deviant enough by the linguist mapping these sets to already existing ISOcat data categories. There might be a stronger drive for harmonization if there was already a standardized core of data categories in ISOcat. However, at the moment all data categories are owned by users and none are standardized and owned by a Thematic Domain Group. Also there is the perceived need to own all data categories which are related to a specific tagset or resource type, and this is indeed a valid concern if the owner of the equivalent data categories is not responding to requests to include the required Data Element Name. The ultimate source of proliferation of data categories in ISOcat is in the end that the registry is based on a grass roots approach, i.e., each user can add data categories that he or she needs. The idea is that in due time, when these data categories have matured, they can be submitted for standardization to one of the Thematic Domain Groups and a coherent, but extensible, core of standardized data categories will appear. However, due to this open approach it is also fairly easy for users to create data categories which are (almost) equivalent to already existing data categories.

RELcat: a Relation Registry
Originally envisioned to be a registry to allow crosswalks between data categories from various registries the role of the Relation Registry has become more important due to the proliferation of data categories within one registry as sketched above. Also it is clear that various (competing) ontological relationships between data categories need to be stored. Driven by the increasing need of storing these relationships the development of a Relation Registry, named RELcat, has started. The natural form of these relationships is a triple of: 1. Subject: the data category the statement is about; 2. Predicate: the relationship between the subject data category and the object data category; 3. Object: the data category the relationship is directed at. As there is the possibility of competing ontological relationships, i.e., private views of specific users, the user should be able to indicate which sets of relationships she wants use. This accommodated by extending the triple to a quad: 4. Set: the set of which this relationship statement is a part. These requirements meet the specifications of a RDF quad store. The current quad store RELcat uses is OpenAnzo 5 which supports queries that indicate which graphs are to be considered. The subject and object URIs in the Relation Registry are obvious: these are the Persistent IDentifiers of the ISOcat data categories or general URLs for concepts from other registries (Windhouwer and Wright, 2012). However, which predicates to use? OWL (W3C, 2009) and SKOS (W3C, 2009) are prime candidates to provide the predicates. But linguistic knowledge resources are already available in various formats, which might not all be readily translated to these formal defined predicates. Also these formal predicates are defined between specific subclasses of RDF resource, e.g., RDF class or RDF property. For example in OWL the predicate owl:equivalentProperty can only be used between to RDF resources of type RDF property. As one of the aims is that 5 See http://www.openanzo.org/ Text="hit" Text="the" Text="ball" Figure 1 The use of the noun phrase concept in a) a feature structure and b) a parse tree RELcat doesn't prescribe design choices for linguistic knowledge bases the taken approach tries to be maximally flexible. It does so by building on the following core taxonomy of relationship types: The same can be done for the other OWL and SKOS predicates, as well as any other existing vocabulary. The benefit of this approach is that 1. existing linguistic knowledge bases can be loaded into or accessed by RELcat as they are, and 2. generic algorithms can use the relationships without intimate knowledge of all these different vocabularies. Relationships between ISOcat data categories and Dublin Core metadata elements can be, for example, specified as follows: @prefix relcat : <http://www.isocat.org/relcat/set/> . @prefix rel : <http://www.isocat.org/relcat/relations#> . @prefix dc : <http://purl.org/dc/elements/1.1/> . @prefix isocat : <http://www.isocat.org/datcat/> . The result will include the Dublin Core metadata elements, but may, when the GOLD relation set is selected, also include GOLD concepts. where relationships are defined in an OWL ontology annotated with ISOcat data category references. But it also still possible to access the original predicates and, for example, combine OWL-based relation sets and feed it into an OWL reasoner. The current alpha version of RELcat supports basic storage and retrieval of relationship sets based on this taxonomy in the OpenAnzo quad store. It also supports SPARQL query templates (so the actual SPARQL is hidden from the user) on (combinations of) relation sets. Sets currently available in RELcat include: 1. relationships between ISOcat metadata categories and Dublin Core metadata elements as created for CMDI (Broeder, Uytvanck et al., 2012), and 2. relationships between ISOcat categories and GOLD concepts as created for the RELISH project (Aristar-Dry, Drude et al., 2012). Although the current version is missing a (web) user interface to add new relationships the collection of sets is already expanding, which indicates the increasing the demand for this kind of registry.

Future work
The current development focuses on supporting and exploiting properties of relation types, like symmetric, transitive and inverses, which will allow the user to specify only a minimal set of relationships. For this RELcat will need a basic reasoner, e.g., based on RDFS Plus (Allemang and Handler, 2008). Depending on the future size of ISOcat and RELcat the use of (almost) same as relationships might require special handling. These (loose) equivalence relationships might possibly lead to a combinatorial explosion. Some commercial triple stores already offer specific support to handle these kinds of large graphs: 1. Oracle Database Semantic Technologies (Oracle, 2005) support so called owl:sameAs cliques which can be consolidated by choosing a clique representative, and 2. OWLIM (ontotext, 2011) supports a similar approach using a so called master node. In RELcat the same kind of approach could be followed: 1) either select one representative from the same-as clique or 2) give the same-as clique its own identifier. A possible complication there is the dynamic combination of relationship sets, which might require that these cliques need to be computed on the fly.

Conclusion
With the increasing usage of ISOcat various causes for proliferation of data categories have appeared. Some of these have fundamental reasons, e.g., due to basic properties of the DCR data model like data category types, and cannot be remedied. To retain semantic interoperability the ontological relationships between these data categories need to be stored in a semantic network to store semantic crosswalks. RELcat is starting to provide the basic facilities to do so, and is specifically aiming at easily integration of already existing linguistic knowledge bases.