text:
"""
Function to obtain series of states out of DAE Solution. It receives the solution, the dynamical system
and a tuple containing the symbol name of the Dynamic Injection device and the symbol of the state.

"""

function get_state_series(sim::DynamicSimulation, ref::NTuple{2, Symbol})
    ix = get(sim.dyn_system.global_state_index[ref[1]], ref[2], nothing)
    return sim.solution.t, [value[ix] for value in sim.solution.u]
end

"""
Function to obtain the voltage magnitude series out of the DAE Solution. It receives the solution, the dynamical system and the bus number.

"""

function get_voltagemag_series(sim::DynamicSimulation, bus_number::Int64)
    bus_size = length(sim.dyn_system.buses)
    return sim.solution.t, [sqrt(value[bus_number]^2 + value[bus_number+bus_size]^2) for value in sim.solution.u]
end

"""
Function to print initial states. It receives the vector of initial states and the dynamical system.
"""

function print_init_states(x0::Vector{Float64}, dyn_system::DynamicSystem)
    for (ix, val_sys) in dyn_system.global_state_index
        println(ix)
        println("====================")
        for (k, val) in val_sys
            print(k, " ", x0[val], "\n")
        end
        println("====================")
    end
end


text:
# ubrk.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ubrk.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
UBRK defines constants for ICU break iterator API:
CHARACTER, TITLE, COUNT, DONE,
WORD, WORD_TAG, WORD_NONE, WORD_NONE_LIMIT, WORD_NUMBER, WORD_NUMBER_LIMIT,
WORD_LETTER, WORD_LETTER_LIMIT, WORD_KANA, WORD_KANA_LIMIT, WORD_IDEO, WORD_IDEO_LIMIT,
LINE, LINE_TAG, LINE_SOFT, LINE_SOFT_LIMIT, LINE_HARD, LINE_HARD_LIMIT,
SENTENCE, SENTENCE_TAG, SENTENCE_TERM, SENTENCE_TERM_LIMIT, SENTENCE_SEP, SENTENCE_SEP_LIMIT
"""
module UBRK

const CHARACTER = Int32(0)
const WORD      = Int32(1)
const LINE      = Int32(2)
const SENTENCE  = Int32(3)
const TITLE     = Int32(4) # Deprecated API (use word break iterator)
const COUNT     = Int32(5)

"""Value indicating all text boundaries have been returned."""
const DONE      = Int32(-1)

"""
    Constants for the word break tags returned by ICU.get_rule_status().
    A range of values is defined for each category of word, to allow for
    further subdivisions of a category in future releases.
    Applications should check for tag values falling within the range,
    rather than for single individual values.
"""
const WORD_TAG = 0

"""
    Tag value for "words" that do not fit into any of other categories.
    Includes spaces and most punctuation.
"""
const WORD_NONE  = Int32(0)
"""
    Upper bound for tags for uncategorized words.
"""
const WORD_NONE_LIMIT = Int32(100)
"""
    Tag value for words that appear to be numbers, lower limit.
"""
const WORD_NUMBER     = Int32(100)
"""
    Tag value for words that appear to be numbers, upper limit.
"""
const WORD_NUMBER_LIMIT   = Int32(200)
"""
    Tag value for words that contain letters, excluding hiragana, katakana or
    ideographic characters, lower limit.
"""
const WORD_LETTER         = Int32(200)
"""
    Tag value for words containing letters, upper limit
"""
const WORD_LETTER_LIMIT   = Int32(300)
"""
    Tag value for words containing kana characters, lower limit
"""
const WORD_KANA           = Int32(300)
"""
    Tag value for words containing kana characters, upper limit
"""
const WORD_KANA_LIMIT     = Int32(400)
"""
    Tag value for words containing ideographic characters, lower limit
"""
const WORD_IDEO           = Int32(400)
"""
    Tag value for words containing ideographic characters, upper limit
"""
const WORD_IDEO_LIMIT     = Int32(500)

"""
    Constants for the line break tags returned by ICU.get_rule_status().
    A range of values is defined for each category of word, to allow for
    further subdivisions of a category in future releases.
    Applications should check for tag values falling within the range,
    rather than for single individual values.
"""
const LINE_TAG = 0

"""
    Tag value for soft line breaks, positions at which a line break
    is acceptable but not required
"""
const LINE_SOFT            = Int32(0)
"""
    Upper bound for soft line breaks.
"""
const LINE_SOFT_LIMIT      = Int32(100)
"""
    Tag value for a hard, or mandatory line break
"""
const LINE_HARD            = Int32(100)
"""
    Upper bound for hard line breaks.
"""
const LINE_HARD_LIMIT      = Int32(200)

"""
    Constants for the sentence break tags returned by ICU.get_rule_status().
    A range of values is defined for each category of sentence, to allow for
    further subdivisions of a category in future releases.
    Applications should check for tag values falling within the range,
    rather than for single individual values.
"""
const SENTENCE_TAG = 0

"""
    Tag value for for sentences  ending with a sentence terminator
    ('.', '?', '!', etc.) character, possibly followed by a
    hard separator (CR, LF, PS, etc.)
"""
const SENTENCE_TERM       = Int32(0)
"""
    Upper bound for tags for sentences ended by sentence terminators.
"""
const SENTENCE_TERM_LIMIT = Int32(100)
"""
    Tag value for for sentences that do not contain an ending
    sentence terminator ('.', '?', '!', etc.) character, but
    are ended only by a hard separator (CR, LF, PS, etc.) or end of input.
"""
const SENTENCE_SEP        = Int32(100)
"""
    Upper bound for tags for sentences ended by a separator.
"""
const SENTENCE_SEP_LIMIT  = Int32(200)
end # module UBRK

export UBRK, UBrk

macro libbrk(s)     ; _libicu(s, iculib,     "ubrk_")     ; end

const UBrkType = Int32

"""
    UBrk is a type along with methods for finding the location of boundaries in text.
    A UBrk maintains a current position and scan over text returning the index of characters
    where boundaries occur.

    Line boundary analysis determines where a text string can be broken when line-wrapping.
    The mechanism correctly handles punctuation and hyphenated words.

    Note: The locale keyword "lb" can be used to modify line break behavior according to
    the CSS level 3 line-break options, see <http://dev.w3.org/csswg/css-text/#line-breaking>.
    For example: "ja@lb=strict", "zh@lb=loose".

    Sentence boundary analysis allows selection with correct interpretation of periods within
    numbers and abbreviations, and trailing punctuation marks such as quotation marks and
    parentheses.

    Note: The locale keyword "ss" can be used to enable use of segmentation suppression data
    (preventing breaks in English after abbreviations such as "Mr." or "Est.", for example),
    as follows: "en@ss=standard".

    Word boundary analysis is used by search and replace functions, as well as within text
    editing applications that allow the user to select words with a double click.
    Word selection provides correct interpretation of punctuation marks within and following
    words. Characters that are not part of a word, such as symbols or punctuation marks,
    have word-breaks on both sides.

    Character boundary analysis identifies the boundaries of "Extended Grapheme Clusters",
    which are groupings of codepoints that should be treated as character-like units for
    many text operations.

    Please see Unicode Standard Annex #29, Unicode Text Segmentation,
    http://www.unicode.org/reports/tr29/ for additional information
    on grapheme clusters and guidelines on their use.

    Title boundary analysis locates all positions, typically starts of words, that should
    be set to Title Case when title casing the text.

    The text boundary positions are found according to the rules described in
    Unicode Standard Annex #29, Text Boundaries, and Unicode Standard Annex #14,
    Line Breaking Properties.  These are available at http://www.unicode.org/reports/tr14/
    and http://www.unicode.org/reports/tr29/.
"""
type UBrk
    p::Ptr{Void}
    s
    r

    function UBrk(kind::Integer, s::Vector{UInt16}, loc::ASCIIStr)
        err = UErrorCode[0]
        p = ccall(@libbrk(open), Ptr{Void},
                  (UBrkType, Cstring, Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  kind, loc, s, length(s), err)
        @assert SUCCESS(err[1])
        # Retain pointer to input vector, otherwise it might be GCed
        self = new(p, s, Void())
        finalizer(self, close)
        self
    end
    function UBrk(kind::Integer, s::Vector{UInt16})
        err = UErrorCode[0]
        p = ccall(@libbrk(open), Ptr{Void},
                  (UBrkType, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
                  kind, C_NULL, s, length(s), err)
        @assert SUCCESS(err[1])
        # Retain pointer to input vector, otherwise it might be GCed
        self = new(p, s, Void())
        finalizer(self, close)
        self
    end
    function UBrk(rules::Vector{UInt16}, s::Vector{UInt16})
        err = UErrorCode[0]
        # Temporary disable UParseError and pass C_NULL
        #p_err = Ref(UParseError())
        p = ccall(@libbrk(openRules), Ptr{Void},
                  (Ptr{UChar}, Int32, Ptr{UChar}, Int32, Ptr{UParseError}, Ptr{UErrorCode}),
                  rules, length(rules), s, length(s), C_NULL, err)
        @assert SUCCESS(err[1])
        # Retain pointer to input vector, otherwise it might be GCed
        self = new(p, s, rules)
        finalizer(self, close)
        self
    end
end
UBrk(kind::Integer, s::UniStr) =
    UBrk(kind, Vector{UInt16}(s)[1:end-1])
UBrk(kind::Integer, s::UniStr, loc::AbstractString) =
    UBrk(kind, Vector{UInt16}(s)[1:end-1], ascii(loc))
UBrk(rules::UniStr, s::UniStr) =
    UBrk(Vector{UInt16}(rules)[1:end-1], Vector{UInt16}(s)[1:end-1])

"""
    Close the Break Iterator and return any resource, if not already closed
"""
close(bi::UBrk) =
    bi.p == C_NULL ||
        (ccall(@libbrk(close), Void, (Ptr{Void},), bi.p) ; bi.p = C_NULL ; bi.s = Void() ; bi.r = Void())

"""
    Determine the most recently-returned text boundary.

    Arguments:
    bi - The break iterator to use.
    Returns the character index most recently returned by next, previous, first, or last.
"""
current(bi::UBrk) = ccall(@libbrk(current), Int32, (Ptr{Void},), bi.p)

"""
    Advance the iterator to the boundary following the current boundary.
    Returns the character index of the next text boundary, or UBRK.DONE
    if all text boundaries have been returned.
"""
next(bi::UBrk) = ccall(@libbrk(next), Int32, (Ptr{Void},), bi.p)

"""
    Set the iterator position to the boundary preceding the current boundary.
    Returns the character index of the preceding text boundary, or UBRK.DONE
    if all text boundaries have been returned.
"""
previous(bi::UBrk) = ccall(@libbrk(previous), Int32, (Ptr{Void},), bi.p)

"""
    Set the iterator position to zero, the start of the text being scanned.
    Returns the new iterator position (zero).
"""
first(bi::UBrk) = ccall(@libbrk(first), Int32, (Ptr{Void},), bi.p)

"""
    Set the iterator position to the index immediately beyond the last character in
    the text being scanned.
    This is not the same as the last character.
    Returns the character offset immediately beyond the last character in the
    text being scanned.
"""
last(bi::UBrk) = ccall(@libbrk(last), Int32, (Ptr{Void},), bi.p)

"""
    Set the iterator position to the first boundary preceding the specified offset.
    The new position is always smaller than offset, or UBRK.DONE

    Arguments:
    bi - The break iterator to use
    offset - The offset to begin scanning.
    Returns the text boundary preceding offset, or UBRK.DONE
