.packageName <- "Malmig"
"Fst" <-
function(rval,N){
  k<-N/sum(N)
  Fst.val<-k%*%diag(rval)
}
"Mal.eq" <-
function(S,P,N){
  phi<-diag(0/N)
  Pt<-t(P)
  x<-0
  repeat{
    x<-x+1
    S1<-mtx.exp(S,x)
    P1<-mtx.exp(P,x)
    Pt1<-mtx.exp(Pt,x)
    D<-(1-phi)/(2*N)
    D<-diag(D)
    D<-diag(D) ## everything till here is similar to a normal phi calculation
    toll<-phi ## I use toll as a comparison mark. toll is phi a n-1 cycles
    toll1<-signif(toll,6) ## optional. I set the number of significant digits to 6
    phi<-phi+(S1%*%Pt1%*%D%*%P1%*%S1) ## that's phi at n cycles
    phi1<-signif(phi,6) ## optional. As for toll
    if (identical(toll1,phi1)){ ## logical condition. If toll (that is, phi for n-1) and phi are identical
      return(x-1) ## return the value of n-1
      break ## and stop, because the Malecot model has reached its asymptot
    }
  }
}
"Phi" <-
function(S,P,N,n){
  if (n < 1){
    return("Number of cycles too low!!!")
  }
  phi<-diag(0/N) ## creating the first phi matrix
  Pt<-t(P)
  x<-0 ## needed for a correct counting cycle
  for (i in 1:n){
    x<-x+1 ## start the counting cycle
    S1<-mtx.exp(S,x) ## powering S
    P1<-mtx.exp(P,x) ## powering P
    Pt1<-mtx.exp(Pt,x) ## powering the transpose of P
    D<-(1-phi)/(2*N) ## calculating the diagonal of the D matrix
    D<-diag(D) ## extracting the diagonal of the above
    D<-diag(D) ## creating the REAL D matix, which is a diagonal matrix
    phi<-phi+(S1%*%Pt1%*%D%*%P1%*%S1) ## Malecot model
  }
}
"R" <-
function(PHI,N){
  k<-N/sum(N) ## the relative population of each k populaion on the total population of the area in study
  rmu<-PHI%*%k ## k is a list, corced to vertical vector. Here I calculate the row wheight phi mean
  mu<-k%*%rmu ## k is now coerced to a linear vector. Here I calculated the overall mean phi
  az<-matrix(rep(rmu,length(rmu)),ncol=length(rmu))
  ax<-az+t(az)
  mu<-as.numeric(mu)
  r.mat<-(PHI+mu-ax)/(1-mu)
}
"col.sto" <-
function(x){
  y<-apply(x,2,sum)
  x1<-t(t(x)/y)
  x1
}
"mtx.exp" <-
function(X,n){
## Function to calculate the n-th power of a matrix X
  if(n != round(n)) {
    n <- round(n)
    warning("rounding exponent `n' to", n)
  }
  phi <- diag(nrow = nrow(X))
  pot <- X # the first power of the matrix.
  while (n > 0)
    {
      if (n %% 2)
        phi <- phi %*% pot
      n <- n %/% 2
      pot <- pot %*% pot
    }
  return(phi)
}
"sym.P" <-
function(x){
  alpha<-x[upper.tri(x)]
  x1<-t(x)
  beta<-x1[upper.tri(x1)]
  gamma<-(alpha+beta)/2
  x[upper.tri(x)]<-gamma
  x2<-t(x)
  x[lower.tri(x)]<-x2[lower.tri(x2)]
  x
}
