defmodule ELA.Matrix do
  alias  ELA.Vector, as: Vector

  @moduledoc"""
  Contains  operations for working with matrices.
  """

  @doc"""
  Returns an identity matrix with the provided dimension.
  """
  @spec identity(number) ::[[number]]
  def identity(n) do
    for i <- 1..n, do: for j <- 1..n, do: (fn
      i, j when i === j -> 1
      _, _ -> 0
    end).(i, j)
  end
  
  @doc"""
  Returns a matrix filled wiht zeroes as with n rows and m columns.
  """
  @spec new(number, number) :: [[number]]
  def new(n, m) do
    for _ <- 1..n, do: for _ <- 1..m, do: 0
  end
  
  @docs"""
  Transposes a matrix.
  """
  @spec transp([[number]]) :: [[number]]
  def transp(a) do
    List.zip(a) |> Enum.map(&Tuple.to_list(&1))
  end

  @doc"""
  Performs elmentwise addition
  """
  @spec add([[number]], [[number]]) :: [[number]]
  def add(a, b) do
    if dim(a) !== dim(b) do
      raise(ArgumentError, "Matrices must have same dimensions for addition.")
    end
    Enum.map(Enum.zip(a, b), fn({a, b}) -> Vector.add(a, b) end)
  end

  @doc"""
  Performs elementwise subtraction
  """
  @spec sub([[number]], [[number]]) :: [[number]]
  def sub(a, b) when length(a) !== length(b),
  do: raise(ArgumentError, "The number of rows in the matrices must match.")
  def sub(a, b) do
    Enum.map(Enum.zip(a, b), fn({a, b}) -> Vector.add(a, Vector.scalar(b, -1)) end)
  end

  @doc"""
  Elementwise mutiplication with a scalar.
  """
  @spec scalar([[number]], number) :: [[number]]
  def scalar(a, s) do
    Enum.map(a, fn(r) -> Vector.scalar(r, s) end)
  end
  
  @doc"""
  Elementwise multiplication with two matrices.
  This is known as the Hadmard product.
  """
  def hadmard(a, b) when length(a) !== length(b),
  do: raise(ArgumentError, "The number of rows in the matrices must match.")
  def hadmard(a, b) do
    Enum.map(Enum.zip(a, b), fn({u, v}) -> Vector.hadmard(u, v) end)
  end

  @doc"""
  Matrix multiplication. Can also multiply matrices with vectors.
  Always returns a matrix.
  """
  def mult(v, b) when is_number(hd(v)) and is_list(hd(b)), do: mult([v], b)
  def mult(a, v) when is_number(hd(v)) and is_list(hd(a)), do: mult(a, [v])
  def mult(a, b) do
    Enum.map(a, fn(r) ->
      Enum.map(transp(b), &Vector.dot(r, &1))
    end)
  end
  
  @doc"""
  Returns a tuple with the matrix dimensions as {rows, cols}.
  """
  @spec dim([[number]]) :: {integer, integer}
  def dim(a) when length(a) === 0, do: 0
  def dim(a) do
    {length(a), length(Enum.at(a, 0))}
  end

  @doc"""
  Pivots them matrix a on the element on row n, column m (zero indexed).
  Pivoting performs row operations to make the 
  pivot element 1 and all others in the same column 0.
  """
  @spec pivot([[number]], number, number) :: [[number]]
  def pivot(a, n, m) do
    pr = Enum.at(a, n)  #Pivot row
    pe = Enum.at(pr, m) #Pivot element

    a
    |> List.delete_at(n)
    |> Enum.map(&Vector.sub(&1, Vector.scalar(pr, Enum.at(&1, m) / pe)))
    |> List.insert_at(n, Vector.scalar(pr, 1 / pe))
  end
  
  @doc"""
  Returns a row equivalent matrix on reduced row echelon form.
  """
  @spec reduce([[number]]) :: [[number]]
  def reduce(a, i \\ 0)  
  def reduce(a, i) do
    r = Enum.at(a, i)
    j = Enum.find_index(r, fn(e) -> e != 0 end)
    a = pivot(a, i, j)
    
    unless j === length(r) - 1 or
           i === length(a) - 1
      do
        reduce(a, i + 1)
      else
	a
    end  end
    
  @doc"""
  Returns the LU-decomposition of a matrix as a tuple {u, l, p}.
  """    
  def lu(a) do
    p = lu_pivot(a)
    a = mult(p, a)
    a = memoize_matrix(a)
    IO.puts "#{inspect map_size(a)}"
    {u, l} = lu_rows(a)
    {u, l, p}
					   
