text:
using Documenter
using RealNeuralNetworks
using RealNeuralNetworks.Neurons
using RealNeuralNetworks.Neurons.Segments
using RealNeuralNetworks.Neurons.Segments.Synapses
using RealNeuralNetworks.SWCs

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

makedocs(
	modules=[RealNeuralNetworks, Neurons, Segments, Synapses, SWCs],
	sitename="RealNeuralNetworks.jl",
	authors="Jingpeng Wu",
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://seung-lab.github.io/RealNeuralNetworks.jl/latest/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-2",
        highlights = ["yaml"],
    ),
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Getting Started" => "man/getting_started.md",
            "man/examples.md",
            "man/references.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            hide("Internals" => "lib/internals.md", Any[
                "lib/internals/Manifests.md",
                "lib/internals/NodeNets.md",
                "lib/internals/PointArrays.md",
            ])
        ],
        "man/contributing.md",
    ],
)

deploydocs(
	repo="github.com/seung-lab/RealNeuralNetworks.jl",
	target="build",
	deps=nothing,
	make=nothing
)


text:
@doc doc"""

```math
dS = μSdt + \sqrt{v}SdW_1
dv = κ(Θ-v)dt + σ\sqrt{v}dW_2
dW_1 dW_2 = ρ dt
```

"""
function HestonProblem(μ,κ,Θ,σ,ρ,u0,tspan;seed=UInt64(0),kwargs...)
  f = function (du,u,p,t)
    du[1] = μ*u[1]
    du[2] = κ*(Θ-u[2])
  end
  g = function (du,u,p,t)
    du[1] = √u[2]*u[1]
    du[2] = Θ*√u[2]
  end
  Γ = [1 ρ;ρ 1] # Covariance Matrix
  noise_rate_prototype = nothing

  if seed == 0
    seed = rand(UInt64)
  end
  noise = CorrelatedWienerProcess!(Γ,tspan[1],zeros(2),zeros(2),
                                   rng = Xorshifts.Xoroshiro128Plus(seed))

  sde_f = SDEFunction{true}(f,g)
  SDEProblem(sde_f,g,u0,tspan,noise=noise,seed=seed,kwargs...)
end

@doc doc"""

``d \ln S(t) = (r(t) - q(t) - \frac{Θ(t,S)^2}{2})dt + σ dW_t``

Solves for ``log S(t)``.
"""
function GeneralizedBlackScholesProblem(r,q,Θ,σ,u0,tspan;kwargs...)
  f = function (u,p,t)
    r(t) - q(t) - Θ(t,exp(u))^2 / 2
  end
  g = function (u,p,t)
    σ
  end
  SDEProblem{false}(f,g,u0,tspan;kwargs...)
end

@doc doc"""

``d \ln S(t) = (r(t) - \frac{Θ(t,S)^2}{2})dt + σ dW_t``

Solves for ``log S(t)``.
"""
mutable struct BlackScholesProblem{uType,tType,tupType,isinplace,NP,F,F2,C,ND,MM} <: DiffEqBase.AbstractSDEProblem{uType,tType,isinplace,ND}
  r::tType
  Θ::tType
  σ::tType
  u0::uType
  tspan::tupType
  p::Nothing
  f::F
  g::F2
  noise::NP
  callback::C
  noise_rate_prototype::ND
  mass_matrix::MM
  seed::UInt64
end

BlackScholesProblem(r,Θ,σ,u0,tspan;callback = CallbackSet(),
                    noise_rate_prototype = nothing,seed=UInt64(0)) =
                    GeneralizedBlackScholesProblem(r,(t)->0,Θ,σ,u0,tspan,callback=callback,seed=seed)

@doc doc"""

``dx = a(b(t)-x)dt + σ dW_t``

"""
function ExtendedOrnsteinUhlenbeckProblem(a,b,σ,u0,tspan;kwargs...)
  f = function (u,p,t)
    a*(b(t)-u)
  end
  g = function (u,p,t)
    σ
  end
  SDEProblem{false}(f,g,u0,tspan;kwargs...)
end

@doc doc"""

``dx = a(r-x)dt + σ dW_t``

"""
function OrnsteinUhlenbeckProblem(a,r,σ,u0,tspan;kwargs...)
  f = function (u,p,t)
    a*(r-u)
  end
  g = function (u,p,t)
    σ
  end
  SDEProblem{false}(f,g,u0,tspan;kwargs...)
end

@doc doc"""

``dx = μ dt + σ dW_t``

"""
function GeometricBrownianMotionProblem(μ,σ,u0,tspan;kwargs...)
  f = function (u,p,t)
    μ
  end
  g = function (u,p,t)
    σ
  end
  SDEProblem{false}(f,g,u0,tspan;kwargs...)
end

@doc doc"""

``dx = σ(t)e^{at} dW_t``

"""
function MfStateProblem(a,σ,u0,tspan;kwargs...)
  f = function (u,p,t)
    0
  end
  g = function (u,p,t)
    σ(t)*exp(a*t)
  end
  SDEProblem{false}(f,g,u0,tspan;kwargs...)
end

ext: jl

max_stars_repo_name: JuliaDiffEq/FinancialModels

max_issues_repo_name: JuliaDiffEq/FinancialModels

max_forks_repo_name: JuliaDiffEq/FinancialModels


text:
#-------------------------------------------------------------------------------
# Types representing the ply data model

"""
    plyname(data)

Return the name that `data` is associated with when serialized in a ply file
"""
function plyname
end

const PropNameList = Union{AbstractVector,Tuple}

#--------------------------------------------------
"""
    ArrayProperty(name, T)

A ply `property \$T \$name`, modelled as an abstract vector, with a name which
can be retrieved using `plyname()`.
"""
mutable struct ArrayProperty{T,Name} <: AbstractVector{T}
    name::Name
    data::Vector{T}
end

#=
# FIXME: Ambiguous constructor
function ArrayProperty(names::PropNameList, data::AbstractVector{T}) where {T}
    if length(names) != length(T)
        error("Number of property names in $names does not match length($T)")
    end
    ArrayProperty(names, data)
end
=#

ArrayProperty(name::AbstractString, ::Type{T}) where {T} = ArrayProperty(String(name), Vector{T}())

Base.summary(prop::ArrayProperty) = "$(length(prop))-element $(typeof(prop)) \"$(plyname(prop))\""

# AbstractArray methods
Base.size(prop::ArrayProperty) = size(prop.data)
Base.getindex(prop::ArrayProperty, i::Int) = prop.data[i]
Base.setindex!(prop::ArrayProperty, v, i::Int) = prop.data[i] = v
Base.IndexStyle(::Type{<:ArrayProperty}) = IndexLinear()

# List methods
Base.resize!(prop::ArrayProperty, len) = resize!(prop.data, len)
Base.push!(prop::ArrayProperty, val) = (push!(prop.data, val); prop)

# Ply methods
plyname(prop::ArrayProperty) = prop.name

#--------------------------------------------------
"""
    ListProperty(name, S, T)
    ListProperty(name, list_of_vectors)

A ply `property list \$S \$T \$name`, modelled as a abstract vector of vectors,
with a name which can be retrieved using `plyname()`.
"""
mutable struct ListProperty{S,T} <: AbstractVector{Vector{T}}
    name::String
    start_inds::Vector{Int}
    data::Vector{T}
end

ListProperty(name, ::Type{S}, ::Type{T}) where {S,T} = ListProperty{S,T}(String(name), ones(Int,1), Vector{T}())

function ListProperty(name::AbstractString, a::AbstractVector)
    # Construct list from an array of arrays
    prop = ListProperty(name, Int32, eltype(a[1]))
    foreach(ai->push!(prop,ai), a)
    prop
end

Base.summary(prop::ListProperty) = "$(length(prop))-element $(typeof(prop)) \"$(plyname(prop))\""

# AbstractArray methods
Base.length(prop::ListProperty) = length(prop.start_inds)-1
Base.size(prop::ListProperty) = (length(prop),)
Base.getindex(prop::ListProperty, i::Int) = prop.data[prop.start_inds[i]:prop.start_inds[i+1]-1]
Base.IndexStyle(::Type{<:ListProperty}) = IndexLinear()
# TODO: Do we need Base.setindex!() ?  Hard to provide with above formulation...

# List methods
function Base.resize!(prop::ListProperty, len)
    resize!(prop.start_inds, len+1)
    prop.start_inds[1] = 1
end
function Base.push!(prop::ListProperty, list)
    push!(prop.start_inds, prop.start_inds[end]+length(list))
    append!(prop.data, list)
    prop
end

# Ply methods
plyname(prop::ListProperty) = prop.name

#--------------------------------------------------
"""
    PlyElement(name, [len | props...])

Construct a ply `element \$name \$len`, containing a list of properties with a
name which can be retrieved using `plyname`.  Properties can be accessed with
the array interface, or looked up by indexing with a string.

The expected length `len` is used if it is set, otherwise the length shared by
the property vectors is used.
"""
mutable struct PlyElement
    name::String
    prior_len::Int  # Length as expected, or as read from file
    properties::Vector
end

PlyElement(name::AbstractString, len::Int=-1) = PlyElement(name, len, Vector{Any}())

function PlyElement(name::AbstractString, props::AbstractVector...)
    PlyElement(name, -1, collect(props))
end

function Base.show(io::IO, elem::PlyElement)
    prop_names = join(["\"$(plyname(prop))\"" for prop in elem.properties], ", ")
    print(io, "PlyElement \"$(plyname(elem))\" of length $(length(elem)) with properties [$prop_names]")
end

# Table-like methods
function Base.length(elem::PlyElement)
    # Check that lengths are consistent and return the length
    if elem.prior_len != -1
        return elem.prior_len
    end
    if isempty(elem.properties)
        return 0
    end
    len = length(elem.properties[1])
    if any(prop->len != length(prop), elem.properties)
        proplens = [length(p) for p in elem.properties]
        throw(ErrorException("Element $(plyname(elem)) has inconsistent property lengths: $proplens"))
    end
    return len
end

function Base.getindex(element::PlyElement, prop_name)
    # Get first property with a matching name
    for prop in element.properties
        if plyname(prop) == prop_name
            return prop
        end
    end
    error("property $prop_name not found in Ply element $(plyname(element))")
end

# List methods
Base.iterate(elem::PlyElement, s...) = iterate(elem.properties, s...)

# Ply methods
plyname(elem::PlyElement) = elem.name

#--------------------------------------------------
"""
    PlyComment(string; [obj_info=false])

A ply comment.

Nonstandard [obj_info header lines](
http://docs.pointclouds.org/1.5.1/structpcl_1_1io_1_1ply_1_1obj__info.html)
may be represented by setting obj_info flag.
"""
struct PlyComment
    comment::String
    obj_info::Bool # Set for comment-like "obj_info" lines
    location::Int # index of previous element (TODO: move this out of the comment)
end

PlyComment(string::AbstractString; obj_info::Bool=false) =
    PlyComment(string, obj_info, -1)

function Base.:(==)(a::PlyComment, b::PlyComment)
    a.comment == b.comment         &&
    a.obj_info == b.obj_info &&
    a.location == b.location
end

#--------------------------------------------------
"""
    Ply()

Container for the contents of a ply file.  This type directly models the
contents of the header.  Ply elements and comments can be added using
`push!()`, elements can be iterated over with the standard iterator
interface, and looked up by indexing with a string.
"""
mutable struct Ply
    elements::Vector{PlyElement}
    comments::Vector{PlyComment}
