text:
#=

AssetsCollection
================
    A container to hold all information regarding the collection of assests to
    be optimized

    Methods:
    -------

Author: Fady Shoukry, Azamat Berdyshev, Shen Wang
Date: 01/23/2016
=#

type AssetsCollection{T1<:Real, T2<:AbstractString}
    # NOTE: All arrays in the collection are expected to have matching indices

    # A list of indicator strings for the assets in the collection
    names::Vector{T2}
    # A list of the expected returns of the assets in the collection
    returns::Vector{T1}
    # A covariance matrix for the assets in the collection
    covariance::Matrix{T1}

    function AssetsCollection(names, returns, covariance)
        # Do some validation work on the covariance matrix
        if !issymmetric(covariance)
        # Make sure it's symmetric
            error("Covariance Matrix must be symmetric")
        elseif !isposdef(covariance)
        # Make sure it's positive semi-definite
            error("Covariance Matrix must be positive semi-definite")
        end
        n_cov = size(covariance, 1)
        n_ret = size(returns, 1)
        n_names = size(names, 1)
        # Make sure the size of the inputs match
        if n_cov != n_ret
            error("Covariance matrix size ($n_cov x $n_cov) and the length of
            returns ($n_ret) must be equal")
        elseif n_ret != n_names
            error("length of names ($n_names) and the length of
            returns ($n_ret) must be equal")
        end

        new(names, returns, covariance)
    end
end

AssetsCollection{T1<:Real, T2<:AbstractString}(names::Vector{T2},
    returns::Vector{T1}, covariance::Matrix{T1}) =
    AssetsCollection{T1, T2}(names, returns, covariance)

function length{T1<:Real,
    T2<:AbstractString}(assets::AssetsCollection{T1, T2})
    """
    Return the number of assets in the collection
    """
    size(assets.returns, 1)
end

function getReturnForAsset{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, name::AbstractString)
    """
    Access the expected return of a given asset by its name
    """
    idx = findfirst(assets.names, name)
    if idx == 0
        error("The asset collection has no asset named $name")
    end
    assets.returns[idx]
end

function setReturnForAsset{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, name::AbstractString, value)
    """
    Set the expected return of a given asset by its name
    """
    # First convert the value to the correct type
    value = convert(T1, value)
    idx = findfirst(assets.names, name)
    if idx == 0
        error("The asset collection has no asset named $name")
    end
    assets.returns[idx] = value
end

function getVarForAsset{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, name::AbstractString)
    """
    Access the variance of a given asset by its name
    """
    idx = findfirst(assets.names, name)
    if idx == 0
        error("The asset collection has no asset named $name")
    end
    assets.covariance[idx, idx]
end

function setVarForAsset{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, name::AbstractString, value)
    """
    Set the variance of a given asset by its name
    """
    # First convert the value to the correct type
    value = convert(T1, value)
    idx = findfirst(assets.names, name)
    if idx == 0
        error("The asset collection has no asset named $name")
    end
    saved_value = assets.covariance[idx, idx]
    assets.covariance[idx, idx] = value
    # Check to see if the positive semi-definitiveness has been violated
    if !isposdef(assets.covariance)
        assets.covariance[idx, idx] = saved_value
        error("The updated covariance matrix is no longer positive
        semi-definite")
    end
end

function getCoVarForAssetPair{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, asset1::AbstractString,
    asset2::AbstractString)
    """
    Access the covariance of 2 given assets by their names
    """
    idx_1 = findfirst(assets.names, asset1)
    idx_2 = findfirst(assets.names, asset2)
    if idx_1 == 0
        error("The asset collection has no asset named $asset1")
    elseif idx_2 == 0
        error("The asset collection has no asset named $asset2")
    end
    assets.covariance[idx_1, idx_2]
end

function setCoVarForAssetPair{T1<:Real, T2<:AbstractString}(
    assets::AssetsCollection{T1, T2}, asset1::AbstractString,
    asset2::AbstractString, value)
    """
    Set the covariance of 2 given assets by their names
    """
    # First convert the value to the correct type
    value = convert(T1, value)
    idx_1 = findfirst(assets.names, asset1)
    idx_2 = findfirst(assets.names, asset2)
    if idx_1 == 0
        error("The asset collection has no asset named $asset1")
    elseif idx_2 == 0
        error("The asset collection has no asset named $asset2")
    end
    saved_value = assets.covariance[idx_1, idx_2]
    assets.covariance[idx_1, idx_2] = value
    assets.covariance[idx_2, idx_1] = value
    # Check to see if the positive semi-definitiveness has been violated
    if !isposdef(assets.covariance)
        assets.covariance[idx_1, idx_2] = saved_value
        assets.covariance[idx_2, idx_1] = saved_value
        error("The updated covariance matrix is no longer positive
        semi-definite")
    end
end

function getCovariance{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2})
    return assets.covariance::Matrix{T1}
end

function setCovariance{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2}, covariance::Matrix{T1})
    # Perform validation of the covariance matrix
    n = length(assets)
    if n != size(covariance, 1)
        # Make sure it's the correct size
        error("Covariance Matrix must have size $n x $n")
    elseif !issymmetric(covariance)
    # Make sure it's symmetric
        error("Covariance Matrix must be symmetric")
    elseif !isposdef(covariance)
    # Make sure it's positive semi-definite
        error("Covariance Matrix must be positive semi-definite")
    end

    assets.covariance = covariance
end

function getReturns{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2})
    return assets.returns::Vector{T1}
end

function setReturns{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2}, returns::Vector{T1})
    assets.returns = returns
end

function getNames{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2})
    return assets.names::Vector{T2}
end

function setNames{T1<:Real, T2<:AbstractString}(assets::AssetsCollection{T1, T2}, names::Vector{T2})
    assets.names = names
end

function Base.show(io::IO, a::AssetsCollection)

    a_length = length(a.names)
    a_df = DataFrames.DataFrame(Names = AbstractString[], Expected_Returns = Real[])

    for i = 1:a_length
        push!(a_df, [a.names[i], a.returns[i]])
    end

    print(io, a_df)
end

function getAssetAndReturnsFromCSV(filepath::AbstractString)
    asset_table = DataFrames.readtable(filepath, header=false)
    return vec(convert(Vector{ASCIIString}, asset_table[1])), vec(convert(Vector{Float64}, asset_table[2]))
end


text:
# Autogenerated wrapper script for ROCmCompilerSupport_jll for x86_64-linux-gnu-cxx11
export libamd_comgr

