text:
using Dates
using Printf
using Libdl
using Sockets
@static if Base.Sys.islinux()
    # RawProcDLLName = joinpath( pwd(),"../Jl-Aux/dll/RawProc.so"  )
	extention =".so"
   # RawDataSaveDLLName= "RawProc"
    # ghDLLRawProc = Libdl.dlopen(RawProcDLLName)
else
	extention =".dll"
    #RawProcDLLName = joinpath( pwd(),"../Jl-Aux/dll/RawProc.dll" )
    # RawDataSaveDLLName= "RawDataSave.dll"
    # ghDLLRawProc = Libdl.dlopen(RawProcDLLName)
end

ghDLLRawProc = nothing
fetchJuliaFilepointer =nothing
function dllLoad()
    println(Base.@__MODULE__,"================" )
    global ghDLLRawProc
    RawProcDLLName = "RawProc"*extention
    @warn "try load dll:",RawProcDLLName ,ENV["LD_LIBRARY_PATH"]

    ghDLLRawProc = Libdl.dlopen(RawProcDLLName)
    @warn RawProcDLLName " dll load ok" ghDLLRawProc
    # Name = ghDLLRawProc
    # baseGlbPath= replace(baseGlbPath,"\\"=>"/")
    eval("Main.ghDLLRawProc=ghDLLRawProc")
    return
end

Base.MainInclude.include("clearFolderlst.jl")

function gettupleFromStr(str::String ,idx )
    len =length( str )
    if idx <= len
        return str[ idx]
    else
        return Char(0x00)
    end
end
# ghDLL = nothing
function Init(needSaveFile,needJuliaSend ,Logfolder,  folderA, folderB,NeedFlowControl )
    println(Base.@__MODULE__ )
    global fetchJuliaFilepointer
    global ghDLLRawProc
    @warn "Init "
    # global ghDLL
    println( "------------------------------1-Init------------------------------")
    # println( "-------------------------------2------------------------------")
    @warn "init 1:::",Logfolder, ",", folderA,",", folderB,needSaveFile,needJuliaSend

    #initSaveData = Libdl.dlsym( Main.ghDLLRawProc,:initSaveData )

    initSaveData = Libdl.dlsym( Main.ghDLLRawProc,:initSaveData )
   @warn    "----- initSaveData:"  initSaveData
    # println( "------------------------------3-------------------------------")
    t = ccall(initSaveData,
    Int32,
    (Cstring ,Cstring, Cstring ),
    Logfolder,  folderA, folderB)
    @warn "init 2:::"

   # initJuliaData = Libdl.dlsym(Main.ghDLLRawProc,:initJuliaData )
    initJuliaData = Libdl.dlsym(Main.ghDLLRawProc,:initJuliaData )
@warn    "----- initJuliaData:"  initJuliaData
    t = ccall( initJuliaData ,
    Cint,
    (Cuchar  ,Cuchar,Cuchar ),
    needSaveFile,  needJuliaSend,NeedFlowControl)
    @warn "init 2:::"
    # initJuliaData( bool needSaveFile, bool needJuliaSend );
    # println( "------------------------------4-------------------------------")

     @warn "..............  ............fetchFiberDatas  Entered.."  ghDLLRawProc   :fetchJuliaFile
    fetchJuliaFilepointer = Libdl.dlsym(ghDLLRawProc,:fetchJuliaFile ) # fetchJuliaFile
	@warn "fetchFiberDatas  leaved.."   fetchJuliaFilepointer

end

function StartupMainThread(  )
    println( "------------------------------1--StartupMainThread-----------------------------")
    global ghDLLRawProc
    startProcessRawData = Libdl.dlsym(ghDLLRawProc,:startProcessRawData )
    t = ccall( startProcessRawData ,
    Cvoid,
    ( ) )
end
function InitIPAddressSvr(IpStr,Port )
    println(Base.@__MODULE__ )
    global ghDLLRawProc
    println( "------------------------------1--InitIPAddressSvr-----------------------------")
    IpStr1 = IpStr
    if IpStr =="0.0.0.0"
        IpStr1= string( Sockets.getipaddr() )
    end
    println( IpStr1 )

    initSvrIP = Libdl.dlsym(ghDLLRawProc,:initSvrIP )
    t = ccall( initSvrIP,
    Cvoid,
    (Cstring ,UInt16  ),
    IpStr1,  Port )
end

function waitAllComplete()
    while true
        sleep(60*100);
        println("Continue...")
    end
end

function ClearExpiredFiles(RAW1MaxSpace, RAW2MaxSpace, RAW3MaxSpace, RAW4MaxSpace,Unit,root1,root2="")
    @warn "ClearExpiredFiles :::"
    # println( "RAW1 MaxSpace:$RAW1MaxSpace Mib,RAW2 MaxSpace:$RAW2MaxSpace Mib,RAW3 MaxSpace:$RAW3MaxSpace Mib,RAW4 MaxSpace:$RAW4MaxSpace Mib" )
    while true
        # println( "start clear expired files,-$root1:$(now() )- ")
        getSizes( root1 ,RAW1MaxSpace, RAW2MaxSpace, RAW3MaxSpace, RAW4MaxSpace,Unit)
        # println( "Complete clear expired files,--$root1:$(now() )-")
        println()
        sleep(10)
        if root2 != ""
            # println( "start clear expired files,-$root2:$(now() )- ")
            getSizes( root2 ,RAW1MaxSpace, RAW2MaxSpace, RAW3MaxSpace, RAW4MaxSpace,Unit)
            # println( "Complete clear expired files,--$root2:$(now() )- ")
        end
        sleep(200)
    end
end

const constbufferbyteSize =512*8192
mutable struct  FiberRawData
    rawDatatype::Ref{Cshort}
    machineName::Array{UInt8,1}
    startDTStr::Array{UInt8,1}
    cfgdata::Array{UInt16,1}
    datas::Array{UInt16,1}
    actualDatabytes::Ref{Cint}
    function FiberRawData()
        rawDatatype = Ref{Cshort}(0)
        machineName=Array{UInt8,1}(undef,20)
        startDTStr=Array{UInt8,1}(undef,14)
        cfgdata=Array{UInt16,1}(undef,7)
        datas=Array{UInt16,1}(undef, constbufferbyteSize )
        actualDatabytes=Ref{Cint}(0)
        new(rawDatatype ,machineName , startDTStr, cfgdata, datas, actualDatabytes   )
    end
end

const constbufferbyteSize =512*8192
mutable struct  fileFiberRawData
    rawDatatype::UInt16
    machineName::Array{UInt8,1}
    startDTStr::Array{UInt8,1}
    cfgdata::Array{UInt16,1}
    datas::Array{UInt16,1}
    actualData::UInt
    function fileFiberRawData(actualData)
        rawDatatype = 0
        machineName=Array{UInt8,1}(undef,20)
        startDTStr=Array{UInt8,1}(undef,14)
        cfgdata=Array{UInt16,1}(undef,7)
        datas=Array{UInt16,1}(undef, actualData )
        # actualDatabytes=Ref{Cint}(0)

        new(rawDatatype ,machineName , startDTStr, cfgdata, datas , actualData   )
    end
end

function ShowContent( fdo::FiberRawData)
    println("-------------------------------")
    println( fdo.rawDatatype[] )
    println( String(fdo.startDTStr))

    mn = String(fdo.machineName )
    j=1
    for i = 1:20
        if !isprint( mn[i] )
            break
        end
        j = i
    end
    # println( j)
    a =mn[1:j]
    println( ( a   ) )

    # println( ( fdo.machineName   ) )
    println( fdo.actualDatabytes[]  )
