text:
# need special handling for showing a string as a textmime
# type, since in that case the string is assumed to be
# raw data unless it is text/plain
israwtext(::MIME, x::AbstractString) = true
israwtext(::MIME"text/plain", x::AbstractString) = false
israwtext(::MIME, x) = false

InlineIOContext(io, KVs::Pair...) = IOContext(
    io,
    :limit=>true, :color=>true, :jupyter=>true,
    KVs...
)

# convert x to a string of type mime, making sure to use an
# IOContext that tells the underlying show function to limit output
function limitstringmime(mime::MIME, x)
    buf = IOBuffer()
    if istextmime(mime)
        if israwtext(mime, x)
            return String(x)
        else
            show(InlineIOContext(buf), mime, x)
        end
    else
        b64 = Base64.Base64EncodePipe(buf)
        if isa(x, Vector{UInt8})
            write(b64, x) # x assumed to be raw binary data
        else
            show(InlineIOContext(b64), mime, x)
        end
        close(b64)
    end
    return String(take!(buf))
end


text:
module TestFoldl

using Test
using LazyGroupBy
using LazyGroupBy: foldxl

@testset "tuple" begin
    @testset for fold in [foldl, foldxl]
        foldl = nothing

        @test fold.(tuple, grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) ==
              Dict(false => ((0, 4), 0), true => ((((((7, 3), 1), 5), 9), 3), 5))

        @test fold.(tuple, grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]); init = -1) ==
              Dict(
            false => (((-1, 0), 4), 0),
            true => (((((((-1, 7), 3), 1), 5), 9), 3), 5),
        )
    end
end

extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))

@testset "extrema" begin
    @test mapfoldl.(
        x -> (x, x),
        extrema_rf,
        grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]),
    ) == Dict(false => (0, 4), true => (1, 9))
end

end  # module

ext: jl


text:
# ---
# title: 1558. Minimum Numbers of Function Calls to Make Target Array
# id: problem1558
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-numbers-of-function-calls-to-make-target-array/description/>
# hidden: true
# ---
#
# ![](https://assets.leetcode.com/uploads/2020/07/10/sample_2_1887.png)
#
# Your task is to form an integer array `nums` from an initial array of zeros
# `arr` that is the same size as `nums`.
#
# Return the minimum number of function calls to make `nums` from `arr`.
#
# The answer is guaranteed to fit in a 32-bit signed integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,5]
#     Output: 5
#     Explanation: Increment by 1 (second element): [0, 0] to get [0, 1] (1 operation).
#     Double all the elements: [0, 1] -> [0, 2] -> [0, 4] (2 operations).
#     Increment by 1 (both elements)  [0, 4] -> [1, 4] -> **[1, 5]** (2 operations).
#     Total of operations: 1 + 2 + 2 = 5.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,2]
#     Output: 3
#     Explanation: Increment by 1 (both elements) [0, 0] -> [0, 1] -> [1, 1] (2 operations).
#     Double all the elements: [1, 1] -> **[2, 2]** (1 operation).
#     Total of operations: 2 + 1 = 3.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [4,2,5]
#     Output: 6
#     Explanation: (initial)[0,0,0] -> [1,0,0] -> [1,0,1] -> [2,0,2] -> [2,1,2] -> [4,2,4] -> **[4,2,5]** (nums).
#
#
# **Example 4:**
#
#
#
#     Input: nums = [3,2,2,4]
#     Output: 7
#
#
# **Example 5:**
#
#
#
#     Input: nums = [2,4,8,16]
#     Output: 8
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `0 <= nums[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
################################################
### Boundary layer conductance formulations ####
################################################

#' Boundary Layer Conductance according to Thom 1972
#'
#' An empirical formulation for the canopy boundary layer conductance
#'              for heat transfer based on a simple ustar dependency.
#'
#' - ustar     Friction velocity (m s-1)
#' - Sc        Optional: Schmidt number of additional quantities to be calculated
#' - Sc_name   Optional: Name of the additional quantities, has to be of same length than
#'                  `Sc_name`
#' - constants k - von-Karman constant
#'                  Sc_CO2 - Schmidt number for CO2
#'                  Pr - Prandtl number (if `Sc` is provided)
#'
#'
#' # Details
 The empirical equation for Rb suggested by Thom 1972 is:
#'
#'      ``Rb = 6.2ustar^-0.67``
#'
#'    Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
#'    Gb for other quantities x is calculated as (Hicks et al. 1987):
#'
#'      ``Gb_x = Gb / (Sc_x / Pr)^0.67``
#'
#'  where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
#'
#' # Value
 a DataFrame with the following columns:
#'  - Gb_h: Boundary layer conductance for heat transfer (m s-1)
#'  - Rb_h: Boundary layer resistance for heat transfer (s m-1)
#'  - kB_h: kB-1 parameter for heat transfer
#'  \item{Gb_Sc_name}{Boundary layer conductance for `Sc_name` (m s-1). Only added if `Sc_name` and
#'                    `Sc_name` are provided}
#'
#' #References
#' Thom, A., 1972: Momentum, mass and heat exchange of vegetation.
#'             Quarterly Journal of the Royal Meteorological Society 98, 124-134.
#'
#'             Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
#'             A preliminary multiple resistance routine for deriving dry deposition velocities
#'             from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
#'
#' #See also
#' [`Gb_Choudhury`](@ref), [`Gb_Su`](@ref), [`aerodynamic_conductance!`](@ref)
#'
#' ```@example; output = false
#' ```
#' Gb_Thom(seq(0.1,1.4,0.1))
#'
#' ## calculate Gb for SO2 as well
#' Gb_Thom(seq(0.1,1.4,0.1),Sc=1.25,Sc_name="SO2")
#'
"""
"""
function Gb_Thom(ustar,Sc=nothing,Sc_name=nothing,constants=BigleafConstants())

  check_input(nothing,ustar)

  Rb_h = 6.2*ustar^-0.667
  Gb_h = 1/Rb_h
  kB_h = Rb_h*constants.k*ustar

  if (!isnothing(Sc) | !isnothing(Sc_name))
    if (length(Sc) != length(Sc_name))
      stop("arguments 'Sc' and 'Sc_name' must have the same length")
end
    if (!is_numeric(Sc))
      stop("argument 'Sc' must be numeric")
end
end

  Sc   = c(constants.Sc_CO2,Sc)
  Gb_x = DataFrame(lapply(Sc,function(x) Gb_h / (x/constants.Pr)^0.67))
  colnames(Gb_x) = paste0("Gb_",c("CO2",Sc_name))

  return(DataFrame(Gb_h,Rb_h,kB_h,Gb_x))
end

#' Boundary Layer Conductance according to Choudhury & Monteith 1988
#'
#' A formulation for the canopy boundary layer conductance
#'              for heat transfer according to Choudhury & Monteith 1988.
#'
#' - data             DataFrame or matrix containing all required variables
#' - Tair             Air temperature (degC)
#' - pressure         Atmospheric pressure (kPa)
#' - wind             Wind speed at sensor height (m s-1)
#' - ustar            Friction velocity (m s-1)
#' - H                Sensible heat flux (W m-2)
#' - leafwidth        Leaf width (m)
#' - LAI              One-sided leaf area index
#' - zh               Canopy height (m)
#' - zr               Instrument (reference) height (m)
#' - d                Zero-plane displacement height (-), can be calculated using `roughness_parameters`
#' - z0m              Roughness length for momentum (m). If not provided, calculated from `roughness_parameters`
#'                         within `wind_profile`
#' - stab_formulation Stability correction function used (If `stab_correction = true`).
#'                         Either `Val(:Dyer_1970)` or `Val(:Businger_1971)`.
#' - Sc               Optional: Schmidt number of additional quantities to be calculated
#' - Sc_name          Optional: Name of the additonal quantities, has to be of same length than
#'                         `Sc_name`
#' - constants        k - von-Karman constant
#'                         Sc_CO2 - Schmidt number for CO2
#'                         Pr - Prandtl number (if `Sc` is provided)
#'
#' # Value
 A data frame with the following columns:
#'  - Gb_h: Boundary layer conductance for heat transfer (m s-1)
#'  - Rb_h: Boundary layer resistance for heat transfer (s m-1)
#'  - kB_h: kB-1 parameter for heat transfer
#'  \item{Gb_Sc_name}{Boundary layer conductance for `Sc_name` (m s-1). Only added if `Sc_name` and
#'                    `Sc_name` are provided}
#'
#' # Details
 Boundary layer conductance according to Choudhury & Monteith 1988 is
#'          given by:
#'
#'            ``Gb_h = LAI((2a/\\alpha)*sqrt(u(h)/w)*(1-exp(-\\alpha/2)))``
#'
#'          where u(zh) is the wind speed at the canopy surface, approximated from
#'          measured wind speed at sensor height zr and a wind extinction coefficient ``\\alpha``:
#'
#'            ``u(zh) = u(zr) / (exp(\\alpha(zr/zh -1)))``.
#'
#'          ``\\alpha`` is modeled as an empirical relation to LAI (McNaughton & van den Hurk 1995):
#'
#'            ``\\alpha = 4.39 - 3.97*exp(-0.258*LAI)``
#'
#'          Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
#'          Gb for other quantities x is calculated as (Hicks et al. 1987):
#'
#'            ``Gb_x = Gb / (Sc_x / Pr)^0.67``
#'
#'          where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
#'
#' # Note
#' If the roughness length for momentum (`z0m`) is not provided as input, it is estimated
#'       from the function `roughness_parameters` within `wind_profile`. This function
#'       estimates a single `z0m` value for the entire time period! If a varying `z0m` value
#'       (e.g. across seasons or years) is required, `z0m` should be provided as input argument.
#'
#' #References
#' Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat
#'             budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
#'
#'             McNaughton, K. G., Van den Hurk, B_J_J_M., 1995: A 'Lagrangian' revision of
#'             the resistors in the two-layer model for calculating the energy budget of a
#'             plant canopy. Boundary-Layer Meteorology 74, 261-288.
#'
#'             Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
#'             A preliminary multiple resistance routine for deriving dry deposition velocities
#'             from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
#'
#' #See also
#' [`Gb_Thom`](@ref), [`Gb_Su`](@ref), [`aerodynamic_conductance!`](@ref)
#'
#' ```@example; output = false
#' ```
#' ## bulk canopy boundary layer resistance for a closed canopy (LAI=5)
#' ## with large leaves (leafwdith=0.1)
#' df = DataFrame(Tair=25,pressure=100,wind=c(3,4,5),ustar=c(0.5,0.6,0.65),H=c(200,230,250))
#' Gb_Choudhury(data=df,leafwidth=0.1,LAI=5,zh=25,d=17.5,zr=40)
#'
#' ## same conditions, but smaller leaves (leafwidth=0.01)
#' Gb_Choudhury(data=df,leafwidth=0.01,LAI=5,zh=25,d=17.5,zr=40)
#'
#' @export
function Gb_Choudhury(data,Tair="Tair",pressure="pressure",wind="wind",ustar="ustar",H="H",
                         leafwidth,LAI,zh,zr,d,z0m=nothing,stab_formulation=c(Val(:Dyer_1970),Val(:Businger_1971)),
                         Sc=nothing,Sc_name=nothing,constants=BigleafConstants())

  stab_formulation = match_arg(stab_formulation)

  check_input(data,list(Tair,pressure,wind,ustar,H))

  alpha   = 4.39 - 3.97*exp(-0.258*LAI)

  if (isnothing(z0m))
    estimate_z0m = true
    z0m = nothing
else
    estimate_z0m = false
end

  wind_zh = wind_profile(data=data,z=zh,Tair=Tair,pressure=pressure,ustar=ustar,H=H,
                          zr=zr,estimate_z0m=estimate_z0m,zh=zh,d=d,z0m=z0m,frac_{z0m}=nothing,
                          stab_correction=true,stab_formulation=stab_formulation)

  ## avoid zero windspeed
  wind_zh = pmax(0.01,wind_zh)

  if (!isnothing(Sc) | !isnothing(Sc_name))
    if (length(Sc) != length(Sc_name))
      stop("arguments 'Sc' and 'Sc_name' must have the same length")
end
    if (!is_numeric(Sc))
      stop("argument 'Sc' must be numeric")
end
end

  Gb_h = LAI*((0.02/alpha)*sqrt(wind_zh/leafwidth)*(1-exp(-alpha/2)))
  Rb_h = 1/Gb_h
  kB_h = Rb_h*constants.k*ustar

  Sc   = c(constants.Sc_CO2,Sc)
  Gb_x = DataFrame(lapply(Sc,function(x) Gb_h / (x/constants.Pr)^0.67))
  colnames(Gb_x) = paste0("Gb_",c("CO2",Sc_name))

  return(DataFrame(Gb_h,Rb_h,kB_h,Gb_x))
end