"""
preceding(bi::UBrk, off) = ccall(@libbrk(preceding), Int32, (Ptr{Void}, Int32), bi.p, off)

"""
    Advance the iterator to the first boundary following the specified offset.
    The value returned is always greater than offset, or UBRK.DONE

    Arguments:
    - bi - The break iterator to use.
    - offset - The offset to begin scanning.
    Returns the text boundary following offset, or UBRK.DONE
"""
following(bi::UBrk, off) = ccall(@libbrk(following), Int32, (Ptr{Void}, Int32), bi.p, off)

"""
    Get a locale for which text breaking information is available.
    A UBrk in a locale returned by this function will perform the correct
    text breaking for the locale.

    Arguments:
    index - The index of the desired locale.

    Returns: A locale for which number text breaking information is available, or 0 if none.
"""
function get_available(index)
    loc = ccall(@libbrk(getAvailable), Ptr{UInt8}, (Int32, ), index)
    loc == C_NULL ? "" : ASCIIStr(loc)
end

"""
    Determine how many locales have text breaking information available.
    This function is most useful as determining the loop ending condition for
    calls to get_available.

    Returns: the number of locales for which text breaking information is available.
"""
count_available() = ccall(@libbrk(countAvailable), Int32, ())

"""
    Returns true if the specfied position is a boundary position.  As a side
    effect, leaves the iterator pointing to the first boundary position at
    or after "offset".

    Arguments:
    - bi - The break iterator to use.
    - offset - The offset to check.

    Returns true if "offset" is a boundary position.
"""
isboundary(bi::UBrk, off) = ccall(@libbrk(isBoundary), Int32, (Ptr{Void}, Int32), bi.p, off) != 0

"""
    Return the status from the break rule that determined the most recently returned break
    position.  The values appear in the rule source within brackets, {123}, for example.
    For rules that do not specify a status, a default value of 0 is returned.
"""
get_rule_status(bi::UBrk) =
    ccall(@libbrk(getRuleStatus), Int32, (Ptr{Void},), bi.p)

"""
    Get the statuses from the break rules that determined the most recently
    returned break position.  The values appear in the rule source within brackets,
    {123}, for example.  The default status value for rules that do not explicitly
    provide one is zero.

    For word break iterators, the possible values are defined module UBRK

    Arguments:
    bi - The break iterator to use
    fillinvec - an array to be filled in with the status values.
    capacity - the length of the supplied vector.  A length of zero causes
               the function to return the number of status values, in the
               normal way, without attemtping to store any values.
    status - receives error codes.

    Returns: The number of rule status values from rules that determined
             the most recent boundary returned by the break iterator.
"""
function get_rule_status_vec(bi::UBrk, fillinvec::Vector{Int32})
    err = UErrorCode[0]
    cnt = ccall(@libbrk(getRuleStatusVec), Int32,
                (Ptr{Void}, Ptr{Int32}, Int32,  Ptr{UErrorCode}),
                bi.p, fillinvec, length(fillinvec), err)
    @assert SUCCESS(err[1])
    cnt
end

function get_rule_status_len(bi::UBrk)
    err = UErrorCode[0]
    cnt = ccall(@libbrk(getRuleStatusVec), Int32,
                (Ptr{Void}, Ptr{Int32}, Int32,  Ptr{UErrorCode}),
                bi.p, C_NULL, 0, err)
    @assert SUCCESS(err[1])
    cnt
end

"""
    Return the locale of the break iterator. You can choose between
    the valid and the actual locale.

    Arguments:
    bi - break iterator
    type - locale type (valid or actual)
    status - error code

    Return: locale string
"""
function get_locale_by_type(bi::UBrk, typ, status)
    # ccall(@libbrk(getLocaleByType), Ptr{UInt8},
    #       (const UBreakIterator *bi, ULocDataLocaleType type, UErrorCode* status);
end

"""
    Sets an existing iterator to point to a new piece of text

    Arguments:
    bi - The iterator to use
    text - The text to be set

    Returns status code
"""
function set! end

set!(bi::UBrk, str::AbstractString) = set!(bi, cvt_utf16(str))
set!(bi::UBrk, str::UniStr) = set!(bi, Vector{UInt16}(str))

function set!(bi::UBrk, v::Vector{UInt16})
    err = UErrorCode[0]
    ccall(@libbrk(setText), Void,
          (Ptr{Void}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
          bi.p, v, length(v)-1, err)
    # Retain pointer so that it won't be GCed
    bi.s = v
    @assert SUCCESS(err[1])
    err[1]
end

"""
    Sets an existing iterator to point to a new piece of text.

    All index positions returned by break iterator functions are
    native indices from the UText. For example, when breaking UTF-8
    encoded text, the break positions returned by next, previous, etc.
    will be UTF-8 string indices, not UTF-16 positions.

    Arguments:
    bi - The iterator to use
    text - The text to be set.
           This function makes a shallow clone of the supplied UText.  This means
           that the caller is free to immediately close or otherwise reuse the
           UText that was passed as a parameter, but that the underlying text itself
           must not be altered while being referenced by the break iterator.

    Returns status code
"""
function set!(bi::UBrk, t::UText)
    err = UErrorCode[0]
    ccall(@libbrk(setUText), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{UErrorCode}),
          bi.p, t.p, err)
    @assert SUCCESS(err[1])
    err[1]
end

ext: jl


text:
using Weave

cd("src/")
weave("index.Rmd", informat="markdown",
      out_path = ".", doctype = "github")
cd("../")


text:
"""
    function dragon_and_cow()

A dragon smiting a cow, possible credit to kube@csua.berkeley.edu

# Example

