text:
# AdaBoost ensemble module
module AdaBoost

	using j4pr.ClassifierCombiner
	using j4pr: countapp, ohenc_integer
	using Distances: colwise, Euclidean
	using MLDataPattern: nobs, getobs, datasubset, undersample, targets
	using StatsBase: sample!,weights
	using StatsFuns: softmax!
	using UnicodePlots: lineplot
	export AdaBoostEnsemble, adaboost_train, adaboost_exec

	abstract type BoostType end
	struct AdaBoostM1 <: BoostType end
	struct AdaBoostM2 <: BoostType end

	"""
	AdaBoost ensemble object (it describes the ensemble).

	# Fields
	 * `weights::Vector{Float64}` is a vector weights associated to each ensemble member
	 * `f_train::Function` is the training function
	 * `f_exec::Function` is the execution function
	 * `members::Vector` is a vector containing the trained classifiers
	 * `boost_type::AdaBoost.BoostType` indicates the type of AdaBoost approach used
	 * `combiner::ClassifierCombiner.AbstractCombiner` the means by which ensemble outputs are combined
	"""
	struct AdaBoostEnsemble{F, T<:Vector{S} where S, B<:BoostType, U<:AbstractCombiner}
		weights::Vector{Float64}	# vector of ensemble weights
		f_train::Function		# training function
		f_exec::F			# execution function
		members::T			# ensemble members (trained)
		boost_type::B			# the type of boosting
		combiner::U			# combiner
	end

	# Aliases
	const AdaBoostEnsembleM1Cell{F<:Void, T, B<:AdaBoostM1, U} = AdaBoostEnsemble{F,T,B,U}			# M1 boost, no execution function
	const AdaBoostEnsembleM1Classic{F<:Function, T, B<:AdaBoostM1, U} = AdaBoostEnsemble{F,T,B,U}		# M1 boost, with execution function
	const AdaBoostEnsembleM2Cell{F<:Void, T, B<:AdaBoostM2, U} = AdaBoostEnsemble{F,T,B,U}			# M2 boost, no execution function
	const AdaBoostEnsembleM2Classic{F<:Function, T, B<:AdaBoostM2, U} = AdaBoostEnsemble{F,T,B,U}		# M2 boost, with execution function

	# Printers
	Base.show(io::IO, m::AdaBoostEnsembleM1Cell) = print(io, "AdaBoost M1 ensemble, $(length(m.weights)) members, Cell version")
	Base.show(io::IO, m::AdaBoostEnsembleM2Cell) = print(io, "AdaBoost M2 ensemble, $(length(m.weights)) members, Cell version")
	Base.show(io::IO, m::AdaBoostEnsembleM1Classic) = print(io, "AdaBoost M1 ensemble, $(length(m.weights)) members")
	Base.show(io::IO, m::AdaBoostEnsembleM2Classic) = print(io, "AdaBoost M2 ensemble, $(length(m.weights)) members")

	"""
		adaboost_train(X, y, L, f_train, f_exec, boost_type)

	Trains an AdaBoost classifier. The function returns a `Adaboost` object.

	# Arguments
	 * `X::AbstractMatrix` is the training data; columns are assumed to be samples, lines variables
	 * `y::AvstractVector` are the labels
	 * `L::Int` is the ensemble size
	 * `f_train::Function` is the training function for the ensemble members
	 * `f_exec::Function` is the execution function for the ensemble members
	 * `boost_type::AdaBoost.BoostType` specifies the training algorithm to use. Supported: `AdaBoostM1` and `AdaBoostM2`;

	It is important to note that the training method calls internally `x->f_train((X,y))` where `X` is the input data.
	Considering a generic training function `g` desired to be used in the ensemble, of the following forms:
	 * `g(X, y, arg1, arg2)`, one should provide as argument something as `f_train = (x)->g_train(x[1],y[2], arg1, arg2) # uses labels`
	 * `g(X, arg1, arg2)`, one should provide as argument something as `f_train = (x)->g_train(x[1],arg1, arg2) # ignores labels`
	 * `g(X, arg1, arg2)` where `X==(data,labels)`, one should provide as argument something as `f_train = (x)->g_train(x, arg1, arg2) # MLDataUtils container case`
	 * `g(args...)::CellFunU`, one can simply provide `g.f` (i.e. the training function of the untrained function cell) since it already contains the training arguments.

	During the training, errors are calculated by calling `f_exec(clf,X)` for a given trained ensemble member `clf`.
	"""
	adaboost_train(X::T where T<:AbstractMatrix, y::S where S<:AbstractVector, L::Int, f_train::Function, f_exec::Union{Function,Void}, ::AdaBoostM1) =
		adaboost_train_m1(X, y, L, f_train, f_exec)
	adaboost_train(X::T where T<:AbstractMatrix, y::S where S<:AbstractVector, L::Int, f_train::Function, f_exec::Union{Function,Void}, ::AdaBoostM2) =
		adaboost_train_m2(X, y, L, f_train, f_exec)

	function adaboost_train_m1(X::T where T<:AbstractMatrix, y::S where S<:AbstractVector, L::Int, f_train::Function, f_exec::Union{Function,Void})

		# Initializations
		n = nobs(X)			# number of observations
		yu = sort(unique(y))		# unique labels
		C = length(yu)			# number of classes
		w = fill(1/n, n) 		# sample weights
		beta = fill(0.0,L)		# classifier weight

		k = 1
		it = 0
		maxit = max(L, 10_000)
		mismatches = zeros(n)
		members = Vector{typeof(f_train(getobs(undersample((X,y)))))}(L)
		idxs = Vector{Int}(n)								# sampled observation indexes
		idx = collect(1:n)								# observation indexes
		while k <= L
			it += 1
			sample!(idx, weights(w), idxs)						# get sampling indices

			# Check for missing labels, if any don't go further, resample
			if isempty(setdiff(yu, unique(y[idxs])))

				members[k] = f_train(getobs(datasubset((X,y), idxs)));		# train classifier

				# Split execution between the cases where the execution function is implicit (for FunctionCells) or specified (generic)
				if f_exec isa Void
					mismatches = float(yu[targets(indmax, members[k](X))] .!= y)	# find mismatches (FunctionCell classifier)
				else
					mismatches = float(f_exec(members[k],X) .!= y)		# find mismatches
				end

				err = sum(mismatches .* w)					# calculate error for member 'k'
				if err == 0.0
					fill!(w, 1.0/n)						# re-initialize weights
					k += 1
				elseif err >= (1.0-1/C)
					fill!(w,1.0/n)						# reinitialize weights but do not change iteration
				else
					beta[k] = err/(1-err)					# beta is directly proportional with the error
					w = w .* (beta[k].^(1.0.-mismatches)) 			# if mismatch keep weight, if not, decrease weight
					w = w./sum(w)						# normalize
					k+=1
				end
			end

			if it == maxit && k<L							# Break if maximum number of iterations reached
				warn("[adaboost] Maximum number of $maxit iterations reached. Training interrupted.")
				members = members[1:k]
				break
			end
		end
		return AdaBoostEnsemble(log.(1 ./(beta.+eps())), f_train, f_exec, members, AdaBoostM1(), WeightedVoteCombiner(L, C, log.(1 ./(beta.+eps()))))
	end

	function adaboost_train_m2(X::T where T<:AbstractMatrix, y::S where S<:AbstractVector{V}, L::Int, f_train::Function, f_exec::Union{Function,Void}) where {V}

		# Labels need to be a matrix of Floats encoded OneOfK to easily compute distances
		# between their representation and the probabilities return by the classifiers (ensemble members)
		yu::Vector{V} = sort(unique(y))
		yenc::Vector{Int} = Int.(ohenc_integer(y,yu))

		# Initializations
		n = nobs(X)			# number of observations
		C = length(yu)			# number of classes
		beta = fill(0.0,L)		# classifier weight

		k = 1
		it = 0
		maxit = max(L, 10_000)
		members = Vector{typeof(f_train(getobs(undersample((X,y)))))}(L)

		# Initialize weight vector (mislabel distribution)
		# (The function does not seem to bring any advantage...)
		function _init_weights_!(w::Vector{Float64}, l::Vector{Int})
			ul::Vector{Int} = unique(l)
			counts::Vector{Float64} = length(l) - countapp(l, ul)
			for i in eachindex(l)
				@inbounds w[i] = 1/(length(l)-counts[find(ul.==l[i])][1])
			end
		end
		w = Vector{Float64}(n)
		fill!(w, 1.0/n)
		#_init_weights_!(w, yenc)

		idxs = Vector{Int}(n)								# sampled observation indexes
		idx = collect(1:n)								# observation indexes
		p = zeros(Float64, C, n) 							# ensemble outputs

		while k <= L
			it += 1
			sample!(idx, weights(w), idxs)						# get sampling indices

			# Check for missing labels, if any don't go further, resample
			if isempty(setdiff(yu, unique(yenc[idxs])))

				members[k] = f_train(getobs(datasubset((X,yenc), idxs)));	# train classifier

				# Split execution between the cases where the execution function is implicit (for FunctionCells) or specified (generic)
				if f_exec isa Void
					p = members[k](X)
				else
					p = f_exec(members[k],X)				# 'f_exec' must return a matrix of class-wise probabilities sized C x n
				end

				# Calculate the pseudo-loss
				err=0.0
				@simd for i in 1:n
					@inbounds err+=w[i]*(1-p[yenc[i],i])
				end

				if err == 0.0
					#_init_weights_!(w,yenc)				# re-initialize weights
					fill!(w, 1.0/n)
					k += 1
				elseif err >= (1-1/C)
					#_init_weights_!(w,yenc)				# reinitialize weights but do not change iteration
					fill!(w, 1.0/n)
				else
					beta[k] = err/(1-err)					# beta is directly proportional with the error
					@simd for i in 1:n
						@inbounds w[i] = w[i]*beta[k]^p[yenc[i],i]
					end
					w = w./sum(w)						# normalize
					k+=1
				end
			end

			if it == maxit && k<L							# Break if maximum number of iterations reached
				warn("[adaboost] Maximum number of $maxit iterations reached. Training interrupted.")
				members = members[1:k]
				break
			end
		end
		return AdaBoostEnsemble(log.(1 ./(beta.+eps())), f_train, f_exec, members, AdaBoostM2(), WeightedMeanCombiner(L,C,1.0,log.(1 ./(beta.+eps()))))
	end

	"""
		adaboost_exec(ensemble, X)

	Executes the AdaBoost ensemble.

	# Arguments
	 * `ensemble::AdaBoost.AdaBoostEnsemble` the AdaBoost ensemble
	 * `X::AbstractMatrix` is the data; columns are assumed to be samples, lines variables
	"""
	function adaboost_exec(ensemble::AdaBoostEnsembleM1Classic, X::T where T<:AbstractMatrix)
		out = [ensemble.f_exec(ensemble.members[i], getobs(X)) for i in eachindex(ensemble.members)]
		return combiner_exec(ensemble.combiner, out)

		# If results have to be of matrix form:
		#return float(MLLabelUtils.convertlabel(MLLabelUtils.LabelEnc.OneOfK(Val{ensemble.combiner.C}), combiner_exec(out, ensemble.combiner)))
	end

	function adaboost_exec(ensemble::AdaBoostEnsembleM2Classic, X::T where T<:AbstractMatrix)
		out = combiner_exec(ensemble.combiner,
			 		[ensemble.f_exec(ensemble.members[i], getobs(X)) for i in eachindex(ensemble.members)])
		for j in 1:size(out,2)
			@inbounds softmax!(view(out,:,j))
		end
		return out
	end

	# j4pr case: no execution function, the ensemble members contain it already
	function adaboost_exec(ensemble::AdaBoostEnsembleM1Cell, X::T where T<:AbstractMatrix)
		out = [targets(indmax, ensemble.members[i](getobs(X))) for i in eachindex(ensemble.members)]  # the 'targets' function is used to get the labels
		return combiner_exec(ensemble.combiner, out)

		# If results have to be of matrix form:
		#return float(MLLabelUtils.convertlabel(MLLabelUtils.LabelEnc.OneOfK(Val{ensemble.combiner.C}), combiner_exec(out, ensemble.combiner)))
	end

	function adaboost_exec(ensemble::AdaBoostEnsembleM2Cell, X::T where T<:AbstractMatrix)
		out = combiner_exec(ensemble.combiner,
		     			[ensemble.members[i](getobs(X)) for i in eachindex(ensemble.members)])
		for j in 1:size(out,2)
			@inbounds softmax!(view(out,:,j))
		end
		return out
	end

end

##########################
# FunctionCell Interface #
##########################
"""
	adaboost(f, L; boost_type=AdaBoostM2())

Creates an untrained AdaBoost ensemble.

# Arguments
 * `f::CellFunU` is a untrained cell corresponding to the classifer/transform that the ensemble will contain
 * `L::Int` is the size of the ensemble

# Keyword arguments
 * `boost_type::AdaBoost.BoostType` specifies whether to use the discrete `AdaBoostM1()` or continuous `AdaBoostM2()` approaches in AdaBoost (default `AdaBoostM1()`)

For more information:
	[1] J. Friedman, T. Hastie, R. Tibshirani "Additive Logistic Regression: A statistical view of boosting" 2000 The Annals of Statistics 28(2): 337-407.

Try `?j4pr.AdaBoost.adaboost_train` and `?j4pr.AdaBoost.adaboost_exec` and for more details.

# Examples
```
julia> using j4pr
```
"""
adaboost(f::F where F<:CellFunU, L::Int; boost_type::AdaBoost.BoostType=AdaBoost.AdaBoostM2()) =
	FunctionCell(adaboost, (f, L), ModelProperties(), "Adaboost ensemble ($(f.tinfo), $L members)"; boost_type=boost_type) # untrained function cell

############################
# DataCell/Array Interface #
############################
"""
	adaboost(x, f, L; boost_type=AdaBoostM2())

Trains an AdaBoost ensemble using the data `x`.
"""
# Training
adaboost(x::T where T<:CellDataL, f::F where F<:CellFunU, L::Int; boost_type::AdaBoost.BoostType=AdaBoost.AdaBoostM2()) =
	adaboost((getx!(x), gety(x)), f, L; boost_type=boost_type)

adaboost(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, f::F where F<:CellFunU, L::Int; boost_type::AdaBoost.BoostType=AdaBoost.AdaBoostM2()) =
	adaboost((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), f, L; boost_type=boost_type)

adaboost(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, f::F where F<:CellFunU, L::Int; boost_type::AdaBoost.BoostType=AdaBoost.AdaBoostM2()) =
begin

	@assert nobs(x[1]) == nobs(x[2]) "[adaboost] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	ensembledata = AdaBoost.adaboost_train(x[1], yenc, L, f.f, nothing, boost_type)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(adaboost, Model(ensembledata, modelprops), "Adaboost ensemble ($(f.tinfo), $L members)")
end

# Execution
adaboost(x::T where T<:CellData, model::Model{<:AdaBoost.AdaBoostEnsemble}) =
	datacell(adaboost(getx!(x), model), gety(x))
adaboost(x::T where T<:AbstractVector, model::Model{<:AdaBoost.AdaBoostEnsemble}) =
	adaboost(mat(x, LearnBase.ObsDim.Constant{2}()), model)
adaboost(x::T where T<:AbstractMatrix, model::Model{<:AdaBoost.AdaBoostEnsemble}) =
	AdaBoost.adaboost_exec(model.data,x)


text:
using Test: @testset, @test, @test_throws, @inferred
import ConvolutionalRBMs as ConvRBMs
using ConvolutionalRBMs: reshape_maybe, expand_tuple, samepad, splitpad

@testset "reshape_maybe" begin
    @test (@inferred reshape_maybe(1, ())) == 1
    @test_throws Exception reshape_maybe(1, 1)
    @test_throws Exception reshape_maybe(1, (1,))

    @test (@inferred reshape_maybe(fill(1), ())) == 1
    @test (@inferred reshape_maybe(fill(1), (1,))) == [1]
    @test (@inferred reshape_maybe(fill(1), (1,1))) == hcat([1])
    @test_throws Exception reshape_maybe(fill(1), (1,2))

    @test (@inferred reshape_maybe([1], ())) == 1
    @test (@inferred reshape_maybe([1], (1,))) == [1]
    @test (@inferred reshape_maybe([1], (1,1))) == hcat([1])
    @test_throws Exception reshape_maybe([1], (1,2))

    A = randn(2,2)
    @test (@inferred reshape_maybe(A, 4)) == reshape(A, 4)
end

@testset "expand_tuple" begin
    @test expand_tuple(Val(2), 1) == (1,1)
    @test expand_tuple(Val(2), 2) == (2,2)
    @test expand_tuple(Val(3), 1) == (1,1,1)
    @test expand_tuple(Val(2), (1,2)) == (1,2)
    @test expand_tuple(Val(3), (1,2,3)) == (1,2,3)
    @test_throws MethodError expand_tuple(Val(3), (1,2))
end

@testset "samepad" begin
    @test expand_tuple(Val(2), 1) == (1,1)
    @test expand_tuple(Val(2), 2) == (2,2)
    @test expand_tuple(Val(3), 1) == (1,1,1)
    @test expand_tuple(Val(2), (1,2)) == (1,2)
    @test expand_tuple(Val(3), (1,2,3)) == (1,2,3)
    @test_throws MethodError expand_tuple(Val(3), (1,2))
end

@testset "splitpad" begin
    @test splitpad((1,2,3,4,5,6)) == (lo = (1,3,5), hi = (2,4,6))
end

ext: jl


text:
mutable struct Permute_Instruct
    funcs::AbstractVector{<:Function}
    weights::AbstractVector{<:AbstractFloat}
    args::AbstractVector{<:AbstractVector{<:Tuple{<:Symbol,<:Any}}}
    model_limit::Integer
    func_limit::Integer
    min_clmps::AbstractVector{<:AbstractFloat}
    max_clmps::AbstractVector{<:AbstractFloat}
    override_time::AbstractFloat
    override_weights::AbstractVector{<:AbstractFloat}
    Permute_Instruct(funcs,
                    weights,
                    model_limit,
                    func_limit;
                    min_clmps=fill(.01,length(funcs)),
                    max_clmps=fill(1.,length(funcs)),
                    override_time=0.,
                    override_weights=zeros(length(funcs)),
                    args=[Vector{Tuple{Symbol,Any}}() for i in 1:length(funcs)])=assert_permute_instruct(funcs,weights,args,model_limit,func_limit,min_clmps,max_clmps,override_time,override_weights) && new(funcs,weights,args,model_limit,func_limit,min_clmps,max_clmps,override_time,override_weights)
end

function assert_permute_instruct(funcs,weights,args,model_limit,func_limit,min_clmps, max_clmps, override_time, override_weights)
    !(length(funcs)==length(args)==length(weights)==length(min_clmps)==length(max_clmps)) && throw(ArgumentError("A valid Permute_Instruct must have as many tuning weights and argument vectors as functions!"))
    model_limit<1 && throw(ArgumentError("Permute_Instruct limit on models to permute must be positive Integer!"))
    func_limit<1 && throw(ArgumentError("Permute_Instruct limit on fuction calls per model permtued must be positive Integer!"))
    any(min_clmps.>=max_clmps) && throw(ArgumentError("Permute_Instruct max_clmps must all be greater than corresponding min_clmps!"))
    sum(min_clmps)>1. && throw(ArgumentError("Sum of minimum clamps must be <1.0 to maintain a valid probability vector!"))
    any(max_clmps.>1.) && throw(ArgumentError("Permute_Tuner maximum clamps cannot be >1.0!"))
    override_time<0. && throw(ArgumentError("Permute_Instruct override_time must be 0 (disabled) or positive float!"))
    override_time>0. && !isprobvec(override_weights) && throw(ArgumentError("Permute_Instruct override_weights must be valid probvec!"))
    return true
end

function permute_IPM(e::IPM_Ensemble, instruction::Permute_Instruct)
    call_report=Vector{Tuple{Int64,Float64,Float64}}()
	for model = 1:instruction.model_limit
		m_record = rand(e.models)
        m = deserialize(m_record.path)

        model_blacklist=copy(m.permute_blacklist)
        filteridxs, filtered_funcs, filtered_weights, filtered_args = filter_permutes(instruction, model_blacklist)

        for call in 1:instruction.func_limit
            start=time()
            funcidx=rand(filtered_weights)
            permute_func=filtered_funcs[funcidx]
            pos_args,kw_args=get_permfunc_args(permute_func,e,m,filtered_args[funcidx])
            new_m=permute_func(pos_args...;kw_args...)
            push!(call_report,(filteridxs[funcidx],time()-start,new_m.log_Li - e.contour))

            if length(new_m.permute_blacklist) > 0 && new_m.log_Li == -Inf #in this case the blacklisted function will not work on this model at all; it should be removed from the functions to use
                vcat(model_blacklist,new_m.permute_blacklist)
                filteridxs, filtered_funcs, filtered_weights, filtered_args = filter_permutes(instruction, model_blacklist)
            end

			dupecheck(new_m,m) && new_m.log_Li > e.contour && return new_m, call_report
		end
	end
	return nothing, call_report
end

function permute_IPM(e::IPM_Ensemble, job_chan::RemoteChannel, models_chan::RemoteChannel, comms_chan::RemoteChannel) #ensemble.models is partially updated on the worker to populate arguments for permute funcs
	persist=true
    id=myid()
    put!(comms_chan,id)
    while persist
        wait(job_chan)
        start=time()
        e.models, e.contour, instruction = fetch(job_chan)
        instruction == "stop" && (persist=false; break)

        call_report=Vector{Tuple{Int64,Float64,Float64}}()
        for model=1:instruction.model_limit
			found::Bool=false
            m_record = rand(e.models)

            remotecall_fetch(isfile,1,m_record.path) ? m = (remotecall_fetch(deserialize,1,m_record.path)) : (break)

            model_blacklist=copy(m.permute_blacklist)
            filteridxs, filtered_funcs, filtered_weights, filtered_args = filter_permutes(instruction, model_blacklist)

            for call in 1:instruction.func_limit
                start=time()
                funcidx=rand(filtered_weights)
                permute_func=filtered_funcs[funcidx]
                pos_args,kw_args=get_permfunc_args(permute_func,e,m,filtered_args[funcidx])
                new_m=permute_func(pos_args...;kw_args...)
                push!(call_report,(filteridxs[funcidx],time()-start,new_m.log_Li - e.contour))

                if length(new_m.permute_blacklist) > 0 && new_m.log_Li == -Inf #in this case the blacklisted function will not work on this model at all; it should be removed from the functions to use
                    vcat(model_blacklist,new_m.permute_blacklist)
                    filteridxs, filtered_funcs, filtered_weights, filtered_args = filter_permutes(instruction, model_blacklist)
                end

                dupecheck(new_m,m) && new_m.log_Li > e.contour && ((put!(models_chan, (new_m ,id, call_report))); found=true; break)

			end
            found==true && break;
            wait(job_chan)
            fetch(job_chan)!=e.models && (break) #if the ensemble has changed during the search, update it
			model==instruction.model_limit && (put!(models_chan, ("quit", id, call_report));persist=false)#worker to put "quit" on channel if it fails to find a model more likely than contour
		end
	end
end
                function filter_permutes(instruction, model_blacklist)
                    filteridxs=findall(p->!in(p,model_blacklist),instruction.funcs)
                    filtered_funcs=instruction.funcs[filteridxs]
                    filtered_weights=filter_weights(instruction.weights, filteridxs)
                    filtered_args=instruction.args[filteridxs]
                    return filteridxs, filtered_funcs, filtered_weights, filtered_args
                end

                function filter_weights(weights, idxs)
                    deplete=0.
                    for (i, weight) in enumerate(weights)
                        !in(i, idxs) && (deplete+=weight)
                    end
                    weights[idxs].+=deplete/length(idxs)
                    new_weights=weights[idxs]./sum(weights[idxs])
                    return Categorical(new_weights)
                end

                function get_permfunc_args(func::Function,e::IPM_Ensemble, m::ICA_PWM_Model, args::Vector{Tuple{Symbol,Any}})
                    pos_args=[]
                    argparts=Base.arg_decl_parts(methods(func).ms[1])
                    argnames=[Symbol(argparts[2][n][1]) for n in 2:length(argparts[2])]
                    for argname in argnames #assemble basic positional arguments from ensemble and model fields
                        if argname == Symbol('m')
                            push!(pos_args,m)
                        elseif argname in fieldnames(IPM_Ensemble)
                            push!(pos_args,getfield(e,argname))
                        elseif argname in fieldnames(ICA_PWM_Model)
                            push!(pos_args,getfield(m,argname))
                        else
                            throw(ArgumentError("Positional argument $argname of $func not available in the ensemble or model!"))
                        end
                    end

                    kw_args = NamedTuple()
                    if length(args) > 0 #if there are any keyword arguments to pass
                        sym=Vector{Symbol}()
                        val=Vector{Any}()
                        for arg in args
                            push!(sym, arg[1]); push!(val, arg[2])
                        end
                        kw_args = (;zip(sym,val)...)
                    end

                    return pos_args, kw_args
                end

                function dupecheck(new_model, model)
                    (new_model.sources==model.sources && new_model.mix_matrix==model.mix_matrix) ? (return false) : (return true)
                end


text:
@testset "one hot encoding" begin

    categories = ["foo", "bar", "baz"]
    ohe = OneHotEncoding(categories)
    @test ohe isa Transform
    @test cardinality(ohe) == OneToMany()

    @testset "Basic" begin

        @testset "specify return type" begin
            x = ["foo", "bar", "baz"]
            expected = [1 0 0; 0 1 0; 0 0 1]
            transformed = FeatureTransforms.apply(x, OneHotEncoding{AbstractFloat}(categories))
            @test transformed == expected
            @test transformed isa AbstractMatrix{AbstractFloat}
        end

        @testset "duplicate values" begin
            x = ["foo", "baz", "bar", "baz"]
            expected = [1 0 0; 0 0 1; 0 1 0; 0 0 1]
            @test FeatureTransforms.apply(x, ohe) == expected

            # Cannot pass duplicate values as the categories for the OHE constructor
            @test_throws ArgumentError OneHotEncoding(x)
        end

        @testset "value does not exist as a category" begin
            x = ["foo", "baz", "bar", "dne"]
            @test_throws KeyError FeatureTransforms.apply(x, ohe)
        end
    end
end


text:
module Totalsizeof

importall Base
export totalsizeof

# Helper: Pointer cache is used to break circular references in objects and avoid double countings
is_seen!(x, ptr_cache) = in(pointer_from_objref(x), ptr_cache) || (push!(ptr_cache, pointer_from_objref(x)); false)

# Helper: Catch types without size method
sizeof_catch(x) = try sizeof(x) catch 0 end

# Composite Types
function totalsizeof(x, ptr_cache = Set())
  is_seen!(x, ptr_cache) && return 0
  result = sizeof_catch(x)

  this_names = nothing

  # Check if type has names method
  try
    this_names = names(x)
  catch
    this_names = []
  end

  el = nothing
  isempty(this_names) ? result :
    result + mapreduce(i ->
      # Check if some elements are not reachable
      (try el = getfield(x,i) catch el = nothing end; totalsizeof(el, ptr_cache))
      , +, this_names)
end

# Tuples
function totalsizeof(x::Tuple, ptr_cache = Set())
  is_seen!(x, ptr_cache) && return 0
  isempty(x) ? 0 : mapreduce(i -> totalsizeof(i, ptr_cache), +, x)
end

# Array
function  totalsizeof{T<:Any}(x::Array{T}, ptr_cache = Set())
  is_seen!(x, ptr_cache) && return 0
  result = sizeof(x)
  el = nothing
  for i in 1:length(x)
    # Check if some elements are undefined
    try el = x[i] catch el = nothing end
    result += totalsizeof(el, ptr_cache)
  end
  result
end

# String, Function, IntrinsicFunction
totalsizeof(x::Union(String, Function, IntrinsicFunction), ptr_cache = Set()) = is_seen!(x, ptr_cache) ? 0 : sizeof(x)

end # Module end


text:
include("main_functions.jl")

#example
instance_name = "case1888rte"
matpower_instance_path = "D:\\repo\\data\\data_Matpower\\matpower\\case14.m"
output_instance_path = "D:\\repo\\data\\data_ROPF\\ROPFu"
output_decomposition_path = "D:\\repo\\data\\data_sdp"
max_time = 300 #1 hour
ROPF = ROPF_infos(instance_name,
matpower_instance_path,
output_instance_path,
"cholesky",
output_decomposition_path)

# typeofconstraint = NoConstraint
typeofconstraint = MAXkshuntsConstraint
# typeofconstraint = MAXkmovesConstraint
nb_max_moves_or_shunts = 4

# UB, LB = optvalue_bounds(ROPF, nb_max_moves_or_shunts, typeofconstraint)

BB_parameters = BB_infos("deepfirst", "1", 0.75, 0)

UB, LB = solve_BandB(ROPF, max_time, BB_parameters, nb_max_moves_or_shunts, typeofconstraint)


text:
export IHTVariable, IHTVariables, use_A2_as_minor_allele, standardize_snpmatrix
export use_A2_as_minor_allele, project_k!, project_group_sparse!, _iht_gradstep, _init_iht_indices
export _iht_omega, _iht_backtrack, std_reciprocal, gIHTResults

"""
Object to contain intermediate variables and temporary arrays. Used for cleaner code in L0_reg
"""
mutable struct IHTVariable{T <: Float, V <: DenseVector}

    #TODO: Consider changing b and b0 to SparseVector
    #TODO: Add in itc
    # itc  :: T             # the intercept
    # itc0 :: T             # old intercept
    b     :: Vector{T}     # the statistical model, most will be 0
    b0    :: Vector{T}     # previous estimated model in the mm step
    xb    :: Vector{T}     # vector that holds x*b
    xb0   :: Vector{T}     # previous xb in the mm step
    xk    :: SnpLike{2}    # the n by k subset of the design matrix x corresponding to non-0 elements of b
    gk    :: Vector{T}     # gk = df[idx] is a temporary array of length `k` that arises as part of the gradient calculations. I avoid doing full gradient calculations since most of `b` is zero.
    xgk   :: Vector{T}     # x * gk also part of the gradient calculation
    idx   :: BitVector     # BitArray indices of nonzeroes in b for A_mul_B
    idx0  :: BitVector     # previous iterate of idx
    r     :: V             # n-vector of residuals
    df    :: V             # the gradient: df = -x' * (y - xb - intercept)
    group :: Vector{Int64} # vector denoting group membership
end

function IHTVariables{T <: Float}(
    x :: SnpLike{2},
    y :: Vector{T},
    J :: Int64,
    k :: Int64
)
    n, p = size(x)
    p += 1 # add 1 for the intercept, need to change this to use itc later

    # itc  = zero(T)
    # itc0 = zero(T)
    b     = zeros(T, p)
    b0    = zeros(T, p)
    xb    = zeros(T, n)
    xb0   = zeros(T, n)
    xk    = SnpArray(n, J*k)
    gk    = zeros(T, J*k)
    xgk   = zeros(T, n)
    idx   = falses(p)
    idx0  = falses(p)
    r     = zeros(T, n)
    df    = zeros(T, p)
    group = ones(Int64, p)

    return IHTVariable{T, typeof(y)}(b, b0, xb, xb0, xk, gk, xgk, idx, idx0, r, df, group)
    # return IHTVariable{T, typeof(y)}(itc, itc0, b, b0, xb, xb0, xk, gk, xgk, idx, idx0, r, df)
end

#"""
#Returns ω, a constant we need to bound the step size μ to guarantee convergence.
#"""
#
#function compute_ω!(v::IHTVariable, snpmatrix::Matrix{Float64})
#    #update v.xb
#    A_mul_B!(v.xb, snpmatrix, v.b)
#
#    #calculate ω efficiently (old b0 and xb0 have been copied before calling iht!)
#    return sqeuclidean(v.b, v.b0) / sqeuclidean(v.xb, v.xb0)
#end

"""
This function is needed for testing purposes only.
Converts a SnpArray to a matrix of float64 using A2 as the minor allele. We want this function
because SnpArrays.jl uses the less frequent allele in each SNP as the minor allele, while PLINK.jl
always uses A2 as the minor allele, and it's nice if we could cross-compare the results.
"""
function use_A2_as_minor_allele(snpmatrix :: SnpArray)
    n, p = size(snpmatrix)
    matrix = zeros(n, p)
    for i in 1:p
        for j in 1:n
            if snpmatrix[j, i] == (0, 0); matrix[j, i] = 0.0; end
            if snpmatrix[j, i] == (0, 1); matrix[j, i] = 1.0; end
            if snpmatrix[j, i] == (1, 1); matrix[j, i] = 2.0; end
            if snpmatrix[j, i] == (1, 0); matrix[j, i] = missing; end
        end
    end
    return matrix
end

"""
This function computes the gradient step v.b = P_k(β + μ∇f(β)), and updates v.idx.
Recall calling axpy! implies v.b = v.b + μ*v.df, but v.df stores an extra negative sign.
"""
function _iht_gradstep{T <: Float}(
    v  :: IHTVariable{T},
    μ  :: T,
    J  :: Int,
    k  :: Int;
)
   # n = length(v.xb)
   # v.itc = v.itc0 + μ*(n*μ - sum(v.r) + n*v.itc0)
   BLAS.axpy!(μ, v.df, v.b) # take the gradient step: v.b = b + μ∇f(b) (which is an addition since df stores X(-1*(Y-Xb)))
   project_group_sparse!(v.b, v.group, J, k) # project to doubly sparse vector
   v.idx .= v.b .!= 0       # Update idx. (find indices of new beta that are nonzero)

   # If the k'th largest component is not unique, warn the user.
   sum(v.idx) <= J*k || warn("More than J*k components of b is non-zero! Need: VERY DANGEROUS DARK SIDE HACK!")
end

"""
When initializing the IHT algorithm, take largest elements of each group of df as nonzero
components of b. This function set v.idx = 1 for those indices.
"""
function _init_iht_indices{T <: Float}(
    v     :: IHTVariable{T},
    J     :: Int,
    k     :: Int
)
    perm = sortperm(v.df, by = abs, rev = true)
    project_group_sparse!(v.df, v.group, J, k)
    v.idx[find(v.df)] = true
    v.gk = zeros(T, sum(v.idx))

    return nothing
end

"""
this function calculates the omega (here a / b) used for determining backtracking
"""
function _iht_omega{T <: Float}(
    v :: IHTVariable{T}
)
    a = sqeuclidean(v.b, v.b0::Vector{T}) :: T
    b = sqeuclidean(v.xb, v.xb0::Vector{T}) :: T
    return a, b
end

"""
this function for determining whether or not to backtrack. True = backtrack
"""
function _iht_backtrack{T <: Float}(
    v       :: IHTVariable{T},
    ot      :: T,
    ob      :: T,
    mu      :: T,
    mu_step :: Int,
    nstep   :: Int
)
    mu*ob > 0.99*ot              &&
    sum(v.idx) != 0              &&
    sum(xor.(v.idx,v.idx0)) != 0 &&
    mu_step < nstep
end

"""
Compute the standard deviation of a SnpArray in place
"""
function std_reciprocal{T <: Float}(A::SnpArray, mean_vec::Vector{T})
    m, n = size(A)
    @assert n == length(mean_vec) "number of columns of snpmatrix doesn't agree with length of mean vector"
    std_vector = zeros(T, n)

    @inbounds for j in 1:n
        @simd for i in 1:m
            (a1, a2) = A[i, j]
            if !isnan(a1, a2) #only add if data not missing
                std_vector[j] += (convert(T, a1 + a2) - mean_vec[j])^2
            end
        end
        std_vector[j] = 1.0 / sqrt(std_vector[j] / (m - 1))
    end
    return std_vector
end

""" Projects the point y onto the set with at most m active groups and at most
n active predictors per group. The variable group encodes group membership. Currently
assumes there are no unknown group membership.
TODO: 1)make this function operate in-place
      2)check if sortperm can be replaced by something that doesn't sort the whole array
"""
function project_group_sparse!{T <: Float}(
    y     :: Vector{T},
    group :: Vector{Int64},
    m     :: Int64,
    n     :: Int64
)
    groups = maximum(group)
    group_count = zeros(Int, groups)         #counts number of predictors in each group
    z = zeros(groups)                        #l2 norm of each group
    perm = zeros(Int64, length(y))           #vector holding the permuation vector after sorting
    sortperm!(perm, y, by = abs, rev = true)

    #calculate the magnitude of each group, where only top predictors contribute
    for i in eachindex(y)
        j = perm[i]
        k = group[j]
        if group_count[k] < n
            z[k] = z[k] + y[j]^2
            group_count[k] = group_count[k] + 1
        end
    end

    #go through the top predictors in order. Set predictor to 0 if criteria not met
    group_rank = zeros(Int64, length(z))
    sortperm!(group_rank, z, rev = true)
    group_rank = invperm(group_rank)
    fill!(group_count, 1)
    for i in eachindex(y)
        j = perm[i]
        k = group[j]
        if (group_rank[k] > m) || (group_count[k] > n)
            y[j] = 0.0
        else
            group_count[k] = group_count[k] + 1
        end
    end
end

"""
an object that houses results returned from a group IHT run
"""
immutable gIHTResults{T <: Float, V <: DenseVector}
    time  :: T
    loss  :: T
    iter  :: Int
    beta  :: V
    J     :: Int64
    k     :: Int64
    group :: Vector{Int64}

    #gIHTResults{T,V}(time::T, loss::T, iter::Int, beta::V) where {T <: Float, V <: DenseVector{T}} = new{T,V}(time, loss, iter, beta)
    gIHTResults{T,V}(time, loss, iter, beta, J, k, group) where {T <: Float, V <: DenseVector{T}} = new{T,V}(time, loss, iter, beta, J, k, group)
end

# strongly typed external constructor for gIHTResults
gIHTResults(time::T, loss::T, iter::Int, beta::V, J::Int, k::Int, group::Vector{Int}) where {T <: Float, V <: DenseVector{T}} = gIHTResults{T, V}(time, loss, iter, beta, J, k, group)

"""
a function to display gIHTResults object
"""
function Base.show(io::IO, x::gIHTResults)
    println(io, "IHT results:")
    println(io, "\nCompute time (sec):     ", x.time)
    println(io, "Final loss:             ", x.loss)
    println(io, "Iterations:             ", x.iter)
    println(io, "Max number of groups:   ", x.J)
    println(io, "Max predictors/group:   ", x.k)
    println(io, "IHT estimated ", countnz(x.beta), " nonzero coefficients.")
    non_zero = find(x.beta)
    print(io, DataFrame(Group=x.group[non_zero], Predictor=non_zero, Estimated_β=x.beta[non_zero]))
    return nothing
end


text:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using State

include("test-unitaire-LSAtTmod.jl")
include("test-unitaire-NLPAtXmod.jl")


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf, LaTeXStrings
# if !(pwd() in LOAD_PATH) push!(LOAD_PATH, pwd()) end
# using TaylorBigF

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 3.3 Extrapolation

p=500;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("1"); n=BigInt(201); iN=convert(Int64,n); dx=BigFloat("2")*L/(n-1); x=-L:dx:L;
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end
V=TaylorBigF.calc_Vandermonde(x,n,n)
xi=x[end]+dx:dx:100x[end]
Vi=TaylorBigF.calc_Vandermonde(xi,length(xi),n)

f=sin.(x); a=V\f;
fi_comp=Vi*a
fi_analyt=sin.(xi)
error1=convert(Array{Float64},fi_comp-fi_analyt)
iok=abs.(error1).<1

# plot(x,f,color=:black,width=5,xlabel= "x" ,ylabel="f(x)",
#     label="given f")
# plot!(xi[iok],fi_analyt[iok],color=:black,width=2,label="exact f")
# plot!(xi[iok],fi_comp[iok],color=:black,width=2,
#         label=string("extrapolated f, p=",p),linestyle=:dash)
# plot!(xi[iok],fi_comp[iok],color=:black,width=2,
#         label=string("extrapolated f, p=",p),linestyle=:dashdot)
# plot!(xi[iok],fi_comp[iok],color=:black,width=2,
#         label=string("extrapolated f, p=",p),linestyle=:dashdotdot)

# savefig(string("extrap.pdf"))

r_root=Array{Float64}(undef,0); for i=1:Int64(n) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
1/r_root[end]
xi[iok][end]


text:
using BSON: @save, @load
using Plots, ArgParse, DelimitedFiles

