text:
if Base.HOME_PROJECT[] !== nothing
    # JuliaLang/julia/pull/28625
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, PDFIO

makedocs(
    format = Documenter.HTML(),
    modules = [PDFIO],
    sitename = "PDFIO",
    pages = [
        "README.md",
        "Architecture and Design" => "arch.md",
        "Encryption in PDF" => "encrypt.md",
        "Digital Signatures" => "digsig.md",
        "API Reference" => "index.md",
        "PDFIO License" => "LICENSE.md"
    ]
)

deploydocs(
    repo   = "github.com/sambitdash/PDFIO.jl.git",
)


text:
function SeisDelay(d; delay=0.1, dt=0.001)

    nt = size(d,1)
    nf = nt
    D = fft(d,1);
    dw = 2*pi/nt/dt
    nw = round(Int,nt/2) + 1
    for iw = 1 : nw
	w = (iw-1)*dw
	D[iw,:] = D[iw,:]*exp(-1im*w*delay)
    end
    # symmetries
    for iw=nw+1:nf
	D[iw,:] = conj(D[nf-iw+2,:])
    end
    d2 = ifft(D,1)
    d2 = real(d2[1:nt,:])
    return d2

end

ext: jl


text:
#=
!NOTE:
1. This sections defines abstract types for structs/functions that need to be defined/overloaded so custom
algorithms(algorithms.jl), kernels inside an Baytes package (kernel.jl) can be used.
2. Super types defined before files are usually arguments in algorithm/kernel.
3. Functions shown in both, algorithms.jl and kernel.jl, need to be defined in either case.
=#
############################################################################################
# Define Abstract Super Types - so we do not need to load ModelWrappers for abstraction and function definitions
abstract type AbstractModelWrapper end
abstract type AbstractObjective end
abstract type AbstractResult end

############################################################################################
# Include files
include("algorithms.jl")
include("kernels.jl")
include("utility.jl")

############################################################################################
# Export
export
    AbstractModelWrapper,
    AbstractObjective,
    AbstractResult


text:
using WeightedLS
using LinearAlgebra
using Test

@testset "normal LS" begin
    A = Float64[1 0; 1 -2]
    B = Float64[32; -4]
    x = WeightedLS.normalLS(A,B)

    @test isapprox(norm(x-Float64[32.0; 18.0]), 0.0; atol = 1e-5)

    A = Float64[10 12;20 22;21 0]
    B = Float64[1;2;56]
    x = WeightedLS.normalLS(A,B)
    @test isapprox(norm(x-Float64[2.662964939354169; -2.285446442736200]), 0.0; atol = 1e-5)
end

@testset "weighted LS" begin

    a = Float64[1; 1; 1; 1; 1]
    m = size(a,1)
    W = Float64[1; Inf; Inf; 1; 1]

    @test WeightedLS.checkConstraint(a,W,m) == 2

    a = Float64[1; 0.5; 1; 1; 1]
    m = size(a,1)
    W = Float64[1; Inf; Inf; 1; 1]

    @test WeightedLS.checkConstraint(a,W,m) == 3

    a = Float64[1; 1; 1; 1; 1]
    m = size(a,1)
    W = Float64[1; 1; 1; 1; 1]

    @test WeightedLS.checkConstraint(a,W,m) == -1

    A = Float64[1 0 0 0 0 1;
                0 0 0 0 1 0;
                0 0 1 1 0 0;
                0 1 0 0 0 0;
                0 0 0 0 0 1]
    B = Float64[0;0;0;0;0]
    W = Float64[0;1;0;1;1]
    W = 1 ./W
    # println(W)
    Rd, outweights = WeightedLS.weightedLS(A,B,W)
    # println(outweights)
    # show(stdout, "text/plain", Rd)

    expected_Rd = Float64[1  0  0  0  0  1  0;
                          0  1  0  0  0  0  0;
                          0  0  1  1  0  0  0;
                          0  0  0  0  1  0  0;
                          0  0  0  0  0  1  0;]
    @test isapprox(norm(Rd-expected_Rd), 0.0; atol = 1e-5)

    A = Float64[    1 0 0;
                    0 1 0;
                    0 0 1;
                   -1 0 1;
                    1 0 0;
                    0 1 0;
                    0 0 1;
                    0 -1 1;
                    1 0 0;
                    0 1 0;
                    0 0 1]
    B = zeros(11,1)
    W = ones(11,1)
    W[4] = Inf; W[8] = Inf;
    # println(W)
    Rd, outweights = WeightedLS.weightedLS(A,dropdims(B,dims=2),dropdims(W,dims=2))
    # println(outweights)
    # show(stdout, "text/plain", Rd)

    expected_Rd = Float64[-1  0 1 0;
                           0 -1 1 0;
                           0  0 1 0;]
    @test isapprox(norm(Rd[1:3,:]-expected_Rd), 0.0; atol = 1e-5)

    A = Float64[1 0; 1 -2]
    B = Float64[32; -4]
    W = Float64[1; Inf]
    Rd, outweights = WeightedLS.weightedLS(A,B,W)
    R = Rd[1:2,1:2]
    d = Rd[1:2,3]
    # println(R)
    # println(d)
    x = R\d

    @test isapprox(norm(x-Float64[32.0; 18.0]), 0.0; atol = 1e-5)
end


text: isarmstrong(n::Int) = n == sum(digits(n).^length(digits(n)))

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
function compute_offset(p::Integer, segment_lo::Integer)
    p² = p * p
    if p² >= segment_lo
        byte = cld(p², 30)

        # Wheel index stores both the index of the prime number mod 30
        # and the index of the active multiple. We start crossing off
        # p * p, so that would be to_idx(p % 30) twice. We combine those values
        # as a number between 1 ... 64.
        wheel = to_idx(p % 30)
        wheel_index = 8 * (wheel - 1) + wheel

        return byte, wheel_index
    else
        # p * q will be the first number to cross off
        q = cld(segment_lo, p)
        q_quot, q_rem = divrem(q, 30)

        remainders = (1, 7, 11, 13, 17, 19, 23, 29, 31)
        i = 1
        while remainders[i] < q_rem
            i += 1
        end

        # maybe wrap around
        q_rem = i == 9 ? 1 : remainders[i]
        q_quot = i == 9 ? q_quot + 1 : q_quot
        i = i == 9 ? 1 : i

        # Our actual first acceptable multiple
        q = 30q_quot + q_rem

        byte = cld(p * q, 30)
        wheel_index = 8 * (to_idx(p % 30) - 1) + i

        return byte, wheel_index
    end
end

struct Siever
    prime_div_30::Int

    # byte_index is the integer range 30(byte_index - 1) up to 30byte_index - 1
    byte_index::Int

    # Stores the next prime number to be crossed off.
    # If `p` is the prime number and `q` the next multiple to be stored
    # Wheel index 8 * to_idx(p % 30) * to_idx(q % 30)
    wheel_index::Int

    function Siever(p::Int, segment_lo::Int)
        byte, wheel = compute_offset(p, segment_lo)
        return new(p ÷ 30, byte, wheel)
    end

    Siever(prime_div_30, byte_index, wheel_index) = new(prime_div_30, byte_index, wheel_index)
end

function Base.show(io::IO, p::Siever)
    print(io, 30p.prime_div_30 + ps[(p.wheel_index - 1) ÷ 8 + 1], " (", p.byte_index, ", ", p.wheel_index, ")")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 332e4ec2-b1c5-47c0-a12e-db0195023e18
# hideall

using Genie, Genie.Router

# ╔═╡ 69bec8f4-56b1-48ee-aa2c-ba4310e79568
# hideall

using Genie.Assets;

# ╔═╡ 32df4cc8-39ac-11ec-2f5e-b176daa6db9a
md"""
# Working with Web Sockets

Genie provides a powerful workflow for client-server communication over websockets. The system hides away the complexity of the network level communication, exposing powerful abstractions which resemble Genie's familiar MVC workflow: the clients and the server exchange messages over `channels` (which are the equivalent of `routes`).

## Registering `channels`

The messages are mapped to a matching channel, where are processed by Genie's `Router` which extracts the payload and invokes the designated handler (controller method or function). For most purposes, the `channels` are the functional equivalents of `routes` and are defined in the same way:
"""

# ╔═╡ a54bd0a5-7ec6-46db-8d42-a4cfc3eb1c2f
md"""
If you are using Julia REPL and Pluto notebook. Execute the below lines:

```julia
julia> import Pkg
julia> Pkg.add("Genie")
```
"""

# ╔═╡ c7ea1c3e-2a61-41d4-ba37-21027f56453a
md"""

```julia
julia> using Genie, Genie.Router

julia> channel("/foo/bar") do
         # process request
       end

```

"""

# ╔═╡ a8b81c5a-6a04-4ef6-8522-82582da4f78a
# hideall

channel("/foo/bar") do
	#process request
end

# ╔═╡ 72010195-6acc-4a4c-bb58-5331f8007192
md"""

```julia
julia> module YourController
        function your_handler
		   #process request
        end
       end
```
"""

# ╔═╡ dc4430bf-4555-4b3a-9c9c-9c14ef4ac806
# hideall

module YourController
	function your_handler
	 	#process request
	end
end

# ╔═╡ bd80dd1c-6e99-4c60-91d3-e470b9f5020c
md"""

```julia
julia> channel("/baz/bax", YourController.your_handler)
```
"""

# ╔═╡ a6ec9f7f-c574-4af0-a755-5d090a00d423
# hideall

channel("/baz/bax", YourController.your_handler)

# ╔═╡ b65c4644-f5f8-4a68-8df7-8064b1274ba5
md"""
The above `channel` definitions will handle websocket messages sent to `/foo/bar` and `/baz/bax`.
"""

# ╔═╡ f082d891-61ba-4a9f-b1ea-943ed2aba815
md"""
## Setting up the client

In order to enable WebSockets communication in the browser we need to load a JavaScript file. This is provided by Genie, through the `Assets` module. Genie makes it extremely easy to setup the WebSockets infrastructure on the client side, by providing the `Assets.channels_support()` method. For instance, if we want to add support for WebSockets to the root page of a web app, all we need is this:
"""

# ╔═╡ 3911bc31-f4ea-493e-95a9-9d743670c4c6
md"""

```julia
julia> using Genie, Genie.Router, Genie.Assets
```
"""

# ╔═╡ 266ce7f3-4364-44ec-90b8-61e06dba5d27
md"""

```julia
julia> route("/") do
		 Assets.channels_support()
	   end
```
"""

# ╔═╡ 70db282b-ec1c-40f7-bd7c-59ed2e94891c
# hideall

route("/") do
    Assets.channels_support()
end

# ╔═╡ 45c63333-2726-400e-a179-996ce88f456a
md"""
Literally, that is all we need in order to be able to push and receive messages between client and server.
"""

# ╔═╡ 15c148af-9708-4cec-a390-17a87a02a8fb
md"""
---

## Try it!

You can follow through by running the following Julia code in a Julia REPL:

```julia
using Genie, Genie.Router, Genie.Assets

Genie.config.websockets_server = true # enable the websockets server

route("/") do
  Assets.channels_support()
end

up() # start the servers
```
"""

# ╔═╡ 78550cf1-a2ef-487e-87a1-589fe3c7c048
# hideall

begin
	Genie.config.websockets_server = true; # enable the websockets server

	route("/") do
  		Assets.channels_support();
	end;

	up();
end;

# ╔═╡ 5350a2e0-2f95-4fb6-a9e5-18b4fc58d84f
md"""

to shutdown server

```julia
julia> down();
```
"""

# ╔═╡ 52576077-67bd-4a24-b81b-be6f475d4192
# hideall

down();

