---
title: "fishSim-vignette"
author: "Shane M Baylis"
date: "`r Sys.Date()`"
# output: rmarkdown::html_vignette
output:
    rmarkdown::pdf_document:
       number_sections: true
vignette: >
  %\VignetteIndexEntry{fishSim-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(mc.cores = 2)
```

fishSim is a package for demographic simulation and exploration of kin
relationships. It includes functions to build populations and move,
mate, kill, age, and switch the sex of individuals within those
populations or defined subpopulations. It _doesn't_ include any
one-line functions that define and run an entire simulation. Instead,
users will have to build their own for-loop over years, with each year
containing mating, mortality, ageing, sampling, and/or movement events
as needed. With a little planning, this makes fishSim _very_
adaptable. Importantly, fishSim is individual-based and retains
parentage information for all individuals, allowing full
reconstruction of kin relationships to arbitrary depth.

In addition to population simulation features, fishSim includes a set
of convenience functions allowing the user to check the scenario
population growth rate, find simulation scenarios with null population
growth (`PoNG()`), selectively archive subsets of the population for
increased computational speed. For a given population, fishSim also
includes functions to look up pairs of animals, compare their sets of
ancestors, and classify them into kin categories based on their
nearest shared ancestor(s), or show patterns of shared ancestors for
members of a pair.

fishSim is developed in response to the demands of close-kin
mark-recapture projects. 'typical' use-cases therefore involve an
understanding of both demography and kin relationships. This vignette
will cover two use scenarios - one relatively simple one, with:

- age-specific (but not sex-specific) fecundity,
- constant adult survival,
- a single population with no subdivisions (i.e., one 'stock'),
- no sex-switching, and
- first-year survival set by `PoNG()`, giving a flat population size through time,
- one-shot sampling of animals from the population

...and a one designed to show off the bells and whistles, with:

- age-specific, sex-specific fecundity with within-season mating exhaustion for females,
- age-specific, sex-specific survival,
- three subpopulations ('stocks') with markovian movement between subpopulations,
- male-to-female sex-switching,
- a growing population,
- lethal sampling, occurring through time during the simulation.

In the 'bells-and-whistles' scenario, we will also cover archiving
tools. For both scenarios, we will take a brief look at relationship
patterns between pairs of animals.

\newpage

# The relatively-simple scenario

## Setup

In order to simulate a population, we need an object to hold that
population's data, and the population will need some founding
members. Our data object will have one line per individual, so let's
call it `indiv`.

```{r}
library(fishSim)
indiv <- makeFounders(stocks = c(1))
head(indiv)
```

Now, what's here? `indiv` has one line per individual, and each
individual has an ID code (`indiv$Me`), a sex (`indiv$Sex`), a
father's ID code (`indiv$Dad` - but which is 'founder' here, because
founding animals can't really be said to have parents), a mother's ID
code (`indiv$Mum` - ditto), a birth year `indiv$BirthY`, a death year
(`indiv$DeathY` - NA for any animals that are still alive), a stock
membership (`indiv$Stock`), an age in years (`indiv$AgeLast`), and a flag
indicating whether the individual has been sampled (`indiv$SampY`).

The age structure, sex ratio, and stock memberships of founders can
all be specified in the `makeFounders()` call. We've only changed one
of the default values - instead of the default 3 stocks, we're going
to have a single panmictic founder population. The other defaults give
us 1000 founders, with an even sex ratio, with an age-structure that
implies 70% annual survival, but with a hard age-limit of 20 years.

## Mating

Having set up our population, we can now get its members to breed, get
older, switch sex and die, using `mate()` (or `altMate()` or
`bondedMate()`), `birthdays()`, `sexSwitch()`, and `mort()`
functions. If we had more than one stock, we could also get
individuals to move between stocks with markovian movement
probabilities, using `move()`. Let's look at those.

```{r}
nrow(indiv)
indiv <- mate(indiv, year = 1)
nrow(indiv) ## 200 newborns
tail(indiv) ## newborns are added to the end
```

Here are some newborns, generated using the default `mate()`
settings. Unlike founders, these newborns all have a father's ID and a
mother's ID, are all aged 0, and all have a birth year of equal to the
`year` argument. There are 200 newborns, because `mate()` generates
new members as a proportion of the extant population size (i.e., we
specify fecundity for the population, not the individuals), and that
proportion (argument `fecundity`) is 0.2 by default. Each individual
was born in a clutch, and each clutch contained a Poisson-distributed
number of clutch-members, the default mean ('batchSize') of
0.5. `mate()` simply kept producing newborns by mating together random
male/female pairs until `nrow(indiv) * 0.2` new offspring were
produced.

Alternatively, we can mate our population using `bondedMate()`, which
stores a mate ID for mature animals and generates offspring within
each pair unless one member dies or moves to another stock, or the
pair divorces. Re-mating occurs among surviving animals, but this
system can give repeated breeding between pair members for many
breeding seasons.

Finally, we can mate our population using `altMate()`, which specifies
individual maturities by age, and a probability distribution of number
of offspring per mature female. Let's generate some new founders and
try that.

```{r}
indiv <- makeFounders(stocks = c(1))
nrow(indiv) ## 1000 founders
indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = 1)
tail(indiv)
```

In this call to `altMate()`, females became sexually mature at 2 years
of age (`firstBreed = 2`), and each sexually-mature female produced a
Poisson-distributed number of offspring with mean set by
`batchSize`. It is possible to specify age-specific or sex-specific
fecundities, set paternity within-clutch to be single or multiple, to
exhaust fathers within each breeding season (females, in this system,
always breed to exhaustion each season), and to set a sex ratio for
offspring by additional arguments to `altMate()`, and these options
can also be set in `mate()` and `bondedMate()`.

## Mortality

Now, let's kill some of the population. Mortality probabilities can be
flat, age-specific or stock-specific (and with a little extra effort,
sex- or age:stock:sex-specific too), or we can randomly kill animals
such that the population is reduced to a certain size, and we can set
an age past which no animal will survive. For now, let's set a flat
20% mortality rate (i.e., the probability of death is equal for all
animals).

```{r}
nrow(indiv)
indiv <- mort(indiv, year = 1, type = "flat", mortRate = 0.2)
nrow(indiv)
head(indiv, n = 15)
```

Mortality does not remove dead animals from 'indiv' - we will need to
refer to them later - it just updates their death year to the value
given in `year`. In general, though, animals with a non-NA value for
death year will not move, mate, switch sex, or have birthdays.

## Birthdays

On the topic of birthdays: that's the final thing that must be done
before we can put this all together and run a full (albeit basic)
demographic simulation: we need to be able to increase the age of our
animals, and that is what the `birthdays()` function is for. It's very
simple: it takes all living members of the population, and increments
their age by 1.

```{r}
tail(indiv)
indiv <- birthdays(indiv)
tail(indiv)
```

## Sampling from the population

It is often useful to keep track of which animals in the population
have been captured, and if capture is lethal, to mark captured animals
as dead through a separate process from the normal 'mort'
functions. The function `capture()` takes a sample of individuals,
marks them as 'captured' by updating their 'captured' value, and, if
capture is lethal, marks them as dead by updating their death year, as
in `mort()`. It is possible to make `capture()` sex-specific,
age-specific, or age:sex specific, so that only one sex may be
captured, only one age-class may be captured, or only one sex at one
age-class may be captured, respectively.

```{r}
## non-lethal sampling of 5 females
indiv <- capture(indiv, n = 5, year = 1, fatal = FALSE, sex = "F")
## lethal sampling of 5 animals, from either sex
indiv <- capture(indiv, n = 5, year = 1, fatal = TRUE)
## non-lethal sampling of 5 one-year-old females
indiv <- capture(indiv, n = 5, year = 1, fatal = FALSE, sex = "F", age = 1)
```

## What's my population doing? Can I make it stay the same size?

A couple more functions deserve a mention here: `check_growthrate()`
and `PoNG()`. `check_growthrate()` tells you how quickly your
population is growing. `PoNG()` tells you what you need first-year
survival to be, in order for your population to maintain a constant
size (within limits - the estimation is Leslie Matrix-based, and there
are ways to fool Leslie Matrices that are out of scope for this
vignette but covered briefly in the `PoNG()` and `check_growthrate()`
documentation). Let's try both of those, using the mating and
mortality settings we have already used.

```{r}
check_growthrate(mateType = "flat", mortType = "flat",
                 batchSize = 0.9, firstBreed = 2, mortRate = 0.2)

