# comes from analysis_quanti.r, multicar identical without file test


split_vals_generic = function(targetFile)
{
    values = read.table(targetFile, sep = " ", header = F)
    #Get pop numbers index
    popIndex = values[,1]
    #Get pops ids
    pops = unique(popIndex)
    #Make a mask to select inds for each pops
    splitter = (lapply(pops, '==', popIndex))
    #Isolate values
    values = values[,3]
    #Separate pops
    values = lapply(splitter, function(x) values[x])
    return(list(values, pops))
}

calc_qst_generic = function(species, trait, step, printV)
{
    targetFile = filesList[[argSpecies]][[trait]][[argStep]]
    print(paste("CSVF",targetFile))
    if (file.info(targetFile)$size > 0)
    {
        vp = split_vals_generic(targetFile)
        values = vp[[1]]
        Vb = var(sapply(values, mean))
        Vw = mean(sapply(values, var, na.rm=TRUE))

        if (printV)
        {
            traitName = substring(argTraits, 1, nchar(argTraits)-6) #CARE, works only if trait has structure "Txy-add"
            print(paste("[F] :",traitName,"-> Vb :", Vb))
            print(paste("[F] :",traitName,"-> Vw :", Vw))
        }
        else
        {
            print(paste("Vb : ", Vb))
            print(paste("Vw : ", Vw))
        }
        Qst = Vb / (Vb + 2 * Vw)
        return(Qst)
    }
    else
    {
        print("File empty")
    }
}

pop_mean_generic = function(species, trait, step)
{
    targetFile = filesList[[argSpecies]][[trait]][[argStep]]
    if (file.info(targetFile)$size > 0)
    {
        vp = split_vals_generic(targetFile)
        values = vp[[1]]
        Vb = sapply(values, mean)
        return(Vb)
    }
}

rat_corr_generic = function(A1, A2)
{
    #return(cov(A1, A2)/(var(A1)*var(A2)))
    rateC <- cov(A1, A2)/(sd(A1)*sd(A2))
    if(is.nan(rateC))
    {
        print("ALERT NAN RAT CORR")
        print("A1")
        print(A1)
        print("A2")
        print(A2)
        print("sdA1")
        print(sd(A1))
        print("sdA2")
        print(sd(A2))
        print("covA1A2")
        print(cov(A1, A2))
    }
    return(rateC)
}

valByloc_generic = function(targetFile, argSelected, focus, interestingLoci)
{
  values = read.table(targetFile, sep = " ", header = F)
  #Get pop numbers index
  popIndex = values[,1]
  #Get pops ids
  pops = unique(popIndex)
  npops = length(pops)
  #Make a mask to select inds for each pops
  splitter = (lapply(pops, '==', popIndex))
  #Separate pops
  nloc = length(argSelected)
  valuesByLoc = vector(mode="list", length=nloc)
  valuesByLocMean = vector(mode="list", length=nloc)


  names(valuesByLoc) <- argSelected
  names(valuesByLocMean) <- argSelected

  for (loc in seq(nloc))
  {
      values2 = values[,2+loc]#c(-1, -2)]
      values2 = lapply(splitter, function(x) values2[x])
      valuesByLoc[[loc]] = values2
      valuesByLocMean[[loc]] = lapply(values2, mean)
  }

  #print(paste("VLB:", valuesByLoc))

  nloclt = length(interestingLoci)
  print(interestingLoci)
  Tinter = matrix(nrow=nloclt, ncol=nloclt)
  Tintra = matrix(nrow=nloclt, ncol=nloclt)

  theta_w = 0
  theta_b = 0
  denom_w = 0
  denom_b = 0

  #test covariance
  for(i in seq(length(interestingLoci)))
  {
      for(j in seq(length(interestingLoci)))
      {
          vI = unlist(valuesByLocMean[[interestingLoci[i]]]) #sapply(seq(npops), function(x) valuesByLocMean[[interestingLoci[i]]][x])
          print(paste(i,j))
          if (j != i)
          {
              vJ = unlist(valuesByLocMean[[interestingLoci[j]]])
              Tinter[i,j] = cov(vI, vJ)
              #TIP: to avoid mismatch between pos and pop name, use names() and call by name instead of pos
              allIntra = lapply(seq(npops), function(x) cov(valuesByLoc[[interestingLoci[i]]][[x]], valuesByLoc[[interestingLoci[i]]][[x]]))
              if (nloclt < 5)
              {
                  print(paste("INTRA? ", i, " ", j))
                  print(unlist(allIntra))
              }
              Tintra[i,j] = mean(unlist(allIntra))

              theta_w = theta_w + Tintra[i,j]
              theta_b = theta_b + Tinter[i,j]
          } else {
              Tinter[i,j] = sd(vI)
              print(i)
              print(pops)
              Tintra[i,j] = mean(unlist(lapply(seq(npops), function(x) sd(valuesByLoc[[interestingLoci[i]]][[x]]))))

              denom_w = denom_w + Tintra[i,j]
              denom_b = denom_b + Tinter[i,j]
          }
      }
  }
  print("Tinter")
  print(Tinter)
  print("Tintra")
  print(Tintra)

  theta_w = theta_w / denom_w
  print(paste("[F] :",focus,"-> theta_w :", theta_w))

  theta_b = theta_b / denom_b
  print(paste("[F] :",focus,"-> theta_b :", theta_b))
  return(theta_w)
}