# ╔═╡ e8137d45-a1c8-480c-80b3-05454f3d241b
md"""
Now if you visit <http://localhost:8000> you'll get a blank page -- which, however, includes all the necessary functionality for WebSockets communication! If you use the browser's developer tools, the Network pane will indicate that a `channels.js` file was loaded and that a WebSockets request was made (Status 101 over GET). Additionally, if you peek at the Console, you will see a `Subscription ready` message. The output in the console should be something like:
"""

# ╔═╡ a26be9e1-42ce-43d6-8c36-d2a08ba75f99
md"""

```text
Subscription ready
channels.js:133 Overwrite window.parse_payload to handle messages from the server
channels.js:134 Subscription: OK
```
"""

# ╔═╡ 0f99c27a-29df-4aad-a387-12799664569f
md"""
**What happened?**

At this point, by invoking `Assets.channels_support()`, Genie has done the following:

* loaded the bundled `channels.js` file which provides a JS API for communicating over WebSockets
* has created two default channels, for subscribing and unsubscribing: `/__/subscribe` and `/__/unsubscribe`
* has invoked `/__/subscribe` and created a WebSockets connection between client and server

### Pushing messages from the server

We are ready to interact with the client. Go to the Julia REPL running the web app and run:
"""

# ╔═╡ fabc9669-ce12-470f-b551-3d8f69cb221b
md"""

```julia
julia> Genie.WebChannels.connected_clients()
1-element Array{Genie.WebChannels.ChannelClient,1}:
 Genie.WebChannels.ChannelClient(HTTP.WebSockets.WebSocket{HTTP.ConnectionPool.Transaction{Sockets.TCPSocket}}(T0  🔁    0↑🔒    0↓🔒 100s 127.0.0.1:8001:8001 ≣16, 0x01, true, UInt8[0x7b, 0x22, 0x63, 0x68, 0x61, 0x6e, 0x6e, 0x65, 0x6c, 0x22  …  0x79, 0x6c, 0x6f, 0x61, 0x64, 0x22, 0x3a, 0x7b, 0x7d, 0x7d], UInt8[], false, false), ["__"])
```
"""

# ╔═╡ 79660882-e6d9-4125-8601-277888c199a5
# hideall

Genie.WebChannels.connected_clients();

# ╔═╡ 229008fe-bad8-455b-853f-57eb39e87665
md"""
We have one connected client to the `__` channel! We can send it a message:
"""

# ╔═╡ b35f5984-9f46-4303-9f0a-af7c79cf7a53
md"""
```julia
julia> Genie.WebChannels.broadcast("__", "Hey!")
true
```
"""

# ╔═╡ 4b6f463c-7453-4f06-a51d-29cc42e66550
md"""
If you look in the browser's console you will see the "Hey!" message! By default, the client side handler simply outputs the message. We're also informed that we can "Overwrite window.parse_payload to handle messages from the server". Let's do it. Run this in the current REPL (it will overwrite our root route handler):
"""

# ╔═╡ 3085ac1c-ebf5-4199-b255-a7c10efafbd2
route("/") do
  Assets.channels_support() *
  """
  <script>
  window.parse_payload = function(payload) {
    console.log('Got this payload: ' + payload);
  }
  </script>
  """
end;

# ╔═╡ 0135778f-63a0-48e1-8e12-e3407d8d5c8b
md"""
Now if you reload the page and broadcast the message, it will be picked up by our custom payload handler. However, chances are you'll also get an error when broadcasting (don't worry though, the error is just logged, it's not breaking the application as it's not critical):
"""

# ╔═╡ e65f8790-8bbe-49ce-9678-23dfa714b6c4
md"""

```julia
julia> Genie.WebChannels.broadcast("__", "Hey!")
┌ Error: Base.IOError("stream is closed or unusable", 0)
└ @ Genie.WebChannels ~/.julia/dev/Genie/src/WebChannels.jl:220
true

```
"""

# ╔═╡ 3c9eb96d-a4e3-48a5-8d90-31c0dea9bd1f
md"""
The error is caused by the fact that, by reloading the page, our previously connected WebSockets client is now unreachable. However, we still keep a reference to it - and when we try to broadcast to it, we find that the stream has been closed. We can fix this by calling
"""

# ╔═╡ 1e16ffe0-f14a-42d9-b05b-ff899763679a
md"""
```julia
julia> Genie.WebChannels.unsubscribe_disconnected_clients()
```
"""

# ╔═╡ 6fa1c839-4dd4-421b-992f-13827883cea6
md"""
The output of `unsubscribe_disconnected_clients()` is the collection of remaining (connected) clients.
"""

# ╔═╡ 3e6e402d-601c-45ea-acad-44cb9ed5219d
md"""
---

**Heads up!**

Although harmless, the error indicates that you have disconnected clients in memory. If you don't need the data, purge them to free memory.

---

At any time, we can check the connected clients with `Genie.WebChannels.connected_clients()` and the disconnected ones with `Genie.WebChannels.disconnected_clients()`.

### Pushing messages from the client

We can also push messages from client to server. As we don't have a UI, we'll use the browser's console and Genie's JavaScript API to send the message. But first, we need to set up the `channel` which will receive our message. Run this in the active Julia REPL:
"""

# ╔═╡ 989907fb-8179-4029-92a4-e0ee44e52ccb
channel("/__/echo") do
  "Received: $(params(:payload))"
end;

# ╔═╡ bf14645c-3fbe-4531-983a-998cc5d629ea
md"""
Now that our endpoint is up, go to the browser's console and run:
```javascript
Genie.WebChannels.sendMessageTo('__', 'echo', 'Hello!')
```
"""

# ╔═╡ f57df6f7-e612-47d8-8856-c2b0eddde564
md"""
The console will immediately display the response from the server:

```text
Received: Hello!  channels.js:74:3
Got this payload: Received: Hello!
```
"""

# ╔═╡ 4e8d5270-116d-4f1c-b410-95724c97d359
md"""
## Summary

This concludes our intro to working with WebSockets in Genie. You now have the knowledge to set up the communication between client and server, send messages from both server and clients, and perform various tasks using the `WebChannels` API.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~3.3.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "b2667530e42347b10c10ba6623cfebc09ac5c7b6"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.2.4"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "d962b5a47b6d191dbcd8ae0db841bc70a05a3f5b"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.13"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "YAML"]
git-tree-sha1 = "4f5526913b239ece648a73d4445c7978f1251935"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "3.3.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "10c95cebcfa37c1f510a726c90886db4745e1238"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.15.11"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "d911b6a12ba974dabe2291c6d450094a7226b372"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.1"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─32df4cc8-39ac-11ec-2f5e-b176daa6db9a
# ╟─a54bd0a5-7ec6-46db-8d42-a4cfc3eb1c2f
# ╟─c7ea1c3e-2a61-41d4-ba37-21027f56453a
# ╠═332e4ec2-b1c5-47c0-a12e-db0195023e18
# ╠═a8b81c5a-6a04-4ef6-8522-82582da4f78a
# ╟─72010195-6acc-4a4c-bb58-5331f8007192
# ╟─dc4430bf-4555-4b3a-9c9c-9c14ef4ac806
# ╟─bd80dd1c-6e99-4c60-91d3-e470b9f5020c
# ╠═a6ec9f7f-c574-4af0-a755-5d090a00d423
# ╟─b65c4644-f5f8-4a68-8df7-8064b1274ba5
# ╟─f082d891-61ba-4a9f-b1ea-943ed2aba815
# ╟─3911bc31-f4ea-493e-95a9-9d743670c4c6
# ╠═69bec8f4-56b1-48ee-aa2c-ba4310e79568
# ╟─266ce7f3-4364-44ec-90b8-61e06dba5d27
# ╠═70db282b-ec1c-40f7-bd7c-59ed2e94891c
# ╟─45c63333-2726-400e-a179-996ce88f456a
# ╟─15c148af-9708-4cec-a390-17a87a02a8fb
# ╠═78550cf1-a2ef-487e-87a1-589fe3c7c048
# ╟─5350a2e0-2f95-4fb6-a9e5-18b4fc58d84f
# ╠═52576077-67bd-4a24-b81b-be6f475d4192
# ╟─e8137d45-a1c8-480c-80b3-05454f3d241b
# ╟─a26be9e1-42ce-43d6-8c36-d2a08ba75f99
# ╟─0f99c27a-29df-4aad-a387-12799664569f
# ╟─fabc9669-ce12-470f-b551-3d8f69cb221b
# ╠═79660882-e6d9-4125-8601-277888c199a5
# ╟─229008fe-bad8-455b-853f-57eb39e87665
# ╟─b35f5984-9f46-4303-9f0a-af7c79cf7a53
# ╟─4b6f463c-7453-4f06-a51d-29cc42e66550
# ╠═3085ac1c-ebf5-4199-b255-a7c10efafbd2
# ╟─0135778f-63a0-48e1-8e12-e3407d8d5c8b
# ╟─e65f8790-8bbe-49ce-9678-23dfa714b6c4
# ╟─3c9eb96d-a4e3-48a5-8d90-31c0dea9bd1f
# ╟─1e16ffe0-f14a-42d9-b05b-ff899763679a
# ╟─6fa1c839-4dd4-421b-992f-13827883cea6
# ╟─3e6e402d-601c-45ea-acad-44cb9ed5219d
# ╠═989907fb-8179-4029-92a4-e0ee44e52ccb
# ╟─bf14645c-3fbe-4531-983a-998cc5d629ea
# ╟─f57df6f7-e612-47d8-8856-c2b0eddde564
# ╟─4e8d5270-116d-4f1c-b410-95724c97d359
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Parameter searching api

"""
Parameter searching Api. It can be run in single core or multiple-processing

