text:
"""
    StateCheck

Module with a minimal set of functions for getting statistics
and basic I/O from ClimateMachine DG state arrays (`MPIStateArray` type).
Created for regression testing and code change tracking and debugging.
StateCheck functions iterate over named variables in an `MPIStateArray`,
calculate and report their statistics and/or write values for all or
some subset of points at a fixed frequency.

# Functions

 - [`sccreate`](@ref) Create a StateCheck call back variable.
 - [`scdocheck`](@ref) Check StateCheck variable values against reference values.
 - [`scprintref`](@ref) Print StateCheck variable in format for creating reference values.
"""
module StateCheck

# Imports from standard Julia packages
using Formatting
using MPI
using PrettyTables
using Printf
using StaticArrays
using Statistics

# Imports from ClimateMachine core
import ..GenericCallbacks: EveryXSimulationSteps
import ..MPIStateArrays: MPIStateArray
import ..VariableTemplates: flattenednames

####
# For testing put a new function signature here!
# Needs to go in src/Utilities/VariableTemplates/var_names.jl
# This handles SMatrix case
flattenednames(::Type{T}; prefix = "") where {T <: SArray} =
    ntuple(i -> "$prefix[$i]", length(T))
####

"""
    VStat

Type for returning variable statistics.
"""
struct VStat
    max
    min
    mean
    std
end

# Global functions to expose
export sccreate # Create a state checker call back
export scdocheck
export scprintref

const nt_freq_def = 10 # default frequency (in time steps) for output.
const prec_def = 15    # default precision used for formatted output table

# TODO: this should use the new callback interface

"""
    sccreate(
        io::IO,
        fields::Array{<:Tuple{<:MPIStateArray, String}},
        nt_freq::Int = $nt_freq_def;
        prec = $prec_def
    )

Create a "state check" call-back for one or more `MPIStateArray` variables
that will report basic statistics for the fields in the array.

  -  `io` an IO stream to use for printed output
  -  `fields` a required first argument that is an array of one or more
                        `MPIStateArray` variable and label string pair tuples.
                        State array statistics will be reported for the named symbols
                        in each `MPIStateArray` labeled with the label string.
  -  `nt_freq` an optional second argument with default value of
                        $nt_freq_def that sets how frequently (in time-step counts) the
                        statistics are reported.
  -  `prec` a named argument that sets number of decimal places to print for
                        statistics, defaults to $prec_def.

# Examples
```julia
using ClimateMachine.VariableTemplates
using StaticArrays
using ClimateMachine.MPIStateArrays
using MPI
MPI.Init()
FT=Float64
F1=@vars begin; ν∇u::SMatrix{3, 2, FT, 6}; κ∇θ::SVector{3, FT}; end
F2=@vars begin; u::SVector{2, FT}; θ::SVector{1, FT}; end
Q1=MPIStateArray{Float32,F1}(MPI.COMM_WORLD,ClimateMachine.array_type(),4,9,8);
Q2=MPIStateArray{Float64,F2}(MPI.COMM_WORLD,ClimateMachine.array_type(),4,6,8);
cb=ClimateMachine.StateCheck.sccreate([(Q1,"My gradients"),(Q2,"My fields")],1; prec=$prec_def);
```
"""
function sccreate(
    io::IO,
    fields::Array{<:Tuple{<:MPIStateArray, String}},
    nt_freq::Int = nt_freq_def;
    prec = prec_def,
)
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        println(io, "# SC Start: creating state check callback")
    end

    ####
    # Print fields that the call back created by this call will query
    ####
    nr = 0
    for f in fields
        print_head = true
        Q = f[1]
        lab = f[2]
        V = typeof(Q).parameters[2]
        slist = typeof(Q).parameters[2].names
        l_s = length(slist)
        nss = 0
        if l_s == 0
            println(
                io,
                "# SC  MPIStateArray labeled \"$lab\" has no named symbols.",
            )
        else
            ns = 0
            for s in slist
                ns = ns + 1
                if MPI.Comm_rank(MPI.COMM_WORLD) == 0
                    if print_head
                        println(
                            io,
                            "# SC Creating state check callback labeled \"$lab\" for symbols",
                        )
                    end
                    print_head = false
                    println(io, "# SC ", s)
                end
                for n in
                    flattenednames(fieldtype(V, ns), prefix = fieldname(V, ns))
                    nss = nss + 1
                end
            end
        end
        nr = nr + nss
    end

    ###
    # Initialize total calls counter for the call back
    ###
    n_cb_calls = 0

    ###
    # Create holder for most recent stats
    ###
    cur_stats_dict = Dict()
    cur_stats_flat = Array{Any}(undef, nr)

    # Save io pointer
    iosave = io

    ######
    # Create the callback
    ######
    cb = EveryXSimulationSteps(nt_freq) do
        # Track which timestep this is
        n_cb_calls = n_cb_calls + 1
        n_step = (n_cb_calls - 1) * nt_freq + 1
        ns_str = @sprintf("%7.7d", n_step - 1)
        io = iosave
        # Obscure trick to do with running in cells in notebook that close Base.stdout
        # between different stages. This affects parsing Literate/Documenter examples.
        if !isopen(io)
            io = Base.stdout
        end

        ## Print header
        nprec = min(max(1, prec), 20)
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            println(
                io,
                "# SC +++++++++++ClimateMachine StateCheck call-back start+++++++++++++++++",
            )
            println(
                io,
                "# SC  Step  |   Label    |  Field   |                            Stats                       ",
            )
        end
        h_var_fmt = "%" * sprintf1("%d", nprec + 8) * "s"
        min_str = sprintf1(h_var_fmt, " min() ")
        max_str = sprintf1(h_var_fmt, " max() ")
        ave_str = sprintf1(h_var_fmt, " mean() ")
        std_str = sprintf1(h_var_fmt, " std() ")
        fmt_str = [" min() ", " max() ", " mean() ", " std() "]
        fmt_str = sprintf1.(Ref(h_var_fmt), fmt_str)
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            print(io, "# SC -------|------------|----------|")
            println(io, join(fmt_str, "|"), "|")
        end

        ## Iterate over the set of MPIStateArrays for this callback
        nr = 0
        for f in fields
            olabel = f[2]
            ol_str = @sprintf("%12.12s", olabel)
            m_array = f[1]

            # Get descriptor for MPIStateArray

            V = typeof(m_array).parameters[2]

            ## Iterate over fields in each MPIStateArray
            #  (use ivar to index individual arrays within the MPIStateArray)
            ivar = 0
            stats_val_dict = Dict()
            for i in 1:length(V.names)
                for n in flattenednames(fieldtype(V, i), prefix = fieldname(V, i))
                    nr = nr + 1
                    ivar = ivar + 1
                    n_str = @sprintf("%9.9s", n)
                    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
                        print(io, "# SC ", ns_str, "|", ol_str, "|", n_str, " |")
                    end
                    stats_string = scstats(m_array, ivar, nprec)
                    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
                        println(io, join(stats_string[1:4], "|"), "|")
                    end
                    stats_val_dict[n] = stats_string[5]
                    cur_stats_flat[nr] = [
                        olabel,
                        n,
                        stats_string[5].max,
                        stats_string[5].min,
                        stats_string[5].mean,
                        stats_string[5].std,
                    ]
                end
            end
            cur_stats_dict[olabel] = stats_val_dict
        end
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            println(
                io,
                "# SC +++++++++++ClimateMachine StateCheck call-back end+++++++++++++++++++",
            )
        end
    end

    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        println(io, "# SC Finish: creating state check callback")
    end

    return cb
end

function scstats(V, ivar, nprec)

    # Get number of MPI procs
    nproc = MPI.Comm_size(V.mpicomm)

    npr = nprec
    fmt = @sprintf("%%%d.%de", npr + 8, npr)

    getByField(V, ivar) = (V.realdata[:, ivar, :])
    Vmcomm = V.mpicomm

    # Get local and global field sizes (degrees of freedom).
    n_size_loc = length(getByField(V, ivar))
    n_size_tot = MPI.Allreduce(n_size_loc, +, Vmcomm)

    # Min
    phi_loc = minimum(getByField(V, ivar))
    phi_min = MPI.Allreduce(phi_loc, min, Vmcomm)
    phi = phi_min
    # minVstr=@sprintf("%23.15e",phi)
    min_v_str = sprintf1(fmt, phi)

    # Max
    phi_loc = maximum(getByField(V, ivar))
    phi_max = MPI.Allreduce(phi_loc, max, Vmcomm)
    phi = phi_max
    # maxVstr=@sprintf("%23.15e",phi)
    max_v_str = sprintf1(fmt, phi)

    # Ave
    phi_loc = mean(getByField(V, ivar))
    phi_loc = phi_loc * n_size_loc / n_size_tot
    phi_sum = MPI.Allreduce(phi_loc, +, Vmcomm)
    phi_mean = phi_sum
    phi = phi_mean
    # aveVstr=@sprintf("%23.15e",phi)
    ave_v_str = sprintf1(fmt, phi)

    # Std
    phi_loc = (getByField(V, ivar) .- phi_mean) .^ 2
    phi_loc = sum(phi_loc)  # Sum local data explicitly since GPU Allreduce
    # does not take arrays yet.
    phi_sum = MPI.Allreduce(phi_loc, +, Vmcomm)
    n_val_sum = n_size_tot
    phi_std = sqrt(sum(phi_sum) / (n_val_sum - 1))
    phi = phi_std
    # stdVstr=@sprintf("%23.15e",phi)
    std_v_str = sprintf1(fmt, phi)

    vals = VStat(phi_min, phi_max, phi_mean, phi_std)

    return min_v_str, max_v_str, ave_v_str, std_v_str, vals
end

function sccreate(
    fields::Array{<:Tuple{<:MPIStateArray, String}},
    nt_freq::Int = nt_freq_def;
    prec = prec_def,
)
    return sccreate(Base.stdout, fields, nt_freq; prec = prec)
end

