Week 8 - Search

var math = require('mathjs')
var Plot = require('plotly')

Problem 15.1:

(a) Plot the Rosenbrock function [Rosenbrock, 1960] $$f = (1 − x)^2 + 100(y − x^2)^2.$$

function rosenbrock(x){
    return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2;
}

function threeDplot(x_data, y_data, z_data, title, range){
//     Display an array assert a graph.
    var NotebookPlot = Plot.createPlot().constructor;
    NotebookPlot.prototype._toHtml = NotebookPlot.prototype.render;
    
    var data = [{x: x_data, y:y_data, z: z_data, type:'surface'}];
    var layout = {title: title, autosize:true, autocolorscale:true} 
    
    var newPlot = Plot.createPlot(data, layout);
    
    $$html$$ = newPlot.render()
}

range = 2

var rosenbrock_x = math.range(-range, range, .1)._data

var rosenbrock_z = [];
for (i=0; i<rosenbrock_x.length; i++){
    var z_row = []
    for (j=0; j<rosenbrock_x.length; j++){
        z_row.push(rosenbrock([rosenbrock_x[i], rosenbrock_x[j]]));
    }
    rosenbrock_z.push(z_row);
}

threeDplot(rosenbrock_x, rosenbrock_x, rosenbrock_z, 'Rosenbrock function', range);

(b) Pick a stopping criterion and use the downhill simplex method to search for its minimum starting from $x = y = −1$, plotting the search path, and comparing the computing time and memory used.

function compare_rosenbrock(a, b){
    return rosenbrock(a) - rosenbrock(b);
}

function insert(arr, func_values, element, element_value){
    var pos = position(element_value, func_values);
    arr.splice(pos, 0, element);
    func_values.splice(pos, 0, element_value);
}

function position(element, array, start, end) {
    start = start || 0;
    end = end || array.length;
    var pivot = parseInt(start + (end - start) / 2, 10);
    
    if (array[pivot] === element) return pivot;
    if (end-start <= 1){
        if (array[pivot] < element){
            return pivot+1;
        } else {
            return pivot;
        }
    } 
    
    if (array[pivot] < element) {
        return position(element, array, pivot, end);
    } else {
        return position(element, array, start, pivot);
    }
}

function shrink(arr, func_values, func){
    for (i=1; i<arr.length; i++){
        arr[i] = math.multiply(.5, math.add(arr[i], arr[0]));
        func_values[i] = func(arr[i]);
    }
}

function plot_simplex(points, func_values){
    var P = math.matrix(points)
    var x = math.flatten(P.subset(math.index(math.range(0, P._size[0]), 0)))
    var y = math.flatten(P.subset(math.index(math.range(0, P._size[0]), 1)))

    var NotebookPlot = Plot.createPlot().constructor;
    NotebookPlot.prototype._toHtml = NotebookPlot.prototype.render;
    
    var data = [{x: rosenbrock_x, y:rosenbrock_x, z: rosenbrock_z, type:'surface'},
               {x: x, y: y, z: func_values, type: 'scatter3d', mode: 'lines+markers'}];
    var layout = {title: 'simplex', autosize:true, autocolorscale:true} 
    
    var newPlot = Plot.createPlot(data, layout);
    
    $$html$$ = newPlot.render()
}

