text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HeatExchangerDetailed.mso 197 2007-03-08 14:31:57Z bicca $
#*------------------------------------------------------------------
type ShellandTubes_LMTD
	ShellandTubes_LMTD()=begin
		new(
			ShellandTubesBasic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"LMTD Correction Factor Model",
				:Valid=>["Bowmann","Fakeri" , "User Specified"],
				:Default=>"Bowmann"
			)),
			fraction(Dict{Symbol,Any}(
				:Default=>0.8
			)),
			LMTD_Basic(),
			positive(Dict{Symbol,Any}(
				:Brief=>" Capacity Ratio for LMTD Correction Fator",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator ",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator when 2 Pass Shell Side",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakeri Equation",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakeri Equation",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakeri Equation when 2 Pass Shell Side",
				:Lower=>1e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Non - Dimensional Variable for LMTD Correction Fator in Fakeri Equationwhen 2 Pass Shell Side",
				:Lower=>1e-6
			)),
			[
				:(_base_1.Details.Q = _base_1.Details.Ud*_base_1.Pi*_base_1.Tubes.TubeOD*_base_1.Tubes.NumberOfTubes*_base_1.Tubes.TubeLength*Method.LMTD*Method.Fc),
				:(Phi*(2*((_base_1.InletShell.T+ _base_1.OutletShell.T)-(_base_1.InletTube.T+ _base_1.OutletTube.T))) = (sqrt(((_base_1.InletShell.T- _base_1.OutletShell.T)*(_base_1.InletShell.T- _base_1.OutletShell.T))+((_base_1.OutletTube.T - _base_1.InletTube.T)*(_base_1.OutletTube.T - _base_1.InletTube.T))))),
				:(R*(_base_1.OutletTube.T - _base_1.InletTube.T ) = (_base_1.InletShell.T-_base_1.OutletShell.T)),
				:(P*(_base_1.InletShell.T- _base_1.InletTube.T)= (_base_1.OutletTube.T-_base_1.InletTube.T)),
				:(Method.DT0 = _base_1.InletShell.T - _base_1.OutletTube.T),
				:(Method.DTL = _base_1.OutletShell.T - _base_1.InletTube.T),
				:(Phi*(2*((_base_1.InletShell.T+ _base_1.OutletShell.T)-(_base_1.InletTube.T+ _base_1.OutletTube.T))) = (sqrt(((_base_1.InletShell.T- _base_1.OutletShell.T)*(_base_1.InletShell.T- _base_1.OutletShell.T))+((_base_1.OutletTube.T - _base_1.InletTube.T)*(_base_1.OutletTube.T - _base_1.InletTube.T))))),
				:(R*(_base_1.OutletShell.T - _base_1.InletShell.T ) = (_base_1.InletTube.T-_base_1.OutletTube.T)),
				:(P*(_base_1.InletTube.T- _base_1.InletShell.T)= (_base_1.OutletShell.T-_base_1.InletShell.T)),
				:(Method.DT0 = _base_1.InletTube.T - _base_1.OutletShell.T),
				:(Method.DTL = _base_1.OutletTube.T - _base_1.InletShell.T),
				:(Method.Fc = FLMTDcorrection),
				:(lambdaN =1),
				:(lambda1 =1),
				:(Rho =1),
				:(Pc = P),
				:(lambdaN =1),
				:(lambda1 =1),
				:(Rho =1),
				:(Pc*(2-P)= P),
				:(Method.Fc= (sqrt(2)*Pc)/((1-Pc)*ln( abs( ( 2-Pc*0.585786)/( 2-Pc*3.414214))))),
				:(Pc = (sqrt(abs(( 1-P*R)/(1-P)))-1)/(sqrt(abs(( 1-P*R)/(1-P)))-R)),
				:(Method.Fc = sqrt(R*R+1)*ln(abs((1-Pc*R)/(1-Pc)))/((1-R)*ln( abs( ( 2-Pc*(R+1-sqrt(R*R+1)))/ ( 2-Pc*(R + 1 + sqrt(R*R+1))))))),
				:(Pc = P),
				:(Rho*(1-P*R) = (1-P)),
				:(lambdaN = 1),
				:(lambda1 = 1),
				:(Method.Fc = (2*Phi )/(ln(abs((1+Phi )/(1-Phi ))))),
				:(lambdaN = (1/ln(sqrt(abs(Rho))))*((2*sqrt(abs(Rho))-2)/(sqrt(abs(Rho))+1))),
				:(lambda1 = (1/ln(abs(Rho)))*((2*Rho-2)/(Rho+1))),
				:(Method.Fc = ((2*Phi *(lambdaN/lambda1))/(ln(abs((1+Phi *(lambdaN/lambda1))/(1-Phi *(lambdaN/lambda1))))))*(1/lambdaN)),
				:(lambdaN =1),
				:(lambda1 =1),
				:(Pc = P),
				:(Method.Fc = FLMTDcorrection),
				:(Rho =1),
				:(Rho = 1),
				:(Method.Fc = (sqrt(2)*P)/((1-P)*ln( abs( ( 2-P*0.585786)/( 2-P*3.414214))))),
				:(Method.Fc = sqrt(R*R+1)*ln(abs((1-P*R)/(1-P)))/((1-R)*ln( abs( ( 2-P*(R+1-sqrt(R*R+1)))/ ( 2-P*(R + 1 + sqrt(R*R+1))))))),
				:(Rho*(1-P*R) = (1-P)),
				:(Method.Fc = (4*Phi)/(ln(abs((1+2*Phi)/(1-2*Phi))))),
				:(Method.Fc = (2*Phi*(Rho+1)*ln(abs(Rho)))/( ln(abs((1+2*Phi)/(1-2*Phi)))*(Rho-1))),
			],
			[
				"Exchange Surface Area","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation ","R: Capacity Ratio for LMTD Correction Fator","P: Non - Dimensional Variable for LMTD Correction Fator","Temperature Difference at Inlet","Temperature Difference at Outlet","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation ","R: Capacity Ratio for LMTD Correction Fator","P: Non - Dimensional Variable for LMTD Correction Fator","Temperature Difference at Inlet","Temperature Difference at Outlet",""," Variable not in use with this equation"," Variable not in use with this equation"," Variable not in use with this equation"," Variable not in use with this equation"," Variable not in use with Bowmann equation"," Variable not in use with Bowmann equation"," Variable not in use with Bowmann equation","Non Dimensional Variable for LMTD Correction Fator when 2 Pass Shell Side","LMTD Correction Fator when 2 Pass Shell Side","Non Dimensional Variable for LMTD Correction Fator when 2 Pass Shell Side","LMTD Correction Fator when 2 Pass Shell Side"," Variable not in use with Fakeri equation","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation"," Variable not in use when Rho = 1"," Variable not in use when Rho = 1","LMTD Correction Fator when 2 Pass Shell Side","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation","LMTD Correction Fator when 2 Pass Shell Side"," Variable not in use when 1 Pass Shell Side"," Variable not in use when 1 Pass Shell Side"," Variable not in use when 1 Pass Shell Side",""," Variable not in use with this equation"," Variable not in use with Bowmann equation","LMTD Correction Fator when 1 Pass Shell Side","LMTD Correction Fator when 1 Pass Shell Side","Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation","LMTD Correction Fator when 1 Pass Shell Side","LMTD Correction Fator when 1 Pass Shell Side",
			],
			[:LMTDcorrection,:FLMTDcorrection,],
			[:Method,:R,:P,:Pc,:Rho,:Phi,:lambdaN,:lambda1,]
		)
	end
	_base_1::ShellandTubesBasic
	LMTDcorrection::DanaSwitcher
	FLMTDcorrection::fraction
	Method::LMTD_Basic
	R::positive
	P::positive
	Pc::positive
	Rho::positive
	Phi::positive
	lambdaN::positive
	lambda1::positive
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ShellandTubes_LMTD
function setEquationFlow(in::ShellandTubes_LMTD)
	addEquation(1)
	let switch=HotSide
		if HotSide==InletTube.T > InletShell.T
			set(switch,"tubes")
		end
		if HotSide==InletTube.T < InletShell.T
			set(switch,"shell")
		end
		if switch=="shell"
			addEquation(2)
			addEquation(3)
			addEquation(4)
			addEquation(5)
			addEquation(6)
		elseif switch=="tubes"
			addEquation(7)
			addEquation(8)
			addEquation(9)
			addEquation(10)
			addEquation(11)
		end
	end
	let switch=ShellType
		if switch=="Fshell"
			let switch=LMTDcorrection
				if switch=="User Specified"
					#Just For Initialization
					addEquation(12)
					addEquation(13)
					addEquation(14)
					#" Variable not in use with this equation"
					#	Phi = 1;
					addEquation(15)
					addEquation(16)
				elseif switch=="Bowmann"
					addEquation(17)
					addEquation(18)
					#" Variable not in use with Bowmann equation"
					#	Phi = 1;
					addEquation(19)
					if R == 1
						addEquation(20)
						addEquation(21)
					else
						addEquation(22)
						addEquation(23)
					end
				elseif switch=="Fakeri"
					addEquation(24)
					addEquation(25)
					#"Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation "
					#	Phi = (sqrt(((Inlet.Hot.T - Outlet.Hot.T)*(Inlet.Hot.T- Outlet.Hot.T))+((Outlet.Cold.T -  Inlet.Cold.T)*(Outlet.Cold.T -  Inlet.Cold.T))))/(2*((Inlet.Hot.T + Outlet.Hot.T)-( Inlet.Cold.T + Outlet.Cold.T)));
					if Rho == 1
						addEquation(26)
						addEquation(27)
						addEquation(28)
					else
						addEquation(29)
						addEquation(30)
						addEquation(31)
					end
				end
			end
		elseif switch=="Eshell"
			addEquation(32)
			addEquation(33)
			addEquation(34)
			let switch=LMTDcorrection
				if switch=="User Specified"
					#Just For Initialization
					addEquation(35)
					#" Variable not in use with this equation"
					#	lambdaN =1;
					#" Variable not in use with this equation"
					#	lambda1 =1;
					#" Variable not in use with this equation"
					#	Phi = 1;
					addEquation(36)
				elseif switch=="Bowmann"
					#" Variable not in use with Bowmann equation"
					#	Phi  = 1;
					addEquation(37)
					if R == 1
						addEquation(38)
					else
						addEquation(39)
					end
				elseif switch=="Fakeri"
					#"Non Dimensional Variable for LMTD Correction Fator in Fakeri Equation "
					#	Phi  = (sqrt(((Inlet.Hot.T- Outlet.Hot.T)*(Inlet.Hot.T- Outlet.Hot.T))+((Outlet.Cold.T - Inlet.Cold.T)*(Outlet.Cold.T - Inlet.Cold.T))))/(2*((Inlet.Hot.T+ Outlet.Hot.T)-(Inlet.Cold.T+ Outlet.Cold.T)));
					addEquation(40)
					if Rho == 1
						addEquation(41)
					else
						addEquation(42)
					end
				end
			end
		end
	end
end
function atributes(in::ShellandTubes_LMTD,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/STHE"
	fields[:Brief]="Shell and Tubes Heat Exchangers"
	fields[:Info]="to be documented."
	drive!(fields,_)
	return fields
end
ShellandTubes_LMTD(_::Dict{Symbol,Any})=begin
	newModel=ShellandTubes_LMTD()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
include("datasets.jl")
include("attributes_from_dataset.jl")
include("models_cross_validation.jl")
include("graphics.jl")

# este archivo tiene funciones para probar los modelos, y devuelve gráficas con
# la std y media

function topology_test(inputs::Array{Float64,2}, targets::Array{Any,1}, topology::Array{Int64,1}, numFolds::Int64)

    @show(topology);
    learningRate = 0.01;
    numMaxEpochs = 1000;
    validationRatio = 0.2;
    maxEpochsVal = 30;
    numRepetitionsAANTraining = 50;

    classes = unique(targets);
    targets = oneHotEncoding(targets, classes);

    μ_testAcc = Array{Float64,1}(undef, numFolds);
    μ_testF1         = Array{Float64,1}(undef, numFolds);

    crossValidationIndices = crossvalidation(size(inputs,1), numFolds);

    for numFold in 1:numFolds

        local trainingInputs, testInputs, trainingTargets, testTargets;
        trainingInputs    = inputs[crossValidationIndices.!=numFold,:];
        testInputs        = inputs[crossValidationIndices.==numFold,:];
        trainingTargets   = targets[crossValidationIndices.!=numFold,:];
        testTargets       = targets[crossValidationIndices.==numFold,:];

        μ_testAccEachRepetition = Array{Float64,1}(undef, numRepetitionsAANTraining);
        μ_testF1EachRepetition         = Array{Float64,1}(undef, numRepetitionsAANTraining);

        for numTraining in 1:numRepetitionsAANTraining

            local trainingIndices, validationIndices;
            (trainingIndices, validationIndices) = holdOut(size(trainingInputs,1), validationRatio*size(trainingInputs,1)/size(inputs,1));

            local ann;
            (ann, training_acc, _, _, _, _, _) = trainClassANN(topology,
                trainingInputs[trainingIndices,:],   trainingTargets[trainingIndices,:],
                trainingInputs[validationIndices,:], trainingTargets[validationIndices,:],
                testInputs,                          testTargets;
                maxEpochs=numMaxEpochs, learningRate=learningRate, maxEpochsVal=maxEpochsVal);
            #=  Sin validación
            local ann;
            ann, = trainClassANN(topology,
                trainingInputs, trainingTargets,
                testInputs,     testTargets;
                maxEpochs=numMaxEpochs, learningRate=learningRate);
                =#
            (acc, _, _, _, _, _, F1, _) = confusionMatrix(collect(ann(testInputs')'), testTargets);

            μ_testAccEachRepetition[numTraining] = acc;
            μ_testF1EachRepetition[numTraining]         = F1;

        end;
        μ_testAcc[numFold] = mean(μ_testAccEachRepetition);
        μ_testF1[numFold]         = mean(μ_testF1EachRepetition);

    end;
    println("Average test accuracy on a ", numFolds, "-fold crossvalidation: ", 100*mean(μ_testAcc), ", with a standard deviation of ", 100*std(μ_testAcc));
    println("Average test F1 on a ", numFolds, "-fold crossvalidation: ", 100*mean(μ_testF1), ", with a standard deviation of ", 100*std(μ_testF1));
    println();
    return (topology, mean(μ_testAcc), std(μ_testAcc), mean(μ_testF1), std(μ_testF1));
end;

function rna_loop_1(inputs::Array{Float64,2}, targets::Array{Any,1}, nnpplayer::Int64, stlynn::Int64, numFolds::Int64)
    μ_acc = [];
    σ_acc = [];
    μ_f1 = [];
    σ_f1 = [];
    topologyarr = [];
    for i in stlynn:nnpplayer
        (topology, μ_testAcc, σ_stdAcc,
            μ_testF1, σ_stdF1) = topology_test(inputs, targets, [i], numFolds);
        push!(μ_acc,μ_testAcc);
        push!(σ_acc,σ_stdAcc);
        push!(μ_f1,μ_testF1);
        push!(σ_f1,σ_stdF1);
        push!(topologyarr, string("[",i,"]"));
    end;
    return (μ_acc, σ_acc, μ_f1, σ_f1, topologyarr);
end;

function rna_loop_2(inputs::Array{Float64,2}, targets::Array{Any,1}, nnpplayer::Int64, stlynn::Int64, numFolds::Int64)
    μ_acc = [];
    σ_acc = [];
    μ_f1 = [];
    σ_f1 = [];
    topologyarr = [];
    for i in stlynn:nnpplayer
        for j in stlynn:nnpplayer
            (topology, μ_testAcc, σ_stdAcc,
                μ_testF1, σ_stdF1) = topology_test(inputs, targets, [i,j], numFolds);
            push!(μ_acc,μ_testAcc);
            push!(σ_acc,σ_stdAcc);
            push!(μ_f1,μ_testF1);
            push!(σ_f1,σ_stdF1);
            push!(topologyarr, string("[",i,",",j,"]"));
        end;
    end;
    return (μ_acc, σ_acc, μ_f1, σ_f1, topologyarr);
end;

function testRNAs(inputs::Array{Float64,2}, targets::Array{Any,1}, parameters::Dict, numFolds::Int64, rep::Symbol)
    layers = parameters["layers"];
    if (layers==0)
        (_, μ_testAcc, stdμ_testAcc, μ_testF1, σ_stdF1) = topology_test(inputs,targets,[0]);
        println("Test accuracies for a topology of 0 hidden layers: " + string(μ_testAcc));
        println("Standard deviation for test accuracies for a topology of 0 hidden layers: " + string(μ_testAcc));
    else
        if (layers==1)
            (μ_acc, σ_acc, μ_f1, σ_f1, topologyarr) = rna_loop_1(inputs, targets, parameters["maxNNxlayer"], parameters["fstNeuron"], numFolds);
        elseif (layers==2)
            (μ_acc, σ_acc, μ_f1, σ_f1, topologyarr) = rna_loop_2(inputs, targets, parameters["maxNNxlayer"], parameters["fstNeuron"], numFolds);
        else
            println("Test para arquitecturas de 3 capas o superiores no implementados.");
        end;
        printAccStdRNA(μ_acc, σ_acc, μ_f1, σ_f1, topologyarr, rep);
    end;
end;

function testSVM(inputs::Array{Float64,2}, targets::Array{Any,1}, parameters::Dict, numFolds::Int64, rep::Symbol)
    kernel = parameters["kernel"];
    if (kernel=="linear")
        modelHyperparameters = Dict();
        modelHyperparameters["kernel"] = kernel;
        modelHyperparameters["kernelDegree"] = 3;
        modelHyperparameters["kernelGamma"] = 2;
        modelHyperparameters["C"] = 1000000;
        (μ_testAcc, σ_stdAcc, μ_testF1, σ_stdF1) = modelCrossValidation(:SVM, modelHyperparameters, inputs, targets, numFolds);
        println(string("Test accuracies for linear kernel: ",μ_testAcc))
        println(string("Standard deviation for test accuracies of linear kernel: ",σ_stdAcc))

    else
        μ_acc = [];
        σ_acc = [];
        μ_f1 = [];
        σ_f1 = [];
        kernel° = parameters["kernelDegree"];
        max_γ = parameters["maxGamma"];
        plot3d = rep==:Plot3D;
        if (kernel=="poly")
            for degree in 1:kernel°
                if !plot3d
                    local_μ_acc = []
                    local_σ_acc = []
                    local_μ_f1 = []
                    local_σ_f1 = []
                end;
                for γ in 1:max_γ
                    modelHyperparameters = Dict();
                    modelHyperparameters["kernel"] = kernel;
                    modelHyperparameters["kernelDegree"] = degree;
                    modelHyperparameters["kernelGamma"] = γ;
                    modelHyperparameters["C"] = 1000000;
                    (μ_testAcc, σ_stdAcc, μ_testF1, σ_stdF1) = modelCrossValidation(:SVM, modelHyperparameters, inputs, targets, numFolds);
                    if plot3d
                        push!(μ_acc,μ_testAcc);
                        push!(σ_acc,σ_stdAcc);
                        push!(μ_f1,μ_testF1);
                        push!(σ_f1,σ_stdF1);
                    else
                        push!(local_μ_acc,μ_testAcc);
                        push!(local_σ_acc,σ_stdAcc);
                        push!(local_μ_f1,μ_testF1);
                        push!(local_σ_f1,σ_stdF1);
                    end;
                end;
                if !plot3d
                    printAccStd(local_μ_acc, local_σ_acc, local_μ_f1, local_σ_f1, max_γ, "Kernel γ", rep);
                end;
            end;
            if plot3d
                pyplot();
                plot(1:kernel°, 1:max_γ, μ_acc, st=:surface, xlabel = "Kernel°", ylabel = "Kernel γ", camera=(22,30));
            end;

        else
            for γ in 1:max_γ
                modelHyperparameters = Dict();
                modelHyperparameters["kernel"] = "rbf";
                modelHyperparameters["kernelDegree"] = kernel°;
                modelHyperparameters["kernelGamma"] = γ;
                modelHyperparameters["C"] = 1000000;
                (μ_testAcc, σ_stdAcc, μ_testF1, σ_stdF1) = modelCrossValidation(:SVM, modelHyperparameters, inputs, targets, numFolds);
                push!(μ_acc,μ_testAcc);
                push!(σ_acc,σ_stdAcc);
                push!(μ_f1,μ_testF1);
                push!(σ_f1,σ_stdF1);
            end;
            printAccStd(μ_acc, σ_acc, μ_f1, σ_f1, max_γ, "Kernel γ", rep);
        end;
    end;

end;

function testDecisionTree(inputs::Array{Float64,2}, targets::Array{Any,1}, maxDepth::Int64, numFolds::Int64, rep::Symbol)
    μ_acc = [];
    σ_acc = [];
    μ_f1 = [];
    σ_f1 = [];
    for depth in 1:maxDepth
        (μ_testAcc, σ_stdAcc, μ_testF1, σ_stdF1) = modelCrossValidation(:DecisionTree, Dict("maxDepth" => depth), inputs, targets, numFolds);
        push!(μ_acc,μ_testAcc);
        push!(σ_acc,σ_stdAcc);
        push!(μ_f1,μ_testF1);
        push!(σ_f1,σ_stdF1);
    end;
    printAccStd(μ_acc, σ_acc, μ_f1, σ_f1, maxDepth, "maxDepth", rep);
end;

function testKNN(inputs::Array{Float64,2}, targets::Array{Any,1}, max_Neigh::Int64, numFolds::Int64, rep::Symbol)
    μ_acc = [];
    σ_acc = [];
    μ_f1 = [];
    σ_f1 = [];
    for numNeighbors in 1:max_Neigh
        (μ_testAcc, σ_stdAcc, μ_testF1, σ_stdF1) = modelCrossValidation(:kNN, Dict("numNeighbors" => numNeighbors), inputs, targets, numFolds);
        push!(μ_acc,μ_testAcc);
        push!(σ_acc,σ_stdAcc);
        push!(μ_f1,μ_testF1);
        push!(σ_f1,σ_stdF1);
    end;
    printAccStd(μ_acc, σ_acc, μ_f1, σ_f1, max_Neigh, "Number of Neighbors", rep);
end;

function testingModels(modelType::Symbol, parameters::Dict, inputs::Array{Float64,2}, targets::Array{Any,1}, numFolds::Int64; rep::Symbol=:All)
    if modelType==:ANN
        testRNAs(inputs,targets,parameters,numFolds,rep);
    elseif modelType==:SVM
        testSVM(inputs,targets,parameters,numFolds,rep);
    elseif modelType==:DecisionTree
        testDecisionTree(inputs,targets,parameters["maxDepth"],numFolds,rep)
    else
        testKNN(inputs,targets,parameters["maxNeighbors"],numFolds,rep)
    end;
end;

ext: jl


text:
using Base.Iterators: repeated

mutable struct Offsets
    word::Int
    token::Int
end

struct IndexedTokenization{T<:AbstractTokenization} <: WrappedTokenization{T}
    base::T
end
IndexedTokenization() = IndexedTokenization(DefaultTokenization())

_offsets(s, w=0, t=0) = (meta = getmeta(s); hasmeta(s) && haskey(getmeta(s), :offsets) ? meta.offsets : Offsets(w,t))

@inline splitting(p::ParentStages, t::IndexedTokenization, s::TokenStages, x) = zip(repeated(_offsets(s)), Iterators.Flatten((true, Iterators.repeated(false))), splitting(p, t.base, s, x))
@inline splitting(p::ParentStages, t::IndexedTokenization, d::DocumentStage, x) = enumerate(splitting(p, t.base, d, x))

@inline wrap(p::ParentStages, t::IndexedTokenization, s::TokenStages, (i, f, x)) = updatemeta(wrap(p, t.base, s, x), (offsets = i, isfirst = f))
@inline wrap(p::ParentStages, t::IndexedTokenization, d::DocumentStage, (i, x)) = updatemeta(wrap(p, t.base, d, x), (sentence_id = i,))

@inline wrap(p::ParentStages, t::IndexedTokenization, s::TokenStages) = wrap(p, t.base, s)

function wrap(p::ParentStages, t::IndexedTokenization, w::SubWordStage)
    meta = getmeta(w)
    if hasmeta(w) && haskey(meta, :offsets)
        offsets = meta.offsets
        word_id = meta.isfirst ? (offsets.word += 1) : offsets.word
    else
        word_id = 1
    end
    return updatemeta(wrap(p, t.base, w), (word_id = word_id,))
end

function wrap(p::ParentStages, t::IndexedTokenization, x::TokenStage)
    x = wrap(p, t.base, x)
    meta = getmeta(x)
    if hasmeta(x) && haskey(meta, :offsets)
        offsets = meta.offsets
        word_id = haskey(meta, :word_id) ? meta.word_id : (offsets.word += 1)
        token_id = offsets.token += 1
        meta = Base.structdiff(meta, NamedTuple{(:offsets, :isfirst)})
    else
        word_id = token_id = 1
    end
    return setmeta(x, updatemeta(meta, (word_id = word_id, token_id = token_id)))
end


text:
using Dates

# using Dates
# import Franklin
# using Weave
#
# function hfun_bar(vname)
#     val = Meta.parse(vname[1])
#     return round(sqrt(val); digits=2)
# end
#
# function hfun_m1fill(vname)
#     var = vname[1]
#     return pagevar("index", var)
# end
#
# function lx_baz(com, _)
#     # keep this first line
#     brace_content = Franklin.content(com.braces[1]) # input string
#     # do whatever you want here
#     return uppercase(brace_content)
# end
#
# function hfun_eval(arg)
#     x = Core.eval(Franklin, Meta.parse(join(arg)))
#     io = IOBuffer()
#     show(io, "text/plain", x)
#     return String(take!(io))
# end
#

function write_posts(rpaths)::String
    sort_posts!(rpaths)
    curyear = Dates.year(pagevar(rpaths[1], :date))
    io = IOBuffer()
    write(io, "
$curyear
")
    write(io, "
")
    for rp in rpaths
        year = Dates.year(pagevar(rp, :date))
        if year < curyear
            write(io, "
$year
")
            curyear = year
        end
        title = pagevar(rp, :title)
        descr = pagevar(rp, :descr)
        pubdate = Dates.format(Date(pagevar(rp, :date)), "U d")
        path = joinpath(splitpath(rp)[1:2]...)
        write(
            io,
            """

                $pubdate
                $title
                $descr

            """,
        )
    end
    write(io, "
")  #= posts =#
    return String(take!(io))
end

function sort_posts!(rpaths)
    sorter(p) = pagevar(p, :date)
    return sort!(rpaths; by=sorter, rev=true)
end

function hfun_allposts()::String
    jp(p) = joinpath("posts", p)
    rpaths = [
        joinpath(jp(p), "index.md") for
        p in readdir("posts") if isdir(jp(p))
    ]
    return write_posts(rpaths)
end

function hfun_alltags()
    # dictionary {id => Xranklin.Tag}
    # Tag has fields
    #   id (e.g. "the_tag")
    #   name (e.g. "The Tag")
    #   locs (list of rpaths with this tag)
    all_tags   = get_all_tags()
    tag_prefix = getgvar(:tags_prefix, "tag")

    io = IOBuffer()
    write(io, "
")

    for tag_id in keys(all_tags)
        tag_name   = all_tags[tag_id].name
        n_with_tag = length(all_tags[tag_id].locs)
        write(io,
            """

              <nobr>

                  $(tag_name)
                  (
                  $n_with_tag
                  )

              </nobr>

            """,
        )
    end

    write(io, "
")
    return String(take!(io))
end

# NOTE: this is now an internal function; you could overwrite it with your own
# stuff but let's wait a bit for now

# function hfun_taglist()
#     tag = Franklin.locvar(:fd_tag)::String
#     rpaths = Franklin.globvar("fd_tag_pages")[tag]
#     return write_posts(rpaths)
# end

#
# function hfun_weave2html(document)
#     f_name = tempname(pwd()) * ".html"
#     weave(first(document); out_path=f_name)
#     text = read(f_name, String)
#     final = x ->
#         replace(x, r"" => "") |> # Removes weave code block syntax
#         x ->
#             replace(x, "" => "") |> # Removes weave code block syntax
#             x ->
#                 replace(
#                     x,
#                     "
\n" => "
", # Replaces weave code block syntax with Franklin's
#                 ) |> x -> replace(x, "
" => "
")("<!DOCTYPE html>\n" *
#                                                                   split(text, "")[2]) # Replaces weave code block syntax with Franklin's
#     rm(f_name)
#     return final
# end
#

function hfun_posttags()
    tags = get_page_tags()
    base = globvar(:tags_prefix)
    return """
""" *
        join((
            """

              $name

            """
            for (id, name) in tags),
            """•"""
        ) *
        """
"""
end

function hfun_socialicons()
    """

    """
end

#
# """
#     newpost(;title::String, descr::String, tags::Vector{String}, code=false)
# """
# function newpost(;title::String, descr::String, tags::Vector{String}, code=false)
#     path = joinpath(@__DIR__, "posts", replace(lowercase(title), " " => "-"))
#     post = joinpath(path, "index.md")
#     mkpath(path)
#     touch(post)
#     y = Dates.year(Dates.today())
#     m = Dates.month(Dates.today())
#     d = Dates.day(Dates.today())
#     open(post, "w") do io
#         write(io, """
#         +++
#         title = "$title"
#         descr = "$descr"
#         rss = "$descr"
#         date = Date($y, $m, $d)
#         hascode = $code
#         tags = $(sort(tags))
#         +++
#
#         {{ posttags }}
#
#         ## $title
#
#         \\toc
#
#         ### Subtitle
#         """)
#     end
# end

max_forks_repo_name: jvaverka/XranklinTest


text:
using Printf
using Test
using IterableTables
using TableTraits
using IteratorInterfaceExtensions

include( "../src/poverty_inequality.jl" )

function dodecile()
    data = zeros( 100_000, 2 )
    inc = 0.0
    pop = 0.0
    total = 0.0
    for i in 1:100_000
        pop = 1.0
        inc += 1.0
        total += inc
        data[i,1] = pop
        data[i,2] = inc
    end
    print( "total income $total ; population $pop\n")
    deciles = binify( data, 10, 1, 2 )
    print( "Deciles $deciles \n")
    quintiles = binify( data, 20, 1, 2 )
    print( "Quintiles $quintiles \n")
    percentiles = binify( data, 100, 1, 2 )
    print( "Percentiles $percentiles \n")
end

dodecile()


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    varshape(d::Distributions.Distribution)::AbstractValueShape

Get the value shape of the variates of distribution `d`.
"""
varshape(d::Distribution{Univariate}) = ScalarShape{Real}()
varshape(d::Distribution{Multivariate}) = ArrayShape{Real}(size(d)...)
varshape(d::Distribution{Matrixvariate}) = ArrayShape{Real}(size(d)...)

@deprecate valshape(d::Distribution) varshape(d)

"""
    vardof(d::Distributions.Distribution)

Get the number of degrees of freedom of the variates of distribution `d`.
Equivalent to `totalndof(varshape(d))`.
"""
vardof(d::Distribution) = totalndof(varshape(d))

"""
    unshaped(d::Distributions.Distribution)

Turns `d` into a `Distributions.Distribution{Multivariate}` based on
`varshape(d)`.
"""
unshaped(d::Distribution{Multivariate}) = d
unshaped(d::MatrixReshaped) = d.d

const PlainVariate = Union{Univariate,Multivariate,Matrixvariate}

struct StructVariate{T} <: VariateForm end

const NamedTupleVariate{names} = StructVariate{NamedTuple{names}}  # ToDo: Use StructVariate{<:NamedTuple{names}} instead?

Random.rand(rng::AbstractRNG, d::Distribution{NamedTupleVariate{names}}) where names = stripscalar(rand(rng, d, ()))

function Random.rand(rng::AbstractRNG, d::Distribution{NamedTupleVariate{names}}, dims::Tuple{}) where names
    shape = varshape(d)
    x = Vector{default_unshaped_eltype(shape)}(undef, totalndof(varshape(d)))
    rand!(rng, unshaped(d), x)
    shape(x)
end

function Random.rand(rng::AbstractRNG, d::Distribution{NamedTupleVariate{names}}, dims::Dims) where names
    shape = varshape(d)
    X_flat = Array{default_unshaped_eltype(shape)}(undef, totalndof(varshape(d)), dims...)
    X = ArrayOfSimilarVectors(X_flat)
    rand!(rng, unshaped(d), X)
    shape.(X)
end

function Random.rand!(rng::AbstractRNG, d::Distribution{NamedTupleVariate{names}}, x::ShapedAsNT{<:NamedTuple{names}}) where names
    valshape(x) >= varshape(d) || throw(ArgumentError("Shapes of variate and value are not compatible"))
    rand!(rng, unshaped(d), unshaped(x))
    x
end

function _aov_rand_impl!(rng::AbstractRNG, d::Distribution{Multivariate}, X::ArrayOfSimilarVectors{<:Real})
    rand!(rng, unshaped(d), flatview(X))
end

# Workaround for current limitations of ArraysOfArrays.unshaped for standard arrays of vectors
function _aov_rand_impl!(rng::AbstractRNG, d::Distribution{Multivariate}, X::AbstractArray{<:AbstractVector{<:Real}})
    rand!.(Ref(rng), Ref(unshaped(d)), X)
end

function Random.rand!(rng::AbstractRNG, d::Distribution{NamedTupleVariate{names}}, X::ShapedAsNTArray{<:NamedTuple{names}}) where names
    elshape(X) >= varshape(d) || throw(ArgumentError("Shapes of variate and value are not compatible"))
    _aov_rand_impl!(rng, unshaped(d), unshaped.(X))
    X
end

function Distributions.logpdf(d::Distribution{NamedTupleVariate{names}}, x::NamedTuple{names}) where names
    logpdf(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.logpdf(d::Distribution{NamedTupleVariate{names}}, x::AbstractArray{<:NamedTuple{names},0}) where names
    logpdf(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.pdf(d::Distribution{NamedTupleVariate{names}}, x::NamedTuple{names}) where names
    pdf(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.pdf(d::Distribution{NamedTupleVariate{names}}, x::AbstractArray{<:NamedTuple{names},0}) where names
    pdf(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.insupport(d::Distribution{NamedTupleVariate{names}}, x::NamedTuple{names}) where names
    insupport(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.insupport(d::Distribution{NamedTupleVariate{names}}, x::AbstractArray{<:NamedTuple{names},0}) where names
    insupport(unshaped(d), unshaped(x, varshape(d)))
end

function Distributions.insupport(d::Distribution{NamedTupleVariate{names}}, X::AbstractArray{<:NamedTuple{names},N}) where {N,names}
    Distributions.insupport!(BitArray(undef, size(X)), d, X)
end

function Distributions.insupport!(r::AbstractArray{Bool,N}, d::Distribution{NamedTupleVariate{names}}, X::AbstractArray{<:NamedTuple{names},N}) where {N,names}
    r .= insupport.(Ref(d), X)
end


text:
typealias RandVector{T} RandArray{T,1}
typealias RandMatrix{T} RandArray{T,2}

# ## Constructors
# ## ============
RandArray{T,N}(xs::Array{RandVar{T},N}) = RandArray{T,N}(xs)
RandArray(T::DataType, nrows::Int64) =
  RandArray{T,1}(Array(RandVar{T},nrows))
RandArray(T::DataType, nrows::Int64, ncols::Int64) =
  RandArray{T,2}(Array(RandVar{T},nrows,ncols))
RandArray{T<:Real}(a::Array{T}) = RandArray(Any[ConstantRandVar(x) for x in a])

RandArray{T<:SymbolicRandVar}(a::Array{T}) = RandArray(convert(Array{Any},a))

# Fall back when type inference fails
function RandArray(xs::Array{Any,2})
  rtype = rangetype(xs[1,1])
  out = Array(RandVar{rtype},size(xs)...)
  for i = 1:size(xs,1)
    for j = 1:size(xs,2)
      @assert isa(xs[i,j],RandVar)
      @assert rangetype(xs[i,j]) == rtype "inconsistent types of random variables in array"
      out[i,j] = xs[i,j]
    end
  end
  RandArray(out)
end

# Fall back when type inference fails
function RandArray(xs::Array{Any,1})
  rtype = rangetype(xs[1])
  out = Array(RandVar{rtype},size(xs)...)
  for i = 1:size(xs,1)
    @assert isa(xs[i],RandVar)
    @assert rangetype(xs[i]) == rtype "inconsistent types of random variables in array"
    out[i] = xs[i]
  end
  RandArray(out)
end

# ## Conversion
# ## ==========
convert{T}(::Type{RandArray{T,2}}, arr::Matrix{T}) =
  RandArray(RandVar{T}[ConstantRandVar{T}(arr[i,j]) for i=1:size(arr,1), j = 1:size(arr,2)])
convert{T}(::Type{RandArray{T,1}}, arr::Vector{T}) = RandArray(RandVar{T}[ConstantRandVar{T}(a) for a in arr])
convert{T<:RandVar,N}(::Type{RandArray{T,N}}, A::Array{T,N}) = RandArray(A)
convert{T,N}(::Type{Array{RandVar{T},N}}, XS::RandArray{T,N}) = XS.array

convert{T, R<:RandVar}(::Type{Array{RandVar{T}}}, arr::Vector{R}) = RandVar{T}[arr...]
convert{T, R<:RandVar}(::Type{Array{RandVar{T}}}, arr::Matrix{R}) =
  RandVar{T}[arr[i,j] for i=1:size(arr,1), j = 1:size(arr,2)]

call{T}(Xs::RandArray{T,1}, ω::Omega) =
  [call(Xs.array[i],ω) for i = 1:size(Xs.array,1)]
call{T}(Xs::RandArray{T,2}, ω::Omega) =
  [call(Xs.array[i,j],ω) for i = 1:size(Xs.array,1), j = 1:size(Xs.array,2)]

# ## Array Access/Updating
# ## =====================
getindex(Xs::RandVector, i::Int64) = Xs.array[i]
getindex(Xs::RandMatrix, i::Int64, j::Int64) = Xs.array[i,j]
getindex(Xs::RandMatrix, i::Int64) = ((c,r) = divrem(i-1,size(Xs,1)); Xs.array[r+1,c+1])

setindex!(X::RandArray, Y, i::Real) = setindex!(X.array, Y, i)
setindex!(X::RandArray, Y::Sigma.RandArray, i::Real) = setindex!(X.array, Y.array, i)
setindex!(X::RandArray, Y::RandArray, is...) = setindex!(X.array, Y.array, is...)
setindex!(X::RandArray, Y, is...) = setindex!(X.array, Y, is...)

# setindex!{T}(X::RandVector,v::T,i::Int64) = X.array[i] = v
# setindex!{T}(X::RandArray,v::T,i::Int64,j::Int64) = X.array[i,j] = v

function setindex!{T}(X::RandArray{T,2}, Y::RandArray{T,2},
                      u::UnitRange{Int64}, i::Int64)
  setindex!(X.array, Y.array, u, i)
end

rangetype{T,N}(XS::RandArray{T,N}) = Array{T,N}
dims(XS::RandArray) = union(map(dims, XS.array)...)
ndims(XS::RandArray) = length(dims(XS))

size(Xs::RandArray, i::Int) = size(Xs.array, i)
size(Xs::RandArray) = size(Xs.array)
endof(Xs::RandArray) = endof(Xs.array)

getindex(Xs::RandMatrix, i::Int64, js::UnitRange{Int64}) = RandArray(Xs.array[i,js])
getindex(Xs::RandMatrix, is::UnitRange{Int64}, j::Int64) = RandArray(Xs.array[is,j])
getindex(Xs::RandMatrix, is::UnitRange{Int64}, js::UnitRange{Int64}) = RandArray(Xs.array[is,js])
getindex(Xs::RandVector, is::UnitRange{Int64}) = RandArray(Xs.array[is])

function getindex(Xs::RandVector, is::StepRange)
  Ys = RandArray(eltype(Xs),length(is))
  j = 1
  for i in is
    Ys[j] = Xs[i]
    j += 1
  end
  Ys
end

## Concats
hcat(Xs::RandArray, Ys::RandArray) =
  RandArray(hcat(Xs.array, Ys.array))

vcat(Xs::RandArray, Ys::RandArray) =
  RandArray(hcat(Xs.array, Ys.array))

# ## Iteration
start(Xs::RandArray) = start(Xs.array)
next(Xs::RandArray, state) = next(Xs.array, state)
done(Xs::RandArray, state) = done(Xs.array, state)

## Arbitrary Array  Functions
@compat similar{T,N}(X::RandArray{T,N}, elem_type, dims::Tuple{Vararg{Int}}) = RandArray(T,dims...)
# @compat similar{T}(X::RandArray{T,1}, elem_type::Type{RandVar{Float64}}, dims::Tuple{Int64}) = RandArray(T)
# similar(::Sigma.RandArray{Float64,1}, ::Type{Sigma.RandVar{Float64}}, ::Tuple{Int64})

# ## Complex Array Functions
# ## =======================

# import Base:all
"Is every element in Xs true, returns Bool-valued RandVar"
function all{N}(Xs::RandArray{Bool,N})
  x = Xs[1]
  for i = 2:length(Xs)
    x &= Xs[i]
  end
  xs
end

# ## Arithmetic
# ## ==========

# Ambiguity Fixes
(.^){N}(::Base.Irrational{:e}, XS::RandArray{Base.Irrational{:e}, N}) = RandArray((.^)(e,XS.array))
(-){N}(x::Bool, Y::RandArray{Bool,N}) = convert(Int, x) - convert(RandVar{Int}, Y) #FIX ME TESTME: THis prolly doesn't work
(+){N}(x::Bool, Y::RandArray{Bool,N}) = convert(Int, x) + convert(RandVar{Int}, Y) #FIX ME TESTME: THis prolly doesn't work
(+){N}(Y::RandArray{Bool,N}, x::Bool) = convert(RandVar{Int}, Y) + convert(Int, x) #FIX ME TESTME: THis prolly doesn't work

# # Here, we extract the arrays of both args and apply op
for op = (:+, :-, :*, :.*, :/, :&, :|, :.^)
  @eval ($op){T<:Real}(XS::RandArray{T}, YS::RandArray{T}) = RandArray(($op)(XS.array,YS.array))
  # Interop with 'normal arrays' promote them to RandArrays
  @eval ($op){T<:Real,N}(XS::RandArray{T,N}, YS::Array{T,N}) = RandArray(($op)(XS.array,YS))
  @eval ($op){T<:Real,N}(YS::Array{T,N}, XS::RandArray{T,N}) = RandArray(($op)(YS,XS.array))

  # Point wise arithmetic against rand variable (first arg)
  @eval ($op){T<:Real,N}(Y::RandVar{T}, XS::RandArray{T,N}) = RandArray(map(x->($op)(Y,x), XS.array))
  @eval ($op){T<:Real,N}(y::T, XS::RandArray{T,N}) = ($op)(ConstantRandVar(y),XS)

  # Point wise arithmetic against rand variable (second arg)
  @eval ($op){T<:Real,N}(XS::RandArray{T,N}, Y::RandVar{T}) = RandArray(map(x->($op)(x,Y), XS.array))
  @eval ($op){T<:Real,N}(XS::RandArray{T,N}, y::T) = ($op)(XS,ConstantRandVar(y))
end

# Inequalities
for (array_op, elem_op) = ((:.>, :>), (:.>=, :>=), (:.<, :<), (:.<=, :<=), (:(.==), :(==)), (:.!=, :!=))
  @eval ($array_op){T,N}(XS::RandArray{T,N}, YS::RandArray{T,N}) =
    RandArray(map((x,y)->($elem_op)(x,y), XS.array,YS.array))
  # Interop with 'normal arrays' promote them to RandArrays
  @eval ($array_op){T,N}(XS::RandArray{T,N}, ys::Array{T,N}) = RandArray(map((x,y)->($elem_op)(x,y), XS.array,ys))
  @eval ($array_op){T,N}(ys::Array{T,N}, XS::RandArray{T,N}) = RandArray(map((x,y)->($elem_op)(y,x), ys,XS.array))
end

for op = (:(==), :!=)
  eval(
  quote
  function ($op){T,N}(XS::RandArray{T,N}, YS::RandArray{T,N})
    if size(XS) == size(YS)
      all([($op)(XS.array[i],YS.array[i]) for i = 1:length(XS)])
    else
      ConstantRandVar(false)
    end
  end

  function ($op){T,N}(XS::RandArray{T,N}, ys::Array{T,N})
    if size(XS) == size(ys)
      all([($op)(XS.array[i],ys[i]) for i = 1:length(XS)])
    else
      ConstantRandVar(false)
    end
  end

  function ($op){T,N}(ys::Array{T,N}, XS::RandArray{T,N})
    if size(XS) == size(ys)
      all([($op)(XS.array[i],ys[i]) for i = 1:length(XS)])
    else
      ConstantRandVar(false)
    end
  end
  end)
end

# Unary Functions
for op = (:transpose, :ctranspose)
  @eval ($op)(X::RandArray) = RandArray($op(X.array))
end

# Unary ElementWise Functions
for op = (:abs, :asin, :sqrt, :exp, :log, :cos, :sin, :tan, :acos, :asin, :atan,
          :cosh, :sinh, :tanh, :acosh, :asinh, :atanh, :abs, :atan2, :max, :min,
          :sign)
  @eval ($op){T,N}(X::RandArray{T,N}) = RandArray{Float64,N}(map($op,X.array))
end

function lambda{T}(XS::RandArray{T,2})
  X_fns = map(lambda,XS.array)
  ω -> [X_fns[i,j](ω) for i = 1:size(XS,1), j = 1:size(XS,2)]
end

function lambda{T}(XS::RandArray{T,1})
  X_fns = map(lambda,XS.array)
  ω -> [X_fns[i](ω) for i = 1:size(XS,1)]
end

function isapprox{T}(XS::RandArray{T}, ys::Array{T},epsilon = DEFAULT_PREC)
  sum(abs(XS - ys)) <= epsilon
end

function isapprox{T}(ys::Array{T}, XS::RandArray{T}, epsilon = DEFAULT_PREC)
  sum(abs(ys - XS)) <= epsilon
end

function isapprox{T}(XS::RandArray{T}, YS::RandArray{T}; epsilon = DEFAULT_PREC)
  sum(abs(XS - YS)) <= epsilon
end

## Ifelse
ifelse{T,N}(A::RandVar{Bool}, B::Array{T,N}, C::Array{T,N}) =
  RandArray(map((b,c)->ifelse(A,b,c),B,C))

ifelse{T,N}(A::RandVar{Bool}, B::RandArray{T,N}, C::Array{T,N}) =
  RandArray(map((b,c)->ifelse(A,b,c),B,C))

ifelse{T,N}(A::RandVar{Bool}, B::Array{T,N}, C::RandArray{T,N}) =
  RandArray(map((b,c)->ifelse(A,b,c),B,C))

print(io::IO, x::RandArray) = print(io, "what")
show(io::IO, x::RandArray) = show(io, "whats")

# print(io::IO, A::RandArray) = (print("YOU"); print(io, typeof(A),"\n",A.array))
# print(A::RandArray) = (print("ME"); print(A.array))
# println(A::RandArray) = (println("MES"); println(A.array))
# println(io::IO, A::RandArray) = (println("MESA"); println(io, A.array))
# show(io::IO, A::RandArray) = show(io, A.array)
# show(A::RandArray) = show(A.array)


text:
module HCurlTests

using Gridap
using LinearAlgebra
using Test

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition) |> simplexify
order = 0
u((x,y)) = 2*VectorValue(-y,x)
reffe = ReferenceFE(nedelec,order)
V = TestFESpace(model,reffe,dirichlet_tags="boundary")
U = TrialFESpace(V,u)
Ω = Triangulation(model)
dΩ = Measure(Ω,2*(order+1))
a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ
l(v) = ∫(u⋅v)dΩ
op = AffineFEOperator(a,l,U,V)
uh = solve(op)
e = u - uh
el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition) |> simplexify
order = 0
u((x,y,z)) = 2*VectorValue(-y,x,0.) - VectorValue(0.,-z,y)
reffe = ReferenceFE(nedelec,order)
V = TestFESpace(model,reffe,dirichlet_tags="boundary")
U = TrialFESpace(V,u)
Ω = Triangulation(model)
dΩ = Measure(Ω,2*(order+1))
a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ
l(v) = ∫(u⋅v)dΩ
op = AffineFEOperator(a,l,U,V)
uh = solve(op)
e = u - uh
el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)
order = 1
u((x,y,z)) = 2*VectorValue(-y,x,0.) - VectorValue(0.,-z,y)
reffe = ReferenceFE(nedelec,order)
V = TestFESpace(model,reffe,dirichlet_tags="boundary")
U = TrialFESpace(V,u)
Ω = Triangulation(model)
dΩ = Measure(Ω,2*(order+1))
a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ
l(v) = ∫(u⋅v)dΩ
op = AffineFEOperator(a,l,U,V)
uh = solve(op)
e = u - uh
el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
@test el2 < 1.0e-10

end # module


text:
using SinkholeDynamics
using Test

@testset "SinkholeDynamics.jl" begin
    # Write your tests here.
end


text:
try
    ver = readchomp(`protoc --version`)
    @info "ProtoBuf compiler found: $ver"
catch
    @warn "No ProtoBuf compiler present."
    exit()
end

using ProtoBuf

proto_path = abspath(joinpath(@__DIR__, "..", "src", "proto"))
plugin_path = abspath(joinpath(dirname(pathof(ProtoBuf)), "..", "plugin",
    Sys.iswindows() ? "protoc-gen-julia_win.bat" : "protoc-gen-julia"))
julia_out = abspath(joinpath(@__DIR__, "julia_pb"))

if !isdir(julia_out)
    mkdir(julia_out)
end

function genproto(filename)
    run(`protoc --proto_path=$proto_path --julia_out=$julia_out --plugin=protoc-gen-julia=$plugin_path $filename`)
end
@info "Generating proto files"
genproto("API.proto")


text:
using Test, Pkg
@testset "Land" begin
    include("test_heat_parameterizations.jl")
    include("test_water_parameterizations.jl")
    include("prescribed_twice.jl")
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
info("Loading HTTP methods...")

"""
Initiate a websocket|client connection to server defined by url. If the server accepts
the connection and the upgrade to websocket, f is called with an open websocket|client

e.g. say hello, close and leave
```julia
import HTTP
using WebSockets
WebSockets.open("ws://127.0.0.1:8000") do ws
    write(ws, "Hello")
    println("that's it")
end;
```
If a server is listening and accepts, "Hello" is sent (as a Vector{UInt8}).

On exit, a closing handshake is started. If the server is not currently reading
(which is a blocking function), this side will reset the underlying connection (ECONNRESET)
after a reasonable amount of time and continue execution.
"""
function open(f::Function, url; verbose=false, subprotocol = "", kw...)

    key = base64encode(rand(UInt8, 16))
    headers = [
        "Upgrade" => "websocket",
        "Connection" => "Upgrade",
        "Sec-WebSocket-Key" => key,
        "Sec-WebSocket-Version" => "13"
    ]
    if subprotocol != ""
        push!(headers, "Sec-WebSocket-Protocol" => subprotocol )
    end

    if in('#', url)
        throw(ArgumentError(" replace '#' with %23 in url: $url"))
    end

    uri = HTTP.URIs.URI(url)
    if uri.scheme != "ws" && uri.scheme != "wss"
        throw(ArgumentError(" bad argument url: Scheme not ws or wss. Input scheme: $(uri.scheme)"))
    end

    try
        HTTP.open("GET", uri, headers;
                reuse_limit=0, verbose=verbose ? 2 : 0, kw...) do http
                    _openstream(f, http, key)
                end
    catch err
        if typeof(err) <: Base.UVError
            throw(WebSocketClosedError(" while open ws|client: $(string(err))"))
        elseif typeof(err) <: HTTP.ExceptionRequest.StatusError
            return err.response
        else
           rethrow(err)
        end
    end
end
"Called by open with a stream connected to a server, after handshake is initiated"
function _openstream(f::Function, http::HTTP.Streams.Stream, key::String)

    HTTP.startread(http)

    status = http.message.status
    if status != 101
        return
    end

    check_upgrade(http)

    if HTTP.header(http, "Sec-WebSocket-Accept") != generate_websocket_key(key)
        throw(WebSocketError(0, "Invalid Sec-WebSocket-Accept\n" *
                                "$(http.message)"))
    end

    io = HTTP.ConnectionPool.getrawstream(http)
    ws = WebSocket(io,false)
    try
        f(ws)
    finally
        close(ws)
    end

end

"""
Used as part of a server definition. Call this if
is_upgrade(http.message) returns true.

Responds to a WebSocket handshake request.
If the connection is acceptable, sends status code 101
and headers according to RFC 6455, then calls
user's handler function f with the connection wrapped in
a WebSocket instance.

f(ws)           is called with the websocket and no client info
f(headers, ws)  also receives a dictionary of request headers for added security measures

On exit from f, a closing handshake is started. If the client is not currently reading
(which is a blocking function), this side will reset the underlying connection (ECONNRESET)
after a reasonable amount of time and continue execution.

If the upgrade is not accepted, responds to client with '400'.

e.g. server with local error handling. Combine with WebSocket.open example.
```julia
import HTTP
using WebSockets

badgatekeeper(reqdict, ws) = sqrt(-2)
handlerequest(req) = HTTP.Response(501)

try
    HTTP.listen("127.0.0.1", UInt16(8000)) do http
        if WebSockets.is_upgrade(http.message)
            WebSockets.upgrade(badgatekeeper, http)
        else
            HTTP.Servers.handle_request(handlerequest, http)
        end
    end
catch err
    showerror(STDERR, err)
    println.(catch_stacktrace()[1:4])
end
```
"""
function upgrade(f::Function, http::HTTP.Stream)
    # Double check the request. is_upgrade should already have been called by user.
    check_upgrade(http)
    if !HTTP.hasheader(http, "Sec-WebSocket-Version", "13")
        HTTP.setheader(http, "Sec-WebSocket-Version" => "13")
        HTTP.setstatus(http, 400)
        HTTP.startwrite(http)
        return
    end
    if HTTP.hasheader(http, "Sec-WebSocket-Protocol")
        requestedprotocol = HTTP.header(http, "Sec-WebSocket-Protocol")
        if !hasprotocol(requestedprotocol)
            HTTP.setheader(http, "Sec-WebSocket-Protocol" => requestedprotocol)
            HTTP.setstatus(http, 400)
            HTTP.startwrite(http)
            return
        else
            HTTP.setheader(http, "Sec-WebSocket-Protocol" => requestedprotocol)
        end
    end
    key = HTTP.header(http, "Sec-WebSocket-Key")
    if length(base64decode(key)) != 16 # Key must be 16 bytes
        HTTP.setstatus(http, 400)
        HTTP.startwrite(http)
        return
    end
    # This upgrade is acceptable. Send the response.
    HTTP.setheader(http, "Sec-WebSocket-Accept" => generate_websocket_key(key))
    HTTP.setheader(http, "Upgrade" => "websocket")
    HTTP.setheader(http, "Connection" => "Upgrade")
    HTTP.setstatus(http, 101)
    HTTP.startwrite(http)
    # Pass the connection on as a WebSocket.
    io = HTTP.ConnectionPool.getrawstream(http)
    ws = WebSocket(io, true)
    # If the callback function f has two methods,
    # prefer the more secure one which takes (request, websocket)
    try
        if applicable(f, http.message, ws)
            f(http.message, ws)
        else
            f(ws)
        end
#    catch err
#        warn("WebSockets.HTTP.upgrade: Caught unhandled error while calling argument function f, the handler / gatekeeper:\n\t")
#        mt = typeof(f).name.mt
#        fnam = splitdir(string(mt.defs.func.file))[2]
#        print_with_color(:yellow, STDERR, "f = ", string(f) * " at " * fnam * ":" * string(mt.defs.func.line) * "\nERROR:\t")
#        showerror(STDERR, err, catch_stacktrace())
    finally
        close(ws)
    end
end

"It is up to the user to call 'is_upgrade' on received messages.
This provides double checking from within the 'upgrade' function."
function check_upgrade(http)
    if !HTTP.hasheader(http, "Upgrade", "websocket")
        throw(WebSocketError(0, "Check upgrade: Expected \"Upgrade => websocket\"!\n$(http.message)"))
    end
    if !(HTTP.hasheader(http, "Connection", "upgrade") || HTTP.hasheader(http, "Connection", "keep-alive, upgrade"))
        throw(WebSocketError(0, "Check upgrade: Expected \"Connection => upgrade or Connection => keep alive, upgrade\"!\n$(http.message)"))
    end
end

"""
Fast checking for websockets vs http requests, performed on all new HTTP requests.
Similar to HttpServer.is_websocket_handshake
"""
function is_upgrade(r::HTTP.Message)
    if (r isa HTTP.Request && r.method == "GET")  || (r isa HTTP.Response && r.status == 101)
        if HTTP.header(r, "Connection", "") != "keep-alive"
            # "Connection => upgrade" for most and "Connection => keep-alive, upgrade" for Firefox.
            if HTTP.hasheader(r, "Connection", "upgrade") || HTTP.hasheader(r, "Connection", "keep-alive, upgrade")
                if lowercase(HTTP.header(r, "Upgrade", "")) == "websocket"
                    return true
                end
            end
        end
    end
    return false
end
# Inline docs in 'WebSockets.jl'
target(req::HTTP.Messages.Request) = req.target
subprotocol(req::HTTP.Messages.Request) = HTTP.header(req, "Sec-WebSocket-Protocol")
origin(req::HTTP.Messages.Request) = HTTP.header(req, "Origin")

"""
WebsocketHandler(f::Function) <: HTTP.Handler

A simple Function-wrapper for HTTP.
The provided argument should be one of the forms
    `f(WebSocket) => nothing`
    `f(HTTP.Request, WebSocket) => nothing`
The latter form is intended for gatekeeping, ref. RFC 6455 section 10.1

f accepts a `WebSocket` and does interesting things with it, like reading, writing and exiting when finished.

Take note of the very similar WebSocketHandler (capital 'S'), which is a subtype of HttpServer, an alternative
to HTTP.
"""
struct WebsocketHandler{F <: Function} <: HTTP.Handler
    func::F # func(ws) or func(request, ws)
end

"""
    WebSockets.ServerWS(::HTTP.HandlerFunction, ::WebSockets.WebsocketHandler)

WebSockets.ServerWS is an argument type for WebSockets.serve. Instances
include .in  and .out channels, see WebSockets.serve.
"""
mutable struct ServerWS{T <: HTTP.Servers.Scheme, H <: HTTP.Handler, W <: WebsocketHandler}
    handler::H
    wshandler::W
    logger::IO
    in::Channel{Any}
    out::Channel{Any}
    options::HTTP.ServerOptions

    ServerWS{T, H, W}(handler::H, wshandler::W, logger::IO = HTTP.compat_stdout(), ch=Channel(1), ch2=Channel(2),
                 options=HTTP.ServerOptions()) where {T, H, W} =
        new{T, H, W}(handler, wshandler, logger, ch, ch2, options)
end

ServerWS(h::Function, w::Function, l::IO=HTTP.compat_stdout();
            cert::String="", key::String="", args...) = ServerWS(HTTP.HandlerFunction(h), WebsocketHandler(w), l;
                                                         cert=cert, key=key, args...)
function ServerWS(handler::H,
                wshandler::W,
                logger::IO = HTTP.compat_stdout();
                cert::String = "",
                key::String = "",
                args...) where {H <: HTTP.Handler, W <: WebsocketHandler}
    if cert != "" && key != ""
        serverws = ServerWS{HTTP.Servers.https, H, W}(handler, wshandler, logger, Channel(1), Channel(2), HTTP.ServerOptions(; sslconfig=HTTP.MbedTLS.SSLConfig(cert, key), args...))
    else
        serverws = ServerWS{HTTP.Servers.http, H, W}(handler, wshandler, logger, Channel(1), Channel(2), HTTP.ServerOptions(; args...))
    end
    return serverws
end

"""
    WebSockets.serve(server::ServerWS, port)
    WebSockets.serve(server::ServerWS, host, port)
    WebSockets.serve(server::ServerWS, host, port, verbose)

A wrapper for HTTP.listen.
Puts any caught error and stacktrace on the server.out channel.
To stop a running server, put HTTP.Servers.KILL on the .in channel.
```julia
    @shedule WebSockets.serve(server, "127.0.0.1", 8080)
```
After a suspected connection task failure:
```julia
    if isready(server.out)
        err = take!(myserver_WS.out)
    end
    if isready(myserver_WS.out)
        stack_trace = take!(server_WS.out)
    end
```
"""
function serve(server::ServerWS{T, H, W}, host, port, verbose) where {T, H, W}

    tcpserver = Ref{HTTP.Sockets.TCPServer}()

    @async begin
        while !isassigned(tcpserver)
            sleep(1)
        end
        while true
            val = take!(server.in)
            val == HTTP.Servers.KILL && close(tcpserver[])
        end
    end

    HTTP.listen(host, port;
           tcpref=tcpserver,
           ssl=(T == HTTP.Servers.https),
           sslconfig = server.options.sslconfig,
           verbose = verbose,
           tcpisvalid = server.options.ratelimit > 0 ? HTTP.Servers.check_rate_limit :
                                                     (tcp; kw...) -> true,
           ratelimits = Dict{IPAddr, HTTP.Servers.RateLimit}(),
           ratelimit = server.options.ratelimit) do stream::HTTP.Stream
                            try
                                if is_upgrade(stream.message)
                                    upgrade(server.wshandler.func, stream)
                                else
                                    HTTP.Servers.handle_request(server.handler.func, stream)
                                end
                            catch err
                                put!(server.out, err)
                                put!(server.out, catch_stacktrace())
                            end
            end
    return
end
serve(server::WebSockets.ServerWS, host, port) =  serve(server, host, port, false)
serve(server::WebSockets.ServerWS, port) =  serve(server, "127.0.0.1", port, false)


text:
module PixelArt

include("functions.jl")

export pixel

end


text:
ceil(x::FastRational{T}) where {T<:Integer} = FastRational{T}(numerator(ceil(T, x.num//x.den)), one(T))
floor(x::FastRational{T}) where {T<:Integer} = FastRational{T}(numerator(floor(T, x.num//x.den)), one(T))
trunc(x::FastRational{T}) where {T<:Integer} = FastRational{T}(numerator(trunc(T, x.num//x.den)), one(T))
round(x::FastRational{T}) where {T<:Integer} = FastRational{T}(numerator(round(T, x.num//x.den)), one(T))

ceil(::Type{I}, x::FastRational{T}) where {I<:Integer, T<:Integer} = ceil(I, x.num//x.den)
floor(::Type{I}, x::FastRational{T}) where {I<:Integer, T<:Integer} = floor(I, x.num//x.den)
trunc(::Type{I}, x::FastRational{T}) where {I<:Integer, T<:Integer} = trunc(I, x.num//x.den)
round(::Type{I}, x::FastRational{T}) where {I<:Integer, T<:Integer} = round(I, x.num//x.den)

ceil(::Type{Integer}, x::FastRational{T}) where {T<:Integer} = ceil(Integer, x.num//x.den)
floor(::Type{Integer}, x::FastRational{T}) where {T<:Integer} = floor(Integer, x.num//x.den)
trunc(::Type{Integer}, x::FastRational{T}) where {T<:Integer} = trunc(Integer, x.num//x.den)
round(::Type{Integer}, x::FastRational{T}) where {T<:Integer} = round(Integer, x.num//x.den)

round(x::FastRational{T}, ::RoundingMode{:ToZero}) where {T<:Integer} = signbit(x) ? ceil(x) : floor(x)
round(x::FastRational{T}, ::RoundingMode{:FromZero}) where {T<:Integer} = signbit(x) ? floor(x) : ceil(x)
round(x::FastRational{T}, ::RoundingMode{:Up}) where { T<:Integer} = FastRational{T}(ceil(T, x.num//x.den), one(T))
round(x::FastRational{T}, ::RoundingMode{:Down}) where {T<:Integer} = FastRational{T}(floor(T, x.num//x.den), one(T))
round(x::FastRational{T}, ::RoundingMode{:Nearest}) where {T<:Integer} = FastRational{T}(round(T, x.num//x.den), one(T))

round(::Type{I}, x::FastRational{T}, ::RoundingMode{:ToZero}) where {I<:Integer, T<:Integer} = signbit(x) ? ceil(I,x) : floor(I,x)
round(::Type{I}, x::FastRational{T}, ::RoundingMode{:FromZero}) where {I<:Integer, T<:Integer} = signbit(x) ? floor(I,x) : ceil(I,x)
round(::Type{I}, x::FastRational{T}, ::RoundingMode{:Up}) where {I<:Integer, T<:Integer} = ceil(I, x.num//x.den)
round(::Type{I}, x::FastRational{T}, ::RoundingMode{:Down}) where {I<:Integer, T<:Integer} = floor(I, x.num//x.den)
round(::Type{I}, x::FastRational{T}, ::RoundingMode{:Nearest}) where {I<:Integer, T<:Integer} = round(I, x.num//x.den)

round(::Type{Integer}, x::FastRational{T}, ::RoundingMode{:ToZero}) where {T<:Integer} = signbit(x) ? ceil(T,x) : floor(T,x)
round(::Type{Integer}, x::FastRational{T}, ::RoundingMode{:FromZero}) where {T<:Integer} = signbit(x) ? floor(T,x) : ceil(T,x)
round(::Type{Integer}, x::FastRational{T}, ::RoundingMode{:Up}) where {T<:Integer} = ceil(T, x.num//x.den)
round(::Type{Integer}, x::FastRational{T}, ::RoundingMode{:Down}) where {T<:Integer} = floor(T, x.num//x.den)
round(::Type{Integer}, x::FastRational{T}, ::RoundingMode{:Nearest}) where {T<:Integer} = round(T, x.num//x.den)


text:
import GeometryTypes

struct Mesh{N, T, PT <: Polyhedron{T}} <: GeometryTypes.GeometryPrimitive{N, T}
    polyhedron::PT
end
function Mesh{N}(polyhedron::Polyhedron{T}) where {N, T}
    return Mesh{N, T, typeof(polyhedron)}(polyhedron)
end
function Mesh(polyhedron::Polyhedron, ::StaticArrays.Size{N}) where N
    return Mesh{N[1]}(polyhedron)
end
function Mesh(polyhedron::Polyhedron, N::Int)
    # This is type unstable but there is no way around that,
    # use polyhedron built from StaticArrays vector to avoid that.
    return Mesh{N}(polyhedron)
end
function Mesh(polyhedron::Polyhedron)
    return Mesh(polyhedron, FullDim(polyhedron))
end

# Creates a scene for the vizualisation to be used to truncate the lines and rays
function scene(vr::VRep, ::Type{T}) where T
    # First compute the smallest rectangle containing the P-representation (i.e. the points).
    (xmin, xmax) = extrema(map((x)->x[1], points(vr)))
    (ymin, ymax) = extrema(map((x)->x[2], points(vr)))
    (zmin, zmax) = extrema(map((x)->x[3], points(vr)))
    width = max(xmax-xmin, ymax-ymin, zmax-zmin)
    if width == zero(T)
        width = 2
    end
    scene = GeometryTypes.HyperRectangle{3, T}([(xmin + xmax) / 2 - width,
                                                (ymin + ymax) / 2 - width,
                                                (zmin + zmax) / 2 - width],
                                               2 * width * ones(T, 3))
    # Intersection of rays with the limits of the scene
    (start, ray) -> begin
        times = max.((Vector(minimum(scene))-start) ./ ray, (Vector(maximum(scene))-start) ./ ray)
        times[ray .== 0] .= Inf # To avoid -Inf with .../(-0)
        time = minimum(times)
        start + time * ray
    end
end

function _isdup(zray, triangles)
    for tri in triangles
        normal = tri[2]
        if isapproxzero(cross(zray, normal)) && dot(zray, normal) > 0 # If A[j,:] is almost 0, it is always true...
            # parallel and equality or inequality and same sense
            return true
        end
    end
    false
end
_isdup(poly, hidx, triangles) = _isdup(get(poly, hidx).a, triangles)

function fulldecompose(poly_geom::Mesh{3}, ::Type{T}) where T
    poly = poly_geom.polyhedron
    exit_point = scene(poly, T)

    triangles = Tuple{Tuple{Vector{T},Vector{T},Vector{T}}, Vector{T}}[]

    function decomposeplane(hidx)
        h = get(poly, hidx)
        # xray should be the rightmost ray
        xray = nothing
        # xray should be the leftmost ray
        yray = nothing
        zray = h.a
        isapproxzero(zray) && return

        # Checking rays
        counterclockwise(a, b) = dot(cross(a, b), zray)
        line = nothing
        lineleft = false
        lineright = false
        function checkleftright(r::Union{Ray, Line})
            cc = counterclockwise(r, line)
            if !isapproxzero(cc)
                if cc < 0 || islin(r)
                    lineleft = true
                end
                if cc > 0 || islin(r)
                    lineright = true
                end
            end
        end
        for l in incidentlines(poly, hidx)
            if !isapproxzero(l)
                if line === nothing
                    line = l
                else
                    checkleftright(l)
                end
            end
        end
        for r in incidentrays(poly, hidx)
            if !isapproxzero(r)
                if line === nothing
                    if xray === nothing || counterclockwise(r, xray) > 0
                        xray = coord(r) # r is more right than xray
                    end
                    if yray === nothing || counterclockwise(r, yray) < 0
                        yray = coord(r) # r is more left than xray
                    end
                else
                    checkleftright(r)
                end
            end
        end

        # Checking vertices
        face_vert = pointtype(poly)[]
        for x in points(poly)
            if _isapprox(dot(x, zray), h.β)
                push!(face_vert, x)
            end
        end

        if line !== nothing
            if isempty(face_vert)
                center = origin(pointtype(poly), 3)
            else
                center = first(face_vert)
            end
            hull = pointtype(poly)[]
            push!(hull, exit_point(center, line))
            if lineleft
                push!(hull, exit_point(center, cross(zray, line)))
            end
            push!(hull, exit_point(center, -line))
            if lineright
                push!(hull, exit_point(center, cross(line, zray)))
            end
            hulls = (hull,)
        else
            #if length(face_vert) < 3 # Wrong, they are also the rays
            #  error("Not enough vertices and rays to form a face, it may be because of numerical rounding. Otherwise, please report this bug.")
            #end
            if length(face_vert) < 3 && (xray == nothing || (length(face_vert) < 2 && (yray == xray || length(face_vert) < 1)))
                return
            end
            if xray == nothing
                sweep_norm = cross(zray, [1,0,0])
                if sum(abs, sweep_norm) == 0
                    sweep_norm = cross(zray, [0,1,0])
                end
            else
                sweep_norm = cross(zray, xray)
            end
            sort!(face_vert, by = x -> dot(x, sweep_norm))
            xtoy_hull = getsemihull(face_vert, 1, counterclockwise, yray)
            if yray == nothing
                ytox_hull = getsemihull(face_vert, -1, counterclockwise, yray)
            else
                ytox_hull = Any[]
                push!(ytox_hull, face_vert[1])
                if last(xtoy_hull) != face_vert[1]
                    push!(ytox_hull, last(xtoy_hull))
                end
                push!(ytox_hull, exit_point(last(xtoy_hull), yray))
                push!(ytox_hull, exit_point(face_vert[1], xray))
            end
            hulls = (xtoy_hull, ytox_hull)
        end
        for hull in hulls
            if length(hull) >= 3
                a = pop!(hull)
                b = pop!(hull)
                while !isempty(hull)
                    c = pop!(hull)
                    push!(triangles, ((a,b,c), zray))
                    b = c
                end
            end
        end
    end

    for hidx in eachindex(hyperplanes(poly))
        decomposeplane(hidx)
    end
    # If there is already a triangle, his normal is an hyperplane and it is the only face
    if isempty(triangles)
        for hidx in eachindex(halfspaces(poly))
            if !_isdup(poly, hidx, triangles)
                decomposeplane(hidx)
            end
        end
    end

    ntri = length(triangles)
    pts = Vector{GeometryTypes.Point{3, T}}(undef, 3ntri)
    faces = Vector{GeometryTypes.Face{3, Int}}(undef, ntri)
    ns = Vector{GeometryTypes.Normal{3, T}}(undef, 3ntri)
    for i in 1:ntri
        tri = pop!(triangles)
        normal = tri[2]
        for j = 1:3
            idx = 3*(i-1)+j
            #ns[idx] = -normal
            ns[idx] = normal
        end
        faces[i] = collect(3*(i-1) .+ (1:3))
        k = 1
        for k = 1:3
            # reverse order of the 3 vertices so that if I compute the
            # normals with the `normals` function, they are in the good
            # sense.
            # I know I don't use the `normals` function but I don't know
            # what is the OpenGL convention so I don't know if it cares
            # about the order of the vertices.
            pts[3*i-k+1] = tri[1][k]
        end
    end
    # If the type of ns is Rational, it also works.
    # The normalized array in in float but then it it recast into Rational
    map!(normalize, ns, ns)
    (pts, faces, ns)
end

fulldecompose(poly::Mesh{N, T}) where {N, T} = fulldecompose(poly, typeof(one(T)/2))

GeometryTypes.isdecomposable(::Type{<:GeometryTypes.Point{3}},  ::Type{<:Mesh{3}}) = true
GeometryTypes.isdecomposable(::Type{<:GeometryTypes.Face{3}},   ::Type{<:Mesh{3}}) = true
GeometryTypes.isdecomposable(::Type{<:GeometryTypes.Normal{3}}, ::Type{<:Mesh{3}}) = true
function GeometryTypes.decompose(PT::Type{<:GeometryTypes.Point}, poly::Mesh)
    points = fulldecompose(poly)[1]
    map(PT, points)
end
function GeometryTypes.decompose(FT::Type{<:GeometryTypes.Face}, poly::Mesh)
    faces = fulldecompose(poly)[2]
    GeometryTypes.decompose(FT, faces)
end
function GeometryTypes.decompose(NT::Type{<:GeometryTypes.Normal}, poly::Mesh)
    ns = fulldecompose(poly)[3]
    map(NT, ns)
end

# In AbstractPlotting, when asking to plot an object, it calls this constructor
# which is only defined for `GeometryTypes.GeometryPrimitive` which is a
# supertype of `Polyhedra.Mesh`
GeometryTypes.GLNormalMesh(p::Polyhedron) = GeometryTypes.GLNormalMesh(Mesh(p))


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

## Definitions of special intervals and associated functions

## Empty interval:
"""`emptyinterval`s are represented as the interval [∞, -∞]; note
that this interval is an exception to the fact that the lower bound is
larger than the upper one."""
emptyinterval(::Type{T}) where T<:Real = Interval{T}(Inf, -Inf)
emptyinterval(x::Interval{T}) where T<:Real = emptyinterval(T)
emptyinterval() = emptyinterval(precision(Interval)[1])
emptyinterval(::Type{<:Integer}) = emptyinterval(Float64)
const ∅ = emptyinterval(Float64)

isempty(x::Interval) = x.lo == Inf && x.hi == -Inf

const ∞ = Inf

## Entire interval:
"""`entireinterval`s represent the whole Real line: [-∞, ∞]."""
entireinterval(::Type{T}) where T<:Real = Interval{T}(-Inf, Inf)
entireinterval(x::Interval{T}) where T<:Real = entireinterval(T)
entireinterval() = entireinterval(precision(Interval)[1])
entireinterval(::Type{<:Integer}) = entireinterval(Float64)

isentire(x::Interval) = x.lo == -Inf && x.hi == Inf
isinf(x::Interval) = isentire(x)
isunbounded(x::Interval) = x.lo == -Inf || x.hi == Inf

isfinite(x::Interval) = isfinite(x.lo) && isfinite(x.hi)

"""
    isthin(x)

Checks if `x` is the set consisting of a single exactly
representable float. Any float which is not exactly representable
does *not* yield a thin interval. Corresponds to `isSingleton` of
the standard.
"""
isthin(x::Interval) = x.lo == x.hi

"""
    iscommon(x)

Checks if `x` is a **common interval**, i.e. a non-empty,
bounded, real interval.
"""
function iscommon(a::Interval)
    (isentire(a) || isempty(a) || isnai(a) || isunbounded(a)) && return false
    true
end

"`widen(x)` widens the lowest and highest bounds of `x` to the previous and next representable floating-point numbers, respectively."
widen(x::Interval{T}) where {T<:AbstractFloat} = Interval(prevfloat(x.lo), nextfloat(x.hi))

"""
    wideinterval(x::AbstractFloat)

Returns the interval Interval( prevfloat(x), nextfloat(x) ).
"""
wideinterval(x::T) where {T<:AbstractFloat} = Interval( prevfloat(x), nextfloat(x) )

"""
    isatomic(x::Interval)

Check whether an interval `x` is *atomic*, i.e. is unable to be split.
This occurs when the interval is empty, or when the upper bound equals the lower bound or the `nextfloat` of the lower bound.
"""
isatomic(x::Interval) = isempty(x) || (x.hi == x.lo) || (x.hi == nextfloat(x.lo))

iszero(x::Interval) = iszero(x.lo) && iszero(x.hi)

contains_zero(X::Interval{T}) where {T} = zero(T) ∈ X

# """
#     widen(x)
#
# Widens the lowest and highest bounds of `x` to the previous
# and next representable floating-point numbers, respectively.
# """
# widen(x::Interval{T}) where T<:AbstractFloat =
#     Interval(prevfloat(x.lo), nextfloat(x.hi))


text:
using LinearAlgebra
using MultivariateSeries

X = @ring x1 x2

n = length(X)
d = 4
r = 3

w0  = rand(r)
Xi0 = rand(Float64, n,r)
println("w0=",w0, "   Xi0=",Xi0)
s0 = series(w0, Xi0, X, d)

w, Xi = decompose(s0)
println("w=",w, "     Xi=",Xi)

eps = 1.e-1
w  += eps*rand(Float64, length(w))
Xi += eps*rand(Float64, size(Xi))
println("w=",w, "     Xi=",Xi)

w, Xi = vdm_newton(w, Xi, s0, monoms(X,d), eps=1.e-10, maxit=25)
println("w=",w, "     Xi=",Xi)

s = series(w, Xi, X, d)
println("Error: ", norm(s-s0,Inf))


text:
export Line, PeriodicLine

## Standard interval

# angle is π*a where a is (false==0) and (true==1)
# or ranges from (-1,1].  We use 1 as 1==true.

"""
    Line{a}(c)

represents the line at angle `a` in the complex plane, centred at `c`.
"""
struct Line{angle,T<:Number} <: SegmentDomain{T}
    center::T
    α::Float64
    β::Float64

    #TODO get this inner constructor working again.
    Line{angle,T}(c,α,β) where {angle,T} = new{angle,T}(c,α,β)
    Line{angle,T}(c) where {angle,T} = new{angle,T}(c,-1.,-1.)
    Line{angle,T}() where {angle,T} = new{angle,T}(zero(T),-1.,-1.)
end

const RealLine{T} = Union{Line{false,T},Line{true,T}}

Line{a}(c,α,β) where {a} = Line{a,typeof(c)}(c,α,β)
Line{a}(c::Number) where {a} = Line{a,typeof(c)}(c)
Line{a}() where {a} = Line{a,Float64}()

angle(d::Line{a}) where {a} = a*π

reverseorientation(d::Line{true}) = Line{false}(d.center,d.β,d.α)
reverseorientation(d::Line{false}) = Line{true}(d.center,d.β,d.α)
reverseorientation(d::Line{a}) where {a} = Line{a-1}(d.center,d.β,d.α)

# ensure the angle is always in (-1,1]
Line(c,a,α,β) = Line{mod(a/π-1,-2)+1,typeof(c)}(c,α,β)
Line(c,a) = Line(c,a,-1.,-1.)
# true is negative orientation, false is positive orientation
# this is because false==0 and we take angle=0
Line(b::Bool) = Line{b}()
Line() = Line(false)

isambiguous(d::Line)=isnan(d.center)
convert(::Type{Domain{T}},d::Line{a}) where {a,T<:Number} = Line{a,T}(d.center, d.α, d.β)
convert(::Type{Line{a,T}},::AnyDomain) where {a,T<:Number} = Line{a,T}(NaN)
convert(::Type{IT},::AnyDomain) where {IT<:Line}=Line(NaN,NaN)

## Map interval

##TODO non-1 alpha,beta

isempty(::Line) = false

function line_tocanonical(α,β,x)
    @assert α==β==-1. || α==β==-.5

    if α==β==-1.
        2x/(1+sqrt(1+4x^2))
    elseif α==β==-.5
        x/sqrt(1 + x^2)
    end
end

function line_tocanonicalD(α,β,x)
    @assert α==β==-1. || α==β==-.5

    if α==β==-1.
        2/(1+4x^2+sqrt(1+4x^2))
    elseif α==β==-0.5
        (1 + x^2)^(-3/2)
    end
end
function line_fromcanonical(α,β,x)
    #TODO: why is this consistent?
    if α==β==-1.
        x/(1-x^2)
    else
        x*(1 + x)^α*(1 - x)^β
    end
end
function line_fromcanonicalD(α,β,x)
    if α==β==-1.
        (1+x^2)/(1-x^2)^2
    else
        (1 - (β-α)x - (β+α+1)x^2)*(1+x)^(α-1)*(1-x)^(β-1)
    end
end

function line_invfromcanonicalD(α,β,x)
    if α==β==-1.
        (1-x^2)^2/(1+x^2)
    else
        1/(1 - (β-α)x - (β+α+1)x^2)*(1+x)^(1-α)*(1-x)^(1-β)
    end
end

tocanonical(d::Line,x) = line_tocanonical(d.α,d.β,cis(-angle(d)).*(x-d.center))
tocanonical(d::Line{false},x) = line_tocanonical(d.α,d.β,x-d.center)
tocanonical(d::Line{true},x) = line_tocanonical(d.α,d.β,d.center-x)

tocanonicalD(d::Line,x) = cis(-angle(d)).*line_tocanonicalD(d.α,d.β,cis(-angle(d)).*(x-d.center))
tocanonicalD(d::Line{false},x) = line_tocanonicalD(d.α,d.β,x-d.center)
tocanonicalD(d::Line{true},x) = -line_tocanonicalD(d.α,d.β,d.center-x)

fromcanonical(d::Line,x) = cis(angle(d))*line_fromcanonical(d.α,d.β,x)+d.center
fromcanonical(d::Line{false},x) = line_fromcanonical(d.α,d.β,x)+d.center
fromcanonical(d::Line{true},x) = -line_fromcanonical(d.α,d.β,x)+d.center

fromcanonicalD(d::Line,x) = cis(angle(d))*line_fromcanonicalD(d.α,d.β,x)
fromcanonicalD(d::Line{false},x) = line_fromcanonicalD(d.α,d.β,x)
fromcanonicalD(d::Line{true},x) = -line_fromcanonicalD(d.α,d.β,x)

invfromcanonicalD(d::Line,x) = cis(-angle(d))*line_invfromcanonicalD(d.α,d.β,x)
invfromcanonicalD(d::Line{false},x) = line_invfromcanonicalD(d.α,d.β,x)
invfromcanonicalD(d::Line{true},x) = -line_invfromcanonicalD(d.α,d.β,x)

==(d::Line{a},m::Line{a}) where {a} = d.center == m.center && d.β == m.β &&d.α == m.α

# algebra
*(c::Real,d::Line{false}) = Line{sign(c)>0 ? false : true}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Real,d::Line{true}) = Line{sign(c)>0 ? true : false}(isapprox(d.center,0) ? d.center : c*d.center,d.α,d.β)
*(c::Number,d::Line) = Line(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c),d.α,d.β)
*(d::Line,c::Number) = c*d
for OP in (:+,:-)
    @eval begin
        $OP(c::Number,d::Line{a}) where {a} = Line{a}($OP(c,d.center),d.α,d.β)
        $OP(d::Line{a},c::Number) where {a} = Line{a}($OP(d.center,c),d.α,d.β)
    end
end

## Periodic line

# angle is (false==0) and π (true==1)
# or ranges from (-1,1]
struct PeriodicLine{angle,T} <: PeriodicDomain{T}
    center::T
    L::Float64
    PeriodicLine{angle,T}(c,L) where {angle,T} = new{angle,T}(c,L)
    PeriodicLine{angle,T}(c) where {angle,T} = new{angle,T}(c,1.)
    PeriodicLine{angle,T}(d::PeriodicLine) where {angle,T} = new{angle,T}(d.center,d.L)
    PeriodicLine{angle,T}() where {angle,T} = new{angle,T}(0.,1.)
end

(::Type{PeriodicLine{a}})(c,L) where {a} = PeriodicLine{a,typeof(c)}(c,L)

PeriodicLine(c,a) = PeriodicLine{a/π,eltype(c)}(c,1.)
PeriodicLine() = PeriodicLine{false,Float64}(0.,1.)
PeriodicLine(b::Bool) = PeriodicLine{b,Float64}()

isambiguous(d::PeriodicLine) = isnan(d.center) && isnan(d.angle)
convert(::Type{Domain{T}},d::PeriodicLine{a}) where {a,T<:Number} = PeriodicLine{a,T}(d.center, d.L)
convert(::Type{PeriodicLine{T,TT}},::AnyDomain) where {T<:Number,TT} = PeriodicLine{T,TT}(NaN,NaN)
convert(::Type{IT},::AnyDomain) where {IT<:PeriodicLine} = PeriodicLine(NaN,NaN)

angle(d::PeriodicLine{a}) where {a}=a*π

reverseorientation(d::PeriodicLine{true})=PeriodicLine{false}(d.center,d.L)
reverseorientation(d::PeriodicLine{false})=PeriodicLine{true}(d.center,d.L)
reverseorientation(d::PeriodicLine{a}) where {a}=PeriodicLine{a-1}(d.center,d.L)

tocanonical(d::PeriodicLine{false},x) = real(2atan((x-d.center)/d.L))
fromcanonical(d::PeriodicLine{false},θ) = d.L*tan(θ/2) + d.center

tocanonical(d::PeriodicLine{a},x) where {a} = tocanonical(PeriodicLine{false,Float64}(0.,d.L),exp(-π*im*a)*(x-d.center))
fromcanonical(d::PeriodicLine{a},v::AbstractArray) where {a} =
    [fromcanonical(d,vk) for vk in v]
fromcanonical(d::PeriodicLine{a},x) where {a} =
    exp(π*im*a)*fromcanonical(PeriodicLine{false,Float64}(0.,d.L),x)+d.center

function invfromcanonicalD(d::PeriodicLine{false})
    @assert d.center==0  && d.L==1.0
    a=Fun(PeriodicSegment(),[1.,0,1])
end

mappoint(a::PeriodicLine{false},b::Circle,x) = b.radius*((a.L*im-(x-a.center))./(a.L*im+(x-a.center)))+b.center
function mappoint(b::Circle,a::PeriodicLine{false},x)
    y=(x-b.center)./b.radius
    a.center+a.L*im*(1-y)./(y+1)
end

# algebra
*(c::Number,d::PeriodicLine)=PeriodicLine(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c))
*(d::PeriodicLine,c::Number)=c*d
for OP in (:+,:-)
    @eval begin
        $OP(c::Number,d::PeriodicLine{a,T}) where {a,T}=PeriodicLine{a,promote_type(eltype(c),T)}($OP(c,d.center),d.L)
        $OP(d::PeriodicLine{a,T},c::Number) where {a,T}=PeriodicLine{a,promote_type(eltype(c),T)}($OP(d.center,c),d.L)
    end
end

arclength(d::Union{Line,PeriodicLine}) = Inf
first(d::Union{Line,PeriodicLine})= -Inf
last(d::Union{Line,PeriodicLine})= Inf
complexlength(d::Union{Line,PeriodicLine})=Inf

## vectorized

for typ in (:Line,:PeriodicLine)
    @eval function convert(::Type{$typ},d::ClosedInterval)
        a,b=d.left,d.right
        @assert abs(a) == abs(b) == Inf

        if isa(a,Real) && isa(b,Real)
            if a==Inf
                @assert b==-Inf
                $typ(true)
            else
                @assert a==-Inf&&b==Inf
                $typ(false)
            end
        elseif abs(real(a)) < Inf
            @assert real(a)==real(b)
            @assert sign(imag(a))==-sign(imag(b))

            $typ(real(b),angle(b))
        elseif isnan(real(a)) && isnan(real(b))  # hack for -im*Inf
            $typ([imag(a)*im,imag(b)*im])
        elseif isnan(real(a))  # hack for -im*Inf
            $typ([real(b)+imag(a)*im,b])
        elseif isnan(real(b))  # hack for -im*Inf
            $typ([a,real(a)+imag(b)*im])
        elseif abs(imag(a)) < Inf
            @assert imag(a)==imag(b)
            @assert sign(real(a))==-sign(real(b))

            $typ(imag(b),angle(b))
        else
            @assert angle(b) == -angle(a)

            $typ(0.,angle(b))
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool <: SwaggerModel

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool(;)
        o = new()
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool = Dict{Symbol,Symbol}()
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrBool }, name::Symbol, val)
end


text:
"""
equivalent to [0; cumsum(x)]
"""
function zero_cumsum(x::Array{T, 1}) where T
    res = Array{T, 1}(undef, length(x) + 1)

    res[1] = zero(T)
    for i in 1:length(x)
        res[i + 1] = res[i] + x[i]
    end

    return res
end

function zero_cumsum(f::Function, x::Array{T, 1}) where T
    res = Array{T, 1}(undef, length(x) + 1)

    res[1] = zero(T)
    for i in 1:length(x)
        res[i + 1] = res[i] + f(x[i])
    end

    return res
end

"""
# Description
Creates a segment cost function assuming times series data has Normal distribution with changing mean
and variance. The returned function takes two indices and calculates twice the negative log-likelihood
this model between these indices.

# Arguments
* `data::Array{Float64}`: Time series

# Returns
* `cost::Function`: Function which takes indices (s, t) where
                    0 ≤ s < t < n and n is the length of the time series
                    and returns the cost of the segment [s+1, ..., t]

# Example
```julia
n = 1000
λ = 100
μ, σ = Normal(0.0, 10.0), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)
cost = NormalMeanSegment(data, σ = 1.0)
cps, cost = BS(cost, n)
```

# See also
NormalVarSegment, NormalMeanVarSegment, ExponentialSegment, PoissonSegment, GammaRateSegment, GammaShapeSegment, Nonparametric, PELT, BS
"""
function NormalMeanSegment(data::Array{Float64}, σ::Real = 1.0)
    cd = zero_cumsum(data)
    cd_2 = zero_cumsum(x -> x^2, data)
    cost(s::Int64, t::Int64) =
        ( cd_2[t+1] - cd_2[s+1] - abs2(cd[t+1] - cd[s+1]) / (t-s) ) / (σ^2)
    return cost
end

"""
# Description
Creates a segment cost function assuming times series data has Normal distribution with changing variance
and fixed mean.

# See also
NormalMeanSegment
"""
function NormalVarSegment(data::Array{Float64}, μ::Real)
    ss = zero_cumsum(x -> (x - μ)^2, data)
    cost(s::Int64, t::Int64) = (t-s) * log( (ss[t+1] - ss[s+1])/(t-s) )
end

"""
# Description
Creates a segment cost function assuming times series data has Normal distribution with changing mean and variance.

# See also
NormalMeanSegment
"""
function NormalMeanVarSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cd_2 = zero_cumsum(x -> x^2, data)
    function cost(s::Int64, t::Int64)
      mu = (cd[t + 1] - cd[s + 1]) / (t - s)
      sig = ( cd_2[t+1] - cd_2[s+1] ) / (t-s) - mu^2
      return (t-s)*log(sig) + (cd_2[t+1] - cd_2[s+1])/sig - 2*(cd[t+1] - cd[s+1])*mu/sig + ((t-s)*mu^2)/sig
    end
    return cost
end

"""
# Description
Creates a segment cost function assuming times series data has Exponential distribution with changing mean

# See also
NormalMeanSegment
"""
function ExponentialSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cost(s::Int64, t::Int64) = -*(t-s) * ( log(t-s) - log(cd[t+1] - cd[s+1]))
    return cost
end

"""
# Description
Creates a segment cost function assuming times series data has Poisson distribution with changing mean

# See also
NormalMeanSegment
"""
function PoissonSegment(data::Array{Float64})
    cd = zero_cumsum(data)
    cost(s::Int64, t::Int64) = -2*(cd[t+1]-cd[s+1]) * ( log(cd[t+1]-cd[s+1]) - log(t-s) - 1 )
    return cost
end

"""
# Description
Creates a segment cost function assuming times series data has Gamma distribution with changing shape parameter

# See also
NormalMeanSegment
"""
function GammaShapeSegment(data::Array{Float64}, beta::Float64)
    lcd = zero_cumsum(log, data)
    function cost(s::Int64, t::Int64)
        alpha_hat = invdigamma( log(beta) + (lcd[t+1] - lcd[s+1])/(t-s) )
        cost = (t-s)*( alpha_hat*log(beta) - lgamma(alpha_hat) ) + (alpha_hat-1)*(lcd[t+1]-lcd[s+1])
        return cost*-2
    end
end

"""
# Description
Creates a segment cost function assuming times series data has Gamma distribution with changing rate parameter

# See also
NormalMeanSegment
"""
function GammaRateSegment(data::Array{Float64}, alpha::Float64)
    cd = zero_cumsum(data)
    function cost(s::Int64, t::Int64)
        beta_hat = ( alpha * (t-s) )/( cd[t+1] - cd[s+1] )
        cost = (t-s)*alpha*log(beta_hat) - beta_hat*( cd[t+1] - cd[s+1] )
        return cost*-2
    end
end

"""
# Description
Creates a segment cost function based on the Nonparametric model

# Arguments
`data::Array{Float64}`: time series
`K::Int`: nonparametric model parameter

# See also
NormalMeanSegment
"""
function NonparametricSegment(data::Array{Float64}, K::Int64)
    n = length(data)
    if (K > n)
        K = n
    end
    x = sort(data)
    yk = -1 .+ (2 .* collect(1:K) ./ K .- 1/K)
    c = -log(2*n-1)
    pK = (1 .+ exp.(c .* yk)) .^ -1.0
    j = floor(Int64, (n-1)*pK[1]+1)
    Q = vcat(0,cumsum(data .< x[j]) .+ 0.5 .* cumsum(data.==x[j]))
    for i in 2:K
        j = convert(Int64,floor((n-1)*pK[i]+1))
        Q =vcat(Q,(vcat(0,cumsum(data .< x[j]) + 0.5*cumsum(data.==x[j]))))
    end

    function cost(s::Int64, t::Int64)
        Fkl = zeros(Float64, K)
        for i in 1:K
            Fkl[i] = ((Q[(t+1)+(i-1)*(n+1)]) -  (Q[(s+1)+(i-1)*(n+1)]))/(t-s)
        end
        cost = 0
        for i in 1:K
            if (Fkl[i] > 0 && Fkl[i] < 1)
                cost  = cost + (t-s)*(Fkl[i]*log(Fkl[i]) + (1-Fkl[i])*log(1-Fkl[i]))
            end
        end
        return 2*c*cost/K
    end
end

"""
# Description
Create a segment cost function for piecewise linear regressions, fitted using OLS (assuming Normally distributed errors).

# Arguments
`data::Array{Float64}`: time series

# See also
NormalMeanSegment
"""
function OLSSegment(data::Array{Float64})
    function cost(s::Int64, t::Int64)
	      if t-s > 1
            y = data[s:t]
	          m = length(y)
	          x = 1:m
	          a = ( 2*(2*m+1)*sum(y) - 6*sum(x.*y) )/( m*(m-1) )
	          b = ( 12*sum(x.*y) - 6*(m+1)*sum(y) )/( m*(m-1)*(m+1) )
	          sig = sum( (y-a-b*x).^2 )/m
	          return m/2 * ( log(sig) + 1 )
	      else
	          return Inf
	      end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PolicyAssignmentProperties <: SwaggerModel
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    policyDefinitionId::Any # spec type: Union{ Nothing, String } # spec name: policyDefinitionId
    scope::Any # spec type: Union{ Nothing, String } # spec name: scope

    function PolicyAssignmentProperties(;displayName=nothing, policyDefinitionId=nothing, scope=nothing)
        o = new()
        validate_property(PolicyAssignmentProperties, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(PolicyAssignmentProperties, Symbol("policyDefinitionId"), policyDefinitionId)
        setfield!(o, Symbol("policyDefinitionId"), policyDefinitionId)
        validate_property(PolicyAssignmentProperties, Symbol("scope"), scope)
        setfield!(o, Symbol("scope"), scope)
        o
    end
end # type PolicyAssignmentProperties

const _property_map_PolicyAssignmentProperties = Dict{Symbol,Symbol}(Symbol("displayName")=>Symbol("displayName"), Symbol("policyDefinitionId")=>Symbol("policyDefinitionId"), Symbol("scope")=>Symbol("scope"))
const _property_types_PolicyAssignmentProperties = Dict{Symbol,String}(Symbol("displayName")=>"String", Symbol("policyDefinitionId")=>"String", Symbol("scope")=>"String")
Base.propertynames(::Type{ PolicyAssignmentProperties }) = collect(keys(_property_map_PolicyAssignmentProperties))
Swagger.property_type(::Type{ PolicyAssignmentProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PolicyAssignmentProperties[name]))}
Swagger.field_name(::Type{ PolicyAssignmentProperties }, property_name::Symbol) =  _property_map_PolicyAssignmentProperties[property_name]

function check_required(o::PolicyAssignmentProperties)
    true
end

function validate_property(::Type{ PolicyAssignmentProperties }, name::Symbol, val)
end


text:
using Plots
import Distributions
using Random
using Turing
using Statistics
using StatsPlots
using FillArrays

d1 = Normal(1,1)
d2 = Normal(2,1)

sample1 = rand(d1,20)
sample2 = rand(d2,20)

histogram(sample1)
histogram!(sample2)

@model function bayesgroups(sample1,sample2)
    m ~ Normal(0,10)
    σ ~ Uniform(0.001,5)
    μ ~ Uniform(-10,10)
    sample1 ~ MvNormal(fill(μ+m/2,length(sample1)),σ)
    sample2 ~ MvNormal(fill(μ-m/2,length(sample2)),σ)
end

model1 = bayesgroups(sample1,sample2)
chain = sample(model1,NUTS(0.65),1000)

plot(chain)


text:
# This file was generated, do not modify it. # hide
ŷ = MLJ.predict(mtm, Xtrain)
cross_entropy(ŷ, ytrain) |> mean

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
Base.@kwdef struct AdjointSolver{T, LS<:AbstractLinearSolver{T}, FS<:FittingSet{T}}
    prob::NonlinearProblem{T}
    targ::TargetFunction{T}
    linalg::LS
    fset::FS
    μ::Vector{T}
    gp::Vector{T}
    g::Vector{T}
end

function AdjointSolver{T}(prob::NonlinearProblem{T}, targ::TargetFunction{T}, linalg::AbstractLinearSolver{T}, fset::FittingSet{T}) where {T}
    Nt = size(prob.C, 2)
    Nx = length(fset.cache.xbuf)
    kwargs = (
        μ = Array{T}(undef, Nt),
        gp = Array{T}(undef, Nt),
        g = Array{T}(undef, Nx),
    )
    AdjointSolver{T, typeof(linalg), typeof(fset)}(; prob, targ, linalg, fset, kwargs...)
end

function solve!(solver::AdjointSolver{T}; verbose=false) where {T}

    @unpack μ, g, gp, prob, linalg, targ, fset = solver

    fill!(μ, zero(T))
    fill!(g, zero(T))
    @inbounds for n = length(prob.pviews):-1:2
        @unpack jac_next = prob.pviews[n]
        jac = linalg.Fbuf[n]

        # Adjoint vector calculation
        grad!(gp, targ, n)
        @turbo for i = 1:length(gp)
            gp[i] = -(jac_next[i] * μ[i] + gp[i])
        end
        # TODO: Somehow faster than 'ldiv!'
        μ .= jac \ gp

        # Gradient calculation with respect to parameters
        grad!(g, fset, prob, targ, μ, n)

        verbose && println("n: $n, mu: $μ")
    end

    # Gradient of the L2 term of the objective function
    grad!(g, targ.terms.L2)
    # Gradient correction according to parameter scaling
    unscaleg!(g, fset.scale)

    return solver
end


text:
function plot_glc_limiting_study(simid, Ds, ϵs, cg, is_glclim_stst)

    ## -----------------------------------------------------
    # simulation
    Dyn_ugs, Dyn_Xs = [], []
    net_Uugs, net_UXs = [], []
    ϵs_ = []

    @showprogress for (ϵi, ϵ) in enumerate(ϵs)

        for D in Ds

            ## ------------------------------------------------------
            # Dynamic
            simparams = (;D, ϵ, cg)
            status = get_status(simid, simparams)
            if (status != STST_SIM_STATUS)
                @warn("Wrong status ", (D, ϵ, cg), status)
                continue
            end
            datfile = Dyn.simdat_file(simid, simparams)
            if !isfile(datfile)
                @warn("File missing ", (D, ϵ, cg), datfile)
                continue
            end
            S = ldat(datfile)
            dyn_ug = S.ug_av
            dyn_X = S.X

            ## ------------------------------------------------------
            # LP
            net = S.V.net
            Dyn.fix!(net, :z, D)
            Dyn.ub!(net, :ug, cg * D / S.X)
            Dyn.lb!(net, :ug, 0.0)
            Dyn.reducebox!(net)
            Lug, Uug = Dyn.fva(net, :ug)
            UX = cg * D / Lug
            if !is_glclim_stst[simparams]
                @warn("Not glc-lim stst ", (D, ϵ, cg))
                continue
            end

            push!(Dyn_ugs, dyn_ug)
            push!(Dyn_Xs, dyn_X)
            push!(net_Uugs, Uug)
            push!(net_UXs, UX)

            push!(ϵs_, ϵ)

        end

    end

    isempty(ϵs_) && error("Not data recorded!!!")
    max_ϵs_ = maximum(ϵs)
    ms0 = 5

    ps = Plots.Plot[]
    for (p, dyn_vs, net_vs) in [
            (plot(;xlabel = _textbf("Dynamic (u_g)"), ylabel = _textbf("c_{g}D/X")), Dyn_ugs, net_Uugs),
            (plot(;xlabel = _textbf("Dynamic (X)"), ylabel = _textbf("Theo. maximum (X)")), Dyn_Xs, net_UXs)
        ]

        all_ = [dyn_vs; net_vs] |> sort!
        lims = (last(all_) - first(all_))
        lims = (first(all_) - (lims * 0.1), last(all_) + (lims * 0.1))
        plot!(p, all_, all_;
            ls = :dash, lw = 3, c = :black, alpha = 0.7, label = "", xlim = lims, ylim = lims
        )

        scatter!(p, dyn_vs, net_vs;
            label = "",
            markerstrokecolor = :white,
            color = :black,
            ms = ms0 .+ (2 * ms0) .* (ϵs_ ./ max_ϵs_),
        )
        push!(ps, p)
    end

    return ps

end


text:
using Test

include("nfa.jl")

M = NFA(Dict(
             ("q0", '0') => ["q0", "q1"],
             ("q0", '1') => ["q0", "q2"],
             ("q1", '0') => ["q3",],
             ("q2", '0') => ["q2", "q3"],
             ("q2", '1') => ["q3",],
             ("q3", '0') => ["q3",],
             ("q3", '1') => ["q3",]
       ), "q0", Set(["q3"]))

@test !evaluate!(M, "")
@test !evaluate!(M, "01")
@test evaluate!(M, "00")
@test evaluate!(M, "1111110000")

# Taken from here:
# https://www.cs.wcupa.edu/rkline/assets/img/FCS/nfas/n2d2.jpg?1379001651
# regex: (a* + (ab)*) b*
M = NFA(Dict(
             ("1",  λ)  => ["2", "3"],
             ("2",  λ)  => ["5"],
             ("2", 'a') => ["2"],
             ("5", 'b') => ["5"],
             ("3", 'a') => ["4"],
             ("3",  λ)  => ["5"],
             ("4", 'b') => ["3"],
       ), "1", Set(["5"]))

@test evaluate!(M, "")
@test evaluate!(M, "b")
@test evaluate!(M, "a")
@test evaluate!(M, "ab")
@test evaluate!(M, "aa")
@test evaluate!(M, "abab")
@test !evaluate!(M, "aba")

max_issues_repo_name: MadhavJivrajani/DFA-Sim

max_forks_repo_name: MadhavJivrajani/DFA-Sim


text:
#' MANUALLY Start `jl`

using Pkg
cd("/home/kar/project/notebook/reports/02__stats/stats_rethinking/StatsRethinking")
Pkg.activate(".")

using RCall
@rlibrary ggplot2
@rlibrary ggdist
@rlibrary theme
reval("suppressMessages(library(theme))");

R"packageVersion('theme')"
@rlibrary utils
packageVersion("theme")

ggplot() + stat_halfeye(aes(randn(100)))


text:
include("preamble.jl")

# -----------------------
# Figures with VTK
# -----------------------

using LazySets, WriteVTK
import Polyhedra

# propagate sets
rot(θ) = [cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1]
H = Hyperrectangle(ones(3), fill(0.1, 3))
adom = range(0, 2pi, length=20)
X = [rot(θ) * H for θ in adom]

# convexify
R = [CH(X[i], X[i+1]) for i in 1:length(X)-1]
push!(R, CH(X[end], X[1]))
RCH = concretize.(R)

# save
writevtk(X, file="rotations")
writevtk(RCH, file="rotations_CH")

using ReachabilityAnalysis

@taylorize function lorenz!(dx, x, p, t)
    local σ = 10.0
    local β = 8.0 / 3.0
    local ρ = 28.0
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
    return dx
end

X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])
prob = @ivp(x' = lorenz!(x), dim=3, x(0) ∈ X0);
alg = TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000);

sol = solve(prob, T=20.0, alg=alg);

solZ = set.(overapproximate(sol, Zonotope))
writevtk(solZ, file="lorenz") # generates lorenz.vtu


text:
using BlsData
using Base.Test

b0 = Bls()
b1 = Bls("test")

data = get_data(b0, "PRS85006092")
data = get_data(b0, ["PRS85006092", "LNS11000000"])
data = get_data(b0, "PRS85006092"; startyear=1990)
data = get_data(b0, ["PRS85006092", "LNS11000000"]; endyear=2014)
data = get_data(b0, ["PRS85006092", "LNS11000000"]; startyear=1980, endyear=2014)
data = get_data(b0, "PRS85006092"; catalog=true)


text:
using SafeTestsets, Test

@safetestset "CanonicalSymplecticForms" begin include("test_CanonicalSymplecticForms.jl") end

@safetestset "Plan Unit Tests" begin
    @safetestset "FirstFinDiffPlans" begin include("test_FirstFinDiffPlans.jl") end
    @safetestset "FirstFourierPlans" begin include("test_FirstFourierPlans.jl") end

    @safetestset "SecondChebyshevPlans" begin include("test_SecondChebyshevPlans.jl") end
    @safetestset "SecondFinDiffPlans" begin include("test_SecondFinDiffPlans.jl") end
end

@safetestset "PoincareInvariants" begin
    include("test_PoincareInvariants.jl")
end


text:
# indicator of the L1 norm ball with given radius

export IndBallL1

"""
**Indicator of a ``L_1`` norm ball**

    IndBallL1(r=1.0)

Returns the indicator function of the set
```math
S = \\left\\{ x : ∑_i |x_i| \\leq r \\right\\}.
```
Parameter `r` must be positive.
"""
struct IndBallL1{R <: Real} <: ProximableFunction
  r::R
  function IndBallL1{R}(r::R) where {R <: Real}
    if r <= 0
      error("parameter r must be positive")
    else
      new(r)
    end
  end
end

is_convex(f::IndBallL1) = true
is_set(f::IndBallL1) = true

IndBallL1(r::R=1.0) where {R <: Real} = IndBallL1{R}(r)

function (f::IndBallL1)(x::AbstractArray{T}) where T <: RealOrComplex
  if norm(x,1) - f.r > 1e-12
    return +Inf
  end
  return zero(T)
end

function prox!(y::AbstractArray{T}, f::IndBallL1, x::AbstractArray{T}, gamma::R=one(R)) where {R<: Real, T <: RealOrComplex{R}}
  # TODO: a faster algorithm
  if norm(x,1) - f.r < 1e-14
    y .= x
    return zero(T)
  else # do a projection of abs(x) onto simplex then recover signs
    n = length(x)
    p = abs.(view(x,:))
    sort!(p, rev=true)
    s = zero(R)
    @inbounds for i = 1:n-1
      s = s + p[i]
      tmax = (s - f.r)/i
      if tmax >= p[i+1]
        @inbounds for j in eachindex(x)
          y[j] = sign(x[j])*max(abs(x[j])-tmax, zero(R))
        end
        return zero(T)
      end
    end
    tmax = (s + p[n] - f.r)/n
    @inbounds for j in eachindex(x)
      y[j] = sign(x[j])*max(abs(x[j])-tmax, zero(R))
    end
    return zero(T)
  end
end

fun_name(f::IndBallL1) = "indicator of an L1 norm ball"
fun_dom(f::IndBallL1) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndBallL1) = "x ↦ 0 if ‖x‖_1 ⩽ r, +∞ otherwise"
fun_params(f::IndBallL1) = "r = $(f.r)"

function prox_naive(f::IndBallL1, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  # do a simple bisection (aka binary search) on λ
  L = 0.0
  U = maximum(abs, x)
  λ = L
  v = 0.0
  maxit = 120
  for iter in 1:maxit
    λ = 0.5*(L + U)
    v = sum(max.(abs.(x) .- λ, 0.0))
    # modify lower or upper bound
    (v < f.r) ? U = λ : L = λ
    # exit condition
    if abs(L - U) < 1e-15
      break
    end
  end
  return sign.(x) .* max.(0.0, abs.(x) .- λ), 0.0
end


text:
module day1

readInput() = parse.(Int, readlines(joinpath(@__DIR__, "input.txt")))

function part1(input)
    for i = 1:length(input), j = i+1:length(input)
        a, b = input[i], input[j]
        a + b == 2020 && return a * b
    end
end

function part2(input)
    for i = 1:length(input), j = i+1:length(input), k = j+1:length(input)
        a, b, c = input[i], input[j], input[k]
        a + b + c == 2020 && return a * b * c
    end
end

end