end

Ply() = Ply(Vector{PlyElement}(), Vector{String}())

function Base.show(io::IO, ply::Ply)
    buf = IOBuffer()
    write_header(ply, buf, true)
    headerstr = String(take!(buf))
    headerstr = replace(strip(headerstr), "\n"=>"\n ")
    print(io, "$Ply with header:\n $headerstr")
end

# List methods
Base.push!(ply::Ply, el::PlyElement) = (push!(ply.elements, el); ply)
function Base.push!(ply::Ply, c::PlyComment)
    push!(ply.comments, PlyComment(c.comment, c.obj_info, length(ply.elements)+1))
    ply
end

# Element search and iteration
function Base.getindex(ply::Ply, elem_name::AbstractString)
    for elem in ply.elements
        if plyname(elem) == elem_name
            return elem
        end
    end
    error("$elem_name not found in Ply element list")
end

Base.length(ply::Ply) = length(ply.elements)
Base.iterate(ply::Ply, s...) = iterate(ply.elements, s...)


text:
# turn a openAI finetune-transformer-lm npy to bson

using JSON
using ZipFile

using Flux: loadparams!

using BytePairEncoding

"""
    openAInpy2bson(path;
                   raw=false,
                   saveto="./",
                   startsym="_start_",
                   delisym="_delimiter_",
                   clfsym="_classify_",
                   unksym="<unk>")

turn openai released gpt format(npy) into BSON file. Set `raw` to `true` to remain the origin data format in bson.
"""
function openAInpy2bson(path;
                        raw=false,
                        saveto="./",
                        startsym="_start_",
                        delisym="_delimiter_",
                        clfsym="_classify_",
                        unksym="<unk>")
  if iszip(path)
    filename = basename(path)[1:end-4]
    data = ZipFile.Reader(path)
  else
    filename = basename(isdirpath(path) ? path[1:end-1] : path)
    data = path
  end

  tokenizer = gpt_tokenizer
  weights, bpe, raw_vocab = readnpzfolder(data)

  iszip(path) && close(data)

  if raw
    bsonname = normpath(joinpath(saveto, filename * ".npbson"))
    BSON.@save bsonname weights bpe raw_vocab tokenizer
  else
    bsonname = normpath(joinpath(saveto, filename * ".bson"))
    vocab = build_vocab(raw_vocab;
                        startsym=startsym,
                        delisym=delisym,
                        clfsym=clfsym,
                        unksym=unksym,
                        )
    gpt_model = load_gpt_from_npbson(weights, length(vocab))
    BSON.@save bsonname gpt_model bpe vocab tokenizer
  end
end

readnpz(path) = error("readnpz require NPZ.jl installed. run `Pkg.add(\"NPZ\"); using NPZ`")

@init @require NPZ="15e1cf62-19b3-5cfa-8e77-841668bca605" begin
  import .NPZ

  function readnpz(dir)
    shapes = JSON.parsefile(joinpath(dir, "params_shapes.json"))
    offsets = accumulate(+, prod.(shapes))
    shapes = map(s -> length(s) > 1 ? (s[end], s[end-1]) : s, shapes)
    params = cat([npzread(joinpath(dir, "params_$(i).npy")) for i = 0:9]..., dims=1)
    params = [collect(reshape(selectdim(params, 1, a+1:b), s...)) for (a, b, s) in zip([0;offsets[1:end-1]], offsets, shapes)]
    params
  end

  function readnpz(z::ZipFile.Reader)
    shapes = JSON.parse(findfile(z, "params_shapes.json"))
    offsets = accumulate(+, prod.(shapes))
    shapes = map(s -> length(s) > 1 ? (s[end], s[end-1]) : s, shapes)
    params = cat([NPZ.npzreadarray(findfile(z, "params_$(i).npy")) for i = 0:9]..., dims=1)
    params = [collect(reshape(selectdim(params, 1, a+1:b), s...)) for (a, b, s) in zip([0;offsets[1:end-1]], offsets, shapes)]
    params
  end
end

function readnpzfolder(dir)
  emp = JSON.parsefile(joinpath(dir, "encoder_bpe_40000.json"))
  vocab = map(first, sort!(collect(emp), by=(x)->x.second))
  bpe = Bpe(joinpath(dir, "vocab_40000.bpe"))
  weights = readnpz(dir)

  weights, bpe, vocab
end

function readnpzfolder(z::ZipFile.Reader)
  emp = JSON.parse(findfile(z, "encoder_bpe_40000.json"))
  vocab = map(first, sort!(collect(emp), by=(x)->x.second))
  bpe = mktemp(
    (fn, f)-> begin
      zf = findfile(z, "vocab_40000.bpe")
      buffer = Vector{UInt8}(undef, zf.uncompressedsize)
      write(f, read!(zf, buffer))
      close(f)
      Bpe(fn)
    end
  )
  weights = readnpz(z)
  weights, bpe, vocab
end

function build_vocab(raw_vocab;
                     startsym="_start_",
                     delisym="_delimiter_",
                     clfsym="_classify_",
                     unksym="<unk>")
    vocab = copy(raw_vocab)
    push!(vocab, startsym)
    push!(vocab, delisym)
    push!(vocab, clfsym)

    vocab = Vocabulary(vocab, unksym)
end

load_gpt_from_npbson(path::AbstractString) = (@assert isnpbson(path); load_gpt_from_npbson(BSON.load(path)))
load_gpt_from_npbson(bson) = load_gpt_from_npbson(bson[:weights], length(bson[:raw_vocab])+3)
function load_gpt_from_npbson(weights, vocab_size; n_special = 3)
    embed = Embed(768, vocab_size)
    pe = PositionEmbedding(768, 512; trainable=true)

    gpt = Gpt(768, 12, 768*4, 12; act=gelu, pdrop=0.1, attn_pdrop=0.1)

    # pm_name = (:pe, :embed, (
    #     (:mh.iqproj.W, :mh.ikproj.W, :mh.ivproj.W),
    #     (:mh.iqproj.b, :mh.ikproj.b, :mh.ivproj.b),
    #     :mh.oproj.W, :mh.oproj.b,
    #     :LN1.α, :LN1.β,
    #     :pw.din.W, :pw.din.b,
    #     :pw.dout.W, :pw.dout.b,
    #     :LN2.α, :LN2.β),
    #            :xN,
    #            )

    loadparams!(embed, [hcat(weights[2],
                             randn(768, n_special) .* 0.02)])
    loadparams!(pe, [weights[1]])
    for i = 1:12
        mhW = weights[12(i-1) + 3]
        mhb = weights[12(i-1) + 4]
        loadparams!(gpt.ts[i].mh.iqproj,[selectdim(mhW, 1, 1:768),
                                         selectdim(mhb, 1, 1:768)])
        loadparams!(gpt.ts[i].mh.ikproj,[selectdim(mhW, 1, 768+1:2*768),
                                         selectdim(mhb, 1, 768+1:2*768)])
        loadparams!(gpt.ts[i].mh.ivproj,[selectdim(mhW, 1, 2*768+1:3*768),
                                         selectdim(mhb, 1, 2*768+1:3*768)])
        loadparams!(gpt.ts[i].mh.oproj,[weights[12(i-1) + 5],
                                        weights[12(i-1) + 6]])
        loadparams!(gpt.ts[i].mhn,[weights[12(i-1) + 7],
                                   weights[12(i-1) + 8]])
        loadparams!(gpt.ts[i].pw.din,[weights[12(i-1) + 9],
                                      weights[12(i-1) + 10]])
        loadparams!(gpt.ts[i].pw.dout,[weights[12(i-1) + 11],
                                       weights[12(i-1) + 12]])
        loadparams!(gpt.ts[i].pwn,[weights[12(i-1) + 13],
                                   weights[12(i-1) + 14]])
    end

    ce = CompositeEmbedding(tok=embed, pe=pe)
    TransformerModel(ce, gpt)
end


text:
using DrillMudsThermalProps
using Documenter

DocMeta.setdocmeta!(DrillMudsThermalProps, :DocTestSetup, :(using DrillMudsThermalProps); recursive=true)

makedocs(;
    modules=[DrillMudsThermalProps],
    authors="Eduardo Alves",
    repo="https://github.com/Eduardo-BDMAlves/DrillMudsThermalProps.jl/blob/{commit}{path}#{line}",
    sitename="DrillMudsThermalProps.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Eduardo-BDMAlves.github.io/DrillMudsThermalProps.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Eduardo-BDMAlves/DrillMudsThermalProps.jl",
    devbranch="Development",
)


text:
using ACTRSimulators, Test, ACTRModels, Random, DataFrames
import ACTRSimulators: start!, press_key!
Random.seed!(8985)
include("task.jl")

scheduler = ACTRScheduler(;model_trace=true, store=true)
task = SimpleTask(;scheduler)
procedural = Procedural()
T = vo_to_chunk() |> typeof
memory = [Chunk(;animal=:dog), Chunk(;animal=:cat)]
declarative = Declarative(;memory)
visicon = VisualObject[]
visual_location = VisualLocation(buffer=T[])
visual = Visual(buffer=T[])

actr = ACTR(;scheduler, procedural, visual_location, visual, declarative, visicon)

function can_attend(actr)
    c1(actr) = !isempty(actr.visual_location.buffer)
    c2(actr) = !actr.visual.state.busy
    return c1, c2
end

function can_stop(actr)
    c1(actr) = !actr.visual.state.empty
    return (c1,)
end

function attend_action(actr, task)
    buffer = actr.visual_location.buffer
    chunk = deepcopy(buffer[1])
    clear_buffer!(actr.visual_location)
    attending!(actr, chunk)
    return nothing
end

function stop(actr, task)
    stop!(actr.scheduler)
end

rule1 = Rule(;conditions=can_attend, action=attend_action, actr, task, name="Attend")
push!(procedural.rules, rule1)

rule2 = Rule(;conditions=can_stop, action=stop, actr, task, name="Stop")
push!(procedural.rules, rule2)
run!(actr, task)
chunk = actr.visual.buffer[1]
@test chunk.slots == (color = :black, text = "hello", x = 300.0, y = 300.0)

observed = map(x->x.description, scheduler.complete_events)
expected = [
    "Starting",
    "Present Stimulus",
    "Selected Attend",
    "Attend",
    "Selected Stop"
]
@test expected == observed


text:
include("BoundaryConditionsCylindrical.jl")
include("BoundaryConditionsCartesian.jl")


text:
###############################################################################
#
#   fmpz_mpoly.jl : Flint multivariate polynomials over fmpz
#
###############################################################################

export FmpzMPolyRing, fmpz_mpoly

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

parent_type{S, N}(::Type{fmpz_mpoly{S, N}}) = FmpzMPolyRing{S, N}

elem_type{S, N}(::FmpzMPolyRing{S, N}) = fmpz_mpoly{S, N}

vars(a::FmpzMPolyRing) = a.S

function gens{S, N}(R::FmpzMPolyRing{S, N})
   A = Array{fmpz_mpoly{S, N}}(R.num_vars)
   for i = 1:R.num_vars
      z = R()
      ccall((:fmpz_mpoly_gen, :libflint), Void,
            (Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}), z, i - 1, R)
      A[i] = z
   end
   return A