function NelderMeadStep(points, func_values, func){
    var point = points.pop();
    var mean = math.mean(points, 0);
    func_values.pop();
    
    // Reflect
    var new_point = math.subtract(math.multiply(2, mean), point);
    var new_value = func(new_point);
    
    if (new_value < func_values[0]){
        // If refelction worked try reflect and grow
        var doubled_point = math.subtract(math.multiply(3, mean), 
                                          math.multiply(2, point));
        var doubled_value = func(doubled_point);
        
        if (doubled_value < new_value){
            // If reflect and grow worked accept
            insert(points, func_values, doubled_point, doubled_value);
        } else {
            // Else just reflect
            insert(points, func_values, new_point, new_value);
        }
    } else if (new_value < func_values[func_values.length-1]){
        // If reflect wasn't worst or best accept and move on
        insert(points, func_values, new_point, new_value);
    } else {
        // Else try reflect and shrink
        new_point = math.subtract(math.multiply(1.5, mean), 
                                  math.multiply(.5, point));
        new_value = func(new_point);
        
        if (new_value < func_values[func_values.length-1]){
            insert(points, func_values, new_point, new_value);
        } else {
            // If reflect and shrink didn't work try just shrinking
            new_point = math.multiply(.5, math.add(mean, point));
            new_value = func(new_point);
            
            if (new_value < func_values[func_values.length-1]){
                insert(points, func_values, new_point, new_value);
            } else {
                // If nothing worked shrink the full simplex
                points.push(point);
                shrink(points, func_values, func);
            }
        }
    }
}

var start_point = [-1, -1];
var stop_criterion = .001;
var history = [];
var max_num_steps = 1000;

var points = [math.add(start_point, [0.1, 0]), math.add(start_point, [0, 0.1]), start_point];
// points = math.sort(points, compare_rosenbrock);

var func_values = []
for (i=0; i<points.length; i++) func_values.push(rosenbrock(points[i]));
// func_values = math.sort(func_values);

while (func_values[0] > stop_criterion && max_num_steps > 0){
//     plot_simplex(points, func_values);
    NelderMeadStep(points, func_values, rosenbrock);
    max_num_steps--;
//     console.log(points)
    console.log(func_values)
}

console.log(points[0])

(c) Repeat with the direction set method.

function line_min(point, dir, func){
    gold_ratio = (1+math.sqrt(5))/2;
    var max_num_steps = 1000;
    
    var points = [point, math.add(point, dir)];
    var func_values = [func(point), func(points[1])];
    
    if (func_values[0] < func_values[1]){
        dir = math.multiply(-1, dir);
        points = [points[1], points[0]];
        func_values = [func_values[1], func_values[0]];
    } 
    
    points.push(math.add(points[1], math.multiply(gold_ratio, dir)));
    func_values.push(func(points[2]));
    
    while (func_values[2] < func_values[1]){
        points.splice(0, 1);
        func_values.splice(0, 1);
        
        points.push(math.add(points[1], 
                             math.multiply(gold_ratio, 
                                           math.subtract(points[1], points[0]))));
        func_values.push(func(points[2]));
    }

    while ((func_values[2] - func_values[1] > .001 ||
            func_values[0] - func_values[1] > .001) &&
           max_num_steps > 0) {
        var new_point = math.add(points[1], 
                                 math.multiply(1/gold_ratio, 
                                              math.subtract(points[1], points[0])));
        var new_value = func(new_point);
        if (new_value > func_values[1]){
            points = [new_point, points[1], points[0]]
            func_values = [new_value, func_values[1], func_values[0]]
        } else {
            points = [points[1], new_point, points[2]]
            func_values = [func_values[1], new_value, func_values[0]]
        }
        max_num_steps--;
    }
    
    return points[1];
}

function remove_biggest(D, D_contibution){
    var idx = 0;
    
    for (i=0; i<D.length; i++){
        if (D_contibution[i] > D_contibution[idx]){
            idx = i;
        }
    }
    
    D.splice(idx, 1);
    D_contibution.splice(idx, 1);
}

function dir_set(start_point, D, stop_criterion, func){
    var D_contibution = math.zeros(D.length)._data;
    var history = [];

    var point = start_point;
    var value = func(point);

    while (value > stop_criterion){
        var new_point = point;
        
        for (i=0; i<D.length; i++){
            new_point = line_min(new_point, D[i], func);
            new_value = func(new_point);
            D_contibution[i] += value - new_value;
            value = new_value;
            
            console.log(D[i], new_point, value)
        }
        
        new_dir = math.subtract(new_point, point);
        remove_biggest(D, D_contibution);
        D.splice(0, 0,new_dir);
        D_contibution.splice(0, 0, 0);
        
        point = new_point;
    }
}

start_point = [-1, -1];
stop_criterion = .001;
max_num_steps = 1000;
func = rosenbrock;

