text:
#TODO: Add Keyword Constructors for all datatypes
#TODO MAYBE: Create + methods to build cells
nca_default() = struct_electrode(
                                  "NCA",                # Chemistry type
                                  0.03,                 # binder weight fraction
                                  0.03,                 # conductive material weight fraction
                                  4.71,                 # AM density
                                  1.77,                 # binder density
                                  1.77,                 # conductive material density
                                  192.0,                # first charge specific capacity
                                  180.0,                # reversible specific capacity
                                  150.0,                # total thickness cathode
                                  0.2,                  # porosity
                                  0.0,                  # extra length
                                  15.0,                 # Pos foil (Al) thickness                      microns
                                  2.7,                  # Positive foil density
                                  7.0,                  # Tab thickness                                microns
                                  6.0,                  # Tab width                                      mm
                                  8.88,                 # Tab density (Ni)                              g/cm3
                                )

nmc622_default() = struct_electrode(
                                    "NMC622",     # Chemistry type
                                    0.03,         # binder weight fraction
                                    0.03,         # Conductive material weight fraction
                                    4.8,          # AM density
                                    1.77,         # binder density
                                    1.77,         # conductive material density
                                    190.0,        # first charge specific capacity
                                    171.0,        # reversible specific capactity
                                    150.0,        # total thickness cathode
                                    0.2,          # porosity
                                    0.0,          # extra length
                                    15.0,         # Pos foil (Al) thickness                      microns
                                    2.7,          # Pos foil density
                                    7.0,          # Tab thickness                                microns
                                    6.0,          # Tab width                                      mm
                                    8.88,         # Tab density
                                  )

lfp_default() = struct_electrode(
                                  "LFP",        # Chemistry type
                                  0.03,         # binder weight fraction
                                  0.03,         # conductive material weight fraction
                                  3.6,          # AM density
                                  1.77,         # binder density
                                  1.77,         # conductive material density
                                  160.0,        # first charge specific capacity
                                  150.0,        # reversible specific capacity
                                  150.0,        # total thickness cathode
                                  0.2,          # porosity
                                  0.0,          # extra length
                                  15.0,         # Pos foil (Al) thickness                      microns
                                  2.7,          # Pos foil density
                                  7.0,          # Tab thickness                                microns
                                  6.0,          # Tab width                                      mm
                                  8.88,         # Tab density
                                )

graphite_default() = struct_electrode(

                                  "Graphite",               # Chemistry type

                                  0.03,                     # binder weight fraction
                                  0.03,                     # conductive weight fraction

                                  2.24,                     # AM density
                                  1.77,                     # binder density
                                  1.77,                     # conductive density

                                  371.93,                   # first charge specific capacity
                                  350.0,                    # reversible specific capacity

                                  150.0,                    # total thickness cathode
                                  0.3,                      # porosity

                                  40.0,                     # Extra length                                   mm

                                  8.0,                      # Neg foil (Cu) thickness                      microns
                                  8.92,                     # Neg foil (Cu) density                         g/cm3

                                  7.0,                      # Tab thickness                                microns
                                  6.0,                      # Tab width                                      mm
                                  8.88,                     # Tab density (Ni)                              g/cm3
                                  )

cell(;form_factor="Cyl", size="18650", dimensional_delta=0.01,
    can_th=0.05, can_rho=7.7, sep_th=20.0, sep_rho=1.2,
    sep_por=0.4, v_nom=3.6, electrolyte_name="EC:DMC",
    electrolyte_rho=1.2, NP=1.15,
    anode=graphite_default(), cathode=nca_default()) = struct_cell_general(
                                      form_factor,              # Form factor                           ["Cyl","Pris","Pou"]
                                       size,                    # Cell size                                    5 digits
                                       dimensional_delta,       # dimensional_delta                               %
                                       can_th,                  # canister thickness                             (mm)
                                       can_rho,                 # Stainless Steel density                        g/cm3

                                       sep_th,                  # Separator thickness                          microns
                                       sep_rho,                 # Separator density                             g/cm3
                                       sep_por,                 # Separator porosity

                                       v_nom,                   # Nominal voltage                                 V

                                       electrolyte_name,        # Electrolyte name
                                       electrolyte_rho,         # Electrolyte density                           g/cm3

                                       NP,                      # N to P ratio

                                       anode,                   # Anode Struct
                                       cathode,                 # Cathode Struct
                                       )

cathode_default_cost() = struct_electrode_costs(
                                                [24.0, 0.9],            # Cost of Positive Active Material                       # USD/kg
                                                [6.6, 1.0],             # Cost of Positive Conductive                            # USD/kg
                                                [9.5, 1.0],             # Cost of Positive Binder                                # USD/kg
                                                [3.1, 1.0],             # Cost of Positive Binder Solvent                        # USD/kg
                                                [10193620.403506],      # Baseline annual positive active material               # kg/year
                                                [13744207.2856261],     # Baseline annual positive binder solvent evaporated     # kg/year
                                                [96.0 / 4.0],           # Multiplier for binder solvent requirement
                                                [0.3, 1.0],             # Cost of Positive Current Collector                     # USD/m2
                                                [0.15, 0.8],            # Cost of Positive Terminal Assembly                     # USD/cell
                                                [0.922],                # Positive Electrode Material (dry) yield
                                                [0.902],                # Positive Current Collector yield
                                              )

anode_default_cost() = struct_electrode_costs(
                                              [12.5, 0.95],                # Cost of Negative Active Material                       # USD/kg
                                              [6.6, 1.0],                  # Cost of Negative Conductive                            # USD/kg
                                              [10.0, 1.0],                 # Cost of Negative Binder                                # USD/kg
                                              ########### Negative binder solvent cost - water
                                              [0.0, 0.9],                  # !!!!!!!!!!! Check !!!!!!!!!!!!!!!!!
                                              [6564082.00851682],          # Baseline annual negative active material               # kg/year
                                              [8291472.01075809],          # Baseline annual negative binder solvent evaporated     # kg/year
                                              [96.0 / 4.0],                # Multiplier for binder solvent requirement
                                              [1.2, 1.0],                  # Cost of Negative Current Collector                     # USD/m2
                                              [0.15, 0.8],                 # Cost of Negative Terminal Assembly                     # USD/cell
                                              [0.922],                     # Negative Electrode Material (dry) yield
                                              [0.902],                     # Negative Current Collector yield
                                            )

cell_costs_default() = struct_cell_costs(
                                          [0.22, 0.8],           # Cost of Cell Canister                                  # USD/cell
                                          [1.1, 1.0],            # Cost of Separator                                      # USD/m2
                                          0.98,                  # Separator yield
                                          [15.0, 1.0],           # Cost of Electrolyte                                    # USD/L
                                          0.94,                  # Electrolyte yield
                                          0.95,                  # Cell yield
                                          cathode_default_cost(),
                                          anode_default_cost(),
                                        )