"""
    scprintref(cb)

Print out a "state check" call-back table of values in a format
suitable for use as a set of reference numbers for CI comparison.

 - `cb` callback variable of type ClimateMachine.GenericCallbacks.Every*
"""
function scprintref(cb)
    sc = cb.callback
    io = sc.iosave
    # Obscure trick to do with running in cells in notebook
    if !isopen(io)
        io = Base.stdout
    end
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        # Get print format lengths for cols 1 and 2 so they are aligned
        # for readability.
        phi = sc.cur_stats_flat
        f = 1
        a1l = maximum(length.(map(
            i -> (phi[i])[f],
            range(1, length = length(phi)),
        )))
        f = 2
        a2l = maximum(length.(String.((map(
            i -> (phi[i])[f],
            range(1, length = length(phi)),
        )))))
        fmt1 = @sprintf("%%%d.%ds", a1l, a1l) # Column 1
        fmt2 = @sprintf("%%%d.%ds", a2l, a2l) # Column 2
        fmt3 = @sprintf("%%28.20e")         # All numbers at full precision
        # Create an string of spaces to be used for formatting
        sp = "                                                                           "

        # Write header
        println(io, "# BEGIN SCPRINT")
        println(io, "# varr - reference values (from reference run)    ")
        println(io, "# parr - digits match precision (hand edit as needed) ")
        println(io, "#")
        println(io, "# [")
        println(
            io,
            "#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],",
        )
        println(
            io,
            "#  [         :                :          :        :      :          :           ],",
        )
        println(io, "# ]")
        #
        # Write tables
        #  Reference value and precision match tables are separate since it is more
        #  common to update reference values occasionally while precision values are
        #  typically changed rarely and the precision values are hand edited from experience.
        #
        # Write table of reference values
        println(io, "varr = [")
        for lv in sc.cur_stats_flat
            s1 = lv[1]
            l1 = length(s1)
            s1 = sp[1:(a1l - l1)] * "\"" * s1 * "\""
            s2 = lv[2]
            if typeof(s2) == String
                l2 = length(s2)
                s2 = sp[1:(a2l - l2)] * "\"" * s2 * "\""
                s22 = ""
            end
            if typeof(s2) == Symbol
                s22 = s2
                l2 = length(String(s2))
                s2 = sp[1:(a2l - l2 + 1)] * ":"
            end
            s3 = sprintf1(fmt3, lv[3])
            s4 = sprintf1(fmt3, lv[4])
            s5 = sprintf1(fmt3, lv[5])
            s6 = sprintf1(fmt3, lv[6])
            println(
                io,
                " [ ",
                s1,
                ", ",
                s2,
                s22,
                ",",
                s3,
                ",",
                s4,
                ",",
                s5,
                ",",
                s6,
                " ],",
            )
        end
        println(io, "]")

        # Write table of reference match precisions using default precision that
        # can be hand updated.
        println(io, "parr = [")
        for lv in sc.cur_stats_flat
            s1 = lv[1]
            l1 = length(s1)
            s1 = sp[1:(a1l - l1)] * "\"" * s1 * "\""
            s2 = lv[2]
            if typeof(s2) == String
                l2 = length(s2)
                s2 = sp[1:(a2l - l2)] * "\"" * s2 * "\""
                s22 = ""
            end
            if typeof(s2) == Symbol
                s22 = s2
                l2 = length(String(s2))
                s2 = sp[1:(a2l - l2 + 1)] * ":"
            end
            s3 = sprintf1(fmt3, lv[3])
            s4 = sprintf1(fmt3, lv[4])
            s5 = sprintf1(fmt3, lv[5])
            s6 = sprintf1(fmt3, lv[6])
            println(
                io,
                " [ ",
                s1,
                ", ",
                s2,
                s22,
                ",",
                "    16,    16,    16,    16 ],",
            )
        end
        println(io, "]")
        println(io, "# END SCPRINT")
    end
end

"""
    scdocheck(cb, ref_dat)

Compare a current State check call-back set of values with a
reference set and match precision table pair.

 - `cb` `StateCheck` call-back variables
 - `ref_dat` an array of reference values and precision to match tables.
"""
function scdocheck(cb, ref_dat)
    sc = cb.callback
    io = sc.iosave
    # Obscure trick to do with running in cells in notebook
    if !isopen(io)
        io = Base.stdout
    end
    MPI.Comm_rank(MPI.COMM_WORLD) == 0 || return true
    println(
        io,
        "# SC +++++++++++ClimateMachine StateCheck ref val check start+++++++++++++++++",
    )
    println(
        io,
        "# SC \"N( )\" bracketing indicates field failed to match      ",
    )
    println(io, "# SC \"P=\"  row pass count      ")
    println(io, "# SC \"F=\"  row failure count   ")
    println(io, "# SC \"NA=\" row not checked count      ")
    println(io, "# SC ")
    println(
        io,
        "# SC        Label         Field      min()      max()     mean()      std() ",
    )
    irow = 1
    i_val = 1
    i_prec = 2
    all_pass = true

    for row in sc.cur_stats_flat
        ## Debugging
        # println(row)
        # println(ref_dat[i_val][irow])
        # println(ref_dat[i_prec][irow])
        row_pass = true
        row_col_pass = 0
        row_col_na = 0
        row_col_fail = 0

        ## Make array copy for reporting
        res_dat = copy(ref_dat[i_prec][irow])

        ## Check MPIStateArrayName
        cval = row[1]
        rval = ref_dat[i_val][irow][1]
        if cval != rval
            all_pass = false
            row_pass = false
            row_col_fail += 1
            res_dat[1] = "N" * "(" * rval * ")"
        else
            res_dat[1] = cval
            row_col_pass += 1
            res_dat[1] = rval
        end

        ## Check term name
        cval = row[2]
        rval = ref_dat[i_val][irow][2]
        if cval != rval
            all_pass = false
            row_pass = false
            if typeof(rval) == String
                res_dat[2] = "N" * "(" * rval * ")"
            else
                res_dat[2] = "N" * "(" * string(rval) * ")"
            end
            row_col_fail += 1
        else
            res_dat[2] = cval
            row_col_pass += 1
        end

        # Check numeric values
        nv = 3
        for nv in [3, 4, 5, 6]
            fmt = @sprintf("%%28.20e")
            lfld = 28
            ndig = 20
            cval = row[nv]
            cvalc = sprintf1(fmt, cval)
            rval = ref_dat[i_val][irow][nv]
            rvalc = sprintf1(fmt, rval)
            pcmp = ref_dat[i_prec][irow][nv]

            # Skip if compare digits set to 0
            if pcmp > 0

                # Check exponent part
                ep1 = cvalc[(lfld - 3):lfld]
                ep2 = rvalc[(lfld - 3):lfld]
                if ep1 != ep2
                    nmatch = 0
                else
                    # Now check individual digits left to right
                    dp1 = cvalc[2:(3 + pcmp + 1)]
                    dp2 = rvalc[2:(3 + pcmp + 1)]
                    nmatch = 0
                    imatch = 1
                    for c in dp1
                        if c == dp2[imatch]
                            nmatch = imatch
                        else
                            break
                        end
                        imatch = imatch + 1
                    end
                end
                # Check for trailing exact 0s number change numerically
                if nmatch < pcmp
                    if rval != 0
                        e_diffl10 = round(log10(abs((rval - cval) / rval)))
                    else
                        e_diffl10 = round(log10(abs(rval - cval)))
                    end
                    if e_diffl10 < -pcmp
                        nmatch = Int(-e_diffl10)
                    end
                end
                if nmatch < pcmp
                    all_pass = false
                    row_pass = false
                    res_dat[nv] = "N(" * string(nmatch) * ")"
                    row_col_fail += 1
                else
                    res_dat[nv] = string(nmatch)
                    row_col_pass += 1
                end
            else
                res_dat[nv] = "0"
                row_col_na += 1
            end
        end

        #
        # println(resDat)
        @printf(
            io,
            "# SC %12.12s, %12.12s, %9.9s, %9.9s, %9.9s, %9.9s",
            res_dat[1],
            res_dat[2],
            res_dat[3],
            res_dat[4],
            res_dat[5],
            res_dat[6]
        )
        @printf(
            io,
            " :: P=%d, F=%d, NA=%d\n",
            row_col_pass,
            row_col_fail,
            row_col_na
        )
        # Next row
        irow = irow + 1
    end
    println(
        io,
        "# SC +++++++++++ClimateMachine StateCheck ref val check end+++++++++++++++++",
    )
    return all_pass

end

end # module

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using Flux: onehotbatch, onecold
using Flux.Losses: logitcrossentropy
using Statistics, Random
using Logging: with_logger
using TensorBoardLogger: TBLogger, tb_overwrite, set_step!, set_step_increment!
using ProgressMeter: @showprogress
import BSON
using CUDA
using CSV, DataFrames

function getimages(filename)
    filepath = pwd() * "/images/" * filename

    mtrx = Matrix(DataFrame(CSV.File(filepath)))
    return mtrx'
end

function getlabels(filename)
    filepath = pwd() * "/images/" * filename
    vec(Matrix(DataFrame(CSV.File(filepath))))
end