var D = [[1, 0], [0, 1]];

dir_set(start_point, D, stop_criterion, func)

(d) Repeat with the conjugate gradient method.

function rosenbrock_grad(x){
    return [-2*(1-x[0]) + -400*(x[1] - x[0]**2)*x[0], 200*(x[1] - x[0]**2)]
}

function conjugate_grad(start_point, stop_criterion, func, grad){
    var history = [];

    var point = start_point;
    var dir = grad(point);

    while (func(point) > stop_criterion){
        var old_gradient = grad(point);
        
        point = line_min(point, dir, func);
        
        var new_gradient = grad(point);
        var gamma = math.dot(new_gradient, math.subtract(new_gradient, old_gradient)) / math.dot(old_gradient, old_gradient);
        
        dir = math.add(new_gradient, math.multiply(gamma, dir));  
        
        console.log(point, dir);
    }
}

start_point = [-1, -1];
stop_criterion = .001;
func = rosenbrock;
grad = rosenbrock_grad;

conjugate_grad(start_point, stop_criterion, func, grad)

(e) Repeat with the Levenberg-Marquardt method (using the gradient and Hessian of the function).

Problem 15.2:

Consider a $1D$ spin glass defined by a set of spins $S_1, S_2, . . . , S_N$ , where each spin can be either $+1$ or $−1$. The energy of the spins is defined by $$E = −\sum_{i=1}^N J_iS_iS_{i+1}, $$ where the $J_i$’s are random variables drawn from a Gaussian distribution with zero mean and unit variance, and $S_{N+1} = S1$ (periodic boundary conditions). Find a low-energy configuration of the spins by simulated annealing. The minimum possible energy is bounded by $$E_{min} = −\sum_{i=1}^N |Ji| ;$$ compare your result with this. At each iteration flip a single randomly chosen spin, and if the energy increases by $∆E$ accept the move with a probability $$p = e^{−β∆E}$$ (always accept a move that decreases the energy). Take $β$ to be proportional to time, $β = αt$ (where $t$ is the number of iterations), and repeat the problem for $α = 0.1, 0.01,$ and $0.001$ for $N = 100$. Choose a single set of values for the $J$’s and the starting values for the spins and use these for each of the cooling rates.

// Standard Normal variate using Box-Muller transform.
function randn_bm(N) {
    N = N || 1;
    
    if (N == 1){
        var u = 1 - Math.random(); // Subtraction to flip [0, 1) to (0, 1].
        var v = 1 - Math.random();
        return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v );
    } else {
        var arr = [];
        for (i=0; i<N; i++){
            var u = 1 - Math.random();
            var v = 1 - Math.random();
            arr.push(Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ));
        }
    }
    return arr;
}

function E(S, J){
    var sum = 0;
    for (i=0; i<S.length; i++){
        sum += J[i]*S[i]*S[(i+1)%S.length];
    }
    return -sum;
}

function plotHistory(history, alpha){
    var NotebookPlot = Plot.createPlot().constructor;
    NotebookPlot.prototype._toHtml = NotebookPlot.prototype.render;
    
    var data = [{x: math.range(0, history.length)._data, y:history, type:'line'}];
    var layout = {title: 'Simulated Annealing - parameter = ' + alpha} 
    
    var newPlot = Plot.createPlot(data, layout);
    
    $$html$$ = newPlot.render()
}

N = 100;
S = math.ones(N)._data;
J = randn_bm(N);
max_num_steps = 1000;
minE = -math.sum(math.abs(J));

console.log(minE);

function simulated_annealing(S, J, E_func, alpha, max_num_steps, min_E){
    var energy = E_func(S, J);
    var t = 0;
    var history = [energy];
    
    while (energy > min_E+.001 && t < max_num_steps){  
        var oldS = math.clone(S);
        var old_energy = energy;
        
        var pos = math.randomInt(S.length);
        S[pos] *= -1;
        
        energy = E_func(S, J);
        if (energy > old_energy) {
            var r = math.random();
            if (r > math.exp(-t*alpha*(energy-old_energy))){
                S = oldS;
                energy = old_energy;
            }
        }
        
        t++;
        
        history.push(energy);
    }
    
    return history;
}