end
# Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8} ,In
function fetchFiberDatas(fdobj)
    # type1=Ref{Cshort}(0)
    # println("---fetchFiberDatas")
    # @warn "fetchFiberDatas  Entered.."  ghDLLRawProc   Main.fetchJuliaFilepointer
    #fetchJuliaFile1 = Libdl.dlsym(ghDLLRawProc,:fetchJuliaFile ) # fetchJuliaFile
	# @warn "fetchFiberDatas  leaved.."   Main.fetchJuliaFilepointer  fetchJuliaFilepointer
	# println("234324324")
    t = ccall( fetchJuliaFilepointer,
        Cint,
        (Ref{Cshort},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        fdobj.rawDatatype,
        fdobj.machineName ,
        fdobj.startDTStr ,
        pointer(reinterpret(UInt8,fdobj.cfgdata)),
        pointer(reinterpret(UInt8,fdobj.datas) ) ,fdobj.actualDatabytes
    );
    #  @warn "fetchFiberDatas  Exit .."
    #println( t )
    if t == 0
        println("--------Exception:Cannot get Socket Data in expected Time ")
        exit(-1)
    else

    end
    # println( fdobj.rawDatatype,"=================================")
    # ShowContent( fdobj)
end

function getCfg( cfgJl::String )
    evs =""
    cfgObj =open( cfgJl,"r")
    skip( cfgObj,3)
    for line in Base.readlines(cfgObj)
        if  strip( line ) != "" && line[1:1] != "#"
            evs = evs *line *"; "
        end
    end

    reps = Base.Meta.parse( evs )
    # println( reps )
    Core.eval(Main,reps )
end

function getStringFromArray( machineName::Array{UInt8,1})
    a=0
    for a1 in machineName
        if a1 == 0x00
            break
        else
            a +=1
			if a>=8
			break
		end
        end
    end
    if a >0
        return String(machineName[ 1:a ])
    else
        return "NO-Name"
    # println("here ",String(machineName[ 1:a ]) )
    end
end

const baseDatetime ="20170211"*"000000"
glbjzdt = Dates.DateTime(2017,02,11)

function DT2UInt( DT::String,format1="yyyymmddHHMMSS"  )
    # return UInt(100)
    # df = Dates.DateFormat(format1);

    # curdt = Dates.DateTime(DT,format1)
    year = parse( UInt,DT[1:4])
    month= parse( UInt,DT[5:6])
    day= parse( UInt,DT[7:8])
    hour= parse( UInt,DT[9:10])
    minute=  parse( UInt,DT[11:12])
    secs= parse( UInt,DT[13:14])
    curdt = Dates.DateTime(year,month,day,hour,minute,secs)
    a =  curdt- glbjzdt

    return  floor(UInt,a.value/1000)
    # return  floor(UInt,curdt.value/1000)
end

function UInt2DTStr( DT::UInt,format1="yyyymmddHHMMSS" )
    # try
        # df = Dates.DateFormat("yyyymmddHHMMSS");
        # jzdt = Dates.DateTime(baseDatetime  ,df )
        # a =  UInt(1)::UInt64
        # a= b1*1 +b2*256 +b3*256*256 +b4*256*256*256
        jzdt1 = glbjzdt + Dates.Second( DT )
        # print( jzdt1," -----" )
        dt =@sprintf( "%04d.%02d.%02d %02d.%02d.%02d",Dates.year(jzdt1 ) ,   Dates.month(jzdt1 ), Dates.day(jzdt1 ), Dates.hour(jzdt1 ), Dates.minute(jzdt1 ), Dates.second(jzdt1 ) )
        return dt
    # catch(e )
    #     println(e)
    #     return "2000.10.10.00.10.00"
    # end
end

function showCurDT2File(lastWriteDT::UInt, curDt::UInt)::UInt
    # return 0::UInt
    if curDt - lastWriteDT < 120
        return  UInt(0);
    end

    info1 = "CurDatatime=" *string(curDt )*"\r\n"
    fileLog_wind.LogEvent(info1);

    a1 = Dates.now()
    # interval = (a1 - lasttime )
    # if interval.value > 1*60*1000
    #    lasttime = a1
    # @warn  "Computer time is :" ,a1
    @warn  "CurDatatime=" *string(curDt )* ","*UInt2DTStr( curDt )   *"" ," ;Computer time is :" ,a1
    return curDt
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
@testset "preserve" begin
    @testset "basic" begin
        @system SPreserve(Controller) begin
            a => 1 ~ track
            b(a) => a + 1 ~ accumulate
            c(b) => b ~ preserve
        end
        s = instance(SPreserve)
        @test s.a' == 1 && s.b' == 0 && s.c' == 0
        update!(s)
        @test s.a' == 1 && s.b' == 2 && s.c' == 0
    end

    @testset "optional" begin
        @system SPreserveOptional(Controller) begin
            a ~ preserve(optional)
            b => 1 ~ preserve(optional)
        end
        s = instance(SPreserveOptional)
        @test isnothing(s.a') && s.b' == 1
        update!(s)
        @test isnothing(s.a') && s.b' == 1
    end

    @testset "parameter" begin
        @system SParameter(Controller) begin
            a => 1 ~ preserve(parameter)
        end
        s = instance(SParameter)
        @test s.a' == 1
        update!(s)
        @test s.a' == 1
    end

    @testset "parameter with config" begin
        @system SParameterConfig(Controller) begin
            a => 1 ~ preserve(parameter)
        end
        o = SParameterConfig => :a => 2
        s = instance(SParameterConfig; config=o)
        @test s.a' == 2
    end

    @testset "parameter with config alias" begin
        @system SParameterConfigAlias(Controller) begin
            a: aa => 1 ~ preserve(parameter)
            bb: b => 1 ~ preserve(parameter)
        end
        o = SParameterConfigAlias => (:a => 2, :b => 2)
        s = instance(SParameterConfigAlias; config=o)
        @test s.a' == 2
        @test s.b' == 2
    end

    @testset "parameter missing" begin
        @system SParameterMissing(Controller) begin
            a => 1 ~ preserve(parameter)
        end
        o = SParameterMissing => :a => missing
        s = instance(SParameterMissing; config=o)
        @test s.a' == 1
    end

    @testset "parameter nothing" begin
        @system SParameterNothing(Controller) begin
            a => 1 ~ preserve(parameter, optional)
        end
        o = SParameterNothing => :a => nothing
        s = instance(SParameterNothing; config=o)
        @test s.a' === nothing
    end

    @testset "minmax" begin
        @system SPreserveMinMax(Controller) begin
            a => 0 ~ preserve(parameter, min=1)
            b => 0 ~ preserve(parameter, max=2)
            c => 0 ~ preserve(parameter, min=1, max=2)
        end
        s = instance(SPreserveMinMax)
        @test s.a' == 1
        @test s.b' == 0
        @test s.c' == 1
    end

    @testset "parameter minmax" begin
        @system SParameterMinMax(Controller) begin
            a => 0 ~ preserve(parameter, min=-1)
            b => 0 ~ preserve(parameter, max=1)
            c => 0 ~ preserve(parameter, min=-1, max=1)
        end
        o1 = SParameterMinMax => (:a => 2, :b => 2, :c => 2)
        s1 = instance(SParameterMinMax; config=o1)
        @test s1.a' == 2
        @test s1.b' == 1
        @test s1.c' == 1
        o2 = SParameterMinMax => (:a => -2, :b => -2, :c => -2)
        s2 = instance(SParameterMinMax; config=o2)
        @test s2.a' == -1
        @test s2.b' == -2
        @test s2.c' == -1
    end

    @testset "round" begin
        @system SPreserveRound(Controller) begin
            a => 1.4 ~ preserve(round)
            b => 1.4 ~ preserve(round=:round)
            c => 1.4 ~ preserve(round=:ceil)
            d => 1.4 ~ preserve(round=:floor)
            e => 1.4 ~ preserve(round=:trunc)
        end
        s = instance(SPreserveRound)
        @test s.a' === s.b'
        @test s.b' === 1.0
        @test s.c' === 2.0
        @test s.d' === 1.0
        @test s.e' === 1.0
    end

    @testset "parameter round" begin
        @system SParameterRound(Controller) begin
            a ~ preserve::int(parameter, round=:floor)
            b ~ preserve::int(parameter, round=:trunc)
        end
        o1 = :0 => (a = 1.4, b = 1.4)
        s1 = instance(SParameterRound; config=o1)
        @test s1.a' === 1
        @test s1.b' === 1
        o2 = :0 => (a = -1.4, b = -1.4)
        s2 = instance(SParameterRound; config=o2)
        @test s2.a' === -2
        @test s2.b' === -1
    end
end

ext: jl


text:
include("GeneticAlgorithm.jl")
using PyPlot

##############################################################################
# Plotting

function plotResults(results)
    figure(figsize=(8,6))
    subplot(221)
    p1 = plot(results["entropy"])
    title("Population Entropy")
    subplot(222)
    p2 = plot(results["numIndividuals"])
    title("Population Size")
    subplot(223)
    p3 = plot(results["cost"])
    title("Cost")
    subplot(224)
    genTime = [results["generationTime"][i]-results["generationTime"][i-1]
        for i=2:length(results["generationTime"])]
    p4 = plot(genTime)
    title("Execution Time")
end

function plotSurvivors(population, survivors)
    indexScore, diversityScore, fitnessScore = fitnessDiversityScore(population,hamming)
    fig = figure()
    title("Survivors Plot")
    xlabel("Fitness")
    ylabel("Diversity")
    G = scatter(fitnessScore[survivors],diversityScore[survivors],color="green", label="Survivors", s=30)
    R = scatter(fitnessScore[.~survivors],diversityScore[.~survivors],color="red", label = "Non-survivors", s=30)
    legend(loc="right");
end

function plotSurvivorsUnnorm(fitness, divContrib, survivors)
    fig = figure()
    title("Survivors Plot - Unnormalized")
    xlabel("Fitness")
    ylabel("Diversity")
    G = scatter(fitness[survivors],divContrib[survivors],color="green", label="Survivors", s=30)
    R = scatter(fitness[.~survivors],divContrib[.~survivors],color="red", label = "Non-survivors", s=30)
    #plt.ylim(69.0, 68.4)
    plt.ylim(maximum(divContrib),minimum(divContrib))
    legend(loc="right")
end

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
import Dolo

model = Dolo.yaml_import("./model_tax.yaml")

model.calibration

m0, s0, x0, p = model.calibration[:exogenous, :states, :controls, :parameters]

Dolo.transition(model, m0, s0, x0, m0, p)

model = Dolo.yaml_import("examples/models/rbc.yaml")

model.calibration

m0, s0, x0, p = model.calibration[:exogenous, :states, :controls, :parameters]

Dolo.transition(model, m0, s0, x0, m0, p)


text:
# Returns a partitioning of `range` into sub-ranges of length `partition_size`.
# If `range` is not divisible by `partition_size`, the last sub-range will be
# of length `range` mod `partition_size`.
function partition_range(range::UnitRange, partition_size::Int)
    n_partitions = ceil(Int, last(range) / partition_size)
    partitions = Vector{UnitRange{Int}}(undef, n_partitions)

    for p in 1:n_partitions
        l::Int = first(range) + (p-1) * partition_size
        r::Int = min(l - 1 + partition_size, last(range))
        partitions[p] = l:r
    end

    return partitions
end

function calculate_sample_weights(targets::T) where {T <: AbstractVector}
    t = targets
    N = length(t)
    ψ₀ = sum(t .!= 1) / N
    ψ₁ = sum(t .== 1) / N
    ψ = [q == 1 ? ψ₁ : ψ₀ for q in t]
end

function gaussian_projection_matrix(::Type{T},
                                    n_neurons::Int,
                                    n_features::Int) where {T <: Number}
    randn(T, n_neurons, n_features) ./ sqrt(n_features)
end

function project(samples::T₁,
                 weights::T₂,
                 activation_functions::Vector{T₃}) where {T₁ <: AbstractMatrix,
                                                          T₂ <: AbstractMatrix,
                                                          T₃ <: ActivationFunction}
    X = samples
    W = weights
    fs = activation_functions
    L = length(fs)

    H = W * X
    for l in 1:L
        H[l,:] .= fs[l].(H[l,:])
    end

    H
end


text:
"Return some planning options for a dist network that also provide flexibility to transmission."
function solve_td_decoupling_distribution(mn_data::Dict{String,Any}; optimizer, setting=Dict{String,Any}(), number_of_candidates)
    _FP.require_dim(mn_data, :sub_nw)
    if _FP.dim_length(mn_data, :sub_nw) > 1
        Memento.error(_LOGGER, "A single distribution network is required ($(dim_length(mn_data, :sub_nw)) found)")
    end
    if number_of_candidates < 1
        Memento.warn(_LOGGER, "The number of requested distribution network planning candidates must be positive: substituting $number_of_candidates with 1.")
        number_of_candidates = 1
    end

    candidates = Dict{String,Any}()
    add_dist_candidate!(candidates, min_investments(mn_data, optimizer, setting)...)
    if number_of_candidates >= 2
        add_dist_candidate!(candidates, max_power(mn_data, optimizer, setting)...)
        add_dist_candidates_cost!(candidates, mn_data)
        if number_of_candidates > 2
            intermediate_costs = collect(range(candidates["min"]["cost"], candidates["max"]["cost"], length = number_of_candidates)[2:end-1])
            for i in 1:number_of_candidates-2
                add_dist_candidate!(candidates, intermediate_investment(mn_data, i, intermediate_costs[i], optimizer, setting)...)
            end
        end
    end
    add_dist_candidates_cost!(candidates, mn_data)
    return candidates
end

"Run a model with usual parameters and model type; error if not solved to optimality."
function run_td_decoupling_model(data::Dict{String,Any}, build_function::Function, optimizer; kwargs...)
    Memento.info(_LOGGER, "running $(String(nameof(build_function)))...")
    result = _PM.run_model(
        data, _FP.BFARadPowerModel, optimizer, build_function;
        ref_extensions = [_FP.ref_add_gen!, _FP.ref_add_storage!, _FP.ref_add_ne_storage!, _FP.ref_add_flex_load!, _PM.ref_add_on_off_va_bounds!, _FP.ref_add_ne_branch_allbranches!, _FP.ref_add_frb_branch!, _FP.ref_add_oltc_branch!],
        solution_processors = [_PM.sol_data_model!, _FP.sol_td_coupling!],
        multinetwork = true,
        kwargs...
    )
    Memento.info(_LOGGER, "solved $(String(nameof(build_function))) in $(round(Int,result["solve_time"])) seconds")
    if result["termination_status"] ∉ (_PM.OPTIMAL, _PM.LOCALLY_SOLVED)
        Memento.error(_LOGGER, "Unable to solve $(String(nameof(build_function))) ($(result["optimizer"]) termination status: $(result["termination_status"]))")
    end
    return result
end

## Distribution candidates

function min_investments(mn_data, optimizer, setting)
    min_base = run_td_decoupling_model(mn_data, _FP.post_flex_tnep, optimizer; setting)

    mn_data_current_investments = deepcopy(mn_data)
    add_ne_branch_indicator!(mn_data_current_investments, min_base)
    add_ne_storage_indicator!(mn_data_current_investments, min_base)

    min_import = run_td_decoupling_model(mn_data_current_investments, build_max_import_with_current_investments, optimizer; setting)
    min_export = run_td_decoupling_model(mn_data_current_investments, build_max_export_with_current_investments, optimizer; setting)

    return ("min", min_import, min_export, min_base)
end

function max_power(mn_data, optimizer, setting)
    max_import_1 = run_td_decoupling_model(mn_data, build_max_import, optimizer; setting)
    mn_data_max_import_2 = deepcopy(mn_data)
    add_td_coupling_power_active!(mn_data_max_import_2, max_import_1)

    max_export_1 = run_td_decoupling_model(mn_data, build_max_export, optimizer; setting)
    mn_data_max_export_2 = deepcopy(mn_data)
    add_td_coupling_power_active!(mn_data_max_export_2, max_export_1)
    _FP.shift_nws!(mn_data_max_export_2)

    mn_data_max_pair = _FP.merge_multinetworks!(mn_data_max_import_2, mn_data_max_export_2, :sub_nw)
    max_pair = run_td_decoupling_model(mn_data_max_pair, build_min_cost_with_same_investments_in_all_sub_nws, optimizer; setting)
    max_import, max_export = split_result(max_pair)

    mn_data_max_base = deepcopy(mn_data)
    add_ne_branch_indicator!(mn_data_max_base, max_import) # Same investments as max_export
    add_ne_storage_indicator!(mn_data_max_base, max_import) # Same investments as max_export
    max_base = run_td_decoupling_model(mn_data_max_base, build_min_cost_with_fixed_investments, optimizer; setting)

    return ("max", max_import, max_export, max_base)
end

function intermediate_investment(mn_data, i, intermediate_cost, optimizer, setting; kwargs...)
    mn_data_1 = deepcopy(mn_data)
    mn_data_2 = deepcopy(mn_data)
    _FP.shift_nws!(mn_data_2)
    mn_data_pair = _FP.merge_multinetworks!(mn_data_1, mn_data_2, :sub_nw)
    add_max_cost!(mn_data_pair, intermediate_cost)
    pair_1 = run_td_decoupling_model(mn_data_pair, build_max_flex_band_with_bounded_cost, optimizer; setting, kwargs...)

    add_td_coupling_power_active!(mn_data_pair, pair_1)
    pair_2 = run_td_decoupling_model(mn_data_pair, build_min_cost_with_same_investments_in_all_sub_nws, optimizer; setting, kwargs...)
    i_import, i_export = split_result(pair_2)

    mn_data_base = deepcopy(mn_data)
    add_ne_branch_indicator!(mn_data_base, i_import) # Same investments as i_export
    add_ne_storage_indicator!(mn_data_base, i_import) # Same investments as i_export
    i_base = run_td_decoupling_model(mn_data_base, build_min_cost_with_fixed_investments, optimizer; setting, kwargs...)

    return ("intermediate_$i", i_import, i_export, i_base)
end

## Results

"""
Generate two results from one result of a run command, by assigning half of nws to each result

Nw ids of the original result are ordered and splitted; the first half is used in both the generated
results. Other fields are copied.
"""
function split_result(res::Dict{String,Any})
    res_1 = Dict{String,Any}()
    res_2 = Dict{String,Any}()
    for (k_res,v_res) in res
        if k_res == "solution"
            res_1["solution"] = Dict{String,Any}()
            res_2["solution"] = Dict{String,Any}()
            for (k_sol, v_sol) in v_res
                if k_sol == "nw"
                    nw_ids = string.(sort(parse.(Int,keys(v_sol))))
                    if isodd(length(nw_ids))
                        Memento.error(_LOGGER, "Attempting to split a result having an odd number of nws.")
                    else
                        nn = length(nw_ids)
                        nw_ids_1 = nw_ids[1:nn÷2]
                        nw_ids_2 = nw_ids[nn÷2+1:nn]
                        res_1["solution"]["nw"] = Dict{String,Any}()
                        res_2["solution"]["nw"] = Dict{String,Any}()
                        for pos in 1:nn÷2
                            res_1["solution"]["nw"][nw_ids_1[pos]] = v_sol[nw_ids_1[pos]]
                            res_2["solution"]["nw"][nw_ids_1[pos]] = v_sol[nw_ids_2[pos]] # In res_2 the same nw_ids of res_1 are used
                        end
                    end
                else
                    res_1["solution"][k_sol] = v_sol
                    res_2["solution"][k_sol] = deepcopy(v_sol)
                end
            end
        else
            res_1[k_res] = v_res
            res_2[k_res] = deepcopy(v_res)
        end
    end
    return res_1, res_2
end

## Distribution candidates handling

function add_dist_candidate!(dist_candidates::Dict{String,Any}, name::String, r_import::Dict{String,Any}, r_export::Dict{String,Any}, r_base::Dict{String,Any} = Dict{String,Any}())

    c = dist_candidates[name] = Dict{String,Any}()

    result = c["result"] = Dict{String,Any}()
    result["export"] = r_export
    result["import"] = r_import
    if !isempty(r_base)
        result["base"] = r_base
    end

    # Results to compare with: r_export is used at first, then r_import and possibly r_base are compared with r_export.
    ctrl_res = isempty(r_base) ? (r_import,) : (r_import, r_base)

    # Store ordered ids of components (nw, branch, etc.) and check that they are the same in each result.
    ids = c["ids"] = Dict{String,Any}()
    ids["nw"] = keys(r_export["solution"]["nw"])
    for res in ctrl_res
        if ids["nw"] ≠ keys(res["solution"]["nw"])
            Memento.error(_LOGGER, "Results of flex candidate \"$name\" have different nw ids.")
        end
    end
    ids["nw"] = string.(sort(parse.(Int, ids["nw"])))
    first_period = ids["nw"][1]
    for comp in ("branch", "ne_branch", "storage", "ne_storage")
        if haskey(r_export["solution"]["nw"][first_period], comp)
            comp_keys = keys(r_export["solution"]["nw"][first_period][comp])
            for res in ctrl_res
                if comp_keys ≠ keys(res["solution"]["nw"][first_period][comp])
                    Memento.error(_LOGGER, "Results of flex candidate \"$name\" have different $comp ids.")
                end
            end
            ids[comp] = string.(sort(parse.(Int, comp_keys)))
        else
            ids[comp] = Vector{String}()
        end
    end

    # Check that investment decisions are the same in each result (only first period is used);
    # create a component key only if component is present in results.
    investment = c["investment"] = Dict{String,Any}()
    function _isbuilt(result, component, built_keyword)
        Bool.(round.(result["solution"]["nw"][first_period][component][k][built_keyword] for k in ids[component]))
    end
    for (comp, built_keyword) in ("ne_branch"=>"built", "ne_storage"=>"isbuilt")
        if !isempty(ids[comp])
            built = _isbuilt(r_export, comp, built_keyword)
            for res in ctrl_res
                if built ≠ _isbuilt(res, comp, built_keyword)
                    Memento.error(_LOGGER, "Results of flex candidate \"$name\" have different $comp investment decisions.")
                end
            end
            investment[comp] = Dict{String,Any}(ids[comp] .=> built)
        end
    end
end

function add_dist_candidates_cost!(dist_candidates::Dict{String,Any}, mn_data::Dict{String,Any})

    # Store investment cost of each network component regardless of whether it is built or not (it depends on dist candidate)
    cost_lookup = Dict{Int,Any}()
    for n in _FP.nw_ids(mn_data; hour=1, scenario=1)
        sn_data = mn_data["nw"]["$n"]
        cl_n = cost_lookup[n] = Dict{String,Any}()
        if haskey(sn_data, "ne_branch")
            cl_n_comp = cl_n["ne_branch"] = Dict{String,Float64}()
            for (i,comp) in sn_data["ne_branch"]
                cl_n_comp[i] = comp["construction_cost"]
            end
        end
        if haskey(sn_data, "ne_storage")
            cl_n_comp = cl_n["ne_storage"] = Dict{String,Float64}()
            for (i,comp) in sn_data["ne_storage"]
                cl_n_comp[i] = comp["eq_cost"] + comp["inst_cost"]
            end
        end
    end

    # Compute the cost of each dist candidate
    for candidate in values(dist_candidates)
        candidate["cost"] = sum( sum( sum( cost_lookup[n][comp_name][id] * built for (id, built) in comp_inv) for (comp_name, comp_inv) in candidate["investment"]) for n in _FP.nw_ids(mn_data; hour=1, scenario=1))
    end
end


text:
# Wrapper for some BLIS utility routines.
#

macro blis_group_lvutil_normv(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          n,        BliDim,
                          x,        Ptr{xType}, incx, BliInc,
                          norm,     Ptr{rType})
    end
end

@blis_group_lvutil_normv norm1v
@blis_group_lvutil_normv normfv
@blis_group_lvutil_normv normiv

macro blis_group_lvutil_normm(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          diaga,    BliDoff,
                          uploa,    BliUpLo,
                          m,        BliDim,
                          n,        BliDim,
                          a,        Ptr{xType}, rsa,  BliInc, csa, BliInc,
                          norm,     Ptr{rType})
    end
end

@blis_group_lvutil_normm norm1m
@blis_group_lvutil_normm normfm
@blis_group_lvutil_normm normim

macro blis_group_lvutil_mkspecial(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          uploa,    BliUpLo,
                          m,        BliDim,
                          a,        Ptr{xType}, rsa,  BliInc, csa, BliInc)
    end
end

@blis_group_lvutil_mkspecial mkherm
@blis_group_lvutil_mkspecial mksymm
@blis_group_lvutil_mkspecial mktrim

macro blis_group_lvutil_randv(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          n,        BliDim,
                          x,        Ptr{xType}, incx, BliInc)
    end
end

@blis_group_lvutil_randv randv
@blis_group_lvutil_randv randnv

macro blis_group_lvutil_randm(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          uploa,    BliUpLo,
                          m,        BliDim,
                          n,        BliDim,
                          a,        Ptr{xType}, rsa,  BliInc, csa, BliInc)
    end
end

@blis_group_lvutil_randm randm
@blis_group_lvutil_randm randnm

@blis_ccall_group(asumv,
                  Cvoid,
                  n,        BliDim,
                  x,        Ptr{xType}, incx, BliInc,
                  asum,     Ptr{rType})
@blis_ccall_group(printv,
                  Cvoid,
                  s1,       Ptr{Cchar},
                  m,        BliDim,
                  x,        Ptr{xType}, incx, BliInc,
                  format,   Ptr{Cchar},
                  s2,       Ptr{Cchar})
@blis_ccall_group(printm,
                  Cvoid,
                  s1,       Ptr{Cchar},
                  m,        BliDim,
                  n,        BliDim,
                  a,        Ptr{xType}, rsa,  BliInc, csa, BliInc,
                  format,   Ptr{Cchar},
                  s2,       Ptr{Cchar})
@blis_ccall_group(sumsqv,
                  Cvoid,
                  n,        BliDim,
                  x,        Ptr{xType}, incx, BliInc,
                  scale,    Ptr{rType},
                  sumsq,    Ptr{rType})


text:
import MultipleShootingPEM
ms = MultipleShootingPEM
using Base.Test
using Calculus

include("test_util.jl")
include("test_sensitivity_equations.jl")
include("test_space_state_model.jl")
include("test_one_shoot_simulation.jl")
include("test_multiple_shooting.jl")
include("test_optimization_problem.jl")


text:
@testset "Rule Matching Tests" begin

    @test ALifeBenchmark._is_match("", "") == true
    @test ALifeBenchmark._is_match("", "1") == true
    @test ALifeBenchmark._is_match("1", "1") == true
    @test ALifeBenchmark._is_match("1", "100") == true
    @test ALifeBenchmark._is_match("1", "0") == false

    @test ALifeBenchmark._is_match("101", "0") == false
    @test ALifeBenchmark._is_match("101", "1") == false
    @test ALifeBenchmark._is_match("101", "101") == true
    @test ALifeBenchmark._is_match("101", "10111") == true

end

@testset "Rule Decoding Tests" begin

    @test ALifeBenchmark._get_rule_at("11011011011111111", 1) == ALifeBenchmark.Rule("1", "1", "1", true, true, true, true, true, true, 1, 17)
    @test ALifeBenchmark._get_rule_at("11011011011111111", 2) == nothing
    @test ALifeBenchmark._get_rule_at("11011011011111111", 17) == nothing

    @test ALifeBenchmark.get_rules("11011011011111111") == [
        ALifeBenchmark.Rule("1", "1", "1", true, true, true, true, true, true, 1, 17)
    ]

    @test ALifeBenchmark._get_rule_at("111011001011100000", 1) == ALifeBenchmark.Rule("", "1", "0", false, true, true, true, false, false, 1, 15)
    @test ALifeBenchmark._get_rule_at("111011001011100000", 2) == ALifeBenchmark.Rule("1", "0", "1", true, false, false, false, false, false, 2, 18)
    @test ALifeBenchmark._get_rule_at("111011001011100000", 3) == nothing

    @test ALifeBenchmark.get_rules("111011001011100000") == [
        ALifeBenchmark.Rule("", "1", "0", false, true, true, true, false, false, 1, 15),
        ALifeBenchmark.Rule("1", "0", "1", true, false, false, false, false, false, 2, 18)
    ]

end

@testset "Rule Filtering Tests" begin

    @test ALifeBenchmark.find_best_matches("0", [""]) == [""]
    @test ALifeBenchmark.find_best_matches("0", ["", "0", "1"]) == ["0"]
    @test ALifeBenchmark.find_best_matches("0", ["", "0", "1", "01"]) == ["0"]
    @test ALifeBenchmark.find_best_matches("000000", ["01", "0100", "110001", "0", "000"]) == ["000"]
    @test ALifeBenchmark.find_best_matches("010001", ["01", "0100", "110001"]) == ["0100"]

    @test ALifeBenchmark.filter_rules([
        ALifeBenchmark.Rule("", "11", "111", false, false, false, false, false, false, 0, 0),
        ALifeBenchmark.Rule("0", "11", "111", false, false, false, false, false, false, 1, 1),
        ALifeBenchmark.Rule("01", "11", "111", false, false, false, false, false, false, 2, 2),
        ALifeBenchmark.Rule("00", "11", "111", false, false, false, false, false, false, 3, 3)
    ]) == [
        ALifeBenchmark.Rule("00", "11", "111", false, false, false, false, false, false, 3, 3),
        ALifeBenchmark.Rule("01", "11", "111", false, false, false, false, false, false, 2, 2),
        ALifeBenchmark.Rule("", "11", "111", false, false, false, false, false, false, 0, 0)
    ]

    @test ALifeBenchmark.filter_rules([
        ALifeBenchmark.Rule("", "11", "111", false, false, false, false, false, false, 0, 0),
        ALifeBenchmark.Rule("0", "11", "111", false, false, false, false, false, false, 1, 1),
        ALifeBenchmark.Rule("01", "11", "0011", false, false, false, false, false, false, 2, 2),
        ALifeBenchmark.Rule("00", "11", "111", false, false, false, false, false, false, 3, 3),
        ALifeBenchmark.Rule("000", "11", "111", false, false, false, false, false, false, 4, 4),
        ALifeBenchmark.Rule("11", "11", "111", false, false, false, false, false, false, 5, 5)
        ]) == Any[
            ALifeBenchmark.Rule("00", "11", "111", false, false, false, false, false, false, 3, 3),
            ALifeBenchmark.Rule("000", "11", "111", false, false, false, false, false, false, 4, 4),
            ALifeBenchmark.Rule("01", "11", "0011", false, false, false, false, false, false, 2, 2),
            ALifeBenchmark.Rule("11", "11", "111", false, false, false, false, false, false, 5, 5)
    ]

    @test ALifeBenchmark.filter_rules([
        ALifeBenchmark.Rule("", "11", "111", false, false, false, false, false, false, 0, 0),
        ALifeBenchmark.Rule("0", "11", "111", false, false, false, false, false, false, 1, 1),
        ALifeBenchmark.Rule("01", "11", "0011", false, false, false, false, false, false, 2, 2),
        ALifeBenchmark.Rule("00", "11", "0011", false, false, false, false, false, false, 1, 3),
        ALifeBenchmark.Rule("000", "11", "111", false, false, false, false, false, false, 4, 4),
        ALifeBenchmark.Rule("11", "11", "111", false, false, false, false, false, false, 5, 5)
        ]) == Any[
        ALifeBenchmark.Rule("00", "11", "0011", false, false, false, false, false, false, 1, 3),
        ALifeBenchmark.Rule("000", "11", "111", false, false, false, false, false, false, 4, 4),
        ALifeBenchmark.Rule("11", "11", "111", false, false, false, false, false, false, 5, 5)
    ]

    @test ALifeBenchmark.filter_rules([
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 5),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 1, 6)
        ]) == Any[
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 5),
    ]

    @test ALifeBenchmark.filter_rules([
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 5),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 10, 16)
        ]) == Any[
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 5),
    ]

