src/backends/simulators/jssim.js

/*
 * Copyright (c) 2018 Isaac Phoenix (tearsofphoenix@icloud.com).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/*
Contains a (slow) JavaScript simulator.

    Please compile the c++ simulator for large-scale simulations.
*/
import assert from 'assert'
import math from 'mathjs'
import {
  matrixDot,
  matrixRangeAssign,
  matrixRangeIndicesAssign,
  zeros
} from '../../libs/util'
import {
  len, setEqual, complexVectorDot
} from '../../libs/polyfill'
import { stringToArray } from '../../ops/qubitoperator'

/**
 * @class JSSimulator
 * @desc
NodeJS implementation of a quantum computer simulator.

    This Simulator can be used as a backup if compiling the c++ simulator is
not an option (for some reason). It has the same features but is much
slower, so please consider building the c++ version for larger experiments.
 */
export default class Simulator {
  /**
   * @constructor
   */
  constructor() {
    // ignore seed
    this._state = math.ones(1)
    this._map = {}
    this._numQubits = 0
  }

  /**
  Return the qubit index to bit location map and the corresponding state
vector.

    This function can be used to measure expectation values more efficiently (emulation).

    @return {Array}
A tuple where the first entry is a dictionary mapping qubit indices
to bit-locations and the second entry is the corresponding state
vector
   */
  cheat() {
    return [this._map, this._state._data.slice(0)]
  }

  /**
  Measure the qubits with IDs ids and return a list of measurement
outcomes (true/false).

    @param {number[]} ids List of qubit IDs to measure.

    @return {boolean[]} List of measurement results (containing either true or false).
   */
  measureQubits(ids) {
    const P = Math.random()
    let val = 0.0
    let i_picked = 0
    while (val < P && i_picked < len(this._state)) {
      val = math.add(val, math.abs(this._getState(i_picked) || math.complex(0, 0)) ** 2)
      i_picked += 1
    }

    i_picked -= 1

    const res = []
    const pos = ids.map((ID) => {
      res.push(false)
      return this._map[ID]
    })

    let mask = 0
    val = 0

    pos.forEach((looper, i) => {
      res[i] = ((i_picked >> looper) & 1) === 1
      mask |= (1 << looper)
      val |= ((res[i] & 1) << looper)
    })

    let nrm = 0.0
    this._state.forEach((looper, _i) => {
      const i = _i[0]
      if ((mask & i) !== val) {
        this._setState(i, 0.0)
      } else {
        const tmp = math.abs(looper)
        nrm = math.add(nrm, math.multiply(tmp, tmp))
      }
    })
    // normalize
    const scale = 1.0 / Math.sqrt(nrm)
    this._state = math.multiply(this._state, scale)
    return res
  }

  /**
  Allocate a qubit.

    @param {number} ID ID of the qubit which is being allocated.
   */
  allocateQubit(ID) {
    this._map[ID] = this._numQubits
    this._numQubits += 1
    this._state.resize([1 << this._numQubits], 0)
  }

  /**
  Return the classical value of a classical bit (i.e., a qubit which has
been measured / uncomputed).

   @param {number} ID ID of the qubit of which to get the classical value.
   @param {number} tolerance Tolerance for numerical errors when determining
whether the qubit is indeed classical.

    @throws {Error} If the qubit is in a superposition, i.e., has not been measured / uncomputed.
*/
  getClassicalValue(ID, tolerance = 1.e-10) {
    const pos = this._map[ID]
    let up = false
    let down = false

    for (let i = 0; i < len(this._state); i += (1 << (pos + 1))) {
      for (let j = 0; j < (1 << pos); ++j) {
        if (math.abs(this._getState(i + j)) > tolerance) {
          up = true
        }

        if (math.abs(this._getState(i + j + (1 << pos)) || 0) > tolerance) {
          down = true
        }

        if (up && down) {
          throw new Error('Qubit has not been measured / '
          + 'uncomputed. Cannot access its '
          + 'classical value and/or deallocate a '
          + 'qubit in superposition!')
        }
      }
    }

    return down
  }

