text:
struct LinearStateSpaceModel{T}
    A::Matrix{T}
    B::Matrix{T}
    Q::Matrix{T}
end

function LinearStateSpaceModel(A::Matrix, Q::Matrix)
    LinearStateSpaceModel(A,zeros(size(A,1),1),Q)
end

struct Simulator <: Transducer
    lssm::LinearStateSpaceModel
end

function Transducers.start(rf::R_{Simulator}, result)
    private_state = (0, zeros(size(xform(rf).lssm.A,1)))
    return wrap(rf, private_state, start(inner(rf), result))
end

function step(l::LinearStateSpaceModel, x, input)
    l.A*x + l.B*input + Gaussian(l.Q)
end

function Transducers.next(rf::R_{Simulator}, result, input)
    wrapping(rf, result) do st, iresult
        (n, x) = st
        n += 1
        new_x = step(xform(rf).lssm,x,input)
        iresult = next(inner(rf), iresult, new_x)
        return (n, new_x), iresult
    end
end


text:
@with_kw struct OrthogonalRegularizer
    β::Float32 = 1f-4 # regularization coefficient
end

function (R::OrthogonalRegularizer)(π)
    reg = 0f0
    dev = Zygote.ignore(()->device(π))
    for l in Zygote.ignore(()->filter((l)->hasproperty(l, :weight), layers(π)))
        W = to2D(l.weight)
        prod = W' * W
        mat = Zygote.ignore(() -> ones(Float32, size(prod)...) .- Matrix{Float32}(I, size(prod)...) |> dev)
        reg += norm(prod .* mat)^2
    end
    R.β*reg
end

ext: jl


text:
module Rhyolite

import Vulkan
const Vk = Vulkan

using ResultTypes

using TimerOutputs
using MLStyle
using Dictionaries
using UnPack
using Accessors

using SPIRV: ImageType
using SPIRV

import glslang_jll
const glslangValidator = glslang_jll.glslangValidator(identity)

const to = TimerOutput()
const Optional{T} = Union{T,Nothing}

const debug_callback_c = Ref{Ptr{Cvoid}}(C_NULL)
const debug_messenger = Ref{Vk.DebugUtilsMessengerEXT}()

function __init__()
    # for debugging in Vulkan
    debug_callback_c[] =
        @cfunction(default_debug_callback, UInt32, (DebugUtilsMessageSeverityFlagEXT, DebugUtilsMessageTypeFlagEXT, Ptr{vk.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))
end

include("utils.jl")
include("handle_wrappers.jl")
include("dispatch.jl")
include("init.jl")
include("memory.jl")
include("commands.jl")
include("frames.jl")

include("shaders/dependencies.jl")
include("shaders/resources.jl")
include("shaders/vertex.jl") # type piracy
include("shaders/formats.jl")
include("shaders/specification.jl")
include("shaders/source.jl")
include("shaders/compilation.jl")

include("descriptor_allocator.jl")
include("descriptor_set_layouts.jl")

include("pipelines/binding.jl")
include("pipelines/cache.jl")

include("draw/render_set.jl")

include("passes.jl")
include("device.jl")

export
        # handles
        AbstractHandle,
        handle,
        Allocated,
        memory,
        Created,
        info,

        # initialization
        init,
        init_debug,
        debug_messenger,

        # setup introspection
        require_extension,
        require_feature,
        require_layer,

        # memory
        download_data,
        upload_data,
        buffer_size,
        MemoryDomain,
        find_memory_type,
        LinearAllocator,
        reset!,

        # descriptor allocator
        DescriptorAllocator,
        find_pool!,
        allocate_descriptor_sets!,
        free_descriptor_sets!,

        # descriptor set layout cache
        DescriptorSetLayoutCache,

        # queue dispatch
        QueueDispatch,
        queue_infos,
        submit,
        present,
        queue_family_indices,

        # frames
        Frame,
        FrameSynchronization,
        FrameState,
        wait_hasrendered,
        next_frame!,
        command_buffers,

        # commands
        @record,
        ThreadedCommandPool,

        # shaders
        ShaderDependencies,
        ShaderResource,
        ShaderLanguage, SPIR_V, GLSL, HLSL,
        ShaderSpecification,
        ShaderCache, find_shader!, find_source!,
        Shader,
        SampledImage,
        StorageBuffer,
        collect_bindings,
        create_descriptor_set_layouts,
        create_descriptor_set_layouts!,
        vertex_input_attribute_descriptions,

        # pipeline
        BindRequirements,
        BindState,
        GraphicsPipelineCache,
        PipelineLayoutCache,

        # render sets
        RenderSet,
        RenderInfo,
        prepare!,

        # passes
        PassType,
        OpaqueShadingPass,
        LightingPass,
        TransparentShadingPass

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct SummarySize
    seen::IdDict{Any,Any}
    frontier_x::Vector{Any}
    frontier_i::Vector{Int}
    exclude::Any
    chargeall::Any
end

"""
    Base.summarysize(obj; exclude=Union{...}, chargeall=Union{...}) -> Int

Compute the amount of memory, in bytes, used by all unique objects reachable from the argument.

# Keyword Arguments
- `exclude`: specifies the types of objects to exclude from the traversal.
- `chargeall`: specifies the types of objects to always charge the size of all of their
  fields, even if those fields would normally be excluded.
"""
function summarysize(obj;
                     exclude = Union{DataType, Core.TypeName, Core.MethodInstance},
                     chargeall = Union{Core.TypeMapEntry, Method})
    @nospecialize obj exclude chargeall
    ss = SummarySize(IdDict(), Any[], Int[], exclude, chargeall)
    size::Int = ss(obj)
    while !isempty(ss.frontier_x)
        # DFS heap traversal of everything without a specialization
        # BFS heap traversal of anything with a specialization
        x = ss.frontier_x[end]
        i = ss.frontier_i[end]
        val = nothing
        if isa(x, SimpleVector)
            nf = length(x)
            if isassigned(x, i)
                val = x[i]
            end
        elseif isa(x, Array)
            nf = length(x)
            if ccall(:jl_array_isassigned, Cint, (Any, UInt), x, i - 1) != 0
                val = x[i]
            end
        else
            nf = nfields(x)
            ft = typeof(x).types
            if !isbitstype(ft[i]) && isdefined(x, i)
                val = getfield(x, i)
            end
        end
        if nf > i
            ss.frontier_i[end] = i + 1
        else
            pop!(ss.frontier_x)
            pop!(ss.frontier_i)
        end
        if val !== nothing && !isa(val, Module) && (!isa(val, ss.exclude) || isa(x, ss.chargeall))
            size += ss(val)::Int
        end
    end
    return size
end

(ss::SummarySize)(@nospecialize obj) = _summarysize(ss, obj)
# define the general case separately to make sure it is not specialized for every type
@noinline function _summarysize(ss::SummarySize, @nospecialize obj)
    isdefined(typeof(obj), :instance) && return 0
    # NOTE: this attempts to discover multiple copies of the same immutable value,
    # and so is somewhat approximate.
    key = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)
    haskey(ss.seen, key) ? (return 0) : (ss.seen[key] = true)
    if nfields(obj) > 0
        push!(ss.frontier_x, obj)
        push!(ss.frontier_i, 1)
    end
    if isa(obj, UnionAll) || isa(obj, Union)
        # black-list of items that don't have a Core.sizeof
        sz = 2 * sizeof(Int)
    else
        sz = Core.sizeof(obj)
    end
    if sz == 0
        # 0-field mutable structs are not unique
        return gc_alignment(0)
    end
    return sz
end

(::SummarySize)(obj::Symbol) = 0
(::SummarySize)(obj::SummarySize) = 0

function (ss::SummarySize)(obj::String)
    key = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)
    haskey(ss.seen, key) ? (return 0) : (ss.seen[key] = true)
    return Core.sizeof(Int) + Core.sizeof(obj)
end

function (ss::SummarySize)(obj::DataType)
    key = pointer_from_objref(obj)
    haskey(ss.seen, key) ? (return 0) : (ss.seen[key] = true)
    size::Int = 7 * Core.sizeof(Int) + 6 * Core.sizeof(Int32)
    size += 4 * nfields(obj) + ifelse(Sys.WORD_SIZE == 64, 4, 0)
    size += ss(obj.parameters)::Int
    if isdefined(obj, :types)
        size += ss(obj.types)::Int
    end
    return size
end

function (ss::SummarySize)(obj::Core.TypeName)
    key = pointer_from_objref(obj)
    haskey(ss.seen, key) ? (return 0) : (ss.seen[key] = true)
    return Core.sizeof(obj) + (isdefined(obj, :mt) ? ss(obj.mt) : 0)
end

function (ss::SummarySize)(obj::Array)
    haskey(ss.seen, obj) ? (return 0) : (ss.seen[obj] = true)
    headersize = 4*sizeof(Int) + 8 + max(0, ndims(obj)-2)*sizeof(Int)
    size::Int = headersize
    datakey = unsafe_convert(Ptr{Cvoid}, obj)
    if !haskey(ss.seen, datakey)
        ss.seen[datakey] = true
        dsize = Core.sizeof(obj)
        if isbitsunion(eltype(obj))
            # add 1 union selector byte for each element
            dsize += length(obj)
        end
        size += dsize
        if !isempty(obj) && !Base.allocatedinline(eltype(obj))
            push!(ss.frontier_x, obj)
            push!(ss.frontier_i, 1)
        end
    end
    return size
end

function (ss::SummarySize)(obj::SimpleVector)
    key = pointer_from_objref(obj)
    haskey(ss.seen, key) ? (return 0) : (ss.seen[key] = true)
    size::Int = Core.sizeof(obj)
    if !isempty(obj)
        push!(ss.frontier_x, obj)
        push!(ss.frontier_i, 1)
    end
    return size
end

function (ss::SummarySize)(obj::Module)
    haskey(ss.seen, obj) ? (return 0) : (ss.seen[obj] = true)
    size::Int = Core.sizeof(obj)
    for binding in names(obj, all = true)
        if isdefined(obj, binding) && !isdeprecated(obj, binding)
            value = getfield(obj, binding)
            if !isa(value, Module) || parentmodule(value) === obj
                size += ss(value)::Int
                if isa(value, UnionAll)
                    value = unwrap_unionall(value)
                end
                if isa(value, DataType) && value.name.module === obj && value.name.name === binding
                    # charge a TypeName to its module (but not to the type)
                    size += ss(value.name)::Int
                end
            end
        end
    end
    return size
end

function (ss::SummarySize)(obj::Task)
    haskey(ss.seen, obj) ? (return 0) : (ss.seen[obj] = true)
    size::Int = Core.sizeof(obj)
    if isdefined(obj, :code)
        size += ss(obj.code)::Int
    end
    size += ss(obj.storage)::Int
    size += ss(obj.donenotify)::Int
    size += ss(obj.result)::Int
    # TODO: add stack size, and possibly traverse stack roots
    return size
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using SymbolicRegression
using DataDrivenDiffEq

"""
$(TYPEDEF)

Options for using SymbolicRegression.jl within the `solve` function.
Automatically creates [`Options`](https://astroautomata.com/SymbolicRegression.jl/stable/api/#Options) with the given specification.
Sorts the operators stored in `functions` into unary and binary operators on conversion.

# Fields
$(FIELDS)

"""
struct EQSearch <: AbstractSymbolicRegression
    "Operators used for symbolic regression"
    functions::AbstractVector{Function}
    "Additionally keyworded arguments passed to SymbolicRegression.Options"
    kwargs
end

function EQSearch(functions = Function[/, +, *, exp, cos]; kwargs...)
    return EQSearch(functions, kwargs)
end

function to_options(x::EQSearch)
    binops = Tuple([fi for fi in x.functions if is_binary(fi)])
    unaops = Tuple([fi for fi in x.functions if is_unary(fi)])
    fnames = fieldnames(SymbolicRegression.Options)
    ks = Dict(
        [k => v for (k,v) in x.kwargs if k ∈ fnames]
    )
    return SymbolicRegression.Options(;
        binary_operators = binops, unary_operators = unaops,
        ks...
    )
end

function DiffEqBase.solve(prob::AbstractDataDrivenProblem, alg::EQSearch;
    max_iter::Int = 10,
    weights = nothing,
    numprocs = nothing, procs = nothing,
    multithreading = false,
    runtests::Bool = true,
    eval_expression = false
    )

    opt = to_options(alg)

    Y = get_target(prob)

    # Inputs
    X̂, _, t, U = get_oop_args(prob)

    t = iszero(t) ? [] : t
    # Cat the inputs
    X = vcat([x for x in (X̂, U, permutedims(t)) if !isempty(x)]...)
    # Use the eqssearch of symbolic regression
    hof = SymbolicRegression.EquationSearch(X, Y, niterations = max_iter, weights = weights, options = opt,
            numprocs = numprocs, procs = procs, multithreading = multithreading,
            runtests = runtests)
    # Sort the paretofront
    doms = map(1:size(Y, 1)) do i
        calculateParetoFrontier(X, Y[i, :], hof[i], opt)
    end

    build_solution(prob, alg, doms; eval_expression = eval_expression)
end