'''
    Enum.each(indexes, &Enum.each(&1, fn
	  {i, j} when i >  j -> u = put_in u[i][j], a[i][j] - sum_upper.(i, j)
	  {i, j} when i == j -> l = put_in l[i][j], 1
	  {i, j} when i <  j -> l = put_in l[i][j], a[i][j] - sum_lower.(i, j)
    end))

    u = Enum.map(u_i, &Enum.map(&1, fn
	{1, j} -> a[1][j]
	{i, j} when i <= j -> a[i][j] - sum_upper.(i, j)
    end))


    indexes = for i <- 1..map_size(a), do:  for j <- 1..map_size(a), do: {i, j}
    a = Enum.map(indexes, &Enum.map(&1, fn
      1, j -> a[1][j]
      i, 1 -> a[i][1]/a[1][1]
      i, j -> a[i][j] - Enum.sum(for k <- 1..i-1, do: a[k][j] * a[i][k])
      i, j -> (a[i][j] - Enum.sum(for k <- 1..j-1, do: a[k][j] * a[i][k]))/a[j][j]
    end))




Börja med första raden och rekursera nedåt
Sätt första element a[i][0]
Generera indexmängd for k <- 1..map_size(a), {i, k}
mappa indexmängden mot rätt värden
modifiera a med en nya raden
passa vidare
avsluta när i == map_size(a)
'''
    # {u, l, p}
  end
  defp lu_rows(a) do lu_rows(a, 1, memoize_matrix(a), memoize_matrix(a)) end
  defp lu_rows(a, i, u, l) when i === map_size(a) + 1, do: {u, l}
  defp lu_rows(a, i, u, l) do
    
  lower_indexes = for j <- 1..i, do: {i, j}
  upper_indexes = for j <- i..map_size(a), do: {i, j}
  IO.puts "lower indexes #{inspect lower_indexes}"
  IO.puts "upper indexes #{inspect upper_indexes}"
  u = memoize_matrix(new(map_size(a), map_size(a)))
  l = memoize_matrix(new(map_size(a), map_size(a)))
    
  l = do_lower(a, lower_indexes, u, l)
  u = do_upper(a, upper_indexes, u, l)

  #IO.puts "lower #{inspect l}"
  IO.puts "upper #{inspect u}"

#    a = Enum.map(a[i], fn
#      {j, v} when i <= j -> sum_upper.(i, j) 
#      {j, v} when i >  j -> sum_lower.(i, j) 
#    end)))

  lu_rows(a, i + 1, u, l)
  end

  defp do_lower(a, [], u, l), do: l
  defp do_lower(a, [{i, j} | t], u, l) do
    IO.puts "doing lower i #{i}, j #{j}"
    IO.puts "a[i][j] #{a[i][j]}"
    #IO.puts "mut #{a[k][j] * a[i][k]}"
    #IO.puts " sum #{Enum.sum(for k <- 1..j-1, do: a[k][j] * a[i][k])/a[j][j]}"
  
    sum_lower = fn
      i, 1 -> a[i][1]/u[1][1]
      i, j -> (a[i][j] - Enum.sum(for k <- 1..j-1, do: u[k][j] * i[i][k]))/u[j][j]
    end
    do_lower(a, t, u, put_in(l[i][j], sum_lower.(i, j)))
  end

  defp do_upper(a, [], u, l), do: u
  defp do_upper(a, [{i, j} | t], u, l) do
    sum_upper = fn
      1, j -> a[1][j]
      i, j -> a[i][j] - Enum.sum(for k <- 1..i-1, do: u[k][j] * i[i][k])
    end
    #IO.puts "doing upper i #{i}, j #{j}"
    do_upper(a, t, put_in(u[i][j], sum_upper.(i, j)), l)
  end
  
  defp lu_decompose_rows(a, u, l, m, i \\ 1)
  defp lu_decompose_rows(a, u, l, m, i) do
    l = put_in l[i][0], u[0][0] / a[i][0]
    m = put_in m[i][0], 0
    {u, l, m} = lu_decompose_row(a, u, l, m, i, 1)
  end

  def lu_decompose_row(a, u, l, m, i, j) do
    IO.puts "doing thing"
    IO.puts "u[#{i}][#{j}] #{u[i][j]}"
    IO.puts "a[#{i}][#{j}] #{a[i][j]}"
    IO.puts "u[#{i} - 1][#{j}] #{u[i - 1][j]}"
    IO.puts "l[#{i}][#{j} - 1}] #{l[i][j - 1]}"
      

    IO.puts "u[#{i}][#{j}] #{u[i][j]}"


    #IO.puts "#{u[j] |> Enum.slice(0, j - 1))}"
    
#    u_k_j = a[i][j] - (u
#      |> Enum.at(j)
#      |> Enum.slice(0, j - 1))

#    IO.puts "ukj #{u_k_j}"
    
#    l_i_k = a[i][j] - (l
#      |> Enum.at(j - 1)
#      |> Enum.slice(0, i))
#    s = Vector.scalar(u_k_j, l_i_k)
    
#    u = put_in u[i][j], a[i][j] - s
    #(u[i - 1][j] * l[i][j - 1]))
#xs    IO.puts "lik #{l_i_k}"

  end

  defp lu_pivot(a), do: lu_pivot(transp(a), identity(length(a)), 0)
  defp lu_pivot(a, p, i) when i === length(a) - 1, do: p
  defp lu_pivot(a, p, i) do
    r = Enum.drop(Enum.at(a, i), i)
    
    j = i + Enum.find_index(r, fn(x) ->
      x === abs(Enum.max(r))
    end)

    p = 
      case {i, j} do
	{i, i} -> p
	_ -> p
	  |> List.update_at(i, &(&1 = Enum.at(p, j)))
	  |> List.update_at(j, &(&1 = Enum.at(p, i)))
      end
    
    lu_pivot(a, p, i + 1)
  end


  defp memoize_matrix(list, map \\ %{}, index \\ 1)
  defp memoize_matrix([], map, _), do: map
  defp memoize_matrix([h|t], map, index) do
    map = Map.put(map, index, memoize_matrix(h))
    memoize_matrix(t, map, index + 1)
  end
  defp memoize_matrix(other, _, _), do: other
end