# Transpose tensor indices



Eval_transpose = ->
  push(cadr(p1))
  Eval()

  # add default params if they
  # have not been passed
  if (cddr(p1) == symbol(NIL))
    push_integer(1)
    push_integer(2)
  else
    push(caddr(p1))
    Eval()
    push(cadddr(p1))
    Eval()
  transpose()

transpose = ->
  i = 0
  j = 0
  k = 0
  l = 0
  m = 0
  ndim = 0
  nelem = 0
  t = 0
  ai = []
  an = []
  for i in [0...MAXDIM]
    ai[i] = 0
    an[i] = 0

  #U **a, **b

  save()

  # by default p3 is 2 and p2 is 1
  p3 = pop() # index to be transposed
  p2 = pop() # other index to be transposed
  p1 = pop() # what needs to be transposed

  # a transposition just goes away when
  # applied to a scalar
  if (isNumericAtom(p1))
    push p1
    restore()
    return

  # transposition goes away for identity matrix
  if ((isplusone(p2) and isplustwo(p3)) or (isplusone(p3) and isplustwo(p2)))
    if isidentitymatrix(p1)
      push p1
      restore()
      return

  # a transposition just goes away when
  # applied to another transposition with
  # the same columns to be switched
  if (istranspose(p1))
    innerTranspSwitch1 = car(cdr(cdr(p1)))
    innerTranspSwitch2 = car(cdr(cdr(cdr(p1))))

    if ( equal(innerTranspSwitch1,p3) and equal(innerTranspSwitch2,p2) ) or
      ( equal(innerTranspSwitch2,p3) and equal(innerTranspSwitch1,p2) ) or
      (( equal(innerTranspSwitch1,symbol(NIL)) and equal(innerTranspSwitch2,symbol(NIL)) ) and ((isplusone(p3) and isplustwo(p2)) or ((isplusone(p2) and isplustwo(p3)))))
        push car(cdr(p1))
        restore()
        return

  # if operand is a sum then distribute
  # (if we are in expanding mode)
  if (expanding && isadd(p1))
    p1 = cdr(p1)
    push(zero)
    while (iscons(p1))
      push(car(p1))
      # add the dimensions to switch but only if
      # they are not the default ones.
      push(p2)
      push(p3)
      transpose()
      add()
      p1 = cdr(p1)
    restore()
    return

  # if operand is a multiplication then distribute
  # (if we are in expanding mode)
  if (expanding && ismultiply(p1))
    p1 = cdr(p1)
    push(one)
    while (iscons(p1))
      push(car(p1))
      # add the dimensions to switch but only if
      # they are not the default ones.
      push(p2)
      push(p3)
      transpose()
      multiply()
      p1 = cdr(p1)
    restore()
    return

  # distribute the transpose of a dot
  # if in expanding mode
  # note that the distribution happens
  # in reverse as per tranpose rules.
  # The dot operator is not
  # commutative, so, it matters.
  if (expanding && isinnerordot(p1))
    p1 = cdr(p1)
    accumulator = []
    while (iscons(p1))
      accumulator.push [car(p1),p2,p3]
      p1 = cdr(p1)

    for eachEntry in [accumulator.length-1..0]
      push(accumulator[eachEntry][0])
      push(accumulator[eachEntry][1])
      push(accumulator[eachEntry][2])
      transpose()
      if eachEntry != accumulator.length-1
        inner()

    restore()
    return


  if (!istensor(p1))
    if (!isZeroAtomOrTensor(p1))
      #stop("transpose: tensor expected, 1st arg is not a tensor")
      push_symbol(TRANSPOSE)
      push(p1)
      # remove the default "dimensions to be switched"
      # parameters
      if (!isplusone(p2) or !isplustwo(p3)) and (!isplusone(p3) or !isplustwo(p2))
        push(p2)
        push(p3)
        list(4)
      else
        list(2)
      restore()
      return
    push(zero)
    restore()
    return

  ndim = p1.tensor.ndim
  nelem = p1.tensor.nelem

  # is it a vector?
  # so here it's something curious - note how vectors are
  # not really special two-dimensional matrices, but rather
  # 1-dimension objects (like tensors can be). So since
  # they have one dimension, transposition has no effect.
  # (as opposed as if they were special two-dimensional
  # matrices)
  # see also Ran Pan, Tensor Transpose and Its Properties. CoRR abs/1411.1503 (2014)

  if (ndim == 1)
    push(p1)
    restore()
    return

  push(p2)
  l = pop_integer()

  push(p3)
  m = pop_integer()

  if (l < 1 || l > ndim || m < 1 || m > ndim)
    stop("transpose: index out of range")

  l--
  m--

  p2 = alloc_tensor(nelem)

  p2.tensor.ndim = ndim

  for i in [0...ndim]
    p2.tensor.dim[i] = p1.tensor.dim[i]

  p2.tensor.dim[l] = p1.tensor.dim[m]
  p2.tensor.dim[m] = p1.tensor.dim[l]

  a = p1.tensor.elem
  b = p2.tensor.elem

  # init tensor index

  for i in [0...ndim]
    ai[i] = 0
    an[i] = p1.tensor.dim[i]

  # copy components from a to b

  for i in [0...nelem]

    # swap indices l and m

    t = ai[l]; ai[l] = ai[m]; ai[m] = t;
    t = an[l]; an[l] = an[m]; an[m] = t;

    # convert tensor index to linear index k

    k = 0
    for j in [0...ndim]
      k = (k * an[j]) + ai[j]

    # swap indices back

    t = ai[l]; ai[l] = ai[m]; ai[m] = t;
    t = an[l]; an[l] = an[m]; an[m] = t;

    # copy one element

    b[k] = a[i]

    # increment tensor index

    # Suppose the tensor dimensions are 2 and 3.
    # Then the tensor index ai increments as follows:
    # 00 -> 01
    # 01 -> 02
    # 02 -> 10
    # 10 -> 11
    # 11 -> 12
    # 12 -> 00

    for j in [(ndim - 1)..0]
      if (++ai[j] < an[j])
        break
      ai[j] = 0

  push(p2)
  restore()