end

@testset "Overlapping Tests" begin

    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 10)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 10)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 20)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 10)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 7)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 20),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 3, 7)
    ) == true
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 10, 20),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 30)
    ) == true

    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 4),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 5, 10)
    ) == false
    @test ALifeBenchmark._overlap(
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 4),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 20, 40)
    ) == false

    rules = [
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 4),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 20, 40),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 50, 60)
    ]
    to_add = [
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 4, 10),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 61, 70),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 65, 80)
    ]
    ALifeBenchmark._append_non_overlapping!(rules, to_add)
    @test rules == [
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 0, 4),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 20, 40),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 50, 60),
        ALifeBenchmark.Rule("", "", "", false, false, false, false, false, false, 61, 70)
    ]

end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using GLAbstraction, GLPlot, Reactive, Colors
GLPlot.init()

function zdata(i, j, t)
    x         = (i - 0.5)
    z         = (j - 0.5)
    radius    = sqrt((x * x) + (z * z))

    r = sin(10.0f0 * radius + t)
    Float32(r + rand(1.0:0.01:1.1))
end
function zcolor(i, j, t)
    x         = (i - 0.5)
    z         = (j - 0.5)
    radius     = sqrt((x * x) + (z * z))

    r = sin(10.0f0 * radius + t)
    g = cos(10.0f0 * radius + t)
    b = radius
    return RGBA{Float32}(r, g, b, 1)
