diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..c1c9f4d
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1 @@
+.git*
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..35e53cc
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,30 @@
+name: CI
+
+on:
+  push:
+    branches:
+      - '**'
+    paths-ignore:
+      - '*.md'
+      - '*.MD'
+      - '*.ignore'
+      - LICENSE
+     
+jobs:
+  
+  circtools:
+    
+    runs-on: ubuntu-latest
+
+    steps:
+
+      - uses: actions/checkout@v2
+
+      # Run everything via the Bioconductor docker image as it has almost all dependencies preinstalled
+      - name: devtools-check-docker
+        run: |
+          bioc_install='BiocManager::install(c("limma", "ggplot2"))'
+          dev_check='devtools::check("/circtools/")'
+          dev_doc='devtools::document("/circtools/")'
+          testthat='testthat::test_file("/circtools/tests/testthat/all_tests.R")'
+          docker run -v "$(pwd)":"/circtools/" bioconductor/bioconductor_docker:RELEASE_3_18 Rscript --vanilla -e "${bioc_install}; ${dev_check}; ${dev_doc}; ${testthat}"
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..83e6752
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,12 @@
+*.RProj
+*.Rproj
+*.rproj
+..R*
+gitHeadInfo.gin
+R/.Rapp.history
+*.DS_*
+.Rproj.user
+.Rhistory
+qmethod.Rproj
+~*
+*.tar.gz
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..279a4dd
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,15 @@
+Package: circtools
+Title: Functions related to analysis of circular data in R
+Version: 0.0.9
+Description: Functions related to analysis of circular data in R.
+Author: Alexander Bender [cut,cre]
+Maintainer: Alexander Bender <atpoint90@gmail.com>
+License: LGPL (>=2)
+Encoding: UTF-8
+Imports: 
+  limma, stats, methods
+Suggests:
+  ggplot2
+URL: https://github.com/atpoint/circtools
+BugReports: https://github.com/atpoint/circtools/issues
+RoxygenNote: 7.2.3
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..107ba35
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,18 @@
+# Generated by roxygen2: do not edit by hand
+
+export(circular_density)
+export(circular_distance)
+export(circular_distance_reference)
+export(circular_intersect)
+export(make_data)
+export(make_design)
+export(rad2period)
+export(run_cosinor)
+importFrom(limma,eBayes)
+importFrom(limma,lmFit)
+importFrom(limma,topTable)
+importFrom(methods,is)
+importFrom(methods,new)
+importFrom(stats,density)
+importFrom(stats,model.matrix)
+importFrom(stats,rnorm)
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..1d5f898
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,328 @@
+#' Detection of rhythmicity using limma with a cosinor model.
+#'
+#' The function detects periodic rhythmicity based on the cosinor model as described in the limma design matrix guide:
+#' \url{https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/designmatrices.html#cyclical-models}
+#' For this it uses the limma framework to first fit per-observation cosinor models with \code{limma::lmFit()} and then
+#' obtains moderated pvalues and FDRs with \code{limma::eBayes()}. It also returns estimates for mesor, amplitude and acrophase
+#' based on the model coefficients. The function assumes that input data are "ready-to-use",  meaning that normalization, prefiltering (if applicable) appropriate 
+#' for this type of data has been done. 
+#'
+#' Since the input data `x` (see argument description) can be both a numeric matrix or an `EList`, the user has a great flexibility over the analysis:
+#' For an ordinary limma analysis (default settings) `x` could be a numeric matrix with data normalized and on log2 scale.
+#' For a limma-trend analysis (\code{eBayes_args=list(trend=TRUE)}) `x` could be a numeric matrix of logCPMs from an RNA-seq experiment. 
+#' For a limma-voom analysis `x` could be the output of \code{limma::voom()}.
+#' Since limma automatically uses **weights** if provided in the EList as x$weights, one can put them there in case weighted regression is desired.
+#' Weights could be estimated using \code{limma::arrayWeights} based on the cosinor model. A `model.matrix()` for the cosinor model can be obtained with \code{make_design()}.
+#' Alternatively, `x` could be the output of \code{limma::voomWithQualityWeights()}. 
+#' 
+#'
+#' @param x A numeric matrix or data.frame with observations in rows and samples in columns. Alternatively, an \code{EList} object.
+#' @param time A numeric vector with the time information for each column of x.
+#' @param period The numeric period of the rhythmicity, default 24.
+#' @param lmFit_args A list with arguments for \code{limma::lmFit()}.
+#' @param eBayes_args A list with arguments for \code{limma::eBayes()}.
+#'
+#' @examples
+#' period <- 24
+#' n_reps <- 1
+#' time <- rep(seq(0, 24, by = 2), each = n_reps)
+#' acrophases <- c(1, 10, 21)
+#' y <- lapply(acrophases, function(ac) make_data(time = time, acrophase = ac, error_sd = .1)$value)
+#'
+#' # x is a numeric matrix
+#' x <- matrix(unlist(y), byrow = TRUE, ncol = length(time))
+#' rownames(x) <- paste0("acrophase_", acrophases)
+#' run_cosinor(x = x, time = time)
+#'
+#' # x is an EList
+#' xe <- new("EList", list(E = x))
+#' run_cosinor(x = xe, time = time)
+#'
+#' # x is an EList with weights
+#' set.seed(1)
+#' xew <- new("EList", list(E = x, weights = rnorm(ncol(x), 1)))
+#' run_cosinor(x = xew, time = time)
+#'
+#' @importFrom limma lmFit eBayes topTable
+#' @importFrom methods is new
+#'
+#' @export
+#'
+run_cosinor <- function(x, time, period = 24, lmFit_args = list(), eBayes_args = list()) {
+  is_class <- is(x, "EList") | is(x, "data.frame") | is(x, "matrix")
+  if (!is_class) stop("x must either be a numeric matrix or data.frame, or a limma EList")
+
+  # Ensure consistent input format and row/colnames are set
+  if (is(x, "EList")) {
+    x$E <- as.matrix(x$E)
+
+    rownames(x$E) <- if (is.null(rownames(x$E))) paste0("row_", 1:nrow(x$E)) else rownames(x$E)
+    colnames(x$E) <- if (is.null(colnames(x$E))) paste0("sample_", 1:ncol(x$E)) else colnames(x$E)
+
+    y <- x$E
+  } else {
+    x <- as.matrix(x)
+
+    rownames(x) <- if (is.null(rownames(x))) paste0("row_", 1:nrow(x)) else rownames(x)
+    colnames(x) <- if (is.null(colnames(x))) paste0("sample_", 1:ncol(x)) else colnames(x)
+
+    y <- x
+  }
+
+  # Check input is numeric and finite
+  is_numeric <- is.numeric(as.vector(y))
+  has_NAs <- sum(!is.finite(as.vector(y))) > 0
+  if (!is_numeric | has_NAs) stop("x contains non-numeric entries and/or NAs/NaNs/infinite values")
+  rm(y)
+
+  # Check that time is numeric and of same length as samples in x
+  if (!is.numeric(time)) stop("time must be numeric")
+  len_time <- length(time)
+  len_x <- ncol(x)
+  if (len_time != len_x) stop("time must have the same length as ncol(x) -- one time entry per column")
+
+  if (!is.numeric(period) | period < 0) stop("period must be numeric and positive")
+  if (!is.list(lmFit_args)) stop("lmFit_args must be a list")
+  if (!is.list(eBayes_args)) stop("eBayes_args must be a list")
+
+  # Fit cosinor model and extract stats
+  design <- make_design(time = time, period = period)
+  fit <- do.call(lmFit, c(list(x, design = design), lmFit_args))
+  fit <- do.call(eBayes, c(list(fit), eBayes_args))
+  coef_pos <- grep("^sinphase$|^cosphase$", colnames(fit$coefficients))
+  tt <- topTable(fit = fit, coef = coef_pos, number = Inf, sort.by = "none")
+  tt <- data.frame(id = rownames(x), tt)
+
+  tt$amplitude <- sqrt(rowSums(tt[, c("sinphase", "cosphase")]^2))
+  tt$acrophase <- (atan2(tt$sinphase, tt$cosphase) / 2 / pi * period + period) %% period
+  tt <- tt[, c("id", "P.Value", "adj.P.Val", "amplitude", "acrophase", "AveExpr", "sinphase", "cosphase", "F")]
+
+  return(tt)
+}
+
+#' Simulate cosinor data
+#'
+#' Accept numeric timepoints and cosinor parameters and returns a data.frame with time and values.
+#'
+#' @param time A numeric vector with time information
+#' @param amplitude The amplitude of the cosinor curve
+#' @param acrophase The acrophase of the cosinor curve
+#' @param period The period of the time series
+#' @param mesor The mesor of the cosinor curve
+#' @param error_sd An error value to simulate noise
+#' @param use_seed A random seed to make error_sd reproducible
+#'
+#' @examples
+#' tt <- seq(0, 23.9, by = .1)
+#' data <- make_data(time = tt, acrophase = 5)
+#' plot(data$time, data$value, pch = 20)
+#' abline(v = 5)
+#'
+#' @importFrom stats rnorm
+#'
+#' @export
+#'
+make_data <- function(time, amplitude = 1, acrophase = 1, period = 24, mesor = 0, error_sd = .1, use_seed = 1) {
+  if (!is.numeric(time) | sum(time < 0) > 0 | length(time) == 1) stop("time must be a numeric vector with non-negative positive values")
+
+  set.seed(use_seed)
+  error <- rnorm(length(time), mean = 0, sd = error_sd)
+  v <- mesor + amplitude * cos(2 * pi / period * (time - acrophase)) + error
+
+  data <- data.frame(time = time, value = v)
+
+  return(data)
+}
+
+#' Make a design matrix using a cosinor model.
+#'
+#' @param time A numeric vector with timepoints
+#' @param period The period of the time series
+#'
+#' @importFrom stats model.matrix
+#' @export
+#'
+make_design <- function(time, period = 24) {
+  if (!is.numeric(time)) stop("time must be numeric")
+  if (!is.numeric(period) | period < 0) stop("period must be numeric and positive")
+
+  d <- data.frame(
+    sinphase = sin(2 * pi * time / period),
+    cosphase = cos(2 * pi * time / period)
+  )
+
+  design <- model.matrix(~ sinphase + cosphase, d)
+
+  return(design)
+}
+
+#' Calculate densities of periodic values such as acrophases between user-defined start and end points, respecting the circular nature of the data.
+#' 
+#' Density calculation is prone to skewed density estimates at the start and end of the data range if the method is not aware of the circular nature of the data.
+#' For that this function appends the data to itself at the front and rear end, and then uses \code{stats::density} with a default bandwidth of 0.5 to calculate densities
+#' that can be used for histogram/ridge-like plots. The user should provide a start point `from` representing the smallest possible value of the periodoc values, and an end point
+#' `to` which represents the period of the circular data / the time series.
+#'
+#' @param x Numeric vector with input value.
+#' @param from,to The left and right-most points of the grid at which the density is to be estimated.
+#'   The defaults are 0 and 24, assuming a 24h clock cycle. The "to" argument is typically the period of the time series.
+#' @param density_args arguments for \code{stats::density}
+#'
+#' @examples
+#'
+#' set.seed(1)
+#' data <- abs(c(rnorm(1000, 12, 5), rnorm(500, 23))) %% 24
+#'
+#' # linear density estimation incorrectly has low density at around x=0, as it does not know that at
+#' # around x=24 density is high
+#' library(ggplot2)
+#' ggplot(data = data.frame(y = data), aes(x = y)) +
+#'   geom_density()
+#'
+#' # circular density estimation is aware that x=0 is almost identical to x=24
+#' data_circ <- circular_density(data)
+#' ggplot(data = data_circ, aes(x = time, y = density)) +
+#'   geom_line()
+#'
+#' @importFrom stats density
+#'
+#' @export
+#'
+circular_density <- function(x, from = 0, to = 24, density_args = list()) {
+  if (!is.numeric(x)) stop("x must be numerical")
+  ft <- c(from, to)
+  if (!is.numeric(ft) | sum(ft < 0) > 0 | sum(is.na(ft)) > 0) stop("from/to must be non-negative & numeric")
+
+  # Set a default bandwidth if none is set via density_args
+  if (!is.list(density_args)) stop("density_args must be a list")
+  if (!"bw" %in% names(density_args)) {
+    density_args$bw <- .5
+  }
+
+  # Append data to itself to respect circular nature
+  y <- c(x - to, x, x + to)
+
+  h <- do.call(density, c(list(y, from = from, to = to), density_args))
+  d <- data.frame(time = h$x, density = h$y)
+  d <- d[d$time >= from & d$time <= to, ]
+
+  return(d)
+}
+
+#' Convert radians hours to period hours.
+#'
+#' @param x a radians time value, from 0 to 2*pi.
+#' @param period The numeric period.
+#'
+#' @export
+#'
+rad2period <- function(x, period = 24) {
+  if (x < 0 | x > (2 * pi)) stop("x must be between 0 and 2*pi")
+  if (!is.numeric(period) | period < 0) stop("period must be numeric and not be negative")
+
+  period_hours <- (x / (2 * pi)) * period
+
+  return(period_hours)
+}
+
+#' Circular distance between two values, given a period.
+#'
+#' Calculates circular distance and direction with respect to a period.
+#' Positive values mean that x is ahead of y in circular space.
+#' For example, x=1 and y=23 gives +2 because the circular distance is 2 and
+#' on a clock with 24h x is two hours later than y. Likewise,
+#' x=23 and y=1 gives negative 2 because distance again is 2 but on a 24h clock x
+#' is earlier than y.
+#'
+#' @param x,y The two values to calculate circular distance for.
+#' @param period THe period of the time series
+#' 
+#' @examples
+#' # result is +2 because on a circular clock x is 2h ahead of y
+#' circular_distance(x=1, y=23, period=24)
+#' 
+#' # result is -5 because on a circular clock x is 5h behind y
+#' circular_distance(x=15, y=20, period=24)
+#'
+#' @export
+#'
+circular_distance <- function(x, y, period = 24) {
+  xyp <- c(x, y, period)
+
+  if (!is.numeric(xyp)) stop("x, y, period must be numeric and >= 0")
+  if (sum(!is.finite(xyp)) > 0) stop("x, y, period must be finite")
+  if (length(xyp) != 3) stop("x, y, period must all have length of 1")
+  if (sum(xyp < 0) > 0) stop("x, y, period must not be negative")
+
+  diff <- abs(x - y)
+  d <- min(diff, period - diff)
+
+  steps1 <- ((x - y) %% period + period) %% period
+  steps2 <- ((y - x) %% period + period) %% period
+
+  z <- if (steps1 <= steps2) d else d * (-1)
+
+  return(z)
+}
+
+#' Calculate absolute circular distance between one query value and a vector with one or many values given a fixed period.
+#'
+#' @param query A numeric value
+#' @param reference A vector with numeric values
+#' @param period A period for the circle
+#'
+#' @export
+#'
+circular_distance_reference <- function(query, reference, period = 24) {
+  qr <- c(query, reference)
+  qrp <- c(query, reference, period)
+
+  if (!is.numeric(qrp) | sum(!is.finite(qrp)) > 0 | sum(qrp < 0) > 0) stop("query, reference and period must all be numeric and finite and >= 0")
+  if (sum(qr < 0) > 0) stop()
+
+  diff <- abs(reference - query)
+  d <- pmin(diff, period - diff)
+
+  return(d)
+}
+
+#' Circular-aware intersection of a query value with a time window.
+#'
+#' Intersect a single timepoint with a timepoint window in a circular-aware fashion given a period,
+#' returning TRUE or FALSE if there is an intersection or not. The lower limit of the window is inclusive,
+#' the upper one is not. Meaning, a query of 2 or a window of 2-4 would return TRUE but for a window 1-2 would
+#' return FALSE.
+#'
+#' @param query The query value.
+#' @param window A vector of two values representing the lower and upper bounds of the time window.
+#'
+#' @examples
+#' query <- c(2)
+#' window1 <- c(0, 3) # TRUE
+#' window2 <- c(23, 3) # TRUE
+#' window3 <- c(3, 5) # FALSE
+#' window4 <- c(23, 2) # FALSE
+#'
+#' circular_intersect(query, window1)
+#' circular_intersect(query, window2)
+#' circular_intersect(query, window3)
+#' circular_intersect(query, window4)
+#'
+#' @export
+#'
+circular_intersect <- function(query, window) {
+  cc <- c(query, window)
+
+  if (!is.numeric(cc) > 0 | sum(!is.finite(cc)) > 0 | sum(cc < 0) > 0) stop("Both query and window must be numeric, finite and non-negative")
+  if (length(query) > 1) stop("query must be a single value")
+  if (length(window) != 2) stop("window must be of length 2")
+  if (length(unique(window)) == 1) stop("The two window values are the same")
+
+  if (window[1] > window[2]) {
+    x <- query >= window[1] | query < window[2]
+  } else {
+    x <- query >= window[1] & query < window[2]
+  }
+
+  return(x)
+}
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..844ac41
--- /dev/null
+++ b/README.md
@@ -0,0 +1,31 @@
+# circtools
+
+A little package with functions related to the analysis of circadian / circular data in R. 
+Note that it is mainly intended for my own work, but others are free to use it if you find it useful.
+
+It includes functions to detect rhythmicity in circadian data using cosinor regression as well as
+adaptations of basic data analysis functions such as intersection, distance and density calculation but modified to
+respect the circular / periodic nature of the input data. For details, please see the function help at `?function_name`.
+
+Functions:  
+  
+- `run_cosinor()` detects circadian rhythmicity using a cosinor model via the [limma](https://bioconductor.org/packages/release/bioc/html/limma.html) framework.
+It scales efficiently to almost any realistic number of samples and observations.  
+- `circular_density()` calculates densities of a vector of circular/periodic values, such as acrophases, respecting the circular nature of the data, thereby avoiding skewed density estimates at the beginning and end of the data range  
+- `circular_distance()` calculates directional circular distance between two timepoints  
+- `circular_distance_reference()` calculates absolute circular distance between a query and one or many reference timepoints  
+- `circular_intersect()` intersects a timepoint with a time window respecting periodicity of data  
+- `make_data()` simulates circadian data based on a cosinor model with adjustable mesor, amplitude and acrophase  
+- `make_design()` creates a design matrix based on the cosinor model for a provided set of timepoints  
+- `rad2period()` converts acrohase values from radians hours to period hours  
+
+## Installation
+
+```r
+# Install Bioconductor
+if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
+BiocManager::install()
+
+# Install limma and circtools
+BiocManager::install(c("limma", "atpoint/circtools"))
+```
diff --git a/man/circular_density.Rd b/man/circular_density.Rd
new file mode 100644
index 0000000..ca7f07f
--- /dev/null
+++ b/man/circular_density.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{circular_density}
+\alias{circular_density}
+\title{Calculate densities of periodic values such as acrophases between user-defined start and end points, respecting the circular nature of the data.}
+\usage{
+circular_density(x, from = 0, to = 24, density_args = list())
+}
+\arguments{
+\item{x}{Numeric vector with input value.}
+
+\item{from, to}{The left and right-most points of the grid at which the density is to be estimated.
+The defaults are 0 and 24, assuming a 24h clock cycle. The "to" argument is typically the period of the time series.}
+
+\item{density_args}{arguments for \code{stats::density}}
+}
+\description{
+Density calculation is prone to skewed density estimates at the start and end of the data range if the method is not aware of the circular nature of the data.
+For that this function appends the data to itself at the front and rear end, and then uses \code{stats::density} with a default bandwidth of 0.5 to calculate densities
+that can be used for histogram/ridge-like plots. The user should provide a start point `from` representing the smallest possible value of the periodoc values, and an end point
+`to` which represents the period of the circular data / the time series.
+}
+\examples{
+
+set.seed(1)
+data <- abs(c(rnorm(1000, 12, 5), rnorm(500, 23))) \%\% 24
+
+# linear density estimation incorrectly has low density at around x=0, as it does not know that at
+# around x=24 density is high
+library(ggplot2)
+ggplot(data = data.frame(y = data), aes(x = y)) +
+  geom_density()
+
+# circular density estimation is aware that x=0 is almost identical to x=24
+data_circ <- circular_density(data)
+ggplot(data = data_circ, aes(x = time, y = density)) +
+  geom_line()
+
+}
diff --git a/man/circular_distance.Rd b/man/circular_distance.Rd
new file mode 100644
index 0000000..47268ce
--- /dev/null
+++ b/man/circular_distance.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{circular_distance}
+\alias{circular_distance}
+\title{Circular distance between two values, given a period.}
+\usage{
+circular_distance(x, y, period = 24)
+}
+\arguments{
+\item{x, y}{The two values to calculate circular distance for.}
+
+\item{period}{THe period of the time series}
+}
+\description{
+Calculates circular distance and direction with respect to a period.
+Positive values mean that x is ahead of y in circular space.
+For example, x=1 and y=23 gives +2 because the circular distance is 2 and
+on a clock with 24h x is two hours later than y. Likewise,
+x=23 and y=1 gives negative 2 because distance again is 2 but on a 24h clock x
+is earlier than y.
+}
+\examples{
+# result is +2 because on a circular clock x is 2h ahead of y
+circular_distance(x=1, y=23, period=24)
+
+# result is -5 because on a circular clock x is 5h behind y
+circular_distance(x=15, y=20, period=24)
+
+}
diff --git a/man/circular_distance_reference.Rd b/man/circular_distance_reference.Rd
new file mode 100644
index 0000000..d555bc5
--- /dev/null
+++ b/man/circular_distance_reference.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{circular_distance_reference}
+\alias{circular_distance_reference}
+\title{Calculate absolute circular distance between one query value and a vector with one or many values given a fixed period.}
+\usage{
+circular_distance_reference(query, reference, period = 24)
+}
+\arguments{
+\item{query}{A numeric value}
+
+\item{reference}{A vector with numeric values}
+
+\item{period}{A period for the circle}
+}
+\description{
+Calculate absolute circular distance between one query value and a vector with one or many values given a fixed period.
+}
diff --git a/man/circular_intersect.Rd b/man/circular_intersect.Rd
new file mode 100644
index 0000000..39a7ec3
--- /dev/null
+++ b/man/circular_intersect.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{circular_intersect}
+\alias{circular_intersect}
+\title{Circular-aware intersection of a query value with a time window.}
+\usage{
+circular_intersect(query, window)
+}
+\arguments{
+\item{query}{The query value.}
+
+\item{window}{A vector of two values representing the lower and upper bounds of the time window.}
+}
+\description{
+Intersect a single timepoint with a timepoint window in a circular-aware fashion given a period,
+returning TRUE or FALSE if there is an intersection or not. The lower limit of the window is inclusive,
+the upper one is not. Meaning, a query of 2 or a window of 2-4 would return TRUE but for a window 1-2 would
+return FALSE.
+}
+\examples{
+query <- c(2)
+window1 <- c(0, 3) # TRUE
+window2 <- c(23, 3) # TRUE
+window3 <- c(3, 5) # FALSE
+window4 <- c(23, 2) # FALSE
+
+circular_intersect(query, window1)
+circular_intersect(query, window2)
+circular_intersect(query, window3)
+circular_intersect(query, window4)
+
+}
diff --git a/man/make_data.Rd b/man/make_data.Rd
new file mode 100644
index 0000000..1021b31
--- /dev/null
+++ b/man/make_data.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{make_data}
+\alias{make_data}
+\title{Simulate cosinor data}
+\usage{
+make_data(
+  time,
+  amplitude = 1,
+  acrophase = 1,
+  period = 24,
+  mesor = 0,
+  error_sd = 0.1,
+  use_seed = 1
+)
+}
+\arguments{
+\item{time}{A numeric vector with time information}
+
+\item{amplitude}{The amplitude of the cosinor curve}
+
+\item{acrophase}{The acrophase of the cosinor curve}
+
+\item{period}{The period of the time series}
+
+\item{mesor}{The mesor of the cosinor curve}
+
+\item{error_sd}{An error value to simulate noise}
+
+\item{use_seed}{A random seed to make error_sd reproducible}
+}
+\description{
+Accept numeric timepoints and cosinor parameters and returns a data.frame with time and values.
+}
+\examples{
+tt <- seq(0, 23.9, by = .1)
+data <- make_data(time = tt, acrophase = 5)
+plot(data$time, data$value, pch = 20)
+abline(v = 5)
+
+}
diff --git a/man/make_design.Rd b/man/make_design.Rd
new file mode 100644
index 0000000..bebe64c
--- /dev/null
+++ b/man/make_design.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{make_design}
+\alias{make_design}
+\title{Make a design matrix using a cosinor model.}
+\usage{
+make_design(time, period = 24)
+}
+\arguments{
+\item{time}{A numeric vector with timepoints}
+
+\item{period}{The period of the time series}
+}
+\description{
+Make a design matrix using a cosinor model.
+}
diff --git a/man/rad2period.Rd b/man/rad2period.Rd
new file mode 100644
index 0000000..d0b13e5
--- /dev/null
+++ b/man/rad2period.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{rad2period}
+\alias{rad2period}
+\title{Convert radians hours to period hours.}
+\usage{
+rad2period(x, period = 24)
+}
+\arguments{
+\item{x}{a radians time value, from 0 to 2*pi.}
+
+\item{period}{The numeric period.}
+}
+\description{
+Convert radians hours to period hours.
+}
diff --git a/man/run_cosinor.Rd b/man/run_cosinor.Rd
new file mode 100644
index 0000000..280b43f
--- /dev/null
+++ b/man/run_cosinor.Rd
@@ -0,0 +1,58 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{run_cosinor}
+\alias{run_cosinor}
+\title{Detection of rhythmicity using limma with a cosinor model.}
+\usage{
+run_cosinor(x, time, period = 24, lmFit_args = list(), eBayes_args = list())
+}
+\arguments{
+\item{x}{A numeric matrix or data.frame with observations in rows and samples in columns. Alternatively, an \code{EList} object.}
+
+\item{time}{A numeric vector with the time information for each column of x.}
+
+\item{period}{The numeric period of the rhythmicity, default 24.}
+
+\item{lmFit_args}{A list with arguments for \code{limma::lmFit()}.}
+
+\item{eBayes_args}{A list with arguments for \code{limma::eBayes()}.}
+}
+\description{
+The function detects periodic rhythmicity based on the cosinor model as described in the limma design matrix guide:
+\url{https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/designmatrices.html#cyclical-models}
+For this it uses the limma framework to first fit per-observation cosinor models with \code{limma::lmFit()} and then
+obtains moderated pvalues and FDRs with \code{limma::eBayes()}. It also returns estimates for mesor, amplitude and acrophase
+based on the model coefficients. The function assumes that input data are "ready-to-use",  meaning that normalization, prefiltering (if applicable) appropriate 
+for this type of data has been done.
+}
+\details{
+Since the input data `x` (see argument description) can be both a numeric matrix or an `EList`, the user has a great flexibility over the analysis:
+For an ordinary limma analysis (default settings) `x` could be a numeric matrix with data normalized and on log2 scale.
+For a limma-trend analysis (\code{eBayes_args=list(trend=TRUE)}) `x` could be a numeric matrix of logCPMs from an RNA-seq experiment. 
+For a limma-voom analysis `x` could be the output of \code{limma::voom()}.
+Since limma automatically uses **weights** if provided in the EList as x$weights, one can put them there in case weighted regression is desired.
+Weights could be estimated using \code{limma::arrayWeights} based on the cosinor model. A `model.matrix()` for the cosinor model can be obtained with \code{make_design()}.
+Alternatively, `x` could be the output of \code{limma::voomWithQualityWeights()}.
+}
+\examples{
+period <- 24
+n_reps <- 1
+time <- rep(seq(0, 24, by = 2), each = n_reps)
+acrophases <- c(1, 10, 21)
+y <- lapply(acrophases, function(ac) make_data(time = time, acrophase = ac, error_sd = .1)$value)
+
+# x is a numeric matrix
+x <- matrix(unlist(y), byrow = TRUE, ncol = length(time))
+rownames(x) <- paste0("acrophase_", acrophases)
+run_cosinor(x = x, time = time)
+
+# x is an EList
+xe <- new("EList", list(E = x))
+run_cosinor(x = xe, time = time)
+
+# x is an EList with weights
+set.seed(1)
+xew <- new("EList", list(E = x, weights = rnorm(ncol(x), 1)))
+run_cosinor(x = xew, time = time)
+
+}
diff --git a/tests/testthat/all_tests.R b/tests/testthat/all_tests.R
new file mode 100644
index 0000000..e046625
--- /dev/null
+++ b/tests/testthat/all_tests.R
@@ -0,0 +1,71 @@
+test_that("run_cosinor works", {
+  period <- 24
+  n_reps <- 1
+  time <- rep(seq(0, 24, by = 2), each = n_reps)
+  acrophases <- c(1, 10, 21)
+  y <- lapply(acrophases, function(ac) make_data(time = time, acrophase = ac, error_sd = .1)$value)
+  x <- matrix(unlist(y), byrow = TRUE, ncol = length(time))
+  expect_true(is(run_cosinor(x = x, time = time), "data.frame"))
+
+  xE <- new("EList", list(E = x))
+  expect_true(is(run_cosinor(x = xE, time = time), "data.frame"))
+
+  xE$weights <- rnorm(nrow(xE), 1, .01)
+  expect_true(is(run_cosinor(x = xE, time = time), "data.frame"))
+
+  xNA <- x
+  xNA[1, 1] <- NA
+  expect_error(run_cosinor(list(1), time = time))
+  expect_error(run_cosinor(x = xNA, time = time))
+  expect_error(run_cosinor(x = xNA, time = time[1:2]))
+  expect_error(run_cosinor(x = xNA, time = as.character(time)))
+  expect_error(run_cosinor(x = x, time = time, period = "24"))
+  expect_error(run_cosinor(x = xNA, time = time, eBayes_args = "1"))
+  expect_error(run_cosinor(x = xNA, time = time, lmFit_args = "1"))
+})
+
+test_that("make_design works", {
+  expect_true(is(make_design(time = c(1, 2, 3), period = 24), "matrix"))
+  expect_error(make_design(time = c(1, 2, "3")))
+  expect_error(make_design(time = c(1, 2, 3), period = "24"))
+  expect_error(make_design(time = c(1, 2, 3), period = -1))
+})
+
+test_that("circular_density works", {
+  data <- abs(c(rnorm(1000, 12, 5), rnorm(500, 23))) %% 24
+  expect_true(is(circular_density(data), "data.frame"))
+  expect_error(circular_density(data, to = -1))
+  expect_error(circular_density(LETTERS[1:10], from = -1))
+  expect_error(circular_density(data, to = -1))
+  expect_error(circular_density(data, to = NA))
+  expect_error(circular_density(data, density_args = 1))
+})
+
+test_that("rad2period works", {
+  expect_type(rad2period(pi), "double")
+  expect_error(rad2period(2 * pi + 1))
+  expect_error(rad2period(1, period = -1))
+})
+
+test_that("circular_distance works", {
+  expect_type(circular_distance(1, 2), "double")
+  expect_error(circular_distance(1, NA))
+  expect_error(circular_distance(1, 1, NA))
+  expect_error(circular_distance(NA, 1, 2))
+  expect_error(circular_distance(1, 2, -3))
+})
+
+test_that("circular_distance_reference works", {
+  expect_type(circular_distance_reference(1, c(1, 2, 3), period = 24), "double")
+  expect_error(circular_distance_reference(1, 1, NA))
+  expect_error(circular_distance_reference(-1, 1, 24))
+  expect_error(circular_distance_reference(1, -1, 24))
+})
+
+test_that("circular_intersect works", {
+  expect_type(circular_intersect(1, c(1, 2)), "logical")
+  expect_error(circular_intersect(-1, c(1, 2)))
+  expect_error(circular_intersect(1, c(1, -2)))
+  expect_error(circular_intersect(1, c(-1, 2)))
+  expect_error(circular_intersect(1, c(1, 1)))
+})