Update of coord computation, lineToPoints

Author: Stefan Fenn

baysianOptimization.html

@@ -14,20 +14,37 @@
     <canvas id="grid" width="600px" height="600px"></canvas>
     <canvas id="fn2" width="600px" height="600px"></canvas>
     <canvas id="gp_mean" width="600px" height="600px"></canvas>
+    <canvas id="gp_conv" width="600px" height="600px"></canvas>
+    <canvas id="gp_acq_pi" width="600px" height="600px"></canvas>
+    <canvas id="gp_acq_ei" width="600px" height="600px"></canvas>
+    <canvas id="gp_acq_scaled_ei" width="600px" height="600px"></canvas>
     <canvas id="observations" width="600px" height="600px"></canvas>
 </div>
+<div style="width:600px;height:600px;padding-left:620px">
+    <label>Kernel Parameter</label><input type="range" min="1" max="30" value="10"
+                                          onchange="window.updateKernelParam(this.value);"></input>
+    <br/>
+    <button onclick="clearObservations()">Clear Observations</button>
+    <br/>
+    <button onclick="addRandomObservationPoint()">Add random Observation Point</button>
+    <br/>
+    <button onclick="addEiSuggestedObservationPoint()">Add EI suggested Observation Point</button>
+    <br/>
+    <button onclick="addScaledEiSuggestedObservationPoint()">Add ScaledEI suggested Observation Point</button>
+</div>
+
 
 <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjs/6.6.4/math.min.js"></script>
 <script src="graphplot.js" type="module"></script>
 <script type="module">
-    import {drawFunction, drawGrid, drawPoints, xPixToCoord, xCoordToPix} from './graphplot.js';
+    import {clear, drawFunction, drawGrid, drawPoints, lineToPoints, xCoordToPix, xPixToCoord, yCoordToPix} from './graphplot.js';
 
     const canvas = document.getElementById("grid");
     const config = {
-        xMin: -0.2,
-        xMax: 3.4,
-        yMin: -0.2,
-        yMax: 4.2,
+        xMin: -1,
+        xMax: 2,
+        yMin: -2.2,
+        yMax: 2.2,
         ctx: canvas.getContext("2d")
     };
 
@@ -43,29 +60,33 @@
     const configGpMean = Object.assign({}, config);
     configGpMean.ctx = canvasGpMean.getContext("2d");
 
+    const canvasGpConv = document.getElementById("gp_conv");
+    const configGpConv = Object.assign({}, config);
+    configGpConv.ctx = canvasGpConv.getContext("2d");
+
+    const canvasGpAcqPi = document.getElementById("gp_acq_pi");
+    const configGpAcqPi = Object.assign({}, config);
+    configGpAcqPi.ctx = canvasGpAcqPi.getContext("2d");
+
+    const canvasGpAcqEi = document.getElementById("gp_acq_ei");
+    const configGpAcqEi = Object.assign({}, config);
+    configGpAcqEi.ctx = canvasGpAcqEi.getContext("2d");
+
+    const canvasGpAcqScaledEi = document.getElementById("gp_acq_scaled_ei");
+    const configGpAcqScaledEi = Object.assign({}, config);
+    configGpAcqScaledEi.ctx = canvasGpAcqScaledEi.getContext("2d");
+
     drawGrid(config);
     const func = x => {
-        if ((x < 0) || (x > Math.PI)) return 0;
-        let result = 0;
-        for (let i = 1; i <= 10; i++) {
-            result += Math.sin(x) * Math.pow(Math.sin(i * x * x / Math.PI), 20);
-        }
+        let result = -x * Math.sin(10 * x);
+        if ((x < -1) || (x > 2)) return -2;
         return result;
     };
-    drawFunction(configF2, "#000000", func);
-
-    let observationsX = [0.2, 0.5, 0.9, 1.2, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.5, 2.8, 3];
-    let observationsY = observationsX.map(x => func(x));
-    let observations = observationsX.map(x => [x, func(x)]);
-    drawPoints(
-        configObs,
-        "#0000FF",
-        observations
-    );
+    drawFunction(configF2, "#606060", func);
 