using hsa_rocr_jll
using ROCmDeviceLibs_jll
JLLWrappers.@generate_wrapper_header("ROCmCompilerSupport")
JLLWrappers.@declare_library_product(libamd_comgr, "libamd_comgr.so.2")
function __init__()
    JLLWrappers.@generate_init_header(hsa_rocr_jll, ROCmDeviceLibs_jll)
    JLLWrappers.@init_library_product(
        libamd_comgr,
        "lib64/libamd_comgr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using Slack
using Test

endpoint = ARGS[1]
token = ARGS[2]

data = "Hello World"
attachment = Dict("attachments" => [Dict("fallback" => "Required plain-text summary of the attachment",
                                        "color" => "#36a64f",
                                        "pretext" => "Optional text that appears above the attachment block",
                                        "author_name" => "Bobby Tables",
                                        "author_link" => "http://flickr.com/bobby/",
                                        "author_icon" => "http://flickr.com/icons/bobby.jpg",
                                        "title" => "Slack API Documentation",
                                        "title_link" => "https://api.slack.com/",
                                        "text" => "Optional text that appears within the attachment",
                                        "fields" => [Dict("title" => "Priority",
                                            "value" => "High",
                                            "short" => false
                                            )]
                                    )])

response = sendtoslack(data, endpoint)
@test response == "ok"

response2 = sendattachmenttoslack(attachment, endpoint)
@test response2 == "ok"

response3 = getchannels(token)
println(response3)


text:
using DataStructures

function main()
    N, M = parse.(Int, split(readline()))
    As, Bs = zeros(Int, M), zeros(Int, M)
    for i in 1:M
        As[i], Bs[i] = parse.(Int, split(readline()))
    end

    ans = 0
    for i in 1:M
        ds = IntDisjointSets(N)
        for j in 1:M
            if j ≠ i
                union!(ds, As[j], Bs[j])
            end
        end

        cnt = 0
        for j in 1:N
            if j == find_root(ds, j)
                cnt += 1
            end
        end

        if cnt > 1
            ans += 1
        end
    end

    println(ans)
end

main()


text:
## String to Char8 Array ##
st = "This is an ASCII string."
ar = c8"This is an ASCII array."
another = char8(st[12:16])

## Macth on Char8 Array ##
m = match(r"ASCII",ar)
dump(m)

## Change characters ##
swap(ar,'a','*')

## Using a IntSet as alphabet ##
vowels = IntSet('A','a','E','e','I','i','O','o','U','u')

idx = in(ar,vowels)
ar[ idx ]
check(ar,vowels)

## Swap using Dict ##
c8dic = (Char8=>Char8)[ char8('A')=>char8('4'), char('I')=>char8('1')]

swap!(ar,c8dic)

dic = ['A'=>1,'S'=>2,'C'=>3,'I'=>4]

intarray = swap(c8"ASCII Array",dic)


text:
using BinDeps
using Compat

vers = "1.8.1"

tagfile = "installed_vers"
target = "libblosc.$(Libdl.dlext)"
url = "https://bintray.com/artifact/download/julialang/generic/"
if !isfile(target) || !isfile(tagfile) || readchomp(tagfile) != "$vers $(Sys.WORD_SIZE)"
    if Compat.KERNEL == :NT
        run(download_cmd(url*"libblosc$(Sys.WORD_SIZE)-$vers.dll", target))
    elseif Compat.KERNEL == :Darwin
        run(download_cmd(url*"libblosc$(Sys.WORD_SIZE)-$vers.dylib", target))
    else
        tarball = "c-blosc-$vers.tar.gz"
        srcdir = "c-blosc-$vers/blosc"
        if !isfile(tarball)
            run(download_cmd("https://github.com/Blosc/c-blosc/archive/v$vers.tar.gz", tarball))
        end
        run(unpack_cmd(tarball, ".", ".gz", ".tar"))
        cd(srcdir) do
            println("Compiling libblosc...")
            run(`make -f ../../make.blosc LIB=../../$target`)
        end
    end
    open(tagfile, "w") do f
        println(f, "$vers $(Sys.WORD_SIZE)")
    end
end


text:
############################
#### Big-Leaf Physiology ###-------------------------------------------------------------------------------
############################

#' Bulk Intercellular CO2 Concentration
#'
#' Bulk canopy intercellular CO2 concentration (Ci) calculated based on Fick's law
#'              given surface conductance (Gs), gross primary productivity (GPP) and
#'              atmospheric CO2 concentration (Ca).
#'
#' - data             Data_Frame or matrix with all required columns
#' - Ca               Atmospheric or surface CO2 concentration (umol mol-1)
#' - GPP              Gross primary productivity (umol CO2 m-2 s-1)
#' - Gs               Surface conductance to water vapor (mol m-2 s-1)
#' - Rleaf            Ecosystem respiration stemming from leaves (umol CO2 m-2 s-1); defaults to 0
#' - missing_Rleaf_as_NA if Rleaf is provided, should missing values be treated as `missing` (`true`)
#'                            or set to 0 (`false`, the default)?
#' - constants        DwDc - Ratio of the molecular diffusivities for water vapor and CO2 (-)
#'
#' # Details
 Bulk intercellular CO2 concentration (Ci) is given by:
#'
#'            ``Ci = Ca - (GPP - Rleaf)/(Gs/1.6)``
#'
#'          where Gs/1.6 (mol m-2 s-1) represents the surface conductance to CO2.
#'          Note that Gs is required in mol m-2 s-1 for water vapor. Gs is converted to
#'          its value for CO2 internally.
#'          Ca can either be atmospheric CO2 concentration (as measured), or surface
#'          CO2 concentration as calculated from [`surface_CO2`](@ref).
#'
#' # Note
#' The equation is based on Fick's law of diffusion and is equivalent to the
#'       often used equation at leaf level (ci = ca - An/gs).
#'       Note that GPP and Gs have a different interpretation than An and gs.
#'       Gs comprises non-physiological contributions (i.e. physical evaporation)
#'       and is confounded by physical factors (e.g. energy balance non-closure).
#'       GPP does not account for dark respiration and is further subject to uncertainties
#'       in the NEE partitioning algorithm used. Leaf respiration can be provided,
#'       but it is usually not known at ecosystem level (as a consequence, Ci is likely to be
#'       slightly underestimated)
#'       This function should be used with care and the resulting Ci might not be
#'       readily comparable to its leaf-level analogue and/or physiological meaningful.
#'
#' # Value
 - Ci -: Bulk canopy intercellular CO2 concentration (umol mol-1)
#'
#' #References
#' Kosugi Y. et al., 2013: Determination of the gas exchange phenology in an
#'             evergreen coniferous forest from 7 years of eddy covariance flux data using
#'             an extended big-leaf analysis. Ecol Res 28, 373-385.
#'
#'             Keenan T., Sabate S., Gracia C., 2010: The importance of mesophyll conductance in
#'             regulating forest ecosystem productivity during drought periods. Global Change Biology
#'             16, 1019-1034.
#'
#' ```@example; output = false
#' ```
#' # calculate bulk canopy Ci of a productive ecosystem
#' intercellular_CO2(Ca=400,GPP=40,Gs=0.7)
#'
#' # note the sign convention for NEE
#'
"""
"""
function intercellular_CO2(data,Ca="Ca",GPP="GPP",Gs="Gs_mol",Rleaf=nothing,
                              missing_Rleaf_as_NA=false,constants=BigleafConstants())

  check_input(data,list(Ca,GPP,Gs))

  if(!isnothing(Rleaf))
    if(!missing_Rleaf_as_NA){Rleaf[ismissing(Rleaf)] = 0 }
else
    cat("Respiration from the leaves is ignored and set to 0.",fill=true)
    Rleaf = 0
end

  Ci = Ca - (GPP - Rleaf)/(Gs/constants.DwDc)

  return(Ci)

end

#' Bulk Canopy Photosynthetic Capacity (Vcmax and Jmax)
#'
#' Bulk canopy maximum carboxylation rate (Vcmax25), and maximum electron
#'              transport rate (Jmax25) at 25 degrees Celsius from bulk intercellular
#'              CO2 concentration using the Farquhar et al. 1980 model for C3 photosynthesis.
#'
#' - data      Data_Frame or matrix with all required columns
#' - C3        C3 vegetation (`true`, the default) or C4 vegetation (`false`)?
#' - Temp      Surface (or air) temperature (degC)
#' - GPP       Gross primary productivity (umol m-2 s-1)
#' - Ci        Bulk canopy intercellular CO2 concentration (umol mol-1)
#' - PPFD      Photosynthetic photon flux density (umol m-2 s-1)
#' - PPFD_j    PPFD threshold, below which the canopy is considered to
#'                  be RuBP regeneration limited. Defaults to 500 umol m-2 s-1.
#' - PPFD_c    PPFD threshold, above which the canopy is considered to
#'                  be Rubisco limited. Defaults to 1000 umol m-2 s-1.
#' - Rleaf     Ecosystem respiration stemming from leaves (umol CO2 m-2 s-1); defaults to 0
#' - Oi        Intercellular O2 concentration (mol mol-1)
#' - Kc25      Michaelis-Menten constant for CO2 at 25 degC (umol mol-1)
#' - Ko25      Michaelis-Menten constant for O2 at 25 degC (mmol mol-1)
#' - Gam25     Photorespiratory CO2 compensation point ('Gamma star')
#'                  at 25 degC (umol mol-1)
#' - Kc_Ha     Activation energy for Kc (kJ mol-1)
#' - Ko_Ha     Activation energy for Ko (kJ mol-1)
#' - Gam_Ha    Activation energy for Gam (kJ mol-1)
#' - Vcmax_Ha  Activation energy for Vcmax (kJ mol-1)
#' - Vcmax_Hd  Deactivation energy for Vcmax (kJ mol-1)
#' - Vcmax_dS  Entropy term for Vcmax (kJ mol-1 K-1)
#' - Jmax_Ha   Activation energy for Jmax (kJ mol-1)
#' - Jmax_Hd   Deactivation energy for Jmax (kJ mol-1)
#' - Jmax_dS   Entropy term for Jmax (kJ mol-1 K-1)
#' - Theta     Curvature term in the light response function of J (-)
#' - alpha_canopy Canopy absorptance (-)
#' - missing_Rleaf_as_NA if Rleaf is provided, should missing values be treated as `missing` (`true`)
#'                            or set to 0 (`false`, the default)?
#' - Ci_C4        intercellular CO2 concentration below which photosynthesis
#'                     is considered to be CO2-limited (umol mol-1), ignored
#'                     if `C3 = true`.
#' - constants    Kelvin - conversion degree Celsius to Kelvin
#'                     Rgas - universal gas constant (J mol-1 K-1)
#'                     kJ2J - conversion kilojoule (kJ) to joule (J)
#'                     J2kJ - conversion joule (J) to kilojoule (kJ)
#'                     se_median - conversion standard error (SE) of the mean to SE of the median
#'
#' # Details
 The maximum carboxylation rate at 25degC (Vcmax25) and the maximum electron
#'          transport rate at 25degC (Jmax25), which characterize photosynthetic capacity,
#'          are calculated as at leaf level.
#'          The required variables Gs and Ci can be calculated from
#'          [`surface_conductance`](@ref) and [`intercellular_CO2`](@ref), respectively.
#'
#'          Gas exchange parameters are taken from Bernacchi et al. 2001 (apparent values, which
#'          assume an infinite mesophyll conductance). Negative and very low Ci values
#'          (the threshold is set to Ci < 80umol mol-1 at the moment) are filtered out.
#'
#'          Vcmax is calculated from the photosynthesis model by Farquhar et al. 1980.
#'          If net photosynthesis is Rubisco-limited (RuBP-saturated carboxylation
#'          rate, i.e. light has to be (near-)saturating):
#'
#'            ``Vcmax = (GPP * (Ci + Kc*(1.0 + Oi/Ko))) / (Ci - Gam)``
#'
#'          where Kc and Ko are the Michaelis-Menten constants for CO2 and O2 (mmol mol-1),
#'          respectively, Oi is the O2 concentration, and Gam is the photorespiratory CO2
#'          compensation point (umol mol-1).
#'          Under low-light conditions, the electron transport rate J is calculated from
#'          the RuBP regeneration-limited photosynthesis rate:
#'
#'            ``J = (GPP * (4.0 * Ci + 8.0 * Gam) / (Ci - Gam)``
#'
#'          In this function, bulk canopy photosynthesis is assumed to be Rubisco/RuBP-regeneration
#'          limited, if incoming PPFD is above/below a specified threshold or range. These ranges
#'          are determined by the parameters `PPFD_j` and `PPFD_c`. If, for example,
#'          `PPFD_j = c(100,400)`, all conditions with a PPFD between 100 and 400 are assumed
#'          to be in the RuBP-regeneration (i.e. light-limited) photosynthesis domain. The electron
#'          transport rate J is then only calculated for periods that meet this criterion.
#'
#'          Jmax is calculated from J and absorbed irradiance:
#'
#'            \deqn{J = (APPFD_PSII + Jmax - sqrt((APPFD_PSII + Jmax)^2 -
#'                     4.0 * Theta * APPFD_PSII * Jmax)) / (2.0 * Theta)
#'                 }
#'
#'          where APPFD_PSII is the absorbed PPFD by photosystem II (PS II),
#'          and Theta is a curvature parameter. APPFD_PSII is calculated as
#'
#'            ``PPFD * alpha_canopy * 0.85 * beta``
#'
#'          where alpha_canopy is canopy-scale absorptance, 0.85 is a correction factor,
#'          and beta is the fraction of photons absorbed by PS II (assumed 0.5).
#'          alpha_canopy accounts for non-absorbing components of the ecosystem such as
#'          stems or soil, and is very likely ecosystem-specific. This parameter is relatively
#'          sensitive for the determination of Jmax25 at some sites.
#'
#'          Vcmax and Jmax at canopy level are assumed to follow the same temperature response
#'          as at leaf level. Hence, the respective parameter k at 25degC (k25) is calculated as
#'          (see e.g. Kattge & Knorr 2007):
#'
#'            \deqn{k25 = k /
#'                        ( exp(Ha * (Temp - Tref) / (Tref * Rgas * Temp)) *
#'                        (1 + exp((Tref * dS - Hd) / (Tref * Rgas))) /
#'                        (1 + exp((Temp * dS - Hd) / (Temp * Rgas)))
#'                        )
#'                  }
#'
#'          where Ha is the activation energy (kJ mol-1), Hd is the deactivation energy (kJ mol-1),
#'          and dS is the entropy term (kJ mol-1 K-1) of the respective parameter. Tref is set
#'          to 298.15 K.
#'
#'          For C4 photosynthesis, the simplified model by von Caemmerer 2000 is used.
#'          For light-saturated photosynthesis, Vcmax is given by:
#'
#'            ``Vcmax = GPP``
#'
#'          Note that in addition to the range `PPFD_c`, the range `Ci_C4`
#'          discards all periods with low Ci, in which photosynthesis is likely to
#'          be CO2-limited (see von Caemmerer 2000 for details).
#'
#'          In the light-limited case, J is calculated as:
#'
#'            ``J = 3 * GPPj / (1 - 0.5) ``
#'
#'          The calculation of Jmax25 and Vcmax25 is identical to C3 photosynthesis
#'          as described above.
#'
#' # Note
#'   The critical assumption is that bulk canopy photosynthesis is limited by
#'         one of the two limitation states. Incoming PPFD is assumed to determine
#'         the limitation states. Note however that the ranges (`PPFD_j` and `PPFD_c`)
#'         are likely ecosystem-specific. E_g. dense canopies presumably require higher
#'         `PPFD_c` thresholds than open canopies. A threshold of 500 umol m-2 s-1 PPFD
#'         for Rubisco-limited photosynthesis was assumed a reasonable working assumption (see Kosugi et al. 2013).
#'         Here, `PPFD_c` defaults to 1000 umol m-2 s-1. Note that even under very high/low irradiances,
#'         not all photosynthetically active plant material of an ecosystem will be in the same
#'         limitation state. Note that parameters describing bulk canopy photosynthetic capacity are not directly
#'         comparable to their leaf-level counterparts, as the former integrate over the entire canopy
#'         depth (i.e. are given per ground area, and not per leaf area).
#'         In general, the function should be used with care!
#'
#' # Value
 a DataFrame with the following columns:
#'         - Vcmax25: maximum bulk canopy carboxylation rate at 25degC (umol m-2 (ground) s-1)
#'         - Jmax25: maximum bulk canopy electron transport rate at 25degC (umol m-2 (ground) s-1)
#'
#' #References
#' Lloyd J. et al., 1995: A simple calibrated model of Amazon rainforest productivity
#'             based on leaf biochemical properties. Plant, Cell and Environment 18, 1129-1145.
#'
#'             Rayment M_B., Loustau D., Jarvis P_G., 2002: Photosynthesis and respiration
#'             of black spruce at three organizational scales: shoot, branch and canopy.
#'             Tree Physiology 22, 219-229.
#'
#'             Kosugi Y. et al., 2013: Determination of the gas exchange phenology in an
#'             evergreen coniferous forest from 7 years of eddy covariance flux data using
#'             an extended big-leaf analysis. Ecol Res 28, 373-385.
#'
#'             Ueyama M. et al, 2016: Optimization of a biochemical model with eddy covariance
#'             measurements in black spruce forests of Alaska for estimating CO2 fertilization
#'             effects. Agricultural and Forest Meteorology 222, 98-111.
#'
#'             Bernacchi C_J., Singsaas E_L., Pimentel C., Portis JR A_R., Long S_P., 2001:
#'             Improved temperature response functions for models of Rubisco-limited
#'             photosynthesis. Plant, Cell and Environment 24, 253-259.
#'
#'             Bernacchi C_J., Pimentel C., Long S_P., 2003: In vivo temperature response
#'             functions of parameters required to model RuBP-limited photosynthesis.
#'             Plant, Cell and Environment 26, 1419-1430.
#'
#'             von Caemmerer, 2000: Biochemical models of leaf photosynthesis. Techniques
#'             in plant sciences No. 2. CSIRO Publishing, Collingwood VIC, Australia.
#'
#' #See also
#' [`intercellular_CO2`](@ref), [`Arrhenius_temp_response`](@ref)
#'
#' ```@example; output = false
#' ```
#' DE_Tha_Jun_2014_2 = filter_data(DE_Tha_Jun_2014,quality_control=false,
#'                                  vars_qc=c("Tair","precip","VPD","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar","LE"),
#'                                  filter_vals_min=c(5,200,0.2,0),
#'                                  filter_vals_max=c(missing,missing,missing,missing),NA_as_invalid=true,
#'                                  quality_ext="_qc",good_quality=c(0,1),
#'                                  missing_qc_as_bad=true,GPP="GPP",doy="doy",
#'                                  year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
#'                                  tprecip=0.1,precip_hours=24,records_per_hour=2)
#'
#' # calculate Ga
#' Ga = aerodynamic_conductance(DE_Tha_Jun_2014_2,Gb_model=Val(:Thom_1972))[,"Ga_h"]
#'
#' # calculate Gs from the the inverted PM equation
#' Gs_PM = surface_conductance(DE_Tha_Jun_2014_2,Tair="Tair",pressure="pressure",
#'                              Rn="Rn",G="G",S=nothing,VPD="VPD",Ga=Ga,
#'                              formulation=Val(:PenmanMonteith))[,"Gs_mol"]
#'
#' # calculate Ci
#' Ci = intercellular_CO2(DE_Tha_Jun_2014_2,Ca="Ca",GPP="GPP",Gs=Gs_PM)
#'
#' # calculate Vcmax25 and Jmax25
#' photosynthetic_capacity(DE_Tha_Jun_2014_2,Temp="Tair",Ci=Ci,PPFD_j=c(200,500),PPFD_c=1000)
#'
#' @importFrom stats optimize
#'
#' @export
function photosynthetic_capacity(data,C3=true,Temp,GPP="GPP",Ci,PPFD="PPFD",PPFD_j=c(200,500),PPFD_c=1000,
                                    Rleaf=nothing,Oi=0.21,Kc25=404.9,Ko25=278.4,Gam25=42.75,
                                    Kc_Ha=79.43,Ko_Ha=36.38,Gam_Ha=37.83,Vcmax_Ha=65.33,Vcmax_Hd=200,
                                    Vcmax_dS=0.635,Jmax_Ha=43.9,Jmax_Hd=200,Jmax_dS=0.640,
                                    Theta=0.7,alpha_canopy=0.8,missing_Rleaf_as_NA=false,Ci_C4=100,
                                    constants=BigleafConstants())

  check_input(data,list(Temp,GPP,Ci,PPFD))

  Temp = Temp + constants.Kelvin
  Tref = 25.0 + constants.Kelvin

  if (C3){  # C3 vegetation
    Kc_Ha    = Kc_Ha * constants.kJ2J
    Ko_Ha    = Ko_Ha * constants.kJ2J
    Gam_Ha   = Gam_Ha * constants.kJ2J

    # Temperature dependencies of photosynthetic parameters
    Kc  = Kc25 * exp(Kc_Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp))
    Ko  = Ko25 * exp(Ko_Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp))
    Gam = Gam25 * exp(Gam_Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp))
    Ko  = Ko * constants.J2kJ

    # basic filtering on Ci
    Ci[Ci < 80 | ismissing(Ci)] = missing

    # Presumed limitation states
    GPPc = GPPj = GPP
    GPPj[PPFD < PPFD_j[1] | PPFD > PPFD_j[2] | ismissing(PPFD)] = missing
    GPPc[PPFD < PPFD_c | ismissing(PPFD)] = missing

    if(!isnothing(Rleaf))
      if(!missing_Rleaf_as_NA){Rleaf[ismissing(Rleaf)] = 0 }