function load_data(args)
    train_data_file = "csvTrainImages.csv"
    test_data_file = "csvTestImages.csv"
    train_label_file = "csvTrainLabel.csv"
    test_label_file = "csvTestLabel.csv"

    train_data = getimages(train_data_file)
    test_data = getimages(test_data_file)
    train_labels = getlabels(train_label_file)
    test_labels = getlabels(test_label_file)

    xtrain = reshape(train_data, 32, 32, 1, :)
    xtest = reshape(test_data, 32, 32, 1, :)

    ytrain, ytest = onehotbatch(train_labels, 1:28), onehotbatch(test_labels, 1:28)

    train_loader = DataLoader((xtrain, ytrain), batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader((xtest, ytest), batchsize=args.batchsize)

    return train_loader, test_loader
end

function LeNet5(; imgsize=(32,32,1), nclasses=28)
    out_conv_size = (imgsize[1]÷4 - 3, imgsize[2]÷4 - 3, 16)

    return Chain(
            Conv((5, 5), imgsize[end]=>6, relu),
            MaxPool((2, 2)),
            Conv((5, 5), 6=>16, relu),
            MaxPool((2, 2)),
            Flux.flatten,
            Dense(prod(out_conv_size), 120, relu),
            Dense(120, 84, relu),
            Dense(84, nclasses)
          )
end

loss(ŷ, y) = logitcrossentropy(ŷ, y)

function eval_loss_accuracy(loader, model, device)
    l = 0f0
    acc = 0
    ntot = 0
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        l += loss(ŷ, y) * size(x)[end]
        acc += sum(onecold(ŷ |> cpu) .== onecold(y |> cpu))
        ntot += size(x)[end]
    end
    return (loss = l/ntot |> round4, acc = acc/ntot*100 |> round4)
end

## utility functions
num_params(model) = sum(length, Flux.params(model))
round4(x) = round(x, digits=4)

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 3e-4             # learning rate
    λ = 0                # L2 regularizer param, implemented as weight decay
    batchsize = 128      # batch size
    epochs = 10          # number of epochs
    seed = 0             # set seed > 0 for reproducibility
    use_cuda = true      # if true use cuda (if available)
    infotime = 1 	     # report every `infotime` epochs
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    tblogger = true      # log training with tensorboard
    savepath = "runs/"    # results path
end

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.use_cuda && CUDA.functional()

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader = load_data(args)
    @info "Dataset MNIST: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = LeNet5() |> device
    @info "LeNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0 # add weight decay, equivalent to L2 regularization
        opt = Optimiser(WeightDecay(args.λ), opt)
    end

    ## LOGGING UTILITIES
    if args.tblogger
        tblogger = TBLogger(args.savepath, tb_overwrite)
        set_step_increment!(tblogger, 0) # 0 auto increment since we manually set_step!
        @info "TensorBoard logging at \"$(args.savepath)\""
    end

    function report(epoch)
        train = eval_loss_accuracy(train_loader, model, device)
        test = eval_loss_accuracy(test_loader, model, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
        if args.tblogger
            set_step!(tblogger, epoch)
            with_logger(tblogger) do
                @info "train" loss=train.loss  acc=train.acc
                @info "test"  loss=test.loss   acc=test.acc
            end
        end
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        @showprogress for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                    ŷ = model(x)
                    loss(ŷ, y)
                end

            Flux.Optimise.update!(opt, ps, gs)
        end

        ## Printing and logging
        epoch % args.infotime == 0 && report(epoch)
        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model.bson")
            let model = cpu(model) #return model to cpu before serialization
                BSON.@save modelpath model epoch
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
end

train()

ext: jl

max_stars_repo_name: whyharsha/arabicdigits

max_issues_repo_name: whyharsha/arabicdigits

max_forks_repo_name: whyharsha/arabicdigits


text:
module FakeModule2
    export f, g

    f(x) = x
    g(x) = x + 1
end


text:
using PyCall
@pyimport matplotlib.pyplot as plt
include("MogTrack.jl")

# Make a MogTrack
mogtrack = MogTrack(50)
from_prior!(mogtrack)

# Define a grid
xmin, xmax, ymin, ymax = -5.0, 5.0, -5.0, 5.0
ni, nj = 201, 201
x = Array(linspace(xmin, xmax, nj))
y = Array(linspace(ymax, ymin, ni))

# Evaluate the MogTrack on the grid
f = Array(Float64, (ni, nj))

plt.ion()
steps = 1000
for(k in 1:steps)
	perturb!(mogtrack)

	for(j in 1:nj)
		for(i in 1:ni)
			f[i, j] = evaluate(mogtrack, x[j], y[i])
		end
	end

	plt.hold(false)
	plt.imshow(-f, interpolation="nearest", cmap="gray")
	plt.title("Model $k/$steps")
	plt.gca()[:set_xticks]([])
	plt.gca()[:set_yticks]([])
	plt.draw()
end

plt.ioff()
plt.show()

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
module InverseProblems

using OrderedCollections
using Suppressor: @suppress

using ..Observations: obs_str, AbstractObservation, SyntheticObservations, initialize_simulation!, FieldTimeSeriesCollector,
    observation_times, observation_names

using ..TurbulenceClosureParameters: free_parameters_str, new_closure_ensemble

using OffsetArrays, Statistics

using Oceananigans: short_show, run!, fields, FieldTimeSeries, CPU
using Oceananigans.OutputReaders: InMemory
using Oceananigans.Fields: interior, location
using Oceananigans.Grids: Flat, Bounded,
                          Face, Center,
                          RectilinearGrid, offset_data,
                          topology, halo_size,
                          interior_parent_indices

using Oceananigans.Models.HydrostaticFreeSurfaceModels: SingleColumnGrid, YZSliceGrid, ColumnEnsembleSize

import ..Observations: normalize!

#####
##### Output maps (maps from simulation output to observation space)
#####

abstract type AbstractOutputMap end

output_map_type(fp) = output_map_str(fp)

struct ConcatenatedOutputMap{T}
    time_indices::T
end

ConcatenatedOutputMap(; time_indices = Colon()) = ConcatenatedOutputMap(time_indices)

output_map_str(::ConcatenatedOutputMap) = "ConcatenatedOutputMap"

"""
    ConcatenatedVectorNormMap()

Forward map transformation of simulation output to a scalar by
taking a naive `norm` of the difference between concatenated vectors of the
observations and simulation output.
"""
struct ConcatenatedVectorNormMap end

output_map_str(::ConcatenatedVectorNormMap) = "ConcatenatedVectorNormMap"

#####
##### InverseProblems
#####

struct InverseProblem{F, O, S, T, P}
    observations :: O
    simulation :: S
    time_series_collector :: T
    free_parameters :: P
    output_map :: F
end

"""
    InverseProblem(observations, simulation, free_parameters; output_map=ConcatenatedOutputMap())

Return an `InverseProblem`.
"""
function InverseProblem(observations,
                        simulation,
                        free_parameters;
                        output_map = ConcatenatedOutputMap(),
                        time_series_collector = nothing)

    if isnothing(time_series_collector) # attempt to construct automagically
        simulation_fields = fields(simulation.model)
        collected_fields = NamedTuple(name => simulation_fields[name] for name in observation_names(observations))
        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))
    end

    return InverseProblem(observations, simulation, time_series_collector, free_parameters, output_map)
end

function Base.show(io::IO, ip::InverseProblem)
    sim_str = "Simulation on $(short_show(ip.simulation.model.grid)) with Δt=$(ip.simulation.Δt)"

    out_map_type = output_map_type(ip.output_map)
    out_map_str = output_map_str(ip.output_map)

    print(io, "InverseProblem{$out_map_type}", '\n',
        "├── observations: $(obs_str(ip.observations))", '\n',
        "├── simulation: $sim_str", '\n',
        "├── free_parameters: $(free_parameters_str(ip.free_parameters))", '\n',
        "└── output map: $out_map_str")

    return nothing
end

tupify_parameters(ip, θ) = NamedTuple{ip.free_parameters.names}(Tuple(θ))
tupify_parameters(ip, θ::NamedTuple) = NamedTuple(name => θ[name] for name in ip.free_parameters.names)

"""
    expand_parameters(ip, θ)

Convert `θ` to `Vector{<:NamedTuple}`, where the elements
correspond to `ip.free_parameters`.

`θ` may represent an ensemble of parameter sets via:

* `θ::Vector{<:Vector}` (caution: parameters must be ordered correctly!)
* `θ::Matrix` (caution: parameters must be ordered correctly!)
* `θ::Vector{<:NamedTuple}`

or a single parameter set if `θ::Vector{<:Number}`.

If `length(θ)` is less the the number of ensemble members in `ip.simulation`, the
last parameter set is copied to fill the parameter set ensemble.
"""
function expand_parameters(ip, θ::Vector)
    Nfewer = n_ensemble(ip) - length(θ)
    Nfewer < 0 && throw(ArgumentError("There are $(-Nfewer) more parameter sets than ensemble members!"))

    θ = [tupify_parameters(ip, θi) for θi in θ]

    # Fill out parameter set ensemble
    Nfewer > 0 && append!(θ, [θ[end] for _ = 1:Nfewer])

    return θ
end

# Expand single parameter set
expand_parameters(ip, θ::Union{NamedTuple,Vector{<:Number}}) = expand_parameters(ip, [θ])
expand_parameters(ip, θ::NamedTuple) = expand_parameters(ip, [θ])

# Convert matrix to vector of vectors
expand_parameters(ip, θ::Matrix) = expand_parameters(ip, [θ[:, i] for i = 1:size(θ, 2)])

#####
##### Forward map evaluation given vector-of-vector (one parameter vector for each ensemble member)
#####

const OneDimensionalEnsembleGrid = RectilinearGrid{<:Any, Flat, Flat, Bounded}
const TwoDimensionalEnsembleGrid = RectilinearGrid{<:Any, Flat, Bounded, Bounded}

n_ensemble(grid::Union{OneDimensionalEnsembleGrid, TwoDimensionalEnsembleGrid}) = grid.Nx
n_observations(grid::OneDimensionalEnsembleGrid) = grid.Ny
n_observations(grid::TwoDimensionalEnsembleGrid) = 1
n_z(grid::Union{OneDimensionalEnsembleGrid, TwoDimensionalEnsembleGrid}) = grid.Nz
n_y(grid::TwoDimensionalEnsembleGrid) = grid.Ny
n_ensemble(ip::InverseProblem) = n_ensemble(ip.simulation.model.grid)

""" Transform and return `ip.observations` appropriate for `ip.output_map`. """
observation_map(ip::InverseProblem) = transform_time_series(ip.output_map, ip.observations)

"""
    forward_run!(ip, parameters)

Initialize `ip.simulation` with `parameters` and run it forward. Output is stored
in `ip.time_series_collector`.
"""
function forward_run!(ip::InverseProblem, parameters)
    observations = ip.observations
    simulation = ip.simulation
    closures = simulation.model.closure

    θ = expand_parameters(ip, parameters)

    simulation.model.closure = new_closure_ensemble(closures, θ)

    initialize_simulation!(simulation, observations, ip.time_series_collector)

    @suppress run!(simulation)

    return nothing
end

"""
    forward_map(ip, parameters)

Run `ip.simulation` forward with `parameters` and return the data,
transformed into an array format expected by `EnsembleKalmanProcesses.jl`.
"""
function forward_map(ip, parameters)

    # Run the simulation forward and populate the time series collector
    # with model data.
    forward_run!(ip, parameters)

    # Transform the model data according to `ip.output_map` into
    # the array format expected by EnsembleKalmanProcesses.jl
    # The result has `size(output) = (output_size, ensemble_capacity)`,
    # where `output_size` is determined by both the `output_map` and the
    # data collection dictated by `ip.observations`.
    output = transform_output(ip.output_map, ip.observations, ip.time_series_collector)

    # (output_size, ensemble_size)
    return output
end

(ip::InverseProblem)(θ) = forward_map(ip, θ)

"""
    transform_time_series(::ConcatenatedOutputMap, time_series::SyntheticObservations)

Concatenates flattened, normalized data for each field in the `time_series`.
"""
function transform_time_series(output_map::ConcatenatedOutputMap, time_series::SyntheticObservations)
    flattened_normalized_data = []

    for field_name in keys(time_series.field_time_serieses)
        field_time_series = time_series.field_time_serieses[field_name]
        field_time_series_data = Array(interior(field_time_series))[:, :, :, output_map.time_indices]

        Nx, Ny, Nz, Nt = size(field_time_series_data)

        field_time_series_data = reshape(field_time_series_data, Nx, Ny * Nz * Nt)

        normalize!(field_time_series_data, time_series.normalization[field_name])

        push!(flattened_normalized_data, field_time_series_data)
    end

    transformed = hcat(flattened_normalized_data...)

    return Matrix(transpose(transformed))
end

