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

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

Cancels an Amazon Braket job.

# Arguments
- `job_arn`: The ARN of the Amazon Braket job to cancel.

"""
function cancel_job(jobArn; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "PUT",
        "/job/$(jobArn)/cancel";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function cancel_job(
    jobArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "PUT",
        "/job/$(jobArn)/cancel",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    cancel_quantum_task(client_token, quantum_task_arn)
    cancel_quantum_task(client_token, quantum_task_arn, params::Dict{String,<:Any})

Cancels the specified task.

# Arguments
- `client_token`: The client token associated with the request.
- `quantum_task_arn`: The ARN of the task to cancel.

"""
function cancel_quantum_task(
    clientToken, quantumTaskArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "PUT",
        "/quantum-task/$(quantumTaskArn)/cancel",
        Dict{String,Any}("clientToken" => clientToken);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function cancel_quantum_task(
    clientToken,
    quantumTaskArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "PUT",
        "/quantum-task/$(quantumTaskArn)/cancel",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => clientToken), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_job(algorithm_specification, client_token, device_config, instance_config, job_name, output_data_config, role_arn)
    create_job(algorithm_specification, client_token, device_config, instance_config, job_name, output_data_config, role_arn, params::Dict{String,<:Any})

Creates an Amazon Braket job.

# Arguments
- `algorithm_specification`: Definition of the Amazon Braket job to be created. Specifies
  the container image the job uses and information about the Python scripts used for entry
  and training.
- `client_token`: A unique token that guarantees that the call to this API is idempotent.
- `device_config`: The quantum processing unit (QPU) or simulator used to create an Amazon
  Braket job.
- `instance_config`: Configuration of the resource instances to use while running the
  hybrid job on Amazon Braket.
- `job_name`: The name of the Amazon Braket job.
- `output_data_config`: The path to the S3 location where you want to store job artifacts
  and the encryption key used to store them.
- `role_arn`: The Amazon Resource Name (ARN) of an IAM role that Amazon Braket can assume
  to perform tasks on behalf of a user. It can access user resources, run an Amazon Braket
  job container on behalf of user, and output resources to the users' s3 buckets.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"checkpointConfig"`: Information about the output locations for job checkpoint data.
- `"hyperParameters"`: Algorithm-specific parameters used by an Amazon Braket job that
  influence the quality of the training job. The values are set with a string of JSON
  key:value pairs, where the key is the name of the hyperparameter and the value is the value
  of th hyperparameter.
- `"inputDataConfig"`: A list of parameters that specify the name and type of input data
  and where it is located.
- `"stoppingCondition"`:  The user-defined criteria that specifies when a job stops running.
- `"tags"`: A tag object that consists of a key and an optional value, used to manage
  metadata for Amazon Braket resources.
"""
function create_job(
    algorithmSpecification,
    clientToken,
    deviceConfig,
    instanceConfig,
    jobName,
    outputDataConfig,
    roleArn;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "POST",
        "/job",
        Dict{String,Any}(
            "algorithmSpecification" => algorithmSpecification,
            "clientToken" => clientToken,
            "deviceConfig" => deviceConfig,
            "instanceConfig" => instanceConfig,
            "jobName" => jobName,
            "outputDataConfig" => outputDataConfig,
            "roleArn" => roleArn,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_job(
    algorithmSpecification,
    clientToken,
    deviceConfig,
    instanceConfig,
    jobName,
    outputDataConfig,
    roleArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "POST",
        "/job",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "algorithmSpecification" => algorithmSpecification,
                    "clientToken" => clientToken,
                    "deviceConfig" => deviceConfig,
                    "instanceConfig" => instanceConfig,
                    "jobName" => jobName,
                    "outputDataConfig" => outputDataConfig,
                    "roleArn" => roleArn,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_quantum_task(action, client_token, device_arn, output_s3_bucket, output_s3_key_prefix, shots)
    create_quantum_task(action, client_token, device_arn, output_s3_bucket, output_s3_key_prefix, shots, params::Dict{String,<:Any})

Creates a quantum task.

# Arguments
- `action`: The action associated with the task.
- `client_token`: The client token associated with the request.
- `device_arn`: The ARN of the device to run the task on.
- `output_s3_bucket`: The S3 bucket to store task result files in.
- `output_s3_key_prefix`: The key prefix for the location in the S3 bucket to store task
  results in.
- `shots`: The number of shots to use for the task.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"deviceParameters"`: The parameters for the device to run the task on.
- `"jobToken"`: The token for an Amazon Braket job that associates it with the quantum task.
- `"tags"`: Tags to be added to the quantum task you're creating.
"""
function create_quantum_task(
    action,
    clientToken,
    deviceArn,
    outputS3Bucket,
    outputS3KeyPrefix,
    shots;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "POST",
        "/quantum-task",
        Dict{String,Any}(
            "action" => action,
            "clientToken" => clientToken,
            "deviceArn" => deviceArn,
            "outputS3Bucket" => outputS3Bucket,
            "outputS3KeyPrefix" => outputS3KeyPrefix,
            "shots" => shots,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_quantum_task(
    action,
    clientToken,
    deviceArn,
    outputS3Bucket,
    outputS3KeyPrefix,
    shots,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "POST",
        "/quantum-task",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "action" => action,
                    "clientToken" => clientToken,
                    "deviceArn" => deviceArn,
                    "outputS3Bucket" => outputS3Bucket,
                    "outputS3KeyPrefix" => outputS3KeyPrefix,
                    "shots" => shots,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the devices available in Amazon Braket.

# Arguments
- `device_arn`: The ARN of the device to retrieve.

"""
function get_device(deviceArn; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "GET",
        "/device/$(deviceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_device(
    deviceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "GET",
        "/device/$(deviceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the specified Amazon Braket job.

# Arguments
- `job_arn`: The ARN of the job to retrieve.

"""
function get_job(jobArn; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "GET", "/job/$(jobArn)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_job(
    jobArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "GET",
        "/job/$(jobArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Retrieves the specified quantum task.

# Arguments
- `quantum_task_arn`: the ARN of the task to retrieve.

"""
function get_quantum_task(quantumTaskArn; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "GET",
        "/quantum-task/$(quantumTaskArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_quantum_task(
    quantumTaskArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "GET",
        "/quantum-task/$(quantumTaskArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Shows the tags associated with this resource.

# Arguments
- `resource_arn`: Specify the resourceArn for the resource whose tags to display.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Searches for devices using the specified filters.

# Arguments
- `filters`: The filter values to use to search for a device.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned in the response. Use the
  token returned from the previous request continue results where the previous request ended.
"""
function search_devices(filters; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "POST",
        "/devices",
        Dict{String,Any}("filters" => filters);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function search_devices(
    filters, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "POST",
        "/devices",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("filters" => filters), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Searches for Amazon Braket jobs that match the specified filter values.

# Arguments
- `filters`: The filter values to use when searching for a job.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned in the response. Use the
  token returned from the previous request to continue results where the previous request
  ended.
"""
function search_jobs(filters; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "POST",
        "/jobs",
        Dict{String,Any}("filters" => filters);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function search_jobs(
    filters, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "POST",
        "/jobs",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("filters" => filters), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Searches for tasks that match the specified filter values.

# Arguments
- `filters`: Array of SearchQuantumTasksFilter objects.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of results to return in the response.
- `"nextToken"`: A token used for pagination of results returned in the response. Use the
  token returned from the previous request continue results where the previous request ended.
"""
function search_quantum_tasks(filters; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "POST",
        "/quantum-tasks",
        Dict{String,Any}("filters" => filters);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function search_quantum_tasks(
    filters, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "POST",
        "/quantum-tasks",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("filters" => filters), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Add a tag to the specified resource.

# Arguments
- `resource_arn`: Specify the resourceArn of the resource to which a tag will be added.
- `tags`: Specify the tags to add to the resource.

"""
function tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return braket(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tags" => tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    resourceArn,
    tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tags" => tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

Remove tags from a resource.

# Arguments
- `resource_arn`: Specify the resourceArn for the resource from which to remove the tags.
- `tag_keys`: Specify the keys for the tags to remove from the resource.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return braket(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return braket(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
#=

  Global constraint assignment in Julia ConstraintSolver.jl

  From the Picat documentation
  """
  assignment(FDVars1,FDVars2): This constraint ensures that FDVars2 is a dual
  assignment of FDVars1, i.e., if the ith element of FDVars1 is j,
  then the jth element of FDVars2 is i.
  """

  Compare with the single array constraint inverse (tested in inverse.jl)

  The number of solutions is - of course - n!

  The number of solutions were X != Y for n=1.. is
  [0,0,2,14,94,644,4808,39556,360260,3619304]
  https://oeis.org/search?q=0%2C0%2C2%2C14%2C94%2C644%2C4808%2C39556%2C360260%2C3619304&sort=&language=english&go=Search
  https://oeis.org/A066052
  """
  A066052: Number of permutations in the symmetric group S_n with order >= 3.
  """

  This model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function assignment_ctr_test(n,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    @variable(model, 1 <= x[1:n] <= n, Int)
    @variable(model, 1 <= y[1:n] <= n, Int)
    @constraint(model, x in CS.AllDifferentSet())
    @constraint(model, y in CS.AllDifferentSet())

    assignment_ctr(model,x, y)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_diffs = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols")

        # if print_solutions
        for sol in 1:num_sols
            # println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x; result=sol))
            yy = convert.(Integer,JuMP.value.(y; result=sol))
            if print_solutions
                println("x:$xx y:$yy")
            end
            # println()
            if xx != yy
            #    println("x:$xx y:$yy")
                num_diffs += 1
            end
        end
        println("num_diffs:$num_diffs")
        # end
        # com = JuMP.backend(model).optimizer.model.inner
        # println("info:$(com.info)")
        # println("start_time:$(com.start_time)")
        # println("solve_time:$(com.solve_time)")
        println()
    else
        println("status:$status")
    end

    return status, num_diffs
end

function test_assignment_ctr_range(from,to)
    print_solutions=true
    all_solutions=true
    num_diffs = []
    for n in from:to
        println("n:$n")
        if n > 4
            print_solutions = false
        end
        if n > 9
            all_solutions=false
            print_solutions=true
        end
        @time status, diffs = assignment_ctr_test(n,print_solutions,all_solutions)
        if status != MOI.OPTIMAL
            break
        else
            push!(num_diffs,diffs)
        end
    end
    println("num_diffs: $num_diffs")
end

test_assignment_ctr_range(1,9)
# assignment_ctr_test(5)

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Euler equations in 2D space
function main(mesh::AbstractGrid,ρ₀::Function,m₀::Function,n₀::Function,E₀::Function,γ;
            CFL=0.04,P_order=2,T=0.0,α=1.1,bdr::Function=+)
quad1D_num=P_order*3+3
quad1D_p,quad1D_w=gausslegendre(quad1D_num)
P_dof=(P_order+1)^2
quad_num=quad1D_num^2
quad_p=Array{Float64,2}(undef,quad_num,2)
quad_w=Array{Float64,1}(undef,quad_num)
for ℓ₂=1:quad1D_num, ℓ₁=1:quad1D_num
    ℓ=(ℓ₂-1)*quad1D_num+ℓ₁
    quad_p[ℓ,1]=quad1D_p[ℓ₁]
    quad_p[ℓ,2]=quad1D_p[ℓ₂]
    quad_w[ℓ]=quad1D_w[ℓ₁]*quad1D_w[ℓ₂]
end
basis1DP=Array{Float64,2}(undef,P_order+1,quad1D_num)
Gradbasis1DP=Array{Float64,2}(undef,P_order+1,quad1D_num)
Grad2basis1DP=Array{Float64,2}(undef,P_order+1,quad1D_num)
for i = 0 : P_order
    for ℓ = 1 : quad1D_num
        basis1DP[i+1,ℓ]     = JacobiP(i,0,0,quad1D_p[ℓ])
        Gradbasis1DP[i+1,ℓ] = GradnJacobiP(i,0,0,quad1D_p[ℓ],1)
        Grad2basis1DP[i+1,ℓ] = GradnJacobiP(i,0,0,quad1D_p[ℓ],2)
    end
end
basisP=Array{Float64,2}(undef,P_dof,quad_num)
GradrbasisP=Array{Float64,2}(undef,P_dof,quad_num)
GradsbasisP=Array{Float64,2}(undef,P_dof,quad_num)
Grad2rrbasisP=Array{Float64,2}(undef,P_dof,quad_num)
Grad2rsbasisP=Array{Float64,2}(undef,P_dof,quad_num)
Grad2ssbasisP=Array{Float64,2}(undef,P_dof,quad_num)
boundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
GradrboundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
GradsboundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
Grad2rrboundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
Grad2rsboundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
Grad2ssboundaryP=Array{Float64,2}(undef,P_dof,quad1D_num*4)
for i=0:P_order, j=0:P_order, ℓ=1:quad_num
    basisP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad_p[ℓ,1])*JacobiP(j,0,0,quad_p[ℓ,2])
    GradrbasisP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad_p[ℓ,1],1)*JacobiP(j,0,0,quad_p[ℓ,2])
    GradsbasisP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad_p[ℓ,1])*GradnJacobiP(j,0,0,quad_p[ℓ,2],1)
    Grad2rrbasisP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad_p[ℓ,1],2)*JacobiP(j,0,0,quad_p[ℓ,2])
    Grad2rsbasisP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad_p[ℓ,1],1)*GradnJacobiP(j,0,0,quad_p[ℓ,2],1)
    Grad2ssbasisP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad_p[ℓ,1])*GradnJacobiP(j,0,0,quad_p[ℓ,2],2)
end
for i=0:P_order, j=0:P_order, ℓ=1:quad1D_num
    boundaryP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad1D_p[ℓ])*JacobiP(j,0,0,-1.)
    boundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=JacobiP(i,0,0,1.)*JacobiP(j,0,0,quad1D_p[ℓ])
    boundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=JacobiP(i,0,0,quad1D_p[ℓ])*JacobiP(j,0,0,1.)
    boundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=JacobiP(i,0,0,-1.)*JacobiP(j,0,0,quad1D_p[ℓ])

    GradrboundaryP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad1D_p[ℓ],1)*JacobiP(j,0,0,-1.)
    GradrboundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=GradnJacobiP(i,0,0,1.,1)*JacobiP(j,0,0,quad1D_p[ℓ])
    GradrboundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=GradnJacobiP(i,0,0,quad1D_p[ℓ],1)*JacobiP(j,0,0,1.)
    GradrboundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=GradnJacobiP(i,0,0,-1.,1)*JacobiP(j,0,0,quad1D_p[ℓ])

    GradsboundaryP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad1D_p[ℓ])*GradnJacobiP(j,0,0,-1.,1)
    GradsboundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=JacobiP(i,0,0,1.)*GradnJacobiP(j,0,0,quad1D_p[ℓ],1)
    GradsboundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=JacobiP(i,0,0,quad1D_p[ℓ])*GradnJacobiP(j,0,0,1.,1)
    GradsboundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=JacobiP(i,0,0,-1.)*GradnJacobiP(j,0,0,quad1D_p[ℓ],1)

    Grad2rrboundaryP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad1D_p[ℓ],2)*JacobiP(j,0,0,-1.)
    Grad2rrboundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=GradnJacobiP(i,0,0,1.,2)*JacobiP(j,0,0,quad1D_p[ℓ])
    Grad2rrboundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=GradnJacobiP(i,0,0,quad1D_p[ℓ],2)*JacobiP(j,0,0,1.)
    Grad2rrboundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=GradnJacobiP(i,0,0,-1.,2)*JacobiP(j,0,0,quad1D_p[ℓ])

    Grad2rsboundaryP[i*(P_order+1)+j+1,ℓ]=GradnJacobiP(i,0,0,quad1D_p[ℓ],1)*GradnJacobiP(j,0,0,-1.,1)
    Grad2rsboundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=GradnJacobiP(i,0,0,1.,1)*GradnJacobiP(j,0,0,quad1D_p[ℓ],1)
    Grad2rsboundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=GradnJacobiP(i,0,0,quad1D_p[ℓ],1)*GradnJacobiP(j,0,0,1.,1)
    Grad2rsboundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=GradnJacobiP(i,0,0,-1.,1)*GradnJacobiP(j,0,0,quad1D_p[ℓ],1)

    Grad2ssboundaryP[i*(P_order+1)+j+1,ℓ]=JacobiP(i,0,0,quad1D_p[ℓ])*GradnJacobiP(j,0,0,-1.,2)
    Grad2ssboundaryP[i*(P_order+1)+j+1,ℓ+quad1D_num]=JacobiP(i,0,0,1.)*GradnJacobiP(j,0,0,quad1D_p[ℓ],2)
    Grad2ssboundaryP[i*(P_order+1)+j+1,ℓ+2*quad1D_num]=JacobiP(i,0,0,quad1D_p[ℓ])*GradnJacobiP(j,0,0,1.,2)
    Grad2ssboundaryP[i*(P_order+1)+j+1,ℓ+3*quad1D_num]=JacobiP(i,0,0,-1.)*GradnJacobiP(j,0,0,quad1D_p[ℓ],2)
end
@.quad1D_p = ( quad1D_p + 1 ) / 2
@.quad1D_w = quad1D_w / 2
@.quad_p = ( quad_p + 1 ) / 2
@.quad_w = quad_w / 4
Coef_ρ = zeros(Float64,P_dof,mesh.NumofE)
Coef_m = zeros(Float64,P_dof,mesh.NumofE)
Coef_n = zeros(Float64,P_dof,mesh.NumofE)
Coef_E = zeros(Float64,P_dof,mesh.NumofE)
x = Array{Float64,2}(undef,quad_num,mesh.NumofE)
y = Array{Float64,2}(undef,quad_num,mesh.NumofE)
@inbounds Threads.@threads for j=1:mesh.NumofE
    x₁=mesh.Coordinates[mesh.Elements[1,j],1]
    x₂=mesh.Coordinates[mesh.Elements[2,j],1]
    y₁=mesh.Coordinates[mesh.Elements[1,j],2]
    y₂=mesh.Coordinates[mesh.Elements[4,j],2]
    for ℓ=1:quad_num
        x[ℓ,j]=(1-quad_p[ℓ,1])*x₁+quad_p[ℓ,1]*x₂
        y[ℓ,j]=(1-quad_p[ℓ,2])*y₁+quad_p[ℓ,2]*y₂
    end
end
ρ0=ρ₀(x,y)
m0=m₀(x,y)
n0=n₀(x,y)
E0=E₀(x,y)
@inbounds for j=1:mesh.NumofE
    for ℓ=1:quad_num
        for mm=1:P_dof
            Coef_ρ[mm,j]+=quad_w[ℓ]*basisP[mm,ℓ]*ρ0[ℓ,j]
            Coef_m[mm,j]+=quad_w[ℓ]*basisP[mm,ℓ]*m0[ℓ,j]
            Coef_n[mm,j]+=quad_w[ℓ]*basisP[mm,ℓ]*n0[ℓ,j]
            Coef_E[mm,j]+=quad_w[ℓ]*basisP[mm,ℓ]*E0[ℓ,j]
        end
    end
    for mm=1:P_dof
        Coef_ρ[mm,j]=Coef_ρ[mm,j]/(1/4)
        Coef_m[mm,j]=Coef_m[mm,j]/(1/4)
        Coef_n[mm,j]=Coef_n[mm,j]/(1/4)
        Coef_E[mm,j]=Coef_E[mm,j]/(1/4)
    end
end
time_init  = 0.0
time_final = T
time_now = time_init

num_iter=0
while time_now < time_final
    Coef_ρₜ,Coef_mₜ,Coef_nₜ,Coef_Eₜ,λₘₐₓ,μₘₐₓ =
        RHS(mesh,time_now,Coef_ρ,Coef_m,Coef_n,Coef_E,basisP,GradrbasisP,GradsbasisP,Grad2rrbasisP,Grad2rsbasisP,Grad2ssbasisP,
            boundaryP,GradrboundaryP,GradsboundaryP,Grad2rrboundaryP,Grad2rsboundaryP,Grad2ssboundaryP,
            quad_num,quad_p,quad_w,quad1D_num,quad1D_p,quad1D_w,γ,α,bdr)
    Δt=(CFL/(λₘₐₓ+μₘₐₓ)*sqrt(minimum(mesh.Tarea)))^((P_order+1)/3)
    Coef_ρ1=Coef_ρ+Δt*Coef_ρₜ
    Coef_m1=Coef_m+Δt*Coef_mₜ
    Coef_n1=Coef_n+Δt*Coef_nₜ
    Coef_E1=Coef_E+Δt*Coef_Eₜ

    Coef_ρ1ₜ,Coef_m1ₜ,Coef_n1ₜ,Coef_E1ₜ,λₘₐₓ,μₘₐₓ =
        RHS(mesh,time_now+Δt,Coef_ρ1,Coef_m1,Coef_n1,Coef_E1,basisP,GradrbasisP,GradsbasisP,Grad2rrbasisP,Grad2rsbasisP,Grad2ssbasisP,
            boundaryP,GradrboundaryP,GradsboundaryP,Grad2rrboundaryP,Grad2rsboundaryP,Grad2ssboundaryP,
            quad_num,quad_p,quad_w,quad1D_num,quad1D_p,quad1D_w,γ,α,bdr)
    Coef_ρ2=3/4*Coef_ρ+1/4*(Coef_ρ1+Δt*Coef_ρ1ₜ)
    Coef_m2=3/4*Coef_m+1/4*(Coef_m1+Δt*Coef_m1ₜ)
    Coef_n2=3/4*Coef_n+1/4*(Coef_n1+Δt*Coef_n1ₜ)
    Coef_E2=3/4*Coef_E+1/4*(Coef_E1+Δt*Coef_E1ₜ)

    Coef_ρ2ₜ,Coef_m2ₜ,Coef_n2ₜ,Coef_E2ₜ,λₘₐₓ,μₘₐₓ =
        RHS(mesh,time_now+Δt/2,Coef_ρ2,Coef_m2,Coef_n2,Coef_E2,basisP,GradrbasisP,GradsbasisP,Grad2rrbasisP,Grad2rsbasisP,Grad2ssbasisP,
            boundaryP,GradrboundaryP,GradsboundaryP,Grad2rrboundaryP,Grad2rsboundaryP,Grad2ssboundaryP,
            quad_num,quad_p,quad_w,quad1D_num,quad1D_p,quad1D_w,γ,α,bdr)
    Coef_ρ=1/3*Coef_ρ+2/3*(Coef_ρ2+Δt*Coef_ρ2ₜ)
    Coef_m=1/3*Coef_m+2/3*(Coef_m2+Δt*Coef_m2ₜ)
    Coef_n=1/3*Coef_n+2/3*(Coef_n2+Δt*Coef_n2ₜ)
    Coef_E=1/3*Coef_E+2/3*(Coef_E2+Δt*Coef_E2ₜ)

    time_now+=Δt
    num_iter+=1
end
return mesh,Coef_ρ,Coef_m,Coef_n,Coef_E,time_now,num_iter,basisP
end

function RHS(mesh,time_now,Coef_ρ,Coef_m,Coef_n,Coef_E,ϕ,∂ᵣϕ,∂ₛϕ,∂ᵣᵣϕ,∂ᵣₛϕ,∂ₛₛϕ,
    ϕᵧ,∂ᵣϕᵧ,∂ₛϕᵧ,∂ᵣᵣϕᵧ,∂ᵣₛϕᵧ,∂ₛₛϕᵧ,Nq,x̂,ω,Nqᵧ,x̂ᵧ,ωᵧ,γ,α,bdr::Function)
    Coef_ρₜ=zeros(eltype(Coef_ρ),size(Coef_ρ))
    Coef_mₜ=zeros(eltype(Coef_ρ),size(Coef_ρ))
    Coef_nₜ=zeros(eltype(Coef_ρ),size(Coef_ρ))
    Coef_Eₜ=zeros(eltype(Coef_ρ),size(Coef_ρ))
    λₘₐₓ=0.0
    μₘₐₓ=0.0
    P_dof=size(ϕ,1)
    nx=[0.,1.,0.,-1.]
    ny=[-1.,0.,1.,0.]
    EdtoVmap=[1 2;2 3;4 3;1 4]
    for j=1:mesh.NumofE
        for ℓ=1:Nq
            ρ=0.0;m=0.0;n=0.0;E=0.0;
            for mm=1:P_dof
                ρ+=Coef_ρ[mm,j]*ϕ[mm,ℓ]
                m+=Coef_m[mm,j]*ϕ[mm,ℓ]
                n+=Coef_n[mm,j]*ϕ[mm,ℓ]
                E+=Coef_E[mm,j]*ϕ[mm,ℓ]
            end
            c=sqrt(γ*(γ-1)*(E-1/2*m^2/ρ-1/2*n^2/ρ)/ρ)
            u=m/ρ
            v=n/ρ
            λₘₐₓ=max(λₘₐₓ,abs(u)+c)
            μₘₐₓ=max(μₘₐₓ,abs(v)+c)
        end
        for i=1:4
            EE=mesh.EtoEmap[i,j]
            if EE == -1
                V₁ = EdtoVmap[i,1]
                V₂ = EdtoVmap[i,2]
                xᵥ₁ = mesh.Coordinates[mesh.Elements[V₁,j],1]
                xᵥ₂ = mesh.Coordinates[mesh.Elements[V₂,j],1]
                yᵥ₁ = mesh.Coordinates[mesh.Elements[V₁,j],2]
                yᵥ₂ = mesh.Coordinates[mesh.Elements[V₂,j],2]
                for ℓ=1:Nqᵧ
                    ℓ⁻=ℓ+(i-1)*Nqᵧ
                    xᵧ=(1-x̂ᵧ[ℓ])*xᵥ₁+x̂ᵧ[ℓ]*xᵥ₂
                    yᵧ=(1-x̂ᵧ[ℓ])*yᵥ₁+x̂ᵧ[ℓ]*yᵥ₂
                    ρ⁻=0.0;m⁻=0.0;n⁻=0.0;E⁻=0.0;
                    for mm=1:P_dof
                        ρ⁻+=Coef_ρ[mm,j]*ϕᵧ[mm,ℓ⁻]
                        m⁻+=Coef_m[mm,j]*ϕᵧ[mm,ℓ⁻]
                        n⁻+=Coef_n[mm,j]*ϕᵧ[mm,ℓ⁻]
                        E⁻+=Coef_E[mm,j]*ϕᵧ[mm,ℓ⁻]
                    end
                    ρ⁺,m⁺,n⁺,E⁺ =
                        bdr(xᵧ,yᵧ,time_now,nx[i],ny[i],ρ⁻,m⁻,n⁻,E⁻)
                    c⁺=sqrt(γ*(γ-1)*(E⁺-1/2*m⁺^2/ρ⁺-1/2*n⁺^2/ρ⁺)/ρ⁺)
                    u⁺=m⁺/ρ⁺
                    v⁺=n⁺/ρ⁺
                    λₘₐₓ=max(λₘₐₓ,abs(u⁺)+c⁺)
                    μₘₐₓ=max(μₘₐₓ,abs(v⁺)+c⁺)
                end
            end
        end
    end

    for j=1:mesh.NumofE
        Aⱼ=mesh.Tarea[j]
        ∂r∂x=mesh.Jacobians[1,1,j]
        ∂s∂y=mesh.Jacobians[2,2,j]
        for ℓ=1:Nq
            ρ=0.0;m=0.0;n=0.0;E=0.0;
            for mm=1:P_dof
                ρ+=Coef_ρ[mm,j]*ϕ[mm,ℓ]
                m+=Coef_m[mm,j]*ϕ[mm,ℓ]
                n+=Coef_n[mm,j]*ϕ[mm,ℓ]
                E+=Coef_E[mm,j]*ϕ[mm,ℓ]
            end
            u=m/ρ
            v=n/ρ
            □₁ₓ=m
            □₁ₙ=n
            □₂ₓ=(γ-1)*E+(3-γ)/2*m*u-(γ-1)/2*n*v
            □₂ₙ=m*v
            □₃ₓ=n*u
            □₃ₙ=(γ-1)*E-(γ-1)/2*m*u+(3-γ)/2*n*v
            □₄ₓ=γ*E*u-(γ-1)/2*m*u^2-(γ-1)/2*m*v^2
            □₄ₙ=γ*E*v-(γ-1)/2*n*u^2-(γ-1)/2*n*v^2
            Ψ₁ₓ=Aⱼ*ω[ℓ]*∂r∂x*□₁ₓ
            Ψ₁ₙ=Aⱼ*ω[ℓ]*∂s∂y*□₁ₙ
            Ψ₂ₓ=Aⱼ*ω[ℓ]*∂r∂x*□₂ₓ
            Ψ₂ₙ=Aⱼ*ω[ℓ]*∂s∂y*□₂ₙ
            Ψ₃ₓ=Aⱼ*ω[ℓ]*∂r∂x*□₃ₓ
            Ψ₃ₙ=Aⱼ*ω[ℓ]*∂s∂y*□₃ₙ
            Ψ₄ₓ=Aⱼ*ω[ℓ]*∂r∂x*□₄ₓ
            Ψ₄ₙ=Aⱼ*ω[ℓ]*∂s∂y*□₄ₙ
            for mm=1:P_dof
                Coef_ρₜ[mm,j]+=∂ᵣϕ[mm,ℓ]*Ψ₁ₓ+∂ₛϕ[mm,ℓ]*Ψ₁ₙ
                Coef_mₜ[mm,j]+=∂ᵣϕ[mm,ℓ]*Ψ₂ₓ+∂ₛϕ[mm,ℓ]*Ψ₂ₙ
                Coef_nₜ[mm,j]+=∂ᵣϕ[mm,ℓ]*Ψ₃ₓ+∂ₛϕ[mm,ℓ]*Ψ₃ₙ
                Coef_Eₜ[mm,j]+=∂ᵣϕ[mm,ℓ]*Ψ₄ₓ+∂ₛϕ[mm,ℓ]*Ψ₄ₙ
            end
        end
        for i=1:4
            hᵢ=mesh.Elength[i,j]
            EE=mesh.EtoEmap[i,j]
            pp=mod(i+2-1,4)+1
            if EE != -1
                for ℓ=1:Nqᵧ
                    ∂r′∂x′=mesh.Jacobians[1,1,EE]
                    ∂s′∂y′=mesh.Jacobians[2,2,EE]
                    ℓ⁻=ℓ+(i-1)*Nqᵧ
                    ℓ⁺=ℓ+(pp-1)*Nqᵧ
                    ρ⁻=0.0;m⁻=0.0;n⁻=0.0;E⁻=0.0;
                    ρ⁺=0.0;m⁺=0.0;n⁺=0.0;E⁺=0.0;
                    for mm=1:P_dof
                        ρ⁻+=Coef_ρ[mm,j]*ϕᵧ[mm,ℓ⁻]
                        m⁻+=Coef_m[mm,j]*ϕᵧ[mm,ℓ⁻]
                        n⁻+=Coef_n[mm,j]*ϕᵧ[mm,ℓ⁻]
                        E⁻+=Coef_E[mm,j]*ϕᵧ[mm,ℓ⁻]

                        ρ⁺+=Coef_ρ[mm,EE]*ϕᵧ[mm,ℓ⁺]
                        m⁺+=Coef_m[mm,EE]*ϕᵧ[mm,ℓ⁺]
                        n⁺+=Coef_n[mm,EE]*ϕᵧ[mm,ℓ⁺]
                        E⁺+=Coef_E[mm,EE]*ϕᵧ[mm,ℓ⁺]
                    end
                    u⁻=m⁻/ρ⁻
                    v⁻=n⁻/ρ⁻
                    u⁺=m⁺/ρ⁺
                    v⁺=n⁺/ρ⁺

                    □₁ₓ⁻=m⁻
                    □₁ₙ⁻=n⁻
                    □₂ₓ⁻=(γ-1)*E⁻+(3-γ)/2*m⁻*u⁻-(γ-1)/2*n⁻*v⁻
                    □₂ₙ⁻=m⁻*v⁻
                    □₃ₓ⁻=n⁻*u⁻
                    □₃ₙ⁻=(γ-1)*E⁻-(γ-1)/2*m⁻*u⁻+(3-γ)/2*n⁻*v⁻
                    □₄ₓ⁻=γ*E⁻*u⁻-(γ-1)/2*m⁻*u⁻^2-(γ-1)/2*m⁻*v⁻^2
                    □₄ₙ⁻=γ*E⁻*v⁻-(γ-1)/2*n⁻*u⁻^2-(γ-1)/2*n⁻*v⁻^2

                    □₁ₓ⁺=m⁺
                    □₁ₙ⁺=n⁺
                    □₂ₓ⁺=(γ-1)*E⁺+(3-γ)/2*m⁺*u⁺-(γ-1)/2*n⁺*v⁺
                    □₂ₙ⁺=m⁺*v⁺
                    □₃ₓ⁺=n⁺*u⁺
                    □₃ₙ⁺=(γ-1)*E⁺-(γ-1)/2*m⁺*u⁺+(3-γ)/2*n⁺*v⁺
                    □₄ₓ⁺=γ*E⁺*u⁺-(γ-1)/2*m⁺*u⁺^2-(γ-1)/2*m⁺*v⁺^2
                    □₄ₙ⁺=γ*E⁺*v⁺-(γ-1)/2*n⁺*u⁺^2-(γ-1)/2*n⁺*v⁺^2

                    □₁ₓ=1/2*(□₁ₓ⁻+□₁ₓ⁺+α*λₘₐₓ*(ρ⁻-ρ⁺)*nx[i])
                    □₁ₙ=1/2*(□₁ₙ⁻+□₁ₙ⁺+α*μₘₐₓ*(ρ⁻-ρ⁺)*ny[i])
                    □₂ₓ=1/2*(□₂ₓ⁻+□₂ₓ⁺+α*λₘₐₓ*(m⁻-m⁺)*nx[i])
                    □₂ₙ=1/2*(□₂ₙ⁻+□₂ₙ⁺+α*μₘₐₓ*(m⁻-m⁺)*ny[i])
                    □₃ₓ=1/2*(□₃ₓ⁻+□₃ₓ⁺+α*λₘₐₓ*(n⁻-n⁺)*nx[i])
                    □₃ₙ=1/2*(□₃ₙ⁻+□₃ₙ⁺+α*μₘₐₓ*(n⁻-n⁺)*ny[i])
                    □₄ₓ=1/2*(□₄ₓ⁻+□₄ₓ⁺+α*λₘₐₓ*(E⁻-E⁺)*nx[i])
                    □₄ₙ=1/2*(□₄ₙ⁻+□₄ₙ⁺+α*μₘₐₓ*(E⁻-E⁺)*ny[i])
                    Ψ₁=-hᵢ*ωᵧ[ℓ]*(□₁ₓ*nx[i]+□₁ₙ*ny[i])
                    Ψ₂=-hᵢ*ωᵧ[ℓ]*(□₂ₓ*nx[i]+□₂ₙ*ny[i])
                    Ψ₃=-hᵢ*ωᵧ[ℓ]*(□₃ₓ*nx[i]+□₃ₙ*ny[i])
                    Ψ₄=-hᵢ*ωᵧ[ℓ]*(□₄ₓ*nx[i]+□₄ₙ*ny[i])
                    for mm=1:P_dof
                        Coef_ρₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₁
                        Coef_mₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₂
                        Coef_nₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₃
                        Coef_Eₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₄
                    end
                end
            else
                V₁ = EdtoVmap[i,1]
                V₂ = EdtoVmap[i,2]
                xᵥ₁ = mesh.Coordinates[mesh.Elements[V₁,j],1]
                xᵥ₂ = mesh.Coordinates[mesh.Elements[V₂,j],1]
                yᵥ₁ = mesh.Coordinates[mesh.Elements[V₁,j],2]
                yᵥ₂ = mesh.Coordinates[mesh.Elements[V₂,j],2]
                for ℓ=1:Nqᵧ
                    ℓ⁻=ℓ+(i-1)*Nqᵧ
                    ℓ⁺=ℓ+(pp-1)*Nqᵧ
                    xᵧ=(1-x̂ᵧ[ℓ])*xᵥ₁+x̂ᵧ[ℓ]*xᵥ₂
                    yᵧ=(1-x̂ᵧ[ℓ])*yᵥ₁+x̂ᵧ[ℓ]*yᵥ₂
                    ρ⁻=0.0;m⁻=0.0;n⁻=0.0;E⁻=0.0;
                    for mm=1:P_dof
                        ρ⁻+=Coef_ρ[mm,j]*ϕᵧ[mm,ℓ⁻]
                        m⁻+=Coef_m[mm,j]*ϕᵧ[mm,ℓ⁻]
                        n⁻+=Coef_n[mm,j]*ϕᵧ[mm,ℓ⁻]
                        E⁻+=Coef_E[mm,j]*ϕᵧ[mm,ℓ⁻]
                    end
                    ρ⁺,m⁺,n⁺,E⁺=bdr(xᵧ,yᵧ,time_now,nx[i],ny[i],ρ⁻,m⁻,n⁻,E⁻)

                    u⁻=m⁻/ρ⁻
                    v⁻=n⁻/ρ⁻
                    u⁺=m⁺/ρ⁺
                    v⁺=n⁺/ρ⁺

                    □₁ₓ⁻=m⁻
                    □₁ₙ⁻=n⁻
                    □₂ₓ⁻=(γ-1)*E⁻+(3-γ)/2*m⁻*u⁻-(γ-1)/2*n⁻*v⁻
                    □₂ₙ⁻=m⁻*v⁻
                    □₃ₓ⁻=n⁻*u⁻
                    □₃ₙ⁻=(γ-1)*E⁻-(γ-1)/2*m⁻*u⁻+(3-γ)/2*n⁻*v⁻
                    □₄ₓ⁻=γ*E⁻*u⁻-(γ-1)/2*m⁻*u⁻^2-(γ-1)/2*m⁻*v⁻^2
                    □₄ₙ⁻=γ*E⁻*v⁻-(γ-1)/2*n⁻*u⁻^2-(γ-1)/2*n⁻*v⁻^2

                    □₁ₓ⁺=m⁺
                    □₁ₙ⁺=n⁺
                    □₂ₓ⁺=(γ-1)*E⁺+(3-γ)/2*m⁺*u⁺-(γ-1)/2*n⁺*v⁺
                    □₂ₙ⁺=m⁺*v⁺
                    □₃ₓ⁺=n⁺*u⁺
                    □₃ₙ⁺=(γ-1)*E⁺-(γ-1)/2*m⁺*u⁺+(3-γ)/2*n⁺*v⁺
                    □₄ₓ⁺=γ*E⁺*u⁺-(γ-1)/2*m⁺*u⁺^2-(γ-1)/2*m⁺*v⁺^2
                    □₄ₙ⁺=γ*E⁺*v⁺-(γ-1)/2*n⁺*u⁺^2-(γ-1)/2*n⁺*v⁺^2

                    □₁ₓ=1/2*(□₁ₓ⁻+□₁ₓ⁺+α*λₘₐₓ*(ρ⁻-ρ⁺)*nx[i])
                    □₁ₙ=1/2*(□₁ₙ⁻+□₁ₙ⁺+α*μₘₐₓ*(ρ⁻-ρ⁺)*ny[i])
                    □₂ₓ=1/2*(□₂ₓ⁻+□₂ₓ⁺+α*λₘₐₓ*(m⁻-m⁺)*nx[i])
                    □₂ₙ=1/2*(□₂ₙ⁻+□₂ₙ⁺+α*μₘₐₓ*(m⁻-m⁺)*ny[i])
                    □₃ₓ=1/2*(□₃ₓ⁻+□₃ₓ⁺+α*λₘₐₓ*(n⁻-n⁺)*nx[i])
                    □₃ₙ=1/2*(□₃ₙ⁻+□₃ₙ⁺+α*μₘₐₓ*(n⁻-n⁺)*ny[i])
                    □₄ₓ=1/2*(□₄ₓ⁻+□₄ₓ⁺+α*λₘₐₓ*(E⁻-E⁺)*nx[i])
                    □₄ₙ=1/2*(□₄ₙ⁻+□₄ₙ⁺+α*μₘₐₓ*(E⁻-E⁺)*ny[i])
                    Ψ₁=-hᵢ*ωᵧ[ℓ]*(□₁ₓ*nx[i]+□₁ₙ*ny[i])
                    Ψ₂=-hᵢ*ωᵧ[ℓ]*(□₂ₓ*nx[i]+□₂ₙ*ny[i])
                    Ψ₃=-hᵢ*ωᵧ[ℓ]*(□₃ₓ*nx[i]+□₃ₙ*ny[i])
                    Ψ₄=-hᵢ*ωᵧ[ℓ]*(□₄ₓ*nx[i]+□₄ₙ*ny[i])
                    for mm=1:P_dof
                        Coef_ρₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₁
                        Coef_mₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₂
                        Coef_nₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₃
                        Coef_Eₜ[mm,j]+=ϕᵧ[mm,ℓ⁻]*Ψ₄
                    end
                end
            end
        end
        for mm=1:P_dof
            Coef_ρₜ[mm,j]=Coef_ρₜ[mm,j]/(1/4*Aⱼ)
            Coef_mₜ[mm,j]=Coef_mₜ[mm,j]/(1/4*Aⱼ)
            Coef_nₜ[mm,j]=Coef_nₜ[mm,j]/(1/4*Aⱼ)
            Coef_Eₜ[mm,j]=Coef_Eₜ[mm,j]/(1/4*Aⱼ)
        end
    end
    return Coef_ρₜ,Coef_mₜ,Coef_nₜ,Coef_Eₜ,λₘₐₓ,μₘₐₓ
end


text:
"""
Wrapper to indicate that a given handler should not be handled within `autorun`
"""
struct NoAutoRun{T}
  value::T
end

"""
    noautorun(handlers...)::Function

Constructs a wrapper function which can be used within @syntax_eff to skip all given handler types within the autorun.

Example
-------
```julia
@syntax_eff noautorun(Vector, Identity) begin
  a = [1,2,3]
  b = Identity(a + 2)
end
```
will actually run no handler at all in the implicit autorun, as both handlers are marked for ignore.
"""
function noautorun(handlers...)
  value -> value isa Union{handlers...} ? NoAutoRun(value) : value
end

"""
    autorun(eff)

special effectrunner which recognizes effecttypes used within `eff` and calls the effects in order, such that the
first effect found will at the end be the most outer container, and the last different effect found will be the most
inner container of the result value.
"""
autorun(eff::Eff) = runlast_ifpossible(_autorun((), eff))
function _autorun(handlers, eff::Eff)
  # for autorun we only deal with simple handlers `handler = typeof(eff.effectful)` where `eff_applies(handler, eff.effectful)`
  handler = eff_autohandler(eff.effectful)

  if handler ∈ handlers || !eff_applies(handler, eff.effectful)
    # If we encounter a handler which we already triggered, we don't want to trigger it again, however
    # we need to make sure that subsequent unseen handlers will be triggered correctly.
    # Hence, we build a continuation with calling `_autorun`.

    # The same we want to happen if we found a handler which cannot be handled automatically.
    # We just skip it and leave it for later handling.
    interpreted_continuation = if isempty(eff.cont)
      Continuation()
    else
      Continuation(x -> _autorun(handlers, eff.cont(x)))
    end
    # unwrap NoAutoRun if found
    value = eff.effectful isa NoAutoRun ? eff.effectful.value : eff.effectful
    Eff(value, interpreted_continuation)
  else
    # If this is a new valid handler, we trigger standard handler interpretation, with the one difference, that before
    # recursing into runhandler on the interpreted_continuation, we want to handle all yet unseen nested handlers.
    # We do this by first calling `_autorun` before calling `runhandler` within the interpreted_continuation.
    handlers_new = (handler, handlers...)
    interpreted_continuation = if isempty(eff.cont)
      Continuation(x -> _eff_pure(handler, x))
    else
      Continuation(x -> runhandler(handler, _autorun(handlers_new, eff.cont(x))))
    end
    _eff_flatmap(handler, interpreted_continuation, eff.effectful)
  end
end

"""
    eff_autohandler(value) = Base.typename(typeof(value)).wrapper

Overwrite this if the default autohandler extraction does not work for your case.
E.g. for `value::Vector` the default would return `Array`, hence we need to overwrite it individually.
"""
eff_autohandler(value) = Base.typename(typeof(value)).wrapper

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
module ArrayFire

typealias AF ArrayFire

export AFArray, AF
# array.jl
export to_host, cast, lookup, available_backends, set_backend, jl_size
export device_ptr, unlock_device_ptr
# image.jl
export unwrap, wrap

@windows? (
begin
  const libaf = Libdl.find_library(["af"])
  #const libaf = abspath(joinpath(dirname(@__FILE__), "..", "deps", "af"))
end : begin
  const libaf = Libdl.find_library(["libaf"])
end)
#!isempty(libaf) || error("ArrayFire library cannot be found.")

const af_fin = [true]
finalizable() = af_fin[1]
set_finalizable(b::Bool) = af_fin[1] = b

include("types.jl")
include("libaf.jl")
include("array.jl")
include("device.jl")
include("math.jl")
include("vector.jl")
include("image.jl")

aftype(::Type{Float32}) = f32
aftype(::Type{Complex{Float32}}) = c32
aftype(::Type{Float64}) = f64
aftype(::Type{Complex{Float64}}) = c64
aftype(::Type{Bool}) = b8
aftype(::Type{Int32}) = s32
aftype(::Type{UInt32}) = u32
aftype(::Type{UInt8}) = u8
aftype(::Type{Int64}) = s64
aftype(::Type{UInt64}) = u64
aftype(::Type{Int16}) = s16
aftype(::Type{UInt16}) = u16

function checkerror(err)
  if err != AF_SUCCESS
    throw(Int(err))
  end
end

end # ArrayFire


text:
export FiniteGroup

export group_order
export period_length
export group_product
export group_inverse

export generate_subgroup
export issubgroup
export isabelian
export minimal_generating_set
export group_multiplication_table

export element, elements
export element_name, element_names

export group_isomorphism
export ishomomorphic

"""
    FiniteGroup

Finite group, with elements {1, 2, 3,..., n}. The identity element is always 1.
Can be constructed using `FiniteGroup(multiplication_table)`

# Fields
* `multiplication_table::Matrix{Int}`: multiplication table
* `period_lengths::Vector{Int}`: period length (order) of every element
* `inverses::Vector{Int}`: inverse of every element
* `conjugacy_classes::Vector{Vector{Int}}`: conjugacy classes

# Examples
```jldoctest
julia> using TightBindingLattice

julia> FiniteGroup([1 2; 2 1])
FiniteGroup([1 2; 2 1], [1, 2], [1, 2], [[1], [2]])
```
"""
struct FiniteGroup <: AbstractGroup
    multiplication_table::Matrix{Int}

    period_lengths::Vector{Int}
    inverses::Vector{Int}
    conjugacy_classes::Vector{Vector{Int}}

    """
        FiniteGroup(multiplication_table)
    """
    function FiniteGroup(mtab::AbstractMatrix{<:Integer})
        n_elem = size(mtab, 1)
        if size(mtab, 2) != n_elem
            throw(ArgumentError("Multiplication table should be a square matrix"))
        end
        # check identity and closure
        if (mtab[:,1] != 1:n_elem) || (mtab[1,:] != 1:n_elem)
            throw(ArgumentError("Element 1 should be identity"))
        end
        elements = BitSet(1:n_elem)
        # check closure and inverse (sudoku)
        for i in 2:n_elem
            if (BitSet(mtab[:,i]) != elements) || (BitSet(mtab[i,:]) != elements)
                throw(ArgumentError("Multiplication not a group"))
            end
        end
        # check associativity
        for i in 1:n_elem, j in 1:n_elem, k in 1:n_elem
            if mtab[mtab[i, j], k] != mtab[i, mtab[j, k]]
                throw(ArgumentError("Multiplication not associative"))
            end
        end
        # compute cycles
        period_lengths = zeros(Int, n_elem)
        for idx in 1:n_elem
            jdx = idx
            for i in 1:n_elem
                if jdx == 1
                    period_lengths[idx] = i
                    break
                end
                jdx = mtab[jdx, idx]
            end
        end
        # compute inverses
        inverses = zeros(Int, n_elem)
        for i in 1:n_elem
            for j in i:n_elem
                if mtab[i,j] == 1
                    inverses[i] = j
                    inverses[j] = i
                    break
                end
            end
        end
        # compute conjugacy classes
        conjugacy_classes = Vector{Int}[]
        let
            adjacency = [Int[] for i in 1:n_elem]
            for i in 1:n_elem, j in 1:n_elem
                k = mtab[ mtab[j,i], inverses[j] ]
                push!(adjacency[i], k)
                push!(adjacency[k], i)
            end
            for i in eachindex(adjacency)
                sort!(adjacency[i])
                unique!(adjacency[i])
            end
            visited = falses(n_elem)
            for i in 1:n_elem
                visited[i] && continue
                visited[adjacency[i]] .= true
                push!(conjugacy_classes, adjacency[i])
            end
        end

        @assert all(
            period_lengths[first(cc)] == period_lengths[i]
            for cc in conjugacy_classes for i in cc
        )

        new(mtab, period_lengths, inverses, conjugacy_classes)
    end
end

function Base.:(==)(lhs::FiniteGroup, rhs::FiniteGroup)
    return (lhs.multiplication_table == rhs.multiplication_table)
end

"""
    element(group, idx)

Return the element of index `idx`. For `FiniteGroup`, this is somewhat meaningless
since the `idx`th element is `idx`. The sole purpose of this function is the bounds checking.
"""
element(group::FiniteGroup, idx) = Base.OneTo(group_order(group))[idx]

"""
    elements(group)

Return the elements of the group.
"""
elements(group::FiniteGroup) = Base.OneTo(group_order(group))

"""
    element_name(group, idx)

Return the name of element at index `idx`, which is just the string of `idx`.
"""
element_name(group::FiniteGroup, idx) = string.(element(group, idx))

"""
    element_names(group)

Return the names of element.
"""
element_names(group::FiniteGroup) = string.(elements(group))

"""
    group_order(group)

Order of group (i.e. number of elements)
"""
group_order(group::FiniteGroup) = size(group.multiplication_table, 1)

"""
    group_order(group, g)

Order of group element (i.e. period length)
"""
group_order(group::FiniteGroup, g) = group.period_lengths[g]

"""
    period_length(group, g)

Order of group element (i.e. period length)
"""
period_length(group::FiniteGroup, g) = group.period_lengths[g]

"""
    group_multiplication_table(group)

Return multiplcation table of the group.
"""
group_multiplication_table(group::FiniteGroup) = group.multiplication_table

"""
    isabelian(group)

Check if the group is abelian.
"""
function isabelian(group::FiniteGroup)
    return group.multiplication_table == transpose(group.multiplication_table)
end

"""
    group_product(group)

Return a function which computes the group product.
"""
function group_product(group::FiniteGroup) # a bit like currying
    function product(lhs::Integer, rhs::Integer)
        return group.multiplication_table[lhs, rhs]
    end
    function product(lhs::Integer, rhs::AbstractSet{<:Integer})
        return BitSet([group.multiplication_table[lhs, x] for x in rhs])
    end
    function product(lhs::AbstractSet{<:Integer}, rhs::Integer)
        return BitSet([group.multiplication_table[x, rhs] for x in lhs])
    end
    function product(lhs::AbstractSet{<:Integer}, rhs::AbstractSet{<:Integer})
        return BitSet([group.multiplication_table[x, y] for x in lhs for y in rhs])
    end
    return product
end

"""
    group_product(group, lhs, rhs)

Return the result of group multiplication of `lhs` and `rhs`.
If both `lhs` and `rhs` are integers, return an integer.
If either of them is a set (`AbstractSet`) of integers, then return a `BitSet`.
"""
function group_product(group::FiniteGroup, lhs::Integer, rhs::Integer)
    return group.multiplication_table[lhs, rhs]
end

function group_product(group::FiniteGroup, lhs::AbstractSet{<:Integer}, rhs::Integer)
    return BitSet([group_product(group, x, rhs) for x in lhs])
end

function group_product(group::FiniteGroup, lhs::Integer, rhs::AbstractSet{<:Integer})
    return BitSet([group_product(group, lhs, x) for x in rhs])
end

function group_product(
    group::FiniteGroup,
    lhs::AbstractSet{<:Integer},
    rhs::AbstractSet{<:Integer},
)
    return BitSet([group_product(group, x, y) for x in lhs for y in rhs])
end

"""
    group_inverse(group)

Get a function which gives inverse.
"""
function group_inverse(group::FiniteGroup)
    inverse(idx::Integer) = group.inverses[idx]
    inverse(idx::AbstractVector{<:Integer}) = group.inverses[idx]
    return inverse
end

"""
    group_inverse(group, g)

Get inverse of element/elements `g`.
"""
group_inverse(group::FiniteGroup, g::Integer) = group.inverses[g]
group_inverse(group::FiniteGroup, g::AbstractVector{<:Integer}) = group.inverses[g]

"""
    conjugacy_class(group::FiniteGroup, i::Integer)

Conjugacy class of the element `i`.
"""
function conjugacy_class(group::FiniteGroup, i::Integer)
    return findfirst(
        c -> let j = searchsortedfirst(c, i)
            j <= length(c) && c[j] == i
        end,
        group.conjugacy_classes
    )
end

"""
    generate_subgroup(group::FiniteGroup, idx::Integer)

subgroup generated by `generators`. ⟨ {`idx`} ⟩
"""
function generate_subgroup(group::FiniteGroup, idx::Integer)
    out = BitSet()
    sizehint!(out, group_order(group, idx))
    jdx = 1
    for i in 1:group_order(group, idx)
        push!(out, jdx)
        jdx = group_product(group, jdx, idx)
    end
    @assert jdx == 1
    return out
end

"""
    generate_subgroup(group::FiniteGroup, generators)

subgroup generated by `generators`. ⟨ S ⟩
"""
function generate_subgroup(
    group::FiniteGroup,
    generators::G
) where {G<:Union{<:AbstractSet{<:Integer}, <:AbstractVector{<:Integer}}}
    change = true
    subgroup = BitSet(generators)
    push!(subgroup, 1)
    while change
        change = false
        for g1 in generators, g2 in subgroup
            g3 = group_product(group, g1, g2)
            if !(g3 in subgroup)
                change = true
                push!(subgroup, g3)
            end
        end
    end
    return subgroup
end

"""
    issubgroup(group, subset)

Check whether the given subset is a subgroup of `group`.
"""
function issubgroup(group::FiniteGroup, subset::AbstractSet{<:Integer})
    return all(group_product(group, x, y) in subset for x in subset for y in subset)
end

"""
    minimal_generating_set(group)

Get minimally generating set of the finite group.
"""
function minimal_generating_set(group::FiniteGroup)
    ord_group::Int = group_order(group)
    element_queue::Vector{Tuple{Int, Int}} = collect(enumerate(group.period_lengths))
    sort!(element_queue, by=item->(-item[2], item[1]))

    queue_begin = 1
    span = BitSet([1])
    generators = Int[]
    while queue_begin <= ord_group && length(span) < ord_group
        next_index = queue_begin
        next_elem = element_queue[queue_begin][1]
        next_span = generate_subgroup(group, group_product(group, span, next_elem))
        for i in (queue_begin+1):ord_group
            (g, _) = element_queue[i]
            new_span = generate_subgroup(group, group_product(group, span, g))
            if length(new_span) > length(next_span)
                next_index = i
                next_elem = g
                next_span = new_span
            end
        end
        queue_begin = next_index + 1
        span = next_span
        push!(generators, next_elem)
    end
    return generators
end

"""
    group_isomorphism(group1, group2)

Find the isomorphism ϕ: G₁ → G₂. Return nothing if not isomorphic.
"""
function group_isomorphism(group1::FiniteGroup, group2::FiniteGroup)
    group_order(group1) != group_order(group2) && return nothing
    sort(group1.period_lengths) != sort(group2.period_lengths) && return nothing
    let cl1 = sort(length.(group1.conjugacy_classes)),
        cl2 = sort(length.(group1.conjugacy_classes))
        cl1 != cl2 && return nothing
    end
    ord_group = group_order(group1)

    pl1_list = group1.period_lengths
    pl2_list = group2.period_lengths
    cci1_list = zeros(Int, ord_group)

    for (cc_index, cc) in enumerate(group1.conjugacy_classes), i in cc
        cci1_list[i] = cc_index
    end

    element_map = zeros(Int, ord_group)
    class_map = zeros(Int, length(group1.conjugacy_classes))

    function suggest(i::Integer)
        # @assert element_map[i] == 0
        cci1 = cci1_list[i]
        cc1 = group1.conjugacy_classes[cci1]
        if class_map[cci1] != 0
            return (j for j in group2.conjugacy_classes[class_map[cci1]]
                    if !(j in element_map) && (pl2_list[j] == pl1_list[i]))
        else
            return (j for (icc2, cc2) in enumerate(group2.conjugacy_classes)
                          if length(cc2) == length(cc1)
                      for j in cc2
                          if !(j in element_map) && (pl2_list[j] == pl1_list[i]))
        end
    end

    function dfs(i::Integer)
        i > ord_group && return true
        element_map[i] != 0 && return dfs(i+1)
        for j in suggest(i)
            cset = class_map[cci1_list[i]] == 0
            element_map[i] = j
            class_map[cci1_list[i]] = conjugacy_class(group2, j)
            if dfs(i+1)
                return true
            else
                element_map[i] = 0
                if cset
                    class_map[cci1_list[i]] = 0
                end
            end
        end
        return false
    end

    if dfs(1)
        return element_map
    else
        return nothing
    end
end

"""
    group_multiplication_table(elements, product=(*))

Generate a multiplication table from elements with product.
"""
function group_multiplication_table(
    elements::AbstractVector{ElementType},
    product::Function=Base.:(*)
) where {ElementType}
    element_lookup = Dict(k=>i for (i, k) in enumerate(elements))
    ord_group = length(elements)
    mtab = zeros(Int, (ord_group, ord_group))
    for i in 1:ord_group, j in 1:ord_group
        mtab[i,j] = element_lookup[ product(elements[i], elements[j]) ]
    end
    return mtab
end

"""
    ishomomorphic(group, representation; product=(*), equal=(==))

Check whether `representation` is homomorphic to `group` under `product` and `equal`,
order preserved.
"""
function ishomomorphic(
    group::FiniteGroup,
    representation::AbstractVector;
    product::Function=Base.:(*),
    equal::Function=Base.:(==)
)
    ord_group = group_order(group)
    if length(representation) != ord_group
        return false
    end
    for i in 1:ord_group, j in 1:ord_group
        if !equal( product(representation[i], representation[j]),
                   representation[ group_product(group, i, j)] )
            return false
        end
    end
    return true
end

# function group_isomorphism_naive(group1::FiniteGroup, group2::FiniteGroup)
#     group_order(group1) != group_order(group2) && return nothing
#     sort(group1.period_lengths) != sort(group2.period_lengths) && return nothing
#
#     ord_group = group_order(group1)
#     element_groups1 = Dict{Tuple{Int, Int}, Vector{Int}}() # group by period lengths and conjugacy class size
#     element_groups2 = Dict{Tuple{Int, Int}, Vector{Int}}() # group by period lengths
#
#     for i in 1:group_order(group1)
#         pl = group1.period_lengths[i]
#         cc = length(group1.conjugacy_classes[conjugacy_class(group1, i)])
#         if !haskey(element_groups1, (pl, cc))
#             element_groups1[(pl, cc)] = Int[]
#         end
#         push!(element_groups1[(pl, cc)], i)
#     end
#
#     for i in 1:group_order(group2)
#         pl = group2.period_lengths[i]
#         cc = length(group2.conjugacy_classes[conjugacy_class(group2, i)])
#         if !haskey(element_groups2, (pl, cc))
#             element_groups2[(pl, cc)] = Int[]
#         end
#         push!(element_groups2[(pl, cc)], i)
#     end
#
#     #q = sort([(pl, i) for (pl, els) in element_groups1 for i in els], rev=true)
#     plccs = sort(collect(keys(element_groups1)))
#     mapping = zeros(Int, ord_group)
#     for perm_set in Iterators.product([
#             permutations(1:length(element_groups1[plcc]), length(element_groups1[plcc]))
#             for plcc in plccs
#         ]...)
#         #@assert length(pls) == length(perm_set)
#         for (ipl, (plcc, perm)) in enumerate(zip(plccs, perm_set))
#             elg1, elg2 = element_groups1[plcc], element_groups2[plcc]
#             for j in 1:length(perm)
#                 mapping[ element_groups1[plcc][j] ] = element_groups2[plcc][ perm[j] ]
#             end
#         end
#         #mtab1p = zeros(Int, (ord_group, ord_group))
#         isiso = true
#         for i in 1:ord_group, j in 1:ord_group
#             #mtab1p[mapping[i], mapping[j]] = mapping[group1.multiplication_table[i, j]]
#             if group2.multiplication_table[mapping[i], mapping[j]] != mapping[group1.multiplication_table[i, j]]
#                 isiso = false
#                 break
#             end
#         end
#         isiso && return mapping
#         #!isiso && continue
#         #mtab1p == group2.multiplication_table && return mapping
#     end
#     return nothing
# end


text:
module Fragmentation
using ..scenario, ..Utils, ..Types, ..ExternalBallistics
using LinearAlgebra
export frag, fragFly!
function frag(proj::Projectile2D, zones::Array{Zone,1}, shapes::Array{FragShapes,1},cf::ShapeFactor, αₑ_bar::Array{Float64,1})
    θ  = asin(-proj.velocity[3]/norm(proj.velocity))
    Ψ = atan(proj.velocity[2]/proj.velocity[1])
    T = [cos(θ)*cos(Ψ) -sin(Ψ) sin(θ)*cos(Ψ); cos(θ)*sin(Ψ) cos(Ψ) sin(θ)*sin(Ψ); -sin(θ) 0.0 cos(θ)]
    frags = []
    for z in zones
        f = createFragment(z,shapes,cf, proj, αₑ_bar)

    #=    for i in f

            v = i.velocity#[i.vx, i.vy, i.vz]
            vrot = rotation(v, T)
            i.velocity = vrot
            #i.vy = vrot[2]
            #i.vz = vrot[3]
        end
        =#
        push!(frags,f)
    end
    return frags
end

function fragFly!(fragments, proj::Projectile2D, w_bar::Array{Float64,1}, target::Target2D, g₀::Float64, ω_bar::Array{Float64,1}, coeffAero;R=6.356766*1e6, tspan=(0.0,1000.0))
    for i=1:length(fragments)
        for j=1:length(fragments[i])

            u0=[proj.position[1], proj.position[2], proj.position[3],fragments[i][j].velocity[1], fragments[i][j].velocity[2], fragments[i][j].velocity[3]]

            p= [w_bar, fragments[i][j].Aₚ, R, g₀,ω_bar, fragments[i][j].mass,target.position,coeffAero,fragments[i][j]]
            ExternalBallistics.Pmm.trajectoryPMM!(u0, tspan, p, fragments[i][j])
        end

    end

end

end


text:
# Grid size
const WIDTH  = 1000
const HEIGHT = 1000

struct Instruction
    instruction::String
    start_x    ::Int16
    end_x      ::Int16
    start_y    ::Int16
    end_y      ::Int16
end

# Executes an instruction and returns the grid
function execute(ins::Instruction, grid)
    # toggle
    if ins.instruction == "toggle"
        grid[ins.start_x:ins.start_y,ins.end_x:ins.end_y] = .!grid[ins.start_x:ins.start_y,ins.end_x:ins.end_y]
    # turn off
    elseif ins.instruction == "turn off"
        grid[ins.start_x:ins.start_y,ins.end_x:ins.end_y] .= 0
    # turn on
    else
         grid[ins.start_x:ins.start_y,ins.end_x:ins.end_y] .= 1
    end

    return grid
end

function solve(input_file)

    # 1000 * 1000 binary grid
    grid = BitArray(zeros(WIDTH, HEIGHT))

    lines = readlines(input_file)
    for line in lines
        instruction_regexp = r"(turn on|toggle|turn off) ([0-9]+),([0-9]+) through ([0-9]+),([0-9]+)"
        tokens = match(instruction_regexp, line)
        instruction = Instruction(tokens[1],
                                  parse(Int16, tokens[2]),
                                  parse(Int16, tokens[3]),
                                  parse(Int16, tokens[4]),
                                  parse(Int16, tokens[5])
        )

        execute(instruction, grid)
    end

    result = sum(grid)
    return result
end

println(solve("input.txt"))


text:
using Soss

N = 10000

m = @model x, λ begin
    σ ~ Exponential(λ)
    β ~ Normal(0,1)
    y ~ For(1:N) do j
        Normal(x[j] * β, σ)
    end
    return y
end

x = randn(N)
λ = 1.0

trace = simulate(m(x=x, λ=1.0)).trace
y = trace.y

post = m(; x, λ) | (; y)
ℓ = symlogdensity(post)


text:
using SparseArrays

function strang_matrix(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      abs(i-j)<=1 && (A[i,j]+=1)
      i==j && (A[i,j]-=3)
  end
  A
end

# tests for full and sparse function
@testset "Full and Sparse functions:" begin
    N = 100
    d_order = 2
    approx_order = 2
    x = collect(1:1.0:N).^2

    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:Dirichlet0,:Dirichlet0)
    mat = convert(Array,A)
    sp_mat = sparse(A)
    @test mat == sp_mat

    @test convert(Array, A, 10) == strang_matrix(10) # Strang Matrix is defined with the center term +ve
    @test convert(Array, A, N) == strang_matrix(N) # Strang Matrix is defined with the center term +ve
    @test convert(Array,A) == sp_mat
    @test opnorm(A, Inf) == opnorm(mat, Inf)

    # testing correctness
    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N).^4 - 2*collect(1:1.0:N).^3 + collect(1:1.0:N).^2;
    y = convert(Array{BigFloat, 1}, y)

    A = DerivativeOperator{BigFloat}(d_order,approx_order,one(BigFloat),N,:Dirichlet0,:Dirichlet0)
    boundary_points = A.boundary_point_count
    mat = convert(Array, A, N)
    sp_mat = sparse(A)
    @test mat == sp_mat

    res = A*y
    @test res[boundary_points[1] + 1: N - boundary_points[2]] ≈ 24.0*ones(N - sum(boundary_points)) atol=10.0^-approx_order
    @test A*y ≈ mat*y atol=10.0^-approx_order
    @test A*y ≈ sp_mat*y atol=10.0^-approx_order
    @test sp_mat*y ≈ mat*y atol=10.0^-approx_order
end

@testset "Indexing tests" begin
    N = 1000
    d_order = 4
    approx_order = 10

    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:Dirichlet0,:Dirichlet0)
    @test A[1,1] ≈ 13.717407 atol=1e-4
    @test A[:,1] == (convert(Array,A))[:,1]
    @test A[10,20] == 0

    for i in 1:N
        @test A[i,i] == A.stencil_coefs[div(A.stencil_length, 2) + 1]
    end

    # Indexing Tests
    N = 1000
    d_order = 2
    approx_order = 2

    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N,:Dirichlet0,:Dirichlet0)
    M = convert(Array,A)

    @test A[1,1] == -2.0
    @test A[1:4,1] == M[1:4,1]
    @test A[5,2:10] == M[5,2:10]
    @test A[60:100,500:600] == M[60:100,500:600]
end

@testset begin "Operations on matrices"
    N = 51
    M = 101
    d_order = 2
    approx_order = 2

    xarr = range(0,stop=1,length=N)
    yarr = range(0,stop=1,length=M)
    dx = xarr[2]-xarr[1]
    dy = yarr[2]-yarr[1]
    F = [x^2+y for x = xarr, y = yarr]

    A = DerivativeOperator{Float64}(d_order,approx_order,dx,length(xarr),:None,:None)
    B = DerivativeOperator{Float64}(d_order,approx_order,dy,length(yarr),:None,:None)

    @test A*F ≈ 2*ones(N,M) atol=1e-2
    @test F*B ≈ 8*ones(N,M) atol=1e-2
    @test A*F*B ≈ zeros(N,M) atol=1e-2

    G = [x^2+y^2 for x = xarr, y = yarr]

    @test A*G ≈ 2*ones(N,M) atol=1e-2
    @test G*B ≈ 8*ones(N,M) atol=1e-2
    @test A*G*B ≈ zeros(N,M) atol=1e-2
end

@testset "Linear combinations of operators" begin
    # Only tests the additional functionality defined in "operator_combination.jl"
    N = 10
    Random.seed!(0); LA = DiffEqArrayOperator(rand(N,N))
    LD = DerivativeOperator{Float64}(2,2,1.0,N,:Dirichlet0,:Dirichlet0)
    @test_broken begin
      L = 1.1*LA - 2.2*LD + 3.3*I
      # Builds convert(L) the brute-force way
      fullL = zeros(N,N)
      v = zeros(N)
      for i = 1:N
          v[i] = 1.0
          fullL[:,i] = L*v
          v[i] = 0.0
      end
      @test convert(L) ≈ fullL
      @test exp(L) ≈ exp(fullL)
      for p in [1,2,Inf]
          @test opnorm(L,p) ≈ opnorm(fullL,p)
          @test opnormbound(L,p) ≈ 1.1*opnorm(LA,p) + 2.2*opnorm(LD,p) + 3.3
      end
  end
end


text:
module BenchAutoMPO

using BenchmarkTools
using ITensors

suite = BenchmarkGroup()

let N = 30
  s = siteinds("S=1/2", N; conserve_qns=false)
  a = AutoMPO()
  for k in 1:N, l in 1:N, m in 1:N, n in 1:N
    a .+= "projDn", k, "projDn", l, "projDn", m, "projDn", n
  end

  # Precompile
  MPO(a, s)

  suite["Quartic Hamiltonian"] = @benchmarkable MPO($a, $s)
end

let N = 30
  s = siteinds("S=1/2", N; conserve_qns=true)
  a = AutoMPO()
  for k in 1:N, l in 1:N, m in 1:N, n in 1:N
    a .+= "projDn", k, "projDn", l, "projDn", m, "projDn", n
  end

  # Precompile
  MPO(a, s)

  suite["Quartic QN Hamiltonian"] = @benchmarkable MPO($a, $s)
end

end

BenchAutoMPO.suite


text:
export PDataMA57, AInv

mutable struct PDataMA57{T} <: PDataFact{T}
    L::SparseMatrixCSC{T,Int64} # sparse L
    D::SparseMatrixCSC{T,Int64} # block diagonal 1X1 and 2X2
    #P::Array{T,2}  # permutation matrix: should be vector
    # LDL' = PHP'
    pp::Array{Int64,1}   # permutation vector LDL' = H[pp,pp]
    # P*v = v[pp]
    # P'*v = v[invperm(pp)]
    s::Array{T,1}  # Scaling matrix S=diagm(s)
    Δ::Array{T,1}  # diagonal, eigenvalues of D
    Q::SparseMatrixCSC{T,Int64}  # orthogonal matrix, eigenvectors of D:  should be sparse
    # QΔQ'  =  D
    g̃::Array{T,1}  # transformed gradient
    λ::T
    success::Bool        # previous iteration was successfull
    OK::Bool             # preprocess success

    PDataMA57() = new{Nothing}()
    PDataMA57(L, D, pp, s, Δ, Q, g̃, l, success, OK) =
        new{eltype(g̃)}(L, D, pp, s, Δ, Q, g̃, l, success, OK)
end

function preprocess(PData::PDataMA57, H, g, n1, n2)
    M = Ma57
    L = SparseMatrixCSC{Float64,Int64}
    D57 = SparseMatrixCSC{Float64,Int64}
    pp = Array{Int64,1}
    s = Array{Float64}
    ρ = Float64
    ncomp = Int64

    H57 = convert(SparseMatrixCSC{Float64,Int64}, H)  #  Hard coded Int
    try
        M = Ma57(H57)#, print_level = -1)
        ma57_factorize(M)
    catch
        println("*******   Problem in MA57_0")
        M = Ma57(H57, print_level = -1)
        ma57_factorize(M)
        res = PDataMA57()
        res.OK = false
        return res
    end
    ###### Tried in an alternative version ####
    # M.control.cntl[1] = 1.0 # diff
    # M.control.icntl[7] = 0 # diff
    ######################
    try
        (L, D57, s, pp) = ma57_get_factors(M)
    catch
        println("*******   Problem after MA57_1")
        # println(" Cond(H) = $(cond(full(H)))")
        res = PDataMA57_0()
        res.OK = false
        return res
    end

    #################  Future object BlockDiag operator?
    vD1 = diag(D57)       # create internal representation for block diagonal D
    vD2 = diag(D57, 1)     #
    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
    vQ2 = zeros(length(vD2))      #
    vQ2m = zeros(length(vD2))     #
    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
    # if D diagonal, nothing more will be computed

    i = 1
    while i < length(vD1)
        if vD2[i] == 0.0
            i += 1
        else
            mA = [vD1[i] vD2[i]; vD2[i] vD1[i+1]] #  2X2 submatrix
            # DiagmA, Qma = eig(mA)                   #  spectral decomposition of mA
            X = eigen(mA)
            DiagmA = X.values
            Qma = X.vectors
            veig[i] = DiagmA[1]
            vQ1[i] = Qma[1, 1]
            vQ2[i] = Qma[1, 2]
            vQ2m[i] = Qma[2, 1]
            vQ1[i+1] = Qma[2, 2]
            veig[i+1] = DiagmA[2]
            i += 2
        end
    end

    Q = sparse(SparseArrays.spdiagm(0 => vQ1, 1 => vQ2m, -1 => vQ2))
    Δ = veig

    l_m, = findmin(Δ)
    sg = s .* g
    L = SparseMatrixCSC{Float64,Int64}(L)  #### very important, the \ command doesn't work with Int32
    ĝ = L \ (sg[pp])
    g̃ = Q' * ĝ

    n_g = norm(g)
    λ = max(-l_m, 0.0)
    return PDataMA57(L, D57, pp, s, Δ, Q, g̃, λ, true, true)
end

function AInv(X::PDataMA57, d̃::Array{T,1}) where {T}
    d̂ = X.Q * d̃
    u = X.L' \ d̂
    return u[invperm(X.pp)] .* X.s
end


text:
export xdp_load, xdp_unload

@enum xdp_action begin
    XDP_ABORTED
    XDP_DROP
    XDP_PASS
    XDP_TX
    XDP_REDIRECT
end

function xdp_load(iface::String, bpf::Vector{UInt8})
    mktemp() do path, io
        write(io, bpf)
        flush(io)
        run(`ip link set $iface xdp obj $path verbose`)
    end
end
function xdp_unload(iface::String)
    run(`ip link set $iface xdp off`)
end


text:
const _clfft_error_codes = Dict{Int,Symbol}(
    cl.CL_INVALID_GLOBAL_WORK_SIZE => :CLFFT_INVALID_GLOBAL_WORK_SIZE,
    cl.CL_INVALID_MIP_LEVEL => :CLFFT_INVALID_MIP_LEVEL,
    cl.CL_INVALID_BUFFER_SIZE => :CLFFT_INVALID_BUFFER_SIZE,
    cl.CL_INVALID_GL_OBJECT => :CLFFT_INVALID_GL_OBJECT,
    cl.CL_INVALID_OPERATION => :CLFFT_INVALID_OPERATION,
    cl.CL_INVALID_EVENT => :CLFFT_INVALID_EVENT,
    cl.CL_INVALID_EVENT_WAIT_LIST => :CLFFT_INVALID_EVENT_WAIT_LIST,
    cl.CL_INVALID_GLOBAL_OFFSET => :CLFFT_INVALID_GLOBAL_OFFSET,
    cl.CL_INVALID_WORK_ITEM_SIZE => :CLFFT_INVALID_WORK_ITEM_SIZE,
    cl.CL_INVALID_WORK_GROUP_SIZE => :CLFFT_INVALID_WORK_GROUP_SIZE,
    cl.CL_INVALID_WORK_DIMENSION => :CLFFT_INVALID_WORK_DIMENSION,
    cl.CL_INVALID_KERNEL_ARGS => :CLFFT_INVALID_KERNEL_ARGS,
    cl.CL_INVALID_ARG_SIZE => :CLFFT_INVALID_ARG_SIZE,
    cl.CL_INVALID_ARG_VALUE => :CLFFT_INVALID_ARG_VALUE,
    cl.CL_INVALID_ARG_INDEX => :CLFFT_INVALID_ARG_INDEX,
    cl.CL_INVALID_KERNEL => :CLFFT_INVALID_KERNEL,
    cl.CL_INVALID_KERNEL_DEFINITION => :CLFFT_INVALID_KERNEL_DEFINITION,
    cl.CL_INVALID_KERNEL_NAME => :CLFFT_INVALID_KERNEL_NAME,
    cl.CL_INVALID_PROGRAM_EXECUTABLE => :CLFFT_INVALID_PROGRAM_EXECUTABLE,
    cl.CL_INVALID_PROGRAM => :CLFFT_INVALID_PROGRAM,
    cl.CL_INVALID_BUILD_OPTIONS => :CLFFT_INVALID_BUILD_OPTIONS,
    cl.CL_INVALID_BINARY => :CLFFT_INVALID_BINARY,
    cl.CL_INVALID_SAMPLER => :CLFFT_INVALID_SAMPLER,
    cl.CL_INVALID_IMAGE_SIZE => :CLFFT_INVALID_IMAGE_SIZE,
    cl.CL_INVALID_IMAGE_FORMAT_DESCRIPTOR => :CLFFT_INVALID_IMAGE_FORMAT_DESCRIPTOR,
    cl.CL_INVALID_MEM_OBJECT => :CLFFT_INVALID_MEM_OBJECT,
    cl.CL_INVALID_HOST_PTR => :CLFFT_INVALID_HOST_PTR,
    cl.CL_INVALID_COMMAND_QUEUE => :CLFFT_INVALID_COMMAND_QUEUE,
    cl.CL_INVALID_QUEUE_PROPERTIES => :CLFFT_INVALID_QUEUE_PROPERTIES,
    cl.CL_INVALID_CONTEXT => :CLFFT_INVALID_CONTEXT,
    cl.CL_INVALID_DEVICE => :CLFFT_INVALID_DEVICE,
    cl.CL_INVALID_PLATFORM => :CLFFT_INVALID_PLATFORM,
    cl.CL_INVALID_DEVICE_TYPE => :CLFFT_INVALID_DEVICE_TYPE,
    cl.CL_INVALID_VALUE => :CLFFT_INVALID_VALUE,
    cl.CL_MAP_FAILURE => :CLFFT_MAP_FAILURE,
    cl.CL_BUILD_PROGRAM_FAILURE => :CLFFT_BUILD_PROGRAM_FAILURE,
    cl.CL_IMAGE_FORMAT_NOT_SUPPORTED => :CLFFT_IMAGE_FORMAT_NOT_SUPPORTED,
    cl.CL_IMAGE_FORMAT_MISMATCH => :CLFFT_IMAGE_FORMAT_MISMATCH,
    cl.CL_MEM_COPY_OVERLAP => :CLFFT_MEM_COPY_OVERLAP,
    cl.CL_PROFILING_INFO_NOT_AVAILABLE => :CLFFT_PROFILING_INFO_NOT_AVAILABLE,
    cl.CL_OUT_OF_HOST_MEMORY => :CLFFT_OUT_OF_HOST_MEMORY,
    cl.CL_OUT_OF_RESOURCES => :CLFFT_OUT_OF_RESOURCES,
    cl.CL_MEM_OBJECT_ALLOCATION_FAILURE => :CLFFT_MEM_OBJECT_ALLOCATION_FAILURE,
    cl.CL_COMPILER_NOT_AVAILABLE => :CLFFT_COMPILER_NOT_AVAILABLE,
    cl.CL_DEVICE_NOT_AVAILABLE => :CLFFT_DEVICE_NOT_AVAILABLE,
    cl.CL_DEVICE_NOT_FOUND => :CLFFT_DEVICE_NOT_FOUND,
    cl.CL_SUCCESS => :CLFFT_SUCCESS,
    4096 => :CLFFT_BUGCHECK,
    4097 => :CLFFT_NOTIMPLEMENTED,
    4098 => :CLFFT_TRANSPOSED_NOTIMPLEMENTED,
    4099 => :CLFFT_FILE_NOT_FOUND,
    4100 => :CLFFT_FILE_CREATE_FAILURE,
    4101 => :CLFFT_VERSION_MISMATCH,
    4102 => :CLFFT_INVALID_PLAN,
    4103 => :CLFFT_DEVICE_NO_DOUBLE,
)


text:
# using StaticArrays

include("typestructs.jl")
include("Vec.jl")
const dtypes = Union{UInt8, Int8, UInt16, Int16, Int32, Float32, Float64}
size_t = UInt64

using CxxWrap
@wrapmodule(joinpath(@__DIR__,"lib","libopencv_julia"), :cv_wrap)
function __init__()
    @initcxx

    if jlopencv_core_get_sizet()==4
        size_t = UInt32
    end
end
const Scalar = Union{Tuple{}, Tuple{Number}, Tuple{Number, Number}, Tuple{Number, Number, Number}, NTuple{4, Number}}

include("Mat.jl")

const InputArray = Union{AbstractArray{T, 3} where {T <: dtypes}, CxxMat}

include("mat_conversion.jl")
include("types_conversion.jl")

function cpp_to_julia(var)
    return var
end
function julia_to_cpp(var)
    return var
end

function cpp_to_julia(var::Tuple)
    ret_arr = Array{Any, 1}()
    for it in var
        push!(ret_arr, cpp_to_julia(it))
    end
    return tuple(ret_arr...)
end

function cpp_to_julia(var::CxxBool)
    return Bool(var)
end

function julia_to_cpp(var::Bool)
    return CxxBool(var)
end

include("cv_cxx_wrap.jl")

include("cv_manual_wrap.jl")


text:
using ITensors,
      Test

@testset "diagITensor (DiagBlockSparse)" begin

  @testset "diagITensor get and set elements" begin
    i = Index(QN(0)=>2,QN(1)=>3; tags="i")

    D = diagITensor(QN(),i,dag(i'))

    for b in eachnzblock(D)
      @test flux(D,b) == QN()
    end

    D[i(1),i'(1)] = 1
    D[i(2),i'(2)] = 2
    D[i(3),i'(3)] = 3
    D[i(4),i'(4)] = 4
    D[i(5),i'(5)] = 5

    @test_throws ErrorException D[i(1),i'(2)] = 2.0

    @test D[i(1),i'(1)] == 1
    @test D[i(2),i'(2)] == 2
    @test D[i(3),i'(3)] == 3
    @test D[i(4),i'(4)] == 4
    @test D[i(5),i'(5)] == 5
  end

  @testset "diagITensor Tuple constructor" begin
    i = Index(QN(0)=>2, QN(1)=>3; tags="i")

    D = diagITensor((i, dag(i')))

    for b in eachnzblock(D)
      @test flux(D, b) == QN()
    end
  end

  @testset "delta" begin
    i = Index(QN(0)=>2,QN(1)=>3; tags="i")
    ĩ = sim(i; tags="i_sim")
    j = Index(QN(0)=>2,QN(1)=>3,QN(2)=>4; tags="j")

    A = randomITensor(QN(), i, dag(j))

    δiĩ = δ(dag(i), ĩ)

    @test storage(δiĩ) isa NDTensors.DiagBlockSparse{ElT,
                                                   ElT} where {ElT<:Number}

    B = A * δiĩ

    A = permute(A, i, j)
    B = permute(B, ĩ, j)

    @test norm(dense(NDTensors.tensor(A)) -
               dense(NDTensors.tensor(B))) ≈ 0
  end

  @testset "delta Tuple constructor" begin
    i = Index(QN(0)=>2, QN(1)=>3; tags="i")
    ĩ = sim(i; tags="i_sim")

    δiĩ = δ((dag(i), ĩ))

    for b in eachnzblock(δiĩ)
      @test flux(δiĩ, b) == QN()
    end
  end

  @testset "Regression test for QN delta contraction bug" begin
    # http://itensor.org/support/2814/block-sparse-itensor-wrong-results-multiplying-delta-tensor
    s = Index([QN(("N",i,1))=>1 for i = 1:2])
    l = dag(addtags(s,"left"))
    r = addtags(s,"right")
    u = addtags(s,"up")
    d = dag(addtags(s,"down"))
    A = emptyITensor(l,r,u,d)
    A[1,1,1,1] = 1.0
    A[1,1,2,2] = 1.0
    A[2,2,1,1] = 1.0
    A[2,2,2,2] = 1.0
    δlr = δ(dag(l), dag(r))
    δud = δ(dag(u), dag(d))
    A1 = A * δlr
    denseA1 = dense(A) * dense(δlr)
    A2 = A1 * δud
    denseA2 = denseA1 * dense(δud)
    @test dense(A1) ≈ denseA1
    @test dense(A2) ≈ denseA2
    @test A2[] ≈ 4
  end

  @testset "Regression test for printing a QN Diag ITensor" begin
    # https://github.com/ITensor/NDTensors.jl/issues/61
    i = Index([QN()=>2])
    A = randomITensor(i', dag(i))
    U, S, V = svd(A, i')
    # Test printing S
    io = IOBuffer()
    show(io, S)
    sS = String(take!(io))
    @test sS isa String
    # Test printing U
    io = IOBuffer()
    show(io, U)
    sU = String(take!(io))
    @test sU isa String
  end
end

nothing


text:
v=:v≖Cross(:A,:B)
e1=[1,0,0];e2=[0,1,0];e3=[0,0,1]
A=:A≖3*:c*Vec(e2)
B=:B≖0.5*Vec(e3)
c=:c≖299792458
vs=v&[A,B,c]
@test vs.rhs==4.49688687e8*Vec(e1)

using LinearAlgebra
v1=Vec([1,2,3]);v2=Vec([3,2,1]);
@test dot(v1,v2)==10


text:
mutable struct Clock{T<:Number}
  time::T
  timevec::Array{T,1}
  iteration::Int
  Nt::Int
  Δt::T
  t_start::Real   # time when the simulation starts
end

function Clock(T::DataType,start_time::Number,iteration::Int,Nt::Int,Δt::Number)
    timevec = Array{T}(start_time:Δt:Δt*Nt)
    t_start = 0.0
    Clock{T}(start_time, timevec, iteration, Nt, Δt, t_start)
end

## Feedback functions
"""Converts time step into percent for feedback."""
function progress(i::Int,Nt::Int)
    if ((i+1)/Nt*100 % 1) < (i/Nt*100 % 1)  # update every 1 percent steps.
        percent = Int(round((i+1)/Nt*100))
        print("\r\u1b[K")   # remove previous percentage
        print("$percent%")  # print new one
    end
end

"""Returns a human readable string representing seconds in terms of days, hours, minutes or seconds."""
function readable_secs(secs::Real)
    days = Int(floor(secs/3600/24))
    hours = Int(floor((secs/3600) % 24))
    minutes = Int(floor((secs/60) % 60))
    seconds = Int(floor(secs%3600%60))

    if days > 0
        return "$(days)d, $(hours)h"
    elseif hours > 0
        return "$(hours)h, $(minutes)min"
    elseif minutes > 0
        return "$(minutes)min, $(seconds)s"
    else
        return "$(seconds)s"
    end
end

"""Estimates the total time the model integration will take."""
function duration_estimate(i::Int,Nt::Int,t_start::Real)
    time_per_step = (time()-t_start)/i
    time_total = Int(round(time_per_step*Nt))
    time_to_go = Int(round(time_per_step*(Nt-i)))

    s1 = "Model integration will take approximately "*readable_secs(time_total)*","
    s2 = "and is hopefully done on "*Dates.format(now() + Dates.Second(time_to_go),Dates.RFC1123Format)

    println(s1)     # print inline
    println(s2)
end

"""Feedback function that calls duration estimate, nan_detection and progress."""
function feedback(clock)

    i = clock.iteration
    Nt = clock.Nt

    if i == 1 # measure time after time loop executed once
        clock.t_start = time()
    elseif i == 100  # estimate time after 10 iterations
        duration_estimate(i,Nt,clock.t_start)
    end

    if i > 100      # show percentage only after duration is estimated
        progress(i,Nt)
    end

    if i == Nt
        s = " Integration done in "*readable_secs(time()-clock.t_start)*"."
        println(s)
    end

    return nothing
end


text:
export VoidActor, void

"""
    VoidActor{D}() where D

Void actor does nothing with input data, error and complete events, can be useful for debugging (e.g. to start side-effects with [`tap`](@ref) operator)

# Examples

```jldoctest
using Rocket

source = from([ 0, 1, 2 ])
actor  = VoidActor{Int}()

subscribe!(source, actor)
;

# output

```

See also: [`Actor`](@ref), [`void`](@ref), [`tap`](@ref)
"""
struct VoidActor{T} <: Actor{T} end

on_next!(actor::VoidActor{T}, data::T) where T = begin end
on_error!(actor::VoidActor, err)               = begin end
on_complete!(actor::VoidActor)                 = begin end

struct VoidActorFactory <: AbstractActorFactory end

create_actor(::Type{L}, factory::VoidActorFactory) where L = VoidActor{L}()

"""
    void()
    void(::Type{T}) where T

Creation operator for the `VoidActor` actor.

# Examples

```jldoctest
using Rocket

actor = void(Int)
actor isa VoidActor{Int}

# output
true

```

See also: [`VoidActor`](@ref), [`AbstractActor`](@ref)
"""
void()                  = VoidActorFactory()
void(::Type{T}) where T = VoidActor{T}()


text:
using StochasticDiffEq, Test, Random, DiffEqNoiseProcess
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear

Random.seed!(100)

for prob in (prob_sde_linear, prob_sde_2Dlinear)
  prob2 = remake(prob; tspan = (1.0, 0.0))

  solEM = solve(prob2, EM(), dt = -1//2^(4), tstops = [0.33])
  @test solEM.t[15] > solEM.t[16]

  solSRIW1 = solve(prob2, SRIW1(), tstops = [0.33])
  @test solSRIW1.t[15] > solSRIW1.t[16]
end

# Test reversing the noise

α=1
β=1
u₀=1/2
f(u,p,t) = α*u
g(u,p,t) = β*u
dt = 1//2^(4)
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,(0.0,1.0))
sol =solve(prob,EulerHeun(),dt=0.01,save_noise=true)
_sol = deepcopy(sol) # to make sure the plot is correct
W3 = NoiseGrid(reverse!(_sol.t),reverse!(_sol.W))
prob3 = SDEProblem(f,g,sol[end],(1.0,0.0),noise=W3)
sol2 = solve(prob3,EulerHeun(),dt=0.01)
@test sol.u ≈ reverse!(sol2.u) atol=1e-1


text:
#   Author: V. Vitale
#   Feb 2022
#   Conctraction utilities

## initial E and F matrices for the left and right vacuum states
function initial_L(W::Array)
    sW=size(W)
    L = ones(1,sW[1],1)#ones(sW[1],1,1)
    return L
end

function initial_R(W::Array)
    sW=size(W)
    R = ones(1,sW[2],1)#ones(sW[2],1,1)
    return R
end

function construct_R(A::MPS, W::MPO)
    R = Dict()
    R[A.N] = initial_R(W.data[A.N])
    for i in A.N:-1:1
        R[i-1] = contract_from_right(R[i], A.data[i], W.data[i])
    end
    return R
end

function construct_L(A::MPS, W::MPO)
    L = Dict()
    L[1] = initial_L(W.data[1])
    return L
end

## tensor contraction from the left hand side
## +-    +--X-
## |     |  |
## L' =  L--W-
## |     |  |
## +-    +--X-

function contract_from_left(E_L::Array,X::Array,W::Array)
    @tensor temp_1[:] := E_L[1,-3,-4]  * X[1,-2,-1]
    @tensor temp_2[:] := temp_1[-1, -2, -5, 1 ]* conj( X[1, -4,-3] )
    @tensor temp[:] := temp_2[ -1, 2, -3, 3, 1 ] * W[1,-2, 2, 3]
    return temp
end

## tensor contraction from the right hand side
##  -+     -Y--+
##   |      |  |
##  -R' =  -W--R
##   |      |  |
##  -+     -Y--+
function contract_from_right(E_R::Array,Y::Array,W::Array)
    @tensor temp_1[:] :=  Y[-1,-2,1]  * E_R[ 1, -3,-4]
    @tensor temp_2[:] := temp_1[ -1, -2, -5,  1] * conj( Y[-3,-4,1] )
    @tensor temp[:] := temp_2[ -1, 1, -3, 2, 6] * W[ -2 6 1 2 ]
    return temp
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
const _HEAD_HASH_OK_REGEX = Regex("(?<hash>[0-9a-f]{40})")
function _HEAD_hash(repodir)
    !_check_gitdir(repodir) && return ""
    out = _run("git -C $(repodir) rev-parse HEAD 2>&1"; verbose = false, ignorestatus = true)
    m = match(_HEAD_HASH_OK_REGEX, out)
    isnothing(m) ? "" : m[:hash]
end
_HEAD_hash(gl::GitLink) = _HEAD_hash(repo_dir(gl))


text:
export ContinuousSpace

struct ContinuousSpace{T,R} <: Space{PiecewiseSegment{T},R}
    domain::PiecewiseSegment{T}
end

ContinuousSpace(d::PiecewiseSegment{T}) where {T} =
    ContinuousSpace{T,real(eltype(T))}(d)

Space(d::PiecewiseSegment) = ContinuousSpace(d)

isperiodic(C::ContinuousSpace) = isperiodic(domain(C))

spacescompatible(a::ContinuousSpace,b::ContinuousSpace) = domainscompatible(a,b)
conversion_rule(a::ContinuousSpace,
                b::PiecewiseSpace{CD,DD,RR}) where {CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},DD,RR<:Real} where {DDD,RRR} = a

plan_transform(sp::ContinuousSpace,vals::AbstractVector) =
    TransformPlan{eltype(vals),typeof(sp),false,Nothing}(sp,nothing)

function *(P::TransformPlan{T,SS,false},vals::AbstractVector{T}) where {T,SS<:ContinuousSpace}
    S = P.space
    n=length(vals)
    d=domain(S)
    K=ncomponents(d)
    k=div(n,K)

    PT=promote_type(prectype(d),eltype(vals))
    if k==0
        vals
    elseif isperiodic(d)
        ret=Array{PT}(undef, max(K,n-K))
        r=n-K*k

        for j=1:r
            cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)),
                          vals[(j-1)*(k+1)+1:j*(k+1)])
            if j==1
                ret[1]=cfs[1]-cfs[2]
                ret[2]=cfs[1]+cfs[2]
            elseif j < K
                ret[j+1]=cfs[1]+cfs[2]
            end
            ret[K+j:K:end]=cfs[3:end]
        end

        for j=r+1:K
            cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)),
                          vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k])
            if length(cfs)==1 && j <K
                ret[j+1]=cfs[1]
            elseif j==1
                ret[1]=cfs[1]-cfs[2]
                ret[2]=cfs[1]+cfs[2]
            elseif j < K
                ret[j+1]=cfs[1]+cfs[2]
            end
            ret[K+j:K:end]=cfs[3:end]
        end

        ret
    else
        ret=Array{PT}(undef, n-K+1)
        r=n-K*k

        for j=1:r
            cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)),
                          vals[(j-1)*(k+1)+1:j*(k+1)])
            if j==1
                ret[1]=cfs[1]-cfs[2]
            end

            ret[j+1]=cfs[1]+cfs[2]
            ret[K+j+1:K:end]=cfs[3:end]
        end

        for j=r+1:K
            cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)),
                          vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k])

            if length(cfs) ≤ 1
                ret .= cfs
            else
                if j==1
                    ret[1]=cfs[1]-cfs[2]
                end
                ret[j+1]=cfs[1]+cfs[2]
                ret[K+j+1:K:end]=cfs[3:end]
            end
        end

        ret
    end
end

components(S::ContinuousSpace) = map(ChebyshevDirichlet{1,1},components(domain(S)))
canonicalspace(S::ContinuousSpace) = PiecewiseSpace(components(S))
convert(::Type{PiecewiseSpace}, S::ContinuousSpace) = canonicalspace(S)

blocklengths(C::ContinuousSpace) = Fill(ncomponents(C.domain),∞)

block(C::ContinuousSpace,k) = Block((k-1)÷ncomponents(C.domain)+1)

## components

components(f::Fun{CS},j::Integer) where {CS<:ContinuousSpace} = components(Fun(f,canonicalspace(f)),j)
components(f::Fun{CS}) where {CS<:ContinuousSpace} = components(Fun(f,canonicalspace(space(f))))

function points(f::Fun{CS}) where {CS<:ContinuousSpace}
    n=ncoefficients(f)
    d=domain(f)
    K=ncomponents(d)

    m=isperiodic(d) ? max(K,n+2K-1) : n+K
    points(f.space,m)
end

## Conversion

coefficients(cfsin::AbstractVector,A::ContinuousSpace,B::PiecewiseSpace) =
    defaultcoefficients(cfsin,A,B)

coefficients(cfsin::AbstractVector,A::PiecewiseSpace,B::ContinuousSpace) =
    default_Fun(Fun(A,cfsin),B).coefficients

coefficients(cfsin::AbstractVector,A::ContinuousSpace,B::ContinuousSpace) =
    default_Fun(Fun(A,cfsin),B).coefficients

# We implemnt conversion between continuous space and PiecewiseSpace with Chebyshev dirichlet
function Conversion(ps::PiecewiseSpace{CD,DD,RR},cs::ContinuousSpace) where {CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},
                                                                    DD,RR<:Real} where {DDD,RRR}
    @assert ps == canonicalspace(cs)
    ConcreteConversion(ps,cs)
end

function Conversion(cs::ContinuousSpace,ps::PiecewiseSpace{CD,DD,RR}) where {CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},
                                                                    DD,RR<:Real} where {DDD,RRR}
    @assert ps == canonicalspace(cs)
    ConcreteConversion(cs,ps)
end

bandwidths(C::ConcreteConversion{PiecewiseSpace{CD,DD,RR},CS}) where {CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},
                                                                                DD,RR<:Real,CS<:ContinuousSpace} where {DDD,RRR} =
    1,ncomponents(domain(rangespace(C)))

function getindex(C::ConcreteConversion{PiecewiseSpace{CD,DD,RR},CS,T},
                  k::Integer,j::Integer) where {T,DD,RR<:Real,CS<:ContinuousSpace,
                                                CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}}}  where {DDD,RRR}
    d=domain(rangespace(C))
    K=ncomponents(d)
    if isperiodic(d)
        if k==j==1
            one(T)
        elseif k==1 && j==K+1
            -one(T)
        elseif 2≤k≤K && (j==k-1 || j==K+k-1)
            one(T)
        elseif K<k && j==k+K
            one(T)
        else
            zero(T)
        end
    else
        if k==j==1
            one(T)
        elseif k==1 && j==K+1
            -one(T)
        elseif 2≤k≤K+1 && (j==k-1 || j==K+k-1)
            one(T)
        elseif K+1<k && j==k+K-1
            one(T)
        else
            zero(T)
        end
    end
end

bandwidths(C::ConcreteConversion{<:ContinuousSpace,
                                     PiecewiseSpace{CD,DD,RR}}) where {CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},
                                              DD,RR<:Real}  where {DDD,RRR} =
            isperiodic(domainspace(C)) ? (2ncomponents(domain(rangespace(C)))-1,1) :
                                         (ncomponents(domain(rangespace(C))),1)

function getindex(C::ConcreteConversion{<:ContinuousSpace,
                                            PiecewiseSpace{CD,DD,RR},T},
                      k::Integer,j::Integer) where {T,CD<:Tuple{Vararg{ChebyshevDirichlet{1,1,DDD,RRR}}},
                                        DD,RR<:Real} where {DDD,RRR}
    d=domain(domainspace(C))
    K=ncomponents(d)
    if isperiodic(d)
        if k < K && (j==k || j==k+1)
            one(T)/2
        elseif k==K && (j==K || j==1)
            one(T)/2
        elseif K+1≤k≤2K && j==k-K
            -one(T)/2
        elseif K+1≤k<2K && j==k-K+1
            one(T)/2
        elseif k==2K && j==1
            one(T)/2
        elseif k>2K && j==k-K
            one(T)
        else
            zero(T)
        end
    else
        if k≤K && (j==k || j==k+1)
            one(T)/2
        elseif K+1≤k≤2K && j==k-K
            -one(T)/2
        elseif K+1≤k≤2K && j==k-K+1
            one(T)/2
        elseif k>2K && j==k-K+1
            one(T)
        else
            zero(T)
        end
    end
end

# Dirichlet for Squares

const TensorChebyshevDirichlet = TensorSpace{<:Tuple{<:ChebyshevDirichlet{1,1,<:IntervalOrSegment},
                                  <:ChebyshevDirichlet{1,1,<:IntervalOrSegment}}}

Dirichlet(S::TensorChebyshevDirichlet,k) = k == 0 ? ConcreteDirichlet(S,0) : tensor_Dirichlet(S,k)

Dirichlet(d::RectDomain) =
    Dirichlet(ChebyshevDirichlet{1,1}(factor(d,1))*ChebyshevDirichlet{1,1}(factor(d,2)))

isblockbanded(::Dirichlet{<:TensorChebyshevDirichlet}) = true

blockbandwidths(::Dirichlet{<:TensorChebyshevDirichlet}) = (0,2)

colstop(B::Dirichlet{<:TensorChebyshevDirichlet}, j::Integer) = j ≤ 3 ? 4 : 4(block(domainspace(B),j).n[1]-1)

function getindex(B::ConcreteDirichlet{<:TensorChebyshevDirichlet}, k::Integer,j::Integer)
    T = eltype(B)
    ds = domainspace(B)
    rs = rangespace(B)
    if j == 1 && k ≤ 4
        one(T)
    elseif j == 2 && k ≤ 2
        -one(T)
    elseif j == 2 && k ≤ 4
        one(T)
    elseif j == 3 && (k == 1 || k == 4)
        -one(T)
    elseif j == 3 && (k == 2 || k == 3)
        one(T)
    elseif j == 5 && (k == 2 || k == 4)
        -one(T)
    elseif j == 5 && (k == 1 || k == 3)
        one(T)
    elseif j == 5 || j ≤ 3
        zero(T)
    else
        K = Int(block(rs,k))
        J = Int(block(ds,j))
        m = mod(k-1,4)
        s,t =  blockstart(ds,J),  blockstop(ds,J)
        if K == J-1 && (m == 1  && j == s ||
                       (m == 0  && j == t))
            one(T)
        elseif K == J-1 && ((m == 3 && j == s) ||
                            (m == 2 && j == t))
            iseven(K) ? one(T) : -one(T)
        elseif K == J-2 && m == 1 && j == s+1
            one(T)
        elseif K == J-2 && m == 2 && j == t-1
            iseven(K) ? one(T) : -one(T)
        elseif K == J-2 && m == 0 && j == t-1
            -one(T)
        elseif K == J-2 && m == 3 && j == s+1
            iseven(K) ? -one(T) : one(T)
        else
            zero(T)
        end
    end
end

function BlockBandedMatrix(S::SubOperator{T,<:ConcreteDirichlet{<:TensorChebyshevDirichlet},
                                Tuple{UnitRange{Int},UnitRange{Int}}}) where {T}
    P=parent(S)
    ret=BlockBandedMatrix(Zeros, S)
    kr,jr=parentindices(S)

    K1=block(rangespace(P),kr[1])
    Kr1=blockstart(rangespace(P),K1)
    J1=block(domainspace(P),jr[1])
    Jr1=blockstart(domainspace(P),J1)

    if ret.rows[1]>0 && ret.cols[1]==1
        view(ret,Block(1),Block(1))[:,1] = 1
    end

    if ret.rows[1] > 0 && ret.cols[2] > 0
        B=view(ret,Block(1),Block(2))

        k_sh = kr[1]-1; j_sh = max(jr[1]-2,0)
        if j_sh == 0
            # first column
            k_sh == 0 && (B[1,1-j_sh]=-1)
            k_sh ≤  1 && (B[2-k_sh,1-j_sh]=-1)
            k_sh ≤ 2 && (B[3-k_sh,1-j_sh]=1)
            B[4-k_sh,1-j_sh]=1
        end
        # second column
        k_sh == 0 && (B[1-k_sh,2-j_sh]=-1)
        k_sh ≤ 1  && (B[2-k_sh,2-j_sh]=1)
        k_sh ≤ 2  && (B[3-k_sh,2-j_sh]=1)

        B[4-k_sh,2-j_sh]=-1
    end

    if ret.rows[1] > 0 && ret.cols[3] > 0
        B=view(ret,Block(1),Block(3))

        k_sh = kr[1]-1; j_sh = max(jr[1]-4,0)
        # second column
        k_sh == 0 && (B[1-k_sh,2-j_sh]=1)
        k_sh ≤ 1 && (B[2-k_sh,2-j_sh]=-1)
        k_sh ≤ 2 && (B[3-k_sh,2-j_sh]=1)
        B[4-k_sh,2-j_sh]=-1
    end
    for K=Block(2):2:Block(min(length(ret.rows),length(ret.cols)-1))
        J = K+1  # super-diagonal block
        N = ret.rows[K.n[1]]
        M = ret.cols[J.n[1]]
        if N ≠ 0 && M ≠ 0
            # calculate shift
            k_sh = K == K1 ? kr[1]-Kr1 : 0
            j_sh = J == J1 ? jr[1]-Jr1 : 0
            B = view(ret,K,J)

            1 ≤ 2-k_sh ≤ N && j_sh == 0 && (B[2-k_sh,1-j_sh]=1)
            1 ≤ 4-k_sh ≤ N && j_sh == 0 && (B[4-k_sh,1-j_sh]=1)
            k_sh == 0 && 1 ≤ J.n[1]-j_sh ≤ M && (B[1-k_sh,J.n[1]-j_sh]=1)
            k_sh ≤ 2 &&  1 ≤ J.n[1]-j_sh ≤ M && (B[3-k_sh,J.n[1]-j_sh]=1)
        end
    end
    for K=Block(3):2:Block(min(length(ret.rows),length(ret.cols)-1))
        J = K+1  # super-diagonal block
        N = ret.rows[K.n[1]]
        M = ret.cols[J.n[1]]
        if N ≠ 0 && M ≠ 0
            # calculate shift
            k_sh = K == K1 ? kr[1]-Kr1 : 0
            j_sh = J == J1 ? jr[1]-Jr1 : 0
            B = view(ret,K,J)

            1 ≤ 2-k_sh ≤ N && j_sh == 0 && (B[2-k_sh,1-j_sh]=1)
            1 ≤ 4-k_sh ≤ N && j_sh == 0 && (B[4-k_sh,1-j_sh]=-1)
            k_sh == 0 && 1 ≤ J.n[1]-j_sh ≤ M && (B[1-k_sh,J.n[1]-j_sh]=1)
            1 ≤ 3-k_sh ≤ N &&  1 ≤ J.n[1]-j_sh ≤ M && (B[3-k_sh,J.n[1]-j_sh]=-1)
        end
    end
    for K=Block(2):2:Block(min(length(ret.rows),length(ret.cols)-2))
        J = K+2  # super-diagonal block
        N = ret.rows[K.n[1]]
        M = ret.cols[J.n[1]]

        if N ≠ 0 && M ≠ 0
            B=view(ret,K,J)
            # calculate shift
            k_sh = K == K1 ? kr[1]-Kr1 : 0
            j_sh = J == J1 ? jr[1]-Jr1 : 0
            B = view(ret,K,J)

            1 ≤ 2-k_sh ≤ N && 1 ≤ 2-j_sh ≤ M && (B[2-k_sh,2-j_sh]=1)
            1 ≤ 4-k_sh ≤ N && 1 ≤ 2-j_sh ≤ M && (B[4-k_sh,2-j_sh]=-1)
            k_sh == 0 && 1 ≤ J.n[1]-j_sh-1 ≤ M && (B[1,J.n[1]-j_sh-1]=-1)
            1 ≤ 3-k_sh ≤ N &&  1 ≤ J.n[1]-j_sh-1 ≤ M && (B[3-k_sh,J.n[1]-j_sh-1]=1)
        end
    end
    for K=Block(3):2:Block(min(length(ret.rows),length(ret.cols)-2))
        J = K+2
        B=view(ret,K,J)
        N = ret.rows[K.n[1]]
        M = ret.cols[J.n[1]]

        if N ≠ 0 && M ≠ 0
            B=view(ret,K,J)
            # calculate shift
            k_sh = K == K1 ? kr[1]-Kr1 : 0
            j_sh = J == J1 ? jr[1]-Jr1 : 0
            B = view(ret,K,J)

            1 ≤ 2-k_sh ≤ N && 1 ≤ 2-j_sh ≤ M && (B[2-k_sh,2-j_sh]=1)
            1 ≤ 4-k_sh ≤ N && 1 ≤ 2-j_sh ≤ M && (B[4-k_sh,2-j_sh]=1)
            k_sh == 0 && 1 ≤ J.n[1]-j_sh-1 ≤ M && (B[1,J.n[1]-j_sh-1]=-1)
            1 ≤ 3-k_sh ≤ N &&  1 ≤ J.n[1]-j_sh-1 ≤ M && (B[3-k_sh,J.n[1]-j_sh-1]=-1)
        end
    end

    ret
end

union_rule(A::PiecewiseSpace, B::ContinuousSpace) = union(A, convert(PiecewiseSpace, B))
union_rule(A::ConstantSpace, B::ContinuousSpace) = B
union_rule(A::ContinuousSpace, B::PolynomialSpace{<:IntervalOrSegment}) =
    Space(domain(A) ∪ domain(B))

function approx_union(a::AbstractVector{T}, b::AbstractVector{V}) where {T,V}
    ret = sort!(union(a,b))
    for k in length(ret)-1:-1:1
        isapprox(ret[k] , ret[k+1]; atol=10eps()) && deleteat!(ret, k+1)
    end
    ret
end

function union_rule(A::ContinuousSpace{<:Real}, B::ContinuousSpace{<:Real})
    p_A,p_B = domain(A).points, domain(B).points
    a,b = minimum(p_A),  maximum(p_A)
    c,d = minimum(p_B),  maximum(p_B)
    @assert !isempty((a..b) ∩ (c..d))
    ContinuousSpace(PiecewiseSegment(approx_union(p_A, p_B)))
end

function integrate(f::Fun{<:ContinuousSpace})
    cs = cumsum.(components(f))
    for k=1:length(cs)-1
        cs[k+1] += last(cs[k])
    end
    Fun(Fun(cs, PiecewiseSpace), space(f))
end

cumsum(f::Fun{<:ContinuousSpace}) = integrate(f)


text:
using MovingWeightedLeastSquares

xs = [0.1 * i for i in 0:100]
fs = hcat([sin(x) for x in xs], [cos(x) for x in xs])
println(fs)

m = mwls(xs, fs, 0.5, (d, e) -> exp(-d^2))

println(m(1))


text:
const INPUT_DATA_CEC2020 = BSON.load(joinpath(@__DIR__, "data.bson"))

const NAME_OF_PROBLEMS_CEC2020 = [
                                      "cec2020_f1", "cec2020_f2",
                                      "cec2020_f3", "cec2020_f4",
                                      "cec2020_f5", "cec2020_f6",
                                      "cec2020_f7", "cec2020_f8",
                                      "cec2020_f9", "cec2020_f10"
                                     ]

function load_data_from_dict(fnum, D)
    fnums = [1,2,3,7,4,16,6,22,24,25]
    # dimensions = [2,5,10,15,20,30,50,100]

    fnum_= fnums[fnum]

    if D==2&&(fnum_==4||fnum_==16||fnum_==6)
        error("function $fnum not implemented for D = $D.")
    end
    fnum = fnum_

    str_fnum = string(fnum)
    str_D = string(D)
    o = INPUT_DATA_CEC2020["o_fnum"*str_fnum*"_D"*str_D]
    M = INPUT_DATA_CEC2020["M_fnum"*str_fnum*"_D"*str_D]
    S = INPUT_DATA_CEC2020["S_fnum"*str_fnum*"_D"*str_D]

    return o, M, S
end

get_cec2020_minimum(fnum) = [100,1100,700,1900,1700,1600,2100,2200,2400,2500][fnum]

"""
    get_cec2020_config(fnum::Int; n = 10)

Return a dictionary containing the information about the problem in CEC2020 benchmark.
`n` is number of variables (2,5,10,15,20,30,50,100)
"""
function get_cec2020_config(fnum::Int; n = 10)
    @assert 1 <= fnum <= 10
    @assert n in [2,5,10,15,20,30,50,100]

    return Dict(:function => NAME_OF_PROBLEMS_CEC2020[fnum],
                :fn    => 1,
                :n     => n,
                :minimum => get_cec2020_minimum(fnum),
                :gn    => 0,
                :hn    => 0,
                :xmin  => fill(-100.0, n),
                :xmax  => fill(100.0, n))
end

"""
    get_cec2020_config(fnum::Int; n = 10)

Return a dictionary containing the information about the problem in CEC2020 benchmark.
`n` is number of variables (2,5,10,15,20,30,50,100)
"""
function get_cec2020_problem(fnum::Int; n = 10)

    !(1 <= fnum <= 10) && error("Problem $fnum is not implemented yet or does not exist.")

    # objetive function x -> f(x)
    f = eval(Symbol(NAME_OF_PROBLEMS_CEC2020[fnum]))

    # problem configuration
    conf = get_cec2020_config(fnum; n = n)

    return f, conf
end