Function discrete alloc

Allocation-functions

syntax

• discrete_alloc(TypeNames, LandUnitDomain, SuitabilityMaps, Partitionings, PartitioningNames, AtomicRegions, AtomicRegionMap, MinClaims, MaxClaims, Threshold, FeasibleSolution)

The last argument (FeasibleSolution) may be omitted. Variants for a single or no regional partitioning (discrete_alloc_sp, discrete_alloc_np) and for large sets of land-use types (discrete_alloc_16) are listed on the Allocation functions page.

definition

The discrete_alloc function is used to let the GeoDms perform discrete allocation: it assigns to each land unit the land use type with the highest suitability, under the constraint of regional minimum and maximum claims per land use type.

#	argument	description	type	value type
1	TypeNames	This attribute contains the land-use type names. Its domain unit defines the set of land-use types. For each name, the SuitabilityMaps, MinClaims and MaxClaims containers must contain a subitem with exactly that name.	attribute	string
2	LandUnitDomain	This domain unit defines the set of land units. It can be of any value type, as long as the total number of land units fits in a uint32; it usually defines a grid domain or a subset thereof.	unit	uint32 or point
3	SuitabilityMaps	A container that contains an attribute (domain: LandUnitDomain) for each land-use type that provides suitability for that land-use type for each land unit. All suitability maps must have the same values unit.	container	int32
4	Partitionings	maps each land-use type to the id of the partitioning that defines its claim regions	attribute	uint8
5	PartitioningNames	maps each partitioning to a partitioning name, such as "Province" or "Municipality". For each name, the AtomicRegions unit must contain an attribute with exactly that name that maps each atomic region to its region in that partitioning.	attribute	string
6	AtomicRegions	Unit that defines the set of atomic regions: the smallest common refinement of the regions of all partitionings, i.e. each atomic region is an intersection of one region from each partitioning. It is usually derived with the overlay function, see AtomicRegions.	unit	uint8, uint16 or uint32
7	AtomicRegionMap	attribute that defines to which atomic region each land unit belongs.	attribute	value type of AtomicRegions
8	MinClaims	container with for each land-use type an attribute (domain: the partitioning of that type) with the minimum number of land units that should be allocated within each region to that land-use type	container	uint32
9	MaxClaims	container with for each land-use type an attribute (domain: the partitioning of that type) with the maximum number of land units that may be allocated within each region to that land-use type	container	uint32
10	Threshold	minimum suitability value that a land use type must have on a land unit to be allocatable there. Land units for which all suitabilities are below the threshold remain unallocated (null). A higher threshold reduces the set of allocation alternatives, which can result in violated claims even when the FeasibilityTest succeeded. A value of null allows any suitability.	parameter	int32
11	FeasibleSolution	arbitrary container, reserved as a feasibility certificate for possible future use. It is never evaluated and may be omitted.	container	(empty)

Result

The result of the discrete_alloc function is a container with:

• attribute landuse (domain): an attribute that gives per land unit the resulting land use type. Land units for which all suitabilities are below the Threshold get the null value.
• parameter statusFlag: a boolean value indicating if the allocation was feasible
• parameter status: a text value that explains the value of statusFlag
• container shadow_prices: a container with for each LandUseType an attribute<Price> %Lut.Name%(%Lut.ClaimRegions%) with the resulting shadow price per claim region.
• container total_allocated: a container with for each LandUseType an attribute<UInt32> %Lut.Name%(%Lut.ClaimRegions%).
• attribute bid_price (domain): an attribute that indicates the augmented suitability of the highest bidder per each land unit, i.e. $S_{ij} + \lambda_j$ for the land use type $j$ that the land unit $i$ was allocated to; null for unallocated land units.

If the FeasibilityTest fails, no allocation is attempted: landuse, bid_price, shadow_prices and total_allocated are all filled with null values and only status and statusFlag are informative. If the allocation itself completes with violated claims (statusFlag false), landuse and bid_price still contain the best allocation found, but shadow_prices and total_allocated are filled with null values.

status

Examples of possible resulting status texts are:

DiscrAlloc completed with a total magnified suitability of 48317397677 over 715179 land units = 67559.866379 per land unit, which is optimal

This text means that an optimal allocation has been generated with an average magnified suitability of 67559.866379. Assuming a magnification factor of 100000, this represents an average suitability (aka transition potential) of 0.67559866379