Performance
------------
```
julia> addprocs(exeflags="--project=$(Base.active_project())");
julia> @time @everywhere using RetroSignalModel
julia> nprocs()
9
julia> @time rs.paramSearching(rs.rtgM4(); num_sim=100, distributed=true, saveall=true)
  2.378919 seconds (78.95 k allocations: 4.964 MiB)
```
"""
function paramSearching(m::RTGmodel;num_sim=10, save_iter=num_sim, distributed=false, saveall=true,batch_size=nprocs(),conditions=getConditions(;df=DataTables.BoolCond), kwags...)
    #create table
    df = get_search_df(m, conditions)

    ran = 1:save_iter:num_sim
    iter = ran.step
    try
        for _ in 1:length(ran)
            # parameter searching
            res =  @showprogress  pmap( x->paramSearching_Once(x[1];conditions=x[2], x[3]...),
                                    fill([m, conditions, kwags], iter),
                                    distributed=distributed,
                                    batch_size=batch_size)

            # Merge result
            map(r-> r.isValid || saveall ? push_sim!(df, r.test_log, r.p, r.u) : nothing, res)
        end
        return df
    catch e
        return df
    end

end

"""
Sampling one parameter set for validation
"""
function paramSearching_Once(m::RTGmodel;Break=false, idx_hill_coefs=[1], K_dist=K_dist, K_N_dist=K_N_dist,conditions=getConditions(;df=DataTables.BoolCond),show_valid=true,ssmethod=SSMETHOD,kwags...)
    # Sampling parameters
    p = init_p(m.model;idx_hill_coefs=idx_hill_coefs, K_dist=K_dist, K_N_dist=K_N_dist)
    m.p[:] = p

    u = getSteadySol(m, m.u; ssmethod=ssmethod).u

    # Test validity
    isValid, test_log = try_conditions(m;conditions=conditions, Break=Break,kwags...)

    # Print valid result
    show_valid && isValid ? @show(p) : nothing

    return (; isValid=isValid, test_log=test_log, u=u, p=p)
end

# Data processing
function get_search_df(m, conditions)
    # Labels
    condNs = ["condition $i" for i in range(1,stop= size(conditions)[1])]
    pNs = string.(keys(init_p(m.model)))
    uNames = string.(keys(init_u(m)))
    # Create empty dataframe
    df_key = [condNs..., pNs...,uNames...]
    trials = DataFrame( [ [] for i in eachindex(df_key) ], df_key)

    return trials
end

function push_sim!(df, conds, ps, us)
    vec = [conds..., ps...,us...]
    push!(df, vec)
end


text:
@add_problem mip function mip_lp_fallback_interface(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable()
    p = minimize(x, x>=4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.3 atol=atol rtol=rtol
    end

    x = Variable(2)
    p = minimize(norm(x,1), x[1]>=4.3; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4.3 atol=atol rtol=rtol
    end
end

@add_problem mip function mip_integer_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(:Int)
    p = minimize(x, x>=4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol=atol rtol=rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x>=4.3; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 10 atol=atol rtol=rtol
    end

    x = Variable(:Int)
    y = Variable()
    p = minimize(sum(x + y), x>=4.3, y>=7; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 12 atol=atol rtol=rtol
    end

    x = Variable(2, :Int)
    p = minimize(norm(x, 1), x[1]>=4.3; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol=atol rtol=rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x[1]>=4.3, x>=0; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 5 atol=atol rtol=rtol
    end

    x = Variable(2, :Int)
    p = minimize(sum(x), x>=.5; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol=atol rtol=rtol
    end
end

@add_problem mip function mip_binary_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2, :Bin)
    p = minimize(sum(x), x>=.5; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol=atol rtol=rtol
    end

    x = Variable(2, :Bin)
    p = minimize(sum(x), x[1]>=.5, x>=0; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol=atol rtol=rtol
    end
end


text:
module BandedMatrices
using Base, FillArrays, LazyArrays, MatrixFactorizations, LinearAlgebra, SparseArrays, Random
using LinearAlgebra.LAPACK
import Base: axes, axes1, getproperty, iterate, tail
import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!,
                        copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym,
                        _chol!
import LinearAlgebra.BLAS: libblas
import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
import LinearAlgebra: cholesky, cholesky!, cholcopy, norm, diag, eigvals!, eigvals, eigen!, eigen,
            qr, qr!, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular,
            chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet,
            svdvals, svdvals!, QRPackedQ, checknonsingular, ipiv2perm, _apply_ipiv!, tril!,
            triu!, Givens
import MatrixFactorizations: ql, ql!, QLPackedQ, reflector!, reflectorApply!
import SparseArrays: sparse

import Base: getindex, setindex!, *, +, -, ==, <, <=, >,
                >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds

import Base: convert, size, view, unsafe_indices,
                first, last, size, length, unsafe_length, step,
                to_indices, to_index, show, fill!, promote_op,
                MultiplicativeInverses, OneTo, ReshapedArray,
                               similar, copy, convert, promote_rule, rand,
                            IndexStyle, real, imag, Slice, pointer, unsafe_convert, copyto!

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                        materialize, materialize!

import LazyArrays: MemoryLayout, @lazymul, @lazylmul, @lazyldiv,
                    AbstractStridedLayout, AbstractColumnMajor, AbstractRowMajor,
                    _copyto!, MatMulVec, MatMulMat, transposelayout, triangulardata,
                    ConjLayout, conjlayout, SymmetricLayout, symmetriclayout, symmetricdata,
                    triangularlayout, MatMulVec, MatLdivVec, TriangularLayout,
                    AbstractBandedLayout, DiagonalLayout, LayoutApplyStyle,
                    ArrayMulArrayStyle, HermitianLayout, hermitianlayout, hermitiandata,
                    MulAdd, materialize!, BlasMatMulMat, BlasMatMulVec, VcatLayout, ZerosLayout,
                    AbstractColumnMajor, MulLayout, colsupport, rowsupport,
                    DenseColumnMajor, DenseRowMajor, ArrayLdivArray
import FillArrays: AbstractFill

export BandedMatrix,
       bandrange,
       brand,
       bandwidth,
       BandError,
       band,
       Band,
       BandRange,
       bandwidths,
       colrange,
       rowrange,
       isbanded,
       Zeros,
       Fill,
       Ones,
       Eye

if VERSION < v"1.2-"
    import Base: has_offset_axes
    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    import Base: require_one_based_indexing
end

include("blas.jl")
include("lapack.jl")

include("generic/AbstractBandedMatrix.jl")
include("generic/broadcast.jl")
include("generic/matmul.jl")
include("generic/Band.jl")
include("generic/utils.jl")
include("generic/indexing.jl")

include("banded/BandedMatrix.jl")
include("banded/BandedLU.jl")
include("banded/bandedqr.jl")
include("banded/bandedql.jl")
include("banded/gbmm.jl")
include("banded/linalg.jl")

include("symbanded/symbanded.jl")
include("symbanded/ldlt.jl")
include("symbanded/BandedCholesky.jl")
include("symbanded/SplitCholesky.jl")
include("symbanded/bandedeigen.jl")

include("tribanded.jl")

include("interfaceimpl.jl")

@deprecate setindex!(A::BandedMatrix, v, b::Band) A[b] .= v
@deprecate setindex!(A::BandedMatrix, v, ::BandRangeType, j::Integer) A[BandRange,j] .= v
@deprecate setindex!(A::BandedMatrix, v, kr::Colon, j::Integer) A[:,j] .= v
@deprecate setindex!(A::BandedMatrix, v, kr::AbstractRange, j::Integer) A[kr,j] .= v
@deprecate setindex!(A::BandedMatrix, v, ::Colon, ::Colon) A[:,:] .= v
@deprecate setindex!(A::BandedMatrix, v, ::Colon) A[:] .= v
@deprecate setindex!(A::BandedMatrix, v, k::Integer, ::BandRangeType) A[k,BandRange] .= v
@deprecate setindex!(A::BandedMatrix, v, k::Integer, jr::AbstractRange) A[k,jr] .= v
@deprecate setindex!(A::BandedMatrix, v, kr::AbstractRange, jr::AbstractRange) A[kr,jr] .= v
@deprecate setindex!(A::BandedMatrix, v, ::BandRangeType) BandedMatrices.bandeddata(A) .= v

end #module


text:
using Revise
using PyCall
using PyPlot
using Bem2d

function plotfunction(els, xcenters, ycenters, xnodes, ynodes, uconst, σconst, tconst, uquad, σquad, tquad, xobs, yobs, uconstfarfield, σconstfarfield, uquadfarfield, σquadfarfield)
    scalestress = 1e6
    contourlevels = 5
    contourlevels = collect(LinRange(4, 7, 6))
    fontsize = 7
    xmin = -5e3
    xmax = 5e3
    cmap = rycroftcmap()

    figure(figsize = (6, 5), dpi=300)
    ax = subplot(3, 3, 1)
    mat = reshape(log10.(abs.(σconstfarfield[:, 1])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{xx}\|$ (Pa) CS", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels, linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    ylabel(L"$y$ (m)", fontsize = fontsize); xlim([xmin, xmax]); xticks([])
    yticks([-2500, 0, 2500]);
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 2)
    mat = reshape(log10.(abs.(σconstfarfield[:, 2])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{yy}\|$ (Pa) CS", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels, linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    xlim([xmin, xmax]); xticks([]); yticks([])
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 3)
    mat = reshape(log10.(abs.(σconstfarfield[:, 3])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{xy}\|$ (Pa) CS", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels, linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    xlim([xmin, xmax]); xticks([]); yticks([])
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 4)
    mat = reshape(log10.(abs.(σquadfarfield[:, 1])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{xx}\|$ (Pa) 3QN", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels, linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    ylabel(L"$y$ (m)", fontsize = fontsize); xlim([xmin, xmax]); xticks([])
    yticks([-2500, 0, 2500]);
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 5)
    mat = reshape(log10.(abs.(σquadfarfield[:, 2])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{yy}\|$ (Pa) 3QN", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels, linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    xlim([xmin, xmax]); xticks([]); yticks([])
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 6)
    mat = reshape(log10.(abs.(σquadfarfield[:, 3])), size(xobs))
    mat[findall(x -> x < minimum(contourlevels), mat)] .= minimum(contourlevels)
    mat[findall(x -> x > maximum(contourlevels), mat)] .= maximum(contourlevels)
    contourf(xobs, yobs, mat, contourlevels, cmap = cmap)
    cbar = colorbar(orientation = "horizontal", fraction = 0.20, pad = 0.05, extend = "both")
    cbar.ax.tick_params(labelsize = fontsize)
    cbar.set_label(label = L"$\log\|\sigma_{xy}\|$ (Pa) 3QN", fontsize = fontsize)
    contour(xobs, yobs, mat, contourlevels,  linewidths = 0.25, colors = "k", linestyles = "-")
    gca().set_aspect("equal")
    xlim([xmin, xmax]); xticks([]); yticks([])
    ax.tick_params(labelsize = fontsize)

    ax = subplot(3, 3, 7)
    plot(xcenters, σconst[1:3:end] ./ scalestress, ".k", markerfacecolor = "blue", markeredgewidth = 0.0, label = L"$σ_{xx}$ CS", zorder = 2)
    plot(xnodes, σquad[1:3:end] ./ scalestress, "ok", markerfacecolor = "red", markeredgewidth = 0.25, label = L"$σ_{xx}$ 3QN", zorder = 1)
    legend(loc = "lower right", fontsize = fontsize, frameon=false)
    xlim([xmin, xmax]); xticks([xmin, 0, xmax])
    xlabel(L"$x$ (m)", fontsize = fontsize)
    ylabel(L"$\sigma$ (MPa)", fontsize = fontsize)
    ax.tick_params(labelsize = fontsize)
    yticks([-10, 0, 10])
    ylim([-1e7 / scalestress, 1e7 / scalestress])

    ax = subplot(3, 3, 8)
    plot(xcenters, σconst[2:3:end] ./ scalestress, ".k", markerfacecolor = "blue", markeredgewidth = 0.0, label = L"$σ_{yy}$ CS", zorder = 2)
    plot(xnodes, σquad[2:3:end] ./ scalestress, "ok", markerfacecolor = "red", markeredgewidth = 0.25, label = L"$σ_{yy}$ 3QN", zorder = 1)
    legend(loc = "lower right", fontsize = fontsize, frameon=false)
    xlim([xmin, xmax]); xticks([xmin, 0, xmax])
    xlabel(L"$x$ (m)", fontsize = fontsize)
    ax.tick_params(labelsize = fontsize)
    ylim([-1e7 / scalestress, 1e7 / scalestress]); yticks([])

    ax = subplot(3, 3, 9)
    plot(xcenters, σconst[3:3:end] ./ scalestress, ".k", markerfacecolor = "blue", markeredgewidth = 0.0, label = L"$σ_{xy}$ CS", zorder = 2)
    plot(xnodes, σquad[3:3:end] ./ scalestress, "ok", markerfacecolor = "red", markeredgewidth = 0.25, label = L"$σ_{xy}$ 3QN", zorder = 1)
    legend(loc = "lower right", fontsize = fontsize, frameon=false)
    xlim([xmin, xmax]); xticks([xmin, 0, xmax])
    xlabel(L"$x$ (m)", fontsize = fontsize)    # ylabel(L"$\sigma$ (MPa)", fontsize = fontsize)
    ax.tick_params(labelsize = fontsize)
    ylim([-1e7 / scalestress, 1e7 / scalestress]); yticks([])
    tight_layout()
    show()
end

function fig_constquadstresscomp()
    # Material and geometric constants
    mu = 3e10
    nu = 0.25
    nels = 20
    els = Elements(Int(1e5))
    L = 10000
    x1, y1, x2, y2 = discretizedline(-L, 0, L, 0, nels)
    x2[11] = 500.0
    x1[12] = 500.0
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Convenience data structures
    idx = getidxdict(els)
    partialsconst = initpartials(els)
    partialsquad = initpartials(els)

    partialsconst["disp"]["fault"]["fault"], partialsconst["stress"]["fault"]["fault"], partialsconst["trac"]["fault"]["fault"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["fault"], mu, nu)
    partialsquad["disp"]["fault"]["fault"], partialsquad["stress"]["fault"]["fault"], partialsquad["trac"]["fault"]["fault"] = partialsquaddispstress(slip2dispstress, els, idx["fault"], idx["fault"], mu, nu)

    # Observation coordinates for plotting
    npts = 300
    obswidth = 5e3
    xobs, yobs = obsgrid(-obswidth, -2.5e3, obswidth, 2.5e3, npts)

    # Evaluation points
    xcenters = els.xcenter[1:els.endidx]
    ycenters = els.ycenter[1:els.endidx]
    xnodes = sort(els.xnodes[1:els.endidx, :][:])
    ynodes = sort(els.ynodes[1:els.endidx, :][:])

    # Slip for volumetric models
    slope = 0.0001
    constxslip = zeros(2, nels)
    constyslip = zeros(2, nels)
    constxslip[1, :] .= 1.0
    constxslip[2, :] = slope .* xcenters
    quadxslip = zeros(2, nels, 3)
    quadyslip = zeros(2, nels, 3)
    quadxslip[1, :, :] .= 1
    quadxslip[2, :, :] = transpose(reshape(slope .* xnodes, 3, nels))

    # Slip for partials
    slipconst = zeros(2, 2 * nels)
    slipquad = zeros(2, 6 * nels)

    # Constant strike-slip
    slipconst[1, 1:2:end] .= 1
    slipquad[1, 1:2:end] .= 1

    # Linear strike-slip
    slope = 0.0001
    slipconst[2, 1:2:end] = slope * xcenters
    slipquad[2, 1:2:end] = slope * xnodes

    # Predict on-fault displacements, stresses, and tractions
    close("all")
    for i in 2:size(slipconst)[1]
        # On fault evaluation
        dispconst = partialsconst["disp"]["fault"]["fault"] * slipconst[i, :]
        stressconst = partialsconst["stress"]["fault"]["fault"] * slipconst[i, :]
        tracconst = partialsconst["trac"]["fault"]["fault"] * slipconst[i, :]
        dispquad = partialsquad["disp"]["fault"]["fault"] * slipquad[i, :]
        stressquad = partialsquad["stress"]["fault"]["fault"] * slipquad[i, :]
        tracquad = partialsquad["trac"]["fault"]["fault"] * slipquad[i, :]

        # Far-field evaluation
        dispconstfarfield, stressconstfarfield = constdispstress(slip2dispstress, xobs, yobs, els, idx["fault"], constxslip[i, :], constyslip[i, :], mu, nu)
        dispquadfarfield, stressquadfarfield = quaddispstress(slip2dispstress, xobs, yobs, els, idx["fault"], quadxslip[i, :, :], quadyslip[i, :, :], mu, nu)

        # Plot 9 panel results
        plotfunction(els, xcenters, ycenters, xnodes, ynodes, dispconst, stressconst, tracconst, dispquad, stressquad, tracquad, reshape(xobs, npts, npts), reshape(yobs, npts, npts), dispconstfarfield, stressconstfarfield, dispquadfarfield, stressquadfarfield)
    end

end
fig_constquadstresscomp()


text:
#
#----------------------------------------------------------------------
export ConvectionDiffusion
#----------------------------------------------------------------------
struct ConvectionDiffusion{T,U} <: Equation
    name::String

    fld::Field{T}
    vx ::Array{T}
    vy ::Array{T}

    ν  ::Array{T}     # viscosity
    f  ::Array{T}     # forcing
    rhs::Array{T}     # RHS
    exH::Array{Array} # storage history of explicit term (convection)

    tstep::TimeStepper{T,U}

    JrVD::Array{T,2}
    JsVD::Array{T,2}

    mshV::Mesh{T} # velocity/scalar mesh
    mshD::Mesh{T} # dealiasing mesh

    set0!::Function # initial condition
    set∂!::Function # dirichlet boundary condition
    setF!::Function # forcing function
    setν!::Function # viscosity
end
#--------------------------------------#
function ConvectionDiffusion(name::String
                            ,fld::Field,vx::Array,vy::Array
                            ,tstep::TimeStepper
                            ,mshD::Mesh
                            ,set0! =fixU!,set∂! =fixU!
                            ,setF! =fixU!,setν! =fixU!)

    ν   = zero(fld.u)
    f   = zero(fld.u)
    rhs = zero(fld.u)

    k = length(fld.uh)
    exH = Array[zero(fld.u) for i in 1:k]

    mshV = fld.msh
    JrVD = interpMat(mshD.zr,mshV.zr)
    JsVD = interpMat(mshD.zs,mshV.zs)

    return ConvectionDiffusion(name
                              ,fld,vx,vy
                              ,ν,f,rhs,exH
                              ,tstep
                              ,JrVD,JsVD
                              ,mshV,mshD
                              ,set0!,set∂!,setF!,setν!)
end
#--------------------------------------#
#function ConvectionDiffusion(name::String,bc::Array{Char,1}
#                            ,mshV::Mesh,mshD::Mesh
#                            ;Ti=0.,Tf=0.,dt=0.,k=3)
#
#    fld = Field(bc,mshV)
#    vx  = zero(fld.u)
#    vy  = zero(fld.u)
#
#    tstep = TimeStepper(Ti,Tf,dt,k)
#
#    return ConvectionDiffusion(name
#                              ,fld,vx,vy
#                              ,tstep
#                              ,mshD)
#end
#----------------------------------------------------------------------
# solve
#----------------------------------------------------------------------
function opLHS(u::Array,cdn::ConvectionDiffusion)
    @unpack fld, mshV, ν = cdn
    @unpack bdfB = cdn.tstep

    Au = hlmz(u,ν,bdfB[1],mshV)

    Au .= gatherScatter(Au,mshV)
    Au .= mask(Au,fld.M)
    return Au
end
#--------------------------------------#
function opPrecond(u::Array,cdn::ConvectionDiffusion)
    @unpack fld, tstep, mshV = cdn
    Mu = u ./ mshV.B ./ tstep.bdfB[1]
    return Mu
end
#--------------------------------------#
function makeRHS!(cdn::ConvectionDiffusion)
    @unpack fld, rhs, ν, f, mshV, mshD = cdn
    @unpack vx,vy,exH,JrVD,JsVD = cdn
    @unpack bdfA, bdfB = cdn.tstep

    rhs  .=      mass(f     ,mshV) # forcing
    rhs .-= ν .* lapl(fld.ub,mshV) # boundary data

    for i=1:length(fld.uh)
        exH[i] .= -advect(fld.uh[i],vx,vy,mshV,mshD,JrVD,JsVD)
        rhs   .-= bdfB[1+i] .* mass(fld.uh[i],mshV)
        rhs   .+= bdfA[i]   .* exH[i]
    end

    rhs .= mask(rhs,fld.M)
    rhs .= gatherScatter(rhs,mshV)
    return
end
#--------------------------------------#
function solve!(cdn::ConvectionDiffusion)
    @unpack rhs, mshV, fld = cdn
    @unpack u,ub = fld

    opL(u) = opLHS(u,cdn)
    opM(u) = opPrecond(u,cdn)

    pcg!(u,rhs,opL;opM=opM,mult=mshV.mult,ifv=false)
    u .+= ub
    return
end
#----------------------------------------------------------------------
# time evolution
#----------------------------------------------------------------------
function updateHist!(cdn::ConvectionDiffusion)

    updateHist!(cdn.fld)

    return
end
#--------------------------------------#
function evolve!(cdn::ConvectionDiffusion)

    @unpack fld, f, ν, mshV, tstep = cdn
    @unpack time = tstep

    cdn.set∂!(fld.ub,mshV.x,mshV.y,time[1])
    cdn.setF!(f     ,mshV.x,mshV.y,time[1])
    cdn.setν!(ν     ,mshV.x,mshV.y,time[1])

    makeRHS!(cdn)
    solve!(cdn)

    return
end
#--------------------------------------#
export step!
#--------------------------------------#
function step!(cdn::ConvectionDiffusion)

    updateHist!(cdn)
    updateHist!(cdn.tstep)
    evolve!(cdn)

    return
end
#--------------------------------------#
export simulate!
#--------------------------------------#
function simulate!(cdn::ConvectionDiffusion,callback!::Function)

    @unpack fld, mshV, tstep = cdn
    @unpack time, istep, dt, Tf = tstep

    cdn.set0!(fld.u,mshV.x,mshV.y,time[1])

    callback!(cdn)
    while time[1] <= Tf[1]

        step!(cdn)
        callback!(cdn)

        if(time[1] < 1e-12) break end

    end

    return
end
#----------------------------------------------------------------------
#


text:
module Test_PathPlanning

import Modia3D
using  ModiaLang
@usingModiaPlot

const ptp_path = Modia3D.PTP_path(["angle1", "angle2", "angle3"],
                                  positions=[0.0 2.0 3.0;
                                             0.5 3.0 4.0;
                                             0.8 1.5 0.3;
                                             0.2 1.5 0.8],
                                  startTime=0.1,
                                  v_max=2*ones(3),
                                  a_max=3*ones(3))

Modia3D.plotPath(ptp_path, plot)

Modia3D.plotPath(ptp_path, plot, names=("angle2", "angle3"), heading="Test of PTP plots",
                 tend=0.9*ptp_path.Tend, figure=2, ntime=99)

Modia3D.plotPath(ptp_path, plot, onlyPositions=false, ntime=1000, figure=3)

Modia3D.plotPath(ptp_path, plot, names=("angle2", "angle3"), heading="Test of PTP plots",
                 tend = 0.9*ptp_path.Tend, onlyPositions=false, ntime=1000, figure=4)

# Test getPosition
angles = zeros(3)
Modia3D.getPosition!(ptp_path, 0.5, angles)
println("... angles(time=0.5) = $angles")

angle2 = Modia3D.getPosition(ptp_path, 2, 0.5)
println("... angle2(time=0.5) = ", angle2)

der_angles  = zeros(3)
der2_angles = zeros(3)
Modia3D.getPosition!(ptp_path, 0.5, angles, der_angles, der2_angles)
println("... angles(          time=0.5) = $angles")
println("... der(angles)(     time=0.5) = $der_angles")
println("... der(der(angles))(time=0.5) = $der2_angles")

end


text:
module SamplingHS3

export mh_step , hr_step
export HRSamplingTrackingData
export NSConfig , nested_sampling

export interp1_lin , bin2d

export LineSegment , create_line , dist2LS , dist2line
export GaussMixture , GaussParam , draw_gaussmixture , mvnpdf

# package code goes here
include("util/interp.jl")
include("util/basics.jl")
include("util/stopcontrol.jl")

include("samplers/hr.jl")
include("samplers/mh.jl")

include("util/dist2line.jl")

include("ns/ns.jl")
include("samplers/mcmc.jl")

include("ns/ns_utils.jl")
include("ns/ns_analysis.jl")

#include("util/plot_utils.jl")

end # module


text:
# This file is a part of RadiationSpectra.jl, licensed under the MIT License (MIT).

__precompile__(true)

module RadiationSpectra

using DelimitedFiles
using LinearAlgebra
using Serialization
using Statistics

using ArraysOfArrays
using Distributions
using ForwardDiff
using IntervalSets
using LinearAlgebra
using LsqFit
using NLSolversBase
using Optim
using RecipesBase
using Requires
using SpecialFunctions
using StatsBase
using Tables
using TypedTables
using ValueShapes
using Random123

import Base: print, println, show
import Statistics: mean
import StatsBase: std, mean_and_std

export FitFunction
export get_ndims, get_nparams
export set_fitranges!
export set_parameter_names!
export set_initial_parameters!
export set_parameter_bounds!
export get_fitted_parameters
export get_initial_parameters
export get_standard_deviations

export llhfit!, lsqfit!, linear_regression!

export peakfinder
export calibrate_spectrum

export mean, std, mean_and_std

export fit_single_peak_histogram
export fit_single_peak_histogram_refined

"""
    get_example_spectrum()::Histogram