else
      cat("Respiration from the leaves is ignored and set to 0.",fill=true)
      Rleaf = 0
end

    # calculate Vcmax and J (electron transport rate)
    Vcmax = (GPPc-Rleaf) * (Ci + Kc*(1.0 + Oi/Ko)) / (Ci - Gam)
    J     = (GPPj-Rleaf) * (4.0 * Ci + 8.0 * Gam) / (Ci - Gam)

else {  # C4 vegetation

    # Presumed limitation states (C4)
    GPPc = GPPj = GPP
    GPPj[PPFD < PPFD_j[1] | PPFD > PPFD_j[2] | ismissing(PPFD) | Ci < 0] = missing
    GPPc[PPFD < PPFD_c | Ci < Ci_C4 | ismissing(PPFD)] = missing

    Vcmax = GPPc
    J     = 3 * GPPj / (1 - 0.5)

end

  # calculate Jmax from J
  APPFD_PSII = PPFD * alpha_canopy * 0.85 * 0.5

  calcJmax  = which(complete_cases(J,APPFD_PSII))
  if (length(calcJmax) > 0)
    Jmax = sapply(calcJmax, function(i) tryCatch(optimize(function(Jmax){abs(J[i] - c((APPFD_PSII[i] + Jmax -
                                                                                sqrt((APPFD_PSII[i] + Jmax)^2 -
                                                                                4.0 * Theta * APPFD_PSII[i] * Jmax)) /
                                                                                (2.0 * Theta)))},
                                                           interval=c(0,1000),tol=1e-02)$minimum,
                                                  error=function(err){missing}
    )
    )
else
    @warn"Not enough observations to calculate Jmax!")
    Jmax = missing
end

  # calculate Vcmax25 and Jmax25
  Vcmax25 = Arrhenius_temp_response(Vcmax,Temp-constants.Kelvin,Ha=Vcmax_Ha,
                                     Hd=Vcmax_Hd,dS=Vcmax_dS,constants=constants)

  Jmax25 = Arrhenius_temp_response(Jmax,Temp[calcJmax]-constants.Kelvin,Ha=Jmax_Ha,
                                    Hd=Jmax_Hd,dS=Jmax_dS,constants=constants)

  # calculate medians and standard errors of the median
  Vcmax25_Median = median(Vcmax25,na_rm=true)
  Vcmax25_SE     = constants.se_median * sd(Vcmax25,na_rm=true)/sqrt((sum(!ismissing(Vcmax25))))
  Jmax25_Median  = median(Jmax25,na_rm=true)
  Jmax25_SE      = constants.se_median * sd(Jmax25,na_rm=true)/sqrt((sum(!ismissing(Jmax25))))

  return(c("Vcmax25"=round(Vcmax25_Median,2),"Vcmax25_SE"=round(Vcmax25_SE,2),
           "Jmax25"=round(Jmax25_Median,2),"Jmax25_SE"=round(Jmax25_SE,2)))

end

#' (Modified) Arrhenius Temperature Response Function
#'
#' (Modified) Arrhenius function describing
#'              the temperature response of biochemical parameters.
#'
#' - param Parameter measured at measurement temperature (umol m-2 s-1)
#' - Temp  Measurement temperature (degC)
#' - Ha    Activation energy for param (kJ mol-1)
#' - Hd    Deactivation energy for param (kJ mol-1)
#' - dS    Entropy term for param (kJ mol-1 K-1)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  Rgas - universal gas constant (J mol-1 K-1)
#'                  kJ2J - conversion kilojoule (kJ) to joule (J)
#'
#' # Details
 The function returns the biochemical rate at a reference
#'          temperature of 25degC given a predefined temperature response function.
#'          This temperature response is given by a modified form of the Arrhenius
#'          function:
#'
#'             \deqn{param25 = param /
#'                         ( exp(Ha * (Temp - Tref) / (Tref*Rgas*Temp)) *
#'                         (1 + exp((Tref*dS - Hd) / (Tref * Rgas))) /
#'                         (1 + exp((Temp*dS - Hd) / (Temp * Rgas)))
#'                         )
#'                  }
#'
#'          where param is the value/rate of the parameter at measurement temperature,
#'          Temp is temperature in K, Tref is reference temperature (298.15K), and Rgas
#'          is the universal gas constant (8.314 J K-1 mol-1). Ha is the activation
#'          energy (kJ mol-1), Hd is the deactivation energy (kJ mol-1), and dS the
#'          entropy term (kJ mol-1 K-1) of the respective parameter.
#'
#'          If either Hd or dS or both are not provided, the equation above reduces
#'          to the first term (i.e. the common Arrhenius equation without the deactivation
#'          term.)
#'
#' # Value
 param25 - value of the input parameter at the reference temperature of 25degC (umol m-2 s-1)
#'
#' #References
#' Johnson F_H., Eyring H., Williams R_W. 1942:
#'             The nature of enzyme inhibitions in bacterial luminescence: sulfanilamide,
#'             urethane, temperature and pressure. Journal of cellular and comparative
#'             physiology 20, 247-268.
#'
#'             Kattge J., Knorr W., 2007: Temperature acclimation in a biochemical
#'             model of photosynthesis: a reanalysis of data from 36 species.
#'             Plant, Cell and Environment 30, 1176-1190.
#'
#' @export
function Arrhenius_temp_response(param,Temp,Ha,Hd,dS,constants=BigleafConstants())

  Temp = Temp + constants.Kelvin
  Tref = 25.0 + constants.Kelvin

  Ha = ifelse(missing(Ha),missing,Ha*constants.kJ2J)
  Hd = ifelse(missing(Hd),missing,Hd*constants.kJ2J)
  dS = ifelse(missing(dS),missing,dS*constants.kJ2J)

  if (ismissing(Ha))

    stop("Activation energy (Ha) has to be provided!")

end

  if (ismissing(Hd) & ismissing(dS))

    param25 = param / exp(Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp))

elseif (!ismissing(Hd) & !ismissing(dS))

    param25 = param /
      ( exp(Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp)) *
        (1 + exp((Tref*dS - Hd) / (Tref * constants.Rgas))) /
        (1 + exp((Temp*dS - Hd) / (Temp * constants.Rgas)))
      )

elseif ((!ismissing(Hd) & ismissing(dS)) | (ismissing(Hd) & !ismissing(dS)) )

    @warn"Both Hd and dS have to be provided for a temperature response
             that considers a temperature optimum and a deactivation term!
             Continue considering activation energy (Ha) only...")

    param25 = param / exp(Ha * (Temp - Tref) / (Tref*constants.Rgas*Temp))

end

  return(param25)

end

#' Stomatal Slope Parameter "g1"
#'
#' Estimation of the intrinsic WUE metric "g1" (stomatal slope)
#'              from nonlinear regression.
#'
#' - data       DataFrame or matrix containing all required columns
#' - Tair       Air (or surface) temperature (deg C)
#' - pressure   Atmospheric pressure (kPa)
#' - GPP        Gross primary productivity (umol CO2 m-2 s-1)
#' - Gs         Surface conductance to water vapor (mol m-2 s-1)
#' - VPD        Vapor pressure deficit (kPa)
#' - Ca         Atmospheric CO2 concentration (air or surface) (umol mol-1)
#' - Rleaf      Ecosystem respiration stemming from leaves (umol CO2 m-2 s-1); defaults to 0
#' - model      Stomatal model used. One of `"USO","Ball&Berry","Leuning"`.
#' - robust_nls Use robust nonlinear regression (`\link[robustbase]{nlrob`})? Default is `false`.
#' - nmin       Minimum number of data required to perform the fit; defaults to 40.
#' - fitg0      Should g0 and g1 be fitted simultaneously?
#' - g0         Minimum stomatal conductance (mol m-2 s-1); ignored if `fitg0 = true`.
#' - fitD0      Should D0 be fitted along with g1 (and g0 if `fitg0 = true`)?; only used if `model = "Leuning"`.
#' - D0         Stomatal sensitivity parameter to VPD; only used if `model = "Leuning"` and `fitD0 = false`.
#' - Gamma      Canopy CO2 compensation point (umol mol-1); only used if `model = "Leuning"`.
#'                   Can be a constant or a variable. Defaults to 50 umol mol-1.
#' - constants  Kelvin - conversion degree Celsius to Kelvin
#'                   Rgas - universal gas constant (J mol-1 K-1)
#'                   DwDc - Ratio of the molecular diffusivities for water vapor and CO2
#' - missing_Rleaf_as_NA if Rleaf is provided, should missing values be treated as `missing` (`true`)
#'                            or set to 0 (`false`, the default)?
#' - ...        Additional arguments to `\link[stats]{nls`} or `\link[robustbase]{nlrob`} if `robust_nls = true`.
#'
#' # Details
 All stomatal models were developed at leaf-level, but its parameters