DiscrAlloc.FeasibilityTest failed: total of the minimum claims is 800000 while there are only 715179 allocatable cells, total minimum claims should be decreased with at least 84821 cells;

This text means that no allocation was attempted because the claims themselves cannot be met, see FeasibilityTest.

ClaimRange(type 0 (Urban), Nuts2 3 (ITC3: Liguria)) = [min 25883, max 26057]; 26062 allocated for price {0, 0};

This indicates that five more land units were allocated than the max claim because no alternatives could be found for allocating urban land use without violating other restrictions.

statusFlag

The statusFlag indicates whether a feasible allocation has been generated: it is true if and only if the FeasibilityTest succeeded and the resulting allocation satisfies all minimum and maximum claims. Land units that remain unallocated due to the Threshold do not violate any restriction themselves, but they can cause minimum claims to remain unmet.

total_allocated

The total_allocated container presents a count of the allocation results per land use type per corresponding claim region.

This is similar to what can be produced with reg_count (Landuse, typeNames, RegionContainer, RegionRefs).

where

RegionContainer is a container with names partitionNames and values equal to the corresponding AtomicRegions/AR[atomicRegionMap] and RegionRefs is defined as partitionNames[partitioning].

In future versions of the GeoDMS, the total_allocated sub-container could be removed or made optional, thus the usage of this part of the results is deprecated since they can easily be obtained by using reg_count explicitly.

applies to

since version

example

container Disc_alloc := 
	discrete_alloc(
		  LandUsePreparation/CBSKlasse/name			// 1 string 	attribute
		, Input/Compacted/ADomain				// 2 uint32	unit
		, Input/Compacted/Suitability				// 3 int32	container
		, const(0[Partitioning], LandUsePreparation/CBSKlasse)	// 4 uint8	attribute
		, Partitioning/name					// 5 string	attribute
		, Input/AggRegio					// 6 uint16	unit
		, Input/Compacted/AtomicRegionMap			// 7 uint16	attribute
		, Input/Claims/minclaims				// 8 uint32	container
		, Input/Claims/maxclaims				// 9 uint32	container
		, threshold						// 10 int32	parameter
		, FeasibleSolution					// 11 (empty)	container
	);

frequent error messages

• Value 255 (a.k.a. null-value) out of range of valid Atomic Regions
• Value 65535 (a.k.a. null-value) out of range of valid Atomic Regions
  • cause: the AtomicRegionMap contains values that are no element of the AtomicRegions unit; most often these are null values (255 for uint8, 65535 for uint16) for land units that are not covered by any region. There might be gaps between regions, which will become null with poly2grid.
  • solution: check for null values in the AtomicRegionMap.
• Claim attribute should contain UInt32 values
  • solution: configure all MinClaims and MaxClaims attributes with value type uint32.
• Land Unit set had N elements, but total nr of elements in tiles is M. Use a land unit set with a completely covering tiling
  • solution: the tiling of the LandUnitDomain must cover the domain completely.

discussion

shadow prices as a splitter

Assuming that there is one solution to the optimization problem, knowing the shadow prices that result from the allocation means that the allocation can be regenerated by simply taking, for each land unit $i$, the land use type $j$ with the highest augmented suitability, thus the $j$ for which $S_{ij} + \lambda_j \ge S_{ik} + \lambda_k$ for all $k$. With regional claims, $\lambda_j$ denotes the shadow price of the claim for type $j$ in the claim region that contains land unit $i$.

The $S_{ij}$ can be regarded as the coordinates of N points in a K-dimensional vector space. If we see this space from the perspective of the point given by the coordinates of $-\lambda_j$, we see that the resulting allocation of land unit $i$ is equal to the direction $j$ in which the point $S_i$ has the highest value above $-\lambda$. The hyperplane given by $x_j + \lambda_j = x_k + \lambda_k$ differentiates between allocation towards $j$ respectively $k$. Therefore, the shadow prices can be regarded as a splitter in the K-dimensional vector space of land unit suitabilities.

Therefore, within the context of the discrete_alloc function, the shadow prices are also known as a splitter.

