text:
export find_library, find_binary,
       find_cuda_library, find_cuda_binary,
       find_toolkit, find_toolkit_version,
       find_libdevice, find_libcudadevrt,
       find_host_compiler, find_toolchain

function resolve(path)
    if islink(path)
        dir = dirname(path)
        resolve(joinpath(dir, readlink(path)))
    else
        path
    end
end

# return a list of valid directories, resolving symlinks and pruning duplicates
function valid_dirs(dirs)
    map!(resolve, dirs, dirs)
    filter(isdir, unique(dirs))
end

const gcc_major_versions = 3:9
const gcc_minor_versions = 0:9

## generic discovery routines

"""
    find_library(names; locations=String[], versions=VersionNumber[], word_size=Sys.WORD_SIZE)

Wrapper for Libdl.find_library, performing a more exhaustive search:

- variants of the library name (including version numbers, platform-specific tags, etc);
- various subdirectories of the `locations` list, and finally system library directories.

Returns the full path to the library.
"""
function find_library(names::Vector{String};
                      locations::Vector{String}=String[],
                      versions::Vector{VersionNumber}=VersionNumber[],
                      word_size::Integer=Sys.WORD_SIZE)
    @trace "Request to look for library $(join(names, ", "))" locations

    # figure out names
    all_names = String[]
    if Sys.iswindows()
        # priority goes to the `names` argument, as per `Libdl.find_library`
        for name in names
            # first look for unversioned libraries, for upgrade resilience
            append!(all_names, ["$(name)$(word_size)", name])
            for version in versions
                append!(all_names, ["$(name)$(word_size)_$(version.major)",
                                    "$(name)$(word_size)_$(version.major)$(version.minor)"])
            end
        end
    elseif Sys.isunix()
        # most UNIX distributions ship versioned libraries (also see JuliaLang/julia#22828)
        for name in names
            # first look for unversioned libraries, for upgrade resilience
            push!(all_names, "lib$(name).$(Libdl.dlext)")
            for version in versions
                append!(all_names, ["lib$(name).$(Libdl.dlext).$(version.major)",
                                    "lib$(name).$(Libdl.dlext).$(version.major).$(version.minor)"])
            end
        end
    else
        # let Libdl do all the work
        all_names = ["lib$name" for name in names]
    end
    unique!(all_names)

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        if Sys.iswindows()
            push!(all_locations, joinpath(location, "bin"))
        else
            push!(all_locations, joinpath(location, "lib"))
            if word_size == 64
                push!(all_locations, joinpath(location, "lib64"))
            end
        end
    end

    @trace "Looking for library $(join(all_names, ", "))" locations=all_locations
    name_found = Libdl.find_library(all_names, all_locations)
    if isempty(name_found)
        return nothing
    end

    # find the full path of the library (which Libdl.find_library doesn't guarantee to return)
    path = Libdl.dlpath(name_found)
    @debug "Found library $(basename(path)) at $(dirname(path))"
    return path
end

"""
    find_binary(names; locations=String[])

Similar to `find_library`, performs an exhaustive search for a binary in various
subdirectories of `locations`, and finally PATH.
"""
function find_binary(names::Vector{String};
                     locations::Vector{String}=String[])
    @trace "Request to look for binary $(join(names, ", "))" locations

    # figure out names
    all_names = String[]
    if Sys.iswindows()
        all_names = ["$name.exe" for name in names]
    else
        all_names = names
    end

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        push!(all_locations, joinpath(location, "bin"))
    end
    let path = ENV["PATH"]
        dirs = split(path, Sys.iswindows() ? ';' : ':')
        filter!(path->!isempty(path), dirs)
        append!(all_locations, dirs)
    end

    @trace "Looking for binary $(join(all_names, ", "))" locations=all_locations
    all_paths = [joinpath(location, name) for name in all_names, location in all_locations]
    paths = String[]
    for path in all_paths
        try
            if ispath(path)
                push!(paths, path)
            end
        catch
            # some system disallow `stat` on certain paths
        end
    end

    if isempty(paths)
        return nothing
    else
        path = first(paths)
        @debug "Found binary $(basename(path)) at $(dirname(path))"
        return path
    end
end

## CUDA-specific discovery routines

# FIXME: CUDA on 32-bit Windows isn't supported

const cuda_names = Dict(
    "cuda"      => Sys.iswindows() ? ["nvcuda"] : ["cuda"],
    "nvml"      => Sys.iswindows() ? ["nvml"]   : ["nvidia-ml"],
    "nvtx"      => ["nvToolsExt"]
)

const cuda_versions = Dict(
    # https://developer.nvidia.com/cuda-toolkit-archive
    "toolkit"   => [v"1.0", v"1.1",
                    v"2.0", v"2.1", v"2.2",
                    v"3.0", v"3.1", v"3.2",
                    v"4.0", v"4.1", v"4.2",
                    v"5.0", v"5.5",
                    v"6.0", v"6.5",
                    v"7.0", v"7.5",
                    v"8.0",
                    v"9.0", v"9.1", v"9.2",
                    v"10.0", v"10.1"],
    # https://developer.nvidia.com/rdp/cudnn-archive
    "cudnn"     => [v"1.0",
                    v"2.0",
                    v"3.0",
                    v"4.0",
                    v"5.0", v"5.1",
                    v"6.0",
                    v"7.0", v"7.1", v"7.2", v"7.3", v"7.4"],
    "nvtx"      => [v"1.0.0"]
)

# simplified find_library/find_binary entry-points,
# looking up name aliases and known version numbers
# and passing the (optional) toolkit dirs as locations.
find_cuda_library(name::String, toolkit_dirs::Vector{String}=String[];
                  versions::Vector{VersionNumber}=reverse(get(cuda_versions, name, cuda_versions["toolkit"])),
                  kwargs...) =
    find_library(get(cuda_names, name, [name]);
                 versions=versions,
                 locations=toolkit_dirs,
                 kwargs...)
find_cuda_binary(name::String, toolkit_dirs::Vector{String}=String[]; kwargs...) =
    find_binary(get(cuda_names, name, [name]);
                locations=toolkit_dirs,
                kwargs...)