#' Boundary Layer Conductance according to Su et al. 2001
#'
#' A physically based formulation for the canopy boundary layer conductance
#'              to heat transfer according to Su et al. 2001.
#'
#' - data      DataFrame or matrix containing all required variables
#' - Tair      Air temperature (degC)
#' - pressure  Atmospheric pressure (kPa)
#' - ustar     Friction velocity (m s-1)
#' - wind      Wind speed (m s-1)
#' - H         Sensible heat flux (W m-2)
#' - zh        Canopy height (m)
#' - zr        Reference height (m)
#' - d         Zero-plane displacement height (-), can be calculated using `roughness_parameters`
#' - z0m       Roughness length for momentum (m). If not provided, calculated from `roughness_parameters`
#'                  within `wind_profile`
#' - Dl        Leaf characteristic dimension (m)
#' - fc        Fractional vegetation cover [0-1] (if not provided, calculated from LAI)
#' - LAI       One-sided leaf area index (-)
#' - N         Number of leaf sides participating in heat exchange (defaults to 2)
#' - Cd        Foliage drag coefficient (-)
#' - hs        Roughness height of the soil (m)
#' - stab_formulation Stability correction function used (If `stab_correction = true`).
#'                         Either `Val(:Dyer_1970)` or `Val(:Businger_1971)`.
#' - Sc        Optional: Schmidt number of additional quantities to be calculated
#' - Sc_name   Optional: Name of the additional quantities, has to be of same length than
#'                  `Sc_name`
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  pressure0 - reference atmospheric pressure at sea level (Pa)
#'                  Tair0 - reference air temperature (K)
#'                  Sc_CO2 - Schmidt number for CO2
#'                  Pr - Prandtl number (if `Sc` is provided)
#'
#' # Value
 A DataFrame with the following columns:
#'  - Gb_h: Boundary layer conductance for heat transfer (m s-1)
#'  - Rb_h: Boundary layer resistance for heat transfer (s m-1)
#'  - kB_h: kB-1 parameter for heat transfer
#'  \item{Gb_Sc_name}{Boundary layer conductance for `Sc_name` (m s-1). Only added if `Sc_name` and
#'                    `Sc_name` are provided}
#'
#' # Details
 The formulation is based on the kB-1 model developed by Massman 1999.
#'          Su et al. 2001 derived the following approximation:
#'
#'            ``kB-1 = (k Cd fc^2) / (4Ct ustar/u(zh)) + kBs-1(1 - fc)^2``
#'
#'          If fc (fractional vegetation cover) is missing, it is estimated from LAI:
#'
#'            ``fc = 1 - exp(-LAI/2)``
#'
#'          The wind speed at the top of the canopy is calculated using function
#'          [`wind_profile`](@ref).
#'
#'          Ct is the heat transfer coefficient of the leaf (Massman 1999):
#'
#'            ``Ct = Pr^-2/3 Reh^-1/2 N``
#'
#'          where Pr is the Prandtl number (set to 0.71), and Reh is the Reynolds number for leaves:
#'
#'            ``Reh = Dl wind(zh) / v``
#'
#'          kBs-1, the kB-1 value for bare soil surface, is calculated according
#'          to Su et al. 2001:
#'
#'            ``kBs^-1 = 2.46(Re)^0.25 - ln(7.4)``
#'
#'          Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.
#'          Gb for other quantities x is calculated as (Hicks et al. 1987):
#'
#'            ``Gb_x = Gb / (Sc_x / Pr)^0.67``
#'
#'          where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).
#'
#' # Note
#' If the roughness length for momentum (`z0m`) is not provided as input, it is estimated
#'       from the function `roughness_parameters` within `wind_profile`. This function
#'       estimates a single `z0m` value for the entire time period! If a varying `z0m` value
#'       (e.g. across seasons or years) is required, `z0m` should be provided as input argument.
#'
#'
#' #References
#' Su, Z., Schmugge, T., Kustas, W. & Massman, W., 2001: An evaluation of
#'             two models for estimation of the roughness height for heat transfer between
#'             the land surface and the atmosphere. Journal of Applied Meteorology 40, 1933-1951.
#'
#'             Massman, W., 1999: A model study of kB H- 1 for vegetated surfaces using
#'            'localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.
#'
#'             Hicks, B_B., Baldocchi, D_D., Meyers, T_P., Hosker, J_R., Matt, D_R., 1987:
#'             A preliminary multiple resistance routine for deriving dry deposition velocities
#'             from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
#'
#' #See also
#' [`Gb_Thom`](@ref), [`Gb_Choudhury`](@ref), [`aerodynamic_conductance!`](@ref)
#'
#' ```@example; output = false
#' ```
#' # Canopy boundary layer resistance (and kB-1 parameter) for a set of meteorological conditions,
#' # a leaf characteristic dimension of 1cm, and an LAI of 5
#' df = DataFrame(Tair=25,pressure=100,wind=c(3,4,5),ustar=c(0.5,0.6,0.65),H=c(200,230,250))
#' Gb_Su(data=df,zh=25,zr=40,d=17.5,Dl=0.01,LAI=5)
#'
#' # the same meteorological conditions, but larger leaves
#' Gb_Su(data=df,zh=25,zr=40,d=17.5,Dl=0.1,LAI=5)
#'
#' # same conditions, large leaves, and sparse canopy cover (LAI = 1.5)
#' Gb_Su(data=df,zh=25,zr=40,d=17.5,Dl=0.1,LAI=1.5)
#'
"""
"""
function Gb_Su(data,Tair="Tair",pressure="pressure",ustar="ustar",wind="wind",
                  H="H",zh,zr,d,z0m=nothing,Dl,fc=nothing,LAI=nothing,N=2,Cd=0.2,hs=0.01,
                  stab_formulation=c(Val(:Dyer_1970),Val(:Businger_1971)),
                  Sc=nothing,Sc_name=nothing,constants=BigleafConstants())

  stab_formulation = match_arg(stab_formulation)

  check_input(data,list(Tair,pressure,ustar,wind,H))

  if (isnothing(fc))
    if (isnothing(LAI))
      stop("one of 'fc' or 'LAI' must be provided",call.=false)
else
      fc = (1-exp(-LAI/2))
end
end

  if (isnothing(z0m))
    estimate_z0m = true
    z0m = nothing
else
    estimate_z0m = false
end

  wind_zh = wind_profile(data=data,z=zh,Tair=Tair,pressure=pressure,ustar=ustar,H=H,
                          zr=zr,estimate_z0m=estimate_z0m,zh=zh,d=d,z0m=z0m,frac_{z0m}=nothing,
                          stab_correction=true,stab_formulation=stab_formulation)

  v   = kinematic_viscosity(Tair,pressure,constants)
  Re  = Reynolds_Number(Tair,pressure,ustar,hs,constants)
  kBs = 2.46 * (Re)^0.25 - log(7.4)
  Reh = Dl * wind_zh / v
  Ct  = 1*constants.Pr^-0.6667*Reh^-0.5*N

  kB_h = (constants.k*Cd)/(4*Ct*ustar/wind_zh)*fc^2 + kBs*(1 - fc)^2
  Rb_h = kB_h/(constants.k*ustar)
  Gb_h = 1/Rb_h

  if (!isnothing(Sc) | !isnothing(Sc_name))
    if (length(Sc) != length(Sc_name))
      stop("arguments 'Sc' and 'Sc_name' must have the same length")
end
    if (!is_numeric(Sc))
      stop("argument 'Sc' must be numeric")
end
end

  Sc   = c(constants.Sc_CO2,Sc)
  Gb_x = DataFrame(lapply(Sc,function(x) Gb_h / (x/constants.Pr)^0.67))
  colnames(Gb_x) = paste0("Gb_",c("CO2",Sc_name))

  return(DataFrame(Gb_h,Rb_h,kB_h,Gb_x))
end