"""
    transform_time_series(map, time_serieses::Vector)

Return the `transform_time_series` of each `time_series` in `time_serieses` vector.
"""
transform_time_series(map, time_serieses::Vector) =
    vcat(Tuple(transform_time_series(map, time_series) for time_series in time_serieses)...)

"""
    observation_map_variance_across_time(map::ConcatenatedOutputMap, observation::SyntheticObservations)

Returns an (Nx, Ny*Nz*Nfields, Ny*Nz*Nfields) array storing the covariance of each element of the observation
map measured across time, for each ensemble member, where `Nx` is the ensemble size, `Ny` is the batch size,
`Nz` is the number of grid elements in the vertical, and `Nfields` is the number of fields in `observation`.
"""
function observation_map_variance_across_time(map::ConcatenatedOutputMap, observation::SyntheticObservations)

    N_fields = length(keys(observation.field_time_serieses))

    a = transform_time_series(map, observation)
    a = transpose(a) # (Nx, Ny*Nz*Nt)

    example_field_time_series = values(observation.field_time_serieses)[1]

    Nx, Ny, Nz, Nt = size(interior(example_field_time_series))

    # Assume all fields have the same size
    b = reshape(a, Nx, Ny * Nz, Nt, N_fields) # (Nx, Ny*Nz, Nt, Nfields)

    c = cat((b[:, :, :, i] for i = 1:N_fields)..., dims = 2) # (Nx, Ny*Nz*Nfields, Nt)

    ds = [reshape(var(c[:, :, 1:t], dims = 3), Nx, Ny * Nz, N_fields) for t = 1:Nt]

    e = cat(ds..., dims = 2)

    replace!(e, NaN => 0) # variance for first time step is zero

    return reshape(e, Nx, Ny * Nz * Nt * N_fields)
end

observation_map_variance_across_time(map::ConcatenatedOutputMap, observations::Vector) =
    hcat(Tuple(observation_map_variance_across_time(map, observation) for observation in observations)...)

observation_map_variance_across_time(ip::InverseProblem) = observation_map_variance_across_time(ip.output_map, ip.observations)

function transform_output(map::ConcatenatedOutputMap,
    observations::Union{SyntheticObservations,Vector{<:SyntheticObservations}},
    time_series_collector)

    # transposed_output isa Vector{SyntheticObservations} where SyntheticObservations is Nx by Nz by Nt
    transposed_output = transpose_model_output(time_series_collector, observations)

    return transform_time_series(map, transposed_output)
end

vectorize(observation) = [observation]
vectorize(observations::Vector) = observations

const YZSliceObservations = SyntheticObservations{<:Any,<:YZSliceGrid}

transpose_model_output(time_series_collector, observations::YZSliceObservations) =
    SyntheticObservations(time_series_collector.field_time_serieses,
                          time_series_collector.grid,
                          time_series_collector.times,
                          nothing,
                          nothing,
                          observations.normalization)

"""
    transpose_model_output(time_series_collector, observations)

Transpose a `NamedTuple` of 4D `FieldTimeSeries` model output collected by `time_series_collector`
into a Vector of `SyntheticObservations` for each member of the observation batch.

Return a 1-vector in the case of singleton observations.
"""
function transpose_model_output(time_series_collector, observations)
    observations = vectorize(observations)
    times = time_series_collector.times

    transposed_output = []

    n_ensemble = time_series_collector.grid.Nx
    n_batch = time_series_collector.grid.Ny
    Nz = time_series_collector.grid.Nz
    Hz = time_series_collector.grid.Hz
    Nt = length(times)

    grid = drop_y_dimension(time_series_collector.grid)

    for j = 1:n_batch
        observation = observations[j]
        time_serieses = OrderedDict{Any,Any}()

        for name in keys(observation.field_time_serieses)
            loc = LX, LY, LZ = location(observation.field_time_serieses[name])
            topo = topology(grid)

            field_time_series = time_series_collector.field_time_serieses[name]

            raw_data = parent(field_time_series.data)
            data = OffsetArray(view(raw_data, :, j:j, :, :), 0, 0, -Hz, 0)

            time_series = FieldTimeSeries{LX, LY, LZ, InMemory}(data, grid, nothing, times)
            time_serieses[name] = time_series
        end

        # Convert to NamedTuple
        time_serieses = NamedTuple(name => time_series for (name, time_series) in time_serieses)

        batch_output = SyntheticObservations(time_serieses,
                                             grid,
                                             times,
                                             nothing,
                                             nothing,
                                             observation.normalization)

        push!(transposed_output, batch_output)
    end

    return transposed_output
end

function drop_y_dimension(grid::RectilinearGrid{<:Any, <:Flat, <:Flat, <:Bounded})
    new_size = ColumnEnsembleSize(Nz=grid.Nz, ensemble=(grid.Nx, 1), Hz=grid.Hz)
    z_domain = (grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz])
    new_grid = RectilinearGrid(size=new_size, z=z_domain, topology=(Flat, Flat, Bounded))
    return new_grid
end

end # module


text:
function init_test_capital_grid()
    kGrid = CapitalGrid(-2.0, 50.0, 300);
    @assert validate_capital_grid(kGrid);
    return kGrid
end

function validate_capital_grid(kGrid :: CapitalGrid)
    isValid = (kGrid.kMax > kGrid.kMin);
    return isValid
end

function validate_capital_grid(kGridV :: AbstractVector{T}) where T
    isValid = all(diff(kGridV) .> 0.0);
    return isValid
end

k_min(kGrid :: CapitalGrid, t) = kGrid.kMin;
k_max(kGrid :: CapitalGrid, t) = kGrid.kMax;
kprime_min(kGrid :: CapitalGrid, t) = k_min(kGrid, t+1);
kprime_max(kGrid :: CapitalGrid, t) = k_max(kGrid, t+1);
make_k_grid(kGrid :: CapitalGrid, t) = LinRange(k_min(kGrid, t), k_max(kGrid, t), kGrid.nk);

# ------------


text:
"""
# Available functions
- `initialize`
- `terminate`
- `get_tc_temperature` (tc = thermocouple)
- `get_voltage`
- `get_current`
- `get_resistance(;wire)` # wire must be set to 2 or 4
- `get_channel` # (some kind of input detection not selection)
"""
struct KeysightDMM34465A <: MultiMeter end

"""
    get_tc_temperature(multimeter)

Perform take a measurement with the probe mode set to thermocouple
"""
function get_tc_temperature(obj::Instr{KeysightDMM34465A})
    units = get_temp_unit(obj)
    return f_query(obj, "MEASURE:TEMPERATURE? TC"; timeout=0) * units
end

"""
    set_tc_type(multimeter; type="K")

# Keywords
- `type`: Can be E, J, K, N, R, T (Defaults to K)
"""
function set_tc_type(obj::Instr{KeysightDMM34465A}; type="K")
    if !(string(type) in ["E", "J", "K", "N", "R", "T"])
        error("$type must be one of [E, J, K, N, R, T]")
    end
    write(obj, "CONFIGURE:TEMPERATURE TC,$type")
end

"""
Returns voltage

# Keywords
- `type`: "DC" | "AC" (Default DC)

"""
function get_voltage(obj::Instr{KeysightDMM34465A}; type="DC")
    !(type in ["AC","DC"]) && error("$type not valid!\nMust be AC or DC")
    f_query(obj, "MEASURE:VOLTAGE:$type?"; timeout=0) * V
end

"""
    get_current(obj::Instr{KeysightDMM34465A}; type="DC")

Returns current

# Keywords
- `type`: "DC" | "AC" (Default DC)

"""
function get_current(obj::Instr{KeysightDMM34465A}; type="DC")
    !(type in ["AC","DC"]) && error("$type not valid!\nMust be AC or DC")
    f_query(obj, "MEASURE:CURRENT:$type?"; timeout=0) * A
end

"""
    get_resistance(multimeter; wire=2)
    get_resistance(multimeter; wire=4)
Returns resistance

# Keywords
- `wire`: 2 | 4 (Required)

"""
function get_resistance(obj::Instr{KeysightDMM34465A}; wire)
    if wire == 2
        f_query(obj, "MEASURE:RESISTANCE?"; timeout=0) * R
    elseif wire == 4
        f_query(obj, "MEASURE:FRESISTANCE?"; timeout=0) * R
    else
        error("wire flag must be 2 or 4 not $wire")
    end
end

"""
    set_temp_unit_celsius(multimeter)

Sets the temperature unit to celcius
"""
set_temp_unit_celsius(obj::Instr{KeysightDMM34465A}) =
    write(obj, "UNIT:TEMPERATURE C")

"""
    set_temp_unit_farenheit(multimeter)

Sets the temperature unit to farenheit
"""
set_temp_unit_farenheit(obj::Instr{KeysightDMM34465A}) =
    write(obj, "UNIT:TEMPERATURE F")

"""
    set_temp_unit_kelvin(multimeter)

Sets the temperature unit to kelvin
"""
set_temp_unit_kelvin(obj::Instr{KeysightDMM34465A}) =
    write(obj, "UNIT:TEMPERATURE K")

"""
    get_temp_unit(multimeter)

Returns C, F or K depending on the set temperature unit
"""
function get_temp_unit(obj::Instr{KeysightDMM34465A})
   units = query(obj, "UNIT:TEMPERATURE?")
   return if units == "C"
       u"C"
   elseif units == "F"
       u"F"
   elseif units == "K"
       u"K"
   else
       error("Expected [C, F, K]. Got: $units")
   end
end

"""
    get_channel(obj::Instr{KeysightDMM34465A})
Indicates which input terminals are selected on the front panel
Front/Rear switch. This switch is not programmable; this query reports
the position of the switch, but cannot change it.

Do not toggle the Front/Rear switch with active signals on the
terminals. This switch is not intended to be used in this way and may be damaged by high voltages or currents, possibly compromising the
instrument's safety features.

# Returns
- "FRON" or "REAR"
"""
get_channel(obj::Instr{KeysightDMM34465A}) =
    query(obj, "ROUTE:TERMINALS?")


text:
using Optim
using GR
using Plots

gr()

f(x) =  x[1]^2+3*x[1]+3+x[1]^4/5

start = [4.0]
iter = 10

GD_res = optimize(f, start, GradientDescent(),
               Optim.Options(g_tol = 1e-20,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))

CGD_res = optimize(f, start, ConjugateGradient(),
               Optim.Options(g_tol = 1e-20,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))

SA_res = optimize(f, start, SimulatedAnnealing(),
               Optim.Options(g_tol = 1e-20,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))
GD_trace = collect(Iterators.flatten(Optim.x_trace(GD_res)))
CGD_trace = collect(Iterators.flatten(Optim.x_trace(CGD_res)))
SA_trace = collect(Iterators.flatten(Optim.x_trace(SA_res)))