Returns an uncalibrated radiation spectrum for testing and demonstrating purpose.
"""
function get_example_spectrum()::Histogram
    weights::Array{Int, 1} = readdlm(joinpath(@__DIR__, "..", "examples", "data", "hpge-spectrum.txt"))[:]
    h = Histogram(0:100:100 * length(weights), :left)
    h.weights = weights
    return h
end

include("Statistics/Statistics.jl")

include("PeakFinder/PeakFinder.jl")

include("Fitting/FitFunction.jl")
include("Fitting/general_model_functions.jl")

include("Fitting/llhfit.jl")

include("Fitting/lsqfit.jl")

include("Fitting/predefined_fit_routines.jl")

include("AutoCalibration/AutoCalibration.jl")

function __init__()
    @require BAT="c0cd4b16-88b7-57fa-983b-ab80aecada7e" include("Fitting/BAT/BAT.jl")
end

end # module


text:
export DrumBoiler

"""
    DrumBoiler

Type representing the IFAC drum boiler benchmark problem.

From E. J. Davison, ‘Benchmark Problems for Control System Design: Report of the IFAC Theory Committee’, IFAC, 1990.
"""
struct DrumBoiler <: ControlBenchmarkProblem end

@addbenchmark DrumBoiler "ifac02" "Drum boiler" "A drum boiler where water is converted to steam" Set([:linear, :unstable, :continuous, :process, :nonminimum_phase, :ifac])

"""
    controlbenchmark(::DrumBoiler) -> sys