function build_solution(prob::AbstractDataDrivenProblem, alg::EQSearch, doms; eval_expression = false)

    opt = to_options(alg)
    @variables x[1:size(prob.X, 1)] u[1:size(prob.U,1)] t
    x = Symbolics.scalarize(x)
    u = Symbolics.scalarize(u)
    x_ = [x;u;t]

    # Build a dict
    subs = Dict([SymbolicUtils.Sym{Number}(Symbol("x$(i)")) => x_[i] for i in 1:size(x_, 1)]...)
    # Create a variable
    eqs = vcat(map(x->node_to_symbolic(x[end].tree, opt), doms))
    eqs = map(x->substitute(x, subs), eqs)
    res_ = Basis(
        eqs, x, iv = t, controls = u, eval_expression = eval_expression
    )

    X = get_target(prob)
    Y = res_(get_oop_args(prob)...)

    # Build the metrics
    complexities = map(x->countNodes(x[end].tree), doms)
    complexity = sum(complexities)
    retcode = :success

    error = norm(X-Y, 2)
    k = free_parameters(res_)
    aic = AICC(k, X, Y)
    errors = zeros(eltype(X), size(Y, 1))
    aiccs = zeros(eltype(X), size(Y, 1))
    j = 1
    for i in 1:size(Y,1)
        errors[i] = norm(X[i,:].-Y[i,:],2)
        aiccs[i] = AICC(k, X[i:i, :], Y[i:i,:])
    end

    metrics = (
        Complexity = complexity,
        Error = error,
        AICC = aic,
        complexities = complexities,
        Errors = errors,
        AICCs = aiccs,
    )

    inputs = (
        Problem = prob,
        Algorithm = opt,
    )

    return DataDrivenSolution(
        res_, retcode, [], alg, doms, inputs, metrics
    )
end


text:
#Tutorial 2
#=
Daniel
=#
# ## is a break on the code

# control+c to stop the process
my_integer=5
my_float=5.
my_range=1:7
my_vector=[1, 3, 5]
my_string="hello"
my_character='o'
my_boolean=true

typeof(my_integer)
typeof(my_float)

typeof(
    my_character
)

# use \ to write simbols like alpha
α=5;
#show the command
@show α + my_integer
print("α = $α \n", "My integer is: $(my_integer)")

#println automatically saparated linea

#column vector
A = [1, 3, 5]
#row vector
A = [1 2 3]

ones(3)
zeros(5)
C=[1 3; 5 "test"]
C=[1 3; 5 9]

A = 1:9
# number with a distance o 0.5
B= 1:0.5:9
# secon element of break
B[3]

typeof(C)
eltype(C)
size(C)
length(C)
ndims(C)

#indexing
C[end]
C[1,2]
@show C
@show C[2,:] #show select row*

using LinearAlgebra

C*I
C^2
C.^2 #by element operation

sin.([1 2 3])

inv(C)
inv(C)*C # too costly
C\C # this notation is better

# to installa a package type ], add pkg name

using Random
rand(1)
rand(3)
randn(3)

Random.seed!(1234) #name of the seed
rand(1)
rand(1)

Random.seed!(1234)
rand(1)

#if/else statement
a=3
if a==3
    print("a is 3")
end

a=true
if a
    print("a is true")
end

a ? print("a is true") : nothing

if a
    print("a is true")
else
    print("a is no true")
end

b=true

if a | b
    print("a or b is true (or both)")
end

if a & b  # double & makes faster the process, just evaluate the firs statement.
    print("a and b are true")
end

#for loops
vals = ["a" "b" "c"]

for val in vals
    println(val)
end

#global
a = 3
for val in vals
    a = a + 1
    return a
end

for i in eachindex(vals)
    println(vals[i])
end

# index and value
for (index, value) in enumerate(vals)
    println("Element $index of vals is $value")
end

# Functions
function cobbdouglas(l, k, α, β = 1 - α; A=1)
    y=A*l^α*k^β
    return y
end

cobbdouglas(1, 3, 0.5)

@show output = cobbdouglas(1, 3, 0.5)

cobbdouglas.([1 2], 3, 0.5)

cobbdouglas(l, k, α) = A*l^α*k^(1-α)
cobbdouglas(1, 3, 0.5)


text:
# Dantzig-Wolfe Decomposition
# Edward J. Xu
# 2019.5.28
module DantzigWolfeOptim
export doDantzigWolfeOptim
########################################################################################################################
using JuMP
# using CPLEX
using Gurobi
using LinearAlgebra
include("FuncSub.jl")
include("FuncMas.jl")
include("FuncOptim.jl")
########################################################################################################################

function doDantzigWolfeOptim(
        vec_c, mat_a, vec_b, matIndexSub, numSub, numXInSub, m_mat_h, vecRowMatD,
        wheBranch, vecWhiSub, vecWhiVar, vecWhiBranch
        )
    timeStart = time()
    epsilon = 0.00001
    gurobi_env = Gurobi.Env()
    vecStrNameVar = getVecStrNameVar(numSub, numXInSub)
    println("#### 1/4,  Set vecModelSub #####################################################") ########################
    vecModelSub = setModelSub(
        numSub, matIndexSub, m_mat_h, vec_c, mat_a, vec_b, vecRowMatD, gurobi_env, numXInSub,
        wheBranch, vecWhiSub, vecWhiVar, vecWhiBranch
        )
    println("#### 2/4,  Set modelMas ########################################################") ########################
    (modMas, vecConsRef, vecConsConvex, vecLambda) = setModelMas(vecModelSub[1].mat_e, vec_b, numSub, gurobi_env)
    println("#### 3/4,  Begin Optim #########################################################") ########################
    (vecLambdaResult, extremePointForSub, extremePoints) = doOptim(
        vecModelSub, modMas, vecConsRef, vecConsConvex, vecLambda, epsilon
        )
    println("#### 4/4,  Print Result ########################################################") ########################
    # printResult(vecLambdaResult, extremePointForSub, extremePoints, vecStrNameVar, matIndexSub)
    println("################################################################################",
            "#### Elapsed time is $(time() - timeStart) seconds.")
    println("#### End #######################################################################") ########################
end

end


text:
include("parsemat.jl")
include("planeconnect.jl")
include("rectangulararray.jl")


text:
module Calmato

using Printf

using JuMP
using McCormick; xlogx
using RecipesBase
import EAGO
import ForwardDiff
import Ipopt
import GLPK
import MathOptInterface
import HTTP
import JSON
import TOML
import PeriodicTable
import ZipFile

MOI = MathOptInterface

export read_db, read_tdb, search_db, select, merge, merge!
export init_system
export equilib
export solidify

const R = 8.3144598 # gas constant
const P = 1.0 # pressure

function McCormick.xlogx(d::T) where T <: ForwardDiff.Dual
    return d * log(d)
end

include("essentials.jl")
include(joinpath("database", "database.jl"))
include(joinpath("database", "utils.jl"))
include(joinpath("database", "tdb.jl"))
include(joinpath("database", "toml.jl"))
include(joinpath("database", "operations.jl"))
include(joinpath("database", "tdbdb.jl"))
include("system.jl")
include("equilib.jl")
include("solidify.jl")

end


text:
module BetaTest

using Test
using ReactiveMP
using Distributions
using Random

import DomainSets

@testset "Generic Functional Distributions" begin

    @testset "ContinuousUnivariateLogPdf" begin

        @testset "Constructor" begin
            f  = (x) -> -x ^ 2
            d1 = ContinuousUnivariateLogPdf(f)
            d2 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), f)

            @test typeof(d1) === typeof(d2)
            @test d1 ≈ d2
        end

        @testset "pdf/logpdf" begin
            d1 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> -x ^ 2)

            f32_points1 = range(Float32(-10.0), Float32(10.0), length = 50)
            f64_points1 = range(-10.0, 10.0, length = 50)
            bf_points1  = range(BigFloat(-10.0), BigFloat(10.0), length = 50)
            points1 = vcat(f32_points1, f64_points1, bf_points1)

            @test all(map(p -> -p^2 == d1(p), points1))
            @test all(map(p -> -p^2 == logpdf(d1, p), points1))
            @test all(map(p -> exp(-p^2) == pdf(d1, p), points1))
            @test all(map(p -> -p^2 == d1([ p ]), points1))
            @test all(map(p -> -p^2 == logpdf(d1, [ p ]), points1))
            @test all(map(p -> exp(-p^2) == pdf(d1, [ p ]), points1))

            d2 = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -x ^ 4)

            f32_points2 = range(Float32(0.0), Float32(10.0), length = 50)
            f64_points2 = range(0.0, 10.0, length = 50)
            bf_points2  = range(BigFloat(0.0), BigFloat(10.0), length = 50)
            points2 = vcat(f32_points2, f64_points2, bf_points2)

            @test all(map(p -> -p^4 == d2(p), points2))
            @test all(map(p -> -p^4 == logpdf(d2, p), points2))
            @test all(map(p -> exp(-p^4) == pdf(d2, p), points2))
            @test all(map(p -> -p^4 == d2([ p ]), points2))
            @test all(map(p -> -p^4 == logpdf(d2, [ p ]), points2))
            @test all(map(p -> exp(-p^4) == pdf(d2, [ p ]), points2))

            @test_throws AssertionError d2(-1.0)
            @test_throws AssertionError logpdf(d2, -1.0)
            @test_throws AssertionError pdf(d2, -1.0)
            @test_throws AssertionError d2([ -1.0 ])
            @test_throws AssertionError logpdf(d2, [ -1.0 ])
            @test_throws AssertionError pdf(d2, [ -1.0 ])

            @test_throws AssertionError d2(Float32(-1.0))
            @test_throws AssertionError logpdf(d2, Float32(-1.0))
            @test_throws AssertionError pdf(d2, Float32(-1.0))
            @test_throws AssertionError d2([ Float32(-1.0) ])
            @test_throws AssertionError logpdf(d2, [ Float32(-1.0) ])
            @test_throws AssertionError pdf(d2, [ Float32(-1.0) ])

            @test_throws AssertionError d2(BigFloat(-1.0))
            @test_throws AssertionError logpdf(d2, BigFloat(-1.0))
            @test_throws AssertionError pdf(d2, BigFloat(-1.0))
            @test_throws AssertionError d2([ BigFloat(-1.0) ])
            @test_throws AssertionError logpdf(d2, [ BigFloat(-1.0) ])
            @test_throws AssertionError pdf(d2, [ BigFloat(-1.0) ])

            d3 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(Float32), (x) -> -x ^ 2)

            @test all(map(p -> -p^2 == d3(p), points1))
            @test all(map(p -> -p^2 == logpdf(d3, p), points1))
            @test all(map(p -> exp(-p^2) == pdf(d3, p), points1))
            @test all(map(p -> -p^2 == d3([ p ]), points1))
            @test all(map(p -> -p^2 == logpdf(d3, [ p ]), points1))
            @test all(map(p -> exp(-p^2) == pdf(d3, [ p ]), points1))

            d4 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(BigFloat), (x) -> -x ^ 2)

            @test all(map(p -> -p^2 == d4(p), points1))
            @test all(map(p -> -p^2 == logpdf(d4, p), points1))
            @test all(map(p -> exp(-p^2) == pdf(d4, p), points1))
            @test all(map(p -> -p^2 == d4([ p ]), points1))
            @test all(map(p -> -p^2 == logpdf(d4, [ p ]), points1))
            @test all(map(p -> exp(-p^2) == pdf(d4, [ p ]), points1))

            d5 = ContinuousUnivariateLogPdf(DomainSets.HalfLine{Float32}(), (x) -> -x ^ 2)

            @test all(map(p -> -p^2 == d5(p), points2))
            @test all(map(p -> -p^2 == logpdf(d5, p), points2))
            @test all(map(p -> exp(-p^2) == pdf(d5, p), points2))
            @test all(map(p -> -p^2 == d5([ p ]), points2))
            @test all(map(p -> -p^2 == logpdf(d5, [ p ]), points2))
            @test all(map(p -> exp(-p^2) == pdf(d5, [ p ]), points2))

            d6 = ContinuousUnivariateLogPdf(DomainSets.HalfLine{BigFloat}(), (x) -> -x ^ 2)

            @test all(map(p -> -p^2 == d6(p), points2))
            @test all(map(p -> -p^2 == logpdf(d6, p), points2))
            @test all(map(p -> exp(-p^2) == pdf(d6, p), points2))
            @test all(map(p -> -p^2 == d6([ p ]), points2))
            @test all(map(p -> -p^2 == logpdf(d6, [ p ]), points2))
            @test all(map(p -> exp(-p^2) == pdf(d6, [ p ]), points2))
        end

        @testset "support" begin
            d1 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> 1.0)
            @test minimum(support(d1)) === -Inf
            @test maximum(support(d1)) === Inf

            d2 = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> 1.0)
            @test minimum(support(d2)) === 0.0
            @test maximum(support(d2)) === Inf
        end

        @testset "vague" begin
            d = vague(ContinuousUnivariateLogPdf)

            @test typeof(d) <: ContinuousUnivariateLogPdf
            @test d ≈ ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> 1.0)
        end

        @testset "prod" begin
            d1 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> 2.0 * -x ^ 2)
            d2 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> 3.0 * -x ^ 2)

            pr1 = prod(ProdAnalytical(), d1, d2)
            pt1 = ContinuousUnivariateLogPdf(DomainSets.FullSpace(), (x) -> logpdf(d1, x) + logpdf(d2, x))

            @test isapprox(pr1, pt1, atol = 1e-12)

            d3 = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> 2.0 * -x ^ 2)
            d4 = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> 3.0 * -x ^ 2)

            pr2 = prod(ProdAnalytical(), d3, d4)
            pt2 = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> logpdf(d3, x) + logpdf(d4, x))

            @test isapprox(pr2, pt2, atol = 1e-12)

            @test !isapprox(pr1, pr2, atol = 1e-12)
        end

        @testset "convert" begin
            d = DomainSets.FullSpace()
            l = (x) -> 1.0

            c = convert(ContinuousUnivariateLogPdf, d, l)
            @test typeof(c) <: ContinuousUnivariateLogPdf
            @test isapprox(c, ContinuousUnivariateLogPdf(d, l), atol = 1e-12)

            c2 = convert(ContinuousUnivariateLogPdf, c)
            @test typeof(c2) <: ContinuousUnivariateLogPdf
            @test isapprox(c2, ContinuousUnivariateLogPdf(d, l), atol = 1e-12)

        end

    end

end

end


text:
# general.jl

export create_froms, count_trans_L, transTens_MLE, forecDist,
    deepcopyFields, timemod!, etr;

function create_froms(K::UInt16, L::UInt16)
  x = UInt16(1):K
  xx = collect(Iterators.product(Iterators.repeated(x,Int(L))...))
  (vec(xx), length(xx))
end
function create_froms(K::Int, L::Int)
  x = 1:K
  xx = collect(Iterators.product(Iterators.repeated(x,L)...))
  (vec(xx), length(xx))
end

"""
    count_trans_L(S, K, L)

Returns a tensor with counts of observed `L`th order transitions in `S`
  for a chain with `K` states. The first dimension is time now, the next
  dimension is lag 1, then lag 2, and so on.

### Example
```julia
  S = [1,2,1,1,1,2,4,1,3,1,1,3,1,2,1,1,2]
  count_trans_L(S, 4, 2)
```
"""
function count_trans_L(S::Vector{Int}, K::Int, L::Int)
    @assert minimum(S) >= 1 && maximum(S) <= K
    N = zeros(Int, fill(K,L+1)...)

    for tt in (L+1):length(S)
      Srev_now = S[range(tt, step=-1, length=(L+1))]
      N[(Srev_now)...] += 1
    end

    N
end

"""
    transTens_MLE(S, K, L)

Returns a tensor with estimated transition probabilities for `L`th order transitions in `S`
  for a chain with `K` states. The first dimension is time now, the next
  dimension is lag 1, then lag 2, and so on.

### Example
```julia
  S = [1,2,1,1,1,2,4,1,3,1,1,3,1,2,1,1,2]
  N = count_trans_L(S, 4, 2)
  transTens_MLE(S, 4, 2)
  transTens_MLD(N)
```
"""
function transTens_MLE(S::Vector{Int}, K::Int, L::Int)
    @assert minimum(S) >= 1 && maximum(S) <= K
    N = count_trans_L(S, K, L)
    N ./ sum(N, dims=1)
end
function transTens_MLE(N::Union{Array{Float64}, Array{Int64}})
    N ./ sum(N, dims=1)
end

"""
    forecDist(Slag::Vector{Int}, λ::Vector{Float64}, Q::Matrix{Float64})

Computes the forecast distribution given lagged values and parameters.

### Example
```julia
    lam = [0.6, 0.3, 0.1]
    Q = reshape([0.7, 0.3, 0.25, 0.75], (2,2))
    Slagrev = [1,1,2]
    forecDist(Slagrev, lam, Q)
```
"""
function forecDist(Slagrev::Vector{Int}, λ::Vector{Float64}, Q::Matrix{Float64})
    K = size(Q)[1]
    L = size(λ)[1]
    @assert size(Slagrev)[1] == L
    w = zeros(Float64, K)
    for k in 1:K
        for ℓ in 1:L
            w[k] += λ[ℓ] * Q[k, Slagrev[ℓ]]
        end
    end
    w
end

## from Arthur Lui
function deepcopyFields(state::T, fields::Vector{Symbol}) where T
  substate = Dict{Symbol, Any}()

  for field in fields
    substate[field] = deepcopy(getfield(state, field))
  end

  return substate
end

function postSimsInit(monitor::Union{Vector{Symbol}, Nothing}, n_keep::Int, init_state::Union{ParamsMTD, ParamsMTDg, ParamsMMTD})

    if typeof(monitor) == Vector{Symbol}
        state = deepcopyFields(init_state, monitor)
    else
        state = Dict{Symbol, Any}()
    end

    state[:llik] = 0.0

    sims = [ deepcopy(state) for i = 1:n_keep ]

    return sims
end

## MCMC timing for benchmarks
function timemod!(n::Int64, model::Union{ModMTD, ModMTDg, ModMMTD}, niter::Int, outfilename::String)
    outfile = open(outfilename, "a+")
    write(outfile, "timing for $(niter) iterations each:\n")
    for i in 1:n
        tinfo = @timed mcmc!(model, niter, save=false, report_filename=outfilename)
        write(outfile, "trial $(i), elapsed: $(tinfo[2]) seconds, allocation: $(tinfo[3]/1.0e6) Megabytes\n")
    end
    close(outfile)
end

## estimate time remaining
function etr(timestart::DateTime, n_keep::Int, thin::Int, outfilename::String)
    timeendburn = now()
    durperiter = (timeendburn - timestart).value / 1.0e5 # in milliseconds
    milsecremaining = durperiter * (n_keep * thin)
    estimatedfinish = now() + Dates.Millisecond(Int64(round(milsecremaining)))
    report_file = open(outfilename, "a+")
    write(report_file, "Completed burn-in at $(durperiter/1.0e3*1000.0) seconds per 1000 iterations \n
      $(durperiter/1.0e3/60.0*1000.0) minutes per 1000 iterations \n
      $(durperiter/1.0e3/60.0/60.0*1000.0) hours per 1000 iterations \n
      estimated completion time $(estimatedfinish)\n\n")
    close(report_file)
end

max_stars_repo_name: mheiner/MTD

max_issues_repo_name: mheiner/MTD

max_forks_repo_name: mheiner/MTD


text:
module Day6

using Debugger

inp = Dict{String,String}()
for line in readlines("day6.e")
    a, b = split(line, ')')
    inp[b] = a # b is in orbit around a
end

function h(input, x)
    return (x in keys(input) ? union(h(input, input[x]), [x]) : [])
end

function solve(input)
    counts = Dict(k => h(input, k) for k in keys(input))

    sum(x -> length(x), values(counts)) # a

    a = counts["YOU"]
    b = counts["SAN"]
    while a[1] == b[1]
        @info a b
        deleteat!(a, 1)
        deleteat!(b, 1)
    end

    length(a) + length(b) - 2
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
#=
    For this particular positions
        r = [1 2 5; 3 5 -1/2]
    the expression below are exact solutions.

    Check sagemath notebook for the proof.
=#
xx = [1im -(3*√(173)/5177717)*(-250*√(173)*im - 26661)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-250*√(173)*im - 26661)*exp(0.5*im*√(173)) 1im]
xy = [0 -(3*√(173)/5177717)*(144*√(173)*im + 3864)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(144*√(173)*im + 3864)*exp(0.5*im*√(173)) 0]
xz = [0 -(3*√(173)/5177717)*(-264*√(173)*im - 7084)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-264*√(173)*im - 7084)*exp(0.5*im*√(173)) 0]

yx = [0 -(3*√(173)/5177717)*(144*√(173)*im + 3864)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(144*√(173)*im + 3864)*exp(0.5*im*√(173)) 0]
yy = [1im -(3*√(173)/5177717)*(-130*√(173)*im - 23441)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-130*√(173)*im - 23441)*exp(0.5*im*√(173)) 1im]
yz = [0 -(3*√(173)/5177717)*(-396*√(173)*im - 10626)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-396*√(173)*im - 10626)*exp(0.5*im*√(173)) 0]

zx = [0 -(3*√(173)/5177717)*(-264*√(173)*im - 7084)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-264*√(173)*im - 7084)*exp(0.5*im*√(173)) 0]
zy = [0 -(3*√(173)/5177717)*(-396*√(173)*im - 10626)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(-396*√(173)*im - 10626)*exp(0.5*im*√(173)) 0]
zz = [1im -(3*√(173)/5177717)*(380*√(173)*im - 9756)*exp(0.5*im*√(173)); -(3*√(173)/5177717)*(380*√(173)*im - 9756)*exp(0.5*im*√(173)) 1im]

K_exact = [xx xy xz; yx yy yz; zx zy zz]

#=
        Creating individual elements
=#
let
    using LinearAlgebra
    N = 2
    r = [1 2 5; 3 5 -1/2]
    δ(x,y) = Int(==(x,y))

    P(x) =  1 - 1/x + 1/x^2
    Q(x) = -1 + 3/x - 3/x^2

    α_range = β_range = [-1, +1, 0]

    for (α_idx, α) ∈ enumerate(α_range)
        for (β_idx, β) ∈ enumerate(β_range)
            for j ∈ 1:N
                for n ∈ 1:N
                    if j≠n
                        r_jn_vec = r[j,:] - r[n,:]
                        r_jn = norm(r_jn_vec)

                        term1 = im.*δ(j,n)*δ(α,β)
                        term2 = (1 - δ(j,n))*(3/2)*(exp(im*r_jn)/r_jn)

                        term3 = P(im*r_jn)*δ(α,β)
                        term4 = Q(im*r_jn)*r_jn_vec[α_idx]*r_jn_vec[β_idx]/r_jn^2
                    else
                        term1 = im.*δ(j,n)*δ(α,β)
                        term2 = term3 = term4 = 0
                    end
                    temp =  term1 + term2*(term3 + term4 )
                    println("α:$(α) ... β:$(β) ... j:$(j) ... n:$(n) ... temp:$(temp)")
                end
            end
            println("---------")
        end
    end
end

#=
        Creating sub matrices but not storing them
=#
let
    using LinearAlgebra
    N = 2
    r = [1 2 5; 3 5 -1/2]
    δ(x,y) = Int(==(x,y))
    Xt, Yt, Zt = r[:,1], r[:,2] ,r[:,3]

    Xjn = Xt*ones(1,N) - ones(N,1)*Xt'
    Yjn = Yt*ones(1,N) - ones(N,1)*Yt'
    Zjn = Zt*ones(1,N) - ones(N,1)*Zt'

    array_XYZ_jn = [Xjn, Yjn, Zjn]
    Rjn = sqrt.(Xjn.^2 + Zjn.^2 + Yjn.^2)

    P(x) =  1 - 1/x + 1/x^2
    Q(x) = -1 + 3/x - 3/x^2

    α_range = β_range = [-1, +1, 0]

    for (α_idx, α) ∈ enumerate(α_range)
        for (β_idx, β) ∈ enumerate(β_range)
            term1 = im*I(N).*δ(α,β)
            term2 = (3/2)*exp.(im*Rjn)./Rjn
            term3 = P.(im*Rjn).*δ(α,β)
            term4 = (array_XYZ_jn[α_idx].*array_XYZ_jn[β_idx]).*Q.(im*Rjn)./Rjn.^2

            term2[findall(isnan.(term2))] .= 0
            term2[findall(isinf.(term2))] .= 0

            term3[findall(isnan.(term3))] .= 0
            term3[findall(isinf.(term3))] .= 0

            term4[findall(isnan.(term4))] .= 0
            term4[findall(isinf.(term4))] .= 0

            K = term1 + term2.*(term3 .+ term4)
            println("------")
            display(K)
        end
    end
end

#=
        Creating matrices generically
=#
let
    using LinearAlgebra
    N = 2
    r = [1 2 5; 3 5 -1/2]
    δ(x,y) = Int(==(x,y))

    Xt, Yt, Zt = r[:,1], r[:,2] ,r[:,3]

    Xjn = Xt*ones(1,N) - ones(N,1)*Xt'
    Yjn = Yt*ones(1,N) - ones(N,1)*Yt'
    Zjn = Zt*ones(1,N) - ones(N,1)*Zt'

    array_XYZ_jn = [Xjn, Yjn, Zjn]
    Rjn = sqrt.(Xjn.^2 + Zjn.^2 + Yjn.^2)

    P(x) =  1 - 1/x + 1/x^2
    Q(x) = -1 + 3/x - 3/x^2

    α_range = β_range = [-1, +1, 0]

    term2 = (3/2)*exp.(im*Rjn)./Rjn
    term2[findall(isnan.(term2))] .= 0
    term2[findall(isinf.(term2))] .= 0

    P_Rjn = P.(im*Rjn)
    P_Rjn[findall(isnan.(P_Rjn))] .= 0

    Q_Rjn_over_Rjn_squared = Q.(im*Rjn)./Rjn.^2
    Q_Rjn_over_Rjn_squared[findall(isnan.(Q_Rjn_over_Rjn_squared))] .= 0

    A = []
    for (α_idx, α) ∈ enumerate(α_range)
        B = []
        for (β_idx, β) ∈ enumerate(β_range)
            term1 = im*I(N).*δ(α,β)
            # term2 = (3/2)*exp.(im*Rjn)./Rjn  ## Defined outside
            term3 = P_Rjn.*δ(α,β)
            term4 = Q_Rjn_over_Rjn_squared.*(array_XYZ_jn[α_idx].*array_XYZ_jn[β_idx])

            K = term1 + term2.*(term3 .+ term4)
            push!(B, K)
        end
        push!(A, vcat(B[1:length(β_range)]...)  )
    end

    K = hcat(A[1:length(α_range)]...)
    all(K .≈ K_exact)
end


text:
# Autogenerated wrapper script for MAGEMin_jll for aarch64-apple-darwin-libgfortran5
export MAGEMin, libMAGEMin