plot(x->f([x,]), lab = "Univariate Function", linecolor = :grey)

plot!(GD_trace,
    [Optim.f_trace(GD_res)],
    lab = "Gradient Descent",
    zcolor = reverse(GD_trace), m = (4, 0.9, :reds, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :red)

plot!(CGD_trace,
    [Optim.f_trace(CGD_res)],
    lab = "Conjugated GD",
    zcolor = reverse(CGD_trace), m = (4, 0.9, :greens, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :green)

plot!(SA_trace,
    [Optim.f_trace(SA_res)],
    lab = "Simulated Annealing",
    zcolor = reverse(SA_trace), m = (4, 0.9, :blues, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :blue)

max_stars_repo_name: lukaskln/Optim-Visualization

max_issues_repo_name: lukaskln/Optim-Visualization

max_forks_repo_name: lukaskln/Optim-Visualization


text:
__precompile__(true)

"""
# Temporal

Time series analysis package for Julia

Type and methods for time series data structures with the `TS` object

Inspired by the pandas DataFrame in python and the xts package in R

"""
module Temporal
    using Dates
    using HTTP
    using JSON
    using Dates
    using Printf
    using RecipesBase
    include("ts.jl")
    include("show.jl")
    include("utils.jl")
    include("indexing/getindex.jl")
    include("indexing/setindex.jl")
    include("indexing/stringrange.jl")
    include("combine.jl")
    include("collapse.jl")
    include("operations.jl")
    include("models.jl")
    include("filter.jl")
    include("ohlc.jl")
    include("convert.jl")
    include("data/utils.jl")
    include("data/yahoo.jl")
    include("data/google.jl")
    include("data/quandl.jl")
    include("data/text.jl")
    include("viz.jl")
    export
        # foundational
        TS, ts,
        # combining
        ojoin, ijoin, ljoin, rjoin, merge, hcat, vcat, head, tail, subset,
        # missingness
        nanrows, nancols, dropnan, fillnan, fillnan!, ffill!, bfill!, linterp!,
        # operations
        numfun, arrfun, operation,
        ones, zeros, trues, falses, isnan, countnz, sign, round,
        sum, mean, maximum, minimum, prod, cumsum, cumprod, diff, lag,
        shift, pct_change,
        rename!, rename,
        # aggregation
        mondays, tuesdays, wednesdays, thursdays, fridays, saturdays, sundays,
        bow, eow, bom, eom, boq, eoq, boy, eoy, collapse, apply,
        # financial
        has_open, has_high, has_low, has_close, has_volume, is_ohlc, is_ohlcv,
        op, hi, lo, cl, vo, ohlc, ohlcv, hlc, hl, hl2, hlc3, ohlc4,
        # utilities
        namefix, namefix!, autoidx, findalphanum,
        SANITIZE_NAMES, RANGE_DELIMITER, set_sanitize_names_option, set_range_delimiter_option,
        # models
        acf,
        # data
        tsread, tswrite, quandl, quandl_auth, quandl_meta, quandl_search, yahoo, google
end


text:
using Dates
using AppliGeneralLedger
using AppliSales # needed for test.jl and runtests.jl

include("./spec.jl")

struct MetaInvoice <: Structure
    _order_id::String
    _training_id::String
    _date::DateTime
    _currency::String
    _currency_ratio::Float64
    # Constructors
    MetaInvoice(order_id, training_id) = new(order_id, training_id, now(), "€", 1.0) #2
    MetaInvoice(order_id, training_id, date, currency, currency_ratio) = new(order_id, training_id, now(), currency, currency_ratio)
end # MetaInvoice

struct Header <: Structure
    _invoice_nbr::String
    _name::String
    _address::String
    _zip::String
    _city::String
    _country::String
    _order_ref::String
    _name_contact::String
    _email_contact::String
end # Header

struct OpentrainingItem <: BodyItem
    _name_training::String
    _date::DateTime
    _price_per_student::Float64
    _students::Array{String, 1}
    _vat_perc::Float64
    # constructors
    OpentrainingItem(name_training, date, price_per_student, students) = new(name_training, date, price_per_student, students, 0.21)
    OpentrainingItem(name_training, date, price_per_student, students, vat_perc) = new(name_training, date, price_per_student, students, vat_perc)
end # OpentrainingItem

struct UnpaidInvoice <: Invoice
    _id::String
    _meta::MetaInvoice
    _header::Header
    _body::OpentrainingItem
end # UnpaidInvoice

struct BankStatement <: Payment
	_date::Date
	_descr::String
	_iban::String
	_amount::Float64
end # BankStatement

struct PaidInvoice <: Invoice
    _id::String
    _meta::MetaInvoice
    _header::Header
    _body::OpentrainingItem
    _stm::BankStatement
end # PaidInvoice

"""
    stm(i::PaidInvoice)

Returns the Bankstatement of a paid invoice.
"""
stm(i::PaidInvoice) = i._stm

# Fields Invoice
meta(i::Invoice)::MetaInvoice = i._meta
header(i::Invoice)::Header = i._header
body(i::Invoice)::BodyItem = i._body
id(i::Invoice)::String = i._id

# fields MetaInvoice
order_id(m::MetaInvoice)::String = m._order_id
training_id(m::MetaInvoice)::String = m._training_id
date(m::MetaInvoice)::DateTime = m._date
currency(m::MetaInvoice)::String = m._currency
currency_ratio(m::MetaInvoice)::Float64 = m._currency_ratio

# fields Header
invoice_nbr(h::Header)::String = h._invoice_nbr
name(h::Header)::String = h._name
address(h::Header)::String = h._address
zip(h::Header)::String = h._zip
city(h::Header)::String = h._city
country(h::Header)::String = h._country
order_ref(h::Header)::String = h._order_ref
name_contact(h::Header)::String = h._name_contact
email_contact(h::Header)::String = h._email_contact

# field OpenTraining
name_training(o::OpentrainingItem)::String = o._name_training
date(o::OpentrainingItem)::DateTime = o._date
price_per_student(o::OpentrainingItem)::Float64 = o._price_per_student
students(o::OpentrainingItem)::Array{String, 1} = o._students
vat_perc(o::OpentrainingItem)::Float64 = o._vat_perc

# gfields BankStatement
date(b::BankStatement)::Date = b._date
descr(b::BankStatement)::String = b._descr
iban(b::BankStatement)::String = b._iban
amount(b::BankStatement)::Float64 = b._amount


text:
"""
    convert2image(array) -> Array{RGB}

Convert the given CIFAR-10 horizontal-major tensor WHCN (or feature
vector/matrix) to a vertical-major HWN `RGB` array.

```julia
julia> CIFAR10.convert2image(CIFAR10.traintensor()) # full training dataset
32×32×50000 Array{RGB{N0f8},3}:
[...]

julia> CIFAR10.convert2image(CIFAR10.traintensor(1)) # first training image
32×32 Array{RGB{N0f8},2}:
[...]
```
"""
function convert2image(array::AbstractArray{T}) where {T<:Number}
    nlast = size(array)[end]
    array = reshape(array, 32, 32, 3, :)
    array = permutedims(array, (3, 2, 1, 4))
    if size(array)[end] == 1 && nlast != 1
        array = dropdims(array, dims=4)
    end
    img = _colorview(RGB, _norm_array(array))

    img
end

_norm_array(array::AbstractArray) = array
_norm_array(array::AbstractArray{<:Integer}) = reinterpret(N0f8, convert(Array{UInt8}, array))


text:
using MPI
using ClimateMachine
using Logging
using Test
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.SystemSolvers
using ClimateMachine.ODESolvers
using LinearAlgebra
using Printf
using Dates
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.VTK: writevtk, writepvtu

if !@isdefined integration_testing
    if length(ARGS) > 0
        const integration_testing = parse(Bool, ARGS[1])
    else
        const integration_testing = parse(
            Bool,
            lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
        )
    end
end

const output = parse(Bool, lowercase(get(ENV, "JULIA_CLIMA_OUTPUT", "false")))

include("advection_diffusion_model.jl")

struct Pseudo1D{n, α, β, μ, δ} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(
    ::Pseudo1D{n, α, β},
    aux::Vars,
    geom::LocalGeometry,
) where {n, α, β}
    # Direction of flow is n with magnitude α
    aux.u = α * n

    # diffusion of strength β in the n direction
    aux.D = β * n * n'
end

function initial_condition!(
    ::Pseudo1D{n, α, β, μ, δ},
    state,
    aux,
    x,
    t,
) where {n, α, β, μ, δ}
    ξn = dot(n, x)
    # ξT = SVector(x) - ξn * n
    state.ρ = exp(-(ξn - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
end
Dirichlet_data!(P::Pseudo1D, x...) = initial_condition!(P, x...)
function Neumann_data!(
    ::Pseudo1D{n, α, β, μ, δ},
    ∇state,
    aux,
    x,
    t,
) where {n, α, β, μ, δ}
    ξn = dot(n, x)
    ∇state.ρ =
        -(
            2n * (ξn - μ - α * t) / (4 * β * (δ + t)) *
            exp(-(ξn - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
        )
end

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model, testname)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state_conservative(model, eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dg, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(pvtuprefix, prefixes, (statenames..., exactnames...))

        @info "Done writing VTK: $pvtuprefix"
    end
end

function run(
    mpicomm,
    ArrayType,
    dim,
    topl,
    N,
    timeend,
    FT,
    dt,
    n,
    α,
    β,
    μ,
    δ,
    vtkdir,
    outputtime,
    linearsolvertype,
    fluxBC,
)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    model = AdvectionDiffusion{dim, fluxBC}(Pseudo1D{n, α, β, μ, δ}())
    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
        direction = EveryDirection(),
    )

    vdg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
        state_auxiliary = dg.state_auxiliary,
        direction = VerticalDirection(),
    )

    Q = init_ode_state(dg, FT(0))

    ode_solver = ARK548L2SA2KennedyCarpenter(
        dg,
        vdg,
        LinearBackwardEulerSolver(linearsolvertype(); isadjustable = false),
        Q;
        dt = dt,
        t0 = 0,
        split_explicit_implicit = false,
    )

    eng0 = norm(Q)
    @info @sprintf """Starting
    norm(Q₀) = %.16e""" eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %.16e
                runtime = %s
                norm(Q) = %.16e""",
                gettime(ode_solver),
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end
    callbacks = (cbinfo,)
    if ~isnothing(vtkdir)
        # create vtk dir
        mkpath(vtkdir)

        vtkstep = 0
        # output initial step
        do_output(
            mpicomm,
            vtkdir,
            vtkstep,
            dg,
            Q,
            Q,
            model,
            "advection_diffusion",
        )

        # setup the output callback
        cbvtk = EveryXSimulationSteps(floor(outputtime / dt)) do
            vtkstep += 1
            Qe = init_ode_state(dg, gettime(ode_solver))
            do_output(
                mpicomm,
                vtkdir,
                vtkstep,
                dg,
                Q,
                Qe,
                model,
                "advection_diffusion",
            )
        end
        callbacks = (callbacks..., cbvtk)
    end

    numberofsteps = convert(Int64, cld(timeend, dt))
    dt = timeend / numberofsteps

    @info "time step" dt numberofsteps dt * numberofsteps timeend

    solve!(
        Q,
        ode_solver;
        numberofsteps = numberofsteps,
        callbacks = callbacks,
        adjustfinalstep = false,
    )

    # Print some end of the simulation information
    engf = norm(Q)
    Qe = init_ode_state(dg, FT(timeend))

    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    base_num_elem = 4

    expected_result = Dict()
    expected_result[2, 1, Float64] = 7.2801198255507391e-02
    expected_result[2, 2, Float64] = 6.8160295851506783e-03
    expected_result[2, 3, Float64] = 1.4439137164205592e-04
    expected_result[2, 4, Float64] = 2.4260727323386998e-06
    expected_result[3, 1, Float64] = 1.0462203776357534e-01
    expected_result[3, 2, Float64] = 1.0280535683502070e-02
    expected_result[3, 3, Float64] = 2.0631857053908848e-04
    expected_result[3, 4, Float64] = 3.3460492914169325e-06
    expected_result[2, 1, Float32] = 7.2801239788532257e-02
    expected_result[2, 2, Float32] = 6.8159680813550949e-03
    expected_result[2, 3, Float32] = 1.4439738879445940e-04
    # This is near roundoff so we will not check it
    # expected_result[2, 4, Float32] = 2.6432753656990826e-06
    expected_result[3, 1, Float32] = 1.0462204366922379e-01
    expected_result[3, 2, Float32] = 1.0280583053827286e-02
    expected_result[3, 3, Float32] = 2.0646647317335010e-04
    expected_result[3, 4, Float32] = 2.0226731066941284e-05

    numlevels = integration_testing ? 4 : 1

    @testset "$(@__FILE__)" begin
        for FT in (Float64, Float32)
            result = zeros(FT, numlevels)
            for dim in 2:3
                for fluxBC in (true, false)
                    for linearsolvertype in (SingleColumnLU, ManyColumnLU)
                        d = dim == 2 ? FT[1, 10, 0] : FT[1, 1, 10]
                        n = SVector{3, FT}(d ./ norm(d))

                        α = FT(1)
                        β = FT(1 // 100)
                        μ = FT(-1 // 2)
                        δ = FT(1 // 10)
                        for l in 1:numlevels
                            Ne = 2^(l - 1) * base_num_elem
                            brickrange = (
                                ntuple(
                                    j -> range(
                                        FT(-1);
                                        length = Ne + 1,
                                        stop = 1,
                                    ),
                                    dim - 1,
                                )...,
                                range(FT(-5); length = 5Ne + 1, stop = 5),
                            )

                            periodicity = ntuple(j -> false, dim)
                            topl = StackedBrickTopology(
                                mpicomm,
                                brickrange;
                                periodicity = periodicity,
                                boundary = (
                                    ntuple(j -> (1, 2), dim - 1)...,
                                    (3, 4),
                                ),
                            )
                            dt = (α / 4) / (Ne * polynomialorder^2)

                            outputtime = 0.01
                            timeend = 0.5

                            @info (
                                ArrayType,
                                FT,
                                dim,
                                linearsolvertype,
                                l,
                                fluxBC,
                            )
                            vtkdir = output ?
                                "vtk_advection" *
                            "_poly$(polynomialorder)" *
                            "_dim$(dim)_$(ArrayType)_$(FT)" *
                            "_$(linearsolvertype)_level$(l)" :
                                nothing
                            result[l] = run(
                                mpicomm,
                                ArrayType,
                                dim,
                                topl,
                                polynomialorder,
                                timeend,
                                FT,
                                dt,
                                n,
                                α,
                                β,
                                μ,
                                δ,
                                vtkdir,
                                outputtime,
                                linearsolvertype,
                                fluxBC,
                            )
                            # test the errors significantly larger than floating point epsilon
                            if !(dim == 2 && l == 4 && FT == Float32)
                                @test result[l] ≈
                                      FT(expected_result[dim, l, FT])
                            end
                        end
                        @info begin
                            msg = ""
                            for l in 1:(numlevels - 1)
                                rate = log2(result[l]) - log2(result[l + 1])
                                msg *= @sprintf(
                                    "\n  rate for level %d = %e\n",
                                    l,
                                    rate
                                )
                            end
                            msg
                        end
                    end
                end
            end
        end
    end
end

nothing

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.7516814527569986 + 0.07211062597860263*1/(1 + exp(-(0.800956367833014 + 0.1575422515944389*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.30181017969165413*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.5024447348308962*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + -0.7463877606567912*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.28099609060437203*1/(1 + exp(-(-0.8241959750660866 + 0.48752067922981546*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + -0.8942226763869914*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + -0.1882956347370066*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.0613423792438379*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + 0.19401527858092393*1/(1 + exp(-(-0.014178711859067938 + -0.6496362145686132*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5744456236018198*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.06990002879133872*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.11738131892224501*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.43991506890783594*1/(1 + exp(-(-0.38559964579408224 + -0.6882459841770032*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5869261088566509*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.7286086900917508*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.9385542676650069*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))))))))) + 1/(1 + exp(-(0.5442086655639207 + -0.859033717055464*1/(1 + exp(-(0.800956367833014 + 0.1575422515944389*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.30181017969165413*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.5024447348308962*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + -0.7463877606567912*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.6994107397450926*1/(1 + exp(-(-0.8241959750660866 + 0.48752067922981546*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + -0.8942226763869914*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + -0.1882956347370066*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.0613423792438379*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.21447765318835765*1/(1 + exp(-(-0.014178711859067938 + -0.6496362145686132*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5744456236018198*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.06990002879133872*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.11738131892224501*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + 0.9247182813449877*1/(1 + exp(-(-0.38559964579408224 + -0.6882459841770032*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5869261088566509*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.7286086900917508*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.9385542676650069*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))))))))) + 1/(1 + exp(-(-0.33457802196328057 + 0.8511359109817849*1/(1 + exp(-(0.800956367833014 + 0.1575422515944389*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.30181017969165413*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.5024447348308962*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + -0.7463877606567912*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + 0.9056512987027685*1/(1 + exp(-(-0.8241959750660866 + 0.48752067922981546*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + -0.8942226763869914*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + -0.1882956347370066*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.0613423792438379*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.901983268277478*1/(1 + exp(-(-0.014178711859067938 + -0.6496362145686132*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5744456236018198*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.06990002879133872*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.11738131892224501*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.02854306019224717*1/(1 + exp(-(-0.38559964579408224 + -0.6882459841770032*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5869261088566509*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.7286086900917508*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.9385542676650069*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))))))))) + 1/(1 + exp(-(0.85223750045835 + -0.7468516096711739*1/(1 + exp(-(0.800956367833014 + 0.1575422515944389*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.30181017969165413*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.5024447348308962*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + -0.7463877606567912*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.1269236474608202*1/(1 + exp(-(-0.8241959750660866 + 0.48752067922981546*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + -0.8942226763869914*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + -0.1882956347370066*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.0613423792438379*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.5407492362268407*1/(1 + exp(-(-0.014178711859067938 + -0.6496362145686132*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5744456236018198*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.06990002879133872*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.11738131892224501*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + -0.8670343803826994*1/(1 + exp(-(-0.38559964579408224 + -0.6882459841770032*1/(1 + exp(-(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])))) + 0.5869261088566509*1/(1 + exp(-(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])))) + 0.7286086900917508*1/(1 + exp(-(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])))) + 0.9385542676650069*1/(1 + exp(-(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""

    u,v,η,sst = run_model()

runs ShallowWaters with default parameters as defined in src/DefaultParameters.jl

# Examples
```jldoc
julia> u,v,η,sst = run_model(Float64,nx=200,output=true)
```
"""
function run_model(::Type{T}=Float32;     # number format
    kwargs...                            # all additional parameters
    ) where {T<:Real}

    P = Parameter(T=T;kwargs...)
    return run_model(T,P)
end

function run_model(P::Parameter)
    @unpack T = P
    return run_model(T,P)
end

function run_model(::Type{T},P::Parameter) where {T<:Real}

    @unpack Tprog = P

    G = Grid{T,Tprog}(P)
    C = Constants{T,Tprog}(P,G)
    F = Forcing{T}(P,G)
    S = ModelSetup{T,Tprog}(P,G,C,F)

    Prog = initial_conditions(Tprog,S)
    Diag = preallocate(T,Tprog,G)

    Prog = time_integration(Prog,Diag,S)
    return Prog
end


text:
@testset "Verified linear solver" begin

    n = 4
    Arat = reshape(1 .//(1:n^2), n, n)
    brat = Arat*fill(1//1, n)
    Afloat = float(Arat)
    bfloat = float(brat)

    x, cert = epsilon_inflation(Afloat, bfloat)

    @test all(ones(n) .∈ x)
    @test cert

    Ain = convert.(IA.Interval{Float64}, IA.Interval.(Arat, Arat))
    bin = convert.(IA.Interval{Float64}, IA.Interval.(brat, brat))

    x, cert = epsilon_inflation(Ain, bin)

    @test all(ones(n) .∈ x)
    @test cert

    # big float test
    Abig = BigFloat.(Arat)
    bbig = BigFloat.(brat)

    x, cert = epsilon_inflation(Abig, bbig)

    @test cert
    @test all(diam.(x) .< 1e-50)
    @test all(ones(n) .∈ x)

    # case when should not be possible to certify
    A = [1..2 1..4;0..1 0..1]
    b = A*[1; 1]

    _, cert = epsilon_inflation(A, b)
    @test !cert
end


text:
using Test
using ImagineHardware
using Unitful
import Unitful: V, s, Hz, m, μm, A, mA, N, g, μF

p = retrieve_positioner("NanoSX800")
@test name(p) == "NanoSX800"
@test vmin_in(p) == -20.0V
@test vmax_in(p) == 130.0V
@test vspan_in(p) == vmax_in(p) - vmin_in(p)
@test ImagineHardware.mon_slope(p, true) == 10.0V / 800.0μm
@test ImagineHardware.mon_slope(p, false) == 10.0V / 900.0μm
@test closed_loop_pad(p) == 50.0μm
@test vmin_out(p, true) == 0.0V
@test vmax_out(p, true) == 10.0V
@test vmin_out(p, true) == 0.0V
@test vmax_out(p, true) == 10.0V
@test vmin_out(p, false) == 0.0V
@test vmax_out(p, false) == 10.0V
@test vspan_out(p, true) == vmax_out(p, true) - vmin_out(p, true)
@test vspan_out(p, false) == vmax_out(p, false) - vmin_out(p, false)
@test pmin(p, true) == 0.0μm
@test pmax(p, true) == 800.0μm
@test pmin(p, false) == -50.0μm
@test pmax(p, false) == 850.0μm
@test pspan(p, true) == pmax(p, true) - pmin(p, true)
@test pspan(p, false) == pmax(p, false) - pmin(p, false)
@test ImagineHardware.blocked_pad(p) == 125μm
@test closed2open(0.0μm, p) == 0.0μm
@test closed2open(800.0μm, p) == 800.0μm
#@test_throws Exception closed2open(-10.0μm, p) #relaxed this requirement because the piezo can go this far in closed loop mode, though it's not supposed to
@test open2closed(0.0μm, p) == 0.0μm
@test open2closed(800.0μm, p) == 800.0μm
@test_throws Exception open2closed(-10.0μm, p)
@test max_displacement(p) == 1150μm
@test capacitance(p) == 7.0μF
@test motion_axis(p) == 1
@test stiffness(p) == stiffness(p, 1) == 0.2N/μm
@test stiffness_all(p) == (0.2N/μm, 2.5N/μm, 2.5N/μm)
@test effective_mass(p) == 77.1917g
@test resonance(p) == resonance(p, 1) == 210.0Hz
@test resonance_all(p) == (210.0Hz, 500.0Hz, 700.0Hz)

@test mon2pos(0.0V, p, true) == 0.0μm
@test mon2pos(0.0V, p, false) == -50.0μm
@test mon2pos(10.0V, p, true) == 800.0μm
@test mon2pos(10.0V, p, false) == 850.0μm
@test_throws Exception mon2pos(10.1V, p, false)
@test_throws Exception mon2pos(10.1V, p, true)
@test pos2mon(0.0μm, p, true) == 0.0V
@test pos2mon(-50.0μm, p, false) == 0.0V
@test pos2mon(800.0μm, p, true) == 10.0V
@test pos2mon(850.0μm, p, false) == 10.0V
@test_throws Exception pos2mon(1000μm, p, false)
#@test_throws Exception pos2mon(808μm, p, true) #relaxed this requirement becase the piezo can go this far in closed loop mode, though it's not supposed to

@test displacement2pos(0μm, p, false) == -50μm #blocked internally
@test displacement2pos(125μm, p, false) == -50μm
@test displacement2pos(125μm, p, false) == -50μm
@test displacement2pos(1025μm, p, false) == 850μm
@test displacement2pos(175μm, p, true) == 0μm
@test displacement2pos(975μm, p, true) == 800μm
@test_throws Exception displacement2pos(0μm, p, true)
@test_throws Exception displacement2pos(125μm, p, true)
@test_throws Exception displacement2pos(1025μm, p, true)
#including outside force
@test displacement2pos(125μm, p, false, 0.0N) == displacement2pos(125μm, p, false)
@test displacement2pos(125μm, p, false, 10.0N) == displacement2pos(125μm + uconvert(μm, 10.0N/stiffness(p)), p, false)

@test_throws Exception pos2displacement(-50μm, p, false)
@test_throws Exception pos2displacement(850μm, p, false)
@test pos2displacement(0μm, p, true) == 175μm
@test pos2displacement(800μm, p, true) == 975μm
@test pos2displacement(400μm, p, true) == 575μm
@test pos2displacement(400μm, p, false) == 575μm
#including outside force
@test pos2displacement(400μm, p, true, 0.0N) == 575μm
@test pos2displacement(400μm, p, false, 0.0N) == 575μm
@test pos2displacement(400μm, p, true, 10.0N) == 575μm - 10.0N/stiffness(p)
@test pos2displacement(400μm, p, true, 10.0N) == 575μm - 10.0N/stiffness(p)

@test input2displacement(-20.0V, p) == 0.0μm
@test input2displacement(130.0V, p) == max_displacement(p)
@test input2displacement(55.0V, p) == max_displacement(p) / 2
@test_throws Exception input2displacement(-30.0V, p)
@test_throws Exception input2displacement(140.0V, p)

@test displacement2input(0.0μm, p) == -20.0V
@test displacement2input(max_displacement(p) / 2, p) == 55.0V
@test displacement2input(max_displacement(p), p) == 130.0V
@test_throws Exception displacement2input(1200μm, p)
@test_throws Exception displacement2input(-10μm, p)

@test input2pos(-20.0V, p, false) == displacement2pos(0.0μm, p, false)
@test input2pos(-20.0V, p, false) == displacement2pos(0.0μm, p, false)
@test input2pos(130.0V, p, false) == displacement2pos(1150.0μm, p, false)
@test input2pos(130.0V, p, false) == displacement2pos(1150.0μm, p, false)
@test input2pos(55.0V, p, false) == displacement2pos(max_displacement(p)/2, p, false)
@test input2pos(55.0V, p, true) == displacement2pos(max_displacement(p)/2, p, true)
@test_throws Exception input2pos(-30.0V, p)
@test_throws Exception input2pos(140.0V, p)
#including outside force
@test input2pos(-20.0V, p, false, 0.0N) == input2pos(-20.0V, p, false)
@test input2pos(130.0V, p, false, 0.0N) == input2pos(130.0V, p, false)
@test input2pos(55.0V, p, false, 0.0N) == input2pos(55.0V, p, false)
@test input2pos(55.0V, p, true, 0.0N) == input2pos(55.0V, p, true)
@test input2pos(-20.0V, p, false, 10.0N) == displacement2pos(0.0μm, p, false, 10.0N)

#including outside force
@test pos2input(0μm, p, true, 0.0N) == pos2input(0μm, p, true)
@test pos2input(800μm, p, true, 0.0N) == pos2input(800μm, p, true)
@test pos2input(0μm, p, true, 10.0N) == displacement2input(pos2displacement(0μm, p, true, 10.0N), p)

#adjusted_displacement - compression = mon_displacement
#adjusted_displacement = mon_displacement + compression
#                  = mon_displacement + f / stiffness(p)
mon2displacement(0.0V, p, true) == pos2displacement(0.0μm, p, true)
@test_throws Exception mon2displacement(0.0V, p, false)
@test mon2displacement(5.0V, p, true) == max_displacement(p)/2
@test mon2displacement(0.0V, p, true, 0.0N) == pos2displacement(0.0μm, p, true, 0.0N)
@test_throws Exception mon2displacement(0.0V, p, false, 0.0N)
@test mon2displacement(0.0V, p, true, 1.0N) == pos2displacement(0.0μm, p, true, 1.0N)
@test_throws Exception mon2displacement(0.0V, p, false, 1.0N)

@test max_force(p, 0μm, 1) == stiffness(p) * max_displacement(p)
@test max_force(p, 800μm, -1) == stiffness(p) * -800μm # - pos2displacement(800.0μm, p, true))
@test max_force(p, 0μm, 800μm, true) == max_force(p, pos2displacement(0.0μm, p, true), 1)
@test max_force(p, 800μm, 0μm, true) == max_force(p, pos2displacement(800.0μm, p, true), -1) == -max_force(p, pos2displacement(0.0μm, p, true), 1)
@test max_force(p, 0μm, 800μm, false) == max_force(p, pos2displacement(0.0μm, p, false), 1)
@test max_force(p, 800μm, 0μm, false) == max_force(p, pos2displacement(800.0μm, p, false), -1) == -max_force(p, pos2displacement(0.0μm, p, false), 1)

@test loaded_resonance(100.0Hz, 70.0g, 100.0g) == uconvert(Unitful.Hz, 100.0Hz * sqrt(70.0g/(70.0g+100.0g)))
@test loaded_resonance(p, 500g) == loaded_resonance(resonance(p), effective_mass(p), 500g)


text:
"""
    reducedForm(qq, qrows, qcols, bcols, condn)

Compute reduced-form coefficient matrix, b.
"""
function reducedForm(qq::Array{Float64,2}, qrows::Int64, qcols::Int64, bcols::Int64, neq::Int64, condn::Float64)

    bb = zeros(qrows, bcols)

    left = 1 : (qcols - qrows)
    right = (qcols - qrows + 1) : qcols
    qtmp = similar(qq)

    nonsing = ( 1 / cond(qq[:, right], 1) ) > condn
     if nonsing
        qtmp[:, left] = -qq[:, right] \ qq[:, left]
        bb = qtmp[1 : neq, 1 : bcols]
    else  # rescale by dividing row by maximal qr element
        # 'inverse condition number small, rescaling '
        themax = maximum(abs.(qtmp[:, right]), 2)
        oneover = diagm(1 ./ themax[:, 1])

        nonsing = ( 1 / cond( oneover*qtmp[:, right], 1) ) > condn
        if nonsing
            qtmp[:, left] = -(oneover*qtmp[:, right]) \ (oneover*(qtmp[:, left]))
            bb = qtmp[1:neq, 1:bcols]
        end
    end

    return (nonsing, bb)

end # reducedForm


text:
"""
Used to plot macro geographic map as vector. From the free data source: https://www.diva-gis.org/data

kwarg "adminlevel" refers to the available shape files in the DIVA dataset:
    0... country outlines,
    1... regional (e.g. home nations in uk case),
    2... county
"""

function DIVA(;
    territory::String="GBR",
    exclShapeBelowLen::Int=200,
    adminLevel=0)

    # Paths
    env = loadPaths()
    subDir = "$(territory)_adm"
    shapeFile = "$(territory)_adm$(adminLevel).shp"
    path = joinpath(env.paths[:basemap_DIVA],subDir,shapeFile)

    # Load shapefile
    shapeSource = open(path, "r") do io
        read(io, Shapefile.Handle)
    end

    # Main loop
    traces = [adminLevelTraces(
        shapeSource.shapes[iShape],
        exclShapeBelowLen,
        territory)
    for iShape in 1:length(shapeSource.shapes)]

    # Merge trace vectors
    traces = vcat(traces...)

    # Layout scale ratio constrained based on midLat
    allLat = vcat([trace[:y] for trace in traces]...)
    midLat = (maximum(allLat)+minimum(allLat))/2
    latLonRatio = cos(midLat*π/180)
    layout = Layout(
        showlegend=false,
        xaxis=attr(zeroline=false,
            scaleanchor="y", scaleratio=latLonRatio),
        hovermode="closest"
    )

    # Plot (suppress in application with semi-colon)
    plt = plot(traces,layout);

    return plt
end

function adminLevelTraces(shapeObj,exclShapeBelowLen,territory)
    shape = GeoInterface.coordinates(shapeObj)

    # Remove shapes with boundaries shorter than exclShapeBelowLen::Int
    shape = [shp for shp in shape if length(shp[1]) > exclShapeBelowLen]

    # Set boundary name
    shapeNames = "$(territory) basemap"

    # Build traces of all polygons in shape (i.e. separate islands)
    coordsets = [hcat(shp[1]...) for shp in shape]
    traces = [scatter(
            x=coordset[1,:][:],
            y=coordset[2,:][:],
            name=shapeNames,mode="lines",line=attr(width=0.1,color="#444444BB"),
            fill="tonexty",fillcolor="#58B19F10",
            hoverinfo="skip",hovertemplate=nothing)
    for coordset in coordsets]

    return traces
end


text:
function triangular_system(dict::Dict{Symbol,Union{Expr,Float64, Int64}})
    context = Dict() # context
    solutions = Dict{Symbol, Float64}()
    finished = false
    N = length(dict)
    n = 0
#     return
    while ~(  (finished) || n>N )
        done_smthg = false
        n += 1
        for k in keys(dict)
            if ~(in(k,keys(solutions)))
                expr = dict[k]
                try
                    sol = eval( :(let $([:($x=$y) for (x, y) in solutions]...); $expr end) )
                    sol = Float64(sol)
                    solutions[k] = sol
                    context[k] = sol
                    done_smthg = true
                catch
                    0
                end
            end
        end
        if done_smthg==false
            finished=true
        end
    end
    if length(solutions)<length(dict)
        error("Not a triangular system")
#         throw(Exception())
    else
        return solutions
    end
end

function string_dict_to_symbol_dict(dict::Dict{String, String})
    res = Dict{Symbol, Union{Expr,Float64}}()
    for k in keys(dict)
        v = dict[k]
        res[Symbol(k)] = parse(v)
    end
    return res
end

function triangular_system(dict::Dict{String, String})
    return triangular_system(string_dict_to_symbol_dict(dict))
end


text:
using BenchmarkTools

function answer1(input = "input0")
    total = 0
    for line in readlines(input)
        code, output = split(line," | ")
        outputs = split(output," ")
        for i in [2,3,4,7]
            total += sum(sizeof.(outputs) .== i)
        end
    end
    total
end

@show answer1("input")

"I don't want to compare entire arrays, thus this function can give 9 or 3 with incorrect input"
function screen(decoded)
    if sum(decoded) == 2
        return 1
    elseif sum(decoded) == 3
        return 7
    elseif sum(decoded) == 4
        return 4
    elseif sum(decoded) == 7
        return 8
    elseif sum(decoded) == 5
        if decoded[2] == 1
            return 5
        elseif decoded[5] == 1
            return 2
        else
            return 3
        end
    elseif sum(decoded) == 6
        if decoded[4] == 0
            return 0
        elseif decoded[3] == 0
            return 6
        else
            return 9
        end
    end
    -1
end
"Takes vector and returns the corresponding array for input to screen."
function todecodedarray(vectortodecode::Vector{Char})
    answer = zeros(Int, 7)
    indexes = Int.(vectortodecode .- 96)
    answer[indexes] .= 1
    answer
end

"Clearly not the most elegant way. I solved it on paper and then translated the operations to Julia.
it is probably easier to keep operating until you have isolated all the numbers, but I got first the
function relating the mixed up input with the actual input"
function answer2(input="input0")
    i = 1
    total = 0
    for line in readlines(input)
        key = repeat(['.'],7)
        real_key = ['a','b','c','d','e','f','g']
        code, output = split(line," | ")
        outputs = split(output," ")
        codes = split.(code," ")
        s1 = only.(split(codes[length.(codes) .== 2][1], ""))
        s3 = [""]
        s4 = only.(split(codes[length.(codes) .== 4][1], ""))
        s6 = [""]
        s7 = only.(split(codes[length.(codes) .== 3][1], ""))
        s8 = only.(split(codes[length.(codes) .== 7][1], ""))

        key[1] = setdiff(s7,s1)[1]
        codes6 = codes[length.(codes) .== 6]
        for i in eachindex(codes6)
            next = mod((i+1),3) + 1
            prev = mod((i+2),3) + 1
            i = mod(i,3) + 1
            code = codes6[i]
            compare = setdiff(s8,only.(split(code, "")))
            if length( compare ) == 1 && compare[1] ∈ s1
                key[3] = compare[1]
                key[6] = setdiff(s1,compare)[1]
                s6 = split(code, "")
                codes6 = [codes6[next],codes6[prev]]
                code
                break
            end
        end
        codes5 = codes[length.(codes) .== 5]

        for i in eachindex(codes5)
            next = mod((i+1),3) + 1
            prev = mod((i+2),3) + 1
            i = mod(i,3) + 1
            codes5[i],codes5[next],codes5[prev]
            joined = union(setdiff(codes5[i],codes5[next]),setdiff(codes5[i],codes5[prev]))
            if sum(in(joined).(s1)) == 2
                s3 = codes5[i]
                codes5 = [codes5[next],codes5[prev]]
                break
            end
        end
        key[4] = setdiff(intersect(s3,s4),s1)[1]
        key[7] = setdiff(s3,key)[1]
        key[2] = setdiff(s4,key)[1]
        key[5] = setdiff(real_key,key)[1]
        decryptDict = Dict(key .=> real_key)
        decrypt(c) = decryptDict[c]
        for (i,encryptedNumber) in enumerate(outputs)
            total += screen(todecodedarray(map(decrypt,only.(split(encryptedNumber,""))))) * 10 ^(4-i)
        end
    end
    total
end

@show answer2("input")


text:
# Arvind Prasadan
# 2019 July
# DMF Package

"""
    permute_IP(Q::Matrix, Q_hat::Matrix)

# Inputs
- Q::Matrix: true vectors
- Q_hat::Matrix: estimate of Q

# Outputs
- Vector of 'best' inner products
"""
function permute_IP(Q, Q_hat)
    k = size(Q, 2)

    @assert(k == size(Q_hat, 2), "Q and Q_hat have differing numbers of vectors")
    @assert(size(Q, 1) == size(Q_hat, 1), "Q and Q_hat are not compatible")

    # Normalize for unit norm
    Q = mapslices(normalize, Q; dims = 1)
    Q_hat = mapslices(normalize, Q_hat; dims = 1)

    # Generate permutations
    pt = permutations(1:k) |> collect

    # Find maximum sum of inner products
    pt_idx = argmax([sum(diag(abs.(Q' * Q_hat[:, x]))) for x in pt])
    pt = pt[pt_idx]

    return diag(Q' * Q_hat[:, pt])
end


text:
using Documenter, RemoteS

makedocs(
    modules = [RemoteS],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true", assets = ["assets/custom.css"]),
    authors = "Joaquim Luis",
    sitename = "RemoteS.jl",
    pages = Any[
        "Gallery"                  => [
            "Aqua orbits"          => "gallery/Aqua_orbits/remotes_sat_tracks.md",
            "Aqua SST"             => "gallery/Aqua_sst/remotes_L2_SST.md",
            "Cloud-Native HLS dat" => "gallery/HLS/cloud-native-hls-data.md",
            "Landsat 8 images"     => "gallery/L8cube_img/remotes_L8_cube_img.md",
            "Landsat 8 NDVI"       => "gallery/L8cube_ndvi/remotes_L8_NDVI.md",
        ],
        "Index"                    => "index.md"],

    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/GenericMappingTools/RemoteS.jl.git",
	target  = "build",
    devbranch = "main",
    devurl = "dev",
    #versions = ["stable" => "v^", "v#.#", devurl => devurl],
    push_preview = true
)

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
using RunTests
using Base.Test

@testmodule ExtraStatisticsTest begin
  using MLLibrary

  function test_meanColumn()
    A = [1 2 3 4;5 6 7 8]
    @test MLLibrary.Statistics.meanColumn(A)==[3.0 4.0 5.0 6.0]
    # @test 1==0
  end

  function test_stdColumn()
    A = [1 2 3 4;5 6 7 8]
    @test MLLibrary.Statistics.stdColumn(A)==[2.8284271247461903 2.8284271247461903 2.8284271247461903 2.8284271247461903]
    # @test 1==0
  end

  function test_sumColumn()
    A = [1 2 3 4;5 6 7 8]
    @test MLLibrary.Statistics.sumColumn(A)==[6 8 10 12]
    # @test 1==0
  end

  function test_getMinIndices()
    A=[1;2;-1;5;3;2;9;-3; 10]
    @test MLLibrary.Statistics.getMinIndices(A)==[8]
    @test MLLibrary.Statistics.getMinIndices(A,2)==[8, 3]
    @test MLLibrary.Statistics.getMinIndices(A,3)==[8, 3, 1]
    @test MLLibrary.Statistics.getMinIndices(A,4)==[8, 3, 1, 2]
  end

  function test_getMajority()
    A=[1;2;-1;5;3;2;9;-3;10]
    @test MLLibrary.Statistics.getMajority(A)==2
    A=[1;2;1;5;3;2;9;-3;10]
    maj=MLLibrary.Statistics.getMajority(A)
    @test maj==1 || maj==2
    A=[1]
    @test MLLibrary.Statistics.getMajority(A)==1

  end
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
# A./*B

FAhadamarddivB(A,B)=(A./B,nothing)

function FAhadamarddivB_inplace(value,aux,A,B)
    copy!(value,A./B)
end

function DAhadamarddivB(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,B)
    if derivativeIDX==1
        axpy!(1.0,grad_c./B,grad_n) # should really use BLAS here
    elseif derivativeIDX==2
        axpy!(-1.0,grad_c.*A./(B.*B),grad_n)
    end
end

if PROC=="GPU"
    function FAhadamarddivB(A::CudaArray,B::CudaArray)
        tmp=CudaArray(Float64,size(A))
        vdiv!(1.0,A,B,tmp)
        return (tmp,nothing)
    end

    function FAhadamarddivB_inplace(value::CudaArray,aux,A::CudaArray,B::CudaArray)
        vdiv!(1.0,A,B,value)
    end

    function DAhadamarddivB(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,B::CudaArray)
    if derivativeIDX==1
        vdivupdate!(1.0,grad_c,B,grad_n)
    elseif derivativeIDX==2
        vAoverBupdate!(1.0,grad_c,A,B,grad_n)
    end
end

end

Derivative[FAhadamarddivB]=DAhadamarddivB
Inplace[FAhadamarddivB]=FAhadamarddivB_inplace

import Base. ./
./(A::ADnode,B::ADnode)=ADnode(FAhadamarddivB,[A B])
export ./

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff