#!/usr/bin/Rscript --silent

args <- commandArgs(TRUE)

USAGE_STRING = "============= ReMix Testing Script   v. June 28, 2014 ========================
Usage:

$ printscore-admix.R <DISTPATH> <CONDSTRING> <reference> <input> <nvec> <cond-pop>

 	Arguments:
 1 DISTPATH  path to the directory containing readmix.r
 2 Condition string - free or equal, the second for equal weights
 3 Reference file - a CSV file with mean population admixture
 4 Input fie
 5 nvec
 6 Initial populations file
 
 If the directory contains a file with name (your_reference_file_name).sdev , 
 it will be used too.
"

NES_OPTS <- 6

if(length(args) < NES_OPTS){
	cat(USAGE_STRING)
	stop()
}

# Argument parsing
optarg <- 1
DISTPATH = args[[optarg]]
optarg <- optarg + 1

CONDSTRING = args[[optarg]]
optarg <- optarg + 1
# Reference file : a csv file with (mean) admixture proportions of known populations
REFERENCE_FILE = args[[optarg]] # 'Harappa1.csv'
optarg <- optarg + 1
#  Additional file: errors for reference:
REF_SD_FILE = paste(REFERENCE_FILE, ".sdev", sep="")

# Infile: file that has the test sample as its second row (the first row is the header naming the ancestral populations)
INFILE = args[[optarg]]
optarg <- optarg + 1
# Populations to begin with

N_REP = as.numeric(args[[optarg]])
optarg <- optarg + 1

keep.init = TRUE
PRESENT_POPS_FILE = args[[optarg]]
optarg <- optarg + 1
CONDITIONAL = TRUE
	
present_pops = read.table(PRESENT_POPS_FILE, stringsAsFactors = FALSE)

if ( ncol(present_pops) != 1 ){
	stop("Error in the populations file: there should be one population per row \n")
}
present_pops = present_pops[[1]]  # convert to vector

comments = grep("^#", present_pops)
if( length(comments) > 0 ){
	present_pops = present_pops[-comments]
}

keep.init = rep(TRUE, length(present_pops) )

tentative = grep("^[?]", present_pops)
if (length(tentative) > 0){
	present_pops[tentative] = gsub( "^[?]", "", present_pops[tentative])
	keep.init[tentative] = FALSE
}

# needed no more:
# ADM_COEF_FILE_NAME = "admixture-coefs.csv"
# Reordering columns and converting from percentages
prep_input <- function(m) {
	m = as.matrix(m)
	rsums <- rowSums(m)
	if (any(rsums <= 0)) stop("Bad input.")
	m <- m/rsums
	m
}

R = read.csv(REFERENCE_FILE, row.names = 1)
R <- prep_input(R)
K <- dim(R)[2]
# Reference population errors (sdev)
Rsd = NULL
if ( file.exists(REF_SD_FILE) ) {
	Rsd = read.csv(REF_SD_FILE, row.names=1)
	Rsd = as.matrix(Rsd)
 	if ( !all(dim(Rsd) == dim(R)) ) {
		cat("Reference stdev table size (", dim(Rsd),  ")does not match that of reference means", dim(R), ".\n")
		stop()
 	}
} else message("Reference stdev table not found")
 
test = read.csv(INFILE, row.names = 1, header = TRUE) 

library("lpSolve")
 
source(paste(DISTPATH, "/R/readmix_deep.r", sep = ""))

test_group <- as.matrix(test[ , , drop = FALSE])
nmembers <- dim(test_group)[1]
colnames(test_group) <- colnames(test)
rownames(test_group) <- rownames(test)
lpo <- group_solution(R, test_group, present_pops)
pr <- group_prediction(R, test_group, present_pops, lpo)
for (nm in seq(nmembers)) {
	cat(lpo[nm,]$maxerr,"\n")
	cat(pr[nm,]$maxerr,"\n")
}

