/*  


This file contains C code for the the infinite island model with no mutation in dioecious populations and 
the finite island model with k-allele mutation in s demes with dioecious populations.


*/


/* INFINITE ISLAND MODEL CODE STARTS HERE */

/* This program simulates the infinite island model with no mutation in dioecious populations */
/* Random number generator functions ignuin(), ranf() and gscgn() from the C library RANLIB (www.netlib.org/random/) are employed */
/* This program produces and prints estimates of Q1, Q2 and Q3 (in the array Qdat) which can then be used to estimate m and N */
/* No warranty, guarantee or technical support is associated with this code, users are responsible for their results */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "ranlib.h"

/* 
m =  migration rate
N = effective population size

reps = number of replicate simulations 
nL = number of loci 
s = number of demes sampled
n =  number of individuals genotyped per deme 
k =  number of alleles at each locus

nt = number of generations
*/

/* Main demographic parameters */
#define m 0.05
#define N 20

/* Sampling design */
#define reps 10
#define nL 20
#define n 50
#define k 8

/* Demographic scenario */
#define nt 1000
#define s 20

/* Number of random number generators in RANLIB */
#define ngen 32

int ikdraw(int,int);
 
int main()
{
  static long genlst[ngen]={0}; long igen;
  int iN1, iN2, ik; 
  int IndM[nL][s][N]={{{0}}}, IndF[nL][s][N]={{{0}}}, IndMo[nL][s][N]={{{0}}}, IndFo[nL][s][N]={{{0}}}, G[nL][s][2*n]={{{0}}}; 
  int qbit; double Q1e, p2s, Q1s, Q2s; 
  double P[s][k]={{0}}, Qa[nL][3]={{0}}, Qdat[reps][5]={{0}};
  int i0, i1, i2, i3, i4, i5, gct;
  FILE * fp;
   
  Q1e=((1-m)*(1-m))/2/N/(1-(1-m)*(1-m)*(1-1/2/N))+(1-(1-m)*(1-m))/(1-(1-m)*(1-m)*(1-1/2/N))/k; /* Equilibrium Q1 */
   
  srand(time(NULL)); /* Seed the generator from the standard C library */
  gct=1; 
  for (i0=1; i0<=ngen; i0++)
    { genlst[i0-1]=i0; }  
   
  /* Repeat for reps number of replicate simulations */
  for (i0=1; i0<=reps; i0++)
    {
	setall(rand(),time(NULL)+i0);  /* Seed the generators of RANLIB for each replicate*/
	
    /* Draw random initial genotypes */
	for (i1=1; i1<=s; i1++)
	  { for (i2=1; i2<=N; i2++)
		  { for (i3=1; i3<=nL; i3++)
		    { 
		    IndM[i3-1][i1-1][i2-1]=ignuin(1,k); IndMo[i3-1][i1-1][i2-1]=IndM[i3-1][i1-1][i2-1]; 
			IndF[i3-1][i1-1][i2-1]=ignuin(1,k); IndFo[i3-1][i1-1][i2-1]=IndF[i3-1][i1-1][i2-1]; 
		    }
		  } 
	  }
	   
    /* Iterate for nt-1 generations */
	if (nt>=2)
      {
	  for (i1=1; i1<=nt-1; i1++)
	    {
		/* Pick a random number generator for each generation */
		if (gct==ngen) {gct=1;} else {gct=gct+1;} 
		igen=*(genlst+gct-1); gscgn(1,&igen); 
		/* Iterate through demes */
		for (i2=1; i2<=s; i2++)
		  {
		  for (i3=1; i3<=N/2; i3++)
			{
			/* Males */
			if (ranf()<m)   /* Male is a migrant, generate genotype using Q1e */
			  { 
			  for (i4=1; i4<=nL; i4++)
				{ IndM[i4-1][i2-1][2*i3-1-1]=ignuin(1,k);
		        if (ranf()<Q1e) { IndM[i4-1][i2-1][2*i3-1]=IndM[i4-1][i2-1][2*i3-1-1]; } else { ik=ikdraw(k,IndM[i4-1][i2-1][2*i3-1-1]); IndM[i4-1][i2-1][2*i3-1]=ik; } 
				}
			  }	
		    else   /* Male is a resident, generate genotype by random mating */
			  { iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndM[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][i2-1][2*iN1-1-1+(rand()%2)]; IndM[i4-1][i2-1][2*i3-1]=IndFo[i4-1][i2-1][2*iN2-1-1+(rand()%2)]; 
				}
			  }
			/* Females */
			if (ranf()<m)   /* Female is a migrant, generate genotype using Q1e */
			  { 
			  for (i4=1; i4<=nL; i4++)
				{ IndF[i4-1][i2-1][2*i3-1-1]=ignuin(1,k);
		        if (ranf()<Q1e) { IndF[i4-1][i2-1][2*i3-1]=IndF[i4-1][i2-1][2*i3-1-1]; } else { ik=ikdraw(k,IndF[i4-1][i2-1][2*i3-1-1]); IndF[i4-1][i2-1][2*i3-1]=ik; } 
				}
			  }	
		    else   /* Female is a resident, generate genotype by random mating */
			  { iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndF[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][i2-1][2*iN1-1-1+(rand()%2)]; IndF[i4-1][i2-1][2*i3-1]=IndFo[i4-1][i2-1][2*iN2-1-1+(rand()%2)]; 
				}
			  }
			}  /* Individuals loop */
          }  /* Demes loop */
		  
		/* Update Indo */
		for (i2=1; i2<=s; i2++)
		  { for (i3=1; i3<=N; i3++)
			  { for (i4=1; i4<=nL; i4++)
				  { IndMo[i4-1][i2-1][i3-1]=IndM[i4-1][i2-1][i3-1]; IndFo[i4-1][i2-1][i3-1]=IndF[i4-1][i2-1][i3-1]; }
			  } 
		  }
		} /* Time loop */		
	  }
    else
	  { }
	 
	/* Sample n diploid genotypes before migration */
    for (i1=1; i1<=s; i1++)
	  { 	 
	  for (i2=1; i2<=n; i2++)
	    {
		/* Individual is a resident, generate genotype by random mating */
		iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
		for (i3=1; i3<=nL; i3++)
		  {
		  G[i3-1][i1-1][2*i2-1-1]=IndM[i3-1][i1-1][2*iN1-1-1+(rand()%2)]; G[i3-1][i1-1][2*i2-1]=IndF[i3-1][i1-1][2*iN2-1-1+(rand()%2)]; 
		  }
		}
      }
	
	/* Initialize Qa to zero */
	  for (i1=1; i1<=nL; i1++)
	    { Qa[i1-1][0]=0; Qa[i1-1][1]=0; Qa[i1-1][2]=0; }
	/* Calculate Q1, Q2, Q3 averaged over demes for each locus  */
    for (i1=1; i1<=nL; i1++)
	  { 	
	  /* Calculate Q1, Q2 averaged over demes */		
	  for (i2=1; i2<=s; i2++)
	    {
	    /* Initialize deme-specific variables to zero */
	    p2s=0; Q1s=0; Q2s=0;
	    for (i3=1; i3<=k; i3++)
	      { P[i2-1][i3-1]=0; }
	    for (i3=1; i3<=n; i3++)
	      {
		  /* Calculate Q1 */
		  if (G[i1-1][i2-1][2*i3-1-1]==G[i1-1][i2-1][2*i3-1]) {qbit=1;} else {qbit=0;} Q1s=Q1s+(double)qbit/n; 
		  /* Calculate P */
		  i4=G[i1-1][i2-1][2*i3-1-1]; i5=G[i1-1][i2-1][2*i3-1];
		  P[i2-1][i4-1]=P[i2-1][i4-1]+(double)1/2/n; P[i2-1][i5-1]=P[i2-1][i5-1]+(double)1/2/n;
	      }
	    for (i3=1; i3<=k; i3++)
	      { p2s=p2s+(double)(P[i2-1][i3-1]*P[i2-1][i3-1]); }
	    Q2s=n*p2s/(n-1)-(1+Q1s)/2/(n-1);
	    Qa[i1-1][0]=Qa[i1-1][0]+Q1s/s; Qa[i1-1][1]=Qa[i1-1][1]+Q2s/s;
	    }
      /* Calculate Q3 averaged over demes */	
	  for (i2=1; i2<=s-1; i2++)
	    {
	    for (i3=i2+1; i3<=s; i3++)
	      {
	      p2s=0;
	      for (i4=1; i4<=k; i4++)
	        {
		    p2s=p2s+(double)(P[i2-1][i4-1]*P[i3-1][i4-1]);
		    }
	      Qa[i1-1][2]=Qa[i1-1][2]+2*p2s/s/(s-1);	 
	      }
	    } 
	  }  /* Loop over loci ends */
   
   /* Average Q's over loci and save */
   Qdat[i0-1][0]=0; Qdat[i0-1][1]=0; Qdat[i0-1][2]=0; Qdat[i0-1][3]=0; Qdat[i0-1][4]=0;
   for (i1=1; i1<=nL; i1++)
     { Qdat[i0-1][2]=Qdat[i0-1][2]+Qa[i1-1][1-1]; Qdat[i0-1][3]=Qdat[i0-1][3]+Qa[i1-1][2-1]; Qdat[i0-1][4]=Qdat[i0-1][4]+Qa[i1-1][3-1]; }
   Qdat[i0-1][2]=Qdat[i0-1][2]/nL; Qdat[i0-1][3]=Qdat[i0-1][3]/nL; Qdat[i0-1][4]=Qdat[i0-1][4]/nL;
   Qdat[i0-1][0]=i0; Qdat[i0-1][1]=n;
   
   } /* Loop over reps ends */
      