#'          can also be estimated at ecosystem level (but be aware of caveats).
#'
#'          The unified stomatal optimization (USO) model is given by (Medlyn et al. 2011):
#'
#'             ``gs = g0 + 1.6*(1.0 + g1/sqrt(VPD)) * An/ca``
#'
#'          The semi-empirical model by Ball et al. 1987 is defined as:
#'
#'             ``gs = g0 + g1* ((An * rH) / ca)``
#'
#'          Leuning 1995 suggested a revised version of the Ball&Berry model:
#'
#'             ``gs = g0 + g1*An / ((ca - \\gamma) * (1 + VPD/D0))``
#'
#'          where ``\\gamma`` is by default assumed to be constant, but likely varies with temperature and among
#'          plant species.
#'          The equations above are valid at leaf-level. At ecosystem level, An is replaced by GPP (or GPP - Rleaf,
#'          where Rleaf is leaf respiration), and gs (stomatal conductance) by Gs (surface conductance).
#'          The parameters in the models are estimated using nonlinear regression (`\link[stats]{nls`}) if
#'          `robust_nls = false` and weighted nonlinear regression if `robust_nls = true`.
#'          The weights are calculated from `\link[robustbase]{nlrob`}, and `\link[stats]{nls`}
#'          is used for the actual fitting.
#'          Alternatively to measured VPD and Ca (i.e. conditions at instrument height), conditions at
#'          the big-leaf surface can be provided. Those can be calculated using [`surface_conditions`](@ref).
#'
#'
#' # Value
 A `nls` model object, containing information on the fitted parameters, their uncertainty range,
#'         model fit, etc.
#'
#' #References
#' Medlyn B_E., et al., 2011: Reconciling the optimal and empirical approaches to
#'             modelling stomatal conductance. Global Change Biology 17, 2134-2144.
#'
#'             Ball T_J., Woodrow I_E., Berry J_A. 1987: A model predicting stomatal conductance
#'             and its contribution to the control of photosynthesis under different environmental conditions.
#'             In: Progress in Photosynthesis Research, edited by J_Biggins, pp. 221-224, Martinus Nijhoff Publishers,
#'             Dordrecht, Netherlands.
#'
#'             Leuning R., 1995: A critical appraisal of a combined stomatal-photosynthesis
#'             model for C3 plants. Plant, Cell and Environment 18, 339-355.
#'
#'             Knauer, J. et al., 2018: Towards physiologically meaningful water-use efficiency estimates
#'             from eddy covariance data. Global Change Biology 24, 694-710.
#'
#' #See also
#' [`surface_conductance`](@ref)
#'
#' ```@example; output = false
#' ```
#' ## filter data to ensure that Gs is a meaningful proxy to canopy conductance (Gc)
#' DE_Tha_Jun_2014_2 = filter_data(DE_Tha_Jun_2014,quality_control=false,
#'                                  vars_qc=c("Tair","precip","VPD","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar","LE"),
#'                                  filter_vals_min=c(5,200,0.2,0),
#'                                  filter_vals_max=c(missing,missing,missing,missing),NA_as_invalid=true,
#'                                  quality_ext="_qc",good_quality=c(0,1),
#'                                  missing_qc_as_bad=true,GPP="GPP",doy="doy",
#'                                  year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
#'                                  tprecip=0.1,precip_hours=24,records_per_hour=2)
#'
#' # calculate Gs from the the inverted PM equation
#' Ga = aerodynamic_conductance(DE_Tha_Jun_2014_2,Gb_model=Val(:Thom_1972))[,"Ga_h"]
#'
#' # if G and/or S are available, don't forget to indicate (they are ignored by default).
#' Gs_PM = surface_conductance(DE_Tha_Jun_2014_2,Tair="Tair",pressure="pressure",
#'                              Rn="Rn",G="G",S=nothing,VPD="VPD",Ga=Ga,
#'                              formulation=Val(:PenmanMonteith))[,"Gs_mol"]
#'
#' ### Estimate the stomatal slope parameter g1 using the USO model
#' mod_USO = stomatal_slope(DE_Tha_Jun_2014_2,model="USO",GPP="GPP",Gs=Gs_PM,
#'                           robust_nls=false,nmin=40,fitg0=false)
#'
#' ### Use robust regression to minimize influence of outliers in Gs
#' mod_USO = stomatal_slope(DE_Tha_Jun_2014_2,model="USO",GPP="GPP",Gs=Gs_PM,
#'                           robust_nls=true,nmin=40,fitg0=false)
#'
#' ### Estimate the same parameter from the Ball&Berry model and prescribe g0
#' mod_BB = stomatal_slope(DE_Tha_Jun_2014_2,model="Ball&Berry",GPP="GPP",
#'                          robust_nls=false,Gs=Gs_PM,g0=0.01,nmin=40,fitg0=false)
#'
#' ## same for the Leuning model, but this time estimate both g1 and g0 (but fix D0)
#' mod_Leu = stomatal_slope(DE_Tha_Jun_2014_2,model="Leuning",GPP="GPP",Gs=Gs_PM,
#'                           robust_nls=false,nmin=40,fitg0=false,D0=1.5,fitD0=false)
#'
#' @importFrom stats nls na_exclude
#' @importFrom robustbase nlrob
#'
#' @export
function stomatal_slope(data,Tair="Tair",pressure="pressure",GPP="GPP",Gs="Gs_mol",
                           VPD="VPD",Ca="Ca",Rleaf=nothing,model=c("USO","Ball&Berry","Leuning"),
                           robust_nls=false,nmin=40,fitg0=false,g0=0,fitD0=false,
                           D0=1.5,Gamma=50,missing_Rleaf_as_NA=false,
                           constants=BigleafConstants(),...)

  model = match_arg(model)

  check_input(data,list(Tair,pressure,GPP,Gs,VPD,Ca))

  df   = DataFrame(Tair,pressure,GPP,Gs,VPD,Ca)
  DwDc = constants.DwDc  # ...to work within nls()

  if (model == "Leuning")
    check_input(data,Gamma)
    df$Gamma = Gamma
end

  if(!isnothing(Rleaf))
    if(!missing_Rleaf_as_NA){Rleaf[ismissing(Rleaf)] = 0 }
else
    cat("Respiration from the leaves is ignored and set to 0.",fill=true)
    Rleaf = 0
end

  GPP = (GPP - Rleaf)

  if (model == "Leuning")
    nr_data = sum(!ismissing(GPP) & !ismissing(Gs) & !ismissing(VPD) & !ismissing(Ca) & !ismissing(Gamma))
else
    nr_data = sum(!ismissing(GPP) & !ismissing(Gs) & !ismissing(VPD) & !ismissing(Ca))
end

  if (nr_data < nmin)
    stop("number of data is less than 'nmin'. g1 is not fitted to the data.")
else

    if (model == "USO")

      if (fitg0)
        if (robust_nls)
          df$DwDc = rep(DwDc,nrow(df))
          mod_weights = nlrob(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,data=df,start=list(g0=0,g1=3),
                               na_action=na_exclude,...)$w
          mod = nls(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,start=list(g0=0,g1=3),weights=mod_weights,...)
else
          mod = nls(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,start=list(g0=0,g1=3),...)
end
else
        if (robust_nls)
          df$g0   = rep(g0,nrow(df))    # g0 as constant does not work in the nlrob function...
          df$DwDc = rep(DwDc,nrow(df))  # same with constants.DwDc
          mod_weights = nlrob(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,data=df,start=list(g1=3),
                               na_action=na_exclude,...)$w
          mod = nls(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,start=list(g1=3),weights=mod_weights,...)
else
          mod = nls(Gs ~ g0 + DwDc*(1.0 + g1/sqrt(VPD))*GPP/Ca,start=list(g1=3),...)
end
end

elseif (model == "Leuning")

      if (fitg0)
        if (fitD0)
          if (robust_nls)
            mod_weights = nlrob(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),data=df,
                                 start=list(g0=0,g1=9,D0=1.5),na_action=na_exclude,...)$w
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g0=0,g1=9,D0=1.5),
                       weights=mod_weights,...)
else
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g0=0,g1=9,D0=1.5),...)
end
else
          if (robust_nls)
            df$D0  = rep(D0,nrow(df))
            mod_weights = nlrob(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),data=df,
                                 start=list(g0=0,g1=9),na_action=na_exclude,...)$w
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g0=0,g1=9),
                       weights=mod_weights,...)
else
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g0=0,g1=9),...)
end
end
else
        if (fitD0)
          if (robust_nls)
            df$g0    = rep(g0,nrow(df))
            mod_weights = nlrob(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),data=df,
                                 start=list(g1=9,D0=1.5),na_action=na_exclude,...)$w
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g1=9,D0=1.5),
                       weights=mod_weights,...)
else
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g1=9,D0=1.5),...)
end
else
          if (robust_nls)
            df$g0  = rep(g0,nrow(df))
            df$D0  = rep(D0,nrow(df))
            mod_weights = nlrob(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),data=df,
                                 start=list(g1=9),na_action=na_exclude,...)$w
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g1=9),
                       weights=mod_weights,...)
else
            mod = nls(Gs ~ g0 + g1*GPP / ((Ca - Gamma) * (1 + VPD/D0)),start=list(g1=9),...)
end
end
end

elseif (model == "Ball&Berry")

      rH = VPD_to_rH(VPD,Tair)
      df$rH = rH

      if (fitg0)
        if (robust_nls)
          mod_weights = nlrob(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g0=0,g1=9),data=df,
                               na_action=na_exclude,...)$w
          mod = nls(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g0=0,g1=9),weights=mod_weights,...)
else
          mod = nls(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g0=0,g1=9),...)
end
else
        if (robust_nls)
          df$g0   = rep(g0,nrow(df))
          mod_weights = nlrob(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g1=9),data=df,
                               na_action=na_exclude,...)$w
          mod = nls(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g1=9),weights=mod_weights,...)
else
          mod = nls(Gs ~ g0 + g1 * (GPP * rH) / Ca,start=list(g1=9),...)
end
end

end

end

  return(mod)

end

#' Ecosystem Light Response
#'
#' calculates GPP_ref at a reference (usually saturating) PPFD and
#'              ecosystem quantum yield (alpha) using a rectangular light response curve.
#'
#' - data      DataFrame or matrix containing all required columns
#' - NEE       Net ecosystem exchange (umol CO2 m-2 s-1)
#' - Reco      Ecosystem respiration (umol CO2 m-2 s-1)
#' - PPFD      Photosynthetic photon flux density (umol m-2 s-1)
#' - PPFD_ref  Reference PPFD (umol m-2 s-1) for which GPP_ref is estimated.
#'                  Default is 2000 umol m-2 s-1.
#' - ...       Additional arguments to `\link[stats]{nls`}
#'
#' # Details
 A rectangular light response curve is fitted to NEE data. The curve
#'          takes the form as described in Falge et al. 2001:
#'
#'             \deqn{-NEE = \\alpha PPFD / (1 - (PPFD / PPFD_ref) + \\alpha
#'                          PPFD / GPP_ref) - Reco}
#'
#'          where ``\\alpha`` is the ecosystem quantum yield (umol CO2 m-2 s-1) (umol quanta m-2 s-1)-1,
#'          and GPP_ref is the GPP at the reference PPFD (usually at saturating light). ``\\alpha``
#'          represents the slope of the light response curve, and is a measure for the light use
#'          efficiency of the canopy.
#'
#'          The advantage of this equation over the standard rectangular light response
#'          curve is that GPP_ref at PPFD_ref is more readily interpretable
#'          as it constitutes a value observed in the ecosystem, in contrast to
#'          GPP_ref (mostly named 'beta') in the standard model that occurs at infinite light.
#'          `PPFD_ref` defaults to 2000 umol m-2 s-1, but other values can be used. For
#'          further details refer to Falge et al. 2001.
#'
#' # Note
#'   Note the sign convention. Negative NEE indicates that carbon is taken up
#'         by the ecosystem. Reco has to be 0 or larger.
#'
#' # Value
 A `nls` model object containing estimates (+/- SE) for alpha and GPP_ref.
