#make reltype = 6+number.generations-currun OR just the specific number
#(PAR$id.vector[number.generations+8]+1):PAR$id.vector[number.generations+9]


#for debugging

 #females.married <-TEMP$females.married
 #reltpe <- TEMP$reltpe
 #children.vector <- TEMP$children.vector
 #corrections <- list()
 #founder.pop <- TRUE

#save(males.effects.mate,females.effects.mate,patgenes,matgenes,TEMP$size.mate,PAR$num.par,a.effects.key,aa.effects.key,
#                           d.effects.key,reltpe,PAR$id.vector,PAR$genenames.row,PAR$parnames,a.mean.correction,
#                           a.sd.correction,aa.mean.correction,aa.sd.correction,d.mean.correction,d.sd.correction,
#                           f.mean.correction,f.sd.correction,s.mean.correction,
#                           s.sd.correction,e.mean.correction,e.sd.correction,age.mean.correction,
#                           age.sd.correction,slope.mean.correction,slope.sd.correction,vt.model,sibling.age.cor,
#                           females.married,children.vector,mvrnormal,effects.mate,beta.matrix,PAR$am.multiplier,file="toMPFunc.RData")

#load("toMPFunc.RData")



make.phenotype <- function(males.effects.mate, females.effects.mate, patgenes, matgenes, females.married, reltpe, children.vector,
                           corrections, twin.age=1, founder.pop=FALSE, center=TRUE,PARAM, A.effects.key, D.effects.key, VARIANCE,
                           BETA.matrix){

                           #PAR$num.par,a.effects.key,aa.effects.key,D.effects.key,
                           #PAR$id.vector=PAR$id.vector,PAR$genenames.row=PAR$genenames.row,PAR$parnames=PAR$parnames,a.mean.correction=a.mean.correction,
                           #a.sd.correction=a.sd.correction,aa.mean.correction=aa.mean.correction,aa.sd.correction=aa.sd.correction,d.mean.correction=d.mean.correction,
                           #d.sd.correction=d.sd.correction,f.mean.correction=f.mean.correction,f.sd.correction=f.sd.correction,s.mean.correction=s.mean.correction,
                           #s.sd.correction=s.sd.correction,e.mean.correction=e.mean.correction,e.sd.correction=e.sd.correction,age.mean.correction=age.mean.correction,
                           #age.sd.correction=age.sd.correction,slope.mean.correction=slope.mean.correction,slope.sd.correction=slope.sd.correction,
                           #vt.model=vt.model,PAR$sibling.age.cor=PAR$sibling.age.cor,females.married=females.married,
                           #beta.matrix,PAR$am.multiplier){


############################################################################
#Create GENERATION CUR PHENOTYPES

#make size.gen
size.gen <- number.children <- dim(patgenes)[2]
number.genes <- dim(patgenes)[1]

#create ID's for each person & relatives two generations back
#note if max.size.id is big, it will be hard to see numbers in these matrices
Subj.ID <- sample((PARAM$id.vector[reltpe+1]+1):PARAM$id.vector[reltpe+2],size=size.gen,replace=FALSE)

#Name the rows & columns of genes.mate
genenames.col <- paste(rep(paste("subj",reltpe,".",sep=""),size.gen),Subj.ID,sep="")

#Stack the genes - one complete genotype per column
genes.stack <- rbind(patgenes,matgenes)
dimnames(genes.stack) <- list(PARAM$genenames.row,genenames.col)

#Make new effects.cur matrix
effects.cur <- matrix(0,nrow=PARAM$num.par,ncol=size.gen,dimnames=list(PARAM$parnames,NULL)) 
effects.cur["female",] <- (runif(size.gen)<.5)*1
if (reltpe==1) effects.cur["female",] <- rep(sample(c(0,1),size=size.gen/2,replace=TRUE),each=2)   #MZ twins are the same sex

Father.ID <- males.effects.mate["Subj.ID",]
Mother.ID <- females.effects.mate["Subj.ID",]
Fathers.Father.ID <- males.effects.mate["Father.ID",]
Fathers.Mother.ID <- males.effects.mate["Mother.ID",]
Mothers.Father.ID <- females.effects.mate["Father.ID",]
Mothers.Mother.ID <- females.effects.mate["Mother.ID",]
Spouse.ID <- rep(0,size.gen)
Relative.type <- rep(reltpe,size.gen) #7=twin/spouse parents; 8 on up are grandparents, greatgrandparents, etc...

effects.cur <- rbind(Subj.ID, Father.ID, Mother.ID, Fathers.Father.ID,Fathers.Mother.ID,
                     Mothers.Father.ID,Mothers.Mother.ID,Spouse.ID,Relative.type,genes.stack,effects.cur)
dimnames(effects.cur)[[2]] <- genenames.col
################################




### ADDITIVE GENETIC EFFECTS 
#Note: we need to correct the mean and SD of these effects with the *same* correction factor from generation 0, making
#all variance components relative to the mean=0 and var=1 in generation 1; if we standardized each generation, this would
#be wrong because the mean=0 and var=1 always, by definition
a.indx <- c(genes.stack)
a.matrix.prelim <- matrix(A.effects.key[a.indx,1],ncol=ncol(genes.stack))
a.prelim <- colSums(a.matrix.prelim)
if (founder.pop){corrections$a.mean <- mean(a.prelim);corrections$a.sd <- 1/sd(a.prelim)}
effects.cur["A",]<- (a.prelim-corrections$a.mean)*corrections$a.sd
################################




### ADDITIVE BY ADDITIVE EPISTATIC GENETIC EFFECTS 
aa.pp <-  genes.stack[1:number.genes,]*rbind(genes.stack[2:number.genes,],genes.stack[1,])
aa.mm <-  genes.stack[(number.genes+1):(number.genes*2),]*rbind(genes.stack[(number.genes+2):(number.genes*2),],genes.stack[(number.genes+1),])
aa.pm <- genes.stack[1:number.genes,]*rbind(genes.stack[(number.genes+2):(number.genes*2),],genes.stack[(number.genes+1),])
aa.mp <- rbind(genes.stack[2:number.genes,],genes.stack[1,])*genes.stack[(number.genes+1):(number.genes*2),]
aa.fin <- rbind(aa.pp,aa.mm,aa.pm,aa.mp)

aa.indx <- c(aa.fin %% 1e4)+1   #the "+1" is necessary to avoid zeros in the index
aa.matrix.prelim <- matrix(A.effects.key[aa.indx],nrow=number.genes*4)
aa.prelim <- colSums(aa.matrix.prelim)
if (founder.pop){corrections$aa.mean <- mean(aa.prelim);corrections$aa.sd <- 1/sd(aa.prelim)}
effects.cur["AA",] <- (aa.prelim-corrections$aa.mean)*corrections$aa.sd 
################################




### DOMINANCE GENETIC EFFECTS 
d.indx <- c((patgenes*matgenes) %% 1e4)+1  
d.matrix.prelim <- matrix(D.effects.key[d.indx],ncol=ncol(genes.stack))
d.prelim <- colSums(d.matrix.prelim)
if (founder.pop){corrections$d.mean <- mean(d.prelim);corrections$d.sd <- 1/sd(d.prelim)}
effects.cur["D",] <- (d.prelim-corrections$d.mean)*corrections$d.sd 
################################




### FAMILIAL ENVIRONMENTAL EFFECTS 
#Vertical Transmission;  transmission of parenting phenotypes to offspring;
f.prelim <- (males.effects.mate["parenting.phenotype",]+females.effects.mate["parenting.phenotype",])*sqrt(.5)
if (founder.pop) {corrections$f.sd <- 1}
corrections$f.mean <- mean(f.prelim)         #NOTE - this doesn't allow F to change means across generations
effects.cur["F",] <- (f.prelim-corrections$f.mean)*corrections$f.sd   
################################




### SIBLING ENVIRONMENT
s.prelim <- sqrt(1-VARIANCE$A.S.cor^2)*males.effects.mate["sib.env",] + VARIANCE$A.S.cor*effects.cur["A",]
if (founder.pop){corrections$s.mean <- mean(s.prelim);corrections$s.sd <- 1/sd(s.prelim)}
effects.cur["S",] <- (s.prelim-corrections$s.mean)*corrections$s.sd 
################################




### CHILDREN'S SIBLING ENVIRONMENT 
#this is a factor that allows the children of this individual to share a sib env. It does not affect the parenting phenotype
sib.prelim <- rnorm(size.gen)
if (founder.pop){corrections$sib.mean <- mean(sib.prelim); corrections$sib.sd <- 1/sd(sib.prelim)}
effects.cur["sib.env",] <- (sib.prelim-corrections$sib.mean)*corrections$sib.sd
################################




### UNIQUE ENVIRONMENTAL EFFECTS
u.prelim <- sqrt(1-VARIANCE$A.U.cor^2)*rnorm(size.gen) + VARIANCE$A.U.cor*effects.cur["A",]
if (founder.pop){corrections$u.mean <- mean(u.prelim);corrections$u.sd <- 1/sd(u.prelim)}
effects.cur["U",] <- (u.prelim-corrections$u.mean)*corrections$u.sd
################################




### MZ & TW ENVIRONMENTAL EFFECTS
if (reltpe==1){
  mz.prelim <- rep(rnorm(length(children.vector)),times=children.vector)
  tw.prelim <- rep(rnorm(length(children.vector)),times=children.vector)}
if (reltpe==2){
  mz.prelim <- rnorm(size.gen)
  tw.prelim <- rep(rnorm(length(children.vector)),times=children.vector)}
if (reltpe>2){
  mz.prelim <- rnorm(size.gen)
  tw.prelim <- rnorm(size.gen)}
if (founder.pop){corrections$mz.mean <- mean(mz.prelim);corrections$mz.sd <- 1/sd(mz.prelim)
               corrections$tw.mean <- mean(tw.prelim);corrections$tw.sd <- 1/sd(tw.prelim)}
effects.cur["MZ",] <- (mz.prelim-corrections$mz.mean)*corrections$mz.sd
effects.cur["TW",] <- (tw.prelim-corrections$tw.mean)*corrections$tw.sd

#the empirical value of MZ or TW (unstandardized); we'll use this later for summarizing (below, Finding Variance Components)
if (reltpe>2){empirical.tw <- 0}
if (reltpe==2){empirical.tw <- var(effects.cur["TW",])} 
if (reltpe==1){empirical.tw <- var(effects.cur["MZ",])} 
################################




### SEX EFFECTS
sex.prelim <- effects.cur["female",]
effects.cur["SEX",] <- (sex.prelim-.5)/.5        #always use population mean (.5) and sd of sex (.5)
if (center==FALSE) effects.cur["SEX",] <- effects.cur["SEX",]+.5    #this uncenters sex (range -.5 to 1.5), making possible scalar interactions with no main effects
################################




### AGE EFFECTS 
#make siblings ages s.t. they are clustered within families according to PARAM$sibling.age.cor
suppressWarnings(varcovar <- matrix(rep(c(1,rep(PARAM$sibling.age.cor,10)),9),nrow=10))
norm.age <- mvrnormal(females.married,mu=rep(0,10),Sigma=varcovar,empirical=FALSE) #normally dist age, s.t. sibs have intraclass correlation
cum.prob.age.spacing <- pnorm(norm.age)          #convert normally dist. ages to cumulative probability; i.e., 0 to 1 uniform
age.new <- qunif(cum.prob.age.spacing,PARAM$range.pop.age[1],PARAM$range.pop.age[2])  #find unif value at these probabilities

#for ancestors, parents, & children of twins
if (reltpe>3){
row.selector <- rep(1:females.married,time=children.vector)
age2 <- age.new[row.selector,]                   #matrix with #rows=number.children, repeated within families
age3 <- as.vector(t(age2))                       #vectorize age
adder1 <- rep(1:10,ceiling(number.children/10))[1:number.children] 
adder2 <- 0:(number.children-1)*10
age.fin <- age3[(adder1+adder2)]                 #choose ages
effects.cur["cur.age",] <- age.fin
effects.cur["age.at.mar",] <-  runif(size.gen,18,36)-(effects.cur["female",]*3) #makes females on avg 3 years younger @ marriage
if (founder.pop){corrections$age.mean <- mean(effects.cur["cur.age",]);corrections$age.sd <- 1/sd(effects.cur["cur.age",])}
effects.cur["AGE",] <- (effects.cur["cur.age",]-corrections$age.mean)*corrections$age.sd
effects.cur["AGE.M",] <- (effects.cur["age.at.mar",]-corrections$age.mean)*corrections$age.sd}

#for sibs
if (reltpe==3){
sibplustwin <- rep(1,females.married)+children.vector
num.tot <- sum(sibplustwin)
row.selector <- rep(1:females.married,time=sibplustwin)
age2 <- age.new[row.selector,]                   #matrix with #rows=number of sibs + twins, repeated within families
age3 <- as.vector(t(age2))                       #vectorize age
adder1 <- rep(1:10,ceiling(num.tot/10))[1:num.tot] 
adder2 <- 0:(num.tot-1)*10
age.fin <- age3[(adder1+adder2)]                 #choose ages for sibs + twins (to be carried forward)
sib.chooser <- rep(which(children.vector != 0),times=children.vector[children.vector != 0])+0:(number.children-1)
sib.age <- age.fin[sib.chooser]
twin.age <- age.fin[-sib.chooser] #twin age gets carried forward, to be input into make.phenotype for MZ and DZ twins
effects.cur["age.at.mar",] <- runif(size.gen,18,36)-(effects.cur["female",]*3) #makes females on avg 3 years younger @ marriage
effects.cur["cur.age",] <- sib.age
effects.cur["AGE",] <- (effects.cur["cur.age",]-corrections$age.mean)*corrections$age.sd
effects.cur["AGE.M",] <- (effects.cur["age.at.mar",]-corrections$age.mean)*corrections$age.sd}

#for dz & mz
if (reltpe<3){
tw.age <- rep(twin.age[which(children.vector == 2)],each=2)
effects.cur["age.at.mar",] <- runif(size.gen,18,36)-(effects.cur["female",]*3) #makes females on avg 3 years younger @ marriage
effects.cur["cur.age",] <-  tw.age
effects.cur["AGE",] <- (effects.cur["cur.age",]-corrections$age.mean)*corrections$age.sd
effects.cur["AGE.M",] <- (effects.cur["age.at.mar",]-corrections$age.mean)*corrections$age.sd}
################################





### SEX by ADDITIVE GENETIC EFFECTS 
sex.slope.prelim <- colSums(matrix(A.effects.key[a.indx,2],ncol=ncol(genes.stack)))   #individual gene effects for A,sex slopes

#correction factors to be used hereafter
if (founder.pop){corrections$sex.slope.mean <- mean(sex.slope.prelim);corrections$sex.slope.sd <- 1/sd(sex.slope.prelim)}
effects.cur["A.slopes.sex",] <- (sex.slope.prelim-corrections$sex.slope.mean)*corrections$sex.slope.sd
effects.cur["A.by.SEX",] <- effects.cur["A.slopes.sex",]*effects.cur["SEX",]
################################





### AGE by ADDITIVE GENETIC EFFECTS 
age.slope.prelim <- colSums(matrix(A.effects.key[a.indx,3],ncol=ncol(genes.stack)))   #individual gene effects for A,age slopes

#correction factors to be used hereafter
if (founder.pop){corrections$age.slope.mean <- mean(age.slope.prelim);corrections$age.slope.sd <- 1/sd(age.slope.prelim)}
effects.cur["A.slopes.age",] <- (age.slope.prelim-corrections$age.slope.mean)*corrections$age.slope.sd
effects.cur["A.by.AGE",] <- effects.cur["A.slopes.age",]*effects.cur["AGE",]
effects.cur["A.by.AGE.M",] <- effects.cur["A.slopes.age",]*effects.cur["AGE.M",]
################################





### S by ADDITIVE GENETIC EFFECTS 
S.slope.prelim <- colSums(matrix(A.effects.key[a.indx,4],ncol=ncol(genes.stack)))   #individual gene effects for A,S slopes

#correction factors to be used hereafter
if (founder.pop){corrections$S.slope.mean <- mean(S.slope.prelim);corrections$S.slope.sd <- 1/sd(S.slope.prelim)}
effects.cur["A.slopes.S",] <- (S.slope.prelim-corrections$S.slope.mean)*corrections$S.slope.sd
effects.cur["A.by.S",] <- effects.cur["A.slopes.S",]*effects.cur["S",]
################################





### U by ADDITIVE GENETIC EFFECTS 
U.slope.prelim <- colSums(matrix(A.effects.key[a.indx,5],ncol=ncol(genes.stack)))   #individual gene effects for A,U slopes

#correction factors to be used hereafter
if (founder.pop){corrections$U.slope.mean <- mean(U.slope.prelim);corrections$U.slope.sd <- 1/sd(U.slope.prelim)}
effects.cur["A.slopes.U",] <- (U.slope.prelim-corrections$U.slope.mean)*corrections$U.slope.sd
effects.cur["A.by.U",] <- effects.cur["A.slopes.U",]*u.prelim
################################





### CREATE PHENOTYPE
effects.cur["cur.phenotype",] <-    BETA.matrix %*%  effects.cur[PARAM$varnames,]
effects.cur["mating.phenotype",] <- BETA.matrix %*% (effects.cur[PARAM$varnames.M,]*PARAM$am.multiplier)
parenting.prelim <- as.vector(BETA.matrix %*% effects.cur[PARAM$varnames.M,])
if (founder.pop){corrections$parenting.mean <- mean(parenting.prelim); corrections$parenting.sd <- 1/sd(parenting.prelim)}
effects.cur["parenting.phenotype",] <- (parenting.prelim-corrections$parenting.mean)*corrections$parenting.sd  #we want parenting phenotype to start with mean=0, sd=1
################################



#make a list of what is returned from the function
list(effects.mate=effects.cur,size.mate=ncol(effects.cur),empirical.tw=empirical.tw,twin.age=twin.age,corrections=corrections)

#END FUNCTION
}

############################################################################