-    const kernel = function (x1, x2) {
+    const kernelWithParameter = function (x1, x2, param) {
         let n = math.norm(math.subtract(x1, x2));
-        return math.exp(-0.5 * n * n);
+        return math.exp(-0.5 * param * n * n);
     }
 
     const K = function (x0, x1, kernel) {
@@ -81,54 +102,237 @@
     }
 
     const posterioriMean = function (x, x_new, sigma, y, kernel) {
-        let scope = {
+        const scope = {
             x_new: x_new,
             x: x,
             xl: x.length,
             sigma: sigma,
             y: y,
             K: (x1, x2) => K(x1, x2, kernel)
-        }
+        };
+        //Seite 8 Conditional of a joint Gaussian
+        //http://cs229.stanford.edu/section/more_on_gaussians.pdf
+        //und
+        //GP-Buch Seite 15, 16
         return math.evaluate('K(x_new,x) * inv(K(x,x) + sigma .* identity(xl)) * transpose(y)', scope);
-    }
+    };
+
+    const posteriorCov = function (x, x_new, sigma, kernel) {
+        const scope = {
+            x: x,
+            xl: x.length,
+            sigma: sigma,
+            K: (x1, x2) => K(x1, x2, kernel)
+        };
 
-    const zip = function (a, b) {
-        let zip = [];
-        for (let i = 0; i < a.length; i++) {
-            zip.push([a[i], b[i]]);
+        const r = [];
+        const invSigma = math.evaluate('inv(K(x,x) + sigma .* identity(xl))', scope).valueOf();
+        for (let x_i_new of x_new) {
+            const scope2 = {
+                x_new: [x_i_new],
+                x: x,
+                xl: x.length,
+                sigma: sigma,
+                K: (x1, x2) => K(x1, x2, kernel),
+                invSigma: invSigma
+            };
+            const stdSqr = math.evaluate('K(x_new, x_new) - K(x_new, x) * invSigma * K(x, x_new)', scope2).valueOf()[0][0];
+            r.push(math.abs(stdSqr));
         }
-        return zip;
-    }
+        return r;
+    };
 
-    const x_new = [];
-    const xCoordStart = xCoordToPix(config, -0.2);
-    const xCoordEnd = xCoordToPix(config, 3.4);
-    for (let i = xCoordStart; i < xCoordEnd; i++) {
-        x_new.push(xPixToCoord(config, i))
-    }
-    const y_new = posterioriMean(observationsX, x_new, 1e-10, observationsY, kernel).valueOf();
-    const xy_new = zip(x_new, y_new);
-    drawPoints(
-        configGpMean,
-        "#00FF00",
-        xy_new
-    );
-    drawFunction(configGpMean, "#00ff00", (x)=> {
-        const w = xCoordToPix(config, x);
-        return y_new[w-xCoordStart];
-    });
+    const cdfNormal = function (x) {
+        return (1 - math.erf(-x / Math.sqrt(2))) / 2;
+    };
+
+    const pdfNormal = function (x) {
+        const m = Math.sqrt(2 * Math.PI);
+        const e = Math.exp(-Math.pow(x, 2) / 2);
+        return e / m;
+    };
+
+
+    const acquisitionPI = function (mean, standardDeviation, fMax) {
+        const pi = [];
+        for (let w = 0; w < mean.length; w++) {
+            pi.push(cdfNormal((mean[w] - fMax) / standardDeviation[w], 0, 1));
+        }
+        return pi;
+    };
+
+    const acquisitionEI = function (mean, standardDeviation, fMax, xi) {
+        const eis = [];
+        let maxValue = 0;
+        let minValue = 1000;
+
+        for (let w = 0; w < mean.length; w++) {
+            if (standardDeviation[w] === 0) {
+                eis.push(0);
+                continue;
+            }
+            const std = standardDeviation[w];
+            let dx = mean[w] - fMax - xi;
+            dx = Math.max(-7, Math.min(7, dx / std)) * std;
+            let z = dx / std;
+            let ei = dx * cdfNormal(z) + std * pdfNormal(z);
+
+            if (isNaN(ei)) {
+                ei.push(0);
+                continue;
+            }
+
+            maxValue = Math.max(maxValue, ei);
+            minValue = Math.min(minValue, ei);
+            eis.push(ei);
+        }
+        return eis.map(v => (v - minValue) / (maxValue - minValue));
+    };
+
+    const acquisitionScaledEI = function (mean, standardDeviation, fMax, xi) {
+        const scaledEis = [];
+        let maxValue = 0;
+        let minValue = 1000;
+
+        for (let w = 0; w < mean.length; w++) {
+            if (standardDeviation[w] === 0) {
+                scaledEis.push(0);
+                continue;
+            }
+            const std = standardDeviation[w];
+            let dx = mean[w] - fMax - xi;
+            dx = Math.max(-5, Math.min(5, dx / std)) * std;
+            let z = dx / std;
+            const ei = dx * cdfNormal(z) + std * pdfNormal(z);
+            const vi = std*std*( (z*z + 1) * cdfNormal(z) + z * pdfNormal(z) ) - ei * ei;
+            const scaledEi = ei / math.sqrt(vi);
+
+            if (isNaN(scaledEi)) {
+                scaledEis.push(0);
+                continue;
+            }
+
+            maxValue = Math.max(maxValue, scaledEi);
+            minValue = Math.min(minValue, scaledEi);
+            scaledEis.push(scaledEi);
+        }
+        return scaledEis.map(v => (v - minValue) / (maxValue - minValue));
+    };
+    
+    let observationsX = [-0.75, 0.1, 1.2, 1.65];
+    let observationsY = observationsX.map(x => func(x));
+    let observations = observationsX.map(x => [x, func(x)]);
+    let kernelParam = 10;
+    let acquisitionEi = [];
+    let acquisitionScaledEi = [];
+
+    const drawAll = function () {
+        drawPoints(
+            configObs,
+            "#0000FF",
+            observations
+        );
+
+        const kernel = (x1, x2) => kernelWithParameter(x1, x2, kernelParam);
+        const sigma = 0.000001;
+        const x_new = [];
+        const xCoordStart = xCoordToPix(config, config.xMin);
+        const xCoordEnd = xCoordToPix(config, config.xMax);
+        for (let i = xCoordStart; i < xCoordEnd; i++) {
+            x_new.push(xPixToCoord(config, i))
+        }
+        const mean = posterioriMean(observationsX, x_new, sigma, observationsY, kernel).valueOf();
+        drawFunction(configGpMean, "#000000", (x) => {
+            const w = xCoordToPix(config, x);
+            return mean[w - xCoordStart];
+        });
+
+        const standardDeviation = posteriorCov(observationsX, x_new, sigma, kernel);
+
+        const xs = [...Array(config.ctx.canvas.width).keys()];
+        const ysU = xs.map(x => {
+            return mean[x] + 1.96 * standardDeviation[x];
+        });
+        const ysL = xs.map(x => {
+            return mean[x] - 1.96 * standardDeviation[x];
+        });
+        const acquisitionPi = acquisitionPI(mean, standardDeviation, Math.max(...observationsY));
+        acquisitionEi = acquisitionEI(mean, standardDeviation, Math.max(...observationsY), 0.01);
+        acquisitionScaledEi = acquisitionScaledEI(mean, standardDeviation, Math.max(...observationsY), 0);
+
+        clear(configGpConv);
+        configGpConv.ctx.fillStyle = "#0000b050";
+        configGpConv.ctx.lineWidth = 3;
+        configGpConv.ctx.beginPath();
+        lineToPoints(configGpConv, xs, ysU);
+        lineToPoints(configGpConv, [...xs].reverse(), [...ysL].reverse());
+        configGpConv.ctx.fill();
+
+
+        drawFunction(configGpAcqEi, "#008000", (x) => {
+            const w = xCoordToPix(config, x);
+            return acquisitionEi[w - xCoordStart]  + config.yMin+0.01;
+        });
+
+/*
+        drawFunction(configGpAcqScaledEi, "#000080", (x) => {
+            const w = xCoordToPix(config, x);
+            return acquisitionScaledEi[w - xCoordStart] + config.yMin+0.01;
+        });
+*/
+    };
+
+    window.updateKernelParam = function (param) {
+        kernelParam = param;
+        drawAll();
+    };
+    window.addRandomObservationPoint = function () {
+        const x = math.random() * (config.xMax - config.xMin) + config.xMin;
+        observationsX.push(x);
+        observationsY = observationsX.map(x => func(x));
+        observations = observationsX.map(x => [x, func(x)]);
+        drawAll();
+    };
+    window.addEiSuggestedObservationPoint = function () {
+        const observationsW = observationsX.map(x => xCoordToPix(config, x));
+        function argMax(array) {
+            return array.map((y, i) => [i, y]).filter(a => !observationsW.includes(a[0])).reduce((r, a) => (a[1] > r[1] ? a : r))[0];
+        }
+        const x = xPixToCoord(config, argMax(acquisitionEi));
+        observationsX.push(x);
+        observationsY = observationsX.map(x => func(x));
+        observations = observationsX.map(x => [x, func(x)]);
+        drawAll();
+    };
+    window.addScaledEiSuggestedObservationPoint = function () {
+        const observationsW = observationsX.map(x => xCoordToPix(config, x));
+        function argMax(array) {
+            return array.map((y, i) => [i, y]).filter(a => !observationsW.includes(a[0])).reduce((r, a) => (a[1] > r[1] ? a : r))[0];
+        }
+        const w = argMax(acquisitionScaledEi);
+        const x = xPixToCoord(config, w);
+        observationsX.push(x);
+        observationsY = observationsX.map(x => func(x));
+        observations = observationsX.map(x => [x, func(x)]);
+        drawAll();
+    };
+    window.clearObservations = function () {
+        observationsX = [];
+        window.addRandomObservationPoint();
+    };
 
+    drawAll();
 </script>
 
 
 <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
 <script id="MathJax-script" async
         src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js">
 </script>
-<div style="width:600px;padding-top:620px">
+<div style="width:600px;padding-top:20px">
     <h3>Plot of function</h3>
 
     <p>
-        \[ f(x) = \sum_{i=1}^{10}\sin(x)\ \sin(\frac{i\ x^2}{\pi})^{20} \]
+        \[ f(x) = - x \sin(10 x) \]
     </p>
 </div>

graphplot.js

@@ -91,13 +91,13 @@ export function drawGrid(c, xFactor, yFactor, axisColor, gridColor) {
 export function xCoordToPix(c, xCoord) {
     const xDiff = c.xMax - c.xMin;
     const xPixPerUnit = c.ctx.canvas.width / xDiff;
-    return (xCoord - c.xMin) * xPixPerUnit;
+    return Math.round((xCoord - c.xMin) * xPixPerUnit);
 }
 
 export function yCoordToPix(c, yCoord) {
     const yDiff = c.yMax - c.yMin;
     const yPixPerUnit = c.ctx.canvas.height / yDiff;
-    return c.ctx.canvas.height - (yCoord - c.yMin) * yPixPerUnit;
+    return Math.round(c.ctx.canvas.height - (yCoord - c.yMin) * yPixPerUnit);
 }
 
 export function xPixToCoord(c, xPix) {
@@ -136,6 +136,15 @@ export function drawFunction(c, strokeStyle, func) {
     c.ctx.stroke();
 }
 
+export function lineToPoints(c, wPoints, yPoints) {
+    for (let i=0; i < wPoints.length; i++) {
+        const w = wPoints[i];
+        const y = yPoints[i];
+        const h = yCoordToPix(c, y);
+        c.ctx.lineTo(w, h);
+    }
+}
+
 export function drawPoints(c, fillStyle, points) {
     clear(c);
     c.ctx.fillStyle = fillStyle;