text:
"""
    RelinKeys

Stores relinearization keys generated by a `KeyGenerator` instance.

See also: [`KeyGenerator`](@ref)
"""
mutable struct RelinKeys <: SEALObject
  handle::Ptr{Cvoid}

  function RelinKeys()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    # RelinKeys are created as KSwitchKeys since they share the same data
    retval = ccall((:KSwitchKeys_Create1, libsealc), Clong,
                   (Ref{Ptr{Cvoid}},),
                   handleref)
    @check_return_value retval
    return RelinKeys(handleref[])
  end

  function RelinKeys(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::RelinKeys)
  if isallocated(object)
    @check_return_value ccall((:KSwitchKeys_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function parms_id(key::RelinKeys)
  parms_id = zeros(UInt64, 4)
  retval = ccall((:KSwitchKeys_GetParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 key, parms_id)
  @check_return_value retval
  return parms_id
end

function save_size(compr_mode, key::RelinKeys)
  result = Ref{Int64}(0)
  retval = ccall((:KSwitchKeys_SaveSize, libsealc), Clong,
                 (Ptr{Cvoid}, UInt8, Ref{Int64}),
                 key, compr_mode, result)
  @check_return_value retval
  return Int(result[])
end
save_size(key::RelinKeys) = save_size(ComprModeType.default, key)

function save!(buffer::DenseVector{UInt8}, length::Integer,
               compr_mode::ComprModeType.ComprModeTypeEnum, key::RelinKeys)
  out_bytes = Ref{Int64}(0)
  retval = ccall((:KSwitchKeys_Save, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, UInt8, Ref{Int64}),
                 key, buffer, length, compr_mode, out_bytes)
  @check_return_value retval
  return Int(out_bytes[])
end
function save!(buffer::DenseVector{UInt8}, length::Integer, key::RelinKeys)
  return save!(buffer, length, ComprModeType.default, key)
end
function save!(buffer::DenseVector{UInt8}, key::RelinKeys)
  return save!(buffer, length(buffer), key)
end

function load!(key::RelinKeys, context::SEALContext, buffer::DenseVector{UInt8}, length)
  in_bytes = Ref{Int64}(0)
  retval = ccall((:KSwitchKeys_Load, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt8}, UInt64, Ref{Int64}),
                 key, context, buffer, length, in_bytes)
  @check_return_value retval
  return Int(in_bytes[])
end
load!(key::RelinKeys, context::SEALContext, buffer::DenseVector{UInt8}) = load!(key,
                                                                                context,
                                                                                buffer,
                                                                                length(buffer))


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Expectations for the parser

"""
    expect_keyword(stream::InputStream, keyword::Symbol)

Read a token from an [`InputStream`](@ref) and check that it is a given [`Keyword`](@ref).
"""
function expect_keyword(stream::InputStream, keyword::Symbol)
    token = read_token(stream)
    isa(token.value, Keyword) || throw(WrongTokenType(token.loc,
                                                      "Expected a keyword instead of '$(typeof(token.value))'\nValid keyword: $keyword",
                                                      token.length))
    token.value.value == keyword || throw(InvalidKeyword(token.loc,
                                                         "Invalid '$(token.value.value)' keyword\nValid keyword: $keyword",
                                                         token.length))
    token
end

"""
    expect_keyword(stream::InputStream, keywords_list::Union{NTuple{N, Symbol} where {N}, AbstractVector{Symbol}})

Read a token from an [`InputStream`](@ref) and check that it is a [`Keyword`](@ref) in `keywords_list`.
"""
function expect_keyword(stream::InputStream, keywords_list::Union{NTuple{N, Symbol} where {N}, AbstractVector{Symbol}})
    token = read_token(stream)
    isa(token.value, Keyword) || throw(WrongTokenType(token.loc,
                                                      "Expected a keyword instead of '$(typeof(token.value))'\nValid keywords:\n\t$(join(keywords_list, "\n\t"))",
                                                      token.length))
    token.value.value ∈ keywords_list || throw(InvalidKeyword(token.loc,
                                                              "Invalid '$(token.value.value)' keyword\nValid keywords:\n\t$(join(keywords_list, "\n\t"))",
                                                              token.length))
    token
end

"""
    expect_command(stream::InputStream)

Read a token from an [`InputStream`](@ref) and check that it is a [`Command`](@ref).
"""
function expect_command(stream::InputStream)
    token = read_token(stream)
    isa(token.value, Command) || throw(WrongTokenType(token.loc,
                                                      "Expected a command instead of '$(typeof(token.value))'",
                                                      token.length))
    token
end

"""
    expect_command(stream::InputStream, command::Command)

Read a token from an [`InputStream`](@ref) and check that it is a given [`Command`](@ref).
"""
function expect_command(stream::InputStream, command::Command)
    token = expect_command(stream)
    token.value == command || throw(InvalidCommand(token.loc,
                                                  "Invalid command '$(token.value)'\nValid command: $command",
                                                  token.length))
    token
end

"""
    expect_command(stream::InputStream, commands::Union{NTuple{N, Command} where {N}, AbstractVector{Command}})

Read a token from an [`InputStream`](@ref) and check that it is a [`Command`](@ref) in the given `commands`.
"""
function expect_command(stream::InputStream, commands::Union{NTuple{N, Command} where {N}, AbstractVector{Command}})
    token = expect_command(stream)
    token.value ∈ commands || throw(InvalidCommand(token.loc,
                                                  "Invalid command '$(token.value)'\nValid commands:\n\t$(join(commands, "\n\t"))",
                                                  token.length))
    token
end

"""
    expect_type(stream::InputStream)

Read a token from an [`InputStream`](@ref) and check that it is a [`LiteralType`](@ref).
"""
function expect_type(stream::InputStream)
    token = read_token(stream)
    isa(token.value, LiteralType) || throw(WrongTokenType(token.loc,
                                                          "Expected a type instead of '$(typeof(token.value))'",
                                                          token.length))
    token
end

"""
    expect_type(stream::InputStream, type::LiteralType)

Read a token from an [`InputStream`](@ref) and check that it is a given [`LiteralType`](@ref).
"""
function expect_type(stream::InputStream, type::LiteralType)
    token = expect_type(stream)
    token.value == type || throw(WrongValueType(token.loc, "Expected type '$type', got: '$(token.value)'", token.length))
    token
end

"""
    expect_type(stream::InputStream, types::Union{NTuple{N, LiteralType} where {N}, AbstractVector{LiteralType}})

Read a token from an [`InputStream`](@ref) and check that it is a [`LiteralType`](@ref) in the given `types`.
"""
function expect_type(stream::InputStream, types::Union{NTuple{N, LiteralType} where {N}, AbstractVector{LiteralType}})
    token = expect_type(stream)
    token.value ∈ types || throw(WrongValueType(token.loc,
                                               "Invalid type '$(token.value)'\nValid types:\n\t$(join(types, "\n\t"))",
                                               token.length))
    token
end

"""
    expect_identifier(stream::InputStream)

Read a token from an [`InputStream`](@ref) and check that it is an [`Identifier`](@ref).
"""
function expect_identifier(stream::InputStream)
    token = read_token(stream)
    isa(token.value, Identifier) || throw(WrongTokenType(token.loc,
                                                         "Got token '$(typeof(token.value))' instead of 'Identifier'",
                                                         token.length))
    token
end

"""
    expect_string(stream::InputStream)

Read a token from an [`InputStream`](@ref) and check that it is a [`LiteralString`](@ref).
"""
function expect_string(stream::InputStream)
    token = read_token(stream)
    isa(token.value, LiteralString) && return token
    throw(WrongTokenType(token.loc,"Got token '$(typeof(token.value))' instead of 'LiteralString'", token.length))
end

"""
    expect_number(stream::InputStream, scene::Scene)

Read a token from an [`InputStream`](@ref) and check that it is a [`LiteralNumber`](@ref).
"""
function expect_number(stream::InputStream, scene::Scene)
    token = read_token(stream)
    isa(token.value, LiteralNumber) && return token
    throw(WrongTokenType(token.loc,
                         "Got '$(typeof(token.value))' instead of 'LiteralNumber'",
                         token.length))
end

"""
    expect_symbol(stream::InputStream, symbol::LiteralSymbol)

Read a token from an [`InputStream`](@ref) and check that it is the requested [`LiteralSymbol`](@ref).
"""
function expect_symbol(stream::InputStream, symbol::Symbol)
    token = read_token(stream)
    isa(token.value, LiteralSymbol) || throw(WrongTokenType(token.loc,
                                                            "Expected a symbol instead of '$(typeof(token.value))'\nValid symbol: $symbol",
                                                            token.length))

    token.value.value == symbol || throw(InvalidSymbol(token.loc,
                                                       "Invalid symbol '$(token.value.value)' \nValid symbol: $symbol",
                                                       token.length))
    token
end

"""
    expect_symbol(stream::InputStream, symbols::Union{Tuple{N, Symbol} where {N}, AbstractVector{Symbol}})

Read a token from an [`InputStream`](@ref) and check that it is one of the requested [`LiteralSymbol`](@ref)s.
"""
function expect_symbol(stream::InputStream, symbols::Union{Tuple{N, Symbol} where {N}, AbstractVector{Symbol}})
    token = read_token(stream)
    isa(token.value, LiteralSymbol) || throw(WrongTokenType(token.loc,
                                                            "Expected a symbol instead of '$(typeof(token.value))'\nValid symbols:\n\t$(join(symbols, "\n\t"))",
                                                            token.length))

    token.value.value ∈ symbols || throw(InvalidSymbol(token.loc,
                                                       "Invalid symbol '$(token.value.value)'\nValid symbols:\n\t$(join(symbols, "\n\t"))",
                                                       token.length))
    token
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

pathToFMU = "https://github.com/ThummeTo/FMI.jl/raw/main/model/" * ENV["EXPORTINGTOOL"] * "/SpringPendulum1D.fmu"

myFMU = fmi2Load(pathToFMU)
comp = fmi2Instantiate!(myFMU; loggingOn=false)
@test comp != 0

@test fmi2SetupExperiment(comp) == 0
@test fmi2EnterInitializationMode(comp) == 0
@test fmi2ExitInitializationMode(comp) == 0

targetValues = [[0.0, -10.0], [1.0, 0.0]]
dir_ders_buffer = zeros(fmi2Real, 2)
sample_ders_buffer = zeros(fmi2Real, 2, 1)
for i in 1:fmi2GetNumberOfStates(myFMU)

    if fmi2ProvidesDirectionalDerivative(myFMU)
        # multi derivatives calls
        sample_ders = fmi2SampleDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi2SampleDirectionalDerivative!(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], sample_ders_buffer)

        @test sum(abs.(sample_ders[:,1] - targetValues[i])) < 1e-3
        @test sum(abs.(sample_ders_buffer[:,1] - targetValues[i])) < 1e-3

        dir_ders = fmi2GetDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmi2GetDirectionalDerivative!(comp, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], dir_ders_buffer)

        @test sum(abs.(dir_ders - targetValues[i])) < 1e-3
        @test sum(abs.(dir_ders_buffer - targetValues[i])) < 1e-3

        # single derivative call
        dir_der = fmi2GetDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences[1], myFMU.modelDescription.stateValueReferences[1])
        @test dir_der == targetValues[1][1]
    else
        @warn "Skipping directional derivative testing, FMU from $(ENV["EXPORTINGTOOL"]) doesn't support directional derivatives."
    end
end