  /**
  Deallocate a qubit (if it has been measured / uncomputed).

   @param {number} ID ID of the qubit to deallocate.

   @throws {Error} If the qubit is in a superposition, i.e., has not been measured / uncomputed.
   */
  deallocateQubit(ID) {
    const pos = this._map[ID]
    const cv = this.getClassicalValue(ID)
    const newstate = math.zeros(1 << (this._numQubits - 1))
    let k = 0
    for (let i = (1 << pos) * cv; i < len(this._state); i += 1 << (pos + 1)) {
      matrixRangeIndicesAssign(newstate, k, k + (1 << pos), this._state, i)
      k += (1 << pos)
    }

    const newmap = {}
    Object.keys(this._map).forEach((key) => {
      const value = this._map[key]
      if (value > pos) {
        newmap[key] = value - 1
      } else if (parseInt(key, 10) !== ID) {
        newmap[key] = value
      }
    })

    this._map = newmap
    this._state = newstate
    this._numQubits -= 1
  }

  /**
  Get control mask from list of control qubit IDs.

    @return {number} A mask which represents the control qubits in binary.
   */
  getControlMask(ctrlids) {
    let mask = 0
    ctrlids.forEach((ctrlid) => {
      const ctrlpos = this._map[ctrlid]
      mask |= (1 << ctrlpos)
    })
    return mask
  }

  /**
  Emulate a math function (e.g., BasicMathGate).

    @param {function} f Function executing the operation to emulate.
    @param {Array.<number[]>} qubitIDs List of lists of qubit IDs to which
        the gate is being applied. Every gate is applied to a tuple of
        quantum registers, which corresponds to this 'list of lists'.
    @param {number[]} ctrlQubitIDs List of control qubit ids.
   */
  emulateMath(f, qubitIDs, ctrlQubitIDs) {
    const mask = this.getControlMask(ctrlQubitIDs)
    // determine qubit locations from their IDs
    const qb_locs = []
    qubitIDs.forEach((qureg) => {
      qb_locs.push([])
      qureg.forEach((qubitID) => {
        qb_locs[qb_locs.length - 1].push(this._map[qubitID])
      })
    })

    const newstate = math.zeros(len(this._state))

    this._state.forEach((looper, _i) => {
      const i = _i[0]
      if ((mask & i) === mask) {
        const argList = zeros(qb_locs.length)
        qb_locs.forEach((qb, qri) => {
          qb.forEach((il, qi) => {
            argList[qri] |= (((i >> il) & 1) << qi)
          })
        })

        const res = f(argList)
        let newI = i

        qb_locs.forEach((qb, qri) => {
          qb.forEach((il, qi) => {
            if (!(((newI >> il) & 1) == ((res[qri] >> qi) & 1))) {
              newI ^= (1 << il)
            }
          })
        })
        newstate.subset(math.index(newI), looper)
      } else {
        newstate.subset(math.index(i), looper)
      }
    })

    this._state = newstate
  }

  /**
  Return the expectation value of a qubit operator w.r.t. qubit ids.

    @param {Array.<Array>} termsArray Operator Array (see QubitOperator.terms)
    @param {number[]} IDs List of qubit ids upon which the operator acts.

    @return Expectation value
   */
  getExpectationValue(termsArray, IDs) {
    let expectation = 0.0
    const current_state = math.clone(this._state)
    termsArray.forEach(([term, coefficient]) => {
      this.applyTerm(term, IDs)
      const tmp = complexVectorDot(current_state, this._state)
      const delta = math.multiply(coefficient, tmp)
      expectation = math.add(expectation, delta)
      this._state = math.clone(current_state)
    })
    if (math.im(expectation) === 0) {
      return math.re(expectation)
    }
    return expectation
  }

  /**
  Apply a (possibly non-unitary) qubit operator to qubits.

    @param {Array.<Array>} termsArray Operator array (see QubitOperator.terms)
    @param {number[]} IDs List of qubit ids upon which the operator acts.
  */
  applyQubitOperator(termsArray, IDs) {
    let new_state = math.zeros(len(this._state))
    const current_state = math.clone(this._state)
    termsArray.forEach(([term, coefficient]) => {
      this.applyTerm(term, IDs)
      const temp = math.multiply(this._state, coefficient)
      new_state = math.add(new_state, temp)
      this._state = math.clone(current_state)
    })
    this._state = new_state
  }