alpha = .1;

his = simulated_annealing(S, J, E, alpha, max_num_steps, minE);
plotHistory(his, alpha);

console.log('Minimal E = ' + his[his.length-1])

alpha = .01;

his = simulated_annealing(S, J, E, alpha, max_num_steps, minE);
plotHistory(his, alpha);

console.log('Minimal E = ' + his[his.length-1])

alpha = .001;

his = simulated_annealing(S, J, E, alpha, max_num_steps, minE);
plotHistory(his, alpha);

console.log('Minimal E = ' + his[his.length-1])

Problem 15.3:

Now solve the same problem with a genetic algorithm (keep the same values for the $J$’s as the previous problem). Start with a population of $100$ randomly drawn sets of spins. At each time step evaluate the energy of each member of the population, and then assign it a probability for reproduction of $$p ∝ e^{−β(E−E_{min})}.$$ Generate $100$ new strings by, for each string, choosing two of the strings from the previous population by drawing from this probability distribution, choosing a random crossover point, taking the bits to the left of the crossover point in the first string and the bits to the right in the second, and then mutating by randomly flipping a bit in the resulting string. Plot the minimum energy in the population as a function of time step for $β = 10, 1, 0.1, 0.01$.

population_size = 100;

var population = [];
for (i=0; i<population_size; i++){
    var new_string = [];
    
    for (j=0; j<N; j++){
        r = math.random();
        if (r < .5){
            new_string.push(-1);
        } else {
            new_string.push(1);
        }
    }
    
    population.push(new_string);
}

function new_element(strings, fitnesses){   
    idx = math.range(0, strings.length)._data;
    fitnesses = math.multiply(fitnesses, 1/math.sum(fitnesses));
    
    parents_pos = math.pickRandom(idx, 2, fitnesses);
    parents = [strings[parents_pos[0]], strings[parents_pos[1]]];
        
    crossover = math.randomInt(0, parents[0].length);
    mutation = math.randomInt(0, parents[0].length);
        
    var child = parents[0].slice(0, crossover).concat(parents[1].slice(crossover, parents[1].length));
    child[mutation] *= -1;
    
    return child;
}

function GA(strings, J, E_func, beta, max_num_steps, min_E){
    var t = 0;
    var history = [];
    
    var energies = [];
    var fitnesses = [];

    for (l=0; l<strings.length; l++){
        var energy = E_func(strings[l], J);
        energies.push(energy);
        fitnesses.push(math.exp(-beta*(energy-min_E)));
    }
        
    var minimal_energy = math.min(energies);
    history.push(minimal_energy);

    while (minimal_energy > min_E + .001 && t < max_num_steps){        
        var new_strings = [];
                    
        for (l=0; l<strings.length; l++){
            new_strings.push(new_element(strings, fitnesses));
        } 

        strings = new_strings;
                    
        for (l=0; l<strings.length; l++){
            var energy = E_func(strings[l], J);
            energies[l] = energy;
            fitnesses[l] = math.exp(-beta*(energy-min_E));
        }
        
        minimal_energy = math.min(energies);        
        history.push(minimal_energy);
        
        t++;
    }
    
    return history;
}

// DO NOT WORK!!!
// beta = 10;

// his = GA(population, J, E, beta, max_num_steps, minE);
// plotHistory(his, beta);

// console.log('Minimal E = ' + his[his.length-1])

beta = 1;

his = GA(population, J, E, beta, max_num_steps, minE);
plotHistory(his, beta);

console.log('Minimal E = ' + his[his.length-1])

beta = .1;

his = GA(population, J, E, beta, max_num_steps, minE);
plotHistory(his, beta);

console.log('Minimal E = ' + his[his.length-1])

beta = .01;

his = GA(population, J, E, beta, max_num_steps, minE);
plotHistory(his, beta);

console.log('Minimal E = ' + his[his.length-1])