text:
#=
This file contains all of the common interface functions.
Many of them simply call corresponding functions in jl.
=#
export generateFor, useLog, domain, solverType, functionSpace, trialSpace, testSpace, finiteVolumeOrder,
        nodeType, timeStepper, setSteps, matrixFree, customOperator, customOperatorFile,
        mesh, exportMesh, variable, coefficient, parameter, testSymbol, index, boundary, addBoundaryID,
        referencePoint, timeInterval, initial, weakForm, fluxAndSource, flux, source, assemblyLoops,
        exportCode, importCode, printLatex,
        solve, cachesimSolve, finalize_femshop, cachesim, output_values,
        morton_nodes, hilbert_nodes, tiled_nodes, morton_elements, hilbert_elements,
        tiled_elements, ef_nodes, random_nodes, random_elements

# Begin configuration setting functions

function generateFor(lang; filename=project_name, header="", params=nothing)
    outputDirPath = pwd()*"/"*uppercasefirst(filename);
    if !isdir(outputDirPath)
        mkdir(outputDirPath);
    end
    framew = 0;
    if !in(lang, [CPP,MATLAB,DENDRO,HOMG])
        # lang should be a filename for a custom target
        # This file must include these three functions:
        # 1. get_external_language_elements() - file extensions, comment chars etc.
        # 2. generate_external_code_layer(var, entities, terms, lorr, vors) - Turns symbolic expressions into code
        # 3. generate_external_files(var, lhs_vol, lhs_surf, rhs_vol, rhs_surf) - Writes all files based on generated code
        include(lang);
        set_custom_gen_target(get_external_language_elements, generate_external_code_layer, generate_external_files, outputDirPath, filename, head=header);
    else
        if lang == DENDRO
            framew = DENDRO;
            lang = CPP;
        elseif lang == HOMG
            framew = HOMG;
            lang = MATLAB;
        else
            framew = 0;
        end
        set_included_gen_target(lang, framew, outputDirPath, filename, head=header);
    end
    if !(params === nothing)
        set_codegen_parameters(params);
    end
end

function useLog(name=project_name; dir=output_dir, level=2)
    init_log(name, dir, level);
end

function domain(dims; shape=SQUARE, grid=UNIFORM_GRID)
    config.dimension = dims;
    config.geometry = shape;
    config.mesh_type = grid;
end

function solverType(type)
    set_solver(type);
end

function functionSpace(;space=LEGENDRE, order=0, orderMin=0, orderMax=0)
    config.trial_function = space;
    config.test_function = space;
    if orderMax > orderMin && orderMin >= 0
        config.p_adaptive = true;
        config.basis_order_min = orderMin;
        config.basis_order_max = orderMax;
    else
        config.basis_order_min = max(order, orderMin);
        config.basis_order_max = max(order, orderMin);
    end
end

function trialSpace(;space=LEGENDRE, order=0, orderMin=0, orderMax=0)
    #TODO
    functionSpace(space=space, order=order, orderMin=orderMin, orderMax=orderMax);
end

function testSpace(;space=LEGENDRE, order=0, orderMin=0, orderMax=0)
    #TODO
    functionSpace(space=space, order=order, orderMin=orderMin, orderMax=orderMax);
end

function nodeType(type)
    config.elemental_nodes = type;
end

function timeStepper(type; cfl=0)
    set_stepper(type, cfl);
end

function setSteps(dt, steps)
    set_specified_steps(dt, steps);
end

function matrixFree(shallwe=true; maxiters=100, tol=1e-6)
    config.linalg_matrixfree = shallwe;
    config.linalg_matfree_max = maxiters;
    config.linalg_matfree_tol = tol;
end

function customOperator(name, handle)
    s = Symbol(name);
    add_custom_op(s, handle);
end

function customOperatorFile(filename)
    log_entry("Adding custom operators from file: "*string(filename), 2);
    add_custom_op_file(filename);
end

# End configuration functions, begin problem definition functions

function mesh(msh; elsperdim=5, bids=1, interval=[0,1])
    if msh == LINEMESH
        log_entry("Building simple line mesh with nx elements, nx="*string(elsperdim));
        meshtime = @elapsed(add_mesh(simple_line_mesh(elsperdim.+1, bids, interval)));
        log_entry("Grid building took "*string(meshtime)*" seconds");

    elseif msh == QUADMESH
        log_entry("Building simple quad mesh with nx*nx elements, nx="*string(elsperdim));
        meshtime = @elapsed(add_mesh(simple_quad_mesh(elsperdim.+1, bids, interval)));
        log_entry("Grid building took "*string(meshtime)*" seconds");

    elseif msh == HEXMESH
        log_entry("Building simple hex mesh with nx*nx*nx elements, nx="*string(elsperdim));
        meshtime = @elapsed(add_mesh(simple_hex_mesh(elsperdim.+1, bids, interval)));
        log_entry("Grid building took "*string(meshtime)*" seconds");

    else # msh should be a mesh file name
        # open the file and read the mesh data
        mfile = open(msh, "r");
        log_entry("Reading mesh file: "*msh);
        meshtime = @elapsed(mshdat=read_mesh(mfile));
        log_entry("Mesh reading took "*string(meshtime)*" seconds");
        add_mesh(mshdat);
        close(mfile);
    end
end

function exportMesh(filename, format=MSH_V2)
    # open the file to write to
    mfile = open(filename, "w");
    log_entry("Writing mesh file: "*filename);
    output_mesh(mfile, format);
    close(mfile);
end

function finiteVolumeOrder(order)
    # Not sure where to put these yet. Let's hand them over to FVSolver.
    (parent, child) = divide_parent_grid(grid_data, order);
    set_parent_and_child(parent, child);
end

function variable(name, type=SCALAR, location=NODAL; index=nothing)
    varind = var_count + 1;
    varsym = Symbol(name);
    # Just make an empty variable with the info known so far.
    var = Variable(varsym, [], varind, type, location, [], index, 0, [], false);
    add_variable(var);
    return var;
end

function coefficient(name, val, type=SCALAR, location=NODAL; index=nothing)
    csym = Symbol(name);
    nfuns = makeFunctions(val); # if val is constant, nfuns will be 0
    return add_coefficient(csym, type, location, val, nfuns);
end

function parameter(name, val, type=SCALAR)
    if length(parameters) == 0
        coefficient("parameterCoefficientForx", "x")
        coefficient("parameterCoefficientFory", "y")
        coefficient("parameterCoefficientForz", "z")
        coefficient("parameterCoefficientFort", "t")
    end
    if typeof(val) <: Number
        newval = [val];
    elseif typeof(val) == String
        # newval will be an array of expressions
        # search for x,y,z,t symbols and replace with special coefficients like parameterCoefficientForx
        newval = [swap_parameter_xyzt(Meta.parse(val))];

    elseif typeof(val) <: Array
        newval = Array{Expr,1}(undef,length(val));
        for i=1:length(val)
            newval[i] = swap_parameter_xyzt(Meta.parse(val[i]));
        end
        newval = reshape(newval,size(val));
    else
        println("Error: use strings to define parameters");
        newval = 0;
    end

    return add_parameter(Symbol(name), type, newval);
end

function testSymbol(symb, type=SCALAR)
    add_test_function(Symbol(symb), type);
end

function index(name; range=[1])
    if length(range) == 2
        range = Array(range[1]:range[2]);
    end
    idx = Indexer(Symbol(name), range, range[1])
    add_indexer(idx);
    return idx;
end

function boundary(var, bid, bc_type, bc_exp=0)
    # The expression may contain variable symbols.
    # Parse it, replace variable symbols with the appropriate parts
    if typeof(bc_exp) <: Array
        for i=1:length(bc_exp)
            if typeof(bc_exp[i]) == String
                ex = Meta.parse(bc_exp[i]);
                ex = replace_var_symbols_with_values(ex);
                bc_exp[i] = string(ex);
            end
        end
    else
        if typeof(bc_exp) == String
            ex = Meta.parse(bc_exp);
            ex = replace_var_symbols_with_values(ex);
            bc_exp = string(ex);
        end
    end
    nfuns = makeFunctions(bc_exp);
    add_boundary_condition(var, bid, bc_type, bc_exp, nfuns);
end

function addBoundaryID(bid, trueOnBdry)
    # trueOnBdry(x, y, z) = something # points with x,y,z on this bdry segment evaluate true here
    if typeof(trueOnBdry) == String
        trueOnBdry = stringToFunction("trueOnBdry", "x,y=0,z=0", trueOnBdry);
    end
    add_boundary_ID_to_grid(bid, trueOnBdry, grid_data);
end

function referencePoint(var, pos, val)
    add_reference_point(var, pos, val);
end

function timeInterval(T)
    prob.time_dependent = true;
    if time_stepper === nothing
        timeStepper(EULER_IMPLICIT);
    end
    prob.end_time = T;
end

function initial(var, ics)
    nfuns = makeFunctions(ics);
    add_initial_condition(var.index, ics, nfuns);
end