function parseCommandLine(; dir=0)

        # initialize the settings (the description is for the help screen)
        # s = ArgParseSettings(description = "Example 1 for argparse.jl: minimal usage.")
		s = ArgParseSettings(description = "Dynamical Structure Factor code.")

        @add_arg_table! s begin
           "--Filename", "-f"         # another option, with short form
                help = "File name"
                arg_type = String
                default = "Ge00Sb00Te100T823K.cpmd"
			"--pathIn", "-i"         # another option, with short form
                 help = "Directory where input file is stored."
                 arg_type = String
                 default = pwd()
			"--pathOut", "-o"         # another option, with short form
                 help = "Directory where everything will be saved. This directory can be created before or during execution. Default is a new directory named after the filename, located where the filename is."
                 arg_type = String
                 default = "temp"
           "--n", "-N"
                    help = "number of particles"
					arg_type = Int
                    default = 300
			"--mass", "-m"
                     help = "mass"
 					arg_type = Float32
                     default = 1.0f0
           "--Lchar", "-L"
                help = "box length"
                arg_type = Float32
                default = 42.2976f0
           "--t_max", "-t"
                help = "t maximum. NO LONGER NEEDED"
                arg_type = Int64
                default = 0
			"--dt", "-s"
                 help = "time step between configurations"
                 arg_type = Float32
                 default = 10.0f0
			 "--dim", "-d"
	 			 help = "dimensions of time step between configurations"
	 			 arg_type = String
	 			 default = "ps"
			"--sparse", "-b"
				help = "Use sparse matrices? Should be faster..."
				arg_type = Int
				default = 0
			"--time_window", "-w"
				help = "List of time window of samples. Should be space separated"
				nargs = '*'
				arg_type = Int64
				default = [100, 500, 1000]
			"--knum", "-k"
				help = "List with kmin and kmax. Should be space separated"
				nargs = '*'
				arg_type = Int64
				default = [1, 10]
			"--test"
				help = "Test code"
				arg_type = Bool
				default = false
	   end

        return parse_args(s) # the result is a Dict{String,Any}
end

function init_dict(kwargs)
	if kwargs["test"] == false
		filename = kwargs["Filename"]
		path_to_file = kwargs["pathIn"] * "/" * kwargs["Filename"]
	    pathout = kwargs["pathOut"] == "temp" ? kwargs["pathIn"] * "/" * split(kwargs["Filename"],".")[1] * "/" : kwargs["pathOut"]

		dict = Dict()
	    dict[:PathIn] = path_to_file
	    dict[:PathOut] = pathout
		dict[:t_max_parsed] = kwargs["t_max"]
		dict[:partitions] = Dict([(:t_samp, []), (:samples, []), (:t_max, [])])
	    dict[:L_char] = kwargs["Lchar"]
		dict[:num_part] = kwargs["n"]
	    dict[:title] = filename
		dict[:dt] = kwargs["dt"]
		dict[:dim] = kwargs["dim"]
		dict[:t_samp_list] = kwargs["time_window"]
		dict[:knum] = kwargs["knum"]
		dict[:mass] = kwargs["mass"]
	else
		kwargs["Filename"] = "Test"
		filename = kwargs["Filename"]
		path_to_file = kwargs["pathIn"] * "/" * kwargs["Filename"]
	    pathout = kwargs["pathIn"] * "/" * kwargs["Filename"] * "/"
		dict = Dict()
	    dict[:PathIn] = path_to_file
	    dict[:PathOut] = pathout
		dict[:t_max_parsed] = kwargs["t_max"]
		dict[:partitions] = Dict([(:t_samp, []), (:samples, []), (:t_max, [])])
	    dict[:L_char] = 10 #kwargs["Lchar"]
		dict[:num_part] = 1000 #kwargs["n"]
	    dict[:title] = filename
		dict[:dt] = 0.01 #kwargs["dt"]
		dict[:dim] = kwargs["dim"]
		dict[:t_samp_list] = Array{Int64}([1] ./ dict[:dt]) #kwargs["time_window"]
		dict[:knum] = kwargs["knum"]
		dict[:mass] = kwargs["mass"]
		dict[:t_max] = 10
	end

    isdir(pathout) || mkdir(pathout)

    @save (dict[:PathOut] * "Dict.bson") dict
    return dict
end

# function init_dict(path_to_file, pathout, t_MAX, L_char, num_part, title, Δt, dim, t_samp_list, kmax, kmin, m)
#     isdir(pathout) || mkdir(pathout)
#     dict = Dict()
#     dict[:PathIn] = path_to_file
#     dict[:PathOut] = pathout
# 	dict[:t_max] = t_MAX
# 	dict[:partitions] = Dict([(:t_samp, []), (:samples, []), (:t_max, [])])
#     dict[:L_char] = L_char
# 	dict[:num_part] = num_part
#     dict[:title] = title
# 	dict[:dt] = Δt
# 	dict[:dim] = dim
# 	dict[:t_samp_list] = t_samp_list
# 	dict[:knum] = [kmin, kmax]
# 	dict[:mass] = m
#
#     @save (dict[:PathOut] * "Dict.bson") dict
#     return dict
# end

function load_dict(path)
	@load (path * "Dict.bson") dict
	return dict
end

function load_data(dict)
	t_max = dict[:t_max_parsed]
	N = dict[:num_part]
	path = dict[:PathIn]
    data1 = readdlm(path)
    x = reshape(data1[:,2],N,:)'
    y = reshape(data1[:,3],N,:)'
    z = reshape(data1[:,4],N,:)'
    vx = reshape(data1[:,5],N,:)'
    vy = reshape(data1[:,6],N,:)'
    vz = reshape(data1[:,7],N,:)'
    # samp = Int(floor(t_max/t_sample))
    # d = [(x[r,:],y[r,:],z[r,:],vx[r,:],vy[r,:],vz[r,:]) for r in Iterators.partition(1:t_max,t_sample)][1:samp]

    # t_max = t_max == 0 ? end : t_max
    d = t_max != 0 ? (x[1:t_max,:],y[1:t_max,:],z[1:t_max,:],vx[1:t_max,:],vy[1:t_max,:],vz[1:t_max,:]) : (x,y,z,vx,vy,vz)
	dict[:t_max] = size(x,1)
	@info "Data loaded. $(dict[:t_max]) configurations detected."
	d
end