end

function gen{S, N}(R::FmpzMPolyRing{S, N}, i::Int)
   z = R()
   ccall((:fmpz_mpoly_gen, :libflint), Void,
            (Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}), z, i - 1, R)
   return z
end

function isgen(a::fmpz_mpoly)
   R = parent(a)
   return Bool(ccall((:fmpz_mpoly_is_gen, :libflint), Cint,
              (Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}), a, a.parent))
end

function isgen(a::fmpz_mpoly, i::Int)
   R = parent(a)
   return Bool(ccall((:fmpz_mpoly_is_gen_i, :libflint), Cint,
              (Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}), a, i - 1, a.parent))
end

function coeff(a::fmpz_mpoly, i::Int)
   z = fmpz()
   ccall((:fmpz_mpoly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
         z, a, i, a.parent)
   return z
end

function deepcopy(a::fmpz_mpoly)
   z = parent(a)()
   ccall((:fmpz_mpoly_set, :libflint), Void,
         (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
         z, a, a.parent)
   return z
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

num_vars(x::fmpz_mpoly) = parent(x).num_vars

function max_degrees{S, N}(a::fmpz_mpoly{S, N})
   R = a.parent
   A = Array(Int, N)
   ccall((:fmpz_mpoly_max_degrees, :libflint), Void,
         (Ref{Int}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}), A, a, R)
   biggest = 0
   for i = 1:N
      if A[i] > biggest
         biggest = A[i]
      end
  end
  return biggest, A
end

###############################################################################
#
#   AbstractString{} I/O
#
###############################################################################

function show(io::IO, x::fmpz_mpoly)
   if length(x) == 0
      print(io, "0")
   else
      cstr = ccall((:fmpz_mpoly_get_str_pretty, :libflint), Ref{UInt8},
          (Ref{fmpz_mpoly}, Ptr{Ref{UInt8}}, Ref{FmpzMPolyRing}),
          x, [string(s) for s in vars(parent(x))], x.parent)
      print(io, unsafe_string(cstr))

      ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
   end
end

function show(io::IO, p::FmpzMPolyRing)
   const max_vars = 5 # largest number of variables to print
   n = p.num_vars
   print(io, "Multivariate Polynomial Ring in ")
   if n > max_vars
      print(io, p.num_vars)
      print(io, " variables ")
   end
   for i = 1:min(n - 1, max_vars - 1)
      print(io, string(p.S[i]), ", ")
   end
   if n > max_vars
      print(io, "..., ")
   end
   print(io, string(p.S[n]))
   print(io, " over ")
   show(io, base_ring(p))
end

###############################################################################
#
#   Basic arithmetic
#
###############################################################################

function -(a::fmpz_mpoly)
   z = parent(a)()
   ccall((:fmpz_mpoly_neg, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, a.parent)
   return z
end

function +{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   ccall((:fmpz_mpoly_add, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent)
   return z
end

function -{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   ccall((:fmpz_mpoly_sub, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent)
   return z
end

function *{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   ccall((:fmpz_mpoly_mul_johnson, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent)
   return z
end

function mul_array{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   ccall((:fmpz_mpoly_mul_array, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent)
   return z
end

###############################################################################
#
#   Ad hoc arithmetic
#
###############################################################################

function *(a::fmpz_mpoly, b::Int)
   r = parent(a)()
   ccall((:fmpz_mpoly_scalar_mul_si, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

function *(a::fmpz_mpoly, b::fmpz)
   r = parent(a)()
   ccall((:fmpz_mpoly_scalar_mul_fmpz, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz}, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

*(a::fmpz, b::fmpz_mpoly) = b*a

*(a::Int, b::fmpz_mpoly) = b*a

function +(a::fmpz_mpoly, b::Int)
   r = parent(a)()
   ccall((:fmpz_mpoly_add_si, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

function +(a::fmpz_mpoly, b::fmpz)
   r = parent(a)()
   ccall((:fmpz_mpoly_add_fmpz, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz}, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

+(a::fmpz, b::fmpz_mpoly) = b + a

+(a::Int, b::fmpz_mpoly) = b + a

function -(a::fmpz_mpoly, b::Int)
   r = parent(a)()
   ccall((:fmpz_mpoly_sub_si, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

function -(a::fmpz_mpoly, b::fmpz)
   r = parent(a)()
   ccall((:fmpz_mpoly_sub_fmpz, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz}, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

-(a::fmpz, b::fmpz_mpoly) = -(b - a)

-(a::Int, b::fmpz_mpoly) = -(b - a)

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpz_mpoly, b::Int)
   b < 0 && throw(DomainError())
   z = parent(a)()
   ccall((:fmpz_mpoly_pow_fps, :libflint), Void,
         (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
         z, a, b, parent(a))
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function =={S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   return Bool(ccall((:fmpz_mpoly_equal, :libflint), Cint,
               (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
               a, b, a.parent))
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

function ==(a::fmpz_mpoly, b::Int)
   return Bool(ccall((:fmpz_mpoly_equal_si, :libflint), Cint,
               (Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
               a, b, a.parent))
end

==(a::Int, b::fmpz_mpoly) = b == a

function ==(a::fmpz_mpoly, b::fmpz)
   return Bool(ccall((:fmpz_mpoly_equal_fmpz, :libflint), Cint,
               (Ref{fmpz_mpoly}, Ref{fmpz}, Ref{FmpzMPolyRing}),
               a, b, a.parent))
end

==(a::fmpz, b::fmpz_mpoly) = b == a

==(a::fmpz_mpoly, b::Integer) = a == fmpz(b)

==(a::Integer, b::fmpz_mpoly) = b == a

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divides_array{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   d = ccall((:fmpz_mpoly_divides_array, :libflint), Cint,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent)
   d == -1 && error("Polynomial too large for divides_array")
   return d == 1, z
end

function divides_monagan_pearce{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   z = parent(a)()
   d = Bool(ccall((:fmpz_mpoly_divides_monagan_pearce, :libflint), Cint,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       z, a, b, a.parent))
   return d, z
end

###############################################################################
#
#   Division with remainder
#
###############################################################################

function div_monagan_pearce{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   q = parent(a)()
   ccall((:fmpz_mpoly_div_monagan_pearce, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
        Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       q, a, b, a.parent)
   return q
end

function divrem_monagan_pearce{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   q = parent(a)()
   r = parent(a)()
   ccall((:fmpz_mpoly_divrem_monagan_pearce, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
        Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       q, r, a, b, a.parent)
   return q, r
end

function divrem_array{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   q = parent(a)()
   r = parent(a)()
   d = ccall((:fmpz_mpoly_divrem_array, :libflint), Cint,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
        Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}),
       q, r, a, b, a.parent)
   d == 0 && error("Polynomials too large for divrem_array")
   return q, r
end

function divrem_monagan_pearce{S, N}(a::fmpz_mpoly{S, N}, b::Array{fmpz_mpoly{S, N}, 1})
   len = length(b)
   q = [parent(a)() for i in 1:len]
   r = parent(a)()
   ccall((:fmpz_mpoly_divrem_ideal, :libflint), Void,
       (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
        Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
       q, r, a, b, len, a.parent)
   return q, r
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(a::fmpz_mpoly, b::Int)
   r = parent(a)()
   ccall((:fmpz_mpoly_scalar_divexact_si, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Int, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

function divexact(a::fmpz_mpoly, b::fmpz)
   r = parent(a)()
   ccall((:fmpz_mpoly_scalar_divexact_fmpz, :libflint), Void,
        (Ref{fmpz_mpoly}, Ref{fmpz_mpoly}, Ref{fmpz}, Ref{FmpzMPolyRing}),
        r, a, b, a.parent)
   return r
end

divexact(a::fmpz_mpoly, b::Integer) = divexact(a, fmpz(b))

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function addeq!{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N})
   ccall((:fmpz_mpoly_add, :libflint), Void,
         (Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
          Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}), a, a, b, a.parent)
   return a
end

function mul!{S, N}(a::fmpz_mpoly{S, N}, b::fmpz_mpoly{S, N}, c::fmpz_mpoly{S, N})
   ccall((:fmpz_mpoly_mul_johnson, :libflint), Void,
         (Ref{fmpz_mpoly}, Ref{fmpz_mpoly},
          Ref{fmpz_mpoly}, Ref{FmpzMPolyRing}), a, b, c, a.parent)
   return a
end

function setcoeff!(a::fmpz_mpoly, i::Int, c::Int)
   ccall((:fmpz_mpoly_set_coeff_si, :libflint), Void,
        (Ref{fmpz_mpoly}, Int, Int, Ref{FmpzMPolyRing}), a, i, c, a.parent)
   return a
end

function setcoeff!(a::fmpz_mpoly, i::Int, c::fmpz)
   ccall((:fmpz_mpoly_set_coeff_fmpz, :libflint), Void,
        (Ref{fmpz_mpoly}, Int, Ref{fmpz}, Ref{FmpzMPolyRing}),
                                                          a, i, c, a.parent)
   return a
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule{V <: Integer}(::Type{fmpz_mpoly}, ::Type{V}) = fmpz_mpoly

promote_rule(::Type{fmpz_mpoly}, ::Type{fmpz}) = fmpz_mpoly

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (R::FmpzMPolyRing{S, N}){S, N}()
   z = fmpz_mpoly(R)
   z.parent = R
   return z
end

function (R::FmpzMPolyRing{S, N}){S, N}(b::Int)
   z = fmpz_mpoly(R, b)
   z.parent = R
   return z
end

function (R::FmpzMPolyRing{S, N}){S, N}(b::fmpz)
   z = fmpz_mpoly(R, b)
   z.parent = R
   return z
end

function (R::FmpzMPolyRing{S, N}){S, N}(b::Integer)
   z = fmpz_mpoly(R, fmpz(b))
   z.parent = R
   return z
end

function (R::FmpzMPolyRing{S, N}){S, N}(a::Array{fmpz, 1}, b::Array{NTuple{N, Int}, 1})
   z = fmpz_mpoly(R, a, b)
   z.parent = R
   return z
end

###############################################################################
#
#   PolynomialRing constructor
#
###############################################################################

function PolynomialRing(R::FlintIntegerRing, s::Array{String, 1}; cached::Bool = true, ordering::Symbol = :lex)
   U = [Symbol(x) for x in s]
   N = (ordering == :deglex || ordering == :degrevlex) ? length(U) + 1 : length(U)
   # default to 8 bit exponent fields
   parent_obj = FmpzMPolyRing{ordering, N}(U, cached)
   return tuple(parent_obj, gens(parent_obj))
end


text:
using ChainRulesCore
struct DiffractorRuleConfig <: RuleConfig{Union{HasReverseMode,HasForwardsMode}} end

@Base.constprop :aggressive accum(a, b) = a + b
@Base.constprop :aggressive accum(a::Tuple, b::Tuple) = map(accum, a, b)
@Base.constprop :aggressive @generated function accum(x::NamedTuple, y::NamedTuple)
    fnames = union(fieldnames(x), fieldnames(y))
    gradx(f) = f in fieldnames(x) ? :(getfield(x, $(quot(f)))) : :(ZeroTangent())
    grady(f) = f in fieldnames(y) ? :(getfield(y, $(quot(f)))) : :(ZeroTangent())
    Expr(:tuple, [:($f=accum($(gradx(f)), $(grady(f)))) for f in fnames]...)
end
@Base.constprop :aggressive accum(a, b, c, args...) = accum(accum(a, b), c, args...)
@Base.constprop :aggressive accum(a::NoTangent, b) = b
@Base.constprop :aggressive accum(a, b::NoTangent) = a
@Base.constprop :aggressive accum(a::NoTangent, b::NoTangent) = NoTangent()


text:
"""
    xlabel!([scene,] xlabel)

Set the x-axis label for the given Scene.
Defaults to using the current Scene.
"""
function xlabel!(scene, xlabel::AbstractString)
    axis = scene[Axis]
    @assert !isnothing(axis)
    scene[Axis][:names][:axisnames][] = (xlabel, scene[Axis][:names][:axisnames][][2:end]...)
    nothing
end
xlabel!(xlabel::AbstractString) = xlabel!(current_scene(), xlabel)

"""
    ylabel!([scene,] ylabel)

Set the y-axis label for the given Scene.
Defaults to using the current Scene.
"""
function ylabel!(scene, ylabel::AbstractString)
    axis = scene[Axis]
    @assert !isnothing(axis)
    if axis isa Axis2D
        scene[Axis][:names][:axisnames][] = (scene[Axis][:names][:axisnames][][1], ylabel)
    elseif axis isa Axis3D
        scene[Axis][:names][:axisnames][] = (scene[Axis][:names][:axisnames][][1], ylabel, scene[Axis][:names][:axisnames][][3])
    else
        @error("Unknown axis type $(typeof(axis)).")
    end
    nothing
end
ylabel!(ylabel::AbstractString) = xlabel!(current_scene(), ylabel)

"""
    zlabel!([scene,] zlabel)

Set the z-axis label for the given Scene.
Defaults to using the current Scene.
!!! warning
    The Scene must have an Axis3D.  If not, then this function will error.
"""
function zlabel!(scene, zlabel::AbstractString)
    axis = scene[Axis]
    @assert !isnothing(axis)
    @assert axis isa Axis3D
    scene[Axis][:names][:axisnames][] = (scene[Axis][:names][:axisnames][][1], scene[Axis][:names][:axisnames][][2], zlabel)
    nothing
end
zlabel!(zlabel::AbstractString) = zlabel!(current_scene(), zlabel)

################################################################################

function setlims!(scene::Scene, lims::NTuple{2, Real}, dim=1)
    ol = scene.limits[]                          # get the Scene's limits as values
    o_origin = ol.origin                         # get the original origin
    o_widths = ol.widths                         # get the original widths
    n_widths = convert(Vector, o_widths)         # convert to mutable form
    n_origin = convert(Vector, o_origin)         # convert to mutable form
    n_origin[dim] = lims[1]                      # set the new origin in dim
    n_widths[dim] = lims[2] - lims[1]            # set the new width in dim
    scene.limits[] = FRect3D(n_origin, n_widths) # set the limits of the scene
    nothing
end

"""
    xlims!(limits::Real...)
    xlims!(limits::NTuple{2, Real})
    xlims!(scene, limits::Real...)
    xlims!(scene, limits::NTuple{2, Real})

Set the x-limits for the given Scene (defaults to current Scene).
"""
xlims!(scene::Scene, lims::NTuple{2, Real}) = setlims!(scene, lims, 1)

"""
    ylims!(limits::Real...)
    ylims!(limits::NTuple{2, Real})
    ylims!(scene, limits::Real...)
    ylims!(scene, limits::NTuple{2, Real})

Set the y-limits for the given Scene (defaults to current Scene).
"""
ylims!(scene::Scene, lims::NTuple{2, Real}) = setlims!(scene, lims, 2)

"""
    zlims!(limits::Real...)
    zlims!(limits::NTuple{2, Real})
    zlims!(scene, limits::Real...)
    zlims!(scene, limits::NTuple{2, Real})

Set the z-limits for the given Scene (defaults to current Scene).
"""
zlims!(scene::Scene, lims::NTuple{2, Real}) = setlims!(scene, lims, 3)

xlims!(scene::Scene, lims::Real...) = xlims!(current_scene(), lims)
ylims!(scene::Scene, lims::Real...) = ylims!(current_scene(), lims)
zlims!(scene::Scene, lims::Real...) = zlims!(current_scene(), lims)

xlims!(lims::NTuple{2, Real}) = xlims!(current_scene(), lims)
ylims!(lims::NTuple{2, Real}) = ylims!(current_scene(), lims)
zlims!(lims::NTuple{2, Real}) = zlims!(current_scene(), lims)

xlims!(lims::Real...) = xlims!(current_scene(), lims)
ylims!(lims::Real...) = ylims!(current_scene(), lims)
zlims!(lims::Real...) = zlims!(current_scene(), lims)


text:
#= DM2_1 avant modification
function GRASP(
    cost, #Une array représentant les couts, de taille n
    matrix, #une matrice de taille m*n representant les contrantes
    n, #nombre de variables
    m, #nombbre de contraintes
    alpha # le alpha que nous allons utiliser
     )

     #intitialisation des listes utilisé
     (desactive_condition, stop1, variables_actives, stop2, util, SOL) = initaliser(m,n)

    #Création de la solution
    while desactive_condition!=stop1 && variables_actives!=stop2

        util = Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2) #Foction d'utilité

        PosCandidat::Int64 = choixcandidat(util, alpha) #Choix parmis les candidats de util

        SOL[PosCandidat] = true #On ajoute le candidat à la solution

        Desactive!(PosCandidat, matrix, desactive_condition, m, variables_actives,n) #Desactive! le candidat selectionné et les conditions où il apparait
    end

    Z = calculz(SOL,cost,n)
    return(SOL, Z, desactive_condition)

end

# =========================================================================== #
function initaliser(m,n)
    #Creating a set of lines that will be avaluated
    desactive_condition= zeros(Bool, m) #si desactive_condition[j]=0 la ligne sera évalué
    stop1= ones(Bool,m) #la une condition d'arret, s'active quand toutes les lignes ont été traité

    #Creation the an array of activated variables
    variables_actives = ones(Bool, n) #Si variables_actives[i]=1, la variable sera traité
    stop2= zeros(Bool , n) #condition d'arret, s'active quand toutes les variables sont traité

    #On crée une liste contenant nos utilités
    util = zeros(Float64, n)

    #Initialisation de la Solution
    SOL = zeros(Bool, n)

    return(desactive_condition, stop1, variables_actives, stop2, util, SOL)
end

 #=================================================================================#

#Détermine l'utilité des éléments basé sur leur nombres d'apparition dans les condition et leur valeur dans les couts
function Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)

    util = zeros(Float64, n) #réinitialisation du vecteur

    for j=1:n #Pour chaque variable
            if variables_actives[j] == 1 #on verifie si la variable est active
                    for i=1:m #Pour chaque condition
                        if desactive_condition[i]==0 #On verifie si la condition doit être évalué
                            K=matrix[i,j]
                            util[j]=util[j]+K
                        end
                    end

                #On divise le le cout de chqaue variable par son nombre d'appartition dans les conditions
                if util[j]!=false #On évide de divisé par 0...
                    util[j] =  cost[j] / util[j]
                else #Si la variable n'a pas été décompté, on la désactive
                    variables_actives[j] = 0
                end
            end
    end
    return util
end
#============================================================#

function Desactive!(PosCandidat, matrix, desactive_condition, m , variables_actives , n) #On désactive les lignes où est le candidat
    for i=1:m
        if matrix[i,PosCandidat] == 1
            desactive_condition[i] = 1
            for j=1:n
                if matrix[i,j] ==1
                    variables_actives[j] = 0
                end
            end
        end
    end
end

#=========================================================================#

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end

#==========================================#

function choixcandidat(util, alpha) #Permet de construire et de choisir un candidat dans la liste de candidat restreint
    max = maximum(util)
    min = minimum(util)
    Candidat_Restreint = Float64[]
    cpt::Int=1
    for i in 1:length(util)
        if util[i]>= min+(alpha*(max-min))
            push!(Candidat_Restreint,i)
        end
    end
    prob = rand(1:length(Candidat_Restreint))
    return(Candidat_Restreint[prob])
end

=#

#DM2_2 avant modification
include("DM1_1.jl")
include("DM2_1.jl")
include("DM1_2_corrige.jl")
#=========================================#
function ReactiveGRASP( #Utilisation du Reactive GRASP
    matrix, #Un matrice de taille n*m qui représente les conditions
    cost, #Une liste des couts de chaque variable
    n, #nombre de variables
    m, #nombre de contraintes
    ite, #on tierera ce nombre avant de recalculer des probabilités de chaque alpha
    coupe, #Le nombre de coupe du segment [0,1], nous permet de construire nos alpha
    temps #La ressource en temps en seconde alloué
    )

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe)

    evol_p=Float64[] #liste qui stokera l'évolution de nos probabilités au fur et à mesure des itérations.
    t=time()
    cpt = 1

    while (time()-t <= temps) #Condtion sur les ressources en temps
        append!(evol_p,p)
        while (cpt <= ite)

            for i in 1:length(p) #On s'assure que chaque probabilité se soit exprimé au moins une fois
                (SOL,z, crts) = GRASP(cost, matrix, n, m, p[i])
                z_cumul[i] += z
            end

            prob=rand(Float64)

            alpha_choisit = choix_alpha(p,prob)
            (SOL, z, crts) = GRASP(cost, matrix, n, m, p[alpha_choisit])

            #Amelioration
            (SOL, z) = exchange1_1(SOL,n,m,cost,crts,matrix)

            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z

            if z > zBest
                zBest = z
            end

            if z < zWorst
                zWorst = z
            end

            cpt = cpt+1 #Compteur du nombre d'itération
        end

        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end
    println("zBest: ", zBest, "    zWorst: ", zWorst)
    println(evol_p) #A commenter si on ne souhaite pas voir la liste d'évolution des probabilités
    return(evol_p)
end

#=======================================#
function intialiser(matrix,cost, n, m, ite, coupe) #initialisationd des variables

    p = Float64[]
    for i in 1:coupe
        push!(p,i/coupe)
    end

    nb_iteration = ones(Int64, length(p))

    z_cumul = zeros(Int64, length(p))

    (SOL,z) =Glouton(cost, matrix, n, m)

    zBest = z

    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end

#=======================================#

function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    q=zeros(Float64,length(p))

    somme_q=0

    for i in 1:length(p)

        moyenne=(z_cumul[i]/nb_iteration[i])

        q[i]=(moyenne-zWorst)/(zBest-zWorst)

        somme_q+=q[i]

    end

    for i in 1:length(q)
        p[i]=q[i]/somme_q
    end

    for i in 2:length(p)-1
        p[i]+=p[i-1]
    end
    p[length(p)]=1

end

#=======================================#

function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end
    end
end

#===#

#=
function exchange1_1(solution,n,m,couts,crts,matrix)
    x_best = Solution(solution,calculz(solution,couts,n))
    # notre boucle k
    for i in 1:n
        if  x_best.x[i]  == 1
            x_prime = deepcopy(x_best)
            x_prime.x[i] = 0
            x_prime.objectif = calculz(x_prime.x,couts,n)
            # notre boucle p
            for j = 1:n
                if x_prime.x[j] == 0 && j != i # j != i car ça sers à rien de remettre une variable enlevée
                    x_seconde = deepcopy(x_prime)
                    x_seconde.x[j] = 1
                    x_seconde.objectif = calculz(x_prime.x,couts,n)

                    #on met à jour notre meilleur Solution
                    if x_seconde.objectif > x_best.objectif
                        x_best = deepcopy(x_seconde)
                    end
                end
            end
        end
    end
    #x,objectif= Solution.x
    return(x_best.x,x_best.objectif)
end
=#

#Reactive grasp final

#=
include("DM1_1.jl")
include("DM2_1.jl")

function ReactiveGRASP(matrix,cost, n, m, ite, coupe,temps)
#    println("^^^^^^^^ Debut R-GRASP ^^^^^^^^")

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe)

#    println("valeurs initiales: ")
#    println("p: ", p)
#    println("nb_iteration: ", nb_iteration)
#    println("z_cumul: ", z_cumul)
#    println("zBest: ", zBest , "  zWorst: ", zWorst)
    evol_p=Float64[]
    t=time()
    while (time()-t <= temps)
        cpt = 1
        append!(evol_p,p)
        while (cpt <= ite)
            for i in 1:length(p)
                (SOL,z) = GRASP(cost, matrix, n, m, p[i])
                z_cumul[i] += z
            end

            cpt = cpt+1

            prob=rand(Float64)

            alpha_choisit = choix_alpha(p,prob)
#            println("alpha_choisit: ", alpha_choisit, "  p: ", p)
            (SOL,z) = GRASP(cost, matrix, n, m, p[alpha_choisit])

            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z

            if z > zBest
                zBest = z
            end

            if z < zWorst
                zWorst = z
            end

#            println(" ============= \n ite: ",cpt)
        end
#        println("zBest: ", zBest, "   zWorst: ", zWorst)
#        println("  z_cumul: ",z_cumul)
#        println( "    nb_iteration : ", nb_iteration , "    p: ", p)
        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end
    println("^^^^^^^^ Fin ReactiveGRASP ^^^^^^^^")
    println("zBest: ", zBest, "    zWorst: ", zWorst)
    println(evol_p)

end

#=======================================#

function intialiser(matrix,cost, n, m, ite, coupe)

    p = Float64[]
    for i in 1:coupe
        push!(p,i/coupe)
    end

    nb_iteration = ones(Int64, length(p))

    z_cumul = zeros(Int64, length(p))

    (SOL,z) =Glouton(cost, matrix, n, m)

    zBest = z

    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end

#=======================================#

function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    q=zeros(Float64,length(p))

    somme_q=0
#    println("q: ", q)
#    println("nb_iteration: ", nb_iteration)

    for i in 1:length(p)
#        println("q: ", q, "   i:", i)
#        if nb_iteration[i]!=0
#            println("z_cumul[i]: ", z_cumul[i],"   nb_iteration[i]: ", nb_iteration[i] )
            moyenne=(z_cumul[i]/nb_iteration[i])
#            println("moyenne: ", moyenne)
#        else
#            moyenne=0
#        end
#        if moyenne>0
            q[i]=(moyenne-zWorst)/(zBest-zWorst)
#            println("q[i]: ", q[i])

#        else
#            q[i]= zWorst/(zBest-zWorst)
#        end
        somme_q+=q[i]

    end
#    println("q: ", q, "      somme_q", somme_q)

    for i in 1:length(q)
        p[i]=q[i]/somme_q
    end

    for i in 2:length(p)-1
        p[i]+=p[i-1]
    end
    p[length(p)]=1
#    println("p: ", p)

end

#=======================================#

function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end
    end
end
=#
#=
#DM2_1
# =========================================================================== #

# =========================================================================== #

#Implementing an x0 viable solution huristic with a gluton method
#receiving an array of lenght()=m , a matrix of size m*n and thoose size
#returing a Bool array of lenght=n
#Asked question:
#Mettre en place une heuristique de construction d’une x0 réalisable.

function RandomGlouton(cost, matrix, n, m, alpha)

    #Creating a set of lines that will be avaluated
    actif= zeros(Bool, m) #if actif[j]=0 the line will be avaluated
    stop1= ones(Bool,m)

    #Creation the an array of activated variables
    varactive = ones(Bool, n)
    stop2= zeros(Bool , n)
#    println("\n varactive: ", varactive, "\n stop2: ", stop2)

    #creating my utility array for my m variables
    util = zeros(Float64, n)

    #initialisation de la Solution
    SOL = zeros(Float64, n)

#    println( "\n Matrix: ", matrix, "\n Cost: ", cost,"\n n= ", n, "   m= ", m,"\n actif: ", actif, "\n stop: ", stop1, "\n util: ", util, "\n SOL: ", SOL, "\n \n")

    ite = 0
    #Création de la solution
    while actif!=stop1 && varactive!=stop2

#      println("\n ++++++++++++++++++++++++++++++++++++++++")
#      println("\n itération: ", ite)
      #println("\n actif : ", actif, "\n stop1 : " , "\n varactive : " , varactive , "\n stop2 : ")
#      println("\n ++++++++++++++++++++++++++++++++++++++++")

        #Foction d'utilité

        util = Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)

        #Choix parmis les candidats de util
#        PosCandidat::Int64
#        PosCandidat = PosMax(util,n)
        PosCandidat::Int64 = choixcandidat(util, alpha)

        #On ajoute le candidat à la solution
        SOL[PosCandidat] = 1

        #Desactive! le candidat selectionné
        Desactive!(PosCandidat, matrix, actif, m, varactive,n)
        ite=ite+1
    end
    #println(calculz(SOL,cost,n))
    #return(SOL)
    Z = calculz(SOL,cost,n)
    println("Solution avant améliration: " , Z)
#    return(SOL, actif, Z)
    return(SOL, Z)

end

# =========================================================================== #

#Detremine an utility  based on an active (actif) set of matrix's lines
function Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)
#    println("==========================")
    util = zeros(Float64, n) #réinitialisation du vecteur

    #On each column of matrix
    for j=1:n
            if varactive[j] == 1
                    #For each line
                    for i=1:m
                        #checking if the line is active
                        if actif[i]==0 #&& i <= m
                            K=matrix[i,j]
                            util[j]=util[j]+K
#                            println("====== \n", "i: ", i ,"    j: ", j , "   k: ", K )
#                            println("\n varactive: ", varactive)
#                            println("\n actif : ", actif)
#                            println("\n util : ", util)
                        end
                    end

                # dividing the number if iterations of the variable j by its cost (not a zero)

                if util[j]!=false #On évide de divisé par 0...
#                    println("\n util avant modif: ", util)
                    util[j] =  cost[j] / util[j]
#                    println("\n util après modif: ", util)

                else
                    varactive[j] = 0
                end
            end
    end
#    println("\n Util: ", util)
#    println("\n \n ===================Utilite fini =============")
    return util
end

# =========================================================================== #
#=
function PosMax(util,n) #On choisie un candidat parmie les différents coûts, par ordre lexicographique
#println("%%%%%%%%%%%%%%%%%%%%")
#    Pos::Int64
    Pos=0
    ValCan=0.0

    for i=1:n

        if util[i] > ValCan

            ValCan = util[i]
            Pos = i
        end
    end
    #println("\n \n %%%%%% Pos max fini %%%%% \n pos: ", Pos)
    return Pos
end
=#

# =========================================================================== #

function Desactive!(PosCandidat, matrix, actif, m , varactive , n) #On désactive les lignes où est le candidat
#println("\n======================= desactive ==============================")
#println("~~~~~~~~~~~~~~~~~~~~~~~~")
    for i=1:m
#        println("\n Boucle i: ", i)
        if matrix[i,PosCandidat] == 1
#            println("\n If matrix[i,PosCandidat] == 1: " , matrix[i,PosCandidat])
            actif[i] = 1
#            println("\n actif[i] = 1: ", actif[i])
            for j=1:n
#                println("\n for j=1:n ; j: ", j)
                if matrix[i,j] ==1
#                    println("\n if matrix[i,j] ==1 :", matrix[i,j])
                    varactive[j] = 0
#                    println("\n varactive[j] = 0: ", varactive)
                end
            end
        end
    end
#    println("\n \n ~~~~~~~~~Desactive fini! ~~~~~~~\n ")
#    println("varactive: ", varactive)
#    println("actif: ", actif)
end

#=========================================================================#

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end

#==========================================#
function choixcandidat(util, alpha)
    max = maximum(util)
    min = minimum(util)
#    Rut = Array{Float64}[]
    Rut = Float64[]
#    println("util: ", util, "\n type de util: ", typeof(util))
#    println("min: ", min, " max: ", max, " Rut: ", Rut)
    cpt::Int=1
    for i in 1:length(util)
#        println("debut for: ", i)
        if util[i]>= min+alpha*(max-min)
            push!(Rut,i)
        end
    end
    prob=rand(1:length(Rut))
#    println("fini, Rut: ", Rut, " Prob: ", prob)
    return(Rut[prob])
end
=#

#==
#Version 0
function choixcandidat(util, alpha)
    max = maximum(util)
    min = minimum(util)
#    Rut = Array{Float64}[]
    Rut = Float64[]
#    println("util: ", util, "\n type de util: ", typeof(util))
#    println("min: ", min, " max: ", max, " Rut: ", Rut)
    cpt::Int=1
    for i in 1:length(util)
#        println("debut for: ", i)
        if util[i]>= min+alpha*(max-min)
            push!(Rut,util[i])
        end
    end
    prob=rand(1:length(Rut))
#    println("fini, Rut: ", Rut, " Prob: ", prob)
    return(Rut[prob])
end

=#
#=
#Test run
z = Array{Float64}(undef, 42)
alpha = 0.9
println("z: ", z, " alpha : ", alpha)
cand=choixcandidat(z,alpha)
println("\n candidat choisit: ",cand)
==#
#=========================#
#=

DM1_1

# =========================================================================== #

# =========================================================================== #

#Implementing an x0 viable solution huristic with a gluton method
#receiving an array of lenght()=m , a matrix of size m*n and thoose size
#returing a Bool array of lenght=n
#Asked question:
#Mettre en place une heuristique de construction d’une x0 réalisable.

function Glouton(cost, matrix, n, m)

    #Creating a set of lines that will be avaluated
    actif= zeros(Bool, m) #if actif[j]=0 the line will be avaluated
    stop1= ones(Bool,m)

    #Creation the an array of activated variables
    varactive = ones(Bool, n)
    stop2= zeros(Bool , n)
#    println("\n varactive: ", varactive, "\n stop2: ", stop2)

    #creating my utility array for my m variables
    util = zeros(Float64, n)

    #initialisation de la Solution
    SOL = zeros(Float64, n)

#    println( "\n Matrix: ", matrix, "\n Cost: ", cost,"\n n= ", n, "   m= ", m,"\n actif: ", actif, "\n stop: ", stop1, "\n util: ", util, "\n SOL: ", SOL, "\n \n")

    ite = 0
    #Création de la solution
    while actif!=stop1 && varactive!=stop2

#      println("\n ++++++++++++++++++++++++++++++++++++++++")
#      println("\n itération: ", ite)
      #println("\n actif : ", actif, "\n stop1 : " , "\n varactive : " , varactive , "\n stop2 : ")
#      println("\n ++++++++++++++++++++++++++++++++++++++++")

        #Foction d'utilité

        util = Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)

        #Choix parmis les candidats de util
#        PosCandidat::Int64
        PosCandidat = PosMax(util,n)

        #On ajoute le candidat à la solution
        SOL[PosCandidat] = 1

        #Desactive! le candidat selectionné
        Desactive!(PosCandidat, matrix, actif, m, varactive,n)
        ite=ite+1
    end
    #println(calculz(SOL,cost,n))
    #return(SOL)
    Z = calculz(SOL,cost,n)
    println("Solution avant améliration: " , Z)
#    return(SOL, actif, Z)
    return(SOL, Z)

end

# =========================================================================== #

#Detremine an utility  based on an active (actif) set of matrix's lines
function Utilite(cost, matrix, actif, n, m , varactive, stop1, stop2)
#    println("==========================")
    util = zeros(Float64, n) #réinitialisation du vecteur

    #On each column of matrix
    for j=1:n
            if varactive[j] == 1
                    #For each line
                    for i=1:m
                        #checking if the line is active
                        if actif[i]==0 #&& i <= m
                            K=matrix[i,j]
                            util[j]=util[j]+K
#                            println("====== \n", "i: ", i ,"    j: ", j , "   k: ", K )
#                            println("\n varactive: ", varactive)
#                            println("\n actif : ", actif)
#                            println("\n util : ", util)
                        end
                    end

                # dividing the number if iterations of the variable j by its cost (not a zero)

                if util[j]!=false #On évide de divisé par 0...
#                    println("\n util avant modif: ", util)
                    util[j] =  cost[j] / util[j]
#                    println("\n util après modif: ", util)

                else
                    varactive[j] = 0
                end
            end
    end
#    println("\n Util: ", util)
#    println("\n \n ===================Utilite fini =============")
    return util
end

# =========================================================================== #

function PosMax(util,n) #On choisie un candidat parmie les différents coûts, par ordre lexicographique
#println("%%%%%%%%%%%%%%%%%%%%")
#    Pos::Int64
    Pos=0
    ValCan=0.0

    for i=1:n

        if util[i] > ValCan

            ValCan = util[i]
            Pos = i
        end
    end
    #println("\n \n %%%%%% Pos max fini %%%%% \n pos: ", Pos)
    return Pos
end

# =========================================================================== #

function Desactive!(PosCandidat, matrix, actif, m , varactive , n) #On désactive les lignes où est le candidat
#println("\n======================= desactive ==============================")
#println("~~~~~~~~~~~~~~~~~~~~~~~~")
    for i=1:m
#        println("\n Boucle i: ", i)
        if matrix[i,PosCandidat] == 1
#            println("\n If matrix[i,PosCandidat] == 1: " , matrix[i,PosCandidat])
            actif[i] = 1
#            println("\n actif[i] = 1: ", actif[i])
            for j=1:n
#                println("\n for j=1:n ; j: ", j)
                if matrix[i,j] ==1
#                    println("\n if matrix[i,j] ==1 :", matrix[i,j])
                    varactive[j] = 0
#                    println("\n varactive[j] = 0: ", varactive)
                end
            end
        end
    end
#    println("\n \n ~~~~~~~~~Desactive fini! ~~~~~~~\n ")
#    println("varactive: ", varactive)
#    println("actif: ", actif)
end

#=========================================================================#

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end
=#

max_stars_repo_path: Archives/Triste vie 10-08/brouillon.jl

max_issues_repo_path: Archives/Triste vie 10-08/brouillon.jl

max_forks_repo_path: Archives/Triste vie 10-08/brouillon.jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Libcroco"
version = v"0.6.13"

# Collection of sources required to build Libcroco
sources = [
    ArchiveSource("https://download.gnome.org/sources/libcroco/$(version.major).$(version.minor)/libcroco-$(version).tar.xz",
                  "767ec234ae7aa684695b3a735548224888132e063f92db585759b422570621d4")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libcroco-*/

FLAGS=()
if [[ "${target}" == *-apple-* ]]; then
    # We purposefully use an old binutils, so we must disable -Bsymbolic
    FLAGS+=(--disable-Bsymbolic)
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-gtk-doc "${FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    # Why must you flaunt the well-accepted ways of versioning your filename, libcroco?!
    # And even worse, why must you do so IN A SYNACTICALLY AMBIGUOUS MANNER?!
    LibraryProduct(["libcroco", "libcroco-$(version.major)", "libcroco-$(version.major).$(version.minor)"], :libcroco),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Glib_jll", v"2.59.0"; compat="2.59"),
    Dependency("XML2_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module SH

import Covering
import Voronoi

function intersectWithHalfSpace(W::Covering.Section, p::Tuple{Real, Real}, q::Tuple{Real, Real})
    pq = Voronoi.Geometry.subVector(q, p)
    v = Voronoi.Geometry.rotateVectorCCW(pq)

    function isVisible(point)
        function dot(a, b)
            return a[1]*b[1] + a[2]*b[2]
        end

        u = Voronoi.Geometry.subVector(point, p)

        return (dot(u, v) >= 0)
    end

    function intersectionPoint(a::Tuple{Real, Real}, b::Tuple{Real, Real})
        p1 = p
        dir1 = pq

        p2 = a
        dir2 = Voronoi.Geometry.subVector(b, a)

        vec = [p2[1] - p1[1], p2[2] - p1[2]]
        M = [dir1[1] dir2[1]; dir1[2] dir2[2]]
        res = M\vec
        return (p1[1] + res[1]*dir1[1], p1[2] + res[1]*dir1[2])
    end

    newW = Covering.Section((0, 0), nothing)
    global current = nothing
    function putPoint(p::Tuple{Real, Real})
        global current
        if current == nothing
            newW.borderHead = Covering.Edge(p, nothing)
            current = newW.borderHead
        else
            current.next = Covering.Edge(p, nothing)
            current = current.next
        end
    end

    w = W.borderHead
    ranOnce = false

    while w != W.borderHead || !ranOnce
        if w == nothing
            break
        end
        a = w.origin
        b = w.next.origin

        if isVisible(a)
            if isVisible(b)
                putPoint(b)
            else
                putPoint(intersectionPoint(a, b))
            end
        elseif isVisible(b)
            putPoint(intersectionPoint(a, b))
            putPoint(b)
        else
            # print none
        end

        w = w.next
        ranOnce = true
    end

    if current != nothing
        current.next = newW.borderHead
    end

    return newW
end

function intersectOne(Vᵢ::Voronoi.Diagram.Region, A::Covering.Section)
    W = Covering.Section(A)

    W = intersectWithHalfSpace(W, Vᵢ.borderHead.origin, Vᵢ.borderHead.next.origin)

    v = Vᵢ.borderHead.next
    while v != nothing && v != Vᵢ.borderHead
        W = intersectWithHalfSpace(W, v.origin, v.next.origin)
        v = v.next
    end

    return W
end

function intersectWithCircle(W::Covering.Section, center::Tuple{Real, Real}, r::Real)
    function isVisible(point)
        dist = Voronoi.Geometry.distanceSquared(center, point)

        return (dist <= r^2)
    end

    newW = Covering.Section(center, nothing)
    global current = nothing
    function putPoint(p::Tuple{Real, Real}, isArc=false)
        global current
        if current == nothing
            if isArc
                newW.borderHead = Covering.Arc(p, nothing)
            else
                newW.borderHead = Covering.Edge(p, nothing)
            end
            current = newW.borderHead
        else
            if isArc
                current.next = Covering.Arc(p, nothing)
            else
                current.next = Covering.Edge(p, nothing)
            end
            current = current.next
        end
    end

    w = W.borderHead
    ranOnce = false

    while w != W.borderHead || !ranOnce
        if w == nothing
            break
        end
        a = w.origin
        b = w.next.origin

        if isVisible(a)
            if isVisible(b)
                putPoint(b)
            else
                putPoint(Covering.segmentArcIntersections(a, b, center, r)[1], true) # put arc
            end
        elseif isVisible(b)
            putPoint(Covering.segmentArcIntersections(a, b, center, r)[1])
            putPoint(b)
        else
            intersections = Covering.segmentArcIntersections(a, b, center, r)

            # arcs with no intersection or just tangential are ignored
            if length(intersections) == 2
                putPoint(intersections[1])
                putPoint(intersections[2], true) # put arc
            end
        end

        w = w.next
        ranOnce = true
    end

    if current != nothing
        current.next = newW.borderHead
    end

    return newW
end

function intersect(V::Voronoi.Diagram.DCEL, A::Covering.Section, r::Real)
    n = length(V.regions)
    sections = Array{Covering.Section, 1}(undef, n)

    for i in 1:n
        sections[i] = intersectWithCircle(intersectOne(V.regions[i], A), V.regions[i].generator, r)
    end

    return Covering.Partition(sections, r)
end

end

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
using BasesAndSamples
using Test

using AbstractAlgebra # for polynomial bases

@testset "BasesAndSamples.jl" begin
    # Write your tests here.
    @testset "Bases" begin
        R, (x,) = PolynomialRing(RealField, ["x"])
        R2, (u,v) = PolynomialRing(RealField, ["u", "v"])
        # This is the intuitive order
        @test basis_monomial(3, x) == [R(1), x, x^2, x^3]
        @test basis_monomial(3, u,v) == [R2(1), u, v, u^2, u*v, v^2, u^3, u^2*v, u*v^2,v^3]

        #this shows that the functions generate a basis, but this doesn't mean that it is entirely correct
        functions = [basis_laguerre, basis_chebyshev, basis_jacobi,basis_gegenbauer]
        d = 3
        params = [[d,3,x],[d,x],[d,3,3,x],[d,3,x]]
        for i=1:length(functions)
            @test total_degree.(functions[i](params[i]...)) == [i for i=0:d]
        end
    end

    @testset "samples" begin
        #This tests whether there are enough points, not whether the points are correct
        #univariate tests:
        d = 5
        @test length(sample_points_chebyshev(d)) == d+1
        @test length(sample_points_chebyshev_mod(d)) == d+1
        @test length(sample_points_simplex(1,d)) == d+1
        @test length(sample_points_rescaled_laguerre(d)) == d+1

        #2 variate
        @test length(sample_points_padua(d)) == binomial(2+d,2)
        @test length(sample_points_simplex(2,d)) == binomial(2+d,2)

        #3 variate
        @test length(sample_points_simplex(3,d)) == binomial(3+d,3)
    end
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Abundant
using NumberTheory, IterTools, Counts

export ModuleAbundant
export isAbundant, is005101, I005101, F005101, L005101, V005101

"""
``n`` is an abundant number if ``σ(n) > 2n``. An abundant number is a number for which the sum of its proper divisors is greater than the number itself.

* isAbundant, is005101, I005101, F005101, L005101, V005101.
"""
const ModuleAbundant = ""

"""
Is ``n`` an abundant number, i.e. is ``σ(n) > 2n`` ?
"""
isAbundant(n) = σ(n) > 2n

"""
Is ``n`` a term of sequence A005101?
"""
is005101(n) = isAbundant(n)

"""
Iterate over the first ``n`` abundant numbers.
"""
I005101(n) = takeFirst(isAbundant, n)

"""
Iterate over the abundant numbers which do not exceed ``n (1 ≤ i ≤ n)``.
"""
F005101(n) = Iterators.filter(isAbundant, 1:n)

"""
Return a list of the  first ``n`` abundant numbers.
"""
L005101(n) = collect(I005101(n))

"""
Return the value of the ``n``-th abundant number.
"""
V005101(n) = nth(I005101(n), n)

#START-TEST-########################################################

using Test, SeqTests, SeqUtils

function test()
    @testset "Abundant" begin
        @test isAbundant(100800) == true
        @test isAbundant(2402400) == true
        @test isAbundant(49008960) == true

        if is_oeis_installed()
            SeqTest([L005101], 'L')
        end
    end
end

function demo()
    println(V005101(15))
    println(L005101(15))

    for a in I005101(15) print(a, ", ") end; println("...")
    for a in F005101(40) print(a, ", ") end; println("...")

    for n in 40:50 println(n, " ↦ ", isAbundant(n)) end
    for n in 1:6   println(n, " ↦ ", L005101(n)) end

    for a in I005101(1000)
        isodd(a) && print(a, " ")
    end
    println()
end

"""
L005101(5000)
    0.018731 seconds (40.38 k allocations: 759.328 KiB)
"""
function perf()
    GC.gc()
    @time L005101(5000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
Test Summary: | Pass  Total
Abundant      |    3      3

72

[12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56, 60, 66, 70, 72]
12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56, 60, 66, 70, 72, ...
12, 18, 20, 24, 30, 36, 40, ...

40 ↦ true
41 ↦ false
42 ↦ true
43 ↦ false
44 ↦ false
45 ↦ false
46 ↦ false
47 ↦ false
48 ↦ true
49 ↦ false
50 ↦ false

1 ↦ [12]
2 ↦ [12, 18]
3 ↦ [12, 18, 20]
4 ↦ [12, 18, 20, 24]
5 ↦ [12, 18, 20, 24, 30]
6 ↦ [12, 18, 20, 24, 30, 36]

945 1575 2205 2835 3465

0.018731 seconds (40.38 k allocations: 759.328 KiB)
=#


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

n = 50
κ = 1000
A = SymTridiagonal(2*ones(n), -ones(n-1))
A = (κ-1)*A/4 + eye(n)
cond(A)
b = zeros(n); b[1] = (κ - 1)/4

Λ, Q = eig(A)

κ = Λ[1]/Λ[n]
α = 2/(Λ[1] + Λ[end])
β = 0.98
#β = 0

# Function from momentum2.jl
Z = cat(2,[ runmomentumcheat(A,b,i)[2] for i = 1:400]...)
heatmap(Z)

plot(b2)
plot!(A\b)

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
#! /usr/bin/julia

# Rosetta Code, Rep-string

function list_reps{T<:String}(r::T)
    n = length(r)
    replst = T[]
    for m in 1:n>>1
        s = r[1:chr2ind(r,m)]
        (s^(div(n,m)+1))[1:chr2ind(r,n)] == r || continue
        push!(replst, s)
    end
    return replst
end

tests = {"1001110011",
         "1110111011",
         "0010010010",
         "1010101010",
         "1111111111",
         "0100101101",
         "0100100",
         "101",
         "11",
         "00",
         "1",
         "\u2200\u2203\u2200\u2203\u2200\u2203\u2200\u2203"}

for r in tests
    replst = list_reps(r)
    rlen = length(replst)
    print(@sprintf("    %s ", r))
    if rlen == 0
        println("is not a rep-string.")
    else
        println("is a rep-string of ", join(replst, ", "), ".")
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
#=

  A coin application in Julia ConstraintSolver.jl

  From "The ECLiPSe Book" pages 99f and 234 ff
  The solution in ECLiPSe is at page 236.

  """
  What is the minimum number of coins that allows one to pay _exactly_
  any amount smaller than one Euro? Recall that there are six different
  euro cents, of denomination 1, 2, 5, 10, 20, 50
  """

  There are 4 optimal solutions (8 coins):
    x=[1,2,1,2,1,1]
    x=[1,2,2,1,1,1]
    x=[2,1,1,2,1,1]
    x=[2,1,2,1,1,1]

  This problem was implemented in coins3.jl

  This model implement a related problem:
  Is there a set of coin denominations which make it possible to change
  into 1..99 with less than 8 coins?

  Answer: Yes, there is. For example the following configuration:

  num_coins:7
  variables:[1,2,3,6,12,25,50]
  x:[1,1,1,1,1,1,1]

  Though, it come with a cost of using 7 different coin denominations.

  Also, we can see that there is a set with just 5 different
  denominations for the 8 coin change. However, it requires that
  more than one coin of some of the denominations

  This is done by comment this constraint:
        NumCoins =< 8,

  Result:
    num_coins:8
    variables:[1,2,4,11,33]
    x:[1,1,2,2,2]

Here's the run. As we see it can be done with:
- 8 coins of 5 (and 6) different denominations
- 7 coins of 7 different denominations.

"""
n:5
num_sols:1

solution #1
num_coins:8
coins:[1, 2, 4, 11, 33]
xs:[1, 1, 2, 2, 2]
 15.108249 seconds (24.47 M allocations: 12.097 GiB, 56.45% gc time)

n:6
num_sols:1

solution #1
num_coins:8
coins:[1, 2, 3, 5, 11, 33]
xs:[1, 1, 1, 1, 2, 2]
 21.561126 seconds (34.97 M allocations: 14.951 GiB, 59.19% gc time)

n:7
num_sols:1

solution #1
num_coins:7
coins:[1, 2, 3, 6, 13, 24, 50]
xs:[1, 1, 1, 1, 1, 1, 1]
 15.875412 seconds (16.94 M allocations: 14.988 GiB, 64.24% gc time)
"""

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function coins3_2(n,print_solutions=true,all_solutions=false,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy" => :BFS,
                                                            # "traverse_strategy" => :DFS,
                                                            # "traverse_strategy" => :DBFS,

                                                            "branch_split" => :Smallest,
                                                            # "branch_split" => :Biggest,
                                                            # "branch_split" => :InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    # Which coins should we use?
    # we assume that there are no coin of a denomination larger than 50.
    @variable(model, 1 <= coins[1:n] <= 50, Int)

    # Total number of coins used. Can we do it in 8 or less than 8 coins?
    @variable(model, 1 <= num_coins <= 8, Int)

    # Use either 1 or 2 coins of each denomination
    @variable(model, 1 <= xs[1:n] <= 2, Int)

    @constraint(model, coins in CS.AllDifferent())
    increasing_strict(model, coins) # symmetry breaking

    # Checks that all changes from 1 to 99 can be made.
    table = resize_matrix([ [i,j,i * j] for i in 0:99, j in 0:99])
    for j in 1:99
       tmp = @variable(model, [1:n], CS.Integers(0:99*99))

       # scalar_product(model,coins,tmp,j) # non-linear constraint
       # Reformulation:
       t = @variable(model, [1:n], CS.Integers(0:99*99)) # the product of coins[i] * tmp[i]
       for i in 1:n
           @constraint(model, [coins[i], tmp[i], t[i]] in CS.TableSet(table))
           # Ensure that we don't use more coins than in x[i]
           @constraint(model, tmp[i] <= xs[i])
       end
       @constraint(model, j == sum(t))

    end
    @constraint(model,num_coins == sum(xs))

    @objective(model,Min,num_coins)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                coins_val = convert.(Integer,JuMP.value.(coins; result=sol))
                xs_val = convert.(Integer,JuMP.value.(xs; result=sol))
                num_coins_val = convert.(Integer,JuMP.value.(num_coins; result=sol))
                println("num_coins:$num_coins_val")
                println("coins:$coins_val")
                println("xs:$xs_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

begin
  for n in 1:7
    println("\nn:$n")
    @time status = coins3_2(n ,true,false,125)
    if status == MOI.TIME_LIMIT
      break
    end
  end
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Dates

mutable struct InterpolatedDiscountCurve{DC <: DayCount, P <: Interpolation, B <: BusinessCalendar} <: InterpolatedCurve{P}
  settlementDays::Int
  referenceDate::Date
  dc::DC
  interp::P
  cal::B
  dates::Vector{Date}
  times::Vector{Float64}
  data::Vector{Float64}
end

function InterpolatedDiscountCurve(dates::Vector{Date},
                                  discounts::Vector{Float64},
                                  dc::DC,
                                  interpolator::P) where {DC <: DayCount, P <: Interpolation}
  idc = InterpolatedDiscountCurve{DC, P, NullCalendar}(0, dates[1], dc, interpolator, NullCalendar(), dates, zeros(length(dates)), discounts)

  initialize!(idc)

  return idc
end

function initialize!(idc::InterpolatedDiscountCurve)
  length(idc.dates) == length(idc.data) || error("dates / data mismatch")
  idc.times[1] = 0.0
  @simd for i = 2:length(idc.dates)
    @inbounds idc.times[i] = year_fraction(idc.dc, idc.dates[1], idc.dates[i])
  end

  # initialize interpolator
  Math.initialize!(idc.interp, idc.times, idc.data)

  Math.update!(idc.interp)

  return idc
end

function discount_impl(idc::InterpolatedDiscountCurve, t::Float64)
  if t <= idc.times[end]
    return Math.value(idc.interp, t)
  end

  # flat fwd extrapolation
  tMax = idc.times[end]
  dMax = idc.data[end]
  instFwdMax = Math.derivative(idc.interp, tMax) / dMax
  return dMax * exp(-instFwdMax * (t - tMax))
end


text:
"""
    aras(decisionMat, weights, fs)
Apply ARAS (Additive Ratio ASsessment) method for a given matrix, weights and, type of criteria.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fs::Array{Function,1}`: m-vector of type of criteria. The benefit criteria shown with "maximum", and the cost criteria shown with "minimum".

 # Description
aras() applies the ARAS method to rank n alternatives subject to m criteria and criteria type vector.

# Output
- `::ARASResult`: ARASResult object that holds multiple outputs including scores and best index.

# Examples
```julia-repl
julia> df = DataFrame(
    :K1 => [105000.0, 120000, 150000, 115000, 135000],
    :K2 => [105.0, 110, 120, 105, 115],
    :K3 => [10.0, 15, 12, 20, 15],
    :K4 => [4.0, 4, 3, 4, 5],
    :K5 => [300.0, 500, 550, 600, 400],
    :K6 => [10.0, 8, 12, 9, 9]
)
julia> df
5×6 DataFrame
 Row │ K1        K2       K3       K4       K5       K6
     │ Float64   Float64  Float64  Float64  Float64  Float64
─────┼───────────────────────────────────────────────────────
   1 │ 105000.0    105.0     10.0      4.0    300.0     10.0
   2 │ 120000.0    110.0     15.0      4.0    500.0      8.0
   3 │ 150000.0    120.0     12.0      3.0    550.0     12.0
   4 │ 115000.0    105.0     20.0      4.0    600.0      9.0
   5 │ 135000.0    115.0     15.0      5.0    400.0      9.0
julia> result = aras(df, w, fs);
julia> result.bestIndex
2
julia> result.scores
5-element Array{Float64,1}:
 0.8142406768388222
 0.8928861957614441
 0.764157900073527
 0.8422546181927358
 0.8654063509472654
```
# References
Zavadskas, E. K., & Turskis, Z. (2010). A new additive ratio assessment (ARAS) method in multicriteria decision‐making. Technological and Economic Development of Economy, 16(2), 159-172.
Yıldırım, B. F. (2015). "Çok Kriterli Karar Verme Problemlerinde ARAS Yöntemi". Kafkas Üniversitesi İktisadi ve İdari Bilimler Fakültesi Dergisi, 6 (9), 285-296. http://dx.doi.org/10.18025/kauiibf.65151
"""
function aras(decisionMat::DataFrame, weights::Array{Float64,1}, fs::Array{Function,1})::ARASResult

    # mat = convert(Matrix, decisionMat)
    mat = Matrix(decisionMat)

    nrows, ncols = size(mat)
    w = unitize(weights)
    referenceRow = apply_columns(fs, mat)
    extendMat = [mat; referenceRow']

    for col in 1:ncols
        if fs[col] == minimum
            extendMat[:, col] = 1 ./ extendMat[:, col]
        end
    end

    normalizedMat = similar(extendMat)

    for col in 1:ncols
        for row in 1:nrows + 1
            normalizedMat[row, col] = extendMat[row, col] ./ sum(extendMat[:, col])
        end
    end

    # optimality = similar(normalizedMat)

    optimality_degrees = zeros(Float64, nrows + 1)
    for i in 1:nrows + 1
        optimality_degrees[i] = w .* normalizedMat[i, :] |> sum
    end

    utility_degrees = zeros(Float64, nrows)
    for i in 1:nrows
        utility_degrees[i] = optimality_degrees[i] /  optimality_degrees[end]
    end

    orderings = sortperm(utility_degrees)
    bestIndex = orderings |> last

    result = ARASResult(
        referenceRow,
        extendMat,
        normalizedMat,
        optimality_degrees,
        utility_degrees,
        orderings,
        bestIndex
    )

    return result
end

"""
    aras(setting)
Apply ARAS (Additive Ratio ASsessment) method for a given matrix, weights and, type of criteria.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

 # Description
aras() applies the ARAS method to rank n alternatives subject to m criteria and criteria type vector.

# Output
- `::ARASResult`: ARASResult object that holds multiple outputs including scores and best index.
"""
function aras(setting::MCDMSetting)::ARASResult
    aras(
        setting.df,
        setting.weights,
        setting.fns
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
@with_polymake @testset "Eichler" begin
  K, a = Hecke.rationals_as_number_field()
  A = Hecke.AlgQuat(K, K(1), K(1))
  A = Hecke.AlgAss(K, A.mult_table)
  OA = maximal_order(A)
  b = @inferred Hecke.principal_generator_eichler(2 * OA)
  @test b * OA == 2 * OA
end


text:
##############################
#      Moral Vector Type     #
##############################

"""
    type MoralVector{K, T <: Real}

Type representing a Moral Vector

### Fields
* `morals` (Vector{T}): normalized `K`-vector representing the moral dimension values (or, moral vector)
"""
struct MoralVector{K, T <: Real}
    morals::Vector{T}
end

"""
    MoralVector()

Construct a MoralVector with a random moral vector with default number of components `KMORAL`.
"""
MoralVector() = MoralVector{KMORAL, Float64}(randSphere(KMORAL))

"""
    MoralVector(k::Int)

Construct a MoralVector with a random moral vector with default number of components `k`.
"""
MoralVector(k::Int) = MoralVector{k, Float64}(randSphere(k))

"""
    MoralVector(m::Vector{T})

Construct a MoralVector with an input moral vector `m` to be normalized.
"""
MoralVector(m::Vector{T}) where T <: Real = MoralVector{length(m), T}(normalize(m))

###############################
# Moral Vector Redefinitions  #
###############################

length(mvector::MoralVector) = length(mvector.morals)
size(mvector::MoralVector)   = size(mvector.morals)

getindex(mvector::MoralVector, i::Int)   = mvector.morals[i]
getindex(mvector::MoralVector, c::Colon) = mvector.morals[:]

start(mvector::MoralVector)        = 1
done(mvector::MoralVector, s::Int) = s > length(mvector)
next(mvector::MoralVector, s::Int) = (mvector[s], s+1)

dot{K,T}(i::MoralVector{K,T}, j::MoralVector{K,T}) = (i[:] ⋅ j[:])

function show(io::IO, mvector::MoralVector)
    N = length(mvector)
    print(io, N, "-dimensional Moral Vector:")
    for i in mvector
        @printf io "\n %.5f" i
    end
end

+(mvector::MoralVector, nvector::MoralVector) = MoralVector(mvector[:] + nvector[:])
-(mvector::MoralVector, nvector::MoralVector) = MoralVector(mvector[:] - nvector[:])
-(mvector::MoralVector) = MoralVector(-mvector[:])

# need a better way to sum weigthed MoralVector's
*(p::T, mvector::MoralVector{K, T}) where {T <: Real, K} = p*mvector[:]

###############################
#  Moral Vector Definitions   #
###############################

"""
    consensus(ags::Vector{MoralVector}, i::Int, j::Int)

Evaluates the consensus of MoralVectors i and j from a Vector of MoralVectors
"""
function consensus(ags::Vector{MoralVector{K, T}}, i::Int, j::Int) where {T <: Real, K}
    N = length(ags)
    if 0 < i <= N && 0 < j <= N
        return ags[i] ⋅ ags[j]
    else
        error("The indices given are out of range")
    end
end

"""
    consensus(ags::::Vector{MoralVector})

Evaluates the consensus of Vector of MoralVectors, that is, the matrix of the correlations between the MoralVectors inner representations
"""
function consensus(ags::Vector{MoralVector{K, T}}; format = :none) where {T <: Real, K}
    N = length(ags)
    ψ = zeros(N, N)

    for i in 1:N for j in (i+1):N
        ψ[i, j] = consensus(ags, i, j)
    end end

    if format == :ordered
        return ψ[ψ .!= 0.]
    end
    return ψ
end


text:
module Healpix
###
export nsideok, nside2pixarea, nside2resol
export Resolution, nside2npix, npix2nside
export ang2pixNest, ang2pixRing, pix2angNest, pix2angRing
export vec2pixNest, vec2pixRing, pix2vecNest, pix2vecRing
export pix2ringpos
export Order, Ring, Nest
export AbstractMap, Map, OrderedMap, MaskedMap
export getcoords, resize, reorder
export ang2vec, vec2ang, ang2pix, pix2ang
export readMapFromFITS, savePixelsToFITS, saveToFITS, conformables
export ringWeightPath, readWeightRing
export pixelWindowPath, readPixelWindowT, readPixelWindowP
export Alm, numberOfAlms, almIndexL0, almIndex, readAlmFromFITS
export getringinfo!, getringinfo, getinterpolRing
export pix2xyfRing, xyf2pixRing, pix2xyfNest, xyf2pixNest
export ring2nest, nest2ring

import FITSIO
import Base: getindex, setindex!

const NSIDE_MAX = 8192

using PyCall; const healpy = PyNULL()
function __init__()
    copy!(healpy, pyimport("healpy"))
end
export healpy

include("math.jl")
include("datatables.jl")
include("resolution.jl")
include("interp.jl")
include("xyf.jl")

include("maps.jl")
include("pixels.jl")
include("projections.jl")
include("alm.jl")
include("mapmaking.jl")
include("fits.jl")

end


text:
using Test
using MarketObservables
using MarketObservables: OrderBook, OrderBookException
using MarketObservables: FiniteDepthOrderBook, InfiniteDepthOrderBook
using MarketObservables: Level, price, bid, ask, spread
using MarketObservables: Volume, null_volume

@testset "Test Level" begin
    l1 = Level(1.34, 100.2)
    l2 = Level(1.35, 110.2)

    @test l1 < l2
end

@testset "Test FiniteDepthOrderBook" begin
    asks = [Level(110.0, 10.0), Level(111.0, 12.0)]  # asks (ascending asks)
    bids = [Level(100.0, 10.0), Level(99.0, 15.0)]  # bids (descending prices)

    ob = OrderBook(asks, bids)
    @test typeof(ob) <: FiniteDepthOrderBook

    @test bid(ob)[1] == 100.0
    @test ask(ob)[1] == 110.0
    @test spread(ob)[1] == 10.0

    @testset "Low volume" begin
        volume = Volume(15.0)
        expected_bid = (100.0 * 10.0 + 99.0 * 5.0) / 15.0
        calc_bid, = bid(ob, volume)
        @test calc_bid == expected_bid
        expected_ask = (110.0 * 10.0 + 111 * 5.0) / 15.0
        @test ask(ob, volume)[1] == expected_ask
        @test spread(ob, volume)[1] == expected_ask - expected_bid
        @test price(ob, volume) == ask(ob, volume)
        @test price(ob, -volume) == bid(ob, volume)
    end

    @testset "Volume too big for depth of orderbook" begin
        volume = Volume(50.0)
        _spread, _remaining_volume_ask, _remaining_volume_bid = spread(ob, volume)
        @test _remaining_volume_ask != 0
        @test _remaining_volume_bid != 0
    end

    @testset "Construction errors" begin
        asks = [Level(111.0, 10.0), Level(110.0, 12.0)]  # asks (NOT ascending asks)
        bids = [Level(100.0, 10.0), Level(99.0, 15.0)]  # bids (descending prices)
        @test_throws OrderBookException OrderBook(asks, bids)

        asks = [Level(110.0, 10.0), Level(111.0, 12.0)]  # asks (ascending asks)
        bids = [Level(99.0, 10.0), Level(100.0, 15.0)]  # bids (NOT descending prices)
        @test_throws OrderBookException OrderBook(asks, bids)

        asks = [Level(111.0, 10.0), Level(111.0, 12.0)]  # asks (NOT unique prices)
        bids = [Level(100.0, 10.0), Level(99.0, 15.0)]  # bids (descending prices)
        @test_throws OrderBookException OrderBook(asks, bids)

        asks = [Level(110.0, 10.0), Level(111.0, 12.0)]  # asks (ascending asks)
        bids = [Level(100.0, 10.0), Level(100.0, 15.0)]  # bids (NOT unique prices)
        @test_throws OrderBookException OrderBook(asks, bids)
    end

end

@testset "Test InfiniteDepthOrderBook" begin
    @testset "define with bid/ask" begin
        _ask, _bid = 110.0, 100.0
        ob = OrderBook(_ask, _bid)  # with bid/ask
        @test typeof(ob) <: InfiniteDepthOrderBook
        @test spread(ob)[1] == 10.0

        volume = Volume(50.0)
        @test bid(ob)[1] == 100.0
        @test ask(ob)[1] == 110.0
        @test bid(ob, volume)[1] == 100.0
        @test ask(ob, volume)[1] == 110.0
        _spread, _remaining_volume_ask, _remaining_volume_bid = spread(ob, volume)
        @test _remaining_volume_ask == 0
        @test _remaining_volume_bid == 0
    end

    @testset "define with price (bid==ask)" begin
        price = 100.0
        ob = OrderBook(price)
        @test typeof(ob) <: InfiniteDepthOrderBook
        volume = Volume(50.0)
        @test bid(ob)[1] == price
        @test ask(ob)[1] == price
        _spread, = spread(ob, volume)
        @test _spread == 0
    end

    @testset "Construction errors" begin
        _ask, _bid = 100.0, 110.0  #  NOT ask >= bid
        @test_throws OrderBookException OrderBook(_ask, _bid)
    end
end