text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: service_quotas

using AWS.Compat
using AWS.UUIDs
"""
    AssociateServiceQuotaTemplate()

Associates the Service Quotas template with your organization so that when new accounts are created in your organization, the template submits increase requests for the specified service quotas. Use the Service Quotas template to request an increase for any adjustable quota value. After you define the Service Quotas template, use this operation to associate, or enable, the template.

"""

associate_service_quota_template(; aws_config::AWSConfig=global_aws_config()) = service_quotas("AssociateServiceQuotaTemplate"; aws_config=aws_config)
associate_service_quota_template(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("AssociateServiceQuotaTemplate", args; aws_config=aws_config)

"""
    DeleteServiceQuotaIncreaseRequestFromTemplate()

Removes a service quota increase request from the Service Quotas template.

# Required Parameters
- `AwsRegion`: Specifies the AWS Region for the quota that you want to delete.
- `QuotaCode`: Specifies the code for the quota that you want to delete.
- `ServiceCode`: Specifies the code for the service that you want to delete.

"""

delete_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("DeleteServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
delete_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("DeleteServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    DisassociateServiceQuotaTemplate()

Disables the Service Quotas template. Once the template is disabled, it does not request quota increases for new accounts in your organization. Disabling the quota template does not apply the quota increase requests from the template.   Related operations    To enable the quota template, call AssociateServiceQuotaTemplate.    To delete a specific service quota from the template, use DeleteServiceQuotaIncreaseRequestFromTemplate.

"""

disassociate_service_quota_template(; aws_config::AWSConfig=global_aws_config()) = service_quotas("DisassociateServiceQuotaTemplate"; aws_config=aws_config)
disassociate_service_quota_template(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("DisassociateServiceQuotaTemplate", args; aws_config=aws_config)

"""
    GetAWSDefaultServiceQuota()

Retrieves the default service quotas values. The Value returned for each quota is the AWS default value, even if the quotas have been increased..

# Required Parameters
- `QuotaCode`: Identifies the service quota you want to select.
- `ServiceCode`: Specifies the service that you want to use.

"""

get_awsdefault_service_quota(QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetAWSDefaultServiceQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_awsdefault_service_quota(QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetAWSDefaultServiceQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    GetAssociationForServiceQuotaTemplate()

Retrieves the ServiceQuotaTemplateAssociationStatus value from the service. Use this action to determine if the Service Quota template is associated, or enabled.

"""

get_association_for_service_quota_template(; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetAssociationForServiceQuotaTemplate"; aws_config=aws_config)
get_association_for_service_quota_template(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetAssociationForServiceQuotaTemplate", args; aws_config=aws_config)

"""
    GetRequestedServiceQuotaChange()

Retrieves the details for a particular increase request.

# Required Parameters
- `RequestId`: Identifies the quota increase request.

"""

get_requested_service_quota_change(RequestId; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetRequestedServiceQuotaChange", Dict{String, Any}("RequestId"=>RequestId); aws_config=aws_config)
get_requested_service_quota_change(RequestId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetRequestedServiceQuotaChange", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequestId"=>RequestId), args)); aws_config=aws_config)

"""
    GetServiceQuota()

Returns the details for the specified service quota. This operation provides a different Value than the GetAWSDefaultServiceQuota operation. This operation returns the applied value for each quota. GetAWSDefaultServiceQuota returns the default AWS value for each quota.

# Required Parameters
- `QuotaCode`: Identifies the service quota you want to select.
- `ServiceCode`: Specifies the service that you want to use.

"""

get_service_quota(QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetServiceQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_service_quota(QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetServiceQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    GetServiceQuotaIncreaseRequestFromTemplate()

Returns the details of the service quota increase request in your template.

# Required Parameters
- `AwsRegion`: Specifies the AWS Region for the quota that you want to use.
- `QuotaCode`: Specifies the quota you want.
- `ServiceCode`: Specifies the service that you want to use.

"""

get_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("GetServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    ListAWSDefaultServiceQuotas()

Lists all default service quotas for the specified AWS service or all AWS services. ListAWSDefaultServiceQuotas is similar to ListServiceQuotas except for the Value object. The Value object returned by ListAWSDefaultServiceQuotas is the default value assigned by AWS. This request returns a list of all service quotas for the specified service. The listing of each you'll see the default values are the values that AWS provides for the quotas.   Always check the NextToken response parameter when calling any of the List* operations. These operations can return an unexpected list of results, even when there are more results available. When this happens, the NextToken response parameter contains a value to pass the next call to the same API to request the next part of the list.

# Required Parameters
- `ServiceCode`: Specifies the service that you want to use.

# Optional Parameters
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
"""

list_awsdefault_service_quotas(ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListAWSDefaultServiceQuotas", Dict{String, Any}("ServiceCode"=>ServiceCode); aws_config=aws_config)
list_awsdefault_service_quotas(ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListAWSDefaultServiceQuotas", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    ListRequestedServiceQuotaChangeHistory()

Requests a list of the changes to quotas for a service.

# Optional Parameters
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from.
- `ServiceCode`: Specifies the service that you want to use.
- `Status`: Specifies the status value of the quota increase request.
"""

list_requested_service_quota_change_history(; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistory"; aws_config=aws_config)
list_requested_service_quota_change_history(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistory", args; aws_config=aws_config)

"""
    ListRequestedServiceQuotaChangeHistoryByQuota()

Requests a list of the changes to specific service quotas. This command provides additional granularity over the ListRequestedServiceQuotaChangeHistory command. Once a quota change request has reached CASE_CLOSED, APPROVED, or DENIED, the history has been kept for 90 days.

# Required Parameters
- `QuotaCode`: Specifies the service quota that you want to use
- `ServiceCode`: Specifies the service that you want to use.

# Optional Parameters
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from.
- `Status`: Specifies the status value of the quota increase request.
"""

list_requested_service_quota_change_history_by_quota(QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistoryByQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
list_requested_service_quota_change_history_by_quota(QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistoryByQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    ListServiceQuotaIncreaseRequestsInTemplate()

Returns a list of the quota increase requests in the template.

# Optional Parameters
- `AwsRegion`: Specifies the AWS Region for the quota that you want to use.
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from.
- `ServiceCode`: The identifier for a service. When performing an operation, use the ServiceCode to specify a particular service.
"""

list_service_quota_increase_requests_in_template(; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServiceQuotaIncreaseRequestsInTemplate"; aws_config=aws_config)
list_service_quota_increase_requests_in_template(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServiceQuotaIncreaseRequestsInTemplate", args; aws_config=aws_config)

"""
    ListServiceQuotas()

Lists all service quotas for the specified AWS service. This request returns a list of the service quotas for the specified service. you'll see the default values are the values that AWS provides for the quotas.   Always check the NextToken response parameter when calling any of the List* operations. These operations can return an unexpected list of results, even when there are more results available. When this happens, the NextToken response parameter contains a value to pass the next call to the same API to request the next part of the list.

# Required Parameters
- `ServiceCode`: The identifier for a service. When performing an operation, use the ServiceCode to specify a particular service.

# Optional Parameters
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from.
"""

list_service_quotas(ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServiceQuotas", Dict{String, Any}("ServiceCode"=>ServiceCode); aws_config=aws_config)
list_service_quotas(ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServiceQuotas", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    ListServices()

Lists the AWS services available in Service Quotas. Not all AWS services are available in Service Quotas. To list the see the list of the service quotas for a specific service, use ListServiceQuotas.

# Optional Parameters
- `MaxResults`: (Optional) Limits the number of results that you want to include in the response. If you don't include this parameter, the response defaults to a value that's specific to the operation. If additional items exist beyond the specified maximum, the NextToken element is present and has a value (isn't null). Include that value as the NextToken request parameter in the call to the operation to get the next part of the results. You should check NextToken after every operation to ensure that you receive all of the results.
- `NextToken`: (Optional) Use this parameter in a request if you receive a NextToken response in a previous request that indicates that there's more output available. In a subsequent call, set it to the value of the previous call's NextToken response to indicate where the output should continue from.
"""

list_services(; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServices"; aws_config=aws_config)
list_services(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("ListServices", args; aws_config=aws_config)

"""
    PutServiceQuotaIncreaseRequestIntoTemplate()

Defines and adds a quota to the service quota template. To add a quota to the template, you must provide the ServiceCode, QuotaCode, AwsRegion, and DesiredValue. Once you add a quota to the template, use ListServiceQuotaIncreaseRequestsInTemplate to see the list of quotas in the template.

# Required Parameters
- `AwsRegion`: Specifies the AWS Region for the quota.
- `DesiredValue`: Specifies the new, increased value for the quota.
- `QuotaCode`: Specifies the service quota that you want to use.
- `ServiceCode`: Specifies the service that you want to use.

"""

put_service_quota_increase_request_into_template(AwsRegion, DesiredValue, QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("PutServiceQuotaIncreaseRequestIntoTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
put_service_quota_increase_request_into_template(AwsRegion, DesiredValue, QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("PutServiceQuotaIncreaseRequestIntoTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)

"""
    RequestServiceQuotaIncrease()

Retrieves the details of a service quota increase request. The response to this command provides the details in the RequestedServiceQuotaChange object.

# Required Parameters
- `DesiredValue`: Specifies the value submitted in the service quota increase request.
- `QuotaCode`: Specifies the service quota that you want to use.
- `ServiceCode`: Specifies the service that you want to use.

"""

request_service_quota_increase(DesiredValue, QuotaCode, ServiceCode; aws_config::AWSConfig=global_aws_config()) = service_quotas("RequestServiceQuotaIncrease", Dict{String, Any}("DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
request_service_quota_increase(DesiredValue, QuotaCode, ServiceCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = service_quotas("RequestServiceQuotaIncrease", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), args)); aws_config=aws_config)


text:
@testset "norm" begin
   @test norm([3.0,4.0].') ≈ 5.0
   @test norm([3.0,4.0].', 1) ≈ 4.0
   @test norm([3.0,4.0].', Inf) ≈ 7.0
end

ext: jl


text:
# Put all the old model related

# model initialization and parameter selections
postagv(m)=m[1]; deprelv(m)=m[2]; lcountv(m)=m[3]; rcountv(m)=m[4]; distancev(m)=m[5]; parserv(m)=m[6];

function makepmodel1(d; GPUFEATS=false)
    m = ([d["postagv"],d["deprelv"],d["lcountv"],d["rcountv"],d["distancev"],d["parserv"]],
         [d["postago"],d["deprelo"],d["lcounto"],d["rcounto"],d["distanceo"],d["parsero"]])

    if gpu() >= 0
        if GPUFEATS
            return map2gpu(m)
        else
            m = map2cpu(m)
            m[1][6] = map2gpu(m[1][6])
            m[2][6] = map2gpu(m[2][6])
            return m
        end
    else
        return map2cpu(m)
    end
end

function makepmodel2(o, s; ftype=Float32, intype=:normal)
    model = Any[]
    dpostag, ddeprel, dcount = o[:embed]
    for (k,n,d) in ((:postag,17,dpostag),(:deprel,37,ddeprel),(:lcount,10,dcount),(:rcount,10,dcount),(:distance,10,dcount))
        if intype == :normal
            push!(model, [ initr(d) for i=1:n ])
            #@msg "Random normal initialization" # dbg purposes
        else
            push!(model, [ initzeros(d) for i=1:n ])
            #@msg "Zero initialization" # dbg purposes
        end
    end
    p = o[:arctype](s)
    f = features([p], o[:feats], model)
    mlpdims = (length(f), o[:hidden]..., p.nmove)
    info("mlpdims=$mlpdims")
    parser = Any[]
    for i=2:length(mlpdims)
        push!(parser, initx(mlpdims[i],mlpdims[i-1]))
        push!(parser, initx(mlpdims[i],1))
    end
    push!(model,parser)
    optim = optimizers(model,o[:optimization])
    return model,optim
end

function splitmodel(pmodel)
    mlpmodel = Any[]
    mlptemp = parserv(pmodel)
    for i=1:length(mlptemp); push!(mlpmodel, mlptemp[i]); end
    featmodel = Array{Any}(5) # Array(Any,5) deprecated
    for k in 1:5 # (:postag,:deprel,:lcount,:rcount,:distance)
        featmodel[k] = Any[]
        pmodel_k = pmodel[k]
        for i in 1:length(pmodel_k)
            push!(featmodel[k], pmodel_k[i])
        end
    end
    return (featmodel,mlpmodel)
end

makepmodel(d, o, s) = (haskey(d, "parserv") ? makepmodel1(d) : makepmodel2(o, s; intype=:normal))


text:
module TestMapfoldl
include("preamble.jl")

@testset "mapfoldl" begin
    xf = opcompose(Filter(isodd), Map(inc))

    ans = @test_deprecated mapfoldl(xf, +, 0:4, init = 0)
    DEPWARN_ERROR || @test ans == 6

    ans = @test_deprecated mapfoldl(xf, +, eduction(x - 1 for x in 1:5), init = 0)
    DEPWARN_ERROR || @test ans == 6
end

@testset "identityof error" begin
    if !DEPWARN_ERROR
        err = @test_deprecated @test_error mapfoldl(Map(identity), +, Any[])
        @test_broken err isa EmptyResultError

        err = @test_deprecated @test_error mapfoldl(Map(identity), +, Vector{Int}[])
        @test_broken err isa EmptyResultError
    end
end

end  # module


text:
module Fuzzy

	export TriangularMF, GaussianMF, BellMF, TrapezoidalMF, SigmoidMF

	export Rule, FISMamdani, FISSugeno

	export eval_fis

	export minimum_value, algebraic_product, bounded_difference, drastic_product, einstein_product, hamacher_product

	export maximum_value, algebraic_sum, bounded_sum, drastic_sum, einstein_sum, hamacher_sum

	# T-Norm
	include("T-norm.jl")

	# S-Norm
	include("S-norm.jl")

	# Membership functions
	include("MF.jl")

	# FIS
	include("FIS.jl")

	# Evaluation functions
	include("Eval.jl")

end


text:
include("../examples/lorenz.jl")
using JLD
using LinearAlgebra
function p_y_g_x(a, σ_o)
        term =  log(1/sqrt(2π)/σ_o)
        return sum(term .- 0.5*a.*a/σ_o/σ_o)
end
function obs(x, σ_o)
    return x .+ σ_o*randn(d)
end
function obs1(x, σ_o)
	return [x[3] + σ_o*randn()]
end

function resample(x, w)
    # multinomial resampling
    cdfw = cumsum(w)
    new_pts = similar(x)
	N_p = size(w)[1]
    for j = 1:N_p
        r = rand()
        for i = 1:N_p
            if cdfw[i] >= r
                new_pts[:,j] .= x[:,i]
                break
            end
        end
    end
    return new_pts
end

function sir(y,x_t,x_ip,Δ,obs_fun,σ_o,σ_d,N_thr=10)
    K = size(y)[2]
    N_p = size(x_ip)[2]

    N_y = K*Δ
	@show N_y
    rmse_trj = zeros(K)
    N_runup = 3000

	w = SharedArray{Float64}(N_p)
	w .= 1.0/N_p
	x = SharedArray{Float64}(d,N_p)
	t = @distributed for i=1:N_p
		x[:,i] = rand(d)
		for k = 1:N_runup
			x[:,i] = next(x[:,i],s)
		end
	end
	wait(t)

	mean_x = zeros(d,K)
    j = 1
	count = 0
    for k = 1:N_y
        if k % Δ == 0
			filter_mean = sum(x,dims=2)/N_p
			rmse_trj[j] = norm(filter_mean - x_t[:,k])/sqrt(d)

			t = @distributed for i = 1: N_p
                logwi = log(w[i]) + p_y_g_x(y[:,j] .- obs_fun(x[:,i],σ_o), σ_o)
                w[i] = exp(logwi)
			end
            w ./= sum(w)
            N_eff = 1.0/sum(w.*w)
            if (N_eff < N_thr)
				count = count + 1
			    x .= resample(x, w)
		    end
			j = j+1
		end
        for i = 1:N_p
            x[:,i] .= next(x[:,i],σ_d)
        end
    end
	@show "Resampling was triggered ", count, " times out of ", N_y
    return rmse_trj
end
"""
    assimilate(K, Np, σ_o, σ_d,
				   Δ, Nth)

Perform `K` data assimilation steps using an SIR algorithm with `Np` particles. Other inputs:

    1. `σ_o`: std of Gaussian observation noise
	2. `σ_d`: std of Gaussian dynamics noise, which is added to every component and at every timestep.
	3. `Δ`: inter-observation number of timesteps, e.g., if `Δ = 5,` observations are assumed available at timestep `0,5,10,...,5(K-1)`.
	4. `Nth`: number of particles below which to resample
	5. `obs`: observation map, e.g., `obs(x) = x[3]`

Output:

`mean_rmse`: Summary scalar performance metric:

# Examples
```julia-repl
julia> rmse_avg = assimilate(500, 1000, 0.1, 0.1, 1, 10, obs)
```
"""
function assimilate(K, Np, σ_o, σ_d,
				   Δ, Nth, obsfun)
	Ny = K*Δ
	ytest = obsfun(rand(d),σ_o)
	dy = size(ytest)[1]
    x_true = ones(d,Ny)
    x0_true = rand(d)
    Nrunup = 2000
    for i = 1:Nrunup
        x0_true = next(x0_true, 0.0)
    end
    x_true[:,1] .= x0_true
    y = zeros(dy, K)
	y[:,1] .= obsfun(x0_true, σ_o)
	k = 1
    for i = 2:Ny
        x_true[:,i] = next(x_true[:,i-1], σ_d)
		if i % Δ == 0
        	y[:,k] .= obsfun(x_true[:,i], σ_o)
			k = k+1
		end
    end
    # store trajectory of w and x.
    rmse_trj = sir(y, x_true, Δ, obsfun, σ_o, σ_d, Nth)
	mean_rmse = sum(rmse_trj)/K
    return mean_rmse
end


text:
# TODO: Build IPython 1.0 dependency? (wait for release?)

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

# print to stderr, since that is where Pkg prints its messages
eprintln(x...) = println(STDERR, x...)

# on OSX, make sure PATH changes in profile are used in detecting ipython
const ipython = @osx ? chomp(readall(`bash -lc "which ipython"`)) : "ipython"

ipyvers = try
    convert(VersionNumber, chomp(readall(`$ipython --version`)))
catch e
    error("IPython is required for IJulia, got error $e")
end

if ipyvers < v"1.0.0-dev"
    error("IPython 1.0 or later is required for IJulia, got $ipyvers instead")
else
    eprintln("Found IPython version $ipyvers ... ok.")
end

#######################################################################
# Create Julia profile for IPython and fix the config options.

# create julia profile (no-op if we already have one)
eprintln("Creating julia profile in IPython...")
run(`$ipython profile create julia`)

juliaprof = chomp(readall(`$ipython locate profile julia`))

# set c.$s in prof file to val, or nothing if it is already set
# unless overwrite is true
function add_config(prof::String, s::String, val, overwrite=false)
    p = joinpath(juliaprof, prof)
    r = Regex(string("^[ \\t]*c\\.", replace(s, r"\.", "\\."), "\\s*=.*\$"), "m")
    if isfile(p)
        c = readall(p)
        if ismatch(r, c)
            m = replace(match(r, c).match, r"\s*$", "")
            if !overwrite || m[search(m,'c'):end] == "c.$s = $val"
                eprintln("(Existing $s setting in $prof is untouched.)")
            else
                eprintln("Changing $s to $val in $prof...")
                open(p, "w") do f
                    print(f, replace(c, r, old -> "# $old"))
                    print(f, """
c.$s = $val
""")
                end
            end
        else
            eprintln("Adding $s = $val to $prof...")
            open(p, "a") do f
                print(f, """

c.$s = $val
""")
            end
        end
    else
        eprintln("Creating $prof with $s = $val...")
        open(p, "w") do f
            print(f, """
c = get_config()
c.$s = $val
""")
        end
    end
end

# add Julia kernel manager if we don't have one yet
add_config("ipython_config.py", "KernelManager.kernel_cmd",
           """["$(escape_string(joinpath(JULIA_HOME,(@windows?"julia.bat":"julia-basic"))))", "$(escape_string(joinpath(Pkg.dir("IJulia"),"src","kernel.jl")))", "{connection_file}"]""",
           true)

# make qtconsole require shift-enter to complete input
add_config("ipython_qtconsole_config.py",
           "IPythonWidget.execute_on_complete_input", "False")

# set Julia notebook to use a different port than IPython's 8888 by default
add_config("ipython_notebook_config.py", "NotebookApp.port", 8998)

#######################################################################
# Copying files into the correct paths in the profile lets us override
# the files of the same name in IPython.

rb(filename::String) = open(readbytes, filename)
eqb(a::Vector{Uint8}, b::Vector{Uint8}) =
    length(a) == length(b) && all(a .== b)

# copy IJulia/deps/src to destpath/destname if it doesn't
# already exist at the destination, or if it has changed (if overwrite=true).
function copy_config(src::String, destpath::String,
                     destname::String=src, overwrite=true)
    dest = joinpath(destpath, destname)
    srcbytes = rb(joinpath(Pkg.dir("IJulia"), "deps", src))
    if !isfile(dest) || (overwrite && !eqb(srcbytes, rb(dest)))
        eprintln("Copying $src to Julia IPython profile.")
        open(dest, "w") do f
            write(f, srcbytes)
        end
    else
        eprintln("(Existing $destname file untouched.)")
    end
end

# copy IJulia icon to profile so that IPython will use it
mkpath(joinpath(juliaprof, "static", "base", "images"))
for T in ("png", "svg")
    copy_config("ijulialogo.$T",
                joinpath(juliaprof, "static", "base", "images"),
                "ipynblogo.$T")
end

# Use our own version of tooltip to handle identifiers ending with !
# (except for line 211, tooltip.js is identical to the IPython version)
# IPython might make his configurable later, at which point the logic
# should be moved to custom.js or a config file.
mkpath(joinpath(juliaprof, "static", "notebook", "js"))
copy_config("tooltip.js", joinpath(juliaprof, "static", "notebook", "js"))

# custom.js can contain custom js login that will be loaded
# with the notebook to add info and/or monkey-patch some javascript
# -- e.g. we use it to add .ipynb metadata that this is a Julia notebook
mkpath(joinpath(juliaprof, "static", "custom"))
copy_config("custom.js", joinpath(juliaprof, "static", "custom"))

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


text:
module GenieAdmin

function install(dest::String; force = false) :: Nothing
    src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

    for f in readdir(src)
        isfile(f) && continue
        isdir(f) || mkpath(joinpath(src, f))

        Genie.Plugins.install(joinpath(src, f), dest, force = force)
    end

    nothing
end

end


text:
module MieSeries

export PECSphere
export excite
export farfield
export scatteredfield

using SpecialFunctions

"""
    Ylm(K, r)

Computes the spherical harmonics up to order `K-1` in direction `r` .
"""
function Ylm(K, r)

    r = r / norm(r)
    x = r[1]; y = r[2]; z = r[3]

    T = eltype(r)
    Y = zeros(Complex{T}, K*K)

    Y[1] = 1/√4π

    facphip =  x + im*y;
    facphim = -x + im*y;
    for m ∈ 1:(K-1)
        Y[(m+1)*(m+1)] = -Y[m*m] * facphip * sqrt((m+0.5)/m)
        Y[m*m+1] = -Y[(m-1)*(m-1)+1] * facphim * sqrt((m+0.5)/m)
    end

    Y[3] = √3 * z * Y[1]
    for l ∈ 2:(K-1)
        Y[(l+1)*(l+1)-1] = sqrt(2l+1) * z * Y[l*l]
        Y[l*l+2] = sqrt(2*l+1) * z * Y[(l-1)*(l-1)+1]
    end

    for l ∈ 2:(K-1)
        for m ∈ (2-l):(l-2)
            Y[l*(l+1)+m+1] =
                sqrt((4*l^2-1)/(l^2-m^2)) * z * Y[l*(l-1)+m+1] -
                sqrt(((2*l+1)*((l-1)^2-m^2))/((l^2-m^2)*(2*l-3))) * Y[(l-1)*(l-2)+m+1]
        end
    end

    Y
end

function Xlm(K, r)

    r = r / norm(r)
    x = r[1]; y = r[2]; z = r[3];

    Y = Ylm(K, r);
    Xx = zeros(Y); Xy = zeros(Y); Xz = zeros(Y);
    Yx = zeros(Y); Yy = zeros(Y); Yz = zeros(Y);
    Zx = zeros(Y); Zy = zeros(Y); Zz = zeros(Y);

    for l = 1:(K-1)
        Ymin = sqrt(2*l) * Y[(l+1)*(l+1)-1];
        Xx[l*(l+1) + l + 1] = 1/sqrt(4*l*(l+1)) * Ymin;
        Xy[l*(l+1) + l + 1] = im/sqrt(4*l*(l+1)) * Ymin;
        Xz[l*(l+1) + l + 1] = l * Y[(l+1)*(l+1)] / sqrt(l*(l+1));

        # m = -l
        Yplus = sqrt(2*l) * Y[l*l + 2];
        Xx[l*(l+1) - l + 1] = 1/sqrt(4*l*(l+1)) * Yplus;
        Xy[l*(l+1) - l + 1] =-im/sqrt(4*l*(l+1)) * Yplus;
        Xz[l*(l+1) - l + 1] = -l * Y[l*l + 1] / sqrt(l*(l+1));

        # other values of m
        for m = (1-l):(l-1)
            Yplus = sqrt((l-m)*(l+m+1)) * Y[l*(l+1) + m+1 + 1];
            Ymin  = sqrt((l+m)*(l-m+1)) * Y[l*(l+1) + m-1 + 1];
            Xx[l*(l+1) + m + 1] = 1/sqrt(4*l*(l+1)) * (Yplus + Ymin);
            Xy[l*(l+1) + m + 1] =-im/sqrt(4*l*(l+1)) * (Yplus - Ymin);
            Xz[l*(l+1) + m + 1] = m * Y[l*(l+1) + m + 1] / sqrt(l*(l+1));
        end
    end

    # Y
    Yx = Y * x;
    Yy = Y * y;
    Yz = Y * z;

    # Z aka n × X
    Zx = y * Xz - z * Xy
    Zy = z * Xx - x * Xz
    Zz = x * Xy - y * Xx

    Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz
end

"""
    MNlm(numterms, position, wavenumber, besselfn)

Compute the multipole solutions to Maxwell's equations up to the orbital angular
momentum `numterms`. The argument `besselfn` can be one of `sbessel` or `shankel`.

The convention is such that the temporal variation of a fasor is

```math
a(t) = Re(a*e^{-iωt})
```

In particular, the Hankel function of the first kind corresponds to outgoing waves.

The multipole solutions are defined by

```math
    M_{l,m}(p) = besself(kr) X_{lm}(Θ,ϕ)
    N_{l,m}(p) = \frac{1}{k} ∇ × (besself(kr) X_{lm}(Θ,ϕ))
```

with ``X_{l,m}`` the vector spherical harmonics as defined in [Jackson, Classical Electrodynamics],
i.e. normalised in ``L^2(S)"`` with ``S`` the unit sphere.
"""
function MNlm(K, r, k, f)

    # Compute the spherical harmonics
    Xx, Xy, Xz, Yx, Yy, Yz, Px, Py, Pz = Xlm(K, r)

    # Argument van de besselfuncties (ook fysisch houden) en golfgetal:
    nor = k*norm(r)
    r = k * r

    # De functies invullen
    C = complex(eltype(r))
    Mx = zeros(C, K*K)
    My = zeros(C, K*K)
    Mz = zeros(C, K*K)
    Nx = zeros(C, K*K)
    Ny = zeros(C, K*K)
    Nz = zeros(C, K*K)

    for l = 1:(K-1)
        bessel = f(l, nor);
        besseldiff = ((l+1)*f(l-1, nor) - l*f(l+1, nor)) / (2*l+1);

        for m = -l:l
            lm = l*(l+1)+m+1
        	Mx[lm] = bessel * Xx[lm]
            My[lm] = bessel * Xy[lm]
            Mz[lm] = bessel * Xz[lm]

            factor = im * sqrt(l*(l+1)) * bessel / nor
            Nx[lm] = factor * Yx[lm] + besseldiff * Px[lm]
            Ny[lm] = factor * Yy[lm] + besseldiff * Py[lm]
            Nz[lm] = factor * Yz[lm] + besseldiff * Pz[lm]
        end
    end
    Mx, My, Mz, Nx, Ny, Nz
end

type PECSphere
    radius
    background_permittivity
    background_permeability
end

type Expansion
    numterms
    omega
    la1
    la2
    ra1
    ra2
end

"""
    excite(object, omega, pol)

Compute the expansion coefficients of the field traces in terms of harmonics.
The polarisation can be choosen from `:left`, `:right`, or ':x' (default).
"""
function excite(obj, omega, pol)

    r = obj.radius
    T = typeof(r)
    C = complex(T)
    K = 20

    # polarisation
    if pol == :left
        pols = [-1]
    elseif pol == :right
        pols = [+1]
    elseif pol == :x
        pols = [-1,1]
    end

    ra1 = zeros(C,K)
    ra2 = zeros(C,K)
    la1 = zeros(C,K)
    la2 = zeros(C,K)

    eps0 = obj.background_permittivity
    mu0  = obj.background_permeability
    k0 = omega * sqrt(eps0 * mu0);

    # Bessels, hankels and derivatives
    ll = collect(1:K)
    jk0   = sbessel.(ll, k0*r)
    djk0 = ((ll+1).*sbessel.(ll-1, k0*r)  - ll.*sbessel.(ll+1, k0*r))  ./ (2*ll+1)
    hk0   = shankel.(ll, k0*r)
    dhk0 = ((ll+1).*shankel.(ll-1, k0*r)  - ll.*shankel.(ll+1, k0*r))  ./ (2*ll+1)

    # compute the expansion coefficients for all requested polarisations
    for p in pols
        for l in ll
            a1 = (1*(-im)^l * sqrt(4*pi*(2*l+1)) * jk0[l]) / (-hk0[l])
            a2 = (p*(-im)^l * sqrt(4*pi*(2*l+1)) * djk0[l])/(-dhk0[l])
            # rescale for linear polarisation
            if pol == :x
                a1 = a1/2
                a2 = a2/2
            end
            # assign
            if p == -1
                la1[l] = a1
                la2[l] = a2
            elseif p == +1
                ra1[l] = a1
                ra2[l] = a2
            end
        end
    end

    Expansion(K, omega, la1, la2, ra1, ra2)
end

function farfield(obj, expansion, punten)

    # allocate space
    T = eltype(punten)
    C = complex(T)

    Ex = zero(C)
    Ey = zero(C)
    Ez = zero(C)

    Hx = zero(C)
    Hy = zero(C)
    Hz = zero(C)

    # basisfuncties in posities punten
    Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz = Xlm(expansion.numterms+1, punten)

    eps0 = obj.background_permittivity
    mu0  = obj.background_permeability
    k = expansion.omega * sqrt(eps0 * mu0);
    #k = expansion.omega * sqrt(obj.eps0*obj.mu0)
    Z = sqrt(mu0 / eps0)

    # build fields
    for p = [-1, 1]
        if p == -1
            a1 = expansion.la1
            a2 = expansion.la2
        else
            a1 = expansion.ra1
            a2 = expansion.ra2
        end
        m = -p
        for l in 1:expansion.numterms

            h =  im^(l+1) / k # shankel(l, k*r);
            dh = h * (-im)    # ((l+1).*shankel(l-1, k*r)  - l.*shankel(l+1, k*r))  ./ (2*l+1);

            # contributie in X_{l,m}
            Ex += a1[l] * h .* Xx[l*(l+1) + m + 1]
            Ey += a1[l] * h .* Xy[l*(l+1) + m + 1]
            Ez += a1[l] * h .* Xz[l*(l+1) + m + 1]
            Hx += im / Z * a2[l] * h .* Xx[l*(l+1) + m + 1]
            Hy += im / Z * a2[l] * h .* Xy[l*(l+1) + m + 1]
            Hz += im / Z * a2[l] * h .* Xz[l*(l+1) + m + 1]

            # contributie in Z_{l,m}
            Ex += a2[l] * dh .* Zx[l*(l+1) + m + 1]
            Ey += a2[l] * dh .* Zy[l*(l+1) + m + 1]
            Ez += a2[l] * dh .* Zz[l*(l+1) + m + 1]
            Hx += im/Z * a1[l] * dh .* Zx[l*(l+1) + m + 1]
            Hy += im/Z * a1[l] * dh .* Zy[l*(l+1) + m + 1]
            Hz += im/Z * a1[l] * dh .* Zz[l*(l+1) + m + 1]

        end
    end

    Ex, Ey, Ez, Hx, Hy, Hz
end

function scatteredfield(obj, expansion, punten)

    # outputarrays
    T = eltype(punten)
    C = complex(T)

    Ex = zero(C)
    Ey = zero(C)
    Ez = zero(C)

    Hx = zero(C)
    Hy = zero(C)
    Hz = zero(C)

    eps0 = obj.background_permittivity
    mu0  = obj.background_permeability
    k = expansion.omega * sqrt(eps0*mu0)
    Z = sqrt(mu0 / eps0)

    # multipole solutions of Maxwell's equations
    Mx, My, Mz, Nx, Ny, Nz = MNlm(expansion.numterms+1, punten, k, shankel)

    # radius
    r = norm(punten)

    # build fields
    for p in [-1, 1]
        if p == -1
            A1 = expansion.la1
            A2 = expansion.la2
        else
            A1 = expansion.ra1
            A2 = expansion.ra2
        end
        m = -p
        for l in 1:expansion.numterms

            # linear index
            lm = l*(l+1)+m+1

            h  = shankel(l, k*r)
            dh = ((l+1)*shankel(l-1, k*r)  - l*shankel(l+1, k*r))  / (2*l+1)

            a1 = A1[l]
            a2 = A2[l]

            # contribution in M_{l,m}
            Ex += a1 * Mx[lm]
            Ey += a1 * My[lm]
            Ez += a1 * Mz[lm]
            Hx += a2 * im / Z * Mx[lm]
            Hy += a2 * im / Z * My[lm]
            Hz += a2 * im / Z * Mz[lm]

            # contributie in N_{l,m}
            Ex += a2 * Nx[lm]
            Ey += a2 * Ny[lm]
            Ez += a2 * Nz[lm]
            Hx += a1 * im / Z * Mx[lm]
            Hy += a1 * im / Z * My[lm]
            Hz += a1 * im / Z * Mz[lm]

        end
    end

    Ex, Ey, Ez, Hx, Hy, Hz
end

"""
    sbessel(nu, x)

returns the spherical Bessel functions jnu(x).
"""
sbessel(nu, x) = sqrt(pi ./(2x)) .* besselj(nu + 0.5, x)

"""
    shankel(nu, x)

Computes the spherical hankel function (of the first kind).
"""
shankel(nu, x) = sqrt(pi ./(2x)) .* besselh(nu + 0.5, x)

end # module


text:
"""
    TPI.Bcast!(data, root, comm)

Broadcasts the content of `data` on thread `root` to the other threads.
"""
function Bcast!(data::A, root, comm::Comm) where {T, A<:DenseArray{T}}
    sz = size(data)

    if Comm_rank(comm) == root
        # pass to the other threads the dst (reduced array)
        comm.passthrough[root] = convert(UInt, pointer(data))

        sync(comm.barrier)
    else
        # wait for thread 1 to compute the sum
        sync(comm.barrier)

        ptr = convert(Ptr{T}, comm.passthrough[root])
        arr = UnsafeArray(ptr, sz)
        data .= arr
    end
    return data
end


text:
export  vdm_jac, vdm_iter!, vdm_newton

#----------------------------------------------------------------------
function vdm_jac(w, Xi, L)
    n = size(Xi,1)
    r = size(Xi,2)
    J = fill(typeof(Xi[1])(0), length(L), r + r*n)

    for k in 1:r
        for i in 1:length(L)
            J[i,k] = 1.;
            for v in 1:n
                J[i,k]*= Xi[v,k]^((L[i].z)[v])
            end
            for l in 1:n
                J[i,r+(k-1)*n+l] = w[k];
                for v in 1:n
                    if v==l
                        if (L[i].z)[v] > 0
                            J[i,r+(k-1)*n+l] *= (L[i].z)[v]*Xi[v,k]^((L[i].z)[v]-1)
                        else
                            J[i,r+(k-1)*n+l] = 0.0
                        end
                    else
                        J[i,r+(k-1)*n+l] *= Xi[v,k]^((L[i].z)[v])
                    end
                end
            end
        end
    end
    J
end

#----------------------------------------------------------------------
function vdm_iter!(w, Xi, s0, L)
    f = fill(typeof(Xi[1])(0), length(L))
    mnt = moment(w,Xi)
    for i in 1:length(L)
        f[i] = mnt(L[i].z)-s0[L[i]]
    end
    print("|ε|=", norm(f))
    J = vdm_jac(w,Xi,L)
    dwXi = J\f
    delta = norm(dwXi)
    println("  δ=", delta)

    r = length(w)
    n = size(Xi,1)
    for i in 1:r
        w[i] -= dwXi[i]
    end
    for i in 1:r
        for j in 1:n
            Xi[j,i] -= dwXi[r+(i-1)*n+j]
        end
    end

    delta
end

#----------------------------------------------------------------------
function vdm_newton(w0, Xi0, s0, L;  args...)
    eps::Float64 = 1.e-5
    cmax::Int64  = 10
    for arg in args
        if arg[1]==:maxit
            cmax=arg[2]
        end
        if arg[1]==:eps
            eps=arg[2]
        end
    end

    w = copy(w0)
    Xi = copy(Xi0)
    c = 0; delta = 1.0
    while delta > eps && c< cmax
        delta = vdm_iter!(w,Xi,s0,L)
        c += 1;
    end
    w, Xi
end

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


text:
function pyramide(r::Robot)
    move_angle_left!(r)
    lenght = moves!(r, Ost)
    hight = moves!(r, Nord)
    move_angle_left!(r)
    j=lenght+1
        for i in 1:hight
            while j-i != 0
                putmarker!(r)
                move!(r,Ost)
                j-=1
            end
            putmarker!(r)
            j=lenght+1
            back(r)
            move!(r,Nord)
        end
end

function move_angle_left!(r::Robot) #перемещает в левый нижний угол
    for side in (Sud, West)
        while isborder(r, side) == false
            move!(r, side)
        end
    end
end

function moves!(r::Robot,side::HorizonSide) #возвращает кол-во шагов в опред направлении
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end


text:
# Tolerance
epsilon = 0.0001

"""
Function used to transform a column with numerical values into one or several binary columns.

Arguments:
 - data: table which contains the column that will be binarized (1 row = 1 individual, 1 column = 1 feature);
 - header: header of the column of data that will be binarized
 - intervals: array of values which delimits the binarization (ex : [2, 4, 6, 8] will lead to 3 columns respectively equal to 1 if the value of column "header" is in [2, 3], [4, 5] and [6, 7])

 Example:
  createColumns(:Age, [1, 17, 50, Inf], data, features) will create 3 binary columns in features named "Age1-16", "Age17-49", "Age50-Inf"
"""
function createColumns(header::Symbol, intervals, data::DataFrames.DataFrame, features::DataFrames.DataFrame)
    for i in 1:size(intervals, 1) - 1
        lb = intervals[i]
        ub = intervals[i+1]
        features[!, Symbol(header, lb, "-", (ub-1))] = ifelse.((data[!, header] .>= lb) .& (data[!, header] .< ub), 1, 0)
    end
end

"""
Create the train and test csv files of a data set

Arguments:
 - dataFolder: folder which contains the data set csv file (ex: "./data")
 - dataSet: name of the data set (ex: "titanic")

Important remark: the first column of the output files must correspond to the class of each individual (and it must be 0 or 1)
"""
function createFeatures(dataFolder::String, dataSet::String)

    # Get the input file path
    rawDataPath = dataFolder * dataSet * ".csv"

    # Test its existence
    if !isfile(rawDataPath)
        println("Error in createFeatures: fichier non trouve: ", rawDataPath)
        return
    end

    # Put the data in a DataFrame variable
    rawData = CSV.read(rawDataPath,  header=true)

    # Output files path
    trainDataPath = dataFolder * dataSet * "_train.csv"
    testDataPath = dataFolder * dataSet * "_test.csv"

    # If the train or the test file do not exist
    if !isfile(trainDataPath) || !isfile(testDataPath)

        println("=== Creating the features")

        # Create the table that will contain the features
        features::DataFrame = DataFrames.DataFrame()

        # Create the features of the titanic data set
        if dataSet == "titanic"

            # Add the column related to the class (always do it first!)
            # Remark: here the rawData already contain 0/1 values so no change needed
            features.Survived = rawData.Survived

            #### First example of the binarization of a column with numerical values
            # Add columns related to the ages
            # -> 3 columns ([0, 16], [17, 49], [50, +infinity[)
            # Ex : if the age is 20, the value of these 3 columns will be 0, 1 and 0, respectively.
            createColumns(:Age, [0, 17, 50, Inf], rawData, features)

            # Add columns related to the fares
            # -> 3 columns ([0, 9], [10, 19], [20, +infinity[)
            createColumns(:Fare,  [0, 10, 20, Inf], rawData, features)

            #### First example of the binarization of a column with categorical values
            # Add 1 column for the sex (female or not)
            # Detailed description of the command:
            # - create in DataFrame "features" a column named "Sex"
            # - for each row of index i of "rawData", if column "Sex" is equal to "female", set the value of column "Sex" in row i of features to 1; otherwise set it to 0
            features.Sex = ifelse.(rawData.Sex .== "female", 1, 0)

            # Add columns related to the passenger class
            # -> 3 columns (class 1, class 2 and class 3)

            # For each existing value in the column "Pclass"
            for a in sort(unique(rawData.Pclass))

                # Create 1 feature column named "Class1", "Class2" or "Class3"
                features[!, Symbol("Class", a)] = ifelse.(rawData.Pclass .<= a, 1, 0)
            end

            # Add a column related  to the number of relatives
            # -> 1 column (0: no relatives, 1: at least one relative)
            features.Relative = ifelse.(rawData[!, Symbol("Siblings/Spouses Aboard")] + rawData[!, Symbol("Parents/Children Aboard")] .> 0, 1, 0)
        end

        if dataSet == "kidney"
            # add class
            features.class = ifelse.(rawData.class .== "ckd", 1, 0)

            #createColumns(:age, [0, 17, 50, Inf], rawData, features)
            createColumns(:age, [0, 45, Inf], rawData, features)
            #createColumns(:bp, [0, 80, Inf], rawData, features)

            features.sg = ifelse.(rawData.sg .==1.020 , 1, 0) .+ ifelse.(rawData.sg .==1.025 , 1, 0)

            #features.al = ifelse.(rawData.al .== 0, 1, 0)

            #features.bgr = ifelse.(rawData.bgr .>= 160, 1, 0)

            #features.pc = ifelse.(rawData.pc .== "normal", 1, 0)

        end

        if dataSet == "flowers"
            features.class = ifelse.(rawData.class .=="Setosa" , 1, 0) .+ ifelse.(rawData.class .=="Versicolor" , 2, 0)
            #createColumns(:sepalLength, [0, 6, Inf], rawData, features)
            #createColumns(:sepalWidth, [0, 3.1, Inf], rawData, features)
            #createColumns(:petalLength, [0, 2.1,5.1, Inf], rawData, features)
            #createColumns(:petalWidth, [0, 1,1.8, Inf], rawData, features)
            createColumns(:sL, [0, 6, Inf], rawData, features)
            createColumns(:sW, [0, 3.1, Inf], rawData, features)
            createColumns(:pL, [0, 2.3,5.1, Inf], rawData, features)
            createColumns(:pW, [0, 0.7,1.8, Inf], rawData, features)

        end

        # Shuffle the individuals
        features = features[shuffle(1:size(features, 1)),:]

        # Split them between train and test
        trainLimit = trunc(Int, size(features, 1) * 2/3)
        train = features[1:trainLimit, :]
        test = features[(trainLimit+1):end, :]

        #CSV.write(trainDataPath, train)
        #CSV.write(testDataPath, test)

    # If the train and test file already exist
    else
        println("=== Warning: Existing features found, features creation skipped")
        println("=== Loading existing features")
        train = CSV.read(trainDataPath)
        test = CSV.read(testDataPath)
    end

    println("=== ... ", size(train, 1), " individuals in the train set")
    println("=== ... ", size(test, 1), " individuals in the test set")
    println("=== ... ", size(train, 2), " features")

    return train, test
end

function InitializeModel(cmax::Int64, RgenX::Float64, RgenB::Float64, n::Int64, d::Int64, S::Set{Int64}, t::DataFrame)
    m = Model(with_optimizer(CPLEX.Optimizer))
    set_parameter(m,"CPX_PARAM_SCRIND",0)

    @variable(m, 0 <= x[1:n] <= 1)
    @variable(m, b[1:d], Bin)

    @objective(m, Max, sum(x[i] for i in S) - RgenX*sum(x[i] for i in 1:n) - RgenB*sum(b[j] for j in 1:d))

    @constraint(m, ct1[i in 1:n, j in 1:d], x[i] <= 1 + (t[i,j] - 1)*b[j])
    @constraint(m, ct2[i in 1:n, j in 1:d], x[i] >= 1 + sum((t[i,j] - 1)*b[j] for j in 1:d))
    @constraint(m, ct3, sum(x[i] for i in 1:n) <= cmax)
    return m,x,b
end

function P(model::AbstractModel, S::Set{Int64},x,b)
    #optimize!(model)
    value = objective_value(model)
    #print("la valeur du modele est $value")
    xSolution=JuMP.value.(x)
    bSolution=JuMP.value.(b)
    #bSolution = getvalue.(b)
    support = sum(xSolution[i] for i in S)
    #println("bSolution : ", bSolution)
    #println("support : $support")
    return support, bSolution
end

function createRules(dataSet::String, resultsFolder::String, train::DataFrames.DataFrame)

    # Output file
    rulesPath = resultsFolder * dataSet * "_rules.csv"
    rules::DataFrame = similar(train,0)

    if !isfile(rulesPath)

        println("=== Generating the rules")

        # Transactions
        t::DataFrame = train[:, 2:end]

        # Class of the transactions
        # Help: to get the class of transaction nb i:
        # - do not use: transactionClass[i]
        # - use: transactionClass[i, 1]
        transactionClass::DataFrame = train[:, 1:1]

        # Number of features
        d::Int64 = size(t, 2)

        # Number of transactions
        n::Int64 = size(t, 1)

        #mincovy::Float64 = 0.90
        #97 pour test rapide, 95 test lent
        mincovy::Float64 = 0.95
        iter_lim::Int64 = 5
        RgenX::Float64 = 0.1 / n
        RgenB::Float64 = 0.1 / (n * d)

        ##################
        # Find the rules for each class
        ##################

        listeClasses=[]
        for i in 1:n
            push!(listeClasses,transactionClass[i,1])
        end
        listeClasses= unique(listeClasses)
        #for y = 0:1
        #for y = 0:2
        for y in listeClasses
            println("-- generating rules for class $y")
            S = Set{Int64}() # transactions of class y
            for i in 1:n
                if (transactionClass[i,1] == y)
                    push!(S, i)
                end
            end

            support=n+1
            iter = 1
            cmax = n
            #cmax=Int(n)
            #println("cmax=$cmax")
            modelY,x,b=InitializeModel(cmax, RgenX, RgenB, n, d, S, t)

            while (cmax >= n*mincovy)
                if iter<iter_lim
                    optimize!(modelY)
                    if termination_status(modelY) == MOI.OPTIMAL
                        supportTemp, rule =P(modelY,S,x,b)
                        if iter == 1
                            support = supportTemp
                        end
                        if supportTemp >= support
                            push!(rules,append!([y],rule))
                            @constraint(modelY,  sum(b[j] for j in 1:d if rule[j]==0)+sum(1-b[j] for j in 1:d if rule[j]==1) >= 1)
                            iter += 1
                        else
                            cmax -= 1
                            iter = 1
                        end
                    else
                        break
                    end
                else
                    cmax -= 1
                    iter = 1
                end
            end

            # Help: Let rule be a rule that you want to add to rules
            # - if it is the first rule, use: rules = rule
            # - if it is not the first rule, use: rules = append!(rules, rule)
        end

        #CSV.write(rulesPath, rules)

    else
        println("=== Warning: Existing rules found, rules creation skipped")
        println("=== Loading the existing rules")
        rules = CSV.read(rulesPath)
    end

    println("========Exititng rule generation Algorithm========")
    return rules
end

"""
Sort the rules

Arguments
  - dataSet: name of the dataset folder
  - resultsFolder: name of the folder in which the results are written
  - train: train data set (1 row = 1 individual)
  - rules: rules which must be sorted (1 row = 1 rule)
  - tilim: maximal running time of CPLEX in seconds
"""
function sortRules(dataSet::String, resultsFolder::String, train::DataFrames.DataFrame, rules::DataFrames.DataFrame, tilim::Int64)

    orderedRulesPath = resultsFolder * dataSet * "_ordered_rules.csv"

    if !isfile(orderedRulesPath)

        println("=== Sorting the rules")

        # Transactions
        t = train[:, 2:end]
        # Class of the transactions
        transactionClass = train[:, 1:1]
        # Number of features
        d = size(t, 2)
        # Number of transactions
        n = size(t, 1)

        # Add the two null rules in first position
        nullrules = similar(rules, 0)
        push!(nullrules, append!([0], zeros(d)))
        push!(nullrules, append!([1], zeros(d)))
        rules = vcat(nullrules, rules)

        # Remove duplicated rules
        rules = unique(rules)

        # Number of rules
        L = size(rules)[1]

        Rrank = 1/L

        ################
        # Compute the v_il and p_il constants
        # p_il = :
        #  0 if rule l does not apply to transaction i
        #  1 if rule l applies to transaction i and   correctly classifies it
        # -1 if rule l applies to transaction i and incorrectly classifies it
        ################
        p = zeros(n, L)

        # For each transaction and each rule
        for i in 1:n
            for l in 1:L

                # If rule l applies to transaction i
                # i.e., if the vector t_i - r_l does not contain any negative value
                if !any(x->(x<-epsilon), [sum(t[i, k]-rules[l, k+1]) for k in 1:d])

                    # If rule l correctly classifies transaction i
                    if transactionClass[i, 1] == rules[l, 1]
                        p[i, l] = 1
                    else
                        p[i, l] = -1
                    end
                end
            end
        end

        v = abs.(p)

        ################
        # Create and solve the model
        ###############
        m = Model(with_optimizer(CPLEX.Optimizer))
        set_parameter(m, "CPX_PARAM_TILIM", tilim)

        # u_il: rule l is the highest which applies to transaction i
        @variable(m, u[1:n, 1:L], Bin)

        # r_l: rank of rule l
        @variable(m, 1 <= r[1:L] <= L, Int)

        # rstar: rank of the highest null rule
        @variable(m, 1 <= rstar <= L)
        @variable(m, 1 <= rB <= L)

        # g_i: rank of the highest rule which applies to transaction i
        @variable(m, 1 <= g[1:n] <= L, Int)

        # s_lk: rule l is assigned to rank k
        @variable(m, s[1:L,1:L], Bin)

        # Rank of null rules
        rA = r[1]
        rB = r[2]

        # rstar == rB?
        @variable(m, alpha, Bin)

        # rstar == rA?
        @variable(m, 0 <= beta <= 1)

        # Maximize the classification accuracy
        @objective(m, Max, sum(p[i, l] * u[i, l] for i in 1:n for l in 1:L)
                   + Rrank * rstar)

        # Only one rule is the highest which applies to transaction i
        @constraint(m, [i in 1:n], sum(u[i, l] for l in 1:L) == 1)

        # g constraints
        @constraint(m, [i in 1:n, l in 1:L], g[i] >= v[i, l] * r[l])
        @constraint(m, [i in 1:n, l in 1:L], g[i] <= v[i, l] * r[l] + L * (1 - u[i, l]))

        # Relaxation improvement
        @constraint(m, [i in 1:n, l in 1:L], u[i, l] >= 1 - g[i] + v[i, l] * r[l])
        @constraint(m, [i in 1:n, l in 1:L], u[i, l] <= v[i, l])

        # r constraints
        @constraint(m, [k in 1:L], sum(s[l, k] for l in 1:L) == 1)
        @constraint(m, [l in 1:L], sum(s[l, k] for k in 1:L) == 1)
        @constraint(m, [l in 1:L], r[l] == sum(k * s[l, k] for k in 1:L))

        # rstar constraints
        @constraint(m, rstar >= rA)
        @constraint(m, rstar >= rB)
        @constraint(m, rstar - rA <= (L-1) * alpha)
        @constraint(m, rA - rstar <= (L-1) * alpha)
        @constraint(m, rstar - rB <= (L-1) * beta)
        @constraint(m, rB - rstar <= (L-1) * beta)
        @constraint(m, alpha + beta == 1)

        # u_il == 0 if rstar > rl (also improve relaxation)
        @constraint(m, [i in 1:n, l in 1:L], u[i, l] <= 1 - (rstar - r[l])/ (L - 1))

        optimize!(m)

        ###############
        # Write the rstar highest ranked rules and their corresponding class
        ###############

        # Number of rules kept in the classifier
        # (all the rules ranked lower than rstar are removed)
        relevantNbOfRules=L-trunc(Int, JuMP.value(rstar))+1

        # Sort the rules and their class by decreasing rank
        rulesOrder = JuMP.value.(r)
        println(rulesOrder)
        orderedRules = rules[sortperm(L.-rulesOrder), :]

        orderedRules = orderedRules[1:relevantNbOfRules, :]

        #CSV.write(orderedRulesPath, orderedRules)

    else
        println("=== Warning: Sorted rules found, sorting of the rules skipped")
        println("=== Loading the sorting rules")
        orderedRules = CSV.read(orderedRulesPath)
    end

    return orderedRules

end

"""
Compute for a given data set the precision and the recall of
- each class
- the whole data set (with and without weight for each class)

Arguments
  - orderedRules: list of rules of the classifier (1st row = 1st rule to test)
  - dataset: the data set (1 row = 1 individual)
"""
function goodShowStatistics(orderedRules::DataFrames.DataFrame, dataSet::DataFrames.DataFrame)

    # Number of transactions
    n = size(dataSet, 1)
    nbr_classes=3

    # Statistics with respect to class 0:
    # - true positive;
    # - true negative;
    # - false positive;
    # - false negative

    confusion_matrix=zeros(Int16,nbr_classes,nbr_classes)

    # For all transaction i in the data set
    for i in 1:n

        # Get the first rule satisfied by transaction i

        ruleId = findfirst(all, collect(eachrow(Array{Float64, 2}(orderedRules[:, 2:end])  .<= Array{Float64, 2}(DataFrame(dataSet[i, 2:end])))))
        # If transaction i is classified correctly (i.e., if it is a true)
        true_y=dataSet[i, 1]+1
        pred_y=orderedRules[ruleId, 1]+1
        confusion_matrix[pred_y,true_y]+=1
    end
    println(confusion_matrix)

    println("Class\tPrec.\tRecall\tSize")
    avg_precision=0
    avg_recall=0
    wavg_precision=0
    wavg_recall=0

    for c in 1:nbr_classes
        denominateur_precision=0
        denominateur_recall=0
        for iterateur in 1:nbr_classes
            denominateur_precision+=confusion_matrix[c,iterateur]
            denominateur_recall+=confusion_matrix[iterateur,c]
        end

        prec=confusion_matrix[c,c]/denominateur_precision
        rec=confusion_matrix[c,c]/denominateur_recall
        avg_precision+=prec
        avg_recall+=rec
        wavg_precision+=prec*denominateur_recall/ size(dataSet, 1)
        wavg_recall+=rec*denominateur_recall/ size(dataSet, 1)

        virgule=2
        prec=round(prec,digits=virgule)
        rec=round(rec,digits=virgule)

        println("$c\t$prec\t$rec\t$denominateur_recall")
    end
    println("avg\t", round(avg_precision/nbr_classes, digits=2), "\t", round(avg_recall/nbr_classes, digits=2))
    println("wavg\t", round(wavg_precision, digits=2), "\t", round(wavg_recall, digits=2))
    println("")

end

function badShowStatistics(orderedRules::DataFrames.DataFrame, dataSet::DataFrames.DataFrame)

    # Number of transactions
    n = size(dataSet, 1)

    # Statistics with respect to class 0:
    # - true positive;
    # - true negative;
    # - false positive;
    # - false negative
    tp::Int = 0
    fp::Int = 0
    fn::Int = 0
    tn::Int = 0

    # Number of individuals in each class
    classSize = Array{Int, 1}([0, 0, 0])

    # For all transaction i in the data set
    for i in 1:n

        # Get the first rule satisfied by transaction i
        ruleId = findfirst(all, collect(eachrow(Array{Float64, 2}(orderedRules[:, 2:end])  .<= Array{Float64, 2}(DataFrame(dataSet[i, 2:end])))))
        # If transaction i is classified correctly (i.e., if it is a true)
        if orderedRules[ruleId, 1] == dataSet[i, 1]
            # If transaction i is of class 0
            if dataSet[i, 1] == 0
                tp += 1
                classSize[1] += 1
            else
                tn += 1
                classSize[2] += 1
            end
            # If it is a negative
        else
            # If transaction i is of class 0
            if dataSet[i, 1] == 0
                fn += 1
                classSize[1] += 1
            else
                fp += 1
                classSize[2] += 1
            end
        end
    end

    precision = Array{Float64, 1}([tp / (tp+fp), tn / (tn+fn)])
    recall = Array{Float64, 1}([tp / (tp + fn), tn / (tn + fp)])

    println("Class\tPrec.\tRecall\tSize")
    println("0\t", round(precision[1], digits=2), "\t", round(recall[1], digits=2), "\t", classSize[1])
    println("1\t", round(precision[2], digits=2), "\t", round(recall[2], digits=2), "\t", classSize[2], "\n")
    println("avg\t", round((precision[1] + precision[2])/2, digits=2), "\t", round((recall[1] + recall[2])/2, digits=2))
    println("w. avg\t", round(precision[1] * classSize[1] / size(dataSet, 1) + precision[2] * classSize[2] / size(dataSet, 1), digits = 2), "\t", round(recall[1] * classSize[1] / size(dataSet, 1) + recall[2] * classSize[2] / size(dataSet, 1), digits = 2), "\n")

end


text:
const MBDValue = MDB_val
const EnvironmentFlags = Unsigned
MDBValue() = MDB_val(zero(Csize_t), C_NULL)
MDBValue(_::Nothing) = MDBValue()
MDBValue(val::String) = MDB_val(Csize_t(sizeof(val)), convert(Ptr{Cvoid},pointer(val)))
function MDBValue(val::T) where {T}
    isbitstype(T) && error("Can not wrap a $T in MDBValue. Use a $T array instead")
    val_size = sizeof(eltype(val))*length(val)
    return MDB_val(Csize_t(val_size), convert(Ptr{Cvoid},pointer(val)))
end

convert(::Type{T}, mdb_val_ref::Ref{MDB_val}) where {T} = _convert(T, mdb_val_ref[])
function _convert(::Type{String}, mdb_val::MDB_val)
    unsafe_string(convert(Ptr{UInt8}, mdb_val.mv_data), mdb_val.mv_size)
end
function _convert(::Type{Vector{T}}, mdb_val::MDB_val) where {T}
    res = unsafe_wrap(Array, convert(Ptr{UInt8}, mdb_val.mv_data), mdb_val.mv_size)
    reinterpret(T,res)
end
function _convert(::Type{T}, mdb_val::MDB_val) where {T}
    unsafe_load(convert(Ptr{T}, mdb_val.mv_data))
end

"""Return the LMDB library version and version information

Function returns tuple `(VersionNumber, String)` that contains a library version and a library version string.
"""
function version()
    major = Cint[0]
    minor = Cint[0]
    patch = Cint[0]
    ver_str = _mdb_version(major, minor, patch)
    return VersionNumber(major[1],minor[1],patch[1]), unsafe_string(ver_str)
end

"""Return a string describing a given error code

Function returns description of the error as a string. It accepts following arguments:
* `err::Int32`: An error code.
"""
function errormsg(err::Cint)
    errstr = _mdb_strerror(err)
    return unsafe_string(errstr)
end

""" Check if binary flag is set in provided value"""
isflagset(value, flag) = (value & flag) == flag


text:
# a rough translation of openlibm/src/e_pow.c under following licence:

## Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
##
## Developed at SunSoft, a Sun Microsystems, Inc. business.
## Permission to use, copy, modify, and distribute this
## software is freely granted, provided that this notice
## is preserved.

# note: this doesn't handle edge cases or negative numbers correctly.

function log2_ext(x)
    sig,n = frexp(x)
    ax = 2sig; n -= 1

    # reduce to
    # k = 0, sqrt(3)/2 = 0.866 < ax < sqrt(3/2) = 1.22
    # k = 1, sqrt(3/2) = 1.22 < ax < sqrt(3) = 1.73

    bp = 1.0
    dp_h = 0.0
    dp_l = 0.0

    if ax < sqrt(3/2)
    elseif ax < sqrt(3)
        bp = 1.5
        # log2(1.5)
        dp_h = 5.84962487220764160156e-01
        dp_l = 1.35003920212974897128e-08
    else
        n += 1
        ax *= 0.5
    end
    # log(x) = 2s + 2/3 s^3 + 2/5 s^5...

    u = ax - bp # usual f, reduced to -0.28 < u < 0.23
    v = 1.0/(ax + bp)
    ss = u*v # ss = (x-1)/(x+1) or (x-1.5)/(x+1.5) < 0.125 = 2^-3
    s_h = trunc32(ss)
    t_h = trunc32(ax+bp) # according to comment
    t_l = ax - (t_h-bp)
    s_l = v*((u-s_h*t_h)-s_h*t_l) # division correction?

    z = ss*ss # < 2^-6
    # log(ax) = 2s + 2/3*(s^3 + r) = 2/3s (3 + s^2 + r)
    r = z*z*@horner(z,
                    5.99999999999994648725e-01,
                    4.28571428578550184252e-01,
                    3.33333329818377432918e-01,
                    2.72728123808534006489e-01,
                    2.30660745775561754067e-01,
                    2.06975017800338417784e-01) # < 2^-12
    r += s_l*(s_h+ss)
    s2 = s_h*s_h
    # t = @split_add(3.0,s2,r)
    t_h = trunc32(3.0 + s2 + r)
    t_l = r+(s2+(3.0-t_h))
    # u+v = ss*(1+...)
    # p = @split t*s
    u = s_h*t_h
    v = s_l*t_h+t_l*ss

    # 2/(3log2)
    cp    =  9.61796693925975554329e-01 # /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
    cp_h  =  9.61796700954437255859e-01 # /* 0x3FEEC709, 0xE0000000 =(float)cp */
    cp_l  = -7.02846165095275826516e-09 # /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/

    # 2/(3log2)*(ss+...)
    p_h = trunc32(u+v)
    p_l = v-(p_h-u)
    z_h = cp_h*p_h
    z_l = cp_l*p_h + p_l*cp + dp_l

    t = float64(n)
    t1 = trunc32(((z_h+z_l) + dp_h) + t)
    t2 = z_l - (((t1 - t) - dp_h) - z_h)
    Double64(t1,t2)
end

function mul_ext(t::Double64,y::Float64)
    # split y, y*(t1+t2)
    y1 = trunc32(y)
    p_l = (y-y1)*t.hi + y*t.lo
    p_h = y1*t.hi
    Double64(p_h,p_l)
end

function exp2(p::Double64)
    z = p.lo + p.hi
    # check for under/overflow

    # 2^(p_l+p_h)

    zr = round(z)
    n = int(zr)
    p_h = p.hi - zr
    p_l = p.lo

    lg2  =  6.93147180559945286227e-01 # /* 0x3FE62E42, 0xFEFA39EF */
    lg2_h  =  6.93147182464599609375e-01 # /* 0x3FE62E43, 0x00000000 */
    lg2_l  = -1.90465429995776804525e-09 # /* 0xBE205C61, 0x0CA86C39 */

    # 1) t1 = 2(e^z-1-z)/(e^z-1) = z - z^2/6 + z^4/360 - z^6/15120 + ...
    # 2) r = z*t1/(t1-2) = 1 - e^z + z = -(z^2/2 + s^3/6 + ...)
    # 3) e^z = 1 - r + z
    t = trunc32(p_l + p_h)
    u = t*lg2_h
    v = (p_l-(t-p_h))*lg2+t*lg2_l
    z = u+v
    w = v-(z-u);
    t  = z*z;
    t1  = z - t*@horner(t,
                        1.66666666666666019037e-01,
                        -2.77777777770155933842e-03,
                        6.61375632143793436117e-05,
                        -1.65339022054652515390e-06,
                        4.13813679705723846039e-08)
    r = (z*t1)/(t1 - 2.0) - (w+z*w)
    z = 1.0-(r-z)
    scalbn(z,n)
end

pow(x,y) = exp2(mul_ext(log2_ext(x),y))


text:
using DataStructures

# create output directory
if !isdir(out_dir * "/" * fig_dir)
    mkpath(out_dir * "/" * fig_dir)
end

# create list of figures to put in list
figures = (
    "",
    "_energy_error",
    "_toroidal_momentum_error"
)

# suffix of figure files
fig_suff = ".png"

# collect simulation parameters
function get_params(Δt, ntime, problem)
    params = OrderedDict()
    params["Δt"] = Δt
    params["nt"] = ntime
    params["problem"] = problem
    params["GeometricIntegrators.jl"] = package_version("GeometricIntegrators")
    params["ChargedParticleDynamics.jl"] = package_version("ChargedParticleDynamics")
    params["MagneticEquilibria.jl"] = package_version("MagneticEquilibria")
    params
end

function compute_error_drift(t, err, lint=100)
    const nint = div(t.n, lint)

    drift = zeros(eltype(err), nint)

    for i in 1:nint
        i1 = lint*(i-1)+1
        i2 = lint*i
        drift[i] = maximum(err[i1:i2])
    end

    drift - drift[1]
end

function compute_energy_error(t, q)
    h = zeros(eltype(q), q.nt+1)
    for i in 1:(q.nt+1)
        h[i] = problem.hamiltonian(t.t[i], q.d[:,i])
    end
    h_error = (h - h[1]) / h[1]
end

function compute_toroidal_momentum_error(t, q)
    P = zeros(eltype(q), q.nt+1)
    for i in 1:(q.nt+1)
        P[i] = problem.toroidal_momentum(t.t[i], q.d[:,i])
    end
    P_error = (P - P[1]) / P[1]
end

# function compute_toroidal_momentum_error(t, q)
#     Q = (q.d[:,1:end-1] .+ q.d[:,2:end]) ./ 2
#     P = zeros(eltype(q), q.nt)
#     for i in 1:(q.nt)
#         P[i] = toroidal_momentum(t.t[i], Q[:,i])
#     end
#     P_error = (P - P[1]) / P[1]
# end

function compute_momentum_error(t, q, p)
    p1_error = zeros(q.nt+1)
    p2_error = zeros(q.nt+1)
    p3_error = zeros(q.nt+1)
    p4_error = zeros(q.nt+1)

    for i in 1:(q.nt+1)
        p1_error[i] = p.d[1,i] - problem.α1(t.t[i], q.d[:,i])
        p2_error[i] = p.d[2,i] - problem.α2(t.t[i], q.d[:,i])
        p3_error[i] = p.d[3,i] - problem.α3(t.t[i], q.d[:,i])
        p4_error[i] = p.d[4,i] - problem.α4(t.t[i], q.d[:,i])
    end

    (p1_error, p2_error, p3_error, p4_error)
end

function compute_one_form(t, q)
    p1 = zeros(q.nt+1)
    p2 = zeros(q.nt+1)
    p3 = zeros(q.nt+1)
    p4 = zeros(q.nt+1)

    for i in 1:(q.nt+1)
        p1[i] = problem.α1(t.t[i], q.d[:,i])
        p2[i] = problem.α2(t.t[i], q.d[:,i])
        p3[i] = problem.α3(t.t[i], q.d[:,i])
        p4[i] = problem.α4(t.t[i], q.d[:,i])
    end

    (p1, p2, p3, p4)
end

function power10ticks(x, pos)
    if x == 0
        return "\$ 0 \$"
    end

    exponent    = @sprintf("%2d",   floor(Int64, log10(x)))
    coefficient = @sprintf("%2.0f", x / (10^floor(log10(x))))

    return "\$ $coefficient \\times 10^\{ $exponent \} \$"
end

xf = matplotlib[:ticker][:FuncFormatter](power10ticks)

yf = matplotlib[:ticker][:ScalarFormatter]()
yf[:set_powerlimits]((-1,+1))
yf[:set_scientific](true)
yf[:set_useOffset](true)

function run_guiding_center_4d(int::Integrator, ode::Equation, Δt, nt, filename)
    sol = Solution(ode, Δt, nt)

    try
        integrate!(int, sol)
    catch ex
        if isa(ex, DomainError)
            warn("DOMAIN ERROR")
        elseif isa(ex, ErrorException)
            warn("Simulation exited early.")
            warn(ex.msg)
        end
    end

    run_diagnostics(sol, Δt, nt, filename)
end

function run_guiding_center_4d(tableau::AbstractTableau, ode::Equation, Δt, nt; integrator=Integrator, filename=nothing)
    int = integrator(ode, tableau, Δt)

    println()
    println("Running ", tableau.name, " (", filename, ") for ", nt, " time steps...")

    run_guiding_center_4d(int, ode, Δt, nt, filename)
end

function run_guiding_center_4d(basis::Basis, quadrature::Quadrature, ode::Equation, Δt, nt; integrator=Integrator, filename=nothing)
    int = IntegratorCGVI(ode, basis, quadrature, Δt)

    println()
    println("Running ", filename, " for ", nt, " time steps of Δt = ", Δt, "...")

    run_guiding_center_4d(int, ode, Δt, nt, filename)
end

function run_diagnostics(sol, Δt, nt, filename)
    if filename ≠ nothing
        ndrift = div(nt, 10)

        if nt > 100000
            npl = nplot
        else
            npl = 1
        end

        if nt ≤ 200
            markersize=5
        else
            markersize=.5
        end

        # plot solution
        fig = figure(figsize=(5,5))
        subplots_adjust(left=0.25, right=0.95, top=0.96, bottom=0.16)
        plot(sol.q.d[1,1:npl:sol.counter+1], sol.q.d[2,1:npl:sol.counter+1], ".", markersize=markersize)
        xlim(XLIM)
        ylim(YLIM)
        xlabel(L"$q_{1} (t)$", labelpad=10, fontsize=20)
        ylabel(L"$q_{2} (t)$", labelpad=10, fontsize=20)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * fig_suff)
        close(fig)

        # plot energy error
        h_error = compute_energy_error(sol.t, sol.q)

        fig = figure(figsize=(8,5))
        subplots_adjust(left=0.16, right=0.95, top=0.95, bottom=0.16)
        plot(sol.t.t[1:npl:sol.counter+1], h_error[1:npl:sol.counter+1], ".", markersize=markersize)
        sol.counter > 0 ? xlim(sol.t[0], sol.t[sol.counter]) : xlim(sol.t[0], sol.t[end])
        xlabel(L"$t$", labelpad=10, fontsize=20)
        ylabel(L"$[H(t) - H(0)] / H(0)$", fontsize=20)
        ax = gca()
        ax[:xaxis][:set_major_formatter](xf)
        ax[:yaxis][:set_major_formatter](yf)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_energy_error" * fig_suff)
        close(fig)

        # plot energy drift
        h_drift = compute_error_drift(sol.t, h_error, ndrift)

        fig = figure(figsize=(8,5))
        subplots_adjust(left=0.16, right=0.95, top=0.95, bottom=0.16)
        plot(sol.t.t[div(ndrift,2)+1:ndrift:end], h_drift, ".", markersize=10)
        xlim(sol.t[0], sol.t[end])
        xlabel(L"$t$", labelpad=10, fontsize=20)
        ylabel(L"$H_{drift}(t)$", fontsize=20)
        ax = gca()
        ax[:xaxis][:set_major_formatter](xf)
        ax[:yaxis][:set_major_formatter](yf)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_energy_drift.png")
        close(fig)

        # plot toroidal momentum error
        P_error = compute_toroidal_momentum_error(sol.t, sol.q)

        fig = figure(figsize=(8,5))
        subplots_adjust(left=0.16, right=0.95, top=0.95, bottom=0.16)
        plot(sol.t.t[1:npl:sol.counter+1], P_error[1:npl:sol.counter+1], ".", markersize=markersize)
        sol.counter > 0 ? xlim(sol.t[0], sol.t[sol.counter]) : xlim(sol.t[0], sol.t[end])
        xlabel(L"$t$", labelpad=10, fontsize=20)
        ylabel(L"$[P(t) - P(0)] / P(0)$", fontsize=20)
        ax = gca()
        ax[:xaxis][:set_major_formatter](xf)
        ax[:yaxis][:set_major_formatter](yf)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_toroidal_momentum_error" * fig_suff)
        close(fig)

        # # plot one-form
        # theta = compute_one_form(sol.t, sol.q)
        #
        # for i in 1:sol.q.nd
        #     fig = figure(figsize=(6,5))
        #     subplots_adjust(left=0.18, right=0.95, top=0.95, bottom=0.16)
        #     plot(sol.t.t[1:npl:end], theta[i][1:npl:end], ".", markersize=markersize)
        #     xlim(sol.t[0], sol.t[end])
        #     if maximum(theta[i]) == 0.
        #         ylim(-2E-16, +2E-16)
        #     end
        #     xlabel("\$ t \$", labelpad=10, fontsize=20)
        #     ylabel("\$ \\vartheta_\{ {i} \} (q(t)) \$", labelpad=6, fontsize=20)
        #     ax = gca()
        #     ax[:xaxis][:set_major_formatter](xf)
        #     ax[:yaxis][:set_major_formatter](yf)
        #     savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_one_form" * string(i) * fig_suff)
        #     close(fig)
        # end

        # plot momentum error
        if isdefined(sol, :p)
            p_error = compute_momentum_error(sol.t, sol.q, sol.p)

            for i in 1:sol.p.nd
                fig = figure(figsize=(6,5))
                subplots_adjust(left=0.18, right=0.95, top=0.95, bottom=0.16)
                plot(sol.t.t[1:npl:sol.counter+1], p_error[i][1:npl:sol.counter+1], ".", markersize=markersize)
                sol.counter > 0 ? xlim(sol.t[0], sol.t[sol.counter]) : xlim(sol.t[0], sol.t[end])
                maximum(p_error[i]) == zero(eltype(p_error[i])) ? ylim(-2E-16, +2E-16) : nothing
                xlabel("\$ t \$", labelpad=10, fontsize=20)
                ylabel("\$ p_\{ {i} \} (t) - \\vartheta_\{ {i} \} (q(t)) \$", labelpad=6, fontsize=20)
                ax = gca()
                ax[:xaxis][:set_major_formatter](xf)
                ax[:yaxis][:set_major_formatter](yf)
                savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_momentum_error" * string(i) * fig_suff)
                close(fig)
            end
        end

        # plot Lagrange multiplier
        if isdefined(sol, :λ)
            for i in 1:sol.λ.nd
                fig = figure(figsize=(6,5))
                subplots_adjust(left=0.16, right=0.95, top=0.95, bottom=0.16)
                plot(sol.t.t[1:npl:sol.counter+1], sol.λ.d[i,1:npl:sol.counter+1], ".", markersize=markersize)
                sol.counter > 0 ? xlim(sol.t[0], sol.t[sol.counter]) : xlim(sol.t[0], sol.t[end])
                maximum(sol.λ.d[i,:]) == zero(eltype(sol.λ.d)) ? ylim(-2E-16, +2E-16) : nothing
                xlabel("\$t\$", labelpad=10, fontsize=20)
                ylabel("\$\\lambda_\{{i}\} (t)\$", labelpad=6, fontsize=20)
                ax = gca()
                ax[:xaxis][:set_major_formatter](xf)
                ax[:yaxis][:set_major_formatter](yf)
                savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_lambda" * string(i) * fig_suff)
                close(fig)
            end
        end
    end
end

function run_tableau_list(tableaus, ode, Δt, ntime)
    for tab in tableaus
        if length(tab) ≥ 4
            nt = tab[4]
        else
            nt = ntime
        end

        if length(tab) ≥ 5
            dt = tab[5]
        else
            dt = Δt
        end

        run_guiding_center_4d(tab[2], ode, dt, nt; integrator=tab[1], filename=tab[3])
    end
end

function run_tableau_list_cgvi(tableaus, ode, Δt, ntime)
    for tab in tableaus
        if length(tab) ≥ 4
            nt = tab[4]
        else
            nt = ntime
        end

        if length(tab) ≥ 5
            dt = tab[5]
        else
            dt = Δt
        end

        run_guiding_center_4d(tab[1], tab[2], ode, dt, nt; filename=tab[3])
    end
end


text:
# # Matrix operators

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
Many solutions of partial differential equations with immersed layers, particularly Poisson problems,
lead to saddle-point problems, in which the Schur complement operator is a matrix
composed from some of the surface-grid operators discussed in [Surface-grid operations](@ref).
The package provides some convenience tools for constructing these matrices. By
their nature, the construction of these matrices is slow, since each column
involves the application of the same set of operations. However, the point
of this construction is to do it once and store it for repeated application.
=#

#=
One common saddle point system is

$$A = \begin{bmatrix} L & R \\ R^T & 0 \end{bmatrix}$$

where $L$ is the discrete Laplacian and $R$ and $R^T$ are the regularization
and interpolation operators ([`regularize!`](@ref) and [`interpolate!`](@ref)),
respectively. This system arises in the solution of the Poisson equation
with Dirichlet boundary conditions on the immersed surface. The Schur complement of this is
$S = - R^T L^{-1} R$. This matrix can be obtained using the function [`create_RTLinvR`](@ref).
=#

#=
Another common saddle point system is

$$A = \begin{bmatrix} L & D_s \\ G_s & R_n^T R_n \end{bmatrix}$$

where $D_s$ and $G_s$ are the surface divergence and gradient operators
([`surface_divergence!`](@ref) and [`surface_grad!`](@ref)),
respectively, and $R_n$ and $R_n^T$ are [`regularize_normal!`](@ref) and [`normal_interpolate!`](@ref).
This system arises in the solution of the Poisson equation
with Neumann boundary conditions on the immersed surface. The Schur complement of this is
$S = R_n^T R_n - G_s L^{-1} D_s$. Each of the matrices in this are individually
provided by the package, by the functions [`create_nRTRn`](@ref) and [`create_GLinvD`](@ref),
respectively. However, it is useful to know that the sum of these two
matrices is exactly the matrix $-C_s L^{-1}C_s^T$, where $C_s$ and $C_s^T$ are
surface curl operators [`surface_curl!`](@ref). This complete matrix is provided by
[`create_CLinvCT`](@ref).
=#

#=
Another helpful matrix operator is the surface filter, given by

$$\tilde{R}^T R$$

where $\tilde{R}^T$ is a modified form of the interpolation operator,
designed to return the regularized field to the surface points while
maintaining the integral value of the original field [^1]. We can
obtain this matrix with [`create_surface_filter`](@ref).
=#

#md # ## Matrix construction functions

#md # ```@docs
#md # create_RTLinvR
#md # create_CLinvCT
#md # create_CL2invCT
#md # create_GLinvD
#md # create_GLinvD_cross
#md # create_nRTRn
#md # create_surface_filter
#md # ```

#md # [^1]: Goza, A., et al., (2016) "Accurate computation of surface stresses and forces with immersed boundary methods," J. Comput. Phys., 321, 860--873.


text:
using CMake, p7zip_jll, Tar
using Clang.Generators

function build_libbitfield_native()
    @info "Building libbitfield binary with native tools. This will break if CMake fails to find a C compiler"
    success = true
    try
        src_dir = joinpath(@__DIR__, "bitfield")
        build_dir = joinpath(@__DIR__, "build")

        config_cmd = `$cmake -B $build_dir -S $src_dir`
        if Sys.WORD_SIZE == 32
            if Sys.iswindows()
                config_cmd = `$config_cmd -A win32`
            elseif Sys.islinux()
                config_cmd = `$config_cmd -D CMAKE_C_FLAGS=-march=native -D CMAKE_CXX_FLAGS=-march=native`
            end
        end
        build_cmd = `$cmake --build $build_dir --config Debug`
        install_cmd = `$cmake --install $build_dir --config Debug --prefix $build_dir`

        run(config_cmd)
        run(build_cmd)
        run(install_cmd)
    catch e
        @warn "Building libbitfield with native tools failed: $e"
        success = false
    end
    return success
end

function build_libbitfield_binarybuilder()
    @info "Building libbitfield binary with BinaryBuilder."
    success = true
    try
        cd(@__DIR__) do
            run(`$(Base.julia_cmd()) --project bitfield/build_tarballs.jl`)
            # from Pkg.download_verify_unpack
            tarball_path = only(readdir("products"))
            dest = "build"
            rm(dest; recursive = true)
            Tar.extract(`$(p7zip_jll.p7zip()) x products/$tarball_path -so`, dest)
        end
    catch e
        @warn "Building libbitfield with BinaryBuilder failed: $e"
        success = false
    end
    return success
end

function build_libbitfield()
    success = true
    try
        # Compile binary
        if !build_libbitfield_binarybuilder() && !build_libbitfield_native()
            error("Could not build libbitfield binary")
        end

        # Generate wrappers
        @info "Building libbitfield wrapper"
        args = get_default_args()
        headers = joinpath(@__DIR__, "build", "include", "bitfield.h")
        options = load_options(joinpath(@__DIR__, "bitfield", "generate.toml"))
        lib_path = joinpath(@__DIR__, "build", "lib", Sys.iswindows() ? "bitfield.dll" : "libbitfield")
        options["general"]["library_name"] = "\"$(escape_string(lib_path))\""
        options["general"]["output_file_path"] = joinpath(@__DIR__, "LibBitField.jl")
        ctx = create_context(headers, args, options)
        build!(ctx)

        # Call a function to ensure build is successful
        include("LibBitField.jl")
        m = Base.@invokelatest LibBitField.Mirror(10, 1.5, 1e6, -4, 7, 3)
        Base.@invokelatest LibBitField.toBitfield(Ref(m))
    catch e
        @warn "Building libbitfield failed: $e"
        success = false
    end
    return success
end

@testset "Bitfield" begin
    if build_libbitfield()
        bf = Ref(LibBitField.BitField(Int8(10), 1.5, Int32(1e6), Int32(-4), Int32(7), UInt32(3)))
        m = Ref(LibBitField.Mirror(10, 1.5, 1e6, -4, 7, 3))
        GC.@preserve bf m begin
            pbf = Ptr{LibBitField.BitField}(pointer_from_objref(bf))
            pm = Ptr{LibBitField.Mirror}(pointer_from_objref(m))
            @test LibBitField.toMirror(bf) == m[]
            @test LibBitField.toBitfield(m).a == bf[].a
            @test LibBitField.toBitfield(m).b == bf[].b
            @test LibBitField.toBitfield(m).c == bf[].c
            @test LibBitField.toBitfield(m).d == bf[].d
            @test LibBitField.toBitfield(m).e == bf[].e
            @test LibBitField.toBitfield(m).f == bf[].f
        end
    end
end