Create the system for the IFAC drum boiler sample problem.

`sys` is the linear statespace system describing the drum boiler.
"""
function controlbenchmark( ::DrumBoiler )

    A = [ -3.93     -3.15e-3     0.00        0.00        0.00        4.03e-5     0.00        0.00        0.00;
           3.68e2   -3.05        3.03        0.00        0.00       -3.77e-3     0.00        0.00        0.00;
           2.74e1    7.87e-2    -5.96e-2     0.00        0.00       -2.81e-4     0.00        0.00        0.00;
          -6.47e-2  -5.2e-5      0.00       -2.55e-1    -3.35e-6     3.60e-7     6.33e-5     1.94e-4     0.00;
           3.85e3    1.73e1     -1.28e1     -1.26e4     -2.91       -1.05e-1     1.27e1      4.31e1      0.00;
           2.24e4    1.80e1      0.00       -3.56e1     -1.04e-4    -4.14e-1     9.00e1      5.69e1      0.00;
           0.00      0.00        2.34e-3     0.00        0.00        2.22e-4    -2.03e-1     0.00        0.00;
           0.00      0.00        0.00       -1.27        1.00e-3     7.86e-5     0.00       -7.17e-2     0.00;
          -2.20     -1.77e-3     0.00       -8.44       -1.11e-4     1.38e-5     1.49e-3     6.02e-3    -1.00e-10 ]

    B = [ 0.00   0.00       0.00;
          0.00   0.00       0.00;
          1.56   0.00       0.00;
          0.00  -5.13e-6    0.00;
          8.28  -1.50       3.95e-2;
          0.00   1.78       0.00;
          2.33   0.00       0.00;
          0.00  -2.45e-2    2.84e-3;
          0.00   2.94e-5    0.00 ]

    C = zeros( Float64, 2, 9 )
    C[1, 6] = 1
    C[2, 9] = 1

    D = zeros( Float64, 2, 3 )

    E = [ -1.00e-2;
           0.00;
           0.00;
           0.00;
           5.20e1;
           0.00;
           0.00;
           0.00;
           0.00 ]

    return (; :sys => StateSpace( A, B, C, D, Continuous() ) )
end


text:
#=
   SEND+MORE=MONEY in Julia + ConstraintSolver.jl

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

=#

using ConstraintSolver, JuMP
const CS = ConstraintSolver
include("constraints_utils.jl")

function send_more_money()
    model = Model(optimizer_with_attributes(CS.Optimizer, "all_solutions"=>true,"logging"=>[]))
    @variable(model, 0 <= x[1:8] <= 9, Int)
    s,e,n,d,m,o,r,y = x

    @constraint(model, x in CS.AllDifferent())
    @constraint(model, s != 0)
    @constraint(model, m != 0)
    @constraint(model,
                              1000*s + 100*e + 10*n + d
                           +  1000*m + 100*o + 10*r + e
                 == 10000*m + 1000*o + 100*n + 10*e + y
        )

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        x = convert.(Integer,JuMP.value.(x))
        println(x)
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
    end
end

send_more_money()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "Conversions" begin

    @testset "Conversions: Weight Conversions" begin
        @test weight_conversion("kilogram", "kilogram", 4) == 4
        @test weight_conversion("kilogram", "gram", 1) == 1000
        @test weight_conversion("kilogram", "milligram", 4) == 4000000
        @test weight_conversion("kilogram", "metric-ton", 4) ≈ 0.004
        @test weight_conversion("kilogram", "long-ton", 3) == 0.0029526219
        @test weight_conversion("kilogram", "short-ton", 1) == 0.0011023122
        @test weight_conversion("kilogram", "pound", 4) == 8.8184976808
        @test weight_conversion("kilogram", "ounce", 4) == 141.095962892
        @test weight_conversion("kilogram", "carrat", 3) == 15000
        @test weight_conversion("kilogram", "atomic-mass-unit", 1) == 6.022136652e+26
        @test weight_conversion("gram", "kilogram", 1) ≈ 0.001
        @test weight_conversion("gram", "gram", 3) ≈ 3.0
        @test weight_conversion("gram", "milligram", 2) ≈ 2000.0
        @test weight_conversion("gram", "metric-ton", 4) ≈ 4e-06
        @test weight_conversion("gram", "long-ton", 3) ≈ 2.9526219e-06
        @test weight_conversion("gram", "short-ton", 3) ≈ 3.3069366000000003e-06
        @test weight_conversion("gram", "pound", 3) ≈ 0.0066138732606
        @test weight_conversion("gram", "ounce", 1) ≈ 0.035273990723
        @test weight_conversion("gram", "carrat", 2) ≈ 10.0
        @test weight_conversion("gram", "atomic-mass-unit", 1) ≈ 6.022136652e+23
        @test weight_conversion("milligram", "kilogram", 1) ≈ 1e-06
        @test weight_conversion("milligram", "gram", 2) ≈ 0.002
        @test weight_conversion("milligram", "milligram", 3) ≈ 3.0
        @test weight_conversion("milligram", "metric-ton", 3) ≈ 3e-09
        @test weight_conversion("milligram", "long-ton", 3) ≈ 2.9526219e-09
        @test weight_conversion("milligram", "short-ton", 1) ≈ 1.1023122e-09
        @test weight_conversion("milligram", "pound", 3) ≈ 6.6138732605999995e-06
        @test weight_conversion("milligram", "ounce", 2) ≈ 7.054798144599999e-05
        @test weight_conversion("milligram", "carrat", 1) ≈ 0.005
        @test weight_conversion("milligram", "atomic-mass-unit", 1) ≈ 6.022136652e+20
        @test weight_conversion("metric-ton", "kilogram", 2) == 2000
        @test weight_conversion("metric-ton", "gram", 2) ≈ 2000000
        @test weight_conversion("metric-ton", "milligram", 3) ≈ 3000000000
        @test weight_conversion("metric-ton", "metric-ton", 2) ≈ 2.0
        @test weight_conversion("metric-ton", "long-ton", 3) ≈ 2.9526219
        @test weight_conversion("metric-ton", "short-ton", 2) ≈ 2.2046244
        @test weight_conversion("metric-ton", "pound", 3) ≈ 6613.8732606
        @test weight_conversion("metric-ton", "ounce", 4) ≈ 141095.96289199998
        @test weight_conversion("metric-ton", "carrat", 4) ≈ 20000000
        @test weight_conversion("metric-ton", "atomic-mass-unit", 1) ≈ 6.022136652e+29
        @test weight_conversion("long-ton", "kilogram", 4) == 4064.18432
        @test weight_conversion("long-ton", "gram", 4) == 4064184.32
        @test weight_conversion("long-ton", "milligram", 3) == 3048138240.0
        @test weight_conversion("long-ton", "metric-ton", 4) ≈ 4.06418432
        @test weight_conversion("long-ton", "long-ton", 3) == 2.999999907217152
        @test weight_conversion("long-ton", "short-ton", 1) == 1.119999989746176
        @test weight_conversion("long-ton", "pound", 3) == 6720.000000049448
        @test weight_conversion("long-ton", "ounce", 1) == 35840.000000060514
        @test weight_conversion("long-ton", "carrat", 4) == 20320921.599999998
        @test weight_conversion("long-ton", "atomic-mass-unit", 4) == 2.4475073353955697e+30
        @test weight_conversion("short-ton", "kilogram", 3) == 2721.5519999999997
        @test weight_conversion("short-ton", "gram", 3) == 2721552.0
        @test weight_conversion("short-ton", "milligram", 1) == 907184000.0
        @test weight_conversion("short-ton", "metric-ton", 4) ≈ 3.628736
        @test weight_conversion("short-ton", "long-ton", 3) == 2.6785713457296
        @test weight_conversion("short-ton", "short-ton", 3) == 2.9999999725344
        @test weight_conversion("short-ton", "pound", 2) == 4000.0000000294335
        @test weight_conversion("short-ton", "ounce", 4) == 128000.00000021611
        @test weight_conversion("short-ton", "carrat", 4) == 18143680.0
        @test weight_conversion("short-ton", "atomic-mass-unit", 1) == 5.463186016507968e+29
        @test weight_conversion("pound", "kilogram", 4) == 1.814368
        @test weight_conversion("pound", "gram", 2) == 907.184
        @test weight_conversion("pound", "milligram", 3) == 1360776.0
        @test weight_conversion("pound", "metric-ton", 3) ≈ 0.001360776
        @test weight_conversion("pound", "long-ton", 2) == 0.0008928571152432
        @test weight_conversion("pound", "short-ton", 1) == 0.0004999999954224
        @test weight_conversion("pound", "pound", 3) == 3.0000000000220752
        @test weight_conversion("pound", "ounce", 1) == 16.000000000027015
        @test weight_conversion("pound", "carrat", 1) == 2267.96
        @test weight_conversion("pound", "atomic-mass-unit", 4) == 1.0926372033015936e+27
        @test weight_conversion("ounce", "kilogram", 3) == 0.0850485
        @test weight_conversion("ounce", "gram", 3) == 85.0485
        @test weight_conversion("ounce", "milligram", 4) == 113398.0
        @test weight_conversion("ounce", "metric-ton", 4) ≈ 0.000113398
        @test weight_conversion("ounce", "long-ton", 4) == 0.0001116071394054
        @test weight_conversion("ounce", "short-ton", 4) == 0.0001249999988556
        @test weight_conversion("ounce", "pound", 1) == 0.0625000000004599
        @test weight_conversion("ounce", "ounce", 2) == 2.000000000003377
        @test weight_conversion("ounce", "carrat", 1) == 141.7475
        @test weight_conversion("ounce", "atomic-mass-unit", 1) == 1.70724563015874e+25
        @test weight_conversion("carrat", "kilogram", 1) == 0.0002
        @test weight_conversion("carrat", "gram", 4) == 0.8
        @test weight_conversion("carrat", "milligram", 2) == 400.0
        @test weight_conversion("carrat", "metric-ton", 2) ≈ 4.0000000000000003e-07
        @test weight_conversion("carrat", "long-ton", 3) == 5.9052438e-07
        @test weight_conversion("carrat", "short-ton", 4) == 8.818497600000002e-07
        @test weight_conversion("carrat", "pound", 1) == 0.00044092488404000004
        @test weight_conversion("carrat", "ounce", 2) == 0.0141095962892
        @test weight_conversion("carrat", "carrat", 4) == 4.0
        @test weight_conversion("carrat", "atomic-mass-unit", 4) == 4.8177093216e+23
        @test weight_conversion("atomic-mass-unit", "kilogram", 4) == 6.642160796e-27
        @test weight_conversion("atomic-mass-unit", "gram", 2) == 3.321080398e-24
        @test weight_conversion("atomic-mass-unit", "milligram", 2) == 3.3210803980000002e-21
        @test weight_conversion("atomic-mass-unit", "metric-ton", 3) ≈ 4.9816205970000004e-30
        @test weight_conversion("atomic-mass-unit", "long-ton", 3) == 4.9029473573977584e-30
        @test weight_conversion("atomic-mass-unit", "short-ton", 1) == 1.830433719948128e-30
        @test weight_conversion("atomic-mass-unit", "pound", 3) == 1.0982602420317504e-26
        @test weight_conversion("atomic-mass-unit", "ounce", 2) == 1.1714775914938915e-25
        @test weight_conversion("atomic-mass-unit", "carrat", 2) == 1.660540199e-23
        @test weight_conversion("atomic-mass-unit", "atomic-mass-unit", 2) == 1.999999998903455
    end

    @testset "Conversions: Temparature Conversions" begin
        @test celsius_to_fahrenheit(273.354, 3) == 524.037
        @test celsius_to_fahrenheit(273.354, 0) == 524.0
        @test celsius_to_fahrenheit(-40.0) == -40.0
        @test celsius_to_fahrenheit(-20.0) == -4.0
        @test celsius_to_fahrenheit(0) == 32.0
        @test celsius_to_fahrenheit(20) == 68.0

        @test celsius_to_kelvin(273.354, 3) == 546.504
        @test celsius_to_kelvin(273.354, 0) == 547.0
        @test celsius_to_kelvin(0.0) == 273.15
        @test celsius_to_kelvin(20.0) == 293.15

        @test fahrenheit_to_celsius(273.354, 3) == 134.086
        @test fahrenheit_to_celsius(273.354, 0) == 134.0
        @test fahrenheit_to_celsius(0.0) == -17.78
        @test fahrenheit_to_celsius(20.0) == -6.67
        @test fahrenheit_to_celsius(40.0) == 4.44
        @test fahrenheit_to_celsius(60.0) == 15.56
        @test fahrenheit_to_celsius(80.0) == 26.67

        @test fahrenheit_to_kelvin(273.354, 3) == 407.236
        @test fahrenheit_to_kelvin(273.354, 0) == 407.0
        @test fahrenheit_to_kelvin(0) == 255.37
        @test fahrenheit_to_kelvin(20.0) == 266.48
        @test fahrenheit_to_kelvin(40.0) == 277.59
        @test fahrenheit_to_kelvin(60.0) == 288.71
        @test fahrenheit_to_kelvin(80.0) == 299.82

        @test kelvin_to_celsius(273.354, 3) == 0.204
        @test kelvin_to_celsius(273.354, 0) == 0.0
        @test kelvin_to_celsius(273.15) == 0.0
        @test kelvin_to_celsius(300) == 26.85

        @test kelvin_to_fahrenheit(273.354, 3) == 32.367
        @test kelvin_to_fahrenheit(273.354, 0) == 32.0
        @test kelvin_to_fahrenheit(273.15) == 32.0
        @test kelvin_to_fahrenheit(300) == 80.33
    end

    @testset "Conversions: Length Conversions" begin
        @test length_conversion(50, "mm", "cm") == 5
        @test length_conversion(500, "mm", "m") == 0.5
        @test length_conversion(4321, "mm", "km") == 0.004321
        @test length_conversion(10, "meter", "in") ≈ 393.701
        @test length_conversion(10, "m", "feet") ≈ 32.8084
        @test length_conversion(10, "m", "yd") ≈ 10.9361
        @test length_conversion(1, "in", "mm") ≈ 25.4
        @test length_conversion(1, "in", "cm") ≈ 2.54
        @test length_conversion(5, "in", "m")  ≈ 0.127
        @test length_conversion(3, "ft", "mm") ≈ 914.4
        @test length_conversion(10, "mi", "km") ≈ 16.0934
    end
end

max_stars_repo_name: standardgalactic/Julia-algorithms


text:
struct PseudoPotential{T,relativistic} <: AbstractPseudoPotential{T}
    name::String
    gst_config::Configuration{Orbital{Int}}
    Q::Int
    Vℓ::Vector{GaussianExpansion{T}} # ℓ ∈ 0:ℓmax
    Vℓ′::Vector{GaussianExpansion{T}} # ℓ′ ∈ 1:ℓmax′
    reference::String
end

Base.:(==)(a::PseudoPotential{<:Any,ra}, b::PseudoPotential{<:Any,rb}) where {ra,rb} =
    ra == rb && a.gst_config == b.gst_config && a.Q == b.Q &&
    a.Vℓ == b.Vℓ && a.Vℓ′ == b.Vℓ′

Base.hash(pp::PseudoPotential{<:Any,relativistic}, h::UInt) where relativistic =
    hash((relativistic, pp.gst_config, pp.Q, pp.Vℓ, pp.Vℓ′), h)

const ScalarSORelativisticPseudoPotential{T} = PseudoPotential{T,true}

AtomicPotentials.charge(pp::PseudoPotential) = num_electrons(pp.gst_config)
AtomicPotentials.ground_state(pp::PseudoPotential) = pp.gst_config

function Base.show(io::IO, pp::PseudoPotential{T,relativistic}) where {T,relativistic}
    if relativistic
        write(io, "Scalar-relativistic + spin–orbit pseudo-potential")
    else
        write(io, "Pseudo-potential")
    end
    write(io, " for $(pp.name) ($(pp.gst_config)), Z = $(charge(pp))")
end

function Base.show(io::IO, ::MIME"text/plain", pp::PseudoPotential)
    show(io, pp)
    println(io)
    ℓmax = length(pp.Vℓ)-1
    ℓmax′ = length(pp.Vℓ′)
    print(io, "Long-range Q = $(pp.Q), ℓ ∈ 0:$(ℓmax)")
    ℓmax′ > 0 && print(io, ", ℓ′ ∈ 1:$(length(pp.Vℓ′))")
    println(io, "\nData from \"$(pp.reference)\"")
    ℓk = map(enumerate(0:ℓmax)) do (i,ℓ)
        nk = length(pp.Vℓ[i])
        [[spectroscopic_label(ℓ);repeat([""],nk-1)] 1:nk]
    end |> v -> vcat(v...)
    ℓ′k = map(enumerate(1:ℓmax′)) do (i,ℓ′)
        nk = length(pp.Vℓ′[i])
        [[spectroscopic_label(ℓ′);repeat([""],nk-1)] 1:nk]
    end |> v -> vcat(v...)
    data = [ℓk vcat([[V.n V.β V.B]
                     for V in pp.Vℓ]...)]
    headers = ["ℓ", "k", "n", "β", "B"]
    if ℓmax′ > 0
        data = [data vcat(repeat([""], length(pp.Vℓ[1]), 5),
                          [ℓ′k vcat([[V.n V.β V.B]
                                     for V in pp.Vℓ′]...)],
                          repeat([""], 1, 5))]
        headers = vcat(headers, ["ℓ′", "k", "n", "β", "B"])
    end
    pretty_table(io, data, headers)
end

#=
\[V_{\mathrm{pp}}(\vec{r})=
-\frac{Q}{r}+
\sum_{\ell jk}
B_{\ell j k}
r^{n_k-2}
\exp(-\beta_{\ell j k} r^2)
\mathcal{P}_{\ell j}\]
=#

function (pp::ScalarSORelativisticPseudoPotential{T})(orb::SpinOrbital, r::AbstractVector{T}) where T
    # This returns the diagonal, average part of the relativistic
    # pseudopotential.
    V = -pp.Q./r

    ℓ = orb.orb.ℓ
    ℓ >= length(pp.Vℓ) && return V

    V += pp.Vℓ[ℓ+1](r)

    V
end

function spin_orbit_potential(pp::ScalarSORelativisticPseudoPotential{T}, r::AbstractVector{T},
                              a::SpinOrbital{<:Orbital}, b::SpinOrbital{<:Orbital}) where T
    # This returns the (off-)diagonal, spin–orbit part of the
    # relativistic pseudopotential.
    s = half(1)
    ℓ,mℓa = jmⱼ(a)
    mas  = spin(a)
    ℓb,mℓb = jmⱼ(b)
    mbs = spin(b)

    @assert ℓ == ℓb

    V = zeros(T, length(r))
    (ℓ > length(pp.Vℓ′) || iszero(ℓ)) && return V

    # Eq. (30) of
    # - Dolg, M. (2016). Relativistic Effective Core Potentials. In (Eds.),
    #   Handbook of Relativistic Quantum Chemistry (pp. 449–478). : Springer
    #   Berlin Heidelberg.
    coeff = ℓ*clebschgordan(
        ℓ,mℓa,s,mas,ℓ+s
    )*clebschgordan(
        ℓ,mℓb,s,mbs,ℓ+s
    )
    if abs(mℓa+mas) < ℓ
        coeff -= (ℓ+1)*clebschgordan(
            ℓ,mℓa,s,mas,ℓ-s
        )*clebschgordan(
            ℓ,mℓb,s,mbs,ℓ-s
        )
    end

    # According to
    # http://www.tc.uni-koeln.de/cgi-bin/pp.pl?language=en,job=getreadme
    # the spin–orbit potentials are pre-multiplied by 2/(2ℓ+1), but we
    # need 1/(2ℓ+1).
    coeff /= 2

    # We can index by ℓ, since Vℓ′ does not contain any entry for ℓ = s
    V += coeff*pp.Vℓ′[ℓ](r)

    V
end

function (pp::PseudoPotential{T,false})(orb::SpinOrbital, r::AbstractVector{T}) where T
    V = -pp.Q./r
    ℓ = orb.orb.ℓ
    ℓ >= length(pp.Vℓ) && return V
    V += pp.Vℓ[ℓ+1](r)

    V
end

(pp::PseudoPotential{T})(orb::AbstractOrbital, r::T) where T = pp(orb, [r])[1]


text:
using Test, YaoBlocks, YaoArrayRegister

struct MockedTag{BT, N} <: TagBlock{BT, N}
    content::BT

    MockedTag(x::BT) where {N, BT <: AbstractBlock{N}} = new{BT, N}(x)
end

@test nqubits(MockedTag(X)) == nqubits(X)
@test nqubits(MockedTag(kron(X, Y))) == nqubits(kron(X, Y))

@test getiparams(MockedTag(phase(0.1))) == ()
@test getiparams(MockedTag(cache(phase(0.1)))) == ()
@test getiparams(MockedTag(Rx(0.1))) == ()

@test parameters(MockedTag(phase(0.1))) == [0.1]
@test parameters(MockedTag(cache(phase(0.1)))) == [0.1]
@test parameters(MockedTag(Rx(0.1))) == [0.1]

@test occupied_locs(MockedTag(chain(3, put(1=>X), put(3=>X)))) == occupied_locs(chain(3, put(1=>X), put(3=>X)))

@testset "scale apply" begin
    # factor
    @test factor(-X) == -1
    @test factor(-2*X) == -2

    # type
    @test -X isa Scale{<:Val}
    @test -im*X isa Scale{<:Val}
    @test -2*X isa Scale{<:Number}

    # apply!
    reg = rand_state(1)
    @test apply!(copy(reg), -X) ≈ -apply!(copy(reg), X)
    @test apply!(copy(reg), (-im*Y)') ≈ apply!(copy(reg), im*Y)
    @test apply!(copy(reg), (-2im*Y)') ≈ apply!(copy(reg), 2im*Y)

    # mat
    @test mat(-X) ≈ -mat(X)
    @test mat((-im*Y)') ≈ (-im*mat(Y))'
    @test mat((-2im*Y)') ≈ (-2im*mat(Y))'

    # other
    @test chsubblocks(-X, Y) == -Y
    @test cache_key(Scale(Val(-1), X)) == cache_key(Scale(-1, X))
    @test copy(-X) == -X
    @test copy(-X) isa Scale{<:Val}
end


text:
xdoc = XMLDocument()

xroot = create_root(xdoc, "States")

xs1 = new_child(xroot, "State")
add_cdata(xdoc, xs1, "Massachusetts")

rtxt = """
<?xml version="1.0" encoding="utf-8"?>
<States>
  <State><![CDATA[Massachusetts]]></State>