if ENV["EXPORTINGTOOL"] != "OpenModelica/v1.17.0"
    jac = fmi2GetJacobian(comp, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
    @test jac ≈ hcat(targetValues...)

    jac = fmi2SampleDirectionalDerivative(comp, myFMU.modelDescription.derivativeValueReferences, myFMU.modelDescription.stateValueReferences)
    @test jac ≈ hcat(targetValues...)
end

fmi2Unload(myFMU)


text:
function AngularRaDec(arg0::GroundStation, arg1::Frame, arg2::AbsoluteDate, arg3::Vector{jdouble}, arg4::Vector{jdouble}, arg5::Vector{jdouble}, arg6::ObservableSatellite)
    return AngularRaDec((GroundStation, Frame, AbsoluteDate, Vector{jdouble}, Vector{jdouble}, Vector{jdouble}, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_reference_frame(obj::AngularRaDec)
    return jcall(obj, "getReferenceFrame", Frame, ())
end

function get_station(obj::AngularRaDec)
    return jcall(obj, "getStation", GroundStation, ())
end


text:
module SelectMixedModels

	using MixedModels, ArrayViews, DataArrays, DataFrames, Distributions, PDMats, StatsBase
    using Base.LinAlg.CHOLMOD: CholmodFactor, CholmodSparse, CholmodSparse!,
        chm_scale, CHOLMOD_SYM, CHOLMOD_L, CHOLMOD_Lt, CHOLMOD_P, CHOLMOD_ROW, solve
	using Base.LinAlg: Cholesky, Ac_ldiv_B!, A_rdiv_Bc!, chksquare, transpose!

export
    conAIC,		# Conditional AIC
    conAICbc,   # Conditional AIC with bias correction
    marAIC,     # marginal AIC
    marBIC,     # marginal BIC
    checkzerovar, # check zero variance and return reduced models
    lmmg,       # fit a linear mixed-effects model (LMM)
                # with PLSGeneral type only
    condll,     # compute conditional likelihood (LMM)
    traceHat    # Caculate Trace of Hat matrix in LMM

    abstract InfoCrit

    include("tracehat.jl")
    include("IC.jl")
    include("marIC.jl")
    include("conAIC.jl")
    include("conAICbc.jl")
    include("util.jl")

end # module

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
module MomentClosure

import Catalyst: species, params, reactions, speciesmap, paramsmap, numspecies,
				 numreactions, numparams, substoichmat, prodstoichmat

using ModelingToolkit
using Symbolics: value, var_from_nested_derivative

using SciMLBase, SciMLBase.EnsembleAnalysis
using DiffEqJump
using Random
using Distributions: Geometric

using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint
using SymbolicUtils: Symbolic, expand, simplify, operation, arguments
using SymbolicUtils: @rule, @acrule, isnotflat, flatten_term, istree

using OrderedCollections: OrderedDict
using Combinatorics: permutations
using TupleTools: sort
using Cumulants
using Latexify

using DocStringExtensions

export generate_central_moment_eqs, generate_raw_moment_eqs, bernoulli_moment_eqs,
       get_S_mat, propensities, ReactionSystemMod,
       species, params, speciesmap, paramsmap, numspecies, numreactions, numparams,
       moment_closure, deterministic_IC, JumpProblem,
	   get_raw_moments, get_central_moments, get_cumulants

# reexporting from ModelingToolkit & Symbolics
# needed for ReactionSystemMod definition
export @parameters, @variables

include("reaction_systems.jl")
include("moment_equations.jl")
include("symbolic.jl")
include("moment_convert.jl")
include("stochastic_stoichiometry.jl")
include("central_moment_equations.jl")
include("raw_moment_equations.jl")
include("closure_methods/closure.jl")
include("bernoulli.jl")
include("utils.jl")
include("closure_methods/zero_closure.jl")
include("closure_methods/normal_closure.jl")
include("closure_methods/log_normal_closure.jl")
include("closure_methods/poisson_closure.jl")
include("closure_methods/gamma_closure.jl")
include("closure_methods/derivative_matching.jl")
include("closure_methods/conditional_gaussian.jl")
include("closure_methods/conditional_derivative_matching.jl")
include("stochastic_simulation.jl")

end


text:
module JaCoP

using JavaCall
import MathOptInterface
import ConstraintProgrammingExtensions

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const CleverDicts = MOIU.CleverDicts
const CP = ConstraintProgrammingExtensions

# Check if the package has been built correctly.
if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include("../deps/deps.jl")
else
    error(
        "JaCoP not properly installed. Please run `Pkg.build(\"JaCoP\")` or `]build JaCoP`",
    )
end

if !@isdefined(libjacopjava)
    error(
        "JaCoP not properly built. There probably was a problem when running `Pkg.build(\"JaCoP\")` or `]build JaCoP`",
    )
end

# Initialise the package by setting the right parameters for Java. This assumes
# no other code uses JavaCall...
function __init__()
    if get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "") != "true"
        jacop_java_init()
    end
end

"""
    jacop_java_init(init_java::Bool=true)

Initialises the JVM to be able to use JaCoP. This function must be called once
per Julia process that uses JaCoP.

By default, this function automatically starts the JVM link (from JavaCall.jl).
If other parts of the application require access to the JVM, this
initialisation can be disabled by setting `init_java` to `false`.
In that case, this function *MUST* be called before `JavaCall.init()`, because
it sets Java's CLASSPATH to include JaCoP.
"""
function jacop_java_init(init_java::Bool=true)
    JavaCall.addClassPath(libjacopjava)
    if init_java
        JavaCall.init()
    end
    return
end

include("java_wrapper.jl")
include("MOI/wrapper.jl")
include("MOI/parse.jl") # Must come after wrapper.jl.
include("MOI/wrapper_constraints_cp.jl")
include("MOI/wrapper_constraints_mo.jl")
include("MOI/wrapper_constraints_singlevar.jl")
include("MOI/wrapper_constraints.jl")
include("MOI/wrapper_variables.jl")

end


text:
#=
Given a string s and an integer k, break up the string into multiple lines such that each line has a length of k or less. You must break it up so that words don't break across lines. Each line has to have the maximum possible amount of words. If there's no way to break the text up, then return null.

You can assume that there are no spaces at the ends of the string and that there is exactly one space between each word.

For example, given the string "the quick brown fox jumps over the lazy dog" and k = 10, you should return: ["the quick", "brown fox", "jumps over", "the lazy", "dog"]. No string in the list has a length of more than 10.
=#

using Test
include("Solutions/problem57_string_into_multiplelines.jl")

@test breakstring("the quick brown fox jumps over the lazy dog", 10) == ["the quick", "brown fox", "jumps over", "the lazy", "dog"]

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
import Pkg
using Test

Pkg.activate(".")

using MatrixOptim

include("./milp.jl")
# include("./benders/milp.jl")
# include("./benders/l-shaped.jl")

# @test test_benders() == 0


text:
module Transform
include("constant.jl")

using .Constant: NONO_SIZE

"""
Transform the row permutations to column format, which start from leftmost to rightmost
Return transformed column
"""
function transform_column!(permRows::Union{SubArray,Array})::Array{Array{Int,1},1}
    cols = [Array{Int,1}(undef, 0) for i in 1:NONO_SIZE]
    accumulator = fill(0, NONO_SIZE)
    for permRow in permRows
        for i in 1:NONO_SIZE
            if (permRow & (1 << (i - 1))) != 0
                accumulator[i] += 1
            elseif accumulator[i] > 0
                push!(cols[i], accumulator[i])
                accumulator[i] = 0
            end
        end
    end
    for (i, acc) in enumerate(accumulator)
        if acc != 0
            push!(cols[i], acc)
        end
    end
    return reverse!(cols)
end
end


text:
#
using StreamModels
using Test
using DataStreams

using StreamModels: @formula, set_schema!, modelmatrow!

source = Data.Table((a = collect(1:10),
                     b = rand(10),
                     c = repeat(["a", "b"], inner=5)))

f = @formula ~ 1+a
set_schema!(f, Data.schema(source))

iter = Data.rows(source)
rows = collect(iter)

modelmatrow!(zeros(2), rows[end], (f.terms...))

modelmatrow!(zeros(2), rows[end], (f.terms...))
@code_warntype modelmatrow!(zeros(2), rows[end], (f.terms...))
@code_typed modelmatrow!(zeros(2), rows[end], (f.terms...))
@code_llvm modelmatrow!(zeros(2), rows[end], (f.terms...))


text:
using Test
using RoundingIntegers
@test isempty(detect_ambiguities(RoundingIntegers))

@testset "Basics" begin
    r16 = RInt16(3)
    @test isa(r16, RInt16)
    r = RInt(5)
    @test isa(r, RInt)
    x = Integer(r)
    @test x === 5
    @test isa(convert(RUInt8, 2), RUInt8)
    @test -7 % RUInt8 === RUInt8(0xf9)
    @test convert(RInteger, 3) === RInt(3)
    @test convert(RInteger, 0x03) === RUInt8(3)

    @test typemin(RInt16) === RInt16(typemin(Int16))
    @test typemin(RInt32) === RInt32(typemin(Int32))
    @test widen(r16) === convert(RInt32, r16)

    @test RInt16(r16) === r16
    @test float(r16) === Float64(r16) === 3.0
    @test convert(Float32, r) === 5.0f0

    @test !(r < r)
    @test r <= r
    @test !(r > r)
    @test r >= r
    @test !isless(r, r)

    @test isa(+r, RInt) && +r ==  5
    @test isa(-r, RInt) && -r == -5
    @test ~r16 === RInt16(-4)
    @test r << 2 === RInt(20)
    @test r >> 1 === RInt(2)
    @test r >>> 2 === RInt(1)
    @test RInt16(-5) >>> 3  === RInt16(8191)

    @test r16+r16 === RInt16(6)
    @test r16+r === 8
    @test r-r === RInt(0)
    @test r-r16 === 2
    @test r*r === RInt(25)
    @test 2*r === 10
    @test r*r16 === 15
    @test r/r === 1.0
    @test inv(r16) === inv(Int16(3))
    @test xor(r, r) === RInt(0)
    @test r & r === r
    @test r | r === r
    @test r ÷ r === RInt(1)
    @test r ÷ 5 === 1
    @test r % r === RInt(0)
    @test r % 5 === 0

    @test fld(r, r) === RInt(1)
    @test cld(r, r) === RInt(1)

    @test isodd(r)
    @test !iseven(r16)
    @test !signbit(r)
    @test copysign(r, r) === r
    @test unsigned(r) === RUInt(5)
    @test string(r, base=16) == "5"

    @test string(RInt(7.2)) == "7"
end

@testset "Rounding" begin
    @test RInt(17.3) == 17
    @test RInt(17.8) == 18
    a = RUInt8[2,3]
    a[1] = 5.7
    @test a[1] === RUInt8(6)
end

nothing


text:
using OpenQuantumBase, Test

c = ConstantCouplings(["ZI", "IZ"])
@test isequal(c.mats[1], 2π*σz⊗σi)
@test c.mats[2] == 2π*σi⊗σz
@test c[2](2.0) == 2π*σi⊗σz
res = c(0.2)
@test isequal(res[1], 2π*σz⊗σi)
@test res[2] == 2π*σi⊗σz

crot = rotate(c, σx⊗σi)
@test crot[1](0) ≈ -2π*σz⊗σi
@test crot[2](0.5) ≈ 2π*σi⊗σz

c = ConstantCouplings(["ZI", "IZ"], unit=:ħ)
@test isequal(c.mats[1], σz⊗σi)
@test isequal(c.mats[2], σi⊗σz)
res = c(0.2)
@test res[1] == σz⊗σi
@test res[2] == σi⊗σz
@test [op(0) for op in c] == [σz⊗σi, σi⊗σz]

c = ConstantCouplings([σz⊗σi, σi⊗σz], unit=:ħ)
@test isequal(c.mats[1], σz⊗σi)
@test isequal(c.mats[2], σi⊗σz)
res = c(0.2)
@test res[1] == σz⊗σi
@test res[2] == σi⊗σz
@test [op(0) for op in c] == [σz⊗σi, σi⊗σz]

c = ConstantCouplings(["ZI", "IZ"], sp=true)
@test isequal(c.mats[1], 2π*spσz⊗spσi)
@test isequal(c.mats[2], 2π*spσi⊗spσz)

c1 = TimeDependentCoupling([(s)->s], [σz], unit=:ħ)
@test c1(0.5) == 0.5σz
c2 = TimeDependentCoupling([(s)->s], [σx], unit=:ħ)
c = TimeDependentCouplings(c1, c2)
@test size(c) == (2, 2)
@test [op for op in c(0.5)] == [c1(0.5), c2(0.5)]
@test [op(0.5) for op in c] == [c1(0.5), c2(0.5)]

c = collective_coupling("Z", 2, unit=:ħ)
@test isequal(c(0.1), [σz⊗σi, σi⊗σz])

test_coupling = [(s)->s*σx, (s)->(1-s)*σz]
coupling = CustomCouplings(test_coupling, unit=:ħ)
@test size(coupling) == (2, 2)
@test coupling(0.5) == 0.5 * [σx, σz]
@test [c(0.2) for c in coupling] == [0.2*σx, 0.8*σz]


text:
using Plots
using Distributions

#=
implement function to perform Gibbs sampling for for a linear simple
linear model y = ax + e
=#

a = 2.0
σ = 5.0
t = range(1, stop = 10, length=20)

f(a, σ, t) = a .* t .+ (rand(Normal(0,σ),20))

y = f(a, σ, t)

scatter(t, y)

function Gibbs_Sampler!(y, σ, t, n, r_a, r_s)

    aₗ, aᵤ = r_a
    r_a = range(aₗ, stop = aᵤ, length = n)

    lₗ, lᵤ = r_s
    r_s = range(lₗ, stop = lᵤ, length = n)

    posterior = zeros(n, n)
    for i = 1:n
        for j = 1:n
            x = r_a[j] .* t
            likelihood = r_s[i]^(-length(t)) .* exp.(-sum((y .- x ).^2,dims=1)/2/r_s[i]^2)
            prior = 1/r_s[i] * pdf(Normal(r_a[j],1.0), 0.0)
            posterior[i,j] = likelihood[1] * prior
            posterior = posterior/sum(posterior)
            return likelihood, prior, posterior
        end
    end
end

n = 100
r_a = 0.0, 5.0
r_s = 0.0, 10

likelihood, prior, posterior = Gibbs_Sampler!(y, σ, t, n, r_a, r_s)

heatmap(posterior)


text:
global mpi = false
global mpi_np = 1
global mpi_tmp = ""
global mpi_disks = ""

function init_mpi(f::Function, np::Int, tmp::AbstractString)
    global mpi = true
    global mpi_np = np
    global mpi_tmp = "$tmp/$(DateTime(now()))"
    ENV["MPI_TMP"] = mpi_tmp
    mkdir(mpi_tmp)
    dirs = []
    for i = 1:np
        mkdir(@sprintf("%s/%03d", mpi_tmp, i-1))
        push!(dirs, "'$mpi_tmp'")
    end

    global mpi_disks = join(dirs, '\n')

    res = f()

    deinit_mpi()

    res
end

function deinit_mpi()
    isdir(mpi_tmp) && rm(mpi_tmp, recursive=true)
    mpi = false
    mpi_np = 1
end

function write_mpi_disks()
    open("disks", "w") do file
        write(file, "'$(pwd())'\n")
        write(file, mpi_disks)
    end
end

function clean_mpi_tmp(patterns...)
    mpi || return
    files = readdir(mpi_tmp)
    for p in patterns
        for f in filter(f -> ismatch(p,f), files)
            rm("$mpi_tmp/$f")
        end
    end
end

function mpi_cmd(cmd::AbstractString)
    # If MPI is enabled, return command for running MPI variant of
    # exe, else, return normal exe
    cmd = (ismatch(r"_mpi$", cmd) ?
           (mpi ?
            `mpirun -np $mpi_np $cmd`
            : cmd[1:end-4])
           : cmd)
    `$cmd`
end

function mpi_run(cmd::Union{Cmd,Base.CmdRedirect})
    mpi && write_mpi_disks()
    run(cmd)
end

export init_mpi, write_mpi_disks, clean_mpi_tmp, mpi_cmd, mpi_run


text:
## Interface to the Metis graph-partitioning algorithms,
##   http://glaros.dtc.umn.edu/gkhome/views/metis
## Version 5.1.0 or later of the metis library is required.
## The library should be compiled with IDXTYPEWIDTH = 32 (the default)

module Metis
    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")      # define libmetis and check it can be dlopen'd
    else
        error("Metis package not properly installed. Please run Pkg.build(\"Metis\")")
    end

    using Graphs,Compat                 # for AdjacencyList types
    using LightGraphs                   # metisform
    export
        nodeND,                         # determine fill-reducing permutation
        vertexSep,                      # single separator
        partGraphKway,
        partGraphRecursive

    include("metis_h.jl")               # define constants
    include("util.jl")                  # metisform and testgraph functions

    const metis_options = -ones(Int32, METIS_NOPTIONS) # defaults indicated by -1

    function nodeND(x,verbose::Integer=0)
        n,xadj,adjncy = metisform(x)
        metis_options[METIS_OPTION_DBGLVL] = verbose
        perm = Array(Cint, n)
        iperm = Array(Cint, n)
        err = ccall((:METIS_NodeND,libmetis), Cint,
                    (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                     Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                    &n, xadj, adjncy, C_NULL, metis_options, perm, iperm)
        err == METIS_OK || error("METIS_NodeND returned error code $err")
        perm .+ one(Cint), iperm .+ one(Cint)
    end

    function vertexSep(x,verbose::Integer=0)
        n,xadj,adjncy = metisform(x)

        metis_options[METIS_OPTION_DBGLVL] = verbose
        sepSize = zeros(Cint, 1)
        part = Array(Cint, n)

        err = ccall((:METIS_ComputeVertexSeparator,libmetis), Cint,
                    (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                     Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                    &n, xadj, adjncy, C_NULL, metis_options,
                    sepSize, part)
        err == METIS_OK || error("METIS_ComputeVertexSeparator returned error code $err")

        sizes = zeros(Cint,3)
        for i=1:n
          sizes[part[i]+1] += 1
        end
        sizes, part
    end

    function partGraphKway(x, nparts::Integer)
        n, xadj, adjncy = metisform(x)
        part = Array(Cint, n)
        objval = zeros(Cint, 1)
        err = ccall((:METIS_PartGraphKway,libmetis), Int32,
                    (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                     Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                     Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),
                    &n, &one(Int32), xadj, adjncy, C_NULL, C_NULL, C_NULL, &convert(Int32,nparts),
                    C_NULL, C_NULL, metis_options, objval, part)
        err == METIS_OK || error("METIS_PartGraphKWay returned error code $err")
        objval[1], part .+ one(Cint)
    end

    function partGraphRecursive(x, nparts::Integer)
        n, xadj, adjncy = metisform(x)
        part = Array(Int32, n)
        objval = zeros(Int32, 1)
        err = ccall((:METIS_PartGraphKway,libmetis), Int32,
                    (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                     Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
                     Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),
                    &n, &one(Int32), xadj, adjncy, C_NULL, C_NULL, C_NULL, &convert(Int32,nparts),
                    C_NULL, C_NULL, metis_options, objval, part)
        err == METIS_OK || error("METIS_PartGraphKWay returned error code $err")
        objval[1], part .+ one(Cint)
    end
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.Execute

qc1 = QCircuit(4)
qc1.x(2)
qc1.h(1)
qc1.cx(3, 2)
qc1.measure([1,2,3], [1,2,3])

qc2 = QCircuit(4)
qc2.x(2)
qc2.h(1)
qc2.cx(3, 2)
qc2.measure(1:3, 1:3)

@test qc1 == qc2

# Create circuit with barrier
qc1 = QCircuit(4)
qc1.x(2)
qc1.h(1)
qc1.barrier()
qc1.cx(3, 2)


text:
module PlacementClients
using gRPCClient
using ProtoBuf:google

include("dapr.jl")
using .dapr.proto.placement.v1

import Base: show

# begin service: dapr.proto.placement.v1.Placement

export PlacementBlockingClient, PlacementClient

struct PlacementBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::PlacementBlockingStub

    function PlacementBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = PlacementBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct PlacementClient
    controller::gRPCController
    channel::gRPCChannel
    stub::PlacementStub

    function PlacementClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = PlacementStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::PlacementBlockingClient) = print(io, "PlacementBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::PlacementClient) = print(io, "PlacementClient(", client.channel.baseurl, ")")

import .dapr.proto.placement.v1: ReportDaprStatus
"""
    ReportDaprStatus

- input: Channel{dapr.proto.placement.v1.Host}
- output: Channel{dapr.proto.placement.v1.PlacementOrder}
"""
ReportDaprStatus(client::PlacementBlockingClient, inp::Channel{dapr.proto.placement.v1.Host}) = ReportDaprStatus(client.stub, client.controller, inp)
ReportDaprStatus(client::PlacementClient, inp::Channel{dapr.proto.placement.v1.Host}, done::Function) = ReportDaprStatus(client.stub, client.controller, inp, done)

# end service: dapr.proto.placement.v1.Placement

end # module DaprClients


text:
using Profile
using ProfileView
using BenchmarkTools

# nn_lattice!(s)
#
# @code_warntype nn_lattice!(s)
#
@btime nn_lattice!(s)
rq = eigen(s.H)

#get_V_au_au(s)
mi = deepcopy(s.H)
aha = LAPACK.syev!('V', 'L', mi)
s.λ .= aha[1]
s.Γ .= aha[2]
#@code_warntype get_V_au_au(s)

#@btime get_V_au_au(s)

# get_e_au_n_neutral(s)
#
# @code_warntype get_e_au_n_neutral(s)
#
# @btime get_e_au_n_neutral(s)

# get_e_au_n_coup(s)
#
# @code_warntype get_e_ion(s)
#
# @btime get_e_all(s)
#
# hp = get_e_all(s)
#
# @profile get_e_all(s)
F_temp = zeros(Float64, N+2, 3)

get_F_all(s)

@benchmark get_F_au_au(s)

@code_warntype get_F_all(s)
#
@profiler  get_F_all(s)

@benchmark get_dm(s)

@code_warntype get_blk_akl!(s)

@benchmark get_blk_akl!(s)

get_dm!(s)

get_occupied_unoccpied_states!(s)

get_bkl!(s);

get_blk_akl!(s)

get_F_au_n_neutral(s)

get_F_au_o_neutral(s)

get_F_n_o_neutral(s)

get_F_au_n_ion(s)
get_F_au_o_ion(s)
get_F_n_o_ion(s)
F_au_n_c .= get_F_au_n_coup(s)
F_au_o_c .= get_F_au_o_coup(s)
F_coupling .= @. F_au_n_c + F_au_o_c

F_au .= get_F_au_au(s)
F_au_n_n .= get_F_au_n_neutral(s)
F_au_o_n .= get_F_au_o_neutral(s)
F_n_o_n_temp .= get_F_n_o_neutral(s)
F_n_o_n .= [F_n_o_n_temp; zeros(Float64, N, 3)]
F_neutral .= @. F_au_n_n + F_au_o_n + F_n_o_n - F_au

psi2 = collect([1.2358616421754511e-2,
  -1.3211748286157800e-3,
    0.99992134240773300 ,
           1.3936110728146633e-3,
              6.3955059627936444e-4,
                 3.3549186782848723e-4 ,
                   2.9694158126455239e-4,
      3.4448005730076524e-4,
     2.9591588060183586e-4,
    2.2393790827156967e-4,
    1.4074068340587945e-4])

vh1 = collect([
8.05440079e-03,
  0.999966562,
         1.22047227e-3,
            5.40201319e-4,
               3.10849166e-4 ,
                 1.74339104e-4,
                    1.57985167e-4,
                       1.89679558e-4,
                          1.68481318e-4,
                             1.30383385e-4,
                                8.28867342e-5
])


text:
using Random
# using PyCall
# import OpenAIGym

export render, close
import Base.close

include("utilities/OpenAIGym.jl")

mutable struct GymEnv <: AbstractEnvironment
    gymEnv::OpenAIGym.AbstractGymEnv
    # ExampleEnvironment() = new(rand(rng), rand(rng, ExampleEnvConst.state_initial_range), false)
    reward::Float64
    seed::Int64
    function GymEnv(name::Symbol, version::Symbol; rng=nothing, seed=0)
        new_gym = new(OpenAIGym.GymEnv(name, version), 0.0, seed)
        println(new_gym.gymEnv.pyenv.seed(PyObject(seed)))
        # pycall(new_gym.gymEnv.pyenv[:seed], PyObject(seed))
        return new_gym
    end
end

function reset!(env::GymEnv; rng = nothing, kwargs...)
    OpenAIGym.reset!(env.gymEnv)
end

function Base.close(env::GymEnv)
    close(env.gymEnv)
end

get_actions(env::GymEnv) = OpenAIGym.actions(env.gymEnv, nothing)

function environment_step!(env::GymEnv, action; rng = Random.GLOBAL_RNG, kwargs...)
    r, s = OpenAIGym.step!(env.gymEnv, action)
    env.reward = r
end

function get_reward(env::GymEnv) # -> determines if the agent_state is terminal
    return env.reward
end

function get_state(env::GymEnv) # -> get state of agent
    return env.gymEnv.state
end

function is_terminal(env::GymEnv) # -> determines if the agent_state is terminal
    return OpenAIGym.finished(env.gymEnv)
end

function get_normalized_state(env::GymEnv)
    return env.gymEnv.state
end

function Base.show(io::IO, env::GymEnv)
    println(io, get_state(env))
end

render(env::GymEnv, args...; kwargs...) =
    pycall(env.gymEnv.pyenv."render", PyAny; kwargs...)

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
using Polynomials
using TransportBasedInference: derivative, vander

@testset "Test probabilistic Hermite polynomials"   begin

x = randn()

P0 = ProPolyHermite(0)
@test P0(x) == 1.0
@test derivative(P0,0).coeffs == (1.0,)
@test derivative(P0,0).coeffs == (1.0,)
@test derivative(P0,1).coeffs == ()
@test derivative(P0,2).coeffs  == ()

@test (derivative(P0,1))(x) == 0.0
@test (derivative(P0,2))(x) == 0.0

P0 = ProPolyHermite(0; scaled = true)
@test P0(x) == 1/sqrt(1.0)
@test Cpro(0) == sqrt(1.0)
@test derivative(P0,0).coeffs == (1/sqrt(1.0),)
@test derivative(P0,0).coeffs == (1/sqrt(1.0),)
@test derivative(P0,1).coeffs == ()
@test derivative(P0,2).coeffs  == ()

P1 = ProPolyHermite(1)
@test P1(x) == x
@test derivative(P1, 0).coeffs == (0.0, 1.0)
@test derivative(P1, 1).coeffs == (1.0,)
@test derivative(P1, 2).coeffs == ()

P1 = ProPolyHermite(1; scaled = true)
@test abs(P1(x) - x/sqrt(factorial(1)))<1e-10
@test derivative(P1, 1).coeffs == (1/sqrt(factorial(1)),)
@test derivative(P1, 2).coeffs == ()

# Verify recurrence relation
# Hn+1(x) = x*Hn(x) - Hn′(x)
# Hn′(x) = n*Hn-1(x)

for i=1:11
    Pmm1 = ProPolyHermite(i-1)
    Pm = ProPolyHermite(i)
    Pmp1 = ProPolyHermite(i+1)
    z = 0.5

    @test abs(Pmp1(z) - (z*Pm(z) - (derivative(Pm,1))(z)))<1e-7

    @test abs(derivative(Pm, 0)(z) - Pm(z))<1e-7
    @test abs(derivative(Pm, 1)(z) - i*Pmm1(z))<1e-7

    @test abs(derivative(Pmp1, 2)(z) - (i+1)*(i)*Pmm1(z))<1e-7

    # Verify normalizing constant:
    @test abs(Cpro(i) - sqrt(factorial(i)))<1e-7
    @test abs(Cpro(Pm)- sqrt(factorial(i)))<1e-7

end

x = randn()

@test norm([ProPolyHermite(10).P.coeffs...] - [-945.0, 0.0, 4725.0, 0.0, -3150.0, 0.0, 630.0, 0.0, -45.0, 0.0, 1.0])<1e-10

P10 = ProPolyHermite(10)

@test abs(P10(-1.25) -  dot([-945.0, 0.0, 4725.0, 0.0, -3150.0, 0.0, 630.0, 0.0, -45.0, 0.0, 1.0], map(i->(-1.25)^i,0:10)))<1e-10

@test norm([derivative(ProPolyHermite(10),1).coeffs...] - [Polynomials.derivative(ProPolyHermite(10).P).coeffs...])<1e-7
@test norm([derivative(ProPolyHermite(10),2).coeffs...] - [Polynomials.derivative(Polynomials.derivative(ProPolyHermite(10).P)).coeffs...])<1e-7

# Compute inner product for normalized Probabilistic Hermite polynomials

for i=1:6
    Pmm1 = ProPolyHermite(i-1; scaled = true)
    Pm = ProPolyHermite(i; scaled = true)
    Pmp1 = ProPolyHermite(i+1; scaled = true)

    # Check derivatives
    z = 0.25
    @test abs(derivative(Pm, 0)(z) - Pm(z))<1e-7
    @test abs(derivative(Pm, 1)(z) - i*(1/Cpro(i))*FamilyProPolyHermite[i](z))<1e-7

    @test abs(derivative(Pmp1, 2)(z) - (i+1)*(i)*(1/Cpro(i+1))*FamilyProPolyHermite[i](z))<1e-7
end

P12 = ProPolyHermite(12;  scaled=true)
@test norm([derivative(P12,1).coeffs...] - [Polynomials.derivative(P12.P).coeffs...])<1e-10
@test norm([derivative(P12,2).coeffs...] - [Polynomials.derivative(Polynomials.derivative(P12.P)).coeffs...])<1e-10

end

@testset "Vandermonde matrix test for Probabilistic Hermite polynomials" begin
    # k = 0: check simple polynomial evaluation
    N = 20
    m = 10
    k = 0
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m), k, x)

    @test size(V)==(N, m+1)

    for i=0:m
        @test norm(V[:,i+1] - (FamilyProPolyHermite[i+1]).(x))<1e-8
    end

    # k = 0: check simple polynomial evaluation with scaled = true
    N = 20
    m = 10
    k = 0
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m;scaled = true), k, x)

    @test size(V)==(N, m+1)

    for i=0:m
        factor = exp(loggamma(i+1) - loggamma(i+1-k))

        @test norm(V[:,i+1] - (FamilyScaledProPolyHermite[i+1]).(x)*sqrt(factor))<1e-8
    end

    # k = 1: check first derivative
    N = 20
    m = 10
    k = 1
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m), k, x)

    @test size(V)==(N, m+1)

    for i=1:m
        @test norm(V[:,i+1] - (FamilyProPolyHermite[i]).(x)*i)<1e-8
    end

    # k = 1: check first derivative scaled = true
    N = 20
    m = 10
    k = 1
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m;scaled = true), k, x)

    @test size(V)==(N, m+1)

    for i=1:m
        factor = exp(loggamma(i+1) - loggamma(i+1-k))

        @test norm(V[:,i+1] - (FamilyProPolyHermite[i]).(x)*sqrt(factor)*1/Cpro(i))<1e-8
    end

    # k = 2: check second derivative
    N = 20
    m = 10
    k = 2
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m; scaled = false), k, x)

    @test size(V)==(N, m+1)

    for i=2:m
        @test norm(V[:,i+1] - (FamilyProPolyHermite[i-1]).(x)*i*(i-1))<1e-8
        @test norm(V[:,i+1] - (FamilyProPolyHermite[i-1]).(x)*factorial(i)/factorial(i-k))<1e-8

    end

    # k = 2: check second derivative scaled = true
    N = 20
    m = 10
    k = 2
    x = 0.2*randn(N)

    V = vander(ProPolyHermite(m;scaled = true), k, x)

    @test size(V)==(N, m+1)

    for i=2:m
        factor = exp(loggamma(i+1) - loggamma(i+1-k))

        @test norm(V[:,i+1] - (FamilyProPolyHermite[i-1]).(x)*sqrt(i*(i-1))*1/Cpro(i))<1e-8
        @test norm(V[:,i+1] - (FamilyProPolyHermite[i-1]).(x)*sqrt(factorial(i)/factorial(i-k))*1/Cpro(i))<1e-8

    end

    # # k = 3: check third derivative
    # N = 20
    # m = 10
    # k = 3
    # x = 0.2*randn(N)
    #
    # V = vander(ProPolyHermite(m; scaled = false), k, x)
    #
    # @test size(V)==(N, m+1)
    #
    # for i=3:m
    #     @test norm(V[:,i+1] - (FamilyProPolyHermite[i-2]).(x)*i*(i-1)*(i-2))<1e-6
    #     @test norm(V[:,i+1] - (FamilyProPolyHermite[i-2]).(x)*factorial(i)/factorial(i-k))<1e-6
    #
    # end
    #
    # # k = 3: check third derivative scaled = true
    # N = 20
    # m = 10
    # k = 3
    # x = 0.2*randn(N)
    #
    # V = vander(ProPolyHermite(m;scaled = true), k, x)
    #
    # @test size(V)==(N, m+1)
    #
    # for i=3:m
    #     factor = exp(loggamma(i+1) - loggamma(i+1-k))
    #
    #     @test norm(V[:,i+1] - (FamilyProPolyHermite[i-2]).(x)*sqrt(i*(i-1)*(i-2))*1/Cpro(i))<1e-8
    #     @test norm(V[:,i+1] - (FamilyProPolyHermite[i-2]).(x)*sqrt(factorial(i)/factorial(i-k))*1/Cpro(i))<1e-8
    #
    # end

end

@testset "Verify relation between PhyHermite and ProHermite" begin
# Hn(x) = 2^(n/2)Hen(√2*x)
# Hen(x) = 2^(-n/2)Hn(x/√2)
    x = -0.5:0.1:0.5

    for i=1:6
        Hi  = PhyPolyHermite(i)
        Hei = ProPolyHermite(i)
        for xi in x
            @test abs(Hi(xi)-(2^(i/2)*Hei(sqrt(2)*xi)))<1e-8
            @test abs(Hei(xi)-(2^(-i/2)*Hi(xi/sqrt(2))))<1e-8
        end
    end

end


text:
using Manifolds, Manopt, Test, Dates

@testset "Cyclic Proximal Point" begin
    @testset "Allocating" begin
        n = 100
        N = PowerManifold(Circle(), n)
        f = artificial_S1_signal(n)
        F(M, x) = costL2TV(M, f, 0.5, x)
        proxes = (
            (N, λ, x) -> prox_distance(N, λ, f, x), (N, λ, x) -> prox_TV(N, 0.5 * λ, x)
        )
        o = cyclic_proximal_point(
            N,
            F,
            proxes,
            f;
            λ=i -> π / (2 * i),
            stopping_criterion=StopAfterIteration(100),
            debug=[
                DebugIterate(), " ", DebugCost(), " ", DebugProximalParameter(), "\n", 10000
            ],
            record=[RecordProximalParameter(), RecordIterate(f), RecordCost()],
            return_options=true,
        )
        fR = get_solver_result(o)
        fR2 = cyclic_proximal_point(
            N, F, proxes, f; λ=i -> π / (2 * i), stopping_criterion=StopAfterIteration(100)
        )
        @test fR == fR2
        rec = get_record(o)
        @test F(N, f) > F(N, fR)
        o = CyclicProximalPointOptions(
            N, f; stopping_criterion=StopAfterIteration(1), λ=i -> π / (2 * i)
        )
        p = ProximalProblem(N, F, proxes, [1, 2])
        @test_throws ErrorException get_proximal_map(p, 1.0, f, 3)
        @test_throws ErrorException ProximalProblem(N, F, proxes, [1, 2, 2])
    end
    @testset "Mutating" begin
        n = 3
        M = Sphere(2)
        N = PowerManifold(M, NestedPowerRepresentation(), n)
        f = artificial_S2_lemniscate([0.0, 0.0, 1.0], n)
        F(N, x) = costL2TV(N, f, 0.5, x)
        proxes! = (
            (N, y, λ, x) -> prox_distance!(N, y, λ, f, x),
            (N, y, λ, x) -> prox_TV!(N, y, 0.5 * λ, x),
        )
        proxes = (
            (N, λ, x) -> prox_distance(N, λ, f, x), (N, λ, x) -> prox_TV(N, 0.5 * λ, x)
        )
        s1 = cyclic_proximal_point(
            N, F, proxes, f; λ=i -> π / (2 * i), stopping_criterion=StopAfterIteration(100)
        )
        s2 = cyclic_proximal_point(
            N,
            F,
            proxes!,
            f;
            λ=i -> π / (2 * i),
            stopping_criterion=StopAfterIteration(100),
            evaluation=MutatingEvaluation(),
        )
        @test isapprox(N, s1, s2)
    end
    @testset "Problem access functions" begin
        n = 3
        M = Sphere(2)
        N = PowerManifold(M, NestedPowerRepresentation(), n)
        f = artificial_S2_lemniscate([0.0, 0.0, 1.0], n)
        F(N, x) = costL2TV(N, f, 0.5, x)
        proxes! = (
            (N, y, λ, x) -> prox_distance!(N, y, λ, f, x),
            (N, y, λ, x) -> prox_TV!(N, y, 0.5 * λ, x),
        )
        proxes = (
            (N, λ, x) -> prox_distance(N, λ, f, x), (N, λ, x) -> prox_TV(N, 0.5 * λ, x)
        )
        for i in 1:2
            p1 = ProximalProblem(N, F, proxes)
            g = deepcopy(f)
            get_proximal_map!(p1, g, 1.0, g, i)
            @test isapprox(N, g, get_proximal_map(p1, 1.0, f, i))
            p2 = ProximalProblem(N, F, proxes!; evaluation=MutatingEvaluation())
            g = deepcopy(f)
            get_proximal_map!(p2, g, 1.0, g, i)
            @test isapprox(N, g, get_proximal_map(p2, 1.0, f, i))
        end
    end