## with the current settings, we expect our population to grow by about 1.8% annually.

PoNG(mateType = "flat", mortType = "flat", batchSize = 0.9, firstBreed = 2,
     mortRate = 0.2)

## if first-year mortality was about 0.445, rather than 0.2, our population
## would have null growth. You can also read off a range of possible growth rates
## from the plot.
```

## Turning processes into a simple simulation

To turn those processes into a full demographic simulation, all that
is needed is to repeat the processes in a loop. Let's do that, but set
age-specific mortality so that the long-run average population growth
rate is zero using the first-year mortality rate we got from `PoNG()`,
so that our population is unlikely to explode or become extinct.

```{r}
indiv <- makeFounders(stocks = c(1))
ageMort <- c(0.445, rep(0.2, 100)) ## age-specific mortality is 0.445 for first-years,
## 0.2 for all older age-classes.

for (y in 1:60) {
    indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y) ## y for year
    indiv <- mort(indiv, year = y, type = "age", ageMort = ageMort)    ## age-specific mort
    indiv <- birthdays(indiv)
}
tail(indiv) ## a population with 60 years of births, deaths, and birthdays
nrow(indiv[is.na(indiv$DeathY),]) ## the currently-alive population size. Note that population
                               ## growth only *averages* zero, and variability occurs!