  /**
  Return the probability of the outcome `bit_string` when measuring
the qubits given by the list of ids.

    @param {boolean[]|number[]} bitString Measurement outcome.
    @param {number[]} IDs List of qubit ids determining the ordering.

    @return Probability of measuring the provided bit string.

    @throws {Error} if an unknown qubit id was provided.
   */
  getProbability(bitString, IDs) {
    const n = IDs.length
    for (let i = 0; i < n; ++i) {
      const id = IDs[i]
      const v = this._map[id]
      if (typeof v === 'undefined') {
        throw new Error('get_probability(): Unknown qubit id. '
        + 'Please make sure you have called '
        + 'eng.flush().')
      }
    }
    let mask = 0
    let bit_str = 0
    for (let i = 0; i < n; ++i) {
      mask |= (1 << this._map[IDs[i]])
      bit_str |= (bitString[i] << this._map[IDs[i]])
    }

    let probability = 0.0

    this._state.forEach((val, _i) => {
      const i = _i[0]
      if ((i & mask) === bit_str) {
        const e = val
        probability += math.re(e) ** 2 + math.im(e) ** 2
      }
    })
    return probability
  }

  /**
   * @ignore
   * @param i
   * @return {*}
   * @private
   */
  _getState(i) {
    return this._state.subset(math.index(i))
  }

  /**
   * @ignore
   * @param i
   * @param value
   * @private
   */
  _setState(i, value) {
    this._state.subset(math.index(i), value)
  }

  /**
  Return the probability amplitude of the supplied `bit_string`.
    The ordering is given by the list of qubit ids.

   @param {boolean[]|number[]} bitString Computational basis state
   @param {number[]} IDs List of qubit ids determining the ordering. Must contain all allocated qubits.

    @return Probability amplitude of the provided bit string.

    @throws {Error} if the second argument is not a permutation of all allocated qubits.
   */
  getAmplitude(bitString, IDs) {
    const s1 = new Set(IDs)
    const s2 = new Set(Object.keys(this._map).map(k => parseInt(k, 10)))
    if (!setEqual(s1, s2)) {
      throw new Error('The second argument to get_amplitude() must'
      + ' be a permutation of all allocated qubits. '
      + 'Please make sure you have called '
      + 'eng.flush().')
    }
    let index = 0
    IDs.forEach((item, i) => {
      item = parseInt(item, 10)
      index |= (bitString[i] << this._map[item])
    })
    const ret = this._getState(index)
    return ret
  }

  /**
Applies exp(-i*time*H) to the wave function, i.e., evolves under
the Hamiltonian H for a given time. The terms in the Hamiltonian
are not required to commute.

    This function computes the action of the matrix exponential using
ideas from Al-Mohy and Higham, 2011.
TODO: Implement better estimates for s.

   @param {Array.<Array>} terms_dict Operator dictionary (see QubitOperator.terms) defining the Hamiltonian.
   @param {number} time Time to evolve for
   @param {number[]} ids A list of qubit IDs to which to apply the evolution.
   @param {number[]} ctrlids A list of control qubit IDs.
  */
  emulateTimeEvolution(terms_dict, time, ids, ctrlids) {
    // Determine the (normalized) trace, which is nonzero only for identity
  // terms:
    let tr = 0
    let sum = 0
    const newTerms = []
    terms_dict.forEach(([t, c]) => {
      if (t.length === 0) {
        tr += c
      } else {
        newTerms.push([t, c])
        sum += Math.abs(c)
      }
    })

    terms_dict = newTerms
    const op_nrm = math.abs(time) * sum
    // rescale the operator by s:
    const s = Math.floor(op_nrm + 1)
    const correction = math.exp(math.complex(0, -time * tr / (s * 1.0)))
    const output_state = math.clone(this._state)
    const mask = this.getControlMask(ctrlids)

    for (let i = 0; i < s; ++i) {
      let j = 0
      let nrm_change = 1.0
      let update
      while (nrm_change > 1.e-12) {
        const coeff = math.divide(math.complex(0, -time), s * (j + 1))
        const current_state = math.clone(this._state)
        update = 0
        terms_dict.forEach(([t, c]) => {
          this.applyTerm(t, ids)
          this._state = math.multiply(this._state, c)

          update = math.add(this._state, update)
          // update += this._state
          this._state = math.clone(current_state)
        })
        update = math.multiply(update, coeff)
        this._state = update

        update.forEach((value, [m]) => {
          if ((m & mask) === mask) {
            const idx = math.index(m)
            const v = math.add(output_state.subset(idx), value)
            output_state.subset(idx, v)
          }
        })
        nrm_change = math.norm(update)
        j += 1
      }

      update.forEach((value, [k]) => {
        if ((k & mask) === mask) {
          const idx = math.index(k)
          const v = math.multiply(output_state.subset(idx), correction)
          output_state.subset(idx, v)
        }
      })
      this._state = math.clone(output_state)
    }
  }