#'
#' #References
#' Falge E., et al. 2001: Gap filling strategies for defensible annual
#'             sums of net ecosystem exchange. Agricultural and Forest Meteorology 107,
#'             43-69.
#'
#'             Gilmanov T_G., et al. 2003: Gross primary production and light response
#'             parameters of four Southern Plains ecosystems estimated using long-term
#'             CO2-flux tower measurements. Global Biogeochemical Cycles 17, 1071.
#'
#'             Reichstein M., Stoy P_C., Desai A_R., Lasslop G., Richardson A. 2012:
#'             Partitioning of net fluxes. In: Eddy Covariance. A practical guide to
#'             measurement and data analysis. Aubinet M., Vesala T., Papale D. (Eds.).
#'             Springer.
#'
#' @importFrom stats nls
#'
"""
"""
function light_response(data,NEE="NEE",Reco="Reco",PPFD="PPFD",PPFD_ref=2000,...)

  check_input(data,list(NEE,Reco,PPFD))

  mod = nls(-NEE ~ alpha * PPFD / (1 - (PPFD / PPFD_ref) + (alpha * PPFD / GPP_ref)) - Reco,
             start=list(alpha=0.05,GPP_ref=30),...)

  return(mod)
end

#' Light-Use Efficiency (LUE)
#'
#' Amount of carbon fixed (GPP) per incoming light.
#'
#' - GPP     Gross ecosystem productivity (umol CO2 m-2 s-1)
#' - PPFD    Photosynthetic photon flux density (umol quanta m-2 s-1)
#'
#' # Details
 Light use efficiency is calculated as
#'
#'             ``LUE = sum(GPP)/sum(PPFD)``
#'
#'          where both GPP and PPFD are in umol m-2 s-1. A more meaningful
#'          (as directly comparable across ecosystems) approach is to take
#'          absorbed PPFD rather than incoming PPFD as used here.
#'
#' # Value
 - LUE -: Light use efficiency (-)
#'
#' #See also
#' [`energy_use_efficiency`](@ref)
#'
#' ```@example; output = false
#' ```
#' light_use_efficiency(GPP=20,PPFD=1500)
#'
#' @export
function light_use_efficiency(GPP,PPFD)

  comp = complete_cases(GPP,PPFD)

  LUE = sum(GPP[comp],na_rm=true)/sum(PPFD[comp],na_rm=true)

  return(c("LUE"=LUE))
end

#' Stomatal Sensitivity to VPD
#'
#' Sensitivity of surface conductance to vapor pressure deficit.
#'
#' - data  DataFrame or matrix containing all required columns
#' - Gs    Surface conductance to water vapor (mol m-2 s-1)
#' - VPD   Vapor pressure deficit (kPa)
#' - ...   Additional arguments to `\link[stats]{nls`}
#'
#' # Details
 The function fits the following equation (Oren et al. 1999):
#'
#'             ``Gs = -m ln(VPD) + b``
#'
#'          where b is the reference surface conductance (Gs) at VPD=1kPa (in mol m-2 s-1),
#'          and m is the sensitivity parameter of Gs to VPD (in mol m-2 s-1 log(kPa)-1).
#'          The two parameters b and m are fitted using `\link[stats]{nls`}.
#'          VPD can be the one directly measured at instrument height, or the
#'          one at the surface, as returned by [`surface_conditions`](@ref).
#'
#' # Value
 A `nls` model object containing (amongst others) estimates for the mean
#'         and standard errors of the parameters m and b.
#'
#' #References
#' Oren R., et al. 1999: Survey and synthesis of intra- and interspecific
#'             variation in stomatal sensitivity to vapour pressure deficit. Plant,
#'             Cell & Environment 22, 1515-1526.
#'
#'             Novick K_A., et al. 2016: The increasing importance of atmospheric demand
#'             for ecosystem water and carbon fluxes. Nature Climate Change 6, 1023 - 1027.
#'
#' #See also
#' [`surface_conductance`](@ref)
#'
#' ```@example; output = false
#' ```
#' ## calculate Ga, Gs, and the stomatal sensitivity to VPD for the site FR-Pue in
#' ## May 2012. Data are filtered for daytime, sufficiently high ustar, etc.
#' FR_Pue_May_2012_2 = filter_data(FR_Pue_May_2012,quality_control=true,
#'                                  vars_qc=c("Tair","precip","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar","VPD"),
#'                                  filter_vals_min=c(5,200,0.2,0.3),
#'                                  filter_vals_max=c(missing,missing,missing,missing),
#'                                  NA_as_invalid=true,quality_ext="_qc",
#'                                  good_quality=c(0,1),missing_qc_as_bad=true,
#'                                  precip="precip",tprecip=0.1,precip_hours=24,
#'                                  records_per_hour=2)
#' Ga = aerodynamic_conductance(FR_Pue_May_2012_2)
#' Gs = surface_conductance(FR_Pue_May_2012_2,Ga=Ga[,"Ga_h"])
#' stomatal_sensitivity(FR_Pue_May_2012_2,Gs=Gs[,"Gs_mol"],VPD="VPD")
#'
#' @importFrom stats nls
#'
"""
"""
function stomatal_sensitivity(data,Gs="Gs_mol",VPD="VPD",...)

  check_input(data,list(Gs,VPD))

  mod = nls(Gs ~ -m * log(VPD) + b,start=list(m=0.05,b=0.2),...)

  return(mod)
end


text:
using LowRankModels

A = sprandn(100,100,.5)
glrm = qpca(A, 3)

fit!(glrm)
A_sampled = sample(glrm)
A_sample_missing = sample_missing(glrm)

# tests
obs = !(A.==0)
@assert (A[obs]==A_sample_missing[obs])


text:
using DualNumbers

import Base: push!, length, show, getindex, setindex!, eachindex, isassigned,
             isapprox, zero, one, lastindex

export Leaf, Tape, Node, Branch, ∇

""" Basic unit on the computational graph."""
abstract type Node{T} end

""" A topologically ordered collection of Nodes. """
struct Tape
    tape::Vector{Any}
    Tape() = new(Vector{Any}())
    Tape(N::Int) = new(Vector{Any}(undef, N))
end
function show(io::IO, t::Tape)
    n = length(t)
    print(io, "Tape with ", n, " element", n == 1 ? "" : "s", n > 0 ? ":" : "")
    for i in eachindex(t)
        print(io, "\n  [", i, "]: ")
        if isassigned(tape(t), i)
            show(io, t[i])
        else
            print(io, "#undef")
        end
    end
end
@inline getindex(t::Tape, n::Int) = getindex(tape(t), n)
@inline getindex(t::Tape, node::Node) = getindex(t, pos(node))
@inline lastindex(t::Tape) = length(t)
@inline setindex!(t::Tape, x, n::Int) = (tape(t)[n] = x; t)
@inline eachindex(t::Tape) = eachindex(tape(t))
@inline length(t::Tape) = length(tape(t))
@inline push!(t::Tape, node::Node) = (push!(tape(t), node); t)
@inline isassigned(t::Tape, n::Int) = isassigned(tape(t), n)
@inline isassigned(t::Tape, node::Node) = isassigned(t, pos(node))

# Make `Tape`s broadcast as scalars without a warning on 0.7
Base.Broadcast.broadcastable(tape::Tape) = Ref(tape)

"""
An element at the 'bottom' of the computational graph.

Fields:
val - the value of the node.
tape - The Tape to which this Leaf is assigned.
pos - the location of this Leaf in the tape to which it is assigned.
"""
struct Leaf{T} <: Node{T}
    val::T
    tape::Tape
    pos::Int
end
function Leaf(tape::Tape, val)
    leaf = Leaf(val, tape, length(tape) + 1)
    push!(tape, leaf)
    return leaf
end
show(io::IO, tape::Leaf{T}) where T = print(io, "Leaf{$T} $(unbox(tape))")
show(io::IO, tape::Leaf{T}) where T<:AbstractArray = print(io, "Leaf{$T} $(size(unbox(tape)))")

"""
A Branch is a Node with parents (args).

Fields:
val - the value of this node produced in the forward pass.
f - the function used to generate this Node.
args - Values indicating which elements in the tape will require updating by this node.
tape - The Tape to which this Branch is assigned.
pos - the location of this Branch in the tape to which it is assigned.
"""
struct Branch{T} <: Node{T}
    val::T
    f
    args::Tuple
    kwargs::NamedTuple
    tape::Tape
    pos::Int
end
function Branch(f, args::Tuple, tape::Tape; kwargs...)
    unboxed = unbox.(args)
    branch = Branch(f(unboxed...; kwargs...), f, args, kwargs.data, tape, length(tape) + 1)
    push!(tape, branch)
    return branch
end
show(io::IO, branch::Branch{T}) where T =
    print(io, "Branch{$T} $(unbox(branch)) f=$(getfield(branch, :f))")
show(io::IO, branch::Branch{T}) where T<:AbstractArray =
    print(io, "Branch{$T} $(size(unbox(branch))) f=$(getfield(branch, :f))")

"""
    tape(x::Node)
    tape(x::Tape)

Retrieve the `Tape` in a `Node`, or the underyling vector in a `Tape`.
"""
tape(x::Union{Node,Tape}) = getfield(x, :tape)

"""
    pos(x::Node)
    pos(x)

Location of Node on tape. -1 if not a Node object.
"""
pos(x::Node) = getfield(x, :pos)
pos(x) = -1

"""
    unbox(x::Node)
    unbox(x)

Get `.val` if `x` is a Node, otherwise is equivalent to `identity`.
"""
unbox(x::Node) = getfield(x, :val)
unbox(x) = x

isapprox(n::Node, f) = unbox(n) ≈ f
isapprox(f, n::Node) = n ≈ f
isapprox(n::Node, f::Node) = unbox(n) ≈ unbox(f)

zero(n::Node) = zero(unbox(n))
one(n::Node) = one(unbox(n))

# Leafs do nothing, Branches compute their own sensitivities and update others.
@inline propagate(y::Leaf, rvs_tape::Tape) = nothing
function propagate(y::Branch, rvs_tape::Tape)
    tape = Nabla.tape(rvs_tape)
    ȳ, f = tape[pos(y)], getfield(y, :f)
    args = getfield(y, :args)
    kwargs = getfield(y, :kwargs)
    xs, xids = map(unbox, args), map(pos, args)
    p = preprocess(f, unbox(y), ȳ, xs...)
    for j in eachindex(xs)
        x, xid = xs[j], xids[j]
        if xid > 0
            tape[xid] = isassigned(tape, xid) ?
                ∇(tape[xid], f, Arg{j}, p, unbox(y), ȳ, xs...; kwargs...) :
                ∇(f, Arg{j}, p, unbox(y), ȳ, xs...; kwargs...)
        end
    end
    return nothing
end

function propagate(fwd_tape::Tape, rvs_tape::Tape)
    for n in eachindex(rvs_tape)
        δ = length(rvs_tape) - n + 1
        isassigned(tape(rvs_tape), δ) && propagate(fwd_tape[δ], rvs_tape)
    end
    return rvs_tape
end

""" Initialise a Tape appropriately for being used as a reverse-tape. """
function reverse_tape(y::Node, ȳ)
    tape = Tape(pos(y))
    tape[end] = ȳ
    return tape
end

""" Used to flag which argument is being specified in x̄. """
struct Arg{N} end

"""
    ∇(y::Node{<:∇Scalar})
    ∇(y::Node{T}, ȳ::T) where T

Return a `Tape` object which can be indexed using `Node`s, each element of which contains
the result of multiplying `ȳ` by the transpose of the Jacobian of the function specified by
the `Tape` object in `y`. If `y` is a scalar and `ȳ = 1` then this is equivalent to
computing the gradient of `y` w.r.t. each of the elements in the `Tape`.

    ∇(f::Function, ::Type{Arg{N}}, p, y, ȳ, x...)
    ∇(x̄, f::Function, ::Type{Arg{N}}, p, y, ȳ, x...)

To implement a new reverse-mode sensitivity for the `N^{th}` argument of function `f`. p
is the output of `preprocess`. `x1`, `x2`,... are the inputs to the function, `y` is its
output and `ȳ` the reverse-mode sensitivity of `y`.
"""
∇(y::Node, ȳ) = propagate(tape(y), reverse_tape(y, ȳ))
@inline ∇(y::Node{<:∇Scalar}) = ∇(y, one(unbox(y)))

# This is a fallback method where we don't necessarily know what we'll be adding and whether
# we can update the value in-place, so we'll try to be clever and dispatch.
@inline function ∇(x̄, f, ::Type{Arg{N}}, args...; kwargs...) where N
    return update!(x̄, ∇(f, Arg{N}, args...; kwargs...))