end
t = GLPlot.play_widget(linspace(1, 50, 200))
N = 128
surf = map(t-> [zdata(i/N, j/N, t) for i = 1:N, j = 1:N], t)
color = map(t-> [zcolor(i/N, j/N, t) for i = 1:N, j = 1:N], t)
xr = linspace(-4, 4, N)
obj = glplot(
    surf, :surface, color = color, ranges = (xr, xr)
)


text:
using TermWin
using Base.Test

# write your own tests here
using Lint

msgs = lintpkg( "TermWin", returnMsgs = true )
println( msgs )
sumseverity = 0
if !isempty( msgs )
    sumseverity = sum( x->x.level, msgs )
end
@test sumseverity == 0

include( "dftests.jl" )
include( "strtests.jl" )


text:
# -*- coding: utf-8 -*-
using LightGraphs
using SemanticModels
using SemanticModels.ModelTools.CategoryTheory
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.GraphModels
using ModelingToolkit
using Petri

# # Example

n = 3
a = smallgraph(:house)
b = StarGraph(n)

f = FinSetMorph([5,3,4])
@assert FinSetMorph([5,4,3,2,1])(model(GraphModel, a)).model |>
    adjacency_matrix |>
    collect == [0 1 1 0 0;
                1 0 1 1 0;
                1 1 0 0 1;
                0 1 0 0 1;
                0 0 1 1 0]

# +
# Create a decorated morphism for a
g = FinSetMorph(1:5, [5,4])
# Add graph a as decoration of morphism f
dec_g = Decorated(g, model(GraphModel, a))
# Create a decorated morphism for b
f = FinSetMorph(1:3, [1,2])
# Add graph b as decoration of morphism g
dec_f = Decorated(f, model(GraphModel, b))
# Create the span of decorated morphisms
s = Span(dec_f,dec_g)
# Solve for the decorated cospan that solves the pushout defined by the span
H = CategoryTheory.pushout(s)

dec = decorations(H, GraphModel)[1]

@show collect(edges(dec.model))

# +
@variables S, I, R, I′

# Create S -> I -> R model
sir = Petri.Model([S, I, R], [(S+I, 2I), (I,R)])
sir_model = model(PetriModel, sir)

# Create Susceptible to 2 diseases model
sii = Petri.Model([S, I, I′], [(S+I,  2I ), (S+I′, 2I′)])
sii_model = model(PetriModel, sii)

# Decorate a finite set with SIR
f = FinSetMorph(1:3, [1, 2])
dec_f = Decorated(f, sir_model)

# Decorate a finite set with SII
g = FinSetMorph(1:3, [1, 2])
dec_g = Decorated(g, sii_model)

# Create a span of decorated morphisms
s = Span(dec_f, dec_g)

# Solve the pushout that combines the two Petri decorations
H = CategoryTheory.pushout(s)

map(petri->println(petri.model.S, "\n", petri.model.Δ), decorations(H, PetriModel))
# -


text:
"Diagonalization for EPEC complex"

include("MILPs/MILP_complex.jl") #The MILP model with complex bidding