using MPICH_jll
using NLopt_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MAGEMin")
JLLWrappers.@declare_library_product(libMAGEMin, "@rpath/libMAGEMin.dylib")
JLLWrappers.@declare_executable_product(MAGEMin)
function __init__()
    JLLWrappers.@generate_init_header(MPICH_jll, NLopt_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libMAGEMin,
        "lib/libMAGEMin.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        MAGEMin,
        "bin/MAGEMin",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#scaling for 0.0 should be 1.0
@test GeneralizedMetropolisHastings.tuneerf(0.0,3.0) == 1.0
@test GeneralizedMetropolisHastings.tunelogistic(0.0,7.0) == 1.0

#smalest scale factor should be around sqrt(2)/2
@test_approx_eq_eps GeneralizedMetropolisHastings.tuneerf(-0.5,3.0) sqrt(2)/2 10-4
@test_approx_eq_eps GeneralizedMetropolisHastings.tunelogistic(-0.5,7.0) sqrt(2)/2 10-4

#largest scale factor should be around sqrt(2)
@test_approx_eq_eps GeneralizedMetropolisHastings.tuneerf(0.5,3.0) sqrt(2) 10-4
@test_approx_eq_eps GeneralizedMetropolisHastings.tunelogistic(0.5,7.0) sqrt(2) 10-4

t1 = tuner(:monitor,10)
t2 = tuner(:scale,4,0.7,:logistic)
t3 = tuner(:scale,4,0.5,:erf)

@test period(t1) == 10
@test period(t2) == 4
@test period(t3) == 4

@test verbose(t1)
@test verbose(t2)
@test verbose(t3)

s1 = tunerstate(t1,100,Float64)
s2 = tunerstate(t2,32,Float64)
s3 = tunerstate(t3,33,Float32)

@test all(isnan(rate(s1))) && numtunesteps(s1) == 10
@test all(isnan(rate(s2))) && numtunesteps(s2) == 8
@test all(isnan(rate(s3))) && numtunesteps(s3) == 9

s1.proposed[s1.index] = 10 ; s1.accepted[s1.index] = 8
s2.proposed[s2.index] = 40 ; s2.accepted[s2.index] = 25
s3.proposed[s3.index] = 20 ; s3.accepted[s3.index] = 11

@test rate(s1)[s1.index] == 8/10 && total(s1) == 10
@test rate(s2)[s2.index] == 25/40 && total(s2) == 40
@test rate(s3)[s3.index] == 11/20 && total(s3) == 20

@test tune(t1,s1) == ()
@test tune(t2,s2) < 1.0
@test tune(t3,s3) > 1.0
@test_approx_eq s2.scalefactor[s2.index] GeneralizedMetropolisHastings.tunelogistic(accepted(s2)[s2.index]/proposed(s2)[s2.index]-t2.targetrate,7.0)
@test_approx_eq s3.scalefactor[s3.index] GeneralizedMetropolisHastings.tuneerf(accepted(s3)[s3.index]/proposed(s3)[s3.index]-t3.targetrate,3.0)

nextindex!(s1)
nextindex!(s2)
nextindex!(s3)

@test index(s1) == 2 && index(s2) == 2 && index(s3) == 2

accepted!(s1,IndicatorStationary{Float64}(zeros(6),[1,1,1,2,3,4]))
accepted!(s2,IndicatorStationary{Float64}(zeros(5),[1,2,3,4,5]))
accepted!(s3,IndicatorStationary{Float64}(zeros(11),ones(Int,11)))

@test proposed(s1)[s1.index] == 5 && accepted(s1)[s1.index] == 3 && total(s1) == 15
@test proposed(s2)[s2.index] == 4 && accepted(s2)[s2.index] == 4 && total(s2) == 44
@test proposed(s3)[s3.index] == 10 && accepted(s3)[s3.index] == 0 && total(s3) == 30

s2.index = 9
println("A warning should be printed below this line")
accepted!(s2,IndicatorStationary{Float64}(zeros(5),[1,2,3,4,5]))

###test the show functions
println()
println("====================")
println("Test show() function")
println("====================")
show(t1)
show(t2)
show(t3)
show(s1)
show(s2)
show(s3)
showstep(t1,s1)
showstep(t2,s2)
showstep(t3,s3)
println("====================")
println("End  show() function")
println("====================")
println()


text:
#!/usr/bin/julia

function err()
   println("Usage: please input a number")
end

function even_odd(n)

    if (n % 2 == 0 )
      return "Even"
    else
      return "Odd"
    end
end

try
    println(even_odd(parse(Int, ARGS[1])))
catch e
    err()
end


text:
__precompile__()

module FLAC

using FileIO

# This `using` is literally only just so that `Ogg.__init__()` gets run.  This
# ensures that `libogg` is loaded into the Julia namcespace, which is necessary
# for `libFLAC` to load properly.  This will not be necessary in the future,
# once https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/194 is solved.
using Ogg

export StreamMetaData,
       InfoMetaData,
       VorbisCommentMetaData,
       PaddingMetaData,
       ApplicationMetaData,
       SeekTableMetaData,
       CueSheetMetaData,
       StreamDecoderPtr,
       StreamEncoderPtr,
       initfile!,
       FLACDecoder,
       seek,
       read,
       size,
       length

import Base: read, seek, size, length

const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("FLAC not properly installed. Please run Pkg.build(\"FLAC\")")
end

include("metadata.jl")
include("format.jl")
include("decoder.jl")
include("encoder.jl")

function __init__()
    init_decoder_cfunctions()
    init_encoder_cfunctions()
end

end # module


text:
function text_record(data, text, hdr=Vector{String}())
	# Create a text record to send to pstext. DATA is the Mx2 coordinates array.
	# TEXT is a string or a cell array

	if (isa(data, Vector))  data = data[:,:]  end 	# Needs to be 2D
	if (!isa(data, Array{Float64}))  data = Float64.(data)  end

	if (isa(text, String))
		T = GMTdataset(data, [text], "", Vector{String}(), "", "", 0)
	elseif (isa(text, Array{String}))
		if (text[1][1] == '>')			# Alternative (but risky) way of setting the header content
			T = GMTdataset(data, text[2:end], text[1], Vector{String}(), "", "", 0)
		else
			T = GMTdataset(data, text, (isempty(hdr) ? "" : hdr), Vector{String}(), "", "", 0)
		end
	elseif (isa(text, Array{Array}) || isa(text, Array{Vector{String}}))
		nl_t = length(text);	nl_d = length(data)
		(nl_d > 0 && nl_d != nl_t) && error("Number of data points is not equal to number of text strings.")
		T = Vector{GMTdataset}(undef,nl_t)
		for k = 1:nl_t
			T[k] = GMTdataset((nl_d == 0 ? data : data[k]), text[k], (isempty(hdr) ? "" : hdr[k]), Vector{String}(), "", "", 0)
		end
	else
		error("Wrong type ($(typeof(text))) for the 'text' argin")
	end
	return T
end
text_record(text) = text_record(Array{Float64,2}(undef,0,0), text)
text_record(text::Array{String}, hdr::String) = text_record(Array{Float64,2}(undef,0,0), text, hdr)

# ---------------------------------------------------------------------------------------------------
"""
    D = mat2ds(mat [,txt]; x=nothing, text=nothing, multi=false, kwargs...)

Take a 2D `mat` array and convert it into a GMTdataset. `x` is an optional coordinates vector (must have the
same number of elements as rows in `mat`). Use `x=:ny` to generate a coords array 1:n_rows of `mat`.
  - `txt`   Return a Text record which is a Dataset with data = Mx2 and text in third column. The ``text``
     can be an array with same size as ``mat``rows or a string (will be reapeated n_rows times.)
  - `x`   An optional vector with the xx coordinates
  - `hdr` optional String vector with either one or n_rows multisegment headers.
  - `color` optional array of strings with color names/values. Its length can be smaller than n_rows, case in
     which colors will be cycled.
  - `linethick`, or `lt` for selecting different line thicknesses. Work alike `color`, but should be
     a vector of numbers, or just a single number that is then appl	ied to all lines.
  - `fill`  Optional string array with color names or array of "patterns"
  - `ls` | `linestyle`  Line style. A string or an array of strings with ``length = size(mat,1)`` with line styles.
  - `lt` | `linethick`  Line thickness.
  - `multi` When number of columns in `mat` > 2, or == 2 and x != nothing, make an multisegment Dataset with
     first column and 2, first and 3, etc. Convenient when want to plot a matrix where each column is a line.
  - `datatype` Keep the original data type of `mat`. Default, converts to Float64
  - `proj` or `proj4`  A proj4 string for dataset SRS
  - `wkt`  A WKT SRS
"""
function mat2ds(mat, txt=Vector{String}(); hdr=Vector{String}(), geom=0, kwargs...)
	d = KW(kwargs)

	(!isempty(txt)) && return text_record(mat, txt,  hdr)
	((text = find_in_dict(d, [:text])[1]) !== nothing) && return text_record(mat, text, hdr)

	val = find_in_dict(d, [:multi :multicol])[1]
	multi = (val === nothing) ? false : ((val) ? true : false)	# Like this it will error if val is not Bool

	if ((x = find_in_dict(d, [:x])[1]) !== nothing)
		n_ds = (multi) ? size(mat, 2) : 1
		xx = (x == :ny || x == "ny") ? collect(1.0:size(mat, 1)) : x
		(length(xx) != size(mat, 1)) && error("Number of X coordinates and MAT number of rows are not equal")
	else
		n_ds = (ndims(mat) == 3) ? size(mat,3) : ((multi) ? size(mat, 2) - 1 : 1)
		xx = Vector{Float64}()
	end

	if (!isempty(hdr) && isa(hdr, String))	# Accept one only but expand to n_ds with the remaining as blanks
		bak = hdr;		hdr = Base.fill("", n_ds);	hdr[1] = bak
	elseif (!isempty(hdr) && length(hdr) != n_ds)
		error("The header vector can only have length = 1 or same number of MAT Y columns")
	end

	if ((color = find_in_dict(d, [:color])[1]) !== nothing)
		_color::Array{String} = isa(color, Array{String}) ? color : ["#0072BD", "#D95319", "#EDB120", "#7E2F8E", "#77AC30", "#4DBEEE", "#A2142F"]
	end
	_fill = helper_ds_fill(d)

	# ---  Here we deal with line colors and line thickness. If not provided we override the GMR defaultb -Wthin ---
	val = find_in_dict(d, [:lt :linethick :linethickness])[1]
	_lt = (val === nothing) ? [0.5] : val
	_lts = Vector{String}(undef, n_ds)
	n_thick = length(_lt)
	[_lts[k] = " -W" * string(_lt[((k % n_thick) != 0) ? k % n_thick : n_thick])  for k = 1:n_ds]

	if (color !== nothing)
		n_colors = length(_color)
		if (isempty(hdr))
			hdr = Vector{String}(undef, n_ds)
			[hdr[k]  = _lts[k] * string(",", _color[((k % n_colors) != 0) ? k % n_colors : n_colors])  for k = 1:n_ds]
		else
			[hdr[k] *= _lts[k] * string(",", _color[((k % n_colors) != 0) ? k % n_colors : n_colors])  for k = 1:n_ds]
		end
	else						# Here we just overriding the GMT -W default that is too thin.
		if (isempty(hdr))
			hdr = Vector{String}(undef, n_ds)
			[hdr[k]  = _lts[k] for k = 1:n_ds]
		else
			[hdr[k] *= _lts[k] for k = 1:n_ds]
		end
	end
	# ----------------------------------------

	if ((ls = find_in_dict(d, [:ls :linestyle])[1]) !== nothing && ls != "")
		if (isa(ls, AbstractString) || isa(ls, Symbol))
			[hdr[k] = string(hdr[k], ',', ls) for k = 1:n_ds]
		else
			[hdr[k] = string(hdr[k], ',', ls[k]) for k = 1:n_ds]
		end
	end

	if (!isempty(_fill))				# Paint the polygons (in case of)
		n_colors = length(_fill)
		if (isempty(hdr))
			hdr = Array{String,1}(undef, n_ds)
			[hdr[k]  = " -G" * _fill[((k % n_colors) != 0) ? k % n_colors : n_colors]  for k = 1:n_ds]
		else
			[hdr[k] *= " -G" * _fill[((k % n_colors) != 0) ? k % n_colors : n_colors]  for k = 1:n_ds]
		end
	end

	prj = ((proj = find_in_dict(d, [:proj :proj4])[1]) !== nothing) ? proj : ""
	(prj == "geo" || prj == "geog") && (prj = prj4WGS84)
	(prj != "" && !startswith(prj, "+proj=")) && (prj = "+proj=" * prj)
	wkt = ((wk = find_in_dict(d, [:wkt])[1]) !== nothing) ? wk : ""

	D = Vector{GMTdataset}(undef, n_ds)

	# By default convert to Doubles, except if instructed to NOT to do it.
	(find_in_dict(d, [:datatype])[1] === nothing) && (eltype(mat) != Float64) && (mat = Float64.(mat))
	geom = (geom == 0 && (2 <= length(mat) <= 3)) ? Gdal.wkbPoint : (geom == 0 ? Gdal.wkbUnknown : UInt32(geom))	# Guess geom
	(multi && geom == 0 && size(mat,1) == 1) && (geom = Gdal.wkbPoint)	# One row with many columns and MULTI => Points
	if (isempty(xx))				# No coordinates transmitted
		if (ndims(mat) == 3)
			[D[k] = GMTdataset(view(mat,:,:,k), String[], (isempty(hdr) ? "" : hdr[k]), String[], prj, wkt, geom) for k = 1:n_ds]
		elseif (!multi)
			D[1] = GMTdataset(mat, String[], (isempty(hdr) ? "" : hdr[1]), String[], prj, wkt, geom)
		else
			[D[k] = GMTdataset(mat[:,[1,k+1]], String[], (isempty(hdr) ? "" : hdr[k]), String[], prj, wkt, geom) for k = 1:n_ds]
		end
	else
		if (!multi)
			D[1] = GMTdataset(hcat(xx,mat), String[], (isempty(hdr) ? "" : hdr[1]), String[], prj, wkt, geom)
		else
			[D[k] = GMTdataset(hcat(xx,mat[:,k]), String[], (isempty(hdr) ? "" : hdr[k]), String[], prj, wkt, geom) for k = 1:n_ds]
		end
	end
	return D
end

# ---------------------------------------------------------------------------------------------------
function ds2ds(D::GMTdataset; kwargs...)::Vector{<:GMTdataset}
	# Take one DS and split it in an array of DS, one for each row and optionally add -G,fill>
	# So far only for internal use but may grow in function of needs
	d = KW(kwargs)

	#multi = "r"		# Default split by rows
	#if ((val = find_in_dict(d, [:multi])[1]) !== nothing)  multi = "c"  end		# Then by columns
	_fill = helper_ds_fill(d)

	if ((val = find_in_dict(d, [:color_wrap])[1]) !== nothing)	# color_wrap is a kind of private option for bar-stack
		n_colors = Int(val)
	end

	n_ds = size(D.data, 1)
	if (!isempty(_fill))				# Paint the polygons (in case of)
		hdr = Vector{String}(undef, n_ds)
		[hdr[k] = " -G" * _fill[((k % n_colors) != 0) ? k % n_colors : n_colors]  for k = 1:n_ds]
		if (D.header != "")  hdr[1] = D.header * hdr[1]  end	# Copy eventual contents of first header
	end

	Dm = Vector{GMTdataset}(undef, n_ds)
	for k = 1:n_ds
		Dm[k] = GMTdataset(D.data[k:k, :], String[], (isempty(_fill) ? "" : hdr[k]), String[], "", "", 0)
	end
	Dm[1].comment = D.comment;	Dm[1].proj4 = D.proj4;	Dm[1].wkt = D.wkt
	(size(D.text) == n_ds) && [Dm.text[k] = D.text[k] for k = 1:n_ds]
	Dm
end

# ------------------------------
function helper_ds_fill(d::Dict)
	# Shared by ds2ds & mat2ds
	if ((fill_val = find_in_dict(d, [:fill :fillcolor])[1]) !== nothing)
		_fill::Array{String} = (isa(fill_val, Array{String}) && !isempty(fill_val)) ? fill_val :
		                       ["#0072BD", "#D95319", "#EDB120", "#7E2F8E", "#77AC30", "#4DBEEE", "#A2142F", "0/255/0"]
		n_colors = length(_fill)
		if ((alpha_val = find_in_dict(d, [:fillalpha])[1]) !== nothing)
			if (eltype(alpha_val) <: AbstractFloat && maximum(alpha_val) <= 1)  alpha_val = collect(alpha_val) .* 100  end
			_alpha = Vector{String}(undef, n_colors)
			na = min(length(alpha_val), n_colors)
			[_alpha[k] = join(string('@',alpha_val[k])) for k = 1:na]
			(na < n_colors) && [_alpha[k] = "" for k = na+1:n_colors]
			[_fill[k] *= _alpha[k] for k = 1:n_colors]		# And finaly apply the transparency
		end
	else
		_fill = Vector{String}()
	end
	return _fill
end

# ---------------------------------------------------------------------------------------------------
"""
    I = mat2img(mat::Array{<:Unsigned}; x=[], y=[], hdr=nothing, proj4="", wkt="", cmap=nothing, kw...)

Take a 2D 'mat' array and a HDR 1x9 [xmin xmax ymin ymax zmin zmax reg xinc yinc] header descriptor
and return a GMTimage type.
Alternatively to HDR, provide a pair of vectors, x & y, with the X and Y coordinates.
Optionaly, the HDR arg may be ommited and it will computed from 'mat' alone, but then x=1:ncol, y=1:nrow
When 'mat' is a 3D UInt16 array we automatically compute a UInt8 RGB image. In that case 'cmap' is ignored.
But if no conversion is wanted use option `noconv=true`

    I = mat2img(mat::Array{UInt16}; x=[], y=[], hdr=nothing, proj4::String="", wkt::String="", kw...)

Take a `mat` array of UInt16 and scale it down to UInt8. Input can be 2D or 3D.
If the kw variable `stretch` is used, we stretch the intervals in `stretch` to [0 255].
Use this option to stretch the image histogram.
If `stretch` is a scalar, scale the values > `stretch` to [0 255]
  - stretch = [v1 v2] scales all values >= v1 && <= v2 to [0 255]
  - stretch = [v1 v2 v3 v4 v5 v6] scales firts band >= v1 && <= v2 to [0 255], second >= v3 && <= v4, same for third
  - stretch = :auto | "auto" | true | 1 will do an automatic stretching from values obtained from histogram thresholds

The `kw...` kwargs search for [:layout :mem_layout], [:names] and [:metadata]
"""
function mat2img(mat::AbstractArray{<:Unsigned}, dumb::Int=0; x=Vector{Float64}(), y=Vector{Float64}(), v=Vector{Float64}(), hdr=nothing, proj4::String="", wkt::String="", cmap=nothing, kw...)
	# Take a 2D array of uint8 and turn it into a GMTimage.
	# Note: if HDR is empty we guess the registration from the sizes of MAT & X,Y
	color_interp = "";		n_colors = 0;
	if (cmap !== nothing)
		have_alpha = !all(cmap.alpha .== 0.0)
		nc = have_alpha ? 4 : 3
		colormap = zeros(Clong, 256 * nc)
		n_colors = 256;			# Because for GDAL we always send 256 even if they are not all filled
		@inbounds for n = 1:3	# Write 'colormap' row-wise
			@inbounds for m = 1:size(cmap.colormap, 1)
				colormap[m + (n-1)*n_colors] = round(Int32, cmap.colormap[m,n] * 255);
			end
		end
		if (have_alpha)						# Have alpha color(s)
			[colormap[m + 3*n_colors] = round(Int32, cmap.colormap[m,4] * 255) for m = 1:size(cmap.colormap, 1)]
			n_colors *= 1000				# Flag that we have alpha colors in an indexed image
		end
	else
		(size(mat,3) == 1) && (color_interp = "Gray")
		colormap = zeros(Clong,3)			# Because we need an array
	end

	nx = size(mat, 2);		ny = size(mat, 1);
	reg = (hdr !== nothing) ? Int(hdr[7]) : (nx == length(x) && ny == length(y)) ? 0 : 1
	x, y, hdr, x_inc, y_inc = grdimg_hdr_xy(mat, reg, hdr, x, y)

	mem_layout = (size(mat,3) == 1) ? "TCBa" : "TCBa"		# Just to have something. Likely wrong for 3D
	d = KW(kw)
	((val = find_in_dict(d, [:layout :mem_layout])[1]) !== nothing) && (mem_layout = string(val))
	_names = ((val = find_in_dict(d, [:names])[1]) !== nothing) ? val : String[]
	_meta  = ((val = find_in_dict(d, [:metadata])[1]) !== nothing) ? val : String[]

	I = GMTimage(proj4, wkt, 0, hdr[:], [x_inc, y_inc], reg, NaN, color_interp, _meta, _names,
	             x,y,v,mat, colormap, n_colors, Array{UInt8,2}(undef,1,1), mem_layout, 0)
end

# ---------------------------------------------------------------------------------------------------
function mat2img(mat::AbstractArray{UInt16}; x=Vector{Float64}(), y=Vector{Float64}(), v=Vector{Float64}(), hdr=nothing, proj4::String="", wkt::String="", img8=Matrix{UInt8}(undef,0,0), kw...)
	# Take an array of UInt16 and scale it down to UInt8. Input can be 2D or 3D.
	# If the kw variable 'stretch' is used, we stretch the intervals in 'stretch' to [0 255].
	# Use this option to stretch the image histogram.
	# If 'stretch' is a scalar, scale the values > 'stretch' to [0 255]
	# stretch = [v1 v2] scales all values >= v1 && <= v2 to [0 255]
	# stretch = [v1 v2 v3 v4 v5 v6] scales firts band >= v1 && <= v2 to [0 255], second >= v3 && <= v4, same for third
	# If the img8 argument is used, it should contain a pre-allocated UInt8 array with the exact same size as MAT.
	# The 'scale_only' kw option makes it return just the scaled array an no GMTimage creation (useful when img8 is a view).
	# Use the keyword NOCONV to return GMTimage UInt16 type. I.e., no conversion to UInt8

	d = KW(kw)
	if ((val = find_in_dict(d, [:noconv])[1]) !== nothing)		# No conversion to UInt8 is wished
		return mat2img(mat, 1; x=x, y=y, v=v, hdr=hdr, proj4=proj4, wkt=wkt, d...)
	end
	img = isempty(img8) ? Array{UInt8}(undef, size(mat)) : img8
	(size(img) != size(mat)) && error("Incoming matrix and image holder have different sizes")
	if ((vals = find_in_dict(d, [:histo_bounds :stretch], false)[1]) !== nothing)
		nz = 1
		isa(mat, Array{UInt16,3}) ? (ny, nx, nz) = size(mat) : (ny, nx) = size(mat)

		(vals == "auto" || vals == :auto || (isa(vals, Bool) && vals) || (isa(vals, Real) && vals == 1)) &&
			(vals = [find_histo_limits(mat)...])	# Out is a tuple, convert to vector
		len = length(vals)

		(len > 2*nz) && error("'stretch' has more elements then allowed by image dimensions")
		(len != 1 && len != 2 && len != 6) &&
			error("Bad 'stretch' argument. It must be a 1, 2 or 6 elements array and not $len")

		val = (len == 1) ? convert(UInt16, vals)::UInt16 : convert(Array{UInt16}, vals)::Array{UInt16}
		if (len == 1)
			sc = 255 / (65535 - val)
			@inbounds @simd for k = 1:length(img)
				img[k] = (mat[k] < val) ? 0 : round(UInt8, (mat[k] - val) * sc)
			end
		elseif (len == 2)
			val = [parse(UInt16, @sprintf("%d", vals[1])) parse(UInt16, @sprintf("%d", vals[2]))]
			sc = 255 / (val[2] - val[1])
			@inbounds @simd for k = 1:length(img)
				img[k] = (mat[k] < val[1]) ? 0 : ((mat[k] > val[2]) ? 255 : round(UInt8, (mat[k]-val[1])*sc))
			end
		else	# len = 6
			nxy = nx * ny
			v1 = [1 3 5];	v2 = [2 4 6]
			sc = [255 / (val[2] - val[1]), 255 / (val[4] - val[3]), 255 / (val[6] - val[5])]
			@inbounds @simd for n = 1:nz
				@inbounds @simd for k = 1+(n-1)*nxy:n*nxy
					img[k] = (mat[k] < val[v1[n]]) ? 0 : ((mat[k] > val[v2[n]]) ? 255 : round(UInt8, (mat[k]-val[v1[n]])*sc[n]))
				end
			end
		end
	else
		sc = 255/65535
		@inbounds @simd for k = 1:length(img)
			img[k] = round(UInt8, mat[k]*sc)
		end
	end
	(haskey(d, :scale_only)) && return img			# Only the scaled array is needed. Alows it to be a view
	mat2img(img; x=x, y=y, v=v, hdr=hdr, proj4=proj4, wkt=wkt, d...)
end

# ---------------------------------------------------------------------------------------------------
function mat2img(img::GMTimage; kw...)
	# Scale a UInt16 GMTimage to UInt8. Return a new object but with all old image parameters
	(!isa(img.image, Array{UInt16}))  && return img		# Nothing to do
	I = mat2img(img.image; kw...)
	I.proj4 = img.proj4;	I.wkt = img.wkt;	I.epsg = img.epsg
	I.range = img.range;	I.inc = img.inc;	I.registration = img.registration
	I.nodata = img.nodata;	I.color_interp = img.color_interp;
	I.names = img.names;	I.metadata = img.metadata
	I.x = img.x;	I.y = img.y;	I.colormap = img.colormap;
	I.n_colors = img.n_colors;		I.alpha = img.alpha;	I.layout = img.layout;
	return I
end

# ---------------------------------------------------------------------------------------------------
# This method creates a new GMTimage but retains all the header data from the IMG object
function mat2img(mat, I::GMTimage; names::Vector{String}=String[], metadata::Vector{String}=String[])
	range = copy(I.range);	(size(mat,3) == 1) && (range[5:6] .= extrema(mat))
	GMTimage(I.proj4, I.wkt, I.epsg, range, copy(I.inc), I.registration, I.nodata, I.color_interp, metadata, names, copy(I.x), copy(I.y), zeros(size(mat,3)), mat, copy(I.colormap), I.n_colors, Array{UInt8,2}(undef,1,1), I.layout, 0)
end
function mat2img(mat, G::GMTgrid; names::Vector{String}=String[], metadata::Vector{String}=String[])
	range = copy(G.range);	range[5:6] .= (size(mat,3) == 1) ? extrema(mat) : [0., 255]
	GMTimage(G.proj4, G.wkt, G.epsg, range, copy(G.inc), G.registration, 0.0, "Gray", metadata, names, copy(G.x), copy(G.y), zeros(size(mat,3)), mat, zeros(Clong,3), 0, Array{UInt8,2}(undef,1,1), G.layout*"a", 0)
end

# ---------------------------------------------------------------------------------------------------
"""
    slicecube(I::GMTimage, layer::Int)

Take a slice of a multylayer GMTimage. Return the result still as a GMTimage. `layer` is the slice number.

### Example
Get the fourth layer of the multi-layered 'I' GMTimage object

```
I = slicecube(I, 4)
```
"""
function slicecube(I::GMTimage, layer::Int)
	(layer < 1 || layer > size(I,3)) && error("Layer number ($layer) is out of bounds of image size ($size(I,3))")
	(size(I,3) == 1) && return I		# There is nothing to slice here, but save the user fro the due deserved insult.
	mat = I.image[:,:,layer]
	range = copy(I.range);	range[5:6] .= extrema(mat)
	names = (!isempty(I.names)) ? [I.names[layer]] : I.names
	GMTimage(I.proj4, I.wkt, I.epsg, range, copy(I.inc), I.registration, I.nodata, "Gray", I.metadata, names, copy(I.x), copy(I.y), [0.], mat, zeros(Clong,3), 0, Array{UInt8,2}(undef,1,1), I.layout, I.pad)
end

# ---------------------------------------------------------------------------------------------------
"""
    stackgrids(names::Vector{String}, v=nothing; zcoord=nothing, zdim_name="time", z_unit="", save="", mirone=false)

Stack a bunch of single grids in a multiband cube like file.

- `names`: A string vector with the names of the grids to stack
- `v`: A vector with the vertical coordinates. If not provided one with 1:length(names) will be generated.
  - If `v` is a TimeType use the `z_unit` keyword to select what to store in file (case insensitive).
    - `decimalyear` or `yeardecimal` converts the DateTime to decimal years (Floa64)
	- `milliseconds` (or just `mil`) will store the DateTime as milliseconds since 0000-01-01T00:00:00 (Float64)
	- `seconds` stores the DateTime as seconds since 0000-01-01T00:00:00 (Float64)
	- `Date` or `DateTime` stores as a string representation of a DateTime.
- `zdim_name`: The name of the vertical axes (default is "time")
- `zcoord`: Keyword same as `v` (may use one or the other).
- `save`: The name of the file to be created.
- `mirone`: Does not create a cube file but instead a file named "automatic_list.txt" (or whaterver `save=xxx`)
   to be used in the Mirone `Empilhador` tool.
"""
function stackgrids(names::Vector{String}, v=nothing; zcoord=nothing, zdim_name::String="time", z_unit::String="",
                    save::String="", mirone::Bool=false)
	(v === nothing && zcoord !== nothing) && (v = zcoord)	# Accept both positional and named var for vertical coordinates

	if (isa(v, Vector{<:TimeType}))
		_z_unit = lowercase(z_unit)
		(mirone && z_unit == "") && (_z_unit = "decimalyear")		# For Mirone the default is DecimalYear
		if (_z_unit == "" || _z_unit == "decimalyear" || _z_unit == "yeardecimal")	# Actually, make it default for all and now
			v = GMT.yeardecimal.(v);				z_unit = "Decimal year"
		elseif (startswith(_z_unit, "mil"))
			v = Dates.datetime2epochms.(v);	z_unit = "Milliseconds since 0000-01-01T00:00:00"
		elseif (_z_unit == "seconds")
			v = Dates.datetime2epochms.(v) / 1000.;	z_unit = "Seconds since 0000-01-01T00:00:00"
		elseif (_z_unit == "date" || _zunit == "datetime")			# Crashes Mirone
			v = string.(v);							z_unit = "ISO Date Time"
		end
	end

	_v = isa(v, Int64) ? Float64.(v) : ((v === nothing) ? Vector{Float64}() : v)	# Int64 Crashes Mirone

	if (mirone)
		(save == "") && (save = "automatic_list.txt")
		fid = open(save, "w")
		if isempty(_v)
			[write(fid, name, "\n") for name in names]
		else
			[write(fid, names[k], "\t", string(_v[k]), "\n") for k = 1:length(names)]
		end
		close(fid)
		return nothing
	end

	G = gmtread(names[1], grid=true)	# reading with GDAL f screws sometimes with the "is not a Latitude/Y dimension."
	x, y, range, inc = G.x, G.y, G.range, G.inc		# So read first with GMT anf keep only the coords info.
	G = gdaltranslate(names[1])
	mat = Array{eltype(G)}(undef, size(G,1), size(G,2), length(names))
	mat[:,:,1] = G.z
	for k = 2:length(names)
		G = gdaltranslate(names[k])
		mat[:,:,k] = G.z
	end
	cube = mat2grid(mat, G)
	cube.x = x;		cube.y = y;		cube.range = range;		cube.inc = inc
	cube.z_unit = z_unit
	(isempty(_v) || eltype(_v) == String) ? append!(cube.range, [0., 1.]) : append!(cube.range, [_v[1], _v[end]])
	cube.names = names;		cube.v = _v
	(save != "") && gdalwrite(cube, save, _v, dim_name=zdim_name)
	return (save != "") ? nothing : cube
end

# ---------------------------------------------------------------------------------------------------
"""
    I = image_alpha!(img::GMTimage; alpha_ind::Integer, alpha_vec::Vector{Integer}, alpha_band::UInt8)

Change the alpha transparency of the GMTimage object `img`. If the image is indexed, one can either
change just the color index that will be made transparent by uing `alpha_ind=n` or provide a vector
of transaparency values in the range [0 255]; This vector can be shorter than the orginal number of colors.
Use `alpha_band` to change, or add, the alpha of true color images (RGB).

    Example1: change to the third color in cmap to represent the new transparent color
        image_alpha!(img, alpha_ind=3)

    Example2: change to the first 6 colors in cmap by assigning them random values
        image_alpha!(img, alpha_vec=round.(Int32,rand(6).*255))
"""
function image_alpha!(img::GMTimage; alpha_ind=nothing, alpha_vec=nothing, alpha_band=nothing)
	# Change the alpha transparency of an image
	n_colors = img.n_colors
	if (n_colors > 100000)  n_colors = Int(floor(n_colors / 1000))  end
	if (alpha_ind !== nothing)			# Change the index of the alpha color
		(alpha_ind < 0 || alpha_ind > 255) && error("Alpha color index must be in the [0 255] interval")
		img.n_colors = n_colors * 1000 + Int32(alpha_ind)
	elseif (alpha_vec !== nothing)		# Replace/add the alpha column of the colormap matrix. Allow also shorter vectors
		@assert(isa(alpha_vec, Vector{<:Integer}))
		(length(alpha_vec) > n_colors) && error("Length of alpha vector is larger than the number of colors")
		n_col = div(length(img.colormap), n_colors)
		vec = convert.(Int32, alpha_vec)
		if (n_col == 4)  img.colormap[(end-length(vec)+1):end] = vec;
		else             img.colormap = [img.colormap; [vec[:]; round.(Int32, ones(n_colors - length(vec)) .* 255)]]
		end
		img.n_colors = n_colors * 1000
	elseif (alpha_band !== nothing)		# Replace the entire alpha band
		@assert(eltype(alpha_band) == UInt8)
		ny1, nx1, = size(img.image)
		ny2, nx2  = size(alpha_band)
		(ny1 != ny2 || nx1 != nx2) && error("alpha channel has wrong dimensions")
		(size(img.image, 3) != 3) ? @warn("Adding alpha band is restricted to true color images (RGB)") :
		                            img.alpha = (isa(alpha_band,GMTimage)) ? alpha_band.image : alpha_band
	end
	return nothing
end

# ---------------------------------------------------------------------------------------------------
"""
    image_cpt!(img::GMTimage, cpt::GMTcpt, clear=false)
or

    image_cpt!(img::GMTimage, cpt::String, clear=false)

Add (or replace) a colormap to a GMTimage object from the colors in the cpt.
This should have effect only if IMG is indexed.
Use `image_cpt!(img, clear=true)` to remove a previously existent `colormap` field in IMG
"""
image_cpt!(img::GMTimage, cpt::String) = image_cpt!(img, gmtread(cpt))
function image_cpt!(img::GMTimage, cpt::GMTcpt)
	# Insert the cpt info in the img.colormap member
	n = 1
	colormap = fill(Int32(255), size(cpt.colormap,1) * 4)
	for k = 1:size(cpt.colormap,1)
		colormap[n:n+2] = round.(Int32, cpt.colormap[k,:] .* 255);	n += 4
	end
	img.colormap = colormap
	img.n_colors = size(cpt.colormap,1)
	return nothing
end
function image_cpt!(img::GMTimage; clear::Bool=true)
	if (clear)
		img.colormap, img.n_colors = fill(Int32(0), 3), 0
	end
	return nothing
end

# ---------------------------------------------------------------------------------------------------
"""
    I = ind2rgb(I)

Convert an indexed image I to RGB. It uses the internal colormap to do the conversion.
"""
function ind2rgb(img::GMTimage)
	# ...
	(size(img.image, 3) >= 3) && return img 	# Image is already RGB(A)
	if (img.n_colors == 0)				# If no cmap just replicate the first layer.
		imgRGB = repeat(img.image, 1, 1, 3)
	else
		imgRGB = zeros(UInt8,size(img.image,1), size(img.image,2), 3)
		n = 1
		for k = 1:length(img.image)
			start_c = img.image[k] * 4
			for c = 1:3
				imgRGB[n] = img.colormap[start_c+c];	n += 1
			end
		end
	end
	mat2img(imgRGB, x=img.x, y=img.y, proj4=img.proj4, wkt=img.wkt, mem_layout=img.layout)
end

# ---------------------------------------------------------------------------------------------------
"""
    G = mat2grid(mat; reg=nothing, x=[], y=[], hdr=nothing, proj4::String="", wkt::String="", tit::String="",
                 rem::String="", cmd::String="", names::Vector{String}=String[], scale::Float32=1f0, offset::Float32=0f0)

Take a 2/3D `mat` array and a HDR 1x9 [xmin xmax ymin ymax zmin zmax reg xinc yinc] header descriptor
and return a grid GMTgrid type.
Alternatively to HDR, provide a pair of vectors, x & y, with the X and Y coordinates.
Optionaly, the HDR arg may be ommited and it will computed from `mat` alone, but then x=1:ncol, y=1:nrow
When HDR is not used, REG == nothing [default] means create a gridline registration grid and REG == 1,
or REG="pixel" a pixel registered grid.

For 3D arrays the `names` option is used to give a description for each layer (also saved to file when using a GDAL function).

The `scale` and `offset` options are used when `mat` is an Integer type and we want to save the grid with an scale/offset.

Other methods of this function do:

    G = mat2grid([val]; hdr=hdr_vec, reg=nothing, proj4::String="", wkt::String="", tit::String="", rem::String="")

Create Float GMTgrid with size, coordinates and increment determined by the contents of the HDR var. This
array, which is now MANDATORY, has either the same meaning as above OR, alternatively, containng only
[xmin xmax ymin ymax xinc yinc]
VAL is the value that will be fill the matrix (default VAL = Float32(0)). To get a Float64 array use, for
example, VAL = 1.0 Ay other non Float64 will be converted to Float32

    Example: mat2grid(1, hdr=[0. 5 0 5 1 1])

    G = mat2grid(f::Function, x, y; reg=nothing, proj4::String="", wkt::String="", epsg::Int=0, tit::String="", rem::String="")

Where F is a function and X,Y the vectors coordinates defining it's domain. Creates a Float32 GMTgrid with
size determined by the sizes of the X & Y vectors.

    Example: f(x,y) = x^2 + y^2;  G = mat2grid(f, x = -2:0.05:2, y = -2:0.05:2)

    G = mat2grid(f::String, x=[], y=[])

Whre F is a pre-set function name. Currently available:
   - "ackley", "eggbox", "sombrero", "parabola" and "rosenbrock"
X,Y are vectors coordinates defining the function's domain, but default values are provided for each function.
creates a Float32 GMTgrid.

    Example: G = mat2grid("sombrero")
"""
function mat2grid(val::Real=Float32(0); reg=nothing, hdr=nothing, proj4::String="", wkt::String="", epsg::Int=0,
	tit::String="", rem::String="", names::Vector{String}=String[])

	(hdr === nothing) && error("When creating grid type with no data the 'hdr' arg cannot be missing")
	(!isa(hdr, Array{Float64})) && (hdr = Float64.(hdr))
	(!isa(val, AbstractFloat)) && (val = Float32(val))		# We only want floats here
	if (length(hdr) == 6)
		hdr = [hdr[1], hdr[2], hdr[3], hdr[4], val, val, reg === nothing ? 0. : 1., hdr[5], hdr[6]]
	end
	mat2grid([nothing val]; reg=reg, hdr=hdr, proj4=proj4, wkt=wkt, epsg=epsg, tit=tit, rem=rem, cmd="", names=names)
end

function mat2grid(mat, xx=Vector{Float64}(), yy=Vector{Float64}(); reg=nothing, x=Vector{Float64}(), y=Vector{Float64}(),
	v=Vector{Float64}(), hdr=nothing, proj4::String="", wkt::String="", epsg::Int=0, tit::String="", rem::String="",
	cmd::String="", names::Vector{String}=String[], scale::Float32=1f0, offset::Float32=0f0)
	# Take a 2/3D array and turn it into a GMTgrid

	!isa(mat[2], Real) && error("input matrix must be of Real numbers")
	reg_ = 0
	if (isa(reg, String) || isa(reg, Symbol))
		t = lowercase(string(reg))
		reg_ = (t != "pixel") ? 0 : 1
	elseif (isa(reg, Real))
		reg_ = (reg == 0) ? 0 : 1
	end
	if (isempty(x) && !isempty(xx))  x = xx  end
	if (isempty(y) && !isempty(yy))  y = yy  end
	x, y, hdr, x_inc, y_inc = grdimg_hdr_xy(mat, reg_, hdr, x, y)

	# Now we still must check if the method with no input MAT was called. In that case mat = [nothing val]
	# and the MAT must be finally computed.
	nx = size(mat, 2);		ny = size(mat, 1);
	if (ny == 1 && nx == 2 && mat[1] === nothing)
		fill_val = mat[2]
		mat = zeros(eltype(fill_val), length(y), length(x))
		(fill_val != 0) && fill!(mat, fill_val)
	end

	GMTgrid(proj4, wkt, epsg, hdr[1:6], [x_inc, y_inc], reg_, NaN, tit, rem, cmd, names, x, y, v, mat, "x", "y", "v", "z", "BCB", scale, offset, 0)
end

# This method creates a new GMTgrid but retains all the header data from the G object
function mat2grid(mat, G::GMTgrid)
	Go = GMTgrid(G.proj4, G.wkt, G.epsg, deepcopy(G.range), deepcopy(G.inc), G.registration, G.nodata, G.title, G.remark, G.command, String[], deepcopy(G.x), deepcopy(G.y), [0.], mat, G.x_unit, G.y_unit, G.v_unit, G.z_unit, G.layout, 1f0, 0f0, G.pad)
	grd_min_max!(Go)		# Also take care of NaNs
	Go
end
function mat2grid(mat, I::GMTimage)
	Go = GMTgrid(I.proj4, I.wkt, I.epsg, I.range, I.inc, I.registration, I.nodata, "", "", "", String[], I.x, I.y, [0.], mat, "", "", "", "", I.layout, 1f0, 0f0, I.pad)
	(length(Go.layout) == 4) && (Go.layout = Go.layout[1:3])	# No space for the a|A
	grd_min_max!(Go)		# Also take care of NaNs
	Go
end

function mat2grid(f::Function, x, y; reg=nothing, proj4::String="", wkt::String="", epsg::Int=0, tit::String="", rem::String="")
	z = Array{Float32,2}(undef,length(y),length(x))
	for i = 1:length(x)
		for j = 1:length(y)
			z[j,i] = f(x[i],y[j])
		end
	end
	mat2grid(z; reg=reg, x=x, y=y, proj4=proj4, wkt=wkt, epsg=epsg, tit=tit, rem=rem)
end

function mat2grid(f::String, x=Vector{Float64}(), y=Vector{Float64}())
	# Something is very wrong here. If I add named vars it annoyingly warns
	#	WARNING: Method definition f2(Any, Any) in module GMT at C:\Users\joaqu\.julia\dev\GMT\src\gmt_main.jl:1556 overwritten on the same line.
	if (startswith(f, "ack"))				# Ackley (inverted) https://en.wikipedia.org/wiki/Ackley_function
		f_ack(x,y) = 20 * exp(-0.2 * sqrt(0.5 * (x^2 + y^2))) + exp(0.5*(cos(2pi*x) + cos(2pi*y))) - 22.718281828459045
		if (isempty(x))  x = -5:0.05:5;	y = -5:0.05:5;  end
		mat2grid(f_ack, x, y)
	elseif (startswith(f, "egg"))
		f_egg(x, y) = (sin(x*10) + cos(y*10)) / 4
		if (isempty(x))  x = -1:0.01:1;	y = -1:0.01:1;  end
		mat2grid(f_egg, x, y)
	elseif (startswith(f, "para"))
		f_parab(x,y) = x^2 + y^2
		if (isempty(x))  x = -2:0.05:2;	y = -2:0.05:2;  end
		mat2grid(f_parab, x, y)
	elseif (startswith(f, "rosen"))			# rosenbrock
		f_rosen(x,y) = (1 - x)^2 + 100 * (y - x^2)^2
		if (isempty(x))  x = -2:0.05:2;	y = -1:0.05:3;  end
		mat2grid(f_rosen, x, y)
	elseif (startswith(f, "somb"))			# sombrero
		f_somb(x,y) = cos(sqrt(x^2 + y^2) * 2pi / 8) * exp(-sqrt(x^2 + y^2) / 10)
		if (isempty(x))  x = -15:0.2:15;	y = -15:0.2:15;  end
		mat2grid(f_somb, x, y)
	else
		@warn("Unknown surface '$f'. Just giving you a parabola.")
		mat2grid("para")
	end
end

# ---------------------------------------------------------------------------------------------------
function grdimg_hdr_xy(mat, reg, hdr, x=Vector{Float64}(), y=Vector{Float64}())
# Generate x,y coords array and compute/update header plus increments for grids/images
	nx = size(mat, 2);		ny = size(mat, 1);

	if (!isempty(x) && !isempty(y))		# But not tested if they are equi-spaced as they MUST be
		if ((length(x) != (nx+reg) || length(y) != (ny+reg)) && (length(x) != 2 || length(y) != 2))
			error("size of x,y vectors incompatible with 2D array size")
		end
		one_or_zero = reg == 0 ? 1 : 0
		if (length(x) != 2)				# Check that REGistration and coords are compatible
			(reg == 1 && round((x[end] - x[1]) / (x[2] - x[1])) != nx) &&		# Gave REG = pix but xx say grid
				(@warn("Gave REGistration = 'pixel' but X coordinates say it's gridline. Keeping later reg."); one_or_zero = 1)
		else
			x = collect(range(x[1], stop=x[2], length=nx+reg))
			y = collect(range(y[1], stop=y[2], length=ny+reg))
		end
		x_inc = (x[end] - x[1]) / (nx - one_or_zero)
		y_inc = (y[end] - y[1]) / (ny - one_or_zero)
		zmin, zmax = extrema_nan(mat)
		hdr = [x[1], x[end], y[1], y[end], zmin, zmax]
	elseif (hdr === nothing)
		zmin, zmax = extrema_nan(mat)
		if (reg == 0)  x = collect(1.0:nx);		y = collect(1.0:ny)
		else           x = collect(0.5:nx+0.5);	y = collect(0.5:ny+0.5)
		end
		hdr = [x[1], x[end], y[1], y[end], zmin, zmax]
		x_inc = 1.0;	y_inc = 1.0
	else
		(length(hdr) != 9) && error("The HDR array must have 9 elements")
		(!isa(hdr, Array{Float64})) && (hdr = Float64.(hdr))
		one_or_zero = (hdr[7] == 0) ? 1 : 0
		if (ny == 1 && nx == 2 && mat[1] === nothing)
			# In this case the 'mat' is a tricked matrix with [nothing val]. Compute nx,ny from header
			# The final matrix will be computed in the main mat2grid method
			nx = Int(round((hdr[2] - hdr[1]) / hdr[8] + one_or_zero))
			ny = Int(round((hdr[4] - hdr[3]) / hdr[9] + one_or_zero))
		end
		x = collect(range(hdr[1], stop=hdr[2], length=(nx+Int(hdr[7])) ))
		y = collect(range(hdr[3], stop=hdr[4], length=(ny+Int(hdr[7])) ))
		# Recompute the x|y_inc to make sure they are right.
		x_inc = (hdr[2] - hdr[1]) / (nx - one_or_zero)
		y_inc = (hdr[4] - hdr[3]) / (ny - one_or_zero)
	end
	if (isa(x, UnitRange))  x = collect(x)  end			# The AbstractArrays are much less forgivable
	if (isa(y, UnitRange))  y = collect(y)  end
	if (!isa(x, Vector{Float64}))  x = Float64.(x)  end
	if (!isa(y, Vector{Float64}))  y = Float64.(y)  end
	return x, y, hdr, x_inc, y_inc
end

# ---------------------------------------------------------------------------------------------------
# Convert the HDR vector from grid to pixel registration or vice-versa
grid2pix(GI; pix=true) = grid2pix([GI.range; GI.registration; GI.inc], pix=pix)
function grid2pix(hdr::Vector{Float64}; pix=true)
	((pix && hdr[7] == 1) || (!pix && hdr[7] == 0)) && return hdr 		# Nothing to do
	if (pix)  hdr[1] -= hdr[8]/2; hdr[2] += hdr[8]/2; hdr[3] -= hdr[9]/2; hdr[4] += hdr[9]/2;	hdr[7] = 1.
	else      hdr[1] += hdr[8]/2; hdr[2] -= hdr[8]/2; hdr[3] += hdr[9]/2; hdr[4] -= hdr[9]/2;	hdr[7] = 0.
	end
	return hdr
end

#= ---------------------------------------------------------------------------------------------------
function mksymbol(f::Function, cmd0::String="", arg1=nothing; kwargs...)
	# Make a fig and convert it to EPS so it can be used as a custom symbol is plot(3)
	d = KW(kwargs)
	t = ((val = find_in_dict(d, [:symbname :symb_name :symbol])[1]) !== nothing) ? string(val) : "GMTsymbol"
	d[:savefig] = t * ".eps"
	f(cmd0, arg1; d...)
end
mksymbol(f::Function, arg1; kw...) = mksymbol(f, "", arg1; kw...)
=#

# ---------------------------------------------------------------------------------------------------
"""
    make_zvals_vec(D, user_ids::Vector{String}, vals::Array{<:Real}, sub_head=0, upper=false, lower=false)

- `user_ids` -> is a string vector with the ids (names in header) of the GMTdataset D
- `vals`     -> is a vector with the the numbers to be used in plot -Z to color the polygons.
- `sub_head` -> Position in header where field is to be found in the comma separated string.
Create a vector with ZVALS to use in plot where length(ZVALS) == length(D)
The elements of ZVALS are made up from the `vals` but it can be larger if there are segments with
no headers. In that case it replicates the previously known value until it finds a new segment ID.

Returns a Vector{Float64} with the same length as the number of segments in D. The content is
made up after the contents of `vals` but repeated such that each polygon of the same family, i.e.
with the same `user_ids`, has the same value.
"""
function make_zvals_vec(D, user_ids::Vector{String}, vals::Array{<:Real}, sub_head::Int=0, case::Int=0)::Vector{Float64}

	n_user_ids = length(user_ids)
	@assert(n_user_ids == length(vals))
	data_ids, ind = get_segment_ids(D, case)
	(ind[1] != 1) && error("This function requires that first segment has a a header with an id")
	n_data_ids = length(data_ids)
	(n_user_ids > n_data_ids) &&
		@warn("Number of segment IDs requested is larger than segments with headers in data")

	if (sub_head != 0)
		[data_ids[k] = split(data_ids[k],',')[sub_head]  for k = 1:length(ind)]
	end

	n_seg = (isa(D, Array)) ? length(D) : 1
	zvals = fill(NaN, n_seg)
	n = 1
	for k = 1:n_data_ids
		for m = 1:n_user_ids
			if startswith(data_ids[k], user_ids[m])			# Find first occurence of user_ids[k] in a segment header
				last = (k < n_data_ids) ? ind[k+1]-1 : n_seg
				[zvals[j] = vals[m] for j = ind[k]:last]		# Repeat the last VAL for segments with no headers
				#println("k = ", k, " m = ",m, " pol_id = ", data_ids[k], ";  usr id = ", user_ids[m], " Racio = ", vals[m], " i1 = ", ind[k], " i2 = ",last)
				n = last + 1					# Prepare for next new VAL
				break
			end
		end
	end
	return zvals
end

# ---------------------------------------------------------------------------------------------------
function edit_segment_headers!(D, vals::Array, opt::String)
	# Add an option OPT to segment headers with a val from VALS. Number of elements of VALS must be
	# equal to the number of segments in D that have a header. If numel(val) == 1 must encapsulate it in []

	ids, ind = get_segment_ids(D)
	if (isa(D, Array))
		[D[ind[k]].header *= string(opt, vals[k])  for k = 1:length(ind)]
	else
		D.header *= string(opt, vals[1])
	end
	return nothing
end

# ---------------------------------------------------------------------------------------------------
"""
    ids, ind = get_segment_ids(D, case=0)::Tuple{Vector{String}, Vector{Int}}

Where D is a GMTdataset of a vector of them, returns the segment ids (first text after the '>') and
the idices of those segments.
"""
function get_segment_ids(D, case::Int=0)::Tuple{Vector{String}, Vector{Int}}
	# Get segment ids (first text after the '>') and the idices of those segments
	# CASE -> If == 1 force return in LOWER case. If == 2 force upper case. Default (case = 0) dosen't touch
	if (isa(D, Array))  n = length(D);	d = Dict(k => D[k].header for k = 1:n)
	else                n = 1;			d = Dict(1 => D.header)
	end
	tf = Vector{Bool}(undef,n)					# pre-allocate
	[tf[k] = (d[k] !== "" && d[k][1] != ' ') ? true : false for k = 1:n];	# Mask of non-empty headers
	ind = 1:n
	ind = ind[tf]			# OK, now we have the indices of the segments with headers != ""
	ids = Vector{String}(undef,length(ind))		# pre-allocate
	if (case == 1)
		[ids[k] = lowercase(d[ind[k]]) for k = 1:length(ind)]	# indices of non-empty segments
	elseif (case == 2)
		[ids[k] = uppercase(d[ind[k]]) for k = 1:length(ind)]
	else
		[ids[k] = d[ind[k]] for k = 1:length(ind)]
	end
	return ids, ind
end


text:
function add_modifier(obj::ObservedMeasurement, arg0::EstimationModifier)
    return jcall(obj, "addModifier", void, (EstimationModifier,), arg0)
end

function estimate(obj::ObservedMeasurement, arg0::jint, arg1::jint, arg2::Vector{SpacecraftState})
    return jcall(obj, "estimate", EstimatedMeasurement, (jint, jint, Vector{SpacecraftState}), arg0, arg1, arg2)
end

function get_base_weight(obj::ObservedMeasurement)
    return jcall(obj, "getBaseWeight", Vector{jdouble}, ())
end

function get_dimension(obj::ObservedMeasurement)
    return jcall(obj, "getDimension", jint, ())
end

function get_modifiers(obj::ObservedMeasurement)
    return jcall(obj, "getModifiers", List, ())
end

function get_parameters_drivers(obj::ObservedMeasurement)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_satellites(obj::ObservedMeasurement)
    return jcall(obj, "getSatellites", List, ())
end

function get_theoretical_standard_deviation(obj::ObservedMeasurement)
    return jcall(obj, "getTheoreticalStandardDeviation", Vector{jdouble}, ())
end

function is_enabled(obj::ObservedMeasurement)
    return jcall(obj, "isEnabled", jboolean, ())
end

function set_enabled(obj::ObservedMeasurement, arg0::jboolean)
    return jcall(obj, "setEnabled", void, (jboolean,), arg0)
end


text: include("TestProblems/TestProblems.jl")


text:
function init_constraint_struct(com, set::AbstractBoolSet{F1,F2}, internals) where {F1,F2}
    f = MOIU.eachscalar(internals.fct)

    lhs_fct = f[1:set.lhs_dimension]
    rhs_fct = f[end-set.rhs_dimension+1:end]

    if F1 <: MOI.ScalarAffineFunction
        lhs_fct = get_saf(lhs_fct)
    end
    if F2 <: MOI.ScalarAffineFunction
        rhs_fct = get_saf(rhs_fct)
    end

    if F1 <: MOI.VectorOfVariables
        lhs_fct = get_vov(lhs_fct)
    end
    if F2 <: MOI.VectorOfVariables
        rhs_fct = get_vov(rhs_fct)
    end

    lhs = get_constraint(com, lhs_fct, set.lhs_set)
    rhs = get_constraint(com, rhs_fct, set.rhs_set)

    bs = bool_constraint(set, com, internals, lhs, rhs)

    return bool_constraint(set, com, internals, lhs, rhs)
end

function bool_constraint(::XorSet, com, internals, lhs, rhs)
    XorConstraint(
        internals,
        BoolConstraintInternals(lhs, rhs),
        get_complement_constraint(com, lhs),
        get_complement_constraint(com, rhs)
    )
end

for (set, bool_data) in BOOL_SET_TO_CONSTRAINT
    get(bool_data, :specific_constraint, false) && continue
    @eval begin
        function bool_constraint(::$set, com, internals, lhs, rhs)
            $(bool_data.constraint)(
                com,
                internals,
                lhs,
                rhs
            )
        end
    end
end

"""
    demorgan_complement_set(::Type{<:AbstractBoolSet})

Return the type of the demorgan complement bool set so AndSet <-> OrSet
"""
demorgan_complement_set(::Type{<:AbstractBoolSet}) = nothing
demorgan_complement_set(::Type{<:AndSet}) = OrSet
demorgan_complement_set(::Type{<:OrSet}) = AndSet

demorgan_complement_constraint_type(::Type{<:AbstractBoolSet}) = nothing
demorgan_complement_constraint_type(::Type{<:AndSet}) = OrConstraint
demorgan_complement_constraint_type(::Type{<:OrSet}) = AndConstraint

"""
    complement_set(::Type{<:AbstractBoolSet})

Return the type of the complement bool set so XorSet => XNorSet
"""
complement_set(::Type{<:AbstractBoolSet}) = nothing
complement_set(::Type{<:XorSet}) = XNorSet
complement_set(::Type{<:XNorSet}) = XorSet

"""
    complement_constraint_type(::Type{<:AbstractBoolSet})

Return the constraint of the complement bool set so XorSet => XNorConstraint
"""
complement_constraint_type(::Type{<:AbstractBoolSet}) = nothing
complement_constraint_type(::Type{<:XorSet}) = XNorConstraint
complement_constraint_type(::Type{<:XNorSet}) = XorConstraint

for (set, bool_data) in BOOL_SET_TO_CONSTRAINT
    res_op = get(bool_data, :res_op, :identity)
    if get(bool_data, :needs_call, false)
        @eval begin
            function apply_bool_operator(::Type{<:$set}, lhs_fct, rhs_fct, args...)
                $(Expr(:call, res_op, Expr(:call, bool_data.op, :(lhs_fct(args...)), :(rhs_fct(args...)))))
            end
        end
    else
        @eval begin
            function apply_bool_operator(::Type{<:$set}, lhs_fct, rhs_fct, args...)
                $(Expr(:call, res_op, Expr(bool_data.op, :(lhs_fct(args...)), :(rhs_fct(args...)))))
            end
        end
    end
end

function init_constraint!(
    com::CS.CoM,
    constraint::BoolConstraint,
    fct,
    set::AbstractBoolSet;
)
   init_lhs_and_rhs!(com, constraint, fct, set)
end

function init_lhs_and_rhs!(
    com::CS.CoM,
    constraint::BoolConstraint,
    fct,
    set::AbstractBoolSet;
)
    lhs_feasible = init_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_feasible = init_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    # check only for && constraint if both are feasible
    if set isa AndSet
        return lhs_feasible && rhs_feasible
    end
    return true
end

"""
    is_constraint_solved(
        constraint::BoolConstraint,
        fct,
        set::AbstractBoolSet,
        values::Vector{Int},
    )

Check if the constraint is solved gived the `values`
"""
function is_constraint_solved(
    constraint::BoolConstraint,
    fct,
    set::AbstractBoolSet,
    values::Vector{Int},
)
    apply_bool_operator(typeof(set), lhs_solved, rhs_solved, constraint, values)
end

function lhs_solved(constraint::BoolConstraint, values::Vector{Int})
    lhs_num_vars = get_num_vars(constraint.lhs.fct)
    return is_constraint_solved(constraint.lhs, constraint.lhs.fct, constraint.lhs.set, values[1:lhs_num_vars])
end

function rhs_solved(constraint::BoolConstraint, values::Vector{Int})
    rhs_num_vars = get_num_vars(constraint.rhs.fct)
    return is_constraint_solved(constraint.rhs, constraint.rhs.fct, constraint.rhs.set, values[end-rhs_num_vars+1:end])
end

function is_lhs_constraint_violated(com, constraint)
    is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
end

function is_rhs_constraint_violated(com, constraint)
    is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
end

"""
    activate_lhs!(com, constraint::BoolConstraint)

Activate the lhs constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_lhs!(com, constraint::BoolConstraint)
    lhs = constraint.lhs
    if !constraint.lhs_activated
        !activate_constraint!(com, lhs, lhs.fct, lhs.set) && return false
        constraint.lhs_activated = true
        constraint.lhs_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    activate_rhs!(com, constraint::BoolConstraint)

Activate the rhs constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_rhs!(com, constraint::BoolConstraint)
    rhs = constraint.rhs
    if !constraint.rhs_activated
        !activate_constraint!(com, rhs, rhs.fct, rhs.set) && return false
        constraint.rhs_activated = true
        constraint.rhs_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

function single_reverse_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::AbstractBoolSet,
    var::Variable,
    backtrack_idx::Int,
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        # the variable must be part of the inner constraint
        # Todo: Speed up the `in` here
        if var.idx in inner_constraint.indices
            single_reverse_pruning_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                var,
                backtrack_idx,
            )
        end
    end
end

function reverse_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::AbstractBoolSet,
    backtrack_id::Int,
) where {
    T<:Real,
}
    # check if inner constraint should be deactived again
    if constraint.lhs_activated && backtrack_id == constraint.lhs_activated_in_backtrack_idx
        constraint.lhs_activated = false
        constraint.lhs_activated_in_backtrack_idx = 0
    end
    if constraint.rhs_activated && backtrack_id == constraint.rhs_activated_in_backtrack_idx
        constraint.rhs_activated = false
        constraint.rhs_activated_in_backtrack_idx = 0
    end

    for inner_constraint in (constraint.lhs, constraint.rhs)
        reverse_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            backtrack_id,
        )
    end
end

function restore_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::AbstractBoolSet,
    prune_steps::Union{Int,Vector{Int}},
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        restore_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            prune_steps,
        )
    end
end

function finished_pruning_constraint!(
    com::CS.CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::AbstractBoolSet,
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        finished_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
        )
    end
end

function changed!(com::CS.CoM, constraint::BoolConstraint, fct, set)
    lhs = constraint.lhs
    changed!(com, lhs, lhs.fct, lhs.set)
    rhs = constraint.rhs
    changed!(com, rhs, rhs.fct, rhs.set)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 38db4d40-8d1b-11eb-37fc-89b3bcb5f9c6
begin
import Pkg
	Pkg.activate(mktempdir())
    Pkg.add("InteractiveUtils")
  using InteractiveUtils
    Pkg.add("PlutoUI")
  using PlutoUI

end

# ╔═╡ 7a6936be-8d16-11eb-0a34-c7e3704e4e77
md"""
# Hello!``~~~`` I am your ``~``🛡 Pluto Notebook!
* Scientists use Notebooks to find new things.
* You can have fun finding and learning new things.
* Sliders and Scrubbers let you try.``~~``Don't worry.
* You do not need to erase,``~~`` Just try again.
* Practice and you will be a Notebook Scientist!
"""

# ╔═╡ e9f3d90a-8c81-11eb-2adc-5fb66a6341e5
word1List = [" ____ " "cow" "   cat   " "dog" "rooster" "pig";
" ____ " "said oink oink" "said moo" "said meow" "said cock-a-doodle-doo" "said bow bow"];

# ╔═╡ 93322ebe-8c85-11eb-2ea7-15d8ec2977aa
md"""
!!! first
#### ``~~~~~~~~~~~~~~``Move the sliders to pick the words that go together.
##### ``~~~~~~~~~~~~~~``These are called word pairs. ``~~`` Think of some more word pairs.
noun
$(@bind word1 Slider(1:6; default=1, show_value=false))
verb + word pair
$(@bind word2 Slider(1:6; default=1, show_value=false))
"""

# ╔═╡ 5c1cff3a-8c83-11eb-3067-4bf2500585d4
md"""
# The $(word1List[1, word1]) $(word1List[2, word2]) .
"""

# ╔═╡ 5600460c-8c92-11eb-0d30-6ffbc6077b55
posList = ["_______" "article" "noun" "pronoun" "verb" "adjective" "adverb" "preposition" "conjunction"  "exclamation"];

# ╔═╡ 017d84c4-8cac-11eb-23ad-53aeacf5028b
md"""
``\mathbf1=``
$(@bind pos1 Scrubbable(1 : 10; default=1)) ``~~~~~~~~~~\mathbf2=``
$(@bind pos2 Scrubbable(1 : 10; default=1)) ``~~~~~~~~~~\mathbf3=``
$(@bind pos3 Scrubbable(1 : 10; default=1)) ``~~~~~~~~~~\mathbf4=``
$(@bind pos4 Scrubbable(1 : 10; default=1)) ``~~~~~~~~~~\mathbf5=``
$(@bind pos5 Scrubbable(1 : 10; default=1)) ``~~~~~~~~~~\mathbf6=``
$(@bind pos6 Scrubbable(1 : 10; default=1))
!!! study
###### Left click and drag the text,``~``  or use arrow keys.``~~~``The mouse wheel zooms.
"""

# ╔═╡ b27de8d8-8c8f-11eb-280f-a7eaad09372d
md"""
!!! next
##### ``~~~~~~~~~``Click and Drag a scrub number to  pick the part of speech.
|_____1______|_____2______|_____3______|_____4______|_____5______|_____6______|
|:----------:|:-------:|:--------:|:---------:|:---------:|:-------:|
| The   | cow | jumped  |over|the|moon.|
| $(posList[pos1])|$(posList[pos2]) |  $(posList[pos3])|$(posList[pos4])|$(posList[pos5]) |  $(posList[pos6])|

"""

# ╔═╡ c4a98619-8ca0-4357-ba03-24976e65762c
html"""

Zoom & Move

"""

# ╔═╡ Cell order:
# ╠═38db4d40-8d1b-11eb-37fc-89b3bcb5f9c6
# ╟─7a6936be-8d16-11eb-0a34-c7e3704e4e77
# ╟─e9f3d90a-8c81-11eb-2adc-5fb66a6341e5
# ╟─93322ebe-8c85-11eb-2ea7-15d8ec2977aa
# ╠═5c1cff3a-8c83-11eb-3067-4bf2500585d4
# ╟─5600460c-8c92-11eb-0d30-6ffbc6077b55
# ╟─b27de8d8-8c8f-11eb-280f-a7eaad09372d
# ╠═017d84c4-8cac-11eb-23ad-53aeacf5028b
# ╟─c4a98619-8ca0-4357-ba03-24976e65762c


text:
@system Area begin
    NU: nodal_units ~ hold
    PD: planting_density ~ hold

    leaf_area(x=NU["*"].leaf.area) => sum(x) ~ track(u"cm^2")

    green_leaf_area(x=NU["*"].leaf.green_area) => sum(x) ~ track(u"cm^2")

    green_leaf_ratio(green_leaf_area, leaf_area) => begin
        iszero(leaf_area) ? 0. : (green_leaf_area / leaf_area)
    end ~ track

    LAI(green_leaf_area, PD): leaf_area_index => begin
        green_leaf_area * PD
    end ~ track

    # actualgreenArea is the green area of leaf growing under carbon limitation
    #SK 8/22/10: There appears to be no distinction between these two variables in the code.
    actual_green_leaf_area(green_leaf_area) => green_leaf_area ~ track(u"cm^2")

    senescent_leaf_area(x=NU["*"].leaf.senescent_area) => sum(x) ~ track(u"cm^2")

    potential_leaf_area(x=NU["*"].leaf.potential_area) => sum(x) ~ track(u"cm^2")

    potential_leaf_area_increase(x=NU["*"].leaf.potential_area_increase) => sum(x) ~ track(u"cm^2")

    # calculate relative area increases for leaves now that they are updated
    #TODO remove if unnecessary
    # relative_leaf_area_increase(x=NU["*"].leaf.relative_area_increase) => sum(x) ~ track(u"cm^2")

    #FIXME it doesn't seem to be 'actual' dropped leaf area
    # calculated dropped leaf area YY
    dropped_leaf_area(potential_leaf_area, green_leaf_area) => (potential_leaf_area - green_leaf_area) ~ track(u"cm^2")
end