  /**
   * leave it empty to keep same API with cpp simulator
   */
  run() {
    //
  }

  /**
  Applies a QubitOperator term to the state vector. (Helper function for time evolution & expectation)

    @param {Array} term One term of QubitOperator.terms
    @param {number[]} ids Term index to Qubit ID mapping
    @param {number[]} controlIDs Control qubit IDs
  */
  applyTerm(term, ids, controlIDs = []) {
    const X = [[0.0, 1.0], [1.0, 0.0]]
    const Y = [[0.0, math.complex(0, -1)], [math.complex(0, 1), 0.0]]
    const Z = [[1.0, 0.0], [0.0, -1.0]]
    const gates = {X, Y, Z}
    term.forEach((local_op) => {
      const qb_id = ids[local_op[0]]
      this.applyControlledGate(gates[local_op[1]], [qb_id], controlIDs)
    })
  }

  /**
  Applies the k-qubit gate matrix m to the qubits with indices ids,
    using ctrlids as control qubits.

    @param {Array.<Array.<number>>} m 2^k x 2^k complex matrix describing the k-qubit gate.
    @param {number[]} ids A list containing the qubit IDs to which to apply the gate.
    @param {number[]} ctrlids A list of control qubit IDs (i.e., the gate is only applied where these qubits are 1).
   */
  applyControlledGate(m, ids, ctrlids) {
    const mask = this.getControlMask(ctrlids)
    if (len(m) === 2) {
      const k = ids[0]
      const pos = this._map[k] || this._map[k.toString()]
      this._singleQubitGate(m, pos, mask)
    } else {
      const pos = ids.map(ID => this._map[ID])
      this._multiQubitGate(m, pos, mask)
    }
  }

  /**
  Applies the single qubit gate matrix m to the qubit at position `pos`
using `mask` to identify control qubits.

   @param {Array.<Array.<number>>} m 2x2 complex matrix describing the single-qubit gate.
    @param {number} pos Bit-position of the qubit.
    @param {number} mask Bit-mask where set bits indicate control qubits.
   */
  _singleQubitGate(m, pos, mask) {
    const kernel = (u, d, m) => {
      const ma = math.add
      const mm = math.multiply
      d = d || math.complex(0, 0)
      u = u || math.complex(0, 0)
      const r1 = ma(mm(u, m[0][0]), mm(d, m[0][1]))
      const r2 = ma(mm(u, m[1][0]), mm(d, m[1][1]))
      return [r1, r2]
    }

    const step = 1 << (pos + 1)
    for (let i = 0; i < len(this._state); i += step) {
      for (let j = 0; j < (1 << pos); ++j) {
        if (((i + j) & mask) === mask) {
          const id1 = i + j
          const id2 = id1 + (1 << pos)
          const [r1, r2] = kernel(this._getState(id1), this._getState(id2), m)
          this._setState(id1, r1)
          this._setState(id2, r2)
        }
      }
    }
  }