```

## Looking up relationships between pairs of animals

One of the key advantages of fishSim is its ability to report on
different relationships between pairs of animals in a sample. Are
these two each other's siblings? Half-siblings? Is one the parent of
the other? The grandparent? Some obscure half-cousin, once removed?

Internally, these relationships are stored in terms of 'shared
ancestors at the n'th generation of each member', starting with the
'self' as generation 1. So a pair might have a shared `pairs$TwoThree`
ancestor - that ancestor is the parent of one member, and the
grandparent of the other. Or the pair might have a shared
`pairs$OneTwo` ancestor - in which case, one member's _self_ is the
other member's parent.

If a pair shares a parent (i.e., `pairs$TwoTwo == 1`), then they also
share at least two grandparents (`pairs$ThreeThree >= 2`), four
great-grandparents (`pairs$FourFour >= 4`), and so on. This pattern
holds generally: a half thiatic pair (`pairs$TwoThree == 1`) will also
share at least two `pairs$ThreeFour >= 2` ancestors and at least four
`pairs$FourFive >= 4` ancestors. It is also possible, for example, for
a pair to share two grandparents, but no parents - that is the case
for full cousin pairs.

```{r}
## mark 200 individuals alive at the end of the simulation as 'captured'.
indiv <- capture(indiv, n = 200, year = 60, fatal = FALSE)

## look up each animal's ancestors and look for shared ancestors between each
## pair of sampled animals:
pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2)

POPs <- pairs[pairs$OneTwo == 1,] ## Parent-Offspring pairs
GGPs <- pairs[pairs$OneThree == 1,] ## Grandparent-Grandoffspring pairs
HSPs <- pairs[pairs$TwoTwo == 1,] ## Half-sibling pairs
FSPs <- pairs[pairs$TwoTwo == 2,] ## Full Sibling pairs (self-comparisons
## are automatically excluded)
FCPs <- pairs[pairs$ThreeThree == 2 & pairs$TwoTwo != 1,] ## Full Cousin pairs

## look at the number of shared ancestors at each ancestral
## generation, for one of the half-sibling pairs.
lookAtPair(HSPs[1,])

relatives <- namedRelatives(pairs) ## shows the number of pairs of each relationship type
relatives
```

\newpage

# The 'bells and whistles' scenario

In this second scenario, we will simulate a metapopulation with:

- age-specific, sex-specific fecundity with within-season mating exhaustion for females,
- stock-specific, age-specific survival,
- three subpopulations ('stocks') with markovian movement between subpopulations,
- male-to-female sex-switching,
- a growing population.
- archiving for speed (?)

We set up the population basically as before, with a couple of slight
tweaks. First, we should make sure the founding population has
multiple stocks (we'll give it three). Second, we will set up an
archive matrix. The archive matrix is intended to hold the records of
dead animals in large, long-running simulations: dead animals do not
take part in any further mating, movement, sex-switching, aging, or
mortality events, but if they stay in `indiv`, `indiv` can become huge
and unweildy, slowing down all of those processes. There is of course
a trade-off, in that writing dead animals to the archive takes system
time, so the optimum may be to only archive dead animals once every
few 'years'.

## Setting up the data objects

```{r}
## set up founders with three stocks: two that each contain 30% of the population,
## and one that contains the remaining 40%.

indiv <- makeFounders(pop = 1000, stocks = c(0.3, 0.3, 0.4))

## set up archive - just a matrix with zero rows and eight columns

archive <- make_archive()

```

## Parameterising movement, survival, and maturity structures

Because we will have inter-stock movement in this sim, we will also
need to set up a matrix giving the probability that an individual will
move into another stock, given its current stock membership. Because
survival will be age- and stock-dependent, we will need a matrix of
survival rates with as many columns as stocks, and as many rows as (at
least) the age of the oldest plausible population-member. Because we
will have age-specific, sex-specific fecundity, we will need a male
age-specific maturity curve and a female age-specific maturity
curve. Let's set those up now.

```{r}
## Markovian movement matrix
stocks <- c(0.3, 0.3, 0.4)
admix.m <- matrix(NA, nrow = length(stocks), ncol = length(stocks))
for(i in 1:nrow(admix.m)) {
    admix.m[i,] <- stocks*stocks[i]
}
## admix.m shows movement proportional to starting population sizes.
admix.m
## let's tweak those numbers so that animals tend to stay where
## they are, and not move around so much.
admix.m <- matrix(c(0.23, 0.03, 0.04, 0.03, 0.23, 0.04, 0.04, 0.04, 0.32),
                  nrow = length(stocks), ncol = length(stocks), byrow = FALSE)
admix.m

## Age- and stock-dependent survival
ageStockMort <- matrix(c(0.47, 0.37, 0.27, rep(0.23, 97),
                         0.45, 0.35, 0.25, 0.20, rep(0.22, 96),
			 0.45, 0.3, 0.3, 0.19, rep(0.2, 96)),
                       ncol = length(stocks), nrow = 100)
head(ageStockMort)

## Sex-specific maturity curves

maleCurve <- c(0,0,0,0.1,0.5,0.8,0.85,0.9,0.95,rep(1, 91))
femaleCurve <- c(0,0,0.5,0.9,0.95,rep(1,95))
## maleCurve and femaleCurve should both be long enough that no individuals
## will outlive the curve.
head(maleCurve)
head(femaleCurve)

```

## Checking population growth rates with the new parameters

Now, what's going to happen with these subpopulations, if we run with those parameters?

```{r}
check_growthrate(mateType = "ageSex", mortType = "ageStock", batchSize = 1.6,
                 femaleCurve = femaleCurve,
                 ageStockMort = ageStockMort)
## Not bad. Two of the three populations are increasing. The third will probably be kept viable
## by immigration from the other two. Note that I fiddled 'batchSize' (which is the mean number
## of offspring per mature female per breeding attempt) a bit.
```

## Running the simulation

Now we just run the population simulator, as we did for the simple
population, but with a couple of extra things to specify. Let's give
it 100 years.

```{r}
for (k in 1:100) {
    ## very rarely, switch some males to females
    indiv <- sexSwitch(indiv = indiv, direction = "MF", prob = 1e-04)
    ## move animals according to the markovian matrix we set up before
    indiv <- move(indiv = indiv, moveMat = admix.m)
    ## mate animals using the age-specific, sex-specific curves we set up before
    indiv <- altMate(indiv = indiv, batchSize = 1.6, type = "ageSex", maleCurve = maleCurve,
                     femaleCurve = femaleCurve, year = k)
    ## kill animals on the basis of their age and stock, as set up before
    indiv <- mort(indiv = indiv, year = k, type = "ageStock", ageStockMort = ageStockMort)
    if(k %% 10 == 0) {
        archive <- archive_dead(indiv = indiv, archive = archive)
	    indiv <- remove_dead(indiv = indiv)
	    cat( sprintf( '\rIteration %i complete', k))
	## occasionally sample from the population and do some clean-up. And report on progress.
    }
    if(k %in% c(94:100)) {
	indiv <- capture(indiv, n = 40, fatal = TRUE, year = k)
    }
    indiv <- birthdays(indiv = indiv)
}
archive <- rbind(archive, indiv)
indiv <- archive  ## merge 'indiv' and 'archive', since they were only separated for speed.

```

## Looking up relationships between pairs of animals

Now that the simulation is finished, examining kin relationships is
the nearly same as the simple scenario. The only difference comes from
the fact that we we integrated multiple years of sample-capture into
the simulation, so we don't need to call `capture()` as a separate
operation.

```{r}
## look up each animal's ancestors and look for shared ancestors between each
## pair of sampled animals:
pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2)

POPs <- pairs[pairs$OneTwo == 1,] ## Parent-Offspring pairs
GGPs <- pairs[pairs$OneThree == 1,] ## Grandparent-Grandoffspring pairs
HSPs <- pairs[pairs$TwoTwo == 1,] ## Half-sibling pairs
FSPs <- pairs[pairs$TwoTwo == 2,] ## Full Sibling pairs (self-comparisons
## are automatically excluded)
FCPs <- pairs[pairs$ThreeThree == 2 & pairs$TwoTwo != 1,] ## Full Cousin pairs

## look at the number of shared ancestors at each ancestral
## generation, for one of the half-sibling pairs.
lookAtPair(HSPs[1,])

relatives <- namedRelatives(pairs) ## shows the number of pairs of each relationship type
relatives
```

\newpage

# Addenda

## A note on mark-dependent survival or movement

It is possible to simulate mark-dependent survival or movement by
working on subsets of the 'indiv' matrix within the main simulation
loop, and applying different mortality or movement operations to each
subset. Note that if `mate()`, `altMate()`, `birthdays()`, etc., are
run on a subset of the population, the operation will only affect
members of that subset. In the case of subsets for marked and unmarked
individuals, this would imply that marked animals only mate with other
marked animals, and _vice versa_. It is important to ensure that a
combined `indiv` object is available for any operations that affect
the whole population, or to repeat each opertation for each subset.

```{r, eval = FALSE}
indiv <- makeFounders(stocks = c(1))
ageMort <- c(0.305, rep(0.2, 100)) ## age-specific mortality is 0.305 for first-years,
## 0.2 for all older age-classes.

for (y in 1:60) {
    indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y)
    ## 'mate' happens first, affecting the whole population

    indiv <- capture(indiv, n = 5, year = y, fatal = FALSE)
    indiv_unmarked <- indiv[is.na(indiv$SampY),] ## all the unmarked animals
    indiv_marked <- indiv[!is.na(indiv$SampY),]  ## all the marked animals

    indiv_unmarked <- mort(indiv_unmarked, year = y, type = "age", ageMort = ageMort)
        ## normal mort for unmarked animals
    indiv_marked <- mort(indiv_marked, year = y, type = "age", ageMort = ageMort*2)
        ## double mort for marked animals
    indiv <- rbind(indiv_unmarked, indiv_marked)
        ## stick them back together and proceed as normal
    rm(indiv_marked, indiv_unmarked) ## clean up

    indiv <- birthdays(indiv)
}
```

## Arbitrary variables to affect the simulation

Capture status isn't the only variable that can influence survival or
movement. It is possible to add arbitrary additional variables to an
`indiv` object in order to keep track of anything you might need to
know about an individual, and it is possible to make simulation steps
depend on those variables. You might want a 'hasBred' column, for
instance, keeping track of which animals have produced offspring. Or,
alongside a custom 'myMate' function (see 'custom functions', below),
you might want to keep track of each animal's pair-bonded mate, so
that animals can form pair bonds that persist across years. Here's an
example of a species with two morphs, where the spiny morph has a
slightly different movement-pattern to the smooth morph:

```{r arbitraryExtras}
stocks <- c(0.3, 0.3, 0.4)
indiv <- makeFounders( stocks = stocks)
indiv[,10] <- sample(c("spiny", "smooth"), nrow(indiv), replace = TRUE)
## an extra column to hold stuff. Extras
## _must_ be in column numbers 10 or greater
colnames(indiv)[10] <- "Morph"