function weakForm(var, wf)
    if typeof(var) <: Array
        # multiple simultaneous variables
        wfvars = [];
        wfex = [];
        if !(length(var) == length(wf))
            printerr("Error in weak form: # of unknowns must equal # of equations. (example: @weakform([a,b,c], [f1,f2,f3]))");
        end
        for vi=1:length(var)
            push!(wfvars, var[vi].symbol);
            push!(wfex, Meta.parse((wf)[vi]));
        end
    else
        wfex = Meta.parse(wf);
        wfvars = var.symbol;
    end

    log_entry("Making weak form for variable(s): "*string(wfvars));
    log_entry("Weak form, input: "*string(wf));

    # This is the parsing step. It goes from an Expr to arrays of Basic
    result_exprs = sp_parse(wfex, wfvars);
    if length(result_exprs) == 4 # has surface terms
        (lhs_expr, rhs_expr, lhs_surf_expr, rhs_surf_expr) = result_exprs;
    else
        (lhs_expr, rhs_expr) = result_exprs;
    end

    if typeof(lhs_expr) <: Tuple && length(lhs_expr) == 2 # has time derivative
        if length(result_exprs) == 4 # has surface terms
            log_entry("Weak form, before modifying for time: Dt("*string(lhs_expr[1])*") + "*string(lhs_expr[2])*" + surface("*string(lhs_surf_expr)*") = "*string(rhs_expr)*" + surface("*string(rhs_surf_expr)*")");

            (newlhs, newrhs, newsurflhs, newsurfrhs) = reformat_for_stepper(lhs_expr, rhs_expr, lhs_surf_expr, rhs_surf_expr, config.stepper);

            log_entry("Weak form, modified for time stepping: "*string(newlhs)*" + surface("*string(newsurflhs)*") = "*string(newrhs)*" + surface("*string(newsurfrhs)*")");

            lhs_expr = newlhs;
            rhs_expr = newrhs;
            lhs_surf_expr = newsurflhs;
            rhs_surf_expr = newsurfrhs;

        else # time derivative, but no surface
            log_entry("Weak form, before modifying for time: Dt("*string(lhs_expr[1])*") + "*string(lhs_expr[2])*" = "*string(rhs_expr));

            (newlhs, newrhs) = reformat_for_stepper(lhs_expr, rhs_expr, config.stepper);

            log_entry("Weak form, modified for time stepping: "*string(newlhs)*" = "*string(newrhs));

            lhs_expr = newlhs;
            rhs_expr = newrhs;
        end

    end

    # Here we build a SymExpression for each of the pieces.
    # This is an Expr tree that is passed to the code generator.
    if length(result_exprs) == 4 # has surface terms
        (lhs_symexpr, rhs_symexpr, lhs_surf_symexpr, rhs_surf_symexpr) = build_symexpressions(wfvars, lhs_expr, rhs_expr, lhs_surf_expr, rhs_surf_expr, remove_zeros=true);
        set_symexpressions(var, lhs_symexpr, LHS, "volume");
        set_symexpressions(var, lhs_surf_symexpr, LHS, "surface");
        set_symexpressions(var, rhs_symexpr, RHS, "volume");
        set_symexpressions(var, rhs_surf_symexpr, RHS, "surface");

        log_entry("lhs volume symexpression:\n\t"*string(lhs_symexpr));
        log_entry("lhs surface symexpression:\n\t"*string(lhs_surf_symexpr));
        log_entry("rhs volume symexpression:\n\t"*string(rhs_symexpr));
        log_entry("rhs surface symexpression:\n\t"*string(rhs_surf_symexpr));

        log_entry("Latex equation:\n\t\t \\int_{K}"*symexpression_to_latex(lhs_symexpr)*
                    " dx + \\int_{\\partial K}"*symexpression_to_latex(lhs_surf_symexpr)*
                    " ds = \\int_{K}"*symexpression_to_latex(rhs_symexpr)*
                    " dx + \\int_{\\partial K}"*symexpression_to_latex(rhs_surf_symexpr)*" ds");
    else
        (lhs_symexpr, rhs_symexpr) = build_symexpressions(wfvars, lhs_expr, rhs_expr, remove_zeros=true);
        set_symexpressions(var, lhs_symexpr, LHS, "volume");
        set_symexpressions(var, rhs_symexpr, RHS, "volume");

        log_entry("lhs symexpression:\n\t"*string(lhs_symexpr));
        log_entry("rhs symexpression:\n\t"*string(rhs_symexpr));

        log_entry("Latex equation:\n\t\t\$ \\int_{K}"*symexpression_to_latex(lhs_symexpr)*
                    " dx = \\int_{K}"*symexpression_to_latex(rhs_symexpr)*" dx");
    end

    # ########## This part simply makes a string for printing #############
    # # make a string for the expression
    # if typeof(lhs_expr[1]) <: Array
    #     lhsstring = "";
    #     rhsstring = "";
    #     for i=1:length(lhs_expr)
    #         lhsstring = lhsstring*"lhs"*string(i)*" = "*string(lhs_expr[i][1]);
    #         rhsstring = rhsstring*"rhs"*string(i)*" = "*string(rhs_expr[i][1]);
    #         for j=2:length(lhs_expr[i])
    #             lhsstring = lhsstring*" + "*string(lhs_expr[i][j]);
    #         end
    #         for j=2:length(rhs_expr[i])
    #             rhsstring = rhsstring*" + "*string(rhs_expr[i][j]);
    #         end
    #         if length(result_exprs) == 4
    #             for j=1:length(lhs_surf_expr)
    #                 lhsstring = lhsstring*" + surface("*string(lhs_surf_expr[j])*")";
    #             end
    #             for j=1:length(rhs_surf_expr)
    #                 rhsstring = rhsstring*" + surface("*string(rhs_surf_expr[j])*")";
    #             end
    #         end
    #         lhsstring = lhsstring*"\n";
    #         rhsstring = rhsstring*"\n";
    #     end
    # else
    #     lhsstring = "lhs = "*string(lhs_expr[1]);
    #     rhsstring = "rhs = "*string(rhs_expr[1]);
    #     for j=2:length(lhs_expr)
    #         lhsstring = lhsstring*" + "*string(lhs_expr[j]);
    #     end
    #     for j=2:length(rhs_expr)
    #         rhsstring = rhsstring*" + "*string(rhs_expr[j]);
    #     end
    #     if length(result_exprs) == 4
    #         for j=1:length(lhs_surf_expr)
    #             lhsstring = lhsstring*" + surface("*string(lhs_surf_expr[j])*")";
    #         end
    #         for j=1:length(rhs_surf_expr)
    #             rhsstring = rhsstring*" + surface("*string(rhs_surf_expr[j])*")";
    #         end
    #     end
    # end
    # log_entry("Weak form, symbolic layer:\n\t"*string(lhsstring)*"\n"*string(rhsstring));
    # ################ string ####################################

    # change symbolic layer into code layer
    (lhs_string, lhs_code) = generate_code_layer(lhs_symexpr, var, LHS, "volume", config.solver_type, language, gen_framework);
    (rhs_string, rhs_code) = generate_code_layer(rhs_symexpr, var, RHS, "volume", config.solver_type, language, gen_framework);
    if length(result_exprs) == 4
        (lhs_surf_string, lhs_surf_code) = generate_code_layer(lhs_surf_symexpr, var, LHS, "surface", config.solver_type, language, gen_framework);
        (rhs_surf_string, rhs_surf_code) = generate_code_layer(rhs_surf_symexpr, var, RHS, "surface", config.solver_type, language, gen_framework);
        log_entry("Weak form, code layer: LHS = \n"*string(lhs_string)*"\nsurfaceLHS = \n"*string(lhs_surf_string)*" \nRHS = \n"*string(rhs_string)*"\nsurfaceRHS = \n"*string(rhs_surf_string));
    else
        log_entry("Weak form, code layer: LHS = \n"*string(lhs_string)*" \n  RHS = \n"*string(rhs_string));
    end

    #log_entry("Julia code Expr: LHS = \n"*string(lhs_code)*" \n  RHS = \n"*string(rhs_code), 3);

    if language == JULIA || language == 0
        args = "args; kwargs...";
        makeFunction(args, lhs_code);
        set_lhs(var, lhs_string);

        makeFunction(args, rhs_code);
        set_rhs(var, rhs_string);

        if length(result_exprs) == 4
            makeFunction(args, string(lhs_surf_code));
            set_lhs_surface(var, lhs_surf_string);

            makeFunction(args, string(rhs_surf_code));
            set_rhs_surface(var, rhs_surf_string);
        end

    else
        set_lhs(var, lhs_code);
        set_rhs(var, rhs_code);
    end
end

function fluxAndSource(var, fex, sex)
    flux(var, fex);
    source(var, sex);
end

function flux(var, fex)
    if typeof(var) <: Array
        # multiple simultaneous variables
        symvars = [];
        symfex = [];
        if !(length(var) == length(fex))
            printerr("Error in flux function: # of unknowns must equal # of equations. (example: flux([a,b,c], [f1,f2,f3]))");
        end
        for vi=1:length(var)
            push!(symvars, var[vi].symbol);
            push!(symfex, Meta.parse((fex)[vi]));
        end
    else
        symfex = Meta.parse(fex);
        symvars = var.symbol;
    end

    log_entry("Making flux for variable(s): "*string(symvars));
    log_entry("flux, input: "*string(fex));

    # The parsing step
    (flhs_expr, frhs_expr) = sp_parse(symfex, symvars);
    # Note that this automatically puts a (-) on all rhs terms. These need to be reversed for FV.

    # Modify the expressions according to time stepper.
    # There is an assumed Dt(u) added which is the only time derivative.
    log_entry("flux, before modifying for time: "*string(flhs_expr)*" - "*string(frhs_expr));

    (newflhs, newfrhs) = reformat_for_stepper_fv_flux(flhs_expr, frhs_expr, config.stepper);

    log_entry("flux, modified for time stepping: "*string(newflhs)*" + "*string(newfrhs));

    flhs_expr = newflhs;
    frhs_expr = newfrhs;

    # Here we build a SymExpression for each of the pieces.
    # This is passed to the code generator.
    (lhs_symexpr, rhs_symexpr) = build_symexpressions(symvars, flhs_expr, frhs_expr, remove_zeros=true);
    set_symexpressions(var, lhs_symexpr, LHS, "surface");
    set_symexpressions(var, rhs_symexpr, RHS, "surface");

    log_entry("flux lhs symexpression:\n\t"*string(lhs_symexpr));
    log_entry("flux rhs symexpression:\n\t"*string(rhs_symexpr));
    log_entry("Latex flux equation:\n\t\t \\int_{\\partial K}"*symexpression_to_latex(lhs_symexpr)*
                    " ds - \\int_{\\partial K}"*symexpression_to_latex(rhs_symexpr)*" ds");

    # change symbolic layer into code layer
    (lhs_string, lhs_code) = generate_code_layer(lhs_symexpr, var, LHS, "surface", FV, language, gen_framework);
    (rhs_string, rhs_code) = generate_code_layer(rhs_symexpr, var, RHS, "surface", FV, language, gen_framework);
    log_entry("flux, code layer: \n  LHS = "*string(lhs_string)*" \n  RHS = "*string(rhs_string));

    if language == JULIA || language == 0
        args = "args; kwargs...";
        makeFunction(args, string(lhs_code));
        set_lhs_surface(var, lhs_string);

        makeFunction(args, string(rhs_code));
        set_rhs_surface(var, rhs_string);

    else
        set_lhs_surface(var, lhs_code);
        set_rhs_surface(var, rhs_code);
    end
end

function source(var, sex)
    if typeof(var) <: Array
        # multiple simultaneous variables
        symvars = [];
        symsex = [];
        if !(length(var) == length(sex))
            printerr("Error in source function: # of unknowns must equal # of equations. (example: source([a,b,c], [s1,s2,s3]))");
        end
        for vi=1:length(var)
            push!(symvars, var[vi].symbol);
            push!(symsex, Meta.parse((sex)[vi]));
        end
    else
        symsex = Meta.parse(sex);
        symvars = var.symbol;
    end

    log_entry("Making source for variable(s): "*string(symvars));
    log_entry("source, input: "*string(sex));

    # The parsing step
    (slhs_expr, srhs_expr) = sp_parse(symsex, symvars);
    # Note that this automatically puts a (-) on all rhs terms. These need to be reversed for FV.

    # Modify the expressions according to time stepper.
    # There is an assumed Dt(u) added which is the only time derivative.
    log_entry("source, before modifying for time: "*string(slhs_expr)*" - "*string(srhs_expr));

    (newslhs, newsrhs) = reformat_for_stepper_fv_source(slhs_expr, srhs_expr, config.stepper);

    log_entry("source, modified for time stepping: "*string(newslhs)*" + "*string(newsrhs));

    slhs_expr = newslhs;
    srhs_expr = newsrhs;

    # Here we build a SymExpression for each of the pieces.
    # This is passed to the code generator.
    (lhs_symexpr, rhs_symexpr) = build_symexpressions(symvars, slhs_expr, srhs_expr, remove_zeros=true);
    set_symexpressions(var, lhs_symexpr, LHS, "volume");
    set_symexpressions(var, rhs_symexpr, RHS, "volume");

    log_entry("source lhs symexpression:\n\t"*string(lhs_symexpr));
    log_entry("source rhs symexpression:\n\t"*string(rhs_symexpr));
    log_entry("Latex source equation:\n\t\t \\int_{K} -"*symexpression_to_latex(lhs_symexpr)*
                    " dx + \\int_{K}"*symexpression_to_latex(rhs_symexpr)*" dx");

    # change symbolic layer into code code_layer
    (lhs_string, lhs_code) = generate_code_layer(lhs_symexpr, var, LHS, "volume", FV, language, gen_framework);
    (rhs_string, rhs_code) = generate_code_layer(rhs_symexpr, var, RHS, "volume", FV, language, gen_framework);
    log_entry("source, code layer: \n  LHS = "*string(lhs_string)*" \n  RHS = "*string(rhs_string));

    if language == JULIA || language == 0
        args = "args; kwargs...";
        makeFunction(args, string(lhs_code));
        set_lhs(var, lhs_string);

        makeFunction(args, string(rhs_code));
        set_rhs(var, rhs_string);

    else
        set_lhs(var, lhs_code);
        set_rhs(var, rhs_code);
    end
end

# Creates an assembly loop function for a given variable that nests the loops in the given order.
function assemblyLoops(var, indices)
    # find the associated indexer objects
    indexer_list = [];
    for i=1:length(indices)
        if typeof(indices[i]) == String
            if indices[i] == "elements" || indices[i] == "cells"
                push!(indexer_list, "elements");
            else
                for j=1:length(indexers)
                    if indices[i] == string(indexers[j].symbol)
                        push!(indexer_list, indexers[j]);
                    end
                end
            end

        elseif typeof(indices[i]) == Indexer
            push!(indexer_list, indices[i]);
        end

    end

    (loop_string, loop_code) = generate_assembly_loops(indexer_list, config.solver_type, language);
    log_entry("assembly loop function: \n\t"*string(loop_string));

    if language == JULIA || language == 0
        if config.solver_type == FV
            args = "var, source_lhs, source_rhs, flux_lhs, flux_rhs, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0";
        elseif config.solver_type == CG
            args = "var, bilinear, linear, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; rhs_only=false";
        end
        makeFunction(args, string(loop_code));
        set_assembly_loops(var);
    else
        set_assembly_loops(var, loop_code);
    end
end