/* Print results to a text file tab-delimited */	  
/* The rows of Qdat contain i0, n, Q1, Q2, Q3 for each replicate simulation */
/* The data in Qdat can be used to estimate m and N */
fp = fopen("diIIM_OUTPUT.txt", "w");  
for (i0=1; i0<=reps; i0++)
  { for (i1=1; i1<=4; i1++)
      { fprintf(fp,"%.6f\t",Qdat[i0-1][i1-1]); }
      fprintf(fp,"%.6f\n",Qdat[i0-1][5-1]);
  }
fclose(fp);			
/* Print results from 1st 10 reps to the screen tab-delimited */ 
for (i0=1; i0<=10; i0++)
  { for (i1=1; i1<=4; i1++)
      { printf("%.6f\t",Qdat[i0-1][i1-1]); }
      printf("%.6f\n",Qdat[i0-1][5-1]);
  }		

return 0;
} /* main ends */

/* Function ikdraw randomly selects an allele that is not inot */
int ikdraw(int numk, int inotk)
{
int j0k, ctk, valk; int nsetk[k-1]={0};
ctk=1;
  for (j0k=1; j0k<=numk; j0k++)
  {
    if (j0k!=inotk)
    {nsetk[ctk-1]=j0k; ctk=ctk+1;}
	else 
	{ }
  }
  valk=nsetk[ignuin(0,numk-2)];
return valk;
}