</States>
"""

@test strip(string(xdoc)) == strip(rtxt)

free(xdoc)


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

pids = addprocs_with_testenv(4; topology="master_slave")

let p1 = pids[1], p2 = pids[2]
    @test_throws RemoteException remotecall_fetch(()->remotecall_fetch(myid, p2), p1)
end

function test_worker_counts()
    # check if the nprocs/nworkers/workers are the same on the remaining workers
    np=nprocs()
    nw=nworkers()
    ws=sort(workers())

    for p in workers()
        @test (true, true, true) == remotecall_fetch(p, np, nw, ws) do x,y,z
            (x==nprocs(), y==nworkers(), z==sort(workers()))
        end
    end
end

function remove_workers_and_test()
    while nworkers() > 0
        rmprocs(workers()[1])
        test_worker_counts()
        if nworkers() == nprocs()
            break
        end
    end
end

remove_workers_and_test()

# connect even pids to other even pids, odd to odd.
mutable struct TopoTestManager <: ClusterManager
    np::Integer
end

function Base.launch(manager::TopoTestManager, params::Dict, launched::Array, c::Condition)
    dir = params[:dir]
    exename = params[:exename]
    exeflags = params[:exeflags]

    cmd = `$exename $exeflags --bind-to $(Base.Distributed.LPROC.bind_addr) --worker`
    cmd = pipeline(detach(setenv(cmd, dir=dir)))
    for i in 1:manager.np
        io = open(cmd, "r+")
        Base.Distributed.write_cookie(io)

        wconfig = WorkerConfig()
        wconfig.process = io
        wconfig.io = io.out
        wconfig.ident = i
        wconfig.connect_idents = collect(i+2:2:manager.np)
        push!(launched, wconfig)
    end

    notify(c)
end

const map_pid_ident=Dict()
function Base.manage(manager::TopoTestManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :register
        map_pid_ident[id] = get(config.ident)
    elseif op == :interrupt
        kill(get(config.process), 2)
    end
end

addprocs_with_testenv(TopoTestManager(8); topology="custom")

while true
    if any(x->get(map_pid_ident, x, 0)==0, workers())
        yield()
    else
        break
    end
end

let p1, p2
for p1 in workers()
    for p2 in workers()
        i1 = map_pid_ident[p1]
        i2 = map_pid_ident[p2]
        if (iseven(i1) && iseven(i2)) || (isodd(i1) && isodd(i2))
            @test p2 == remotecall_fetch(p->remotecall_fetch(myid, p), p1, p2)
        else
            @test_throws RemoteException remotecall_fetch(p->remotecall_fetch(myid, p), p1, p2)
        end
    end
end
end

remove_workers_and_test()

# test `lazy` connection setup
function def_count_conn()
    @everywhere function count_connected_workers()
        count(x -> isa(x, Base.Distributed.Worker) && isdefined(x, :r_stream) && isopen(x.r_stream),
                Base.Distributed.PGRP.workers)
    end
end

addprocs_with_testenv(8)
def_count_conn()

# Test for 10 random combinations
wl = workers()
combinations = []
while length(combinations) < 10
    from = rand(wl)
    to = rand(wl)
    if from == to || ((from,to) in combinations) || ((to,from) in combinations)
        continue
    else
        push!(combinations, (from,to))
    end
end

# Initially only master-slave connections ought to be setup
expected_num_conns = 8
let num_conns = sum(asyncmap(p->remotecall_fetch(count_connected_workers,p), workers()))
    @test num_conns == expected_num_conns
end

for (i, (from,to)) in enumerate(combinations)
    remotecall_wait(topid->remotecall_fetch(myid, topid), from, to)
    global expected_num_conns += 2    # one connection endpoint on both from and to
    let num_conns = sum(asyncmap(p->remotecall_fetch(count_connected_workers,p), workers()))
        @test num_conns == expected_num_conns
    end
end

# With lazy=false, all connections ought to be setup during `addprocs`
rmprocs(workers())
addprocs_with_testenv(8; lazy=false)
def_count_conn()
@test sum(asyncmap(p->remotecall_fetch(count_connected_workers,p), workers())) == 64

# Cannot add more workers with a different `lazy` value
@test_throws ArgumentError addprocs_with_testenv(1; lazy=true)

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Wasmer
using Test

@testset "" begin
    engine = WasmEngine()
    store = WasmStore(engine)
    code = wat"""
    (module
        (func $add (param $lhs i32) (param $rhs i32) (result i32)
            local.get $lhs
            local.get $rhs
            i32.add)
        (export "add" (func $add))
    )
    """
    modu = WasmModule(store, code)

    instance = WasmInstance(store, modu)
    add = exports(instance).add

    res = add(Int32(1), Int32(2))

    @test res == Int32(3)
end


text:
module SudokuSolver
    using Statistics
    export show_sudoku, sudoku_solver, block, fixed_values, nr_errors, fill_full, fixed_blocks, start_temp, decide

    """
        show_sudoku(sudoku)

    Shows sudoku with grid lines.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> show_sudoku(sudoku)
    741|59.|382
    3.6|247|195
    295|183|6.7
    -----------
    972|..8|.14
    ..4|.5.|8..
    15.|4..|736
    -----------
    5.3|972|468
    469|835|2.1
    827|.14|953
    ```
    """
    function show_sudoku(sudoku)
    	for i in 1:9
	    	line = ""

		    if i == 4 || i == 7
			    println("-----------")
    		end
	    	for j in 1:9
		    	if j == 4 || j == 7
			    	line = line * "|"
    			end
                if sudoku[i,j] == 0
                    line = line * '.'
                else
	    		    line = line * string(sudoku[i,j])
                end
		    end
    		println(line)
	    end
    end

    """
        block(row, column, sudoku)

    Returns the block a certain position is located in. With blocks as the 3x3 subsections of the sudoku between the gridlines.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> block(1, 1, sudoku)
    3x3 Array{Int64,2}:
    741
    306
    295
    ```
    """
    function block(row, column, sudoku)
        # Michiel: adding the following checks gives errers, there is something wrong here!
        #@assert 1 ≤ row ≤ 3 BoundsError()
        #@assert 1 ≤ column ≤ 3 BoundsError()
        # Michiel: you can do this much more concisely, e.g. `range(3(row-1)+1, 3row, step=1)`
        # if you use `@view` you also don't take a copy of the block (more memory efficient)
	    if row in [1, 2, 3]
    		block = sudoku[1:3, :]
    	elseif row in [4 5 6]
	    	block = sudoku[4:6, :]
	    else
		    block = sudoku[7:9, :]
	    end

    	if column in [1, 2, 3]
	    	block = block[:, 1:3]
    	elseif column in [4 5 6]
	    	block = block[:, 4:6]
    	else
	    	block = block[:, 7:9]
    	end

	    return block
    end

    """
        freq_table(row)

    Makes a frequency table of a given row, column or block. Row can be any array.

    ## Examples

    ```julia-repl
    julia> row = [9 7 2 0 0 8 0 1 4]

    julia> freq_table(row)
    Dict{Any, Any}(
    7 => 1
    9 => 1
    0 => 3
    4 => 1
    2 => 1
    8 => 1
    1 => 1
    )

    ```
    """
    function freq_table(row)
    	freqs = Dict()
	    for value in row
		    if haskey(freqs, value)
			    freqs[value] += 1
    		else
	    		freqs[value] = 1
		    end
    	end
	    return freqs
    end

    """
        duplicates(row)

    Searches for duplicates in a row, column or block. 0 values are ignored. Returns true if a duplicate is found.

    ## Examples

    ```julia-repl
    julia> row = [9 7 2 0 0 8 0 1 4]

    julia> duplicates(row)
    false

    julia> row2 = [9 7 2 0 9 8 0 1 4]

    julia> duplicates(row2)
    true
    ```
    """
    function duplicates(row)
        table = freq_table(row)
        for i in 1:9
            if haskey(table, i)
                if table[i] > 1
                    return true
                end
            end
        end
        return false
    end

    """
        valid(sudoku)

    Checks if a given array is a valid sudoku. Sudoku should be a 9x9 array, blanks are filled with zeros and rows, columns and blocks do not have duplicate values other than zero.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
    3 0 6 2 4 7 1 9 5
    2 9 5 1 8 3 6 0 7
    9 7 2 0 0 8 0 1 4
    0 0 4 0 5 0 8 0 0
    1 5 0 4 0 0 7 3 6
    5 0 3 9 7 2 4 6 8
    4 6 9 8 3 5 2 0 1
    8 2 7 0 1 4 9 5 3]

    julia> valid(sudoku)
    true

    julia> valid(3)
    false
    ```
    """
    function valid(sudoku)
        if size(sudoku) != (9, 9) # 9x9 array
            return false
        end

        for entry in keys(freq_table(sudoku)) # no values other than 0:9
            if !(entry in 0:9)
                return false
            end
        end

        for i in 1:9 # no repeated numbers in rows and columns
            if duplicates(sudoku[i,:])
                return false
            elseif duplicates(sudoku[:,i])
                return false
            end
        end

        for i in [1, 4, 7] # no repeats in blocks
            for j in [1, 4, 7]
                if duplicates(block(i,j,sudoku))
                    return false
                end
            end
        end

        return true
    end

    """
        fixed_values(sudoku)

    Returns an array indicating which values the sudoku started out whith (fixed values).
    Fixed positions are indicated with a 1, non-fixed values with a 0.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> fixed_values(sudoku)
    9×9 Array{Int64,2}:
     1  1  1  1  1  0  1  1  1
     1  0  1  1  1  1  1  1  1
     1  1  1  1  1  1  1  0  1
     1  1  1  0  0  1  0  1  1
     0  0  1  0  1  0  1  0  0
     1  1  0  1  0  0  1  1  1
     1  0  1  1  1  1  1  1  1
     1  1  1  1  1  1  1  0  1
     1  1  1  0  1  1  1  1  1
    ```
    """
    function fixed_values(sudoku)
    	fixed = zero(sudoku)
	    for i in 1:9
		    for j in 1:9
			    if sudoku[i, j] != 0
				    fixed[i, j] = 1
    			end
	    	end
    	end
	    return fixed
    end

    """
        nr_errors(sudoku)

    Calculates the number of errors (repeated values in a row or column) are in a sudoku.
    Assumes that blocks have no repeat values.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7  4  1  5  9  6  3  8  2
        3  8  6  2  4  7  1  9  5
        2  9  5  1  8  3  6  4  7
        9  7  2  6  3  8  2  1  4
        3  6  4  1  5  7  8  5  9
        1  5  8  4  2  9  7  3  6
        5  1  3  9  7  2  4  6  8
        4  6  9  8  3  5  2  7  1
        8  2  7  6  1  4  9  5  3]

    julia> nr_errors(sudoku)
    10
    ```
    """
    function nr_errors(sudoku)
	    errors = 0
    	for i in 1:9
	    	errors += (9 - length(freq_table(sudoku[i,:])))
		    errors += (9 - length(freq_table(sudoku[:,i])))
    	end
	    return errors
    end

    """
        missing_nr(block)

    Finds numbers (1-9) not yet in block.

    ## Examples

    ```julia-repl
    julia> block =  [7  4  1
    3  0  6
    2  9  5]

    julia> missing_nr(block)
    1-element Array{Any,1}:
    8
    ```
    """
    function missing_nr(block)
    	freqs = freq_table(block)
	    values =  []
    	for value in 1:9
	    	if !haskey(freqs, value)
		    	append!(values, value)
    		end
	    end
	    return values
    end

    """
        pick(list)

    Returns a random value from list(value) and the input list without the picked value (list).

    ## Examples

    ```julia-repl
    julia> list = [9, 7, 2, 0, 9, 8, 0, 1, 4]

    julia> pick(list)
    (2, [9, 7, 0, 9, 8, 0, 1, 4])
    ```
    """
    function pick(list)
		i = rand(1:length(list))
    	value = popat!(list, i)
	    return value, list
    end

    """
    fill_block(block)

    Fills a block randomly with non-repeated numbers from 1 to 9.

    ## Examples

    ```julia-repl
    julia> block =  [7  4  1
    3  0  6
    2  9  5]

    julia> fill_block(block)
    3×3 Array{Int64,2}:
     7  4  1
     3  8  6
     2  9  5
    ```
    """
    function fill_block(block)
    	values = missing_nr(block)
	    for (i, v) in enumerate(block)
		    if v == 0
			    nr, values = pick(values)
    			block[i] = nr
	    	end
    	end
    	return block
    end

    """
        fill_full(sudoku)

    Fills all empty spaces of a sudoku at random (numbers between 1 and 9). Numbers are not repeated within a block, but there are repeats in rows and columns.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> fill_full(sudoku)
    9×9 Array{Int64,2}:
     7  4  1  5  9  6  3  8  2
     3  8  6  2  4  7  1  9  5
     2  9  5  1  8  3  6  4  7
     9  7  2  6  3  8  2  1  4
     3  6  4  1  5  7  8  5  9
     1  5  8  4  2  9  7  3  6
     5  1  3  9  7  2  4  6  8
     4  6  9  8  3  5  2  7  1
     8  2  7  6  1  4  9  5  3
    ```
    """
    function fill_full(sudoku)
    	filled = zero(sudoku)
	    for i in [1, 4, 7]
		    for j in [1, 4, 7]
			    fullblock = fill_block(block(i,j,sudoku))
    			filled[i:i+2, j:j+2] = fullblock
	    	end
    	end
	    return filled
    end

    """
        not_fixed(fixed)

    Returns an array of positions with non-fixed values. Non-fixed values are the blank (0) values in the start sudoku. Starts from the result of the fixed_values function.

    ## Examples

    ```julia-repl
    julia> fixed = [ 1  1  1  1  1  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  0  1  0  1  1
        0  0  1  0  1  0  1  0  0
        1  1  0  1  0  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  1  1  1  1  1]

    julia> not_fixed(fixed)
    18-element Array{Any,1}:
     5
     11
     14
     16
     24
     31
     32
     36
     40
     42
     46
     50
     51
     58
     66
     68
     71
     77
    ```
    """
    function not_fixed(fixed)
        nonfixed = []
	    for (i, v) in enumerate(fixed)
		    if v == 0
			    append!(nonfixed, i)
    		end
	    end
    	return nonfixed
    end

    """
        block_nr_ind(nr)

    Blocks are numbered from 1 to 9, from top to bottom and left to right.
    This function returns an index within that block. Row and column will be 3, 6 or 9.

    ## Examples

    ```julia-repl
    julia> block_nr_ind(3)
    (9, 3)
    ```
    """
    function block_nr_ind(nr)
        row = (mod(nr-1, 3) + 1) * 3
        column = (div(nr-1, 3) + 1) * 3
        return row, column
    end

    """
        fixed_blocks(fixed)

    Lists blocks that have more than 2 non-fixed positions. Blocks are numbered from 1 to 9, from top to bottom and left to right.

    ## Examples

    ```julia-repl
    julia> fixed = [ 1  1  1  1  1  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  0  1  0  1  1
        0  0  1  0  1  0  1  0  0
        1  1  0  1  0  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  1  1  1  1  1]

    julia> fixed_blocks(fixed)
    Any
    2
    5
    8
    ```
    """
    function fixed_blocks(fixed)
        blocks = Int[]
        for i in 1:9
            row, col = block_nr_ind(i)
            # Michiel: you might make this faster using the function `count`
            nonfixed = not_fixed(block(row, col, fixed))
            if length(nonfixed) > 2
                append!(blocks, i)
            end
        end
        return blocks
    end

    """
        possible_swap(fixed)

    Generates a list of blocks with more than two non-fixed positions and a dictionary of
    possible positions to swap. Keys are the block numbers used in fixed_blocks.

    ## Examples

    ```julia-repl
    julia> fixed = [ 1  1  1  1  1  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  0  1  0  1  1
        0  0  1  0  1  0  1  0  0
        1  1  0  1  0  0  1  1  1
        1  0  1  1  1  1  1  1  1
        1  1  1  1  1  1  1  0  1
        1  1  1  0  1  1  1  1  1]

    julia> possible_swap(fixed)
    (Any[2, 5, 8], Dict{Any,Any}(2 => Any[2, 5, 9],5 => Any[1, 2, 4, 6, 8, 9],8 => Any[1, 5, 8]))
    ```
    """
    function possible_swap(fixed)
        blocks = fixed_blocks(fixed)
        possible = Dict{Int, Vector{Int}}()  #michiel: Type stability
        for i in blocks
            row, col = block_nr_ind(i)
            possible[i] = not_fixed(block(row, col, fixed))
        end
        return blocks, possible
    end

    """
        swap(sudoku, blocks, possible)

    Randomly swaps the values of two non-fixed positions within the same (random) block. Generates a new state.
    ```
    """
    function swap(sudoku, blocks, possible)
    	blocknr = rand(blocks) # random block where positions can be swapped
        row, col = block_nr_ind(blocknr)
        nrblock = block(row, col, sudoku)
        nonfixed = copy(possible[blocknr])
    	pos1, nonfixed = pick(nonfixed) # Select random positions
	    pos2, nonfixed = pick(nonfixed)
    	temp = nrblock[pos1] # Swap values
	    nrblock[pos1] = nrblock[pos2]
    	nrblock[pos2] = temp
	    sudoku[row-2:row, col-2:col] = nrblock # Replace selected block in sudoku
    	return sudoku
    end

    """
        start_temp(sudoku)

    Chooses a starting temperature by taking the standard deviation of the errors of 50 starting states.
    Returns the best of these stating states as a starting point.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> start_temp(sudoku)
    (
    3.0445
    9×9 Array{Int64,2}:
    9  5  1  2  7  3  8  4  3
    3  8  4  1  5  8  2  7  6
    6  7  2  6  4  9  5  1  9
    1  3  9  7  6  8  4  2  5
    2  8  4  3  1  5  9  6  7
    5  6  7  4  9  2  1  3  8
    4  9  3  7  5  8  6  5  2
    8  2  6  6  3  4  7  9  1
    7  1  5  9  2  1  3  8  4
    )
    ```
    """
    function start_temp(sudoku)
    	errorscores = []
        bestscore = 80
        best = []
    	for i in 1:50
    		state = fill_full(sudoku)
            score = nr_errors(state)
            if score < bestscore
                bestscore = score
               best = state
            end
		    append!(errorscores, nr_errors(state))
	    end
	    return std(errorscores), best
    end

    """
        decide(diff, T)

    Decides whether a new state is used or discarded. Diff is the difference in amount of errors between the new and old states.
    """
    decide(diff, T) = rand() < exp(-diff/T)

    """
        sudoku_solver(sudoku)

    Solves a given sudoku using Simulated Annealing. Sudoku should be a 9x9 array, blanks are filled with zeros and rows, columns and blocks do not have duplicate values other than zero.

    ## Examples

    ```julia-repl
    julia> sudoku =  [7 4 1 5 9 0 3 8 2
	    3 0 6 2 4 7 1 9 5
    	2 9 5 1 8 3 6 0 7
    	9 7 2 0 0 8 0 1 4
    	0 0 4 0 5 0 8 0 0
    	1 5 0 4 0 0 7 3 6
    	5 0 3 9 7 2 4 6 8
    	4 6 9 8 3 5 2 0 1
    	8 2 7 0 1 4 9 5 3]

    julia> sudoku_solver(sudoku)
    741|59.|382
    3.6|247|195
    295|183|6.7
    -----------
    972|.08|.14
    ..4|.5.|8..
    15.|4..|736
    -----------
    5.3|972|468
    469|835|2.1
    827|.14|953
    741|596|382
    386|247|195
    295|183|647
    -----------
    972|368|514
    634|751|829
    158|429|736
    -----------
    513|972|468
    469|835|271
    827|614|953
    ```
    """
    function sudoku_solver(sudoku, cooling = 0.99)
        if !valid(sudoku)
            error("Input is not a valid sudoku")
        end
    	solved = 0
    	stuck = 0
        restart = 0
    	show_sudoku(sudoku)
	    fixed = fixed_values(sudoku)
    	T0, filled = start_temp(sudoku)
        T = T0
        blocks, possible = possible_swap(fixed)
	    errors = nr_errors(filled)
    	iterations = 81 - sum(fixed)
	    state = filled
        best = filled
        bestscore = 90

    	if errors == 0
	    	solved = 1
	    end

    	while solved == 0
	    	lastscore = errors
		    for i in 1:iterations
			    newstate = swap(state, blocks, possible)
    			newerrors = nr_errors(newstate)

	    		if newerrors == 0
		    		solved = 1
			    	break
    			end

	    		diff = newerrors - errors
		    	if diff <= 0
			    	state = newstate
				    errors = newerrors
                    if errors < bestscore
                        best = state
                        bestscore = errors
                    end
    			else
	    			kept = decide(diff, T)
		    		if kept
			    		state = newstate
				    	errors = newerrors
    				end
	    		end
		    end

    		T *= cooling
	    	if errors == 0
		    	solved = 1
			    break
    		end

	    	if errors >= lastscore
			    stuck += 1
		    else
    			stuck = 0
	    	end

            if stuck > 1000 # When running too long:
                restart +=1
                if restart > 20 # Return to best solution
                    state = best
                    errors = bestscore
                    restart = 0
                    stuck = 0
                else
                    T += T0/3 # Or increase temperature
                    stuck = 0
                end
            end
    	end

    	show_sudoku(state)
	    return state
    end
end