## a movement matrix for spiny individuals
admix.spiny <- matrix(NA, nrow = length(stocks), ncol = length(stocks))
for(i in 1:nrow(admix.spiny)) {
    admix.spiny[i,] <- stocks*stocks[i]
}
## admix.spiny shows movement proportional to starting population sizes.
admix.spiny

## a movement matrix for smooth individuals - let's make them more sedentary,
## with 90% chance of staying where they are and 5% chance of moving to each
## other stock
admix.smooth <- matrix(c(0.05), nrow = length(stocks), ncol = length(stocks), byrow = FALSE)
diag(admix.smooth) <- 0.9
admix.smooth <- admix.smooth / sum(admix.smooth) ## standardise so the matrix sums to 1
admix.smooth

for (y in 1:20) {
    indiv <- altMate(indiv, firstBreed = 2, batchSize = 0.9, year = y)
    ## 'mate' happens first, affecting the whole population

    ## newborns need a 'smooth' or 'spiny' status, because altMate won't give one
    indiv$Morph[ is.na( indiv$Morph)] <- sample(c("spiny", "smooth"),
                                                nrow(indiv[ is.na( indiv$Morph),]), replace = TRUE)
    
    indiv <- mort(indiv, year = y, type = "flat", mortRate = 0.2)
    ## then 'mort', also affecting the whole population
    
    indiv_smooth <- indiv[indiv$Morph == "smooth",] ## all the smooth animals
    indiv_spiny <- indiv[indiv$Morph == "spiny",]  ## all the spiny animals

    indiv_smooth <- move(indiv_smooth, moveMat = admix.smooth)
        ## one movement pattern for smooth animals
    indiv_spiny <- move(indiv_spiny, moveMat = admix.spiny)
        ## another movement pattern for spiny animals
    indiv <- rbind(indiv_smooth, indiv_spiny)
        ## stick them back together and proceed as normal
    rm(indiv_smooth, indiv_spiny) ## clean up

    indiv <- birthdays(indiv)
}
```


## Custom functions in fishSim

The objects used in `fishSim` are deliberately kept simple, to
encourage hacking wherever hacking is useful. To that end, most
functions in `fishSim` both input and output a `data.frame` that we
have called `indiv` throughout this vignette.

If your species of interest has a life-history not readily captured by
the existing `fishSim` functions, you may wish to write a custom
function for just that step in the simulation. For instance, you may
have a species in which each female mates with a handful of males in
the year preceding each breeding season, and then uses stored sperm
from these mates to fertilise a clutch of eggs just before giving
birth. Let's call this the 'sperm storage' scenario. The sperm storage
scenario can't be coded using the built-in `mate`, `altMate`, or
`bondedMate` functions - in those functions, each mating event either
has single paternity (i.e., all offspring produced by each female have
the same father) or completely random paternity (i.e., the father of
each offspring is randomly selected from among all mature males in the
same stock), or, using `mate`, within-clutch single paternity but
potentially multiple clutches (each with a different father). In the
situations handled by `mate` and `altMate`, the expected number of
offspring for a female is linearly related to her number of mates,
which doesn't fit with our 'sperm storage' scenario. In the situation
handled by `bondedMate`, one male and one female are socially-bonded
for at least one breeding season, and a female's whole clutch are
fathered by her mate.

To get around the linear relationship between fecundity and number of
mates, we will need to define a custom `mate`-like function, that
inputs an `indiv`-like `data.frame`, generates new offspring the way
the scenario demands, and outputs an `indiv`-like `data.frame` to be
used by the other functions in `fishSim`. Such a function could look
like this:

```{r}
# A custom myMate function
#
# parameter 'indiv' is a standard indiv frame, as from makeFounders()
#
# parameter 'batchSize' sets the number of offspring per breeding female
# such that the number of offspring ~Poisson(lambda = batchSize)
#
# parameter 'breedProb' is the probability that a mature female will
# breed this in this mating event
#
# parameter 'maturityAge' is the knife-edge age at maturity, for both
# males and females
#
# parameter 'year' is the year. It is used to set the $BirthY for new
# recruits
#
# myMate uses a zero-truncated Poisson distribution to determine the
# number of mates for each female that breeds. Parameter 'meanMates'
# sets the T parameter for that distribution