  /**
  Applies the k-qubit gate matrix m to the qubits at `pos`
using `mask` to identify control qubits.

   @param {Array.<number[]>} m 2^k x 2^k complex matrix describing the k-qubit gate.
   @param {number[]} pos List of bit-positions of the qubits.
   @param {number} mask Bit-mask where set bits indicate control qubits.
   see follows the description in https://arxiv.org/abs/1704.01127
   */
  _multiQubitGate(m, pos, mask) {
    const inactive = Object.keys(this._map).map(k => parseInt(k, 10)).filter(p => !pos.includes(p))

    const matrix = math.matrix(m)
    const subvec = zeros(1 << pos.length)
    const subvec_idx = zeros(subvec.length)
    for (let c = 0; c < (1 << inactive.length); ++c) {
      // determine base index (state of inactive qubits)
      let base = 0
      for (let i = 0; i < inactive.length; ++i) {
        base |= ((c >> i) & 1) << inactive[i]
      }

      // check the control mask
      if (mask !== (base & mask)) {
        continue
      }
      // now gather all elements involved in mat-vec mul
      for (let x = 0; x < subvec_idx.length; ++x) {
        let offset = 0
        for (let i = 0; i < pos.length; ++i) {
          offset |= ((x >> i) & 1) << pos[i]
        }
        subvec_idx[x] = base | offset
        subvec[x] = this._getState(subvec_idx[x]) || math.complex(0, 0)
      }

      // perform mat-vec mul
      matrixRangeAssign(this._state, subvec_idx, matrixDot(matrix, subvec))
    }
  }

  /**
  Set wavefunction and qubit ordering.

    @param {Complex[]} wavefunction Array of complex amplitudes describing the wavefunction (must be normalized).
    @param {Array} ordering List of ids describing the new ordering of qubits
(i.e., the ordering of the provided wavefunction).
   */
  setWavefunction(wavefunction, ordering) {
    // wavefunction contains 2^n values for n qubits
    assert(wavefunction.length === (1 << ordering.length))

    // all qubits must have been allocated before
    const f1 = ordering.filter((Id) => {
      const v = this._map[Id]
      return typeof v !== 'undefined'
    }).length === ordering.length
    const f2 = len(this._map) === ordering.length
    if (!f1 || !f2) {
      throw new Error('set_wavefunction(): Invalid mapping provided.'
      + ' Please make sure all qubits have been '
      + 'allocated previously (call eng.flush()).')
    }

    this._state = math.matrix(wavefunction)
    const map = {}
    for (let i = 0; i < ordering.length; ++i) {
      map[ordering[i]] = i
    }
    this._map = map
  }

  /**
  Collapse a quantum register onto a classical basis state.

    @param {number[]} ids Qubit IDs to collapse.
    @param {boolean[]} values Measurement outcome for each of the qubit IDs in `ids`.
    @throws {Error} If probability of outcome is ~0 or unknown qubits are provided.
   */
  collapseWavefunction(ids, values) {
    assert(ids.length === values.length)

    // all qubits must have been allocated before
    const f1 = ids.filter(Id => typeof this._map[Id] !== 'undefined').length === ids.length
    if (!f1) {
      throw new Error('collapse_wavefunction(): Unknown qubit id(s)'
        + ' provided. Try calling eng.flush() before '
        + 'invoking this function.')
    }

    let mask = 0
    let val = 0
    ids.forEach((looper, i) => {
      const pos = this._map[looper]
      mask |= (1 << pos)
      val |= (Math.floor(values[i]) << pos)
    })

    let nrm = 0.0
    this._state.forEach((looper, _i) => {
      const i = _i[0]
      if ((mask & i) === val) {
        nrm += math.abs(this._getState(i)) ** 2
      }
    })

    if (nrm < 1.e-12) {
      throw new Error('collapse_wavefunction(): Invalid collapse! Probability is ~0.')
    }
    const inv_nrm = 1.0 / math.sqrt(nrm)
    this._state.forEach((looper, _i) => {
      const i = _i[0]
      if ((mask & i) !== val) {
        this._setState(i, 0)
      } else {
        this._setState(i, math.multiply(looper, inv_nrm))
      }
    })
  }
}