"""
    find_toolkit()::Vector{String}

Look for directories where (parts of) the CUDA toolkit might be installed. This returns a
(possibly empty) list of paths that can be used as an argument to other discovery functions.

The behavior of this function can be overridden by defining the `CUDA_PATH`, `CUDA_HOME` or
`CUDA_ROOT` environment variables, which should point to the root of the CUDA toolkit.
"""
function find_toolkit()
    dirs = String[]

    # NVTX library (special case for Windows)
    if Sys.iswindows()
        var = "NVTOOLSEXT_PATH"
        basedir = get(ENV, var, nothing)
        if basedir !== nothing && isdir(basedir)
            @trace "Looking for NVTX library via environment variable" var
            suffix = Sys.WORD_SIZE == 64 ? "x64" : "Win32"
            dir = joinpath(basedir, "bin", suffix)
            isdir(dir) && push!(dirs, dir)
        else
            program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
            basedir = joinpath(program_files, "NVIDIA Corporation", "NvToolsExt")
            @trace "Looking for NVTX library in the default directory" basedir
            suffix = Sys.WORD_SIZE == 64 ? "x64" : "Win32"
            dir = joinpath(basedir, "bin", suffix)
            isdir(dir) && push!(dirs, dir)
        end
    end

    # look for environment variables to override discovery
    envvars = ["CUDA_PATH", "CUDA_HOME", "CUDA_ROOT"]
    envdict = Dict(Symbol(var) => ENV[var] for var in envvars if haskey(ENV, var))
    if length(envdict) > 0
        if length(unique(values(envdict))) > 1
            @warn "Multiple CUDA environment variables set to different values" envdict...
        end

        @trace "Looking for CUDA toolkit via environment variables" envdict...
        append!(dirs, collect(values(envdict)))
        return dirs
    end

    # look for the compiler binary (in the case PATH points to the installation)
    nvcc_path = find_cuda_binary("nvcc")
    if nvcc_path !== nothing
        nvcc_dir = dirname(nvcc_path)
        if occursin(r"^bin(32|64)?$", basename(nvcc_dir))
            nvcc_dir = dirname(nvcc_dir)
        end

        @trace "Looking for CUDA toolkit via nvcc binary" path=nvcc_path dir=nvcc_dir
        push!(dirs, nvcc_dir)
    end

    # look for the runtime library (in the case LD_LIBRARY_PATH points to the installation)
    libcudart_path = find_cuda_library("cudart")
    if libcudart_path !== nothing
        libcudart_dir = dirname(libcudart_path)
        if occursin(r"^(lib|bin)(32|64)?$", basename(libcudart_dir))
            libcudart_dir = dirname(libcudart_dir)
        end

        @trace "Looking for CUDA toolkit via CUDA runtime library" path=libcudart_path dir=libcudart_dir
        push!(dirs, libcudart_dir)
    end

    # look in default installation directories
    default_dirs = String[]
    if Sys.iswindows()
        # CUDA versions are installed in separate directories under a single base dir
        program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
        basedir = joinpath(program_files, "NVIDIA GPU Computing Toolkit", "CUDA")
        if isdir(basedir)
            entries = map(dir -> joinpath(basedir, dir), readdir(basedir))
            reverse!(entries) # we want to search starting from the newest CUDA version
            append!(default_dirs, entries)
        end
    else
        # CUDA versions are installed in unversioned dirs, or suffixed with the version
        basedirs = ["/usr/local/cuda", "/opt/cuda"]
        for dir in basedirs
            append!(default_dirs, "$dir-$(ver.major).$(ver.minor)" for ver in cuda_versions["toolkit"])
        end
        append!(default_dirs, basedirs)
        push!(default_dirs, "/usr/lib/nvidia-cuda-toolkit")
        push!(default_dirs, "/usr/share/cuda")
    end
    default_dirs = valid_dirs(default_dirs)
    if !isempty(default_dirs)
        @trace "Looking for CUDA toolkit via default installation directories" dirs=default_dirs
        append!(dirs, default_dirs)
    end

    # filter
    dirs = valid_dirs(dirs)
    @debug "Found CUDA toolkit at $(join(dirs, ", "))"
    return dirs
end

# figure out the CUDA toolkit version (by looking at the `nvcc --version` output)
function find_toolkit_version(toolkit_dirs)
    nvcc_path = find_cuda_binary("nvcc", toolkit_dirs)
    if nvcc_path === nothing
        error("CUDA toolkit at $(join(toolkit_dirs, ", ")) doesn't contain nvcc")
    end

    # parse the nvcc version string
    verstr = withenv("LANG"=>"C") do
        read(`$nvcc_path --version`, String)
    end
    m = match(r"\bV(?<major>\d+).(?<minor>\d+).(?<patch>\d+)\b", verstr)
    m !== nothing || error("could not parse NVCC version info (\"$verstr\")")

    version = VersionNumber(parse(Int, m[:major]),
                            parse(Int, m[:minor]),
                            parse(Int, m[:patch]))
    @debug "CUDA toolkit identified as $version"
    return version
end

"""
    find_libdevice(targets::Vector{VersionNumber}, toolkit_dirs::Vector{String})

Look for the CUDA device library supporting `targets` in any of the CUDA toolkit directories
`toolkit_dirs`. On CUDA >= 9.0, a single library unified library is discovered and returned
as a string. On older toolkits, individual libraries for each of the targets are returned as
a vector of strings.
"""
function find_libdevice(targets::Vector{VersionNumber}, toolkit_dirs)
    @trace "Request to look for libdevice $(join(targets, ", "))" locations=toolkit_dirs

    # figure out locations
    dirs = String[]
    for toolkit_dir in toolkit_dirs
        push!(dirs, toolkit_dir)
        push!(dirs, joinpath(toolkit_dir, "libdevice"))
        push!(dirs, joinpath(toolkit_dir, "nvvm", "libdevice"))
    end

    # filter
    dirs = valid_dirs(dirs)
    @trace "Look for libdevice $(join(targets, ", "))" locations=dirs

    for dir in dirs
        # parse filenames
        libraries = Dict{VersionNumber,String}()
        for target in targets
            path = joinpath(dir, "libdevice.compute_$(target.major)$(target.minor).10.bc")
            if isfile(path)
                libraries[target] = path
            end
        end
        library = nothing
        let path = joinpath(dir, "libdevice.10.bc")
            if isfile(path)
                library = path
            end
        end

        # select
        if library !== nothing
            @debug "Found unified device library at $library"
            return library
        elseif !isempty(libraries)
            @debug "Found split device libraries at $(join(libraries, ", "))"
            return libraries
        end
    end

    return nothing
end

"""
    find_libcudadevrt(toolkit_dirs::Vector{String})

Look for the CUDA device runtime library in any of the CUDA toolkit directories
`toolkit_dirs`.
"""
function find_libcudadevrt(toolkit_dirs)
    locations = toolkit_dirs
    @trace "Request to look for libcudadevrt " locations

    name = nothing
    if Sys.isunix()
        name = "libcudadevrt.a"
    elseif Sys.iswindows()
        name = "cudadevrt.lib"
    else
        error("No support for discovering the CUDA device runtime library on your platform, please file an issue.")
    end

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        if Sys.iswindows()
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib", "x64"))
            elseif Sys.WORD_SIZE == 32
                push!(all_locations, joinpath(location, "lib", "Win32"))
            end
        else
            push!(all_locations, joinpath(location, "lib"))
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib64"))
            end
        end
    end

    @trace "Looking for CUDA device runtime library $name" locations=all_locations
    paths = filter(isfile, map(location->joinpath(location, name), all_locations))

    if isempty(paths)
        return nothing
    else
        path = first(paths)
        @debug "Found CUDA device runtime library $(basename(path)) at $(dirname(path))"
        return path
    end
end

function find_host_compiler(toolkit_version=nothing)
    if !(Sys.iswindows() || Sys.isapple())
        # Unix-like platforms: find compatible GCC binary
        # enumerate possible names for the gcc binary
        # NOTE: this is coarse, and might list invalid, non-existing versions
        gcc_names = [ "gcc", "cuda-gcc" ]
        for major in gcc_major_versions
            push!(gcc_names, "gcc-$major")
            for minor in gcc_minor_versions
                push!(gcc_names, "gcc-$major.$minor")
                push!(gcc_names, "gcc$major$minor")
            end
        end

        # find the binary
        gcc_possibilities = []
        for gcc_name in gcc_names
            # check if the binary exists
            gcc_path = find_binary([gcc_name])
            if gcc_path === nothing
                continue
            end

            # parse the GCC version string
            verstr = withenv("LANG" => "C") do
                readlines(`$gcc_path --version`)[1]
            end
            m = match(Regex("^$(basename(gcc_path)) \\(.*\\) ([0-9.]+)"), verstr)
            if m === nothing
                @warn "Could not parse GCC version info (\"$verstr\"), skipping this compiler."
                continue
            end
            gcc_ver = VersionNumber(m.captures[1])
            @trace "Found GCC $gcc_ver at $gcc_path"

            if toolkit_version === nothing || gcc_supported(gcc_ver, toolkit_version)
                push!(gcc_possibilities, (gcc_path, gcc_ver))
            elseif toolkit_version !== nothing
                @warn "Ignoring $gcc_path v$gcc_ver which isn't supported by CUDA $toolkit_version"
            end
        end

        # select the most recent compiler
        if length(gcc_possibilities) == 0
            error("Could not find a suitable GCC")
        end
        sort!(gcc_possibilities; rev=true, lt=(a, b) -> a[2]<b[2])
        host_compiler, host_version = gcc_possibilities[1]
    elseif Sys.iswindows()
        # Windows: search for MSVC in default locations

        # discover Visual Studio installations
        msvc_paths = String[]
        program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles(x86)" : "ProgramFiles"]
        ## locate ≥ VS2017
        vswhere_dist = joinpath(program_files, "Microsoft Visual Studio", "Installer", "vswhere.exe")
        vswhere_url = "https://github.com/Microsoft/vswhere/releases/download/2.6.7/vswhere.exe"
        let vswhere = Sys.which("vswhere.exe") !== nothing ? "vswhere.exe" : isfile(vswhere_dist) ? vswhere_dist : download(vswhere_url)
            msvc_cmd_tools_dir = readchomp(`$vswhere -latest -products \* -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`)
            if !isempty(msvc_cmd_tools_dir)
                msvc_build_ver = readchomp(joinpath(msvc_cmd_tools_dir, "VC\\Auxiliary\\Build\\Microsoft.VCToolsVersion.default.txt"))
                arch = Sys.WORD_SIZE == 64 ? "x64" : "x86"
                msvc_path = joinpath(msvc_cmd_tools_dir,"VC\\Tools\\MSVC\\$msvc_build_ver\\bin\\Host$arch\\$arch\\cl.exe")
                @trace "Considering MSVC at $msvc_path located with vswhere"
                push!(msvc_paths, msvc_path)
            end
        end
        ## locate VS2012 to 2014
        let envvars = ["VS140COMNTOOLS", "VS120COMNTOOLS", "VS110COMNTOOLS", "VS100COMNTOOLS"]
            envvars_set = filter(var -> haskey(ENV, var), envvars)
            for var in envvars_set
                val = ENV[var]
                arch = Sys.WORD_SIZE == 64 ? "amd64" : "x86"
                msvc_path = joinpath(dirname(dirname(dirname(val))), "VC", "bin", arch, "cl.exe")
                if isfile(msvc_path)
                    @trace "Considering MSVC at $msvc_path located with environment variable $var"
                    push!(msvc_paths, msvc_path)
                end
            end
        end
        ## look in PATH as well
        let msvc_path = Sys.which("cl.exe")
            if msvc_path !== nothing
                push!(msvc_paths, msvc_path)
            end
        end
        isempty(msvc_paths) && error("No Visual Studio installation found")

        # find MSVC versions
        msvc_list = Dict{VersionNumber,String}()
        for path in msvc_paths
            tmpfile = tempname() # TODO: do this with a pipe
            if !success(pipeline(`$path`, stdout=devnull, stderr=tmpfile))
                @warn "Could not execute $path"
                continue
            end
            verstr = read(tmpfile, String)
            m = match(r"Version\s+(\d+(\.\d+)?(\.\d+)?)"i, verstr)
            if m === nothing
                # depending on the locale, this regex might not match
                m = match(r"\b(\d+(\.\d+)?(\.\d+)?)\b"i, verstr)
            end
            if m === nothing
                @warn "Could not parse Visual Studio version info (\"$verstr\"), skipping this compiler."
                continue
            end
            msvc_ver = VersionNumber(m.captures[1])
            msvc_list[msvc_ver] = path
        end

        # check compiler compatibility
        msvc_path, msvc_ver = nothing, nothing
        if toolkit_version !== nothing
            for ver in sort(collect(keys(msvc_list)), rev=true) # search the highest version first
                if msvc_supported(ver, toolkit_version)
                    msvc_path, msvc_ver = msvc_list[ver], ver
                    break
                else
                    @warn "Ignoring $msvc_path v$msvc_ver which isn't supported by CUDA $toolkit_version"
                end
            end
            if msvc_ver === nothing
                error("None of the available Visual Studio C++ compilers ($(join(keys(msvc_list), ", "))) are compatible with CUDA $toolkit_version")
            end
        else
            # take the first found host, which will be the highest version found
            host_pair = first(sort(collect(msvc_list), by=x->x[1], rev=true))
            msvc_path, msvc_ver = last(host_pair), first(host_pair)
        end

        host_compiler, host_version = msvc_path, msvc_ver
    elseif Sys.isapple()
        # GCC is no longer supported on MacOS so let's just use clang
        # TODO: discovery of all compilers, and version matching against the toolkit
        clang_path = find_binary(["clang"])
        if clang_path === nothing
            error("Could not find clang")
        end
        verstr = read(`$clang_path --version`, String)
        m = match(r"version\s+(\d+(\.\d+)?(\.\d+)?)"i, verstr)
        if m === nothing
            error("Could not parse Clang version info (\"$verstr\")")
        end
        clang_ver = VersionNumber(m[1])

        host_compiler, host_version = clang_path, clang_ver
    end

    @debug "Selected host compiler version $host_version at $host_compiler"
    return host_compiler, host_version
end

mutable struct Toolchain
    cuda_compiler::String
    cuda_version::VersionNumber

    host_compiler::String
    host_version::VersionNumber
end
function find_toolchain(toolkit_dirs, toolkit_version=find_toolkit_version(toolkit_dirs))
    # find the CUDA compiler
    nvcc_path = find_cuda_binary("nvcc", toolkit_dirs)
    if nvcc_path === nothing
        error("CUDA toolkit at $(join(toolkit_dirs, ", ")) doesn't contain nvcc")
    end
    nvcc_version = toolkit_version

    # find a suitable host compiler
    host_compiler, host_version = find_host_compiler(toolkit_version)

    return Toolchain(nvcc_path, nvcc_version,
                     host_compiler, host_version)
end


text:
module SimpleLogger

using Printf
using POLO: AbstractLogger
using POLO.Utility: LogLevel, None, Value, Decision, Gradient, Full
import POLO.log

struct Logger{L <: LogLevel} <: AbstractLogger
    function (::Type{Logger})()
        return new{None}()
    end

    function (::Type{Logger})(::Type{L}) where L <: LogLevel
        return new{L}()
    end
end
None() = Logger()
Value() = Logger(Value)
Decision() = Logger(Decision)
Gradient() = Logger(Gradient)
Full() = Logger(Full)

function log(logger::Logger{None},k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    nothing
end

function log(logger::Logger{Value},k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    if k == 1
        println("k    fval")
    end
    println(@sprintf("%d    %.2f",k,fval))
end

function log(logger::Logger{Decision},k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    if k == 1
        println("k    fval    x")
    end
    print(@sprintf("%d    %.2f   ",k,fval))
    println(x')
end

function log(logger::Logger{Gradient},k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    if k == 1
        println("k    fval    g")
    end
    print(@sprintf("%d    %.2f   ",k,fval))
    println(g')
end

end

ext: jl


text:
"""
This is an example of how to extend CloudScat with new geometrical shapes.
Here we will build a cloud geometry consisting in concentric rings.  This
illustrates both the required extensions to define new geometrical shapes
and also how to include shapes defined implicitly by a mathematical function
(this is S = {x: f(x) < 0}).  This involves repeatedly solving the equation
f(x) = 0 along observation paths and therefore one should expect longer
simulation times.

To use this code run julia and type

julia> include("rings.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia rings.jl

- - -
2020 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat
using Parameters, Roots
const co = CloudScat.constants

# We need to import these names besides those exported by CloudScat
import CloudScat: interfaces!, inside, accum!, Accumulator

"""
Geometrical shape that encapsulates concentric rings of period `L`, amplitude
`A` and base altitude `z0`.  The top of the cloud is z = z0 + A * cos(2πρ/L),
with ρ = √(x²+y²).
"""
struct ConcentricRings <: Shape
    "Period (m)"
    L::Float64

    "Amplitude (m)"
    A::Float64

    "Base altitude (m)"
    z0::Float64
end

# This function provides a way to find intersections between the shape and
# straight lines between two points a and b.  We must find values of s such that
# (1 - s)a + sb is an interface point (0 < s < 1).  For each of these points
# we must call accum!(v, s).  Here we use the find_zeros function from the julia
# Roots package.  Note that in some cases some roots may not be found so use this
# with care.
function interfaces!(v::Accumulator, rings::ConcentricRings, a::Point, b::Point)
    function f(s)
        @unpack A, z0, L = rings
        p = @. (1 - s) * a + s * b
        return p[3] - (A * cos(2π * sqrt(p[1]^2 + p[2]^2) / L) + z0)
    end

    for s in find_zeros(f, 0.0, 1.0)
        accum!(v, s)
    end
end

# This function just checks whether a point is inside the shape.
function inside(rings::ConcentricRings, p::Point)
    @unpack A, z0, L = rings
    p[3] < (A * cos(2π * sqrt(p[1]^2 + p[2]^2) / L) + z0)
end

#
# Once a new shape is defined (here, ConcentricRings) one can use it and combine
# with other shapes to define a cloud geometry.
#

function run()
    params = Params(
        # Number of simulated photons
        N = 400000,

        # The wavelength of the scattered light
        λ = 337 * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 10 * co.kilo],
        source_b = [0, 0, 10 * co.kilo],

        # Minimum altitude for observations
        zmin = 7 * co.kilo)

    # A cylinder is defined as Cylinder(zbottom, ztop, xcenter, ycenter, radius)
    c1 = Cylinder(7 * co.kilo,  15 * co.kilo, 0, 0, 20 * co.kilo)

    rings = ConcentricRings(2 * co.kilo, 2 * co.kilo, 10 * co.kilo)

    # We limit the extension of the concentric rings by intersecting the shape
    # with a 20-km-radius cylinder.
    cloud = intersect(c1, rings)

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Cylinder(7 * co.kilo, 60 * co.kilo, 0, 0, 200 * co.kilo)

    # Set a homogeneous cloud droplet density and radius.
    # See variable_droplet_radius.jl for how to deal with inhomogeneous clouds.
    # Radius of the scattering particles
    radius = 10e-6

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, radius)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    # Observers record a time curve and an image of received light.  You can
    # set up as many observers as you wish but they have a significant impact
    # on performance.

    observers = [Observer(
        # Location of the observer
        position = [0, 0, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                   saveto="$(basename).h5")
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()


text:
using StaticThreadSchedules

v = zeros(12)
# suppose this operation scales like log(i)
@threads :static (i->i^2) for i in 1:12
    v[i] = Threads.threadid()
end
print(v)


text:
using ModelTES
using Test

Tc=0.107; Tbath = 0.065
R0 = 1.55e-3; Rl = 0.35e-3; Rn = 10.3e-3; Rpara=0.0
n=3.3; k = 3.50e-9; C = 0.72e-12
alpha = 92.0;β = 1.4; L=300e-9
model1 = ModelTES.ShankRIT(alpha, β, n, Tc, Tbath, k, R0, Rn);
tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model1)
bt = BiasedTES(tes_param, R0)
ihtes = IrwinHiltonTES(bt)

@test Tbath == ihtes.Tbath
@test isapprox(R0, ihtes.R0; rtol=1e-5)
@test Rl == ihtes.Rl
@test C == ihtes.C0
@test isapprox(alpha, ihtes.alpha; rtol=1e-5)
@test isapprox(β, ihtes.beta; rtol=1e-5)
@test L == ihtes.L

f = 10.0 .^ range(0, stop=6, length=100);
n,n1,n2,n3,n4 = ModelTES.noisePSD(ihtes, f);
z = ModelTES.Z(ihtes,f)
zcircuit = ModelTES.Zcircuit(ihtes,f)

zr, zi = real.(z), imag.(z)

# Check issue #4: is the normalization correct on PSD and autocorrelation?
# Estimate the total power by a weird log-space trapezoid rule.
amplifier_noise=1e-35
frequencies = collect(10 .^ range(0, stop=6, length=1501));
psdRef = ModelTES.noisePSD(bt, frequencies, amplifier_noise)[1];
totalpower = psdRef[1]*frequencies[1];
for i=2:length(frequencies)
    global totalpower += (frequencies[i]-frequencies[i-1])*0.5*(psdRef[i]+psdRef[i-1]);
end

sampleTimeA = 1e-8
modelRef = NoiseModel(bt, sampleTimeA, amplifier_noise)
for sampleTimeB in 10 .^ range(-7.5, stop=-6, length=4)
    model = NoiseModel(bt, sampleTimeB, amplifier_noise)
    @test isapprox(model.covarIV[1], modelRef.covarIV[1], rtol=0.01)
    @test isapprox(model.covarIV[1], totalpower, rtol=0.2)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    pix_write_ps_i(
        pix::Pix,
        box::Union{PixBox, Nothing},
        ppi::Integer,
        scale::AbstractFloat
    )::Tuple{Ptr{UInt8}, Csize_t}

Write an image to memory as a PostScript file. If there is an error `(C_NULL, 0)` is returned.

__Arguments:__

| T | Name     | Default  | Description
|:--| :------- | :------- | :----------
| R | pix      |          | The image to write to memory.
| R | box      |          | Optional location to have the image appear on the page.
| R | ppi      |          | The resolution to use for the image.
| R | scale    |          | Scale the image on the page.

__Restrictions:__

  * `ppi` - Must be greater than 0.
  * `scale` - Must be greater than or equal to `0.0`.

__Details:__

**This method is meant for internal use only.**  It returns the memory pointer and number of bytes
allocated by the Leptonica library to write out the image.  The caller is responsible for freeing
this memory with a call to `lept_free()`.

If you don't want the image scaled set the `scale` to `1.0`.  If you want to position the image on
the page you can use a `box` however set `w` and `h` to `0`, then the image will be positioned
where you want it and not scaled.

If you want the image to be a specific size set the scale to `0.0` and use a box setting `w` and
`h` to the width and height you want in thousandth of an inch.  So a width of `5000` and a height of
`2000` will give you an image that is 5 inches across and 2 inches high.

If you want the image scaled set the `scale` to the factor you desire.  Value of `0.5` will cause
the image to be halfed, while a value of `2.0` will double the size of the image.  If you want to
also position the image on the page use a `box` but set the `w` and `h` values to `0`.
"""
function pix_write_ps_i(
            pix::Pix,
            box::Union{PixBox, Nothing},
            ppi::Integer,
            scale::AbstractFloat
        )::Tuple{Ptr{UInt8}, Csize_t}
    local temp   = Ref(Ptr{UInt8}(C_NULL))
    local size   = Ref(Csize_t(0))
    local boxPtr = C_NULL

    # ---------------------------------------------------------------------------------------------
    # Verify the parameters.
    if is_valid(pix) == false
        @error "Pix has been freed."
        return (C_NULL, 0)
    end

    if ppi <= 0
        @error "PPI must be positive."
        return (C_NULL, 0)
    end

    if scale <= Float32(0.0)
        @error "Scale must be positive."
        return (C_NULL, 0)
    end

    # ---------------------------------------------------------------------------------------------
    # If a box was provided allocate a box.
    if box != nothing
        boxPtr = ccall(
            (:boxCreate, LEPTONICA),
            Ptr{Cvoid},
            (Cint, Cint, Cint, Cint),
            box.x, box.y, box.w, box.h
        )
    end

    if box != nothing && boxPtr == C_NULL
        return (C_NULL, 0)
    end

    # ---------------------------------------------------------------------------------------------
    # Call the write function.
    local strPtr = ccall(
        (:pixWriteStringPS, LEPTONICA),
        Ptr{UInt8},
        (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cfloat),
        pix,
        boxPtr,
        ppi,
        scale
    )

    # ---------------------------------------------------------------------------------------------
    # Release the box if used.
    if boxPtr != C_NULL
        ccall(
            (:boxDestroy, LEPTONICA),
            Cvoid,
            (Ptr{Cvoid},),
            Ref(boxPtr)
        )
    end

    if strPtr == C_NULL
        return (C_NULL, 0)
    end

    # ---------------------------------------------------------------------------------------------
    # Figure out how long the string is.
    local size = 0
    while unsafe_load(strPtr, size+1) != 0x00
        size += 1
    end

    return (strPtr, size)
end

# =================================================================================================
"""
    pix_write_ps(
        filename::AbstractString,
        pix::Pix;
        box::Union{PixBox, Nothing} = nothing,
        ppi::Integer = Int32(300),
        scale::AbstractFloat = Float32(1.0)
    )::Bool

Write an image to a PostScript file.  If there is an error `false` is returned.

__Arguments:__

| T | Name     | Default        | Description
|:--| :------- | :------------- | :----------
| R | filename |                | The name of the file to write to.
| R | pix      |                | The image to write to the file.
| O | box      | `nothing`      | Location to have the image appear on the page.
| O | ppi      | `Int32(300)`   | The resolution to use for the image.
| O | scale    | `Float32(1.0)` | Scale the image on the page.

__Restrictions:__

  * `ppi` - Must be greater than 0.
  * `scale` - Must be greater than or equal to `0.0`.

__Details:__

If the file exists it will be overwritten.

If you don't want the image scaled set the `scale` to `1.0`.  If you want to position the image on
the page you can use a `box` however set `w` and `h` to `0`, then the image will be positioned
where you want it and not scaled.

If you want the image to be a specific size set the scale to `0.0` and use a box setting `w` and
`h` to the width and height you want in thousandth of an inch.  So a width of `5000` and a height of
`2000` will give you an image that is 5 inches across and 2 inches high.

If you want the image scaled set the `scale` to the factor you desire.  Value of `0.5` will cause
the image to be halfed, while a value of `2.0` will double the size of the image.  If you want to
also position the image on the page use a `box` but set the `w` and `h` values to `0`.
"""
function pix_write_ps(
            filename,
            pix::Pix;
            box::Union{PixBox, Nothing} = nothing,
            ppi::Integer                = Int32(300),
            scale::AbstractFloat        = Float32(1.0)
        )::Bool
    local file   = nothing
    local result = false
    local data, size = pix_write_ps_i(pix, box, ppi, scale)

    if data != C_NULL
        try
            if size > 0
                file = open(filename; write=true, create=true)
                unsafe_write(file, data, size)
                result = true
            end
        finally
            lept_free(data)

            if file != nothing
                close(file)
            end
        end
    end

    return result
end

# =================================================================================================
"""
    pix_write_ps(
        stream::IO,
        pix::Pix;
        box::Union{PixBox, Nothing} = nothing,
        ppi::Integer = Int32(300),
        scale::AbstractFloat = Float32(1.0)
    )::Bool

Write an image to a PostScript file.  Returns 'false' if there was an error.

__Arguments:__

| T | Name     | Default        | Description
|:--| :------- | :------------- | :----------
| R | stream   |                | The stream to write the image to.
| R | pix      |                | The image to write to the stream.
| O | box      | `nothing`      | Location to have the image appear on the page.
| O | ppi      | `Int32(300)`   | The resolution to use for the image.
| O | scale    | `Float32(1.0)` | Scale the image on the page.

__Restrictions:__

  * `ppi` - Must be greater than 0.
  * `scale` - Must be greater than or equal to `0.0`.

__Details:__

If you don't want the image scaled set the `scale` to `1.0`.  If you want to position the image on
the page you can use a `box` however set `w` and `h` to `0`, then the image will be positioned
where you want it and not scaled.

If you want the image to be a specific size set the scale to `0.0` and use a box setting `w` and
`h` to the width and height you want in thousandth of an inch.  So a width of `5000` and a height of
`2000` will give you an image that is 5 inches across and 2 inches high.

If you want the image scaled set the `scale` to the factor you desire.  Value of `0.5` will cause
the image to be halfed, while a value of `2.0` will double the size of the image.  If you want to
also position the image on the page use a `box` but set the `w` and `h` values to `0`.
"""
function pix_write_ps(
            stream::IO,
            pix::Pix;
            box::Union{PixBox, Nothing} = nothing,
            ppi::Integer                = Int32(300),
            scale::AbstractFloat        = Float32(1.0)
        )::Bool
    local result = false
    local data, size = pix_write_ps_i(pix, box, ppi, scale)

    if data != C_NULL
        try
            if size > 0
                unsafe_write(stream, data, size)
                result = true
            end
        finally
            lept_free(data)
        end
    end

    return result
end

# =================================================================================================
"""
    pix_write_ps(
        pix::Pix;
        box::Union{PixBox, Nothing} = nothing,
        ppi::Integer = Int32(300),
        scale::AbstractFloat = Float32(1.0)
    )::Union{Vector{UInt8}, Nothing}

Write an image to a byte array.  Returns 'nothing' if there was an error.

__Arguments:__

| T | Name     | Default        | Description
|:--| :------- | :------------- | :----------
| R | pix      |                | The image to write to a byte array.
| O | box      | `nothing`      | Location to have the image appear on the page.
| O | ppi      | `Int32(300)`   | The resolution to use for the image.
| O | scale    | `Float32(1.0)` | Scale the image on the page.

__Restrictions:__

  * `ppi` - Must be greater than 0.
  * `scale` - Must be greater than or equal to `0.0`.

__Details:__

If you don't want the image scaled set the `scale` to `1.0`.  If you want to position the image on
the page you can use a `box` however set `w` and `h` to `0`, then the image will be positioned
where you want it and not scaled.

If you want the image to be a specific size set the scale to `0.0` and use a box setting `w` and
`h` to the width and height you want in thousandth of an inch.  So a width of `5000` and a height of
`2000` will give you an image that is 5 inches across and 2 inches high.

If you want the image scaled set the `scale` to the factor you desire.  Value of `0.5` will cause
the image to be halfed, while a value of `2.0` will double the size of the image.  If you want to
also position the image on the page use a `box` but set the `w` and `h` values to `0`.
"""
function pix_write_ps(
            pix::Pix;
            box::Union{PixBox, Nothing} = nothing,
            ppi::Integer                = Int32(300),
            scale::AbstractFloat        = Float32(1.0)
        )::Union{Vector{UInt8}, Nothing}
    local result = nothing
    local data, size = pix_write_ps_i(pix, box, ppi, scale)

    if data != C_NULL
        try
            if size > 0
                result = Vector{UInt8}(undef, size)
                unsafe_copyto!(pointer(result), data, size)
            end
        finally
            lept_free(data)
        end
    end

    return result
end


text:
abstract type AbstractFeatureScaling end

#####
##### Zero-mean unit-variance feature scaling
#####

struct ZeroMeanUnitVarianceScaling{T} <: AbstractFeatureScaling
    μ :: T
    σ :: T
end

"""
    ZeroMeanUnitVarianceScaling(data)

Returns a feature scaler for `data` with zero mean and unit variance.
"""
function ZeroMeanUnitVarianceScaling(data)
    μ, σ = mean(data), std(data)
    return ZeroMeanUnitVarianceScaling(μ, σ)
end

scale(x, s::ZeroMeanUnitVarianceScaling) = (x .- s.μ) / s.σ
unscale(y, s::ZeroMeanUnitVarianceScaling) = s.σ * y .+ s.μ

#####
##### Min-max feature scaling
#####

struct MinMaxScaling{T} <: AbstractFeatureScaling
           a :: T
           b :: T
    data_min :: T
    data_max :: T
end

"""
    MinMaxScaling(data; a=0, b=1)

Returns a feature scaler for `data` with minimum `a` and `maximum `b`.
"""
function MinMaxScaling(data; a=0, b=1)
    data_min, data_max = extrema(data)
    return MinMaxScaling{typeof(data_min)}(a, b, data_min, data_max)
end

scale(x, s::MinMaxScaling) = s.a + (x - s.data_min) * (s.b - s.a) / (s.data_max - s.data_min)
unscale(y, s::MinMaxScaling) = s.data_min .+ (y .- s.a) * (s.data_max - s.data_min) / (s.b - s.a)

#####
##### Convenience functions
#####

(s::AbstractFeatureScaling)(x) = scale(x, s)
Base.inv(s::AbstractFeatureScaling) = y -> unscale(y, s)


text:
function measureXmag(A::fPEPS,
                     Ls::Vector{Environments},
                     Rs::Vector{Environments},
                     col;
                     kwargs...)
    s = Index(2, "Site,SpinInd")
    X = ITensor(s, s')
    is_cu         = is_gpu(A)
    X[s(1), s'(2)] = 0.5
    X[s(2), s'(1)] = 0.5
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    measuredX = zeros(Ny)
    op = is_cu ? cuITensor(X) : X
    Xs = [Operator([row=>col], [op], s, Field) for row in 1:Ny]
    tR = col == Nx ? dummyEnv : Rs[col+1]
    tL = col == 1  ? dummyEnv : Ls[col-1]
    AI = makeAncillaryIs(A, tL, tR, col)
    AF = makeAncillaryFs(A, tL, tR, Xs, col)
    fT = fieldTerms(A, tL, tR, (above=AI,), (above=AF,), Xs, 1, col, A[1, col])
    N  = buildN(A, tL, tR, (above=AI,), 1, col, A[1, col])
    for row in 1:Ny
        measuredX[row] = scalar(fT[row] * dag(A[1, col])')/scalar(N * dag(A[1, col]'))
    end
    return measuredX
end

function measureZmag(A::fPEPS,
                     Ls::Vector{Environments},
                     Rs::Vector{Environments},
                     col;
                     kwargs...)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    is_cu         = is_gpu(A)
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    Nx, Ny = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    measuredZ = zeros(Ny)
    op = is_cu ? cuITensor(Z) : Z
    Zs = [Operator([row=>col], [op], s, Field) for row in 1:Ny]
    #A  = intraColumnGauge(A, col; kwargs...)
    tR = col == Nx ? dummyEnv : Rs[col+1]
    tL = col == 1  ? dummyEnv : Ls[col-1]
    AI = makeAncillaryIs(A, tL, tR, col)
    AF = makeAncillaryFs(A, tL, tR, Zs, col)
    fT = fieldTerms(A, tL, tR, (above=AI,), (above=AF,), Zs, 1, col, A[1, col])
    N  = buildN(A, tL, tR, (above=AI,), 1, col, A[1, col])
    for row in 1:Ny
        measuredZ[row] = scalar(fT[row] * dag(A[1, col]'))/scalar(N * dag(A[1, col]'))
    end
    return measuredZ
end

function measureSmagVertical(A::fPEPS,
                             Ls::Vector{Environments},
                             Rs::Vector{Environments},
                             col;
                             kwargs...)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    Nx, Ny     = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    measuredSV = zeros(Ny)
    is_cu     = is_gpu(A)
    Z = is_cu ? cuITensor(Z) : Z
    P = is_cu ? cuITensor(P) : P
    M = is_cu ? cuITensor(M) : M
    SVs = Operator[]
    for row in 1:Ny-1
        push!(SVs, Operator([row=>col, row+1=>col], [0.5*P, M], s, Vertical))
        push!(SVs, Operator([row=>col, row+1=>col], [0.5*M, P], s, Vertical))
        push!(SVs, Operator([row=>col, row+1=>col], [Z, Z], s, Vertical))
    end
    tR = col == Nx ? dummyEnv : Rs[col+1]
    tL = col == 1  ? dummyEnv : Ls[col-1]
    AI = makeAncillaryIs(A, tL, tR, col)
    AV = makeAncillaryVs(A, tL, tR, SVs, col)
    vTs = verticalTerms(A, tL, tR, (above=AI,), (above=AV,), SVs, 1, col, A[1, col])
    N   = buildN(A, tL, tR, (above=AI,), 1, col, A[1, col])
    nrm = scalar(N * dag(A[1, col]'))
    for (vi, vT) in enumerate(vTs)
        row = SVs[vi].sites[1][1]
        measuredSV[row] += scalar(vT * dag(A[1, col]'))/nrm
    end
    return measuredSV
end

function measureSmagHorizontal(A::fPEPS,
                               Ls::Vector{Environments},
                               Rs::Vector{Environments},
                               col;
                               kwargs...)
    s = Index(2, "Site,SpinInd")
    Z = ITensor(s, s')
    Z[s(1), s'(1)] = 0.5
    Z[s(2), s'(2)] = -0.5
    P = ITensor(s, s')
    M = ITensor(s, s')
    P[s(1), s'(2)] = 1.0
    M[s(2), s'(1)] = 1.0
    is_cu     = is_gpu(A)
    Z = is_cu ? cuITensor(Z) : Z
    P = is_cu ? cuITensor(P) : P
    M = is_cu ? cuITensor(M) : M
    Nx, Ny = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    measuredSH = zeros(Ny)
    SHs = Operator[]
    for row in 1:Ny
        push!(SHs, Operator([row=>col, row=>col+1], [0.5*P, M], s, Horizontal))
        push!(SHs, Operator([row=>col, row=>col+1], [0.5*M, P], s, Horizontal))
        push!(SHs, Operator([row=>col, row=>col+1], [Z, Z], s, Horizontal))
    end
    tR = Rs[col+1]
    tL = col == 1 ? dummyEnv : Ls[col-1]
    AI = makeAncillaryIs(A, tL, tR, col)
    AS = makeAncillarySide(A, tR, tL, SHs, col, :right)
    hTs = connectRightTerms(A, tL, tR, (above=AI,), (above=AS,), SHs, 1, col, A[1, col])
    N  = buildN(A, tL, tR, (above=AI,), 1, col, A[1, col])
    nrm = scalar(N * dag(A[1, col]'))
    for (hi, hT) in enumerate(hTs)
        row = SHs[hi].sites[1][1]
        measuredSH[row] += scalar(hT * dag(A[1, col]'))/nrm
    end
    return measuredSH
end

function measure_correlators_heisenberg(A::fPEPS, H, Ls, Rs, sweep::Int; kwargs...)
    Ny, Nx = size(A)
    x_mag = zeros(Ny, Nx)
    z_mag = zeros(Ny, Nx)
    v_mag = zeros(Ny, Nx)
    h_mag = zeros(Ny, Nx)
    env_maxdim::Int = get(kwargs, :env_maxdim, 1)
    maxdim::Int     = get(kwargs, :maxdim, 1)
    mindim::Int     = get(kwargs, :mindim, 1)
    cutoff::Float64 = get(kwargs, :cutoff, 0.0)
    prefix::String  = get(kwargs, :prefix, "")
    max_gauge_iter::Int = get(kwargs, :max_gauge_iter, 50)
    if iseven(sweep)
        #R_s = buildRs(A, H; mindim=mindim, maxdim=maxdim, env_maxdim=env_maxdim)
        R_s = Vector{PEPS.Environments}(undef, Nx)
        for col in reverse(1:Nx)
            A  = intraColumnGauge(A, col; mindim=mindim, maxdim=maxdim)
            x_mag[:, col] = measureXmag(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            z_mag[:, col] = measureZmag(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            v_mag[:, col] = measureSmagVertical(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            if col < Nx
                h_mag[:, col] = measureSmagHorizontal(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            end
            if col > 1
                A  = gaugeColumn(A, col, :left; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, overlap_cutoff=0.999, max_gauge_iter=max_gauge_iter)
                if col < Nx
                    R_s[col] = buildNextEnvironment(A, R_s[col+1], H, :right, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                else
                    right_H_terms = getDirectional(H[Nx-1], Horizontal)
                    R_s[col] = buildEdgeEnvironment(A, H, right_H_terms, :right, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                end
            end
        end
    else
        #L_s = buildLs(A, H; mindim=mindim, maxdim=maxdim, env_maxdim=env_maxdim)
        L_s = Vector{PEPS.Environments}(undef, Nx)
        for col in 1:Nx
            A  = intraColumnGauge(A, col; mindim=mindim, maxdim=maxdim)
            x_mag[:, col] = measureXmag(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            z_mag[:, col] = measureZmag(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            v_mag[:, col] = measureSmagVertical(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            if col < Nx
                h_mag[:, col] = measureSmagHorizontal(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            end
            if col < Nx
                A  = gaugeColumn(A, col, :right; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, overlap_cutoff=0.999, max_gauge_iter=max_gauge_iter)
                if col > 1
                    L_s[col] = buildNextEnvironment(A, L_s[col-1], H, :left, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                else
                    left_H_terms = getDirectional(H[1], Horizontal)
                    L_s[col] = buildEdgeEnvironment(A, H, left_H_terms, :left, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                end
            end
        end
    end
    writedlm(prefix*"_$(sweep)_x", x_mag)
    writedlm(prefix*"_$(sweep)_z", z_mag)
    writedlm(prefix*"_$(sweep)_v", v_mag)
    writedlm(prefix*"_$(sweep)_h", h_mag)
end

function measure_correlators_ising(A::fPEPS, H, Ls, Rs, sweep::Int; kwargs...)
    Ny, Nx = size(A)
    x_mag = zeros(Ny, Nx)
    z_mag = zeros(Ny, Nx)
    env_maxdim::Int = get(kwargs, :env_maxdim, 1)
    maxdim::Int     = get(kwargs, :maxdim, 1)
    mindim::Int     = get(kwargs, :mindim, 1)
    cutoff::Float64 = get(kwargs, :cutoff, 0.0)
    prefix::String  = get(kwargs, :prefix, "")
    max_gauge_iter::Int = get(kwargs, :max_gauge_iter, 50)
    if iseven(sweep)
        R_s = Vector{PEPS.Environments}(undef, Nx)
        for col in reverse(1:Nx)
            A  = intraColumnGauge(A, col; mindim=mindim, maxdim=maxdim)
            x_mag[:, col] = measureXmag(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            z_mag[:, col] = measureZmag(A, Ls, R_s, col; mindim=mindim, maxdim=maxdim)
            if col > 1
                A  = gaugeColumn(A, col, :left; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, overlap_cutoff=0.999, max_gauge_iter=max_gauge_iter)
                if col < Nx
                    R_s[col] = buildNextEnvironment(A, R_s[col+1], H, :right, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                else
                    right_H_terms = getDirectional(vcat(H[:, Nx-1]), Horizontal)
                    R_s[col] = buildEdgeEnvironment(A, H, right_H_terms, :right, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                end
            end
        end
    else
        L_s = Vector{PEPS.Environments}(undef, Nx)
        for col in 1:Nx
            A  = intraColumnGauge(A, col; mindim=mindim, maxdim=maxdim)
            x_mag[:, col] = measureXmag(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            z_mag[:, col] = measureZmag(A, L_s, Rs, col; mindim=mindim, maxdim=maxdim)
            if col < Nx
                A  = gaugeColumn(A, col, :right; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, overlap_cutoff=0.999, max_gauge_iter=max_gauge_iter)
                if col > 1
                    L_s[col] = buildNextEnvironment(A, L_s[col-1], H, :left, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                else
                    left_H_terms = getDirectional(vcat(H[:, 1]), Horizontal)
                    L_s[col] = buildEdgeEnvironment(A, H, left_H_terms, :left, col; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
                end
            end
        end
    end
    writedlm(prefix*"_$(sweep)_x", x_mag)
    writedlm(prefix*"_$(sweep)_z", z_mag)
end


text:
module DistStat

import MPI
import MPI: COMM_WORLD
using Requires

function __init__()
    MPI.Initialized() || MPI.Init()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
        include("cuda.jl")
        set_device!()
        #CuArrays.allowscalar(false)
    end

end

@inline function Size()
    MPI.Comm_size(COMM_WORLD)
end

@inline function Rank()
    MPI.Comm_rank(COMM_WORLD)
end

include("distdirectives.jl")
include("distarray.jl")
include("distlinalg.jl")
include("reduce.jl")
include("accumulate.jl")
include("broadcast.jl")
include("arrayfunctions.jl")
include("utils.jl")

end # module


text:
export CbcOptimizer

using MathOptInterface
using Cbc.CbcCInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const CbcCI = CbcCInterface

mutable struct CbcOptimizer <: MOI.AbstractOptimizer
    inner::CbcModel
    CbcOptimizer() = new(CbcModel()) # Initializes with an empty model
end

struct CbcModelFormat
    num_rows::Int
    num_cols::Int
    row_idx::Vector{Int}
    col_idx::Vector{Int}
    values::Vector{Float64}
    col_lb::Vector{Float64}
    col_ub::Vector{Float64}
    obj::Vector{Float64}
    row_lb::Vector{Float64}
    row_ub::Vector{Float64}
    function CbcModelFormat(num_rows::Int, num_cols::Int)
        obj = fill(0.0, num_cols)
        row_idx = Int[]
        col_idx = Int[]
        values = Float64[]
        col_lb = fill(-Inf, num_cols)
        col_ub = fill(Inf, num_cols)
        row_lb = fill(-Inf, num_rows)
        row_ub = fill(Inf, num_rows)
        constraint_matrix = Tuple{Int,Int,Float64}[]
        new(num_rows, num_cols, row_idx, col_idx, values, col_lb, col_ub, obj, row_lb, row_ub)
    end
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.SingleVariable, s::MOI.EqualTo)
    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.value
    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.value
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.SingleVariable, s::MOI.LessThan)
    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.upper
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.SingleVariable, s::MOI.GreaterThan)
    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.lower
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.SingleVariable, s::MOI.Interval)
    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.lower
    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.upper
end

function push_terms(row_idx::Vector{Int}, col_idx::Vector{Int}, values::Vector{Float64},
    ci::MOI.ConstraintIndex{F,S}, terms::Vector{MOI.ScalarAffineTerm{Float64}},
    mapping::MOIU.IndexMap) where {F,S}
    for term in terms
        push!(row_idx, mapping.conmap[ci].value)
        push!(col_idx, mapping.varmap[term.variable_index].value)
        push!(values, term.coefficient)
    end
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat,
    mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction, s::MOI.EqualTo)
    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.value - f.constant
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.value - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat,
    mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction, s::MOI.GreaterThan)
    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.lower - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.ScalarAffineFunction, s::MOI.LessThan)
    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.upper - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex, cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.ScalarAffineFunction, s::MOI.Interval)
    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.upper - f.constant
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.lower - f.constant
end

function load_obj(cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.ScalarAffineFunction)
    # We need to increment values of objective function with += to handle cases like $x_1 + x_2 + x_1$
    # This is safe becasue objective function is initialized with zeros in the constructor
    for term in f.terms
        cbc_model_format.obj[mapping.varmap[term.variable_index].value] += term.coefficient
    end
end

function copy_constraints!(cbc_model_format::CbcModelFormat, user_optimizer::MOI.ModelLike,
    mapping::MOIU.IndexMap)
    for (F,S) in MOI.get(user_optimizer, MOI.ListOfConstraints())
        if !(S <: Union{MOI.ZeroOne, MOI.Integer})
            for ci in MOI.get(user_optimizer, MOI.ListOfConstraintIndices{F,S}())
                f = MOI.get(user_optimizer, MOI.ConstraintFunction(), ci)
                s = MOI.get(user_optimizer,  MOI.ConstraintSet(), ci)
                load_constraint(ci, cbc_model_format, mapping, f, s)
            end
        end
    end
end

function update_bounds_for_binary_vars!(col_lb::Vector{Float64}, col_ub::Vector{Float64}, zero_one_indices::Vector{Int})
    for idx in zero_one_indices
        if col_lb[idx] < 0.0
            col_lb[idx] = 0.0
        end
        if col_ub[idx] > 1.0
            col_ub[idx] = 1.0
        end
    end
end

function update_zeroone_indices(user_optimizer::MOI.ModelLike, mapping::MOIU.IndexMap,
    ci::Vector{MOI.ConstraintIndex{F,S}}, zero_one_indices::Vector{Int}) where {F, S}
    for i in 1:length(ci)
        f = MOI.get(user_optimizer, MOI.ConstraintFunction(), ci[i])
        push!(zero_one_indices, mapping.varmap[f.variable].value)
    end
end

function update_integer_indices(user_optimizer::MOI.ModelLike, mapping::MOIU.IndexMap,
    ci::Vector{MOI.ConstraintIndex{F,S}}, integer_indices::Vector{Int}) where {F, S}
    for i in 1:length(ci)
        f = MOI.get(user_optimizer, MOI.ConstraintFunction(), ci[i])
        push!(integer_indices, mapping.varmap[f.variable].value)
    end
end

"""
    function copy!(cbc_optimizer, user_optimizer; copynames=false)

Receive a cbc_optimizer which contains the pointer to the cbc C object and instantiate the object cbc_model_format::CbcModelFormat based on user_optimizer::AbstractModel (also provided by the user).
Function loadProblem of CbcCInterface requires all information stored in cbc_model_format.
"""
function MOI.copy!(cbc_optimizer::CbcOptimizer,
    user_optimizer::MOI.ModelLike; copynames=false)

    mapping = MOIU.IndexMap()

    num_cols = MOI.get(user_optimizer, MOI.NumberOfVariables())
    var_index = MOI.get(user_optimizer, MOI.ListOfVariableIndices())
    for i in 1:num_cols
        mapping.varmap[var_index[i]] = MOI.VariableIndex(i)
    end

    zero_one_indices = Int[]
    integer_indices = Int[]
    list_of_constraints = MOI.get(user_optimizer, MOI.ListOfConstraints())
    num_rows = 0
    for (F,S) in list_of_constraints
        if !(MOI.supportsconstraint(cbc_optimizer, F, S))
            return MOI.CopyResult(MOI.CopyUnsupportedConstraint,
            "Cbc MOI Interface does not support constraints of type " * (F,S) * ".", nothing)
        end

        ci = MOI.get(user_optimizer, MOI.ListOfConstraintIndices{F,S}())

        if F == MOI.SingleVariable
            if S == MOI.ZeroOne
                update_zeroone_indices(user_optimizer, mapping, ci, zero_one_indices)
            elseif S == MOI.Integer
                update_integer_indices(user_optimizer, mapping, ci, integer_indices)
            end
        else
            ## Update conmap for (F,S) for F != MOI.SingleVariable
            ## Single variables are treated by bounds in Cbc, so no
            ## need to add a row
            for i in 1:length(ci)
                mapping.conmap[ci[i]] = MOI.ConstraintIndex{F,S}(num_rows + i)
            end
            num_rows += MOI.get(user_optimizer, MOI.NumberOfConstraints{F,S}())
        end
    end

    cbc_model_format = CbcModelFormat(num_rows, num_cols)

    copy_constraints!(cbc_model_format, user_optimizer, mapping)
    update_bounds_for_binary_vars!(cbc_model_format.col_lb, cbc_model_format.col_ub, zero_one_indices)

    ## Copy objective function
    objF = MOI.get(user_optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    load_obj(cbc_model_format, mapping, objF)
    sense = MOI.get(user_optimizer, MOI.ObjectiveSense())
    MOI.set!(cbc_optimizer, MOI.ObjectiveSense(), sense)

    ## Load the problem to Cbc
    CbcCI.loadProblem(cbc_optimizer.inner, sparse(cbc_model_format.row_idx, cbc_model_format.col_idx,
                                                  cbc_model_format.values, cbc_model_format.num_rows,
                                                  cbc_model_format.num_cols),
                      cbc_model_format.col_lb, cbc_model_format.col_ub, cbc_model_format.obj,
                      cbc_model_format.row_lb, cbc_model_format.row_ub)

    empty!(cbc_model_format.row_idx)
    empty!(cbc_model_format.col_idx)
    empty!(cbc_model_format.values)

    ## Set integer variables
    for idx in vcat(integer_indices, zero_one_indices)
        CbcCI.setInteger(cbc_optimizer.inner, idx-1)
    end

    return MOI.CopyResult(MOI.CopySuccess, "Model was copied succefully.", MOIU.IndexMap(mapping.varmap, mapping.conmap))
end

function MOI.optimize!(cbc_optimizer::CbcOptimizer)
    # Call solve function
    CbcCI.solve(cbc_optimizer.inner)
end

## canadd, canset, canget functions

function MOI.canaddvariable(cbc_optimizer::CbcOptimizer)
    return false
end

## supports constraints

MOI.supportsconstraint(::CbcOptimizer, ::Type{<:Union{MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable}},
::Type{<:Union{MOI.EqualTo{Float64}, MOI.Interval{Float64}, MOI.LessThan{Float64},
MOI.GreaterThan{Float64}, MOI.ZeroOne, MOI.Integer}}) = true

MOI.supports(cbc_optimizer::CbcOptimizer, object::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true

## Set functions

function MOI.write(cbc_optimizer::CbcOptimizer, filename::String)
    if !endswith("filename", "mps")
        error("CbcOptimizer only supports writing .mps files")
    else
        writeMps(cbc_optimizer.inner, filename)
    end
end

# empty!
function MOI.empty!(cbc_optimizer::CbcOptimizer)
    cbc_optimizer.inner = CbcModel()
end

function MOI.set!(cbc_optimizer::CbcOptimizer, object::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MaxSense
        CbcCI.setObjSense(cbc_optimizer.inner, -1)
    else ## Other senses are set as minimization (cbc default)
        CbcCI.setObjSense(cbc_optimizer.inner, 1)
    end
end

## Get functions

function MOI.canget(cbc_optimizer::CbcOptimizer, object::MOI.PrimalStatus)
    if object.N != 1
        return false
    end
    return MOI.get(cbc_optimizer, MOI.ResultCount()) == 1
end

MOI.canget(cbc_optimizer::CbcOptimizer, object::Union{MOI.NodeCount, MOI.ResultCount,
MOI.TerminationStatus, MOI.ObjectiveSense, MOI.ObjectiveValue, MOI.ObjectiveBound, MOI.NumberOfVariables}) = true

MOI.canget(cbc_optimizer::CbcOptimizer, object::MOI.VariablePrimal, indexTypeOrObject::Type{MOI.VariableIndex}) = true

function MOI.isempty(cbc_optimizer::CbcOptimizer)
    return (CbcCI.getNumCols(cbc_optimizer.inner) == 0 && CbcCI.getNumRows(cbc_optimizer.inner) == 0)
end

MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.NumberOfVariables) = getNumCols(cbc_optimizer.inner)

MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.ObjectiveBound) = CbcCI.getBestPossibleObjValue(cbc_optimizer.inner)

MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.NodeCount) = CbcCI.getNodeCount(cbc_optimizer.inner)

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.ObjectiveValue)
    return CbcCI.getObjValue(cbc_optimizer.inner)
end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.VariablePrimal, ref::MOI.VariableIndex)
    variablePrimals = CbcCI.getColSolution(cbc_optimizer.inner)
    return variablePrimals[ref.value]
end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.VariablePrimal, ref::Vector{MOI.VariableIndex})
    variablePrimals = CbcCI.getColSolution(cbc_optimizer.inner)
    return [variablePrimals[vi.value] for vi in ref]
end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.ResultCount)
    if (isProvenInfeasible(cbc_optimizer.inner) || isContinuousUnbounded(cbc_optimizer.inner)
        || isAbandoned(cbc_optimizer.inner) || CbcCI.getObjValue(cbc_optimizer.inner) >= 1e300)
        return 0
    end
    return 1
end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.ObjectiveSense)
    CbcCI.getObjSense(cbc_optimizer.inner) == 1 && return MOI.MinSense
    CbcCI.getObjSense(cbc_optimizer.inner) == -1 && return MOI.MaxSense
end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.TerminationStatus)

    if isProvenInfeasible(cbc_optimizer.inner)
        return MOI.InfeasibleNoResult
    elseif isContinuousUnbounded(cbc_optimizer.inner)
        return MOI.InfeasibleOrUnbounded
    elseif isNodeLimitReached(cbc_optimizer.inner)
        return MOI.NodeLimit
    elseif isSecondsLimitReached(cbc_optimizer.inner)
        return MOI.TimeLimit
    elseif isSolutionLimitReached(cbc_optimizer.inner)
        return MOI.SolutionLimit
    elseif (isProvenOptimal(cbc_optimizer.inner) || isInitialSolveProvenOptimal(cbc_optimizer.inner)
        || MOI.get(cbc_optimizer, MOI.ResultCount()) == 1)
        return MOI.Success
    elseif isAbandoned(cbc_optimizer.inner)
        return MOI.Interrupted
    else
        error("Internal error: Unrecognized solution status")
    end

end

function MOI.get(cbc_optimizer::CbcOptimizer, object::MOI.PrimalStatus)
    if isProvenOptimal(cbc_optimizer.inner) || isInitialSolveProvenOptimal(cbc_optimizer.inner)
        return MOI.FeasiblePoint
    elseif isProvenInfeasible(cbc_optimizer.inner)
        return MOI.InfeasiblePoint
    end
end