```jldoctest
julia> cowsay("Roast beef: it's what's for dinner", cow=Cowsay.dragon_and_cow)
 ____________________________________
< Roast beef: it's what's for dinner >
 ------------------------------------
                       \\                    ^    /^
                        \\                  / \\  // \\
                         \\   |\\___/|      /   \\//  .\\
                          \\  /O  O  \\__  /    //  | \\ \\           *----*
                            /     /  \\/_/    //   |  \\  \\          \\   |
                            @___@`    \\/_   //    |   \\   \\         \\/\\ \\
                           0/0/|       \\/_ //     |    \\    \\         \\  \\
                       0/0/0/0/|        \\///      |     \\     \\       |  |
                    0/0/0/0/0/_|_ /   (  //       |      \\     _\\     |  /
                 0/0/0/0/0/0/`/,_ _ _/  ) ; -.    |    _ _\\.-~       /   /
                             ,-}        _      *-.|.-~-.           .~    ~
            \\     \\__/        `/\\      /                 ~-. _ .-~      /
             \\____(oo)           *.   }            {                   /
             (    (--)          .----~-.\\        \\-`                 .~
             //__\\\\  \\__ Ack!   ///.----..<        \\             _ -~
            //    \\\\               ///-._ _ _ _ _ _ _{^ - - - - ~
```
"""
function dragon_and_cow(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
                       $thoughts                    ^    /^
                        $thoughts                  / \\  // \\
                         $thoughts   |\\___/|      /   \\//  .\\
                          $thoughts  /O  O  \\__  /    //  | \\ \\           *----*
                            /     /  \\/_/    //   |  \\  \\          \\   |
                            @___@`    \\/_   //    |   \\   \\         \\/\\ \\
                           0/0/|       \\/_ //     |    \\    \\         \\  \\
                       0/0/0/0/|        \\///      |     \\     \\       |  |
                    0/0/0/0/0/_|_ /   (  //       |      \\     _\\     |  /
                 0/0/0/0/0/0/\`/,_ _ _/  ) ; -.    |    _ _\\.-~       /   /
                             ,-}        _      *-.|.-~-.           .~    ~
            \\     \\__/        \`/\\      /                 ~-. _ .-~      /
             \\____($eyes)           *.   }            {                   /
             (    (--)          .----~-.\\        \\-\`                 .~
             //__\\\\  \\__ Ack!   ///.----..<        \\             _ -~
            //    \\\\               ///-._ _ _ _ _ _ _{^ - - - - ~
"""
return the_cow
end


text:
using Documenter, MixedLayerThermoclineDynamics

makedocs(modules = [MixedLayerThermoclineDynamics],
        sitename = "MixedLayerThermoclineDynamics.jl",
         doctest = true,
          strict = :doctest)

deploydocs(        repo = "github.com/ClimateFluidPhysics-ANU/MixedLayerThermoclineDynamics.jl.git",
               versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
           push_preview = true
)


text:
using Pkg
packages = [("BSON","0.2.6"),
            ("Calculus","0.5.1"),
            ("CategoricalArrays","0.8.1"),
            #("Clustering","0.14.1"),
            ("Combinatorics","1.0.2"),
            ("CSV","0.6.2"),
            ("DataFrames","0.21.2"),
            #("DataStructures","0.17.18"),
            #("DecisionTree","0.10.3"),
            ("DifferentialEquations","6.14.0"),
            ("Distributions","0.23.4"),
            ("Flux","0.10.4"),
            ("GLM","1.3.9"),
            #("HCubature","1.4.0"),
            ("HypothesisTests","0.10.0"),
            #("HTTP","0.8.15"),
            # ("IJulia",),
            #("Images","0.22.2"),
            #("JSON","0.21.0"),
            # ("Juno",),
            ("KernelDensity","0.5.1"),
            ("LaTeXStrings","1.1.0"),
            #("LIBSVM","0.4.0"),
            #("LightGraphs","1.3.3"),
            ("Measures","0.3.1"),
            #("MLDatasets","0.5.2"),
            ("MultivariateStats","0.7.0"),
            #("NLsolve","4.4.0"),
            ("Plots","1.4.3"),
            ("PyCall","1.91.4"),
            ("PyPlot","2.9.0"),
            ("QuadGK","2.3.1"),
            ("RCall","0.13.7"),
            ("RDatasets","0.6.9"),
            ("Roots","1.0.2"),
            ("SpecialFunctions","0.10.3"),
            ("StatsBase","0.33.0"),
            ("StatsModels","0.6.11"),
            ("StatsPlots","0.14.6"),
            #("TimeSeries","0.18.0")
             ]

for (pN,pV) in packages
    Pkg.add(Pkg.PackageSpec(;name=pN,version=pV))
end

@info "Precompiling packages - this will take a while"
Pkg.precompile()

#@info "Packages not precompiled.\n If you wish to precompile (on Julia V1.3 or above) use `Pkg.precompile()`."

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Base
using Base.Test
using CRF

import CRF.Sequence
import CRF.logsumexp

# Macro hiding the differences between test_throws in Julia v0.2 and v0.3.
macro test_error(error, expression)
    if VERSION < v"0.3.0-prerelease"
        quote
            @test_throws $(expression)
        end
    else
        quote
            @test_throws $(error) $(expression)
        end
    end
end

function feature_vector(yt::Bool, x::Array{Array{Float64,1},1}, t::Int)
    feature_vector(false, yt, x, t)
end

function feature_vector(yp::Bool, yt::Bool, x::Array{Array{Float64,1},1}, t::Int)
    v = Features(4+20+20)

    for ypv in [true, false], ytv in [true, false]
        @append! v ((yp ==ypv) & (yt == ytv))
    end

    for ytv in [true, false], range = 0.0:0.1:0.9
        @append! v ((yt == ytv) & (range <= x[t][1] < (range + 0.1)))
        @append! v ((yt == ytv) & (range <= x[t][2] < (range + 0.1)))
    end

    @test sum(v.x) != 0
    @test v.i == length(v.x) + 1

    return convert(Array{Float64,1}, v)
end

X = [ rand(Float64, 2) for i = 1:100 ]
Y = [ (x[1] < x[2] < 0.5) for x in X ]

s = Sequence(X, Y, feature_vector)

@test s.x == X
@test s.y == Y
@test s.Y == unique(Y)
@test s.f == feature_vector

@test sum(s.F) > 0
@test loglikelihood(s) < 0.0

@test loglikelihood(s) == loglikelihood(s)
@test loglikelihood(s) != loglikelihood(s, Θ=rand(s.k))

@test size(loglikelihood_gradient(s)) == size(s.F)
@test size(label(s)) == size(s.y)

for l in label(s)
    @test l in (true, false)
end

s = Sequence(X, feature_vector; labels=unique(Y))
@test s.x == X
@test s.Y == unique(Y)
@test isempty(s.y)
@test sum(s.F) == 0

# Unlabeld data without label alphabet shouldn't work
@test_error UndefVarError s = Sequence(X, feature_vector)

# Observation / label mismatch shouldn't work
@test_error ArgumentError s = Sequence(X[2:end], Y, feature_vector)

# Feature / weight mismatch shouldn't work
@test_error ArgumentError s = Sequence(X, Y, feature_vector; Θ=rand(2))

# Test logsumexp function
@test logsumexp([0.0:9.0]) == logsumexp([0:9]) == 9.4586297444267107


text:
function turbulence(mtime)

  global Θv,KAm,KAt,γcq,γct,γcm
  ########################################################################
  # Atmospheric vertical diffusion
  ########################################################################
  Θv=ΘA.*(1 .+humid_fac*qA);
  KAm,KAt,γcq,γct,γcm=holtslag(Θv,UA,UA.*0,LH,SH,ustar,ρA[1],ZAF);
  #println("----",maximum(γct),"----")

end


text:
"""
    SimpleSDMPredictor(::Type{EarthEnv}, ::Type{Topography}, layer::Integer=1; resolution::Int64=100, kwargs...)

Allowed resolutions (in km) are 1, 5, 10, 50 (default), and 100

Allowed sources is currently *only* "GMTED"

Allowed aggregations are "mean" (default), "median", "minimum", "maximum", and "std"
"""
function SimpleSDMPredictor(::Type{EarthEnv}, ::Type{Topography}, layer::Integer=1; resolution::Int64=100, source::String="GMTED", aggregation::String="mean", kwargs...)
    @assert resolution in [1, 5, 10, 50, 100]
    @assert source in ["GMTED"]
    @assert aggregation in ["mean", "median", "minimum", "maximum", "std"]
    file = _get_raster(EarthEnv, Topography, layer, resolution, source, aggregation)
    return geotiff(SimpleSDMPredictor, file; kwargs...)
end

function SimpleSDMPredictor(::Type{EarthEnv}, ::Type{Topography}, layers::AbstractArray; kwargs...)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(EarthEnv, Topography, l; kwargs...) for l in layers]
end


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ 7ccb2b9c-1ecd-11eb-205e-11296dd5bb76
using Pkg, DrWatson

# ╔═╡ 855951b0-1ecd-11eb-0108-71f29ba2082f
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using StatsBase, Plots
end

# ╔═╡ 54e90298-1ecd-11eb-01d6-714cdcc0e877
md"## Listing 2.4"

# ╔═╡ 97360536-1ecd-11eb-36d6-038ba8a75dc5
function proportionFished(gF, sF, n, N, withReplacement = false)

    function fishing()
        fishInPond = [ones(Int64,gF); zeros(Int64,sF)]
        fishCaught = Int64[]

        for fish in 1:n
            fished = rand(fishInPond)
            push!(fishCaught,fished)
            if withReplacement == false
                deleteat!(fishInPond, findfirst(x->x==fished, fishInPond))
            end
        end
        sum(fishCaught)
    end

    simulations = [fishing() for _ in 1:N]
    proportions = counts(simulations, 0:n)/N

    if withReplacement
        plot!(0:n, proportions,
             line=:stem, marker=:circle, c=:blue, ms=6, msw=0,
             label="With replacement",
             xlabel="n",
             ylims=(0, 0.6), ylabel="Probability")
    else
        plot!(0:n, proportions,
            line=:stem, marker=:xcross, c=:red, ms=6, msw=0,
            label="Without replacement")
    end
end

# ╔═╡ f505b35c-1ecd-11eb-1eef-6fe09900d4be
begin
	N = 10^6
	goldFish, silverFish, n = 3, 4, 3

	p = plot()
	proportionFished(goldFish, silverFish, n, N)
	proportionFished(goldFish, silverFish, n, N, true)
	p
end

# ╔═╡ f506085a-1ecd-11eb-1e55-db1a82f7bf85
md"## End of listing 2.4"

# ╔═╡ Cell order:
# ╟─54e90298-1ecd-11eb-01d6-714cdcc0e877
# ╠═7ccb2b9c-1ecd-11eb-205e-11296dd5bb76
# ╠═855951b0-1ecd-11eb-0108-71f29ba2082f
# ╠═97360536-1ecd-11eb-36d6-038ba8a75dc5
# ╠═f505b35c-1ecd-11eb-1eef-6fe09900d4be
# ╟─f506085a-1ecd-11eb-1e55-db1a82f7bf85


text:
# measurements.jl
#
# Defines estimators and provides measurements

function sample(qmc_state::BinaryQMCState)
    operator_list = qmc_state.operator_list

    M = length(operator_list) ÷ 2
    spin_prop = copy(qmc_state.left_config)

    @inbounds for op in operator_list[1:M] #propagate half the list only (to the middle)
        if issiteoperator(op) && !isdiagonal(op)
            spin_prop[op[2]] ⊻= 1 #spinflip
        end
    end
    return spin_prop
end

function simulation_cell(qmc_state::BinaryQMCState)
    operator_list = qmc_state.operator_list

    cell = falses(length(qmc_state.left_config), length(operator_list))
    spin_prop = copy(qmc_state.left_config)

    @inbounds for (n, op) in enumerate(operator_list)
        if issiteoperator(op) && !isdiagonal(op)
            spin_prop[op[2]] ⊻= 1 #spinflip
        end
        copy!(cell[:, n], spin_prop)
    end
    return cell
end

magnetization(spin_prop) = mean(x -> 2x - 1, spin_prop)

num_single_site_diag(operator_list) = mean(x -> issiteoperator(x) && isdiagonal(x), operator_list)
num_single_site_offdiag(operator_list) = mean(x -> issiteoperator(x) && !isdiagonal(x), operator_list)
num_single_site(operator_list) = mean(issiteoperator, operator_list)

function autocorrelation(m::Vector)
    N = length(m)

    m′ = m .- mean(m)
    m′ = vcat(m′, zeros(N))
    mw = fft(m′)
    s = abs2.(mw)

    @inbounds chi = real(ifft(s)[1:N])

    @inbounds for i in 1:N
        chi[i] /= (2*N)  # normalize FFT
        chi[i] /= (N - i - 1)
    end
    return chi
end

# use method explained by Sokal to estimate correlation time
# https://pdfs.semanticscholar.org/0bfe/9e3db30605fe2d4d26e1a288a5e2997e7225.pdf
function correlation_time(m::Vector)
    ac = autocorrelation(m)
    ac_0 = ac[1]

    corr_time = 0.0
    @inbounds for M in axes(ac, 1)
        corr_time += (ac[M] / ac_0)
        if M >= 10*corr_time
            break
        end
    end

    return corr_time
end


text:
## -------------  Messages during long calibration

export log_msg

"""
	$(SIGNATURES)

Display a message during calibration. Messages differ in their 'logLevel' (an Integer ≥ 1).
"""
function log_msg(logLevel :: Integer, io,  msgV :: AbstractVector{String})
    top_separator(logLevel, io);
    log_lines(logLevel, io, msgV);
    bottom_separator(logLevel, io);
end

log_msg(logLevel :: Integer, io, msg :: AbstractString) =
    log_msg(logLevel, io, [msg]);

function top_separator(logLevel :: Integer, io)
    if logLevel == 1
        println(io, "\n");
    end
    separator_line(logLevel, io);
end

function separator_line(logLevel :: Integer, io)
    nRepeat = separator_length(logLevel);
    if nRepeat > 0
        println(io, repeat('-', nRepeat));
    end
    return nothing
end

function separator_length(logLevel :: Integer)
    if logLevel == 1
        nRepeat = 40;
    else
        nRepeat = 0;
    end
    return nRepeat
end

function bottom_separator(logLevel :: Integer, io)
    separator_line(logLevel, io);
    if logLevel == 1
        println(io, " ");
    end
end

function log_lines(logLevel :: Integer, io, msgV :: AbstractVector{String})
    indentStr = log_indent(logLevel);
    indentBlank = repeat(' ', length(indentStr));
    for (j, msg) in enumerate(msgV)
        if j == 1
            print(io, indentStr);
        else
            print(io, indentBlank);
        end
        println(io, msg);
    end
end

function log_indent(logLevel :: Integer)
    if logLevel == 1
        return ""
    else
        nSpaces = clamp(2 * (logLevel - 2), 0, 8);
        nDashes = 10 - nSpaces;
        return repeat(' ', nSpaces) * repeat('-', nDashes) * "  ";
    end
end

# ----------


text:
@testset "CIR approx smoothing tests" begin

    Random.seed!(1)

    δ = 3.
    γ = 2.5
    σ = 4.
    Nobs = 2
    dt = 0.011
    Nsteps = 10
    λ = 1.

    α = δ/2
    β = γ/σ^2

    function rec_rcCIR(Dts, x, δ, γ, σ)
        x_new = DualOptimalFiltering.rCIR(1, Dts[1], x[end], δ, γ, σ)
        if length(Dts) == 1
            return Float64[x; x_new]
        else
            return Float64[x; rec_rcCIR(Dts[2:end], x_new, δ, γ, σ)]
        end
    end

    function generate_CIR_trajectory2(times, x0, δ, γ, σ)
        θ1 = δ*σ^2
        θ2 = 2*γ
        θ3 = 2*σ
        Dts = diff(times)
        return rec_rcCIR(Dts, [x0], δ, γ, σ)
    end

    time_grid = [k*dt for k in 0:(Nsteps-1)]
    X = generate_CIR_trajectory2(time_grid, 3, δ*1.2, γ/1.2, σ*0.7)
    Y = map(λ -> rand(Poisson(λ), Nobs), X);
    data = zip(time_grid, Y) |> Dict;

    @test_nowarn DualOptimalFiltering.smooth_CIR_keep_fixed_number(δ, γ, σ, λ, data, 10; silence = false)

    # @test_nowarn DualOptimalFiltering.smooth_CIR_keep_above_threshold(δ, γ, σ, λ, data, 0.01, 10^(-20); silence = false)

    @test_nowarn DualOptimalFiltering.smooth_CIR_keep_fixed_fraction(δ, γ, σ, λ, data, 0.9; silence = false)

end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
function label(colors::AbstractArray)

end


text:
using Documenter, ClearStacktrace

makedocs(;
    modules=[ClearStacktrace],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jkrumbiegel/ClearStacktrace.jl/blob/{commit}{path}#L{line}",
    sitename="ClearStacktrace.jl",
    authors="Julius Krumbiegel",
    assets=String[],
)

deploydocs(;
    repo="github.com/jkrumbiegel/ClearStacktrace.jl",
)


text:
# packaged version of align.jl. Run align_prep.jl or prep_full.jl first.

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

import NLopt
import PlotlyJS
using Plots; plotly()

include("helper.jl")

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

save_BSON_flag = false
save_plot_flag = false

N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

shift_lb = -ones(N_d)
shift_ub = ones(N_d)

Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )

a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)
#

#LS_inds = 1:length(U_cost)
loop_range = 1:length(cost_inds_set)

println("Timing:")
@time obj_funcs, minfs, minxs, rets, ws = NMRCalibrate.alignquantificationcompound(y,
    U_y,
    P_y,
    As,
    Bs,
    Es,
    fs,
    SW,
    Δsys_cs,
    a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
    shift_lb,
    shift_ub,
    cost_inds_set;
    loop_range = loop_range,
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-9,
    maxeval = 50, # 2, # 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    w_lb_default = 1e-1,
    w_ub_default = 100.0,
    β_max_iters = 500, # 2, # 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf);

dummy = 1

#minxs[3][1] = 0.0

#@assert 99==44

if save_BSON_flag
    save_path = joinpath(project_folder, "quantify_results.bson")
    BSON.bson(save_path, region_min_dist = region_min_dist,
    minfs = minfs,
    minxs = minxs,
    rets = rets,
    ws = ws)
end

#### visualize.
# minxs[1][1] = 0.000
# minxs[2][1] = 0.000
# display_reduction_factor = 100
# display_threshold_factor = 0.05/10
display_reduction_factor = 1
display_threshold_factor = 0.001/10
if "L-Isoleucine" in molecule_names
    display_reduction_factor = 1
    display_threshold_factor = 0.001/10
end
display_flag = true

# save data plot.
file_name = "data.html"
plots_save_path = joinpath(project_folder, file_name)

plotquantificationresults(As, Es, ws, project_folder,
    P_y, y,
    region_min_dist,
    obj_funcs, minfs, minxs, rets,
    display_reduction_factor, display_threshold_factor,
    cost_inds_set, loop_range;
    canvas_size = (1000, 400),
    display_flag = display_flag,
    save_plot_flag = save_plot_flag,
    N_viz = 50000)

for r in loop_range
    println("region $(r):")
    println("objective: $(minfs[r]), return status: $(rets)")
    println("shift variable:")
    display(minxs)
    println()
end


text:
function to_xml(doc::Doc2011)
    xml = XMLDocumentNode("1.0")
    ddr = addelement!(xml, "documentoDDR")
    link!(ddr, AttributeNode("cnpj", doc.cnpj))
    link!(ddr, AttributeNode("dataBase", Dates.format(doc.data, "yyyy-mm-dd")))
    link!(ddr, AttributeNode("codigoDocumento", "2011"))
    link!(ddr, AttributeNode("tipoEnvio", encode(doc.tipo)))

    parametros = addelement!(ddr, "parametros")
    add_parametro!(parametros, "31", doc.responsavel.nome)
    add_parametro!(parametros, "32", doc.responsavel.telefone)
    add_parametro!(parametros, "33", doc.responsavel.email)

    contas = addelement!(ddr, "contas")
    for conta in doc.contas
        conta.valor >= 0.01 && add_conta!(contas, conta)
    end

    return xml
end

function add_parametro!(node::EzXML.Node, codigo::String, valor::String)
    par = addelement!(node, "parametro")
    link!(par, AttributeNode("codigoParametro", codigo))
    link!(par, AttributeNode("valorParametro", valor))
    return par
end

function add_conta!(node::EzXML.Node, conta::Conta)
    elm = addelement!(node, "conta")
    link!(elm, AttributeNode("codigoConta", conta.codigo))
    str_valor = @sprintf("%.2f", round(conta.valor, digits=2))
    link!(elm, AttributeNode("valorConta", str_valor))

    if length(conta.detalhes) > 0
        detalhes = addelement!(elm, "detalhamentosDDR")
        for detalhe in conta.detalhes
            detalhe.valor >= 0.01 && add_detalhe!(detalhes, detalhe)
        end
    end

    return elm
end

function add_detalhe!(node::EzXML.Node, detalhe::DetalheConta)
    dd = addelement!(node, "detalhamentoDDR")
    str_valor = @sprintf("%.2f", round(detalhe.valor, digits=2))
    link!(dd, AttributeNode("valorDetalhe", str_valor))

    for el in detalhe.elementos
        node = addelement!(dd, "detalhe")
        link!(node, AttributeNode("codigoElemento", get_codigo_elemento(el)))
        link!(node, AttributeNode("valorElemento", get_valor_elemento(el)))
    end

    return dd
end

write_xml(io::IO, doc::Doc2011) = prettyprint(io, to_xml(doc))
write_xml(path::String, doc::Doc2011) = open(path, "w+") do io write_xml(io, doc) end


text:
#push!(Base.LOAD_PATH,joinpath(dirname(@__FILE__),"..",".."))

using Xpress, Compat.Test, MathOptInterface, MathOptInterface.Test

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test

@testset "MathOptInterface" begin
    @testset "Unit Tests" begin
        config = MOIT.TestConfig()
        solver = Xpress.Optimizer()
        MOIT.basic_constraint_tests(solver, config;
            exclude = [
                (MOI.SingleVariable, MOI.Integer),
                (MOI.SingleVariable, MOI.EqualTo{Float64}),
                (MOI.SingleVariable, MOI.Interval{Float64}),
                (MOI.SingleVariable, MOI.GreaterThan{Float64}),
                (MOI.SingleVariable, MOI.LessThan{Float64}),
            ]
        )
        MOIT.unittest(solver, config, [
            "solve_qp_edge_cases",    # tested below
            "solve_qcp_edge_cases"    # tested below
        ])
        @testset "solve_qp_edge_cases" begin
                MOIT.solve_qp_edge_cases(solver,
                MOIT.TestConfig(atol=1e-2,rtol = 1e-2)
            )
        end
        @testset "solve_qcp_edge_cases" begin
                MOIT.solve_qcp_edge_cases(solver,
                MOIT.TestConfig(atol=1e-2,rtol = 1e-2)
            )
        end
    end

    @testset "Linear tests" begin
        linconfig = MOIT.TestConfig()
        solver = Xpress.Optimizer()
        MOIT.contlineartest(solver, linconfig, [
            "linear12","linear8a","linear8b","linear8c",
            # VariablePrimalStart not implemented.
            "partial_start"
        ])

        solver_nopresolve = Xpress.Optimizer(PRESOLVE = 0)
        MOIT.contlineartest(solver_nopresolve, linconfig, [
            "linear12",
            # VariablePrimalStart not implemented.
            "partial_start"
        ])

        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        MOIT.linear12test(solver, linconfig_nocertificate)
    end

    @testset "Quadratic tests" begin
        quadconfig = MOIT.TestConfig(atol=1e-5, rtol=1e-5, duals=false)
        solver = Xpress.Optimizer()
        MOIT.contquadratictest(solver, quadconfig)
    end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig()
        solver = Xpress.Optimizer()
        MOIT.lintest(solver, linconfig, ["lin3","lin4"])

        solver_nopresolve = Xpress.Optimizer(PRESOLVE=0)
        MOIT.lintest(solver_nopresolve, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig()
        solver = Xpress.Optimizer()
        MOIT.intlineartest(solver, intconfig, ["int2"])
    end

    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig()
        solver = Xpress.Optimizer()
        MOIT.nametest(solver)
        @testset "validtest" begin
            MOIT.validtest(solver)
        end
        @testset "emptytest" begin
            MOIT.emptytest(solver)
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(solver)
        end
        @testset "copytest" begin
            solver2 = Xpress.Optimizer()
            MOIT.copytest(solver,solver2)
        end
    end

    @testset "IIS tests" begin
        @testset "Variable bounds (ScalarAffine)" begin
            model = Xpress.Optimizer()
            x = MOI.add_variable(model)
            c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(2.0))
            c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(1.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
            @test_throws ErrorException MOI.get(model, Xpress.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            Xpress.compute_conflict(model)
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMAL
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c1) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c2) == true
        end

        @testset "Two conflicting constraints (GreaterThan, LessThan)" begin
            model = Xpress.Optimizer()
            x = MOI.add_variable(model)
            y = MOI.add_variable(model)
            b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
            b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
            cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
            c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
            cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
            c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
            @test_throws ErrorException MOI.get(model, Xpress.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            # Two possible IISes: b1, b2, c1 OR b2, c1, c2
            Xpress.compute_conflict(model)
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMAL
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b1) in [true, false]
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b2) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c1) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c2) in [true, false]
        end

        @testset "Two conflicting constraints (EqualTo)" begin
            model = Xpress.Optimizer()
            x = MOI.add_variable(model)
            y = MOI.add_variable(model)
            b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
            b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
            cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
            c1 = MOI.add_constraint(model, cf1, MOI.EqualTo(-1.0))
            cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
            c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
            @test_throws ErrorException MOI.get(model, Xpress.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            # Two possible IISes: b1, b2, c1 OR b2, c1, c2
            Xpress.compute_conflict(model)
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMAL
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b1) in [true, false]
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b2) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c1) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c2) in [true, false]
        end

        @testset "Variables outside conflict" begin
            model = Xpress.Optimizer()
            x = MOI.add_variable(model)
            y = MOI.add_variable(model)
            z = MOI.add_variable(model)
            b1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(0.0))
            b2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [y]), 0.0), MOI.GreaterThan(0.0))
            b3 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [z]), 0.0), MOI.GreaterThan(0.0))
            cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
            c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
            cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 1.0], [x, y, z]), 0.0)
            c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
            @test_throws ErrorException MOI.get(model, Xpress.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            Xpress.compute_conflict(model)
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMAL
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b1) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b2) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), b3) == false
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c1) == true
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c2) == false
        end

        @testset "No conflict" begin
            model = Xpress.Optimizer()
            x = MOI.add_variable(model)
            c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(1.0))
            c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(2.0))

            # Getting the results before the conflict refiner has been called must return an error.
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.OPTIMIZE_NOT_CALLED
            @test_throws ErrorException MOI.get(model, Xpress.ConstraintConflictStatus(), c1)

            # Once it's called, no problem.
            Xpress.compute_conflict(model)
            @test MOI.get(model, Xpress.ConflictStatus()) == MOI.INFEASIBLE
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c1) == false
            @test MOI.get(model, Xpress.ConstraintConflictStatus(), c2) == false
        end
    end
end


text:
# models for plots
# to do
function plot_timecoor(file::String, fit::EasyFit.SingleExponential)

end


text:
using Gadfly, LightGraphs, GraphPlot

gx = kronecker(5, 12, 0.57, 0.19, 0.19)
gplot(gx)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function getmembers(s;kwargs...)
	dict = Dict{Symbol, Any}()
	for symb in fieldnames(typeof(s))
		dict[symb] = getfield(s, symb)
	end
	return dict
end
function modify(s;kwargs...)
	dict = getmembers(s)
	for (k, v) in kwargs
		if haskey(dict, k)
			dict[k] = v
		end
	end
	return typeof(s)(;dict...)
end
keytosymbol(x) = Dict(Symbol(k) => v for (k, v) in pairs(x))
function convertdict(S, fname::String)
	dict = JSON.parsefile(fname)
	res = Dict()
	for (k, v) in pairs(dict)
		T = fieldtype(S, Symbol(k))
		if T != String
			res[k] = eltype(T).(v)
		else
			res[k] = (T).(v)
		end
	end
	return S(InputFile=fname;keytosymbol(res)...)