myMate <- function(indiv, batchSize, breedProb, maturityAge, year, meanMates) {

    library(ids)

    ## select mature females and males at current time-step
    matureFemales <- indiv[indiv$Sex == "F" & indiv$AgeLast >= maturityAge & is.na(indiv$DeathY) ,]
    matureMales <- indiv[indiv$Sex == "M" & indiv$AgeLast >= maturityAge & is.na(indiv$DeathY),]
    ## select subset of mature females that will breed (given breedProb)
    iBreed_F <- matureFemales[sample(1:nrow(matureFemales), size = nrow(matureFemales)*breedProb),]
    iBreed_M <- matureMales

    ## if no mature females or males available, don't add new individuals to the population
    if(nrow(iBreed_F) == 0 | nrow(iBreed_M) == 0) { return(indiv) }

    clutchSize <- rpois(nrow(iBreed_F), batchSize) # how many offsprings per breeding females? 
    nMates <- rTruncPoisson(nrow(iBreed_F), T = meanMates) # how many mates per breeding females?
    mates <- lapply(nMates, function(i) sample(iBreed_M$Me, i, replace = FALSE)) # draw mates for each female based on nMates
    # and then assign fathers to each new recruit based on mates identified previously
    # (note each item in the list corresponds to a unique female)
    dads <- lapply(1:length(mates), function(i) sample(mates[[i]], clutchSize[[i]], replace = TRUE))

    # now set-up table of recruits assembling all variables
    recruits <- data.frame(Me=uuid(sum(clutchSize), drop_hyphens = TRUE), # unique ID
                           Sex=sample(c("M","F"), sum(clutchSize), prob = c(0.5,0.5), replace = TRUE),
                           Mum=rep(iBreed_F$Me, clutchSize), # add mum vector that match litter sizes assigned earlier
                           Dad=unlist(dads), # unlist dads (length should match new recruits number)
                           BirthY=year, # rest of meta-data is simple for recruits. birth year
                           DeathY=NA, # not dead yet
                           Stock=1, # single stock for now
                           AgeLast=0, # recruits start at age 0
                           SampY=NA) # not sampled yet

    outs <- rbind(indiv, recruits)
    return(outs)
}
```

With `myMate` defined, we can go on to write a simulation
using the other `fishSim` tools, as in this stripped-down example:

```{r}
indiv <- makeFounders(stocks = c(1))
ageMort <- c(0.305, rep(0.2, 100)) ## age-specific mortality as in the
## first, 'simple' example
for (y in 1:10) {
    indiv <- myMate(indiv = indiv, batchSize = 3, breedProb = 1,
                      maturityAge = 4, year = y, meanMates = 3)
    indiv <- mort(indiv, year = y, type = "age", ageMort = ageMort)
    indiv <- birthdays(indiv)
    if(y %in% c(6:10)) {
	    indiv <- capture(indiv, n = 40, fatal = FALSE, year = y)
    }
}
pairs <- findRelativesPar(indiv = indiv, sampled = TRUE, nCores = 2)
namedRelatives(pairs)
```

### Pair-bonded mating
