text:
module straight_duct_examples
using FinEtools
using FinEtoolsAcoustics
using UnicodePlots

function straight_duct_Q8_example()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    println("""

    Straight duct with anechoic termination.
    Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    Both real and imaginary components of the pressure should have amplitude of
    rho*c = $(rho*c).

    Two-dimensional model. Quadratic qquadrilateral mesh.
    """)

    fens,fes  =  Q8block(Lx,Ly,n,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes,facing = true, direction = [-1.0 0.0])
    L10 = selectelem(fens,bfes,facing = true, direction = [+1.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(ComplexF64,size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegDomain(fes, GaussRule(2, 3)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegDomain(subset(bfes,L10),GaussRule(1, 3)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegDomain(subset(bfes,L0), GaussRule(1,  3)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    println("Pressure amplitude bounds: ")
    println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")

    println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct_Q8.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.Q8;  scalars = [("Pressure", scalars)])
    @async run(`"paraview.exe" $File`)

    ix = sortperm(geom.values[nLx,1])
    # @pgf a = Axis({
    #           xlabel = "x",
    #           ylabel = L"$\mathrm{Re}P$, $\mathrm{Im}P$",
    #           title = "Steady-state pressure"
    #       },
    #       Plot(Table([:x => geom.values[nLx,1][ix], :y => real(P.values)[nLx][ix]])), LegendEntry("real"),
    #       Plot(Table([:x => geom.values[nLx,1][ix], :y => imag(P.values)[nLx][ix]])), LegendEntry("imag"))
    plt = lineplot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], canvas = DotCanvas, title = "Steady-state pressure", name = "real", xlabel = "x", ylabel = "P")
    plt = lineplot!(plt, geom.values[nLx,1][ix], imag(P.values)[nLx][ix], name = "imag")
    display(plt)

    true

end # straight_duct_Q8_example

function straight_duct_T3_example()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    println("""

    Straight duct with anechoic termination.
    Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    Both real and imaginary components of the pressure should have amplitude of
    rho*c = $(rho*c).

    Triangle mesh.
    """)

    fens,fes  =  T3block(Lx,Ly,n,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes,facing = true, direction = [-1.0 0.0])
    L10 = selectelem(fens,bfes,facing = true, direction = [+1.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(ComplexF64,size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegDomain(fes, TriRule(1)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegDomain(subset(bfes,L10),GaussRule(1, 2)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegDomain(subset(bfes,L0), GaussRule(1,  2)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    println("Pressure amplitude bounds: ")
    println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")

    println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.T3; scalars = [("Pressure", scalars)])
    @async run(`"paraview.exe" $File`)

    ix = sortperm(geom.values[nLx,1])
    # @pgf a = Axis({
    #           xlabel = "x",
    #           ylabel = L"$\mathrm{Re}P$, $\mathrm{Im}P$",
    #           title = "Steady-state pressure"
    #       },
    #       Plot(Table([:x => geom.values[nLx,1][ix], :y => real(P.values)[nLx][ix]])), LegendEntry("real"),
    #       Plot(Table([:x => geom.values[nLx,1][ix], :y => imag(P.values)[nLx][ix]])), LegendEntry("imag"))
    # display(a)
    plt = lineplot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], canvas = DotCanvas, title = "Steady-state pressure", name = "real", xlabel = "x", ylabel = "P")
    plt = lineplot!(plt, geom.values[nLx,1][ix], imag(P.values)[nLx][ix], name = "imag")
    display(plt)

    true

end # straight_duct_T3_example

function allrun()
    println("#####################################################")
    println("# straight_duct_Q8_example ")
    straight_duct_Q8_example()
    println("#####################################################")
    println("# straight_duct_T3_example ")
    straight_duct_T3_example()
end # function allrun

end # module straight_duct_examples


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events
using AWS.Compat
using AWS.UUIDs

"""
    create_detector_model(detector_model_definition, detector_model_name, role_arn)
    create_detector_model(detector_model_definition, detector_model_name, role_arn, params::Dict{String,<:Any})

Creates a detector model.

# Arguments
- `detector_model_definition`: Information that defines how the detectors operate.
- `detector_model_name`: The name of the detector model.
- `role_arn`: The ARN of the role that grants permission to AWS IoT Events to perform its
  operations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"detectorModelDescription"`: A brief description of the detector model.
- `"evaluationMethod"`: Information about the order in which events are evaluated and how
  actions are executed.
- `"key"`: The input attribute key used to identify a device or system to create a detector
  (an instance of the detector model) and then to route each input received to the
  appropriate detector (instance). This parameter uses a JSON-path expression in the message
  payload of each input to specify the attribute-value pair that is used to identify the
  device associated with the input.
- `"tags"`: Metadata that can be used to manage the detector model.
"""
create_detector_model(detectorModelDefinition, detectorModelName, roleArn; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/detector-models", Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "detectorModelName"=>detectorModelName, "roleArn"=>roleArn); aws_config=aws_config)
create_detector_model(detectorModelDefinition, detectorModelName, roleArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/detector-models", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "detectorModelName"=>detectorModelName, "roleArn"=>roleArn), params)); aws_config=aws_config)

"""
    create_input(input_definition, input_name)
    create_input(input_definition, input_name, params::Dict{String,<:Any})

Creates an input.

# Arguments
- `input_definition`: The definition of the input.
- `input_name`: The name you want to give to the input.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"inputDescription"`: A brief description of the input.
- `"tags"`: Metadata that can be used to manage the input.
"""
create_input(inputDefinition, inputName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/inputs", Dict{String, Any}("inputDefinition"=>inputDefinition, "inputName"=>inputName); aws_config=aws_config)
create_input(inputDefinition, inputName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/inputs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputDefinition"=>inputDefinition, "inputName"=>inputName), params)); aws_config=aws_config)

"""
    delete_detector_model(detector_model_name)
    delete_detector_model(detector_model_name, params::Dict{String,<:Any})

Deletes a detector model. Any active instances of the detector model are also deleted.

# Arguments
- `detector_model_name`: The name of the detector model to be deleted.

"""
delete_detector_model(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/detector-models/$(detectorModelName)"; aws_config=aws_config)
delete_detector_model(detectorModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/detector-models/$(detectorModelName)", params; aws_config=aws_config)

"""
    delete_input(input_name)
    delete_input(input_name, params::Dict{String,<:Any})

Deletes an input.

# Arguments
- `input_name`: The name of the input to delete.

"""
delete_input(inputName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/inputs/$(inputName)"; aws_config=aws_config)
delete_input(inputName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/inputs/$(inputName)", params; aws_config=aws_config)

"""
    describe_detector_model(detector_model_name)
    describe_detector_model(detector_model_name, params::Dict{String,<:Any})

Describes a detector model. If the version parameter is not specified, information about
the latest version is returned.

# Arguments
- `detector_model_name`: The name of the detector model.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"version"`: The version of the detector model.
"""
describe_detector_model(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)"; aws_config=aws_config)
describe_detector_model(detectorModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)", params; aws_config=aws_config)

"""
    describe_detector_model_analysis(analysis_id)
    describe_detector_model_analysis(analysis_id, params::Dict{String,<:Any})

Retrieves execution information about a detector model analysis

# Arguments
- `analysis_id`: The ID of the analysis result that you want to retrieve.

"""
describe_detector_model_analysis(analysisId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/analysis/detector-models/$(analysisId)"; aws_config=aws_config)
describe_detector_model_analysis(analysisId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/analysis/detector-models/$(analysisId)", params; aws_config=aws_config)

"""
    describe_input(input_name)
    describe_input(input_name, params::Dict{String,<:Any})

Describes an input.

# Arguments
- `input_name`: The name of the input.

"""
describe_input(inputName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/inputs/$(inputName)"; aws_config=aws_config)
describe_input(inputName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/inputs/$(inputName)", params; aws_config=aws_config)

"""
    describe_logging_options()
    describe_logging_options(params::Dict{String,<:Any})

Retrieves the current settings of the AWS IoT Events logging options.

"""
describe_logging_options(; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/logging"; aws_config=aws_config)
describe_logging_options(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/logging", params; aws_config=aws_config)

"""
    get_detector_model_analysis_results(analysis_id)
    get_detector_model_analysis_results(analysis_id, params::Dict{String,<:Any})

Retrieves one or more analysis results of the detector model.

# Arguments
- `analysis_id`: The ID of the analysis result that you want to retrieve.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
get_detector_model_analysis_results(analysisId; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/analysis/detector-models/$(analysisId)/results"; aws_config=aws_config)
get_detector_model_analysis_results(analysisId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/analysis/detector-models/$(analysisId)/results", params; aws_config=aws_config)

"""
    list_detector_model_versions(detector_model_name)
    list_detector_model_versions(detector_model_name, params::Dict{String,<:Any})

Lists all the versions of a detector model. Only the metadata associated with each detector
model version is returned.

# Arguments
- `detector_model_name`: The name of the detector model whose versions are returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
list_detector_model_versions(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)/versions"; aws_config=aws_config)
list_detector_model_versions(detectorModelName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models/$(detectorModelName)/versions", params; aws_config=aws_config)

"""
    list_detector_models()
    list_detector_models(params::Dict{String,<:Any})

Lists the detector models you have created. Only the metadata associated with each detector
model is returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
list_detector_models(; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models"; aws_config=aws_config)
list_detector_models(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/detector-models", params; aws_config=aws_config)

"""
    list_inputs()
    list_inputs(params::Dict{String,<:Any})

Lists the inputs you have created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to be returned per request.
- `"nextToken"`: The token that you can use to return the next set of results.
"""
list_inputs(; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/inputs"; aws_config=aws_config)
list_inputs(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/inputs", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags (metadata) you have assigned to the resource.

# Arguments
- `resource_arn`: The ARN of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/tags", Dict{String, Any}("resourceArn"=>resourceArn); aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("GET", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn), params)); aws_config=aws_config)

"""
    put_logging_options(logging_options)
    put_logging_options(logging_options, params::Dict{String,<:Any})

Sets or updates the AWS IoT Events logging options. If you update the value of any
loggingOptions field, it takes up to one minute for the change to take effect. If you
change the policy attached to the role you specified in the roleArn field (for example, to
correct an invalid policy), it takes up to five minutes for that change to take effect.

# Arguments
- `logging_options`: The new values of the AWS IoT Events logging options.

"""
put_logging_options(loggingOptions; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("PUT", "/logging", Dict{String, Any}("loggingOptions"=>loggingOptions); aws_config=aws_config)
put_logging_options(loggingOptions, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("PUT", "/logging", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("loggingOptions"=>loggingOptions), params)); aws_config=aws_config)

"""
    start_detector_model_analysis(detector_model_definition)
    start_detector_model_analysis(detector_model_definition, params::Dict{String,<:Any})

Performs an analysis of your detector model. For more information, see Running detector
model analyses in the AWS IoT Events Developer Guide.

# Arguments
- `detector_model_definition`:

"""
start_detector_model_analysis(detectorModelDefinition; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/analysis/detector-models/", Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition); aws_config=aws_config)
start_detector_model_analysis(detectorModelDefinition, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/analysis/detector-models/", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds to or modifies the tags of the given resource. Tags are metadata that can be used to
manage a resource.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tags`: The new or modified tags for the resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/tags", Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the given tags (metadata) from the resource.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tag_keys`: A list of the keys of the tags to be removed from the resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/tags", Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("DELETE", "/tags", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_detector_model(detector_model_definition, detector_model_name, role_arn)
    update_detector_model(detector_model_definition, detector_model_name, role_arn, params::Dict{String,<:Any})

Updates a detector model. Detectors (instances) spawned by the previous version are deleted
and then re-created as new inputs arrive.

# Arguments
- `detector_model_definition`: Information that defines how a detector operates.
- `detector_model_name`: The name of the detector model that is updated.
- `role_arn`: The ARN of the role that grants permission to AWS IoT Events to perform its
  operations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"detectorModelDescription"`: A brief description of the detector model.
- `"evaluationMethod"`: Information about the order in which events are evaluated and how
  actions are executed.
"""
update_detector_model(detectorModelDefinition, detectorModelName, roleArn; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/detector-models/$(detectorModelName)", Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "roleArn"=>roleArn); aws_config=aws_config)
update_detector_model(detectorModelDefinition, detectorModelName, roleArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("POST", "/detector-models/$(detectorModelName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectorModelDefinition"=>detectorModelDefinition, "roleArn"=>roleArn), params)); aws_config=aws_config)

"""
    update_input(input_definition, input_name)
    update_input(input_definition, input_name, params::Dict{String,<:Any})

Updates an input.

# Arguments
- `input_definition`: The definition of the input.
- `input_name`: The name of the input you want to update.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"inputDescription"`: A brief description of the input.
"""
update_input(inputDefinition, inputName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("PUT", "/inputs/$(inputName)", Dict{String, Any}("inputDefinition"=>inputDefinition); aws_config=aws_config)
update_input(inputDefinition, inputName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events("PUT", "/inputs/$(inputName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("inputDefinition"=>inputDefinition), params)); aws_config=aws_config)

ext: jl


text:
################################################################################
#
#   DEFINITION OF CONCRETE HAMILTONIAN <: AbstractHamiltonian
#   - label type L (of bond labels)
#   - unitcell type U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
#   - bond hamiltonian type HB<:AbstractBondHamiltonian{L,NS} where {NS}
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
struct Hamiltonian{L,UC,HB} <: AbstractHamiltonian{L,UC,HB}

    # the unitcell
    unitcell :: UC

    # the bond hamiltonian
    h_bond   :: HB

end

# export the type
export Hamiltonian

################################################################################
#
#	INTERFACING / ACCESSING HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# Function to obtain the unitcell
function unitcell(
            h :: Hamiltonian{L,UC,HB}
        ) :: UC where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS}}

    # return the unitcell
    return h.unitcell
end

# Function to obtain the bond hamiltonian
function bondHamiltonian(
            h :: Hamiltonian{L,UC,HB}
        ) :: HB where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS}}

    # return the bond hamiltonian
    return h.h_bond
end

# easy constructor
function Hamiltonian(
            unitcell :: UC,
            h_bond   :: HB
        ) :: Hamiltonian{L,UC,HB} where {L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS}}

    # return a newly created hamiltonian
    return Hamiltonian{L,UC,HB}(unitcell, h_bond)
end

function saveHamiltonian(
        h  :: H,
        fn :: AbstractString,
        group :: AbstractString = "hamiltonian"
        ;
        append :: Bool = false
    ) where {L,UC<:AbstractUnitcell,N,HB<:AbstractBondHamiltonian{L,N}, H<:Hamiltonian{L,UC,HB}}

    # determine the mode based on if one wants to append stuff
    if append
        mode = "r+"
    else
        mode = "w"
    end

    # group for bond hamiltonian
    group_hb = group*"/bond_hamiltonian"
    # group for unitcell
    group_uc = group*"/unitcell"

    # open the file in mode
    h5open(fn, mode) do file
        # create the group in which the bonds are saved
        group_h = g_create(file, group)
        # save the parameters
        attrs(group_h)["N"] = N
        attrs(group_h)["L"] = string(L)
        # save the groups
        attrs(group_h)["bond_hamiltonian"] = group_hb
        attrs(group_h)["unitcell"]         = group_uc
    end

    # save bond hamiltonian and unitcell in respective groups (append!)
    saveBondHamiltonian(bondHamiltonian(h), fn, group_hb, append=true)
    saveUnitcell(unitcell(h), fn, group_uc, append=true)

    # return nothing
    return nothing
end

function loadHamiltonian(
        ::Type{H},
        fn :: AbstractString,
        group :: AbstractString = "hamiltonian"
    ) where {L,UC<:AbstractUnitcell,N,HB<:AbstractBondHamiltonian{L,N}, H<:Union{Hamiltonian{L,UC,HB}, Hamiltonian}}

    # read attribute data
    attr_data = h5readattr(fn, group)

    # load bond hamiltonian
    hb = loadBondHamiltonian(fn, attr_data["bond_hamiltonian"])
    # load unitcell
    uc = loadUnitcell(fn, attr_data["unitcell"])

    # return the new hamiltonian
    return Hamiltonian(uc, hb)
end

function loadHamiltonian(
        fn :: AbstractString,
        group :: AbstractString = "hamiltonian"
    )

    return loadHamiltonian(Hamiltonian, fn, group)
end


text:
module MasslessChargedParticlePlots

    using GeometricBase
    using GeometricEquations
    using GeometricProblems.Diagnostics

    using LaTeXStrings
    using Measures: mm
    using RecipesBase

    """
    Plots the solution of a massless charged particle together with the energy error.

    Arguments:
    * `sol <: Solution`
    * `params <: NamedTuple`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `xlims=:auto`: xlims for solution plot
    * `ylims=:auto`: ylims for solution plot
    """
    @userplot Plot_Massless_Charged_Particle
    @recipe function f(p::Plot_Massless_Charged_Particle; nplot=1, nt=:auto, xlims=:auto, ylims=:auto, elims=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Massless charged particle plots should be given two arguments: a solution and a parameter tuple. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        H, ΔH = compute_invariant_error(sol.t, sol.q, (t,q) -> equ.invariants[:h](t,q,params))

        size   := (800,300)
        layout := (1,2)
        legend := :none

        right_margin := 10mm
        bottom_margin := 10mm

        guidefontsize := 14
        tickfontsize  := 10

        @series begin
            subplot := 1

            if sol.nt ≤ 200
                markersize := 5
            else
                markersize  := .5
                markercolor := 1
                linecolor   := 1
                markerstrokewidth := .5
                markerstrokecolor := 1
            end
            if nplot > 1
                seriestype := :scatter
            end

            if latex
                xguide := L"x_1"
                yguide := L"x_2"
            else
            # if backend() == Plots.PlotlyJSBackend()
                xguide := "x₁"
                yguide := "x₂"
            end
            xlims  := xlims
            ylims  := ylims
            aspect_ratio := 1
            sol.q[1,0:nplot:end], sol.q[2,0:nplot:end]
        end

        @series begin
            subplot := 2
            if latex
                xguide := L"t"
                yguide := L"[H(t) - H(0)] / H(0)"
            else
                xguide := "t"
                yguide := "[H(t) - H(0)] / H(0)"
            end
            xlims  := (sol.t[0], Inf)
            ylims  := elims
            yformatter := :scientific
            sol.t[0:nplot:nt], ΔH[0:nplot:nt]
        end
    end

    """
    Plots the solution of a massless charged particle.

    Arguments:
    * `sol <: Solution`
    * `params <: NamedTuple`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `xlims=:auto`: xlims for solution plot
    * `ylims=:auto`: ylims for solution plot
    """
    @userplot Plot_Massless_Charged_Particle_Solution
    @recipe function f(p::Plot_Massless_Charged_Particle_Solution; nplot=1, nt=:auto, xlims=:auto, ylims=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Massless charged particle plots should be given two arguments: a solution and a parameter tuple. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        if ntime(sol) ≤ 200
            markersize := 5
        else
            markersize  := .5
            markercolor := 1
            linecolor   := 1
            markerstrokewidth := .5
            markerstrokecolor := 1
        end
        if nplot > 1
            seriestype := :scatter
        end

        legend := :none
        size := (400,400)

        # if backend() == Plots.PlotlyJSBackend() || backend() == Plots.GRBackend()
            guidefontsize := 18
            tickfontsize  := 12
        # else
        #     guidefontsize := 14
        #     tickfontsize  := 10
        # end

        # solution
        @series begin
            if latex
                xguide := L"x_1"
                yguide := L"x_2"
            else
            # if backend() == Plots.PlotlyJSBackend()
                xguide := "x₁"
                yguide := "x₂"
            end

            xlims  := xlims
            ylims  := ylims
            aspect_ratio := 1
            sol.q[1,0:nplot:nt], sol.q[2,0:nplot:nt]
        end
    end

    """
    Plots time traces of the energy error of a massless charged particle.

    Arguments:
    * `sol <: Solution`
    * `params <: NamedTuple`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    """
    @userplot Plot_Massless_Charged_Particle_Energy_Error
    @recipe function f(p::Plot_Massless_Charged_Particle_Energy_Error; nplot=1, nt=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Massless charged particle plots should be given two arguments: a solution and a parameter tuple. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        H, ΔH = compute_invariant_error(sol.t, sol.q, (t,q) -> equ.invariants[:h](t,q,params))

        size   := (800,200)
        legend := :none
        right_margin := 10mm

        # if backend() == Plots.GRBackend()
            left_margin := 10mm
            bottom_margin := 10mm
        # end

        guidefontsize := 12
        tickfontsize  := 10

        @series begin
            if latex
                xguide := L"t"
                yguide := L"[H(t) - H(0)] / H(0)"
            # if backend() == Plots.PlotlyJSBackend()
            else
                xguide := "t"
                yguide := "[H(t) - H(0)] / H(0)"
            end

            xlims  := (sol.t[0], Inf)
            yformatter := :scientific
            sol.t[0:nplot:nt], ΔH[0:nplot:nt]
        end
    end

    """
    Plots time traces of the momentum error of a massless charged particle.

    Arguments:
    * `sol <: Solution`
    * `params <: NamedTuple`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    * `k=0`: index of momentum component (0 plots all components)
    """
    @userplot Plot_Massless_Charged_Particle_Momentum_Error
    @recipe function f(p::Plot_Massless_Charged_Particle_Momentum_Error; k=0, nplot=1, nt=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Massless charged particle plots should be given two arguments: a solution and a parameter tuple. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        Δp = compute_momentum_error(sol.t, sol.q, sol.p, params)

        ntrace = (k == 0 ?    Δp.nd  :    1 )
        trange = (k == 0 ? (1:Δp.nd) : (k:k))

        size   := (800,200*ntrace)
        legend := :none
        layout := (ntrace,1)

        right_margin := 10mm

        if ntrace == 1 #&& backend() == Plots.GRBackend()
            left_margin := 10mm
        end

        guidefontsize := 12
        tickfontsize  := 10

        for i in trange
            @series begin
                if k == 0
                    subplot := i
                end
                if i == Δp.nd || k ≠ 0
                    if latex
                        xguide := L"t"
                    # if backend() == Plots.PlotlyJSBackend()
                    else
                        xguide := "t"
                    end
                else
                    xaxis := false
                    bottom_margin := 10mm
                end
                # if i < Δp.nd && k == 0
                #     if backend() == Plots.PGFPlotsXBackend()
                #         bottom_margin := -6mm
                #     elseif backend() == Plots.GRBackend()
                #         bottom_margin := -3mm
                #     end
                # end
                # if i > 1 && k == 0
                #     if backend() == Plots.PGFPlotsXBackend()
                #         top_margin := -6mm
                #     elseif backend() == Plots.GRBackend()
                #         top_margin := -3mm
                #     end
                # end
                if latex
                    yguide := latexstring("p_$i (t) - \\vartheta_$i (t)")
                # if backend() == Plots.PlotlyJSBackend()
                else
                    yguide := "p" * subscript(i) * "(t) - ϑ" * subscript(i) * "(t)"
                end
                xlims  := (sol.t[0], Inf)
                yformatter := :scientific
                sol.t[0:nplot:nt], Δp[i,0:nplot:nt]
            end
        end
    end

    """
    Plots time traces of the solution of a massless charged particle trajectory and its energy error.

    Arguments:
    * `sol <: Solution`
    * `params <: NamedTuple`

    Keyword aguments:
    * `nplot=1`: plot every `nplot`th time step
    """
    @userplot Plot_Massless_Charged_Particle_Traces
    @recipe function f(p::Plot_Massless_Charged_Particle_Traces; nplot=1, nt=:auto, xlims=:auto, ylims=:auto, elims=:auto, latex=true)
        if length(p.args) != 2 || !(typeof(p.args[1]) <: Solution) || !(typeof(p.args[2]) <: Equation)
            error("Massless charged particle plots should be given two arguments: a solution and a parameter tuple. Got: $(typeof(p.args))")
        end
        sol = p.args[1]
        equ = p.args[2]
        params = equ.parameters

        if nt == :auto
            nt = ntime(sol)
        end

        if nt > ntime(sol)
            nt = ntime(sol)
        end

        H, ΔH = compute_invariant_error(sol.t, sol.q, (t,q) -> equ.invariants[:h](t,q,params))

        size   := (800,600)
        legend := :none
        right_margin := 10mm

        # traces
        layout := (3,1)

        # if backend() == Plots.PlotlyJSBackend() || backend() == Plots.GRBackend()
            guidefontsize := 14
            tickfontsize  := 10
        # else
        #     guidefontsize := 12
        #     tickfontsize  := 10
        # end

        if latex
            ylabels = (L"x_1", L"x_2")
        # if backend() == Plots.PlotlyJSBackend()
        else
            ylabels = ("x₁", "x₂")
        end

        lims = (xlims, ylims, elims)

        for i in 1:2
            @series begin
                subplot := i
                yguide := ylabels[i]
                xlims  := (sol.t[0], Inf)
                ylims  := lims[i]
                xaxis  := false

                # if i == 1 && backend() == Plots.PGFPlotsXBackend()
                #     bottom_margin := -4mm
                # elseif i == 1 && backend() == Plots.GRBackend()
                #     bottom_margin := -2mm
                # end
                # if i == 2 && backend() == Plots.PGFPlotsXBackend()
                #     top_margin := -2mm
                #     bottom_margin := -2mm
                # elseif i == 2 && backend() == Plots.GRBackend()
                #     top_margin := -1mm
                #     bottom_margin := -1mm
                # end

                sol.t[0:nplot:nt], sol.q[i,0:nplot:nt]
            end
        end

        @series begin
            subplot := 3
            # if backend() == Plots.PGFPlotsXBackend()
            #     top_margin := -4mm
            # elseif backend() == Plots.GRBackend()
            #     top_margin := -2mm
            # end
            if latex
                xguide := L"t"
                yguide := L"[H(t) - H(0)] / H(0)"
            # if backend() == Plots.PlotlyJSBackend()
            else
                xguide := "t"
                yguide := "[H(t) - H(0)] / H(0)"
            end
            xlims := (sol.t[0], Inf)
            ylims := elims
            yformatter := :scientific
            sol.t[0:nplot:nt], ΔH[0:nplot:nt]
        end
    end

end


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
using Pkg, DrWatson

# ╔═╡ e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using SpecialFunctions
end

# ╔═╡ ac3ceabe-1ee4-11eb-3935-313e96aafd62
md"## Listing 2.11"

# ╔═╡ e9fa266e-1eeb-11eb-2630-59930c474ed3
n = 2000

# ╔═╡ d595cf1a-1f62-11eb-3503-753f5260ea5c
begin
	probAgivenB(k) = 1- 1/(k+1)
	probB(k) = 6/(pi*(k+1))^2
end

# ╔═╡ d5960d36-1f62-11eb-01dc-b59f88c34a97
begin
	numerical= sum([probAgivenB(k)*probB(k) for k in 0:n])
	analytic = 1 - 6*zeta(3)/pi^2
	Text("Analytic: $(analytic) Numerical: $(numerical)")
end

# ╔═╡ ea086874-1eeb-11eb-2cb0-4186697d4294
md"## End of listing 2.11"

# ╔═╡ Cell order:
# ╟─ac3ceabe-1ee4-11eb-3935-313e96aafd62
# ╠═de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
# ╠═e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
# ╠═e9fa266e-1eeb-11eb-2630-59930c474ed3
# ╠═d595cf1a-1f62-11eb-3503-753f5260ea5c
# ╠═d5960d36-1f62-11eb-01dc-b59f88c34a97
# ╟─ea086874-1eeb-11eb-2cb0-4186697d4294


text:
"""
x[1] is dx/dt (the velocity of the car)
x[2] is x (the position of the car)
x[3] is d phi/dt (the angular velocity of the pendulum)
x[4] is phi (the angle of the pendulum, 0[rad] is at the top)
u    is Vm (the input voltage for the motor)
"""
function inv_pendulum_rk4(x_t::Vector{T}, u, dt::Number) where {T}
    x_t1 = copy(x_t)
    inv_pendulum_rk4!(x_t1, u, dt)
    return x_t1
end

"""
x[1] is dx/dt (the velocity of the car)
x[2] is x (the position of the car)
x[3] is d phi/dt (the angular velocity of the pendulum)
x[4] is phi (the angle of the pendulum, 0[rad] is at the top)
u    is Vm (the input voltage for the motor)
"""
function inv_pendulum_rk4!(x_t::AbstractVector, u, dt::Number)

    temp = zeros(4)

    #get partial results dy1-dy4
    dy1 = dxdt_invpendulum(x_t, u);
    temp = dy1.*(dt/2.0) .+ x_t;

    dy2 = dxdt_invpendulum(temp, u);
    temp = dy2.*(dt/2.0) .+ x_t;

    dy3 = dxdt_invpendulum(temp, u);
    temp = dy3.*dt .+ x_t;

    dy4 = dxdt_invpendulum(temp, u);

    #update x
    x_t .+= dt.*(dy1.+2.0.*(dy2.+dy3).+dy4)./6.0;

    #x_t[4] = mod(x_t[4]+pi, 2*pi) - pi;
end

"""
x[1] is dx/dt (the velocity of the car)
x[2] is x (the position of the car)
x[3] is d phi/dt (the angular velocity of the pendulum)
x[4] is phi (the angle of the pendulum, 0[rad] is at the top)
u    is Vm (the input voltage for the motor)
"""
function dxdt_invpendulum(x::AbstractVector, u)
    # Parameters
    K_m = 0.00767;  # Gegeninduktionskoeffizient
    K_t = 0.00767;  # Motordrehmomentkonstante
    l_p = 0.3302;   # 0.1778;# Pendellaenge von Drehpunkt bis Schwerpunkt
    M_p = 0.23;     # Masse des Pendels
    M_c = 1.073;    # 0.7031;  # 0.94;# Masse des Wagens
    g = 9.81;
    J_p = 0.0079;    # Traegheitsmoment des Pendels bzgl. seines Schwerpunkts
    B_p = 0.00725;   # 0.0024;  # Viskoser Daempungskoeffizient aus Sciht der Pendelachse
    n_g = 1.0;       # Wirkungsgrad Planetengetriebe
    K_g = 3.7;       # Uerbersetzung Planetengetriebe
    n_m = 1.0;       # Wirkungsgrad des Motors
    R_m = 2.6;       # Motorwiderstand
    r_mp = 6.35e-3;  # Radius des Motorritzels
    B_eq = 5.4;      # viskoser Dämpfungskoeffizient

    F_c = -(n_g*K_g*K_g*n_m*K_t*K_m*x[1])/(R_m*r_mp*r_mp)+(n_g*K_g*n_m*K_t*u[1])/(R_m*r_mp);
    nen = (M_c +M_p)*J_p + M_c*M_p*l_p*l_p + M_p*M_p*l_p*l_p*sin(x[4])*sin(x[4]);

    xdot = zeros(4)
    xdot[1] = (-(J_p + M_p*(l_p*l_p))*B_eq*(x[1]) - ((M_p*M_p)*(l_p*l_p*l_p) + J_p*M_p*l_p)*((x[3])*(x[3]))*sin(x[4])- M_p*l_p*B_p*(x[3])*cos(x[4]) + (J_p + M_p*(l_p*l_p))*F_c + (M_p*M_p)*(l_p*l_p)*g*cos(x[4])*sin(x[4]))/nen;
    xdot[2] = x[1];
    xdot[3] = (-M_p*l_p*B_eq*(x[1])*cos(x[4]) -(M_p*M_p)*(l_p*l_p)*(x[3])*(x[3])*sin(x[4])*cos(x[4]) + (M_c + M_p)*M_p*g*l_p*sin(x[4])-(M_c + M_p)*B_p*(x[3]) + F_c*M_p*l_p*cos(x[4]))/nen;
    xdot[4] = x[3];

    return xdot
end

"""
x[1] is dx/dt (the velocity of the car)
x[2] is x (the position of the car)
x[3] is d phi/dt (the angular velocity of the pendulum)
x[4] is phi (the angle of the pendulum, 0[rad] is at the top)
u    is Vm (the input voltage for the motor)
"""
#Returns the received reward for a certain Environment
function inv_pendulum_cost(x::AbstractVector, u)
    #We want to reward the upright position of the pendulum and a position of x=0
    cost = (x[4]^2 + 0.1*x[3]^2 + 0.001*(u'u) + 0.1*x[2]^2)

    #Use the terminal statement to break the trajectory:
    if x[4]^2 > 2.47
        terminal = 1
        cost *= 100
    else
        terminal = 0
    end

    return cost, terminal
end

#ALTERNATIVE For Balancing:
"""
x[1] is dx/dt (the velocity of the car)
x[2] is x (the position of the car)
x[3] is d phi/dt (the angular velocity of the pendulum)
x[4] is phi (the angle of the pendulum, 0[rad] is at the top)
u    is Vm (the input voltage for the motor)
"""
#Returns the received reward for a certain Environment
function inv_pendulum_cost_swingup(x::AbstractVector, u)

    terminal = 0
    #Give reward according to inverse of pendulum
    cost = -5*cos(x[4]) + 0.2*x[2]^2 - 0.5*x[3]^2 + 0.1*(u'u)

    if abs(x[2]) > 5
        cost += 100
        terminal = 1
    end

    #if not fallen over ;)
    #if x[4]^2 < 2.47
    if cos(x[4]) > 0.5
        #cost = (x[4]^2 + 2*x[3]^2 + 0.01*(u'u) + 0.1*x[2]^2) - 50
        cost = (x[4]^2 + 2*x[3]^2 + 0.01*(u'u) + 0.1*x[2]^2) - 50
    end

    return cost,terminal
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
"""
$(DocStringExtensions.README)
"""
module DiffEqGPU

using DocStringExtensions
using KernelAbstractions, CUDA, DiffEqBase, LinearAlgebra, Distributed
using CUDA: CuPtr, CU_NULL, Mem, CuDefaultStream
using CUDA: CUBLAS
using ForwardDiff
import ZygoteRules
using RecursiveArrayTools
import Base.Threads

@kernel function gpu_kernel(@Const(f),du,@Const(u),@Const(p),@Const(t))
    i = @index(Global, Linear)
    @views @inbounds f(du[:,i],u[:,i],p[:,i],t)
end

@kernel function jac_kernel(@Const(f),J,@Const(u),@Const(p),@Const(t))
    i = @index(Global, Linear)-1
    section = 1 + (i*size(u,1)) : ((i+1)*size(u,1))
    @views @inbounds f(J[section,section],u[:,i+1],p[:,i+1],t)
end

@kernel function discrete_condition_kernel(@Const(condition),cur,@Const(u),@Const(t),@Const(p))
    i = @index(Global, Linear)
    @views @inbounds cur[i] = condition(u[:,i],t,FakeIntegrator(u[:,i],t,p[:,i]))
end

@kernel function discrete_affect!_kernel(@Const(affect!),cur,u,t,p)
    i = @index(Global, Linear)
    @views @inbounds cur[i] && affect!(FakeIntegrator(u[:,i],t,p[:,i]))
end

@kernel function continuous_condition_kernel(@Const(condition),out,@Const(u),@Const(t),@Const(p))
    i = @index(Global, Linear)
    @views @inbounds out[i] = condition(u[:,i],t,FakeIntegrator(u[:,i],t,p[:,i]))
end

@kernel function continuous_affect!_kernel(affect!,event_idx,u,t,p)
    for i in event_idx
        @views @inbounds affect!(FakeIntegrator(u[:,i],t,p[:,i]))
    end
end

maxthreads(::CPU) = 1024
maxthreads(::CUDADevice) = 256

function workgroupsize(backend, n)
    min(maxthreads(backend),n)
end

@kernel function W_kernel(@Const(jac), W, @Const(u), @Const(p), @Const(gamma), @Const(t))
    i = @index(Global, Linear)
    len = size(u,1)
    _W = @inbounds @view(W[:, :, i])
    @views @inbounds jac(_W,u[:,i],p[:,i],t)
    @inbounds for i in eachindex(_W)
        _W[i] = gamma*_W[i]
    end
    _one = one(eltype(_W))
    @inbounds for i in 1:len
        _W[i, i] = _W[i, i] - _one
    end
end

@kernel function Wt_kernel(@Const(jac), W, @Const(u), @Const(p), @Const(gamma), @Const(t))
    i = @index(Global, Linear)
    len = size(u,1)
    _W = @inbounds @view(W[:, :, i])
    @views @inbounds jac(_W,u[:,i],p[:,i],t)
    @inbounds for i in 1:len
        _W[i, i] = -inv(gamma) + _W[i, i]
    end
end

function cuda_lufact!(W)
    CUBLAS.getrf_strided_batched!(W, false)
    return nothing
end

function cpu_lufact!(W)
    len = size(W, 1)
    for i in 1:size(W,3)
        _W = @inbounds @view(W[:, :, i])
        generic_lufact!(_W, len)
    end
    return nothing
end

struct FakeIntegrator{uType,tType,P}
    u::uType
    t::tType
    p::P
end

abstract type EnsembleArrayAlgorithm <: DiffEqBase.EnsembleAlgorithm end
struct EnsembleCPUArray <: EnsembleArrayAlgorithm end
struct EnsembleGPUArray <: EnsembleArrayAlgorithm
    cpu_offload::Float64
end

# Work around the fact that Zygote cannot handle the task system
# Work around the fact that Zygote isderiving fails with constants?
function EnsembleGPUArray()
    EnsembleGPUArray(0.2)
end

ZygoteRules.@adjoint function EnsembleGPUArray()
    EnsembleGPUArray(0.0), _ -> nothing
end

function DiffEqBase.__solve(ensembleprob::DiffEqBase.AbstractEnsembleProblem,
                 alg::Union{DiffEqBase.DEAlgorithm,Nothing},
                 ensemblealg::EnsembleArrayAlgorithm;
                 trajectories, batch_size = trajectories,
                 unstable_check = (dt,u,p,t)->false,
                 kwargs...)

    if trajectories == 1
        return DiffEqBase.__solve(ensembleprob,alg,EnsembleSerial();trajectories=1,kwargs...)
    end

    cpu_trajectories = (ensemblealg isa EnsembleGPUArray && ensembleprob.reduction === DiffEqBase.DEFAULT_REDUCTION) ? round(Int,trajectories * ensemblealg.cpu_offload) : 0
    gpu_trajectories = trajectories - cpu_trajectories

    num_batches = gpu_trajectories ÷ batch_size
    num_batches * batch_size != gpu_trajectories && (num_batches += 1)

    if cpu_trajectories != 0 && ensembleprob.reduction === DiffEqBase.DEFAULT_REDUCTION

        cpu_II = (gpu_trajectories+1):trajectories
        function f()
            solve_batch(ensembleprob,alg,EnsembleThreads(),cpu_II,nothing;kwargs...)
        end

        cpu_sols = Channel{Core.Compiler.return_type(f,Tuple{})}(1)
        t = @task begin
            put!(cpu_sols,f())
        end
        schedule(t)
    end

    if num_batches == 1 && ensembleprob.reduction === DiffEqBase.DEFAULT_REDUCTION
       time = @elapsed sol = batch_solve(ensembleprob,alg,ensemblealg,1:gpu_trajectories;unstable_check=unstable_check,kwargs...)
       if cpu_trajectories != 0
         wait(t)
         sol = vcat(sol,take!(cpu_sols))
       end
       return DiffEqBase.EnsembleSolution(sol,time,true)
    end

    converged::Bool = false
    u = ensembleprob.u_init === nothing ? similar(batch_solve(ensembleprob,alg,ensemblealg,1:batch_size;unstable_check=unstable_check,kwargs...), 0) : ensembleprob.u_init

    if nprocs() == 1
        # While pmap works, this makes much better error messages.
        time = @elapsed begin
            sols = map(1:num_batches) do i
                if i == num_batches
                  I = (batch_size*(i-1)+1):gpu_trajectories
                else
                  I = (batch_size*(i-1)+1):batch_size*i
                end
                batch_data = batch_solve(ensembleprob,alg,ensemblealg,I;unstable_check=unstable_check,kwargs...)
                if ensembleprob.reduction !== DiffEqBase.DEFAULT_REDUCTION
                  u, _ = ensembleprob.reduction(u,batch_data,I)
                  return u
                else
                  batch_data
                end
            end
        end
    else
        time = @elapsed begin
            sols = pmap(1:num_batches) do i
                if i == num_batches
                  I = (batch_size*(i-1)+1):gpu_trajectories
                else
                  I = (batch_size*(i-1)+1):batch_size*i
                end
                x = batch_solve(ensembleprob,alg,ensemblealg,I;unstable_check=unstable_check,kwargs...)
                yield()
                if ensembleprob.reduction !== DiffEqBase.DEFAULT_REDUCTION
                  u, _ = ensembleprob.reduction(u,x,I)
                else
                  x
                end

            end
        end
    end

    if ensembleprob.reduction === DiffEqBase.DEFAULT_REDUCTION
        if cpu_trajectories != 0
          wait(t)
          sols = vcat(reduce(vcat,vec.(sols)),take!(cpu_sols))
        else
          sols = reduce(vcat,sols)
        end
        DiffEqBase.EnsembleSolution(sols,time,true)
    else
        DiffEqBase.EnsembleSolution(sols[end], time, true)
    end
end

diffeqgpunorm(u::AbstractArray,t) = sqrt.(sum(abs2, u)./length(u))
diffeqgpunorm(u::Union{AbstractFloat,Complex},t) = abs(u)
diffeqgpunorm(u::AbstractArray{<:ForwardDiff.Dual},t) = sqrt.(sum(abs2∘ForwardDiff.value, u)./length(u))
diffeqgpunorm(u::ForwardDiff.Dual,t) = abs(ForwardDiff.value(u))

function batch_solve(ensembleprob,alg,ensemblealg::EnsembleArrayAlgorithm,I;kwargs...)
    if ensembleprob.safetycopy
        probs = map(I) do i
            ensembleprob.prob_func(deepcopy(ensembleprob.prob),i,1)
        end
    else
        probs = map(I) do i
            ensembleprob.prob_func(ensembleprob.prob,i,1)
        end
    end
    @assert all(p->p.tspan == probs[1].tspan,probs)
    @assert !isempty(I)
    #@assert all(p->p.f === probs[1].f,probs)

    u0 = reduce(hcat,probs[i].u0 for i in 1:length(I))
    p  = reduce(hcat,probs[i].p  for i in 1:length(I))
    sol, solus = batch_solve_up(ensembleprob,probs,alg,ensemblealg,I,u0,p;kwargs...)
    [ensembleprob.output_func(DiffEqBase.build_solution(probs[i],alg,sol.t,solus[i],destats=sol.destats,retcode=sol.retcode),i)[1] for i in 1:length(probs)]
end

function batch_solve_up(ensembleprob,probs,alg,ensemblealg,I,u0,p;kwargs...)
    if ensemblealg isa EnsembleGPUArray
        u0 = CuArray(u0)
        p  = CuArray(p)
    end

    len = length(probs[1].u0)

    if DiffEqBase.has_jac(probs[1].f)
        jac_prototype = ensemblealg isa EnsembleGPUArray ?
                        cu(zeros(Float32,len,len,length(I))) :
                        zeros(Float32,len,len,length(I))
        if probs[1].f.colorvec !== nothing
            colorvec = repeat(probs[1].f.colorvec,length(I))
        else
            colorvec = repeat(1:length(probs[1].u0),length(I))
        end
    else
        jac_prototype = nothing
        colorvec = nothing
    end

    _callback = generate_callback(probs[1],length(I),ensemblealg; kwargs...)
    prob = generate_problem(probs[1],u0,p,jac_prototype,colorvec)

    if hasproperty(alg, :linsolve)
        _alg = remake(alg,linsolve = LinSolveGPUSplitFactorize())
    else
        _alg = alg
    end

    sol  = solve(prob,_alg; kwargs..., callback = _callback,merge_callbacks = false,
                 internalnorm=diffeqgpunorm)

    us = Array.(sol.u)
    solus = [[@view(us[i][:,j]) for i in 1:length(us)] for j in 1:length(probs)]
    (sol,solus)
end

function seed_duals(x::Matrix{V},::Type{T},
                    ::ForwardDiff.Chunk{N} = ForwardDiff.Chunk(@view(x[:,1]),typemax(Int64)),
                    ) where {V,T,N}
  seeds = ForwardDiff.construct_seeds(ForwardDiff.Partials{N,V})
  duals = [ForwardDiff.Dual{T}(x[i,j],seeds[i]) for i in 1:size(x,1), j in 1:size(x,2)]
end

function extract_dus(us)
  jsize = size(us[1],1), ForwardDiff.npartials(us[1][1])
  utype = typeof(ForwardDiff.value(us[1][1]))
  map(1:size(us[1],2)) do k
      map(us) do u
        du_i = zeros(utype, jsize)
        for i in size(u,1)
          du_i[i, :] = ForwardDiff.partials(u[i,k])
        end
        du_i
      end
  end
end

struct DiffEqGPUAdjTag end

ZygoteRules.@adjoint function batch_solve_up(ensembleprob,probs,alg,ensemblealg,I,u0,p;kwargs...)
    pdual = seed_duals(p,DiffEqGPUAdjTag)
    u0 = convert.(eltype(pdual),u0)

    if ensemblealg isa EnsembleGPUArray
        u0     = CuArray(u0)
        pdual  = CuArray(pdual)
    end

    len = length(probs[1].u0)

    if DiffEqBase.has_jac(probs[1].f)
        jac_prototype = ensemblealg isa EnsembleGPUArray ?
                        cu(zeros(Float32,len,len,length(I))) :
                        zeros(Float32,len,len,length(I))
        if probs[1].f.colorvec !== nothing
            colorvec = repeat(probs[1].f.colorvec,length(I))
        else
            colorvec = repeat(1:length(probs[1].u0),length(I))
        end
    else
        jac_prototype = nothing
        colorvec = nothing
    end

    _callback = generate_callback(probs[1],length(I),ensemblealg)
    prob = generate_problem(probs[1],u0,pdual,jac_prototype,colorvec)

    if hasproperty(alg, :linsolve)
        _alg = remake(alg,linsolve = LinSolveGPUSplitFactorize())
    else
        _alg = alg
    end

    sol  = solve(prob,_alg; kwargs..., callback = _callback,merge_callbacks = false,
                 internalnorm=diffeqgpunorm)

    us = Array.(sol.u)
    solus = [[ForwardDiff.value.(@view(us[i][:,j])) for i in 1:length(us)] for j in 1:length(probs)]

    function batch_solve_up_adjoint(Δ)
        dus = extract_dus(us)
        _Δ = Δ[2]
        adj = map(eachindex(dus)) do j
            sum(eachindex(dus[j])) do i
                J = dus[j][i]
                if _Δ[j] isa AbstractVector
                    v = _Δ[j][i]
                else
                    v = @view _Δ[j][i]
                end
                J'v
            end
        end
        (ntuple(_->nothing, 6)...,Array(VectorOfArray(adj)))
    end
    (sol,solus),batch_solve_up_adjoint
end

function generate_problem(prob::ODEProblem,u0,p,jac_prototype,colorvec)
    _f = let f=prob.f.f
        function (du,u,p,t)
            version = u isa CuArray ? CUDADevice() : CPU()
            wgs = workgroupsize(version,size(u,2))
            wait(version, gpu_kernel(version)(f,du,u,p,t;ndrange=size(u,2),
                                              dependencies=Event(version),
                                              workgroupsize=wgs))
        end
    end

    if DiffEqBase.has_jac(prob.f)
        _Wfact! = let jac=prob.f.jac
            function (W,u,p,gamma,t)
                iscuda = u isa CuArray
                version = iscuda ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, W_kernel(version)(jac, W, u, p, gamma, t;
                                                ndrange=size(u,2),
                                                dependencies=Event(version),
                                                workgroupsize=wgs))
                iscuda ? cuda_lufact!(W) : cpu_lufact!(W)
            end
        end
        _Wfact!_t = let jac=prob.f.jac
            function (W,u,p,gamma,t)
                iscuda = u isa CuArray
                version = iscuda ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, Wt_kernel(version)(jac, W, u, p, gamma, t;
                                                 ndrange=size(u,2),
                                                 dependencies=Event(version),
                                                 workgroupsize=wgs))
                iscuda ? cuda_lufact!(W) : cpu_lufact!(W)
            end
        end
    else
        _Wfact! = nothing
        _Wfact!_t = nothing
    end

    if DiffEqBase.has_tgrad(prob.f)
        _tgrad = let tgrad=prob.f.tgrad
            function (J,u,p,t)
                version = u isa CuArray ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, gpu_kernel(version)(tgrad,J,u,p,t;
                                                  ndrange=size(u,2),
                                                  dependencies=Event(version),
                                                  workgroupsize=wgs))
            end
        end
    else
        _tgrad = nothing
    end

    f_func = ODEFunction(_f,Wfact = _Wfact!,
                        Wfact_t = _Wfact!_t,
                        #colorvec=colorvec,
                        jac_prototype = jac_prototype,
                        tgrad=_tgrad)
    prob = ODEProblem(f_func,u0,prob.tspan,p;
                      prob.kwargs...)
end

function generate_problem(prob::SDEProblem,u0,p,jac_prototype,colorvec)
    _f = let f=prob.f.f
        function (du,u,p,t)
            version = u isa CuArray ? CUDADevice() : CPU()
            wgs = workgroupsize(version,size(u,2))
            wait(version, gpu_kernel(version)(f,du,u,p,t;
                                              ndrange=size(u,2),
                                              dependencies=Event(version),
                                              workgroupsize=wgs))
        end
    end

    _g = let f=prob.f.g
        function (du,u,p,t)
            version = u isa CuArray ? CUDADevice() : CPU()
            wgs = workgroupsize(version,size(u,2))
            wait(version, gpu_kernel(version)(f,du,u,p,t;
                                              ndrange=size(u,2),
                                              dependencies=Event(version),
                                              workgroupsize=wgs))
        end
    end

    if DiffEqBase.has_jac(prob.f)
        _Wfact! = let jac=prob.f.jac
            function (W,u,p,gamma,t)
                iscuda = u isa CuArray
                version = iscuda ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, W_kernel(version)(jac, W, u, p, gamma, t;
                                                ndrange=size(u,2),
                                                dependencies=Event(version),
                                                workgroupsize=wgs))
                iscuda ? cuda_lufact!(W) : cpu_lufact!(W)
            end
        end
        _Wfact!_t = let jac=prob.f.jac
            function (W,u,p,gamma,t)
                iscuda = u isa CuArray
                version = iscuda ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, Wt_kernel(version)(jac, W, u, p, gamma, t;
                                                 ndrange=size(u,2),
                                                 dependencies=Event(version),
                                                 workgroupsize=wgs))
                iscuda ? cuda_lufact!(W) : cpu_lufact!(W)
            end
        end
    else
        _Wfact! = nothing
        _Wfact!_t = nothing
    end

    if DiffEqBase.has_tgrad(prob.f)
        _tgrad = let tgrad=prob.f.tgrad
            function (J,u,p,t)
                version = u isa CuArray ? CUDADevice() : CPU()
                wgs = workgroupsize(version,size(u,2))
                wait(version, gpu_kernel(version)(tgrad,J,u,p,t;
                                                  ndrange=size(u,2),
                                                  dependencies=Event(version),
                                                  workgroupsize=wgs))
            end
        end
    else
        _tgrad = nothing
    end

    f_func = SDEFunction(_f,_g,Wfact = _Wfact!,
                        Wfact_t = _Wfact!_t,
                        #colorvec=colorvec,
                        jac_prototype = jac_prototype,
                        tgrad=_tgrad)
    prob = SDEProblem(f_func,_g,u0,prob.tspan,p;
                      prob.kwargs...)
end

function generate_callback(callback::DiscreteCallback,I,ensemblealg)
    if ensemblealg isa EnsembleGPUArray
        cur = CuArray([false for i in 1:I])
    else
        cur = [false for i in 1:I]
    end
    _condition = callback.condition
    _affect!   = callback.affect!

    condition = function (u,t,integrator)
        version = u isa CuArray ? CUDADevice() : CPU()
        wgs = workgroupsize(version,size(u,2))
        wait(version, discrete_condition_kernel(version)(_condition,cur,u,t,integrator.p;
                                                         ndrange=size(u,2),
                                                         dependencies=Event(version),
                                                         workgroupsize=wgs))
        any(cur)
    end

    affect! = function (integrator)
        version = integrator.u isa CuArray ? CUDADevice() : CPU()
        wgs = workgroupsize(version,size(integrator.u,2))
        wait(version, discrete_affect!_kernel(version)(_affect!,cur,integrator.u,integrator.t,integrator.p;
                                                       ndrange=size(integrator.u,2),
                                                       dependencies=Event(version),
                                                       workgroupsize=wgs))
    end

    return DiscreteCallback(condition,affect!,save_positions=callback.save_positions)
end

function generate_callback(callback::ContinuousCallback,I,ensemblealg)
    _condition   = callback.condition
    _affect!     = callback.affect!
    _affect_neg! = callback.affect_neg!

    condition = function (out,u,t,integrator)
        version = u isa CuArray ? CUDADevice() : CPU()
        wgs = workgroupsize(version,size(u,2))
        wait(version, continuous_condition_kernel(version)(_condition,out,u,t,integrator.p;
                                                           ndrange=size(u,2),
                                                           dependencies=Event(version),
                                                           workgroupsize=wgs))
        nothing
    end

    affect! = function (integrator,event_idx)
        version = integrator.u isa CuArray ? CUDADevice() : CPU()
        wgs = workgroupsize(version,size(integrator.u,2))
        wait(version, continuous_affect!_kernel(version)(_affect!,event_idx,integrator.u,integrator.t,integrator.p;
                                                         ndrange=size(integrator.u,2),
                                                         dependencies=Event(version),
                                                         workgroupsize=wgs))
    end

    affect_neg! = function (integrator,event_idx)
        version = integrator.u isa CuArray ? CUDADevice() : CPU()
        wgs = workgroupsize(version,size(integrator.u,2))
        wait(version, continuous_affect!_kernel(version)(_affect_neg!,event_idx,integrator.u,integrator.t,integrator.p;
                                                         ndrange=size(integrator.u,2),
                                                         dependencies=Event(version),
                                                         workgroupsize=wgs))
    end

    return VectorContinuousCallback(condition,affect!,affect_neg!,I,save_positions=callback.save_positions)
end

function generate_callback(callback::CallbackSet,I,ensemblealg)
    return CallbackSet(map(cb->generate_callback(cb,I,ensemblealg),
            (callback.continuous_callbacks..., callback.discrete_callbacks...))...)
end

generate_callback(::Tuple{},I,ensemblealg) = nothing

function generate_callback(x)
    # will catch any VectorContinuousCallbacks
    error("Callback unsupported")
end

function generate_callback(prob,I,ensemblealg; kwargs...)
    prob_cb = get(prob.kwargs, :callback, ())
    kwarg_cb = get(kwargs, :merge_callbacks, false) ? get(kwargs, :callback, ()) : ()

    if isempty(prob_cb) && isempty(kwarg_cb)
        return nothing
    else
        return CallbackSet(generate_callback(prob_cb,I,ensemblealg),
                           generate_callback(kwarg_cb,I,ensemblealg))
    end
end

### GPU Factorization

struct LinSolveGPUSplitFactorize
    len::Int
    nfacts::Int
end
LinSolveGPUSplitFactorize() = LinSolveGPUSplitFactorize(0, 0)

function (p::LinSolveGPUSplitFactorize)(x,A,b,update_matrix=false;kwargs...)
    version = b isa CuArray ? CUDADevice() : CPU()
    copyto!(x, b)
    wgs = workgroupsize(version,p.nfacts)
    wait(version, ldiv!_kernel(version)(A,x,p.len,p.nfacts;
                                        ndrange=p.nfacts,
                                        dependencies=Event(version),
                                        workgroupsize=wgs))
    return nothing
end

function (p::LinSolveGPUSplitFactorize)(::Type{Val{:init}},f,u0_prototype)
    LinSolveGPUSplitFactorize(size(u0_prototype)...)
end

@kernel function ldiv!_kernel(W,u,@Const(len),@Const(nfacts))
    i = @index(Global, Linear)
    section = 1 + ((i-1)*len) : (i*len)
    _W = @inbounds @view(W[:, :, i])
    _u = @inbounds @view u[section]
    naivesolve!(_W, _u, len)
end

function __printjac(A, ii)
    @cuprintf "[%d, %d]\n" ii.offset1 ii.stride2
    @cuprintf "%d %d %d\n%d %d %d\n%d %d %d\n" ii[1, 1] ii[1, 2] ii[1, 3] ii[2, 1] ii[2, 2] ii[2, 3] ii[3, 1] ii[3, 2] ii[3, 3]
    @cuprintf "%2.2f %2.2f %2.2f\n%2.2f %2.2f %2.2f\n%2.2f %2.2f %2.2f" A[ii[1, 1]] A[ii[1, 2]] A[ii[1, 3]] A[ii[2, 1]] A[ii[2, 2]] A[ii[2, 3]] A[ii[3, 1]] A[ii[3, 2]] A[ii[3, 3]]
end

function generic_lufact!(A::AbstractMatrix{T}, minmn) where {T}
    m = n = minmn
    #@cuprintf "\n\nbefore lufact!\n"
    #__printjac(A, ii)
    #@cuprintf "\n"
    @inbounds for k = 1:minmn
        #@cuprintf "inner factorization loop\n"
        # Scale first column
        Akkinv = inv(A[k,k])
        for i = k+1:m
            #@cuprintf "L\n"
            A[i,k] *= Akkinv
        end
        # Update the rest
        for j = k+1:n, i = k+1:m
            #@cuprintf "U\n"
            A[i,j] -= A[i,k]*A[k,j]
        end
    end
    #@cuprintf "after lufact!"
    #__printjac(A, ii)
    #@cuprintf "\n\n\n"
    return nothing
end

struct MyL{T} # UnitLowerTriangular
    data::T
end
struct MyU{T} # UpperTriangular
    data::T
end

function naivesub!(A::MyU, b::AbstractVector, n)
    x = b
    @inbounds for j in n:-1:1
        xj = x[j] = A.data[j,j] \ b[j]
        for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
            b[i] -= A.data[i,j] * xj
        end
    end
    return nothing
end
function naivesub!(A::MyL, b::AbstractVector, n)
    x = b
    @inbounds for j in 1:n
        xj = x[j]
        for i in j+1:n
            b[i] -= A.data[i,j] * xj
        end
    end
    return nothing
end

function naivesolve!(A::AbstractMatrix, x::AbstractVector, n)
    naivesub!(MyL(A), x, n)
    naivesub!(MyU(A), x, n)
    return nothing
end

function solve_batch(prob,alg,ensemblealg::EnsembleThreads,II,pmap_batch_size;kwargs...)

  if length(II) == 1 || Threads.nthreads() == 1
    return DiffEqBase.solve_batch(prob,alg,EnsembleSerial(),II,pmap_batch_size;kwargs...)
  end

  if typeof(prob.prob) <: DiffEqBase.AbstractJumpProblem && length(II) != 1
    probs = [deepcopy(prob.prob) for i in 1:Threads.nthreads()]
  else
    probs = prob.prob
  end

  #
  batch_size = length(II)÷(Threads.nthreads()-1)

  batch_data = tmap(1:(Threads.nthreads()-1)) do i
    if i == Threads.nthreads()-1
      I_local = II[(batch_size*(i-1)+1):end]
    else
      I_local = II[(batch_size*(i-1)+1):(batch_size*i)]
    end
    DiffEqBase.solve_batch(prob,alg,EnsembleSerial(),I_local,pmap_batch_size;kwargs...)
  end
  DiffEqBase.tighten_container_eltype(batch_data)
end

function tmap(f,args...)
  batch_data = Vector{Core.Compiler.return_type(f,Tuple{typeof.(getindex.(args,1))...})}(undef,length(args[1]))
  Threads.@threads for i in 1:length(args[1])
      batch_data[i] = f(getindex.(args,i)...)
  end
  reduce(vcat,batch_data)
end

export EnsembleCPUArray, EnsembleGPUArray, LinSolveGPUSplitFactorize

end # module


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

using FiniteStateTransducers

export get_L
# builds the L transducer
function get_L(lexicon::Dict{S,Vector{S}}) where {S<:AbstractString}
  phones = sort!(unique!(vcat(values(lexicon)...)))
  words = sort!([keys(lexicon)...])
  psym = Dict(p => i for (i,p) in enumerate(phones))
  wsym = Dict(w => i for (i,w) in enumerate(words ))

  L = WFST(psym,wsym)
  add_states!(L,2)
  initial!(L,1)
  final!(L,2)
  ϵ = get_eps(S)
  c = 3
  for w in keys(lexicon)
    pron=lexicon[w]
    for (i,p) in enumerate(pron)
      if i == 1 && (length(pron) != 1)
        add_arc!(L,1,c,p,w,1)
      elseif i == length(pron)
        if i == 1
          add_arc!(L,1,2,p,w)
        else
          add_arc!(L,c,2,p,ϵ)
          c += 1
        end
      else
        add_arc!(L,c,c+1,p,ϵ)
        c += 1
      end
    end
  end
  add_arc!(L,2,1,ϵ,ϵ)
  return L
end

export get_H
# builds the H transucer, 2 state phone per HMM
function get_H(psym; selfloop_prob=0.4)
  isym = Dict{String,Int}()
  c=1
  for p in sort([keys(psym)...])
    if p == "<SIL>"
      isym["$(p)1"] = c
      isym["$(p)2"] = c+1
      c +=2
    else
      for i=1:2
        isym["$p$i"] = c
        c+=1
      end
    end
  end
  Ns = length(psym)

  H = WFST(isym, psym)
  add_states!(H,Ns+1)
  initial!(H,1)
  ϵ = get_eps(String)
  for p in sort([keys(psym)...])
    if p == "<SIL>"
      # from initial state, assume silence
      add_arc!(H, 1  ,  isym["<SIL>1"]+1, "$(p)1", p)
    end
    # this avoids trivial solution of always staying in the same state
    add_arc!(H, isym["$(p)1"]+1, isym["$(p)1"]+1, "$(p)1", ϵ,-log(selfloop_prob))
    # prob of transistion to other state unknown, set to 1
    add_arc!(H, isym["$(p)1"]+1, isym["$(p)2"]+1, "$(p)2", ϵ)
    # to final state
    final!(H,isym["$(p)2"]+1)
  end
  for s in keys(get_final(H))
    for p in keys(psym)
      # prob of transistion to other phone unknown, set to 1
      add_arc!(H, s, isym["$(p)1"]+1, "$(p)1", p) # emitting state
    end
  end
  return H
end

export Hfst2trans
# convert the H transducer into transition matrix
function Hfst2trans(H::WFST)
  Ns = length(get_isym(H))
  A = zeros(Float32,Ns,Ns)
  state2outtr=Dict(i => (get_ilabel.(s),get_weight.(s)) for (i,s) in enumerate(H))
  for (p,s,n,d,e,a) in FiniteStateTransducers.DFS(H,1)
    if d
      intr = get_ilabel(a)
      outtr,w = state2outtr[n]
      for i in eachindex(outtr)
        A[intr,outtr[i]] = exp(-get(w[i]))
      end
    end
  end
  return A
end

export get_lexicon
function get_lexicon()
  lexicon = Dict(
                 "<SIL>" => ["<SIL>"],
                 "OH"  => ["OW"],
                 "ZERO" =>  ["Z", "IH", "R", "OW"],
                 "ONE" =>  ["W", "AH", "N"],
                 "TWO" =>  ["T", "UW"],
                 "THREE" => ["TH", "R", "IY"],
                 "FOUR"  => ["F", "AO", "R"],
                 "FIVE" =>  ["F", "AY", "V"],
                 "SIX" => ["S", "IH", "KS"],
                 "SEVEN"  => ["S", "EH", "V", "AH", "N"],
                 "EIGHT" => ["EY", "T"],
                 "NINE" =>  ["N", "AY", "N"]
                )
  ilexicon = Dict(lexicon[w] => w for w in keys(lexicon))
  return lexicon, ilexicon
end

export get_HL
function get_HL(lexicon)
  L = get_L(lexicon)
  H = get_H(get_isym(L))
  return H,L
end

export get_aA
function get_aA(H; use_log=true)
  A = Hfst2trans(H)
  Ns = size(A,1)
  a = zeros(Float32,Ns) # initial state prob
  a[H.isym["<SIL>1"]] = one(Float32)

  if use_log
    A .= log.(A)
    a .= log.(a)
  end
  return a,A
end


text:
abstract type AbstractParametersVPRK{DT,TT,ET,D,S} <: Parameters{DT,TT} end
abstract type AbstractIntegratorVPRK{DT,TT} <: DeterministicIntegrator{DT,TT} end
abstract type AbstractIntegratorVPRKwProjection{DT,TT} <: AbstractIntegratorVPRK{DT,TT} end

@inline equation(integrator::AbstractIntegratorVPRK) = integrator.params.equ
@inline timestep(integrator::AbstractIntegratorVPRK) = integrator.params.Δt
@inline tableau(integrator::AbstractIntegratorVPRK) = integrator.params.tab

@inline nstages(integrator::AbstractIntegratorVPRK)  = nstages(tableau(integrator))
@inline eachstage(integrator::AbstractIntegratorVPRK) = 1:nstages(integrator)

function update_params!(params::AbstractParametersVPRK, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t̅  = sol.t
    params.q̅ .= sol.q
    params.p̅ .= sol.p
end

function update_solution!(int::AbstractIntegratorVPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))
end

function project_solution!(int::AbstractIntegratorVPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, R::Vector{TT}) where {DT,TT}
    update_solution!(sol.q, sol.q̃, int.cache.U, R, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.G, R, timestep(int))
end

function project_solution!(int::AbstractIntegratorVPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, RU::Vector{TT}, RG::Vector{TT}) where {DT,TT}
    update_solution!(sol.q, sol.q̃, int.cache.U, RU, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.G, RG, timestep(int))
end

function initialize!(int::AbstractIntegratorVPRK{DT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end


text:
using BenchmarkTools
SUITE = BenchmarkGroup()
for file in readdir(@__DIR__)
    if startswith(file, "bench_") && endswith(file, ".jl")
        SUITE[file[length("bench_") + 1:end - length(".jl")]] =
            include(file)
    end
end


text:
module Chessboard

using Mapper2

# For testing, we create a "chessboard" type architecture. The layout of the
# architecture will either be rectangular 2D / 3D, or hexagonal 2D. About half
# of the tiles will be "black_squares", while the other half will be "white_squares"
#
# We will have 3 classes of tasks:
# "black" tasks - can only be mapped to black squares
# "white" tasks - can only be mapped to white squares
# "gray" tasks - can be mapped to both squares.
#
# 1/4 of tasks will be black, 1/4 will be white, and the remaining tasks will
# be gray.
#
# This will help test various neighbor and move generation capabilities of the
# Mapper.
@enum SquareColor White Black Gray

struct Chess <: RuleSet end

export Chess,
    # Colors
    SquareColor,
    White,
    Black,
    Gray,
    # Architecture
    architecture,
    Rectangle2D,
    Rectangle3D,
    Hexagonal2D,
    ChessboardColor,
    HashColor,
    # Taskgraph
    AllGray,
    OddEven,
    Quarters,
    taskgraph,
    linegraph

include("Taskgraph.jl")
include("Architecture.jl")

# Extend "ismappable"
Mapper2.ismappable(::Chess, x::Component) = haskey(x.metadata, "color")

# Two tasks are equivalent if they have the same color.
function Mapper2.isequivalent(::Chess, a::TaskgraphNode, b::TaskgraphNode)
    return a.metadata["color"] == b.metadata["color"]
end

function Mapper2.canmap(::Chess, t::TaskgraphNode, c::Component)
    haskey(c.metadata, "color") || return false

    task_color = t.metadata["color"]
    component_color = c.metadata["color"]

    # Match by color.
    return task_color == Gray || (task_color == component_color)
end

Mapper2.getcapacity(::Chess, args...) = 5

################################################################################
# Map
################################################################################
function build_map(dim, ntasks, nedges)
    arch = architecture(dim, Val{2}(), ChessboardColor())
    taskgraph = build_taskgraph(ntasks, nedges)
    return Map(Chess(), arch, taskgraph)
end

end


text:
import Distributions: MvNormal, pdf, logpdf
import Plots: scatter, contour, contour!, xlabel!, ylabel!
using ColorSchemes
using LaTeXStrings

function prob_cluster(x, y, μ, Σ)
    "Probability under Gaussian"

    # Ensure matrix is Hermitian
    Σ = enforce_Hermitian!(Σ)

    # Ensure matrix is non-negative
    Σ = enforce_nonneg!(Σ)

    # Return weighted probability
    return pdf(MvNormal(μ, Σ), [x, y])
end

function enforce_nonneg!(M)
    "Enforce non-negative elements"
    M = max.(0, M)
end

function enforce_Hermitian!(M)
    "Enforce Hermitian"
    M = (M + M')/2
end

function plot_clusters(observations;
                       μ=[],
                       Σ=[[1. 0.;0. 1.], [1. 0.;0. 1.]],
                       x1=range(-5, stop=3),
                       x2=range(-5, stop=4),
                       colorlist=["reds", "blues", "greens"],
                       size=(500, 500))
    "Scatter data and overlay Gaussian clusters."

    # Show the data
    plt = scatter(observations[:,1],
                  observations[:,2],
                  label="",
                  color="black",
                  xlims=[x1[1] x1[end]],
                  ylims=[x2[1] x2[end]],
                  xlabel=L"$x_1$",
                  ylabel=L"$x_2$",
                  grid=false,
                  size=size)

    # Loop over clusters
    for k = 1:length(μ)

        # Overlay the posterior probabilities for k-th cluster
        contour!(plt, x1, x2,
                 (x, y) -> prob_cluster(x, y, μ[k], Σ[k]),
                 color=colorlist[k],
                 cbar=nothing)
    end
    plot(plt)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
__precompile__()

module Gencode

using SymPy

#export syminit, gencode, compilecode
export syminit, gencode, gencodeall, compilecode, runcode, checkcompilers, setcompilers, tring2cmd

include("syminit.jl");
include("contains.jl");
include("varsassign.jl");
include("sympyassign.jl");
include("gencodebou.jl");
include("gencodeelem.jl");
include("gencodeelem2.jl");
include("gencodeelem3.jl");
include("gencodeface.jl");
include("gencodeface2.jl");
include("nocodeelem.jl");
include("nocodeelem2.jl");
include("nocodeelem3.jl");
include("nocodeface.jl");
include("nocodeface2.jl");
include("gencodeelemface.jl");
include("gencodeall.jl");
include("string2cmd.jl");
include("genlib.jl");
include("checkcompilers.jl");
include("setcompilers.jl");
include("compilecode.jl");
include("runcode.jl");

function gencode(app)

print("generate code...\n");
if !isdir("app")
    mkdir("app");
else
    if isfile("app/opuApp.a")
        rm("app/opuApp.a")
    end
    if isfile("app/cpuApp.a")
        rm("app/cpuApp.a")
    end
    if isfile("app/gpuApp.a")
        rm("app/gpuApp.a")
    end
end

xdg, udg, udg1, udg2, wdg, wdg1, wdg2, odg, odg1, odg2, uhg, nlg, tau, uinf, param, time = syminit(app);

if app.modelnumber==0
    strn = "";
else
    strn = string(app.modelnumber);
end

ncu = app.ncu;
u = udg[1:ncu];
u1 = udg1[1:ncu];
u2 = udg2[1:ncu];
if app.nc>app.ncu
    q = udg[ncu+1:end];
    q1 = udg1[ncu+1:end];
    q2 = udg2[ncu+1:end];
else
    q = [];
    q1 = [];
    q2 = [];
end

if isdefined(Main, Symbol(app.modelfile))
    pdemodel = getfield(Main, Symbol(app.modelfile))
else
    pdemodel = getfield(Main, Symbol("Main"))
end

if isdefined(pdemodel, Symbol("flux"))
    #f = pdemodel.flux(xdg, udg, odg, wdg, uinf, param, time);
    f = pdemodel.flux(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem("Flux" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    error("app.Flux is empty");
end
if isdefined(pdemodel, Symbol("source"))
    #f = pdemodel.source(xdg, udg, odg, wdg, uinf, param, time);
    f = pdemodel.source(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem("Source" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    nocodeelem("Source" * strn);
end
if isdefined(pdemodel, Symbol("sourcew"))
    f = pdemodel.sourcew(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem2("Sourcew" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    nocodeelem2("Sourcew" * strn);
end
if isdefined(pdemodel, Symbol("mass"))
    #f = pdemodel.mass(xdg, udg, odg, wdg, uinf, param, time);
    f = pdemodel.mass(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem("Tdfunc" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    if app.model=="ModelW" || app.model == "modelW" || app.tdep==1
        error("pde.mass is not defined");
    else
        nocodeelem("Tdfunc" * strn);
    end
end
if isdefined(pdemodel, Symbol("avifield"))
    #f = pdemodel.avfield(xdg, udg, odg, wdg, uinf, param, time);
    f = pdemodel.avfield(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem2("Avfield" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    nocodeelem2("Avfield" * strn);
end
if isdefined(pdemodel, Symbol("output"))
    #f = pdemodel.output(xdg, udg, odg, wdg, uinf, param, time);
    f = pdemodel.output(u, q, wdg, odg, xdg, time, param, uinf);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeelem2("Output" * strn, f, xdg, udg, odg, wdg, uinf, param, time);
else
    nocodeelem2("Output" * strn);
end
if isdefined(pdemodel, Symbol("fbou"))
    #f = pdemodel.fbou(xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time);
    f = pdemodel.fbou(u, q, wdg, odg, xdg, time, param, uinf, uhg, nlg, tau);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    f = reshape(f,app.ncu,Int(length(f)/app.ncu));
    gencodeface("Fbou" * strn, f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time);
else
    error("app.Fbou is empty");
end
if isdefined(pdemodel, Symbol("ubou"))
    #f = pdemodel.ubou(xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time);
    f = pdemodel.ubou(u, q, wdg, odg, xdg, time, param, uinf, uhg, nlg, tau);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    f = reshape(f,app.ncu,Int(length(f)/app.ncu));
    gencodeface("Ubou" * strn, f, xdg, udg, odg, wdg, uhg, nlg, tau, uinf, param, time);
else
    nocodeface("Ubou" * strn);
end
if isdefined(pdemodel, Symbol("Fhat"))
    #f = pdemodel.fhat(xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
    f = pdemodel.fhat(u1, q1, wdg1, odg1, xdg, time, param, uinf, uhg, nlg, tau, u2, q2, wdg2, odg2);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeface2("Fhat" * strn, f, xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
else
    nocodeface2("Fhat" * strn);
end
if isdefined(pdemodel, Symbol("uhat"))
    #f = pdemodel.uhat(xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
    f = pdemodel.uhat(u1, q1, wdg1, odg1, xdg, time, param, uinf, uhg, nlg, tau, u2, q2, wdg2, odg2);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeface2("Uhat" * strn, f, xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
else
    nocodeface2("Uhat" * strn);
end
if isdefined(pdemodel, Symbol("stab"))
    #f = pdemodel.stab(xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
    f = pdemodel.stab(u1, q1, wdg1, odg1, xdg, time, param, uinf, uhg, nlg, tau, u2, q2, wdg2, odg2);
    if length(f)==1
        f = reshape([f],1,1);
    end
    f = f[:];
    gencodeface2("Stab" * strn, f, xdg, udg1, udg2, odg1, odg2, wdg1, wdg2, uhg, nlg, tau, uinf, param, time);
else
    nocodeface2("Stab" * strn);
end
if isdefined(pdemodel, Symbol("initu"))
    udg = pdemodel.initu(xdg, param, uinf);
    if length(udg)==1
        udg = reshape([udg],1,1);
    end
    udg = udg[:];
    gencodeelem3("Initu" * strn, udg, xdg, uinf, param);
else
    error("initu is not defined");
end
if isdefined(pdemodel, Symbol("initw"))
    wdg = pdemodel.initw(xdg, param, uinf);
    if length(wdg)==1
        wdg = reshape([wdg],1,1);
    end
    wdg = wdg[:];
    gencodeelem3("Initwdg" * strn, wdg, xdg, uinf, param);
else
    nocodeelem3("Initwdg" * strn);
end
if isdefined(pdemodel, Symbol("initv"))
    odg = pdemodel.initv(xdg, param, uinf);
    if length(odg)==1
        odg = reshape([odg],1,1);
    end
    odg = odg[:];
    gencodeelem3("Initodg" * strn, odg, xdg, uinf, param);
else
    nocodeelem3("Initodg" * strn);
end
if isdefined(pdemodel, Symbol("initq"))
    qdg = pdemodel.initq(xdg, param, uinf);
    if length(qdg)==1
        qdg = reshape([qdg],1,1);
    end
    qdg = qdg[:];
    gencodeelem3("Initq" * strn, qdg, xdg, uinf, param);

    udg = pdemodel.initu(xdg, param, uinf);
    if length(udg)==1
        udg = reshape([udg],1,1);
    end
    udg = udg[:];

    udg = [udg; qdg];
    gencodeelem3("Initudg" * strn, udg, xdg, uinf, param);
else
    if app.model == "ModelW" || app.model == "modelW"
        error("initq is not defined");
    else
        nocodeelem3("Initq" * strn);
        nocodeelem3("Initudg" * strn);
    end
end

# if isdefined(pdemodel, Symbol("initq"))
#     qdg = pdemodel.initq(xdg, uinf, param);
#     if length(qdg)==1
#         qdg = reshape([qdg],1,1);
#     end
#     qdg = qdg[:];
#     gencodeelem3("Initq", qdg, xdg, uinf, param);
# else
#     nocodeelem3("Initq");
# end
# if isdefined(pdemodel, Symbol("inituq"))
#     udg = pdemodel.inituq(xdg, uinf, param);
#     if length(udg)==1
#         udg = reshape([udg],1,1);
#     end
#     udg = udg[:];
#     gencodeelem3("Initudg", udg, xdg, uinf, param);
# else
#     if app.model == "ModelW" || app.model == "modelW"
#         error("inituq is not defined");
#     else
#         nocodeelem3("Initudg");
#     end
# end

end

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
#Came from GLAbstraction/GLUtils.jl
#question: Used?
function close_to_square(n::Real)
    # a cannot be greater than the square root of n
    # b cannot be smaller than the square root of n
    # we get the maximum allowed value of a
    amax = floor(Int, sqrt(n))
    if 0 == rem(n, amax)
        # special case where n is a square number
        return (amax, div(n, amax))
    end
    # Get its prime factors of n
    primeFactors = factor(n)
    # Start with a factor 1 in the list of candidates for a
    candidates = Int[1]
    for (f, _) in primeFactors
        # Add new candidates which are obtained by multiplying
        # existing candidates with the new prime factor f
        # Set union ensures that duplicate candidates are removed
        candidates = union(candidates, f .* candidates)
        # throw out candidates which are larger than amax
        filter!(x -> x <= amax, candidates)
    end
    # Take the largest factor in the list d
    return (candidates[end], div(n, candidates[end]))
end

#Came from GLAbstraction/GLCamera.jl
w_component(::Point{N,T}) where {N,T} = T(1)
w_component(::Vec{N,T}) where {N,T} = T(0)

# @component struct AABB
#     origin  ::Vec3f0
#     diagonal::Vec3f0
# end

# function aabb_ray_intersect(aabb::AABB, entity_pos::Point3f0, ray_origin::Point3f0, ray_direction::Vec3f0)
#     dirfrac = 1.0f0 ./ ray_direction
#     right   = aabb.origin + aabb.diagonal

#     real_origin = aabb.origin + entity_pos
#     real_right  = right + entity_pos
#     t1 = (real_origin - ray_origin) .* dirfrac
#     t2 = (real_right - ray_origin)  .* dirfrac
#     tsmaller = min.(t1, t2)
#     tbigger  = max.(t1,t2)

#     tmin = maximum(tsmaller)
#     tmax = minimum(tbigger)

#     return tmin < tmax
# end


text:
using ..ClangCompiler: libclangex
using LLVM.API: LLVMModuleRef, LLVMOpaqueModule
using LLVM.API: LLVMOpaqueContext, LLVMContextRef
using LLVM.API: LLVMMemoryBufferRef, LLVMGenericValueRef

const time_t = Clong


text:
using DiffEqPDEBase
using Base.Test


text:
"""
    AnatomicalStructureList
"""
struct AnatomicalVector{T<:AbstractAnatomicalStructure} <: AbstractVector{T}
    values::Vector{T}
end

Base.getindex(av::AnatomicalVector, i::Int) = av.values[i]

Base.length(av::AnatomicalVector) = length(av.values)

Base.size(av::AnatomicalVector) = size(av.values)

Base.names(av::AnatomicalVector) = names.(av.values)

function splitbilaterals(av::AnatomicalVector{T}) where {T}
    out = T[]
    for a_i in av
        if isbilateral(a_i)
            push!(out, left(a_i))
            push!(out, right(a_i))
        else
            push!(out, a_i)
        end
    end
    return AnatomicalVector(out)
end

function Base.getindex(a::AnatomicalVector, i::AbstractAnatomicalStructure)
    getindex(a, Base.to_index(a, i))
end

function Base.to_index(a::AnatomicalVector, i::AbstractAnatomicalStructure)
    searchsortedfirst(names(a), names(i))
end

function Base.show(io::IO, a::AnatomicalVector)
    print(io, summary(a))
    for a_i in a
        print(io, summary(a_i))
    end
end


text:
__precompile__(true)
module Visualize

using Compat, FileIO, FieldTraits

using FieldTraits: @reactivecomposed, @field, @composed, @needs
using FieldTraits: Composable, ComposableLike, Field, UsageError, Fields, Links
using FieldTraits: ReactiveComposable, Partial
import FieldTraits: on, default, convertfor

using Colors, ColorVectorSpace, StaticArrays
using GeometryTypes, Quaternions
using Interpolations

import GLAbstraction, ColorVectorSpace
import Transpiler: gli
import Transpiler: mix, smoothstep
import Transpiler.gli: dFdx, dFdy

import Base: scale!

"""
Replacement of Pkg.dir("Visualize") --> Visualize.dir,
returning the correct path
"""
dir(dirs...) = joinpath(normpath(dirname(@__FILE__), ".."), dirs...)

"""
returns path relative to the assets folder
"""
assetpath(folders...) = dir("assets", folders...)

"""
Loads a file from the asset folder
"""
function loadasset(folders...; kw_args...)
    path = assetpath(folders...)
    isfile(path) || isdir(path) || error("Could not locate file at $path")
    load(path; kw_args...)
end

export assetpath, loadasset

include("math.jl")
include("utils.jl")
include("base.jl")
include("events.jl")

include("camera.jl")
include("windowbase.jl")

include("shader/base.jl")
include("images.jl")
include("text/textbase.jl")

include("backends/julia/base.jl")
include("backends/cairo/base.jl")
include("backends/webgl/base.jl")
include("backends/opengl/base.jl")
include("visualize.jl")

using .GLRasterization
using .JLRasterization

export orthographicprojection, perspectiveprojection, lookat
export normalmesh, uvmesh, JLCanvas, JLRasterizer, Area, Framebuffer
export add!, TranslationSpeed, LookAt, EyePosition, ProjectionView, EyePosition
export RotationSpeed, Translation, Rotation, Keyboard, WindowEvents, Window
export Mouse, Pan, View, Projection, IRect, Sampler, isopen
export Resolution, Open, Visible, renderloop, Light, Shading
export glwindow, ShadingFunction, Mesh, Vertices, visualize
export Model, Color, Position, LineVertex, LineSegments, Camera, BasicCamera, PerspectiveCamera

end # module


text:
using FEASTSolver
using LinearAlgebra
using DelimitedFiles
using SparseArrays

A = spdiagm(-1 => fill(-1.0, 99), 0 => fill(2.0, 100), 1 => fill(-1.0, 99))
# A = diagm(-1 => fill(-1.0, 99), 0 => fill(2.0, 100), 1 => fill(-1.0, 99))
B = sprand(100, 100, 0.2)
B = B + B'
C, R = complex(0.05, 0.0), 0.05 # contour
# C, R = complex(0.1, 0.0), 0.1 # contour
# X = rand(ComplexF64, 100, 20) # initial subspace
X = sprand(ComplexF64, 100, 20, 0.9)

contour = circular_contour_trapezoidal(C, R, 8)
# contour = rectangular_contour_gauss(0.0 - R*im, 2*R + R*im, 8)
# contour = rectangular_contour_trapezoidal(0.0 - R*im, 2*R + R*im, 32)
# contour = rectangular_contour_gauss(0.9 - R*im, 1.1 + R*im, 16)

# z = 0.025+0.0im
# println(rational_func(z, contour))

e, v, res = gen_feast!(X, A, B, contour, debug=true, ϵ=10e-15)
# @timev e, v, res = feast!(X, A, contour, debug=false, ϵ=10e-15, iter=1)

convergence_info(e,v,res,contour)
display(e)


text:
using QXTns
using QXTools
using QXContexts

using YaoBlocks
using Reexport

@reexport using YaoAPI

export QXReg, ghz, qx_intrinsic_nodes

"""Register for QX"""
mutable struct QXReg{B} <: AbstractRegister{B}
    qubits::Int
    tnc::TensorNetworkCircuit
    cg::ComputeGraph
    QXReg(n, tnc) = new{1}(n, tnc)
end

function QXReg(n::Int)
    tnc = TensorNetworkCircuit(n)
    QXReg(n, tnc)
end

"""
    YaoAPI.apply!(reg::QXReg, block::AbstractBlock)

Apply the given block to the register
"""
function YaoAPI.apply!(reg::QXReg, block::AbstractBlock)
    for (locs, matrix) in qx_intrinsic_nodes(block)
        push!(reg.tnc, collect(locs), Matrix(matrix))
    end
    return reg
end

YaoAPI.nactive(reg::QXReg) = reg.qubits
YaoAPI.nqubits(reg::QXReg) = reg.qubits

"""Prepare a ghz circuit block"""
function ghz(n)
    chain(put(1=>H), chain(map(x -> control(n, x, x+1=>X), 1:n-1)...))
end

qx_intrinsic_nodes(blk::AbstractBlock) = qx_intrinsic_nodes!([], blk)

function qx_intrinsic_nodes!(list, blk::AbstractBlock)
    for each in subblocks(blk)
        if each isa ControlBlock && length(each.ctrl_locs) < 2 && length(each.locs) == 1
            if each.ctrl_locs[1] < each.locs[1]
                gate = control(2, 1, 2=>each.content)
            else
                gate = control(2, 2, 1=>each.content)
            end

            push!(list, occupied_locs(each) => mat(gate))
        elseif each isa PutBlock
            push!(list, occupied_locs(each)=>mat(first(subblocks(each))))
        elseif each isa PrimitiveBlock
            push!(list, 1=>mat(each))
        else
            qx_intrinsic_nodes!(list, each)
        end
    end
    return list
end

find_primitive_blocks(blk::AbstractBlock) = find_primitive_blocks!([], blk)

function find_primitive_blocks!(list, blk::AbstractBlock)
    for each in subblocks(blk)
        if each isa PrimitiveBlock
            push!(list, each)
        else
            find_primitive_blocks!(list, each)
        end
    end
    return list
end

"""
    YaoAPI.measure(reg::QXReg; nshots=10)

Measure the register and return the measurements. Note that currently this returns
an array of bitstrings and an array of corresponsing amplitudes. This should be replaced
by a rejection sampler which should be able to return measurements distributed according to
the expected output distribution.
"""
function YaoAPI.measure(reg::QXReg; nshots=10)
    if !isdefined(reg, :cg)
        add_input!(reg.tnc, "0"^nqubits(reg))
        add_output!(reg.tnc, "0"^nqubits(reg))
        # must be outputs on network currently to find a plan
        bond_groups, plan, _ = contraction_scheme(reg.tnc.tn, 2;
                                                  seed=rand(Int),
                                                  time=10)

        reg.cg = build_compute_graph(reg.tnc, plan, bond_groups)
    end
    ctx = QXContext(reg.cg)
    sampler_params = Dict(:method => "List",
                          :params =>
                     Dict(:num_samples => nshots,
                          :bitstrings => collect(amplitudes_uniform(nqubits(reg), nothing, nshots))))
    sampler = ListSampler(ctx; sampler_params[:params]...)
    sampler()
end

"""
    QXTools._convert_to_tnc(block::AbstractBlock; kwargs...)

Convert the given YaoBlocks circuit to a TensorNetworkCircuit
"""
function QXTools._convert_to_tnc(block::AbstractBlock; kwargs...)
    tnc = TensorNetworkCircuit(nqubits(block))
    for (locs, matrix) in qx_intrinsic_nodes(block)
        push!(tnc, collect(locs), Matrix(matrix); kwargs...)
    end
    return tnc
end

"""
    QXTools.convert_to_tnc(circ::AbstractBlock;
                           input::Union{String, Nothing}=nothing,
                           output::Union{String, Nothing}=nothing,
                           no_input::Bool=false,
                           no_output::Bool=false,
                           kwargs...)

Function to convert a YaoBlocks circuit to a QXTools tensor network circuit
"""
function QXTools.convert_to_tnc(circ::AbstractBlock;
                        input::Union{String, Nothing}=nothing,
                        output::Union{String, Nothing}=nothing,
                        no_input::Bool=false,
                        no_output::Bool=false,
                        kwargs...)
    tnc = QXTools._convert_to_tnc(circ; kwargs...)
    if !no_input add_input!(tnc, input) end
    if !no_output add_output!(tnc, output) end
    tnc
end