function exportCode(filename)
    # For now, only do this for Julia code because others are already output in code files.
    if language == JULIA || language == 0
        file = open(filename*".jl", "w");
        println(file, "#=\nGenerated functions for "*project_name*"\n=#\n");
        for LorR in [LHS, RHS]
            if LorR == LHS
                codevol = code_strings[1];
                codesurf = code_strings[3];
            else
                codevol = code_strings[2];
                codesurf = code_strings[4];
            end
            for i=1:length(variables)
                var = string(variables[i].symbol);
                if !(codevol[i] === "")
                    func_name = LorR*"_volume_function_for_"*var;
                    println(file, "function "*func_name*"(args; kwargs...)");
                    println(file, codevol[i]);
                    println(file, "end #"*func_name*"\n");
                else
                    println(file, "# No "*LorR*" volume set for "*var*"\n");
                end

                if !(codesurf[i] === "")
                    func_name = LorR*"_surface_function_for_"*var;
                    println(file, "function "*func_name*"(args; kwargs...)");
                    println(file, codesurf[i]);
                    println(file, "end #"*func_name*"\n");
                else
                    println(file, "# No "*LorR*" surface set for "*var*"\n");
                end

                if LorR == RHS
                    if !(assembly_loops[i] == "" || assembly_loops[i] === nothing)
                        func_name = "assembly_function_for_"*var;
                        println(file, "function "*func_name*"(args; kwargs...)");
                        println(file, assembly_loops[i]);
                        println(file, "end #"*func_name*"\n");
                    else
                        println(file, "# No assembly function set for "*var*"\n");
                    end
                end
            end
        end

        close(file);
        log_entry("Exported code to "*filename*".jl");
    else
        # Should we export for other targets?
    end
end

function importCode(filename)
    # For now, only do this for Julia code because others are already output in code files.
    if language == JULIA || language == 0
        file = open(filename*".jl", "r");
        lines = readlines(file, keep=true);
        for LorR in [LHS, RHS]
            if LorR == LHS
                codevol = bilinears;
                codesurf = face_bilinears;
            else
                codevol = linears;
                codesurf = face_linears;
            end

            # Loop over variables and check to see if a matching function is present.
            for i=1:length(variables)
                # Scan the file for a pattern like
                #   function LHS_volume_function_for_u
                #       ...
                #   end #LHS_volume_function_for_u
                #
                # Set LHS/RHS, volume/surface, and the variable name
                var = string(variables[i].symbol);
                vfunc_name = LorR*"_volume_function_for_"*var;
                vfunc_string = "";
                sfunc_name = LorR*"_surface_function_for_"*var;
                sfunc_string = "";
                afunc_name = "assembly_function_for_"*var;
                afunc_string = "";
                for st=1:length(lines)
                    if occursin("function "*vfunc_name, lines[st])
                        # s is the start of the function
                        for en=(st+1):length(lines)
                            if occursin("end #"*vfunc_name, lines[en])
                                # en is the end of the function
                                st = en; # update st
                                break;
                            else
                                vfunc_string *= lines[en];
                            end
                        end
                    elseif occursin("function "*sfunc_name, lines[st])
                        # s is the start of the function
                        for en=(st+1):length(lines)
                            if occursin("end #"*sfunc_name, lines[en])
                                # en is the end of the function
                                st = en; # update st
                                break;
                            else
                                sfunc_string *= lines[en];
                            end
                        end
                    elseif LorR == RHS && occursin("function "*afunc_name, lines[st])
                        # s is the start of the function
                        for en=(st+1):length(lines)
                            if occursin("end #"*afunc_name, lines[en])
                                # en is the end of the function
                                st = en; # update st
                                break;
                            else
                                afunc_string *= lines[en];
                            end
                        end
                    end
                end # lines loop

                # Generate the functions and set them in the right places
                if vfunc_string == ""
                    log_entry("Warning: While importing, no "*LorR*" volume function was found for "*var);
                else
                    makeFunction("args; kwargs...", CodeGenerator.code_string_to_expr(vfunc_string));
                    if LorR == LHS
                        set_lhs(variables[i]);
                    else
                        set_rhs(variables[i]);
                    end
                end
                if sfunc_string == ""
                    log_entry("Warning: While importing, no "*LorR*" surface function was found for "*var);
                else
                    makeFunction("args; kwargs...", CodeGenerator.code_string_to_expr(sfunc_string));
                    if LorR == LHS
                        set_lhs_surface(variables[i]);
                    else
                        set_rhs_surface(variables[i]);
                    end
                end
                if afunc_string == ""
                    log_entry("Warning: While importing, no assembly function was found for "*var*" (using default)");
                else
                    makeFunction("args; kwargs...", CodeGenerator.code_string_to_expr(afunc_string));
                    set_assembly_loops(variables[i]);
                end

            end # vars loop
        end

        close(file);
        log_entry("Imported code from "*filename*".jl");

    else
        # TODO non-julia
    end
end

# Prints a Latex string for the equation for a variable
function printLatex(var)
    if typeof(var) <: Array
        varname = "["*string(var[1].symbol);
        for i=2:length(var)
            varname *= ", "*string(var[i].symbol);
        end
        varname *= "]";
        var = var[1];
    else
        varname = string(var.symbol);
    end
    println("Latex equation for "*varname*":");

    if config.solver_type == FV
        result = raw"$" * " \\frac{d}{dt}\\bar{"*string(var.symbol)*"}";
        if !(symexpressions[1][var.index] === nothing)
            result *= " + \\int_{K}"*symexpression_to_latex(symexpressions[1][var.index])*" dx";
        end
        if !(symexpressions[2][var.index] === nothing)
            result *= " + \\int_{\\partial K}"*symexpression_to_latex(symexpressions[2][var.index])*" ds";
        end
        result *= " = ";
        if !(symexpressions[3][var.index] === nothing)
            result *= "\\int_{K}"*symexpression_to_latex(symexpressions[3][var.index])*" dx";
        end
        if !(symexpressions[4][var.index] === nothing)
            if !(symexpressions[3][var.index] === nothing)
                result *= " + ";
            end
            result *= "\\int_{\\partial K}"*symexpression_to_latex(symexpressions[4][var.index])*" ds";
        end
        if symexpressions[3][var.index] === nothing && symexpressions[4][var.index] === nothing
            result *= "0";
        end
        result *= raw"$";
    else
        if prob.time_dependent
            result = raw"$ \left[";
            if !(symexpressions[1][var.index] === nothing)
                result *= "\\int_{K}"*symexpression_to_latex(symexpressions[1][var.index])*" dx";
            end
            if !(symexpressions[2][var.index] === nothing)
                if !(symexpressions[1][var.index] === nothing)
                    result *= " + ";
                end
                result *= "\\int_{\\partial K}"*symexpression_to_latex(symexpressions[2][var.index])*" ds";
            end
            if symexpressions[1][var.index] === nothing && symexpressions[2][var.index] === nothing
                result *= "0";
            end
            result *= "\\right]_{t+dt} = \\left[";
            if !(symexpressions[3][var.index] === nothing)
                result *= "\\int_{K}"*symexpression_to_latex(symexpressions[3][var.index])*" dx";
            end
            if !(symexpressions[4][var.index] === nothing)
                if !(symexpressions[3][var.index] === nothing)
                    result *= " + ";
                end
                result *= "\\int_{\\partial K}"*symexpression_to_latex(symexpressions[4][var.index])*" ds";
            end
            if symexpressions[3][var.index] === nothing && symexpressions[4][var.index] === nothing
                result *= "0";
            end
            result *= raw"\right]_{t} $";

        else
            result = raw"$";
            if !(symexpressions[1][var.index] === nothing)
                result *= "\\int_{K}"*symexpression_to_latex(symexpressions[1][var.index])*" dx";
            end
            if !(symexpressions[2][var.index] === nothing)
                if !(symexpressions[1][var.index] === nothing)
                    result *= " + ";
                end
                result *= "\\int_{\\partial K}"*symexpression_to_latex(symexpressions[2][var.index])*" ds";
            end
            if symexpressions[1][var.index] === nothing && symexpressions[2][var.index] === nothing
                result *= "0";
            end
            result *= " = ";
            if !(symexpressions[3][var.index] === nothing)
                result *= "\\int_{K}"*symexpression_to_latex(symexpressions[3][var.index])*" dx";
            end
            if !(symexpressions[4][var.index] === nothing)
                if !(symexpressions[3][var.index] === nothing)
                    result *= " + ";
                end
                result *= "\\int_{\\partial K}"*symexpression_to_latex(symexpressions[4][var.index])*" ds";
            end
            if symexpressions[3][var.index] === nothing && symexpressions[4][var.index] === nothing
                result *= "0";
            end
            result *= raw"$";
        end
    end
    println(result);
    println("");

    return result;
end