For an elaborate description of the splitter perspective, the sampling and scaling, and the complexity analysis, see Tokuyama, T. and Nakano, J. (1995) Efficient algorithms for the Hitchcock transportation problem. SIAM Journal on Computing 24(3): 563-578; the implementation (DiscrAlloc.cpp) documents an expected time complexity of $O(N*K)$. See also paragraph 1.3.3 of Koomen, E., Hilferink, M. and Borsboom-van Beurden, J. (2011), Introducing Land Use Scanner, Springer.

Degeneration and SoS

Degenerate cases are cases where different allocations have the same total suitability, which means that there is a continuum of solutions to the optimization problem when $X_{ij}$ are allowed to have any real value between 0 and 1.

The discrete_alloc uses virtual perturbation, also known as (Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms, Edelsbrunner et al.) to decide on degenerate cases. Within the discrete_alloc function, each suitability $S\prime_{ij}$ is not taken as $S_{ij}$ but as $S_{ij} + \epsilon * i * j$ where $\epsilon$ is a positive value that is infinitely smaller than the smallest unit of $S_{ij}$. Since $\epsilon$ is infinitesimal, $S_{ij} < S_{ik}$ implies $S\prime_{ij} < S\prime_{ik}$; only when $S_{ij} = S_{ik}$ while $j \neq k$, the simulated perturbation makes a difference of $\epsilon * i * (k-j)$ that decides the comparison. Internally, perturbed values are administered as a pair (value, perturbation coefficient) that is compared lexicographically.

Also, the shadow prices internally have an explicit $\epsilon$ component, which can make the difference between two otherwise equal land units.

scaling

The discrete_alloc function uses a scaling method to find estimates for the shadow_prices before all data is processed. The land units are processed in the order of a systematic sample: with step size $m$, every $m$-th land unit is visited (starting near the middle of the first stride). The initial step size is the smallest power of the step factor 4 for which the sample contains at most 1000 land units. After each pass, the step size is divided by 4 and the next pass visits the land units of the refined sample that were not visited in earlier passes, so each pass approximately quadruples the total number of processed land units; the final pass, with step size 1, processes all remaining land units. Since this is a deterministic systematic sample and not a random sample, the samples are only representative if the order of the land units does not correlate with the sampling pattern; this can be assured by providing the land units (suitabilities and atomic region map) in a (pseudo) random order.

The shadow prices are initially set to 0. At each pass, each claim is scaled to the sample processed so far: the original claim is multiplied by the sampled fraction of the land units of its claim region (counted per atomic region), where scaled minimum claims are rounded down and scaled maximum claims are rounded up. Scaled claims never decrease from one pass to the next and never exceed the original claims, which are restored exactly at the final pass. All visited land units are allocated to the land use type with the highest augmented suitability, which initially equals the original unaugmented suitability.

Once an allocation of land use type $j$ to land unit $i$ results in an overflow of the claim for $j$ (its count exceeds the maximum claim, or exceeds the minimum claim while its shadow price is positive), a procedure is started to maintain the feasibility of the maximum claims, called UpdateSplitterDown, that tries to move the splitter given by the coordinates of the negative shadow price $-\lambda$ in such a way that one land unit less is allocated to $j$ without violating other maximum claims. To achieve this, a heap is maintained for each facet $(j,k)$ between the set of suitability vectors allocatable to $j$ and the set of suitability vectors allocatable to $k$ that contains all land units $i$ that have been allocated to $j$ but could be allocated to $k$. This heap is ordered on $(S\prime_{ij} - S\prime_{ik})$ where the ordering takes the virtual epsilon component into account. After each allocation of land unit $i$ to type $j$, $i$ is pushed into the heaps of all facets of $j$ to other allocation alternatives $k$ (unless the suitability for $k$ is below the Threshold) so that these alternatives are recorded in an accessible way for future reallocations. UpdateSplitterDown then searches the shortest path through the network with claims as nodes and facets as links, where the cost of the link of facet $(j,k)$ is $(S\prime_{ij} + \lambda_j) - (S\prime_{ik} + \lambda_k)$ for the land unit $i$ on top of the facet's heap, until a claim with spare capacity is reached. Note that the ordered heaps give constant time complexity access to the cheapest reallocation alternative from $j$ to $k$, that these link costs are never negative, and that therefore Dijkstra's algorithm can be and is used. The shadow prices of the claims on the found path are then adjusted and the land units on top of the heaps along the path are reallocated. Reallocated land units are added to the heaps of the facets of their new type directly, but their stale entries in the heaps of their old type are only removed when they surface (lazy deletion); shadow price adjustments do not affect the relative order of the land units within each heap, so the heap order is preserved.

After processing all land units of a pass while maintaining the feasibility of the maximum claim restrictions, a mirrored procedure called UpdateSplitterUp is called repeatedly to get sufficient land allocated to meet the (scaled) minimum claims.

Memory usage

overview

The usage of RAM or Virtual Address Space by the discrete_alloc function for allocating $N$ land units (with a non-tiled domain) to $K$ land use types is $O(K*N)$, to be more precise (assuming a uint16 AtomicRegionMap and the uint8 landuse result of discrete_alloc), it is as follows:

• administration of land use types, partitionings and claim region names: insignificant
• administration of claims (with their shadow prices) and facets: proportional to the number of claims resp. facets (bounded by $K*(K-1)$ per atomic region); normally insignificant compared to the $N$-sized arrays below.
• suitabilityMaps: $K$ page-files of $N*sizeof(Int32)$ are mapped into memory.
• atomicRegionMap: 1 page-file of $N*sizeof(UInt16)$ is mapped into memory.
• resulting landuse: 1 page-file of $N*sizeof(UInt8)$ is mapped into memory.
• resulting bid_price: 1 page-file of $N*sizeof(Int32)$, mapped for writing before the allocation starts and filled when it completes.
• ordered heaps with reallocation options: used part bounded by $(K-1)Nsizeof(UInt32)$, but up to a double amount can be taken from the dynamic memory pool.

This totals $(4K+7)*N$ bytes distributed over $K+3$ memory-mapped page files and up to 2 times $(4K-4)*N$ bytes for the ordered heaps.

Historical note: in a Win32 process, where approximately 1.5 GiB was available as heap addresses, the discrete_allocation without selective threshold was safe as long as $(4K+3)N + 2(4K-4)*N = (12K-5)*N < 1.5$ GiB. With $K=13$ as in the EuClueScanner, this simplified to $151N &lt; 1.5$ GiB or roughly $N &lt; 10$ million allocatable land units. Current GeoDMS builds are 64-bit, so the address space no longer imposes a practical limit; the amounts of available RAM and page file remain relevant for performance.

memory usage of the ordered heaps

For each claim j that has an overlapping region with claim $k$, there is a facet with an ordered heap that registers reallocation options from $j$ to $k$. There is one facet for each ordered pair of claims for different land use types that relate to overlapping land units, so the number of facets is bounded by $K*(K-1)$ times the number of atomic regions.

The number of reallocation alternatives that are stored in these heaps is bounded by $(K-1)*N$, except that for performance reasons, the execution of a reallocation from $j$ to $k$ does not directly remove the other alternatives that became obsolete from the different facets of $j$.

The growth strategy of the ordered heaps is that when the required size overflows the reserved memory, a memory reallocation with double this size is requested from the dynamic memory pool. This can result in double the amount of memory reserved for the heaps that have been in use during any step of the discrete_alloc algorithm. Furthermore, the dynamic aspect of moving land unit indices from one facet to the other due to reallocations has not been thoroughly analyzed.

A strict threshold will reduce the amount of (re)allocation options and, therefore, the sizes of these ordered heaps.

Support for tiled data

The domain of the suitability maps and the atomicRegionMap may be tiled (also known as segmented), provided that the tiling covers the domain completely. The current implementation reads the data of all tiles and solves the allocation for the entire domain in one run; the systematic sample of the scaling method spans all tiles.

Previous (GeoDMS 7) versions processed a tiled domain tile by tile, keeping only the data of one tile at a time in memory, with claims scaled per tile and the allocations of earlier tiles frozen: the Tile By Tile Heuristic. An unfortunate permutation of the domain could then result in choices in the first tile(s) that are not optimal or even result in an Infeasible Solution when taking the alternatives of later tiles into account. This heuristic is no longer applied (at least since GeoDMS 8.035); some status texts still refer to (numbers of) tiles.

see also

• concept overview: Discrete Allocation
• overview of Allocation functions
• configuration example: Discrete Allocation Example
• AtomicRegions
• AtomicRegionMap
• FeasibilityTest