end

# Update regular arrays in-place. Structured array types should not be updated in-place,
# even though it technically "works" (https://github.com/JuliaLang/julia/issues/31674),
# so we'll only permit mutating addition for `Array`s, e.g. `Vector` and `Matrix`.
# Mixed array and scalar adds should not occur, as sensitivities should always have the
# same shape, so we won't bother allowing e.g. updating an array with a scalar on the RHS.
update!(x̄::Array{T,N}, y::AbstractArray{S,N}) where {T,S,N} = x̄ .+= y

# Fall back to using regular addition
update!(x̄, y) = x̄ + y

"""
    ∇(f; get_output::Bool=false)

Returns a function which, when evaluated with arguments that are accepted by `f`, will
return the gradient w.r.t. each of the arguments. If `get_output` is `true`, the result
of calling `f` on the given arguments is also returned.
"""
function ∇(f; get_output::Bool=false)
    return function(args...; kwargs...)
        args_ = Leaf.(Tape(), args)
        y = f(args_...; kwargs...)
        if y isa Node
            ∇f = ∇(y)
            ∇args = map(args_, args) do arg_, arg
                isassigned(∇f, arg_) ? ∇f[arg_] : zero(arg)
            end
        else
            ∇args = zero.(args)
        end
        return get_output ? (y, ∇args) : ∇args
    end
end

# """
#     ∇(f::Function)

# Returns a function which, when evaluated with arguments that are accepted by `f` (`x`),
# will return a Tuple, the first element of which is the output of the function `f` and then
# second element of which is (yet another) function `g`. `g` can either be evaluated with no
# arguments, in which case it will return the gradient of `f` evaluated at `x`.
# Alternatively, it can be evaluated with arguments of the same type and shape as the output
# of `f(x)`, in which case it is equivalent to multiplying them 'from the left' by the
# Jacobian ∂(f(x)) / ∂x.
# """
# function ∇(f::Function)
#     return function(args...)
#         args_ = Leaf.(Tape(), args)
#         y = f(args_...)
#         ∇fx = (ȳ)->∇

#     end
# end

# A collection of methods for initialising nested indexable containers to zero.
for (f_name, scalar_init, array_init) in
    zip((:zerod_container, :oned_container, :randned_container),
        (:zero, :one, nothing),
        (:zeros, :ones, nothing))
    if scalar_init !== nothing
        @eval @inline $f_name(x::Number) = $scalar_init(x)
    end
    if array_init !== nothing
        @eval @inline $f_name(x::AbstractArray{<:Real}) = $array_init(eltype(x), size(x))
    end
    eval(quote
        @inline $f_name(x::Tuple) = map($f_name, x)
        @inline function $f_name(x)
            y = Base.copy(x)
            for n in eachindex(y)
                @inbounds y[n] = $f_name(y[n])
            end
            return y
        end
        $f_name(x::Ref) = Ref($f_name(x[]))
    end)
end
@inline randned_container(x::Number) = randn(typeof(x))
@inline randned_container(x::AbstractArray{<:Real}) = randn(eltype(x), size(x)...)
for T in (:Diagonal, :UpperTriangular, :LowerTriangular)
    @eval @inline randned_container(x::$T{<:Real}) = $T(randn(eltype(x), size(x)...))
end

# Bare-bones FMAD implementation based on DualNumbers. Accepts a Tuple of args and returns
# a Tuple of gradients. Currently scales almost exactly linearly with the number of inputs.
# The coefficient of this scaling could be improved by implementing a version of DualNumbers
# which computes from multiple seeds at the same time.
function dual_call_expr(f, x::Type{<:Tuple}, ::Type{Type{Val{n}}}) where n
    dual_call = Expr(:call, :f)
    for m in 1:Base.length(x.parameters)
        push!(dual_call.args, n == m ? :(Dual(x[$m], 1)) : :(x[$m]))
    end
    return :(dualpart($dual_call))
end
@generated fmad(f, x, n) = dual_call_expr(f, x, n)
function fmad_expr(f, x::Type{<:Tuple})
    body = Expr(:tuple)
    for n in 1:Base.length(x.parameters)
        push!(body.args, dual_call_expr(f, x, Type{Val{n}}))
    end
    return body
end
@generated fmad(f, x) = fmad_expr(f, x)


text:
module PetriNetUtilities

using LightGraphs, GraphPlot, DiffEqBiological

mutable struct PetriNet{GRAPH}
    g::GRAPH
    nodelabel::Vector{Symbol}
    edgelabel::Vector{String}
end

function make_petri_net(rs)
    # extract nodes
    species_nodes  = species(rs)
    reaction_nodes = [Symbol("R$j") for j in eachindex(rs.reactions)]

    number_species = length(species_nodes)

    # initialize edge set
    edge_set = Tuple{Int, Int}[]

    # labels
    nodelabel = vcat(species_nodes, reaction_nodes)
    edgelabel = String[]

    dep_graph = rxtospecies_depgraph(rs)

    for j in eachindex(dep_graph)
        # reaction -> species
        for (k, stoich) in productstoich(rs, j)
            push!(edge_set, (number_species + j, k))
            if stoich > 1
                push!(edgelabel, string(stoich))
            else
                push!(edgelabel, "")
            end
        end

        # species -> reaction
        for (k, stoich) in substratestoich(rs, j)
            push!(edge_set, (k, number_species + j))
            if stoich > 1
                push!(edgelabel, string(stoich))
            else
                push!(edgelabel, "")
            end
        end
    end

    number_nodes = length(species_nodes) + length(reaction_nodes)

    g = SimpleDiGraph(number_nodes)

    for edge in edge_set
        add_edge!(g, edge)
    end

    # return gplot(g,
    #     nodelabel = nodelabel,
    #     edgelabel = edgelabel,
    #     layout    = spectral_layout,
    #     edgelabeldistx=1.0)
    return PetriNet(g, nodelabel, edgelabel)
end

import GraphPlot: gplot

function gplot(petri_net::PetriNet, locs_x_in::Vector{R}, locs_y_in::Vector{R}; keyargs...) where {R <: Real}
    gplot(petri_net.g,
        locs_x_in, locs_y_in;
        nodelabel = petri_net.nodelabel,
        edgelabel = petri_net.edgelabel,
        keyargs...)
end

function gplot(petri_net::PetriNet; layout::Function=spring_layout, keyargs...)
    gplot(petri_net, layout(petri_net.g)...; keyargs...)
end

export make_petri_net, gplot

end # module


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module DerivTreeMDPs

using ExprSearch, LinearDerivTrees
using Reexport
using RLESUtils, Observers
import RLESTypes.SymbolTable
using POMDPs
using Iterators

export DerivTreeMDPParams, DerivTreeState, DerivTreeAction, DerivTreeMDP,
    DerivTreeStateSpace, DerivTreeActionSpace
export DerivTreeTransitionDistr, discount, n_actions, actions, iterator, reward,
    sync!, step!, get_fitness, get_expr
export create_state, create_action, create_transition_distribution

import LinearDerivTrees: step!, get_expr
import Base: ==, hash, rand!, copy!, push!
import ExprSearch: ExprProblem, get_fitness

type DerivTreeMDPParams
  grammar::Grammar
  max_neg_reward::Float64
  step_reward::Float64
  discount::Float64
  observer::Observer
end

function DerivTreeMDPParams(grammar::Grammar, max_neg_reward::Float64, step_reward::Float64,
                            discount::Float64=1.0; observer::Observer=Observer())
  return DerivTreeMDPParams(grammar, max_neg_reward, step_reward, discount, observer)
end

type DerivTreeAction <: Action
  action_id::Int64
end
DerivTreeAction() = DerivTreeAction(-1)

type DerivTreeMDP <: POMDP
  params::DerivTreeMDPParams
  tree::LinearDerivTree #true state of the sim
  problem::ExprProblem #used in calling get_fitness
  statehash::UInt64 #hash for current state for sync'ing purposes
  all_actions::Vector{DerivTreeAction}
  userargs::SymbolTable

  function DerivTreeMDP(p::DerivTreeMDPParams, tree::LinearDerivTree, problem::ExprProblem,
      all_actions::Vector{DerivTreeAction}, userargs::SymbolTable=SymbolTable())
    return new(p, tree, zero(UInt64), problem, all_actions, userargs)
  end

  function DerivTreeMDP(p::DerivTreeMDPParams, tree::LinearDerivTree, problem::ExprProblem,
      userargs::SymbolTable=SymbolTable())
    mdp = new()
    mdp.params = p
    mdp.tree = tree
    mdp.problem = problem
    mdp.statehash = zero(UInt64)
    mdp.all_actions = generate_all_actions(p.grammar)
    mdp.userargs = userargs
    return mdp
  end
end

#TODO: put mdp or tree inside state so that sync! only depends on s, i.e., s fully contains state
type DerivTreeState <: State
  mdp::DerivTreeMDP #sneak this in to have fully contained state
  past_actions::Vector{Int64} #actions taken since initialize!
end
DerivTreeState(mdp::DerivTreeMDP) = DerivTreeState(mdp, Int64[])

type DerivTreeStateSpace <: AbstractSpace
  #mdp::DerivTreeMDP
  action_id::Int64
end

type DerivTreeActionSpace <: AbstractSpace
  mdp::DerivTreeMDP
  action_ids::Vector{Int64}
end

type DerivTreeTransitionDistr <: AbstractDistribution
  nextstate::DerivTreeState
  DerivTreeTransitionDistr(mdp) = new(DerivTreeState(mdp))
  DerivTreeTransitionDistr(state::DerivTreeState) = new(state)
end

function generate_all_actions(grammar::Grammar)
  return [DerivTreeAction(i) for i = 1:maxlength(grammar)]
end

# returns the discount factor
POMDPs.discount(mdp::DerivTreeMDP) = mdp.params.discount

# returns the number of actions in the problem
function POMDPs.n_actions(mdp::DerivTreeMDP, s::DerivTreeState)
  p = mdp.params
  @notify_observer(p.observer, "debug2", ["n_actions called"])
  sync!(s)
  return length(actionspace(s.mdp.tree))
end

# returns the action space for the entire problem
POMDPs.actions(mdp::DerivTreeMDP) = DerivTreeActionSpace(mdp, collect(1:length(mdp.all_actions)))

# fills the action space with the actions availiable from state s
function POMDPs.actions(mdp::DerivTreeMDP, s::DerivTreeState, as::DerivTreeActionSpace)
  p = mdp.params
  @notify_observer(p.observer, "debug2", ["actions called"])
  sync!(s)
  as.action_ids = collect(actionspace(s.mdp.tree)) #reachable actions
  return as
end

function POMDPs.iterator(action_space::DerivTreeActionSpace)
  p = action_space.mdp.params
  @notify_observer(p.observer, "debug2", ["iterator called"])
  all_actions = action_space.mdp.all_actions
  ids = action_space.action_ids
  return imap(id -> all_actions[id], ids) #iterator avoids allocation
end

# fills d with neighboring states reachable from the s,a pair
function POMDPs.transition(mdp::DerivTreeMDP, s::DerivTreeState, a::DerivTreeAction, d::DerivTreeTransitionDistr)
  p = s.mdp.params
  @notify_observer(p.observer, "debug2", ["transition called"])
  step!(s, d.nextstate, a)
  return d
end

# fills sp with random sample from distribution d
function POMDPs.rand!(rng::AbstractRNG, sp::DerivTreeState, d::DerivTreeTransitionDistr)
  return copy!(sp, d.nextstate)
end

function step!(s::DerivTreeState, sp::DerivTreeState, a::DerivTreeAction)
  p = s.mdp.params
  @notify_observer(p.observer, "debug1", ["mdp step! called"])
  sync!(s)
  step!(s.mdp.tree, a.action_id) #deterministic transition
  copy!(sp, s)
  push!(sp, a)
  s.mdp.statehash = hash(sp) #update sync obj
  return sp
end

get_fitness() = error("fitness not defined")

# returns the immediate reward of being in state s and performing action a
POMDPs.reward(mdp::DerivTreeMDP, s::DerivTreeState, a::DerivTreeAction) = POMDPs.reward(mdp, s)