end


text:
module ITensorTDVP

using ITensors
using KrylovKit: exponentiate, eigsolve
using Printf
using TimerOutputs
using Observers

using ITensors:
  AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, set_nsite!

# Compatibility of ITensor observer and Observers
include("update_observer.jl")

# Utilities for making it easier
# to define solvers (like ODE solvers)
# for TDVP
include("solver_utils.jl")

include("applyexp.jl")
include("tdvporder.jl")
include("tdvpinfo.jl")
include("tdvp_step.jl")
include("tdvp_generic.jl")
include("tdvp.jl")
include("dmrg.jl")
include("dmrg_x.jl")

export tdvp, dmrg_x, to_vec, TimeDependentSum

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots
p = heatmap(rand(100, 100);
    colorbar_title = "colorbar title",
    thickness_scaling = 1.6,
)

# %%
using Plots
p = heatmap(rand(100, 100);
    colorbar_title = " \ncolorbar title",
    thickness_scaling = 1.6,
    right_margin = 2Plots.mm,
)

# %%

max_stars_repo_path: 0017/colorbar_title example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0017/colorbar_title example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0017/colorbar_title example.jl

max_forks_repo_name: genkuroki/public


text:
using Plots

struct EMACS_DISPLAY <: AbstractDisplay end

import Base: display

display(d::EMACS_DISPLAY, obj::CLEAR) = println(stdout, "\n<emacs-clear></emacs-clear>")

# function display(d::EMACS_DISPLAY, p::Plots.Plot)
#     # Is this really needed?
#     Base.invokelatest(show, stdout, MIME"image/svg+xml"(), p)
# end

images_tempdir = nothing
function ImagesDir()
    global images_tempdir
    if images_tempdir === nothing
        images_tempdir = mktempdir()
    end
    images_tempdir
end


text:
###################################################
#Auxitial functions flat boundaries
###################################################
export boundStopMin, boundStopMax, boundBounceMin, boundBounceMax, boundReflect

boundStopMin(x,min) = max(x,min)
boundStopMax(x,max) = min(x,max)

boundBounceMin(x,min) = min + (min-x)
boundBounceMax(x,max) = max - (x-max)

boundReflect(x) = -x

function checkInitiallyInBoundBoundaryFlat(simulationBox,localV,boundaries)
    for (i,b) in enumerate(boundaries)
        if b in [Periodic,Bounded]
            if any(localV[:,i] .<= simulationBox[i,1]) || any(localV[:,1] .> simulationBox[i,2])
                error("Particles have to be initialised inside the simulation boundary min<X=<max for flat boundaries that are not of type Free.",
                    " Some particles of axis ", i ," with declared boundary ", b ," have been initialised outside it."
                )
            end
        end
    end
end

###################################################
#Return code boundaries
###################################################

#Free
function returnBound_(s::Symbol,pos::Int,b::Free,p::Program_)

    return quote end
end

#Periodic
function returnBound_(s::Symbol,pos::Int,b::Periodic,p::Program_)

    posS = p.update["Local"][s]

    up = quote localVCopy[ic1_,$posS] += simulationBox[$pos,2]-simulationBox[$pos,1] end
    for i in b.addSymbols["add"]
        posm = p.update["Local"][i]
        push!(up.args,:(localVCopy[ic1_,$posm] += simulationBox[$pos,2]-simulationBox[$pos,1]))
    end

    down = quote localVCopy[ic1_,$posS] -= simulationBox[$pos,2]-simulationBox[$pos,1] end
    for i in b.addSymbols["add"]
        posm = p.update["Local"][i]
        push!(down.args,:(localVCopy[ic1_,$posm] -= simulationBox[$pos,2]-simulationBox[$pos,1]))
    end

    code = quote
                if localVCopy[ic1_,$pos] < simulationBox[$pos,1]
                    $up
                elseif localVCopy[ic1_,$pos] >= simulationBox[$pos,2]
                    $down
                end
           end

    return code
end

#Bounded
function returnBound_(s::Symbol,pos::Int,b::Bounded,p::Program_)

    #Add main symbol if not present
    symbolPresent = false
    for i in keys(b.addSymbols)
        if s in b.addSymbols[i]
            symbolPresent = true
            break
        end
    end
    if !symbolPresent
        push!(b.addSymbols["stop"],s)
    end

    min = quote end
    for i in b.addSymbols["stop"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundStopMin(localVCopy[ic1_,$posm],simulationBox[$pos,1])))
    end
    for i in b.addSymbols["stopMin"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundStopMin(localVCopy[ic1_,$posm],simulationBox[$pos,1])))
    end
    for i in b.addSymbols["bounce"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundBounceMin(localVCopy[ic1_,$posm],simulationBox[$pos,1])))
    end
    for i in b.addSymbols["bounceMin"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundBounceMin(localVCopy[ic1_,$posm],simulationBox[$pos,1])))
    end
    for i in b.addSymbols["reflect"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundReflect(localVCopy[ic1_,$posm])))
    end
    for i in b.addSymbols["reflectMin"]
        posm = p.update["Local"][i]
        push!(min.args,:(localVCopy[ic1_,$posm] = boundReflect(localVCopy[ic1_,$posm])))
    end

    max = quote end
    for i in b.addSymbols["stop"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundStopMax(localVCopy[ic1_,$posm],simulationBox[$pos,2])))
    end
    for i in b.addSymbols["stopMax"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundStopMax(localVCopy[ic1_,$posm],simulationBox[$pos,2])))
    end
    for i in b.addSymbols["bounce"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundBounceMax(localVCopy[ic1_,$posm],simulationBox[$pos,2])))
    end
    for i in b.addSymbols["bounceMax"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundBounceMax(localVCopy[ic1_,$posm],simulationBox[$pos,2])))
    end
    for i in b.addSymbols["reflect"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundReflect(localVCopy[ic1_,$posm])))
    end
    for i in b.addSymbols["reflectMax"]
        posm = p.update["Local"][i]
        push!(max.args,:(localVCopy[ic1_,$posm] = boundReflect(localVCopy[ic1_,$posm])))
    end

    posm = p.update["Local"][s]

    if !emptyquote_(min) && !emptyquote_(max)

        code = quote
            if localVCopy[ic1_,$posm] < simulationBox[$pos,1]
                $min
            elseif localVCopy[ic1_,$posm] >= simulationBox[$pos,2]
                $max
            end
       end

       return code
    elseif !emptyquote_(min)
        code = quote
            if localVCopy[ic1_,$posm] < simulationBox[$pos,1]
                $min
            end
       end

        return code
    elseif !emptyquote_(max)
        code = quote
            if localVCopy[ic1_,$posm] >= simulationBox[$pos,2]
                $max
            end
       end

        return code
    else

        return quote end
    end

end

function returnBound_(b::BoundaryFlat,p::Program_)

    code = quote end

    for (i,bound) in enumerate(b.boundaries)
        s = POSITIONSYMBOLS[i]
        append!(code.args,returnBound_(s,i,bound,p).args)
    end

    return code
end