function run_diagonaluization(EPEC)

    function compute_profits(model)
        lambda = value.(model[:λ][:])
        #NON-STRATEGIC GENERATORS
        gen_1 = value.(model[:g]["GEN1",:])
        gen_2 = value.(model[:g]["GEN2",:])
        gen_3 = value.(model[:g]["GEN3",:])
        gen_5 = value.(model[:g]["GEN5",:])
        gen_6 = value.(model[:g]["GEN6",:])
        gen_7 = value.(model[:g]["GEN7",:])
        gen_8 = value.(model[:g]["GEN8",:])

        #DEMAND
        d = value.(model[:d])
        #NON-STRATEGIC RES
        res_PV = value.(model[:w][:PV,:])
        #DISPATCHED QUANTITIES FROMS STR AGENTS
        dch = value.(model[:dch_st])
        ch = value.(model[:ch_st])
        gen = value.(model[:g_st])
        res = value.(model[:w_st])
        #CALCULATING PROFITS
        earned_ESS = sum(dch[t] * lambda[t] for t in t_steps)
        paid_ESS = sum(ch[t] * lambda[t] for t in t_steps)
        profit_ESS =  earned_ESS - paid_ESS
        # println("PROFIT_ESS: ", profit_ESS)

        earned_GEN = sum(gen[t] * lambda[t] for t in t_steps)
        paid_GEN = sum(dat_str.GEN_STR.MC_g[t] * gen[t] for t in t_steps)
        profit_GEN =  earned_GEN - paid_GEN
        # println("PROFIT_GEN: ", profit_GEN)

        earned_RES = sum(res[t] * lambda[t] for t in t_steps)
        paid_RES = sum(dat_str.WIND_STR.MC_g[t] * res[t] for t in t_steps)
        profit_RES =  earned_RES - paid_RES
        # println("PROFIT_RES: ", profit_RES)

        return profit_ESS, profit_GEN, profit_RES
    end

    ###########READING THE SOLUTIONS OF THE EPEC#############
    mod_dat_str = deepcopy(dat_str)
    mod_dat_str.GEN_STR.G[:] = value.(EPEC[:G_ST])
    mod_dat_str.WIND_STR.G[:] = value.(EPEC[:W_ST])
    mod_dat_str.ESS_STR.DCH[:] = value.(EPEC[:DCH_ST])
    mod_dat_str.ESS_STR.CH[:] = value.(EPEC[:CH_ST])
    #CALCULATING PROFITS OF THE EPEC
    global profit_ESS_EPEC, profit_GEN_EPEC, profit_WIND_EPEC = compute_profits(EPEC)
    global iter_diag = 0 #iteration counter
    global tol_diag = 1.0 #We consider any deviation smaller than this as negligible
    global max_iter_diag = 10

    while iter_diag <= max_iter_diag
        println("...Starting the diagonalization loop...")
        #####INIT ESS MILP######
        current = "ESS_STR"
        objective_type = "dual"
        optimizer = with_optimizer(Gurobi.Optimizer, OutputFlag=0)

        println("SOLUTION METHOD:
                    ....Agent: $current
                    ....Objective_type: $objective_type")

        #INIT MODEL
        if bidding == "complex"
            MILP_ESS = build_MILP_complex(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        elseif bidding ==  "quantity"
            MILP_ESS = build_MILP_quantity(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        else
            println("ERROR: Pick 'complex' or 'quantity' for bidding.")
        end

        #UPDATING OTHER UL AGENT DECISIONS
        for t in t_steps
            unfix(MILP_ESS[:G_ST][t])
            unfix(MILP_ESS[:W_ST][t])
            fix(MILP_ESS[:G_ST][t], mod_dat_str.GEN_STR.G[t]; force = true)
            fix(MILP_ESS[:W_ST][t], mod_dat_str.WIND_STR.G[t]; force = true)
        end

        #OPTIMIZE
        optimize!(MILP_ESS)
        objective_value(MILP_ESS)

        #READ MILP RESULT
        global profit_ESS_MILP_ESS, profit_GEN_MILP_ESS, profit_RES_MILP_ESS = compute_profits(MILP_ESS)
        println("ESS_PROFIT:$profit_ESS_MILP_ESS")

        #####INIT GEN MILP######
        current = "GEN_STR"
        println("SOLUTION METHOD:
                    ....Agent: $current
                    ....Objective_type: $objective_type")

        #INIT MODEL
        if bidding == "complex"
            MILP_GEN = build_MILP_complex(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        elseif bidding ==  "quantity"
            MILP_GEN = build_MILP_quantity(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        else
            println("ERROR: Pick 'complex' or 'quantity' for bidding.")
        end

        #UPDATING OTHER UL AGENT DECISIONS
        for t in t_steps
            unfix(MILP_GEN[:DCH_ST][t])
            unfix(MILP_GEN[:CH_ST][t])
            unfix(MILP_GEN[:W_ST][t])
            fix(MILP_GEN[:DCH_ST][t], mod_dat_str.ESS_STR.DCH[t]; force = true)
            fix(MILP_GEN[:CH_ST][t], mod_dat_str.ESS_STR.CH[t]; force = true)
            fix(MILP_GEN[:W_ST][t], mod_dat_str.WIND_STR.G[t]; force = true)
        end

        #OPTIMIZE
        optimize!(MILP_GEN)
        objective_value(MILP_GEN)

        #READ MILP RESULT
        global profit_ESS_MILP_GEN, profit_GEN_MILP_GEN, profit_RES_MILP_GEN = compute_profits(MILP_GEN)
        println("GEN_PROFIT:$profit_GEN_MILP_GEN")

        #####INIT WIND MILP######
        current = "WIND_STR"
        println("SOLUTION METHOD:
                    ....Agent: $current
                    ....Objective_type: $objective_type")

        #INIT MODEL
        if bidding == "complex"
            MILP_WIND = build_MILP_complex(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        elseif bidding ==  "quantity"
            MILP_WIND = build_MILP_quantity(t_steps, non_str_gens,
            non_str_res, current, optimizer, objective_type)
        else
            println("ERROR: Pick 'complex' or 'quantity' for bidding.")
        end

        #UPDATING OTHER UL AGENT DECISIONS
        for t in t_steps
            unfix(MILP_WIND[:DCH_ST][t])
            unfix(MILP_WIND[:CH_ST][t])
            unfix(MILP_WIND[:G_ST][t])
            fix(MILP_WIND[:DCH_ST][t], mod_dat_str.ESS_STR.DCH[t]; force = true)
            fix(MILP_WIND[:CH_ST][t], mod_dat_str.ESS_STR.CH[t]; force = true)
            fix(MILP_WIND[:G_ST][t], mod_dat_str.GEN_STR.G[t]; force = true)
        end

        #OPTIMIZE
        optimize!(MILP_WIND)
        objective_value(MILP_WIND)

        #READ MILP RESULT
        global profit_ESS_MILP_WIND, profit_GEN_MILP_WIND, profit_RES_MILP_WIND = compute_profits(MILP_WIND)
        println("WIND_PROFIT:$profit_RES_MILP_WIND")
        ######SAVE RESULTS#########
        #USE DF HERE
        println("FINAL_PROFITS:
                        PROFIT_ESS: $profit_ESS_MILP_ESS,
                        PROFIT_GEN: $profit_GEN_MILP_GEN,
                        PROFIT_RES: $profit_RES_MILP_WIND")
        ####UPDATING DECISIONS#########
        #Update mod_dat_str
        mod_dat_str.GEN_STR.G[:] = value.(MILP_GEN[:G_ST])
        mod_dat_str.WIND_STR.G[:] = value.(MILP_WIND[:W_ST])
        mod_dat_str.ESS_STR.DCH[:] = value.(MILP_ESS[:DCH_ST])
        mod_dat_str.ESS_STR.CH[:] = value.(MILP_ESS[:CH_ST])

        #####CALCULATING ERRORS########
        #W.R.T. thair dec. vars
        #NOTE: THIS MAY BE MORE REPRESENTATIVE
        # global err_bids
        # value.(MILP_GEN[:G_ST]) - value.(EPEC[:G_ST])
        global err = 0
        if iter_diag == 0
            err_ESS = profit_ESS_MILP_ESS - profit_ESS_EPEC
            err_GEN = profit_GEN_MILP_GEN - profit_GEN_EPEC
            err_WIND = profit_RES_MILP_WIND - profit_WIND_EPEC
            if err_ESS > 0
                err += err_ESS
            end
            if err_GEN > 0
                err += err_GEN
            end
            if err_WIND > 0
                err += err_WIND
            end
        else
            err_ESS = profit_ESS_MILP_ESS - profit_ESS_MILP_ESS_0
            err_GEN = profit_GEN_MILP_GEN - profit_GEN_MILP_GEN_0
            err_WIND = profit_RES_MILP_WIND - profit_RES_MILP_WIND_0
            if err_ESS > 0
                err += err_ESS
            end
            if err_GEN > 0
                err += err_GEN
            end
            if err_WIND > 0
                err += err_WIND
            end
        end
        #IF TOLERANCE HAS REACHED -> BREAK
        if err <= tol_diag
            println("Diagonalization has converged to a solution after the $iter_diag-th iteration!")
            if iter_diag == 0
                println("
                @@@HURRAY: EPEC SOLUTION WAS NASH-EQUILIBRIUM@@@
                ")
            end
            println("PROFITS:
                                PROFIT_ESS: $profit_ESS_MILP_ESS,
                                PROFIT_GEN: $profit_GEN_MILP_GEN,
                                PROFIT_RES: $profit_RES_MILP_WIND")
            println("DEVIATIONS:
                            DEV_ESS: $err_ESS,
                            DEV_GEN: $err_GEN,
                            DEV_RES: $err_WIND")
            break

        end

        println("DEVIATIONS:
                        DEV_ESS: $err_ESS,
                        DEV_GEN: $err_GEN,
                        DEV_RES: $err_WIND")

        global profit_ESS_MILP_ESS_0 = profit_ESS_MILP_ESS
        global profit_GEN_MILP_GEN_0 = profit_GEN_MILP_GEN
        global profit_RES_MILP_WIND_0 = profit_RES_MILP_WIND
        iter_diag += 1
    end
    return iter_diag, profit_ESS_MILP_ESS, profit_GEN_MILP_GEN, profit_RES_MILP_WIND
end


text:
struct MountainCarAgnosticSysIDAgent
    mountaincar::MountainCar
    model::MountainCar
    model_class::Array{Tuple{MountainCar, MountainCarParameters}}
    horizon::Int64
end

function MountainCarAgnosticSysIDAgent(
    mountaincar::MountainCar,
    model::MountainCar,
    horizon::Int64
)
    model_class = generate_model_class()
    MountainCarAgnosticSysIDAgent(mountaincar, model, model_class, horizon)
end

function generate_model_class()::Array{Tuple{MountainCar, MountainCarParameters}}
    # theta1s = -0.0024:-0.00001:-0.0026
    # theta2s = 2.9:0.01:3.1
    # models = []
    # for theta1 in theta1s
    #     for theta2 in theta2s
    #         push!(models, MountainCarParameters(theta1, theta2))
    #     end
    # end
    # The second parameter below is good for 0.03 rock_c
    models = [
        (MountainCar(0.0), MountainCarParameters(-0.0025, 3)),
        (MountainCar(0.03), MountainCarParameters(-0.003, 3.5)),
    ]
    models
end

function run(
    agent::MountainCarAgnosticSysIDAgent;
    debug = false,
    value_aware = true,
)
    rng = MersenneTwister(0)
    # Choose an initial model
    model, params = agent.model_class[1]
    # Choose an initial policy
    policy, values, _ = value_iteration(model, params)
    # Initialize dataset
    all_transitions::Array{Array{MountainCarContTransition}} = []
    # Initialize value functions
    all_values::Array{Array{Float64}} = []
    # Initialize losses
    losses = [0.0 for _ in 1:length(agent.model_class)]

    # Start iterations
    n = 1000  # number of iterations
    m = 10  # number of rollouts using policy
    p = 10  # number of rollouts using random policy with random start
    for i in 1:n
        push!(all_transitions, [])
        # Do rollouts in the real world using policy to collect data
        for _ in 1:m
            episode = simulate_episode(
                agent.mountaincar,
                true_params,
                agent.horizon;
                policy = policy,
                rng = rng,
            )
            all_transitions[i] = vcat(all_transitions[i], episode)
        end
        # Do rollouts in the real world using random policy to collect data
        for _ in 1:p
            episode = simulate_episode(
                agent.mountaincar,
                true_params,
                agent.horizon;
                rng = rng,
            )
            all_transitions[i] = vcat(all_transitions[i], episode)
        end
        # Store value function
        push!(all_values, values)
        # Collected data, now update model
        for j in 1:length(agent.model_class)
            # Compute loss
            loss = 0.0
            if value_aware
                loss = compute_model_advantage_loss(
                    agent,
                    all_transitions[i],
                    all_values[i],
                    agent.model_class[j],
                )
            else
                loss = compute_l2_loss(
                    all_transitions[i],
                    agent.model_class[j],
                )
            end
            loss = loss / (m + p)
            # println(vec(model_params), " ", loss)
            losses[j] += loss
        end
        # Find best model
        # TODO: Doing FTL now, but need to do better
        model, params = agent.model_class[argmin(losses)]
        # println(vec(params), " ", minimum(losses))
        println([(agent.model_class[j][1].rock_c, vec(agent.model_class[j][2]), losses[j]) for j in 1:length(agent.model_class)], "\t",
                argmin(losses))
        # Compute corresponding policy and values
        policy, values, _ = value_iteration(model, params)
    end
end

function compute_model_advantage_loss(
    agent::MountainCarAgnosticSysIDAgent,
    transitions::Array{MountainCarContTransition},
    values::Array{Float64},
    model_params::Tuple{MountainCar, MountainCarParameters},
)
    # _, values, _ = value_iteration(agent.mountaincar, true_params)
    loss = 0.0
    model, params = model_params
    for transition in transitions
        predicted_state, _ = step(
            model,
            transition.initial_state,
            transition.action,
            params,
        )
        loss += abs(
            values[
                cont_state_to_idx(
                    model,
                    transition.final_state,
                )
            ] -
            values[
                cont_state_to_idx(
                    model,
                    predicted_state,
                )
            ]
        )
    end
    loss
end

function compute_l2_loss(
    transitions::Array{MountainCarContTransition},
    model_params::Tuple{MountainCar, MountainCarParameters},
)
    model, params = model_params
    loss = 0.0
    for transition in transitions
        predicted_state, _ = step(
            model,
            transition.initial_state,
            transition.action,
            params,
        )
        loss += (predicted_state.position - transition.final_state.position)^2 + (predicted_state.speed - transition.final_state.speed)^2
    end
    loss
end


text:
import Hyperopt
using Documenter, Nonconvex

makedocs(
    sitename="Nonconvex.jl",
    pages = [
        "Getting started" => "index.md",
        "Problem definition" => "problem.md",
        "Algorithms" => [
            "Overview" => "algorithms/algorithms.md",
            "algorithms/mma.md",
            "algorithms/ipopt.md",
            "algorithms/nlopt.md",
            "algorithms/auglag.md",
            "algorithms/minlp.md",
            "algorithms/hyperopt.md",
            "algorithms/surrogate.md",
        ],
        "Gradients, Jacobians and Hessians" => "gradients.md",
    ],
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/mohamed82008/Nonconvex.jl.git",
        push_preview=true,
    )
end


text:
"""
Seasonal Hierarchical Hidden Markov Models for Julia.

Extansion of the [HMMBase](https://maxmouchet.github.io/HMMBase.jl/stable/) package.
"""
module SHHMM

# TODO!  Package SHHMM does not have CyclicArrays in its dependencies
using Distributed
using ArgCheck
using Distributions
using LinearAlgebra
using JuMP
using Ipopt
using LsqFit
using SpecialFunctions
using Base: OneTo
using Random: AbstractRNG, GLOBAL_RNG
using StatsFuns: logsumexp
using StatsBase: Weights
using ShiftedArrays: lag, lead
using CyclicArrays: CyclicArray
using Dates: now

# Extended functions
import Base: ==, copy, rand, size
import Distributions: MixtureModel, fit_mle, loglikelihood

export
    # shhmm.jl
    AbstractHMM,
    HMM,
    copy,
    rand,
    size,
    nparams,
    permute,
    statdists,
    istransmat,
    istransmats,
    likelihoods,
    # messages.jl
    forward,
    backward,
    posteriors,
    loglikelihood,
    # likelihoods.jl
    loglikelihoods,
    # mle.jl
    fit_mle,
    # viterbi.jl
    viterbi,
    # utilities.jl,
    gettransmat,
    randtransmat,
    remapseq,
    # PeriodicHMM
    PeriodicHMM,
    # # HierarchicalPeriodicHMM
    HierarchicalPeriodicHMM,
    bin2digit,
    dayx,
    conditional_to,
    # Slice fit for initialization
    sort_wrt_ref!,
    fit_mle_all_slices,
    # Trig
    polynomial_trigo,
    fit_Y!,
    fit_Q!,
    HierarchicalPeriodicHMM_trig

include("hmm.jl")
include("mle.jl")
include("mle_init.jl")
include("messages.jl")
include("viterbi.jl")
include("likelihoods.jl")
include("utilities.jl")
include("experimental.jl")
include("periodichmm.jl")
include("periodichmm_leap.jl")
include("periodichmm_leap_hierarchical.jl")
include("fit_mle_slice.jl")
# include("parametric_BW.jl")
include("parametric_BW_distributed.jl")

include("aux_func.jl")
include("mixture.jl")

end # module


text:
import Plots
using RecipesBase

@recipe function f(qt::QuadTreeMeshes.QuadTree)
  linecolor --> :black
  linewidth --> :1
  seriestype := :path
  grid := false

  x = zeros(0)
  y = zeros(0)

  for qte in qt.elements
    bb = QuadTreeMeshes.get_element_bounding_box(qt, qte)
    xe = [bb.x, bb.x + bb.w, bb.x + bb.w, bb.x, bb.x, NaN]
    ye = [bb.y, bb.y, bb.y + bb.h, bb.y + bb.h, bb.y, NaN]
    append!(x, xe)
    append!(y, ye)
  end

  x, y
end

@recipe function f(mesh::QuadTreeMeshes.QuadTreeMesh; boundaries_only=false)
  linecolor --> (boundaries_only ? :red : :blue)
  linewidth --> :2
  seriestype := :path
  grid := false

  x = zeros(0)
  y = zeros(0)

  qt = mesh.quadtree

  if boundaries_only
    for element in qt.values
      if !isnull(element)
        for bv in get(element).boundaries
          v1, v2 = bv
          v1pos, v2pos = qt.vertices[v1.vertex], qt.vertices[v2.vertex]

          xe = [v1pos[1], v2pos[1], NaN]
          ye = [v1pos[2], v2pos[2], NaN]

          append!(x, xe)
          append!(y, ye)
        end
      end
    end
  else
    for triangle in mesh.triangles
      vi = triangle.vertex_indices
        v1 = qt.vertices[vi[1]]
        v2 = qt.vertices[vi[2]]
        v3 = qt.vertices[vi[3]]

        xe = [v1[1], v2[1], v3[1], v1[1], NaN]
        ye = [v1[2], v2[2], v3[2], v1[2], NaN]

        append!(x, xe)
        append!(y, ye)
    end
  end

  x, y
end


text:
#=
Challenge #260 [Easy] Garage Door Opener

Description

You just got a new garage door installed by the Automata™ Garage Door Company.
You are having a lot of fun playing with the remote clicker, opening and closing
the door, scaring your pets and annoying the neighbors.

The clicker is a one-button remote that works like this:
If the door is OPEN or CLOSED, clicking the button will cause the door to move,
until it completes the cycle of opening or closing.

Door: Closed -> Button clicked -> Door: Opening -> Cycle complete -> Door: Open.

If the door is currently opening or closing, clicking the button will make the
door stop where it is. When clicked again, the door will go the opposite direction,
until complete or the button is clicked again.
We will assume the initial state is CLOSED.

Formal Inputs & Outputs

Input description

Input will be a series of commands (can be hard coded, no need to parse):

button_clicked
cycle_complete
button_clicked
button_clicked
button_clicked
button_clicked
button_clicked
cycle_complete

Output description

Output should be the state of the door and the input commands, such as:

Door: CLOSED
> Button clicked.
Door: OPENING
> Cycle complete.
Door: OPEN
> Button clicked.
Door: CLOSING
> Button clicked.
Door: STOPPED_WHILE_CLOSING
> Button clicked.
Door: OPENING
> Button clicked.
Door: STOPPED_WHILE_OPENING
> Button clicked.
Door: CLOSING
> Cycle complete.
Door: CLOSED
Notes/Hints

This is an example of a simple Finite State Machine with 6 States and 2 inputs.

Bonus

Bonus challenge - The door has an infrared beam near the bottom, and if something
is breaking the beam, (your car, your cat, or a baby in a stroller) the door will
be BLOCKED and will add the following rules:
If the door is currently CLOSING, it will reverse to OPENING until completely OPEN.
It will remain BLOCKED, however, until the input BLOCK_CLEARED is called.
Any other state, it will remain in that position, until the input BLOCK_CLEARED
is called, and then it will revert back to it's prior state before it was blocked.
Button clicks will be discarded. If the door was already in the process of opening,
it will continue to OPEN until CYCLE_COMPLETE is called.

Bonus Challenge Input

button_clicked
cycle_complete
button_clicked
block_detected
button_clicked
cycle_complete
button_clicked
block_cleared
button_clicked
cycle_complete

Bonus Challenge output:

Door: CLOSED
> Button clicked
Door: OPENING
> Cycle complete
Door: OPEN
> Button Clicked
Door: CLOSING
> Block detected!
Door: EMERGENCY_OPENING
> Button clicked.
Door: EMERGENCY_OPENING
> Cycle complete.
Door: OPEN_BLOCKED
> Button clicked
Door: OPEN_BLOCKED
> Block cleared
Door: OPEN
> Button clicked
Door: CLOSING
> Cycle complete
Door: CLOSED

Finally

Have a good challenge idea?
Consider submitting it to /r/dailyprogrammer_ideas
Thanks to /u/Philboyd_Studge for this challenge idea.

( source: https://redd.it/4cb7eh )
=#

function main( )
    local state::AbstractString = "CLOSED"
    local line::AbstractString = "nothing"
    local finish::Bool = true

    @time while finish
        print( "Door: ", state, "\n> " )

        if ( finish = !eof( STDIN ) ) != false
            line = chomp( readline( STDIN ) )

            if "CLOSED" == state || "OPENED" == state
                state = ( "CLOSED" == state ) ? "OPENING" : "CLOSING"
            elseif "OPENING" == state || "CLOSING" == state
                if "button_clicked" == line
                    state = ( "CLOSING" == state ) ? "STOPPED_WHILE_CLOSING" : "STOPPED_WHILE_OPENING"
                elseif "block_detected" == line
                    state = "EMERGENCY_OPENING"
                else
                    state = ( "CLOSING" == state ) ? "CLOSED" : "OPENED"
                end
            elseif "EMERGENCY_OPENING" == state && "cycle_complete" == line
                state = "OPEN_BLOCKED"
            elseif "OPEN_BLOCKED" == state && "block_cleared" == line
                state = "OPENED"
            elseif "STOPPED_WHILE_CLOSING" == state || "STOPPED_WHILE_OPENING" == state
                state = ( "STOPPED_WHILE_CLOSING" == state ) ? "OPENING" : "CLOSING"
            end
        end
    end
end

main( )

max_stars_repo_name: farmfirst/dailyprogrammer

max_issues_repo_name: farmfirst/dailyprogrammer

max_forks_repo_name: farmfirst/dailyprogrammer


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    KBallSearch(domain, k, ball)

A method that searches `k` nearest neighbors and then filters
these neighbors using a norm `ball`.
"""
struct KBallSearch{D,B,T} <: BoundedNeighborSearchMethod
  # input fields
  domain::D
  k::Int
  ball::B

  # state fields
  tree::T
end

function KBallSearch(domain::D, k::Int, ball::B) where {D,B}
  m  = metric(ball)
  xs = [coordinates(centroid(domain, i)) for i in 1:nelements(domain)]
  tree = m isa MinkowskiMetric ? KDTree(xs, m) : BallTree(xs, m)
  KBallSearch{D,B,typeof(tree)}(domain, k, ball, tree)
end

maxneighbors(method::KBallSearch) = method.k

function search!(neighbors, pₒ::Point, method::KBallSearch; mask=nothing)
  k = method.k
  r = method.ball isa IsotropicBall ? first(radii(method.ball)) : one(coordtype(pₒ))

  inds, dists = knn(method.tree, coordinates(pₒ), k, true)

  # keep neighbors inside ball
  keep = dists .≤ r

  # possibly mask some of the neighbors
  isnothing(mask) || (keep .*= mask[inds])

  nneigh = 0
  @inbounds for i in 1:k
    if keep[i]
      nneigh += 1
      neighbors[nneigh] = inds[i]
    end
  end

  nneigh
end


text:
function brinco1()
    r = rand()
    if r < 0.5
        1
    else
        -1
    end
end

function primera_llegada_a_0(L::Int=100)
    pos = 1
    t = 0

    while pos > 0

        pos = pos + brinco1()

        if pos > L
            pos = L
        end

        t += 1
    end

    t
end

max_stars_repo_name: dpsanders/metodos-monte-carlo

max_issues_repo_name: dpsanders/metodos-monte-carlo

max_forks_repo_name: dpsanders/metodos-monte-carlo


text:
"""
    Mechanism{T}

    Multi-rigid-body system.

    origin: global reference frame represented with Origin
    joints: list of JointConstraint objects
    bodies: list of Body objects
    contacts: list of ContactConstraint objects
    system: graph-based representation for mechanism
    residual_entries: containt entries for linear system residual
    matrix_entries: contains entries for linear system matrix
    diagonal_inverses: contains inverted matrices computing during LU factorization
	data_matrix: contains parameter information that is fixed during simulation
	root_to_leaves: list of node connections traversing from root node to leaves
    timestep: time discretization
    gravity: force vector resulting from gravitational potential
    μ: complementarity violation (contact softness)
"""
mutable struct Mechanism{T,Nn,Ne,Nb,Ni}
    origin::Origin{T}
    joints::Vector{<:JointConstraint{T}}
    bodies::Vector{Body{T}}
    contacts::Vector{<:ContactConstraint{T}}

    system::System{Nn}
    residual_entries::Vector{Entry}
    matrix_entries::SparseMatrixCSC{Entry,Int64}
    diagonal_inverses::Vector{Entry}

	data_matrix::SparseMatrixCSC{Entry,Int64}
	root_to_leaves::Vector{Int64}

    timestep::T
    gravity::SVector{3,T}
    μ::T
end

function Mechanism(origin::Origin{T}, bodies::Vector{Body{T}}, joints::Vector{<:JointConstraint{T}}, contacts::Vector{<:ContactConstraint{T}};
    spring=0.0,
    damper=0.0,
    timestep=0.01,
    gravity=[0.0; 0.0;-9.81]) where T

    # reset ids
    resetGlobalID()

    # check body inertia parameters
    check_body.(bodies)

    # dimensions
    Ne = length(joints)
    Nb = length(bodies)
    Ni = length(contacts)
    Nn = Ne + Nb + Ni

    # nodes
    nodes = [joints; bodies; contacts]

    # set IDs
    global_id!(nodes)
    origin.id = 0

    # graph system
    system = create_system(origin, joints, bodies, contacts)
    residual_entries = deepcopy(system.vector_entries)
    matrix_entries = deepcopy(system.matrix_entries)
    diagonal_inverses = deepcopy(system.diagonal_inverses)

	# data gradient system
	data_matrix = create_data_matrix(joints, bodies, contacts)

	# Node ordering from root to leaves, loop joints at the end
	loop_joints = get_loop_joints(bodies, joints)
	nodes = [origin; bodies; joints; contacts]
	root_to_leaves = root_to_leaves_ordering(nodes, loop_joints,
		    exclude_origin=true,
            exclude_loop_joints=false)

    # springs and dampers
    (minimum(spring) > 0.0 && minimum(damper) > 0.0) && set_spring_damper_values!(joints, spring, damper)

    Mechanism{T,Nn,Ne,Nb,Ni}(origin, joints, bodies, contacts, system, residual_entries,
		matrix_entries, diagonal_inverses, data_matrix, root_to_leaves, timestep, get_gravity(gravity), 0.0)
end

Mechanism(origin::Origin{T}, bodies::Vector{Body{T}}, joints::Vector{<:JointConstraint{T}}; kwargs...) where T =
	Mechanism(origin, bodies, joints, ContactConstraint{T}[]; kwargs...)

function Mechanism(filename::String,
    floating::Bool=false,
    T=Float64;
    kwargs...)
    # parse urdf
    origin, links, joints, loopjoints = parse_urdf(filename, floating, T)

    # create mechanism
    mechanism = Mechanism(origin, links, [joints; loopjoints]; kwargs...)

    # initialize mechanism
    set_parsed_values!(mechanism, loopjoints)

    return mechanism
end

Base.length(mechanism::Mechanism{T,N}) where {T,N} = N


text:
#	GridWhale module
#	Copyright (c) 2020 Kronosaur Productions, LLC. All Rights Reserved.
#
#	This module provide functions for interaction with the GridWhale system.

module GridWhale

export Canvas

include("Types.jl")
include("Canvas.jl")
include("GoogleDrive.jl")

end


text:
test_folder = "bifurcation"
cd(test_folder)

openBF.runSimulation("$test_folder.yml", verbose=true, out_files=false)

p = readdlm("$test_folder"*"_results/d1_P.last")

@test isapprox(minimum(p[:,end]), 8827.0, atol=1)
@test isapprox(mean(p[:,end]), 12622.0, atol=1)
@test isapprox(maximum(p[:,end]), 17695.0, atol=1)

rm("$test_folder"*"_results/", recursive=true)
cd("..")


text:
# Generic function wrappers

# Functions taking one argument
@generated function llvmwrap(::Val{Op}, v1::Vec{N,T1}, ::Type{R} = T1) where {Op,N,T1,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    flags = [""]
    if Op in FASTOPS# && T1 <: FloatingTypes
        push!(flags, fastflags(promote_type(T1,R)))
    end
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1)")
        push!(instrs, "%res = call" * join(flags, " ") * " $vtypr $ins($vtyp1 %0)")
    else
        if Op === :~
            @assert T1 <: IntegerTypes
            otherval = -1
        elseif Op === :inv
            @assert T1 <: FloatingTypes
            otherval = 1.0
        else
            otherval = 0
        end
        otherarg = llvmconst(N, T1, otherval)
        push!(instrs, "%res = $ins" * join(flags, " ") * " $vtyp1 $otherarg, %0")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,R}, Tuple{Vec{N,T1}}, v1)
    end
end

# Functions taking one Bool argument
@generated function llvmwrap(::Val{Op}, v1::Vec{N,Bool}, ::Type{Bool} = Bool) where {Op,N}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    ins = llvmins(Op, N, Bool)
    decls = String[]
    instrs = String[]
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$N x i1>")
    otherarg = llvmconst(N, Bool, true)
    push!(instrs, "%res = $ins <$N x i1> $otherarg, %arg1")
    push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,Bool}, Tuple{Vec{N,Bool}}, v1)
    end
end

# Functions taking two arguments
@generated function llvmwrap(
    ::Val{Op}, v1::Vec{N,T1}, v2::Vec{N,T2}, ::Type{R} = T1
) where {Op,N,T1,T2,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    flags = [""]
    T = promote_type(T1,T2,R)
    if Op in FASTOPS# && T <: FloatingTypes
        push!(flags, fastflags(T))
    end
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2)")
        push!(instrs, "%res = call" * join(flags, " ") * " $vtypr $ins($vtyp1 %0, $vtyp2 %1)")
    else
        push!(instrs, "%res = $ins" * join(flags, " ") * " $vtyp1 %0, %1")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,R}, Tuple{Vec{N,T1}, Vec{N,T2}},
            v1, v2)
    end
end
@generated function llvmwrap_notfast(::Val{Op}, v1::Vec{N,T1},
        v2::Vec{N,T2}, ::Type{R} = T1) where {Op,N,T1,T2,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2)")
        push!(instrs, "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,R}, Tuple{Vec{N,T1}, Vec{N,T2}},
            v1, v2)
    end
end

# Functions taking two arguments, returning Bool
@generated function llvmwrap(::Val{Op}, v1::Vec{N,T1},
        v2::Vec{N,T2}, ::Type{Bool}) where {Op,N,T1,T2}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    abtyp = "[$N x $btyp]"
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    atyp1 = "[$N x $typ1]"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    atyp2 = "[$N x $typ2]"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    # if false && N == 1
        # append!(instrs, array2vector("%arg1", N, typ1, "%0", "%arg1arr"))
        # append!(instrs, array2vector("%arg2", N, typ2, "%1", "%arg2arr"))
        # push!(instrs, "%cond = $ins $vtyp1 %arg1, %arg2")
        # push!(instrs, "%res = zext <$N x i1> %cond to $vbtyp")
        # append!(instrs, vector2array("%resarr", N, btyp, "%res"))
        # push!(instrs, "ret $abtyp %resarr")
    # else
    push!(instrs, "%res = $ins $vtyp1 %0, %1")
    push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    # end
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,Bool}, Tuple{Vec{N,T1}, Vec{N,T2}},
            v1, v2)
    end
end
# Functions taking two arguments, returning a bitmask
@generated function llvmwrap_bitmask(::Val{Op}, v1::Vec{N,T1}, v2::Vec{N,T2}) where {Op,N,T1,T2}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    abtyp = "[$N x $btyp]"
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    atyp1 = "[$N x $typ1]"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    atyp2 = "[$N x $typ2]"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    # if T1 <: FloatingTypes && T2 <: FloatingTypes
        # push!(instrs, "%res = $ins reassoc $vtyp1 %0, %1")
    # else
    push!(instrs, "%res = $ins $vtyp1 %0, %1")
    # end
    # push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")

    maskbits = max(8, VectorizationBase.nextpow2(N))
    bitcastname = maskbits == N ? "resu" : "resutrunc"
    push!(instrs, "%$(bitcastname) = bitcast <$N x i1> %res to i$N")
    if maskbits != N
        push!(instrs, "%resu = zext i$N %$(bitcastname) to i$maskbits")
    end
    push!(instrs, "ret i$maskbits %resu")
    julia_mask_type = VectorizationBase.mask_type(maskbits)
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            $julia_mask_type, Tuple{Vec{N,T1}, Vec{N,T2}},
            v1, v2)
    end
end

# Functions taking two Bool arguments, returning Bool
@generated function llvmwrap(::Val{Op}, v1::Vec{N,Bool},
        v2::Vec{N,Bool}, ::Type{Bool} = Bool) where {Op,N}
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$N x $btyp>"
    ins = llvmins(Op, N, Bool)
    decls = String[]
    instrs = String[]
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$N x i1>")
    push!(instrs, "%arg2 = trunc $vbtyp %1 to <$N x i1>")
    push!(instrs, "%res = $ins <$N x i1> %arg1, %arg2")
    push!(instrs, "%resb = zext <$N x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,Bool}, Tuple{Vec{N,Bool}, Vec{N,Bool}},
            v1, v2)
    end
end

# Functions taking three arguments
@generated function llvmwrap(::Val{Op}, v1::Vec{N,T1},
        v2::Vec{N,T2}, v3::Vec{N,T3}, ::Type{R} = T1) where {Op,N,T1,T2,T3,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typ3 = llvmtype(T3)
    vtyp3 = "<$N x $typ3>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2, $vtyp3)")
        push!(instrs,
            "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1, $vtyp3 %2)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1, %2")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,R},
            Tuple{Vec{N,T1}, Vec{N,T2}, Vec{N,T3}},
            v1, v2, v3)
    end
end
@generated function llvmwrap_fast(
    ::Val{Op}, v1::Vec{N,T1}, v2::Vec{N,T2}, v3::Vec{N,T3}, ::Type{R} = T1
) where {Op,N,T1,T2,T3,R}
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$N x $typ2>"
    typ3 = llvmtype(T3)
    vtyp3 = "<$N x $typ3>"
    typr = llvmtype(R)
    vtypr = "<$N x $typr>"
    ins = llvmins(Op, N, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2, $vtyp3)")
        push!(instrs,
            "%res = call fast $vtypr $ins($vtyp1 %0, $vtyp2 %1, $vtyp3 %2)")
    else
        push!(instrs, "%res = $ins fast $vtyp1 %0, %1, %2")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,R},
            Tuple{Vec{N,T1}, Vec{N,T2}, Vec{N,T3}},
            v1, v2, v3)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::Vec{N,T}, ::Val{I}
) where {Op,N,T,I}
    @assert isa(Op, Symbol)
    if I >= 0
        op = Op
        i = I
    else
        if Op === :>> || Op === :>>>
            op = :<<
        else
            @assert Op === :<<
            if T <: Unsigned
                op = :>>>
            else
                op = :>>
            end
        end
        i = -I
    end
    @assert op in (:<<, :>>, :>>>)
    @assert i >= 0
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(op, N, T)
    decls = String[]
    instrs = String[]
    nbits = 8*sizeof(T)
    if (op === :>> && T <: IntTypes) || i < nbits
        count = llvmconst(N, T, min(nbits-1, i))
        push!(instrs, "%res = $ins $vtyp %0, $count")
        push!(instrs, "ret $vtyp %res")
    else
        push!(instrs, "return $vtyp zeroinitializer")
    end
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,T}, Tuple{Vec{N,T}}, v1)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::Vec{N,T}, x2::Unsigned
) where {Op,N,T}
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(Op, N, T)
    decls = String[]
    instrs = String[]
    append!(instrs, scalar2vector("%count", N, typ, "%1"))
    nbits = 8*sizeof(T)
    push!(instrs, "%res = $ins $vtyp %0, %count")
    # push!(instrs, "%inbounds = icmp ult $typ %1, $nbits")
    # if Op === :>> && T <: IntTypes
    #     nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
    #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
    #     push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp %limit")
    # else
    #     push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
    # end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        # Note that this function might be called with out-of-bounds
        # values for x2, assuming that the results are then ignored
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,T}, Tuple{Vec{N,T}, T}, v1, x2 % T)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::Vec{N,T}, x2::Integer
) where {Op,N,T}
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}()
    ValNegOp = Val{NegOp}()
    quote
        $(Expr(:meta, :inline))
        x2 < 0 ? llvmwrapshift($ValNegOp, v1, unsigned(-x2)) : llvmwrapshift($ValOp, v1, unsigned(x2))
    end
end

@generated function llvmwrapshift(
    ::Val{Op},
    v1::Vec{N,T},
    v2::Vec{N,U}
) where {Op,N,T,U<:UIntTypes}
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$N x $typ>"
    ins = llvmins(Op, N, T)
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = $ins $vtyp %0, %1")
    # push!(instrs, "%tmp = $ins $vtyp %0, %1")
    # nbits = llvmconst(N, T, 8*sizeof(T))
    # push!(instrs, "%inbounds = icmp ult $vtyp %1, $nbits")
    # if Op === :>> && T <: IntTypes
    #     nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
    #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
    #     push!(instrs,
    #         "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp %limit")
    # else
    #     push!(instrs,
    #         "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
    # end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{N,T}, Tuple{Vec{N,T}, Vec{N,T}},
            v1, vrem(v2, Vec{N,T}))
    end
end

@generated function llvmwrapshift(
    ::Val{Op},
    v1::Vec{N,T},
    v2::Vec{N,U}
) where {Op,N,T,U<:IntegerTypes}
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}()
    ValNegOp = Val{NegOp}()
    quote
        $(Expr(:meta, :inline))
        vifelse(
            vgreater_or_equal(v2, 0),
            llvmwrapshift($ValOp, v1, vrem(v2, Vec{N,unsigned(U)})),
            llvmwrapshift($ValNegOp, v1, vrem(vsub(v2), Vec{N,unsigned(U)}))
        )
    end
end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::Vec{N,T}, ::Val{I}
# ) where {Op,N,T,I}
#     @assert isa(Op, Symbol)
#     if I >= 0
#         op = Op
#         i = I
#     else
#         if Op === :>> || Op === :>>>
#             op = :<<
#         else
#             @assert Op === :<<
#             if T <: Unsigned
#                 op = :>>>
#             else
#                 op = :>>
#             end
#         end
#         i = -I
#     end
#     @assert op in (:<<, :>>, :>>>)
#     @assert i >= 0
#     typ = llvmtype(T)
#     vtyp = "<$N x $typ>"
#     ins = llvmins(op, N, T)
#     decls = String[]
#     instrs = String[]
#     nbits = 8*sizeof(T)
#     if (op === :>> && T <: IntTypes) || i < nbits
#         count = llvmconst(N, T, min(nbits-1, i))
#         push!(instrs, "%res = $ins $vtyp %0, $count")
#         push!(instrs, "ret $vtyp %res")
#     else
#         push!(instrs, "return $vtyp zeroinitializer")
#     end
#     quote
#         $(Expr(:meta, :inline))
#         Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{N,T}, Tuple{Vec{N,T}}, v1)
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::Vec{N,T}, x2::Unsigned
# ) where {Op,N,T}
#     @assert isa(Op, Symbol)
#     typ = llvmtype(T)
#     vtyp = "<$N x $typ>"
#     ins = llvmins(Op, N, T)
#     decls = String[]
#     instrs = String[]
#     append!(instrs, scalar2vector("%count", N, typ, "%1"))
#     nbits = 8*sizeof(T)
#     push!(instrs, "%tmp = $ins $vtyp %0, %count")
#     push!(instrs, "%inbounds = icmp ult $typ %1, $nbits")
#     if Op === :>> && T <: IntTypes
#         nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
#         push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
#         push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp %limit")
#     else
#         push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
#     end
#     push!(instrs, "ret $vtyp %res")
#     quote
#         $(Expr(:meta, :inline))
#         # Note that this function might be called with out-of-bounds
#         # values for x2, assuming that the results are then ignored
#         Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{N,T}, Tuple{Vec{N,T}, T}, v1, x2 % T)
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::Vec{N,T}, x2::Integer
# ) where {Op,N,T}
#     if Op === :>> || Op === :>>>
#         NegOp = :<<
#     else
#         @assert Op === :<<
#         if T <: Unsigned
#             NegOp = :>>>
#         else
#             NegOp = :>>
#         end
#     end
#     ValOp = Val{Op}()
#     ValNegOp = Val{NegOp}()
#     quote
#         $(Expr(:meta, :inline))
#         x2 < 0 ? llvmwrapshift($ValNegOp, v1, unsigned(-x2)) : llvmwrapshift($ValOp, v1, unsigned(x2))
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op},
#     v1::Vec{N,T},
#     v2::Vec{N,U}
# ) where {Op,N,T,U<:UIntTypes}
#     @assert isa(Op, Symbol)
#     typ = llvmtype(T)
#     vtyp = "<$N x $typ>"
#     ins = llvmins(Op, N, T)
#     decls = String[]
#     instrs = String[]
#     push!(instrs, "%res = $ins $vtyp %0, %1")
#     # push!(instrs, "%tmp = $ins $vtyp %0, %1")
#     # nbits = llvmconst(N, T, 8*sizeof(T))
#     # push!(instrs, "%inbounds = icmp ult $vtyp %1, $nbits")
#     # if Op === :>> && T <: IntTypes
#     #     nbits1 = llvmconst(N, T, 8*sizeof(T)-1)
#     #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
#     #     push!(instrs,
#     #         "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp %limit")
#     # else
#     #     push!(instrs,
#     #         "%res = select <$N x i1> %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
#     # end
#     push!(instrs, "ret $vtyp %res")
#     quote
#         $(Expr(:meta, :inline))
#         Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{N,T}, Tuple{Vec{N,T}, Vec{N,T}},
#             v1, vrem(v2, Vec{N,T}))
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op},
#     v1::Vec{N,T},
#     v2::Vec{N,U}
# ) where {Op,N,T,U<:IntegerTypes}
#     if Op === :>> || Op === :>>>
#         NegOp = :<<
#     else
#         @assert Op === :<<
#         if T <: Unsigned
#             NegOp = :>>>
#         else
#             NegOp = :>>
#         end
#     end
#     ValOp = Val{Op}()
#     ValNegOp = Val{NegOp}()
#     quote
#         $(Expr(:meta, :inline))
#         vifelse(
#             vgreater_or_equal(v2, 0),
#             llvmwrapshift($ValOp, v1, vrem(v2, Vec{N,unsigned(U)})),
#             llvmwrapshift($ValNegOp, v1, vrem(vsub(v2), Vec{N,unsigned(U)}))
#         )
#     end
# end