/* INFINITE ISLAND MODEL CODE ENDS HERE */




/* FINITE ISLAND MODEL CODE STARTS HERE */

/* This program simulates the finite island model with k-allele mutation in s demes with dioecious populations */
/* Random number generator functions ignuin(), ranf() and gscgn() from the C library RANLIB (www.netlib.org/random/) are employed */
/* This program produces and prints estimates of Q1, Q2 and Q3 (in the array Qdat) which can then be used to estimate u, m and N */
/* No warranty, guarantee or technical support is associated with this code, users are responsible for their results */
	
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "ranlib.h"
 
/* 
u = mutation rate
m =  migration rate
N = effective population size 

reps = number of replicate simulations 
nL = number of loci 
n =  number of individuals genotyped per deme 
k =  number of alleles at each locus 

nt = number of generations
s = number of demes 
*/

/* Main demographic parameters */
#define u 0.0005
#define m 0.05
#define N 20

/* Sampling design */
#define reps 10
#define nL 30
#define n 100
#define k 8

/* Demographic scenario */
#define nt 6000
#define s 20

/* Number of generators */
#define ngen 32

int isdraw(int,int); int ikdraw(int,int);
 
int main()
{
  static long genlst[ngen]={0}; long igen;
  int iN1, iN2, ik, is; 
  int IndM[nL][s][N]={{{0}}}, IndF[nL][s][N]={{{0}}}, IndMo[nL][s][N]={{{0}}}, IndFo[nL][s][N]={{{0}}}, G[nL][s][2*n]={{{0}}}; 
  int qbit; double p2s, Q1s, Q2s; 
  double P[s][k]={{0}}, Qa[nL][3]={{0}}, Qdat[reps][5]={{0}};
  int i0, i1, i2, i3, i4, i5, gct;
  FILE * fp;
   
  srand(time(NULL)); /* Seed the generator from the standard C library */
  gct=1; 
  for (i0=1; i0<=ngen; i0++)
    { genlst[i0-1]=i0; }  
   
  /* Repeat for reps number of replicate simulations */
  for (i0=1; i0<=reps; i0++)
    {
	setall(rand(),time(NULL)+i0);  /* Seed the generators of RANLIB for each replicate*/
	
    /* Draw random initial genotypes */
	for (i1=1; i1<=s; i1++)
	  { for (i2=1; i2<=N; i2++)
		  { for (i3=1; i3<=nL; i3++)
		    { 
		    IndM[i3-1][i1-1][i2-1]=ignuin(1,k); IndMo[i3-1][i1-1][i2-1]=IndM[i3-1][i1-1][i2-1]; 
			IndF[i3-1][i1-1][i2-1]=ignuin(1,k); IndFo[i3-1][i1-1][i2-1]=IndF[i3-1][i1-1][i2-1]; 
		    }
		  } 
	  }
	   
    /* Iterate for nt-1 generations */
	if (nt>=2)
      {
	  for (i1=1; i1<=nt-1; i1++)
	    {
		/* Pick a random number generator for each generation */
		if (gct==ngen) {gct=1;} else {gct=gct+1;} 
		igen=*(genlst+gct-1); gscgn(1,&igen); 
		/* Iterate through demes */
		for (i2=1; i2<=s; i2++)
		  {
		  for (i3=1; i3<=N/2; i3++)
			{
			/* Males */
			if (ranf()<m)   /* Male is a migrant, generate genotype by random mating */
			  { is=isdraw(s,i2);	
				iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndM[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][is-1][2*iN1-1-1+(rand()%2)]; IndM[i4-1][i2-1][2*i3-1]=IndFo[i4-1][is-1][2*iN2-1-1+(rand()%2)];
				/* 1st gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndM[i4-1][i2-1][2*i3-1-1]); IndM[i4-1][i2-1][2*i3-1-1]=ik; }
			    else			 
			      { }
			    /* 2nd gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndM[i4-1][i2-1][2*i3-1]); IndM[i4-1][i2-1][2*i3-1]=ik; }
			    else			 
			      { }  
				}
			  }	
		    else   /* Male is a resident, generate genotype by random mating */
			  { iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndM[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][i2-1][2*iN1-1-1+(rand()%2)]; IndM[i4-1][i2-1][2*i3-1]=IndFo[i4-1][i2-1][2*iN2-1-1+(rand()%2)]; 
				/* 1st gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndM[i4-1][i2-1][2*i3-1-1]); IndM[i4-1][i2-1][2*i3-1-1]=ik; }
			    else			 
			      { }
			    /* 2nd gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndM[i4-1][i2-1][2*i3-1]); IndM[i4-1][i2-1][2*i3-1]=ik; }
			    else			 
			      { }  
				}
			  }
			/* Females */
			if (ranf()<m)   /* Female is a migrant, generate genotype by random mating */
			  { is=isdraw(s,i2);	
				iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndF[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][is-1][2*iN1-1-1+(rand()%2)]; IndF[i4-1][i2-1][2*i3-1]=IndFo[i4-1][is-1][2*iN2-1-1+(rand()%2)]; 
				/* 1st gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndF[i4-1][i2-1][2*i3-1-1]); IndF[i4-1][i2-1][2*i3-1-1]=ik; }
			    else			 
			      { }
			    /* 2nd gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndF[i4-1][i2-1][2*i3-1]); IndF[i4-1][i2-1][2*i3-1]=ik; }
				else			 
			      { }
				}
			  }	
		    else   /* Female is a resident, generate genotype by random mating */
			  { iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
			  for (i4=1; i4<=nL; i4++)
				{
				IndF[i4-1][i2-1][2*i3-1-1]=IndMo[i4-1][i2-1][2*iN1-1-1+(rand()%2)]; IndF[i4-1][i2-1][2*i3-1]=IndFo[i4-1][i2-1][2*iN2-1-1+(rand()%2)]; 
				/* 1st gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndF[i4-1][i2-1][2*i3-1-1]); IndF[i4-1][i2-1][2*i3-1-1]=ik; }
			    else			 
			      { }
			    /* 2nd gene might mutate */
			    if (ranf()<u)  
			      { ik=ikdraw(k,IndF[i4-1][i2-1][2*i3-1]); IndF[i4-1][i2-1][2*i3-1]=ik; }
				else			 
			      { } 
				}
			  }
			}  /* Individuals loop */
          }  /* Demes loop */
		  
		/* Update Indo */
		for (i2=1; i2<=s; i2++)
		  { for (i3=1; i3<=N; i3++)
			  { for (i4=1; i4<=nL; i4++)
				  { IndMo[i4-1][i2-1][i3-1]=IndM[i4-1][i2-1][i3-1]; IndFo[i4-1][i2-1][i3-1]=IndF[i4-1][i2-1][i3-1];}
			  } 
		  }
		} /* Time loop */		
	  }
    else
	  { }
	 
	/* Sample n diploid genotypes before migration */
    for (i1=1; i1<=s; i1++)
	  { 	 
	  for (i2=1; i2<=n; i2++)
	    {
		/* Individual is a resident, generate genotype by random mating */
		iN1=ignuin(1,N/2); iN2=ignuin(1,N/2);
		for (i3=1; i3<=nL; i3++)
		  {
		  G[i3-1][i1-1][2*i2-1-1]=IndM[i3-1][i1-1][2*iN1-1-1+(rand()%2)]; G[i3-1][i1-1][2*i2-1]=IndF[i3-1][i1-1][2*iN2-1-1+(rand()%2)]; 
		  /* 1st gene might mutate */
		  if (ranf()<u)  
		    { ik=ikdraw(k,G[i3-1][i1-1][2*i2-1-1]); G[i3-1][i1-1][2*i2-1-1]=ik; }
		  else			 
			{ }
		  /* 2nd gene might mutate */
		  if (ranf()<u)  
		    { ik=ikdraw(k,G[i3-1][i1-1][2*i2-1]); G[i3-1][i1-1][2*i2-1]=ik; }
		  else			 
			{ } 
		  }
		}
      }
	
	/* Initialize Qa to zero */
	for (i1=1; i1<=nL; i1++)
	  { Qa[i1-1][0]=0; Qa[i1-1][1]=0; Qa[i1-1][2]=0; }
	/* Calculate Q1, Q2, Q3 averaged over demes for each locus  */
    for (i1=1; i1<=nL; i1++)
	  { 	
	  /* Calculate Q1, Q2 averaged over demes */		
	  for (i2=1; i2<=s; i2++)
	    {
	    /* Initialize deme-specific variables to zero */
	    p2s=0; Q1s=0; Q2s=0;
	    for (i3=1; i3<=k; i3++)
	      { P[i2-1][i3-1]=0; }
	    for (i3=1; i3<=n; i3++)
	      {
		  /* Calculate Q1 */
		  if (G[i1-1][i2-1][2*i3-1-1]==G[i1-1][i2-1][2*i3-1]) {qbit=1;} else {qbit=0;} Q1s=Q1s+(double)qbit/n; 
		  /* Calculate P */
		  i4=G[i1-1][i2-1][2*i3-1-1]; i5=G[i1-1][i2-1][2*i3-1];
		  P[i2-1][i4-1]=P[i2-1][i4-1]+(double)1/2/n; P[i2-1][i5-1]=P[i2-1][i5-1]+(double)1/2/n;
	      }
	    for (i3=1; i3<=k; i3++)
	      { p2s=p2s+(double)(P[i2-1][i3-1]*P[i2-1][i3-1]); }
	    Q2s=n*p2s/(n-1)-(1+Q1s)/2/(n-1);
	    Qa[i1-1][0]=Qa[i1-1][0]+Q1s/s; Qa[i1-1][1]=Qa[i1-1][1]+Q2s/s;
	    }
      /* Calculate Q3 averaged over demes */	
	  for (i2=1; i2<=s-1; i2++)
	    {
	    for (i3=i2+1; i3<=s; i3++)
	      {
	      p2s=0;
	      for (i4=1; i4<=k; i4++)
	        {
		    p2s=p2s+(double)(P[i2-1][i4-1]*P[i3-1][i4-1]);
		    }
	      Qa[i1-1][2]=Qa[i1-1][2]+2*p2s/s/(s-1);	 
	      }
	    } 
	  }  /* Loop over loci ends */
   
    /* Average Q's over loci and save */
    Qdat[i0-1][0]=0; Qdat[i0-1][1]=0; Qdat[i0-1][2]=0; Qdat[i0-1][3]=0; Qdat[i0-1][4]=0;
    for (i1=1; i1<=nL; i1++)
      { Qdat[i0-1][2]=Qdat[i0-1][2]+Qa[i1-1][1-1]; Qdat[i0-1][3]=Qdat[i0-1][3]+Qa[i1-1][2-1]; Qdat[i0-1][4]=Qdat[i0-1][4]+Qa[i1-1][3-1]; }
    Qdat[i0-1][2]=Qdat[i0-1][2]/nL; Qdat[i0-1][3]=Qdat[i0-1][3]/nL; Qdat[i0-1][4]=Qdat[i0-1][4]/nL;
    Qdat[i0-1][0]=i0; Qdat[i0-1][1]=n;
   
    } /* Loop over reps ends */
      
/* Print results to a text file tab-delimited */	  
/* The rows of Qdat contain i0, n, Q1, Q2, Q3 for each replicate simulation */
/* The data in Qdat can be used to estimate m and N */
fp = fopen("diFIM_OUTPUT.txt", "w");  
for (i0=1; i0<=reps; i0++)
  { for (i1=1; i1<=4; i1++)
      { fprintf(fp,"%.6f\t",Qdat[i0-1][i1-1]); }
      fprintf(fp,"%.6f\n",Qdat[i0-1][5-1]);
  }
fclose(fp);			
/* Print results from 1st 10 reps to the screen tab-delimited */ 
for (i0=1; i0<=10; i0++)
  { for (i1=1; i1<=4; i1++)
      { printf("%.6f\t",Qdat[i0-1][i1-1]); }
      printf("%.6f\n",Qdat[i0-1][5-1]);
  }		

return 0;
} /* main ends */

/* Function isdraw randomly selects a deme that is not inot */
int isdraw(int nums, int inots)
{
int j0s, cts, vals; int nsets[s-1]={0};
cts=1;
  for (j0s=1; j0s<=nums; j0s++)
  {
    if (j0s!=inots)
    {nsets[cts-1]=j0s; cts=cts+1;}
	else 
	{ }
  }
  vals=nsets[ignuin(0,nums-2)];
return vals;
}
/* Function ikdraw randomly selects an allele that is not inot */
int ikdraw(int numk, int inotk)
{
int j0k, ctk, valk; int nsetk[k-1]={0};
ctk=1;
  for (j0k=1; j0k<=numk; j0k++)
  {
    if (j0k!=inotk)
    {nsetk[ctk-1]=j0k; ctk=ctk+1;}
	else 
	{ }
  }
  valk=nsetk[ignuin(0,numk-2)];
return valk;
}

/* FINITE ISLAND MODEL CODE ENDS HERE */