function saveData(ω_list, freq_list, CT, CT_arr, dict; flag="T")
	fig_name = split(dict[:title],".")[1]
	writedlm(dict[:PathOut] * "/" * fig_name * "_modes_$(flag).dat", ω_list')
	writedlm(dict[:PathOut] * "/" * fig_name * "_freq_$(flag).dat", freq_list)

	isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed")
	for k in 1:size(CT,1)
	    writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed/C$(flag)_k_$k.dat", CT[k,:,:])
	end

	isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData")
	for k in 1:size(CT_arr,2), s in 1:size(CT_arr,4)
		tw = dict[:t_samp_list][s]
	    tw2 = Int(floor(tw/2))
	    writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationData/C$(flag)_k_$(k)_TW_$tw.dat", [CT_arr[1:tw2,k,2,s] CT_arr[1:tw2,k,1,s]])
	end
end

function saveData(CT_arr, dict; flag="T")
	fig_name = split(dict[:title],".")[1]
	# isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed")
	# for k in 1:size(CT,1)
	#     writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed/C$(flag)_k_$k.dat", CT[k,:,:])
	# end

	isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData")
	for k in 1:size(CT_arr,2), s in 1:size(CT_arr,4)
		tw = dict[:t_samp_list][s]
	    tw2 = Int(floor(tw/2))
	    writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationData/C$(flag)_k_$(k)_TW_$tw.dat", [CT_arr[1:tw2,k,2,s] CT_arr[1:tw2,k,1,s]])
	end
end

function savePlots(ω_list, freq_list, CT, CT_arr, dict; flag="T", mixedWindows=true)
	dim = dict[:dim]
	#Mean frequency over different time windows.
	m = reshape(mean(freq_list, dims=2),:)
	st = reshape(std(freq_list, dims=2),:)
	fig = plot(1:10, m, ribbon=st, fillalpha=0.2, label="mean", lw=4,
	    	frame=:box, xlabel="k (wave number)", ylabel="freq (1/$dim)", legend=:topleft,
	    	title=dict[:title], margin = 10Plots.mm, ms=5, markershapes = :circle, markerstrokewidth=0)

	# Correlation in mixed windows
	figs_CT = []
    for k in 1:size(CT,1)
		wmaxBound = minimum([200, size(CT,2)])
        f_CT = plot(CT[k,1:wmaxBound,1],CT[k,1:wmaxBound,2], frame=:box, xlabel="freq (1/$dim)", ylabel="C$(flag)", label="k = $k", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
        push!(figs_CT, f_CT)
    end

	# Max freq. Max modes. Correlation per time window.
    figs_freq = []
    figs_ω = []
    figs_CT_arr = []
    for s in 1:size(ω_list,1)
		tw = dict[:t_samp_list][s]
        f_ω = plot(ω_list[s,:], frame=:box, xlabel="k", ylabel="w max", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
        push!(figs_ω, f_ω)

        f_freq = plot(freq_list[:,s], frame=:box, xlabel="k", ylabel="freq max (1/$dim)", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
        push!(figs_freq, f_freq)

        for k in 1:size(CT_arr,2)
			tw2 = Int(floor(tw/2))
            f_CT_arr = plot(CT_arr[1:tw2,k,2,s], CT_arr[1:tw2,k,1,s], frame=:box, xlabel="freq (1/$dim)", ylabel="C$(flag)", label="k = $k, Time Window $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
            push!(figs_CT_arr, f_CT_arr)
        end
    end

	fig_name = split(dict[:title],".")[1]
	savefig(fig, dict[:PathOut] * "/" * fig_name * "_freqMean.png")
	savefig(plot(figs_CT..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_C$(flag).png")
	savefig(plot(figs_freq..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_freq_$(flag).png")
	savefig(plot(figs_ω..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_modes_$(flag).png")
	savefig(plot(figs_CT_arr..., size=(2500,1900)), dict[:PathOut] * "/" * fig_name * "_C$(flag)_arr.png")
end

function savePlots(CT_arr, dict; flag="T_im", mixedWindows=true)
	dim = dict[:dim]

	# # Correlation in mixed windows
	# figs_CT = []
    # for k in 1:size(CT,1)
    #     f_CT = plot(CT[k,1:200,1],CT[k,1:200,2], frame=:box, xlabel="time ($dim)", ylabel="C$(flag)", label="k = $k", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
    #     push!(figs_CT, f_CT)
    # end

	# Correlation per time window.
    figs_CT_arr = []
    for s in 1:size(CT_arr,4)
		tw = dict[:t_samp_list][s]
        # f_ω = plot(ω_list[s,:], frame=:box, xlabel="k", ylabel="w max", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
        # push!(figs_ω, f_ω)

        # f_freq = plot(freq_list[:,s], frame=:box, xlabel="k", ylabel="freq max (1/$dim)", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
        # push!(figs_freq, f_freq)

        for k in 1:size(CT_arr,2)
			tw2 = Int(floor(tw/2))
            f_CT_arr = plot(CT_arr[1:tw2,k,2,s], CT_arr[1:tw2,k,1,s], frame=:box, xlabel="time ($dim)", ylabel="C$(flag)", label="k = $k, Time Window $tw", margin = 5Plots.mm, ms=4, markershapes = :circle, markerstrokewidth=0)
            push!(figs_CT_arr, f_CT_arr)
        end
    end

	fig_name = split(dict[:title],".")[1]
	# savefig(plot(figs_CT..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_C$(flag).png")

	savefig(plot(figs_CT_arr..., size=(2500,1900)), dict[:PathOut] * "/" * fig_name * "_C$(flag)_arr.png")
end

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text:
export FourRooms, FourRoomsCont, MarkovChain, StarProblem

export step!, start!, get_actions, is_terminal, get_reward, get_state
import JuliaRL: AbstractEnvironment, step!, start!, get_actions, is_terminal, get_reward, get_state

include("env/FourRooms.jl")
include("env/FourRoomsCont.jl")
include("env/MarkovChain.jl")


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using CSV
using PyPlot

using HTTP
using StableRNGs

req = HTTP.get("https://raw.githubusercontent.com/rupakc/UCI-Data-Analysis/master/Airfoil%20Dataset/airfoil_self_noise.dat");

df = CSV.read(req.body; header=[
                  "Frequency","Attack_Angle","Chord+Length",
                  "Free_Velocity","Suction_Side","Scaled_Sound"
              ]
       );
df[1:5, :] |> pretty

schema(df)

y, X = unpack(df, ==(:Scaled_Sound), col -> true);

X = transform(fit!(machine(Standardizer(), X)), X);

train, test = partition(eachindex(y), 0.7, shuffle=true, rng=StableRNG(612));

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

coerce!(X, :Frequency=>Continuous)

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

dcr = @load DecisionTreeRegressor pkg=DecisionTree

dcrm = machine(dcr, X, y)

fit!(dcrm, rows=train)
pred_dcrm = MLJ.predict(dcrm, rows=test);

rms(pred_dcrm, y[test])

rfr = @load RandomForestRegressor pkg=DecisionTree

rfr_m = machine(rfr, X, y);

fit!(rfr_m, rows=train);

pred_rfr = MLJ.predict(rfr_m, rows=test);
rms(pred_rfr, y[test])

r_maxD = range(rfr, :n_trees, lower=9, upper=15)
r_samF = range(rfr, :sampling_fraction, lower=0.6, upper=0.8)
r = [r_maxD, r_samF];

tuning = Grid(resolution=7)
resampling = CV(nfolds=6)

tm = TunedModel(model=rfr, tuning=tuning,
                resampling=resampling, ranges=r, measure=rms)

rfr_tm = machine(tm, X, y);

fit!(rfr_tm, rows=train);

pred_rfr_tm = MLJ.predict(rfr_tm, rows=test);
rms(pred_rfr_tm, y[test])

fitted_params(rfr_tm).best_model

r = report(rfr_tm)
res = r.plotting

md = res.parameter_values[:,1]
mcw = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, res.measurements)

xlabel("Number of trees", fontsize=14)
ylabel("Sampling fraction", fontsize=14)
xticks(9:1:15, fontsize=12)
yticks(fontsize=12)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
function test_randomgrids()
    println("Random grids:")
    @testset begin
        d = UnitDisk()
        g = randomgrid(d, 10)
        @test length(g) == 10
        @test length(eltype(g)) == dimension(d)
        @test reduce(&, [x ∈ d for x in g])

        g2 = randomgrid(UnitDisk{BigFloat}(), 10)
        @test eltype(g2[1]) == BigFloat

        g3 = randomgrid(0.0..1.0, 10)
        @test length(g3) == 10
        # 1D is a special case where we don't use vectors of length 1
        @test eltype(g3) == Float64

        box1 = UnitInterval()^2
        p1 = randompoint(box1)
        @test length(p1) == dimension(box1)
        @test p1 ∈ box1
        box2 = 0.0..1.0
        p2 = randompoint(box2)
        @test typeof(p2) == Float64
        @test p2 ∈ box2

        g3 = randompoint(UnionDomain(0.0..1.5,1.5..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 0.0..2.0

        g3 = randompoint(IntersectDomain(0.0..1.5,1.0..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 1.0..1.5

        g3 = randompoint(SetdiffDomain(0.0..1.5,1.0..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 0.0..1.0
    end
end

test_randomgrids()


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 02c15be0-defd-11ea-004b-45c4e09aa541
using Plots

# ╔═╡ 8aeae030-811f-11eb-0075-272d383103ef
using PlutoUI

# ╔═╡ ba6bf300-8122-11eb-0a3e-2d30a150250c
using SymEngine

# ╔═╡ 6cbb32f0-d665-11ea-3f24-41cbe1910b80
using LinearAlgebra

# ╔═╡ 28fccf20-d63c-11ea-3311-411f89a63431
md"# Curvas de nivel"

# ╔═╡ 309959f0-defd-11ea-34c3-71eeda7b523b
md"Curvas cerradas:

$$V(x)\to \infty$$ cuando $$x\to \infty$$"

# ╔═╡ cd88fade-defd-11ea-2097-25b2d1eac285
V1(x1,x2)=x1^2+x2^2

# ╔═╡ ce5fe800-deff-11ea-0b19-7fc23b14bd06
x=y=-1:0.1:1

# ╔═╡ c6c49460-deff-11ea-1c6d-673d863ea07b
begin
	plotly() #Back-end interactivo de Plots
	surface(x,y,V1)
end

# ╔═╡ 5b26b050-d7eb-11ea-2488-b7da296fed24
@bind L1 Slider(0:0.05:2)

# ╔═╡ 2e580002-defe-11ea-1b3d-d5ec461e0715
begin
	gr() #Back-End por defecto
	contour(x, y, V1, fill = true)
	contour!(x,y,V1,levels=[L1],linewidth=3,color=:red, ratio=1)
end

# ╔═╡ 2ec228f0-defd-11ea-329f-53576ab14004
md"Curvas abiertas

$$V(x) \nrightarrow \infty$$ cuando $$x\to \infty$$"

# ╔═╡ ca2fea72-3655-11eb-0984-a532e1f9efcc
md"""Por ejemplo

$V_2(x_1,x_2)=\frac{x_1^2}{1+x_1^2}+x_2^2$

Ya que si $x_1\to \pm\infty$, $V_2 \to 1+x_2^2 \ne   \infty$"""

# ╔═╡ 204e6710-defe-11ea-1fb0-43c7ef863786
V2(x1,x2)=x1^2/(1+x1^2)+x2^2;

# ╔═╡ f3bcc450-df00-11ea-0080-03b54212b9ed
begin
	x2=-5:0.1:5
	y2=-1:0.1:1
	plotly()
	surface(x2,y2,V2,ratio=1)
end

# ╔═╡ ca99149e-df02-11ea-39a0-179e6a7d1344
@bind L2 Slider(0:0.05:2)

# ╔═╡ eb00be70-df00-11ea-0ce9-3b0a79a0f4a3
begin
	gr()
	contour(x2, y2, V2, fill = true)
	contour!(x2,y2,V2,levels=[L2],linewidth=3,color=:red) #, ratio=1)
end

# ╔═╡ d8678810-defc-11ea-05b2-51373f754e61
md"# Semi-Definida Vs Definida
Pensemos **antes** de pintarlo, para ver si hay *sorpresas*"

# ╔═╡ 34a9ca60-df03-11ea-018c-399409d7842e
begin
	f(x1,x2)=x1^2 + x2^2
	#f(x1,x2)=(x1+x2)^2
	#f(x1,x2)=-x1^2
	#f(x1,x2)=-2x1^2 - x2^2 - 2x1*x2
	#f(x1,x2)=-x1^2 + x2^2
end

# ╔═╡ f94dece0-defc-11ea-0546-39a094bbd2b6
begin
	gr()
	contour(x, y, f, fill = true)
	contour!(x,y,f,levels=[0.0],linewidth=3,color=:red, ratio=1)
end

# ╔═╡ 440c2010-df04-11ea-2550-6f679efd6b07
md"""# La Salle y región de atracción
Partimos del ejemplo

$$\dot x_1=(x_2-x_1)(1+x_1^2)^2$$

$$\dot x_2=-x_1+x_1^2$$

"""

# ╔═╡ c3aa90e0-df04-11ea-0736-5176e619708a
function f_ejemplo1(x,p,t)
	#p y t no se usan pero es para usar la función flechas del tema anterior
	x1,x2=x
	dx1=(x2-x1)*(1+x1^2)^2
    dx2=-x1+x1^2
	return [dx1;dx2]
end

# ╔═╡ 116e9d42-812f-11eb-327c-e107809f9603
md"## vamos a empezar a hacer cosas simbólicas (ahora que ya sabemos hacerlas a mano)"

# ╔═╡ d4547220-8121-11eb-1591-ebb8d7742fcc
function jacobiano(f,variables)
	N=length(variables)
	fun=f(variables,[],0)
	M=length(fun)
	J=[diff(fun[i],variables[j]) for i in 1:M, j in 1:N]
end

# ╔═╡ 5eebc400-8123-11eb-029e-8bd1bef4612e
x_1,x_2=symbols("x1,x2") #permite operar simbólicamente

# ╔═╡ f9dcadc0-812e-11eb-005f-d751ee820296
J=jacobiano(f_ejemplo1,[x_1,x_2])

# ╔═╡ f13ed140-81bd-11eb-200f-9fc7aa549363
f_ejemplo1([x_1,x_2],[],0)

# ╔═╡ 6eb00ba0-812b-11eb-06b8-9d782689846a
function linealizar(f,variables,puntos)
	#evalua el jacobiano en un punto
	J=jacobiano(f_ejemplo1,variables)
	for i in 1:length(variables)
		for j in 1:length(J)
		   J[j]=subs(J[j],variables[i],puntos[i])
		end
	end
	return J
end

# ╔═╡ 83941610-df41-11ea-0f0b-636bcab6c386
A=linealizar(f_ejemplo1,[x_1,x_2],[0,0])

# ╔═╡ c4a86fb0-df42-11ea-2ee5-63a5bdd03376
begin
	autovalores=eigvals(Float64.(A))
	md"Los autovalores son $autovalores"
end

# ╔═╡ 3cbc07a0-df43-11ea-0df8-5ba5311f2f43
md"""Sabemos que es estable (**localmente**)
para saber más derivamos $V_2$ del ejemplo anterior:

$$V_2=\frac{x_1^2}{1+x_1^2}+x_2^2$$

$$\dot V_2=\frac{2x_1}{(1+x_1^2)^2}\dot x_1+2x_2\dot x_2

=2x_1(x_2-x_1)+2x_2(-x_1+x_1^2)$$

$$=-2x_1^2(1-x_2)$$

De momento a mano en el tema 4 veremos cómo decirle a la máquina que lo haga

"""

# ╔═╡ c82134f0-8123-11eb-279a-2b03ca565783
md"V no aumenta luego es **Estable**, pero ¿Podemos saber algo más?"

# ╔═╡ a7b42e60-dfcb-11ea-2ae9-87311393dbef
md"""# Usemos LaSalle para ver la región de atracción
## Paso 1 Elegir L
Si elijo bien L V no aumenta en la región V(x)<L (las flechas no apuntan hacia fuera)"""

# ╔═╡ dd573c40-dfc3-11ea-3ec8-7d0b939123c1
dV2(x1,x2)=-2x1^2*(1-x2)

# ╔═╡ 9b580610-dfc5-11ea-2cfb-bb5db8322e86
@bind L3 Slider(0:0.05:2)

# ╔═╡ 5b862210-dfc5-11ea-3958-d5f618a37108
begin
	x3=-2:0.01:2
	y3=x3
	gr()
	contour(x3, y3, dV2, fill = false, linestyle=:dot)
	contour!(x3, y3, V2, fill = false)
	contour!(x3,y3,V2,levels=[L3],linewidth=3,color=:red)
	contour!(x3,y3,dV2,levels=[0],linewidth=3,color=:green)
	plot!(title="L=$(L3)")
end

# ╔═╡ d64c1460-f9d2-11ea-307b-c7d59f18b084
md"""**Ojo no siempre se detectan bien las curvas de nivel**

¿Echais algo en falta?"""

# ╔═╡ 94fc4d7e-f9d3-11ea-0b41-a1d22cd9376b
begin
	plotly()
	surface(x3,y3,dV2)
end

# ╔═╡ 698f68b0-dfcc-11ea-2a38-1fd314c52157
md"""Si elijo L=1 V no aumenta en V(x)<L.

Luego si empiezo en $V(x_0)<1$ entonces $V(x(t))<1$ para todo $t$
estamos atrapados en V<L.

V No aumenta y de hecho disminuye en algunos puntos ¿Dónde deja de disminuir? donde $$\dot V=0)$$"""

# ╔═╡ f2777960-dfcc-11ea-089a-af9cb98c26e7
md"""## Paso 2 calcular la región R
R es simplemente la región en la que $V<L$ y además $\dot V=0$
es decir:
$$\dot V_1=-2x_1^2(1-x_2)=0$$
de modo que:

A) $x_1=0$

B) $1-x_2=0$ esta **NO** puede ser ya que si $x_2=1$ entonces $V=\frac{x_1^2}{1+x_1^2}+1^2 \ge 1$ que está **fuera** de $V<1$

Sabemos entonces que $x_1 \to 0$
"""

# ╔═╡ 9370bbb0-dfcd-11ea-3079-d3d74eee20ea
md"""## Paso 3 calcular la región M
M es el conjunto de todas las solucines que viven en R.
Es decir tenemos que hacer $x_1(t)=0$ (y por tanto $\dot x_1=0$ )y ver que pasa al meterlo en las ecuaciones:

$$\dot x_1=$$

$$0=(x_2-0)(1+0^2)^2 \to x_2=0$$

$$\dot x_2=-0+0^2$$

Solo hay un punto en el que las soluciones *se quedan* en la región R que es (0,0) las otras soluciones "pasan de largo" (ver el dibujo)

"""

# ╔═╡ f9362140-dfcf-11ea-2882-77e3c638d875
md"""## Hemos demostrado que si empezamos en $V<1$ entonces $x \to 0$

O lo que es lo mismo que $V<1$ es una región de atración del origen"""

# ╔═╡ 503ad610-dfd6-11ea-1021-5fb2a45bc262
md"""# ¿Por qué ese V?
Podríamos haber usado $V_1$, veamos que sale:
$$V_1=x_1^2+x_2^2$$

$$\dot V_1=2x_1\dot x_1+2x_2\dot x_2

=2x_1(x_2-x_1)(1+x_1^2)^2+2x_2(-x_1+x_1^2)$$

V no aumenta luego es **Estable**, pero ¿Podemos saber algo más?

"""

# ╔═╡ bf615090-dfd7-11ea-12bb-11d922c37139
dV1(x1,x2)=2x1*(x2-x1)*(1+x1^2)^2+2x2*(-x1+x1^2)

# ╔═╡ f0dd5060-dfd7-11ea-03d8-63269e0c3394
@bind L4 Slider(0:0.05:2)

# ╔═╡ a4399020-dfd7-11ea-2e36-ddc90140e91d
begin
	x4=-2:0.01:2
	y4=x3
	gr()
	contour(x4, y4, dV1, fill = false)
	contour!(x4, y4, V1, fill = false)
	contour!(x4,y4,V1,levels=[L4],linewidth=3,color=:red)
	contour!(x4,y4,dV1,levels=[0],linewidth=3,color=:green)
	plot!(title="L=$(L4)",ratio=1)
end

# ╔═╡ 0edbfe90-dfd8-11ea-038a-97c39cddb672
md"""con L=0.6, V no aumenta

R es la región en la que \dot v=0 que al elegir el l adecuado se reduce a $x_1=0$ el resto es igual que el caso anterior.
"""

# ╔═╡ a96a7e30-d639-11ea-08df-5d727910849c
function flechas!(figura,f,p,rangos;N=10)
	xmin,xmax,ymin,ymax=rangos
	longitud=max(abs(xmax-xmin),abs(ymax-ymin))/N

	x1s=xmin:((xmax-xmin)/N):xmax
	x2s=ymin:((ymax-ymin)/N):ymax

	X1=[];
	X2=[];
	U=[];
	V=[];
	for x1 in x1s #barrido para calcular las flechas
		for x2 in x2s
			push!(X1,x1); #Push añade un componente al final de un vector
			push!(X2,x2);
			F=f([x1,x2],p,0)
			X=0.1*F/(0.1+sqrt(F[1]^2+F[2]^2)) #Saturo el tamaño de la flecha
			push!(U,X[1]);
			push!(V,X[2]);
		end
	end
	quiver!(figura,X1,X2,quiver=(U,V), arrow=arrow(0.1*longitud, 0.1*longitud))
	return figura
end

# ╔═╡ 93e7ebb0-dfc6-11ea-26f1-f7e281edf001
begin
	gr()
	p1=plot()
	flechas!(p1,f_ejemplo1,[],[-2,2,-2,2],N=20)
	contour!(x3,y3,V2,levels=[L3],linewidth=3,color=:red)
	contour!(x3,y3,dV2,levels=[0],linewidth=3,color=:green)
end

# ╔═╡ 8f8b3390-8130-11eb-38df-99e46c4051e0
function trayectoria(x0,f,p,rangos;N=10,pasosMax=100)
	xmin,xmax,ymin,ymax=rangos
	long=max(abs(xmax-xmin),abs(ymax-ymin))/N
	x,y=x0
	pasos=0

	X1=[x]
	Y1=[y]
	X2=[x] #parte de la trayectoria hacia atrás en el tiempo
	Y2=[y]

	while (xmin<x<xmax) & (ymin<y<ymax) & (pasos<pasosMax)
		derivada=f([x,y],p,0)
		dt=0.05*long/(0.1+sqrt(derivada[1]^2+derivada[2]^2))
		Xn=[x;y]+derivada*dt
		x,y=Xn
		pasos=pasos+1
		push!(X1,x)
		push!(Y1,y)
	end
	pasos=0
	x,y=x0
	while (xmin<x<xmax) & (ymin<y<ymax) & (pasos<pasosMax)
		derivada=f([x,y],p,0)
		dt=-0.05*long/(0.1+sqrt(derivada[1]^2+derivada[2]^2))
		Xn=[x;y]+derivada*dt
		x,y=Xn
		pasos=pasos+1
		push!(X2,x)
		push!(Y2,y)
	end
	X=append!(reverse(X2),X1) #junto ambas partes
	Y=append!(reverse(Y2),Y1)
	return (X,Y)
end

# ╔═╡ e4e5a27e-8130-11eb-0c8a-91d509d7e897
function trayectorias!(figura,f,p,rangos;N=10)

	xmin,xmax,ymin,ymax=rangos
	N=round(N/2);
	x1s=xmin:((xmax-xmin)/N):xmax
	x2s=ymin:((ymax-ymin)/N):ymax
	T=10;
	for x1 in x1s
		for x2 in x2s
			X,Y=trayectoria([x1;x2],f,p,rangos,N=N,pasosMax=100)
			figura=plot!(figura,X,Y,legend=false)
		end
	end
	return figura
end

# ╔═╡ 42ab643e-dfd8-11ea-0ed7-392d94c7904e
begin
	gr()
	p=plot()
	flechas!(p,f_ejemplo1,[],[-2,2,-2,2],N=20)
	trayectorias!(p,f_ejemplo1,[],[-2,2,-2,2],N=30)
	contour!(x3,y3,V2,levels=[L3],linewidth=3,color=:red)
	contour!(x3,y3,V1,levels=[L4],linewidth=3,color=:red)
	plot!(title="Comparación de regiones",ratio=1)
end

# ╔═╡ Cell order:
# ╠═02c15be0-defd-11ea-004b-45c4e09aa541
# ╠═8aeae030-811f-11eb-0075-272d383103ef
# ╟─28fccf20-d63c-11ea-3311-411f89a63431
# ╟─309959f0-defd-11ea-34c3-71eeda7b523b
# ╠═cd88fade-defd-11ea-2097-25b2d1eac285
# ╠═ce5fe800-deff-11ea-0b19-7fc23b14bd06
# ╠═c6c49460-deff-11ea-1c6d-673d863ea07b
# ╠═5b26b050-d7eb-11ea-2488-b7da296fed24
# ╠═2e580002-defe-11ea-1b3d-d5ec461e0715
# ╟─2ec228f0-defd-11ea-329f-53576ab14004
# ╟─ca2fea72-3655-11eb-0984-a532e1f9efcc
# ╠═204e6710-defe-11ea-1fb0-43c7ef863786
# ╠═f3bcc450-df00-11ea-0080-03b54212b9ed
# ╠═ca99149e-df02-11ea-39a0-179e6a7d1344
# ╠═eb00be70-df00-11ea-0ce9-3b0a79a0f4a3
# ╟─d8678810-defc-11ea-05b2-51373f754e61
# ╠═34a9ca60-df03-11ea-018c-399409d7842e
# ╠═f94dece0-defc-11ea-0546-39a094bbd2b6
# ╟─440c2010-df04-11ea-2550-6f679efd6b07
# ╠═c3aa90e0-df04-11ea-0736-5176e619708a
# ╟─116e9d42-812f-11eb-327c-e107809f9603
# ╠═ba6bf300-8122-11eb-0a3e-2d30a150250c
# ╠═d4547220-8121-11eb-1591-ebb8d7742fcc
# ╠═5eebc400-8123-11eb-029e-8bd1bef4612e
# ╠═f9dcadc0-812e-11eb-005f-d751ee820296
# ╠═f13ed140-81bd-11eb-200f-9fc7aa549363
# ╠═6eb00ba0-812b-11eb-06b8-9d782689846a
# ╠═83941610-df41-11ea-0f0b-636bcab6c386
# ╠═6cbb32f0-d665-11ea-3f24-41cbe1910b80
# ╠═c4a86fb0-df42-11ea-2ee5-63a5bdd03376
# ╟─3cbc07a0-df43-11ea-0df8-5ba5311f2f43
# ╟─c82134f0-8123-11eb-279a-2b03ca565783
# ╟─a7b42e60-dfcb-11ea-2ae9-87311393dbef
# ╠═dd573c40-dfc3-11ea-3ec8-7d0b939123c1
# ╠═5b862210-dfc5-11ea-3958-d5f618a37108
# ╠═9b580610-dfc5-11ea-2cfb-bb5db8322e86
# ╟─d64c1460-f9d2-11ea-307b-c7d59f18b084
# ╠═94fc4d7e-f9d3-11ea-0b41-a1d22cd9376b
# ╠═93e7ebb0-dfc6-11ea-26f1-f7e281edf001
# ╟─698f68b0-dfcc-11ea-2a38-1fd314c52157
# ╟─f2777960-dfcc-11ea-089a-af9cb98c26e7
# ╟─9370bbb0-dfcd-11ea-3079-d3d74eee20ea
# ╟─f9362140-dfcf-11ea-2882-77e3c638d875
# ╟─503ad610-dfd6-11ea-1021-5fb2a45bc262
# ╠═bf615090-dfd7-11ea-12bb-11d922c37139
# ╠═f0dd5060-dfd7-11ea-03d8-63269e0c3394
# ╠═a4399020-dfd7-11ea-2e36-ddc90140e91d
# ╟─0edbfe90-dfd8-11ea-038a-97c39cddb672
# ╟─a96a7e30-d639-11ea-08df-5d727910849c
# ╟─8f8b3390-8130-11eb-38df-99e46c4051e0
# ╟─e4e5a27e-8130-11eb-0c8a-91d509d7e897
# ╠═42ab643e-dfd8-11ea-0ed7-392d94c7904e

max_stars_repo_path: Tema 2.jl

max_stars_repo_name: Dictino/CNL

max_issues_repo_path: Tema 2.jl

max_issues_repo_name: Dictino/CNL

max_forks_repo_path: Tema 2.jl

max_forks_repo_name: Dictino/CNL


text:
# ---
# title: 12. Integer to Roman
# id: problem12
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, String
# link: <https://leetcode.com/problems/integer-to-roman/description/>
# hidden: true
# ---
#
# Roman numerals are represented by seven different symbols: `I`, `V`, `X`, `L`,
# `C`, `D` and `M`.
#
#
#
#     **Symbol**       **Value**
#     I             1
#     V             5
#     X             10
#     L             50
#     C             100
#     D             500
#     M             1000
#
# For example, `2` is written as `II` in Roman numeral, just two one's added
# together. `12` is written as `XII`, which is simply `X + II`. The number `27`
# is written as `XXVII`, which is `XX + V + II`.
#
# Roman numerals are usually written largest to smallest from left to right.
# However, the numeral for four is not `IIII`. Instead, the number four is
# written as `IV`. Because the one is before the five we subtract it making
# four. The same principle applies to the number nine, which is written as `IX`.
# There are six instances where subtraction is used:
#
#   * `I` can be placed before `V` (5) and `X` (10) to make 4 and 9.
#   * `X` can be placed before `L` (50) and `C` (100) to make 40 and 90.
#   * `C` can be placed before `D` (500) and `M` (1000) to make 400 and 900.
#
# Given an integer, convert it to a roman numeral.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 3
#     Output: "III"
#
#
# **Example 2:**
#
#
#
#     Input: num = 4
#     Output: "IV"
#
#
# **Example 3:**
#
#
#
#     Input: num = 9
#     Output: "IX"
#
#
# **Example 4:**
#
#
#
#     Input: num = 58
#     Output: "LVIII"
#     Explanation: L = 50, V = 5, III = 3.
#
#
# **Example 5:**
#
#
#
#     Input: num = 1994
#     Output: "MCMXCIV"
#     Explanation: M = 1000, CM = 900, XC = 90 and IV = 4.
#
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 3999`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function optimize!(model::Model)
  @assert is_valid(model.env)
  stat = (if model.has_int
    if VERSION < v"0.7.0-DEV.3382"
        @cpx_ccall_intercept(model, mipopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    else
        @cpx_ccall(mipopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    end
  elseif model.has_qc
    if VERSION < v"0.7.0-DEV.3382"
        @cpx_ccall_intercept(model, qpopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    else
        @cpx_ccall(qpopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    end
  else
    if VERSION < v"0.7.0-DEV.3382"
        @cpx_ccall_intercept(model, lpopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    else
        @cpx_ccall(lpopt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp)
    end
  end)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
end

set_branching_priority(model::Model, priority) =
    set_branching_priority(model, Cint[1:num_var(model)], priority)

set_branching_priority(model::Model, indices, priority) =
    set_branching_priority(model, indices, priority, C_NULL)

set_branching_priority(model, indices, priority, direction) =
    set_branching_priority(model, convert(Vector{Cint},indices), convert(Vector{Cint},priority), direction)

function set_branching_priority(model::Model, indices::Vector{Cint}, priority::Vector{Cint}, direction)
    @assert (cnt = length(indices)) == length(priority)
    isa(direction,Vector) && @assert cnt == length(direction)
    stat = @cpx_ccall(copyorder, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cint}),
                      model.env.ptr, model.lp, cnt, indices .- Cint(1), priority, direction)
    stat == 0 || throw(CplexError(model.env, stat))
    return nothing
end

function newlongannotation(model::Model, name::String, defval::Clong)
    stat = @cpx_ccall(newlongannotation, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cchar},
                      Clong),
                      model.env.ptr, model.lp, name, defval)
    stat == 0 || throw(CplexError(model.env, stat))

    return nothing
end

function setlongannotations(model::Model, idx::Cint, objtype::Cint, cnt::Cint, indexArr::Array{Cint},
                valArr::Array{Clong})
    stat = @cpx_ccall(setlongannotations, Cint, (
                Ptr{Cvoid},
                Ptr{Cvoid},
                Cint,
                Cint,
                Cint,
                Ptr{Cvoid},
                Ptr{Cvoid}),
                model.env.ptr, model.lp, idx, objtype, cnt, indexArr, valArr)
    stat == 0 || throw(CplexError(model.env, stat))
    return nothing
end

export setlongannotations, newlongannotation

function c_api_getobjval(model::Model)
  objval = Vector{Cdouble}(undef, 1)
  stat = @cpx_ccall(getobjval, Cint, (
                    Ptr{Cvoid},
                    Ptr{Cvoid},
                    Ptr{Cdouble}
                    ),
                    model.env.ptr, model.lp, objval)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  return objval[1]
end
get_objval(model::Model) = c_api_getobjval(model)

function c_api_solninfo(model::Model)
  solnmethod_p = Ref{Cint}()
  solntype_p = Ref{Cint}()
  pfeasind_p = Ref{Cint}()
  dfeasind_p = Ref{Cint}()
  stat = @cpx_ccall(solninfo, Cint, (
                    Ptr{Cvoid},
                    Ptr{Cvoid},
                    Ptr{Cint},
                    Ptr{Cint},
                    Ptr{Cint},
                    Ptr{Cint}
                    ),
                    model.env.ptr, model.lp, solnmethod_p, solntype_p,
                    pfeasind_p, dfeasind_p)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  return (solnmethod_p[], solntype_p[], pfeasind_p[], dfeasind_p[])
end

function c_api_getx(model::Model, x::FVec)
  nvars = num_var(model)
  stat = @cpx_ccall(getx, Cint, (
                    Ptr{Cvoid},
                    Ptr{Cvoid},
                    Ptr{Cdouble},
                    Cint,
                    Cint
                    ),
                    model.env.ptr, model.lp, x, 0, nvars-Cint(1))
  if stat != 0
    throw(CplexError(model.env, stat))
  end
end

function get_solution(model::Model)
  nvars = num_var(model)
  x = Vector{Cdouble}(undef, nvars)
  c_api_getx(model, x)
  return x
end

function c_api_getdj(model::Model, p::FVec)
    nvars = num_var(model)
    stat = @cpx_ccall(getdj, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, p, 0, nvars-1)
    if stat != 0
       throw(CplexError(model.env, stat))
    end
end

function get_reduced_costs(model::Model)
    nvars = num_var(model)
    p = Vector{Cdouble}(undef, nvars)
    c_api_getdj(model, p)
    return p
end

function c_api_getpi(model::Model, p::FVec)
    ncons = num_constr(model)
    stat = @cpx_ccall(getpi, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, p, 0, ncons-1)
    if stat != 0
       throw(CplexError(model.env, stat))
    end
end

function c_api_getqconstrslack(model::Model, dest::Vector{Float64})
    ncons = num_qconstr(model)
    stat = @cpx_ccall(getqconstrslack, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, dest, 0, ncons-1)
    if stat != 0
       throw(CplexError(model.env, stat))
    end
end

function get_constr_duals(model::Model)
    ncons = num_constr(model)
    p = Vector{Cdouble}(undef, ncons)
    c_api_getpi(model, p)
    return p
end

function c_api_getax(model::Model, Ax::FVec)
    ncons = num_constr(model)
    stat = @cpx_ccall(getax, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, Ax, 0, ncons-1)
    if stat != 0
      throw(CplexError(model.env, stat))
    end
end

function c_api_getxqxax(model::Model, dest::Vector{Float64})
    ncons = num_qconstr(model)
    stat = @cpx_ccall(getxqxax, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, dest, 0, ncons-1)
    if stat != 0
      throw(CplexError(model.env, stat))
    end
end

function get_constr_solution(model::Model)
    ncons = num_constr(model)
    Ax = Vector{Cdouble}(undef, ncons)
    c_api_getax(model, Ax)
    return Ax
end

function get_infeasibility_ray(model::Model)
  ncons = num_constr(model)
  y = Vector{Cdouble}(undef, ncons)
  proof_p = Vector{Cdouble}(undef, 1)
  stat = @cpx_ccall(dualfarkas, Cint, (
                    Ptr{Cvoid},
                    Ptr{Cvoid},
                    Ptr{Cdouble},
                    Ptr{Cdouble}
                    ),
                    model.env.ptr, model.lp, y, proof_p)
  if stat != 0
    error("CPLEX is unable to grab infeasible ray; consider resolving with presolve turned off")
    throw(CplexError(model.env, stat))
  end
  return y
end

function get_unbounded_ray(model::Model)
  solve_stat = get_status(model)
  if solve_stat == :CPX_STAT_UNBOUNDED
    n = num_var(model)
    z = Vector{Cdouble}(undef, n)
    stat = @cpx_ccall(getray, Cint, (
                      Ptr{Cvoid},
                      Ptr{Cvoid},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, z)
    if stat != 0
      throw(CplexError(model.env, stat))
    end
    return z
  else
    error("CPLEX is unable to grab unbounded ray; consider resolving with presolve turned off")
  end
end

const varmap = Dict(
    0 => :NonbasicAtLower,
    1 => :Basic,
    2 => :NonbasicAtUpper,
    3 => :Free
)

const conmap = Dict(
    0 => :NonbasicAtLower,
    1 => :Basic,
    2 => :NonbasicAtUpper
)

function get_basis(model::Model)
    cval = Vector{Cint}(undef, num_var(model))
    rval = Vector{Cint}(undef, num_constr(model))
    stat = @cpx_ccall(getbase, Cint, (Ptr{Cvoid},Ptr{Cvoid},Ptr{Cint},Ptr{Cint}),
                      model.env.ptr, model.lp, cval, rval)
    stat != 0 && throw(CplexError(model.env, stat))

    csense = get_constr_senses(model)
    cbasis = Vector{Symbol}(undef, num_var(model))
    rbasis = Vector{Symbol}(undef, num_constr(model))
    for it in 1:num_var(model)
        cbasis[it] = varmap[cval[it]]
    end
    for it in 1:num_constr(model)
        rbasis[it] = conmap[rval[it]]
        if (rbasis[it] == :NonbasicAtLower) && (csense[it] == convert(Cchar,'L'))
            rbasis[it] = :NonbasicAtUpper
        end
    end
    return cbasis, rbasis
end

get_node_count(model::Model) = @cpx_ccall(getnodecnt, Cint, (Ptr{Cvoid},Ptr{Cvoid}), model.env.ptr, model.lp)

function get_rel_gap(model::Model)
  ret = Vector{Cdouble}(undef, 1)
  stat = @cpx_ccall(getmiprelgap, Cint, (Ptr{Cvoid},Ptr{Cvoid},Ptr{Cdouble}), model.env.ptr, model.lp, ret)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  ret[1]
end

function get_num_cuts(model::Model,cuttype)
    cutcount = Vector{Cint}(undef, 1)

    stat = @cpx_ccall(getnumcuts, Cint, (Ptr{Cvoid},Ptr{Cvoid},Cint,Ptr{Cvoid}), model.env.ptr , model.lp, cuttype, cutcount)
    if stat != 0
        error(CplexError(model.inner.env, stat).msg)
    end
    return cutcount[1]
end

const status_symbols = Dict(
    1   => :CPX_STAT_OPTIMAL,
    2   => :CPX_STAT_UNBOUNDED,
    3   => :CPX_STAT_INFEASIBLE,
    4   => :CPX_STAT_INForUNBD,
    5   => :CPX_STAT_OPTIMAL_INFEAS,
    6   => :CPX_STAT_NUM_BEST,
    7   => :CPX_STAT_FEASIBLE_RELAXED,
    8   => :CPX_STAT_OPTIMAL_RELAXED,
    10  => :CPX_STAT_ABORT_IT_LIM,
    11  => :CPX_STAT_ABORT_TIME_LIM,
    12  => :CPX_STAT_ABORT_OBJ_LIM,
    13  => :CPX_STAT_ABORT_USER,
    20  => :CPX_STAT_OPTIMAL_FACE_UNBOUNDED,
    21  => :CPX_STAT_ABORT_PRIM_OBJ_LIM,
    22  => :CPX_STAT_ABORT_DUAL_OBJ_LIM,
    101 => :CPXMIP_OPTIMAL,
    102 => :CPXMIP_OPTIMAL_TOL,
    103 => :CPXMIP_INFEASIBLE,
    104 => :CPXMIP_SOL_LIM,
    105 => :CPXMIP_NODE_LIM_FEAS,
    106 => :CPXMIP_NODE_LIM_INFEAS,
    107 => :CPXMIP_TIME_LIM_FEAS,
    108 => :CPXMIP_TIME_LIM_INFEAS,
    109 => :CPXMIP_FAIL_FEAS,
    110 => :CPXMIP_FAIL_INFEAS,
    111 => :CPXMIP_MEM_LIM_FEAS,
    112 => :CPXMIP_MEM_LIM_INFEAS,
    113 => :CPXMIP_ABORT_FEAS,
    114 => :CPXMIP_ABORT_INFEAS,
    115 => :CPXMIP_OPTIMAL_INFEAS,
    116 => :CPXMIP_FAIL_FEAS_NO_TREE,
    117 => :CPXMIP_FAIL_INFEAS_NO_TREE,
    118 => :CPXMIP_UNBOUNDED,
    119 => :CPXMIP_INForUNBD,
    120 => :CPXMIP_FEASIBLE_RELAXED,
    121 => :CPXMIP_OPTIMAL_RELAXED
)

get_status(model::Model) = status_symbols[Int(get_status_code(model))]::Symbol

function c_api_getstat(model::Model)
    return @cpx_ccall(getstat, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                      model.env.ptr, model.lp)
end
get_status_code(model::Model) = c_api_getstat(model)

function c_api_getstatstring(model::Model, statind::Cint)
  buffer_str = Vector{Cchar}(undef, 1024)
  @cpx_ccall(
    getstatstring,
    Ptr{Cvoid},
    (Ptr{Cvoid}, Cint, Ptr{Cchar}),
    model.env.ptr, statind, buffer_str
  )
  return unsafe_string(pointer(buffer_str))
end

function c_api_getitcnt(model::Model)
    return @cpx_ccall(
      getitcnt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp
    )
end

function c_api_getbaritcnt(model::Model)
  return @cpx_ccall(
    getbaritcnt, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), model.env.ptr, model.lp
  )
end


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

include("formatting.jl")

const margin = 2
cols(io) = displaysize(io)[2]

function term(io::IO, content::Vector, cols)
    isempty(content) && return
    for md in content[1:end-1]
        term(io, md, cols)
        println(io)
    end
    term(io, content[end], cols)
end

term(io::IO, md::MD, columns = cols(io)) = term(io, md.content, columns)

function term(io::IO, md::Paragraph, columns)
    print(io, " "^margin)
    print_wrapped(io, width = columns-2margin, pre = " "^margin) do io
        terminline(io, md.content)
    end
end

function term(io::IO, md::BlockQuote, columns)
    s = sprint(term, md.content, columns - 10)
    for line in split(rstrip(s), "\n")
        println(io, " "^margin, "|", line)
    end
end

function term(io::IO, md::Admonition, columns)
    print(io, " "^margin, "| ")
    with_output_format(:bold, print, io, isempty(md.title) ? md.category : md.title)
    println(io, "\n", " "^margin, "|")
    s = sprint(term, md.content, columns - 10)
    for line in split(rstrip(s), "\n")
        println(io, " "^margin, "|", line)
    end
end

function term(io::IO, f::Footnote, columns)
    print(io, " "^margin, "| ")
    with_output_format(:bold, print, io, "[^$(f.id)]")
    println(io, "\n", " "^margin, "|")
    s = sprint(term, f.text, columns - 10)
    for line in split(rstrip(s), "\n")
        println(io, " "^margin, "|", line)
    end
end

function term(io::IO, md::List, columns)
    for (i, point) in enumerate(md.items)
        print(io, " "^2margin, isordered(md) ? "$(i + md.ordered - 1). " : "•  ")
        print_wrapped(io, width = columns-(4margin+2), pre = " "^(2margin+2),
                          i = 2margin+2) do io
            term(io, point, columns - 10)
        end
    end
end

function _term_header(io::IO, md, char, columns)
    text = terminline(md.text)
    with_output_format(:bold, io) do io
        print(io, " "^(margin))
        line_no, lastline_width = print_wrapped(io, text,
                                                width=columns - 4margin; pre=" ")
        line_width = min(1 + lastline_width, columns)
        if line_no > 1
            line_width = max(line_width, div(columns, 3))
        end
        char != ' ' && println(io, " "^(margin), string(char) ^ line_width)
    end
end

const _header_underlines = collect("≡=–-⋅ ")
# TODO settle on another option with unicode e.g. "≡=≃–∼⋅" ?

function term(io::IO, md::Header{l}, columns) where l
    underline = _header_underlines[l]
    _term_header(io, md, underline, columns)
end

function term(io::IO, md::Code, columns)
    with_output_format(:cyan, io) do io
        for line in lines(md.code)
            print(io, " "^margin)
            println(io, line)
        end
    end
end

function term(io::IO, br::LineBreak, columns)
   println(io)
end

function term(io::IO, br::HorizontalRule, columns)
   println(io, " " ^ margin, "-" ^ (columns - 2margin))
end

term(io::IO, x, _) = show(io, MIME"text/plain"(), x)

# Inline Content

terminline(md) = sprint(terminline, md)

function terminline(io::IO, content::Vector)
    for md in content
        terminline(io, md)
    end
end

function terminline(io::IO, md::AbstractString)
    print(io, replace(md, r"[\s\t\n]+", " "))
end

function terminline(io::IO, md::Bold)
    with_output_format(:bold, terminline, io, md.text)
end

function terminline(io::IO, md::Italic)
    with_output_format(:underline, terminline, io, md.text)
end

function terminline(io::IO, md::LineBreak)
    println(io)
end

function terminline(io::IO, md::Image)
    terminline(io, "(Image: $(md.alt))")
end

terminline(io::IO, f::Footnote) = with_output_format(:bold, terminline, io, "[^$(f.id)]")

function terminline(io::IO, md::Link)
    terminline(io, "$(terminline(md.text)) ($(md.url))")
end

function terminline(io::IO, code::Code)
    print_with_format(:cyan, io, code.code)
end

terminline(io::IO, x) = show(io, MIME"text/plain"(), x)

# Show in terminal

Base.display(d::Base.REPL.REPLDisplay, md::MD) = term(Base.REPL.outstream(d.repl), md)

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
#!/usr/bin/env julia
# Rendering slices and an isosurface of MRI data.

import GR
const gr3 = GR.gr3

function draw(mesh::Cint; x::Union{Real, Void}=nothing, y::Union{Real, Void}=nothing, z::Union{Real, Void}=nothing)
    gr3.clear()
    gr3.drawmesh(mesh, 1, (0,0,0), (0,0,1), (0,1,0), (1,1,1), (1,1,1))
    gr3.drawslicemeshes(data, x=x, y=y, z=z)
    GR.clearws()
    gr3.drawimage(0, 1, 0, 1, 500, 500, gr3.DRAWABLE_GKS)
    GR.updatews()
end

data = open(stream -> read!(stream, Array{UInt16}(93, 64, 64)), "mri.raw")
data = min.(data, 2000) / 2000.0 * typemax(UInt16)
data = convert(Array{UInt16, 3}, floor.(data))
data = permutedims(data, [3, 2, 1])

GR.setviewport(0, 1, 0, 1)
gr3.cameralookat(-3, 2, -2, 0, 0, 0, 0, 0, -1)

mesh = gr3.createisosurfacemesh(data, (2.0/63, 2.0/63, 2.0/92), (-1.0, -1.0, -1.0), 40000)

GR.setcolormap(1)
for z in linspace(0, 1, 300)
    draw(mesh, x=0.9, z=z)
end
for y in linspace(1, 0.5, 300)
    draw(mesh, x=0.9, y=y, z=1)
end
GR.setcolormap(19)
for x in linspace(0.9, 0, 300)
    draw(mesh, x=x, y=0.5, z=1)
end
for x in linspace(0, 0.9, 300)
    draw(mesh, x=x, z=1)
end
gr3.terminate()


text:
using Base.Test
using JuMP
using MIPVerify
isdefined(:TestHelpers) || include("../../TestHelpers.jl")
using TestHelpers: get_new_model

@testset "flatten.jl" begin

    @testset "Flatten" begin
        @testset "default permutation" begin
            p = Flatten(5)
            @test p.n_dim == 5
            @test all(p.perm .== [5, 4, 3, 2, 1])
        end
        @testset "specified permutation, legal" begin
            p = Flatten([3, 1, 2, 4])
            @test p.n_dim == 4
            @test all(p.perm .== [3, 1, 2, 4])
        end
        @testset "specified permutation, illegal" begin
            @test_throws DomainError Flatten([3, 1, 4])
        end
        @testset "Base.show" begin
            p = Flatten([3, 1, 2, 4])
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) == "Flatten(): flattens 4 dimensional input, with dimensions permuted according to the order [3, 1, 2, 4]"
        end
    end

    @testset "flatten" begin
        srand(31415)
        xs = rand(1:5, (2, 2, 2, 2))
        p1 = Flatten([3, 1, 2, 4])
        @test all(p1(xs) .== [4, 3, 3, 1, 4, 5, 3, 5, 2, 4, 5, 5, 4, 2, 5, 4])
        p2 = Flatten([1, 3, 4, 2])
        @test p2(xs) == [4, 3, 3, 1, 2, 5, 4, 5, 4, 3, 5, 5, 4, 5, 2, 4]
    end

end


text:
using Test, YaoBlocks, YaoArrayRegister

@testset "basic" begin
    reg = rand_state(10)
    block = kron(4, 2 => X)
    c = subroutine(10, block, [1, 3, 9, 2])
    @test nqubits(c) == 10
    @test nactive(c) == 4
    @test isunitary(c) == true
    @test isreflexive(c) == true
    @test ishermitian(c) == true
    blk = kron(4, 2 => Rx(0.3))
    @test chsubblocks(c, [blk]) |> subblocks |> first == blk
    @test occupied_locs(c) == (3,)

    @test apply!(copy(reg), c) == apply!(copy(reg), kron(10, 3 => X))
    @test apply!(rand_state(12, nbatch = 10) |> focus!(Tuple(1:10)...), c) |> nactive == 10
end

@testset "test repeat" begin
    c = subroutine(8, repeat(5, H), 1:5)
    r = rand_state(8)
    r1 = copy(r) |> c
    r2 = copy(r) |> repeat(8, H, 1:5)
    @test r1 ≈ r2
end

@testset "mat" begin
    cc = subroutine(5, kron(X, X), (1, 3))
    @test applymatrix(cc) ≈ mat(cc)
end

@testset "push coverage" begin
    @test subroutine(put(3, 2=>X), (3,4,5))(5) isa Subroutine
    @test subroutine(put(2=>X), (3,4,5))(5) isa Subroutine
    s = subroutine(put(2=>X), (3,4,5))(5)
    q = subroutine(put(3=>Y), (3,4,5))(5)
    z = subroutine(put(3=>Val(im)*X), (3,4,5))(5)
    @test s' == s
    @test iscommute(s', s)
    @test iscommute(s', z)
    @test_throws LocationConflictError subroutine(5, put(3,2=>X), (1,4,2,5))
    @test_throws ErrorException subroutine(5, put(3,2=>X), (7,4,5))
    @test YaoBlocks.PropertyTrait(s) == YaoBlocks.PreserveAll()
end


text:
module Stiffness

using NPZ
using Glob
using SuperStiff.PeriodizeSC

###########################################Functions and structures###########################################
global pwd_ = pwd()
global verbose = 0

macro assertion(ex, text)
    :($ex ? nothing : error("Assertion failed: ", $(text)))
end

struct File_Input_Error <: Exception
    var::Symbol
end

# Exception constructor
Base.showerror(io::IO, e::File_Input_Error) = println(io, e.var, " Keys are not balanced (See function read_data_loop)")

# Static type object containing all important parameters
mutable struct StiffnessArray{T, N, S <: AbstractArray, Y <: Associative{String,Any}} <: AbstractArray{T, N}
    # Contains the data array
    data_::S
    # Can have maximum three non-zero hopping parameters
    hopping_val_::Base.RefValue{Tuple{T, T, T, T}}
    # Contains informations about the calculations to do
    vals_::Y
end

# Outer constructor:
StiffnessArray{T,N}(data_::AbstractArray{Complex{T},N}, hopping_val_::Tuple{T,T,T,T}, vals_::Associative{String,Any}) = StiffnessArray{T,N,typeof(data_),typeof(vals_)}(data_, Ref(hopping_val_), vals_)
ts(t::StiffnessArray) = t.hopping_val_[] # Getter
Base.size(A::StiffnessArray) = size(A.data_)
Base.IndexStyle{T<:StiffnessArray}(::Type{T}) = Base.IndexLinear()

function Base.getindex(A::StiffnessArray, I...)
    checkbounds(Bool, A.data_, I...) && return A.data_[I...]
    I... < 1 ? throw(ErrorException("Out of bounds (lower bound)")) : throw(ErrorException("Out of bounds (upper bound)"))
end

"""
Function used to build main dictionnary containing file-related information.

#Argument(s):

- block_params::Dict{String,Any}: Dict-valued argument produced form JSON parsing.
- print_mu_dop::Int64: Binary-valued argument (1 or 0) to decide to print (or not) a file describing chemical potential vs doping
- pattern::String: String-valued argument representing the extension common to the python binary files (*.npy)

#Returns:

- super_datap: 2d array containing chemical potential on first column and doping on second one (Array{Float64,2})
- stiffness.dat: If print_mu_dop = 1, returns a file and exits the program
- list_of_files: list of files from which the self-energies are to be extracted (Array{String,1})
"""
function read_data_loop(block_params::Dict{String,Any}, print_mu_dop::Int64, pattern::String)
    dict_data_file = Dict{String,Any}()

    data_file = open(readdlm, block_params["data_loop"])
    data_file_header = data_file[1,:]

    indmu = find(x->x=="mu",data_file_header)
    inddop = find(x->x=="ave_mu",data_file_header)
    indM = find(x->x=="ave_M",data_file_header)

    data_file_mu = data_file[:,indmu]
    data_file_dop = data_file[:,inddop]
    data_file_M = data_file[:,indM]
    length(data_file_mu) != length(data_file_dop) && throw(ErrorException("Length of chemical potential list must be the same length as that of the doping!"))

    data_file_mu = filter(x->x!="mu",data_file_mu)
    data_file_dop = filter(x->x!="ave_mu",data_file_dop)
    data_file_M = filter(x->x!="ave_M",data_file_M)
    datap_h = hcat(data_file_mu,data_file_dop)
    datap_h = convert(Array{Float64,2},datap_h)
    data_file_M = convert(Array{Float64,1},data_file_M)

    list_of_files = glob(string(block_params["path_to_files"],pattern))
    # Exits the program if print_mu_dop = 1
    if print_mu_dop == 1
        writedlm("stiffness.dat", datap_h, "\t\t")
        print_with_color(:red, "Printed stiffness.dat for later use (SE_G_fig_producer.py)\n")
        exit(0)
    else
        nothing
    end
    #println(length(data_file_M), "\n", length(list_datap_h), "\n", length(list_of_files))
    return data_file_M, datap_h, list_of_files
end

"""
Function extracting the last number of the binary filename.

#Argument(s):

- list_of_files::Array{String,1}: 1d array-valued argument containing the filenames (Array{String,1})

#Returns:

- list_num: 1d array-valued output containing the last number of the filenames (Array{Int64,1})
"""
function gen_file_num(list_of_files::Array{String,1})
    list_num = Array{Int64,1}(length(list_of_files))
    for (i,files) in enumerate(list_of_files)
        m = [x for x in eachmatch(r"\d+",files)][end]
        list_num[i] = parse(Int64,m.match)
    end
    if verbose > 0
        println("Length list_num: ", length(list_num))
    end
    return list_num
end

"""
Function producing the arguments defining any instance (to initiate the instance) of type PeriodizeSC.ModelVector.

#Argument(s):

- block_params::Dict{String,Any}: Directory containing a specific set of the input parameters extracted from params.json (Dict{String,Any})
- list_of_files::Array{String,1}: 1d array-valued argument containing the filenames (Array{String,1})
- data_mu::Array{Float64,1}: 1d array containing the chemical potentials (Array{Float64,1})
- params::Dict{String,Any}: Directory containing all of the input parameters extracted from params.json (Dict{String,Any})

#Returns:

- list_modulevec: 1d array containing all the ModelVector objects holding the informations read out of each of the files
"""
function gen_modulevec_args(block_params::Dict{String,Any},list_of_files::Array{String,1},data_mu::Array{Float64,1},params::Dict{String,Any})
    list_modulevec = Array{StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}},1}() #; list_t_mu = Array{Tuple{Float64,Float64,Float64,Float64},1}()
    list_num = gen_file_num(list_of_files)
    zip_num_files = collect(zip(list_num,list_of_files))
    sort!(zip_num_files)
    for (num,files) in zip_num_files
        sEvec_c = npzread(files)
        if verbose > 0
            println(num, " : ", (params["t"],params["tp"],params["tpp"],data_mu[num+1]))
        end
        println("sEvec_c: ", typeof(sEvec_c))
        println("Types of variables: ", typeof(StiffnessArray(sEvec_c, (params["t"],params["tp"],params["tpp"],data_mu[num+1]), block_params)))
        push!(list_modulevec, StiffnessArray(sEvec_c, (params["t"],params["tp"],params["tpp"],data_mu[num+1]), block_params)) ## Pushing the data structure here !!!!!
    end
    return list_modulevec
end

"""
Function calculating the superfluid stiffness calling the module SuperStiff.PeriodizeSC in the regime where there is coexistence.
Function to use (most suitable) when computing the c-axis superfluid stiffness. Can also compute the a- and b-axis superfluid stiffness, but with
performance drawbacks.

#Argument(s):

- super_data_M_el::Float64: Value of the AFM order parameter amplitude (Float64)
- modelvec::PeriodizeSC.ModelVector: ModelVector instance of a particular set of informations read out from a file (PeriodizeSC.ModelVector)
- modulevec_el::Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}}: Data structure to handle the different input parameters extracted from params.json (Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}})
- M_tol: Tolerance of the AFM order parameter amplitude above which AFM is considered to have set in (For debugging purposes only)

#Returns:

- list_stiff: 2d array containing the values of the superfluid stiffness before summing over Matsubara frequencies (Matrix{Complex{Float64}})
"""
function calc_stiff_funct_COEX(super_data_M_el::Float64, modelvec::PeriodizeSC.ModelVector, modulevec_el::Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}}, M_tol::Float64)
    list_stiff = Matrix{Complex{Float64}}(0,0)
    cumulant = modulevec_el.vals_["cumulant"]
    println("c-axis computations")
    if modulevec_el.vals_["Periodization"] == 1
        println("Periodization option set to 1")
        if cumulant == 0
            cond = abs(super_data_M_el) > M_tol
            if cond
                list_stiff = PeriodizeSC.calcintegral_RBZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_test)
            else
                list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_SC)
            end
            #list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_test) ## Temporary test <---------------------------------########
        elseif cumulant == 1
            cond = abs(super_data_M_el) > M_tol
            if cond
                list_stiff = PeriodizeSC.calcintegral_RBZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_cum_AFM_SC)
            else
                list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_cum_SC)
            end
            #list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_cum_AFM_SC)
        end
    elseif modulevec_el.vals_["Periodization"] == 0
        println("Periodization option set to 0")
        #list_stiff = PeriodizeSC.calcintegral_BZ(modelvec, PeriodizeSC.make_stiffness_cluster_G_kintegrand)
        list_stiff = PeriodizeSC.calcintegral_RBZ(modelvec, PeriodizeSC.make_stiffness_cluster_G_kintegrand_fourth)
    end
    return list_stiff
end

"""
Function calculating the superfluid stiffness calling the module SuperStiff.PeriodizeSC in the regime where there is NO coexistence.
Function to use (most suitable) when computing the c-axis superfluid stiffness. Can also compute the a- and b-axis superfluid stiffness, but with
performance drawbacks.

#Argument(s):

- modelvec::PeriodizeSC.ModelVector: ModelVector instance of a particular set of informations read out from a file (PeriodizeSC.ModelVector)
- modulevec_el::Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}}: Data structure to handle the different input parameters extracted from params.json (Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}})

#Returns:

- list_stiff: 2d array containing the values of the superfluid stiffness before summing over Matsubara frequencies (Matrix{Complex{Float64}})
"""
function calc_stiff_funct_NOCOEX(modelvec::PeriodizeSC.ModelVector, modulevec_el::Stiffness.StiffnessArray{Float64,4,Array{Complex{Float64},4},Dict{String,Any}})
    list_stiff = Matrix{Complex{Float64}}(0,0)
    cumulant = modulevec_el.vals_["cumulant"]
    AFM_SC = modulevec_el.vals_["AFM_SC_NOCOEX"]
    println("abc option set to \"c\"")
    if modulevec_el.vals_["Periodization"] == 1
        println("Periodization option set to 1")
        if AFM_SC == 0 && cumulant == 0
            list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_SC)
        elseif AFM_SC == 0 && cumulant == 1
            println("CUM")
            list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_cum_SC)
        elseif AFM_SC == 1 && cumulant == 0
            println("Per AFM_SC_1")
            list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_test)
        elseif AFM_SC == 1 && cumulant == 1
            println("CUM AFM_SC_1")
            list_stiff = PeriodizeSC.calcintegral_BZ(modelvec,PeriodizeSC.make_stiffness_kintegrand_cum_AFM_SC)
        end
    elseif modulevec_el.vals_["Periodization"] == 0
        println("Tr.")
        println("Periodization option set to 0")
        if (AFM_SC == 0 || AFM_SC == 1)
            list_stiff = PeriodizeSC.calcintegral_BZ(modelvec, PeriodizeSC.make_stiffness_trace_G_kintegrand)
        end
    end
    return list_stiff
end

end ## end of module Stiffness scope

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff