text:
const tol = 10*eps()

"""
ItoIntegral
A struct detailing an ito integral. It contains a UnivariateFunction detailing the integrand as well as a symbol detailing an id of the integral's processes.

Usual (and most general) contructor is:
    ItoIntegral(brownian_id_::Symbol, f_::UnivariateFunction)
Convenience constructor for ItoIntegrals where the integrand is a flat function is:
    ItoIntegral(brownian_id_::Symbol, variance_::Real)
"""
struct ItoIntegral
    brownian_id_::Symbol
    f_::UnivariateFunction
    function ItoIntegral(brownian_id_::Symbol, variance_::Real)
        return new(brownian_id_, PE_Function(variance_, 0.0, 0.0, 0))
    end
    function ItoIntegral(brownian_id_::Symbol, f_::UnivariateFunction)
        return new(brownian_id_, f_)
    end
end

"""
    get_variance(ito::ItoIntegral, from::Real, to::Real)
    get_variance(ito::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime})
Get the variance of an ItoIntegral from one point of time to another.
"""
function get_variance(ito::ItoIntegral, from::Real, to::Real)
    if from + eps() > to
        return 0.0;
    end
    return evaluate_integral(ito.f_^2, from, to)
end
function get_variance(ito::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime})
    from_fl = years_between(from, base)
    to_fl   = years_between(to, base)
    return get_variance(ito, from_fl, to_fl)
end

"""
    get_volatility(ito::ItoIntegral, on::Union{Date,DateTime})
Get the volatility of an ItoIntegral on a certain date.
"""
function get_volatility(ito::ItoIntegral, on::Union{Date,DateTime})
    return evaluate(ito.f_, on)
end

"""
    get_covariance(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    get_covariance(ito1::ItoIntegral,ito2::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
Get the covariance of two ItoIntegrals over a certain period given the underlying Brownian processes have a correlation of gaussian_correlation.
"""
function get_covariance(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    if from + eps() >= to
        return 0.0;
    end
    return gaussian_correlation * evaluate_integral(ito1.f_ * ito2.f_, from, to)
end
function get_covariance(ito1::ItoIntegral,ito2::ItoIntegral, base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
    from_fl = years_between(from, base)
    to_fl   = years_between(to, base)
    return get_covariance(ito1, ito2, from_fl, to_fl, gaussian_correlation)
end

"""
    get_correlation(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    get_correlation(ito1::ItoIntegral,ito2::ItoIntegral,  base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
Get the correlation of two ItoIntegrals over a certain period given the underlying Brownian processes have a correlation of gaussian_correlation.
"""
function get_correlation(ito1::ItoIntegral,ito2::ItoIntegral, from::Real, to::Real, gaussian_correlation::Real)
    cov =  covar(ito1,ito2, base, from, to, gaussian_correlation)
    var1 = var(ito1, from, to)
    var2 = var(ito2, from, to)
    return gaussian_correlation * (cov / (var1 * var2))
end
function get_correlation(ito1::ItoIntegral,ito2::ItoIntegral,  base::Union{Date,DateTime}, from::Union{Date,DateTime}, to::Union{Date,DateTime}, gaussian_correlation::Real)
    from_fl = years_between(from, base)
    to_fl   = years_between(to, base)
    return get_correlation(ito1, ito2, from_fl, to_fl, gaussian_correlation)
end

"""
    brownians_in_use(itos::Array{ItoIntegral,1}, brownians::Array{Symbol,1})
Determine which Browninan processes are used in an array of ItoIntegrals.
"""
function brownians_in_use(itos::Dict{Symbol,ItoIntegral}, brownians::Array{Symbol,1})
    all_brownians_in_use = unique(map(x -> x.brownian_id_ , values(itos)))
    indices_in_use   = unique(findall(map(x -> x in all_brownians_in_use , brownians)))
    reduced_brownian_list = brownians[indices_in_use]
    return all_brownians_in_use, indices_in_use, reduced_brownian_list
end

"""
    ItoSet
Creates an ItoSet. This contains :
* A correlation matrix of brownian motions.
* A vector giving the axis labels for this correlation matrix.
* A dict of ItoInterals. Here the keys should be ids for the ito integrals and the values should be ItoIntegrals.
Determine which Brownian processes are used in an array of ItoIntegrals.
"""
struct ItoSet{T<:Real}
    brownian_correlation_matrix_::Symmetric{T}
    brownian_ids_::Array{Symbol,1}
    ito_integrals_::Dict{Symbol,ItoIntegral}
    function ItoSet(brownian_corr_matrix::Symmetric{T}, brownian_ids::Array{Symbol,1}, ito_integrals::Dict{Symbol,ItoIntegral}) where T<:Real
        if (size(brownian_ids)[1] != size(brownian_corr_matrix)[1])
            error("The shape of brownian_ids_ must match the number of rows/columns of brownian_correlation_matrix_")
        end
          all_brownians_in_use, used_brownian_indices, brown_ids = brownians_in_use(ito_integrals, brownian_ids)
          if length(setdiff(all_brownians_in_use, brownian_ids)) > 0
              error(string("In creating an ItoSet there are some brownian motions referenced by ito integrals for which there are no corresponding entries in the correlation matrix for brownian motions. Thus an ItoSet cannot be built. These include ", setdiff(all_brownians_in_use, brownian_ids)))
          end
          brownian_corr_matrix_subset      = Symmetric(brownian_corr_matrix[used_brownian_indices,used_brownian_indices])
       return new{T}(brownian_corr_matrix_subset, brown_ids, ito_integrals)
    end
end

"""
    get_correlation(ito::ItoSet, index1::Integer, index2::Integer)
    get_correlation(ito::ItoSet, brownian_id1::Symbol, brownian_id2::Symbol)
Get correlation between brownian motions in an ItoSet.
"""
function get_correlation(ito::ItoSet, index1::Integer, index2::Integer)
    return ito.brownian_correlation_matrix_[index1, index2]
end
function get_correlation(ito::ItoSet, brownian_id1::Symbol, brownian_id2::Symbol)
    index1 = findall(brownian_id1 .== ito.brownian_ids_)[1]
    index2 = findall(brownian_id2 .== ito.brownian_ids_)[1]
    return get_correlation(ito, index1, index2)
end

"""
    get_volatility(ito::ItoSet, index::Integer, on::Union{Date,DateTime})
    get_volatility(ito::ItoSet, ito_integral_id::Symbol, on::Union{Date,DateTime})
Get volatility of an ito_integral on a date.
"""
function get_volatility(ito::ItoSet, index::Integer, on::Union{Date,DateTime})
    return get_volatility(ito.ito_integrals_[index], on)
end

function get_volatility(ito::ItoSet, ito_integral_id::Symbol, on::Union{Date,DateTime})
    ito_integral = ito.ito_integrals_[ito_integral_id]
    return get_volatility(ito_integral, on)
end

"""
    make_covariance_matrix(ito_set_::ItoSet, from::Real, to::Real)
Make a covariance matrix given an ItoSet and a period of time.
"""
function make_covariance_matrix(ito_set_::ItoSet{T}, from::Real, to::Real) where T<:Real
    number_of_itos = length(ito_set_.ito_integrals_)
    ito_ids = collect(keys(ito_set_.ito_integrals_))
    cov = Array{T,2}(undef, number_of_itos,number_of_itos)
    for r in 1:number_of_itos
        rito = ito_set_.ito_integrals_[ito_ids[r]]
        for c in r:number_of_itos
            #if c < r
            #    cov[r,c] = 0.0 # Since at the end we use the Symmetric thing, this is discarded so we don't bother computing it.
            #end
            cito = ito_set_.ito_integrals_[ito_ids[c]]
            cr_correlation = get_correlation(ito_set_, rito.brownian_id_, cito.brownian_id_)
            cov[r,c] = get_covariance(rito, cito, from, to, cr_correlation)
        end
    end
    return Symmetric(cov), ito_ids
end

"""
    CovarianceAtDate
Creates an CovarianceAtDate object. This contains :
* An Itoset
* Time From
* Time To
And in the constructor the following items are generated and stored in the object:
* A covariance matrix
* Labels for the covariance matrix.
* The cholesky decomposition of the covariance matrix.
* The inverse of the covariance matrix.
* The determinant of the covariance matrix.
"""
struct CovarianceAtDate
    ito_set_::ItoSet
    from_::Real
    to_::Real
    covariance_::Symmetric
    covariance_labels_::Array{Symbol,1}
    chol_::Union{Missing,LowerTriangular}
    inverse_::Union{Missing,Symmetric}
    determinant_::Union{Missing,Real}
    """
    CovarianceAtDate(ito_set_::ItoSet, from_::Real, to_::Real;
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
    CovarianceAtDate(ito_set_::ItoSet, from::Union{Date,DateTime}, to::Union{Date,DateTime})
    CovarianceAtDate(old_CovarianceAtDate::CovarianceAtDate, from::Real, to::Real)
    CovarianceAtDate(old_CovarianceAtDate::CovarianceAtDate, from::Union{Date,DateTime}, to::Union{Date,DateTime})
        These are constructors for a CovarianceAtDate struct.
    """
    function CovarianceAtDate(ito_set_::ItoSet, from_::Real, to_::Real;
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
        covariance_, covariance_labels_ = make_covariance_matrix(ito_set_, from_, to_)
        chol_ = missing
        inverse_ = missing
        determinant_ = missing
        if calculate_chol chol_ = LowerTriangular(cholesky(covariance_).L) end
        if calculate_inverse inverse_           = Symmetric(inv(covariance_)) end
        if calculate_determinant determinant_       = det(covariance_) end
        return new(ito_set_, from_, to_, covariance_, covariance_labels_, chol_, inverse_, determinant_)
    end
    function CovarianceAtDate(ito_set_::ItoSet, from::Union{Date,DateTime}, to::Union{Date,DateTime};
             calculate_chol::Bool = true, calculate_inverse::Bool = true, calculate_determinant::Bool = true)
        from_ = years_from_global_base(from)
        to_   = years_from_global_base(to)
        return CovarianceAtDate(ito_set_, from_, to_; calculate_chol = calculate_chol,
                  calculate_inverse = calculate_inverse, calculate_determinant = calculate_determinant)
    end
    function CovarianceAtDate(old_CovarianceAtDate::CovarianceAtDate, from::Real, to::Real)
        return CovarianceAtDate(old_CovarianceAtDate.ito_set_, from, to)
    end
    function CovarianceAtDate(old_CovarianceAtDate::CovarianceAtDate, from::Date, to::Date)
        return CovarianceAtDate(old_CovarianceAtDate.ito_set_, from, to)
    end
end

"""
    get_volatility(covar::CovarianceAtDate, index::Integer, on::Union{Date,DateTime})
    get_volatility(covar::CovarianceAtDate, id::Symbol, on::Union{Date,DateTime})
Get the volatility of an ItoIntegral on a date..
"""
function get_volatility(covar::CovarianceAtDate, index::Integer, on::Union{Date,DateTime})
    return get_volatility(covar.ito_set_, index, on)
end
function get_volatility(covar::CovarianceAtDate, id::Symbol, on::Union{Date,DateTime})
    return get_volatility(covar.ito_set_, id, on)
end

"""
    get_variance(covar::CovarianceAtDate, id::Symbol)
    get_variance(covar::CovarianceAtDate, index::Integer)
Get the variance of an ItoIntegral over a period.
"""
function get_variance(covar::CovarianceAtDate, id::Symbol)
        index = findall(id .== covar.covariance_labels_)[1]
        return get_variance(covar, index)
end
function get_variance(covar::CovarianceAtDate, index::Integer)
    return covar.covariance_[index,index]
end

"""
    get_covariance(covar::CovarianceAtDate, index_1::Integer, index_2::Integer)
    get_covariance(covar::CovarianceAtDate, id1::Symbol, id2::Symbol)
Get the covariance of two ItoIntegrals over a period.
"""
function get_covariance(covar::CovarianceAtDate, index_1::Integer, index_2::Integer)
    return covar.covariance_[index_1,index_2]
end
function get_covariance(covar::CovarianceAtDate, id1::Symbol, id2::Symbol)
    index_1 = findall(id1 .== covar.covariance_labels_)[1]
    index_2 = findall(id2 .== covar.covariance_labels_)[1]
    return get_covariance(covar, index_1, index_2)
end

"""
    get_correlation(covar::CovarianceAtDate, index_1::Integer, index_2::Integer)
    get_correlation(covar::CovarianceAtDate, id1::Symbol, id2::Symbol)
Get the correlation of two ItoIntegrals over a period.
"""
function get_correlation(covar::CovarianceAtDate, index_1::Integer, index_2::Integer)
    covariance = get_covariance(covar, index_1, index_2)
    var1  = get_variance(covar, index_1)
    var2  = get_variance(covar, index_2)
    return covariance/sqrt(var1 * var2)
end
function get_correlation(covar::CovarianceAtDate, id1::Symbol, id2::Symbol)
    index_1 = findall(id1 .== covar.covariance_labels_)[1]
    index_2 = findall(id2 .== covar.covariance_labels_)[1]
    return get_correlation(covar, index_1, index_2)
end

## Random draws
"""
    get_draws(covar::CovarianceAtDate; uniform_draw::Array{T,1} = rand(length(covar.covariance_labels_))) where T<:Real
    get_draws(covar::CovarianceAtDate, num::Integer; twister::MersenneTwister = MersenneTwister(1234), antithetic_variates = false)
get pseudorandom draws from a CovarianceAtDate struct. Other schemes (like quasirandom) can be done by inserting quasirandom
numbers in as the uniform_draw.
If the antithetic_variates control is set to true then every second set of draws will be antithetic to the previous.
"""
function get_draws(covar::CovarianceAtDate; uniform_draw::Array{T,1} = rand(length(covar.covariance_labels_))) where T<:Real
    number_of_itos = length(covar.covariance_labels_)
    normal_draw = quantile.(Ref(Normal()), uniform_draw)
    scaled_draw = covar.chol_ * normal_draw
    first_set_of_draws = Dict{Symbol,Real}(covar.covariance_labels_ .=> scaled_draw)
    return first_set_of_draws
end
function get_draws(covar::CovarianceAtDate, num::Integer; number_generator::NumberGenerator = Mersenne(MersenneTwister(1234), length(covar.covariance_labels_)), antithetic_variates = false)
    if antithetic_variates
        half_num = convert(Int, round(num/2))
        array_of_dicts = Array{Dict{Symbol,Real}}(undef, half_num*2)
        number_of_itos = length(covar.covariance_labels_)
        for i in 1:half_num
            first_entry = 2*(i-1)+1
            draw = next!(number_generator)
            array_of_dicts[first_entry] = get_draws(covar; uniform_draw = draw)
            array_of_dicts[first_entry + 1] = get_draws(covar; uniform_draw = 1.0 .- draw)
        end
        return array_of_dicts
    else
        array_of_dicts = Array{Dict{Symbol,Real}}(undef, num)
        number_of_itos = length(covar.covariance_labels_)
        for i in 1:num
            draw =  next!(number_generator)
            array_of_dicts[i] = get_draws(covar; uniform_draw = draw)
        end
        return array_of_dicts
    end
end

# This is most likely useful for bug hunting.
"""
    get_zero_draws(covar::CovarianceAtDate)
get a draw of zero for all ito_integrals. May be handy for bug hunting.
"""
function get_zero_draws(covar::CovarianceAtDate)
    return Dict{Symbol,Real}(covar.covariance_labels_ .=> 0.0)
end
# This is most likely useful for bug hunting.
"""
    get_zero_draws(covar::CovarianceAtDate, num::Integer)
get an array of zero draws for all ito_integrals. May be handy for bug hunting.
"""
function get_zero_draws(covar::CovarianceAtDate, num::Integer)
    array_of_dicts = Array{Dict{Symbol,Real}}(undef, num)
    for i in 1:num
        array_of_dicts[i] = get_zero_draws(covar)
    end
    return array_of_dicts
end

"""
    pdf(covar::CovarianceAtDate, coordinates::Dict{Symbol,Real})
get the value of the pdf at some coordinates.
"""
function pdf(covar::CovarianceAtDate, coordinates::Dict{Symbol,Real})
    # The pdf is det(2\pi\Sigma)^{-0.5}\exp(-0.5(x - \mu)^\prime \Sigma^{-1} (x - \mu))
    # Where Sigma is covariance matrix, \mu is means (0 in this case) and x is the coordinates.
    rank_of_matrix = length(covar.covariance_labels_)
    x = get.(Ref(coordinates), covar.covariance_labels_, 0)
    one_on_sqrt_of_det_two_pi_covar     = 1/(sqrt(covar.determinant_) * (2*pi)^(rank_of_matrix/2))
    return one_on_sqrt_of_det_two_pi_covar * exp(-0.5 * x' * covar.inverse_ * x)
end

"""
    log_likelihood(covar::CovarianceAtDate, coordinates::Dict{Symbol,Real})
get the log likelihood at some coordinates.
"""
function log_likelihood(covar::CovarianceAtDate, coordinates::Dict{Symbol,Real})
    # The pdf is det(2\pi\Sigma)^{-0.5}\exp(-0.5(x - \mu)^\prime \Sigma(x - \mu))
    # Where Sigma is covariance matrix, \mu is means (0 in this case) and x is the coordinates.
    rank_of_matrix = length(covar.covariance_labels_)
    x = get.(Ref(coordinates), covar.covariance_labels_, 0)
    one_on_sqrt_of_det_two_pi_covar     = -0.5*log(covar.determinant_)  +  (rank_of_matrix/2)*log(2*pi)
    return one_on_sqrt_of_det_two_pi_covar + (-0.5 * x' * covar.inverse_ * x)
end


text:
function lock(name; wait_for_semaphore="", path=metadatadir())
    lock_path =  joinpath(path,"$name.lck")
    for _ in 1:max_lock_retries
        if wait_for_semaphore == "" || semaphore_status(wait_for_semaphore) == 0
            try
                mkdir(lock_path)
                return
            catch e
                sleep(0.1)
            end
        end
    end
    error("Could not retriev lock $name")
end

function unlock(name; path=metadatadir())
    lock_path =  joinpath(path,"$name.lck")
    try
        mkdir(lock_path)
    catch e
        rm(lock_path)
        return
    end
    rm(lock_path)
    error("$name is currently unlocked.")
end

function semaphore_status(name; path=metadatadir())
    sem_path =  joinpath(path,"$name.sem")
    lock(name, path=path)
    if isfile(sem_path)
        n = parse(Int,read(sem_path,String))
    else
        n = 0
    end
    unlock(name, path=path)
    return n
end

function semaphore_enter(name; path=metadatadir())
    sem_path =  joinpath(path,"$name.sem")
    lock(name, path=path)
    if isfile(sem_path)
        n = parse(Int,read(sem_path,String))
    else
        n = 0
    end
    open(sem_path,"w") do f
        write(f, string(n+1))
    end
    unlock(name, path=path)
end

function semaphore_exit(name; path=metadatadir())
    sem_path =  joinpath(path,"$name.sem")
    lock(name, path=path)
    if isfile(sem_path)
        n = parse(Int,read(sem_path,String))
        if n == 1
            rm(sem_path)
        else
            open(sem_path,"w") do f
                write(f, string(n-1))
            end
        end
    else
        unlock(name, path=path)
        error("Semaphore $name is out of balance. Expected a file but there is none")
    end
    unlock(name, path=path)
end

ext: jl

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
export read_file, write_file
export add_int_var, add_float_var, remove_var, add_empty_record, set_number_of_records, remove_record
export find_variable, get_number_of_records, get_number_of_variables
export get_int, get_float, set_int, set_float, set_missing, get_missing_int, get_missing_float
export is_missing, is_discrete, get_int_data, get_float_data, get_id, set_id
export get_num_states, get_state_names, set_state_names
export has_missing_data, is_constant
export discretize, discretize_with_info
export disc_alg_string
export DSL_DISCRETIZE_HIERARCHICAL, DSL_DISCRETIZE_UNIFORMWIDTH, DSL_DISCRETIZE_UNIFORMCOUNT
export DSL_MISSING_INT, DSL_MISSING_FLOAT

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

const DSL_DISCRETIZE_HIERARCHICAL = unsafe_load(cglobal((:EDSL_DISCRETIZE_HIERARCHICAL, LIB_SMILE), Cint))
const DSL_DISCRETIZE_UNIFORMWIDTH = unsafe_load(cglobal((:EDSL_DISCRETIZE_UNIFORMWIDTH, LIB_SMILE), Cint))
const DSL_DISCRETIZE_UNIFORMCOUNT = unsafe_load(cglobal((:EDSL_DISCRETIZE_UNIFORMCOUNT, LIB_SMILE), Cint))

const DSL_MISSING_INT             = unsafe_load(cglobal((:EDSL_MISSING_INT,             LIB_SMILE), Cint))
const DSL_MISSING_FLOAT           = unsafe_load(cglobal((:EDSL_MISSING_FLOAT,           LIB_SMILE), Cfloat))

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

function read_file( dset::Dataset, filename::AbstractString )
	retval = ccall( (:dataset_readFile, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8}), dset.ptr, bytestring(filename) )
	return retval == 0 # true if success
end
function read_file( dset::Dataset, filename::AbstractString, pp::DatasetParseParams )
	retval = ccall( (:dataset_readFile_withparams, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8},Ptr{Void}), dset.ptr, bytestring(filename), pp.ptr )
	return retval == 0 # true if success
end
function write_file( dset::Dataset, filename::AbstractString )
	retval = ccall( (:dataset_writeFile, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8}), dset.ptr, bytestring(filename) )
	return retval == 0 # true if success
end
function write_file( dset::Dataset, filename::AbstractString, wp::DatasetWriteParams )
	retval = ccall( (:dataset_writeFile_withparams, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8},Ptr{Void}), dset.ptr, bytestring(filename), wp.ptr )
	return retval == 0 # true if success
end

function add_int_var( dset::Dataset, id::AbstractString, missingValue::Integer = DSL_MISSING_INT )

	ccall( (:dataset_addIntVar, LIB_SMILE), Cint, (Ptr{Void},Ptr{UInt8},Cint), dset.ptr, bytestring(id), missingValue )
end
function add_float_var( dset::Dataset, varName::AbstractString, missingValue::Cfloat = DSL_MISSING_FLOAT )

	ccall( (:dataset_addFloatVar, LIB_SMILE), Cint, (Ptr{Void},Ptr{UInt8},Cfloat), dset.ptr, bytestring(varName), missingValue )
end
function remove_var( dset::Dataset, var::Integer )

	ccall( (:dataset_removeVar, LIB_SMILE), Void, (Ptr{Void},Int32), dset.ptr, var )
end
function add_empty_record( dset::Dataset )

	ccall( (:dataset_addEmptyRecord, LIB_SMILE), Void, (Ptr{Void},), dset.ptr )
end
function set_number_of_records( dset::Dataset, numRecords::Integer )

	ccall( (:dataset_setNumberOfRecords, LIB_SMILE), Void, (Ptr{Void},Int32), dset.ptr, numRecords )
end
function remove_record( dset::Dataset, rec::Integer )

	ccall( (:dataset_removeRecord, LIB_SMILE), Cint, (Ptr{Void},Cint), dset.ptr, rec )
end

function find_variable( dset::Dataset, id::AbstractString )

	ccall( (:dataset_findVariable, LIB_SMILE), Cint, (Ptr{Void},Ptr{UInt8}), dset.ptr, bytestring(id) )
end
function get_number_of_variables( dset::Dataset )

	ccall( (:dataset_getNumberOfVariables, LIB_SMILE), Cint, (Ptr{Void},), dset.ptr )
end
function get_number_of_records( dset::Dataset )

	ccall( (:dataset_getNumberOfRecords, LIB_SMILE), Cint, (Ptr{Void},), dset.ptr )
end

function get_int( dset::Dataset, var::Integer, rec::Integer )
	@assert( is_discrete(dset, var) )
	ccall( (:dataset_getInt, LIB_SMILE), Cint, (Ptr{Void},Cint,Cint), dset.ptr, var, rec )
end
function get_float( dset::Dataset, var::Integer, rec::Integer )
	@assert( !is_discrete(dset, var) )
	ccall( (:dataset_getFloat, LIB_SMILE), Cfloat, (Ptr{Void},Cint,Cint), dset.ptr, var, rec )
end
function set_int( dset::Dataset, var::Integer, rec::Integer, value::Integer )
	@assert( is_discrete(dset, var) )
	ccall( (:dataset_setInt, LIB_SMILE), Void, (Ptr{Void},Cint,Cint,Cint), dset.ptr, var, rec, value )
end
function set_float( dset::Dataset, var::Integer, rec::Integer, value::AbstractFloat )

	ccall( (:dataset_setFloat, LIB_SMILE), Void, (Ptr{Void},Cint,Cint,Cfloat), dset.ptr, var, rec, value )
end
set_missing( dset::Dataset, var::Integer, rec::Integer ) =
	ccall( (:dataset_setMissing, LIB_SMILE), Void, (Ptr{Void},Cint,Cint), dset.ptr, var, rec )
get_missing_int( dset::Dataset, var::Integer ) =
	ccall( (:dataset_getMissingInt, LIB_SMILE), Cint, (Ptr{Void},Cint), dset.ptr, var )
function get_missing_float( dset::Dataset, var::Integer )

	ccall( (:dataset_getMissingFloat, LIB_SMILE), Cfloat, (Ptr{Void},Cint), dset.ptr, var )
end
function is_missing( dset::Dataset, var::Integer, rec::Integer )

	ccall( (:dataset_isMissing, LIB_SMILE), Bool, (Ptr{Void},Cint,Cint), dset.ptr, var, rec )
end
function is_discrete( dset::Dataset, var::Integer )

	ccall( (:dataset_isDiscrete, LIB_SMILE), Bool, (Ptr{Void},Cint), dset.ptr, var )
end
function get_int_data( dset::Dataset, var::Integer )
	n = get_number_of_records(dset)
	retval = Array(Cint, n)
	ccall( (:dataset_getIntData, LIB_SMILE), Void, (Ptr{Void},Cint,Ptr{Cint}), dset.ptr, var, retval )
	retval
end
function get_float_data( dset::Dataset, var::Integer )
	n = get_number_of_records(dset)
	retval = Array(Cfloat, n)
	ccall( (:dataset_getFloatData, LIB_SMILE), Void, (Ptr{Void},Cint,Ptr{Cfloat}), dset.ptr, var, retval )
	retval
end
function get_id( dset::Dataset, var::Integer )
	str_buf = Array(UInt8, 128)
	ccall( (:dataset_getId, LIB_SMILE), Void, (Ptr{Void},Cint,Ptr{UInt8}), dset.ptr, var, str_buf )
	bytestring(pointer(str_buf))
end
function set_id( dset::Dataset, var::Integer, newId::AbstractString )

	ccall( (:dataset_setId, LIB_SMILE), Cint, (Ptr{Void},Cint,Ptr{UInt8}), dset.ptr, var, bytestring(newId) )
end

function get_num_states( dset::Dataset, var::Integer )

	ccall( (:dataset_getNumStates, LIB_SMILE), Cint, (Ptr{Void},Cint), dset.ptr, var )
end
function get_state_names( dset::Dataset, var::Integer )
	n = get_num_states(dset, var)
	retval = Array(ASCIIString, n)
	ccall( (:dataset_getStateNames, LIB_SMILE), Void, (Ptr{Void},Cint,Ptr{Ptr{UInt8}}), dset.ptr, var, retval )
	retval
end
function set_state_names( dset::Dataset, var::Integer, newNames::Vector{ASCIIString} )
	@assert( length(newNames) == get_num_states(dset, var) )
	ccall( (:dataset_setStateNames, LIB_SMILE), Cint, (Ptr{Void},Cint,Ptr{Ptr{UInt8}}), dset.ptr, var, newNames )
end

function getindex( dset::Dataset, var::Integer, rec::Integer )
	if is_discrete(dset, var)
		return get_int(dset, var, rec)
	else
		return get_float(dset, var, rec)
	end
end
function getindex{I<:Integer}( dset::Dataset, var::Integer, rng::UnitRange{I} )
	if is_discrete(dset, var)
		return map(i->get_int(dset, var, i), rng)
	else
		return map(i->get_float(dset, var, i), rng)
	end
end
function getindex( dset::Dataset, var::Integer )
	m = get_number_of_records(dset)
	if is_discrete(dset, var)
		return map(i->get_int(dset, var, i), 0:(m-1))
	else
		return map(i->get_float(dset, var, i), 0:(m-1))
	end
end
function setindex!( dset::Dataset, val, var::Integer, rec::Integer )
	if is_discrete( dset, var )
		set_int(dset, var, rec, val)
	else
		set_float(dset, var, rec, val)
	end
end
function setindex!{I<:Integer}( dset::Dataset, val::Vector, var::Integer, rng::UnitRange{I} )
	@assert(length(val) == length(rng))
	if is_discrete( dset, var )
		for (i,ind) in enumerate(rng)
			set_int(dset, var, ind, val[i])
		end

	else
		for (i,ind) in enumerate(rng)
			set_float(dset, var, ind, val[i])
		end
	end
end

function has_missing_data( dset::Dataset, var::Integer )

	ccall( (:dataset_hasMissingData, LIB_SMILE), Bool, (Ptr{Void},Int32), dset.ptr, var )
end
function is_constant( dset::Dataset, var::Integer )

	ccall( (:dataset_isConstant, LIB_SMILE), Bool, (Ptr{Void},Int32), dset.ptr, var )
end

function discretize{T<:Real}(
	dset        :: Dataset,
	var         :: Integer,
	nBins       :: Integer;
 	statePrefix :: AbstractString    = "s",
 	edges       :: Vector{T} = Float64[], # seems to require nBins+1 edges but only first nBins-1 edges matter, and edge[0] should be the larger bound of the first bin
 	algorithm   :: Cint      = DSL_DISCRETIZE_HIERARCHICAL
 	)

	@assert( nBins > 1 )
	@assert( algorithm == DSL_DISCRETIZE_HIERARCHICAL || algorithm == DSL_DISCRETIZE_UNIFORMWIDTH || algorithm == DSL_DISCRETIZE_UNIFORMCOUNT )

	if isempty(edges)
		ccall( (:dataset_discretize, LIB_SMILE), Void, (Ptr{Void},Cint,Cint,Ptr{UInt8},Cint), dset.ptr, var, nBins, bytestring(statePrefix), algorithm )
	else
		@assert(length(edges) == nBins+1)
		ccall( (:dataset_discretize_withedges, LIB_SMILE), Void,
			(Ptr{Void},Cint,Cint,Ptr{UInt8},Ptr{Float64}),
			dset.ptr, var, nBins, bytestring(statePrefix), edges )
	end
end
function discretize_with_info(
	dset        :: Dataset,
	var         :: Integer,
	nBins       :: Integer;
	statePrefix :: AbstractString = "s",
	algorithm   :: Cint   = DSL_DISCRETIZE_HIERARCHICAL
	)

	# the edge list contains the boundaries between edges
	# |edges| = nBins-1

	@assert( nBins > 1 )
	@assert( algorithm == DSL_DISCRETIZE_HIERARCHICAL || algorithm == DSL_DISCRETIZE_UNIFORMWIDTH || algorithm == DSL_DISCRETIZE_UNIFORMCOUNT )

	binEdges = Array(Float64, nBins-1)
	nActualBins = ccall( (:dataset_discretize_getedges, LIB_SMILE), UInt32,
		(Ptr{Void},Cint,Cint,Ptr{UInt8},Cint,Ptr{Float64}),
		dset.ptr, var, nBins, bytestring(statePrefix), algorithm, binEdges )
	binEdges[1:nActualBins]
end
function discretize{T<:Real}(
	data  :: Vector{T},
	nBins :: Integer;
	statePrefix :: AbstractString = "s",
	algorithm   :: Cint   = DSL_DISCRETIZE_HIERARCHICAL
	)

	n = length(data)
	@assert(n > 1)

	dset = Dataset()
	add_float_var(dset, "A")
	set_number_of_records(dset, n)
	dset[0,0:n-1] = float32(data)

	minv, maxv = extrema(data)
	binedges = discretize_with_info(dset, 0, nBins, statePrefix=statePrefix, algorithm=algorithm)
	discretized = dset[0,0:n-1]
	(discretized, [minv, binedges, maxv])
end

function disc_alg_string( alg::Cint )
	if alg == DSL_DISCRETIZE_HIERARCHICAL
		return "HIERARCHICAL"
	elseif alg == DSL_DISCRETIZE_UNIFORMWIDTH
		return "UNIFORMWIDTH"
	elseif alg == DSL_DISCRETIZE_UNIFORMCOUNT
		return "UNIFORMCOUNT"
	end
	error("alg not found")
end

# DatasetParseParams - a struct within DSL_dataset
##################################################
function getindex( pp::DatasetParseParams, entry::Symbol )
	if entry == :missingValueToken
		str_buf = Array(UInt8, 128)
		ccall( (:datasetparseparams_getMissingValueToken, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, str_buf )
		return bytestring(pointer(str_buf))
	elseif entry == :missingInt
		return ccall( (:datasetparseparams_getMissingInt, LIB_SMILE), Int32, (Ptr{Void},), pp.ptr )
	elseif entry == :missingFloat
		return ccall( (:datasetparseparams_getMissingFloat, LIB_SMILE), Float32, (Ptr{Void},), pp.ptr )
	elseif entry == :columnIdsPresent
		return ccall( (:datasetparseparams_getColumnIdsPresent, LIB_SMILE), Bool, (Ptr{Void},), pp.ptr )
	end
	error("Index not found")
end
function setindex!( pp::DatasetParseParams, value, entry::Symbol)
	if entry == :missingValueToken
		@assert(isa(value, AbstractString))
		ccall( (:datasetparseparams_setMissingValueToken, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, bytestring(value) )
	elseif entry == :missingInt
		return ccall( (:datasetparseparams_setMissingInt, LIB_SMILE), Void, (Ptr{Void},Int32), pp.ptr, value )
	elseif entry == :missingFloat
		return ccall( (:datasetparseparams_setMissingFloat, LIB_SMILE), Void, (Ptr{Void},Float32), pp.ptr, value )
	elseif entry == :columnIdsPresent
		return ccall( (:datasetparseparams_setColumnIdsPresent, LIB_SMILE), Void, (Ptr{Void},Bool), pp.ptr, value )
	else
		error("Index not found")
	end
end

# DatasetWriteParams - a struct within DSL_dataset
##################################################
function getindex( pp::DatasetWriteParams, entry::Symbol )
	if entry == :missingValueToken
		str_buf = Array(UInt8, 128)
		ccall( (:datasetwriteparams_getMissingValueToken, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, str_buf )
		return bytestring(pointer(str_buf))
	elseif entry == :columnIdsPresent
		return ccall( (:datasetwriteparams_getColumnIdsPresent, LIB_SMILE), Bool, (Ptr{Void},), pp.ptr )
	elseif entry == :useStateIndices
		return ccall( (:datasetwriteparams_getUseStateIndices, LIB_SMILE), Bool, (Ptr{Void},), pp.ptr )
	elseif entry == :separator
		return Char(ccall( (:datasetwriteparams_getSeparator, LIB_SMILE), UInt8, (Ptr{Void},), pp.ptr ))
	elseif entry == :floatFormat
		str_buf = Array(UInt8, 128)
		ccall( (:datasetwriteparams_getFloatFormat, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, str_buf )
		return bytestring(pointer(str_buf))
	end
	error("Index not found")
end
function setindex!( pp::DatasetWriteParams, value, entry::Symbol)
	if entry == :missingValueToken
		@assert(isa(value, AbstractString))
		ccall( (:datasetwriteparams_setMissingValueToken, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, bytestring(value) )
	elseif entry == :columnIdsPresent
		ccall( (:datasetwriteparams_setColumnIdsPresent, LIB_SMILE), Void, (Ptr{Void},Bool), pp.ptr, value )
	elseif entry == :useStateIndices
		ccall( (:datasetwriteparams_setUseStateIndices, LIB_SMILE), Void, (Ptr{Void},Bool), pp.ptr, value )
	elseif entry == :separator
		ccall( (:datasetwriteparams_setSeparator, LIB_SMILE), Void, (Ptr{Void},UInt8), pp.ptr, value )
	elseif entry == :floatFormat
		@assert(isa(value, AbstractString))
		ccall( (:datasetwriteparams_setFloatFormat, LIB_SMILE), Void, (Ptr{Void},Ptr{UInt8}), pp.ptr, bytestring(value) )
	else
		error("Index not found")
	end
end


text:
module FTFSpectrum

using Optim
using BlackBoxOptim
using Statistics
using ..CommonStructures: ModelFit, BoundariesFit, PlaneWave, LayerTMMO, FitProcedure
using ..CommonStructures: _get_beam_parameters
using ..Utils: build_interpolation, average_polarisation, flatten_arrays, array_arrays, _oscillators_input
using ..TransferMatrixMethod: transfer_matrix
using ..RI

export space_solution_ema,
       theoretical_spectrum,
       normalize_reflectance,
       fit_tmm_optics,
       Reflectance,
       Transmittance,
       Ellipsometry,
       NoAlpha,
       UseAlpha,
       MeanAbs,
       SumAbs,
       SumMeanAbs

## Definition of subtypes
mutable struct Ellipsometry{T1,T2} <: FitProcedure where {T1<:Real, T2<:Real}
    X::Array{T1,2}
    Xexp::Array{T2,2}
    function Ellipsometry(
        X::Array{T1,2};
        Xexp::Array{T2,2}=zeros(size(X)),
    ) where {T1<:Real, T2<:Real}
      return new{T1,T2}(X, Xexp)
    end
end
Ellipsometry() = Ellipsometry([]; Xexp=zeros(1,2))

struct Reflectance{T1} <: FitProcedure where {T1<:Real}
    Xexp::Array{T1,1}
    function Reflectance(Xexp::Array{T1,1}) where {T1<:Real}
      return new{T1}(Xexp)
    end
end
Reflectance() = Reflectance([0.0])

struct Transmittance{T1} <: FitProcedure where {T1<:Real}
    Xexp::Array{T1,1}
    function Transmittance(Xexp::Array{T1,1}) where {T1<:Real}
      return new{T1}(Xexp)
    end
end
Transmittance() = Transmittance([0.0])

struct MeanAbsObjFun <: FitProcedure end
MeanAbs() = MeanAbsObjFun()
struct SumAbsObjFun <: FitProcedure end
SumAbs() = SumAbsObjFun()
struct SumMeanAbsObjFun <: FitProcedure end
SumMeanAbs() = SumMeanAbsObjFun()

# Non-exported
struct DoNotAlpha <: FitProcedure end
NoAlpha() = DoNotAlpha()
struct DoAlpha <: FitProcedure end
UseAlpha() = DoAlpha()

# Results of the solution space parameters
struct SpaceSolutionEMA{T1,T2} <: FitProcedure where {T1<:Float64, T2<:PlaneWave}
    solSpace::Array{T1}
    objfunMin::T1
    optOD::T1
    optParams::T1
    optThickness::T1
    spectrumFit::Array
    spectrumExp::Array
    od
    p
    beam::T2
end

# Solution of the optimization with the Optim solver for multilayers
struct FitTMMOptics{T1,T2,T3} <: FitProcedure where {T1<:Float64, T2<:PlaneWave, T3<:NamedTuple}
    spectrumFit::Array{T1}
    spectrumExp::Array{T1}
    optParams::Array{Array{T1,1},1}
    objfunMin::T1
    beam::T2
    layers::Array
    fitOptions::T3
end

"""

    Normalize the experimental reflectance dividing by the reference reflectance and
    multiplying by the theoretical one. The experimental and the reference reflectances
    must have the same scale.

        R = normalize_reflectance(λ, Rexp, Rthe, Rref)

            λ: wavelength of interest [nm]
            Rexp: 2-columns array with wavelength in first one and experimental
                reflectance in the second
            Rthe: 2-columns array with wavelength in first one and theoretical
                reflectance in the second
            Rref: 2-columns array with wavelength in first one and reference reflectance
                in the second

    The output has the same scale as the theoretical reflectance with the same length as λ.

"""
function normalize_reflectance(
    λ::AbstractArray{T0,1}, Rexp::Array{T1,2}, Rthe::Array{T1,2}, Rref::Array{T1,2},
) where {T0<:Real, T1<:Float64}
    # Interpolate for the λ range
    spl_exp = build_interpolation(Rexp)
    spl_ref = build_interpolation(Rref)
    spl_the = build_interpolation(Rthe)
    exp2ref = spl_exp.(λ) ./ spl_ref.(λ)
    # Find out if the sample is overshooting the reference
    (maximum(exp2ref) < 1.0) || return (exp2ref./maximum(exp2ref).*spl_the.(λ))
    return exp2ref.*spl_the.(λ)
end

"""

    Returns the polarisation averaged theoretical spectrum calculated with the transfer
    matrix method. It consiers light hitting the first medium (incidentRI).

        X = theoretical_spectrum(specType, beam, incidentRI, emergentRI)

            specType: type of spectrum to compute, Reflectance() or Transmittance().
            beam: a PlaneWave structure.
            incidentRI: Array containing the index of refraction of the incident material
                for several wavelengths. For instance, RIdb.air(beam.λ).
            emergentRI: Array containing the index of refraction of the emergent material
                for several wavelengths. For instance, RIdb.silicon(beam.λ).
            X: averaged spectrum = beam.p*Xp + (1.0 - beam.p)*Xs

"""
function theoretical_spectrum(
    specType::T0, beam::T1, N1::Array{T2,1}, N2::Array{T2,1},
) where {T0<:FitProcedure, T1<:PlaneWave, T2<:ComplexF64}
    beam_, _ = _get_beam_parameters(beam)
    X = _compute_transfer_matrix(specType, vec([0. 100.0 0.]), [N1 N2 N2], beam_)
    return X
end

"""

    Computes the solution space for a window range of two parameters, the optical thickness
    and porosity. Only works with these two parameters and for a single layer between
    incident and emergent media.

        sol = space_solution_ema(specType, b, beam, Xexp, layers; oftype=MeanAbs())

            specType: Type of spectrum to fit, with accepted values:
                Reflectance(Xexp), Transmittance(Xexp) or Ellipsometry(Xexp).
                For the Reflectance() and Transmittance() you must pass the spectrum array
                Xexp (between 0 and 1) to be fitted. For the Ellipsometry() you have to
                pass the Ψ (first column) and Δ (second column) spectra in degrees.
            b: lower and upper boundaries for the optical thickness and porosity.
               (e.g. [odl, odu, pl pu; Nod, Np]), where Nod and Np indicates the number of
               points to search in the solution space as optional parameters.
            beam: structure from PlaneWave
            layers: columnwise array of three LayerTMMO and ModelFit with information
                about the layers. For space_solution_ema: length(vec(layers)) = 3. Notice that
                the type of the first and third layers is LayerTMMO while the type of the
                second layers is ModelFit.
                oftype: Metric for the objective function to optimise.
                    Can take values: MeanAbs() (default), SumAbs() and SumMeanAbs().
            sol: structure with the solution
                solSpace: matrix with objective function values for the whole space
                objfunMin: minimum value of the objective function found
                optOD: optimum value of optical thickness
                optParams: optimum value of the other parameter
                optThickness: optimum value of thickness
                spectrumFit: optimum spectrum obtained with dopt and xopt
                spectrumExp: experimental spectrum
                od: range of optical thicknesses explored
                p: range of porosities explored
                beam: information of the light used

"""
function space_solution_ema(
    specType::T0, b::T1, beam::T2, layers::Array;
    oftype::T4=MeanAbs(),
) where{T0<:FitProcedure, T1<:BoundariesFit, T2<:PlaneWave, T3<:Float64, T4<:FitProcedure}
    # Generate grids
    _d = LinRange(b.odlo, b.odup, b.Nod)
    _p = LinRange(b.plo, b.pup, b.Np)
    d = similar(_d)
    sol = zeros(Float64, b.Nod, b.Np)
    beam_, _ = _get_beam_parameters(beam)
    # If it is ellipsometry, then fit the norm of ρ
    if isa(specType, Ellipsometry)
        ρexp = @. tan(deg2rad(specType.X[:,1]))*exp(im*deg2rad(specType.X[:,2]))
        specType.Xexp = hcat(real.(ρexp), imag.(ρexp))
        # specType.Xexp = @. real(ρexp*conj(ρexp))
    end
    @inbounds for j in eachindex(_p), i in eachindex(_d)
        neff = _refractive_index_MR(layers[2], [_p[j]], beam_.λ)
        d[i] = _d[i] / mean(real.(neff))
        X = _compute_transfer_matrix(specType, vec([0. d[i] 0.]), [layers[1].n neff layers[3].n], beam_)
        sol[i, j] = _objective_function(oftype, X, specType.Xexp)
    end
    smin = findmin(sol)
    X = _compute_transfer_matrix(
            specType,
            vec([0. d[smin[2][1]] 0.]),
            [layers[1].n _refractive_index_MR(layers[2], [_p[smin[2][2]]], beam_.λ) layers[3].n],
            beam_,
    )
    sol_ = SpaceSolutionEMA(
            sol,
            smin[1],
            _d[smin[2][1]],
            _p[smin[2][2]],
            d[smin[2][1]],
            X,
            specType.Xexp,
            _d,
            _p,
            beam_,
    )
    return sol_
end

"""

    Optimize the parameters of the thin films in a multilayer stack using the selected models and the experimental spectrum.

        sol = fit_tmm_optics(
            specType, xinit, beam, Rexp, layers;
            order=1:length(layers), options=Optim.Options(), alg=SAMIN(),
            lb=0.5.*xinit, ub=1.5.*xinit, σ=ones(size(Xexp)), alpha=false,
            oftype=MeanAbs(),
        )

            specType: Type of spectrum to fit, with accepted values:
                Reflectance(Xexp), Transmittance(Xexp) or Ellipsometry(Xexp).
                For the Reflectance() and Transmittance() you must pass the spectrum array
                Xexp (0<=Xexp<=1) to be fitted. For the Ellipsometry() you have to
                pass the Ψ (first column) and Δ (second column) spectra in degrees, i.e.,
                Xexp = [Ψ Δ].
            xinit: Array with the initial parameters for optimization. For instance, to
                optimise two layers in the system you need to wrap the seeds for the
                algorithm as follow: xinit = [seed1, seed2], even if it is only one
                single layer to fit, xinit = [seed1]. The first parameter of each seed
                must be always the thickness, the rest are the parameters for the selected
                model. The number of seeds must match that of the number of ModelFit layers.
            beam: Structure from PlaneWave.
            layers: Columnwise array of LayerTMMO and ModelFit with information about the
                layers.
                order: Set the order of the layers (build the system) with the layers.
                    This parameter is directly bound to the (previously defined) layers
                    parameter. Order basically contains the indices of each layer inside
                    the parameter layers. Notice that the first and last indices in
                    order cannot be associated to ModelFit layers, as they are the
                    outter media.
                options: Optional Optim.Options structure.
                alg: Optional algorithm method selected, by default takes SAMIN().
                    You can pass the options inside as well, for instance, SAMIN(rt=0.1).
                    Right now, LBFGS(), NelderMead() and SAMIN() are supported from
                    Optim.jl. If you select SAMIN() you need to input lb and ub, otherwise
                    will be set as 0.5 and 1.5 times the xinit argument, respectively.
                    The LBFGS algorithm requires the jacobian for the optimisation but it
                    can be omitted and the algorithm will calculate the differentials
                    numerically.
                lb: Lower bounds for the optimisation variables, by default lb=0.5.*xinit
                ub: Upper bounds for the optimisation variables, by default ub=1.5.*xinit
                σ: Array with the standard deviation of the spectrum for each wavelength,
                   by default is ones(size(Xexp)).
                alpha: Tells the procedure to either use (alpha=true) a linear decrease in
                    the thicknesses of the modelling layers or not (alpha=false). The
                    default is false, do not use alpha.
                oftype: Metric for the objective function to optimise.
                    Can take values: MeanAbs() (default), SumAbs() and SumMeanAbs().

            sol: FitTMMOptics structure with results, with fields
                OptimSolverInfo: Status info from the solver
                spectrumFit: Spectrum obtained with the model and the optimal parameters
                spectrumExp: Input experimental spectrum
                optParams: Array of arrays with optimal parameters
                objfunMin: Optimum value of objective function MSE
                beam: Structure of the light used
                fitOptions: NamedTuple with information of the fitting procedure

"""
function fit_tmm_optics(
    specType::T0,
    xinit,
    beam::T1,
    layers::Array;
    order::AbstractArray{T3,N3}=1:length(layers),
    options=Optim.Options(),
    alg=SAMIN(),
    lb=0.01.*xinit,
    ub=10.0.*xinit,
    σ::Array{T4}=ones(size(specType.Xexp)),
    alpha::T5=false,
    oftype::T7=MeanAbs(),
    SearchRange=(-1.0, 1.0),
    NumDimensions=2,
    FitnessScheme=ScalarFitnessScheme{true}(),
    PopulationSize=50,
    MaxTime=0.0,
    Method=:adaptive_de_rand_1_bin_radiuslimited,
    MaxNumStepsWithoutFuncEvals=100,
    RngSeed=1234,
    MaxFuncEvals=0,
    SaveTrace=false,
    SaveFitnessTraceToCsv=false,
    CallbackInterval=-1.0,
    TargetFitness=nothing,
    TraceMode=:compact,
    MinDeltaFitnessTolerance=1.0e-50,
    FitnessTolerance=1.0e-8,
    TraceInterval=0.5,
    MaxStepsWithoutProgress=10000,
    # CallbackFunction=##80#81(),
    MaxSteps=10000,
    SaveParameters=false,
    SearchSpace=false,
    NumRepetitions=1,
    RandomizeRngSeed=true,
) where{T0<:FitProcedure, T1<:PlaneWave, T3<:Int64, N3, T4<:Real, T5<:Bool, T6<:Symbol, T7<:FitProcedure}
    isa(layers[order[1]], LayerTMMO) || throw("The first layer of the system cannot be ModelFit type. Check the order parameter.")
    isa(layers[order[end]], LayerTMMO) || throw("The last layer of the system cannot be ModelFit type. Check the order parameter.")
    if sum(isa.(layers, ModelFit)) < 1
        error("There should be at least one layer to modelling inside the system (ModelFit type).")
    end
    # Get the parameters for beam
    beam_, _ = _get_beam_parameters(beam)
    # Warm-up
    N = Array{ComplexF64,2}(undef, length(beam_.λ), length(vec(order)))
    d = Array{Float64,1}(undef, length(vec(order)))
    # Check alpha input
    α = alpha ? UseAlpha() : NoAlpha()
    # Build the indices vector for the seeds
    idxSeed = _indicesSeed(layers, vec(order))
    # If it is ellipsometry, then fit norm of ρ
    if isa(specType, Ellipsometry)
        ρexp = @. tan(deg2rad(specType.X[:,1]))*exp(im*deg2rad(specType.X[:,2]))
        specType.Xexp = hcat(real.(ρexp), imag.(ρexp))
        # specType.Xexp = @. real(ρexp*conj(ρexp))
    end
    # Run the solver depending on the algorithm
    if isequal(alg, :BBO)
        optbbo=(
            SearchRange=SearchRange,
            NumDimensions=NumDimensions,
            FitnessScheme=FitnessScheme,
            PopulationSize=PopulationSize,
            MaxTime=MaxTime,
            Method=Method,
            MaxNumStepsWithoutFuncEvals=MaxNumStepsWithoutFuncEvals,
            RngSeed=RngSeed,
            MaxFuncEvals=MaxFuncEvals,
            SaveTrace=SaveTrace,
            SaveFitnessTraceToCsv=SaveFitnessTraceToCsv,
            CallbackInterval=CallbackInterval,
            TargetFitness=TargetFitness,
            TraceMode=TraceMode,
            MinDeltaFitnessTolerance=MinDeltaFitnessTolerance,
            FitnessTolerance=FitnessTolerance,
            TraceInterval=TraceInterval,
            MaxStepsWithoutProgress=MaxStepsWithoutProgress,
            # CallbackFunction=##80#81(),
            MaxSteps=MaxSteps,
            SaveParameters=SaveParameters,
            SearchSpace=SearchSpace,
            NumRepetitions=NumRepetitions,
            RandomizeRngSeed=RandomizeRngSeed,
        )
        solution = _run_fitting_procedure_bbo(
            specType,
            float.(xinit),
            beam_,
            float.(σ),
            layers,
            vec(collect(order)),
            d,
            N,
            idxSeed,
            α,
            oftype,
            optbbo,
        )
    else
        solution = _run_fitting_procedure(
            alg,
            specType,
            float.(xinit),
            beam_,
            float.(σ),
            layers,
            vec(collect(order)),
            options,
            float.(lb),
            float.(ub),
            d,
            N,
            idxSeed,
            α,
            oftype,
        )
    end
    solution_ = FitTMMOptics(
                    solution.X,
                    specType.Xexp,
                    solution.xfinal,
                    solution.solmin,
                    beam_,
                    layers,
                    (
                     fitType=specType,
                     optimAlgorithm=Symbol(summary(solution.sol)),
                     optimSolverInfo=solution.sol,
                    ),
    )
    return solution_
end

## Run the fitting procedure for Optim's solvers depending on input.
function _run_fitting_procedure end

# SAMIN
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:SAMIN, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoNotAlpha, T8<:FitProcedure}
    _checkInput(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(lb),
            flatten_arrays(ub),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer, xinit)
    _multilayer_parameters!(N, d, xfinal, beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# SAMIN, alpha
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:SAMIN, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoAlpha, T8<:FitProcedure}
    _checkInputAlpha(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(lb),
            flatten_arrays(ub),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer[1:end-1], xinit[1:end-1])
    push!(xfinal, [sol.minimizer[end]])
    _multilayer_parameters!(N, d, array_arrays(sol.minimizer[1:end-1], xfinal[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, sol.minimizer[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# NelderMead
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:NelderMead, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoNotAlpha, T8<:FitProcedure}
    _checkInput(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer, xinit)
    _multilayer_parameters!(N, d, xfinal, beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# NelderMead, alpha
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:NelderMead, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoAlpha, T8<:FitProcedure}
    _checkInputAlpha(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer[1:end-1], xinit[1:end-1])
    push!(xfinal, [sol.minimizer[end]])
    _multilayer_parameters!(N, d, array_arrays(sol.minimizer[1:end-1], xfinal[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, sol.minimizer[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# LBFGS
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:LBFGS, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoNotAlpha, T8<:FitProcedure}
    _checkInput(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer, xinit)
    _multilayer_parameters!(N, d, xfinal, beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# LBFGS, alpha
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:LBFGS, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoAlpha, T8<:FitProcedure}
    _checkInputAlpha(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer[1:end-1], xinit[1:end-1])
    push!(xfinal, [sol.minimizer[end]])
    _multilayer_parameters!(N, d, array_arrays(sol.minimizer[1:end-1], xfinal[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, sol.minimizer[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# NewtonTrustRegion
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:NewtonTrustRegion, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoNotAlpha, T8<:FitProcedure}
    _checkInput(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer, xinit)
    _multilayer_parameters!(N, d, xfinal, beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# NewtonTrustRegion, alpha
function _run_fitting_procedure(
    alg::T1,
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    options,
    lb::Array{Array{T3,1},1},
    ub::Array{Array{T3,1},1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
) where {T1<:NewtonTrustRegion, T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoAlpha, T8<:FitProcedure}
    _checkInputAlpha(xinit, layers)
    sol = optimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype),
            flatten_arrays(xinit),
            alg,
            options,
    )
    xfinal = array_arrays(sol.minimizer[1:end-1], xinit[1:end-1])
    push!(xfinal, [sol.minimizer[end]])
    _multilayer_parameters!(N, d, array_arrays(sol.minimizer[1:end-1], xfinal[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, sol.minimizer[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.minimum)
    return sol_
end

# BBO
function _run_fitting_procedure_bbo(
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
    optbbo,
) where {T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoNotAlpha, T8<:FitProcedure}
    _checkInput(xinit, layers)
    sol = bboptimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype);
            SearchRange=optbbo.SearchRange,
            NumDimensions=optbbo.NumDimensions,
            FitnessScheme=optbbo.FitnessScheme,
            PopulationSize=optbbo.PopulationSize,
            MaxTime=optbbo.MaxTime,
            Method=optbbo.Method,
            MaxNumStepsWithoutFuncEvals=optbbo.MaxNumStepsWithoutFuncEvals,
            RngSeed=optbbo.RngSeed,
            MaxFuncEvals=optbbo.MaxFuncEvals,
            SaveTrace=optbbo.SaveTrace,
            SaveFitnessTraceToCsv=optbbo.SaveFitnessTraceToCsv,
            CallbackInterval=optbbo.CallbackInterval,
            TargetFitness=optbbo.TargetFitness,
            TraceMode=optbbo.TraceMode,
            MinDeltaFitnessTolerance=optbbo.MinDeltaFitnessTolerance,
            FitnessTolerance=optbbo.FitnessTolerance,
            TraceInterval=optbbo.TraceInterval,
            MaxStepsWithoutProgress=optbbo.MaxStepsWithoutProgress,
            # CallbackFunction=##80#81(),
            MaxSteps=optbbo.MaxSteps,
            SaveParameters=optbbo.SaveParameters,
            SearchSpace=optbbo.SearchSpace,
            NumRepetitions=optbbo.NumRepetitions,
            RandomizeRngSeed=optbbo.RandomizeRngSeed,
    )
    xfinal = array_arrays(sol.archive_output.best_candidate, xinit)
    _multilayer_parameters!(N, d, xfinal, beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.archive_output.best_fitness)
    return sol_
end

# BBO, alpha
function _run_fitting_procedure_bbo(
    specType::T2,
    xinit::Array{Array{T3,1},1},
    beam::T4,
    σ::Array{T3},
    layers::Array,
    order::Array{T5,1},
    d::Array{T3,1},
    N::Array{T6,2},
    idxSeed::Array{T5,1},
    alpha::T7,
    oftype::T8,
    optbbo,
) where {T2<:FitProcedure, T3<:Float64, T4<:PlaneWave, T5<:Int64, T6<:ComplexF64, T7<:DoAlpha, T8<:FitProcedure}
    _checkInputAlpha(xinit, layers)
    sol = bboptimize(
            x->_fit_objfun(x, specType, beam, σ, layers, order, xinit, d, N, idxSeed, alpha, oftype);
            SearchRange=optbbo.SearchRange,
            NumDimensions=optbbo.NumDimensions,
            FitnessScheme=optbbo.FitnessScheme,
            PopulationSize=optbbo.PopulationSize,
            MaxTime=optbbo.MaxTime,
            Method=optbbo.Method,
            MaxNumStepsWithoutFuncEvals=optbbo.MaxNumStepsWithoutFuncEvals,
            RngSeed=optbbo.RngSeed,
            MaxFuncEvals=optbbo.MaxFuncEvals,
            SaveTrace=optbbo.SaveTrace,
            SaveFitnessTraceToCsv=optbbo.SaveFitnessTraceToCsv,
            CallbackInterval=optbbo.CallbackInterval,
            TargetFitness=optbbo.TargetFitness,
            TraceMode=optbbo.TraceMode,
            MinDeltaFitnessTolerance=optbbo.MinDeltaFitnessTolerance,
            FitnessTolerance=optbbo.FitnessTolerance,
            TraceInterval=optbbo.TraceInterval,
            MaxStepsWithoutProgress=optbbo.MaxStepsWithoutProgress,
            # CallbackFunction=##80#81(),
            MaxSteps=optbbo.MaxSteps,
            SaveParameters=optbbo.SaveParameters,
            SearchSpace=optbbo.SearchSpace,
            NumRepetitions=optbbo.NumRepetitions,
            RandomizeRngSeed=optbbo.RandomizeRngSeed,
    )
    xfinal = array_arrays(sol.archive_output.best_candidate[1:end-1], xinit[1:end-1])
    push!(xfinal, [sol.archive_output.best_candidate[end]])
    _multilayer_parameters!(N, d, array_arrays(sol.archive_output.best_candidate[1:end-1], xfinal[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, sol.archive_output.best_candidate[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    sol_ = (sol=sol, X=X, xfinal=xfinal, solmin=sol.archive_output.best_fitness)
    return sol_
end

# Check the input sizes and the seed with alpha.
function _checkInput(xinit, layers)
    length(xinit) == length(findall(isa.(layers, ModelFit))) || throw("There should be seeds for all active layers (ModelFit).")
    return nothing
end

function _checkInputAlpha(xinit, layers)
    length(xinit) == length(findall(isa.(layers, ModelFit)))+1 || throw("There should be seeds for all active layers (ModelFit) plus the last element with the alpha.")
    return nothing
end

## Fit objective function used for the optimization process.
function _fit_objfun end

# No alpha
function _fit_objfun(
    x::Array{T1,1},
    specType::T2,
    beam::T3,
    σ::Array{T1},
    layers::Array,
    order::Array{T4,1},
    xinit::Array{Array{T1,1},1},
    d::Array{T1,1},
    N::Array{T5,2},
    idxSeed::Array{T4,1},
    alpha::T6,
    oftype::T7,
) where {T1<:Float64, T2<:FitProcedure, T3<:PlaneWave, T4<:Int64, T5<:ComplexF64, T6<:DoNotAlpha, T7<:FitProcedure}
    _multilayer_parameters!(N, d, array_arrays(x, xinit), beam, layers, order, idxSeed)
    X = _compute_transfer_matrix(specType, d, N, beam)
    mse = _objective_function(oftype, X, specType.Xexp; σ=σ)
    return mse
end

# Alpha
function _fit_objfun(
    x::Array{T1,1},
    specType::T2,
    beam::T3,
    σ::Array{T1},
    layers::Array,
    order::Array{T4,1},
    xinit::Array{Array{T1,1},1},
    d::Array{T1,1},
    N::Array{T5,2},
    idxSeed::Array{T4,1},
    alpha::T6,
    oftype::T7,
) where {T1<:Float64, T2<:FitProcedure, T3<:PlaneWave, T4<:Int64, T5<:ComplexF64, T6<:DoAlpha, T7<:FitProcedure}
    _multilayer_parameters!(N, d, array_arrays(x[1:end-1], xinit[1:end-1]), beam, layers, order, idxSeed)
    _monotonic_lin!(d, x[end])
    X = _compute_transfer_matrix(specType, d, N, beam)
    mse = _objective_function(oftype, X, specType.Xexp; σ=σ)
    return mse
end

## Compute the transfer matrix and return the specified spectra.
function _compute_transfer_matrix end

# Transmittance
function _compute_transfer_matrix(
    specType::T0, d::Array{T1,1}, N::Array{T2,2}, beam::T3,
) where {T0<:Transmittance, T1<:Float64, T2<:ComplexF64, T3<:PlaneWave}
    spectra = transfer_matrix(N, d, beam)[1]
    avg = average_polarisation(beam.p, spectra.Tp, spectra.Ts)
    return avg
end

# Reflectance
function _compute_transfer_matrix(
    specType::T0, d::Array{T1,1}, N::Array{T2,2}, beam::T3,
) where {T0<:Reflectance, T1<:Float64, T2<:ComplexF64, T3<:PlaneWave}
    spectra = transfer_matrix(N, d, beam)[1]
    avg = average_polarisation(beam.p, spectra.Rp, spectra.Rs)
    return avg
end

# Ellipsometry
function _compute_transfer_matrix(
    specType::T0, d::Array{T1,1}, N::Array{T2,2}, beam::T3,
) where {T0<:Ellipsometry, T1<:Float64, T2<:ComplexF64, T3<:PlaneWave}
    spectra = transfer_matrix(N, d, beam)[1]
    ρ::Array{ComplexF64,1} = vec(spectra.ρp./spectra.ρs)
    return hcat(-real.(ρ), imag.(ρ))
    # ρnorm = @. real(ρ*conj(ρ))
    # return ρnorm
end

## Metrics (objective function) for the optimisation process.

# MeanAbs() = MeanAbsObjFun
function _objective_function(
    oftype::T0, X::Array{T1}, Xexp::Array{T1};
    σ::Array{T1}=ones(size(Xexp)),
) where {T0<:MeanAbsObjFun, T1<:Float64}
    _absq = abs2.(X .- Xexp)./σ
    _mse = mean(_absq, dims=1)
    objfun = sum(_mse)
    return objfun
end

# SumAbs() = SumAbsObjFun
function _objective_function(
    oftype::T0, X::Array{T1}, Xexp::Array{T1};
    σ::Array{T1}=ones(size(Xexp)),
) where {T0<:SumAbsObjFun, T1<:Float64}
    _absq = abs2.(X .- Xexp)./σ
    _mse = sum(_absq, dims=1)
    objfun = sum(_mse)
    return objfun
end

# SumMeanAbs() = SumMeanAbsObjFun
function _objective_function(
    oftype::T0, X::Array{T1}, Xexp::Array{T1};
    σ::Array{T1}=ones(size(Xexp)),
) where {T0<:SumMeanAbsObjFun, T1<:Float64}
    _absq = abs2.(X .- Xexp)./σ
    _mse = mean(_absq, dims=1)
    _sse = sum(_absq, dims=1)
    # objfun = sqrt(sum(_mse + _sse)/size(Xexp, 1))
    objfun = sum(_mse + _sse)/size(Xexp, 1)
    return objfun
end

## Returns the index of refraction for the selected model. RI is imported from RefractiveIndicesModels.jl.
function _refractive_index_MR(
    layer::T1, argin::AbstractArray{T2,1}, λ::Array{T2,1},
) where {T1<:ModelFit, T2<:Float64}
    if layer.model == :bruggeman
        _errorParametersEMA(layer)
        return RI.bruggeman(argin[1], [layer.N.ninc layer.N.nhost]; df=layer.df)
    elseif layer.model == :looyenga
        _errorParametersEMA(layer)
        return RI.looyenga([argin[1] 1.0-argin[1]], [layer.N.ninc layer.N.nhost]; df=layer.df)
    elseif layer.model == :monecke
        _errorParametersEMA(layer)
        return RI.monecke(argin[1], [layer.N.ninc layer.N.nhost]; dfm=layer.dfm, c=layer.c)
    elseif layer.model == :maxwell
        _errorParametersEMA(layer)
        return RI.maxwell([argin[1] 1.0-argin[1]], [layer.N.ninc layer.N.nhost]; df=layer.df)
    elseif layer.model == :gedf
        _errorParametersEMA(layer)
        return RI.gedf(argin, [layer.N.ninc layer.N.nhost])
    elseif layer.model == :gem
        _errorParametersEMA(layer)
        return RI.gem(argin, [layer.N.ninc layer.N.nhost])
    elseif layer.model == :lorentzlorenz
        _errorParametersEMA(layer)
        return RI.lorentz_lorenz(argin[1], [layer.N.ninc layer.N.nhost])
    elseif layer.model == :sellmeier
        return RI.sellmeier(argin, λ./1e3) # Sellmeier eq takes λ in μm
    elseif layer.model == :cauchyurbach
    	return RI.cauchy_urbach(argin, λ)
    elseif layer.model == :drudelorentz
        argin_ = [argin[1], argin[2:end]]
        x = _oscillators_input(argin, 3)
        return RI.drude_lorentz(x, 1240.0./λ)
    elseif layer.model == :tauclorentz
        argin_ = [argin[1:2], argin[3:end]]
        x = _oscillators_input(argin_, 3)
    	return RI.tauc_lorentz(x, 1240.0./λ)
    elseif layer.model == :forouhibloomer
        argin_ = [argin[1:2], argin[3:end]]
        x = _oscillators_input(argin_, 3)
        return RI.forouhi_bloomer(x, 1240.0./λ)
    elseif layer.model == :lorentzplasmon
        argin_ = [argin[1:3], argin[4:end]]
        x = _oscillators_input(argin_, 3)
        return RI.lorentz_plasmon(x, 1240.0./λ)
    elseif layer.model == :codylorentz
        argin_ = [argin[1:5], argin[6:end]]
        x = _oscillators_input(argin_, 3)
        return RI.cody_lorentz(x, 1240.0./λ; cme=layer.cme)
    end
end

## Error if EMA models do not have the two columns indices of refraction.
function _errorParametersEMA(layer::ModelFit)
    !|(isempty(layer.N[1]), isempty(layer.N[2])) || throw("To use the $(layer.model) model you need to specify the two indices of refraction of the components in ModelFit.")
    return nothing
end

## Construct the multilayers parameters to optimise depending on their type (LayerTMMO, ModelFit).
function _multilayer_parameters!(
    N::Array{T1,2},
    d::Array{T2,1},
    x::Array{Array{T2,1},1},
    beam::T3,
    layers::Array,
    order::Array{T4,1},
    idxSeed::Array{T4,1},
) where{T1<:ComplexF64, T2<:Float64, T3<:PlaneWave, T4<:Int64}
    k::Int64 = 1
    for i in eachindex(order)
        idx = order[i]
        if isa(layers[idx], LayerTMMO)
            d[i] = layers[idx].d
            @views N[:,i] = layers[idx].n
        elseif isa(layers[idx], ModelFit)
            d[i] = x[idxSeed[k]][1] # the thickness is defined first
            @views N[:,i] = _refractive_index_MR(layers[idx], x[idxSeed[k]][2:end], beam.λ)
            k += 1
        end
    end
    return nothing
end

## Build the indices vector for the seeds
function _indicesSeed(layers, order)
    mflayers = findall(isa.(layers, ModelFit))
    idx = 1:length(mflayers)
    orderSeed = []
    for i in eachindex(order), j in eachindex(mflayers)
        if order[i] == mflayers[j]
            push!(orderSeed, idx[j])
        end
    end
    return Int.(orderSeed)
end

## Linear decrease of thickness.
function _monotonic_lin!(d::Array{T1,1}, α::T1) where {T1<:Float64}
    # d .*= α.^(1:length(d))
    # Only for the second active layer until the one before the substrate
    d .*= [1.0; 1.0; α.^(3:length(d)-1); 1.0]
    return nothing
end

end # module


text:
#module ModelParamsTest

#using ModelParams
#using Random, Test

include("object_id_test.jl")
include("parameters_test.jl")
include("param_vector_test.jl")
include("deviation_test.jl")

"""
## Test model
More extensive testing in `SampleModel`
"""
mutable struct Obj1 <: ModelObject
    objId :: ObjectId
    x :: Float64
    y :: Vector{Float64}
    z :: Array{Float64,2}
    pvec :: ParamVector
end

# function Obj1(x, y, z)
#     return Obj1(x, y, z, ParamVector(ObjectId(:pv1)))
# end

function init_obj1()
    objId = ObjectId(:obj1);
    px = Param(:x, "x obj1", "x1", 11.1, 9.9, 1.1, 99.9, true);
    valueY = [1.1, 2.2];
    py = Param(:y, "y obj1", "y1", valueY, valueY .+ 1.0,
        valueY .- 5.0, valueY .+ 5.0, true);
    valueZ = [3.3 4.4; 5.5 7.6];
    pz = Param(:z, "z obj1", "z1", valueZ, valueZ .+ 1.0,
        valueZ .- 5.0, valueZ .+ 5.0, false);
    pvector = ParamVector(objId, [px, py, pz])
    o1 = Obj1(objId, px.value, py.value, pz.value, pvector);
    ModelParams.set_own_values_from_pvec!(o1, pvector, true);
    return o1
end

mutable struct Obj2 <: ModelObject
    objId :: ObjectId
    a :: Float64
    y :: Float64
    b :: Array{Float64,2}
    pvec :: ParamVector
end

# function Obj2(a, y, b)
#     return Obj2(a, y, b, ParamVector(ObjectId(:pv1)))
# end

function init_obj2()
    objId = ObjectId(:obj2);
    pa = Param(:a, "a obj2", "a2", 12.1, 7.9, -1.1, 49.9, true);
    valueY = 9.4;
    py = Param(:y, "y obj2", "y2", valueY, valueY .+ 1.0,
        valueY .- 5.0, valueY .+ 5.0, true);
    valueB = 2.0 .+ [3.3 4.4; 5.5 7.6];
    pb = Param(:b, "b obj2", "b2", valueB, valueB .+ 1.0,
        valueB .- 5.0, valueB .+ 5.0, true);
    pvector = ParamVector(objId, [pa, py, pb]);
    o2 = Obj2(objId, pa.value, py.value, pb.value, pvector);
    ModelParams.set_own_values_from_pvec!(o2, pvector, true);
    return o2
end

mutable struct TestModel <: ModelObject
    objId :: ObjectId
    o1 :: Obj1
    # pvec1 :: ParamVector
    o2 :: Obj2
    # pvec2 :: ParamVector
    a :: Float64
    y :: Float64
end

function init_test_model()
    o1 = init_obj1();
    o2 = init_obj2();
    return TestModel(ObjectId(:testModel), o1, o2, 9.87, 87.73)
end

"""
## Simulates the workflow for calibration

More extensive testing in `SampleModel`
"""
function modelTest()
    m = init_test_model()

    ModelParams.report_params(m, true);
    ModelParams.report_params(m, false);
    nParam, nElem = ModelParams.n_calibrated_params(m; isCalibrated = true);
    @test nParam > 1
    @test nElem > nParam
    nParam, nElem = ModelParams.n_calibrated_params(m; isCalibrated = false);
    @test nParam >= 1
    @test nElem > nParam

    isCalibrated = true;

    # Sync calibrated model values with param vector
    ModelParams.set_own_values_from_pvec!(m.o1, m.o1.pvec, isCalibrated);
    # also sync the non-calibrated default values
    ModelParams.set_default_values!(m.o1, false);
    @test ModelParams.check_own_calibrated_params(m.o1, m.o1.pvec);
    @test ModelParams.check_own_fixed_params(m.o1, m.o1.pvec);

    # The same in one step
    ModelParams.sync_own_values!(m.o2, m.o2.pvec);
    @test ModelParams.check_own_calibrated_params(m.o2, m.o2.pvec);
    @test ModelParams.check_own_fixed_params(m.o2, m.o2.pvec);

    # This step just for testing
    # These are the values that we expect to get back in the end
    d1 = make_dict(m.o1.pvec, isCalibrated);
    d2 = make_dict(m.o2.pvec, isCalibrated);

    # For each model object: make vector of param values
    vv1 = make_guess(m.o1.pvec);
    @test isa(ModelParams.values(vv1), Vector{Float64})
    vv2 = make_guess(m.o2.pvec);
    @test isa(ModelParams.values(vv2), Vector{Float64})

    # This is passed to optimizer as guess
    vAll = [ModelParams.values(vv1); ModelParams.values(vv2)];
    @test isa(vAll, Vector{Float64})

    # Same in one step for all param vectors
    vv = ModelParams.make_guess([m.o1.pvec, m.o2.pvec]);
    @test vAll ≈ ModelParams.values(vv)
    @test ModelParams.lb(vv) ≈ [ModelParams.lb(vv1); ModelParams.lb(vv2)]
    @test ModelParams.ub(vv) ≈ [ModelParams.ub(vv1); ModelParams.ub(vv2)]

    # Now we run the optimizer, which changes `vAll`

    # In the objective function: the guess is reassembled
    # into dicts which are then put into the objects

    # Using in a single convenience method
    vOut = ModelParams.sync_from_vector!([m.o1, m.o2], vAll);
    @test isempty(vOut);
    @test ModelParams.check_own_calibrated_params(m.o1, m.o1.pvec);
    @test ModelParams.check_own_fixed_params(m.o1, m.o1.pvec);
    @test ModelParams.check_own_calibrated_params(m.o2, m.o2.pvec);
    @test ModelParams.check_own_fixed_params(m.o2, m.o2.pvec);

    # The same step-by-step. Only needed for testing
    d11, nUsed1 = vector_to_dict(m.o1.pvec, vAll, isCalibrated);
    @test d11[:x] == d1[:x]
    @test d11[:y] == d1[:y]
    # copy into param vector; then sync with model object
    ModelParams.set_values_from_dict!(m.o1.pvec, d11);
    ModelParams.set_own_values_from_pvec!(m.o1, m.o1.pvec, isCalibrated);
    @test m.o1.x ≈ d1[:x]
    @test m.o1.y ≈ d1[:y]

    # The same in a single convenience function, one object at a time
    # Just for testing
    nUsed11 = ModelParams.sync_from_vector!(m.o1, vAll);
    @test nUsed11 == nUsed1
    @test m.o1.x ≈ d1[:x]
    @test m.o1.y ≈ d1[:y]
    deleteat!(vAll, 1 : nUsed1);

    d22, nUsed2 = ModelParams.vector_to_dict(m.o2.pvec, vAll, isCalibrated);
    @test d22[:a] == d2[:a]
    @test d22[:b] == d2[:b]
    ModelParams.set_values_from_dict!(m.o2.pvec, d22);
    ModelParams.set_own_values_from_pvec!(m.o2, m.o2.pvec, isCalibrated);
    @test m.o2.a ≈ d2[:a]
    @test m.o2.b ≈ d2[:b]
    # Last object: everything should be used up
    @test length(vAll) == nUsed2

    # Test changing parameters
    d22[:a] = 59.34;
    ModelParams.set_own_values_from_dict!(m.o2, d22);
    @test m.o2.a ≈ d22[:a]

    d22[:b] .+= 3.8;
    ModelParams.set_values_from_dict!(m.o2.pvec, d22);
    ModelParams.set_own_values_from_pvec!(m.o2, m.o2.pvec, isCalibrated);
    @test m.o2.b ≈ d22[:b]

    return true
end

"""
    merge_object_arrays!
"""

mutable struct MoaTest
    x
    y
    z
    a
end

function merge_object_arrays_test()
    n = 17;
    Random.seed!(123);
    x = rand(n, 2);
    y = rand(n, 3);
    z = rand(n, 3, 2);
    a = rand(n);
    oTg = MoaTest(x, y, z, a);

    idxV = [3,4,8];
    nSrc = length(idxV);
    x2 = rand(nSrc, 2);
    y2 = 9.3;
    z2 = rand(nSrc, 3, 2);
    a2 = rand(nSrc);
    oSrc = MoaTest(x2, y2, z2, a2);

    ModelParams.merge_object_arrays!(oSrc, oTg, idxV, false, dbg = true);

    @test oTg.x[idxV,:] ≈ x2
    @test oTg.y ≈ y  # Skipped field; not array in source
    @test oTg.z[idxV,:,:] ≈ z2
    @test oTg.a[idxV] ≈ a2
    return true
end

end # module


text:
using Base.Test

# choose the next block for redistribution
# @return the value and index of the block to distribute (1 based)
function choose_block(x)
  findmax(x)
end

# redistribute the max block
# @return the redistributed memory
function redistribute(x)
  val, index = choose_block(x)
  result = x
  result[index] = 0 # zero redistributed block
  len = size(x)[1]
  # we distribute one-by-one until nothing is left
  for i = 1:val
    index = (index == len) ? 1 : (index + 1) # wrap around at last index
    result[index] += 1
  end
  result
end

# count redistributions until cycle is detected
# @return the number of redistributions overall, and the cycle length
function count_until_cycle(x)
  seen = Dict()
  n = 0
  next = x
  # we use the string representation of the array as key to avoid aliasing problems
  while !haskey(seen, "$next")
    seen["$next"] = n
    next = redistribute(next)
    n += 1
  end
  (n, n - seen["$next"])
end

# tests for choose_block
@test choose_block([0, 2, 7, 0]) == (7,3)
# choose_block should choose the first one if multiple max-value blocks exist
@test choose_block([0, 7, 0, 7]) == (7,2)

# tests for redistribution
@test redistribute([0, 0, 0, 1]) == [1, 0, 0, 0]
@test redistribute([0, 2, 7, 0]) == [2, 4, 1, 2]

# tests for detecting cycles
@test count_until_cycle([0, 2, 7, 0]) == (5, 4)

# solve it
input = [11,11,13,7,0,15,5,5,4,4,1,1,7,1,15,11]
println(count_until_cycle(input))


text:
using StatsPlots

# ------------------------------------------------------------------------------------------------ #
#                 Compare the distributions of features in a feature matrices                      #
# ------------------------------------------------------------------------------------------------ #

@userplot FeatureViolin
@recipe function f(H::FeatureViolin; normalise=:feature, normtype=:mix, rightannotations=nothing)
    # if normalise is on, normalise distributions to the range of the first feature matrix but centre uniquely
    F = H.args
    if length(F) == 1
        F = F[1]
        F[isinf.(F)] .= 0
        F2 = nothing
    elseif length(F) == 2
        F, F2 = F
        F[isinf.(F)] .= 0
        F2[isinf.(F2)] .= 0
        if F isa AbstractFeatureArray && F2 isa AbstractFeatureArray
            if size(F, 1) > size(F2, 1)
                (F2, F) = intersectFeatures(F2, F)
            else
                (F, F2) = intersectFeatures(F, F2)
            end
        end
    end
    if normtype == :zscore
        f1 = x -> StatsBase.std(x)
        f2 = x -> mean(x)
    elseif normtype == :unit
        f1 = x -> abs(-(extrema(x)...))
        f2 = x -> mean(extrema(x))
    elseif normtype == :mix
        f1 = x -> abs(-(extrema(x)...))
        f2 = x -> mean(x)
    elseif normtype == :none
        f1 = x -> 1.0
        f2 = x -> 0.0
    end
    if (normalise == :feature) && !isnothing(F2)
        S = mapslices(f1, hcat(F, F2), dims=2)
        C = mapslices(f2, hcat(F, F2), dims=2)
        S[S .== 0] .= 1.0
        F = (F.-C)./S
        F2 = (F2.-C)./S
        yaxis --> nothing
    elseif isnothing(F2) || normalise
        S = mapslices(f1, F, dims=2)
        C = mapslices(f2, F, dims=2)
        S[S .== 0] .= 1.0
        F = (F.-C)./S
        if ~isnothing(F2)
            F2 = (F2.- mapslices(x -> mean(x), F2, dims=2))./S
        end
    end

    if F isa AbstractFeatureArray && F2 isa AbstractFeatureArray
        fnames = replace.(String.(Catch22.featureDims(F)),  '_'=>"\\_")
    else
        fnames = 1:size(F, 1)
    end

    if backend() !== Plots.PyPlotBackend()
        rotation --> 90
        bottom_margin --> 10Plots.cm
    else
        tickfontrotation --> 90
    end
    size --> (700, 900)
    legend := false

    if ~isnothing(F2)
        @series begin
            seriestype := :violin
            side := :right
            seriescolor := :crimson
            ((1:size(F2, 1))', Array(F2)')
        end
    end

    @series begin
        seriestype := :violin
        if ~isnothing(rightannotations)
            top_margin --> 40Plots.mm
            annotations := [(i, 1.1, rightannotations[i]) for i ∈ 1:length(rightannotations)]
        end
        #bottom_margin --> 70Plots.mm
        #top_margin --> 70Plots.mm
        tickfontvalign := :top
        framestyle --> :box
        ymirror --> true
        if ~isnothing(F2); side := :left; seriescolor := :cornflowerblue; end
        xticks --> (1:size(F, 1), fnames)
        ((1:size(F, 1))', Array(F)')
    end

end

# ------------------------------------------------------------------------------------------------ #
#                                Plot the baseline variance mapping                                #
# ------------------------------------------------------------------------------------------------ #
@userplot VarianceMapping
@recipe function f(P::VarianceMapping; interval=rampInterval, names=false)
    F = P.args[1]
    Fₗ = P.args[2]
    Fₕ = P.args[3]
    legend --> false
    #link := :x
    grid --> true
    yticks --> [0.0, 1.0]
    xguide --> "Input Variance"
    yguide --> "Output Variance"
    left_margin --> 10Plots.mm

    if any(Catch22.featureDims(Fₗ) .!= Catch22.featureDims(Fₕ))
        error("High and low dimensional baselines do not have the same features")
    end

    𝛔ₗ, 𝛔ₕ = std(Fₗ, dims=2), std(Fₕ, dims=2)
    𝛔ₕ[𝛔ₕ .< 𝛔ₗ] .= Inf
    fnames = Catch22.featureDims(F)
    𝐟 = interval.(Fₗ, Fₕ, F)
    𝐟 = Catch22.featureVector(𝐟, Catch22.featureDims(Fₗ))
    (F, 𝐟) = intersectFeatures(F, 𝐟)

    xx = -0.1:0.01:1.1
    f₁ = interval(0.0, 1.0)
    yy = f₁.(xx)
    𝛔ᵣ = [𝐟[s, 1](F[s, :]) for s ∈ fnames]

    @series begin
        seriescolor --> :black
        linewidth --> 3
        (xx, yy)
    end

    @series begin
        #clims := extrema(vec(𝛔ᵣ))
        marker_z --> vec(𝛔ᵣ)
        markersize --> 10
        tickfontsize --> 12
        seriescolor --> cgrad(:RdGy_4, rev=true)
        seriestype --> :scatter
        xticks --> ((0.0, 1.0), ("σₗ", "σₕ"))
        (d, 𝛔ₗ′) = intersectFeatures(𝛔, 𝛔ₗ)
        (d, 𝛔ₕ′) = intersectFeatures(𝛔, 𝛔ₕ)
        x = ((vec(𝛔) .- vec(𝛔ₗ′))./(vec(𝛔ₕ′) .- vec(𝛔ₗ′)))
        y = vec(𝛔ᵣ)
        if names
            ftsz = [1+Int(round(y[i]./max(y...))).*8 for i ∈ 1:lastindex(fnames)]
            annotations := [(x[i]+0.05, y[i], text("$(fnames[i])", :black, :left, ftsz[i])) for i ∈ 1:lastindex(fnames)]
        end
        (x, y)
    end

end

@userplot IntervalScaling
@recipe function f(P::IntervalScaling; interval=rampInterval, reftoramp=true)
    if length(P.args) > 0 # It's evolving, but backwards
        σₗ = P.args[1]
    else
        σₗ = 0.0
    end
    if length(P.args) > 1
        σₕ = P.args[2]
    else
        σₕ = π
    end
    fl = rand(1, 1000).*σₗ
    fh = rand(1, 1000).*σₕ

    𝛔 = LinRange(σₗ, σₕ, 10000)
    f = [rand(1, 20).*𝛔[i] for i ∈ 1:length(𝛔)]
    x = StatsBase.std.(f)
    y = [interval(fl, fh, f[i])[1](f[i]) for i ∈ 1:length(x)]
    if reftoramp
        @series begin
            seriestype := :line
            label --> nothing#"Ramp Interval"
            linewidth --> 2.5
            seriescolor --> :black
            (x, [rampInterval(fl, fh, f[i])[1](f[i]) for i ∈ 1:length(x)])
        end
    end
    @series begin
        seriestype := :line
        seriescolor := crimson
        linewidth --> 2.5
        label --> nothing
        xguide --> "σ"
        yguide --> "σ̂"
        (x, y)
    end
end


text:
#=

This file specifies how to create surface layers, given the surface type, and related info

=#

"""
    $(FUNCTIONNAME)(lambertian::LambertianSurfaceScalar{FT})

Computes (in place) surface optical properties for a (scalar) lambertian albedo as [`AddedLayer`](@ref)

    - `lambertian` a [`LambertianSurfaceScalar`](@ref) struct that defines albedo as scalar
    - `SFI` bool if SFI is used
    - `m` Fourier moment (starting at 0)
    - `pol_type` Polarization type struct
    - `quad_points` Quadrature points struct
    - `τ_sum` total optical thickness from TOA to the surface
    - `architecture` Compute architecture (GPU,CPU)
"""
function create_surface_layer!(lambertian::LambertianSurfaceScalar{FT},
                               added_layer::AddedLayer,
                               SFI,
                               m::Int,
                               pol_type,
                               quad_points,
                               τ_sum,
                               architecture) where {FT}

    if m == 0
        @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart, iμ₀, μ₀ = quad_points
        arr_type = array_type(architecture)
        # Albedo normalized by π (and factor 2 for 0th Fourier Moment)
        ρ = 2lambertian.albedo/FT(π)
        # Get size of added layer
        dim = size(added_layer.r⁻⁺)
        Nquad = dim[1] ÷ pol_type.n
        # Ensure matrices are zero:
        tmp = zeros(Real, pol_type.n)
        tmp[1] = ρ
        R_surf = Array(Diagonal(vcat(ρ, zeros(pol_type.n-1))))
        R_surf = repeat(R_surf',Nquad)
        R_surf = repeat(R_surf',Nquad)
        # Move to architecture:
        R_surf = arr_type(R_surf)
        if SFI
            I₀_NquadN = similar(qp_μN);
            I₀_NquadN[:] .=0;
            i_start   = pol_type.n*(iμ₀-1) + 1
            i_end     = pol_type.n*iμ₀;
            I₀_NquadN[i_start:i_end] = pol_type.I₀;

            added_layer.J₀⁺[:] .= 0
            added_layer.J₀⁻[:,1,:] = μ₀*(R_surf*I₀_NquadN) .* exp.(-τ_sum/μ₀)';
            #repeat(repeat(tmp .* exp.(-τ_sum/qp_μN[i_start]), Nquad), 1,1,dim[3])
        end
        R_surf = R_surf * Diagonal(qp_μN.*wt_μN)
        tmp    = ones(pol_type.n*Nquad)
        T_surf = arr_type(Diagonal(tmp))

        added_layer.r⁻⁺ .= R_surf;
        added_layer.r⁺⁻ .= 0;
        added_layer.t⁺⁺ .= T_surf;
        added_layer.t⁻⁻ .= T_surf;

    else
        added_layer.r⁻⁺[:] .= 0;
        added_layer.r⁻⁺[:] .= 0;
        added_layer.t⁺⁺[:] .= 0;
        added_layer.t⁻⁻[:] .= 0;
        added_layer.J₀⁺[:] .= 0;
        added_layer.J₀⁻[:] .= 0;
    end
end


text:
using PLext

#for i=1:20
#    sv = PLext.StateVector(i)
#    PLext.applyPauliX(sv, i, false)
#    PLext.applyPauliX(sv, i, true)
#    PLext.applyHadamard(sv, i, false)
#    PLext.applyHadamard(sv, i, true)
#end
#

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
module CellArrayReindexTests

using Test
using Gridap
using Gridap.CellValuesGallery
using ..CellValuesMocks

a = [[1,2,4,5],[2,4,6,7],[4,3,5,1],[2,3]]
ca0 = CellValueFromArray(a)
ca = -ca0
b = [3,1,2]
cb = CellValueFromArray(b)

cc = reindex(ca,cb)

d = -a[b]
test_index_cell_array(cc,d)

z = [2,3,1]
cz = CellValueFromArray(z)

cc2 = reindex(ca,cz)

d2 = -a[z]

# This tests that ci and ci2 are properly cached
for (k,(ci,ci2)) in enumerate(zip(cc,cc2))
  @test ci == d[k]
  @test ci2 == d2[k]
end

b = [3,3,3,3]
cb = TestIterCellValue(3,4)

cc = reindex(ca,cb)
d = -a[b]
test_iter_cell_array(cc,d)

end # module