2d:null 2e:["$","title",null,{"children":"Hello Linear Programming | Jo Wood | Observable"}] 33:T7c7,What is the objective function? Well, our goal is to select the set of 3 decision variables from the 9 available that minimise the geographic distortion involved in fitting to a grid. How we quantify that distortion is up to us, but one simple way would be to measure the squared distance between each point in its original and gridded positions. The sum of these squared distances gives us a 'cost' of allocation that we wish to minimise. Assuming the original points are scaled to fit within a unit square and both the original and gridded coordinates have the origin in the bottom-left, point ${tex`a`} is located at ${tex`(1,0.5)`} and so its cost to allocate to the top-left grid cell ${tex`(0,1)`} would be ${tex`(1-0)^2 + (0.5-1)^2 = 1.25`} We can perform a similar distance cost calculation for the other 8 variables and sum them to give us our objective function: * ${tex`\min z = 1.25x_0 + 0.25x_1 + 1.25x_2 + x_4 + x_5 + 1.09x_6 + 1.49x_7 + 0.09x_8`} What constraints does this operate under? By declaring our decision variables as `Boolean` we already constrain them to be either 0 or 1. But additionally we need to ensure every point is allocated to a cell and that each cell has at most, one point allocated to it: * ${tex`x_0 + x_1 + x_2 = 1`} (${tex`a`} is allocated to a single cell) * ${tex`x_3 + x_4 + x_5 = 1`} (${tex`b`} is allocated to a single cell) * ${tex`x_6 + x_7 + x_8 = 1`} (${tex`c`} is allocated to a single cell) * ${tex`x_0 + x_3 + x_6 \leq 1`} (top-left grid cell is allocated no more than one point) * ${tex`x_1 + x_4 + x_7 \leq 1`} (top-right grid cell is allocated no more than one point) * ${tex`x_2 + x_5 + x_8 \leq 1`} (bottom-left grid cell is allocated no more than one point) When we solve, we should have a set of true (1) and false (0) decision variable values. The variables with true values give us the cell allocations (refer to the diagram of the 9 decision variables above to see which variables correspond to which allocations):34:T45c,{ // Map data const ukMap = await d3 .json("https://gicentre.github.io/data/uk/ukRegions.json") .then((feats) => topojson.feature(feats, feats.objects.regions)); // Grid data const pts = ukRegions.map((o) => [o.lon, o.lat]); const spacers = [ [3, 3], [0, 2], [3, 2], [0, 0] ]; // This does the grid allocation const grid = await pointsToGrid(pts, 4, 4, spacers); // Show map and gridded labels const r = new Renderer(600, 400) .fill("rgba(220,140,95,0.3") .stroke("white") .strokeWidth(1) .translate(-190, 0) .setProjection(d3.geoTransverseMercator().rotate([2, 0]).center([0, 49])) .geoShape(ukMap) .fill("rgb(75,25,17)") .stroke() .pointSize(2) .points(ukRegions.map((row) => Object.values(row).slice(2, 4))) .textSize(9) .textAlign("middle") .textBaseline("bottom"); ukRegions.forEach((row) => r.text(row.short, row.lon, row.lat)); r.fit([ [0, 0], [8, 8] ]) .flipY() .textSize(12) .translate(480, -50); grid.forEach((p, i) => r.text(ukRegions[i].short, ...p)); display(r.render()); }35:T5a1,{ const london = await d3 .json("https://gicentre.github.io/data/geoTutorials/londonBoroughs.json") .then((feats) => topojson.feature(feats, feats.objects.boroughs)); // Find the centroid locations and names directly from the topoJSON. const centroids = london.features.map((borough) => [ borough.id.split(" ")[0], ...getCentroid(borough) ]); // Grid allocation const pts = centroids.map((b) => b.slice(1, 3)); const spacers = [ [0, 0], [1, 0], [2, 0], [4, 0], [5, 0], [6, 0], [7, 0], [0, 1], [1, 1], [6, 1], [7, 1], [0, 2], [0, 5], [1, 5], [6, 5], [7, 5], [0, 6], [1, 6], [2, 6], [3, 6], [5, 6], [6, 6], [7, 6] ]; const grid = await pointsToGrid(pts, 8, 7, spacers); const r = new Renderer(900, 350) // Draw the map .fill("rgba(220,140,95,0.3") .stroke("white") .strokeWidth(1) .translate(-220, 0) .setProjection(d3.geoTransverseMercator().rotate([2, 0]).center([0, 49])) .geoShape(london); // Map labels r.fill("rgb(75,25,17)") .textAlign("middle") .textBaseline("middle") .textSize(6) .stroke(); centroids.forEach((borough) => r.text(...borough)); // Draw grid r.fit([ [0, 0], [8, 8] ]) .flipY() .textSize(12) .translate(710, -50); grid.forEach((p, i) => r.text(centroids[i][0].substring(0, 3), ...p)); display(r.render()); }36:[{"id":"d2e0a749accf512a14a72fe8110be2fea77037d07414236bdd4aee8328b0bb4d6e04e462c86902f45f53188869719c8b905d24764caba5c1dc5ae13290603c0c","href":"https://static.observableusercontent.com/files/d2e0a749accf512a14a72fe8110be2fea77037d07414236bdd4aee8328b0bb4d6e04e462c86902f45f53188869719c8b905d24764caba5c1dc5ae13290603c0c","name":"gridAllocation.png","type":"image/png","lastModified":1628760625629.564,"size":18601},{"id":"1e34df83f75cac05574f4f439296740ed50a9fbea284cd9171120d64f1233b749ef3ed0eb7b3577584e5dabadc3360903872de3d433db2f845f8594ed78209bd","href":"https://static.observableusercontent.com/files/1e34df83f75cac05574f4f439296740ed50a9fbea284cd9171120d64f1233b749ef3ed0eb7b3577584e5dabadc3360903872de3d433db2f845f8594ed78209bd","name":"simpleGrid.png","type":"image/png","lastModified":1628761353374.853,"size":14285}] 31:{"id":"bf6ec1f80070d845","slug":"hello-linear-programming","title":"Hello Linear Programming","thumbnail":"9213b8cecfe0344d8cf68142e77601578f9fca978e239cf136c14c42e6f1dcd1","access_level":"public","create_time":"$D2021-08-11T06:37:11.919Z","update_time":"$D2026-09-08T19:53:13.513Z","view_time":null,"first_public_time":"$D2021-08-13T17:01:36.119Z","has_unpublished_edits":null,"fork":false,"stars":13,"starred":false,"owner":{"type":"team","id":"9461f58b209074d5","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":"starter_2024"},"version":1629,"publish_version":1629,"role":null,"authored":false,"authors":[{"type":"individual","id":"e2432cc13dc447d3","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":null}],"members":[{"type":"individual","id":"e2432cc13dc447d3","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":2,"collection":{"id":"ae8896a1274d5f06","slug":"hello","title":"Hello","owner":{"id":"9461f58b209074d5","login":"jwolondon","name":"Jo Wood"}},"type":"notebook","body":{"title":"Hello Linear Programming","cells":[{"id":0,"value":"# Hello Linear Programming\n\n[Linear Programming](https://en.wikipedia.org/wiki/Linear_programming) provides a way of solving numeric problems that can be expressed as a series of constraints and some objective. Not all optimisation problems can be expressed in this way, especially as a series of linear equations, but it can provide a concise way of solving a large class of them.\n\nThis notebook uses [Tom Larkworthy's](https://observablehq.com/@tomlarkworthy) original [mip parser](https://observablehq.com/@tomlarkworthy/mip) to parse input and access the [JavaScript functions](https://github.com/jvail/glpk.js) that do the solving.\n\nHere are a few simple examples to show how to solve optimisation problems with LP/MIP in Observable.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"## Example 1: Wheat and Barley\n\n*(adapted from this [tutorial example](https://www.analyticsvidhya.com/blog/2017/02/lintroductory-guide-on-linear-programming-explained-in-simple-english/))*\n\nA farmer wishes to plant up some arable fields with wheat and barley. Each crop has a different cost, a different amount of labour to farm it and a different profit when sold:\n\n| Crop | Cost | Labour |Profit |\n| - | - | - | - |\n| Wheat | £100/hectare| 10 days/hectare| £50/hectare|\n| Barley | £200/hectare| 30 days/hectare| £120/hectare|\n\nThe farmer has 110 hectares of land, a budget of £10000 and 1200 person-days of labour available.*** How many hectares of barley and how many of wheat should be grown to maximise profit?***","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":26,"value":"The first task is to express the goal as an *objective function*. Let ${tex`x`} be the number of hectares of wheat planted, ${tex`y`} the number of barley hectares; that is, the variables that will hold the answer to the problem we are trying to solve. These are the *decision variables*. Let ${tex`z`} be the total profit, which given by the question, we wish to maximise:\n\n* ${tex`\\max z = 50x + 120y`}\n\nNext we need to specify the constraints:\n\n* Costs: ${tex`100x + 200y \\leq 10000`}\n* Labour: ${tex`10x + 30y \\leq 1200`}\n* Land area: ${tex`x + y \\leq 110`}\n\nWe can represent the objective function and constraints with Tom's *mip* library as follows:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":38,"value":"const arablePlanting = await mip()\n .objective(Math.max, \"50x + 120y\")\n .subjectTo(\"100x + 200y <= 10000\")\n .subjectTo(\"10x + 30y <= 1200\")\n .subjectTo(\"x + y <= 110\")\n .solve();","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":45,"value":"Inputs.table(\n Object.entries(arablePlanting.vars).map(([k, v]) => ({\n crop: k,\n hectares: v\n })),\n {width: 200}\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":124,"value":"Which gives us our solution: `${arablePlanting.vars.x.toFixed(1)}` hectares of wheat fields and `${arablePlanting.vars.y.toFixed(1)}` hectares of barley. Note that the optimum solution did not involve planting all fields so even in this simple case, the answer may not be immediately intuitive (labour availability prevented more fields from being planted before the constraint of land area kicked in).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":65,"value":"## Example 2: Buttons and Zips\n\n*(adapted from this [Brunel exam question from 1997](http://people.brunel.ac.uk/~mastjjb/jeb/or/morelp.html))*\n\nA company makes two products – *buttons* and *zips* – using two machines, the *stampulator* and the *rotatulator*. Each kilogram of buttons that is produced requires 50 minutes processing time in the stampulator and 30 minutes processing time in the rotatulator. Each kilogram of zips produced requires 24 minutes in the stampulator and 33 minutes of rotatulator time. The factory can run the stampulator for up to 40 hours per week, and the rotatulator up to 35 hours per week.\n\nAt the start of the week the company has 30kg of buttons and 90kg of zips in stock. It will need to supply 75kg of buttons and 95kg of zips to customers but would like to maximise the stock of both by the end of the week. ***How many kgs of buttons and zips should it produce in the week?***","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":83,"value":"The objective function is to maximise the stock of buttons (**x**) and zips (**y**):\n\n* ${tex`\\max z = x + y`}\n\nThe constraints:\n\n* Stampulator time: ${tex`50x + 24y \\leq 2400`} (where 2400 minutes is 40 hours)\n* Rotatulator time: ${tex`30x + 33y \\leq 2100`} (where 2100 minutes is 35 hours)\n* Minimum button production: ${tex`x \\geq 45`} (75kg to supply minus stock of 30kg)\n* Minimum zip production: ${tex`y \\geq 5`} (95kg to supply minus stock of 90kg)","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":63,"value":"const buttonsAndZips = await mip()\n .objective(Math.max, \"x + y\")\n .subjectTo(\"50x + 24y <= 2400\")\n .subjectTo(\"30x + 33y <= 2100\")\n .subjectTo(\"x >= 45\")\n .subjectTo(\"y >= 5\")\n .solve();","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":99,"value":"Inputs.table(\n Object.entries(buttonsAndZips.vars).map(([k, v]) => ({\n Product: k,\n Stock: v\n })),\n {width: 200}\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":234,"value":"## Example 3: LP vs MIP\n\nSo far we have considered cases where the variables we are wishing to find are continuous (area of a crop, weight of some buttons and zips). In the first example the answer happened to be integer values, but they need not have been (try changing the wheat price from £100 per hectare to £101 per hectare to demonstrate this).\n\nLet's consider a more abstract example with four variables and a set of constraints, again using Linear Programming with continuous variables. ***What are the values a, b, c and d that satisfy the following objective function and constraints?***\n\n* ${tex`\\max z = a + b + 2c - 2d`}\n\nConstraints:\n\n* ${tex`a + 2c \\leq 100`}\n* ${tex`2b - 8d \\leq 0`}\n* ${tex`b - 2c + d \\geq 1`}\n* ${tex`a + b + c + d = 10`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":240,"value":"const lp3 = await mip()\n .objective(Math.max, \"a + b + 2c - 2d\")\n .subjectTo(\"a + 2c <= 100\")\n .subjectTo(\"2b - 8d <= 0\")\n .subjectTo(\"b - 2c + d >= 1\")\n .subjectTo(\"a + b + c + d == 10\")\n .solve();","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":250,"value":"Inputs.table(\n Object.entries(lp3.vars)\n .map(([k, v]) => ({\n variable: k,\n value: v\n }))\n .sort((a, b) => a.variable.localeCompare(b.var)),\n {width: 180}\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":352,"value":"giving ${tex`z=${(lp3.vars.a+lp3.vars.b+2*lp3.vars.c-2*lp3.vars.d).toFixed(1)}`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":348,"value":"Some problems might be better represented as integers, such as counts of things. We can add an integer constraint and therefore apply *mixed integer programming (MIP)* to the problem. We specify that our variables are whole numbers by declaring them as `BigInt`. Here is the same objective function and constraints, but this time solved with MIP:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":295,"value":"const mip3 = await mip()\n .objective(Math.max, \"a + b + 2c - 2d\")\n .subjectTo(\"a + 2c <= 100\")\n .subjectTo(\"2b - 8d <= 0\")\n .subjectTo(\"b - 2c + d >= 1\")\n .subjectTo(\"a + b + c + d == 10\")\n .var(\"a\", BigInt)\n .var(\"b\", BigInt)\n .var(\"c\", BigInt)\n .var(\"d\", BigInt)\n .solve()","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":297,"value":"Inputs.table(\n Object.entries(mip3.vars)\n .map(([k, v]) => ({\n variable: k,\n value: v\n }))\n .sort((a, b) => a.variable.localeCompare(b.var)),\n {width: 180}\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":378,"value":"giving ${tex`z=${(mip3.vars.a+mip3.vars.b+2*mip3.vars.c-2*mip3.vars.d)}`}, which although it happens to be a rounded version of the LP solution, has quite different values for the variables ${tex`a`} to ${tex`d`}.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":117,"value":"## Example 4: Simple grid map allocation\n\nSuppose we have a set of geographic locations as points in 2d space and we wish to represent them as grid cells in a matrix with as few cells as possible. We would like to allocate each location to the grid cell that best represents its geographic position. So we might expect more northerly locations to be in the top rows of the grid, westerly locations on the left etc.\n\nWhile it may not seem quite as obvious as the previous examples, this allocation can also be expressed as an objective function under a set of constraints solvable via LP/MIP.\n\nConsider a simple case of 3 points ${tex`a`}, ${tex`b`} and ${tex`c`} that we wish to allocate to three cells in a 2x2 grid. For this example, we will just consider the upper left three cells as candidates. The fourth unused cell in the bottom right we will refer to as a *spacer*.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":464,"value":" $\"Simple$ ","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":470,"value":"Each point should be allocated to one of three places, giving nine possible allocations in total from which we need to select the optimal three. We can assign a Boolean decision variable to each of these 9 allocations:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":456,"value":" $\"9$ ","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":474,"value":"$33","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":427,"value":"const mip3a = await mip()\n // Minimise cost of allocation\n .objective(\n Math.min,\n \"1.25x0 + 0.25x1 + 1.25x2 + x4 + x5 + 1.09x6 + 1.49x7 + 0.09x8\"\n )\n // Point must be allocated to one cell in grid:\n .subjectTo(\"x0+x1+x2 == 1\")\n .subjectTo(\"x3+x4+x5 == 1\")\n .subjectTo(\"x6+x7+x8 == 1\")\n\n // Grid cell can have at most one point allocated to it\n .subjectTo(\"x0+x3+x6 <= 1\")\n .subjectTo(\"x1+x4+x7 <= 1\")\n .subjectTo(\"x2+x5+x8 <= 1\")\n\n .var(\"x0\", Boolean)\n .var(\"x1\", Boolean)\n .var(\"x2\", Boolean)\n .var(\"x3\", Boolean)\n .var(\"x4\", Boolean)\n .var(\"x5\", Boolean)\n .var(\"x6\", Boolean)\n .var(\"x7\", Boolean)\n .var(\"x8\", Boolean)\n\n .solve();","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":445,"value":"Inputs.table(\n Object.entries(mip3a.vars)\n .map(([k, v]) => ({\n variable: k,\n value: v\n }))\n .sort((a, b) => a.variable.localeCompare(b.var)),\n {width: 180}\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1058,"value":"## Example 5: Larger grid map allocation\n\nThe examples considered so far are all relatively simple ones; simple enough that we could express the objective function and constraint expressions directly in algebraic form. Tom Larkworthy's [mip package](https://observablehq.com/@tomlarkworthy/mip) did the work of translating those algebraic expressions into a form usable by the underlying glpk solver library. For more complex problems we can skip the algebraic formulation and provide the LP directly to [glpk](https://github.com/jvail/glpk.js) as a JSON expression which we generate programmatically.\n\nLet's generalise the grid allocation LP so that it works with a set of points of any size. The first step is to create an array of grid locations with the option of excluding some as 'spacers':","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":878,"value":"function gridLocations(numCols, numRows, spacers = []) {\n const isSpacer = (c, r) => {\n let matched = false;\n spacers.forEach(([sc, sr]) => {\n if (sr == r && sc == c) matched = true;\n });\n return matched;\n };\n\n const grd = [];\n for (let row = 0; row < numRows; row++) {\n for (let col = 0; col < numCols; col++) {\n if (!isSpacer(col, numRows - 1 - row)) {\n grd.push([col, numRows - 1 - row]); // Make grid origin bottom-left.\n }\n }\n }\n return grd;\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1075,"value":"The JSON file we need to create for the LP solver has three parts, the objective function, the constraints and the types of decision variables (all binary in our grid example):\n``` \n{\n objective: {...}\n subjectTo: [{constraint1}, {{constraint2}, ...}]\n types: [v1, v2, ...]\n}\n```\nWe can create functions to generate the necessary objects based on the input points:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1094,"value":"function objective(pts, grd) {\n // Scale points to unit square\n const xDomain = d3.extent(pts.map((p) => p[0]));\n const yDomain = d3.extent(pts.map((p) => p[1]));\n const xNorm = d3.scaleLinear().domain(xDomain).range([0, 1]);\n const yNorm = d3.scaleLinear().domain(yDomain).range([0, 1]);\n const ptsNorm = pts.map(([x, y]) => [xNorm(x), yNorm(y)]);\n\n // This is how costly it is to move from point (x,y) to grid position (c,r)\n const cost = ([x, y], [c, r]) => (x - c) * (x - c) + (y - r) * (y - r);\n\n // Build the expression from the coefficients\n const vars = [];\n let i = 0;\n for (let a = 0; a < ptsNorm.length; a++) {\n for (let b = 0; b < ptsNorm.length; b++) {\n const coeff = cost(ptsNorm[a], grd[b]);\n if (coeff != 0) {\n vars.push({name: \"x\" + i, coef: coeff});\n }\n i++;\n }\n }\n return vars;\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1119,"value":"function pointToOneCell(n) {\n const constraint = (a) =>\n Array.from(Array(n)).map((_, i) => ({name: \"x\" + (a + i), coef: 1}));\n\n return Array.from(Array(n)).map((_, i) => ({\n name: \"pointToOneCell\" + (i + 1),\n vars: constraint(i * n),\n bnds: {type: glpk.GLP_FX, ub: 1.0, lb: 1.0}\n }));\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1150,"value":"function cellMaxOnePoint(n) {\n const constraint = (a) =>\n Array.from(Array(n)).map((_, i) => ({name: \"x\" + (a + i * n), coef: 1}));\n\n return Array.from(Array(n)).map((_, i) => ({\n name: \"cellMaxOnePoint\" + (i + 1),\n vars: constraint(i),\n bnds: {type: glpk.GLP_UP, ub: 1.0, lb: 0}\n }));\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1162,"value":"function vars(n) {\n return Array.from(Array(n * n)).map((_, i) => \"x\" + i);\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1241,"value":"And then combine them to create the LP JSON, send it to the solver and use the solution to allocate grid positions in the same order as input points.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1106,"value":"async function pointsToGrid(pts, numCols, numRows, spacers = []) {\n const grid = gridLocations(numCols, numRows, spacers);\n const lp = {\n name: \"gridMapLP\",\n objective: {\n direction: glpk.GLP_MIN,\n name: \"cost\",\n vars: objective(pts, grid)\n },\n subjectTo: pointToOneCell(pts.length).concat(cellMaxOnePoint(pts.length)),\n binaries: vars(pts.length)\n };\n\n const solution = await glpk.solve(lp, glpk.GLP_MSG_ALL);\n\n const gPos = Object.entries(solution.result.vars)\n .filter((v) => v[1] == 1) // Select true variables\n .map((v) => Number(v[0].slice(1))) // The x-number corresponds to grid allocation.\n .sort((a, b) => a - b); // Keep grid positions in original point order\n\n return gPos.map((v) => grid[v % gPos.length]);\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":935,"value":"Let's apply this gridding to a real example. Here are the point locations representing the standard regions of the UK:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":159,"value":"const ukRegions = [\n { name: \"North East\", short: \"NE\", lon: -1.85, lat: 55.05 },\n { name: \"North West\", short: \"NW\", lon: -2.77, lat: 54.45 },\n { name: \"Yorkshire and The Humber\", short: \"YH\", lon: -1.29, lat: 53.93 },\n { name: \"East Midlands\", short: \"EM\", lon: -0.85, lat: 52.81 },\n { name: \"West Midlands\", short: \"WM\", lon: -2.21, lat: 52.56 },\n { name: \"East\", short: \"Ea\", lon: 0.51, lat: 52.24 },\n { name: \"London\", short: \"Ln\", lon: -0.15, lat: 51.49 },\n { name: \"South East\", short: \"SE\", lon: -0.81, lat: 51.21 },\n { name: \"South West\", short: \"SW\", lon: -3.63, lat: 50.81 },\n { name: \"Wales\", short: \"Wa\", lon: -3.99, lat: 52.07 },\n { name: \"Scotland\", short: \"Sc\", lon: -3.97, lat: 56.18 },\n { name: \"Northern Ireland\", short: \"NI\", lon: -6.45, lat: 54.61 }\n]","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":165,"value":"Inputs.table(ukRegions, { rows: 13, layout: \"auto\", width: 360 })","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":135,"value":"$34","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1379,"value":"And here are the boroughs of London with their grid layout. In this example we insert many spacers around the edge of the region (for more details on approaches to inserting spacers, see [Meulemans et al 2017](https://openaccess.city.ac.uk/id/eprint/15167/)):","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1383,"value":"$35","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1510,"value":"To see a more thorough example of using LP to allocate regions to grid cells, see this [Grid map allocation page](https://observablehq.com/@jwolondon/gridmap-allocation).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":11,"value":"---\n\n## Appendix\n\n*The mixed integer programming library, noting the GPL licence for the solver backend:*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":9,"value":"import {mip, glpk} from \"observable:@tomlarkworthy/mip\";","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1397,"value":"function getCentroid(feat) {\n if (feat.geometry.type === \"MultiPolygon\") {\n const maxPoly = ([aMax, coords], cs) => {\n const area = d3.polygonArea(cs.flat());\n return area > aMax ? [area, cs] : [aMax, coords];\n };\n const [area, coords] = feat.geometry.coordinates.reduce(maxPoly, [0, []]);\n return d3.polygonCentroid(coords.flat());\n } else {\n return d3.polygonCentroid(feat.geometry.coordinates.flat());\n }\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":147,"value":"import { Renderer } from \"observable:@jwolondon/renderer\"","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W36","stdlib":"2","titling":"auto"}} 32:"$31:body"