# Evaluate all of the initial conditions
function eval_initial_conditions()
    dim = config.dimension;

    # build initial conditions
    for vind=1:length(variables)
        if !(variables[vind].ready) && vind <= length(prob.initial)
            if !(prob.initial[vind] === nothing)
                if variables[vind].location == CELL && !(fv_grid === nothing)
                    # Need to use the fv_grid instead of grid_data
                    this_grid_data = fv_grid;
                    this_geo_factors = fv_geo_factors;
                    this_refel = fv_refel;
                else
                    this_grid_data = grid_data;
                    this_geo_factors = geo_factors;
                    this_refel = refel;
                end

                if typeof(prob.initial[vind]) <: Array
                    # Evaluate at nodes
                    nodal_values = zeros(length(prob.initial[vind]), size(this_grid_data.allnodes,2));
                    for ci=1:length(prob.initial[vind])
                        for ni=1:size(this_grid_data.allnodes,2)
                            if typeof(prob.initial[vind][ci]) <: Number
                                nodal_values[ci,ni] = prob.initial[vind][ci];
                            elseif dim == 1
                                nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[ni],0,0,0);
                            elseif dim == 2
                                nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],0,0);
                            elseif dim == 3
                                nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],this_grid_data.allnodes[3,ni],0);
                            end
                        end
                    end

                    # compute cell averages using nodal values if needed
                    if variables[vind].location == CELL
                        nel = size(this_grid_data.loc2glb, 2);
                        for ei=1:nel
                            e = elemental_order[ei];
                            glb = this_grid_data.loc2glb[:,e];
                            vol = this_geo_factors.volume[e];
                            detj = this_geo_factors.detJ[e];

                            for ci=1:length(prob.initial[vind])
                                variables[vind].values[ci,e] = detj / vol * (this_refel.wg' * this_refel.Q * (nodal_values[ci,glb][:]))[1];
                            end
                        end
                    else
                        variables[vind].values = nodal_values;
                    end

                else # scalar
                    nodal_values = zeros(1, size(this_grid_data.allnodes,2));
                    for ni=1:size(this_grid_data.allnodes,2)
                        if typeof(prob.initial[vind]) <: Number
                            nodal_values[ni] = prob.initial[vind];
                        elseif dim == 1
                            nodal_values[ni] = prob.initial[vind].func(this_grid_data.allnodes[ni],0,0,0);
                        elseif dim == 2
                            nodal_values[ni] = prob.initial[vind].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],0,0);
                        elseif dim == 3
                            nodal_values[ni] = prob.initial[vind].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],this_grid_data.allnodes[3,ni],0);
                        end
                    end

                    # compute cell averages using nodal values if needed
                    if variables[vind].location == CELL
                        nel = size(this_grid_data.loc2glb, 2);
                        for e=1:nel
                            glb = this_grid_data.loc2glb[:,e];
                            vol = this_geo_factors.volume[e];
                            detj = this_geo_factors.detJ[e];

                            variables[vind].values[e] = detj / vol * (this_refel.wg' * this_refel.Q * nodal_values[glb])[1];
                        end

                    else
                        variables[vind].values = nodal_values;
                    end

                end

                variables[vind].ready = true;
                log_entry("Built initial conditions for: "*string(variables[vind].symbol));
            end
        end
    end
end

# This will either solve the problem or generate the code for an external target.
function solve(var, nlvar=nothing; nonlinear=false)
    if use_cachesim
        return cachesimSolve(var);
    end

    global time_stepper; # This should not be necessary. It will go away eventually

    # Generate files or solve directly
    if !(!generate_external && (language == JULIA || language == 0)) # if an external code gen target is ready
        if typeof(var) <: Array
            varind = var[1].index;
        else
            varind = var.index;
        end
        generate_all_files(var, bilinears[varind], face_bilinears[varind], linears[varind], face_linears[varind]);

    else
        if typeof(var) <: Array
            varnames = "["*string(var[1].symbol);
            for vi=2:length(var)
                varnames = varnames*", "*string(var[vi].symbol);
            end
            varnames = varnames*"]";
            varind = var[1].index;
        else
            varnames = string(var.symbol);
            varind = var.index;
        end

        # Evaluate initial conditions if not already done
        eval_initial_conditions();

        # If any of the variables are indexed types, generate assembly loops
        var_indexers = [];
        if typeof(var) <: Array
            for vi=1:length(var)
                if !(var[vi].indexer === nothing)
                    push!(var_indexers, var[vi].indexer);
                end
            end
        else
            if !(var.indexer === nothing)
                var_indexers = [var.indexer];
            end
        end
        if length(var_indexers)>0 && assembly_loops[varind] === nothing
            log_entry("Indexed variables detected, but no assembly loops specified. Using default.")
            assemblyLoops(var, ["elements"; var_indexers]);
        end

        # Use the appropriate solver
        if config.solver_type == CG
            lhs = bilinears[varind];
            rhs = linears[varind];
            loop_func = assembly_loops[varind];

            if prob.time_dependent
                time_stepper = init_stepper(grid_data.allnodes, time_stepper);
                if use_specified_steps
                    time_stepper.dt = specified_dt;
				    time_stepper.Nsteps = specified_Nsteps;
                end
                if (nonlinear)
                    if time_stepper.type == EULER_EXPLICIT || time_stepper.type == LSRK4
                        printerr("Warning: Use implicit stepper for nonlinear problem. cancelling solve.");
                        return;
                    end
                	t = @elapsed(result = CGSolver.nonlinear_solve(var, nlvar, lhs, rhs, time_stepper, assemble_func=loop_func));
				else
                	t = @elapsed(result = CGSolver.linear_solve(var, lhs, rhs, time_stepper, assemble_func=loop_func));
				end
                # result is already stored in variables
            else
                # solve it!
				if (nonlinear)
                	t = @elapsed(result = CGSolver.nonlinear_solve(var, nlvar, lhs, rhs, assemble_func=loop_func));
                else
                    t = @elapsed(result = CGSolver.linear_solve(var, lhs, rhs, assemble_func=loop_func));
				end

                # place the values in the variable value arrays
                if typeof(var) <: Array && length(result) > 1
                    tmp = 0;
                    totalcomponents = 0;
                    for vi=1:length(var)
                        totalcomponents = totalcomponents + length(var[vi].symvar);
                    end
                    for vi=1:length(var)
                        components = length(var[vi].symvar);
                        for compi=1:components
                            #println("putting result "*string(compi+tmp)*":"*string(totalcomponents)*":end in var["*string(vi)*"].values["*string(compi)*"]")
                            var[vi].values[compi,:] = result[(compi+tmp):totalcomponents:end];
                        end
                        tmp = tmp + components;
                    end
                elseif length(result) > 1
                    components = length(var.symvar);
                    for compi=1:components
                        var.values[compi,:] = result[compi:components:end];
                    end
                end
            end

            log_entry("Solved for "*varnames*".(took "*string(t)*" seconds)", 1);

        elseif config.solver_type == DG
            lhs = bilinears[varind];
            rhs = linears[varind];
            slhs = face_bilinears[varind];
            srhs = face_linears[varind];

            if prob.time_dependent
                time_stepper = init_stepper(grid_data.allnodes, time_stepper);
                if use_specified_steps
                    time_stepper.dt = specified_dt;
				    time_stepper.Nsteps = specified_Nsteps;
                end
				if (nonlinear)
                	# t = @elapsed(result = DGSolver.nonlinear_solve(var, nlvar, lhs, rhs, slhs, srhs, time_stepper));
                    printerr("Nonlinear solver not ready for DG");
                    return;
				else
                	t = @elapsed(result = DGSolver.linear_solve(var, lhs, rhs, slhs, srhs, time_stepper));
				end
                # result is already stored in variables
            else
                # solve it!
				if (nonlinear)
                	t = @elapsed(result = DGSolver.nonlinear_solve(var, nlvar, lhs, rhs, slhs, srhs));
                else
                    t = @elapsed(result = DGSolver.linear_solve(var, lhs, rhs, slhs, srhs));
				end

                # place the values in the variable value arrays
                if typeof(var) <: Array && length(result) > 1
                    tmp = 0;
                    totalcomponents = 0;
                    for vi=1:length(var)
                        totalcomponents = totalcomponents + length(var[vi].symvar);
                    end
                    for vi=1:length(var)
                        components = length(var[vi].symvar);
                        for compi=1:components
                            var[vi].values[compi,:] = result[(compi+tmp):totalcomponents:end];
                            tmp = tmp + 1;
                        end
                    end
                elseif length(result) > 1
                    components = length(var.symvar);
                    for compi=1:components
                        var.values[compi,:] = result[compi:components:end];
                    end
                end
            end

            log_entry("Solved for "*varnames*".(took "*string(t)*" seconds)", 1);

        elseif config.solver_type == FV
            slhs = bilinears[varind];
            srhs = linears[varind];
            flhs = face_bilinears[varind];
            frhs = face_linears[varind];
            loop_func = assembly_loops[varind];

            if prob.time_dependent
                if fv_grid === nothing
                    time_stepper = init_stepper(grid_data.allnodes, time_stepper);
                else
                    time_stepper = init_stepper(fv_grid.allnodes, time_stepper);
                end
                if use_specified_steps
                    time_stepper.dt = specified_dt;
				    time_stepper.Nsteps = specified_Nsteps;
                end
				if (nonlinear)
                	printerr("Nonlinear solver not ready for FV");
                    return;
				else
                	t = @elapsed(result = FVSolver.linear_solve(var, slhs, srhs, flhs, frhs, time_stepper, loop_func));
				end
                # result is already stored in variables
                log_entry("Solved for "*varnames*".(took "*string(t)*" seconds)", 1);
            else
                # does this make sense?
                printerr("FV assumes time dependence. Set initial conditions etc.");
            end
        end
    end

    # At this point all of the final values for the variables in var should be placed in var.values.
end

# When using cachesim, this will be used to simulate the solve.
function cachesimSolve(var, nlvar=nothing; nonlinear=false)
    if !(!generate_external && (language == JULIA || language == 0))
        printerr("Cachesim solve is only ready for Julia direct solve");
    else
        if typeof(var) <: Array
            varnames = "["*string(var[1].symbol);
            for vi=2:length(var)
                varnames = varnames*", "*string(var[vi].symbol);
            end
            varnames = varnames*"]";
            varind = var[1].index;
        else
            varnames = string(var.symbol);
            varind = var.index;
        end

        lhs = bilinears[varind];
        rhs = linears[varind];

        t = @elapsed(result = linear_solve_cachesim(var, lhs, rhs));
        log_entry("Generated cachesim ouput for "*varnames*".(took "*string(t)*" seconds)", 1);
    end
end

# Writes the values for each variable to filename in the given format.
function output_values(vars, filename; format="raw", ascii=false)
    available_formats = ["raw", "csv", "vtk"];
    if format == "vtk"
        output_values_vtk(vars, filename, ascii)
        log_entry("Writing values to file: "*filename*".vtu");

    else
        filename *= "."*format;
        file = open(filename, "w");
        log_entry("Writing values to file: "*filename);

        if format == "raw"
            output_values_raw(vars, file)
        elseif format == "csv"
            output_values_csv(vars, file)
        else
            println("Unknown output file format("*string(format)*"). Choose from: "*string(available_formats));
        end
        log_entry("Finished writing to "*filename);

        close(file);
    end
end

function finalize_femshop()
    # Finalize generation
    finalize_code_generator();
    if use_cachesim
        CachesimOut.finalize();
    end
    # anything else
    close_log();
    println("Femshop has completed.");
end

### Other specialized functions ###

function cachesim(use)
    log_entry("Using cachesim - Only cachesim output will be generated.", 1);
    global use_cachesim = use;
end

function morton_nodes(griddim)
    t = @elapsed(global grid_data = reorder_grid_recursive(grid_data, griddim, MORTON_ORDERING));
    log_entry("Reordered nodes to Morton. Took "*string(t)*" sec.", 2);
end

function morton_elements(griddim)
    global elemental_order = get_recursive_order(MORTON_ORDERING, config.dimension, griddim);
    log_entry("Reordered elements to Morton.", 2);
end

function hilbert_nodes(griddim)
    t = @elapsed(global grid_data = reorder_grid_recursive(grid_data, griddim, HILBERT_ORDERING));
    log_entry("Reordered nodes to Hilbert. Took "*string(t)*" sec.", 2);
end

function hilbert_elements(griddim)
    global elemental_order = get_recursive_order(HILBERT_ORDERING, config.dimension, griddim);
    log_entry("Reordered elements to Hilbert.", 2);
end

function tiled_nodes(griddim, tiledim)
    t = @elapsed(global grid_data = reorder_grid_tiled(grid_data, griddim, tiledim));
    log_entry("Reordered nodes to tiled. Took "*string(t)*" sec.", 2);
end

function tiled_elements(griddim, tiledim)
    global elemental_order = get_tiled_order(config.dimension, griddim, tiledim, true);
    log_entry("Reordered elements to tiled("*string(tiledim)*").", 2);
end

function ef_nodes()
    t = @elapsed(global grid_data = reorder_grid_element_first(grid_data, config.basis_order_min, elemental_order));
    log_entry("Reordered nodes to EF. Took "*string(t)*" sec.", 2);
end

function random_nodes(seed = 17)
    t = @elapsed(global grid_data = reorder_grid_random(grid_data, seed));
    log_entry("Reordered nodes to random. Took "*string(t)*" sec.", 2);
end

function random_elements(seed = 17)
    global elemental_order = random_order(mesh_data.nel, seed);
    log_entry("Reordered elements to random.", 2);
end


text:
@testset "filter" begin
    r, c = (100, 500)
    data = rand(0:10, r, c)
    obs = DataFrame(A=rand(c), B=rand(c))
    var = DataFrame(C=rand(r), D=rand(r))
    layer_a = rand(r, c)
    layer_b = rand(c, r)

    prof = Profile(data, var, obs)
    prof.layers[:a] = copy(layer_a)
    prof.layers[:b] = copy(layer_b)

    min_genes = 90
    expressed_genes = vec(sum(data .!= 0, dims=1))
    filter_cells!(prof; min_genes=min_genes)
    @test prof.data == data[:, min_genes .<= expressed_genes]
    @test prof.obs == obs[min_genes .<= expressed_genes, :]
    @test prof.var == var
    @test prof.layers[:a] == layer_a[:, min_genes .<= expressed_genes]
    @test prof.layers[:b] == layer_b
    @test_throws ArgumentError filter_cells!(prof; min_genes=-1)
    @test_throws ArgumentError filter_cells!(prof; max_genes=-1)
    @test_throws ArgumentError filter_cells!(prof; min_counts=-1)
    @test_throws ArgumentError filter_cells!(prof; max_counts=-1)

    prof = Profile(data, var, obs)
    prof.layers[:a] = copy(layer_a)
    prof.layers[:b] = copy(layer_b)

    min_cells = 300
    expressed_cells = vec(sum(data .!= 0, dims=2))
    filter_genes!(prof, min_cells=min_cells)
    @test prof.data == data[min_cells .<= expressed_cells, :]
    @test prof.obs == obs
    @test prof.var == var[min_cells .<= expressed_cells, :]
    @test prof.layers[:a] == layer_a[min_cells .<= expressed_cells, :]
    @test prof.layers[:b] == layer_b
    @test_throws ArgumentError filter_genes!(prof; min_cells=-1)
    @test_throws ArgumentError filter_genes!(prof; max_cells=-1)
    @test_throws ArgumentError filter_genes!(prof; min_counts=-1)
    @test_throws ArgumentError filter_genes!(prof; max_counts=-1)
end

ext: jl


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

using Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2

@testset for eltya in (Float32, Float64, Complex64, Complex128, Int)
    aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    asym = aa'+aa                  # symmetric indefinite
    apd  = aa'*aa                 # symmetric positive-definite
    for (a, asym, apd) in ((aa, asym, apd),
                           (view(aa, 1:n, 1:n),
                            view(asym, 1:n, 1:n),
                            view(apd, 1:n, 1:n)))
        ε = εa = eps(abs(float(one(eltya))))

        α = rand(eltya)
        β = rand(eltya)
        eab = eig(α,β)
        @test eab[1] == eigvals(fill(α,1,1),fill(β,1,1))
        @test eab[2] == eigvecs(fill(α,1,1),fill(β,1,1))

        @testset "non-symmetric eigen decomposition" begin
            d, v = eig(a)
            for i in 1:size(a,2)
                @test a*v[:,i] ≈ d[i]*v[:,i]
            end
            f = eigfact(a)
            @test det(a) ≈ det(f)
            @test inv(a) ≈ inv(f)
            @test isposdef(a) == isposdef(f)
            @test eigvals(f) === f[:values]
            @test eigvecs(f) === f[:vectors]
            @test Array(f) ≈ a

            num_fact = eigfact(one(eltya))
            @test num_fact.values[1] == one(eltya)
            h = asym
            @test minimum(eigvals(h)) ≈ eigmin(h)
            @test maximum(eigvals(h)) ≈ eigmax(h)
            @test_throws DomainError eigmin(a - a')
            @test_throws DomainError eigmax(a - a')
        end
        @testset "symmetric generalized eigenproblem" begin
            if isa(a, Array)
                asym_sg = asym[1:n1, 1:n1]
                a_sg = a[:,n1+1:n2]
            else
                asym_sg = view(asym, 1:n1, 1:n1)
                a_sg = view(a, 1:n, n1+1:n2)
            end
            f = eigfact(asym_sg, a_sg'a_sg)
            @test asym_sg*f[:vectors] ≈ (a_sg'a_sg*f[:vectors]) * Diagonal(f[:values])
            @test f[:values] ≈ eigvals(asym_sg, a_sg'a_sg)
            @test prod(f[:values]) ≈ prod(eigvals(asym_sg/(a_sg'a_sg))) atol=200ε
            @test eigvecs(asym_sg, a_sg'a_sg) == f[:vectors]
            @test eigvals(f) === f[:values]
            @test eigvecs(f) === f[:vectors]
            @test_throws KeyError f[:Z]

            d,v = eig(asym_sg, a_sg'a_sg)
            @test d == f[:values]
            @test v == f[:vectors]
        end
        @testset "Non-symmetric generalized eigenproblem" begin
            if isa(a, Array)
                a1_nsg = a[1:n1, 1:n1]
                a2_nsg = a[n1+1:n2, n1+1:n2]
            else
                a1_nsg = view(a, 1:n1, 1:n1)
                a2_nsg = view(a, n1+1:n2, n1+1:n2)
            end
            f = eigfact(a1_nsg, a2_nsg)
            @test a1_nsg*f[:vectors] ≈ (a2_nsg*f[:vectors]) * Diagonal(f[:values])
            @test f[:values] ≈ eigvals(a1_nsg, a2_nsg)
            @test prod(f[:values]) ≈ prod(eigvals(a1_nsg/a2_nsg)) atol=50000ε
            @test eigvecs(a1_nsg, a2_nsg) == f[:vectors]
            @test_throws KeyError f[:Z]

            d,v = eig(a1_nsg, a2_nsg)
            @test d == f[:values]
            @test v == f[:vectors]
        end
    end
end

@testset "eigenvalue computations with NaNs" begin
    for eltya in (NaN16, NaN32, NaN)
        @test_throws(ArgumentError, eig(fill(eltya, 1, 1)))
        @test_throws(ArgumentError, eig(fill(eltya, 2, 2)))
        test_matrix = rand(typeof(eltya),3,3)
        test_matrix[2,2] = eltya
        @test_throws(ArgumentError, eig(test_matrix))
    end
end

# test a matrix larger than 140-by-140 for #14174
let aa = rand(200, 200)
    for a in (aa, view(aa, 1:n, 1:n))
        f = eigfact(a)
        @test a ≈ f[:vectors] * Diagonal(f[:values]) / f[:vectors]
    end
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
const celldm_1 = Symbol("celldm(1)")
const QE_EXECS = [
    "pw.x",
    "projwfc.x",
    "pp.x",
    "ld1.x",
    "ph.x",
    "pw2wannier90.x"
]
#REVIEW: Should we make the flag name a String?
#QE calls these flags
struct QEFlagInfo{T}
    name::Symbol
    # default::Union{T,Nothing,Symbol}  #check again v0.7 Some
    description::String
end
QEFlagInfo() = QEFlagInfo{Nothing}(:error, "")
Base.eltype(x::QEFlagInfo{T}) where T = T

abstract type AbstractBlockInfo end

struct QEControlBlockInfo <: AbstractBlockInfo
    name::Symbol
    flags::Vector{<:QEFlagInfo}
end

#TODO rewrite this this is garbage
function qe_flaginfo(block::AbstractBlockInfo, variable_name::Symbol)
    varstr1 = string(variable_name)
    for var in block.flags
        varstr2 = string(var.name)
        if occursin(varstr2, varstr1) && (length(varstr1) == length(varstr2) || length(varstr1) == length(varstr2) + 2)
            return var
        end
    end
    return QEFlagInfo()
end

struct QEDataBlockInfo <: AbstractBlockInfo
    name                ::Symbol
    description         ::String
    options             ::Vector{Symbol}
    options_description ::String
    flags               ::Vector{<:QEFlagInfo}
end

struct QEInputInfo
    exec::String
    control::Vector{QEControlBlockInfo}
    data::Vector{QEDataBlockInfo}
end

allflags(info::QEInputInfo) = flatten([[i.flags for i in info.control]; [i.flags for i in info.data]])

include(joinpath(depsdir, "qeflags.jl"))
const QEInputInfos = _QEINPUTINFOS()
push!(QEInputInfos, QEInputInfo("pw2wannier90.x", [QEControlBlockInfo(:inputpp,[QEFlagInfo{String}(:outdir        , "location of temporary output files"),
                          QEFlagInfo{String}(:prefix        , "pwscf filename prefix"),
                          QEFlagInfo{String}(:seedname      , "wannier90 input/output filename prefix"),
                          QEFlagInfo{String}(:wan_mode      , "'standalone' or 'library'"),
                          QEFlagInfo{String}(:spin_component, "'none', 'up' or 'down'"),
                          QEFlagInfo{Bool}(:write_mmn     , "compute M_mn matrix"),
                          QEFlagInfo{Bool}(:write_amn     , "compute A_mn matrix"),
                          QEFlagInfo{Bool}(:write_unk     , "write wavefunctions to file"),
                          QEFlagInfo{Bool}(:write_uHu     , "write the hamiltonian elements between different k-values"),
                          QEFlagInfo{Bool}(:wvfn_formatted, "formatted or unformatted output for wavefunctions"),
                          QEFlagInfo{Bool}(:reduce_unk    , "output wavefunctions on a coarse grid to save memory")])], QEDataBlockInfo[]))

qe_input_info(input::DFInput{QE}) = getfirst(x-> occursin(x.exec, input.exec), QEInputInfos)
qe_input_info(exec::AbstractString) = getfirst(x-> occursin(x.exec, exec), QEInputInfos)
qe_input_flags(exec::AbstractString) = allflags(qe_input_info(exec))

function qe_flaginfo(input_info::QEInputInfo, variable_name::Symbol)
    for block in vcat(input_info.control, input_info.data)
        var = qe_flaginfo(block, variable_name)
        if eltype(var) != Nothing
            return var
        end
    end
    return QEFlagInfo()
end

function qe_flaginfo(variable_name::Symbol)
    for info in QEInputInfos
        var = qe_flaginfo(info, variable_name)
        if eltype(var) != Nothing
            return var
        end
    end
    return QEFlagInfo()
end

function qe_block_variable(input_info::QEInputInfo, variable_name)
    for block in vcat(input_info.control, input_info.data)
        var = qe_flaginfo(block, variable_name)
        if eltype(var) != Nothing
            return block, var
        end
    end
    return :error, QEFlagInfo()
end

function qe_flaginfo(exec::Exec, varname)
    for input_info in QEInputInfos
        if occursin(input_info.exec, exec.exec)
            return qe_flaginfo(input_info, varname)
        end
    end
    return QEFlagInfo()
end

function qe_block_info(block_name::Symbol)
    for input_info in QEInputInfos
        for block in [input_info.control;input_info.data]
            if block.name == block_name
                return block
            end
        end
    end
end

qe_all_block_flags(input::DFInput{QE}, block_name) = getfirst(x -> x.name == block, qe_input_info(input).control).flags
qe_all_block_flags(exec::AbstractString, block_name) = getfirst(x -> x.name == block_name, qe_input_info(exec).control).flags

function qe_block_variable(exec::AbstractString, flagname)
    for input_info in QEInputInfos
        if occursin(input_info.exec, exec)
            return qe_block_variable(input_info, flagname)
        end
    end
    return :error, QEFlagInfo()
end

function qe_exec(input::DFInput{QE})
    exec = getfirst(x -> x.exec ∈ QE_EXECS, execs(input))
    if exec === nothing
        error("Input $input does not have a valid QE executable, please set it first.")
    end
    return exec
end

qe_block_variable(input::DFInput, flagname) = qe_block_variable(qe_exec(input).exec, flagname)

flagtype(input::DFInput{QE}, flag) = eltype(qe_flaginfo(qe_exec(input), flag))
flagtype(::Type{QE}, exec, flag) = eltype(qe_flaginfo(exec, flag))


text:
# This is the source code for the "Loading Data with Julia" vignette:
#
using Pkg

Pkg.activate("covid19_env")
Pkg.add(PackageSpec(name = "PlotlyJS", version = "0.18.8"))
Pkg.add(PackageSpec(name = "CSV", version = "0.9.6"))
Pkg.add(PackageSpec(name = "Chain", version = "0.4.8"))
Pkg.add(PackageSpec(name = "DataFrames", version = "1.2.2"))

using  CSV, DataFrames, Chain, PlotlyJS

Pkg.status()

# Loading the data
url = "https://raw.githubusercontent.com/RamiKrispin/coronavirus/master/csv/coronavirus.csv"
file = CSV.File(download(url), missingstring= "NA")

df = DataFrame(file)
df

describe(df)

@chain df begin
    filter(:type => ==("confirmed"),_)
    groupby([:combined_key, :type])
    combine([:cases] .=> sum)
    sort!([:cases_sum], rev = true)
end

@chain df begin
    filter(:type => ==("death"),_)
    groupby([:combined_key, :type])
    combine([:cases] .=> sum)
    sort!([:cases_sum], rev = true)
end

# Visualize the data
df_brazil_confirmed = filter(row -> row.country == "Brazil" && row.type == "confirmed", df)
df_brazil_death = filter(row -> row.country == "Brazil" && row.type == "death", df)

p = make_subplots(rows=2,
                  cols=1,
                  shared_xaxes=true,
                  x_title = "Source: Johns Hopkins University Center for Systems Science and Engineering",
                  vertical_spacing=0.02)

add_trace!(p, scatter(df_brazil_confirmed,x=:date, y=:cases, name = "Confirmed"), row=1, col=1)
add_trace!(p, scatter(df_brazil_death,x=:date, y=:cases, name = "Death"), row=2, col=1)

relayout!(p, title_text="Brazil - Daily Confirmed and Death Cases")
p


text:
include("../src/Curvelet.jl")
#using Curvelet
using PyPlot

PyPlot.close("all")

# Load sample data

input = imread("test/peppers.gif")
#input = imread("test/lena.jpg")
#input = input[:,:,1]
sample_data = input*1.0
N = size(sample_data,1)

# Display unmodified data
figure(1)
gray()
imshow(sample_data)
title("Unmodified Data")

println("Performing Conversion...")
plan        = Curvelet.planCurveletTransform(N, Curvelet.C1)
tmp         = Curvelet.curveletTransform(sample_data, plan)
sample_hat  = Curvelet.inverseCurveletTransform(tmp, plan)
println("Done")
println()
println("(Small) Error from transform was: ", norm(sample_data-sample_hat))
println()

# Display a subset of coefficients
combo = zeros(size(tmp.HP[1]))
N = length(tmp.HP)
for i=1:N
    combo += abs.(tmp.HP[i])
end
figure(10)
title("High Frequency Curves")
imshow(log.(combo))

figure(11)
title("A Single Curvelet Band")
imshow(log.(abs.(tmp.HP[5])))

#
##Create an approximation
#
function vectorize_curvelet(cc::Matrix{Complex{Float64}})
    vec(cc)
end

function vectorize_curvelet(cc::Curvelet.CurveletCoeff)
    data = vectorize_curvelet(cc.LP)
    for i = cc.HP
        data = [data,vec(i)]
    end
    data
end

function build_sparse(image, plan, sparsity)
    tmp = Curvelet.curveletTransform(sample_data, plan)
    data = vectorize_curvelet(tmp)
    ordered = reverse(sort(abs.(data)))
    curvelet_threshold(tmp, ordered[sparsity])
    Curvelet.inverseCurveletTransform(tmp, plan)
end
data = Float64[0.0]
for i=tmp.HP
    data = vcat(data, vec(i))
end
#
data = vec(data)
ordered = reverse(sort(abs.(data)))

#threshold = ordered[30000]
threshold = ordered[300]

function curvelet_threshold(cc::Matrix{Complex{Float64}}, threshold)
    #cc[abs(cc) .< threshold] = 0
end

function curvelet_threshold(cc::Curvelet.CurveletCoeff, threshold)
    curvelet_threshold(cc.LP, threshold)
    for i = 1:length(cc.HP)
        cc.HP[i][abs.(cc.HP[i]) .< threshold] = 0
    end
end

function reconstruction_step(observation, sparsity)
    tmp = Curvelet.curveletTransform(observation, plan)
    data = Float64[0.0]
    for i=tmp.HP
        data = vcat(data, vec(i))
    end
    #
    data = vec(data)
    ordered = reverse(sort(abs.(data)))
    threshold = ordered[sparsity]

    curvelet_threshold(tmp, threshold)
    sample_hat  = Curvelet.inverseCurveletTransform(tmp, plan)
end

begin# reconstruction_example()
    #attempt a reconstruction problem
    #grab around 30% of the points
    px_rate = 0.30

    N = size(sample_data,1)
    good_px = zeros(size(sample_data))
    for i=1:N, j=1:N
        good_px[i,j] = (rand() < px_rate) ? 1 : 0
    end

    observation = sample_data.*good_px
    sample_hat  = copy(observation)

    figure(100)
    title("Approximation Under Missing Data - Missing")
    gray()
    imshow(abs.(observation))

    for iters=1:40
        tr = 1000

        if(iters > 10)
            tr = 10000
        elseif(iters > 20)
            tr = 40000
        end

        smp_next = reconstruction_step(sample_hat, tr)
        println("Appx error  was: ", norm(sample_data-smp_next), " @ iter ",
        iters)
        sample_hat[good_px.==0] = smp_next[good_px.==0]
        if(iters in [1, 5, 10])
            figure(iters+100)
            title("Intermediate Denoising Stage")
            imshow(sample_hat)
        end
    end

    println("Denoising Missing Data Example")
    println()
    println("Obs  error  was: ", norm(sample_data-observation))
    println("Appx error  was: ", norm(sample_data-sample_hat))
    println("Orig energy was: ", norm(sample_data))

    figure(151)
    title("Approximation Under Missing Data - Reconstruction")
    gray()
    imshow(abs.(sample_hat))
    figure(152)
    title("Approximation Under Missing Data - Diff")
    gray()
    imshow(abs.(sample_hat).-sample_data)
end

tmp = Curvelet.curveletTransform(sample_data, plan)
K = round(Int,length(sample_data)*0.10)
threshold = ordered[K]
curvelet_threshold(tmp,threshold)

sample_hat  = Curvelet.inverseCurveletTransform(tmp, plan)
println()
println()
println("Creating a sparse represntation")
println()
println("Sparsity:       ", 100*K/length(sample_data), "%")
println("total energy:   ", norm(sample_data))
println("Appx error was: ", norm(sample_data-sample_hat))

figure(200)
title("Approximation Under Sparsity - Original")
gray()
imshow(abs.(sample_data))
figure(201)
title(string("Approximation Under Sparsity - ",K," sparse"))
gray()
imshow(abs.(sample_hat))


text:
function BoxAndSolarArraySpacecraft(arg0::Vector{Facet}, arg1::PVCoordinatesProvider, arg10::jdouble, arg2::jdouble, arg3::Vector3D, arg4::AbsoluteDate, arg5::Vector3D, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BoxAndSolarArraySpacecraft((Vector{Facet}, PVCoordinatesProvider, jdouble, jdouble, Vector3D, AbsoluteDate, Vector3D, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function BoxAndSolarArraySpacecraft(arg0::Vector{Facet}, arg1::PVCoordinatesProvider, arg2::jdouble, arg3::Vector3D, arg4::AbsoluteDate, arg5::Vector3D, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BoxAndSolarArraySpacecraft((Vector{Facet}, PVCoordinatesProvider, jdouble, Vector3D, AbsoluteDate, Vector3D, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function BoxAndSolarArraySpacecraft(arg0::Vector{Facet}, arg1::PVCoordinatesProvider, arg2::jdouble, arg3::Vector3D, arg4::jdouble, arg5::jdouble, arg6::jdouble)
    return BoxAndSolarArraySpacecraft((Vector{Facet}, PVCoordinatesProvider, jdouble, Vector3D, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function BoxAndSolarArraySpacecraft(arg0::Vector{Facet}, arg1::PVCoordinatesProvider, arg2::jdouble, arg3::Vector3D, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble)
    return BoxAndSolarArraySpacecraft((Vector{Facet}, PVCoordinatesProvider, jdouble, Vector3D, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function BoxAndSolarArraySpacecraft(arg0::jdouble, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg2::jdouble, arg3::PVCoordinatesProvider, arg4::jdouble, arg5::Vector3D, arg6::AbsoluteDate, arg7::Vector3D, arg8::jdouble, arg9::jdouble)
    return BoxAndSolarArraySpacecraft((jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, PVCoordinatesProvider, jdouble, Vector3D, AbsoluteDate, Vector3D, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function BoxAndSolarArraySpacecraft(arg0::jdouble, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg2::jdouble, arg3::PVCoordinatesProvider, arg4::jdouble, arg5::Vector3D, arg6::AbsoluteDate, arg7::Vector3D, arg8::jdouble, arg9::jdouble)
    return BoxAndSolarArraySpacecraft((jdouble, jdouble, jdouble, jdouble, jdouble, PVCoordinatesProvider, jdouble, Vector3D, AbsoluteDate, Vector3D, jdouble, jdouble), arg0, arg1, arg10, arg11, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function BoxAndSolarArraySpacecraft(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::PVCoordinatesProvider, arg4::jdouble, arg5::Vector3D, arg6::jdouble, arg7::jdouble, arg8::jdouble)
    return BoxAndSolarArraySpacecraft((jdouble, jdouble, jdouble, PVCoordinatesProvider, jdouble, Vector3D, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function BoxAndSolarArraySpacecraft(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::PVCoordinatesProvider, arg4::jdouble, arg5::Vector3D, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BoxAndSolarArraySpacecraft((jdouble, jdouble, jdouble, PVCoordinatesProvider, jdouble, Vector3D, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function drag_acceleration(obj::BoxAndSolarArraySpacecraft, arg0::AbsoluteDate, arg1::Frame, arg2::Vector3D, arg3::Rotation, arg4::jdouble, arg5::jdouble, arg6::Vector3D, arg7::Vector{jdouble})
    return jcall(obj, "dragAcceleration", Vector3D, (AbsoluteDate, Frame, Vector3D, Rotation, jdouble, jdouble, Vector3D, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function drag_acceleration(obj::BoxAndSolarArraySpacecraft, arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldVector3D, arg3::FieldRotation, arg4::RealFieldElement, arg5::RealFieldElement, arg6::FieldVector3D, arg7::Vector{RealFieldElement})
    return jcall(obj, "dragAcceleration", FieldVector3D, (FieldAbsoluteDate, Frame, FieldVector3D, FieldRotation, RealFieldElement, RealFieldElement, FieldVector3D, Vector{RealFieldElement}), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function get_drag_parameters_drivers(obj::BoxAndSolarArraySpacecraft)
    return jcall(obj, "getDragParametersDrivers", Vector{ParameterDriver}, ())
end

function get_normal(obj::BoxAndSolarArraySpacecraft, arg0::AbsoluteDate, arg1::Frame, arg2::Vector3D, arg3::Rotation)
    return jcall(obj, "getNormal", Vector3D, (AbsoluteDate, Frame, Vector3D, Rotation), arg0, arg1, arg2, arg3)
end

function get_normal(obj::BoxAndSolarArraySpacecraft, arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldVector3D, arg3::FieldRotation)
    return jcall(obj, "getNormal", FieldVector3D, (FieldAbsoluteDate, Frame, FieldVector3D, FieldRotation), arg0, arg1, arg2, arg3)
end

function get_radiation_parameters_drivers(obj::BoxAndSolarArraySpacecraft)
    return jcall(obj, "getRadiationParametersDrivers", Vector{ParameterDriver}, ())
end

function radiation_pressure_acceleration(obj::BoxAndSolarArraySpacecraft, arg0::AbsoluteDate, arg1::Frame, arg2::Vector3D, arg3::Rotation, arg4::jdouble, arg5::Vector3D, arg6::Vector{jdouble})
    return jcall(obj, "radiationPressureAcceleration", Vector3D, (AbsoluteDate, Frame, Vector3D, Rotation, jdouble, Vector3D, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function radiation_pressure_acceleration(obj::BoxAndSolarArraySpacecraft, arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldVector3D, arg3::FieldRotation, arg4::RealFieldElement, arg5::FieldVector3D, arg6::Vector{RealFieldElement})
    return jcall(obj, "radiationPressureAcceleration", FieldVector3D, (FieldAbsoluteDate, Frame, FieldVector3D, FieldRotation, RealFieldElement, FieldVector3D, Vector{RealFieldElement}), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end


text:
Base.@kwdef struct NativeCompilerTarget <: GPUCompiler.AbstractCompilerTarget
    cpu::String=(LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUName())
    features::String=(LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUFeatures())
end

GPUCompiler.llvm_triple(::NativeCompilerTarget) = Sys.MACHINE

function GPUCompiler.llvm_machine(target::NativeCompilerTarget)
    triple = GPUCompiler.llvm_triple(target)

    t = LLVM.Target(triple=triple)

    tm = LLVM.TargetMachine(t, triple, target.cpu, target.features, reloc=LLVM.API.LLVMRelocPIC)
    GPUCompiler.asm_verbosity!(tm, true)

    return tm
end

GPUCompiler.runtime_slug(job::GPUCompiler.CompilerJob{NativeCompilerTarget}) = "native_$(job.target.cpu)-$(hash(job.target.features))"

module StaticRuntime
    # the runtime library
    signal_exception() = return
    malloc(sz) = ccall("extern malloc", llvmcall, Csize_t, (Csize_t,), sz)
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end

struct StaticCompilerParams <: GPUCompiler.AbstractCompilerParams end

GPUCompiler.runtime_module(::GPUCompiler.CompilerJob{<:Any,StaticCompilerParams}) = StaticRuntime
GPUCompiler.runtime_module(::GPUCompiler.CompilerJob{NativeCompilerTarget}) = StaticRuntime
GPUCompiler.runtime_module(::GPUCompiler.CompilerJob{NativeCompilerTarget, StaticCompilerParams}) = StaticRuntime

GPUCompiler.can_throw(job::GPUCompiler.CompilerJob{<:Any,StaticCompilerParams}) = true
GPUCompiler.can_throw(job::GPUCompiler.CompilerJob{NativeCompilerTarget, StaticCompilerParams}) = true
GPUCompiler.can_throw(job::GPUCompiler.CompilerJob{NativeCompilerTarget}) = true

function native_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, name=GPUCompiler.safe_name(repr(func)), kwargs...)
    source = GPUCompiler.FunctionSpec(func, Base.to_tuple_type(types), kernel, name)
    target = NativeCompilerTarget()
    params = StaticCompilerParams()
    GPUCompiler.CompilerJob(target, source, params), kwargs
end


text:
# This file provides support for the AbstractTrees interface.

using AbstractTrees

# Traits
Base.eltype(::Type{<:TreeIterator{Node{T}}}) where T = Node{T}
Base.IteratorEltype(::Type{<:TreeIterator{Node{T}}}) where T = Base.HasEltype()
Base.parent(node::Node, state::Node) = state.parent

AbstractTrees.parentlinks(::Type{Node{T}}) where T = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{Node{T}}) where T = AbstractTrees.StoredSiblings()
AbstractTrees.rootstate(node::Node) = node
AbstractTrees.children(node::Node) = node
function AbstractTrees.nextsibling(tree::Node, state::Node)
    ns = state.sibling
    return state === ns ? nothing : ns
end
AbstractTrees.printnode(io::IO, node::Node) = print(io, node.data)


text:
function cswap!(M::Matrix{T}, j1::T,j2::T) where T <: Integer
    j1 == j2 && return M
    M[:,[j1,j2]] = M[:,[j2,j1]]
    return M
end
function rswap!(M::Matrix{T}, i1::T,i2::T) where T <: Integer
    i1 == i2 && return M
    M[[i1,i2],:] = M[[i2,i1],:]
    return M
end
function pullfirst!(M::Matrix{T}, idx) where T <: Integer
    rswap!(M, 1, idx[1])
    cswap!(M, 1, idx[2])
end

function eratosthenes(n::T) where T <: Integer
    if n == 1 return [[1]], [] end
    if n == 2 return [[1], [2]], [2] end
    primes = [2]
    factorized = [[1], [2]]
    for k ∈ 3:n
        m = k
        for p ∈ primes
            if m % p == 0
                m ÷= p
                temp = [p; factorized[m]]
                push!(factorized, temp)
                break
            end
        end
        if length(factorized) != k
            push!(primes, k)
            push!(factorized, [k])
        end
    end
    return factorized, primes
end

function eratosthenes!(factorized::Vector{Vector{T}}, primes::Vector{T}, n::T) where T <: Integer
    old_n = length(factorized)
    if n ≤ old_n return factorized, primes end
    for k ∈ (old_n + 1):n
        m = k
        for p ∈ primes
            if m % p == 0
                m ÷= p
                temp = [p; factorized[m]]
                push!(factorized, temp)
                break
            end
        end
        if length(factorized) != k
            push!(primes, k)
            push!(factorized, [k])
        end
    end
    return factorized, primes
end

F, P = eratosthenes(10)
eratosthenes!(F, P, 100)
static_nfactor = length.(F)

function nfactor(k::T) where T <: Integer
    if k < 0 k = abs(k) end
    if k == 0 return typemax(T) end
    if k == 1 return 0 end
    return static_nfactor[k]
end

M = [
   -6 3 2 2 5
    2 2 0 2 7
   -2 5 2 4 3
   -7 3 7 8 9
]

idx = argmin(nfactor.(M))
pullfirst!(M, idx)

nfactor.(M[begin,:])
nfactor.(M[:,begin])


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    SatelliteToolbox orbit propagator API: step
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export step!

"""
    step!(orbp, Δt::Number)

Propagate the orbit in `orbp` by `Δt` [s] using the algorithm of `orbp`. The new
parameters will be written in `orbp`.

# Returns

* The Keplerian elements represented in the inertial frame after the step (see
  `Orbit`) [SI units].
* The position vector represented in the inertial frame after the step [m].
* The velocity vector represented in the inertial frame after the step [m].

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. If the orbit parameters are obtained
from a TLE, then the inertial frame will be TEME. Notice, however, that the
perturbation theory requires an inertial frame with true equator.

"""
function step!(orbp::OrbitPropagatorJ2, Δt::Number)
    # Auxiliary variables.
    orb = orbp.orb
    j2d = orbp.j2d

    # Propagate the orbit.
    (r_i, v_i) = j2!(j2d, j2d.Δt + Δt)

    # Update the elements in the `orb` structure.
    orb.t += Δt/86400
    orb.a  = j2d.al_k*j2d.j2_gc.R0
    orb.e  = j2d.e_k
    orb.i  = j2d.i_k
    orb.Ω  = j2d.Ω_k
    orb.ω  = j2d.ω_k
    orb.f  = j2d.f_k

    # Return the information about the step.
    (copy(orbp.orb), r_i, v_i)
end

function step!(orbp::OrbitPropagatorJ4, Δt::Number)
    # Auxiliary variables.
    orb = orbp.orb
    j4d = orbp.j4d

    # Propagate the orbit.
    (r_i, v_i) = j4!(j4d, j4d.Δt + Δt)

    # Update the elements in the `orb` structure.
    orb.t += Δt/86400
    orb.a  = j4d.al_k*j4d.j4_gc.R0
    orb.e  = j4d.e_k
    orb.i  = j4d.i_k
    orb.Ω  = j4d.Ω_k
    orb.ω  = j4d.ω_k
    orb.f  = j4d.f_k

    # Return the information about the step.
    (copy(orbp.orb), r_i, v_i)
end

function step!(orbp::OrbitPropagatorSGP4{T}, Δt::Number) where T
    # Auxiliary variables.
    orb     = orbp.orb
    sgp4d   = orbp.sgp4d
    sgp4_gc = orbp.sgp4_gc

    # Propagate the orbit.
    (r_teme, v_teme) = sgp4!(sgp4d, sgp4d.Δt + Δt/60)

    # Convert km to m.
    r_teme *= 1000
    v_teme *= 1000

    # Update the elements in the `orb` structure.
    orb.t = sgp4d.epoch + sgp4d.Δt*60/86400
    orb.a = sgp4d.a_k*sgp4_gc.R0*1000
    orb.e = sgp4d.e_k
    orb.i = sgp4d.i_k
    orb.Ω = sgp4d.Ω_k
    orb.ω = sgp4d.ω_k
    orb.f = M_to_f(sgp4d.e_k, sgp4d.M_k)

    # Return the information about the step.
    (copy(orbp.orb), r_teme, v_teme)
end

function step!(orbp::OrbitPropagatorTwoBody, Δt::Number)
    # Auxiliary variables.
    orb = orbp.orb
    tbd = orbp.tbd

    # Propagate the orbit.
    (r_i, v_i) = twobody!(tbd, tbd.Δt + Δt)

    # Update the elements in the `orb` structure.
    orb.t += Δt/86400
    orb.f  = tbd.f_k

    # Return the information about the step.
    (copy(orbp.orb), r_i, v_i)
end


text:
# Representation of MOI functions using SparseArrays
# Might be able to replace in the future by a function in MOI, see
# https://github.com/jump-dev/MathOptInterface.jl/pull/1238
function sparse_array_representation(terms::Vector{MOI.ScalarAffineTerm{T}}, num_variables, index_map) where {T}
    n = length(terms)
    indices = Vector{Int64}(undef, n)
    coefficients = Vector{T}(undef, n)
    for i in eachindex(terms)
        term = terms[i]
        indices[i] = index_map[term.variable].value
        coefficients[i] = term.coefficient
    end
    return sparsevec(indices, coefficients, num_variables)
end
function sparse_array_representation(terms::Vector{MOI.ScalarQuadraticTerm{T}}, num_variables, index_map) where {T}
    n = length(terms)
    I = Vector{Int64}(undef, n)
    J = Vector{Int64}(undef, n)
    V = Vector{T}(undef, n)
    for k in eachindex(terms)
        term = terms[k]
        i = index_map[term.variable_1].value
        j = index_map[term.variable_2].value
        I[k] = i
        J[k] = j
        V[k] = term.coefficient
        if i != j
            push!(I, j)
            push!(J, i)
            push!(V, term.coefficient)
        end
    end
    return sparse(I, J, V, num_variables, num_variables)
end
struct SparseScalarAffineFunction{T}
    terms::SparseVector{T,Int64}
    constant::T
end
function sparse_array_representation(func::MOI.ScalarAffineFunction, num_variables, index_map)
    return SparseScalarAffineFunction(
        sparse_array_representation(func.terms, num_variables, index_map),
        func.constant,
    )
end
struct SparseScalarQuadraticFunction{T}
    quadratic_terms::SparseMatrixCSC{T,Int64}
    affine_terms::SparseVector{T,Int64}
    constant::T
end
function sparse_array_representation(func::MOI.ScalarQuadraticFunction, num_variables, index_map)
    return SparseScalarQuadraticFunction(
        sparse_array_representation(func.quadratic_terms, num_variables, index_map),
        sparse_array_representation(func.affine_terms, num_variables, index_map),
        func.constant,
    )
end
_convert(::Type{F}, ::Nothing) where {F} = zero(F)
_convert(::Type{F}, obj) where {F} = convert(F, obj)

struct IdentityMap end
Base.getindex(::IdentityMap, index) = index
function sparse_array_representation(func::MOI.AbstractFunction, num_variables)
    return sparse_array_representation(func, num_variables, IdentityMap())
end


text:
using PatModules
@auxmodule two ("utils",) begin

    function say_bye()
        utils.say("bye")
    end
end


text:
function construct_network!(canonical::Canonical, sys::PSY.System,
                            system_formulation::Type{CopperPlatePowerModel};
                            kwargs...)

    buses = PSY.get_components(PSY.Bus, sys)
    bus_count = length(buses)

    copper_plate(canonical, :nodal_balance_active, bus_count)

    return
end

function construct_network!(canonical::Canonical, sys::PSY.System,
                            system_formulation::Type{StandardPTDFModel};
                            kwargs...)

    if :PTDF in keys(kwargs)
        buses = PSY.get_components(PSY.Bus, sys)
        ac_branches = PSY.get_components(PSY.ACBranch, sys)
        ptdf_networkflow(canonical,
                         ac_branches,
                         buses,
                         :nodal_balance_active,
                         kwargs[:PTDF])

        dc_branches = PSY.get_components(PSY.DCBranch, sys)
        dc_branch_types = typeof.(dc_branches)
        for btype in Set(dc_branch_types)
            typed_dc_branches = IS.FlattenIteratorWrapper(btype, Vector([[b for b in dc_branches if typeof(b) == btype]]))
            flow_variables!(canonical,
                           StandardPTDFModel,
                           typed_dc_branches)
        end

    else
        throw(ArgumentError("no PTDF matrix supplied"))
    end

    return

end

function construct_network!(canonical::Canonical,
                            sys::PSY.System,
                            ::Type{T};
                            kwargs...) where {T<:PM.AbstractPowerModel}

    incompat_list = [PM.SDPWRMPowerModel,
                     PM.SparseSDPWRMPowerModel,
                     PM.SOCBFPowerModel,
                     PM.SOCBFConicPowerModel]

    if T in incompat_list
       throw(ArgumentError("$(T) formulation is not currently supported in PowerSimulations"))
    end

    powermodels_network!(canonical, T, sys)
    add_pm_var_refs!(canonical, T, sys)

    return

end


text:
using RandomFerns
using Base.Test

test = ["regression_test.jl"]

println("Running tests...")
for curtest in test
  println("Test: $curtest")
  include(curtest)
end


text:
using Econometrics, Plots
function GridExample(points, doprint=false)
    # plot the line
    x = range(-pi,stop=pi/2.0,length=1000)
    obj = theta-> 0.5*theta*sin(theta^2.0)
    y = obj.(x)
    plot(x, y, legend=false)
    # plot the grid points
    x = range(-pi,stop=pi/2.0,length=points)
    y = obj.(x)
    scatter!(x, y)
    # plot the best point found
    scatter!([x[argmin(y)]], [y[argmin(y)]], markersize=10)
    p = plot!(legend=false)
    if doprint savefig("gridsearch.png") end
    return p
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
@testset "Node pop plots" begin

excelPopPlot( mpSim, "excelrequest/plotRequest" )

end  # @testset "Node pop plots"


text:
using Documenter
using Pushover

makedocs(
    format = :html,
    sitename = "Pushover.jl",
    pages = [
        "index.md",
    ]
)

deploydocs(
    repo = "github.com/scls19fr/Pushover.jl.git"
)


text:
function LTT(tree::Tree)
    function aux(subtree, node_times, events, current_depth)
        if subtree.n_nodes < 2
            pushfirst!(node_times, current_depth + subtree.branch_length)
            if subtree.extant
                pushfirst!(events, 0)
            else
                pushfirst!(events, -1)
            end
        else
            aux(subtree.offsprings[2],node_times,events, current_depth + subtree.branch_length)
            aux(subtree.offsprings[1],node_times,events, current_depth + subtree.branch_length)
            pushfirst!(node_times,current_depth+ subtree.branch_length)
            pushfirst!(events,1)
        end
    end

    node_times = Array{Float64,1}(undef,0)
    events = Array{Int64,1}(undef,0)
    aux(tree, node_times, events, 0.)
    indices = sortperm(node_times)
    events = events[indices]
    node_times = node_times[indices]
    if length(events) > 1
        while events[end-1] == 0
            pop!(events)
            pop!(node_times)
        end
    end
    events[1] += 1
    events = cumsum(events, dims = 1)

    return node_times, events
end

function LTT(tree::Tree, times::Array{Float64,1})
    n = length(times)

    function aux(subtree::Tree, from_root::Float64, id_time::Int64, ltt::Array{Int64,1})
        if id_time <= n
            if (subtree.branch_length + from_root) <= times[id_time]
                if subtree.n_nodes < 2
                    if !subtree.extant
                        ltt[id_time] += -1
                    end
                else
                    ltt[id_time] += 1
                    aux(subtree.offsprings[2],from_root+ subtree.branch_length, id_time, ltt)
                    aux(subtree.offsprings[1],from_root+ subtree.branch_length, id_time, ltt)
                end
            else
                aux(subtree,from_root, id_time+1, ltt)
            end
        end

    end

    ltt = fill(0, length(times))
    ltt[1] = 1
    aux(tree, 0., 1, ltt)
    ltt = cumsum(ltt, dims = 1)

    return times, ltt

end

function LTT(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, times::Array{Float64,1})
    n = length(times)

    function aux(subtree::Tree, from_root::Float64, id_time::Int64, ltt::Array{Int64,1}, scale::Float64)
        if id_time <= n
            if (subtree.branch_length/scale + from_root) <= times[id_time]
                if subtree.n_nodes < 2
                    if !subtree.extant
                        ltt[id_time] += -1
                    end
                else
                    ltt[id_time] += 1
                    aux(subtree.offsprings[2],from_root+ subtree.branch_length/scale, id_time, ltt, scale)
                    aux(subtree.offsprings[1],from_root+ subtree.branch_length/scale, id_time, ltt, scale)
                end
            else
                aux(subtree,from_root, id_time+1, ltt, scale)
            end
        end

    end

    ltt = fill(0, length(times))
    ltt[1] = 1
    aux(tree, 0., 1, ltt, 1.)
    #println("1 $ltt ")
    i=0
    md = edge_trees[1].stem_depth
    for et in edge_trees# in 1:length(edge_trees)
        if et.tree.n_nodes > 1
            aux(et.tree, md - et.stem_depth, 1, ltt,et.scale)
        end
    end
    #println("2 $ltt ")
    ltt = cumsum(ltt, dims = 1)
    #println("3 $ltt ")

    return times, ltt
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
#=
read_json_url_julia:
- Julia version: 1.1.0 (2019-01-21)
- Author: arturogonzalez
- Date: 2019-04-22
=#
#=
using Pkg
Pkg.add("HTTP")
Pkg.add("JSON")
=#
using HTTP, JSON

function main()
    """
    :return: List of places and magnitudes, where magnitude is greater than 1.0.
    """
    resp = HTTP.get("https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.geojson")
    str = String(resp.body)
    data = JSON.Parser.parse(str)
    features_data = data["features"]
    # JSON.print(features_data, 2)
    for i in features_data
        place = i["properties"]["place"]
        magnitude = i["properties"]["mag"]
        if magnitude > 1.0
            println(place, " | ", magnitude)
            end
    end
end;
@time main()


text:
export Triangle

# ------------ #
# - Triangle - #
# ------------ #

"""
    Triangle

Triangle is a primitive object. Any complex object can be represented as a mesh
of Triangles.

### Fields:

* `v1`       - Vertex 1
* `v2`       - Vertex 2
* `v3`       - Vertex 3
* `material` - Material of the Triangle
"""
struct Triangle{V, P, Q, R, S, T, U} <: Object
    v1::Vec3{V}
    v2::Vec3{V}
    v3::Vec3{V}
    material::Material{P, Q, R, S, T, U}
end

show(io::IO, t::Triangle) =
    print(io, "Triangle Object:\n    Vertex 1 - ", t.v1, "\n    Vertex 2 - ", t.v2,
          "\n    Vertex 3 - ", t.v3, "\n    ", t.material)

@diffops Triangle

get_intersections_triangle(a, b, c, d) =
    (a > 0 && b > 0 && c > 0 && d > 0) ? a : bigmul(a)

Zygote.@adjoint function get_intersections_triangle(a::T, b::T, c::T, d::T) where {T}
    res = get_intersections_triangle(a, b, c, d)
    function ∇get_intersections_triangle(Δ)
        if a > 0 && b > 0 && c > 0 && d > 0
            return (Δ, T(0), T(0), T(0))
        else
            return (T(0), T(0), T(0), T(0))
        end
    end
    return res, ∇get_intersections_triangle
end

function intersect(t::Triangle, origin, direction)
    normal = normalize(cross(t.v2 - t.v1, t.v3 - t.v1))
    h = (-dot(normal, origin) .+ dot(normal, t.v1)) ./ dot(normal, direction)
    pt = origin + direction * h
    edge1 = t.v2 - t.v1
    edge2 = t.v3 - t.v2
    edge3 = t.v1 - t.v3
    c₁ = pt - t.v1
    c₂ = pt - t.v2
    c₃ = pt - t.v3
    val1 = dot(normal, cross(edge1, c₁))
    val2 = dot(normal, cross(edge2, c₂))
    val3 = dot(normal, cross(edge3, c₃))
    result = broadcast(get_intersections_triangle, h, val1, val2, val3)
    return result
end

function get_normal(t::Triangle, pt, dir)
    # normal not expanded
    normal_nexp = normalize(cross(t.v2 - t.v1, t.v3 - t.v1))
    direction = -sign.(dot(normal_nexp, dir))
    normal = Vec3(repeat(normal_nexp.x, inner = size(pt.x)),
                  repeat(normal_nexp.y, inner = size(pt.y)),
                  repeat(normal_nexp.z, inner = size(pt.z)))
    return normal * direction
end


text:
function lintrig2quadtrig(mesh::GeoMesh)
    geomesh = GeoMesh()
    n_nodes = length(mesh.nodes)

    nn_dict = Dict{(Int, Int), Int}()

    for node in mesh.nodes
        push!(geomesh, node)
    end
    combs = [[1 2]; [2 3]; [3 1]]

    n_eles = 0
    for lintrig in mesh.elements
        vs = lintrig.vertices
        vs_quad = [vs[1],vs[2],vs[3]]
        for i in 1:size(combs,1)
            c = combs[i,:]
            n1 = mesh.nodes[vs[c[1]]]
            n2 = mesh.nodes[vs[c[2]]]
            node_numbs = sort([n1.n, n2.n])

            if !haskey(nn_dict, (node_numbs[1], node_numbs[2]))
                n_nodes += 1
                nn_dict[(node_numbs[1], node_numbs[2])] = n_nodes
                c1 = get_coord(n1)
                c2 = get_coord(n2)

                node = GeoNode2(n_nodes, Point2((c1.x+c2.x)/2, (c1.y+c2.y)/2))
                push!(geomesh, node)
            end
            node_n =  nn_dict[(node_numbs[1], node_numbs[2])]
            push!(vs_quad, node_n)
        end
        n_eles += 1
        push!(geomesh, GeoQTrig(n_eles, vs_quad))
    end
    return geomesh
end


text:
export @xfail

immutable XFail
  fn::Function
end

macro xfail(test::Expr)
  const name = test.args[1].args[1]
  const body = test.args[2]
  quote
    const $(esc(name)) = XFail() do
                           $(esc(body))
                         end
  end
end

function push_test!(tests::Vector{Tuple{AbstractString, Function}}, name::AbstractString, test::XFail)
  function run_test()
    try
      test.fn()
      return false, :green, "$name XPASSED"
    catch
      return true, :red, "$name XFAILED"
    end
  end
  push!(tests, (name, run_test))
end


text:
test1 = ex1()
#nlp = convert(AbstractNLPModel, test1)
nlp = NLMPCC(test1)
x0 = [1.5, 1.5]

@test nlp.mod.meta.ncc  == 1
@test nlp.mod.meta.ncon == 1
@test obj(nlp, nlp.meta.x0) == 0.0
@test grad(nlp, nlp.meta.x0) == grad(test1,[Inf,Inf])
@test hess(nlp, nlp.meta.x0) == sparse(zeros(2,2))

nlp_at_x = NLPAtX(x0, zeros(4))
stop_2 = NLPStopping(nlp, nlp_at_x, optimality_check = KKT)

printstyled("1st scenario:\n")
solveIpopt(stop_2)
@show stop_2.current_state.x, status(stop_2)

stop_s = MPCCStopping(test1, MPCCAtX(x0, zeros(1)), optimality_check = SStat)

printstyled("2nd scenario:\n")
solveIpopt(stop_s)
@show stop_s.current_state.x, status(stop_s)

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

test2 = ex3()
#nlp = convert(AbstractNLPModel, test1)
nlp = NLMPCC(test2)
x0 = [1.5, 1.5, 1.5]

@test nlp.mod.meta.ncc  == 1
@test nlp.mod.meta.ncon == 2

nlp_at_x = NLPAtX(x0, zeros(5))
stop_3 = NLPStopping(nlp, nlp_at_x, optimality_check = KKT)

printstyled("1st scenario:\n")
solveIpopt(stop_3)
@show stop_3.current_state.x, status(stop_3)

stop_s = MPCCStopping(test2, MPCCAtX(x0, zeros(2)), optimality_check = MStat)

printstyled("2nd scenario:\n")
solveIpopt(stop_s)
@show stop_s.current_state.x, status(stop_s)


text:
using Pkg
Pkg.activate(".")
using SignalBroadcastingSim
using Random
using OpenStreetMapX
using DelimitedFiles, CSV
using LightGraphs
using DataFrames
using Statistics
using CSV
using Plots

m = OpenStreetMapX.get_map_data("maps", "torontoF.osm")

p = ModelParams(n_agents=0, discretize_m=50)
s = Simulation(p,m)
SignalBroadcastingSim.plot_sim(s, "road_network.html")

Random.seed!(20)
p = ModelParams(n_agents=500, discretize_m=50)
s = Simulation(p,m)
init!(s)
for i in 1:100
    SignalBroadcastingSim.step!(s)
end
#SignalBroadcastingSim.
plot_sim(s, "sample_200agents_discretize50m_after100steps.html")

Random.seed!(20)
p = ModelParams(n_agents=10000, discretize_m=50)
s = Simulation(p,m)
init!(s)
SignalBroadcastingSim.simulate_steps!(s;max_idle=50000)

steps = length(s.step_infected_agents_count)

pgfplotsx()

plt = plot(lab="", xlim=[0,100],
    legend=:bottomright,
    xlabel = "Simulation step number",
    ylabel="Percentage of agents that received the message" )

median_x = findfirst(>=(0.5), s.step_infected_agents_count./p.n_agents)
plot!(plt, 1:steps,s.step_infected_agents_count./p.n_agents,
    lab="Agents that received the message", color=:blue)

bar!(1:(median_x-1),(s.step_infected_agents_count./p.n_agents)[1:median_x-1],color=:blue,linecolor=nothing, lab="", fillalpha=0.1)
bar!(median_x:steps,ones(steps-median_x+1),color=:red,linecolor=nothing, lab="",fillalpha=0.1)
bar!(median_x:steps,(s.step_infected_agents_count./p.n_agents)[median_x:steps],color=:white,linecolor=:white, lab="")

plot!(plt, 1:steps,s.step_infected_agents_count./p.n_agents,lab="", color=:blue)

annotate!(plt, median_x+24, 0.5, text("Median value (here it is at the step $median_x)", 10, :black))

annotate!(plt, median_x*0.6, 0.1, text("Left tail", 8, :blue))

annotate!(plt, median_x*1.1, .9, text("Upper part of the right tail", 8, :red))

scatter!(plt, [median_x-0.75], [0.5],
    markershape=:star5,
    markerstrokecolor=:black,
lab="")

function savefig2(plot,filename; twidth=raw"\columnwidth", replace_y_lab_pos=nothing)
    savefig(plot,filename)
    dat = read(filename, String)
    open(filename, "w") do f
        println(f,raw"\resizebox{"*twidth*"}{!}{%")
        if replace_y_lab_pos != nothing
            dat = replace(dat, "ylabel style={"=>"ylabel style={at={(axis description cs:$(replace_y_lab_pos))}")
        end

        println(f, dat)
        println(f,raw"}")
    end
end

savefig2(plt,"../5d63f3b5c3791641ba87e19f/sample_simulation_run.tex";twidth=raw"0.9\columnwidth")

plt = plot(lab="", xlim=[0,100],
    legend=:bottomright,
    xlabel = "Simulation step number",
    ylabel="Number of agents that received the message",
    size=(400,300) );

plot!(plt, 1:steps,s.step_infected_agents_count,
    lab="Agents that received the message", color=:blue);

savefig(plt,"../5d63f3b5c3791641ba87e19f/sample_simulation_run_simple.pdf")