#######################################################################################
#Return computation over boundaries
#######################################################################################
function boundariesFunctionDefinition(b::BoundaryFlat,p::Program_, platform::String)

    #Make function to compute the boundaries
    code = quote end

    lUpdate = length(p.agent.declaredSymbols["Medium"])
    for i in 1:p.agent.dims #Adapt boundary updates depending on the boundary type
        ic = Meta.parse(string("ic",i,"_"))
        nm = Meta.parse(string("N",["x","y","z"][i],"_"))
        subcode = quote end

        if p.agent.dims == 1
            v = :(mediumVCopy[ic1_,ic4_])
            v2 = :(mediumV[ic1_,ic4_])
        elseif p.agent.dims == 2
            v = :(mediumVCopy[ic1_,ic2_,ic4_])
            v2 = :(mediumV[ic1_,ic2_,ic4_])
        elseif p.agent.dims == 3
            v = :(mediumVCopy[ic1_,ic2_,ic3_,ic4_])
            v2 = :(mediumV[ic1_,ic2_,ic3_,ic4_])
        end

        if b.boundaries[i].medium == PeriodicBoundaryCondition

            #Lower boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)
            vAss = postwalk(x->@capture(x, t_) && t == ic ? :($nm-1) : x, v)

            push!(subcode.args, :($vUp=$vAss))

            #Upper boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? :($nm) : x, v)
            vAss = postwalk(x->@capture(x, t_) && t == ic ? 2 : x, v)

            push!(subcode.args, :($vUp=$vAss))

        end

        xMin = [:(simulationBox[1,1]) :(ic2_*(simulationBox[2,2]-simulationBox[2,1])/Ny_+simulationBox[2,1]) :(ic3_*(simulationBox[3,2]-simulationBox[3,1])/Nz_+simulationBox[3,1]);
        :(ic1_*(simulationBox[1,2]-simulationBox[1,1])/Nx_+simulationBox[1,1]) :(simulationBox[2,1]) :(ic3_*(simulationBox[3,2]-simulationBox[3,1])/Nz_+simulationBox[3,1]);
        :(ic1_*(simulationBox[1,2]-simulationBox[1,1])/Nx_+simulationBox[1,1]) :(ic2_*(simulationBox[2,2]-simulationBox[2,1])/Ny_+simulationBox[2,1]) :(simulationBox[3,1])]

        xMax = [:(simulationBox[1,2]) :(ic2_*(simulationBox[2,2]-simulationBox[2,1])/Ny_+simulationBox[2,1]) :(ic3_*(simulationBox[3,2]-simulationBox[3,1])/Nz_+simulationBox[3,1]);
        :(ic1_*(simulationBox[1,2]-simulationBox[1,1])/Nx_+simulationBox[1,1]) :(simulationBox[2,2]) :(ic3_*(simulationBox[3,2]-simulationBox[3,1])/Nz_+simulationBox[3,1]);
        :(ic1_*(simulationBox[1,2]-simulationBox[1,1])/Nx_+simulationBox[1,1]) :(ic2_*(simulationBox[2,2]-simulationBox[2,1])/Ny_+simulationBox[2,1]) :(simulationBox[3,2])]

        if typeof(b.boundaries[i].medium) in [NewmannBoundaryCondition,NewmannBoundaryCondition_DirichletBoundaryCondition]

            #Lower boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)
            vUp2 = postwalk(x->@capture(x, t_) && t == ic ? 2 : x, v2) #Original medium variable

            x = xMin[i,1:p.agent.dims]
            push!(subcode.args, :($vUp=$(b.boundaries[i].medium.fMin)($(x...),t)*$vUp2+$vUp2))

        end

        if typeof(b.boundaries[i].medium) in [NewmannBoundaryCondition,DirichletBoundaryCondition_NewmannBoundaryCondition]

            #Upper boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? nm : x, v)
            vUp2 = postwalk(x->@capture(x, t_) && t == ic ? :($nm-1) : x, v2) #Original medium variable

            x = xMin[i,1:p.agent.dims]
            push!(subcode.args, :($vUp=$(b.boundaries[i].medium.fMax)($(x...),t)*$vUp2+$vUp2))

        end

        if typeof(b.boundaries[i].medium) in [DirichletBoundaryCondition,DirichletBoundaryCondition_NewmannBoundaryCondition]

            #Lower boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)

            x = xMin[i,1:p.agent.dims]
            push!(subcode.args, :($vUp=$(b.boundaries[i].medium.fMin)($(x...),t)))

        end

        if typeof(b.boundaries[i].medium) in [DirichletBoundaryCondition,NewmannBoundaryCondition_DirichletBoundaryCondition]

            #Upper boundary
            vUp = postwalk(x->@capture(x, t_) && t == ic ? nm : x, v)

            x = xMax[i,1:p.agent.dims]
            update = :($vUp=$(b.boundaries[i].medium.fMax)($(x...),t))

            push!(subcode.args, update)

        end

        if !(typeof(b.boundaries[i]) in [Bounded,Periodic])
            error("Boundary has to be a bounded boundary but some was set to a non-valid type: ", [typeof(i) for i in b.boundaries])
        end

        subcode = :(for ic4_ in 1:$lUpdate
                        $subcode
                    end
                )

        subcode = simpleGridLoop_(platform,subcode,p.agent.dims-1, indexes = [1,2,3][findall([1,2,3].!=i)])

        push!(code.args,subcode)

    end

    f = wrapInFunction_(:mediumBoundaryStep_!,code)

    push!(p.declareF.args,f)
end


text:
@testset "Testing modeltools.jl" begin

    @testset "Log Utility" begin
    # Test log utility
    u = LogUtility(2.0)
    @test isapprox(u(1.0), 0.0)
    @test isapprox(2.0*log(2.5), u(2.5))
    @test isapprox(u.(ones(5)), zeros(5))
    @test u(2.0) > u(1.0)  # Ensure it is increasing
    @test isapprox(derivative(u, 1.0), 2.0)
    @test isapprox(derivative(u, 2.0), 1.0)

    # Test "extrapolation evaluations"
    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
    @test u(-0.5) < u(-0.1)  # Make sure it doesn't fail with negative values

    @test isapprox(LogUtility().ξ, 1.0)  # Test default constructor

    end

    @testset "CRRA Utility" begin
    # Test CRRA utility
    u = CRRAUtility(2.0)
    @test isapprox(u(1.0), 0.0)
    @test isapprox((2.5^(-1.0) - 1.0) / (-1.0), u(2.5))
    @test isapprox(u.(ones(5)), zeros(5))
    @test u(5.0) > u(3.0)  # Ensure it is increasing
    @test isapprox(derivative(u, 1.0), 1.0)
    @test isapprox(derivative(u, 2.0), 0.25)

    # Test "extrapolation evaluations"
    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
    @test u(-0.5) < u(-0.1)  # Make sure it doesn't fail with negative values

    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at γ=1.0
    end

    @testset "CFE Utility" begin
    # Test CFE Utility
    v = CFEUtility(2.0)
    @test isapprox(v(1.0), -2.0/3.0)
    @test isapprox(v(0.5), -v.ξ * 0.5^(1.0 + 1.0/v.ϕ) / (1.0 + 1.0/v.ϕ))
    @test v(0.5) > v(0.85)  # Ensure it is decreasing
    @test isapprox(derivative(v, 0.25), -0.5)
    @test isapprox(derivative(v, 1.0), -1.0)

    # Test "extrapolation evaluations"
    @test isapprox(v(0.5e-10), v(1e-10) + derivative(v, 1e-10)*(0.5e-10 - 1e-10))
    @test v(-0.5) > v(-0.1)  # Make sure it doesn't fail with negative values

    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at ϕ=1.0
    end

    @testset "Elliptical Utility" begin
    # Test Elliptical Utility
    u = EllipticalUtility(1.0, 2.0)
    @test isapprox(u(sqrt(0.75)), 0.5)
    @test isapprox(u(sqrt(0.51)), 0.7)
    @test u(0.5) > u(0.95)  # Ensure it is decreasing
    @test isapprox(derivative(u, 0.0), 0.0)
    @test isapprox(derivative(u, sqrt(0.5)), -1.0)

    # Test default values
    @test isapprox(EllipticalUtility().b, 0.5223)
    @test isapprox(EllipticalUtility().μ, 2.2926)

    end

end


text:
# Line plots
import .LibCImPlot: PlotLine

function PlotLine(x::AbstractArray{T}, y::AbstractArray{T};
                  count::Integer = min(length(x), length(y)), offset::Integer = 0,
                  stride::Integer = 1, label_id::String = "") where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotLine(x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real, T2<:Real}
    PlotLine(promote(x,y)...; kwargs...)
end

function PlotLine(y::AbstractArray{T}; label_id::String="", count::Integer=length(y),
                  xscale::Real = 1.0, x0::Real = 0.0, offset::Integer=0,
                  stride::Integer=1) where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotLine(x::UnitRange{<:Integer}, y::AbstractArray{T}; xscale::Real = 1.0,
                  x0::Real = 0.0, label_id::String="") where {T<:ImPlotData}

    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride::Cint = sizeof(T)
    LibCImPlot.PlotLine(label_id, y, count,  xscale, x0, offset, stride)
end

function PlotLine(x::StepRange, y::AbstractArray{T}; xscale::Real = 1.0, x0::Real = 0.0,
                  label_id::String= "") where {T<:ImPlotData}

    x.stop < 1 && throw("Range out of bounds")
    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(x.step * sizeof(T))
    LibCImPlot.PlotLine(label_id, y, count, xscale, x0, offset, stride)
end

PlotLine(x::UnitRange{<:Integer}, y::AbstractArray{<:Real}; kwargs...) = PlotLine(x, Float64.(y); kwargs...)
PlotLine(x::StepRange, y::AbstractArray{<:Real}; kwargs...) = PlotLine(x, Float64.(y); kwargs...)


text:
module TestHistogramMSD
using LiterateTest.AssertAsTest: @assert
test() = include(joinpath(@__DIR__, "../examples/histogram_msd.jl"))
end  # module


text:
#
# XML support.
#

include("xml/libexpat.jl")
include("xml/driver.jl")
include("xml/combinator.jl")

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module GfxBase

export Vec3, Vec2, Ray, Vec4

using Images
using StaticArrays
import LinearAlgebra.cross
import LinearAlgebra.normalize

export RGB

# some type aliases:
const Vec3 = SVector{3, Float64} # 3-vector of floats
const Vec4 = SVector{4, Float64}
const Vec2 = SVector{2, Float64} # 2-vector of floats

struct Ray
    origin::Vec3
    direction::Vec3
end

end # module GfxBase


text:
using LaTeXStrings
using PyPlot
using Printf

function dsdemo1()
    println("NTF Synthesis-- 5th-order modulator")

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("NTF Poles and Zeros")

    order = 5
    OSR = 32
    opt = 0
    H = synthesizeNTF(order, OSR, opt)

    plotPZ(H)

    print("paused")
    readline()

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("NTF Magnitude Response")
    fig.subplots_adjust(hspace=0.5)

    subplot(211)

    f = [range(0, 0.75/OSR, length=100); range(0.75/OSR, 0.5, length=100)]
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))

    plot(f, magH, linewidth=1)
    figureMagic([0, 0.5], 0.05, 2, [-100, 10], 10, 2)
    xlabel(string("Normalize frequency ", L"(1 \rightarrow f_s)"))
    ylabel("dB")

    subplot(212)

    fstart = 0.01
    f = range(fstart, 1.2, length=200) / (2 * OSR)
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))
    sigma_H = dbv(rmsGain(H, 0, 0.5/OSR))

    semilogx(f*2*OSR, magH, linewidth=1)
    semilogx([fstart, 1], sigma_H*[1, 1], linewidth=1, marker="o", markerfacecolor="none", markeredgewidth=0.5)
    text(0.15, sigma_H+2, @sprintf("rms gain = %5.0fdB", sigma_H))
    axis([fstart, 1.2, -100, -30])
    grid()
    grid(which="minor", linestyle="dotted")
    xlabel(string("Normalize frequency ", L"(1 \rightarrow f_B)"))
    ylabel("dB")

    println("paused")
    readline()

    println("Optimized zeros")

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("NTF Poles and Optimized Zeros")

    opt = 1
    H = synthesizeNTF(order, OSR, opt)

    plotPZ(H)

    fig.canvas.manager.window.raise_()

    print("paused")
    readline()

    fig = figure(2)
    axs = fig.get_axes()

    f = [range(0, 0.75/OSR, length=100); range(0.75/OSR, 0.5, length=100)]
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))

    axs[1].plot(f, magH, linestyle="--", linewidth=1)

    f = range(fstart, 1.2, length=200) / (2 * OSR)
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))

    axs[2].semilogx(f*2*OSR, magH, linestyle="--", linewidth=1)

    sigma_H = dbv(rmsGain(H, 0, 0.5/OSR))

    axs[2].plot([fstart, 1], sigma_H*[1, 1], linestyle="--", linewidth=1, marker="o", markerfacecolor="none", markeredgewidth=0.5)
    axs[2].text(0.15, sigma_H+2, @sprintf("rms gain = %5.0fdB", sigma_H))

    fig.canvas.manager.window.raise_()

    println("paused")
    readline()

    println("NTF Synthesis-- Bandpass Modulator")

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("Bandpass NTF Poles and Zeros")

    order = 8
    OSR = 64
    opt = 2
    f0 = 0.125
    H = synthesizeNTF(order, OSR, opt, 1.5, f0)

    plotPZ(H)

    fig.canvas.manager.window.raise_()

    print("paused")
    readline()

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("Bandpass NTF/STF Magnitude Response")

    subplot(211)

    f = [range(0, f0-1/(2*OSR), length=50); range(f0-1/(2*OSR), f0+1/(2*OSR), length=100); range(f0+1/(2*OSR), 0.5, length=50)]
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))

    G = zpk(zeros(order÷2), H.matrix[1].p, 1, 1)
    G *= 1 / abs(evalTF(G, exp.(2im* pi * f0)))
    magG = dbv(evalTF(G, z))

    plot(f, magH, linewidth=1)
    plot(f, magG, color="r", linewidth=1)
    axis([0, 0.5, -100, 10])
    grid()
    xlabel(string("Normalize frequency ", L"(1 \rightarrow f_s)"))
    ylabel("dB")

    subplot(212)

    f = range(f0-0.3/OSR, f0+0.3/OSR, length=100)
    z = exp.(2im * pi * f)
    magH = dbv(evalTF(H, z))
    sigma_H = dbv(rmsGain(H, f0-0.25/OSR, f0+0.25/OSR))

    plot(2*OSR*(f .- f0), magH, linewidth=1)
    plot([-0.5, 0.5], sigma_H*[1, 1], marker="o", markerfacecolor="none", markeredgewidth=0.5)
    text(-0.2, sigma_H+2, @sprintf("rms gain = %5.0fdB", sigma_H))
    axis([-0.6, 0.6, -100, -60])
    grid()
    xlabel("Normalized frequency offset")
    ylabel("dB")

    fig.canvas.manager.window.raise_()

    return nothing
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
export User

type User <: AbstractModel
  _table_name::String
  _id::String

  id::Nullable{SearchLight.DbId}
  name::String
  email::String
  password::String
  role_id::Nullable{SearchLight.DbId}

  belongs_to::Vector{SearchLight.SQLRelation}

  User(;
    id = Nullable{SearchLight.DbId}(),
    name = "",
    email = "",
    password = "",
    role_id = Nullable{SearchLight.DbId}(),
    updated_at = DateTime(),

    belongs_to = [SearchLight.SQLRelation(Role)]
  ) = new("users", "id", id, name, email, password, role_id, belongs_to)
end

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
@diffusion_process SineAux{K,R,S} begin
    :parameters
    (a, b, c, σ) --> K

    :auxiliary_info
    t0 --> Float64
    T --> Float64
    v0 --> R
    vT --> S

    :additional
    constdiff --> true
    linear --> true
end

DiffusionDefinition.B(t, P::SineAux) = 0.0 + (t/P.T)/5.0
DiffusionDefinition.β(t, P::SineAux) = (P.vT-P.v0)/P.T
DiffusionDefinition.σ(t, x, P::SineAux) = P.σ

DiffusionDefinition.default_type(::SineAux) = Float64
DiffusionDefinition.default_wiener_type(::SineAux) = Float64


text:
module DH32Utils

export     block_average_spread,
    draw_mamba_model_graph,
    drawSVG,
    drawSVGandPNG,
    hclustplot,
    loess_fit,
    midpoints,
    save_script_and_fig_to_unique_filenames,
    screenareaPNG

using ImageMagick, FileIO, Loess, DataFrames, Clustering, StatPlots

"""
Takes a Mamba model (see Mamba.jl) and a string (stem name of files), and
draws an acyclic graph corresponding to the model.

NOTE: needs the GraphViz dot program!

Symbols in the graph:

- ellipses are stochastic variables

- squares are logic variables

- white symbols = monitored variables

- gray symbols = variables not monitored

- arrows = "generates"

"""

function draw_mamba_model_graph(model::Any, filestem::String)
    filename=filestem*".dot"
    draw(model,filename=filename)

    #we replace the "." in the model name because it causes
    #an error:
    m = replace(readstring(filename),"Mamba.Model","MambaModel")
    write(filename,m)

    #we use the original dot program because it plots with a better
    #layout:
    display("image/svg+xml",open(readstring,`dot $filename -Tsvg`))
end

"""
Plot recipe for cluster dendrogram, from
https://gist.github.com/kescobo/ -> hclustrecipe.jl
"""
function hclustplot(hc::Hclust, useheight::Bool)
    o = indexmap(hc.order)
    n = [x for x in 1:length(o)]

    pos = treepositions(hc, useheight)

    xs = []
    ys = []
    for i in 1: size(hc.merge, 1)
        x1 = pos[hc.merge[i,1]][1]
        x2 = pos[hc.merge[i,2]][1]
        append!(xs, [x1,x1,x2,x2])

        y1 = pos[hc.merge[i,1]][2]
        y2 = pos[hc.merge[i,2]][2]
        useheight ? h = hc.height[i] : h = 1
        newy = maximum([y1,y2]) + h
        append!(ys, [y1,newy,newy,y2])
    end
    return (reshape(xs, 4, size(hc.merge, 1)), reshape(ys, 4, size(hc.merge, 1)))
end

"""
Belongs to plot recipe for cluster dendrogram, from
https://gist.github.com/kescobo/ -> hclustrecipe.jl
"""
function treepositions(hc::Hclust, useheight::Bool)
    order = indexmap(hc.order)
    positions = Dict{}()
    for (k,v) in order
        positions[-k] = (v, 0)
    end
    for i in 1:size(hc.merge,1)
        xpos = mean([positions[hc.merge[i,1]][1], positions[hc.merge[i,2]][1]])
        if hc.merge[i,1] < 0 && hc.merge[i,2] < 0
            useheight ? ypos = hc.height[i] : ypos = 1
        else
            useheight ? h = hc.height[i] : h = 1
            ypos = maximum([positions[hc.merge[i,1]][2], positions[hc.merge[i,2]][2]]) + h
        end

        positions[i] = (xpos, ypos)
    end
    return positions
end

function drawSVG(figstem::String; edit::Bool=false)

    if (stat("$figstem.svg").size == 0) #generate new figure
        cp(joinpath(Pkg.dir("DH32Utils","data"),"template-drawing.svg"),"$figstem.svg")
        run(`inkscape $figstem.svg --with-gui -D -l $figstem.svg`)
    elseif (edit == true) #edit existing figure
        run(`inkscape $figstem.svg --with-gui -D -l $figstem.svg`)
    end

    #show figure
    graph=open(readall,"$figstem.svg")
    display("image/svg+xml",graph)

end

function drawSVGandPNG(figstem::String; edit::Bool=false)

    if (stat("$figstem.svg").size == 0) #generate new figure
        cp(joinpath(Pkg.dir("DH32Utils","data"),"template-drawing.svg"),"$figstem.svg")
        run(`inkscape $figstem.svg --with-gui -D -l $figstem.svg`)
    elseif (edit == true) #edit existing figure
        run(`inkscape $figstem.svg --with-gui -D -l $figstem.svg -e $figstem.png`)
    end

    run(`inkscape $figstem.svg -e $figstem.png`)

    #show figure
    load("$figstem.png")

end

function screenareaPNG(figstem::String; edit::Bool=false)

    if (stat("$figstem.png").size == 0 || edit == true)
        run(`scrot "$figstem.png" -s`) #generate new screenshot
    end

    #show image
    load("$figstem.png")

end

"""
Produces a smooth fit (Loess) for a set of x and y values.

Input:

    - x : input x values

    - y : input y values

    - nsteps (default=100) : number of steps between min. and max. of x

    - normalize (default: true) [see Loess package]

    - span (default: 0.75) [see Loess package]

    - degree (default: 2) [see Loess package]

Output:

    - x steps

    - fitted y values at stepped x positions
"""
function loess_fit{T <: AbstractFloat}(
    x::AbstractVector{T},
    y::AbstractVector{T};
    nsteps::Int64=100,
    normalize::Bool=true, span::T=0.75, degree::Int=2
)
    if (nsteps < 2)
        error("loess_fit: need at least 2 steps")
    end
    if (length(x) != length(y))
        error("loess_fit: x and y input arrays must have same lengths")
    end
    maxx = maximum(x)
    minx = minimum(x)
    if (minx >= maxx)
        error("loess_fit: maximum x value not greater than minimum x value")
    end
    delta = (maxx-minx)/nsteps
    xsteps = collect(minx:delta:maxx)

    model = loess(x, y; normalize=normalize, span=span, degree=degree)
    ymodel = Loess.predict(model, xsteps)

    xsteps, ymodel
end

"""
Computes block averages of a float vector and determines quantiles (default 0.05 and 0.95 quantiles) of the distribution of these averages for each block size.

The vector elements can be weighted with an additional vector of the same length. In this case the weights in each block are normalized to sum up to one.

Input:

- x: Float64 array

- optionally the following named parameters can be given:

    - dimin (default 10): the minimum block interval (= min. block size)

    - dimax (default 0, will be automatically calculated): the maximum block interval (= max. block size)

    - ddi (default 10): increment in block size

    - w (default [1.] = no weights used): vector of weights, if not given will be set to [1., 1., ...., 1.]

    - qmin (default 0.05): lower quantile

    - qmax (defaults 0.95): upper quantile

Output:

    - qmins: array of lower quantile of block averages

    - qmaxs: array of upper quantile of block averages
"""
function block_average_spread(
            x::Array{Float64,1}; #quantity to be averaged
            dimin::Int64=10, #minimum width of blocks
            dimax::Int64=0, #maximum width of blocks
            ddi::Int64=10, #block increment
            w::Array{Float64,1}=[1.], #weights
            qmin::Float64=0.05, #minimum quantile to be reported
            qmax::Float64=0.95 #maximum quantile to be reported
        )
    nx = length(x)
    if w != [1.] #weighted
        if length(w) != nx
            error("error block_averages: length(w)!=length(x)")
        end
        if sum(w .> 0.) != nx
            error("error block_averages: some weights not positive")
        end
    else
        w = ones(nx)
    end
    if dimax == 0 #the last block average is between first and second half of data
        dimax = div(nx,2)
    end
    if dimax < dimin
        error("error block_averages: dimax < dimin")
    end

    ba = zeros(div(nx,dimin)) #array for analysis of block averages

    nout = length(collect(dimin:ddi:dimax))
    qmins = zeros(nout) #output of lower quantiles of block averages
    qmaxs = zeros(nout) #output of upper quantiles of bl. av.

    jout=1
    for di in dimin:ddi:dimax #width of blocks
        ni = div(nx,di) # number of blocks
        bstart = 1 #index of block start
        bend = di #index of block end
        for i in 1:ni
            ba[i]=dot(x[bstart:bend],w[bstart:bend])/sum(w[bstart:bend])
            bstart += di
            bend += di
        end
        qmins[jout] = quantile(ba[1:ni],qmin)
        qmaxs[jout] = quantile(ba[1:ni],qmax)
        jout += 1
    end
    qmins, qmaxs
end

"""

Compute midpoints of given range or vector.

"""
function midpoints(r::Range)
    r[1:length(r)-1] + 0.5*step(r)
end

function midpoints(v::AbstractVector)
    [0.5*(v[i] + v[i+1]) for i in 1:length(v)-1]
end

"""
Example:

    For given file "bla.jl" script the call

    save_script_and_fig_to_unique_filenames("bla")

will generate a copy of that file in "bla_RsTzzu.jl" and save the last figure to "bla_RsTzzu.pdf".

Input:

    - source_stem: stem name of source file, string; this should be the name of the existing script file (without extension_script) that will be written in first output

    - extension_script (optional, default ".jl"): ending of script file

    - extension_fig (optional, default ".pdf"): ending and type of figure

    - separator (optional, defaul "_"): string that separates source stem and unique component of file names

Output:

    - writes script file with name source*unique*extension_script in current directory

    - saves last figure to file named source*unique*extension_fig in current directory

"unique" in both outputs is the same unique filename string.

"""

function save_script_and_fig_to_unique_filenames(source_stem::String="";
                                                 extension_script::String=".jl",
                                                 extension_fig::String=".pdf",
                                                 separator::String="_")

    if extension_script == extension_fig
        error("extension_script and extension_fig must be different")
    end

    if source_stem==""
        error("source file stem name required as input")
    end

    if !isfile(source_stem*extension_script)
        error("source file "*source_stem*extension_script*" not present")
    end

    written=false
    trials=0
    while (!written && trials<1000)
        #propose unique file names
        proposed_unique=randstring(10)
        outfile_script = source_stem*separator*proposed_unique*extension_script
        outfile_fig = source_stem*separator*proposed_unique*extension_fig

        if !isfile(outfile_script) && !isfile(outfile_fig)
            #no files with the proposed names present
            cp(source_stem*extension_script, outfile_script)
            savefig(outfile_fig)
            written=true
            break
        end
        trials += 1
    end
end

end # module