function POMDPs.reward(mdp::DerivTreeMDP, s::DerivTreeState)
  p = mdp.params
  @notify_observer(p.observer, "debug2", ["reward called"])
  sync!(s)

  tree = mdp.tree
  if iscomplete(tree)
    reward = -get_fitness(mdp.problem, get_derivtree(tree), mdp.userargs)
  elseif isdone(tree) #not-compilable
    reward = p.max_neg_reward
  else #each step
    reward = p.step_reward
  end
  return reward
end

# returns a boolean indicating if state s is terminal
function POMDPs.isterminal(mdp::DerivTreeMDP, s::DerivTreeState)
  p = mdp.params
  @notify_observer(p.observer, "debug2", ["isterminal called"])
  sync!(s)
  return isdone(s.mdp.tree)
end

# initializes a model state
POMDPs.create_state(mdp::DerivTreeMDP) = DerivTreeState(mdp)

# initializes a model action
POMDPs.create_action(mdp::DerivTreeMDP) = DerivTreeAction()

# initializes a distirbution over states
POMDPs.create_transition_distribution(mdp::DerivTreeMDP) = DerivTreeTransitionDistr(mdp)

# if you want to use a random policy you need to implement an action space sampling function
# fills action a with a random action form action_space
function POMDPs.rand!(rng::AbstractRNG, a::DerivTreeAction, action_space::DerivTreeActionSpace)
  a.action_id = rand(rng, action_space.action_ids)
  return a
end

#equate all fields
function ==(s1::DerivTreeState, s2::DerivTreeState)
  s1.past_actions == s2.past_actions
end

#hash all fields
function hash(s::DerivTreeState)
  h = hash(DerivTreeState)
  h = hash(s.mdp, h)
  for a in s.past_actions
    h = hash(a, h)
  end
  return h
end

function sync!(s::DerivTreeState)
  p = s.mdp.params
  @notify_observer(p.observer, "debug2", ["sync! called current=$(s.mdp.statehash), s=$(s.past_actions)"])
  #if sync'd, we're done
  s.mdp.statehash == hash(s) && return
  @notify_observer(p.observer, "debug1", ["sync'ing"])
  #resync
  initialize!(s.mdp.tree)
  for a in s.past_actions
    step!(s.mdp.tree, a)
  end
  s.mdp.statehash = hash(s)
end

reset!(s::DerivTreeState) = empty!(s.past_actions)

function copy!(dst::DerivTreeState, src::DerivTreeState)
  dst.mdp = src.mdp
  resize!(dst.past_actions, length(src.past_actions))
  copy!(dst.past_actions, src.past_actions)
  return dst
end

push!(s::DerivTreeState, a::DerivTreeAction) = push!(s.past_actions, a.action_id)

function get_expr(state::DerivTreeState)
  sync!(state) #sync to state
  return get_expr(state.mdp.tree)
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# uses the defaults for each dist.
# otherwise, you can manually construct expectations and glob them together
function expectation(m::UnivariateMixture; kwargs...)
    expectations = [expectation(d; kwargs...) for d in components(m)]
    return MixtureExpectation(expectations, probs(m))
end

function (e::MixtureExpectation)(f::Function; kwargs...)
    return dot(e.mixtureweights, [E(f; kwargs...) for E in e.expectations])
end

weights(e::MixtureExpectation) = e.mixtureweights
expectations(e::MixtureExpectation) = e.expectations
expectation(f::Function, m::UnivariateMixture; kwargs...) = dot(probs(m), [expectation(f, dist; kwargs...) for dist in components(m)])
*(r::Real, e::MixtureExpectation) = MixtureExpectation(r*expectations(e), weights(e))
# *(e::MixtureExpectation, h::AbstractArray) = dot(weights(e), [E*h for E in expectations(e)])
import Base.+
+(expectations::Expectation...) = MixtureExpectation(expectations, ones(length(expectations)))


text:
"""
    Logfile(id, [path = pwd()])

Create a new logfile with the given network `id` at `path`.

Note that the `id` does not have to be equal (or similar for
that matter) to the file name in `path`.
"""
struct Logfile
    id::String
    path::String

    function Logfile(id::String, path::String = pwd(); append = false)
        path = endswith(path, "/") || isdir(path) ? joinpath(path, _make_filename(path, id, append)) : path
        path = endswith(path, ".h5") ? path : path * ".h5"
        mkpath(dirname(path))
        if !append || !isfile(path)
            h5open(path, "w") do root
                # attrs(root)[ATTR_NETWORKNAME] = id
            end
        end
        new(id, path)
    end
end

function _make_filename(dir, id, append = false)
    basename = replace(id, r"\W", "")
    fname = basename
    i = 1
    append && return fname
    while isfile(joinpath(dir, fname * ".h5"))
        fname = basename * "_" * string(i)
        i += 1
    end
    fname
end


text:
using NCDatasets
using ArgParse
using JSON
using Formatting

include("./interpolate.jl")

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--input-file"
            help = "Input file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--input-zdomain-file"
            help = "Resolution keywords. Currently accept: POP2, Standard."
            arg_type = String
            required = true

        "--input-zdomain-varname"
            help = "Resolution keywords. Currently accept: POP2, Standard."
            arg_type = String
            required = true

        "--output-zdomain-file"
            help = "Resolution keywords. Currently accept: POP2, Standard."
            arg_type = String
            required = true

        "--output-zdomain-varname"
            help = "Resolution keywords. Currently accept: POP2, Standard."
            arg_type = String
            required = true

        "--varname"
            help = "Variable name"
            arg_type = String
            required = true

        "--input-zdomain-type"
            help = "Currently accept: endpoints, midpoints"
            arg_type = String
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["input-zdomain-file"], "r") do ds
    global zs_i = ds[parsed["input-zdomain-varname"]][:] |> nomissing
end

Dataset(parsed["output-zdomain-file"], "r") do ds
    global zs_o = ds[parsed["output-zdomain-varname"]][:] |> nomissing
end

if parsed["input-zdomain-type"] == "endpoints"
    println("Input zdomain use endpoints. Take average.")
    zs_i_mid = (zs_i[2:end] + zs_i[1:end-1] ) / 2.0
elseif parsed["input-zdomain-type"] == "midpoints"
    println("Input zdomain use midpoints. Do nothing.")
    zs_i_mid = copy(zs_i)
    # do nothing
else
    throw(ErrorException("Unknown zdomain type: " * string(parsed["input-zdomain-type"])))
end

zs_o_mid = (zs_o[2:end] + zs_o[1:end-1] ) / 2.0

missing_value = 1e20

ds_i    = Dataset(parsed["input-file"], "r")
ds_o    = Dataset(parsed["output-file"], "c")

Nx = ds_i.dim["Nx"]
Ny = ds_i.dim["Ny"]
Nz_o = length(zs_o) - 1

defDim(ds_o, "Nx", Nx)
defDim(ds_o, "Ny", Ny)
defDim(ds_o, "Nz", Nz_o)
defDim(ds_o, "zs", length(zs_o))
defDim(ds_o, "time", Inf)

for (varname, vardata, dims) in (
    ("zs", zs_o, ("zs",)),
)
    println("varname: ", varname)
    v = defVar(ds_o, varname, Float64, dims)
    v[:] = vardata
end

old_data = replace(ds_i[parsed["varname"]][:], missing=>NaN)
new_data = zeros(Float64, Nx, Ny, Nz_o)

println("size of old_data: ", size(old_data))
println("size of new_data: ", size(new_data))

new_data .= NaN

# interpolation function needs
# monotonic increasing function
depth_i = - zs_i_mid
depth_o = - zs_o_mid

for i=1:Nx, j=1:Ny

    # detect rng
    valid_data_cnt = sum(isfinite.(old_data[i, j, :]))

    if valid_data_cnt == 0

        continue

    elseif valid_data_cnt == 1  # no interpolation at all

        println(format("Special case happens at (i, j) = ({},{})", i, j))
        new_data[i, j, :] .= old_data[i, j, 1]

    else

        new_data[i, j, :] = interpolate(
            depth_i[1:valid_data_cnt],
            old_data[i, j, 1:valid_data_cnt],
            depth_o;
            left_copy=true,
            right_copy=true,
        )

    end

end

new_var = defVar(ds_o, parsed["varname"], Float64, ("Nx", "Ny", "Nz", "time"))
new_var.attrib["_FillValue"] = missing_value
new_var[:, :, :, 1] = new_data

close(ds_i)
close(ds_o)

println("Total data cnt: ", length(new_data))
println("Valid data cnt: ", isfinite.(new_data) |> sum)
println("Data holes cnt: ", isnan.(new_data) |> sum)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module AoC2021

using ReTest
using DocStringExtensions

TEST_DATA_DIRECTORY = joinpath(dirname(@__DIR__), "test_data")

"""Information associated with a particular puzzle."""
abstract type PuzzleInfo end

"""
    $(FUNCTIONNAME)(puzzleinfo)

Compute the solution to an Advent of Code puzzle described by `puzzleinfo`.
"""
function solve(::PuzzleInfo) end

for srcfile in readdir(@__DIR__)
    if occursin(r"^Day[0-9]{2}\.jl", srcfile)
        include(srcfile)
    end
end

end # module


text:
function coarsen_grid(grd::Grid)

    _, ks, _ = findnz(grd.Proj)

    _Nxc = Int64(ceil(grd._Nx/2))
    _Nyc = Int64(ceil(grd._Ny/2))
    Nbuffer = Int64(ceil(grd._Nbuffer/2))

    dr = 2*grd.dr
    if iseven(grd._Nbuffer)
        r0 = grd.r0
        nan_outside_boundaries = grd._nan_outside_boundaries[1:2:end,1:2:end]
    else
        r0 = grd.r0 + [1,1]*grd.dr
        nan_outside_boundaries = grd._nan_outside_boundaries[2:2:end,2:2:end]
    end
    Nxc, Nyc = size(nan_outside_boundaries)
    r1 = r0 + dr*[Nxc, Nyc]

    points = zeros(2, grd.Nk)
    projc_js = zeros(grd.Nk)
    Nkc = 0
    for k=1:grd.Nk
        kc = ks[k]
        j = Int(ceil(kc / grd._Nx))
        i = kc - (j - 1) * grd._Nx
        if isodd(i) && isodd(j)
            Nkc += 1
            points[:,Nkc] = grd.points[:,k]

            ic = (i + 1)/2
            jc = (j + 1)/2
            projc_js[Nkc] = ic + (jc-1)*_Nxc
        end
    end

    Projc = sparse([1:Nkc...], projc_js[1:Nkc], ones(Nkc), Nkc, _Nxc*_Nyc)
    return Grid(r0, r1, points[:,1:Nkc], Projc, dr, Nkc, _Nxc, _Nyc, Nbuffer, nan_outside_boundaries)
end

function restriction_and_interpolation(xb, finegrid, coarsegrid, extr_c, extr_f)

    Interp_cart, Restr_cart = cartesian_2d_crs(coarsegrid._Nx, coarsegrid._Ny, 4, 4)
    Interp_i = finegrid.Proj*Interp_cart*transpose(coarsegrid.Proj)
    Restr_i = coarsegrid.Proj*Restr_cart*transpose(finegrid.Proj)
    xb_c = Restr_i*xb

    Interp_mat = extr_f.Proj*Interp_i*extr_c.R*transpose(extr_c.Proj)
    Interp_f = extr_f.Proj*Interp_i*(I-extr_c.R)*xb_c

    Restr_mat = extr_c.Proj*Restr_i*extr_f.R*transpose(extr_f.Proj)
    Restr_f = extr_c.Proj*Restr_i*(I-extr_f.R)*xb

    Interp = AffineMap{Float64}(Interp_mat, Interp_f)
    Restr = AffineMap{Float64}(Restr_mat, Restr_f)
    return Interp, Restr, xb_c
end


text:
using ContinuousDPs
using QuantEcon
using BasisMatrices
using Test

include("test_cdp.jl")


text:
# https://docs.julialang.org/en/v1/manual/modules/

using Test

module Bar
    x = 13.0

    function foo()
        return x
    end
end

#import .Bar # . relative import path

@test 13 == Bar.foo()

module A
    x = 17.0
    function foo()
        return x
    end

    module AA
        x = 19.0
        function foo()
            return x
        end
    end
    module AB
        x = 23.0
        function foo()
            return x
        end
    end

    module AAB
        import ..AA
        import ..AB

        function foo()
            return AA.foo() + AB.foo()
        end
    end
end

#import .A

@test 17 == A.foo()

@test 19 == A.AA.foo()
@test 23 == A.AB.foo()