baseline_cell_assembly_data_default() = struct_manufacturing(
                                                              # Receiving
                                                              [
                                                                14400.0 0.4           # Direct Labor
                                                                3.6 0.6               # Capital Equipment
                                                                900.0 0.5             # Plant Area
                                                                1.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # pos_mat_prep
                                                              [
                                                                14400.0 0.5           # Direct Labor
                                                                2.0 0.7               # Capital Equipment
                                                                600.0 0.6             # Plant Area
                                                                2.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # neg_mat_prep
                                                              [
                                                                14400.0 0.5           # Direct Labor
                                                                2.0 0.7               # Capital Equipment
                                                                600.0 0.6             # Plant Area
                                                                3.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # pos_electrode_coat
                                                              [
                                                                28800.0 0.5           # Direct Labor
                                                                8.0 0.8               # Capital Equipment
                                                                750.0 0.8             # Plant Area
                                                                4.0 1.0               # [volume ratio, special case]
                                                              ],

                                                              # neg_electrode_coat
                                                              [
                                                                28800.0 0.5           # Direct Labor
                                                                8.0 0.8               # Capital Equipment
                                                                750.0 0.8             # Plant Area
                                                                4.0 1.0               # [volume ratio, special case]
                                                              ],

                                                              # bind_slvnt_recov
                                                              [
                                                                14400.0 0.4           # Direct Labor
                                                                3.0 0.6               # Capital Equipment
                                                                225.0 0.6             # Plant Area
                                                                5.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # pos_calendr
                                                              [
                                                                14400.0 0.5           # Direct Labor
                                                                1.0 0.7               # Capital Equipment
                                                                225.0 0.6             # Plant Area
                                                                4.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # neg_calendr
                                                              [
                                                                7200.0 0.5            # Direct Labor
                                                                1.0 0.7               # Capital Equipment
                                                                225.0 0.6             # Plant Area
                                                                4.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # intr_process_mat_handling
                                                              [
                                                                28800.0 0.7           # Direct Labor
                                                                1.5 0.7               # Capital Equipment
                                                                900.0 0.6             # Plant Area
                                                                4.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # electrode_slitting
                                                              [
                                                                28800.0 0.5           # Direct Labor
                                                                2.0 0.7               # Capital Equipment
                                                                300.0 0.6             # Plant Area
                                                                4.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # vacuum_drying
                                                              [
                                                                14400.0 0.5           # Direct Labor
                                                                1.6 0.7               # Capital Equipment
                                                                300.0 0.6             # Plant Area
                                                                4.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # control_lab
                                                              [
                                                                28800.0 0.5           # Direct Labor
                                                                1.5 0.7               # Capital Equipment
                                                                300.0 0.6             # Plant Area
                                                                1.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # cell_winding
                                                              [
                                                                36000.0 0.7           # Direct Labor
                                                                4.0 0.8               # Capital Equipment
                                                                600.0 0.8             # Plant Area
                                                                6.0 2.0               # [volume ratio, special case]
                                                              ],

                                                              # CC_welding
                                                              [
                                                                36000.0 0.7           # Direct Labor
                                                                4.0 0.8               # Capital Equipment
                                                                600.0 0.8             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # cell_insertion_container
                                                              [
                                                                21600.0 0.5           # Direct Labor
                                                                3.0 0.7               # Capital Equipment
                                                                600.0 0.6             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # electrolyte_filling_and_cell_sealing
                                                              [
                                                                36000.0 0.5           # Direct Labor
                                                                5.0 0.7               # Capital Equipment
                                                                900.0 0.6             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # dry_room_control
                                                              [
                                                                14400.0 0.4           # Direct Labor
                                                                20.0 0.6              # Capital Equipment
                                                                100.0 0.4             # Plant Area
                                                                7.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # formation_cycling
                                                              [
                                                                57600.0 0.7           # Direct Labor
                                                                30.0 0.8              # Capital Equipment
                                                                2200.0 0.8            # Plant Area
                                                                6.0 2.0               # [volume ratio, special case]
                                                              ],

                                                              # cell_sealing
                                                              [
                                                                14400.0 0.5           # Direct Labor
                                                                2.0 0.7               # Capital Equipment
                                                                450.0 0.6             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # charge_retention_testing
                                                              [
                                                                21600.0 0.4           # Direct Labor
                                                                4.75 0.7              # Capital Equipment
                                                                900.0 0.6             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # rejected_cell_scrap_recycle_assembly
                                                              [
                                                                36000.0 0.7           # Direct Labor
                                                                2.5 0.7               # Capital Equipment
                                                                600.0 0.6             # Plant Area
                                                                6.0 0.0               # [volume ratio, special case]
                                                              ],

                                                              # Shipping
                                                              [
                                                                28800.0 0.5           # Direct Labor
                                                                5.0 0.7               # Capital Equipment
                                                                900.0 0.6             # Plant Area
                                                                1.0 0.0               # [volume ratio, special case]
                                                              ],
                                                            )

costs_data_default() = struct_general_costs(
                                            715.0,                       # Pack Integration Cost                                  # USD/Pack
                                            3000.0,                      # Building cost                                          #  USD/m2
                                            18.0,                        # Labor cost                                             # $/hr
                                            150.0,                       # No of units mfg per energy_kwh_per_year                Million
                                            0.995,                       # Positive Binder Solvent Recovery percentage
                                            baseline_cell_assembly_data_default(),
                                            1.1,                         # formation_cycling_mult
                                            1.1,                         # cell_stacking_mult
                                            0.2,                         # pos_electrode_coating_solvent_SR
                                            0.2,                         # neg_electrode_coating_solvent_SR
                                          )

baseline_data_default() = struct_baseline(
                                          6000000.0,               # Baseline annual energy                                   # kWh/year
                                          25263158.0,              # Baseline annual cells adjusted for yield                 # cells/year
                                          67669329.6866343,        # Baseline annual electrode area                           # m2/year
                                          6100.0,
                                        )

battery_pack_data_default() = struct_pack_level(240.0,)

ovhd_rate_data_default() = struct_ovhd_rate(
                                            5.0 / 100.0,                      # Launch cost percent direct annual materials
                                            10.0 / 100.0,                     # Launch cost percent direct labor plus variable overhead
                                            16.666667 / 100.0,                # Depreciation rate capital equipment                          # 6-year rate
                                            5.0 / 100.0,                      # Depreciation rate building investment                        # 20 year rate
                                            25.0 / 100.0,                     # General, Sales and Administration rate direct labor plus variable overhead
                                            25.0 / 100.0,                     # General, Sales and Administration rate depreciation
                                            40.0 / 100.0,                     # R & D rate depreciation
                                            40.0 / 100.0,                     # Variable overhead rate direct labor
                                            20.0 / 100.0,                     # Variable overhead rate depreciation
                                            5.0 / 100.0,                      # Profit rate after taxes                                      # % of investment
                                            15.0 / 100.0,                     # Working capital rate                                         # % of annual variable cost
                                            5.6 / 100.0,                      # Warranty cost added to the price rate
                                          )

cost_default() = struct_cost(
  cell_costs_default(),
  baseline_cell_assembly_data_default(),
  costs_data_default(),
  baseline_data_default(),
  ovhd_rate_data_default(),
)

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
macro gpu(_ex); if gpu()>=0; esc(_ex); end; end

macro cuda1(lib,fun,x...)       # return -1 if library missing, error code if run
    if Libdl.find_library(["lib$lib"], []) != ""
        Expr(:ccall, ("$fun","lib$lib"), :UInt32, x...)
    else
        -1
    end
end

macro cuda(lib,fun,x...)        # give an error if library missing, warning if return!=0
    if Libdl.find_library(["lib$lib"], []) != ""
        fx = Expr(:ccall, ("$fun","lib$lib"), :UInt32, x...)
        msg = "$lib.$fun error "
        err = gensym()
        esc(:(if ($err=$fx) != 0; warn($msg, $err); Base.show_backtrace(STDOUT, backtrace()); end))
    else
        Expr(:call,:error,"Cannot find lib$lib, please install it and rerun Pkg.build(\"Knet\").")
    end
end

typealias Cptr Ptr{Void}

"""

`gpu()` returns the id of the active GPU device or -1 if none are
active.

`gpu(true)` resets all GPU devices and activates the one with the most
available memory.

`gpu(false)` resets and deactivates all GPU devices.

`gpu(d::Int)` activates the GPU device `d` if `0 <= d < gpuCount()`,
otherwise deactivates devices.

`gpu(true/false)` resets all devices.  If there are any allocated
KnetArrays their pointers will be left dangling.  Thus
`gpu(true/false)` should only be used during startup.  If you want to
suspend GPU use temporarily, use `gpu(-1)`.

`gpu(d::Int)` does not reset the devices.  You can select a previous
device and find allocated memory preserved.  However trying to operate
on arrays of an inactive device will result in error.

"""
function gpu end

let GPU=-1, GPUCNT=-1, CUBLAS=nothing, CUDNN=nothing
    global gpu, gpuCount, cublashandle, cudnnhandle, cudaRuntimeVersion, cudaDriverVersion

    gpu()=GPU

    function gpuCount() # should not bomb when there is no gpu or nvidia libs
        if GPUCNT == -1
            GPUCNT = try
	        p=Cuint[0]
                # @cuda does not stay quiet so we use @cuda1 here
                # This code is only run once if successful, so nvmlInit here is ok
                @cuda1("nvidia-ml",nvmlInit,())
                @cuda1("nvidia-ml",nvmlDeviceGetCount,(Ptr{Cuint},),p)
                # Let us keep nvml initialized for future ops such as meminfo
                # @cuda1("nvidia-ml",nvmlShutdown,())
	        Int(p[1])
            catch
	        0
            end
        end
        return GPUCNT
    end

    function gpu(i::Int)
        (GPU == i) && return i
        if 0 <= i < gpuCount()
            @cuda(cudart,cudaSetDevice, (Cint,), i)
            cudaRuntimeVersion = (p=Cint[0];@cuda(cudart,cudaRuntimeGetVersion,(Ptr{Cint},),p);Int(p[1]))
            cudaDriverVersion  = (p=Cint[0];@cuda(cudart,cudaDriverGetVersion, (Ptr{Cint},),p);Int(p[1]))
        else
            i = -1
            # @cuda(cudart,cudaDeviceReset,()) # may still go back and use arrays allocated in a previous gpu
        end
        return (GPU = i)
    end

    function gpu(usegpu::Bool)
        if usegpu && gpuCount() > 0
            pick = free = same = -1
            for i=0:gpuCount()-1
                mem = nvmlDeviceGetMemoryInfo(i)
                if mem[2] > free
                    pick = i
                    free = mem[2]
                    same = 1
                elseif mem[2] == free
                    # pick one of equal devices randomly
                    rand(1:(same+=1)) == 1 && (pick = i)
                end
            end
            gpu(pick)
        else
            for i=0:gpuCount()-1
                @cuda(cudart,cudaDeviceReset,())
            end
            gpu(-1)
        end
    end

    function cublashandle(dev=gpu())
        if dev==-1; error("No cublashandle for CPU"); end
        i = dev+2
        if CUBLAS == nothing; CUBLAS=Array(Any,gpuCount()+1); end
        if !isassigned(CUBLAS,i); CUBLAS[i]=cublasCreate(); end
        return CUBLAS[i]
    end

    function cudnnhandle(dev=gpu())
        if dev==-1; error("No cudnnhandle for CPU"); end
        i = dev+2
        if CUDNN == nothing; CUDNN=Array(Any,gpuCount()+1); end
        if !isassigned(CUDNN,i); CUDNN[i]=cudnnCreate(); end
        return CUDNN[i]
    end
end

# cudaGetDeviceCount is deprecated, use gpuCount instead:
cudaGetDeviceCount()=(try; p=Cint[0]; eval(:(ccall(("cudaGetDeviceCount","libcudart"),UInt32,(Ptr{Cint},),$p))); p[1]; catch; 0; end) # will not bomb when there is no gpu
cudaGetDevice()=(d=Cint[-1];@cuda(cudart,cudaGetDevice,(Ptr{Cint},),d);d[1])
cudaGetMemInfo()=(f=Csize_t[0];m=Csize_t[0]; @cuda(cudart,cudaMemGetInfo,(Ptr{Csize_t},Ptr{Csize_t}),f,m); (Int(f[1]),Int(m[1])))
cudaDeviceSynchronize()=@cuda(cudart,cudaDeviceSynchronize,())

function nvmlDeviceGetMemoryInfo(i=gpu())
    0 <= i < gpuCount() || return nothing
    dev = Cptr[0]
    mem = Array(Culonglong,3)
    @cuda("nvidia-ml","nvmlDeviceGetHandleByIndex",(Cuint,Ptr{Cptr}),i,dev)
    @cuda("nvidia-ml","nvmlDeviceGetMemoryInfo",(Cptr,Ptr{Culonglong}),dev[1],mem)
    ntuple(i->Int(mem[i]),length(mem))
end

function cublasCreate()
    handleP = Cptr[0]
    @cuda(cublas,cublasCreate_v2, (Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cuda(cublas,cublasDestroy_v2, (Cptr,), handle))
    global cublasVersion = (p=Cint[0];@cuda(cublas,cublasGetVersion_v2,(Cptr,Ptr{Cint}),handle,p);Int(p[1]))
    return handle
end

function cudnnCreate()
    handleP = Cptr[0]
    @cuda(cudnn,cudnnCreate,(Ptr{Cptr},), handleP)
    handle = handleP[1]
    atexit(()->@cuda(cudnn,cudnnDestroy,(Cptr,), handle))
    global cudnnVersion = Int(ccall((:cudnnGetVersion,:libcudnn),Csize_t,()))
    return handle
end

ext: jl


text:
@ms using JuMP
@ms using MathOptFormat

# import Julia library corresponding to the chosen solver
if args["external_solver"] == "cbc"
    @ms using Cbc
elseif args["external_solver"] == "glpk"
    @ms using GLPK
else
    @ms using CPLEX
end

"""
Creates an empty JuMP optimization model
"""
function new_model(; solver::String = "auto")
    if solver == "auto"
        solver = args["external_solver"]
    end
    if solver == "glpk"
        model = Model(with_optimizer(GLPK.Optimizer))
    elseif solver == "cplex"
        model = Model(with_optimizer(CPLEX.Optimizer))
    elseif solver == "cbc"
        model = Model(with_optimizer(Cbc.Optimizer))
    end
    return model
end

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Store time series file for a TreeMesh with a DG solver
function save_time_series_file(time_series_callback, mesh::TreeMesh, equations, dg::DG)
  @unpack (interval, solution_variables, variable_names,
           output_directory, filename, point_coordinates,
           point_data, time, step, time_series_cache) = time_series_callback
  n_points = length(point_data)

  h5open(joinpath(output_directory, filename), "w") do file
    # Add context information as attributes
    n_variables = length(variable_names)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["equations"] = get_name(equations)
    attributes(file)["polydeg"] = polydeg(dg)
    attributes(file)["n_vars"] = n_variables
    attributes(file)["n_points"] = n_points
    attributes(file)["interval"] = interval
    attributes(file)["variable_names"] = collect(variable_names)

    file["time"] = time
    file["timestep"] = step
    file["point_coordinates"] = point_coordinates
    for p in 1:n_points
      # Store data as 2D array for convenience
      file["point_data_$p"] = reshape(point_data[p], n_variables, length(time))
    end
  end
end

end # @muladd


text:
gc()

#Tests des packages
include("../ALAS/include.jl")

using MPCCmod
using OutputRelaxationmod

using CUTEst
using NLPModels

problems_unconstrained = CUTEst.select(contype="unc")
problems_linconstrained = CUTEst.select(contype="linear",min_var=1, max_var=100)
problems_boundconstrained = CUTEst.select(contype="bounds")
problems_quadconstrained = CUTEst.select(contype="quadratic")
problems_genconstrained = CUTEst.select(contype="general")
#problems_linconstrained[50] NaN ?
nlp = CUTEstModel(problems_linconstrained[1])
print(nlp)

#déclare le mpcc :
@printf("Create an MPCC \n")
@time exemple_nlp=MPCCmod.MPCC(nlp)

@printf("%i %i %i %i %i %i %i %i \n",nlp.counters.neval_obj,
        nlp.counters.neval_cons,nlp.counters.neval_grad,
        nlp.counters.neval_hess,exemple_nlp.G.counters.neval_cons,
        exemple_nlp.G.counters.neval_jac,exemple_nlp.H.counters.neval_cons
        ,exemple_nlp.H.counters.neval_jac)

@printf("Butterfly method:\n")
#résolution avec ALAS Butterfly
@time xb,fb,orb,nb_eval = MPCCsolve.solve(exemple_nlp)
try
 oa=orb.inner_output_alas[1]
 @printf("%s	",orb.solve_message)
end
@show nb_eval

finalize(nlp)


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    addstereohydrogens, removestereohydrogens,
    setdiastereo!, setdiastereo,
    setstereocenter!, setstereocenter,
    optimizewedges!, optimizewedges

"""
    addstereohydrogens(mol::GraphMol) -> GraphMol

Add stereospecific hydrogen nodes to the molecule.

If the molecule was genereted from SDFile, run `setstereocenter!(mol)` in advance to set stereocenter information.
"""
function addstereohydrogens(mol::GraphMol)
    """
    [C@@H](C)(C)C -> [C@@]([H])(C)(C)C
    C[C@@H](C)C -> C[C@@]([H])(C)C
    """
    atoms = nodeattrs(mol)
    hydrogenconnected_ = hydrogenconnected(mol)
    newmol = graphmol(mol)
    mapper = Dict{Int,Int}()
    offset = 0
    for i in 1:nodecount(mol)
        mapper[i] = i + offset
        atoms[i].stereo === :unspecified && continue
        degree(mol, i) == 3 || continue
        hydrogenconnected_[i] == 1 || continue
        n = addnode!(newmol, nodeattrtype(mol)(:H))
        addedge!(newmol, i, n, edgeattrtype(mol)())
        offset += 1
        mapper[n] = i + offset
    end
    return remapnodes(newmol, mapper)
end

"""
    removestereohydrogens(mol::GraphMol) -> GraphMol

Return new molecule with stereospecific hydrogen nodes removed.

If the molecule was genereted from SDFile, run `setstereocenter!(mol)` in advance to set stereocenter information. This will keep stereochemistry of the molecule object. `addstereohydrogens!` can be used to revert the change.
"""
function removestereohydrogens(mol::GraphMol)
    """
    [C@@]([H])(C)(N)O -> [C@@H](C)(N)O
    ([H])[C@@](C)(N)O -> [C@@H](C)(N)O
    C[C@@]([H])(N)O -> C[C@@H](N)O
    C[C@@](N)([H])O -> C[C@H](N)O (reverse direction)
    C[C@@](N)(O)[H] -> C[C@@H](N)O
    """
    atoms = nodeattrs(mol)
    nodes_to_rm = Int[]
    nodes_to_rev = Int[]
    for c in 1:nodecount(mol)
        atoms[c].stereo === :unspecified && continue
        adjs = adjacencies(mol, c)
        length(adjs) == 4 || continue
        hs = Int[]
        rev = false
        for (i, adj) in enumerate(sort(collect(adjs)))
            if atoms[adj].symbol === :H
                if i == 3
                    rev = true
                end
                push!(hs, adj)
            end
        end
        length(hs) == 1 || continue
        append!(nodes_to_rm, hs)
        rev && push!(nodes_to_rev, c)
    end
    newmol = graphmol(mol)
    for c in nodes_to_rev
        a = setstereo(atoms[c],
            atoms[c].stereo === :clockwise ? :anticlockwise : :clockwise)
        setnodeattr!(newmol, c, a)
    end
    ns = setdiff(nodeset(mol), nodes_to_rm)
    return graphmol(nodesubgraph(newmol, ns))
end

"""
    setdiastereo!(mol::Union{SMILES,SDFile}) -> Nothing

Set diastereomerism flags to `Bond.stereo` fields of double bonds.

`Bond.stereo::Symbol` indicates whether the lower index nodes attached to each sides of the double bond is in `:cis`, `:trans` or `:unspecified` configuration.
"""
function setdiastereo!(mol::SMILES)
    edges = edgeattrs(mol)
    for i in 1:edgecount(mol)
        edges[i].order == 2 || continue
        ds = Symbol[]
        n1, n2 = mol.edges[i]
        for (u, v) in [(n1, n2), (n2, n1)]
            adjs = setdiff(adjacencies(mol, u), [v])
            if length(adjs) == 2
                f, s = sort(collect(adjs))
                fe = findedgekey(mol, u, f)
                se = findedgekey(mol, u, s)
                if edges[fe].direction !== :unspecified
                    if f < u
                        push!(ds, edges[fe].direction === :up ? :down : :up)
                    else
                        push!(ds, edges[fe].direction)
                    end
                elseif edges[se].direction !== :unspecified
                    push!(ds, edges[se].direction === :up ? :down : :up)
                end
            elseif length(adjs) == 1
                f = collect(adjs)[1]
                fe = findedgekey(mol, u, f)
                if edges[fe].direction !== :unspecified
                    if f < u
                        push!(ds, edges[fe].direction === :up ? :down : :up)
                    else
                        push!(ds, edges[fe].direction)
                    end
                end
            end
            # TODO: axial chirality
        end
        if length(ds) == 2
            stereo = ds[1] === ds[2] ? :cis : :trans
            bond = setstereo(edges[i], stereo)
            setedgeattr!(mol, i, bond)
        end
    end
    clearcache!(mol)
end

function setdiastereo!(mol::SDFile)
    nodes = nodeattrs(mol)
    edges = edgeattrs(mol)
    coords_ = sdfcoords2d(mol)
    for i in 1:edgecount(mol)
        edges[i].order == 2 || continue
        edges[i].notation == 3 && continue  # stereochem unspecified
        # Get lower indexed edges connected to each side of the bond
        ns = Int[]
        for n in mol.edges[i]
            incs = incidences(mol, n)
            length(incs) in (2, 3) || continue
            f = sort([e for e in incs if e != i])[1]
            push!(ns, neighbors(mol, n)[f])
        end
        length(ns) == 2 || continue
        # Check coordinates
        d1, d2 = [Point2D(coords_, n) for n in mol.edges[i]]
        n1, n2 = [Point2D(coords_, n) for n in ns]
        cond(a, b) = (d1.x - d2.x) * (b - d1.y) + (d1.y - d2.y) * (d1.x - a)
        n1p = cond(n1.x, n1.y)
        n2p = cond(n2.x, n2.y)
        if n1p * n2p < 0
            stereo = :trans
        elseif n1p * n2p > 0
            stereo = :cis
        else
            stereo = :unspecified
        end
        bond = setstereo(edges[i], stereo)
        setedgeattr!(mol, i, bond)
    end
    clearcache!(mol)
end

"""
    setdiastereo(mol::GraphMol) -> GraphMol

Return new molecule with diastereomeric information set.

See [`setdiastereo!`](@ref))
"""
function setdiastereo(mol::GraphMol)
    newmol = graphmol(mol)
    setdiastereo!(newmol)
    return newmol
end

function angeval(u::Point2D, v::Point2D)
    # 0deg -> 1, 90deg -> 0, 180deg -> -1, 270deg-> -2, 360deg -> -3
    uv = dot(u, v) / (norm(u) * norm(v))
    return cross2d(u, v) >= 0 ? uv : -2 - uv
end

function anglesort(coords, center, ref, vertices)
    c = Point2D(coords, center)
    r = Point2D(coords, ref)
    ps = [Point2D(coords, v) for v in vertices]
    vs = [p - c for p in ps]
    per = sortperm([angeval(r, v) for v in vs])
    return vertices[per]
end

"""
    setstereocenter!(mol::SDFile) -> Nothing

Set stereocenter information to Atom.stereo (`:unspecified`, `:clockwise`, `:anticlockwise` or `:atypical`). Clockwise/anticlockwise means the configuration of 2-4th nodes in index label order when we see the chiral center from the node labeled by the lowest index. If there is an implicit hydrogen, its index label will be regarded as the same as the stereocenter atom.
"""
function setstereocenter!(mol::SDFile)
    nodes = nodeattrs(mol)
    edges = edgeattrs(mol)
    coords_ = sdfcoords2d(mol)
    for i in 1:nodecount(mol)
        nbrs = neighbors(mol, i)
        length(nbrs) < 3 && continue # atoms attached to the chiral center
        adjs = Int[]
        upadjs = Int[]
        downadjs = Int[]
        for (inc, adj) in nbrs
            push!(adjs, adj)
            if edges[inc].notation == 1 && mol.edges[inc][1] == i
                push!(upadjs, adj)
            elseif edges[inc].notation == 6 && mol.edges[inc][1] == i
                push!(downadjs, adj)
            end
        end
        (isempty(upadjs) && isempty(downadjs)) && continue  # unspecified
        if (length(adjs) > 4 ||
            length(adjs) - length(upadjs) - length(downadjs) < 2)
            # Hypervalent, not sp3, axial chirality or wrong wedges
            a = setstereo(nodes[i], :atypical)
            setnodeattr!(mol, i, a)
            continue
        end
        if length(adjs) == 4
            # Select a pivot
            pivot = isempty(upadjs) ? downadjs[1] : upadjs[1]
            tri = setdiff(adjs, [pivot])
            quad = adjs
            rev = isempty(upadjs) ? true : false
            # Check wedges
            if length(downadjs) == 2 || length(upadjs) == 2
                ordered = anglesort(coords_, i, pivot, tri)
                sec = isempty(upadjs) ? downadjs[2] : upadjs[2]
                if findfirst(isequal(sec), ordered) != 2
                    a = setstereo(nodes[i], :atypical)
                    setnodeattr!(mol, i, a)
                    continue
                end
            elseif length(downadjs) == 1 && length(upadjs) == 1
                ordered = anglesort(coords_, i, pivot, tri)
                if findfirst(isequal(downadjs[1]), ordered) == 2
                    a = setstereo(nodes[i], :atypical)
                    setnodeattr!(mol, i, a)
                    continue
                end
            end
        elseif length(adjs) == 3
            # Implicit hydrogen is the pivot
            pivot = i
            tri = adjs
            quad = union(adjs, [i])
            rev = isempty(upadjs) ? false : true
        end
        cw = isclockwise(toarray(coords_, sort(tri)))
        if rev
            cw = !cw
        end
        # Arrange the configuration so that the lowest index node is the pivot.
        pividx = findfirst(isequal(pivot), sort(quad))
        if pividx in (2, 4)
            cw = !cw
        end
        stereo = cw ? :clockwise : :anticlockwise
        a = setstereo(nodes[i], stereo)
        setnodeattr!(mol, i, a)
    end
    clearcache!(mol)
end

"""
    setstereocenter(mol::GraphMol) -> GraphMol

Return new molecule with stereocenter information set.

See [`setstereocenter!`](@ref))
"""
function setstereocenter(mol::GraphMol)
    newmol = graphmol(mol)
    setstereocenter!(newmol)
    return newmol
end

"""
    optimizewedges!(mol::SDFile) -> Nothing

Optimize dashes and wedges representations. Typical stereocenters can be drawn as four bonds including only a wedge and/or a dash, so if there are too many dashes and wedges, some of them will be converted to normal single bond without changing any stereochemistry.
"""
function optimizewedges!(mol::SDFile)
    # TODO: no longer used
    edges = edgeattrs(mol)
    coords_ = sdfcoords2d(mol)
    for i in 1:nodecount(mol)
        incs = sort(collect(incidences(mol, i)))
        length(incs) in (3, 4) || continue
        upincs = Int[]
        downincs = Int[]
        for inc in incs
            if edges[inc].notation == 1
                push!(mol.edges[inc][1] == i ? upincs : downincs, inc)
            elseif edges[inc].notation == 6
                push!(downincs, inc)
            end
        end
        (isempty(upincs) && isempty(downincs)) && continue  # unspecified
        newbonds = Tuple{Int,Int}[]
        if length(incs) == 4
            if length(upincs) == 3
                downb = (pop!(setdiff(incs, upincs)), 6)
                others = [(inc, 0) for inc in upincs]
                push!(newbonds, downb)
                append!(newbonds, others)
            elseif length(downincs) == 3
                upb = (pop!(setdiff(incs, downincs)), 1)
                others = [(inc, 0) for inc in downincs]
                push!(newbonds, upb)
                append!(newbonds, others)
            elseif length(upincs) == 2 && length(downincs) == 1
                append!(newbonds, [(inc, 0) for inc in upincs])
            elseif length(downincs) == 2 && length(upincs) == 1
                append!(newbonds, [(inc, 0) for inc in downincs])
            else
                nbrs = neighbors(mol, i)
                badjs = [nbrs[b] for b in setdiff(incs, upincs, downincs)]
                upadjs = [nbrs[b] for b in upincs]
                downadjs = [nbrs[b] for b in downincs]
                if length(upincs) == 2 && isempty(downincs)
                    ordered = anglesort(
                        coords_, i, upadjs[1], union(badjs, [upadjs[2]]))
                    if findfirst(isequal(upadjs[2]), ordered) == 2
                        push!(newbonds,(upincs[2], 0))
                    end
                elseif length(upincs) == 2 && length(downincs) == 2
                    ordered = anglesort(
                        coords_, i, upadjs[1], union(downadjs, [upadjs[2]]))
                    if findfirst(isequal(upadjs[2]), ordered) == 2
                        push!(newbonds, (upincs[2], 0))
                    end
                elseif isempty(upincs) && length(downincs) == 2
                    ordered = anglesort(
                        coords_, i, downadjs[1], union(badjs, [downadjs[2]]))
                    if findfirst(isequal(downadjs[2]), ordered) == 2
                        push!(newbonds, (downincs[2], 0))
                    end
                elseif length(upincs) == 1 && length(downincs) == 1
                    ordered = anglesort(
                        coords_, i, badjs[1],
                        [badjs[2], upadjs[1], downadjs[1]]
                    )
                    u = Point2D(coords_, badjs[1])
                    v = Point2D(coords_, badjs[2])
                    d = cross2d(u, v)
                    if findfirst(isequal(badjs[2]), ordered) == 2
                        # if d == 0, stereochem is ambigious
                        if d > 0
                            push!(newbonds, (ordered[1], 0))
                        elseif d < 0
                            push!(newbonds, (ordered[3], 0))
                        end
                    end
                end
            end
        elseif length(incs) == 3
            if length(upincs) == 3
                append!(newbonds, [(inc, 0) for inc in upincs[1:2]])
            elseif length(downincs) == 3
                append!(newbonds, [(inc, 0) for inc in downincs[1:2]])
            elseif length(upincs) == 2
                downb = (pop!(setdiff(incs, upincs)), 6)
                others = [(inc, 0) for inc in upincs]
                push!(newbonds, downb)
                append!(newbonds, others)
            elseif length(downincs) == 2
                upb = (pop!(setdiff(incs, downincs)), 1)
                others = [(inc, 0) for inc in downincs]
                push!(newbonds, upb)
                append!(newbonds, others)
            end
        end
        for (inc, notation) in newbonds
            if mol.edges[inc][2] == i
                mol.edges[inc] = (i, mol.edges[inc][1])  # Reverse edge
            end
            b = setnotation(edges[inc], notation)
            setedgeattr!(mol, inc, b)
        end
    end
    clearcache!(mol)
end

"""
    optimizewedges(mol::GraphMol) -> GraphMol

Return new molecule with wedge configurations optimized.

See [`optimizewedges!`](@ref))
"""
function optimizewedges(mol::GraphMol)
    newmol = graphmol(mol)
    optimizewedges!(newmol)
    return newmol
end


text:
using CuArrays
using CUDAnative
using CUDAdrv

"""
    empiricalbandwidth(nbytes=2*1024^3; devicenumber=0, numtests=10)

Compute the emperical bandwidth in GB/s of a CUDA device using `nbytes` of
memory.  The device to test can be slected with `devicenumber` and the
bandwidth is an average of `ntests`.
"""
function empiricalbandwidth(nbytes=1024^2; devicenumber=0, ntests=10, pinned=false)
    device!(devicenumber)
    stm = CuStream(CUDAdrv.STREAM_NON_BLOCKING)

    nchars = nbytes÷sizeof(Char)

    # Can't use the following because it poisons the CuArrays memory pool.
    #
    # a = CuArray{Char,1}(Mem.alloc(Mem.Host, nbytes), (nchars,))

    a = pinned ? Mem.alloc(Mem.Host, nbytes) : pointer(rand(Char, nbytes))
    b = CuArray{Char}(undef, nchars)

    Mem.copy!(a, b.buf, nbytes, async=true, stream=stm)
    Mem.copy!(b.buf, a, nbytes, async=true, stream=stm)
    synchronize(stm)

    t_dtoh = CUDAdrv.@elapsed stm for n = 1:ntests
        Mem.copy!(a, b.buf, nbytes, async=true, stream=stm)
    end
    bandwidth_dtoh = nbytes*ntests/(t_dtoh*1e9)

    t_htod = CUDAdrv.@elapsed stm for n = 1:ntests
        Mem.copy!(b.buf, a, nbytes, async=true, stream=stm)
    end
    bandwidth_htod = nbytes*ntests/(t_htod*1e9)

    pinned && Mem.free(a)

    (bandwidth_dtoh, bandwidth_htod)
end


text:
using FvCFD

### Freestream Conditions ###
P = 7669
T = 180
U = [ 0, 0, 539 ]
Cp = 1005
R = 287.05
gamma = 1.4
fluid = Fluid(Cp, R, gamma)

### Choose a mesh ###
meshPath = "test/OFdeltaMesh"

boundaryConditions = [
    supersonicInletBoundary, [P, T, U...],
    zeroGradientBoundary, [],
    zeroGradientBoundary, [],
    symmetryBoundary, [],
    zeroGradientBoundary, [],
    wallBoundary, []
]

### Solve ###
mesh = OpenFOAMMesh(meshPath)
cellPrimitives = initializeUniformSolution3D(mesh, P, T, U...)
solve(mesh, meshPath, cellPrimitives, boundaryConditions, LTSEuler, initDt=0.00000001, endTime=3000, outputInterval=10, targetCFL=0.3, silent=false, restart=false, fluid=fluid)

max_stars_repo_name: henrystoldt/fvCFD


text:
using Backtester
using StandardMarketData
using HDF5Utils
using Statistics
using Dates
using Random
using Test

import Backtester: simulate, transition
mutable struct SignalSimulator{S}
    sgnl::S
end
simulate(sim::SignalSimulator, data) = (sim.sgnl, sim.sgnl)
transition(sim::SignalSimulator) = (s, a) -> a

cd(mktempdir())

F, N, T = 2, 5, 100
Random.seed!(1234)

特征名 = idxmap(string.("f", 1:F))
特征 = randn(Float32, F, N, T)
涨幅 = reshape(mean(特征, dims = 1) / 100, N, T)
涨幅 = circshift(涨幅, (0, 1))
ti, Δt = DateTime(2019, 1, 1), Hour(1)
时间戳 = reshape(range(ti, step = Δt, length = T), 1, :)
时间戳 = datetime2unix.(repeat(时间戳, N, 1))
代码 = MLString{8}[string(n) for n in 1:N, t in 1:T]
买1价 = 卖1价 = 最新价 = 1 .+ cumsum(涨幅, dims = 2)
手续费率 = fill(0f0, N, T)
涨停 = 跌停 = zeros(Float32, N, T)
交易池 = ones(Float32, N, T)
data = Data(特征名, 特征, 涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 涨停, 跌停, 交易池)
sgnl = reshape(sum(特征, dims = 1), N, T)
sgnl = @. ifelse(sgnl > 0, 1f0, -1f0)
pnl′ = sum(mean(abs, 涨幅[:, 2:end] , dims = 1)) / ndays(data) * 240
strat = Strategy(sim = SignalSimulator(sgnl))
pnl = backtest(strat, data, mode = "train")
@test pnl ≈ pnl′

时间戳 = range(ti, step = Δt, length = T ÷ 2)
时间戳 = datetime2unix.(repeat(reshape(时间戳, 1, :), N, 2))
代码 = MLString{8}[string(2t <= T ? n : N + n) for n in 1:N, t in 1:T]
时间戳[:, end - 10:end] .= 0
代码[:, end - 10:end] .= MLString{8}("")
买1价 = 最新价 .- 0.001
卖1价 = 最新价 .+ 0.001
data.手续费率 = fill(1f-4, N, T)
data = Data(特征名, 特征, 涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 涨停, 跌停, 交易池)

sgnl = @. ifelse(rand() > 0.05, sgnl, -sgnl)
strat = Strategy(sim = SignalSimulator(sgnl))
pnl = backtest(strat, data, mode = "train")
pnl′ = sum(mean(abs, 涨幅[:, 2:end] , dims = 1)) / ndays(data) * 240
@test pnl > pnl′ * 0.6
backtest(strat, data, mode = "test")


text:
module Zalgo

export zalgo, upsidedown, fraktur, script

const updc = vcat('\u0300':'\u0315', '\u031A', '\u033D',
    '\u034A':'\u034C', '\u0350':'\u0352',
    '\u0357':'\u0358', '\u035B', '\u035D', '\u035E',
    '\u0360':'\u0361', '\u0363':'\u036F',
    '\uFE20':'\uFE22')

const middledc = vcat('\u031B', '\u0334':'\u0338')

const downdc = vcat('\u0316':'\u0319', '\u031C':'\u0333',
    '\u0339':'\u033C', '\u0347':'\u0349',
    '\u034D':'\u034E', '\u0353':'\u0356', '\u0359',
    '\u035A', '\u035C', '\u035F','\u0362')

const upsidedowndict = Dict{String, String}("a" => "ɐ",
    "b" => "q", "c" => "ɔ", "d" => "p", "e" => "ǝ", "f" => "ɟ", "g" => "ƃ",
    "h" => "ɥ", "i" => "ı", "j" => "ɾ", "k" => "ʞ", "l" => "ן", "m" => "ɯ",
    "n" => "u", "o" => "o", "p" => "d", "q" => "b", "r" => "ɹ", "s" => "s",
    "t" => "ʇ", "u" => "n", "v" => "ʌ", "w" => "ʍ", "x" => "x", "y" => "ʎ",
    "z" => "z", "A" => "∀", "B" => "𐐒", "C" => "Ɔ", "D" => "◖", "E" => "Ǝ",
    "F" => "Ⅎ", "G" => "⅁", "H" => "H", "I" => "I", "J" => "ſ", "K" => "⋊",
    "L" => "˥", "M" => "W", "N" => "N", "O" => "O", "P" => "Ԁ", "Q" => "Ό",
    "R" => "ᴚ", "S" => "S", "T" => "⊥", "U" => "∩", "V" => "Λ", "W" => "M",
    "X" => "X", "Y" => "⅄", "Z" => "Z", "0" => "0", "1" => "Ɩ", "2" => "ᄅ",
    "3" => "Ɛ", "4" => "ㄣ", "5" => "ϛ", "6" => "9", "7" => "ㄥ", "8" => "8",
    "9" => "6")

# build the fraktur dict
frakturdict = Dict{String, Char}()
[frakturdict[string(Char(i + 64))] = vcat('\U1D56C':'\U1D585')[i] for i = 1:26]
[frakturdict[string(Char(i + 96))] = vcat('\U1D586':'\U1D59F')[i] for i = 1:26]

# build the Script dict
scriptdict = Dict{String, Char}()
[scriptdict[string(Char(i + 64))] = vcat('\U1D49C':'\U1D4B5')[i] for i = 1:26]
[scriptdict[string(Char(i + 96))] = vcat('\U1D4EA':'\U1D503')[i] for i = 1:26]

function adddc(letter, dc)
    return string(letter) * string(dc[rand(1:end)])
end

"""
    zalgo(str::String, upmarks = 1:4, middlemarks = 1:4,
        downmarks = 1:4, maxmarks = 6)

Randomly add up to `maxmarks` diacritic marks to each letter of `str`. The `upmarks`,
`middlemarks`, and `downmarks` ranges determine the minimum and maximum number of
diacritic marks added to the letter at that position.
"""
function zalgo(text::String;
        upmarks = 1:4,
        middlemarks = 1:4,
        downmarks = 1:4,
        maxmarks = 6)
    letters = split(text, "")
    zalgostring = String[]
    for letter in letters
        # can't add a diacritic mark to some letters
        if !all(isletter, letter)
            push!(zalgostring, letter)
            continue
        end
        upmarks_added     = rand(upmarks.start:upmarks.stop)
        downmarks_added   = rand(downmarks.start:downmarks.stop)
        middlemarks_added = rand(middlemarks.start:middlemarks.stop)
        newletter = letter
        for i in 1:maxmarks
            randint = rand(1:2)
            if randint == 1
                if upmarks_added > 0
                    newletter = adddc(newletter, updc)
                    upmarks_added -= 1
                end
            elseif randint == 2
                if downmarks_added > 0
                    newletter = adddc(newletter, downdc)
                    downmarks_added -= 1
                end
            else
                if middlemarks_added > 0
                    newletter = adddc(newletter, middledc)
                    middlemarks_added -= 1
                end
            end
        end
        push!(zalgostring, newletter)
    end
    return join(zalgostring)
end

"""
    upsidedown(str)

Return a version of string `str` with upside down letters from the Unicode table.
"""
function upsidedown(str)
    asciistr = filter!(c -> haskey(upsidedowndict, c), split(str, ""))
    return join(map(c -> upsidedowndict[c], asciistr))
end

"""
    fraktur(str)

Return a version of string `str` with Fraktur letters from the Unicode table.
"""
function fraktur(str)
    asciistr = filter!(c -> haskey(frakturdict, c), split(str, ""))
    return join(map(c -> frakturdict[c], asciistr))
end

"""
    upsidedown(str)

Return a version of string `str` with script letters from the Unicode table.
"""
function script(str)
    asciistr = filter!(c -> haskey(scriptdict, c), split(str, ""))
    return join(map(c -> scriptdict[c], asciistr))
end

end # module


text:
export dtwplot

"""
    dtwplot(seq1, seq2, [dist=SqEuclidean()]; transportcost=1, diagonal=false)
    dtwplot(seq1, seq2, D, i1, i2; transportcost=1, diagonal=false)

Given two sequences, perform dynamic time warping and plot
the results. If alignment has already been computed, pass
the indices `i1` and `i2` to make the plot.

`diagonal = true` plots a diagonal marker as visual aid.
"""
dtwplot

function handleargs(seq1, seq2, dist::SemiMetric = SqEuclidean(); kwargs...)
    D = dtw_cost_matrix(seq1, seq2, dist; kwargs...)
    cost, i1, i2 = DynamicAxisWarping.trackback(D)
    seq1, seq2, D, i1, i2
end

function handleargs(seq1, seq2, dist, i2min, i2max; kwargs...)
    D = dtw_cost_matrix(seq1, seq2, dist, i2min, i2max; kwargs...)
    cost, i1, i2 = DynamicAxisWarping.trackback(D)
    seq1, seq2, D, i1, i2
end

handleargs(h; kwargs...) = handleargs(h.args...; kwargs...)

@userplot DTWPlot

@recipe function f(h::DTWPlot; transportcost=1, diagonal=false, postprocess=nothing)
    seq1, seq2, D, i1, i2 = handleargs(h; transportcost=transportcost, postprocess=postprocess)

    n1, n2 = lastlength(seq1), lastlength(seq2)

    all = ndims(seq1) ∈ (1,2)
    # set up the subplots
    legend --> false
    link := :both
    grid --> false
    if all
        layout --> @layout [
            a b{0.8w,0.8h}
            _ c
        ]
    else
        layout --> 1
    end

    left_margin --> 0mm
    bottom_margin --> 0mm
    top_margin --> 0mm
    right_margin --> 0mm

    # heatmap
    @series begin
        clims --> (0, 3 * D[end, end])
        seriestype := :heatmap
        formatter --> (z) -> ""
        subplot := (all ? 2 : 1)
        D
    end

    # the rest of the plots are paths

    # main plot
    s1 = @series begin
        seriestype := :path
        linecolor --> :auto
        linewidth --> 3
        subplot := (all ? 2 : 1)
        formatter --> (z) -> ""
        i1, i2
    end

    if diagonal
        m2 = max(n1, n2)
        m1 = min(n1, n2)
        d = m2-m1
        seriestype := :path
        subplot := (all ? 2 : 1)
        imi, ima = radiuslimits(d,seq1, seq2)
        if d == 0
            @series 1:n1
        else
            @series begin
                [imi ima]
            end
        end
    end

    if all
        if ndims(seq1) == 1
            # left line plot
            @series begin
                subplot := 1
                seq2, 1:n2
            end

            # bottom line plot
            @series begin
                subplot := 3
                1:n1, seq1
            end
        else
            # left line plot
            @series begin
                seriestype := :heatmap
                subplot := 1
                seq2'
            end

            # bottom line plot
            @series begin
                seriestype := :heatmap
                subplot := 3
                seq1
            end

        end
    end
end

@userplot MatchPlot

znorm(x) = (x = x.- mean(x); x ./= std(x))
using Statistics
@recipe function f(h::MatchPlot; transportcost=1, separation=2, ds=1, postprocess=nothing)
    x, y, D, i1, i2 = handleargs(h; transportcost=transportcost, postprocess=postprocess)
    x,y = znorm.((x,y))
    s1 = x .- separation
    s2 = y .+ separation

    @series begin
        s1
    end
    @series begin
        s2
    end
    @series begin
        primary := false
        linecolor --> :black
        seriesalpha --> 0.2
        i = fill(Inf, 1, length(i1))
        vec([i1'; i2'; i][:,1:ds:end]), vec([s1[i1]'; s2[i2]'; i][:,1:ds:end])
    end
end

@recipe function plot(r::DTWSearchResult)
    title --> "DTW-NN Search result"
    yguide --> "Distance"
    label --> round(r.cost, sigdigits=4)
    if length(r.dists) == 1
        @series begin
            seriestype := :scatter
            makersize --> 15
            markershape --> :x
            group := 1
            r.dists
        end
        @series begin
            seriestype := :hline
            linestyle := :dash
            primary := false
            group := 1
            r.dists
        end
    else
        @series r.dists
    end

    @series begin
        seriestype := :vline
        linestyle := :dash
        linecolor := :black
        primary := false
        [r.loc]
    end

end


text:
export LinearMapDom, LinearMap

# External packages.
import .LinearMaps
using .LinearMaps: LinearMap

const LMs = LinearMaps

@auto_hash_equals struct LinearMapDom
  N::Int
end

@instance LinearFunctions{LinearMapDom, LinearMap} begin
  @import adjoint, +

  dom(f::LinearMap) = LinearMapDom(size(f,2))
  codom(f::LinearMap) = LinearMapDom(size(f,1))

  compose(f::LinearMap, g::LinearMap) = g*f
  id(V::LinearMapDom) = LMs.UniformScalingMap(1, V.N)

  oplus(V::LinearMapDom, W::LinearMapDom) = LinearMapDom(V.N + W.N)
  oplus(f::LinearMap, g::LinearMap) = LMs.BlockDiagonalMap(f, g)
  mzero(::Type{LinearMapDom}) = LinearMapDom(0)
  swap(V::LinearMapDom, W::LinearMapDom) =
    LinearMap(swap_lm(V.N), swap_lm(W.N), W.N+V.N, V.N+W.N)

  mcopy(V::LinearMapDom) = LinearMap(mcopy_lm, plus_lm, 2*V.N, V.N)
  delete(V::LinearMapDom) = LinearMap(delete_lm, zero_lm(V.N), 0, V.N)
  plus(V::LinearMapDom) = LinearMap(plus_lm, mcopy_lm, V.N, 2*V.N)
  zero(V::LinearMapDom) = LinearMap(zero_lm(V.N), delete_lm, V.N, 0)

  plus(f::LinearMap, g::LinearMap) = f+g
  scalar(V::LinearMapDom, c::Number) = LMs.UniformScalingMap(c, V.N)
  antipode(V::LinearMapDom) = LMs.UniformScalingMap(-1, V.N)

  pair(f::LinearMap, g::LinearMap) = mcopy(dom(f)) ⋅ (f ⊕ g)
  copair(f::LinearMap, g::LinearMap) = (f ⊕ g) ⋅ plus(codom(f))
  proj1(A::LinearMapDom, B::LinearMapDom) = id(A) ⊕ delete(B)
  proj2(A::LinearMapDom, B::LinearMapDom) = delete(A) ⊕ id(B)
  coproj1(A::LinearMapDom, B::LinearMapDom) = id(A) ⊕ zero(B)
  coproj2(A::LinearMapDom, B::LinearMapDom) = zero(A) ⊕ id(B)
end

swap_lm(n::Int) = x::AbstractVector -> vcat(x[n+1:end], x[1:n])
mcopy_lm(x::AbstractVector) = vcat(x, x)
delete_lm(x::AbstractVector) = eltype(x)[]
plus_lm(x::AbstractVector) = begin
  n = length(x) ÷ 2
  x[1:n] + x[n+1:end]
end
zero_lm(n::Int) = x::AbstractVector -> zeros(eltype(x), n)


text:
#
# Make linked list structure for p and data
#

function makell!(first_atom, next_atom, p, data)
    for iat = 1:data.N
        for ix = 1:2
            # p[iat][ix] = p[iat][ix] > 0 ? p[iat][ix] % data.bdim : p[iat][ix] % data.bdim + data.bdim
            if p[iat][ix] > 0
                p[iat][ix] = p[iat][ix] % data.bdim
            else
                if p[iat][ix] % data.bdim == 0
                    p[iat][ix] = p[iat][ix] % data.bdim
                else
                    p[iat][ix] = p[iat][ix] % data.bdim + data.bdim
                end
            end
        end
        icell = trunc(Int64, p[iat][1]/data.cutoff) + 1
        jcell = trunc(Int64, p[iat][2]/data.cutoff) + 1
        next_atom[iat] = first_atom[icell,jcell]
        first_atom[icell,jcell] = iat
    end
end


text:
module TestProcesses
include("preamble.jl")

@testset "mapfoldl" begin
    # https://clojure.org/reference/transducers#_transduce
    xf = Filter(isodd) |> Map(inc)
    @testset for xs in iterator_variants(0:4)
        @test mapfoldl(xf, +, xs, init=0) == 6
        @test mapfoldl(xf, +, xs, init=100) == 106
    end

    # https://clojuredocs.org/clojure.core/transduce
    xf = Filter(isodd) |> Take(10)
    @testset "$(typeof(xs))" for xs in iterator_variants(0:1000)
        @test mapfoldl(xf, push!, xs, init=Int[]) == 1:2:19
        @test mapfoldl(xf, +, xs, init=0) == 100
        @test mapfoldl(xf, +, xs, init=17) == 117
        @test mapfoldl(xf, string, xs, init="") == "135791113151719"

        @test transduce(xf, push!, Int[], xs) == Reduced(1:2:19)
        @test transduce(xf, +, 0, xs) == Reduced(100)
        @test transduce(xf, +, 17, xs) == Reduced(117)
        @test transduce(xf, string, "", xs) == Reduced("135791113151719")
    end

    @testset "empty" begin
        xf = Filter(_ -> false)
        iter = 1:4
        @test_throws EmptyResultError foldl(+, xf, iter)
        @test foldl(+, xf, iter, init=32.) === 32.

        ed = eduction(xf, iter)
        @test_throws EmptyResultError foldl(+, ed)
        @test foldl(+, ed, init=32.) === 32.

        # @test_throws EmptyResultError mapfoldl(xf, +, iter)  # broken
        @test mapfoldl(xf, +, iter, init=32.) === 32.

        nested_xf = Drop(10^9) |> FlagFirst() |> Map(x -> 32)
        @test_throws EmptyResultError foldl(+, nested_xf, iter)
    end

    @testset "zip-of-arrays" begin
        @testset for arrays in [
                (0:3,),
                (0:3, 1:2:8),
                (0:3, 1:2:8, 10:3:19),
                (reshape([1:9;], (3, 3)), [10:3:34;]),
                ]
            xs = zip(arrays...)
            VERSION >= v"1.1-" && @test xs isa Iterators.Zip
            @test mapfoldl(MapSplat(*), +, xs) == sum(map(*, arrays...))
        end
    end

    @testset "product-of-iterators" begin
        iterator_prototypes = [
            (1, 2),
            0:3,
            [-1, 2, 1],
        ]
        @testset for iterators in Iterators.product(iterator_prototypes,
                                                    iterator_prototypes,
                                                    iterator_prototypes),
                n in 1:length(iterators)
            xs = Iterators.product(iterators[1:n]...)
            @test collect(Map(identity), xs) == collect(xs)[:]
            @test mapfoldl(MapSplat(*), +, xs) == sum(Base.splat(*), xs)
        end
    end

    @testset "broadcast" begin
        @testset for xs in iterator_variants(1:3)
            ys = @~ xs.^2
            @test collect(Map(identity), ys) == copy(ys)
            @test mapfoldl(Filter(isodd), +, ys) == 10
            @test mapfoldl(Filter(isodd), +, ys; init=0) == 10
        end
    end
end

@testset "eduction" begin
    # https://clojuredocs.org/clojure.core/eduction
    xf = Filter(isodd) |> Take(5)
    @testset "$(typeof(xs))" for xs in iterator_variants(0:1000)
        @test collect(xf, xs) == 1:2:9
        @test collect(eduction(xf, xs)) == 1:2:9
        @test collect(Map(x -> 2x), eduction(xf, xs)) == 2(1:2:9)
    end

    @testset for xs in iterator_variants(1:5)
        @test collect(xf, xs) == 1:2:5
        @test collect(eduction(xf, xs)) == 1:2:5
        @test collect(Map(x -> 2x), eduction(xf, xs)) == 2(1:2:5)
        @test collect(eduction(AbortIf(iseven) |> TakeLast(1), xs)) == 1:1
    end

    ed = eduction(xf, 1:5)

    @testset "inference" begin
        xf = Zip(Count(), Map(identity), Map(x -> 2x)) |> MapSplat(*)
        ed = eduction(xf, 1:10)
        @test_broken (@inferred foldl(+, ed)) === 6050
    end

    @testset "foreach" begin
        ed = eduction(Map(identity), 1:10)
        result = []
        @test foreach(ed) do x
            x > 5 && return reduced()
            push!(result, x)
        end === reduced()
    end
end

@testset "setinput" begin
    xf = Zip(Count(), Map(identity), Map(x -> 2x)) |> MapSplat(*)
    ed = eduction(xf, 1:10)

    @testset for xs in iterator_variants(1:10)
        @test setinput(ed, xs).coll isa typeof(xs)
        @test collect(setinput(ed, xs)) == collect(ed)
    end

    @testset "inference" begin
        @test (@inferred setinput(ed, [0])).coll isa Vector{Int}
        @test (@inferred setinput(ed, view([0], 1:1))).coll isa SubArray{Int}
        @test_broken (@inferred setinput(ed, Float64[])).coll isa Vector{Float64}
    end
end

# Custom container without unary push! method
struct MinimalContainer
    internal
end

Base.push!(coll::MinimalContainer, x) =
    MinimalContainer(push!(coll.internal, x))

@testset "append!" begin
    # https://clojuredocs.org/clojure.core/into#example-57294b20e4b050526f331420
    xf = Map(x -> x + 2) |> Filter(isodd)
    @testset for xs in iterator_variants(0:9)
        @test append!(xf, [-1, -2], xs) == [-1, -2, 3, 5, 7, 9, 11]
        @test append!(xf, MinimalContainer([-1, -2]), xs).internal ==
            [-1, -2, 3, 5, 7, 9, 11]
    end
end

@testset "map!" begin
    xs0 = 1:10
    @testset for src in [xs0, collect(xs0)]
        dest = similar(xs0)
        @test map!(Scan(+), dest, src) === dest
        @test dest == cumsum(xs0)
    end

    xs0 = 1:5
    @testset for src in [xs0, collect(xs0)]
        dest = zero(src)
        @test map!(Filter(isodd) |> Scan(+), dest, src) == [1, 0, 4, 0, 9]
        @test map!(ScanEmit(tuple, 0), dest, src) == 0:4
    end

    @testset "Expansive transducers are not allowed" begin
        @test_throws Exception map!(Cat(), [0], 1:1)
        @test_throws Exception map!(ScanEmit(tuple, 0, identity), [0], 1:1)
    end
end

@testset "copy!" begin
    xs0 = 1:10
    @testset for src in [xs0, collect(xs0)]
    # @testset for src in iterator_variants(xs0)  # TODO: fix
        dest = similar(xs0)
        @test copy!(Scan(+), dest, src) === dest
        @test dest == cumsum(src)
    end

    xs0 = 1:5
    @testset for src in [xs0, collect(xs0)]
    # @testset for src in iterator_variants(xs0)  # TODO: fix
        dest = zero(xs0)
        @test copy!(Filter(isodd) |> Scan(+), dest, src) == [1, 4, 9]
    end
end

@testset "simple_transduce" begin
    xf = PartitionBy(x -> x > 0) |>
        Filter(xs -> mean(abs, xs) < 1.0) |>
        Map(prod)

    for _ in 1:100
        xs = randn(100)
        @test transduce(xf, +, 0.0, xs) == simple_transduce(xf, +, 0.0, xs)
    end
end

@testset "Transducer(::Reduction)" begin
    # testing Transducer(::Eduction) which calls Transducer(::Reduction)
    @testset for xf in [
            Map(sin),
            TeeZip(Filter(isfinite) |> Map(tan)),
            Map(sin) |> TeeZip(Filter(isfinite) |> Map(tan)) |>
                Map(cos),
            TeeZip(Map(sin) |> TeeZip(Map(tan))),
            TeeZip(TeeZip(Map(tan)) |> Map(identity)),
            TeeZip(Map(sin) |> TeeZip(Map(tan)) |> Map(identity)),
            Map(cos) |> TeeZip(Map(sin) |> TeeZip(Map(tan)) |> Map(identity)),
            Map(cos) |>
                TeeZip(Map(sin) |> TeeZip(Map(tan)) |> Map(identity)) |>
                Map(first),
            TeeZip(TeeZip(Map(tan))),
            ]
        @test Transducer(eduction(xf, 1:1)) === xf
    end
end

@testset "isexpansive" begin
    expansives = [
        Cat()
        Cat() |> Filter(isodd)
        Filter(isodd) |> Cat()
    ]
    @testset for xf in expansives
        @test isexpansive(xf)
        @test isexpansive(TeeZip(xf))
    end
    nonexpansives = [
        Map(identity)
        Map(identity) |> Count()
    ]
    @testset for xf in nonexpansives
        @test !isexpansive(xf)
        @test !isexpansive(TeeZip(xf))
    end
end

@testset "Non-executable transducers" begin
    @testset for ex in @expressions begin
            # mapfoldl(Map(error), +, Int[])  # broken
            foldl(+, Map(error), Int[])
            foldl(+, eduction(Map(error), Int[]))
        end

        err = @eval @test_error $ex
        @test occursin("EmptyResultError:", sprint(showerror, err))
    end

    @testset for ex in @expressions begin
            mapfoldl(Map(error), +, ["hello"]; init=0)
            foldl(+, Map(error), ["hello"]; init=0)
            foldl(+, eduction(Map(error), ["hello"]); init=0)
            foreach(tuple, Map(error), ["hello"])
            foreach(tuple, eduction(Map(error), ["hello"]))
            collect(Map(error), ["hello"])
            collect(eduction(Map(error), ["hello"]))
            append!(Map(error), [], ["hello"])
            copy!(Map(error), [], ["hello"])
        end

        @eval @test_throws ErrorException("hello") $ex
    end
end

@testset "IdentityNotDefinedError" begin
    @test_throws(
        IdentityNotDefinedError,
        foldl((x, y) -> x + y, Map(identity), 1:1; init=Init),
    )
end

@testset "identityof error" begin
    @test_throws EmptyResultError mapfoldl(Map(identity), right, Any[])
    err = @test_error mapfoldl(Map(identity), +, Any[])
    @test_broken err isa EmptyResultError

    err = @test_error mapfoldl(Map(identity), +, Vector{Int}[])
    @test_broken err isa EmptyResultError
end

@testset "AdHocFoldable" begin
    counting = AdHocFoldable() do rf, acc, _
        i = 0
        while true
            i += 1
            acc = @next(rf, acc, i)
        end
    end

    @test collect(Take(3), counting) == 1:3

    @test foreach(counting) do i
        i == 3 && return reduced()
    end === reduced()
end

@testset "reduce" begin
    # Test that InitialValues.jl can handle right identity
    @test reduce(+, Filter(x -> x == 1), 1:2; basesize=1) == 1
end

@testset "copy" begin
    @test copy(Map(x -> x^2), Vector, 1:3) == (1:3) .^ 2
    @test copy(Map(x -> x^2), [1:3;]) == (1:3) .^ 2
    @test copy(Map(x -> x^2), Set, 1:3) == Set((1:3) .^ 2)
    @test copy(Map(x -> x => x^2), Dict, 2:2) == Dict(2 => 4)
end

end  # module


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

main_paths = ["."]
for main_path in main_paths
    format(
        main_path;
        whitespace_ops_in_indices = true,
        remove_extra_newlines = true,
        verbose = true,
        always_for_in = true,
        whitespace_typedefs = true,
        # always_use_return = true. # Disabled since it throws a lot of false positives
    )
end


text:
using ModelingToolkit, Test

@testset "Parsing Test" begin include("variable_parsing.jl") end
@testset "Basic Variables and Operations" begin include("basic_variables_and_operations.jl") end
@testset "Differentiation Test" begin include("derivatives.jl") end
@testset "Internal Test" begin include("internal.jl") end
@testset "Simplify Test" begin include("simplify.jl") end
@testset "Ambiguity Test" begin include("ambiguity.jl") end
@testset "Components Test" begin include("components.jl") end
@testset "System Construction Test" begin include("system_construction.jl") end


text:
export
	benchmark_solve!,
	shift_fill!

function interp_rows(N::Int,tf::Float64,X::AbstractMatrix)::Matrix
    n,N1 = size(X)
    t1 = range(0,stop=tf,length=N1)
    t2 = collect(range(0,stop=tf,length=N))
    X2 = zeros(n,N)
    for i = 1:n
        interp_cubic = CubicSplineInterpolation(t1, X[i,:])
        X2[i,:] = interp_cubic(t2)
    end
    return X2
end

function convertInf!(A::VecOrMat{Float64},infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

function benchmark_solve!(solver; samples=10, evals=10)
    Z0 = deepcopy(get_trajectory(solver))
    v0 = solver.opts.verbose
    s0 = solver.opts.show_summary
    solver.opts.verbose = 0
    solver.opts.show_summary = false
    b = @benchmark begin
        TO.initial_trajectory!($solver,$Z0)
        solve!($solver)
    end samples=samples evals=evals
    solver.opts.verbose = v0
    solver.opts.show_summary = s0
    return b
end

function benchmark_solve!(solver, data::Dict; samples=10, evals=10)
	b = benchmark_solve!(solver, samples=samples, evals=evals)

   # Run stats
   push!(data[:time], time(median(b))*1e-6)  # ms
   push!(data[:iterations], iterations(solver))
   push!(data[:cost], cost(solver))
   return b
end

function shift_fill!(A::Vector, n=1, mode=:copy)
    N = length(A)
    N <= 1 && return nothing
	@inbounds for k = n+1:N
		A[k-n] = A[k]
    end
    if mode == :copy
        a_last = A[N-n]
        @inbounds for k = N-n:N
            A[k] = copy(a_last)
        end
    else mode == :zero
        @inbounds for k = N-n:N
            A[k] .= zero(A[k])
        end
    end
	return nothing
end

"""
    triukkt

Form an upper-triangular KKT matrix with regularization on the dual variables and
an active set `a`. Takes a matrix `A`, a sparse matrix of size `(Np, Np + Nd)` where `Np`
is the number of primal variables and `Nd` is the number of dual variables. The vector
`a` is a vector of Booleans of size `(Np + Nd,)` denoting which columns to remove from `A`.

The CSC data `colptr`, `rowval`, and `nzval` need to be large enough to store the
resulting information, which should always be the case if `length(colptr) >= n + 1`,
`length(rowval) >= nnz(A) + Nd` and `length(nzval) >= nnz(A) + Nd`. These vectors
are updated in-place.

The function returns `nnz_new`, the new number of non-zero entries. A new sparse matrix can
be created from the updated data:

    n_new = sum(a)
    SparseMatrixCSC(n_new, n_new, colptr[1:n_new+1], rowval[1:nnz_new], nzval[1:nnz_new])

To do this without allocations, you can use the `resize!` command.
"""
function triukkt!(A::SparseMatrixCSC{Tv,Ti}, a::AbstractVector{Bool},
                    colptr, rowval, nzval; reg=zero(Tv)) where {Tv,Ti}
    n,m = size(A)
    m_new = sum(a)
    nnz_new = 0
    for j = 1:m
        nvals = A.colptr[j+1] - A.colptr[j]
        nnz_new += nvals * a[j]
    end

    colptr[1] = 1
    col = 1
    nzcount = 0
    for j = 1:size(A,2)
        nvals = A.colptr[j+1] - A.colptr[j]  # number of elements in the column

        # Keep column if a[j] is true
        if a[j]
            # Assign next column pointer
            colptr[col+1] = colptr[col] + nvals

            # Loop over elements in the column, copying into new data
            for i = 1:nvals
                rowval[nzcount + i] = A.rowval[A.colptr[j] + i - 1]
                nzval[nzcount +  i] = A.nzval[A.colptr[j] + i - 1]
            end
            if col > n
                rowval[nzcount + nvals + 1] = col
                nzval[nzcount + nvals + 1] = -reg
                colptr[col+1] += 1
                nzcount += 1
            end
            nzcount += nvals
            col += 1
        end
    end
    m2 = col - 1 - n
    nnz_new += m2
    return nnz_new
end


text:
# # Heat PDE
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Heat3D.ipynb)
#
#md # !!! note "Overview"
#md #     System type: Affine system\
#md #     State dimension: 48\
#md #     Application domain: Mechanical Engineering
#
# ## Model description
#
# The system is described by the linear differential equations:
#
# ```math
#   \begin{array}{lcl}
#   \dot{x}(t) &=& Ax(t) + Bu(t),\qquad u(t) \in \mathcal{U} \\
#   y(t) &=& C x(t)
#   \end{array}
# ```
#
# There are two versions of this benchmark:
#
# - *(time-varying inputs):* The inputs can change arbitrarily over time: $\forall t: u(t)\in \mathcal{U}$.
#
# - *(constant inputs):* The inputs are uncertain only in their initial value, and
#    constant over time: ``u(0)\in \mathcal{U}``, ``\dot u (t)= 0``.

using ReachabilityAnalysis, SparseArrays, JLD2

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
HEAT01_path = joinpath(examples_dir, "Heat3D", "HEAT01.jld2")


text:
module OptionalArgChecks

using IRTools: @dynamo, IR, recurse!, block, branches, branch!
using MacroTools: postwalk

export @argcheck, @skipargcheck

"""
    @argcheck ex

Marks `ex` as an optional argument check, so when a function is called via
[`@skipargcheck`](@ref), `ex` will be omitted.

```jldoctest
julia> function half(x::Integer)
           @argcheck iseven(x) || throw(DomainError(x, "x has to be an even number"))
           return x ÷ 2
       end
half (generic function with 1 method)

julia> half(4)
2

julia> half(3)
ERROR: DomainError with 3:
x has to be an even number
[...]

julia> @skipargcheck half(3)
1
```
"""
macro argcheck(ex)
    return Expr(:block, Expr(:meta, :begin_argcheck), esc(ex), Expr(:meta, :end_argcheck))
end

@dynamo function skipargcheck(x...)
    ir = IR(x...)
    ir === nothing && return
    next = iterate(ir)
    while next !== nothing
        (x, st), state = next
        if Meta.isexpr(st.expr, :meta) && st.expr.args[1] === :begin_argcheck
            orig = block(ir, x)
            delete!(ir, x)

            (x, st), state = iterate(ir, state)
            while !(Meta.isexpr(st.expr, :meta) && st.expr.args[1] === :end_argcheck)
                delete!(ir, x)
                (x, st), state = iterate(ir, state)
            end

            dest = block(ir, x)
            if orig != dest
                empty!(branches(orig))
                branch!(orig, dest)
            end
            delete!(ir, x)
        end
        next = iterate(ir, state)
    end
    recurse!(ir)
    return ir
end

"""
    @skipargcheck ex

For every function call in `ex`, expressions wrapped in [`@argcheck`](@ref) get omitted
recursively.
"""
macro skipargcheck(ex)
    ex = postwalk(ex) do x
        if Meta.isexpr(x, :call)
            pushfirst!(x.args, GlobalRef(@__MODULE__, :skipargcheck))
        end
        return x
    end
    return esc(ex)
end

end


text:
module DataDrivenLinearModels

# Write your package code here.

end


text:
module TestLearningNetworkMachines

using Test
using ..Models
using ..TestUtilities
using MLJBase
using Tables
using StableRNGs
using Serialization
rng = StableRNG(616161)

# A dummy clustering model:
mutable struct DummyClusterer <: Unsupervised
    n::Int
end
DummyClusterer(; n=3) = DummyClusterer(n)
function MLJBase.fit(model::DummyClusterer, verbosity::Int, X)
    Xmatrix = Tables.matrix(X)
    n = min(size(Xmatrix, 2), model.n)
    centres = Xmatrix[1:n, :]
    levels = categorical(1:n)
    report = (centres=centres,)
    fitresult = levels
    return fitresult, nothing, report
end
MLJBase.transform(model::DummyClusterer, fitresult, Xnew) =
    selectcols(Xnew, 1:length(fitresult))
MLJBase.predict(model::DummyClusterer, fitresult, Xnew) =
    [fill(fitresult[1], nrows(Xnew))...]

N = 20
X = (a = rand(N), b = categorical(rand("FM", N)))

@testset "signature helpers" begin
    @test MLJBase._call(NamedTuple()) == NamedTuple()
    a = source(:a)
    b = source(:b)
    W = source(:W)
    yhat = source(:yhat)
    s = (transform=W,
         report=(a=a, b=b),
         predict=yhat)
    @test MLJBase._report_part(s) == (a=a, b=b)
    @test MLJBase._operation_part(s) == (transform=W, predict=yhat)
    @test MLJBase._nodes(s) == (W, yhat, a, b)
    @test MLJBase._operations(s) == (:transform, :predict)
    R = MLJBase._call(MLJBase._report_part(s))
    @test R.a == :a
    @test R.b == :b
end

@testset "wrapping a learning network in a machine" begin

    # unsupervised:
    Xs = source(X)
    W = transform(machine(OneHotEncoder(), Xs), Xs)
    clust = DummyClusterer(n=2)
    m = machine(clust, W)
    yhat = predict(m, W)
    Wout = transform(m, W)
    rnode = source(:stuff)

    # test of `fitted_params(::NamedTuple)':
    fit!(Wout, verbosity=0)

    @test_throws(MLJBase.ERR_BAD_SIGNATURE,
                 machine(Unsupervised(),
                         predict=yhat,
                         fitted_params=rnode))
    @test_throws(MLJBase.ERR_EXPECTED_NODE_IN_SIGNATURE,
                 machine(Unsupervised(),
                         predict=42))
    @test_throws(MLJBase.ERR_EXPECTED_NODE_IN_SIGNATURE,
                 machine(Unsupervised(), Xs;
                         predict=yhat,
                         transform=Wout,
                         report=(some_stuff=42,)))
    mach = machine(Unsupervised(), Xs;
                   predict=yhat,
                   transform=Wout,
                   report=(some_stuff=rnode,))
    @test mach.args == (Xs, )
    @test mach.args[1] == Xs
    fit!(mach, force=true, verbosity=0)
    Θ = mach.fitresult
    @test Θ.predict == yhat
    @test Θ.transform == Wout
    Θ.report.some_stuff == rnode
    @test report(mach).some_stuff == :stuff

    @test report(mach).machines == fitted_params(mach).machines

    # supervised
    y = rand("ab", N) |> categorical;
    ys = source(y)
    mm = machine(ConstantClassifier(), W, ys)
    yhat = predict(mm, W)
    e = @node auc(yhat, ys)

    @test_throws Exception machine(predict=yhat)
    mach = machine(Probabilistic(), Xs, ys;
                   predict=yhat,
                   report=(training_auc=e,))
    @test mach.model isa Probabilistic
    @test_throws ArgumentError machine(Probabilistic(), Xs, ys)
    @test_throws ArgumentError machine(Probabilistic(), Xs, ys;
                                       report=(training_auc=e,))

    # test extra report items coming from `training_auc=e` above
    fit!(mach, verbosity=0)
    err = auc(yhat(), y)
    @test report(mach).training_auc ≈ err

    # supervised - predict_mode
    @test predict_mode(mach, X) == mode.(predict(mach, X))
    predict_mode(mach, rows=1:2) == predict_mode(mach, rows=:)[1:2]

    # evaluate a learning machine
    evaluate!(mach, measure=LogLoss(), verbosity=0)

    # supervised - predict_median, predict_mean
    X1, y1 = make_regression(20)

    Xs = source(X1); ys = source(y1)
    mm = machine(ConstantRegressor(), Xs, ys)
    yhat = predict(mm, Xs)
    mach = fit!(machine(Probabilistic(), Xs, ys; predict=yhat), verbosity=0)
    @test predict_mean(mach, X1) ≈ mean.(predict(mach, X1))
    @test predict_median(mach, X1) ≈ median.(predict(mach, X1))

end

# build a learning network:
x1 = map(n -> mod(n,3), rand(rng, UInt8, 100)) |> categorical;
x2 = randn(rng, 100);
X = (x1=x1, x2=x2);
y = x2.^2;

Xs = source(X)
ys = source(y)
z = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)
enode = @node mae(ys, yhat)

@testset "replace method for learning network machines, acceleration: $(typeof(accel))" for accel in (CPU1(), CPUThreads())

    fit!(yhat, verbosity=0, acceleration=accel)

    # test nested reporting:
    r = MLJBase.report(yhat)
    d = r.report_given_machine
    ms = machines(yhat)
    @test ms == r.machines |> reverse
    @test all(mach -> report(mach) == d[mach], ms)

    hot2 = deepcopy(hot)
    knn2 = deepcopy(knn)

    # duplicate a learning network machine:
    mach  = machine(Deterministic(), Xs, ys;
                    predict=yhat,
                    report=(mae=enode,))
    mach2 = replace(mach, hot=>hot2, knn=>knn2,
                    ys=>source(ys.data);
                    empty_unspecified_sources=true)
    ss = sources(glb(mach2))
    @test isempty(ss[1])
    mach2 = @test_logs((:warn, r"No replacement"),
                       replace(mach, hot=>hot2, knn=>knn2,
                               ys=>source(ys.data)))
    yhat2 = mach2.fitresult.predict
    fit!(mach, verbosity=0)
    fit!(mach2, verbosity=0)
    @test predict(mach, X) ≈ predict(mach2, X)
    @test report(mach).mae ≈ report(mach2).mae

    @test mach2.args[1]() == Xs()
    @test mach2.args[2]() == ys()

    ## EXTRA TESTS FOR TRAINING SEQUENCE

    # pickout the newly created machines:
    standM2 = machines(yhat2, stand) |> first
    oakM2 = machines(yhat2, oak) |> first
    knnM2 = machines(yhat2, knn) |> first
    hotM2 = machines(yhat2, hot) |> first

    @test_mach_sequence(fit!(yhat2, force=true, acceleration=accel),
                        [(:train, standM2), (:train, hotM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:train, hotM2), (:train, standM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:train, standM2), (:train, hotM2),
                         (:train, oakM2), (:train, knnM2)],
                        [(:train, hotM2), (:train, standM2),
                         (:train, oakM2), (:train, knnM2)])

    @test length(MLJBase.machines(yhat)) == length(MLJBase.machines(yhat2))
    @test MLJBase.models(yhat) == MLJBase.models(yhat2)
    @test sources(yhat) == sources(yhat2)
    @test MLJBase.tree(yhat) == MLJBase.tree(yhat2)
    @test yhat() ≈ yhat2()

    # this change should trigger retraining of all machines except the
    # univariate standardizer:
    hot2.drop_last = true
    @test_mach_sequence(fit!(yhat2, acceleration=accel),
                        [(:skip, standM2), (:update, hotM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:update, hotM2), (:skip, standM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:skip, standM2), (:update, hotM2),
                         (:train, oakM2), (:train, knnM2)],
                        [(:update, hotM2), (:skip, standM2),
                         (:train, oakM2), (:train, knnM2)])
end

@testset "Test serializable of pipeline" begin
    filename = "pipe_mach.jls"
    X, y = make_regression(100, 1)
    pipe = (X -> coerce(X, :x₁=>Continuous)) |> DecisionTreeRegressor()
    mach = machine(pipe, X, y)
    fit!(mach, verbosity=0)

    # Check serializable function
    smach = MLJBase.serializable(mach)
    @test_logs((:warn, MLJBase.warn_serializable_mach(predict),
                (:warn, MLJBase.warn_serializable_mach(predict)),
                fit!(mach, verbosity=0)))
    TestUtilities.generic_tests(mach, smach)
    @test keys(fitted_params(smach)) == keys(fitted_params(mach))
    @test keys(report(smach)) == keys(report(mach))
    # Check data has been wiped out from models at the first level of composition
    @test length(machines(glb(smach))) == length(machines(glb(mach)))
    for submach in machines(glb(smach))
        TestUtilities.test_data(submach)
    end

    # End to end
    MLJBase.save(filename, mach)
    smach = machine(filename)
    @test predict(smach, X) == predict(mach, X)

    rm(filename)
end

@testset "Test serializable of composite machines" begin
    filename = "stack_mach.jls"
    X, y = make_regression(100, 1)
    model = Stack(
        metalearner = DecisionTreeRegressor(),
        tree1 = DecisionTreeRegressor(min_samples_split=3),
        tree2 = DecisionTreeRegressor(),
        measures=rmse)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)

    # Check serializable function
    smach = MLJBase.serializable(mach)
    TestUtilities.generic_tests(mach, smach)
    # Check data has been wiped out from models at the first level of composition
    @test length(machines(glb(smach))) == length(machines(glb(mach)))
    for submach in machines(glb(smach))
        TestUtilities.test_data(submach)
    end

    # Testing extra report field : it is a deepcopy
    @test smach.report.cv_report === mach.report.cv_report

    @test smach.fitresult isa MLJBase.CompositeFitresult

    Serialization.serialize(filename, smach)
    smach = Serialization.deserialize(filename)
    MLJBase.restore!(smach)

    @test MLJBase.predict(smach, X) == MLJBase.predict(mach, X)
    @test keys(fitted_params(smach)) == keys(fitted_params(mach))
    @test keys(report(smach)) == keys(report(mach))

    rm(filename)

    # End to end
    MLJBase.save(filename, mach)
    smach = machine(filename)
    @test predict(smach, X) == predict(mach, X)

    rm(filename)
end

@testset "Test serializable of nested composite machines" begin
    filename = "nested_stack_mach.jls"
    X, y = make_regression(100, 1)

    pipe = (X -> coerce(X, :x₁=>Continuous)) |> DecisionTreeRegressor()
    model = Stack(
        metalearner = DecisionTreeRegressor(),
        pipe = pipe)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)

    MLJBase.save(filename, mach)
    smach = machine(filename)

    @test predict(smach, X) == predict(mach, X)

    # Test data as been erased at the first and second level of composition
    for submach in machines(glb(smach))
        TestUtilities.test_data(submach)
        if submach isa Machine{<:Composite}
            for subsubmach in machines(glb(submach))
                TestUtilities.test_data(subsubmach)
            end
        end
    end

    rm(filename)

end

struct DummyRangeCV
    cv
end

torange(x::UnitRange) = x
torange(x) = minimum(x):maximum(x)

function MLJBase.train_test_pairs(dcv::DummyRangeCV, rows, X, y)
    ttp = MLJBase.train_test_pairs(dcv.cv, rows)
    return [(torange(t),torange(e)) for (t,e) in ttp]
end

@testset "Test serialized filesize does not increase with datasize" begin
    # At the moment it is necessary to resort to a custom resampling strategy,
    # for this test. This is because partial functions working on nodes,
    # such as `selectrows`in learning networks store data.
    # A basic CV would store vectors which would grow in size as the dataset grows.

    dcv = DummyRangeCV(CV(nfolds=3))
    model = Stack(
        metalearner = FooBarRegressor(lambda=1.),
        resampling = dcv,
        model_1 = DeterministicConstantRegressor(),
        model_2=ConstantRegressor())

    filesizes = []
    for n in [100, 500, 1000]
        filename = "serialized_temp_$n.jls"
        X, y = make_regression(n, 1)
        mach = machine(model, X, y)
        fit!(mach, verbosity=0)
        MLJBase.save(filename, mach)
        push!(filesizes, filesize(filename))
        rm(filename)
    end
    @test all(x==filesizes[1] for x in filesizes)
    # What if no serializable procedure had happened
    filename = "full_of_data.jls"
    X, y = make_regression(1000, 1)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)
    serialize(filename, mach)
    @test filesize(filename) > filesizes[1]

    @test_logs (:warn, MLJBase.warn_bad_deserialization(mach.state)) machine(filename)

    rm(filename)
end

mutable struct DummyComposite <: DeterministicComposite
    stand1
    stand2
end

@testset "issue 377" begin
    stand = Standardizer()
    model = DummyComposite(stand, stand)

    Xs = source()
    mach1 = machine(model.stand1, Xs)
    X1 = transform(mach1, Xs)
    mach2 = machine(model.stand2, X1)
    X2 = transform(mach2, X1)

    mach = machine(Unsupervised(), Xs; transform=X2)
    @test_logs((:error, r"The hyper"),
               @test_throws(ArgumentError,
                            MLJBase.network_model_names(model, mach)))
end

end

true


text:
using LinearAlgebra

A = [2 3 0;
     1 2 4;
     0 2 1]

println("norm(A): ", norm(A))