end


text:
################################################################################
#
# Routines for the the evaluation and gradient computation of FWI
#
################################################################################

export objFWI, objFWI!, objFWI_rawinput

### Optim.jl-adapted wrapper

function objFWI!(
	F, Gm, m::Array{R, 2},
	n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector,
	dat::judiVector;
	opt = Options(space_order=SPACE_ORDER),
	dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	if Gm == nothing
		mode = "eval"
	else
		mode = "grad"
	end
	argout = objFWI(m, n, d, o, wav, dat; opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	if Gm != nothing
		fval = argout[1]
		Gm .= argout[2]
	else
		fval = argout
	end
	return fval

end

function objFWI!(
	F, Gm, m::Array{R, 2},
	model::Modelall,
	wav::judiVector,
	dat::judiVector;
	opt = Options(space_order=SPACE_ORDER),
	dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	if Gm == nothing
		mode = "eval"
	else
		mode = "grad"
	end
	model.m = m
	argout = objFWI(m, model, wav, dat; opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
	if Gm != nothing
		fval = argout[1]
		Gm .= argout[2]
	else
		fval = argout
	end
	return fval

end

### Objective wrapper for JUDI data types

function objFWI(
	m::Modelall,
	wav::judiVector,
	dat::judiVector;
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)
	return objFWI_rawinput(m, wav.geometry, dat.geometry, wav.data, dat.data, 1:length(dat.data); opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)
end

### Objective wrapper for array input

function objFWI(
	m::Array{R, 2}, n::NTuple{2, Int64}, d::NTuple{2, R}, o::NTuple{2, R},
	wav::judiVector, dat::judiVector;
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	return objFWI(Model(n, d, o, m), wav, dat; opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

function objFWI(
	m::Array{R, 2}, model::Modelall,
	wav::judiVector, dat::judiVector;
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)
	model.m = m
	return objFWI(model, wav, dat; opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

### Objective: serial/parallel implementation

## Serial implementation

function objFWI_rawinput(
	model_full::Modelall,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2}, dat_data::Array{R, 2};
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, R} = nothing,
	Filter::Union{Nothing, Array{C, 1}} = nothing,
	gradmprec_fun = g->g)

	# Setting pre-defined absorbing layer size
	model_full.nb = NB

	# Load full geometry for out-of-core geometry containers
	typeof(rcv_geom) == GeometryOOC && (rcv_geom = Geometry(rcv_geom))
	typeof(src_geom) == GeometryOOC && (src_geom = Geometry(src_geom))
	length(model_full.n) == 3 ? dims = [3, 2, 1] : dims = [2, 1] # model dimensions for Python are (z,y,x) and (z,x)

	# Limit model to area with sources/receivers
	if opt.limit_m == true
		model = deepcopy(model_full)
		model = limit_model_to_receiver_area(src_geom, rcv_geom, model, opt.buffer_size)
	else
		model = model_full
	end

	# Set up Python model structure
	model_py = devito_model(model, "F", 1, opt, 0)

	# Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
	rcv_geom = remove_out_of_bounds_receivers(rcv_geom, model)

	# Call to objective with Julia/Devito interface function
	return objFWI_jldevito(model_py, model.o, src_geom, rcv_geom, src_data, dat_data; opt = opt, mode = mode, dt_comp = dt_comp, Filter = Filter, gradmprec_fun = gradmprec_fun)

end

## Parallel implementation

function objFWI_rawinput(
	model::Modelall,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data, dat_data,
	src_idx::UnitRange{Int64};
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, Array{Any, 1}} = nothing,
	Filter::Union{Nothing, Array{C, 1}, Array{Array{C, 1}, 1}} = nothing,
	gradmprec_fun = g->g)

	# Setup parallelization
	p = default_worker_pool()
	objFWI_rawinput_par = remote(TWRIdual.objFWI_rawinput)
	objFWI_rawinput = retry(objFWI_rawinput_par)

	# Initialize output
	nsrc = length(src_idx)
	results = Array{Any}(undef, nsrc)

	# Process shots from source channel asynchronously
	@sync begin
		for j = 1:nsrc

			# Local geometry for current position
			opt_loc = subsample(opt, j)
			src_geom_loc = subsample(src_geom, j)
			rcv_geom_loc = subsample(rcv_geom, j)

			# Selecting variables for current shot index
			src_data_loc = src_data[j]
			dat_data_loc = dat_data[j]
			if isa(Filter, Array{Array{C, 1}, 1})
				Filter_loc = Filter[j]
			else
				Filter_loc = Filter
			end
			if dt_comp == nothing
				dt_comp_loc = nothing
			else
				dt_comp_loc = dt_comp[j]
			end

			# Local result
			results[j] = @spawn objFWI_rawinput(model, src_geom_loc, rcv_geom_loc, src_data_loc, dat_data_loc, opt = opt_loc; mode = mode, dt_comp = dt_comp_loc, Filter = Filter_loc, gradmprec_fun = gradmprec_fun)

		end
	end

	# Aggregating results
	if mode == "eval"
		argout = fetch(results[1])
		for j = 2:nsrc
			argout += fetch(results[j])
		end
	elseif mode == "grad"
		arg = fetch(results[1])
		argout1 = arg[1]
		argout2 = arg[2]
		for j = 2:nsrc
			arg = fetch(results[j])
			argout1 += arg[1]
			argout2 += arg[2]
		end
		argout = (argout1, argout2)
	end
	return argout./nsrc

end

### Julia/devito call

function objFWI_jldevito(
	model_py::PyCall.PyObject,
	origin,
	src_geom::Geometry, rcv_geom::Geometry,
	src_data::Array{R, 2}, dat_data::Array{R, 2};
	opt = Options(space_order=SPACE_ORDER),
	mode = "eval", dt_comp::Union{Nothing, R} = nothing,
	Filter::Union{Nothing, Array{C, 1}} = nothing,
	gradmprec_fun = g->g)
	# Interpolate input data to computational grid
	if dt_comp == nothing
		dt_comp = R(model_py.critical_dt)
	end
	q_in = time_resample(src_data, src_geom, dt_comp)[1]
	dat_in = time_resample(dat_data, rcv_geom, dt_comp)[1]
	nt_comp = size(q_in, 1)
	nt_rcv = Int(trunc(rcv_geom.t[1]/dt_comp+1))

	# Set up coordinates with devito dimensions
	src_coords = setup_grid(src_geom, model_py.shape)
	rcv_coords = setup_grid(rcv_geom, model_py.shape)

	fwi_ops = load_codegen_FWI(model_py)
	# Computing output
	if mode == "eval"
		argout = pycall(fwi_ops.fwi_fun, R, src_coords, rcv_coords, q_in, dat_in, Filter, mode)
	elseif mode == "grad"
		argout = pycall(fwi_ops.fwi_fun, Tuple{R, Array{R, 2}}, src_coords, rcv_coords, q_in, dat_in, Filter, mode)
	end

	# Output post-processing
	if mode == "grad"
		gm = remove_padding(argout[2], model_py.nbl, true_adjoint = opt.sum_padding) # remove AB layers
		argout = (argout[1], gradmprec_fun(gm))
	end

	return argout

end

# Loader for python devito modeling
@cache function load_codegen_FWI(model)

	pushfirst!(PyVector(pyimport("sys")."path"), string(pwd(), "/src/TWRI_py/"))
	int = pyimport("wri_interface")
	return int.WRIInterface(model, SPACE_ORDER)
end


text:
using GitHub, MassInstallAction
import HTTP

# This demo adds or updates workflow files (TagBot, Documenter, and CI) for a set of repositories
# in an organization. In general you should copy/paste these lines after editing them
# to suit your needs. Note that you'll need to specify your GitHub personal access token,
# specify the organization/user, and select repositories interactively.

# authenticate
auth = authenticate(#= your-github-personal-access-token =#)   # don't ever post this publicly!

# Get all the repositories in the JuliaDebug organization
rs = repos(#= organization/user goes here =#, true; auth=auth)[1]   # change true -> false for users (true = "org")

# Optional: filter on various properties, e.g.,
# filter!(r -> !r.fork, rs)    # eliminate all repos that are forks of other repos

# Manually pick a subset of the repositories, using an interactive menu
using REPL.TerminalMenus
menu = MultiSelectMenu([rp.name for rp in rs])
sel = request(menu)   # hit arrows to move through list, Enter to select, and 'd' when done
rs = rs[collect(sel)]

## Key/secret creation

# Create a `DOCUMENTER_KEY` secret and corresponding `Documenter` deploy key, unless the repository
# already has a secret of that name.
# This key is used by TagBot as well as Documenter, so do this even if you don't have any `docs/` in your repository.
for r in rs
    has_secret = false
    for s in secrets(r; auth=auth)[1]
        if s.name == "DOCUMENTER_KEY"
            has_secret = true
            break
        end
    end
    has_secret && continue
    pubkey, privkey = GitHub.genkeys()
    create_deploykey(r; auth=auth, params=Dict("key"=>pubkey, "title"=>"Documenter", "read_only"=>false))
    create_secret(r, "DOCUMENTER_KEY"; auth=auth, value=privkey)
end

## TagBot

# Add or modify the TagBot workflow to use `DOCUMENTER_KEY`. Here we use the templates in MassInstallAction directly,
# but alternatively we could create a workflow script here and use `Workflow("TagBot", "TagBot.yml" => script)`.
#
# Related: https://discourse.julialang.org/t/ann-required-updates-to-tagbot-yml/49249
workflow = MassInstallAction.tag_bot()
for r in rs
    MassInstallAction.install(workflow, r; auth=auth, commit_message="Set up TagBot workflow")
end

## Documenter

script = """
name: Documenter
on:
  push:
    branches: [master]
    tags: [v*]
  pull_request:

jobs:
  Documenter:
    name: Documentation
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v2
      - uses: julia-actions/julia-buildpkg@latest
      - uses: julia-actions/julia-docdeploy@latest
        env:
          GITHUB_TOKEN: \${{ secrets.GITHUB_TOKEN }}
          DOCUMENTER_KEY: \${{ secrets.DOCUMENTER_KEY }}
"""
workflow = MassInstallAction.Workflow("Documenter", "Documenter.yml" => script)
for r in rs
    # Check whether the repo has docs, if not skip it
    filesdirs = directory(r, "."; auth=auth)[1]
    idx = findfirst(fd -> fd.typ == "dir" && fd.name == "docs", filesdirs)
    if idx !== nothing
        MassInstallAction.install(workflow, r; auth=auth, commit_message="Build docs on GitHub Actions")
    end
end

## CI

# Default CI workflow that tests on the same versions as the current .travis.yml

ci_pre = """
name: CI
on:
  pull_request:
  push:
    branches:
      - master
    tags: '*'
jobs:
  test:
    name: Julia \${{ matrix.version }} - \${{ matrix.os }} - \${{ matrix.arch }} - \${{ github.event_name }}
    runs-on: \${{ matrix.os }}
    strategy:
      fail-fast: false
      matrix:
        version:
"""
ci_post = """
        os:
          - ubuntu-latest
          - macOS-latest
          - windows-latest
        arch:
          - x64
    steps:
      - uses: actions/checkout@v2
      - uses: julia-actions/setup-julia@v1
        with:
          version: \${{ matrix.version }}
          arch: \${{ matrix.arch }}
      - uses: actions/cache@v1
        env:
          cache-name: cache-artifacts
        with:
          path: ~/.julia/artifacts
          key: \${{ runner.os }}-test-\${{ env.cache-name }}-\${{ hashFiles('**/Project.toml') }}
          restore-keys: |
            \${{ runner.os }}-test-\${{ env.cache-name }}-
            \${{ runner.os }}-test-
            \${{ runner.os }}-
      - uses: julia-actions/julia-buildpkg@v1
      - uses: julia-actions/julia-runtest@v1
      - uses: julia-actions/julia-processcoverage@v1
      - uses: codecov/codecov-action@v1
        with:
          file: lcov.info
"""
for r in rs
    # Check whether the repo has a workflow called "CI.yml", and if so skip it
    try
        filesdirs = directory(r, ".github/workflows"; auth=auth)[1]
        idx = findfirst(fd -> lowercase(fd.name) ∈ ("ci.yml", "ci.yaml"), filesdirs)
        if idx !== nothing
            println("\n\n$(r.name) skipped due to presence of $(filesdirs[idx].name)")
            continue
        end
    catch
    end
    try
        # Parse the .travis.yml file and extract the versions
        url = file(r, ".travis.yml"; auth=auth).download_url
        tscript = String(HTTP.get(url).body)
        lines = split(chomp(tscript), '\n')
        idx = findfirst(isequal("julia:"), lines)
        vs = String[]
        while idx < length(lines)
            idx += 1
            m = match(r"^\s*-\s*(nightly|[\d\.]*)", lines[idx])
            m === nothing && break
            push!(vs, m.captures[1])
        end
        # optional: delete "nightly"
        if !isempty(vs)
            # Build the replacement CI.yml
            script = ci_pre
            for v in vs
                line = " "^10 * "- '" * v * "'\n"
                # optional: comment out nightly
                if v == "nightly"
                    line = "#" * line
                end
                script *= line
            end
            script *= ci_post
            todelete = ["../../.travis.yml"]
            # Find out if there's an appveyor and set it up for deletion
            filesdirs = directory(r, "."; auth=auth)[1]
            idx = findfirst(fd -> fd.typ == "file" && fd.name ∈ ("appveyor.yml", ".appveyor.yml"), filesdirs)
            if idx !== nothing
                push!(todelete, "../../" * filesdirs[idx].name)
            end
            workflow = MassInstallAction.Workflow("CI", "CI.yml" => script, todelete)
            println("\n\n$(r.name) had .travis.yml script:\n", tscript, "\nInstalling the following CI.yml script:\n", script)
            println("OK to submit pull request? (y/n)")
            resp = readline(; keep=false)
            if resp ∈ ("y", "Y", "yes", "Yes", "YES")
                MassInstallAction.install(workflow, r; auth=auth, commit_message="Switch to GitHub Actions for CI")
            end
        end
    catch
        println("\n\nNo .travis.yml found for $(r.name)")
    end
end


text:
module MutationFunctionsModule

import ..CoreModule: CONST_TYPE, Node, copy_node, Options
import ..EquationUtilsModule: count_nodes, count_constants, count_operators, count_depth

# Return a random node from the tree
function random_node(tree::Node)::Node
    if tree.degree == 0
        return tree
    end
    b = 0
    c = 0
    if tree.degree >= 1
        b = count_nodes(tree.l)
    end
    if tree.degree == 2
        c = count_nodes(tree.r)
    end

    i = rand(1:(1 + b + c))
    if i <= b
        return random_node(tree.l)
    elseif i == b + 1
        return tree
    end

    return random_node(tree.r)
end

# Randomly convert an operator into another one (binary->binary;
# unary->unary)
function mutate_operator(tree::Node, options::Options)::Node
    if count_operators(tree) == 0
        return tree
    end
    node = random_node(tree)
    while node.degree == 0
        node = random_node(tree)
    end
    if node.degree == 1
        node.op = rand(1:(options.nuna))
    else
        node.op = rand(1:(options.nbin))
    end
    return tree
end

# Randomly perturb a constant
function mutate_constant(tree::Node, temperature::T, options::Options)::Node where {T<:Real}
    # T is between 0 and 1.

    if count_constants(tree) == 0
        return tree
    end
    node = random_node(tree)
    while node.degree != 0 || node.constant == false
        node = random_node(tree)
    end

    bottom = convert(T, 1//10)
    maxChange = options.perturbationFactor * temperature + convert(T, 1) + bottom
    factor = maxChange^Float32(rand())
    makeConstBigger = rand() > 0.5

    if makeConstBigger
        node.val *= convert(CONST_TYPE, factor)
    else
        node.val /= convert(CONST_TYPE, factor)
    end

    if rand() > options.probNegate
        node.val *= -1
    end

    return tree
end

# Add a random unary/binary operation to the end of a tree
function append_random_op(
    tree::Node, options::Options, nfeatures::Int; makeNewBinOp::Union{Bool,Nothing}=nothing
)::Node
    node = random_node(tree)
    while node.degree != 0
        node = random_node(tree)
    end

    if makeNewBinOp === nothing
        choice = rand()
        makeNewBinOp = choice < options.nbin / (options.nuna + options.nbin)
    end

    if makeNewBinOp
        newnode = Node(
            rand(1:(options.nbin)), make_random_leaf(nfeatures), make_random_leaf(nfeatures)
        )
    else
        newnode = Node(rand(1:(options.nuna)), make_random_leaf(nfeatures))
    end

    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant

    return tree
end

# Insert random node
function insert_random_op(tree::Node, options::Options, nfeatures::Int)::Node
    node = random_node(tree)
    choice = rand()
    makeNewBinOp = choice < options.nbin / (options.nuna + options.nbin)
    left = copy_node(node)

    if makeNewBinOp
        right = make_random_leaf(nfeatures)
        newnode = Node(rand(1:(options.nbin)), left, right)
    else
        newnode = Node(rand(1:(options.nuna)), left)
    end
    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant
    return tree
end

# Add random node to the top of a tree
function prepend_random_op(tree::Node, options::Options, nfeatures::Int)::Node
    node = tree
    choice = rand()
    makeNewBinOp = choice < options.nbin / (options.nuna + options.nbin)
    left = copy_node(tree)

    if makeNewBinOp
        right = make_random_leaf(nfeatures)
        newnode = Node(rand(1:(options.nbin)), left, right)
    else
        newnode = Node(rand(1:(options.nuna)), left)
    end
    if newnode.degree == 2
        node.r = newnode.r
    end
    node.l = newnode.l
    node.op = newnode.op
    node.degree = newnode.degree
    node.val = newnode.val
    node.feature = newnode.feature
    node.constant = newnode.constant
    return node
end

function make_random_leaf(nfeatures::Int)::Node
    if rand() > 0.5
        return Node(randn(CONST_TYPE))
    else
        return Node(rand(1:nfeatures))
    end
end

# Return a random node from the tree with parent, and side ('n' for no parent)
function random_node_and_parent(
    tree::Node, parent::Union{Node,Nothing}; side::Char
)::Tuple{Node,Union{Node,Nothing},Char}
    if tree.degree == 0
        return tree, parent, side
    end
    b = 0
    c = 0
    if tree.degree >= 1
        b = count_nodes(tree.l)
    end
    if tree.degree == 2
        c = count_nodes(tree.r)
    end

    i = rand(1:(1 + b + c))
    if i <= b
        return random_node_and_parent(tree.l, tree; side='l')
    elseif i == b + 1
        return tree, parent, side
    end

    return random_node_and_parent(tree.r, tree; side='r')
end

function random_node_and_parent(tree::Node)::Tuple{Node,Union{Node,Nothing},Char}
    return random_node_and_parent(tree, nothing; side='n')
end

# Select a random node, and replace it an the subtree
# with a variable or constant
function delete_random_op(tree::Node, options::Options, nfeatures::Int)::Node
    node, parent, side = random_node_and_parent(tree)
    isroot = (parent === nothing)

    if node.degree == 0
        # Replace with new constant
        newnode = make_random_leaf(nfeatures)
        node.degree = newnode.degree
        node.val = newnode.val
        node.constant = newnode.constant
        if !newnode.constant
            node.feature = newnode.feature
        end
    elseif node.degree == 1
        # Join one of the children with the parent
        if isroot
            return node.l
        elseif parent.l == node
            parent.l = node.l
        else
            parent.r = node.l
        end
    else
        # Join one of the children with the parent
        if rand() < 0.5
            if isroot
                return node.l
            elseif parent.l == node
                parent.l = node.l
            else
                parent.r = node.l
            end
        else
            if isroot
                return node.r
            elseif parent.l == node
                parent.l = node.r
            else
                parent.r = node.r
            end
        end
    end
    return tree
end

# Create a random equation by appending random operators
function gen_random_tree(length::Int, options::Options, nfeatures::Int)::Node
    # Note that this base tree is just a placeholder; it will be replaced.
    tree = Node(convert(CONST_TYPE, 1))
    for i in 1:length
        # TODO: This can be larger number of nodes than length.
        tree = append_random_op(tree, options, nfeatures)
    end
    return tree
end

function gen_random_tree_fixed_size(node_count::Int, options::Options, nfeatures::Int)::Node
    tree = make_random_leaf(nfeatures)
    cur_size = count_nodes(tree)
    while cur_size < node_count
        if cur_size == node_count - 1  # only unary operator allowed.
            options.nuna == 0 && break # We will go over the requested amount, so we must break.
            tree = append_random_op(tree, options, nfeatures; makeNewBinOp=false)
        else
            tree = append_random_op(tree, options, nfeatures)
        end
        cur_size = count_nodes(tree)
    end
    return tree
end

"""Crossover between two expressions"""
function crossover_trees(tree1::Node, tree2::Node)::Tuple{Node,Node}
    tree1 = copy_node(tree1)
    tree2 = copy_node(tree2)

    node1, parent1, side1 = random_node_and_parent(tree1)
    node2, parent2, side2 = random_node_and_parent(tree2)

    node1 = copy_node(node1)

    if side1 == 'l'
        parent1.l = copy_node(node2)
        # tree1 now contains this.
    elseif side1 == 'r'
        parent1.r = copy_node(node2)
        # tree1 now contains this.
    else # 'n'
        # This means that there is no parent2.
        tree1 = copy_node(node2)
    end

    if side2 == 'l'
        parent2.l = node1
    elseif side2 == 'r'
        parent2.r = node1
    else # 'n'
        tree2 = node1
    end
    return tree1, tree2
end

end


text:
    # error?
function Base.getindex(x::ParseMatch{<:Union{Transformation}},a...)
    println("implement lazy getindex(::$(typeof(x))")
    getindex(get(x), a...)
end

getstateindex(state::Integer, a...) = MatchState()

function Base.getindex(x::ParseMatch{<:Sequence}, i)
    mi = x.parsings
    sequence = mi.sequence
    inner =
    outer_state = x.state
    inner_state = getindex(outer_state, i)
    l = lastindex(mi.parser.parts)
    start = stop = x.stop
    ## walk back
    while l>=i
        stop = start
        start = leftof(sequence,stop,mi.parser.parts[l],outer_state[l]) ##state[end][1]
        l-=1
    end
    pm = ParseMatch(MatchesIterator(inner,mi.sequence, mi.from, mi.till),
                    start, stop, outer_state[i])
end

function Base.getindex(x::ParseMatch{<:Repeat}, i)
    mi = x.parsings
    sequence = mi.sequence
    inner = mi.parser.parser
    outer_state = x.state
    l = state_length(mi.parser,outer_state)
    local inner_state::state_type(inner)
    start = stop = x.stop
    ## walk back
    while l>=i
        stop = start
        # outer_state::Union{<:Vector, Int}
        inner_state, outer_state = poplast!(outer_state)
        start = leftof(sequence,stop,inner,inner_state) ##state[end][1]
        l-=1
    end
    pm = ParseMatch(MatchesIterator(inner,mi.sequence, mi.from, mi.till),
                    start, stop, inner_state)
end

function Base.getindex(x::ParseMatch{<:MatchedSubSequence},a...)
    getindex(x.parsings.sequence, ( nextind(x.parsings.sequence,x.start,i-1) for i in a)...)
end

function Base.getindex(x::ParseMatch{<:WrappedParser},a...)
    mi = x.parsings
    pm = ParseMatch(MatchesIterator(mi.parser.parser,mi.sequence,mi.from, mi.till),
                    x.start,x.stop,x.state)
    getindex(pm,a...)
end


text:
# This file is part of Kpax3. License is MIT.

function test_support_mcmc_constructor()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"
  maxclust = 1

  settings = Kpax3.KSettings(ifile, ofile, maxclust=maxclust)

  x = Kpax3.AminoAcidData(settings)

  (m, n) = size(x.data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = Kpax3.AminoAcidState(x.data, [1; 1; 1; 2; 2; 3], priorR, priorC, settings)

  g = 0
  lp = zeros(Float64, 4, 3, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end

  support = Kpax3.MCMCSupport(state, priorC)

  @test support.m == m
  @test support.n == n

  @test support.u == Int[a for a in 1:n]
  @test support.t == zeros(Float64, n)
  @test support.lp == lp
  @test support.lq == zeros(Float64, 4, 3, m)
  @test support.lr == zeros(Float64, 3, 3, m)
  @test support.vi == 0
  @test support.ni == zeros(Float64, m)
  @test support.ui == zeros(Int, n)
  @test isa(support.wi, Kpax3.KWeight)
  @test support.wi.c == zeros(Float64, m)
  @test support.wi.w == zeros(Float64, 4, m)
  @test support.wi.z == zeros(Float64, 4, m)
  @test support.lpi == zeros(Float64, 4, m)

  @test support.vj == 0
  @test support.nj == zeros(Float64, m)
  @test support.uj == zeros(Int, n)
  @test isa(support.wj, Kpax3.KWeight)
  @test support.wj.c == zeros(Float64, m)
  @test support.wj.w == zeros(Float64, 4, m)
  @test support.wj.z == zeros(Float64, 4, m)
  @test support.lpj == zeros(Float64, 4, m)

  @test support.tmp == zeros(Float64, 4)

  @test support.cl == zeros(Int, n)
  @test support.k == 0

  @test support.lograR == 0.0

  logmlik = 0.0
  tmp = zeros(Float64, 3)
  for b in 1:m
    tmp[1] = priorC.logγ[1]
    tmp[2] = priorC.logγ[2]
    tmp[3] = priorC.logγ[3]

    for l in 1:state.k
      g = state.cl[l]

      tmp[1] += lp[1, g, b]
      tmp[2] += lp[2, g, b]
      tmp[3] += log(exp(priorC.logω[state.k][1] + lp[3, g, b]) + exp(priorC.logω[state.k][2] + lp[4, g, b]))
    end

    logmlik += log(exp(tmp[1]) + exp(tmp[2]) + exp(tmp[3]))
  end

  @test isapprox(support.logmlik, logmlik, atol=ε)

  @test support.logmlikcandidate == 0.0

  nothing
end

test_support_mcmc_constructor()

function test_support_mcmc_resize()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"
  maxclust = 1

  settings = Kpax3.KSettings(ifile, ofile, maxclust=maxclust)

  x = Kpax3.AminoAcidData(settings)

  (m, n) = size(x.data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = Kpax3.AminoAcidState(x.data, [1; 1; 1; 2; 2; 3], priorR, priorC, settings)

  g = 0
  lp = zeros(Float64, 4, 3, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end

  support = Kpax3.MCMCSupport(state, priorC)

  # do nothing
  Kpax3.resizesupport!(support, 3)

  @test support.lp == lp
  @test support.lq == zeros(Float64, 4, 3, m)
  @test support.lr == zeros(Float64, 3, 3, m)

  Kpax3.resizesupport!(support, 6)

  g = 0
  lp = zeros(Float64, 4, 6, m)
  for b in 1:m, l in 1:state.k
    g = state.cl[l]

    lp[1, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[1, b], priorC.B[1, b])
    lp[2, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[2, b], priorC.B[2, b])
    lp[3, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[3, b], priorC.B[3, b])
    lp[4, g, b] = Kpax3.logmarglik(state.n1s[g, b], state.v[g], priorC.A[4, b], priorC.B[4, b])
  end

  @test support.lp == lp
  @test support.lq == zeros(Float64, 4, 6, m)
  @test support.lr == zeros(Float64, 3, 6, m)

  nothing
end

test_support_mcmc_resize()

function test_support_ga_constructor()
  m = 18
  n = 6

  support = Kpax3.GASupport(m, n)

  @test support.m == m
  @test support.n == n

  @test isa(support.oi, Kpax3.KOffspring)
  @test support.oi.R == zeros(Int, n)
  @test support.oi.v == zeros(Int, n)

  @test isa(support.oj, Kpax3.KOffspring)
  @test support.oj.R == zeros(Int, n)
  @test support.oj.v == zeros(Int, n)

  nothing
end

test_support_ga_constructor()


text:
@enum MetisScheme begin
    metis_scheme_recursive_bisection = Metis.METIS_PTYPE_RB
    metis_scheme_direct_k_way = Metis.METIS_PTYPE_KWAY
end

@enum MetisObjective begin
    metis_objective_edge_cut = Metis.METIS_OBJTYPE_CUT
    metis_objective_communication_volume = Metis.METIS_OBJTYPE_VOL
end

@enum MetisCoarsen begin
    metis_coarsen_random = Metis.METIS_CTYPE_RM
    metis_coarsen_sorted_heavy_edge = Metis.METIS_CTYPE_SHEM
end

@enum MetisInitialize begin
    metis_initialize_grow = Metis.METIS_IPTYPE_GROW
    metis_initialize_random = Metis.METIS_IPTYPE_RANDOM
    metis_initialize_edge = Metis.METIS_IPTYPE_EDGE
    metis_initialize_node = Metis.METIS_IPTYPE_NODE
end

@enum MetisRefine begin
    metis_refine_fm = Metis.METIS_RTYPE_FM
    metis_refine_greedy = Metis.METIS_RTYPE_GREEDY
    metis_refine_fm_two_sided = Metis.METIS_RTYPE_SEP2SIDED
    metis_refine_fm_one_sided = Metis.METIS_RTYPE_SEP1SIDED
end

struct MetisPartitioner
    assumesymmetry::Bool
    checksymmetry::Bool
    scheme::Union{Nothing, MetisScheme}
    objective::Union{Nothing, MetisObjective}
    coarsen::Union{Nothing, MetisCoarsen}
    initialize::Union{Nothing, MetisInitialize}
    refine::Union{Nothing, MetisRefine}
    ncuts::Union{Nothing, Int}
    niterations::Union{Nothing, Int}
    seed::Union{Nothing, Int}
    minimizeconnectivity::Union{Nothing, Bool}
    maketwohops::Union{Nothing, Bool}
    forcecontiguous::Union{Nothing, Bool}
    compressgraph::Union{Nothing, Bool}
    imbalance::Union{Nothing, Float64}
    weight::Any
    verbose_info::Bool
    verbose_time::Bool
    verbose_coarsen::Bool
    verbose_refine::Bool
    verbose_initialize::Bool
    verbose_moves::Bool
    verbose_separators::Bool
    verbose_connectivity::Bool
    verbose_contiguous::Bool
    verbose_memory::Bool
end

MetisPartitioner(;
    assumesymmetry = false,
    checksymmetry = true,
    scheme = nothing,
    objective = nothing,
    coarsen = nothing,
    initialize = nothing,
    refine = nothing,
    ncuts = nothing,
    niterations = nothing,
    seed = nothing,
    minimizeconnectivity = nothing,
    maketwohops = nothing,
    forcecontiguous = nothing,
    compressgraph = nothing,
    imbalance = nothing,
    weight = VertexCount(),
    verbose_info = false,
    verbose_time = false,
    verbose_coarsen = false,
    verbose_refine = false,
    verbose_initialize = false,
    verbose_moves = false,
    verbose_separators = false,
    verbose_connectivity = false,
    verbose_contiguous = false,
    verbose_memory = false,
) = MetisPartitioner(
    assumesymmetry,
    checksymmetry,
    scheme,
    objective,
    coarsen,
    initialize,
    refine,
    ncuts,
    niterations,
    seed,
    minimizeconnectivity,
    maketwohops,
    forcecontiguous,
    compressgraph,
    imbalance,
    weight,
    verbose_info,
    verbose_time,
    verbose_coarsen,
    verbose_refine,
    verbose_initialize,
    verbose_moves,
    verbose_separators,
    verbose_connectivity,
    verbose_contiguous,
    verbose_memory,
)

function metis_partition(A::SparseMatrixCSC, wgt, K, method::MetisPartitioner)
    old_options = copy(Metis.options)

    fill!(Metis.options, Cint(-1))
    Metis.options[Metis.METIS_OPTION_NUMBERING] = 1

    if method.scheme === nothing
        if K > 8
            alg = :KWAY
        else
            alg = :RECURSIVE
        end
    else
        if method.scheme == metis_scheme_direct_k_way
            alg = :KWAY
        elseif method.scheme == metis_scheme_recursive_bisection
            alg = :RECURSIVE
        end
    end

    if method.objective !== nothing
        Metis.options[Metis.METIS_OPTION_OBJTYPE] = Cint(method.objective)
    end
    if method.coarsen !== nothing
        Metis.options[Metis.METIS_OPTION_CTYPE] = Cint(method.coarsen)
    end
    if method.initialize !== nothing
        Metis.options[Metis.METIS_OPTION_IPTYPE] = Cint(method.initialize)
    end
    if method.refine !== nothing
        Metis.options[Metis.METIS_OPTION_RTYPE] = Cint(method.refine)
    end
    if method.ncuts !== nothing
        @assert method.ncuts >= 1
        Metis.options[Metis.METIS_OPTION_NCUTS] = Cint(method.ncuts)
    end
    if method.niterations !== nothing
        @assert method.niterations >= 0
        Metis.options[Metis.METIS_OPTION_NITER] = Cint(method.niterations)
    end
    if method.seed !== nothing
        Metis.options[Metis.METIS_OPTION_SEED] = Cint(method.seed)
    end
    if method.minimizeconnectivity !== nothing
        Metis.options[Metis.METIS_OPTION_MINCONN] = Cint(method.minimizeconnectivity)
    end
    if method.maketwohops !== nothing
        Metis.options[Metis.METIS_OPTION_NO2HOP] = Cint(!method.maketwohops)
    end
    if method.forcecontiguous !== nothing
        Metis.options[Metis.METIS_OPTION_CONTIG] = Cint(method.forcecontiguous)
    end
    if method.compressgraph !== nothing
        Metis.options[Metis.METIS_OPTION_COMPRESS] = Cint(method.compressgraph)
    end
    if method.imbalance !== nothing
        imbalance = ceil(Cint, method.imbalance * 1000)
        @assert imbalance > 0
        Metis.options[Metis.METIS_OPTION_UFACTOR] = imbalance
    end

    verbose = Cint(0)
    if method.verbose_info verbose |= Metis.METIS_DBG_INFO end
    if method.verbose_time verbose |= Metis.METIS_DBG_TIME end
    if method.verbose_coarsen verbose |= Metis.METIS_DBG_COARSEN end
    if method.verbose_refine verbose |= Metis.METIS_DBG_REFINE end
    if method.verbose_initialize verbose |= Metis.METIS_DBG_IPART end
    if method.verbose_moves verbose |= Metis.METIS_DBG_MOVEINFO end
    if method.verbose_separators verbose |= Metis.METIS_DBG_SEPINFO end
    if method.verbose_connectivity verbose |= Metis.METIS_DBG_CONNINFO end
    if method.verbose_contiguous verbose |= Metis.METIS_DBG_CONTIGINFO end
    if method.verbose_memory verbose |= Metis.METIS_DBG_MEMORY end

    g = Metis.graph(A)
    if wgt !== nothing
        g = Metis.Graph(g.nvtxs, g.xadj, g.adjncy, convert(Vector{typeof(g.nvtxs)}, wgt))
    end

    if verbose == 0
        asg = @suppress Metis.partition(A, K, alg=alg)
    else
        asg = Metis.partition(A, K, alg=alg)
    end

    Metis.options .= old_options

    return asg
end

function partition_stripe(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, K, method::MetisPartitioner; kwargs...) where {Tv, Ti}
    wgt = method.weight isa VertexCount ? nothing : compute_weight(A, method.weight)
    asg = metis_partition(pattern(SparseMatrixCSC(A)), wgt, K, method)
    return MapPartition(K, asg)
end

function partition_stripe(A::SparseMatrixCSC, K, method::MetisPartitioner; kwargs...)
    return partition_plaid(A, K, method)[end]
end

function partition_plaid(A::Symmetric{Tv, SparseMatrixCSC{Tv, Ti}}, K, method::MetisPartitioner; kwargs...) where {Tv, Ti}
    Φ = partition_stripe(A, K, method, kwargs...)
    return (Φ, Φ)
end

function partition_plaid(A::SparseMatrixCSC, K, method::MetisPartitioner; kwargs...)
    if method.assumesymmetry || (method.checksymmetry && issymmetric(pattern(A)))
        wgt = method.weight isa VertexCount ? nothing : compute_weight(A, method.weight)
        asg = metis_partition(pattern(A), wgt, K, method)
        return (MapPartition(K, asg), MapPartition(K, asg))
    else
        m, n = size(A)
        AP = pattern(A)
        B = hvcat((2, 2), spzeros(Bool, m, m), AP, adjointpattern(AP), spzeros(Bool, n, n))
        wgt = compute_weight(A, method.weight)
        wgt = vcat(zeros(eltype(wgt), m), wgt)
        asg = metis_partition(B, wgt, K, method)
        return (MapPartition(K, asg[1:m]), MapPartition(K, asg[m + 1:end]))
    end
end


text:
"""
    pyiter(x)

Equivalent to `iter(x)` in Python.
"""
pyiter(x) = pynew(errcheck(@autopy x C.PyObject_GetIter(getptr(x_))))
export pyiter

"""
    pynext(x)

Equivalent to `next(x)` in Python.
"""
pynext(x) = pybuiltins.next(x)
export pynext

"""
    unsafe_pynext(x)

Return the next item in the iterator `x`. When there are no more items, return NULL.
"""
unsafe_pynext(x::Py) = GC.@preserve x pynew(errcheck_ambig(C.PyIter_Next(getptr(x))))


text:
module Optimize

include("orbits.jl")
include("observables.jl")

using Schwarz.Orbits: OrbitLibrary, BinnedLibrary
using Schwarz.Observables: KinematicData, binedges

function objective(lib::BinnedLibrary, data::KinematicData)
    #=
    From an orbit library and kinematic data, return the objective function,
    its gradient, constraints, etc.
    =#
    # lib = BinnedLibrary(lib, data.edges)
    # m = fractional mass times sigma^2
    # m_lib is n_orbit rows by n_bins columns matrix
    m_lib = [lib.orbits[i].mass[j] * lib.orbits[i].dispersion[j]
             for i in 1:length(lib), j in 1:length(data)]
    mass_lib = [lib.orbits[i].mass[j]
                for i in 1:length(lib), j in 1:length(data)]
    m_data = [data.mass[j] * data.dispersion[j] for j in 1:length(data)]
    dm_data = [sqrt((data.unc_mass[j] * data.dispersion[j])^2
                    + (data.unc_dispersion[j] * data.mass[j])^2)
               for j in 1:length(data)]
    u_data = 1 ./ dm_data .^ 2
    Q = Array{Float64, 2}(length(lib), length(lib))
    for i in 1:length(lib), j in 1:i
        if i == j
            Q[i, j] = 2 * sum(u_data .* m_lib[i, :].^2)
        else
            x = 4 * sum(u_data .* m_lib[i, :] .* m_lib[j, :])
            Q[i, j] = x
            Q[j, i] = x
        end
    end
    c = -2 * [sum(u_data .* m_data .* m_lib[i, :]) for i in 1:length(lib)]

    # objective
    function eval_f(weights::Vector{Float64})
        0.5 * dot(weights, Q * weights) + dot(c, weights)
    end

    # gradient of objective
    function eval_grad_f(weights::Vector{Float64}, grad_f::Vector{Float64})
        grad_f[:] = (Q * weights) .+ c
    end

    # constraints
    function eval_g(weights::Vector{Float64}, g::Vector{Float64})
        weighted_mass_lib = [dot(weights, mass_lib[:, j]) for j in 1:length(data)]
        # upper bound of 1
        g[1:(end - 1)] = abs(weighted_mass_lib .- data.mass) ./ data.unc_mass
        # upper and lower bound of 1
        g[end] = sum(weights)
    end

    # jacobian of constraints
    function eval_jac_g(weights::Vector{Float64}, mode,
                        rows::Vector{Int32}, cols::Vector{Int32},
                        values::Vector{Float64})
        #=
        rows: n_bins + 1 constraints, cols: n_orbit variables
        =#
        weighted_mass_lib = [dot(weights, mass_lib[:, j]) for j in 1:length(data)]
        sgn = map(j -> sign(data.mass[j] - weighted_mass_lib[j]),
                  1:length(data))
        if mode == :Structure
            idx = 1
            for row in 1:(length(data) + 1)
                for col in 1:length(lib)
                    rows[idx] = row
                    cols[idx] = col
                    idx += 1
                end
            end
        else
            count = 1
            for i in 1:length(lib)
                values[count:(count - 1 + length(data))] =
                    [sgn[j] * mass_lib[i, j] / data.unc_mass[j]
                     for j in 1:length(data)]
                values[count + length(data)] = 1.0
                count += length(data) + 1
            end
        end
    end

    function eval_h(weights::Vector{Float64}, mode,
                    rows::Vector{Int32}, cols::Vector{Int32},
                    obj_factor::Float64, lambda::Vector{Float64},
                    values::Vector{Float64})
        if mode == :Structure
            idx = 1
            for row in 1:(length(lib))
                for col in 1:row
                    rows[idx] = row
                    cols[idx] += 1
                end
            end
        else
            idx = 1
            for row in 1:(length(lib))
                for col in 1:row
                    values[idx] = obj_factor * Q[row, col]
                    idx += 1
                end
            end
        end
    end

    return eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h
end

function fit(lib::BinnedLibrary, data::KinematicData)
    nbins = length(data)
    n = length(lib) # number of variables
    x_L = zeros(n)
    x_U = ones(n)
    m = nbins + 1 # number of constraints
    g_L = vcat(zeros(nbins), 1)
    g_U = ones(m)
    nele_jac = m * n
    nele_hess = sum(1:n)
    eval_f, eval_grad_f, eval_g, eval_jac_g, eval_h = objective(lib, data)

    println("Creating problem")
    prob = createProblem(n, x_L, x_U, m, g_L, g_U, nele_jac, nele_hess,
                         eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
    prob.x = ones(n) / n
    println("Solving problem")
    status = solveProblem(prob)
    println(Ipopt.ApplicationReturnStatus[status])
    println("Weights: ", prob.x)
    println("Objective: ", prob.obj_val)
    return prob
end

end


text:
export  Update,
        addcolumnupdate!,
        update

"""
```
Update( name::Union{Var, AbstractString} )
```
This mutable struct represents a `UPDATE` SQL statement where `name` is the name of the affected table.
"""
mutable struct Update <: SQLStatement
    name::Var
    columns::Vector{Var}
    values::Vector
    wherec::SQLCondition
    returning::Vector{String}

    Update( name::SVar ) = new( name |> Var, Var[], Any[], nocond, String[] )
end

function Base.show( io::IO, up::Update )
    isempty(up.columns) && return
    print( io, "UPDATE ", up.name, "\n    SET ", join( string.( up.columns, " = ", processentry.(up.values) ), ", " ) )
    isempty(up.wherec) || print( io, "\n    WHERE ", up.wherec )
    isempty(up.returning) || print( io, "\n    RETURNING ", join( up.returning, ", " ) )
end

"""
```
addcolumnupdate!( colname::Union{Var, AbstractString}, val )
```
This function permits the user to add a command that the column `colname` with the value `val` to a `Update` struct.

The output of this function is a function that takes a `Update` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> addcolumnupdate!( "var1", "val1" ) |> addcolumnupdate!( "var2", "val2" )
```
sets the commands to update column `var1` to `val1` and `var2` to `val2` in the `Update` struct.
"""
addcolumnupdate!( colname::SVar, val ) = up::Update -> begin
    push!( up.columns, colname |> Var )
    push!( up.values, val )
    up
end

"""
```
update( conn::Conn, df::Bool=true )
```
This function permits the user to execute the SQL equivalent of a `Update` struct on the database that `conn` connects to. The flag `df` determines the format of the output (`DataFrame` or `NamedTuple`) if the statement will return an output.

The output of this function is a function that takes a `Update` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> update(conn)
```
"""
update( conn::Conn, df::Bool=true ) = up::Update -> begin
    res = LibPQ.execute( conn, string( up, ";" ) )
    isempty(up.returning) && return res
    res |> makeunique! |> tableresult(df)
end


text:
# The unit convention throughout the package and other packages in the CESMIX Julia Suite
# is to assume all unspecified units to be atomic units as defined in UnitfulAtomic.jl.

# Here we provide constants for the atomic units we use to make the code more readable.

MASS_UNIT = u"me_au"
LENGTH_UNIT = u"bohr"
ENERGY_UNIT = u"hartree"
TIME_UNIT = u"ħ_au / hartree"
VELOCITY_UNIT = u"bohr * hartree / ħ_au"
TEMPERATURE_UNIT = u"hartree / k_au"
FORCE_UNIT = u"hartree / bohr"

MASS_TYPE = typeof(zeros() * MASS_UNIT)
LENGTH_TYPE = typeof(zeros() * LENGTH_UNIT)
ENERGY_TYPE = typeof(zeros() * ENERGY_UNIT)
TIME_TYPE = typeof(zeros() * TIME_UNIT)
VELOCITY_TYPE = typeof(zeros() * VELOCITY_UNIT)
TEMPERATURE_TYPE = typeof(zeros() * TEMPERATURE_UNIT)
FORCE_TYPE = typeof(zeros() * FORCE_UNIT)