include("Permutations.jl")
using Permutations
using PyPlot


function simulate(N::Int64, M::Int64, K::Int64)

  E1 = zeros(Float64,N)
  Eq1 = zeros(Float64,N)
  V1 = zeros(Float64,N)
  D1 = zeros(Float64,N)

  E2 = zeros(Float64,N)
  Eq2 = zeros(Float64,N)
  V2 = zeros(Float64,N)
  D2 = zeros(Float64,N)

  E3 = zeros(Float64,N)
  Eq3 = zeros(Float64,N)
  V3 = zeros(Float64,N)
  D3 = zeros(Float64,N)

  for n = 1:N
    for step = 1:M
      perm = Permutations.RandomPermutation(n)
      tmp1 = getNoOfFixedPoints(perm)
      E1[n] += tmp1
      Eq1[n] += tmp1^2

      tmp2 = 0
      if K != 0
        tmp2 = getKSizeCycles(K, perm)
        E2[n] += tmp2
      else
        tmp2 = length(Permutations.cycles(perm))
        E2[n] += tmp2
      end
      Eq2[n] += tmp2^2

      tmp3 = getNumberOfRecords(perm)
      E3[n] += tmp3
      Eq3[n] += tmp3^2
    end

    E1[n] /= M
    Eq1[n] /= M
    V1[n] = Eq1[n] - (E1[n])^2
    D1[n] = sqrt(V1[n])

    E2[n] /= M
    Eq2[n] /= M
    V2[n] = Eq2[n] - (E2[n])^2
    D2[n] = sqrt(V2[n])

    E3[n] /= M
    Eq3[n] /= M
    V3[n] = Eq3[n] - (E3[n])^2
    D3[n] = sqrt(V3[n])


  end

  printDiagram(E1, V1, D1, 0)
  printDiagram(E2, V2, D2, 1)
  printDiagram(E3, V3, D3, 2)

end

function getNoOfFixedPoints(perm::Permutations.Permutation)
  noOfFixedPoints = 0

  for i = 1:length(perm)
    if perm[i] == i
      noOfFixedPoints +=1
    end
  end

  return noOfFixedPoints
end

function getKSizeCycles(k::Int64, perm::Permutations.Permutation)
  cycles = Permutations.cycles(perm)
  noOfKSizeCycles = 0

  for i = 1:length(cycles)
    if length(cycles[i]) == k
      noOfKSizeCycles += 1
    end
  end

  return noOfKSizeCycles
end

function getNumberOfRecords(perm::Permutations.Permutation)
  noOfRecords = 1
  max = perm[1]

  for element = 2:length(perm)
    if element > max
      max = element
      noOfRecords += 1
    end
  end

  return noOfRecords
end

function Var(E::Array, Eq::Array)
  var = 0
  for i = 1:length(E)
    var = Eq[i] - (E[i])^2
  end
  return var
end


function Deviation(Var::Float64)
    return sqrt(Var)
end

function printDiagram(E::Array, V::Array, D::Array, case::Int64)

  D1 = zeros(Float64,N)
  D2 = zeros(Float64,N)
  for i = 1:length(D)
    D1[i] = D[i] + 3
    D2[i] = D[i] - 3
  end
  msg = ""
  if case == 0
    msg = " fixed points number"
  elseif case == 1
    msg = " cycles number"
  elseif case == 2
    msg = " records number"
  end
  x = 1:N
  fig, ax = subplots()
  suptitle("Permutations")

  ax[:plot](x, E, linewidth=2, alpha=0.6, label="Average value of$msg")
  ax[:plot](x, V, linewidth=2, alpha=0.6, label="Variation of$msg")
  ax[:plot](x, D1, linewidth=2, alpha=0.6, label="(+)Deviation of$msg")
  ax[:plot](x, D2, linewidth=2, alpha=0.6, label="(-)Deviation of$msg")
  ax[:legend]()
  fileName = string("permutations_$msg.png")
  savefig(fileName, dpi=72)

end



N = 100
M = 100000
K = 0

simulate(N, M, K)