#import .A.AAB
@test 42 == A.AAB.foo()


text:
using LinearResponseVariationalBayes
using HessianReparameterization
using Base.Test

println("Testing HessianReparameterization.")

function test_lgamma_hessian_reparameterization()
	# An example function and transform to check the lgamma and digamma functions.

	K = 2
	a = 1.5

	function f{T <: Number}(y::Array{T})
		y[1]^3 + y[2]^2 + a * y[1] * y[2]
	end

	function y_to_x{T <: Number}(y::Array{T})
		[ lgamma(0.1 * y[1] + 0.9 * y[2]),
		  digamma(0.9 * y[1] + 0.1 * y[2]) ]
	end

	dx_dy_func = get_dx_dy_func(y_to_x, K)
	d2x_dy2_funcs = get_d2x_dy2_funcs(y_to_x, K)
	df_dy_func = get_df_dy_func(f, K)
	d2f_dy2_func = get_d2f_dy2_func(f, K)

	y = Float64[1.2, 2.1]

	df_dy = df_dy_func(y)
	d2f_dy2 = d2f_dy2_func(y)
	dx_dy = dx_dy_func(y)
	d2x_dy2 = get_d2x_dy2(d2x_dy2_funcs, y)

	d2f_dx2 = transform_hessian(dx_dy, d2x_dy2, df_dy, d2f_dy2)

	# Test with numeric derivatives.
	function df_dx_func{T <: Number}(y::Array{T})
		dx_dy_func(y)' \ df_dy_func(y)
	end

	eps = 1e-6
	y_a = y + eps * Float64[1, 0]
	y_b = y + eps * Float64[0, 1]

	x = y_to_x(y)
	x_a = y_to_x(y_a) - x
	x_b = y_to_x(y_b) - x

	f_a = f(y_a) - f(y)
	f_b = f(y_b) - f(y)

	# Test the first derivative.
	df_dx_est = vcat(x_a', x_b') \ Float64[f_a, f_b]
	df_dx_func(y)

	@test_approx_eq_eps df_dx_est df_dx_func(y) 1e-5

	# Test the second derivative.
	g_a = df_dx_func(y_a) - df_dx_func(y)
	g_b = df_dx_func(y_b) - df_dx_func(y)

	d2f_dx2_est =  vcat(x_a', x_b') \ vcat(g_a', g_b')
	@test_approx_eq_eps d2f_dx2_est d2f_dx2 1e-4
end

function test_hessian_reparameterization()
	# An example function and transform.
	K = 2
	a = 1.5

	function f{T <: Number}(y::Array{T})
		y[1]^3 + y[2]^2 + a * y[1] * y[2]
	end

	function y_to_x{T <: Number}(y::Array{T})
		T[ y[1] * (y[2] ^ 2),
		   y[1] + y[2]]
	end

	dx_dy_func = get_dx_dy_func(y_to_x, K)
	d2x_dy2_funcs = get_d2x_dy2_funcs(y_to_x, K)
	df_dy_func = get_df_dy_func(f, K)
	d2f_dy2_func = get_d2f_dy2_func(f, K)

	y = Float64[1.7, 3]

	df_dy = df_dy_func(y)
	d2f_dy2 = d2f_dy2_func(y)
	dx_dy = dx_dy_func(y) # Actually the transpose of the Jacobian
	d2x_dy2 = get_d2x_dy2(d2x_dy2_funcs, y)

	d2f_dx2 = transform_hessian(dx_dy, d2x_dy2, df_dy, d2f_dy2)

	# Test with numeric derivatives.
	function df_dx_func{T <: Number}(y::Array{T})
		dx_dy_func(y) \ df_dy_func(y)
	end

	function df_dx_func{T <: Number}(y::Array{T})
		dx_dy_func(y)' \ df_dy_func(y)
	end

	eps = 1e-8
	y_a = y + eps * Float64[1, 0]
	y_b = y + eps * Float64[0, 1]

	x = y_to_x(y)
	x_a = y_to_x(y_a) - x
	x_b = y_to_x(y_b) - x

	f_a = f(y_a) - f(y)
	f_b = f(y_b) - f(y)

	# Test the first derivative.
	df_dx_est = vcat(x_a', x_b') \ Float64[f_a, f_b]
	df_dx_func(y)

	@test_approx_eq_eps df_dx_est df_dx_func(y) 1e-5

	# Test the second derivative.
	g_a = df_dx_func(y_a) - df_dx_func(y)
	g_b = df_dx_func(y_b) - df_dx_func(y)

	d2f_dx2_est =  vcat(x_a', x_b') \ vcat(g_a', g_b')
	@test_approx_eq_eps d2f_dx2_est d2f_dx2 1e-4
end

test_hessian_reparameterization()

# This will be failing until ForwardDiff.jl releases v0.0.3 to METADATA.
test_lgamma_hessian_reparameterization()


text:
"""
    generate_constraints(
        ; mono::Vector{Vector{S}}=Vector{String}[],
        not_mono::Vector{Vector{S}}=Vector{String}[],
        exc::Vector{Tuple{Vector{S}, Vector{S}}}=Vector{Tuple{Vector{String}, Vector{String}}}[]
    )::ConstraintDict where S<:AbstractString

Generate a dictionary of constraints based on the given arguments.
Mono for all monophyletic groups. not_mono for leafs that are not allowed
to form a monophletic group. And exc for a partial constraint, where one or more
leafs are not allowed to be part of a specific clade.

e.g.: generate_constraints(mono=[["a", "b"], ["a", "c"]],
                           exc=[(["b", "c"],["d"]])
"""
function generate_constraints(
    ; mono::Vector{Vector{S}}=Vector{String}[],
    not_mono::Vector{Vector{S}}=Vector{String}[],
    exc::Vector{Tuple{Vector{S}, Vector{S}}}=Tuple{Vector{String}, Vector{String}}[]
)::ConstraintDict where S<:AbstractString

    constraints_dict = ConstraintDict()
    lengths = [length(mono), length(not_mono), length(exc)]
    # filter out trivial constraints...
    filter!(x -> length(x) > 1, mono)
    filter!(x -> length(x) > 1, not_mono)
    filter!(x -> length(x[1]) > 1, exc)
    filter!(x -> length(x[2]) > 0, exc)
    # ... and inform the user about it
    if (length(mono) < lengths[1]) || (length(not_mono) < lengths[2])
        @warn "Some trivial 'mono' / 'not_mono' type constraints were removed.
         A valid 'mono' / 'not_mono' constraint needs at least 2 elements."
    end # if
    if length(exc) < lengths[3]
        @warn "Some trivial 'exc' type constraints were removed.
         A non-trivial 'exc' constraints needs at least 2 elements in the first, and at least 1 in the second part of the tuple"
    end # end
    constraints_dict[:mono] = unique(mono)
    constraints_dict[:not_mono] = unique(not_mono)
    constraints_dict[:exc] = unique(exc)
    return constraints_dict
end # generate_constraints

"""
    generate_constraints!(
        constraints::Dict;
        mono::Vector{Vector{S}}=Vector{String}[],
        not_mono::Vector{Vector{S}}=Vector{String}[],
        exc::Vector{Tuple{Vector{S}, Vector{S}}}=Vector{Tuple{Vector{String}, Vector{String}}}[]
    )::ConstraintDict where S<:AbstractString

Function that adds further constraints to an existing dictionary of constraints.
See basic generate_constraints function for more info.
"""
function generate_constraints!(constraints::ConstraintDict;
                               mono::Vector{Vector{S}}=Vector{String}[],
                               not_mono::Vector{Vector{S}}=Vector{String}[],
                               exc::Vector{Tuple{Vector{S}, Vector{S}}}=Tuple{Vector{String}, Vector{String}}[]
                               )::ConstraintDict where S<:AbstractString

    constraints2 = generate_constraints(mono=mono, not_mono=not_mono, exc=exc)
    # make sure entries are unique for each constraint category
    for key in keys(constraints)
        constraints[key] = union(constraints[key], constraints2[key])
    end # for
    return constraints
end # generate_constraints

"""
    generate_constraints(filename::String)

Function that creates a dictionary of constraints, based on a txt file with a
specific format, i.e.:

# lines starting with '#' will be ignored
# multiple constraints per line allowed
mono:        A,B; C,D;
mono:        E,F;
# final semicolon is optional
not_mono:    C,D,E
exc :        E,F!G;

See basic generate_constraints function for more info.
"""
function generate_constraints(filename::String)
    mono, not_mono, exc = parse_constraints(filename)
    generate_constraints(mono=mono, not_mono=not_mono, exc=exc)
end # generate_constraints

"""
    generate_constraints!(constraints::ConstraintDict, filename::String)

Function that adds further constraints to an existing dictionary of constraints,
based on a txt file with a specific format (each line one constraint). See basic
generate_constraints function on a file for more info.
"""
function generate_constraints!(constraints::ConstraintDict, filename::String)
    mono, not_mono, exc = parse_constraints(filename)
    generate_constraints!(constraints; mono=mono, not_mono=not_mono, exc=exc)
end # generate_constraints!

"""
    parse_constraints(filename::S
                     )::Tuple{
                        Vector{Vector{S}},
                        Vector{Vector{S}},
                        Vector{Tuple{Vector{S}, Vector{S}}}
                     } where S<:AbstractString

--- INTERNAL ---
Helper function that reads constraints out of a file
"""
function parse_constraints(filename::S)::Tuple{Vector{Vector{S}}, Vector{Vector{S}}, Vector{Tuple{Vector{S}, Vector{S}}}} where S<:AbstractString

    mono, not_mono = [Vector{String}[] for _ = 1:2]
    exc = Vector{Tuple{Vector{AbstractString}, Vector{AbstractString}}}()
    for line in eachline(filename)
        # remove whitespace
        line = filter(x -> !isspace(x), line)
        # skip comment lines
        startswith(line, "#") && continue
        split_l::Vector{AbstractString} = split(line,":")
        length(split_l) != 2 && throw(FileSyntaxError("There should be exactly one colon in each non-comment line."))
        constraints::Vector{AbstractString} = split(split_l[2], ";")
        if split_l[1] == "mono"
            for constraint in constraints
                # makes sure that parsing works with & without trailing semicolon
                constraint == "" && continue
                push!(mono, split(constraint, ","))
            end # for
        elseif split_l[1] == "not_mono"
            for constraint in constraints
                constraint == "" && continue
                push!(not_mono, split(constraint, ","))
            end # for
        elseif split_l[1] == "exc"
            for constraint in constraints
                constraint == "" && continue
                exc_constraint::Vector{AbstractString} = split(constraint, "!")
                length(exc_constraint) != 2 && throw(FileSyntaxError("There should be exactly one exclamation mark in each 'exc' type constraint."))
                push!(exc, (split(exc_constraint[1], ","),
                            split(exc_constraint[2], ",")))
            end # for
        else
            @warn "Skipped line with unsupported constraint type '$split_l[1]'.
             Allowed types are 'mono', 'not_mono' and 'exc'"
        end # if/else
    end # for
    return (mono, not_mono, exc)
end # parse_constraints

function topological(constraints::ConstraintDict, tree::GeneralNode)::Bool
    isempty(constraints) && return true
    for leaves in constraints[:mono]
        lca = find_lca(tree, leaves)
        lca.root && return false
        for child in get_leaves(lca)
            !(child.name in leaves) && return false
        end # for
    end # for
    for leaves in constraints[:not_mono]
        lca = find_lca(tree, leaves)
        length(get_leaves(lca)) == length(leaves) && return false
    end # for
    for leaves in constraints[:exc]
        lca = find_lca(tree, leaves[1])
        lca.root && return false
        for leaf in leaves[2]
            leaf = find_name(tree, leaf)
            (lca.binary == leaf.binary[1:length(lca.binary)]) && return false
        end # for
    end # for
    true
end # topological

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
@testset "Example 2" begin
    Mod.@polyvar x[1:3]
    p = sum(x .* x)
    @test p == x[1]^2 + x[2]^2 + x[3]^2
    @test subs(p, x[1]=>2, x[3]=>3) == x[2]^2 + 13
    A = [1.0 2 3; 4 5 6; 7 8 9]
    @test p(x=>A*vec(x)) == (x[1] + 2x[2] + 3x[3])^2 + (4x[1] + 5x[2] + 6x[3])^2 + (7x[1] + 8x[2] + 9x[3])^2
end