text:
"""
    parseGPExperimentInfo(GPExperimentInfoDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseCGPInfo(GPInfoDict::Dict, experiment::GeneticJuliaStructure)
    nodesFile = get(GPInfoDict, "nodesFile", "")
    if !(typeof(nodesFile) <: String)
        error("nodesFile field must be a string with the path of the file containing the information about the nodes")
    end

    maxTreeDepth = get(GPInfoDict, "maxTreeDepth", false)
    if experiment._experimentInfo._individualType <: GAGenotype
        if maxTreeDepth != false
            @warn "Maximum tree depth set in a non Genetic Programming problem, maximum depth field will be dismissed"
        end
    elseif experiment._experimentInfo._individualType <: GPGenotype
        if maxTreeDepth != false
            if !(typeof(maxTreeDepth) <: Integer)
                error("Maximum tree depth must be an integer number greater than 0")
            end
        else
            @warn "Maximum tree depth not set in a Genetic Programming problem, maximum depth will be set to 3 by default"
            maxTreeDepth = 3
        end
    end

    setCGPInfo(nodesFile=nodesFile, maxTreeDepth=maxTreeDepth, genj=experiment)
end

"""
    parseGPExperimentInfo(GPExperimentInfoDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseSTGPInfo(GPInfoDict::Dict, experiment::GeneticJuliaStructure)
    nodesFile = get(GPInfoDict, "nodesFile", "")
    if !(typeof(nodesFile) <: String)
        error("nodesFile field must be a string with the path of the file containing the information about the nodes")
    end

    maxTreeDepth = get(GPInfoDict, "maxTreeDepth", false)
    if experiment._experimentInfo._individualType <: GAGenotype
        if maxTreeDepth != false
            @warn "Maximum tree depth set in a non Genetic Programming problem, maximum depth field will be dismissed"
        end
    elseif experiment._experimentInfo._individualType <: GPGenotype
        if maxTreeDepth != false
            if !(typeof(maxTreeDepth) <: Integer)
                error("Maximum tree depth must be an integer number greater than 0")
            end
        else
            @warn "Maximum tree depth not set in a Genetic Programming problem, maximum depth will be set to 3 by default"
            maxTreeDepth = 3
        end
    end

    setSTGPInfo(nodesFile=nodesFile, maxTreeDepth=maxTreeDepth, genj=experiment)
end

"""
    parseGPExperimentInfo(GPExperimentInfoDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseGEPInfo(GEPInfoDict::Dict, experiment::GeneticJuliaStructure)
    nodesFile = get(GEPInfoDict, "nodesFile", "")
    if !(typeof(nodesFile) <: String)
        error("nodesFile field must be a string with the path of the file containing the information about the nodes")
    end

    head = get(GEPInfoDict, "head", false)
    if experiment._experimentInfo._individualType <: GAGenotype
        if head != false
            @warn "head length set in a non Genetic Programming problem, head field will be dismissed"
        end
    elseif experiment._experimentInfo._individualType <: GPGenotype
        if head != false
            if !(typeof(head) <: Integer)
                error("head field must be an integer number greater than 0")
            end
        else
            @warn "head field not set in a GEP problem, head length will be set to 7 by default"
            head = 7
        end
    end

    setGEPInfo(nodesFile=nodesFile, head=head, genj=experiment)
end

"""
    parseGPExperimentInfo(GPExperimentInfoDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseGPExperimentInfo(GPExperimentInfoDict::Dict, experiment::GeneticJuliaStructure)
    if experiment._experimentInfo._individualType <: GAGenotype
        @warn "Information about Genetic Programming problem has been given for a non-GP problem ($individualType), this information will be dismissed"
    elseif experiment._experimentInfo._individualType <: CGPGenotype
        parseCGPInfo(GPExperimentInfoDict, experiment)
    elseif experiment._experimentInfo._individualType <: STGPGenotype
        parseSTGPInfo(GPExperimentInfoDict, experiment)
    elseif experiment._experimentInfo._individualType <: GEPGenotype
        parseGEPInfo(GPExperimentInfoDict, experiment)
    #elseif experiment._experimentInfo._individualType <: GEGenotype
        #parseGEInfo(GPExperimentInfoDict, experiment)
    end
end # function

"""
    parseExperimentInfo(experimentInfo::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseExperimentInfo(experimentInfoDict::Dict, experiment::GeneticJuliaStructure)
    individualType = get(experimentInfoDict, "individualType", false)
    if individualType != false
        individualType = eval(Meta.parse(individualType))
        setIndividualType(individualType, genj=experiment)
    else
        error("individualType field in ExperimentInfo is mandatory and must be set")
    end

    gpExperimentInfoDict = get(experimentInfoDict, "GPExperimentInfo", false)
    if individualType <: GPGenotype && gpExperimentInfoDict == false
        @warn "Information about Genetic Programming problem has not been given for a GP problem ($individualType), this information must be specified in the 'GPExperimentInfo' field, default information will be used for the experiment (see setGPExperimentInfo documentation)"
        setGPExperimentInfo(genj=experiment)
    elseif individualType <: GPGenotype && gpExperimentInfoDict != false
        parseGPExperimentInfo(gpExperimentInfoDict, experiment)
    elseif gpExperimentInfoDict != false
        @warn "Information about Genetic Programming problem has been given for a non-GP problem ($individualType), this information will be dismissed"
    end

    randomSeed = get(experimentInfoDict, "randomSeed", false)
    if randomSeed != false
        if !(typeof(randomSeed) <: Integer) || randomSeed < 0
            error("Random seed must be an integer number greater or equal than 0")
        end
        setRandomSeed(randomSeed, genj=experiment)
    end

    rng = get(experimentInfoDict, "rng", false)
    if rng != false
        try
            aux = eval(Meta.parse(rng))
        catch e
            error("The following random number generator is not defined: '$rng'")
        end
        if typeof(aux) <: Random.AbstractRNG
            rng = aux
            setRNG(rng, genj=experiment)
        else
            error("The following random number generator is not a random number generator: '$rng'")
        end
    end

    if rng != false && randomSeed != false
        @warn "Random number generator and random seed have been provided, random seed will be dismissed"
    elseif rng == false && randomSeed == false
        setRNG(Random.GLOBAL_RNG, genj=experiment)
    end

    algorithm = get(experimentInfoDict, "algorithm", false)

    if algorithm == false
        algorithm = basicExperiment
    else
        aux = 0
        try
            aux = eval(Meta.parse(algorithm))
        catch e
            error("The following function of ExperimentInfo is not defined: '$algorithm'")
        end
        if typeof(aux) <: Function
            algorithm = aux
        else
            error("The following function of ExperimentInfo is not a function: '$algorithm'")
        end
    end

    varArgs = get(experimentInfoDict, "arguments", [])

    if !(typeof(varArgs) <: Array)
        varArgs = [varArgs]
    end

    setAlgorithm(algorithm, varArgs..., genj=experiment)
end # function

"""
    parseFitnessFunction(fitnessFunctionDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseFitnessFunction(fitnessFunctionDict::Dict, index::Integer)
    aux = 0

    fitnessFunction = get(fitnessFunctionDict, "function", false)
    if fitnessFunction == false
        error("No function specified for FitnessFunction dictionary number $index")
    end

    weight = get(fitnessFunctionDict, "weight", 1.0)
    if !(typeof(weight) <: Real)
        error("weight field in fitness function number $index must be a real number")
    end

    varArgs = get(fitnessFunctionDict, "arguments", [])

    if !(typeof(varArgs) <: Array)
        varArgs = [varArgs]
    end
    for i=1:length(varArgs)
        try
            arg = Meta.parse(varArgs[i])
            varArgs[i] = eval(arg)
        catch e
            # Empty
        end
    end

    try
        aux = eval(Meta.parse(fitnessFunction))
    catch e
        error("The following function of Evaluator is not defined in FitnessFunction dictionary number $index: '$(fitnessFunction)'")
    end
    if typeof(aux) <: Function
        fitnessFunction = aux
    else
        error("The following function of Evaluator is not a function in FitnessFunction dictionary number $index: '$(fitnessFunction)'")
    end

    return FitnessFunction(fitnessFunction, varArgs..., weight = weight)
end # function

"""
    parseEvaluator(evaluatorDict::Dict)

documentation
"""
function parseEvaluator(evaluatorDict::Dict, experiment::GeneticJuliaStructure)

    aux = 0

    globalFitnessFunction = get(evaluatorDict, "globalFitnessFunction", false)
    globalFitnessFunctionArgs = get(evaluatorDict, "globalFitnessFunctionArgs", [])

    if !(typeof(globalFitnessFunctionArgs) <: Array)
        globalFitnessFunctionArgs = [globalFitnessFunctionArgs]
    end

    if globalFitnessFunction != false
        try
            aux = eval(Meta.parse(globalFitnessFunction))
        catch e
            error("The following function of Evaluator is not defined: '$(globalFitnessFunction)'")
        end
        if typeof(aux) <: Function
            globalFitnessFunction = FitnessFunction(aux, globalFitnessFunctionArgs...)
        else
            error("The following function of Evaluator is not a function: '$(globalFitnessFunction)'")
        end
    else
        globalFitnessFunction = FitnessFunction(noFunc)
    end

    fitnessComparisonMode = get(evaluatorDict, "fitnessComparisonMode", "rawW")
    if !(typeof(fitnessComparisonMode) <: String)
        error("fitnessComparisonMode field in Evaluator must be a string: ($fitnessComparisonMode)")
    end

    fitnessFunctionsStrings = get(evaluatorDict, "fitnessFunctions", [])

    if !(typeof(fitnessFunctionsStrings) <: Array)
        fitnessFunctionsStrings = [fitnessFunctionsStrings]
    end

    lenFitnessFunctions = length(fitnessFunctionsStrings)
    if lenFitnessFunctions == 0
        error("A fitness function must be provided in JSON-> GeneticJuliaStructure-> Evaluator")
    end

    if typeof(fitnessFunctionsStrings[1]) <: Dict
        nFitnessFunctions = length(fitnessFunctionsStrings)
        fitnessFunctions = Array{FitnessFunction}(undef, nFitnessFunctions)
        for i=1:nFitnessFunctions
            fitnessFunctions[i] = parseFitnessFunction(fitnessFunctionsStrings[i]["FitnessFunction"], i)
        end

        setEvaluator(fitnessFunctions, genj=experiment, globalFitnessFunction=globalFitnessFunction, compareFunctionArgs=fitnessComparisonMode)
    else
        weights = Array{Real}(undef,1)
        weightsNumber = get(evaluatorDict, "weights", [])

        if !(typeof(weightsNumber) <: Array)
            weightsNumber = [weightsNumber]
        end
        try
            weights = convert(Array{Real}, weightsNumber)
        catch e
            if isa(e, MethodError)
                error("weights field must be an array of numbers")
            end
        end

        fitnessFunctions = Array{Function}(undef, lenFitnessFunctions)

        for i=1:lenFitnessFunctions
            try
                aux = eval(Meta.parse(fitnessFunctionsStrings[i]))
            catch e
                error("The following function of Evaluator is not defined: '$(fitnessFunctionsStrings[i])'")
            end
            if typeof(aux) <: Function
                fitnessFunctions[i] = aux
            else
                error("The following function of Evaluator is not a function: '$(fitnessFunctionsStrings[i])'")
            end
        end

        lenWeights = length(weights)

        if lenWeights != 0 && lenWeights != lenFitnessFunctions
            error("If weights are provided, they must be in the same amount as fitness functions")
        end

        if lenWeights == 0
            setEvaluator(fitnessFunctions, genj=experiment,
                         globalFitnessFunction=globalFitnessFunction,
                         compareFunctionArgs=fitnessComparisonMode)
        else
            setEvaluator(fitnessFunctions, weights, genj=experiment,
                         globalFitnessFunction=globalFitnessFunction,
                         compareFunctionArgs=fitnessComparisonMode)
        end
    end
end # function

"""
    parseGenerator(generatorDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseGenerator(generatorDict::Dict, experiment::GeneticJuliaStructure)
    popSize = get(generatorDict, "popSize", 100)
    if !(typeof(popSize) <: Integer)
        error("popSize field ($popSize) must be an integer number")
    end

    method = get(generatorDict, "method", false)
    if method == false
        error("Generation method must always be provided")
    end

    generateOneByOne = get(generatorDict, "generateOneByOne", true)
    if generateOneByOne == 1
        generateOneByOne = true
    elseif generateOneByOne == 0
        generateOneByOne = false
    end
    if !(typeof(generateOneByOne) <: Bool)
        error("generateOneByOne field must be boolean, true if generation method creates individuals one by one or false if it generates the whole population")
    end

    aux = 0
    try
        aux = eval(Meta.parse(method))
    catch e
        error("The following function of Generator: '$method', is not defined")
    end
    if typeof(aux) <: Function
        method = aux
    else
        error("The following function of Generator: '$method', is not a function")
    end

    varArgs = get(generatorDict, "arguments", [])

    if !(typeof(varArgs) <: Array)
        varArgs = [varArgs]
    end

    setGenerator(method, varArgs..., genj=experiment, popSize=popSize, generateOneByOne=generateOneByOne)
end # function

"""
    parseSelector(selectorDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseSelector(selectorDict::Dict, experiment::GeneticJuliaStructure)

    method = get(selectorDict, "method", false)
    if method == false
        aux = getDefaultSelector()
        method = aux[1]
        varArgs = aux[2]
        @warn "Method for selection has not been provided, functionArgs or individualMode introduced will be dismissed, '$method' will be the default selection method"

        nSelected = get(selectorDict, "nSelected", 1.0)

        if !(typeof(nSelected) <: Real)
            error("nSelected field ($nSelected) must be a number")
        end

        samplingWithRep = get(selectorDict, "samplingWithRep", true)

        if samplingWithRep == 1
            samplingWithRep = true
        elseif samplingWithRep == 0
            samplingWithRep = false
        end

        if !(typeof(samplingWithRep) <: Bool)
            error("samplingWithRep field in Selector must be boolean")
        end

        setSelector(method,varArgs..., genj=experiment, nSelected=nSelected,
                    samplingWithRep=samplingWithRep)
    else
        aux = 0
        try
            aux = eval(Meta.parse(method))
        catch e
            error("The following function of Selector: '$method', is not defined")
        end
        if typeof(aux) <: Function
            method = aux
        else
            error("The following function of Selector: '$method', is not a function")
        end

        nSelected = get(selectorDict, "nSelected", 1.0)

        if !(typeof(nSelected) <: Real)
            error("nSelected field ($nSelected) must be a number")
        end

        needsComparison = get(selectorDict, "needsComparison", true)

        if needsComparison == 1
            needsComparison = true
        elseif needsComparison == 0
            needsComparison = false
        end

        if !(typeof(needsComparison) <: Bool)
            error("needsComparison field in Selector must be boolean")
        end

        individualMode = get(selectorDict, "individualMode", "fitness")

        if !(typeof(individualMode) <: String)
            error("individualMode field in Selector ($individualMode) must be a string")
        end

        samplingWithRep = get(selectorDict, "samplingWithRep", true)

        if samplingWithRep == 1
            samplingWithRep = true
        elseif samplingWithRep == 0
            samplingWithRep = false
        end

        if !(typeof(samplingWithRep) <: Bool)
            error("samplingWithRep field in Selector must be boolean")
        end

        varArgs = get(selectorDict, "arguments", [])

        if !(typeof(varArgs) <: Array)
            varArgs = [varArgs]
        end

        setSelector(method, varArgs..., genj=experiment, nSelected=nSelected,
                    needsComparison=needsComparison, samplingWithRep=samplingWithRep)
    end
end # function

"""
    parseCrossover(crossoverDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseCrossover(crossoverDict::Dict, experiment::GeneticJuliaStructure)

    method = get(crossoverDict, "method", false)
    if method == false
        type = experiment._experimentInfo._individualType
        if !isdefined(experiment._experimentInfo, :_individualType) || !isPredefinedType(type)
            error("Either a crossover method or one of the predefined types(BinaryGA, IntegerGA, CanonicalGP) must be provided. The predefined types can be set in individualType field.")
        end
        aux = getDefaultCrossoverOp(type)
        method = aux[1]
        varArgs = aux[2]
        @warn "Method for crossover has not been provided, functionArgs introduced will be dismissed, '$method' will be the default crossover method"

        probability = get(crossoverDict, "probability", 0.7)

        if !(typeof(probability) <: AbstractFloat)
            error("probability field ($probability) must be a floating point number")
        end

        nParents = get(crossoverDict, "nParents", 2)

        if !(typeof(nParents) <: Integer)
            error("nParents field ($nParents) must be an integer number")
        end

        nChildren = get(crossoverDict, "nChildren", -1)

        if !(typeof(nChildren) <: Integer)
            error("nChildren field ($nChildren) must be an integer number")
        end

        setCrossoverOperator(method, varArgs..., genj=experiment, probability=probability, nParents=nParents, nChildren=nChildren)

    else
        aux = 0
        try
            aux = eval(Meta.parse(method))
        catch e
            error("The following function of Crossover: '$method', is not defined")
        end
        if typeof(aux) <: Function
            method = aux
        else
            error("The following function of Crossover: '$method', is not a function")
        end

        probability = get(crossoverDict, "probability", 0.7)

        if !(typeof(probability) <: AbstractFloat)
            error("probability field ($probability) must be a floating point number")
        end

        nParents = get(crossoverDict, "nParents", 2)

        if !(typeof(nParents) <: Integer)
            error("nParents field ($nParents) must be an integer number")
        end

        nChildren = get(crossoverDict, "nChildren", -1)

        if !(typeof(nChildren) <: Integer)
            error("nChildren field ($nChildren) must be an integer number")
        end

        varArgs = get(crossoverDict, "arguments", [])

        if !(typeof(varArgs) <: Array)
            varArgs = [varArgs]
        end

        setCrossoverOperator(method, varArgs..., genj=experiment, probability=probability, nParents=nParents, nChildren=nChildren)

    end
end # function

"""
    parseMutation(mutationDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseMutation(mutationDict::Dict, experiment::GeneticJuliaStructure)

    method = get(mutationDict, "method", false)
    if method == false
        type = experiment._experimentInfo._individualType
        if !isdefined(experiment._experimentInfo, :_individualType) || !isPredefinedType(type)
            error("Either a mutation method or one of the predefined types(BinaryGA, IntegerGA, CanonicalGP) must be provided. The predefined types can be set in individualType field.")
        end
        aux = getDefaultMutationOp(type)
        method = aux[1]
        varArgs = aux[2]
        @warn "Method for mutation has not been provided, functionArgs introduced will be dismissed, '$method' will be the default mutation method"

        probability = get(mutationDict, "probability", 0.7)

        if !(typeof(probability) <: AbstractFloat)
            error("probability field ($probability) must be a floating point number")
        end

        setMutationOperator(method, varArgs..., genj=experiment, probability=probability)

    else
        aux = 0
        try
            aux = eval(Meta.parse(method))
        catch e
            error("The following function of Mutation: '$method', is not defined")
        end
        if typeof(aux) <: Function
            method = aux
        else
            error("The following function of Mutation: '$method', is not a function")
        end

        probability = get(mutationDict, "probability", 0.7)

        if !(typeof(probability) <: AbstractFloat)
            error("probability field ($probability) must be a floating point number")
        end

        varArgs = get(mutationDict, "arguments", [])

        if !(typeof(varArgs) <: Array)
            varArgs = [varArgs]
        end

        setMutationOperator(method, varArgs..., genj=experiment, probability=probability)

    end
end # function

"""
    parseReplacement(replacementDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseReplacement(replacementDict::Dict, experiment::GeneticJuliaStructure)

    method = get(replacementDict, "method", false)
    if method == false
        type = experiment._experimentInfo._individualType
        if !isdefined(experiment._experimentInfo, :_individualType) || !isPredefinedType(type)
            error("Either a replacement method or one of the predefined types(BinaryGA, IntegerGA, CanonicalGP) must be provided. The predefined types can be set in individualType field.")
        end
        aux = getDefaultMutationOp(type)
        method = aux[1]
        varArgs = aux[2]
        @warn "Method for replacement has not been provided, functionArgs introduced will be dismissed, '$method' will be the default replacement method"

        setReplacementOperator(method, varArgs..., genj=experiment)

    else
        aux = 0
        try
            aux = eval(Meta.parse(method))
        catch e
            error("The following function of Replacement: '$method', is not defined")
        end
        if typeof(aux) <: Function
            method = aux
        else
            error("The following function of Replacement: '$method', is not a function")
        end

        needsComparison = get(replacementDict, "needsComparison", false)

        if needsComparison == 1
            needsComparison = true
        elseif needsComparison == 0
            needsComparison = false
        end

        if !(typeof(needsComparison) <: Bool)
            error("needsComparison field in Replacement must be boolean")
        end

        eliteSize = get(replacementDict, "eliteSize", 0)

        if !(typeof(eliteSize) <: Integer)
            error("eliteSize field in Replacement must be an integer number")
        end

        individualMode = get(replacementDict, "individualMode", "fitness")

        if !(typeof(individualMode) <: String)
            error("individualMode field in Replacement ($individualMode) must be a string")
        end

        varArgs = get(replacementDict, "arguments", [])

        if !(typeof(varArgs) <: Array)
            varArgs = [varArgs]
        end

        setReplacementOperator(method, varArgs..., genj=experiment, needsComparison=needsComparison, eliteSize=eliteSize)
    end
end # function

"""
    parseStopCondition(stopConditionDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseStopCondition(stopConditionDict::Dict, experiment::GeneticJuliaStructure)

    maxEvaluations = get(stopConditionDict, "maxEvaluations", -1)

    if !(typeof(maxEvaluations) <: Integer)
        error("maxEvaluations field ($maxEvaluations) must be an integer number")
    end

    maxIterations = get(stopConditionDict, "maxIterations", -1)

    if !(typeof(maxIterations) <: Integer)
        error("maxIterations field ($maxIterations) must be an integer number")
    end

    maxIterNotImproving = get(stopConditionDict, "maxIterNotImproving", -1)

    if !(typeof(maxIterNotImproving) <: Integer)
        error("maxIterNotImproving field ($maxIterNotImproving) must be an integer number")
    end

    maxTime = get(stopConditionDict, "maxTime", Inf)

    if !(typeof(maxTime) <: AbstractFloat)
        error("maxTime field ($maxTime) must be a number")
    end

    setStopCondition(genj=experiment, maxEvaluations=maxEvaluations, maxIterations=maxIterations, maxIterNotImproving=maxIterNotImproving, maxTime=maxTime)
end # function

"""
    parseExperimentSummary(summaryDict::Dict, experiment::GeneticJuliaStructure)

documentation
"""
function parseExperimentSummary(summaryDict::Dict, experiment::GeneticJuliaStructure)
    outputFile = get(summaryDict, "outputFile", "")
    if !(typeof(outputFile) <: String)
        error("outputFile field in Experiment Summary must be a string")
    end

    batchSize = get(summaryDict, "batchSize", 0)
    if !(typeof(batchSize) <: Integer)
        error("batchSize field in Experiment Summary must be an integer number")
    end

    printDuringExperiment = get(summaryDict, "printDuringExperiment", false)
    if !(typeof(printDuringExperiment) <: Bool)
        error("printDuringExperiment field in Experiment Summary must be either true or false")
    end

    displayFitness = get(summaryDict, "displayFitness", false)
    if !(typeof(displayFitness) <: Bool)
        error("displayFitness field in Experiment Summary must be either true or false")
    end

    displayBestFitness = get(summaryDict, "displayBestFitness", false)
    if !(typeof(displayBestFitness) <: Bool)
        error("displayBestFitness field in Experiment Summary must be either true or false")
    end

    displayFitnessMean = get(summaryDict, "displayFitnessMean", false)
    if !(typeof(displayFitnessMean) <: Bool)
        error("displayFitnessMean field in Experiment Summary must be either true or false")
    end

    displayFitnessSTD = get(summaryDict, "displayFitnessSTD", false)
    if !(typeof(displayFitnessSTD) <: Bool)
        error("displayFitnessSTD field in Experiment Summary must be either true or false")
    end

    setExperimentSummary(experiment, outputFile=outputFile, batchSize=batchSize,
                         printDuringExperiment=printDuringExperiment, displayFitness=displayFitness,
                         displayBestFitness=displayBestFitness, displayFitnessMean=displayFitnessMean,
                         displayFitnessSTD=displayFitnessSTD)
end # function

"""
    generateMainStructure(genJ::GeneticJuliaStructure, jsonFile::String)

documentation
"""
function generateMainStructure(jsonFile::String; verbose::Bool = true)
    file=open(jsonFile)
    dictionary = JSON.parse(file)
    close(file)

    nExperiments = length(dictionary["Experiments"])
    experiments = Array{GeneticJuliaStructure}(undef, nExperiments)

    if verbose
        println("\n")
        println("Reading JSON file")
        println("##########################")
    end
    for i=1:nExperiments
        if verbose
            print("\t· Experiment $i ")
        end
        experiments[i] = GeneticJuliaStructure()

        experiments[i]._experimentInfo = ExperimentInfo()

        experimentInfoDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "ExperimentInfo", false)
        if experimentInfoDict != false
            parseExperimentInfo(experimentInfoDict, experiments[i])
        end

        evaluator = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Evaluator", false)
        if evaluator != false
            parseEvaluator(evaluator, experiments[i])
        else
            error("Evaluator is a mandatory field in GeneticJuliaStructure")
        end

        generatorDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Generator", false)
        if generatorDict != false
            parseGenerator(generatorDict, experiments[i])
        else
            error("Generator is a mandatory field in GeneticJuliaStructure")
        end

        selectorDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Selector", false)
        if selectorDict != false
            parseSelector(selectorDict, experiments[i])
        end

        crossoverDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Crossover", false)
        if crossoverDict != false
            parseCrossover(crossoverDict, experiments[i])
        end

        mutationDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Mutation", false)
        if mutationDict != false
            parseMutation(mutationDict, experiments[i])
        end

        replacementDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "Replacement", false)
        if replacementDict != false
            parseReplacement(replacementDict, experiments[i])
        end

        stopConditionDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "StopConditions", false)
        if stopConditionDict != false
            parseStopCondition(stopConditionDict, experiments[i])
        end

        summaryDict = get(dictionary["Experiments"][i]["GeneticJuliaStructure"], "ExperimentSummary", false)
        if summaryDict != false
            parseExperimentSummary(summaryDict, experiments[i])
        end

        if verbose
            println("→ Experiment correctly read ✔")
        end
    end
    if verbose
        println("##########################\n")
    end

    return experiments
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using ArgParse

s = ArgParseSettings(description="Run simulations of oceanic free convection.")

@add_arg_table s begin
    "--resolution", "-N"
        arg_type=Int
        required=true
        help="Number of grid points in each dimension (Nx, Ny, Nz) = (N, N, N)."
    "--heat-flux"
        arg_type=Float64
        required=true
        help="Heat flux to impose at the surface [W/m²]. Negative values imply a cooling flux."
    "--dTdz"
        arg_type=Float64
        required=true
        help="Temperature gradient to impose at the bottom [K/m]."
    "--diffusivity"
        arg_type=Float64
        required=true
        help="Diffusivity κ [m²/s]."
    "--dt"
        arg_type=Float64
        required=true
        help="Time step in seconds."
    "--days"
        arg_type=Float64
        required=true
        help="Number of simulation days to run the model."
    "--output-dir", "-d"
        arg_type=AbstractString
        required=true
        help="Base directory to save output to."
end

parsed_args = parse_args(s)
N = parsed_args["resolution"]
Q = parsed_args["heat-flux"]
dTdz = parsed_args["dTdz"]
κ = parsed_args["diffusivity"]

dt, days = parsed_args["dt"], parsed_args["days"]
dt = isinteger(dt) ? Int(dt) : dt
days = isinteger(days) ? Int(days) : days

base_dir = parsed_args["output-dir"]

filename_prefix = "free_convection_N" * string(N) * "_Q" * string(Q) *
                  "_dTdz" * string(dTdz) * "_k" * string(κ) * "_dt" * string(dt) *
                  "_days" * string(days)
output_dir = joinpath(base_dir, filename_prefix)

if !isdir(output_dir)
    println("Creating directory: $output_dir")
    mkpath(output_dir)
end

# Adding a second worker/proccessor for asynchronous NetCDF output writing.
using Distributed
addprocs(1)

# We need to activate the Oceananigans environment on all workers if executing
# this script from the development repository using "julia --project".
@everywhere using Pkg
@everywhere Pkg.activate(".");

@everywhere using Oceananigans
@everywhere using CuArrays
@everywhere using Printf

# Physical constants.
ρ₀ = 1027    # Density of seawater [kg/m³]
cₚ = 4181.3  # Specific heat capacity of seawater at constant pressure [J/(kg·K)]

# Set more simulation parameters.
Nx, Ny, Nz = N, N, N
Lx, Ly, Lz = 100, 100, 100
Δt = dt
Nt = Int(days*60*60*24/dt)
ν = κ  # This is also implicitly setting the Prandtl number Pr = 1.

# To impose a flux boundary condition, the top flux imposed should be negative
# for a heating flux and positive for a cooling flux, thus the minus sign.
top_flux = -Q / (ρ₀ * cₚ)

# Create the model.
model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=ν, κ=κ, arch=GPU(), float_type=Float64)

# Set boundary conditions.
model.boundary_conditions.T.z.left = BoundaryCondition(Flux, top_flux)
model.boundary_conditions.T.z.right = BoundaryCondition(Gradient, dTdz)

# Set initial conditions.
T_prof = 273.15 .+ 20 .+ dTdz .* model.grid.zC  # Initial temperature profile.
T_3d = repeat(reshape(T_prof, 1, 1, Nz), Nx, Ny, 1)  # Convert to a 3D array.

# Add small normally distributed random noise to the top half of the domain to
# facilitate numerical convection.
@. T_3d[:, :, 1:round(Int, Nz/2)] += 0.01*randn()

model.tracers.T.data .= CuArray(T_3d)

# Add a NaN checker diagnostic that will check for NaN values in the vertical
# velocity and temperature fields every 1,000 time steps and abort the simulation
# if NaN values are detected.
nan_checker = NaNChecker(1000, [model.velocities.w, model.tracers.T], ["w", "T"])
push!(model.diagnostics, nan_checker)

# Write full output to disk enough times that we end up with 32 outputs.
n_outputs = 32
output_freq = floor(Int, Nt / n_outputs)
netcdf_writer = NetCDFOutputWriter(dir=output_dir, prefix=filename_prefix * "_",
                                   padding=0, naming_scheme=:file_number,
                                   frequency=output_freq, async=true)
push!(model.output_writers, netcdf_writer)

# With asynchronous output writing we need to wrap our time-stepping using
# an @sync block so that Julia does not just quit if the model finishes time
# stepping while there are still output files that need to be written to disk
# on another worker/processor.
@sync begin
    # Take Ni "intermediate" time steps at a time and print out the current time
    # and average wall clock time per time step.
    Ni = 100
    for i = 1:ceil(Int, Nt/Ni)
        progress = 100 * (model.clock.iteration / Nt)  # Progress %
        @printf("[%06.2f%%] Time: %.1f / %.1f...", progress, model.clock.time, Nt*Δt)

        tic = time_ns()
        time_step!(model, Ni, Δt)

        @printf("   average wall clock time per iteration: %s\n", prettytime((time_ns() - tic) / Ni))
    end
end

ext: jl


text:
export mpi_bcast, mpi_init, mpi_recv, mpi_send,
    mpi_sendrecv, mpi_sum,  mpi_finalize, mpi_initialized, mpi_halo_exchange, mpi_halo_exchange2,
    mpi_finalized, mpi_rank, mpi_size, mpi_sync!, mpi_gather, mpi_SparseTensor, require_mpi

"""
    mpi_init()

Initialized the MPI session. `mpi_init` must be called before any `run(sess, ...)`.
"""
function mpi_init()
    if mpi_initialized()
        @warn "MPI has already been initialized"
        return
    end
    @eval ccall((:mpi_init, $LIBADCME), Cvoid, ())
end

"""
    mpi_finalize()

Finalize the MPI call.
"""
function mpi_finalize()
    mpi_check()
    @eval ccall((:mpi_finalize, $LIBADCME), Cvoid, ())
end

"""
    mpi_rank()

Returns the rank of current MPI process (rank 0 based).
"""
function mpi_rank()
    mpi_check()
    out = @eval ccall((:mpi_rank, $LIBADCME), Cint, ())
    Int64(out)
end

"""
    mpi_size()

Returns the size of MPI world.
"""
function mpi_size()
    mpi_check()
    @eval ccall((:mpi_size, $LIBADCME), Cint, ())
end

"""
    require_mpi()

Throws an error if `mpi_init()` has not been called.
"""
function require_mpi()
    if !mpi_initialized()
        error("MPI has not been initialized. Run `mpi_init()` to initialize MPI first.")
    end
end

"""
    mpi_finalized()

Returns a boolean indicating whether the current MPI session is finalized.
"""
function mpi_finalized()
    require_mpi()
    Bool(@eval ccall((:mpi_finalized, $LIBADCME), Cuchar, ()))
end

"""
    mpi_initialized()

Returns a boolean indicating whether the current MPI session is initialized.
"""
function mpi_initialized()
    Bool(@eval ccall((:mpi_initialized, $LIBADCME), Cuchar, ()))
end

"""
    mpi_sync!(message::Array{Int64,1}, root::Int64 = 0)
    mpi_sync!(message::Array{Float64,1}, root::Int64 = 0)

Sync `message` across all MPI processors.
"""
function mpi_sync!(message::Array{Int64,1}, root::Int64 = 0)
    mpi_check()
    @eval ccall((:mpi_sync, $LIBADCME), Cvoid, (Ptr{Clonglong}, Cint, Cint), $message, Int32(length($message)), Int32($root))
end

function mpi_sync!(message::Array{Float64,1}, root::Int64 = 0)
    mpi_check()
    @eval ccall((:mpi_sync_double, $LIBADCME), Cvoid, (Ptr{Cdouble}, Cint, Cint), $message, Int32(length($message)), Int32($root))
end

function mpi_check()
    require_mpi()
    if mpi_finalized()
        error("MPI has been finalized.")
    end
end

function _mpi_sum(a,root::Int64=0)
    mpi_check()
    mpisum_ = load_system_op("mpisum")
    a,root = convert_to_tensor(Any[a,root], [Float64,Int64])
    out = mpisum_(a,root)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_sum(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Sum `a` on the MPI processor `root`.
"""
function mpi_sum(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_sum(a, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_sum(a, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_sum is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_sum(a, root)
        return reshape(out, s)
    end
end

function _mpi_bcast(a,root::Int64=0)
    mpi_check()
    mpibcast_ = load_system_op("mpibcast")
    a,root = convert_to_tensor(Any[a,root], [Float64,Int64])
    out = mpibcast_(a,root)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_bcast(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Broadcast `a` from processor `root` to all other processors.
"""
function mpi_bcast(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_bcast(a, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_bcast(a, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_bcast is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_bcast(a, root)
        return reshape(out, s)
    end
end

function _mpi_send(a,dest::Int64,tag::Int64=0)
    mpisend_ = load_system_op("mpisend")
    a,dest,tag = convert_to_tensor(Any[a,dest,tag], [Float64,Int64,Int64])
    out = mpisend_(a,dest,tag)
    if !isnothing(length(a))
        return set_shape(out, (length(a),))
    else
        return out
    end
end

"""
    mpi_send(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64,root::Int64 = 0)

Sends `a` to processor `dest`. `a` itself is returned so that the send action can be added to the computational graph.
"""
function mpi_send(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64,root::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_send(a,dest, root)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_send(a, dest, root)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_send is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_send(a, dest, root)
        return reshape(out, s)
    end
end

function _mpi_recv(a,src::Int64,tag::Int64=0)
    mpirecv_ = load_system_op("mpirecv")
    a,src,tag = convert_to_tensor(Any[a,src,tag], [Float64,Int64,Int64])
    mpirecv_(a,src,tag)
end

"""
    mpi_recv(a::Union{Array{Float64}, Float64, PyObject}, src::Int64, tag::Int64 = 0)

Receives an array from processor `src`. `mpi_recv` requires an input for gradient backpropagation.
Typically we can write

```julia
r = mpi_rank()
a = constant(Float64(r))
if r==1
    a = mpi_send(a, 0)
end
if r==0
    a = mpi_recv(a, 1)
end
```

Then `a=1` on both processor 0 and processor 1.
"""
function mpi_recv(a::Union{Array{Float64}, Float64, PyObject}, src::Int64, tag::Int64 = 0)
    a = convert_to_tensor(a, dtype = Float64)
    if length(size(a))==0
        a = reshape(a, (1,))
        out = _mpi_recv(a, src, tag)
        return squeeze(out)
    elseif length(size(a))==1
        return _mpi_recv(a, src, tag)
    elseif nothing in size(a)
        error("The shape of 1st input $(size(a)) contains nothing. mpi_recv is not able to determine
        the output shape.")
    else
        s = size(a)
        a = reshape(a, (-1,))
        out = _mpi_recv(a, src, tag)
        return reshape(out, s)
    end
end

"""
    mpi_sendrecv(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64, src::Int64, tag::Int64=0)

A convenient wrapper for `mpi_send` followed by `mpi_recv`.
"""
function mpi_sendrecv(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64, src::Int64, tag::Int64=0)
    r = mpi_rank()
    @assert src != dest
    if r==src
        a = mpi_send(a, dest, tag)
    elseif r==dest
        a = mpi_recv(a, src, tag)
    end
    a
end

"""
    mpi_gather(u::Union{Array{Float64, 1}, PyObject}, deps::Union{Missing, PyObject} = missing)

Gathers all the vectors from different processes to the root process. The function returns
a long vector which concatenates of local vectors in the order of process IDs.
"""
function mpi_gather(u::Union{Array{Float64, 1}, PyObject}, deps::Union{Missing, PyObject} = missing)
    mpigather_ = load_system_op("mpigather")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    out = mpigather_(u)
    set_shape(out, (mpi_size()*length(u)))
end

function load_plugin_MPITensor()
    if !isfile("$(ADCME.LIBDIR)/libHYPRE.so")
        scriptpath = joinpath(@__DIR__, "..", "deps", "install_hypre.jl")
        include(scriptpath)
    end
    if Sys.isapple()
        oplibpath = joinpath(@__DIR__, "..", "deps", "Plugin", "MPITensor", "build", "libMPITensor.dylib")
    elseif Sys.iswindows()
        oplibpath = joinpath(@__DIR__, "..", "deps", "Plugin", "MPITensor", "build", "MPITensor.dll")
    else
        oplibpath = joinpath(@__DIR__, "..", "deps", "Plugin", "MPITensor", "build", "libMPITensor.so")
    end
    if !isfile(oplibpath)
        PWD = pwd()
        cd(joinpath(@__DIR__, "..", "deps", "Plugin", "MPITensor"))
        if !isdir("build")
        mkdir("build")
        end
        cd("build")
        cmake()
        make()
        cd(PWD)
    end
    oplibpath
end

function load_plugin_MPIHaloExchange(reach::Int64=1)
    if reach==1
        DIR = joinpath(@__DIR__, "..", "deps", "Plugin", "MPIHaloExchange")
        LIB = "HaloExchangeTwoD"
    elseif reach==2
        DIR = joinpath(@__DIR__, "..", "deps", "Plugin", "MPIHaloExchange2")
        LIB = "HaloExchangeNeighborTwo"
    end
    if Sys.isapple()
        oplibpath = joinpath(DIR, "build", "lib$LIB.dylib")
    elseif Sys.iswindows()
        oplibpath = joinpath(DIR, "build", "$LIB.dll")
    else
        oplibpath = joinpath(DIR, "build", "lib$LIB.so")
    end
    if !isfile(oplibpath)
        PWD = pwd()
        cd(DIR)
        if !isdir("build")
        mkdir("build")
        end
        cd("build")
        cmake()
        make()
        cd(PWD)
    end
    oplibpath
end

function mpi_halo_exchange_(oplibpath, u,fill_value,m,n, tag, deps)
    halo_exchange_two_d_ = load_op_and_grad(oplibpath,"halo_exchange_two_d")
    u,fill_value,m,n,tag = convert_to_tensor(Any[u,fill_value,m,n,tag], [Float64,Float64,Int64,Int64,Int64])
    deps = coalesce(deps, u[1,1])
    out = halo_exchange_two_d_(u,fill_value,m,n,tag, deps)
    set_shape(out, (size(u,1)+2, size(u,2)+2))
end

function mpi_halo_exchange2_(oplibpath, u,fill_value,m,n,tag,w)
    halo_exchange_neighbor_two_ = load_op_and_grad(oplibpath,"halo_exchange_neighbor_two")
    w = coalesce(w, u[1,1])
    u,fill_value,m,n,tag,w = convert_to_tensor(Any[u,fill_value,m,n,tag,w], [Float64,Float64,Int64,Int64,Int64,Float64])
    out = halo_exchange_neighbor_two_(u,fill_value,m,n,tag,w)
    set_shape(out, (size(u,1)+4, size(u,2)+4))
end

function mpi_create_matrix(oplibpath, indices,values,ilower,iupper)
    mpi_create_matrix_ = load_op_and_grad(oplibpath,"mpi_create_matrix", multiple=true)
    indices,values,ilower,iupper = convert_to_tensor(Any[indices,values,ilower,iupper], [Int64,Float64,Int64,Int64])
    mpi_create_matrix_(indices,values,ilower,iupper)
end

function mpi_get_matrix(oplibpath, rows,ncols,cols,ilower,iupper,values, N)
    mpi_get_matrix_ = load_op_and_grad(oplibpath,"mpi_get_matrix", multiple=true)
    rows,ncols,cols,ilower_,iupper_,values = convert_to_tensor(Any[rows,ncols,cols,ilower,iupper,values], [Int32,Int32,Int32,Int64,Int64,Float64])
    indices, vals = mpi_get_matrix_(rows,ncols,cols,ilower_,iupper_,values)
    SparseTensor(tf.SparseTensor(indices, vals, (iupper-ilower+1, N)), false)
end

function mpi_tensor_solve(oplibpath, rows,ncols,cols,values,rhs,ilower,iupper,solver,printlevel)
    mpi_tensor_solve_ = load_op_and_grad(oplibpath,"mpi_tensor_solve")
    rows,ncols,cols,values,rhs,ilower,iupper,printlevel = convert_to_tensor(Any[rows,ncols,cols,values,rhs,ilower,iupper,printlevel], [Int32,Int32,Int32,Float64,Float64,Int64,Int64,Int64])
    mpi_tensor_solve_(rows,ncols,cols,values,rhs,ilower,iupper,solver,printlevel)
end

@doc raw"""
    mutable struct mpi_SparseTensor
        rows::PyObject
        ncols::PyObject
        cols::PyObject
        values::PyObject
        ilower::Int64
        iupper::Int64
        N::Int64
        oplibpath::String
    end

A structure to hold local data of a sparse matrix. The global matrix is assumed to be a $M\times N$ square matrix.
The current processor owns rows from `ilower` to `iupper` (inclusive). The data is specified by

- `rows`: an array indicating the rows that contain nonzero values. Note `rows ≥ ilower`.
- `ncols`: an array indicating the number of nonzero values for each row in `rows`.
- `cols`: the column indices for nonzero values. Its length is $\sum_{i=1}^{\mathrm{ncols}} \mathrm{ncols}_i$
- `vals`: the nonzero values corresponding to each column index in `cols`
- `oplibpath`: the backend library (returned by `ADCME.load_plugin_MPITensor`)

All data structure are 0-based. Note if we work with a linear solver, $M=N$.

For example, consider the sparse matrix

```
[  1 0 0 1  ]
[  0 1 2 1  ]
```

We have

```julia
rows = Int32[0;1]
ncols = Int32[2;3]
cols = Int32[0;3;1,2,3]
values = [1.;1.;1.;2.;1.]
iupper = ilower + 2
```
"""
mutable struct mpi_SparseTensor
    rows::PyObject
    ncols::PyObject
    cols::PyObject
    values::PyObject
    ilower::Int64
    iupper::Int64
    N::Int64
    oplibpath::String
end

function Base.:show(io::IO, sp::mpi_SparseTensor)
    if isnothing(length(sp.values))
        len = "?"
    else
        len = length(sp.values)
    end
    print("mpi_SparseTensor($(sp.iupper - sp.ilower + 1), $(sp.N)), range = [$(sp.ilower), $(sp.iupper)], nnz = $(len)")
end

function mpi_SparseTensor(indices::PyObject, values::PyObject, ilower::Int64, iupper::Int64, N::Int64)
    @assert ilower >=0
    @assert ilower <= iupper
    @assert iupper <= N
    oplibpath = load_plugin_MPITensor()
    rows, ncols, cols, out = mpi_create_matrix(oplibpath, indices,values,ilower,iupper)
    mpi_SparseTensor(rows, ncols, cols, out, ilower, iupper, N, oplibpath)
end

@doc raw"""
    mpi_SparseTensor(rows::Union{Array{Int32,1}, PyObject}, ncols::Union{Array{Int32,1}, PyObject}, cols::Union{Array{Int32,1}, PyObject},
        vals::Union{Array{Float64,1}, PyObject}, ilower::Int64, iupper::Int64, N::Int64)

Create a $N\times N$ distributed sparse tensor `A` for the current MPI processor. The current MPI processor owns rows with indices `[ilower, iupper]`.
The submatrix is specified using the CSR format.

- `rows`: an array indicating the rows that contain nonzero values. Note `rows ≥ ilower`.
- `ncols`: an array indicating the number of nonzero values for each row in `rows`.
- `cols`: the column indices for nonzero values. Its length is $\sum_{i=1}^{\mathrm{ncols}} \mathrm{ncols}_i$
- `vals`: the nonzero values corresponding to each column index in `cols`

Note that by default the indices are zero-based.
"""
function mpi_SparseTensor(rows::Union{Array{Int32,1}, PyObject}, ncols::Union{Array{Int32,1}, PyObject}, cols::Union{Array{Int32,1}, PyObject},
     vals::Union{Array{Float64,1}, PyObject}, ilower::Int64, iupper::Int64, N::Int64)
    @assert ilower >=0
    @assert ilower <= iupper
    @assert iupper <= N
    oplibpath = load_plugin_MPITensor()
    rows, ncols, cols, vals = convert_to_tensor(
        [rows, ncols, cols, vals],
        [Int32, Int32, Int32, Float64]
    )
    mpi_SparseTensor(rows, ncols, cols, vals, ilower, iupper, N, oplibpath)
end

"""
    mpi_SparseTensor(sp::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}},
        ilower::Union{Int64, Missing} = missing,
        iupper::Union{Int64, Missing} = missing)

Constructing `mpi_SparseTensor` from a `SparseTensor` or a sparse Array.
"""
function mpi_SparseTensor(sp::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}},
    ilower::Union{Int64, Missing} = missing,
    iupper::Union{Int64, Missing} = missing)
    sp = constant(sp)
    ilower = coalesce(ilower, 0)
    iupper = coalesce(iupper, size(sp, 1)-1)
    N = size(sp, 2)
    mpi_SparseTensor(sp.o.indices, sp.o.values, ilower, iupper, N)
end

function LinearAlgebra.:\(sp::mpi_SparseTensor, b::Union{Array{Float64, 1}, PyObject})
    b = constant(b)
    @assert length(b)==sp.iupper - sp.ilower + 1
    out = mpi_tensor_solve(sp.oplibpath, sp.rows,sp.ncols,
                sp.cols,sp.values,b,
                sp.ilower,sp.iupper,options.mpi.solver, options.mpi.printlevel)
    set_shape(out, (sp.iupper-sp.ilower+1,))
end

function SparseTensor(sp::mpi_SparseTensor)
    mpi_get_matrix(sp.oplibpath, sp.rows,sp.ncols,sp.cols,sp.ilower,sp.iupper,sp.values, sp.N)
end

function Base.:run(sess::PyObject, sp::mpi_SparseTensor)
    run(sess, SparseTensor(sp))
end

@doc raw"""
    mpi_halo_exchange(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
    fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)

Perform Halo exchnage on `u` (a $k \times k$ matrix). The output has a shape $(k+2)\times (k+2)$

- `fill_value`: value used for the boundaries
- `tag`: message tag
- `deps`: a **scalar** tensor; it can be used to serialize the MPI calls
"""
function mpi_halo_exchange(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
    fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)
    @assert size(u,1)==size(u,2)
    oplibpath = load_plugin_MPIHaloExchange()
    mpi_halo_exchange_(oplibpath, u, fill_value, m, n, tag, deps)
end

@doc raw"""
    mpi_halo_exchange2(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
    fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)

Similar to [`mpi_halo_exchange`](@ref), but the reach is 2, i.e., for a $N\times N$ matrix $u$, the output will be a
$(N+4)\times (N+4)$ matrix.
"""
function mpi_halo_exchange2(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
    fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)
    @assert size(u,1)==size(u,2)
    oplibpath = load_plugin_MPIHaloExchange(2)
    mpi_halo_exchange2_(oplibpath, u, fill_value, m, n, tag, deps)
end

function mpi_tensor_transpose(oplibpath, row,col,ncol,val,n,rank,nt)
    require_mpi()
    mpi_tensor_transpose_ = load_op_and_grad(oplibpath,"mpi_tensor_transpose", multiple=true)
    row,col,ncol,val,n,rank,nt = convert_to_tensor(Any[row,col,ncol,val,n,rank,nt], [Int32,Int32,Int32,Float64,Int64,Int64,Int64])
    indices, vals = mpi_tensor_transpose_(row,col,ncol,val,n,rank,nt)
end

"""
    adjoint(A::mpi_SparseTensor)

Returns the adjoint of `A`, i.e., `A'`. Each MPI rank owns the same number of rows.
"""
function Base.:adjoint(A::mpi_SparseTensor)
    oplibpath = load_plugin_MPITensor()
    n = A.iupper - A.ilower + 1
    indices, vals = mpi_tensor_transpose(oplibpath, A.rows, A.cols, A.ncols, A.values, n, mpi_rank(), A.N)
    sp = RawSparseTensor(indices, vals, n, A.N)
    mpi_SparseTensor(sp, A.ilower, A.iupper)
end


text:
"""
    lefttransfer(HL₀::AbstractMatrix, hL::Union{Nothing, AbstractPiecewise},
                    Ψ₁::FiniteCMPS, Ψ₂::FiniteCMPS = Ψ₁;
                    Kmax = 50, tol = defaulttol(Ψ₁))

Solve the non-homogeneous (or homogeneous if hL === nothing) Lindblad equation

```math
∂ₓ(HL) = Q₂' * HL + HL * Q₁ + R₂' * HL * R₁ + hL
```

with `Q₁` and `R₁` the matrices of the finite cMPS `Ψ₁` (similar for `Q₂` and `R₂`, where
`Ψ₂ = Ψ₁` by default), over the whole domain `(a,b) = domain(Ψ₁)`, starting from the initial
condition `HL(a) = HL₀`, all the way up to `b`. The states `Ψ₁` and `Ψ₂` are assumed to have
equal domain and grid; some checks are in place.

The result is given by

```math
⟨HL(x)| = ⟨HL₀|exp(∫ₐˣ T(z) dz) + ∫ₐˣ [ ⟨hL(u)|exp(∫ᵤˣ T(z) dz) ] du
```

and returned as `HL, converged` with `HL::Piecewise` and where `converged::Bool` is `true`
or `false` depending on convergence up to the requested tolerance.
"""
function lefttransfer(HL₀::AbstractMatrix, hL::Union{Nothing, AbstractPiecewise},
                        Ψ₁::FiniteCMPS, Ψ₂::FiniteCMPS = Ψ₁;
                        Kmax = 50, tol = defaulttol(Ψ₁))

    Q₁ = Ψ₁.Q
    R₁s = Ψ₁.Rs
    Q₂ = Ψ₂.Q
    R₂s = Ψ₂.Rs
    grid = nodes(Q₁)
    N = length(grid) - 1
    a = first(grid)
    b = last(grid)
    if Ψ₁ !== Ψ₂
        grid == nodes(Q₂) || throw(DomainMismatch())
    end
    if !isnothing(hL)
        grid == nodes(hL) || throw(DomainMismatch())
    end
    HLᵢ = TaylorSeries([HL₀], a)
    HLs = Vector{typeof(HLᵢ)}(undef, N)
    HLs[1] = HLᵢ
    infoL = true
    for i = 1:N
        HLᵢ = HLs[i]
        xᵢ = grid[i]
        xⱼ = grid[i+1]
        Δxᵢ = grid[i+1] - grid[i]
        Q₁ᵢ = shift!(Q₁[i], xᵢ)
        R₁ᵢs = shift!.(getindex.(R₁s, i), xᵢ)
        Q₂ᵢ = shift!(Q₂[i], xᵢ)
        R₂ᵢs = shift!.(getindex.(R₂s, i), xᵢ)
        hLᵢ = isnothing(hL) ? nothing : shift(hL[i], xᵢ)

        # build Taylor coefficients (i.e. solve triangular problem)
        HLᵢ, converged = _lefttransfer!(HLᵢ, hLᵢ, Q₁ᵢ, R₁ᵢs, Δxᵢ, Q₂ᵢ, R₂ᵢs;
                                        Kmax = Kmax, tol = tol)
        infoL &= converged
        shift!(HLᵢ, (xᵢ+xⱼ)/2)

        # initialize next ρ element
        if i < N
            HLs[i+1] = TaylorSeries([HLᵢ(xⱼ)], xⱼ)
        end
    end
    HL = Piecewise(grid, HLs)
    return HL, infoL
end

"""
    righttransfer(HR₀::AbstractMatrix, hR::Union{Nothing, AbstractPiecewise},
                    Ψ₁::FiniteCMPS, Ψ₂::FiniteCMPS = Ψ₁;
                    Kmax = 50, tol = defaulttol(Ψ₁))

Solve the non-homogeneous (or homogeneous if hR === nothing) Lindblad equation

```math
-∂ₓ(HR) = Q₁ * HR + HR * Q₂' + R₁ * HR * R₂ + hR
```

with `Q₁` and `R₁` the matrices of the finite cMPS `Ψ₁` (similar for `Q₂` and `R₂`, where
`Ψ₂ = Ψ₁` by default), over the whole domain `(a,b) = domain(Ψ₁)`, starting from the initial
condition `HR(b) = HR₀`, all the way down to `a`. The states `Ψ₁` and `Ψ₂` are assumed to
have equal domain and grid; some checks are in place.

The result is given by

```math
|HR(x)⟩ = exp(∫ₓᵇ T(z) dz)|HR₀⟩ + ∫ₓᵇ exp(∫ₓᵘ T(z) dz)|hR(u)⟩ du
```

and returned as `HL, converged` with `HL::Piecewise` and where `converged::Bool` is `true`
or `false` depending on convergence up to the requested tolerance.
"""
function righttransfer(HR₀::AbstractMatrix, hR::Union{Nothing, AbstractPiecewise},
                        Ψ₁::FiniteCMPS, Ψ₂::FiniteCMPS = Ψ₁;
                        Kmax = 50, tol = defaulttol(Ψ₁))

    Q₁ = Ψ₁.Q
    R₁s = Ψ₁.Rs
    Q₂ = Ψ₂.Q
    R₂s = Ψ₂.Rs
    grid = nodes(Q₁)
    N = length(grid) - 1
    a = first(grid)
    b = last(grid)
    if Ψ₁ !== Ψ₂
        grid == nodes(Q₂) || throw(DomainMismatch())
    end
    if !isnothing(hR)
        grid == nodes(hR) || throw(DomainMismatch())
    end
    HRᵢ = TaylorSeries([HR₀], b)
    HRs = Vector{typeof(HRᵢ)}(undef, N)
    HRs[N] = HRᵢ
    infoR = true
    for i = N:-1:1
        HRᵢ = HRs[i]
        xᵢ = grid[i]
        xⱼ = grid[i+1]
        Δxᵢ = grid[i+1] - grid[i]
        Q₁ᵢ = shift!(Q₁[i], xⱼ)
        R₁ᵢs = shift!.(getindex.(R₁s, i), xⱼ)
        Q₂ᵢ = shift!(Q₂[i], xⱼ)
        R₂ᵢs = shift!.(getindex.(R₂s, i), xⱼ)
        hRᵢ = isnothing(hR) ? nothing : shift(hR[i], xⱼ)

        # build Taylor coefficients (i.e. solve triangular problem)
        HRᵢ, converged = _righttransfer!(HRᵢ, hRᵢ, Q₁ᵢ, R₁ᵢs, Δxᵢ, Q₂ᵢ, R₂ᵢs;
                                            Kmax = Kmax, tol = tol)
        infoR &= converged
        shift!(HRᵢ, (xᵢ+xⱼ)/2)

        # initialize next ρ element
        if i > 1
            HRs[i-1] = TaylorSeries([HRᵢ(xᵢ)], xᵢ)
        end
    end
    HR = Piecewise(grid, HRs)
    return HR, infoR
end

"""
    _lefttransfer!(ρ::TaylorSeries, h::TaylorSeries,
                    Q₁::TaylorSeries, R₁s::Tuple{Vararg{TaylorSeries}}, Δx,
                    Q₂::TaylorSeries = Q₁, R₂s::Tuple{Vararg{TaylorSeries}} = R₁s;
                    Kmax = 50, tol = defaulttol(ρ))

Solve the non-homogeneous (or homogeneous if h === nothing) Lindblad equation

```math
∂ₓρ = Q₂'*ρ + ρ*Q₁ + R₂'*ρ*R₁ + h
```

in place (i.e. the Taylor coefficients of `ρ` will be inserted), starting from the initial
conditions `ρ[0]` at position `x = offset(ρ)`. The solution will have `degree(ρ) <= Kmax`
and is found up to a tolerance `tol` in the interval `Δx`, i.e. the last few Taylor
coefficients of `ρ` and `h` should satisfy `ρ[k]*Δx^k < tol`. The result is returned as `ρ,
converged`, where `converged::Bool` is `true` or `false` depending on convergence up to the
given tolerance.

Note that the coefficients of `h` are destroyed in the process, if you want to preserve
those you should take a copy beforehand.
"""
function _lefttransfer!(ρ::TaylorSeries, h::Union{Nothing, TaylorSeries},
                            Q₁::TaylorSeries, R₁s::Tuple{Vararg{TaylorSeries}}, Δx,
                            Q₂::TaylorSeries = Q₁, R₂s::Tuple{Vararg{TaylorSeries}} = R₁s;
                            Kmax = 50, tol = defaulttol(ρ))

    Kmin = max(degree(Q₁), degree(Q₂), maximum(degree, R₁s), maximum(degree, R₂s)) + 1
    temp = zero(ρ[0])
    T = eltype(ρ[0])
    for k = 1:Kmax
        ρᵏ = isnothing(h) ? zero(ρ[0]) : h[k-1]
        # solve triangular system:
        for l = 0:min(k-1, degree(Q₁))
            mul!(ρᵏ, ρ[k-1-l], Q₁[l], one(T), one(T))
        end
        for l = 0:min(k-1, degree(Q₂))
            mul!(ρᵏ, Q₂[l]', ρ[k-1-l], one(T), one(T))
        end
        for (R₁, R₂) in zip(R₁s, R₂s)
            for l = 0:min(k-1, degree(R₂))
                for m = 0:min(k-1-l, degree(R₁))
                    mul!(temp, ρ[k-1-l-m], R₁[m], one(T), zero(T))
                    mul!(ρᵏ, R₂[l]', temp, one(T), one(T))
                end
            end
        end
        rmul!(ρᵏ, 1/k)
        ρ[k] = ρᵏ
        # check for convergence
        if k > Kmin
            converged = true
            for l = k-Kmin+1:k
                if norm(ρ[l]) * (Δx^l) > tol
                    converged = false
                end
            end
            if converged
                return ρ, true
            end
        end
    end
    return ρ, false
end

"""
    _righttransfer!(ρ::TaylorSeries, h::TaylorSeries,
                    Q₁::TaylorSeries, R₁s::Tuple{Vararg{TaylorSeries}}, Δx,
                    Q₂::TaylorSeries = Q₁, R₂s::Tuple{Vararg{TaylorSeries}} = R₁s;
                    Kmax = 50, tol = defaulttol(ρ))

Solve the non-homogeneous (or homogeneous if h === nothing) Lindblad equation

```math
∂ₓρ = -(Q₂'*ρ + ρ*Q₁ + R₂'*ρ*R₁) - h
```

in place (i.e. the Taylor coefficients of `ρ` will be inserted), starting from the initial
conditions `ρ[0]` at position `x = offset(ρ)`. The solution will have `degree(ρ) <= Kmax`
and is found up to a tolerance `tol` in the interval `Δx`, i.e. the last few Taylor
coefficients of `ρ` and `h` should satisfy `ρ[k]*Δx^k < tol`. The result is returned as `ρ,
converged`, where `converged::Bool` is `true` or `false` depending on convergence up to the
given tolerance.

Note that the coefficients of `h` are destroyed in the process, if you want to preserve
those you should take a copy beforehand. Also note that all TaylorSeries are assumed to have
the same offset.
"""
function _righttransfer!(ρ::TaylorSeries, h::Union{Nothing, TaylorSeries},
                            Q₁::TaylorSeries, R₁s::Tuple{Vararg{TaylorSeries}}, Δx,
                            Q₂::TaylorSeries = Q₁, R₂s::Tuple{Vararg{TaylorSeries}} = R₁s;
                            Kmax = 50, tol = defaulttol(ρ))

    Kmin = max(degree(Q₁), degree(Q₂), maximum(degree, R₁s), maximum(degree, R₂s)) + 1
    temp = zero(ρ[0])
    T = eltype(ρ[0])
    for k = 1:Kmax
        ρᵏ = isnothing(h) ? zero(ρ[0]) : rmul!(h[k-1], -1)
        # solve triangular system:
        for l = 0:min(k-1, degree(Q₁))
            mul!(ρᵏ, Q₁[l], ρ[k-1-l], -one(T), one(T))
        end
        for l = 0:min(k-1, degree(Q₂))
            mul!(ρᵏ, ρ[k-1-l], Q₂[l]', -one(T), one(T))
        end
        for (R₁, R₂) in zip(R₁s, R₂s)
            for l = 0:min(k-1, degree(R₁))
                for m = 0:min(k-1-l, degree(R₂))
                    mul!(temp, ρ[k-1-l-m], R₂[m]', one(T), zero(T))
                    mul!(ρᵏ, R₁[l], temp, -one(T), one(T))
                end
            end
        end
        rmul!(ρᵏ, 1/k)
        ρ[k] = ρᵏ
        # check for convergence
        if k > Kmin
            converged = true
            for l = k-Kmin+1:k
                if norm(ρ[l]) * (Δx^l) > tol
                    converged = false
                end
            end
            if converged
                return ρ, true
            end
        end
    end
    # not converged:
    return ρ, false
end


text:
"""
This struct contains the collection of all markets modeled
for price forcasts using CEM and endogeneous Economic Dispatch.
"""

struct MarketCollection{Z, T}

    capacity_market::CapacityMarket
    energy_market::EnergyMarket
    reserveup_market::ReserveUpMarket
    reservedown_market::ReserveDownMarket
    rec_market::RECMarket

    function MarketCollection(
                              c::CapacityMarket,
                              e::EnergyMarket{Z, T},
                              r::ReserveUpMarket,
                              l::ReserveDownMarket{Z, T},
                              rec::RECMarket
    ) where {Z, T}
        new{Z, T}(c, e, r, l, rec)
    end
end


text:
import PotentialFlow

function compute_ẋ!(ẋ, x, t)
    plate, ambient_sys = x
    motion = ẋ[1]
    motion.ċ, motion.c̈, motion.α̇ = motion.kin(t)

    Plates.enforce_no_flow_through!(plate, motion, ambient_sys, t)
    reset_velocity!(ẋ, x)
    self_induce_velocity!(ẋ, x, t)
    display(ẋ)
end

function shed_new_vorticity!(blobs, plate, motion, t, lesp = 0.0, tesp = 0.0)
    z₊ = (blobs[end-1].z + 2plate.zs[end])/3
    z₋ = (blobs[end].z + 2plate.zs[1])/3

    blob₊ = PotentialFlow.Vortex.Blob(z₊, 1.0, δ)
    blob₋ = PotentialFlow.Vortex.Blob(z₋, 1.0, δ)
    PotentialFlow.Plates.enforce_no_flow_through!(plate, motion, blobs, t)

    Γ₊, Γ₋, _, _ = PotentialFlow.Plates.vorticity_flux!(plate, blob₊, blob₋, t, lesp, tesp);

    push!(blobs, PotentialFlow.Vortex.Blob(z₊, Γ₊, blobs[1].δ), PotentialFlow.Vortex.Blob(z₋, Γ₋, blobs[1].δ))
end

L = 1.0   # chord length
N = 128   # number of plate control points (distributed along a extrema Chebyshev grid)

ċ = -1.0L  # translation velocity
α = -π/3   # angle of attack

plate = PotentialFlow.Plate(N, L, zero(ComplexF64), α)
motion = PotentialFlow.Plates.RigidBodyMotion(ċ, 0.0);

Δt = 5e-2

δ = 0.01
lesp = 0.0
tesp = 0.0

Δz₀ = im*3Δt*exp(im*plate.α) # vectors perpendicular to the plate
z₋, z₊ = plate.zs[[1,N]] # LE and TE

blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], 1.0, δ) # First two point vortices are placed close to the LE and TE with unit strength

PotentialFlow.Plates.enforce_no_flow_through!(plate, motion, (), 0)
Γ₊, Γ₋, _, _ = PotentialFlow.Plates.vorticity_flux!(plate, blobs[1], blobs[2], 0.0, lesp, tesp); # Determine strength of first two vortices

blobs = PotentialFlow.Vortex.Blob.(Δz₀ .+ [z₊, z₋], [Γ₊, Γ₋], δ) # Create first two point vortices now with calculated strengths

sys₀ = (plate, blobs)

sys = deepcopy(sys₀)
sys₊ = deepcopy(sys₀) # Used for storage during time-marching
ẋs = (motion, PotentialFlow.allocate_velocity(blobs))

forces = ComplexF64[];


text:
import CLIMAParameters

@testset "Override defaults" begin
  struct EarthParameterSet <: AbstractEarthParameterSet end

  ps = EarthParameterSet()
  CLIMAParameters.Planet.grav(::EarthParameterSet) = 2.0
  @test grav(ps) ≈ 2.0

end


text:
# prompt user for input and return it or the default
function user_input(prompt::String, default::String)
    print(prompt)
    response = readline(stdin)
    if response == ""
        return default
    else
        return response
    end
end

# internal function that prompts user to select a deidentification method
function deid_col!(d::OrderedDict, col_nm::String)
    deid_types = ["Nothing", "Hash", "Hash & Salt", "Date Shift", "Drop"]
    col_types = ["", "hash_cols", "salt_cols", "dateshift_cols", "drop_cols"]

    menu = RadioMenu(deid_types, pagesize=5)
    choice = request("Deidentification Method:", menu)

    if choice == 1
        return nothing
    elseif choice == -1
        println("Menu canceled.")
        return nothing
    else
        push!(d[col_types[choice]], col_nm)
        return nothing
    end
end

# initialize dictionary for each dataset - ensures order consistency
get_ds_dict(name::String) = get_ds_dict(name, "")

# initialize dictionary for each dataset - ensures order consistency
function get_ds_dict(name::String, filename::String)
    d = OrderedDict()
    d["name"] = name
    d["filename"] = filename

    col_types = ["rename_cols", "hash_cols", "salt_cols", "dateshift_cols", "drop_cols"]
    for col in col_types
        d[col] = []
    end

    return d
end

# delete unused deid types in the dictionary
function tidy_up!(d)
    col_types = ["rename_cols", "hash_cols", "salt_cols", "dateshift_cols", "drop_cols"]
    for col in col_types
        if d[col] == []
            delete!(d, col)
        end
    end
end

# recursively print yaml file
function print_yaml(io, yml::AbstractArray, indent::Int)
    firstval = true
    for item in yml
        if typeof(item) <: AbstractDict
            print_yaml(io, item, indent)
        else
            write(io, repeat(' ', indent), "- ", string(item), '\n')
        end
    end
end

# recursively print yaml file
function print_yaml(io, yml::AbstractDict, indent::Int)
    firstval = indent > 0
    for (k, v) in yml
        if typeof(v) <: AbstractArray
            write(io, repeat(' ', indent), k, ":", '\n')
            print_yaml(io, v, indent + 2)
        else
            if firstval
                write(io, repeat(' ', indent), "- ")
                write(io, string(k), ": ", string(v), '\n')

                firstval = false
                indent += 2
            else
                write(io, repeat(' ', indent))
                write(io, string(k), ": ", string(v), '\n')
            end
        end
    end
end

"""
    write_yaml(file::String, yml::AbstractDict)

Recursively writes YAML object to file. A YAML object is a dictionary, which can contain
arrays of YAML objects.  See YAML.jl for more on format.
"""
function write_yaml(file::String, yml::AbstractDict)
    open(file, "w") do io
        print_yaml(io, yml, 0)
    end
end

function remove_spaces(str::String; sub::String = "")
    replace(str, " " => sub)
end

"""
    build_config_from_csv(project_name::String, file::String)

Generates a configuration YAML file from a CSV file that defines the mappings.
The CSV file needs to have at least three named columns, one called **Source Table**
which defines the name of the CSV file the data will be read from, a second called
**Field** which defines the name of the field in the data source and a final column
called **Method** which contains the method to apply (one of **Hash - Research ID**, **Hash**,
**Hash & Salt**, **Date Shift**, or **Drop**).

Any column renames and pre- or post-processing will need to be added manually to
the file.
"""
function build_config_from_csv(project_name::String, file::String;
        max_dateshift_days::Int=30, dateshift_years::Int=0, log_path::String="./logs",
        output_path::String="./output", date_format::String="y-m-dTH:M:S.s",
        seed::Union{Int, Nothing}=nothing, config_file::Union{String, Nothing}=nothing)
    if !isfile(file)
        error("$file does not exist!")
    end

    if config_file === nothing
        config_file = "$project_name.yml"
    end

    primary_id = nothing

    yml = OrderedDict()
    yml["project"] = project_name
    yml["seed"] = string(seed === nothing ? make_seed()[1] : seed)
    yml["max_dateshift_days"] = string(max_dateshift_days)
    yml["dateshift_years"] = string(dateshift_years)
    yml["log_path"] = log_path
    yml["output_path"] = output_path
    yml["date_format"] = date_format

    data_sets = OrderedDict()

    for row in CSV.Rows(file)
        # ignore rows without a method
        if ismissing(row.Method)
            continue
        end

        data_set_name = remove_spaces(getproperty(row, Symbol("Source Table")))
        if !haskey(data_sets, data_set_name)
            data_sets[data_set_name] = get_ds_dict(data_set_name)
        end

        data_set = data_sets[data_set_name]
        col_name = remove_spaces(row.Field, sub="_")
        if row.Method == "Hash"
            push!(data_set["hash_cols"], col_name)
        elseif row.Method == "Hash & Salt"
            push!(data_set["salt_cols"], col_name)
        elseif row.Method == "Date Shift"
            push!(data_set["dateshift_cols"], col_name)
        elseif row.Method == "Drop"
            push!(data_set["drop_cols"], col_name)
        elseif row.Method == "Hash - Research ID"
            push!(data_set["hash_cols"], col_name)
            if primary_id === nothing
                yml["primary_id"] = primary_id = col_name
            elseif primary_id != col_name
                error("Primary ID in $data_set_name, $col_name is inconsistent with the rest of the data set")
            end
        else
            @warn("Unknown method $(row.Method) chosen for $data_set_name")
        end
    end

    yml["datasets"] = []
    for p in pairs(data_sets)
        # p.second is the dictionary
        tidy_up!(p.second)
        push!(yml["datasets"], p.second)
    end

    write_yaml(config_file, yml)
    return yml
end

"""
    build_config(data_dir::String, config_file::String)

Interactively guides user through writing a configuration YAML file for DeIdentification.
The `data_dir` should contain one of each type of dataset you expect to deidentify
(e.g. the data directory `./test/data'` contains `pat.csv`, `med.csv`, and `dx.csv`).
The config builder reads the headers of each CSV file and iteratively asks about the
output name and deidentification type of each column. The results are written to `config_file`.
"""
function build_config(data_dir::String, config_file::String)
    if !isfile(config_file)
        touch(config_file)
    end

    if !isdir(data_dir)
        error("data_dir must be a directory containing datasets to be de-identified")
    end

    println("")
    println("DeIdentification Config Builder")
    println("===============================")
    println("Follow the prompts to build a draft of your config file using the datasets.")
    println("The prompts are all written as 'Prompt [default]: '. If there is no default")
    println("the field is required.")
    println("NOTE: this builder will not ask about pre- or post-processing, add after if needed")
    if lowercase(user_input("Ready to get started? [y] ", "y"))[1] != 'y'
        return nothing
    else
        println("Great! Here we go...")
        println()
    end

    yml = OrderedDict()

    println("Let's start with the project level info")
    println("---------------------------------------")

    yml["project"] = user_input("Project name [$(basename(data_dir))]: ", basename(data_dir))
    seed = make_seed()[1]
    yml["project_seed"] = user_input("Project seed [$seed]: (used for reproducibility) ", string(seed))
    yml["max_dateshift_days"] = user_input("Maximum Date Shift Days [30]: ", "30")
    yml["dateshift_years"] = user_input("Years to add to all dates [0]: ", "0")
    yml["log_path"] = user_input("Path for logs [./logs]: ", "./logs")
    yml["output_path"] = user_input("Path for output files [./output]: ", "./output")
    yml["date_format"] = user_input("Input date format [y-m-dTH:M:S.s]: ", "y-m-dTH:M:S.s")

    yml["primary_id"] = ""
    while yml["primary_id"] == ""
        yml["primary_id"] = user_input("Primary ID Column Name: (REQUIRED - must be present in all datasets) ", "")
    end

    println()
    println("Now let's look at the data sets")
    println("-------------------------------")

    yml["datasets"] = []
    for file in Glob.glob(joinpath(data_dir, "*.csv"))
        bnm = basename(file)
        nm = user_input("Dataset Name [$(bnm[1:end-4])]: ", bnm[1:end-4]) # without '.csv'
        fnm = normpath(file)

        d = get_ds_dict(nm, fnm)

        f = CSV.File(fnm, dateformat = yml["date_format"])

        for i in 1:length(f.names)
            println("")
            orig_nm = string(f.names[i])
            println("[  ", orig_nm, " - ", string(f.types[i]), "  ]")

            # rename col?
            col_nm = user_input("Column Name [$(orig_nm)]: ", orig_nm)
            if col_nm != orig_nm
                push!(d["rename_cols"], OrderedDict("in"=>orig_nm, "out"=>col_nm))
            end

            # all others
            deid_col!(d, col_nm)
        end

        tidy_up!(d)
        push!(yml["datasets"], d)
        println("")
    end

    println('\n', "All set! Writing your config file to ", config_file)

    write_yaml(config_file, yml)

    println("Your file is ready - please review it and add any pre- or post-processing steps as needed.")

    return yml

end


text:
using PyCall

py"""
import sys
import subprocess
try:
    import qiskit
except ImportError as e:
    subprocess.check_call([sys.executable, "-m", "pip", "install", "qiskit"])
"""


text:
# https://github.com/sisl/PGFPlots.jl

# significant contributions by: @pkofod

supported_args(::PGFPlotsBackend) = merge_with_base_supported([
    # :annotations,
    # :background_color_legend,
    :background_color_inside,
    # :background_color_outside,
    # :foreground_color_legend, :foreground_color_grid, :foreground_color_axis,
    #     :foreground_color_text, :foreground_color_border,
    :label,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha, :markerstrokestyle,
    :fillrange, :fillcolor, :fillalpha,
    :bins,
    # :bar_width, :bar_edges,
    :title,
    # :window_title,
    :guide, :lims, :ticks, :scale, :flip, :rotation,
    :tickfont, :guidefont, :legendfont,
    :grid, :legend,
    # :colorbar,
    # :marker_z, :levels,
    # :ribbon, :quiver, :arrow,
    # :orientation,
    # :overwrite_figure,
    # :polar,
    # :normalize, :weights, :contours,
    :aspect_ratio,
    # :match_dimensions,
  ])
supported_types(::PGFPlotsBackend) = [:path, :path3d, :scatter, :steppre, :stepmid, :steppost, :histogram2d, :ysticks, :xsticks, :contour]
supported_styles(::PGFPlotsBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
supported_markers(::PGFPlotsBackend) = [:none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5, :pentagon] #vcat(_allMarkers, Shape)
supported_scales(::PGFPlotsBackend) = [:identity, :ln, :log2, :log10]
is_subplot_supported(::PGFPlotsBackend) = false

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

function _initialize_backend(::PGFPlotsBackend; kw...)
    @eval begin
        import PGFPlots
        export PGFPlots
    end
end

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

const _pgfplots_linestyles = KW(
    :solid => "solid",
    :dash => "dashed",
    :dot => "dotted",
    :dashdot => "dashdotted",
    :dashdotdot => "dashdotdotted",
)

const _pgfplots_markers = KW(
    :none => "none",
    :cross => "+",
    :xcross => "x",
    :utriangle => "triangle*",
    :dtriangle => "triangle*",
    :circle => "*",
    :rect => "square*",
    :star5 => "star",
    :star6 => "asterisk",
    :diamond => "diamond*",
    :pentagon => "pentagon*",
)

const _pgfplots_legend_pos = KW(
    :bottomleft => "south west",
    :bottomright => "south east",
    :topright => "north east",
    :topleft => "north west",
)

const _pgf_series_extrastyle = KW(
    :steppre => "const plot mark right",
    :stepmid => "const plot mark mid",
    :steppost => "const plot",
    :sticks => "ycomb",
    :ysticks => "ycomb",
    :xsticks => "xcomb",
)

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

# takes in color,alpha, and returns color and alpha appropriate for pgf style
function pgf_color(c, a = nothing)
    c = getColor(c)
    cstr = @sprintf("{rgb,1:red,%.8f;green,%.8f;blue,%.8f}", red(c), green(c), blue(c))
    a = float(a == nothing ? alpha(c) : a)
    cstr, a
end

function pgf_fillstyle(d::KW)
    cstr,a = pgf_color(d[:fillcolor], d[:fillalpha])
    "fill = $cstr, fill opacity=$a"
end

function pgf_linestyle(d::KW)
    cstr,a = pgf_color(d[:linecolor], d[:linealpha])
    """
    color = $cstr,
    draw opacity=$a,
    line width=$(d[:linewidth]),
    $(get(_pgfplots_linestyles, d[:linestyle], "solid"))"""
end

function pgf_marker(d::KW)
    shape = d[:markershape]
    cstr, a = pgf_color(d[:markercolor], d[:markeralpha])
    cstr_stroke, a_stroke = pgf_color(d[:markerstrokecolor], d[:markerstrokealpha])
    """
    mark = $(get(_pgfplots_markers, shape, "*")),
    mark size = $(0.5 * d[:markersize]),
    mark options = {
        color = $cstr_stroke, draw opacity = $a_stroke,
        fill = $cstr, fill opacity = $a,
        line width = $(d[:markerstrokewidth]),
        rotate = $(shape == :dtriangle ? 180 : 0),
        $(get(_pgfplots_linestyles, d[:markerstrokestyle], "solid"))
    }"""
end

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

function pgf_series(sp::Subplot, series::Series)
    d = series.d
    st = d[:seriestype]
    style = []
    kw = KW()

    push!(style, pgf_linestyle(d))
    push!(style, pgf_marker(d))

    if d[:fillrange] != nothing
        push!(style, pgf_fillstyle(d))
    end

    # add to legend?
    if sp[:legend] != :none && should_add_to_legend(series)
        kw[:legendentry] = d[:label]
    end

    # function args
    args = if st  == :contour
        d[:z].surf, d[:x], d[:y]
    elseif is3d(st)
        d[:x], d[:y], d[:z]
    else
        d[:x], d[:y]
    end

    # PGFPlots can't handle non-Vector?
    args = map(a -> if typeof(a) <: AbstractVector && typeof(a) != Vector
            collect(a)
        else
            a
        end, args)
    # for (i,a) in enumerate(args)
    #     if typeof(a) <: AbstractVector && typeof(a) != Vector
    #         args[i] = collect(a)
    #     end
    # end

    # include additional style, then add to the kw
    if haskey(_pgf_series_extrastyle, st)
        push!(style, _pgf_series_extrastyle[st])
    end
    kw[:style] = join(style, ',')

    # build/return the series object
    func = if st == :path3d
        PGFPlots.Linear3
    elseif st == :scatter
        PGFPlots.Scatter
    elseif st == :histogram2d
        PGFPlots.Histogram2
    elseif st == :contour
        PGFPlots.Contour
    else
        PGFPlots.Linear
    end
    func(args...; kw...)
end

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

function pgf_axis(sp::Subplot, letter)
    axis = sp[Symbol(letter,:axis)]
    style = []
    kw = KW()

    # axis guide
    kw[Symbol(letter,:label)] = axis[:guide]

    # flip/reverse?
    axis[:flip] && push!(style, "$letter dir=reverse")

    # scale
    scale = axis[:scale]
    if scale in (:log2, :ln, :log10)
        kw[Symbol(letter,:mode)] = "log"
        scale == :ln || push!(style, "log basis $letter=$(scale == :log2 ? 2 : 10)")
    end

    # limits
    # TODO: support zlims
    if letter != :z
        lims = axis_limits(axis)
        kw[Symbol(letter,:min)] = lims[1]
        kw[Symbol(letter,:max)] = lims[2]
    end

    # return the style list and KW args
    style, kw
end

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

function _make_pgf_plot!(plt::Plot)
    plt.o = PGFPlots.Axis[]
    for sp in plt.subplots
        # first build the PGFPlots.Axis object
        style = ["unbounded coords=jump"]
        kw = KW()

        # add to style/kw for each axis
        for letter in (:x, :y, :z)
            if letter != :z || is3d(sp)
                axisstyle, axiskw = pgf_axis(sp, letter)
                merge!(kw, axiskw)
            end
        end

        # bounding box values are in mm
        # note: bb origin is top-left, pgf is bottom-left
        bb = bbox(sp)
        push!(style, """
            xshift = $(left(bb).value)mm,
            yshift = $((height(bb) - (bottom(bb))).value)mm,
            width = $(width(bb).value)mm,
            height = $(height(bb).value)mm,
            axis background/.style={fill=$(pgf_color(sp[:background_color_inside])[1])}
        """)

        if sp[:title] != ""
            push!(style, "title = $(sp[:title])")
        end

        sp[:grid] && push!(style, "grid = major")
        if sp[:aspect_ratio] in (1, :equal)
            kw[:axisEqual] = "true"
        end

        legpos = sp[:legend]
        if haskey(_pgfplots_legend_pos, legpos)
            kw[:legendPos] = _pgfplots_legend_pos[legpos]
        end

        o = PGFPlots.Axis(; style = style, kw...)

        # add the series object to the PGFPlots.Axis
        for series in series_list(sp)
            push!(o, pgf_series(sp, series))
        end

        # add the PGFPlots.Axis to the list
        push!(plt.o, o)
    end
end

function _writemime(io::IO, mime::MIME"image/svg+xml", plt::Plot{PGFPlotsBackend})
    _make_pgf_plot!(plt)
    writemime(io, mime, plt.o)
end

function _writemime(io::IO, mime::MIME"application/pdf", plt::Plot{PGFPlotsBackend})
    _make_pgf_plot!(plt)

    # prepare the object
    pgfplt = PGFPlots.plot(plt.o)

    # save a pdf
    fn = tempname()*".pdf"
    PGFPlots.save(PGFPlots.PDF(fn), pgfplt)

    # read it into io
    write(io, readall(open(fn)))

    # cleanup
    PGFPlots.cleanup(plt.o)
end

function _display(plt::Plot{PGFPlotsBackend})
    # prepare the object
    _make_pgf_plot!(plt)
    pgfplt = PGFPlots.plot(plt.o)

    # save an svg
    fn = string(tempname(), ".svg")
    PGFPlots.save(PGFPlots.SVG(fn), pgfplt)

    # show it
    open_browser_window(fn)

    # cleanup
    PGFPlots.cleanup(plt.o)
end


text:
using JuLIP, Test, Printf
using JuLIP.Testing

h0("   JuLIP Tests   ")

@info("preparing the tests...")

verbose=true
isCI = haskey(ENV, "CI")
notCI = !isCI
eam_W4 = nothing

## ===== some prototype potentials ======
@info("Loading some interatomic potentials . .")
data = joinpath(dirname(pathof(JuLIP)), "..", "data") * "/"
eam_Fe = JuLIP.Potentials.EAM(data * "pfe.plt", data * "ffe.plt", data * "F_fe.plt")
print(" .")
eam_W = JuLIP.Potentials.FinnisSinclair(data*"W-pair-Wang-2014.plt", data*"W-e-dens-Wang-2014.plt")
print(" .")
global eam_W4
try
   global eam_W4 = JuLIP.Potentials.EAM(data * "w_eam4.fs")
catch
   global eam_W4 = nothing
end
println(" done.")

##

julip_tests = [
   ("testaux.jl", "Miscellaneous"),
   ("test_atoms.jl", "Atoms"),
   ("test_build.jl", "Build"),
   ("test_fio.jl", "File IO"),
   ("testanalyticpotential.jl", "Analytic Potential"),
   ("testpotentials.jl", "Potentials"),
   ("test_ad.jl", "AD Potentials"),
   ("testvarcell.jl", "Variable Cell"),
   ("testhessian.jl", "Hessian"),
   ("test_onebody.jl", "One-Body"),
   ]

# add solver tests if not on travis
if !isCI
   push!(julip_tests, ("testsolve.jl", "Solve"))
else
   @info("on CI : don't run solver tests")
end

# "testexpvarcell.jl";  # USE THIS TO WORK ON EXPCELL IMPLEMENTATION

@testset "JuLIP" begin
   for (testfile, testid) in julip_tests
      h1("Testset $(testid)")
      @testset "$(testid)" begin include(testfile); end
   end
end

# TODO:
# - fix performance regression in morseold / testanalyticpotential
# - some other analytic EAM potentials
# - finish proper EAM implementation
#   do this via ASE.jl and importing the data from ASE!!
# - maybe simplify the hessian assembly by skipping the step via
#   hessian_pos


text:
using Test, Astro

print("47.a Moon position")

# Calculate the geocentric longitude, latitude, distance, and equatorial horizontal parallax
# of the moon for 1992 April 12 at 0h TD

L, B, R = moon_dimension3(cal_to_jd(1992, 4, 12))

# longitude should be 133d.152655
# beta -3d.229126
# delta 368409.7 km
# pi 0.991990

@test isapprox(rad2deg(L), 133.162, atol = .1) # 0.1 degree OK?
@test isapprox(rad2deg(B), -3.229126, atol = 1e-6) # degrees
@test isapprox(R, 368409.7, atol = 0.1) # km

# the other routine:

L = moon_dimension(2448724.5, "L")
@test isapprox(rad2deg(L), 133.162655, atol=0.1)

B = moon_dimension(2448724.5, "B")
@test isapprox(rad2deg(B), -3.229126, atol= 0.1) # degrees

R = moon_dimension(2448724.5, "R")
@test isapprox(R, 368409.7, atol= 0.1) # km

println(" passed")

print("47 Moon ascending node longitude")
# times when the mean ascending or descending node of the lunar orbit coincides with the vernal equinox,
# ie when it is equal to 0 or 180.
# It is 180 on 2015 Oct 10, and 0 on 2025 Jan 29
L = rad2deg(moon_mean_ascending_node_longitude(cal_to_jd(2015, 10, 10)))
@test isapprox(L, 180, atol=10e-2)
L = rad2deg(moon_mean_ascending_node_longitude(cal_to_jd(2025, 1, 29)))
@test isapprox(L, 0, atol=10e-2)
println(" passed")

print("48.a (1) Illuminated fraction of moon's disk (high precision)")
k, khi = moon_illuminated_fraction_high(cal_to_jd(1992, 4, 12))
@test isapprox(k, 0.678, atol= 0.001) # percent
@test isapprox(rad2deg(khi), 285, atol= 0.1) # degrees
println(" passed")

print("48.a (2) Illuminated fraction of moon's disk (low precision)")
k = moon_illuminated_fraction_low(cal_to_jd(1992, 4, 12))
@test isapprox(k, 0.680, atol=0.001)
println(" passed")

print("49.a New moon")
jd = moon_nearest_phase(cal_to_jd(1977, 2, 13), 0)
yr, mo, dy = jd_to_cal(jd)
@test isapprox(yr, 1977, atol=1)
@test isapprox(mo, 2, atol=0)
@test isapprox(dy, 18.15, atol=0.01)
println(" passed")

print("49.b Moon last quarter")
jd = moon_nearest_phase(cal_to_jd(2044, 1, 19), 3)
yr, mo, dy = jd_to_cal(jd)
@test isapprox(yr, 2044, atol=0.1)
@test isapprox(mo, 1, atol=0.1)
@test isapprox(dy, 21.99, atol=0.01)
println(" passed")

print("50 Moon apogee")
jd, parallax = moon_apogee_perigee_time_low(cal_to_jd(1988, 10, 1), 1)
# Apogee time
@test jd * seconds_per_day/60 ≈  2447442.3537 * seconds_per_day / 60 atol= 1# minutes
# parallax
@test parallax ≈  3240.679 atol= 1# seconds
println(" passed")

print("50 Moon perigee (p 361)")
jd, para = moon_apogee_perigee_time_low(cal_to_jd(1997, 12, 9), 0)
# Perigee time
@test isapprox(jd, 2450792.2059, atol=0.1)
jd, para = moon_apogee_perigee_time_low(cal_to_jd(1998, 1, 3), 0)
# Perigee time"
@test isapprox(jd, 2450816.8549, atol=0.1)
jd, para = moon_apogee_perigee_time_low(cal_to_jd(1990, 12, 2), 0)
# Perigee time
@test isapprox(jd, 2448227.9505, atol=0.1)
jd, para = moon_apogee_perigee_time_low(cal_to_jd(1990, 12, 30), 0)
# Perigee time
@test isapprox(jd, 2448256.4941, atol=0.1)
println(" passed")

print("51.a Moon ascending node time")
# node()
#    	- jd: julian date
#    	- desc_not_asc: 0 for the ascending node, 1 for the descending node
# returns jd of the closer passage
jd = moon_node(cal_to_jd(1987, 5, 15), 0)
@test isapprox(jd, 2446938.76, atol=0.01)
println(" passed")

print("Lunation")
jd = cal_to_jd(2014,10,14)
lunationnumber = lunation(jd, "brown")
@test isapprox(lunationnumber, 1135.0, atol=.1)
println(" passed")

print("Nearest moon phase to a specific day")
# NASA, at http://eclipse.gsfc.nasa.gov/phase/phase2001gmt.html,
# give New: Oct 23  21:57  FirstQ: Oct 31  02:48  Full:Nov  6  22:23  Last Quarter: Nov 14  15:16
jd = cal_to_jd(2014, 10, 15)
nearest_new = moon_nearest_phase(jd, 0)
nearest_first = moon_nearest_phase(jd, 1)
nearest_full  = moon_nearest_phase(jd + 10, 2) # get next, not nearest
nearest_last  = moon_nearest_phase(jd + 10, 3)

@test isapprox(nearest_new,   date_to_jd(2014,10,23,21,57, 0), atol=.001) # within about 2 minutes
@test isapprox(nearest_first, date_to_jd(2014,10,31,2,48, 0), atol=.001)
@test isapprox(nearest_full,  date_to_jd(2014,11,6,22,23, 0), atol=.001)
@test isapprox(nearest_last,  date_to_jd(2014,11,14,15,16, 0), atol=.001)
println(" passed")

print("Moon position")
# Duffett-Smith: what was the position of the Moon on Feb 26 1979 at 16h 00 m UT
# Should be ra 22 33 29s decl -8 02 42s
jd = date_to_jd(1979,2,26,16,0,0)
geoecl_long, geoecl_lat, rad = moon_dimension3(jd)
ra_rad, dec_rad = ecl_to_equ(geoecl_long, geoecl_lat, obliquity(jd))
ra_deg = ra_rad / (2 * pi) |> fday_to_hms
@test(ra_deg[1] == 22)
@test(ra_deg[2] == 33)
@test isapprox(ra_deg[3], 28, atol = 30) # within 30 seconds OK?
#
# what about testing dec_rad?
#
# long and lat do the same job as moon_dimension3
@test isapprox(moon_longitude(jd), moon_dimension3(jd)[1], atol=0.01)
@test isapprox(moon_latitude(jd), moon_dimension3(jd)[2], atol=0.01)
println(" passed")

println("Moon rise and moon set")
# Duffett-Smith what were the times of moonrise and moonset on Sept 6 1979
# as observed from longitude 0N and latitude 52N
# ephemeris says 18h46m and 04h58m
set_latitude(deg2rad(52))
set_longitude(deg2rad(0))
jd = cal_to_jd(1979,9,6)
r,s = moon_riseset(jd)
d1, d2 = map(jd_to_date, (r,s))
hr, mn, sc = fday_to_hms(r - jd)
@test (hr == 18)
@test isapprox(mn, 46, atol = 30) # TODO moon_riseset test within 30 minutes???

hr, mn, sc = fday_to_hms(s - jd)
@test isapprox(hr, 4, atol = 1) #
@test isapprox(mn, 58, atol = 300) # within 240 minutes (4 hours) ???

println()

"""
 TODO moon-tests.jl more still to be tested:
moon_horizontal_parallax(jd::Float64)
moon_mean_ascending_node_longitude(jd)
moon_mean_perigee_longitude(jd::Float64)
moon_node(jd::Float64, desc_not_asc)
moon_radius(jd::Float64)
moon_true_ascending_node_longitude(jd::Float64)

"""