​

let P = []

let Q = []

let N = 128

let A = []

const norm = 50

let svd

let bases

let s_vectors

​

const red = '#fb8072'

const green = '#b3de69'

const orange = '#fdb462'

const blue = '#80b1d3'

const blue_orange = [blue, orange]

​

function setup() {

  const w = Math.min(window.innerWidth, 600)

  const h = Math.min(window.innerHeight, 400)

  // console.log(w, h)

  createCanvas(w, h);

  // Centre coordinates

  cx = width / 2;

  cy = height / 2;

  // Set default stroke and colour mode

  colorMode(HSB)

  for (let n = 0; n < N; n++) {

    P[n] = [];

    P[n][0] = norm * cos(n * 2 * PI / N);

    P[n][1] = norm * sin(n * 2 * PI / N);

    Q[n] = [...P[n]];

  }

  A[0] = [1, 0]

  A[1] = [0, 1]

  // Create reset button

  button = createButton('Reset')

  button.position(0, height-20)

  button.mousePressed(reset)

  button.style('background-color', 'black')

  button.style('color', '#bbbbbb')

  button.style('border', 'none')

  bases = [

    new Draggable(norm + cx, cy, red),

    new Draggable(cx, cy - norm, green)

  ]

  s_vectors = [

    new Draggable(0, 0, blue),

    new Draggable(0, 0, orange)

  ]

  compute_singular_vectors()

}

​

function draw() {

  background("black")

​

  // let d = 5

  // Draw input circle, compute output points

  for (let n = 0; n < N; n++) {

    strokeWeight(1)

    draw_segment(P[n], P[(n+1)%N], n * 360/N)

    Q[n] = mult(P[n])    

  }

  // Compute and draw output circle

  for (let n = 0; n < N; n++) {

    strokeWeight(2)

    draw_segment(Q[n], Q[(n+1)%N], n * 360/N)

  }

  // Draw i and j

  noStroke()

  fill(red); ellipse(cx + norm, cy, 10)

  fill(green); ellipse(cx, cy - norm, 10)

  // Display right singular vectors

  for (let i = 0; i < 2; i++) {

    let vx = svd.V.data[0][i] * norm + cx

    let vy = -svd.V.data[1][i] * norm + cy

    fill(blue_orange[i])

    ellipse(vx, vy, 10)

  }

  let over_bases = false

  // Check if i or j have been moved

  for (let i = 0; i < 2; i++) {

    let b = bases[i]

    let xy = b.update()

    if (b.dragging) {

      A[0][i] =  (xy[0] - cx)/norm

      A[1][i] = -(xy[1] - cy)/norm

      compute_singular_vectors()

    }

    over_bases = b.over() || over_bases

    // b.show()

  }

  let over_signular_vectors = false

  // Check if the left singular vectors have been moved

  for (let i = 0; i < 2; i++) {

    let s = s_vectors[i]

    // If not under the bases, then allow them to be dragged

    if (!over_bases) s.update()

    // Draw a line between +v to -v

    let xy = s.xy()

    stroke(s.c)

    line(xy[0], xy[1], 2 * cx - xy[0], 2 * cy - xy[1])

    if (s.dragging & ! over_bases) {

      let ux = svd.U.data[0][i]

      let uy = -svd.U.data[1][i]

      let x = s.x - cx

      let y = s.y - cy

      let len = sqrt((x)**2 + (y)**2)

      let theta = atan2(x * uy - y * ux, x * ux + y * uy)

      // print(theta)

      svd.s[i] = len / norm

      let rot = new mlMatrix.Matrix([

        [cos(theta), -sin(theta)],

        [sin(theta), cos(theta)]

      ])

      let new_U = rot.mmul(svd.U)

      let R = new_U.mmul(mlMatrix.Matrix.diag(svd.s)).mmul(svd.V.transpose())

      A = R.to2DArray()

      update_bases()

      let sx, sy

      sx = new_U.data[0][1-i] * svd.s[1-i] * norm + cx

      sy = -new_U.data[1][1-i] * svd.s[1-i] * norm + cy

      s_vectors[1-i].x = sx

      s_vectors[1-i].y = sy

    }

    over_signular_vectors = s.over() || over_signular_vectors

    s.show()

  }

  // Set cursor to grab, if over buttons, default otherwise

  let over = over_bases || over_signular_vectors

  if (over) cursor('grab')

  else cursor('default')

  // Show them at last, so they appear on top

  for (let b of bases) b.show()

}

​

function draw_segment(r, s, c) {

  // Move point to the centre

  let x1 = r[0] + cx

  let y1 = -r[1] + cy

  let x2 = s[0] + cx

  let y2 = -s[1] + cy

  // We're using the HSB colour space

  stroke(c, 100, 100)

  line(x1, y1, x2, y2)

}

​

// Perhaps I should use the library

function mult(p) {

  let qx, qy

  qx = A[0][0] * p[0] + A[0][1] * p[1]

  qy = A[1][0] * p[0] + A[1][1] * p[1]

  return [qx, qy]

}

​

function mousePressed() {

  bases.forEach(b => b.pressed())  

  s_vectors.forEach(s => s.pressed())  

}

​

function mouseReleased() {

  bases.forEach(b => b.released())

  s_vectors.forEach(s => s.released())

}

​

function mouseDragged() {

}

function compute_singular_vectors() {

  svd = new mlMatrix.SVD(A)

  let S = mlMatrix.Matrix.diag(svd.s)

​

  // Update s'vectors

  for (let i = 0; i < 2; i++) {

    let sx, sy

    sx = svd.U.data[0][i] * svd.s[i] * norm + cx

    sy = -svd.U.data[1][i] * svd.s[i] * norm + cy

    s_vectors[i].x = sx

    s_vectors[i].y = sy

  }

}

​

function update_bases() {

  let xy = mult([norm, 0])

  bases[0].x = xy[0] + cx

  bases[0].y = -xy[1] + cy

  xy = mult([0, norm])

  bases[1].x = xy[0] + cx

  bases[1].y = -xy[1] + cy

}

function reset() {

  A[0] = [1, 0]

  A[1] = [0, 1]

  compute_singular_vectors()

  update_bases()

}