text:
# # HK (Hegselmann and Krause) opinion dynamics model

# This is an implementation of a simple version of the
# [Hegselmann and Krause (2002)](http://jasss.soc.surrey.ac.uk/5/3/2.html) model,
# which also features synchronous updating of Agent properties.

# It is a model of opinion formation with the question: which
# parameters' values lead to consensus, polarization or fragmentation?
# It models interacting **groups** of agents (as opposed to interacting pairs, typical in
# the literature) in which it is assumed that if an agent disagrees too much with
# the opinion of a source of influence, the source can no longer influence the
# agent’s opinion. There is then a "bound of confidence". The model shows that the
# systemic configuration is heavily dependent on this parameter's value.

# We implement it as an example of how to implement a [Synchronous update schedule](http://jmckalex.org/compass/syn-and-asynch-expl.html) .
# In a Synchronous update schedule changes made to an agent are not seen by
# other agents until the next clock tick — that is,
# all agents update simultaneously ([Wilensky 2015, p.286](https://mitpress.mit.edu/books/introduction-agent-based-modeling)).

# The model has the following components:

# - A set of n Agents with opinions xᵢ in the range [0,1] as attribute;
# - A bound ϵ in also in the range [0,1] (actually, the range of interesting results is
#   approximately (0, 0.3]);
# - The update rule: at each step every agent adopts the mean of the opinions which are within
#   the confidence bound ( |xᵢ - xⱼ| ≤ ϵ).

# ## Core structures
# We start by defining the Agent type and initializing the model.
# The Agent type has two fields so that we can implement the synchronous update.
using Agents
using Statistics: mean

mutable struct HKAgent{T <: AbstractFloat} <: AbstractAgent
    id::Int
    old_opinion::T
    new_opinion::T
end

function hk_model(;numagents = 100, ϵ = 0.4)
    model = ABM(HKAgent, scheduler = fastest,
                properties = Dict(:ϵ => ϵ))
    for i in 1:numagents
        o = rand()
        add_agent!(model, o, o)
    end
    return model
end

# And some helper functions for the update rule. As there is a filter in
# the rule we implement it outside the `agent_step!` method. Notice that the filter
# is applied to the `:old_opinion` field.
get_old_opinion(agent)::Float64 = agent.old_opinion

function boundfilter(agent,model)
    filter(j->abs(get_old_opinion(agent) - j) < model.properties[:ϵ],
     get_old_opinion.(values(model.agents)))
end

# Now we implement the `agent_step!` and `model_step!` methods.
function agent_step!(agent, model)
    agent.new_opinion = mean(boundfilter(agent,model))
end

function updateold(a)
    a.old_opinion = a.new_opinion
    return a
end

function model_step!(model)
    for i in keys(model.agents)
        agent = id2agent(i, model)
        updateold(agent)
    end
end

# From this implementation we see that to implement synchronous scheduling
# we define an Agent type with an `old` and `new` fields for attributes that
# are changed through synchronous updating. In `agent_step!` we use the `old` field
# and after updating all the agents `new` field we use the `model_step!`
# to update the model  for the next iteration.
# ## Running the model
# Now we can define a method for our simulation run.
# The parameter of interest is the `:new_opinion` field so we assign
# it to variable `agent_properties` and pass it to the `step!` method
# to be collected in a DataFrame.
function model_run(; numagents = 100, iterations = 50, ϵ= 0.05)
    model = hk_model(numagents = numagents, ϵ = ϵ)
    when = 0:5:iterations
    agent_properties = [:new_opinion]
    data = step!(
            model,
            agent_step!,
            model_step!,
            iterations,
            agent_properties,
            when = when
            )
    return(data)
end

data = model_run(numagents = 10, iterations = 20)
data[end-19:end, :]

# Finally we run three scenarios, collect the data and plot it.
using Plots

plotsim(data, ϵ) = plot(
                        data[!, :step],
                        data[!, :new_opinion],
                        leg= false,
                        group = data[!, :id],
                        title = "epsilon = $(ϵ)"
                        )

plt001,plt015,plt03 = map(
                          e -> (model_run(ϵ= e), e) |>
                          t -> plotsim(t[1], t[2]),
                          [0.05, 0.15, 0.3]
                          )

plot(plt001, plt015, plt03, layout = (3,1))


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

using Markdown
using InteractiveUtils

# ╔═╡ 48aa5a70-b65d-11eb-1f4c-ed447341ce4f
begin
	using FFTW, DeconvOptim, TestImages, Napari, Noise
	using IndexFunArrays, FourierTools, Colors
	v = napari.view_image
end

# ╔═╡ ced8e77f-035e-4aae-b236-b5196028861f
md"## Load Packages"

# ╔═╡ 6e88b04a-dcb5-4bc9-a4f9-f5b7dfd9965e
FFTW.set_num_threads(1)

# ╔═╡ bf93f8bd-144b-4c90-892e-01df05956327
begin
	img = Float32.(testimage("resolution_test_512"))
end

# ╔═╡ 35262a15-628e-4639-a0ca-2b53ac7fbecc
md"## Simulate a simple PSF"

# ╔═╡ 594eb15c-7f38-4890-8f26-d35a9fe128ad
function simulate_psf(shape, radius, dtype=Float32)
	arr_ft = ones(dtype, shape)
	aperture = rr(shape) .< radius
	arr = abs2.(ift(arr_ft .* aperture))
	arr ./= sum(arr)
	return arr
end

# ╔═╡ ba837952-d782-46eb-846c-7893401a4d4b
psf = ifftshift(simulate_psf(size(img), 30));

# ╔═╡ ed5e859b-aa48-4d43-8378-367be437f17e
v(psf)

# ╔═╡ 3cf013a8-221d-40ca-91b6-0a48774d9317
md"### Create blurry measured image"

# ╔═╡ 085402ed-d42d-40f5-9f70-6c6c0483aa39
measured = poisson(FourierTools.conv(img, psf), 300)

# ╔═╡ eb362894-9ea4-4a0d-8723-7acd39abd567
v(measured)

# ╔═╡ 4464fa21-69d4-4447-ba27-7299d7e8d454
md"## DeconvOptim.jl"

# ╔═╡ 5f7222dc-c2c5-4fd5-879f-0dddc99c9c6e
reg = TV()

# ╔═╡ 16d79931-5ae0-4515-b1f3-f1b694da80c6
res, o = deconvolution(measured, psf, regularizer=reg, iterations=20, λ=0.02);

# ╔═╡ 01a50c2a-bb57-4625-980d-fd4c9d07044e
v(res)

# ╔═╡ Cell order:
# ╟─ced8e77f-035e-4aae-b236-b5196028861f
# ╠═48aa5a70-b65d-11eb-1f4c-ed447341ce4f
# ╠═6e88b04a-dcb5-4bc9-a4f9-f5b7dfd9965e
# ╠═bf93f8bd-144b-4c90-892e-01df05956327
# ╟─35262a15-628e-4639-a0ca-2b53ac7fbecc
# ╠═594eb15c-7f38-4890-8f26-d35a9fe128ad
# ╠═ba837952-d782-46eb-846c-7893401a4d4b
# ╠═ed5e859b-aa48-4d43-8378-367be437f17e
# ╟─3cf013a8-221d-40ca-91b6-0a48774d9317
# ╠═085402ed-d42d-40f5-9f70-6c6c0483aa39
# ╠═eb362894-9ea4-4a0d-8723-7acd39abd567
# ╟─4464fa21-69d4-4447-ba27-7299d7e8d454
# ╠═5f7222dc-c2c5-4fd5-879f-0dddc99c9c6e
# ╠═16d79931-5ae0-4515-b1f3-f1b694da80c6
# ╠═01a50c2a-bb57-4625-980d-fd4c9d07044e

ext: jl


text:
using MVApp
if length(ARGS) >= 4
    MVApp.julia_main()
else
    println("Input must be start timepoint, end timepoint, and cell IDs")
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Gherkin
using Behavior.Gherkin: ScenarioStep, Background
using Behavior
using Behavior: StepDefinitionContext, StepDefinition, StepDefinitionLocation
using Behavior: Executor, StepExecutionResult, QuietRealTimePresenter, executefeature, ScenarioResult
import Behavior: present

mutable struct FakeRealTimePresenter <: Behavior.RealTimePresenter
    scenarios::Vector{Scenario}
    scenarioresults::Vector{ScenarioResult}
    steps::Vector{ScenarioStep}
    results::Dict{ScenarioStep, StepExecutionResult}
    features::Vector{Feature}

    FakeRealTimePresenter() = new([], [], [], Dict(), [])
end

present(p::FakeRealTimePresenter, scenario::Scenario) = push!(p.scenarios, scenario)
present(p::FakeRealTimePresenter, _scenario::Scenario, result::ScenarioResult) = push!(p.scenarioresults, result)
present(p::FakeRealTimePresenter, step::ScenarioStep) = push!(p.steps, step)
present(p::FakeRealTimePresenter, step::ScenarioStep, result::StepExecutionResult) = p.results[step] = result
present(p::FakeRealTimePresenter, feature::Feature) = push!(p.features, feature)

stepresult(p::FakeRealTimePresenter, step::ScenarioStep) = p.results[step]

@testset "Executor Presentation" begin
    @testset "Execution presentation; Scenario is executed; Scenario is presented" begin
        presenter = FakeRealTimePresenter()
        matcher = FakeStepDefinitionMatcher(Dict())
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[])
        Behavior.executescenario(executor,  Background(),scenario)

        @test presenter.scenarios[1] == scenario
    end

    @testset "Execution presentation; Scenario is executed; ScenarioResult is presented" begin
        presenter = FakeRealTimePresenter()
        matcher = FakeStepDefinitionMatcher(Dict())
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[])
        Behavior.executescenario(executor,  Background(),scenario)

        @test presenter.scenarioresults[1].scenario == scenario
    end

    @testset "Execution presentation; Scenario has on Given; Given is presented" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[given])
        Behavior.executescenario(executor,  Background(),scenario)

        @test presenter.steps[1] == given
    end

    @testset "Execution presentation; Scenario step is successful; Step is presented as successful" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[given])
        Behavior.executescenario(executor,  Background(),scenario)

        @test stepresult(presenter, given) == Behavior.SuccessfulStepExecution()
    end

    @testset "Execution presentation; Scenario step fails; Next is also presented" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        when = When("some action")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                 when => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[given, when])
        Behavior.executescenario(executor,  Background(),scenario)

        @test presenter.steps[2] == when
    end

    @testset "Execution presentation; Scenario step fails; Next has result Skipped" begin
        presenter = FakeRealTimePresenter()
        given = Given("some precondition")
        when = When("some action")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition,
                                                 when => successful_step_definition))
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[given, when])
        Behavior.executescenario(executor,  Background(),scenario)

        @test stepresult(presenter, when) == Behavior.SkippedStep()
    end

    @testset "Execution presentation; Feature is executed; Feature is presented" begin
        presenter = FakeRealTimePresenter()
        matcher = FakeStepDefinitionMatcher(Dict())
        executor = Executor(matcher, presenter)

        scenario = Scenario("Some scenario", String[], ScenarioStep[])
        feature = Feature(FeatureHeader("", [], []), [scenario])
        Behavior.executefeature(executor, feature)

        @test presenter.features[1] == feature
    end
end


text:
using Transformers
using Transformers.Basic
using Test
using Random

import Flux
using Flux: gradient

const tests = [
    "transformer",
    "nntopo",
    "embed",
    "basic",
    "gpt",
    "bert",
]

Random.seed!(0)

if v"1.0.0" <= VERSION  < v"1.4.0"
    @test isempty(detect_ambiguities(Transformers, Transformers.GenerativePreTrain, Transformers.Basic, Transformers.Datasets))
end

if haskey(ENV, "TEST_TRANSFORMERS_PRETRAIN")
    push!(tests, "pretrain")
end

@testset "Transformers" begin
  if Flux.use_cuda[]
      @info "Test CUDA"
    include("test_cuda.jl")
  end

  for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    @info "Test $(uppercase(t))"
    include(fp)
  end
end


text:
export compute

"""
`compute(path::String)`

Call the `compute` function on the configuration file.

Inputs:
======

* path::String - Path to configuration file

"""
function compute(path::String)
    cfg = parse_config(path)
    update_logging!(cfg)
    write_config(cfg)
    T = cfg["precision"] in SINGLE ? Float32 : Float64
    if T == Float32 && (cfg["solver"] in CHOLMOD || cfg["solver"] in MKLPARDISO)
        cswarn("Cholmod & MKLPardiso solver modes work only in double precision. Switching precision to double.")
        T = Float64
    end
    V = cfg["use_64bit_indexing"] in TRUELIST ? Int64 : Int32
    csinfo("Precision used: $(cfg["precision"])", cfg["suppress_messages"] in TRUELIST)
    is_parallel = cfg["parallelize"] in TRUELIST
    if is_parallel
        n = parse(Int, cfg["max_parallel"])
        csinfo("Starting up Circuitscape to use $n processes in parallel", cfg["suppress_messages"] in TRUELIST)
        myaddprocs(n)
    end
    t = @elapsed r = _compute(T, V, cfg)

    csinfo("Time taken to complete job = $t", cfg["suppress_messages"] in TRUELIST)
    is_parallel && rmprocs(workers())
    r
end

function _compute(T, V, cfg)
    is_raster = cfg["data_type"] in RASTER
    scenario = cfg["scenario"]
    if is_raster
        if scenario in PAIRWISE
            raster_pairwise(T, V, cfg)
        elseif scenario in ADVANCED
            raster_advanced(T, V, cfg)
        elseif scenario in ONETOALL
            raster_one_to_all(T, V, cfg)
        else
            raster_one_to_all(T, V, cfg)
        end
    else
        if scenario in PAIRWISE
            network_pairwise(T, V, cfg)
        else
            network_advanced(T, V, cfg)
        end
    end
end

function compute(dict)
    cfg = init_config()
    update!(cfg, dict)
    update_logging!(cfg)
    write_config(cfg)
    T = cfg["precision"] in SINGLE ? Float32 : Float64
    V = cfg["use_64bit_indexing"] in TRUELIST ? Int64 : Int32
    csinfo("Precision used: $(cfg["precision"])", cfg["suppress_messages"] in TRUELIST)
    is_parallel = cfg["parallelize"] in TRUELIST
    if is_parallel
        n = parse(Int, cfg["max_parallel"])
        csinfo("Starting up Circuitscape to use $n processes in parallel", cfg["suppress_messages"] in TRUELIST)
        myaddprocs(n)
    end
    t = @elapsed r = _compute(T, V, cfg)

    csinfo("Time taken to complete job = $t", cfg["suppress_messages"] in TRUELIST)
    is_parallel && rmprocs(workers())

    r
end


text:
module Repos

include("structure.jl")
include("internal.jl")
include("construction.jl")
include("class.jl")
include("id.jl")
include("interfaces.jl")
include("print.jl")

export Repo,RepoClass,setclass!,getclass,class,print_repo,repo_id

end


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v01.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")

include("structure_factor.jl")
include("gen_dr.jl")
include("gen_rho.jl")
include("calc_ewald.jl")

function test_main()
    Ns = [64,64,64]
    LatVecs = gen_lattice_sc(16.0)
    pw_grid = PWGrid(Ns,LatVecs)
    # Atomic positions and nuclear charge
    Xpos = reshape( [0.0, 0.0, 0.0,
                     1.5, 0.0, 0.0], (3,2) )
    println("Xpos =")
    println(Xpos)

    Z = 1.0
    # Calculate structure factor
    Sf = structure_factor( Xpos, pw_grid.G )

    print("sum(Sf)="); println(sum(Sf))

    E_nn = calc_ewald( pw_grid, Xpos, Sf )
    @printf("E_nn = %18.10f\n", E_nn)
end

@time test_main()


text:
"""
Discrete Laplacian in d dimensions and n equidistant discretization points in each direction
"""
function Lap(n::Integer,d::Integer) #returns the tensor of the discrete Laplacian in a box [0,1]^d with n equidistant discretization points in each direction
    A = zeros(n^d,n^d)
    for j in 0:n^d-1
        J = digits(j,base=n,pad=d) #corresponding index as d-tuples
        for k in 0:d-1
            if J[k+1]==0
                A[j+1,j+1] = 2*d
                A[j+1,j+n^k+1] = -1
            elseif J[k+1]==n-1
                A[j+1,j+1] = 2*d
                A[j+1,j-n^k+1] = -1
            else
                A[j+1,j+1] = 2*d
                A[j+1,j-n^k+1] = -1
                A[j+1,j+n^k+1] = -1
            end
        end
    end
    return A
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
###################################################################################################
# Helpers function for Tensorflow.jl to make building models easier
#

import TensorFlow: tanh
import TensorFlow.nn: sigmoid, zero_state
import TensorFlow.nn.rnn_cell: get_input_dim, LSTMStateTuple

########### General Helpers ###################################

"""
create a summary object to be passed to the summary writer to logg value under name tag
"""
function logg_scalar(value::Real, tag::String)
    return tf.tensorflow.Summary(value=[tf.tensorflow.Summary_Value(simple_value = value, tag=tag)])
end

"""
Reset the graph and open a new session
"""
function init_session(g = Graph())
    tf.set_def_graph(g)
    sess = Session(g)
    return sess
end

"""
returns a list of trainable variables whose name contains name
"""
function get_train_vars_by_name(name::String)
    return [var for var in get_def_graph().collections[:TrainableVariables]
            if contains(tf.get_name(var.var_node), name)]
end

function Base.ndims(A::AbstractTensor)
    length(get_shape(A).dims)
end

"""
Compute the Huber Loss
"""
function huber_loss(x, δ::Float64=1.0)
    mask = abs(x) .< δ
    return mask.*0.5.*x.^2 + (1-mask).*δ.*(abs(x) - 0.5*δ)
end

############ MODEL BUILDING ###################################

"""
    compute_flat_dim(dims::Vector{Nullable{Int64}})
Multiply all the known dimensions
"""
function compute_flat_dim(dims::Vector{Nullable{Int64}})
    res = 1
    for i in dims
        if isnull(i)
            continue
        end
        res *= get(i)
    end
    return res
end

"""
flatten an array or a tensor and keep the batch size
    flatten(x::Array{Float64}; batch_dim::Int64=1)
    flatten(x::Tensor; batch_dim::Int64=1)
"""
function flatten(x::Array{Float64}; batch_dim::Int64=1)
    x_dims = size(x)
    N = ndims(x)
    batch_size = x_dims[batch_dim]
    if batch_dim != 1 #permute dims such that it is in the first axis
        perm = [mod1(batch_dim + i, N) for i=0:N-1]
        x = permutedims(x, perm)
    end
    flat_dim = div(length(x), batch_size)
    @assert batch_size == size(x)[1]
    x = reshape(x, (batch_size, flat_dim))
    return x
end

function flatten(x::Tensor; batch_dim::Int64=1)
    shape = get_shape(x)
    N = length(shape.dims)
    flat_dim = compute_flat_dim(shape.dims)
    if batch_dim != 1 #permute dims such that it is in the first axis
        perm = [mod1(batch_dim + i, N) for i=0:N-1]
        x = permutedims(x, perm)
    end
    x = reshape(x, (-1, flat_dim))
    return x
end

"""
    Build a fully connected layer followed by an activation function
    returns the output of the activation layer as a tensor
    assumes that the input is of shape (batch size, n)
    uses Xavier initialization
    dense(input::Tensor, hidden_units; activation=identity, scope="")
"""
function dense(input::Tensor, hidden_units; activation=identity, scope="fc", reuse=false)
    input_shape = get(get_shape(input).dims[end])
    weight_shape = (input_shape, hidden_units)
    var = 2/(input_shape + hidden_units) # Xavier initialization
    fc_W = variable_scope(scope, initializer= Normal(0.0, var),  reuse=reuse) do
        get_variable("weight", [weight_shape...], Float32)
    end
    fc_b = variable_scope(scope, initializer=tf.ConstantInitializer(0.),  reuse=reuse) do
        get_variable("bias", [hidden_units], Float32)
    end
    a = variable_scope(scope,  reuse=reuse) do
        if activation == identity #XXX hack because identity returns a tensor of shape unknown
            a = input*fc_W + fc_b
        else
            a = activation(input*fc_W + fc_b)
        end
        a
    end
    return a
end

"""
    Build a 2d convolutional layer
    returns the output of the activation as a tensor
    Assumes that the input is of shape (batch, height, width, channels)
    uses Xavier initialization
    conv2d(x, W)

    This version also adds an activation layer
    conv2d(inputs::Tensor, num_filters::Int64, kernel_size::Vector{Int64}, activation=identity; stride::Int64=1, padding::String="SAME")

"""
function conv2d(inputs::Tensor, num_filters::Int64, kernel_size::Vector{Int64};
                activation=identity, stride::Int64=1, padding::String="SAME", scope="conv", reuse::Bool=false)
    # assume inputs is of shape batch, height, width, channels
    input_channels = get(get_shape(inputs).dims[end])
    weight_shape = (kernel_size[1], kernel_size[2], input_channels, num_filters)
    var = 1/(kernel_size[1]*kernel_size[2]*input_channels) # Xavier initialization
    conv_W =  variable_scope(scope, initializer=Normal(0.0, var), reuse=reuse) do
        get_variable("weight", [weight_shape...], Float32)
    end
    conv_b = variable_scope(scope, initializer=tf.ConstantInitializer(0.),  reuse=reuse) do
        get_variable("bias", [num_filters], Float32)
    end
    a = variable_scope(scope,  reuse=reuse) do
        z = nn.conv2d(inputs, conv_W, Int64[1, stride, stride, 1], padding) + conv_b
        a = activation(z)
        a = cast(a, Float32)
    end
    return a
end

function conv2d(x, W)
    nn.conv2d(x, W, [1, 1, 1, 1], "SAME")
end

"""
    Build a max pooling layer, divides the height and width of the input by 2
"""
function max_pool_2x2(x)
    nn.max_pool(x, [1, 2, 2, 1], [1, 2, 2, 1], "SAME")
end

"""
    Build a multi-layer perceptron (mlp) model given the list of nodes in each layer and the input tensor
"""
function mlp(input::Tensor,
            hiddens::Vector{Int64},
            num_output::Int64;
            final_activation=identity,
            scope="mlp",
            reuse = false
            )
    a = flatten(input, batch_dim=1)
    for (i,h) in enumerate(hiddens)
        a = dense(a, h, activation=nn.relu, scope=scope*"/fc_$i", reuse=reuse)
    end
    a = dense(a, num_output, activation=final_activation, scope=scope*"/fc_out", reuse=reuse)
    return a
end

"""
    Builds a model starting with convolutional layers followed by a multi-layer perceptron
- convs: [(int, (int,int), int)]
        list of convolutional layers in form of
        (num_outputs, kernel_size, stride)
- hiddens: [int]
        list of sizes of hidden layers

"""

function cnn_to_mlp(inputs, convs, hiddens, num_output;
    final_activation=identity, scope="conv2mlp", reuse=false, dueling=false)
    out = inputs
    for (i,(nfilters, kernel_size, stride)) in enumerate(convs)
        out = conv2d(out, nfilters, kernel_size,
                     activation=nn.relu, stride=stride, scope=scope*"/conv_$i", reuse=reuse)
    end
    out = variable_scope(scope,  reuse=reuse) do
        flatten(out, batch_dim=1)
    end
    if dueling
        state_val_scope = scope*"/state_value"
        state_out = out
        for (i,h) in enumerate(hiddens)
            state_out = dense(state_out, h, activation=nn.relu, scope=state_val_scope*"/fc_$i", reuse=reuse)
        end
        state_out = dense(state_out, 1, activation=final_activation, scope=state_val_scope*"/fc_out", reuse=reuse)

        action_val_scope = scope*"/action_value"
        action_out = out
        for (i,h) in enumerate(hiddens)
            action_out = dense(action_out, h, activation=nn.relu, scope=action_val_scope*"/fc_$i", reuse=reuse)
        end
        action_out = dense(action_out, num_output, activation=final_activation, scope=action_val_scope*"/fc_out", reuse=reuse)

        actions_mean = Ops.expand_dims(mean(action_out, 2),2) # shape bs x 1

        actions_scaled = action_out - actions_mean # broadcast bs x n_actions - bs x 1
        out = state_out + actions_scaled
    else
        for (i,h) in enumerate(hiddens)
            out = dense(out, h, activation=nn.relu, scope=scope*"/fc_$i", reuse=reuse)
        end
        if num_output > 0
            out = dense(out, num_output, activation=final_activation, scope=scope*"/fc_out", reuse=reuse)
        end
    end
    return out
end

# Modify the methods for LSTM and dynamic_rnn to reuse weights.
# this might break if there are changes in tensorflow.jl , maybe it should be a PR

function (cell::nn.rnn_cell.LSTMCell)(input, state, input_dim=-1; reuse=false, scope="")
    N = get_input_dim(input, input_dim) + cell.hidden_size
    T = eltype(state)
    input = Tensor(input)
    X = [input state.h]
    var = 2/N
    local Wi, Wf, Wo, Wg
    tf.variable_scope(scope*"/Weights", initializer=Normal(0.0, var), reuse=reuse) do
        Wi = get_variable("Wi", [N, cell.hidden_size], T)
        Wf = get_variable("Wf", [N, cell.hidden_size], T)
        Wo = get_variable("Wo", [N, cell.hidden_size], T)
        Wg = get_variable("Wg", [N, cell.hidden_size], T)
    end

    local Bi, Bf, Bo, Bg
    tf.variable_scope(scope*"/Bias", initializer=tf.ConstantInitializer(0.0), reuse=reuse) do
        Bi = get_variable("Bi", [cell.hidden_size], T)
        Bf = get_variable("Bf", [cell.hidden_size], T)
        Bo = get_variable("Bo", [cell.hidden_size], T)
        Bg = get_variable("Bg", [cell.hidden_size], T)
    end

    # TODO make this all one multiply
    I = sigmoid(X*Wi + Bi)
    F = sigmoid(X*Wf + Bf)
    O = sigmoid(X*Wo + Bo)
    G = tanh(X*Wg + Bg)
    C = state.c.*F + G.*I
    S = tanh(C).*O
    return (S, LSTMStateTuple(C, S))
end

# modify the function from tensorflow.jl to support parameter reuse and scoping
function tf.nn.rnn(cell, inputs::Vector, sequence_length=nothing; initial_state=nothing, dtype=nothing, scope="", reuse=false)
    if initial_state === nothing
        initial_state = zero_state(cell, first(inputs), dtype)
    end
    outputs = tf.Tensor[]
    state = initial_state
    rnn_scope = scope*"/RNN"
    for (time_step, input) in enumerate(inputs)
        loop_reuse = time_step>1 || reuse
        tf.variable_scope(scope; reuse=loop_reuse) do
            new_output, new_state = cell(input, state, reuse=loop_reuse, scope=rnn_scope)
        end
        if sequence_length!==nothing && time_step > 1 # This should be removed by the julia lowering process
            # Only update output and state for rows that are not yet passed their ends
            have_passed_end = sequence_length .< time_step
            new_output = tf.select(have_passed_end, outputs[end], new_output)
            new_state = tf.select(have_passed_end, state, new_state)
        end
        state = new_state
        push!(outputs, new_output)
    end
    return outputs, state
end

function tf.nn.rnn(cell, inputs::tf.Tensor, sequence_length=nothing; time_major=false, kwargs...)
    input_list = tf.unstack(inputs; axis = (time_major ? 1 : 2))
    tf.nn.rnn(cell, input_list, sequence_length; kwargs...)
end

# functions from tensorflow.jl tutorial, no control over the initialization
function weight_variable(shape)
    initial = map(Float32, rand(Normal(0, .001), shape...))
    return Variable(initial)
end

function bias_variable(shape)
    initial = zeros(shape...)
    return Variable(initial)
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
module GitHub

import Main, ..Git, ..Dir

const AUTH_DATA = {
    "scopes" => ["repo"],
    "note" => "Julia Package Manager",
    "note_url" => "http://docs.julialang.org/en/latest/manual/packages/",
}

function user()
    if !success(`git config --global github.user`)
        error("""
        no GitHub user name configured; please configure it with:

            git config --global github.user USERNAME

        where USERNAME is replaced with your GitHub user name.
        """)
    end
    readchomp(`git config --global github.user`)
end

function json()
    isdefined(:JSON) || try require("JSON")
    catch err
        warn(err)
        error("using the GitHub API requires having the JSON package installed ")
    end
    Main.JSON
end

function curl(url::String, opts::Cmd=``)
    success(`curl --version`) || error("using the GitHub API requires having `curl` installed")
    out, proc = open(`curl -i -s -S $opts $url`,"r")
    head = readline(out)
    status = int(split(head,r"\s+",3)[2])
    for line in eachline(out)
        ismatch(r"^\s*$",line) || continue
        wait(proc); return status, readall(out)
    end
    error("strangely formatted HTTP response")
end
curl(url::String, data::Nothing, opts::Cmd=``) = curl(url,opts)
curl(url::String, data, opts::Cmd=``) =
    curl(url,`--data $(sprint(io->json().print(io,data))) $opts`)

function token(user::String=user())
    tokfile = Dir.path(".github","token")
    isfile(tokfile) && return strip(readchomp(tokfile))
    status, content = curl("https://api.github.com/authorizations",AUTH_DATA,`-u $user`)
    (status != 401 && status != 403) || error("$status: $(json().parse(content)["message"])")
    tok = json().parse(content)["token"]
    mkpath(dirname(tokfile))
    open(io->println(io,tok),tokfile,"w")
    return tok
end

function req(resource::String, data, opts::Cmd=``)
    url = "https://api.github.com/$resource"
    status, content = curl(url,data,`-u $(token()):x-oauth-basic $opts`)
    response = json().parse(content)
    status, response
end

GET(resource::String, data, opts::Cmd=``) = req(resource,data,opts)
HEAD(resource::String, data, opts::Cmd=``) = req(resource,data,`-I $opts`)
PUT(resource::String, data, opts::Cmd=``) = req(resource,data,`-X PUT $opts`)
POST(resource::String, data, opts::Cmd=``) = req(resource,data,`-X POST $opts`)
PATCH(resource::String, data, opts::Cmd=``) = req(resource,data,`-X PATCH $opts`)
DELETE(resource::String, data, opts::Cmd=``) = req(resource,data,`-X DELETE $opts`)

for m in (:GET,:HEAD,:PUT,:POST,:PATCH,:DELETE)
    @eval $m(resource::String, opts::Cmd=``) = $m(resource,nothing,opts)
end

function pushable(owner::String, repo::String, user::String=user())
    status, response = HEAD("repos/$owner/$repo")
    status == 404 && error("repo $owner/$repo does not exist")
    status, response = GET("repos/$owner/$repo/collaborators/$user")
    status == 204 && return true
    status == 404 && return false
    error("unexpected API status code: $status – $(response["message"])")
end

function fork(owner::String, repo::String)
    status, response = POST("repos/$owner/$repo/forks")
    status == 202 || error("forking $owner/$repo failed: $(response["message"])")
    return response
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function mysum_checked(x)
    s = zero(eltype(x))
    @avx check_empty = true for i ∈ eachindex(x)
        s += x[i]
    end
    s
end
function mysum_unchecked(x)
    s = zero(eltype(x))
    @avx for i ∈ eachindex(x)
        s += x[i]
    end
    s
end

@testset "Check Empty" begin
	@show @__LINE__
    x = fill(9999, 100, 10, 10);
    xv = view(x, :, 1:0, :);
    @test mysum_checked(x) == mysum_unchecked(x) == sum(x)
    @test iszero(mysum_checked(xv))
    @test !iszero(mysum_unchecked(xv))

end


text:
using Test
using MPI
using ClimateMachine.Mesh.BrickMesh
using Random

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    Nelemtotal = 113
    Random.seed!(1234)
    globalcode = randperm(Nelemtotal)

    @assert csize > 1

    bs = [
        (i == 1) ? (1:0) :
            BrickMesh.linearpartition(Nelemtotal, i - 1, csize - 1)
        for i in 1:csize
    ]
    as = [BrickMesh.linearpartition(Nelemtotal, i, csize) for i in 1:csize]

    codeb = globalcode[bs[crank + 1]]

    (so, ss, rs) = BrickMesh.getpartition(comm, codeb)

    codeb = codeb[so]
    codec = []
    for r in 0:(csize - 1)
        sendrange = ss[r + 1]:(ss[r + 2] - 1)
        rcounts = MPI.Allgather(Cint(length(sendrange)), comm)
        c = MPI.Gatherv(view(codeb, sendrange), rcounts, r, comm)
        if r == crank
            codec = c
        end
    end

    codea = (1:Nelemtotal)[as[crank + 1]]

    @test sort(codec) == codea
end

main()


text:
using Knet, Test, Base.Iterators, IterTools, Random # , LinearAlgebra, StatsBase

include("common/debug.jl")
include("common/vocab.jl")
include("common/textreader.jl")
#include("common/embed.jl")
#include("common/linear.jl")
include("common/mask.jl")
include("common/mtdata.jl")

function trmor_word_tokenizer(line)
    (a,b) = split(line)
    return split(a, "")
end

function trmor_tag_tokenizer(line)
    (a,b) = split(line)
    c = split(b, ['+','^'])
    w = popfirst!(c)
    w = split(w, "")
    return append!(w, c)
end

trmor_file = "trmor2018.filtered"

trmor_word_vocab = Vocab(trmor_file; tokenizer=trmor_word_tokenizer)
trmor_tag_vocab = Vocab(trmor_file; tokenizer=trmor_tag_tokenizer)

trmor_word_reader = TextReader(trmor_file, trmor_word_vocab)
trmor_tag_reader = TextReader(trmor_file, trmor_tag_vocab)

# buckets based on src length, so make words src for clean uniform word length batches
trmor_data = MTData(trmor_word_reader, trmor_tag_reader; bucketwidth=1)

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
using REPL
using Scratch

const plotly_local_file_path = Ref{Union{Nothing, String}}(nothing)

function _plots_defaults()
    if isdefined(Main, :PLOTS_DEFAULTS)
        copy(Dict{Symbol,Any}(Main.PLOTS_DEFAULTS))
    else
        Dict{Symbol,Any}()
    end
end

function __init__()
    user_defaults = _plots_defaults()
    if haskey(user_defaults, :theme)
        theme(pop!(user_defaults, :theme); user_defaults...)
    else
        default(; user_defaults...)
    end

    insert!(Base.Multimedia.displays, findlast(x -> x isa Base.TextDisplay || x isa REPL.REPLDisplay, Base.Multimedia.displays) + 1, PlotsDisplay())

    atreplinit(i -> begin
        while PlotsDisplay() in Base.Multimedia.displays
            popdisplay(PlotsDisplay())
        end
        insert!(Base.Multimedia.displays, findlast(x -> x isa REPL.REPLDisplay, Base.Multimedia.displays) + 1, PlotsDisplay())
    end)

    @require HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" begin
        fn = joinpath(@__DIR__, "backends", "hdf5.jl")
        include(fn)
    end

    @require InspectDR = "d0351b0e-4b05-5898-87b3-e2a8edfddd1d" begin
        fn = joinpath(@__DIR__, "backends", "inspectdr.jl")
        include(fn)
    end

    @require PGFPlots = "3b7a836e-365b-5785-a47d-02c71176b4aa" begin
        fn = joinpath(@__DIR__, "backends", "deprecated", "pgfplots.jl")
        include(fn)
    end

    @require PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" begin
        fn = joinpath(@__DIR__, "backends", "plotlybase.jl")
        include(fn)
    end

    @require PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925" begin
        fn = joinpath(@__DIR__, "backends", "pgfplotsx.jl")
        include(fn)
    end

    @require PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" begin
        fn = joinpath(@__DIR__, "backends", "plotlyjs.jl")
        include(fn)
    end

    @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" begin
        fn = joinpath(@__DIR__, "backends", "pyplot.jl")
        include(fn)
    end

    @require UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" begin
        fn = joinpath(@__DIR__, "backends", "unicodeplots.jl")
        include(fn)
    end

    @require IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" begin
        if IJulia.inited
            _init_ijulia_plotting()

            IJulia.display_dict(plt::Plot) = _ijulia_display_dict(plt)
        end
    end

    if get(ENV, "PLOTS_HOST_DEPENDENCY_LOCAL", "false") == "true"
        global plotly_local_file_path[] = joinpath(@get_scratch!("plotly"), "plotly-1.54.2.min.js")
        if !isfile(plotly_local_file_path[])
            download("https://cdn.plot.ly/plotly-1.54.2.min.js", plotly_local_file_path[])
        end

        use_local_plotlyjs[] = true
    end

    use_local_dependencies[] = use_local_plotlyjs[]

    @require FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" begin
        _show(io::IO, mime::MIME"image/png", plt::Plot{<:PDFBackends}) = _show_pdfbackends(io, mime, plt)
    end
end


text:
using Hydra, Test
using Hydra: shuffle

@test lane() == lanewidth() == 1

@test @spmd(4, lane()) == 1

@test @spmd(4, lanewidth()) == 4

@test @spmd(4, lanesum(5)) == 20

@test @spmd(4, shuffle(lane(), 4+1-lane())) == 4

@test @spmd(4, shuffle(2lane(), 3)) == 6

function if_statement(x)
  a = x + 1
  if a > 10
    a += 1000
  elseif a > 5
    a += 100
  else
    a += 10
  end

  if x > 0
    a += 10
  else
    a -= 10
  end
  return a
end

function while_loop(x)
  a = x + 1
  i = 0
  while x > i
    if x >= 10
      a += 3
    elseif x >= 5
      a += 2
    else
      a += 1
    end
    i += 1
  end
  a
end

function for_loop(x)
  a = x + 1
  for i=1:x
    if x > 10
      a += x
    else
      a -= x
    end
  end
  a
end

W = 1:60000 |> collect |> x->reshape(x,(200,300))
b = 1:200 |> collect

function rnn(sequence, ht)
  i = 1
  while i <= length(sequence)
    word = sequence[i]
    ht = W * word + b + ht
    i += 1
  end
  ht
end

mask = Hydra.vect(true,true,true,true)

input = Vector{Int32}([1,6,11,-1])
@test Hydra.spmd(mask, if_statement, Hydra.vect(input...)) == Hydra.vect(map(if_statement, input)...)

input = Vector{Int32}([5,5,5,5])
@test Hydra.spmd(mask, while_loop, Hydra.vect(input...)) == Hydra.vect(map(while_loop, input)...)

input = Vector{Int32}([1,5,12,-3])
@test Hydra.spmd(mask, while_loop, Hydra.vect(input...)) == Hydra.vect(map(while_loop, input)...)

input = Vector{Int32}([2,32,-1,21])
@test Hydra.spmd(mask, for_loop, Hydra.vect(input...)) == Hydra.vect(map(for_loop, input)...)

mask = Hydra.vect(true,true)

input = (
  [collect(1:300),collect(301:600),collect(601:900)],
  [collect(901:1200), collect(1201:1500)])
@test Hydra.spmd(mask, rnn, Hydra.vect(input...), zeros(Int64, 200)).batch == Hydra.vect(map(x->rnn(x,zeros(Int64, 200)), input)...).batch


text:
using Documenter
using MLStyle

makedocs(
    modules = [
        MLStyle.MatchCore,
        MLStyle.DataType,
        MLStyle.Err,
        MLStyle.Pervasives,
        MLStyle.Record,

        # standard patterns
        MLStyle.Active,
        MLStyle.LambdaCases,
        MLStyle.WhenCases,
    ],
    clean = false,
    sitename = "MLStyle.jl",
    linkcheck = false,
    pages = [
        "Home" => "index.md",
        "Syntax" => Any[
            "syntax/adt.md",
            "syntax/switch.md",
            "syntax/pattern.md",
            "syntax/pattern-function.md",
            "syntax/when.md",
        ],
        "Tutorials" => Any["tutorials/capture.md", "tutorials/query-lang.md"],
        "Modules" => Any["modules/ast.md", "modules/cond.md"],
    ],
)

deploydocs(
    repo = "github.com/thautwarm/MLStyle.jl",
    target = "build",
    deps = nothing,
    make = nothing,
    julia = "1.4",
)


text:
# Interaction kernels defined here (they are unit types!)

@inline function interaction_kernel_fourier(ω::Real,popstate::PopulationState)
  return interaction_kernel_fourier(ω,popstate.unittype)
end

# a cluncky sharp step !
struct KernelStep <: UnitType
  τ::Float64
end
@inline function interaction_kernel(t::R,ker::KernelStep) where R<:Real
  return (t < zero(R)) || (t > ker.τ) ? zero(R) :  inv(ker.τ)
end
interaction_kernel_upper(t::Real,pop::KernelStep) = interaction_kernel(t,pop) + eps(100.0)

@inline function interaction_kernel_fourier(ω::R,pop::KernelStep) where R<:Real
  if iszero(ω)
    return one(R)
  end
  cc = 2π*pop.τ*ω
  return (sin(cc)-2*im*sin(0.5cc)^2)/cc
end

# Negative exponential
struct KernelExp <: UnitType
  τ::Float64
end
@inline function interaction_kernel(t::R,ker::KernelExp) where R<:Real
  return t < zero(R) ? zero(R) : exp(-t/ker.τ) / ker.τ
end
interaction_kernel_upper(t::Real,ker::KernelExp) = interaction_kernel(t,ker) + eps(100.0)
interaction_kernel_fourier(ω::Real,ker::KernelExp) =  inv( 1 + im*2*π*ω*ker.τ)

# Alpha-shape with delay

struct KernelAlphaDelay <: UnitType
  τ::Float64
  τdelay::Float64
end

@inline function interaction_kernel(t::Real,ker::KernelAlphaDelay)
  Δt = t - ker.τdelay
  if Δt <= 0.0
    return 0.0
  else
    return Δt/ker.τ^2 * exp(-Δt/ker.τ)
  end
end
@inline function interaction_kernel_upper(t::Real,ker::KernelAlphaDelay)
  Δt = t - ker.τdelay
  if Δt <= ker.τ
    return  inv(ker.τ*ℯ) + eps(100.)
  else
    return  interaction_kernel(t,ker) + eps(100.0)
  end
end
@inline function interaction_kernel_fourier(ω::Real,ker::KernelAlphaDelay)
  return exp(-im*2*π*ω*ker.τdelay) / (1+im*2*π*ω*ker.τ)^2
end


text:
using DashBootstrapComponents, DashHtmlComponents

button_groups = html_div([
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size = "lg",
        className = "me-1",
    ),
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size = "md",
        className = "me-1",
    ),
    dbc_buttongroup(
        [dbc_button("Left"), dbc_button("Middle"), dbc_button("Right")],
        size = "sm",
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
## General

function treeadd_edges!(str,ops,node,M,labels)
    labels[node]=str

    nd=node+1
    for op in ops
        Main.Graphs.add_edge!(M,node,nd)
        add_edges!(op,nd,M,labels)
        nd+=treecount(op)
    end
end

function introspect(A::Union{Operator,Space})
    @eval import TikzGraphs
    m=treecount(A)

    M=Main.Graphs.simple_graph(m)
    labels=Array(ASCIIString,m)

    add_edges!(A,1,M,labels)

    Main.TikzGraphs.plot(M,labels)
end

## Spaces

texname(FS)=string(typeof(FS))

texname(::Line)="\$(-\\infty,\\infty)\$"
texname(d::Interval)="\$["*string(d.a)*","*string(d.b)*"]\$"

texname(CS::CosSpace)="Cos"
texname(CS::SinSpace)="Sin"
texname(CS::Chebyshev)="\$T\$"  #*string(first(domain(CS)))*","*string(last(domain(CS)))*"]"
texname(CS::Ultraspherical)="\$U\^{"*string(order(CS))*"}\$"
texname(CS::Taylor)="\$H\^+\$"
texname(CS::Hardy{false})="\$H\^-\$"
texname(::SumSpace)="\$\\oplus\$"
texname(::PiecewiseSpace)="\$\\bigcup\$"
texname{T}(A::ArraySpace{T,1})="["*string(length(A))*"]"
texname{T}(A::ArraySpace{T,2})="["*string(size(A,1))*"\$\\times\$"*string(size(A,2))*"]"
texname(J::JacobiWeight)="\$(1+x)\^{"*string(J.α)*"}(1-x)\^{"*string(J.β)*"}\$"

treecount(S::Union{DirectSumSpace,PiecewiseSpace})=1+mapreduce(treecount,+,S.spaces)
treecount(S::Union{ArraySpace,JacobiWeight})=1+treecount(S.space)
treecount(::Space)=1

add_edges!(FS::Space,nd,M,labels)=(labels[nd]=string(nd)*":"*texname(FS))

for (OP) in (:DirectSumSpace,:PiecewiseSpace)
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),A.spaces,nd,M,labels)
end

for (OP) in (:ArraySpace,:JacobiWeight)
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),[A.space],nd,M,labels)
end

## Operators

treecount(::Operator) = 1

treecount(::Derivative) = 1
treecount(::ConstantOperator) = 1
treecount(::Conversion) = 1
treecount(::ConcreteMultiplication) = 1

treecount(M::Union{MultiplicationWrapper,DerivativeWrapper,IntegralWrapper,LeftIntegralWrapper,RightIntegralWrapper,
                   ConversionWrapper,SpaceOperator,DiagonalArrayOperator,ConstantTimesOperator})=1+treecount(M.op)
treecount(A::Union{PlusOperator,TimesOperator,InterlaceOperator,
                   DiagonalInterlaceOperator,
                   KroneckerOperator})=1+mapreduce(treecount,+,A.ops)

domainrangestr(A)=texname(domainspace(A))*"\$\\rightarrow\$"*texname(rangespace(A))

texname(A::Operator)=string(typeof(A))
texname(D::ConcreteDerivative)=(D.order==1?"\$D":"\$D\^"*string(D.order))*"\$:"*domainrangestr(D)
texname(A::ConstantOperator)=string(A.c)*"I"
texname(C::ConcreteConversion)="C:"*domainrangestr(C)
texname(A::ConcreteMultiplication)=space(A.f)==domainspace(A)==rangespace(A)?"M["*texname(space(A.f))*"]":"M["*texname(space(A.f))*"]:"domainrangestr(A)

texname(D::ConcreteIntegral)=(D.order==1?"\$Q":"\$Q\^"*string(D.order))*"\$:"*domainrangestr(D)
texname(D::ConcreteLeftIntegral)=(D.order==1?"\$Q_{$(first(domain(D)))}":"\$Q_{$(first(domain(D)))}\^"*string(D.order))*"\$:"*domainrangestr(D)

texname(D::DerivativeWrapper)=(D.order==1?"\$(D":"\$(D\^"*string(D.order))*")\$"
texname(A::SpaceOperator)="("*domainrangestr(A)*")"
texname(A::ConstantTimesOperator)=string(A.c)

add_edges!(A::Operator,nd,M,labels)=(labels[nd]=string(nd)*":"*texname(A))

for (WRAP,STR) in ((:MultiplicationWrapper,:"(M)"),(:ConversionWrapper,"(C)"),(:DiagonalArrayOperator,:"DiagArray"))
    @eval add_edges!(A::$WRAP,nd,M,labels)=treeadd_edges!(string(nd)*":"*$STR,[A.op],nd,M,labels)
end

@eval add_edges!(A::Union{SpaceOperator,DerivativeWrapper},nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),[A.op],nd,M,labels)

for (OP,STR) in ((:PlusOperator,:"+"),(:TimesOperator,:"*"),(:InterlaceOperator,:"Interlace"),
    (:DiagonalInterlaceOperator,:"\$\\oplus\$"),
                 (:KroneckerOperator,:"\$\\otimes\$"))
    @eval add_edges!(A::$OP,nd,M,labels)=treeadd_edges!(string(nd)*":"*$STR,A.ops,nd,M,labels)
end

add_edges!(A::ConstantTimesOperator,nd,M,labels)=treeadd_edges!(string(nd)*":"*string(A.c),[A.op],nd,M,labels)


text:
using DataDepsGenerators
using Test

using ReferenceTests

@testset "DataCite test" begin
    @test_reference "references/DataCite Ecology.txt" remove_cite_version(generate(DataCite(), "10.5061/dryad.74699"))

    @test_reference "references/DataCite Forage Fish.txt" remove_cite_version(generate(DataCite(), "10.5063/F1HT2M7Q"))

    @test_reference "references/DataCite Fire Patch.txt" remove_cite_version(generate(DataCite(), "https://search.datacite.org/works/10.15148/0e999ffc-e220-41ac-ac85-76e92ecd0320"))

    @test_reference "references/DataCite Ceramic.txt" remove_cite_version(generate(DataCite(), "10.5281/zenodo.1147572"))
end


text:
using Base.Test
using MIPVerify
using MIPVerify: UnrestrictedPerturbationFamily, BlurringPerturbationFamily
using MIPVerify: LInfNormBoundedPerturbationFamily
isdefined(:TestHelpers) || include("../../../TestHelpers.jl")
using TestHelpers: batch_test_adversarial_example

@testset "Conv + Conv + FC + Softmax" begin

### Parameters for neural net
batch = 1
in1_height = 8
in1_width = 8
stride1_height = 2
stride1_width = 2
strides1 = (1, stride1_height, stride1_width, 1)
pooled1_height = round(Int, in1_height/stride1_height, RoundUp)
pooled1_width = round(Int, in1_width/stride1_width, RoundUp)
in1_channels = 1
filter1_height = 2
filter1_width = 2
out1_channels = 2

in2_height = pooled1_height
in2_width = pooled1_width
stride2_height = 1
stride2_width = 1
strides2 = (1, stride2_height, stride2_width, 1)
pooled2_height = round(Int, in2_height/stride2_height, RoundUp)
pooled2_width = round(Int, in2_width/stride2_width, RoundUp)
in2_channels = out1_channels
filter2_height = 2
filter2_width = 2
out2_channels = 2

A_height = 5
A_width = pooled2_height*pooled2_width*out2_channels

B_height = 3
B_width = A_height

### Choosing data to be used
srand(5)
input_size = (batch, in1_height, in1_width, in1_channels)
x0 = rand(input_size)

kernelc1 = rand(filter1_height, filter1_width, in1_channels, out1_channels)*2-1
biasc1 = rand(out1_channels)*2-1

kernelc2 = rand(filter2_height, filter2_width, in2_channels, out2_channels)*2-1
biasc2 = rand(out2_channels)*2-1

kernelf1 = rand(A_width, A_height)*2-1
biasf1 = rand(A_height)*2-1

kernelf2 = rand(B_width, B_height)*2-1
biasf2 = rand(B_height)*2-1

nn = Sequential(
    [
        Conv2d(kernelc1, biasc1), MaxPool(strides1), ReLU(),
        Conv2d(kernelc2, biasc2), MaxPool(strides2), ReLU(),
        Flatten(4),
        Linear(kernelf1, biasf1), ReLU(),
        Linear(kernelf2, biasf2)
    ],
    "tests.integration.generated_weights.conv+conv+fc+softmax"
)

pp_blur = BlurringPerturbationFamily((5, 5))
pp_unrestricted = UnrestrictedPerturbationFamily()

expected_objective_values = Dict(
    (2, pp_unrestricted, 1, 0.1) => 0.0,
    (2, pp_unrestricted, 1, 1) => 0.991616,
    (2, pp_unrestricted, 1, 2) => NaN,
    (2, pp_unrestricted, Inf, 0.1) => 0.0,
    (2, pp_unrestricted, Inf, 1) => 0.0953527,
    (2, LInfNormBoundedPerturbationFamily(0.09), Inf, 1) => NaN,
    (2, LInfNormBoundedPerturbationFamily(0.0954), Inf, 1) => 0.0953527,
    (2, LInfNormBoundedPerturbationFamily(0.096), Inf, 1) => 0.0953527,
    (2, pp_blur, 1, 0.6) => 0.0,
    (2, pp_blur, 1, 0.625) => 1.40955,
    (2, pp_blur, 1, 0.65) => NaN,
    (2, pp_blur, Inf, 0.6) => 0.0,
    (2, pp_blur, Inf, 0.625) => 0.0570452
)

batch_test_adversarial_example(nn, x0, expected_objective_values)

end


text:
#to plot uncertainities in neural networks
function plot_predictive(data, trajectories, xs; μ=0, σ=0, title=["Plot"], legend = false)

	lt = length(trajectories)
	if lt < 1
		throw("Err: No data")
	elseif lt == 1
		μ = mean(trajectories["1"], dims=2)
		σ = std(trajectories["1"], dims=2)
		(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
		f_axes.scatter(data[:,1],data[:,2], c="red", marker=".")
		f_axes.plot(xs,vec(μ), c="blue")
		f_axes.fill_between(xs, vec(μ+3σ), vec(μ-3σ), alpha=0.5)

		f_axes.set_title(title[1])
		fig.show()
	else
		nrows = Int(ceil(lt/2))
		fig, f_axes = PyPlot.subplots(ncols=2, nrows=nrows)
		for i in 1:lt
			μ = mean(trajectories["$i"], dims=2)
			σ = std(trajectories["$i"], dims=2)
			f_axes[i].scatter(data[:,1],data[:,2], c="red", marker=".")
			f_axes[i].plot(xs,vec(μ), c="blue")
			f_axes[i].fill_between(xs, vec(μ+3σ), vec(μ-3σ), alpha=0.5)
			f_axes[i].set_title(title[i])
		end
		fig.show()
	end
end
#To plot uncertainites of Neural ODE for each variables
function plot_node(t, trajectories, data, data_w_noise, n_vars, datasize, title; mean_fit = nothing)
	lt = length(trajectories)

	start_loc = 1
	end_loc = datasize

	for nv in 1:n_vars
		mx_data_n = maximum(data_w_noise[start_loc:end_loc,:], dims=2)
		mn_data_n = minimum(data_w_noise[start_loc:end_loc,:], dims=2)
		if lt < 1
			throw("no data")
		elseif lt == 1
			mx = maximum(trajectories[1][start_loc:end_loc,:], dims=2)
			mn = minimum(trajectories[1][start_loc:end_loc,:], dims=2)
			(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
			f_axes.fill_between(t, vec(mx_data_n), vec(mn_data_n), alpha=0.3, label ="data", color="red")
			if !(mean_fit isa Nothing)
				f_axes.plot(t,vec(mean_fit[start_loc:end_loc]), c="green", label = "Mean fit")
			end
			f_axes.plot(t,vec(data[start_loc:end_loc,1]), c="red", label = "data without noise")
			f_axes.plot(t,vec(mx), c="blue") #maximum values
			f_axes.plot(t,vec(mn), c="blue") #minimum values
			f_axes.fill_between(t, vec(mx), vec(mn), alpha=0.3, label ="uncertainty", color="blue")
			f_axes.set_title(title[1])
			f_axes.set_xlabel("time")
			f_axes.set_ylabel("solution")
			f_axes.legend()
			fig.show()
			fig.suptitle("Variable :$nv")
		else
			nrows = Int(ceil(lt/2))
			fig, f_axes = PyPlot.subplots(ncols=2, nrows=nrows)
			for ti in 1:lt
				mx = maximum(trajectories[ti][start_loc:end_loc,:], dims=2)
				mn = minimum(trajectories[ti][start_loc:end_loc,:], dims=2)
				f_axes[ti].fill_between(t, vec(mx_data_n), vec(mn_data_n), alpha=0.3, label ="Data", color="red")
				if !(mean_fit isa Nothing)
					f_axes[ti].plot(t,vec(mean_fit[start_loc:end_loc]), c="green", label = "Mean fit")
				end
				f_axes[ti].plot(t,vec(data[start_loc:end_loc,1]), c="red", label = "Data without noise")
				f_axes[ti].plot(t,vec(mx), c="blue") #maximum values
				f_axes[ti].plot(t,vec(mn), c="blue") #minimum values
				f_axes[ti].fill_between(t, vec(mx), vec(mn), alpha=0.3, label ="Uncertainty", color="blue")
				f_axes[ti].set_title(title[ti])
				f_axes[ti].set_xlabel("time")
				f_axes[ti].set_ylabel("solution")
				f_axes[ti].legend(loc=1)
			end
			fig.suptitle("Variable :$nv")
			fig.show()
		end

		start_loc = end_loc+1
		end_loc = start_loc+datasize-1
	end
end
function get_cmap(n, name="hsv")
	return PyPlot.cm.get_cmap(name, n)
end
function plot_single_node(trajectories, data_train, data_forcat, data_train_w_noise,
	data_train_w_noise_f, n_vars, datasize_train, datasize_forcast, t_data, t_forecast)
	(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
	cmap =  get_cmap(10)
	start_loc = 1
	end_loc = datasize_train
	start_locf = 1
	end_locf = datasize_forcast
	locs =  findall(x->x>=t_data[end], t_forecast)
	prepend!(locs,[locs[1]-1])
	new_locs = locs

	un_start_loc = 1
	un_end_loc = datasize_train

	color_set = ["red", "green", "purple", "orange", "blue", "cyan","yellow","pink"]

	for nv in 1:n_vars

		#training data without noise
		f_axes.plot(t_data,vec(data_train[start_loc:end_loc,1]), c=color_set[nv], label = "var $(nv)")
		f_axes.plot(t_forecast[locs],vec(data_forcat[new_locs,1]), c=color_set[nv+2], label = "var $(nv) forecast")

		mx_data_n = maximum(data_train_w_noise[start_loc:end_loc,:], dims=2)
		mn_data_n = minimum(data_train_w_noise[start_loc:end_loc,:], dims=2)
		f_axes.fill_between(t_data, vec(mx_data_n), vec(mn_data_n), alpha=0.3, label ="var $(nv) with noise", color=color_set[nv])

		mx_data_f = maximum(data_train_w_noise_f[new_locs,:], dims=2)
		mn_data_f = minimum(data_train_w_noise_f[new_locs,:], dims=2)
		f_axes.fill_between(t_forecast[locs], vec(mx_data_f), vec(mn_data_f), alpha=0.3, label ="Forecasted var $(nv) with", color=color_set[nv+2])

		mx = maximum(trajectories[un_start_loc:un_end_loc,:], dims=2)
		mn = minimum(trajectories[un_start_loc:un_end_loc,:], dims=2)

		f_axes.fill_between(t_data, vec(mx), vec(mn), alpha=0.3, label ="uncertainty $(nv)", color=cmap(nv+4))

		mx_f = maximum(trajectories[new_locs,:], dims=2)
		mn_f = minimum(trajectories[new_locs,:], dims=2)

		f_axes.fill_between(t_forecast[locs], vec(mx_f), vec(mn_f), alpha=0.3, label ="uncertainty $(nv) forecast", color=cmap(nv+6))

		start_loc = end_loc+1
		end_loc = start_loc+datasize_train-1

		un_start_loc = new_locs[end]+1
		un_end_loc = un_start_loc+datasize_train-1
		new_locs = new_locs .+ datasize_forcast
	end
	f_axes.legend()
	fig.show()
end
function variableplot_single_node(trajectories, data_train, data_forcast, data_train_w_noise,
	data_train_w_noise_f, n_vars, datasize_train, datasize_forcast, t_data, t_forecast)

	locs =  findall(x->x>=t_data[end], t_forecast)
	prepend!(locs,[locs[1]-1])
	new_locs = locs
	tm, tn = size(trajectories)
	un_start_loc = new_locs[end]+1
	un_end_loc = un_start_loc+datasize_train-1

	start_loc = 1
	end_loc = datasize_train

	(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
	f_axes.plot(vec(data_train[start_loc:end_loc,1]),vec(data_train[start_loc+datasize_train:end,1]),
		c="red", label="data")

	mx_data_n = maximum(data_train_w_noise, dims=2)
	mn_data_n = minimum(data_train_w_noise, dims=2)

	# f_axes.fill_between(vec(data_train[start_loc:end_loc]), vec(mx_data_n[start_loc+datasize_train:end]), vec(mn_data_n[start_loc+datasize_train:end]),
	# 	alpha=0.3, color="red", label="noisy data")
	for t in 1:tn
		f_axes.plot(vec(trajectories[start_loc:end_loc,t]),vec(trajectories[un_start_loc:un_end_loc,t]),
		alpha = 0.2, color="blue")
	end
	# mx = maximum(trajectories[un_start_loc:un_end_loc,:], dims=2)
	# mn = minimum(trajectories[un_start_loc:un_end_loc,:], dims=2)

	# f_axes.fill_between(vec(data_train[start_loc:end_loc,1]), vec(mx), vec(mn), alpha=0.3, color="blue", label="uncertainty")

	f_axes.legend()
	f_axes.set_xlabel("var 1")
	f_axes.set_ylabel("var 2")
	fig.show()
end
#To plot predictions
function plot_pred(t, n_ode, ode_data_bkp, datasize, u0)
	(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
	pred = n_ode(vec(u0[:,1])) # Get the prediction using the correct initial condition
	f_axes.plot(t,vec(ode_data_bkp[1:datasize,1]), c="red", marker=".", label = "data")
	f_axes.plot(t,vec(pred[1,:]), c="green", marker=".", label ="prediction")
	f_axes.legend()
	fig.show()
end
#To plot NeuralODE based forecasting
function plot_forecast(u0, t, nt, tspan, new_node, ode_data_bkp,
	ode_data_bkp_f, datasize)
	locs =  findall(x->x>=tspan[2], nt)
	(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
	pred = new_node(vec(u0[:,1])) # Get the prediction using the correct initial condition
	f_axes.plot(t,vec(ode_data_bkp[1:datasize,1]), c="red", marker=".", label = "data")
	f_axes.plot(nt[locs],vec(ode_data_bkp_f[locs,1]), c="purple", marker=".", label = "data")
	f_axes.plot(nt,vec(pred[1,:]), c="green", marker=".", label ="prediction")
	f_axes.legend()
	fig.show()
end

function plot_noise_data(t, ode_data, ode_data_bkp, datasize, len)
	(fig, f_axes) = PyPlot.subplots(ncols=1, nrows=1)
	for i in 1:len
		f_axes.scatter(t,vec(ode_data[1:1:datasize,i]), c="green", marker="*", label ="data with noise")
	end
	f_axes.plot(t,vec(ode_data_bkp[1:1:datasize,1]), c="red", marker=".", label = "data")
	f_axes.legend()
	fig.show()
end


text:
julia> a = [4,2,3,1]
4-element Int32 Array:
 4
 2
 3
 1
julia> sort(a) #out-of-place/non-mutating sort
4-element Int32 Array:
 1
 2
 3
 4

julia> a
4-element Int32 Array:
 4
 2
 3
 1

julia> sort!(a) # in-place/mutating sort
4-element Int32 Array:
 1
 2
 3
 4

julia> a
4-element Int32 Array:
 1
 2
 3
 4

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Different versions of variadic function defined in pub_expr.h
#
#   SCIP_RETCODE SCIPexprCreate(BMS_BLKMEM* blkmem, SCIP_EXPR** expr,
#                               SCIP_EXPROP op, ...);
#
# Prototype without extra arguments:
#
#   function SCIPexprCreate(blkmem, expr, op, ...)
#       ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
#             (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, ...),
#             blkmem, expr, op, ...)
#   end
#
#
# TODO: Fix warnings below by creating thin, non-variadic C wrappers for the
#       following calls.

# SCIP_EXPR_VARIDX
function SCIPexprCreate(blkmem, expr, op, varpos::Cint)
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Cint...),
          blkmem, expr, op, varpos)
end

# SCIP_EXPR_CONST
function SCIPexprCreate(blkmem, expr, op, value::Cdouble)
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Cdouble...),
          blkmem, expr, op, value)
end

# SCIP_EXPR_SQRT, SCIP_EXPR_EXP, SCIP_EXPR_LOG (unary op)
function SCIPexprCreate(blkmem, expr, op, child::Ptr{SCIP_EXPR})
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Ptr{SCIP_EXPR}...),
          blkmem, expr, op, child)
end

# SCIP_EXPR_MINUS, SCIP_EXPR_DIV (binary op)
function SCIPexprCreate(blkmem, expr, op, left::Ptr{SCIP_EXPR}, right::Ptr{SCIP_EXPR})
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Ptr{SCIP_EXPR}...),
          blkmem, expr, op, left, right)
end

# SCIP_EXPR_SUM, SCIP_EXPR_PRODUCT (n-ary op)
function SCIPexprCreate(blkmem, expr, op, nchildren::Cint, children::Vector{Ptr{SCIP_EXPR}})
    # WARNING: In the actual C function, the last two arguments given here are
    # part of the variadic arguments. But since Julia only supports variadic
    # arguments of the same type, we specify the first one (nchildren::Cint)
    # directly, and pretend that only the last one is variadic.
    @assert nchildren == length(children)
    GC.@preserve children begin
        ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
              (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Cint, Ptr{Ptr{SCIP_EXPR}}...),
              blkmem, expr, op, nchildren, pointer(children))
    end
end

# SCIP_EXPR_REALPOWER (mixed binary op)
function SCIPexprCreate(blkmem, expr, op, base::Ptr{SCIP_EXPR}, exponent::Cdouble)
    # WARNING: In the actual C function, the last two arguments given here are
    # part of the variadic arguments. But since Julia only supports variadic
    # arguments of the same type, we specify the first one (base::Ptr{SCIP_EXPR})
    # directly, and pretend that only the last one is variadic.
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Ptr{SCIP_EXPR}, Cdouble...),
          blkmem, expr, op, base, exponent)
end

# SCIP_EXPR_INTPOWER (mixed binary op)
function SCIPexprCreate(blkmem, expr, op, base::Ptr{SCIP_EXPR}, exponent::Cint)
    # WARNING: In the actual C function, the last two arguments given here are
    # part of the variadic arguments. But since Julia only supports variadic
    # arguments of the same type, we specify the first one (base::Ptr{SCIP_EXPR})
    # directly, and pretend that only the last one is variadic.
    ccall((:SCIPexprCreate, libscip), SCIP_RETCODE,
          (Ptr{BMS_BLKMEM}, Ptr{Ptr{SCIP_EXPR}}, SCIP_EXPROP, Ptr{SCIP_EXPR}, Cint...),
          blkmem, expr, op, base, exponent)
end


text:
# File that contains macros used throughout the code :)
#macro debugprint(pos,msg)
#	:(j4pr_debug == true && print_with_color(:cyan, "  -->  [ J4PR ][$(uppercase($pos))] $($msg)"))
#end


text:
#################################################################################
# Comments
#
# - Ideally the net_injection variables would be bounded.  This can be done using an adhoc data model extention
#
#################################################################################
# Model Definitions

""
function instantiate_nip_expr_model(data::Dict{String, Any}, model_constructor; kwargs...)
    return PM.instantiate_model(data, model_constructor, instantiate_nip_expr; kwargs...)
end

""
function instantiate_nip_expr(pm::PM.AbstractPowerModel)
    for (n, network) in PM.nws(pm)
        @assert !PM.ismulticonductor(pm, nw = n)
        PM.variable_bus_voltage(pm, nw = n)
        PM.variable_branch_power(pm, nw = n; bounded = false)
        PM.variable_dcline_power(pm, nw = n)

        PM.constraint_model_voltage(pm, nw = n)

        for i in PM.ids(pm, :ref_buses, nw = n)
            PM.constraint_theta_ref(pm, i, nw = n)
        end

        for i in PM.ids(pm, :bus, nw = n)
            constraint_power_balance_ni_expr(pm, i, nw = n)
        end

        for i in PM.ids(pm, :branch, nw = n)
            PM.constraint_ohms_yt_from(pm, i, nw = n)
            PM.constraint_ohms_yt_to(pm, i, nw = n)

            PM.constraint_voltage_angle_difference(pm, i, nw = n)
        end

        for i in PM.ids(pm, :dcline)
            PM.constraint_dcline_power_losses(pm, i, nw = n)
        end
    end

    return

end

function instantiate_bfp_expr_model(data::Dict{String, Any}, model_constructor; kwargs...)
    return PM.instantiate_model(data, model_constructor, instantiate_bfp_expr; kwargs...)
end

""
function instantiate_bfp_expr(pm::PM.AbstractPowerModel)
    for (n, network) in PM.nws(pm)
        @assert !PM.ismulticonductor(pm, nw = n)
        PM.variable_bus_voltage(pm, nw = n)
        PM.variable_branch_power(pm, nw = n; bounded = false)
        PM.variable_dcline_power(pm, nw = n)

        PM.constraint_model_current(pm, nw = n)

        for i in PM.ids(pm, :ref_buses, nw = n)
            PM.constraint_theta_ref(pm, i, nw = n)
        end

        for i in PM.ids(pm, :bus, nw = n)
            constraint_power_balance_ni_expr(pm, i, nw = n)
        end

        for i in PM.ids(pm, :branch, nw = n)
            PM.constraint_power_losses(pm, i, nw = n)
            PM.constraint_voltage_magnitude_difference(pm, i, nw = n)

            PM.constraint_voltage_angle_difference(pm, i, nw = n)
        end

        for i in PM.ids(pm, :dcline)
            PM.constraint_dcline_power_losses(pm, i, nw = n)
        end
    end

    return

end

function instantiate_vip_expr_model(data::Dict{String, Any}, model_constructor; kwargs...)
    throw(error("VI Models not currently supported"))
end

#################################################################################
# Model Extention Functions

""
function constraint_power_balance_ni_expr(
    pm::PM.AbstractPowerModel,
    i::Int;
    nw::Int = pm.cnw,
)
    if !haskey(PM.con(pm, nw), :power_balance_p)
        PM.con(pm, nw)[:power_balance_p] = Dict{Int, JuMP.ConstraintRef}()
    end
    if !haskey(PM.con(pm, nw), :power_balance_q)
        PM.con(pm, nw)[:power_balance_q] = Dict{Int, JuMP.ConstraintRef}()
    end

    bus = PM.ref(pm, nw, :bus, i)
    bus_arcs = PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc = PM.ref(pm, nw, :bus_arcs_dc, i)

    pni_expr = PM.ref(pm, nw, :bus, i, "pni")
    qni_expr = PM.ref(pm, nw, :bus, i, "qni")

    constraint_power_balance_ni_expr(pm, nw, i, bus_arcs, bus_arcs_dc, pni_expr, qni_expr)

    return

end

""
function constraint_power_balance_ni_expr(
    pm::PM.AbstractPowerModel,
    n::Int,
    i::Int,
    bus_arcs,
    bus_arcs_dc,
    pni_expr,
    qni_expr,
)
    p = PM.var(pm, n, :p)
    q = PM.var(pm, n, :q)
    p_dc = PM.var(pm, n, :p_dc)
    q_dc = PM.var(pm, n, :q_dc)

    PM.con(pm, n, :power_balance_p)[i] = JuMP.@constraint(
        pm.model,
        sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == pni_expr
    )
    PM.con(pm, n, :power_balance_q)[i] = JuMP.@constraint(
        pm.model,
        sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == qni_expr
    )

    return

end

""
function constraint_current_balance_ni_expr(
    pm::PM.AbstractPowerModel,
    i::Int;
    nw::Int = pm.cnw,
)
    if !haskey(PM.con(pm, nw), :kcl_cr)
        PM.con(pm, nw)[:kcl_cr] = Dict{Int, JuMP.ConstraintRef}()
    end
    if !haskey(PM.con(pm, nw), :kcl_ci)
        PM.con(pm, nw)[:kcl_ci] = Dict{Int, JuMP.ConstraintRef}()
    end

    bus = PM.ref(pm, nw, :bus, i)
    bus_arcs = PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc = PM.ref(pm, nw, :bus_arcs_dc, i)

    pni_expr = PM.ref(pm, nw, :bus, i, "pni")
    qni_expr = PM.ref(pm, nw, :bus, i, "qni")

    constraint_current_balance_ni_expr(pm, nw, i, bus_arcs, bus_arcs_dc, pni_expr, qni_expr)

    return

end

""
function constraint_current_balance_ni_expr(
    pm::PM.AbstractPowerModel,
    n::Int,
    i::Int,
    bus_arcs,
    bus_arcs_dc,
    pni_expr,
    qni_expr,
)
    p = PM.var(pm, n, :p)
    q = PM.var(pm, n, :q)
    p_dc = PM.var(pm, n, :p_dc)
    q_dc = PM.var(pm, n, :q_dc)

    PM.con(pm, n, :power_balance_p)[i] = JuMP.@constraint(
        pm.model,
        sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == pni_expr
    )
    PM.con(pm, n, :power_balance_q)[i] = JuMP.@constraint(
        pm.model,
        sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == qni_expr
    )

    return

end

"active power only models ignore reactive power variables"
function variable_reactive_net_injection(pm::PM.AbstractActivePowerModel; kwargs...)
    return
end

""
function constraint_power_balance_ni_expr(
    pm::PM.AbstractActivePowerModel,
    n::Int,
    i::Int,
    bus_arcs,
    bus_arcs_dc,
    pni_expr,
    qni_expr,
)
    p = PM.var(pm, n, :p)
    p_dc = PM.var(pm, n, :p_dc)

    PM.con(pm, n, :power_balance_p)[i] = JuMP.@constraint(
        pm.model,
        sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == pni_expr
    )

    return

end

""
function powermodels_network!(
    psi_container::PSIContainer,
    system_formulation::Type{S},
    sys::PSY.System,
    instantiate_model = instantiate_nip_expr_model,
) where {S <: PM.AbstractPowerModel}

    time_steps = model_time_steps(psi_container)
    pm_data, PM_map = pass_to_pm(sys, time_steps[end])
    buses = PSY.get_components(PSY.Bus, sys)

    remove_undef!(psi_container.expressions[:nodal_balance_active])
    remove_undef!(psi_container.expressions[:nodal_balance_reactive])

    for t in time_steps, bus in buses
        pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["pni"] =
            psi_container.expressions[:nodal_balance_active][bus.number, t]
        pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["qni"] =
            psi_container.expressions[:nodal_balance_reactive][bus.number, t]
    end

    psi_container.pm =
        instantiate_model(pm_data, system_formulation, jump_model = psi_container.JuMPmodel)
    psi_container.pm.ext[:PMmap] = PM_map

    return

end

""
function powermodels_network!(
    psi_container::PSIContainer,
    system_formulation::Type{S},
    sys::PSY.System,
    instantiate_model = instantiate_nip_expr_model,
) where {S <: PM.AbstractActivePowerModel}

    time_steps = model_time_steps(psi_container)
    pm_data, PM_map = pass_to_pm(sys, time_steps[end])
    buses = PSY.get_components(PSY.Bus, sys)

    remove_undef!(psi_container.expressions[:nodal_balance_active])

    for t in time_steps, bus in buses
        pm_data["nw"]["$(t)"]["bus"]["$(PSY.get_number(bus))"]["pni"] =
            psi_container.expressions[:nodal_balance_active][PSY.get_number(bus), t]
        #pm_data["nw"]["$(t)"]["bus"]["$(bus.number)"]["qni"] = 0.0
    end

    psi_container.pm =
        instantiate_model(pm_data, system_formulation, jump_model = psi_container.JuMPmodel)
    psi_container.pm.ext[:PMmap] = PM_map

    return

end

#### PM accessor functions ########

function PMvarmap(system_formulation::Type{S}) where {S <: PM.AbstractDCPModel}
    pm_var_map = Dict{Type, Dict{Symbol, Union{String, NamedTuple}}}()

    pm_var_map[PSY.Bus] = Dict(:va => THETA)
    pm_var_map[PSY.ACBranch] = Dict(:p => (from_to = FLOW_ACTIVE_POWER, to_from = nothing))
    pm_var_map[PSY.DCBranch] =
        Dict(:p_dc => (from_to = FLOW_ACTIVE_POWER, to_from = nothing))

    return pm_var_map
end

function PMvarmap(system_formulation::Type{S}) where {S <: PM.AbstractActivePowerModel}
    pm_var_map = Dict{Type, Dict{Symbol, Union{String, NamedTuple}}}()

    pm_var_map[PSY.Bus] = Dict(:va => THETA)
    pm_var_map[PSY.ACBranch] = Dict(
        :p =>
            (from_to = FLOW_ACTIVE_POWER_FROM_TO, to_from = FLOW_ACTIVE_POWER_TO_FROM),
    )
    pm_var_map[PSY.DCBranch] = Dict(
        :p_dc =>
            (from_to = FLOW_ACTIVE_POWER_FROM_TO, to_from = FLOW_ACTIVE_POWER_TO_FROM),
    )

    return pm_var_map
end

function PMvarmap(system_formulation::Type{S}) where {S <: PM.AbstractPowerModel}
    pm_var_map = Dict{Type, Dict{Symbol, Union{String, NamedTuple}}}()

    pm_var_map[PSY.Bus] = Dict(:va => THETA, :vm => VM)
    pm_var_map[PSY.ACBranch] = Dict(
        :p =>
            (from_to = FLOW_ACTIVE_POWER_FROM_TO, to_from = FLOW_ACTIVE_POWER_TO_FROM),
        :q => (
            from_to = FLOW_REACTIVE_POWER_FROM_TO,
            to_from = FLOW_REACTIVE_POWER_TO_FROM,
        ),
    )
    pm_var_map[PSY.DCBranch] = Dict(
        :p_dc =>
            (from_to = FLOW_ACTIVE_POWER_FROM_TO, to_from = FLOW_ACTIVE_POWER_TO_FROM),
        :q_dc => (
            from_to = FLOW_REACTIVE_POWER_FROM_TO,
            to_from = FLOW_REACTIVE_POWER_TO_FROM,
        ),
    )

    return pm_var_map
end

function PMconmap(system_formulation::Type{S}) where {S <: PM.AbstractActivePowerModel}
    pm_con_map = Dict{Type, Dict{Symbol, String}}()

    pm_con_map[PSY.Bus] = Dict(:power_balance_p => NODAL_BALANCE_ACTIVE)
    return pm_con_map
end

function PMconmap(system_formulation::Type{S}) where {S <: PM.AbstractPowerModel}
    pm_con_map = Dict{Type, Dict{Symbol, String}}()

    pm_con_map[PSY.Bus] = Dict(
        :power_balance_p => NODAL_BALANCE_ACTIVE,
        :power_balance_q => NODAL_BALANCE_REACTIVE,
    )
    return pm_con_map
end

function add_pm_var_refs!(
    psi_container::PSIContainer,
    system_formulation::Type{S},
    sys::PSY.System,
) where {S <: PM.AbstractPowerModel}

    time_steps = model_time_steps(psi_container)
    bus_dict = psi_container.pm.ext[:PMmap].bus
    ACbranch_dict = psi_container.pm.ext[:PMmap].arcs
    ACbranch_types = typeof.(values(ACbranch_dict))
    DCbranch_dict = psi_container.pm.ext[:PMmap].arcs_dc
    DCbranch_types = typeof.(values(DCbranch_dict))

    pm_var_names = keys(psi_container.pm.var[:nw][1])

    pm_var_map = PMvarmap(system_formulation)

    for (pm_v, ps_v) in pm_var_map[PSY.Bus]
        if pm_v in pm_var_names
            container = PSI.container_spec(
                psi_container.JuMPmodel,
                [PSY.get_name(b) for b in values(bus_dict)],
                time_steps,
            )
            assign_variable!(psi_container, ps_v, PSY.Bus, container)
            for t in time_steps, (pm_bus, bus) in bus_dict
                name = PSY.get_name(bus)
                container[name, t] = PM.var(psi_container.pm, t, pm_v)[pm_bus] #pm_vars[pm_v][pm_bus]
            end
        end
    end

    add_pm_var_refs!(
        psi_container,
        PSY.ACBranch,
        ACbranch_types,
        ACbranch_dict,
        pm_var_map,
        pm_var_names,
        time_steps,
    )
    add_pm_var_refs!(
        psi_container,
        PSY.DCBranch,
        DCbranch_types,
        DCbranch_dict,
        pm_var_map,
        pm_var_names,
        time_steps,
    )

end

function add_pm_var_refs!(
    psi_container::PSIContainer,
    d_class::Type,
    device_types::Vector,
    pm_map::Dict,
    pm_var_map::Dict,
    pm_var_names::Base.KeySet,
    time_steps::UnitRange{Int},
)

    for d_type in Set(device_types)
        devices = [d for d in pm_map if typeof(d[2]) == d_type]
        for (pm_v, ps_v) in pm_var_map[d_class]
            if pm_v in pm_var_names
                for dir in fieldnames(typeof(ps_v))
                    isnothing(getfield(ps_v, dir)) && continue
                    # TODO: make a better mapping with the var names in the definitions file
                    var_name = getfield(ps_v, dir)
                    container = PSI.container_spec(
                        psi_container.JuMPmodel,
                        [PSY.get_name(d[2]) for d in devices],
                        time_steps,
                    )
                    assign_variable!(psi_container, var_name, d_type, container)
                    for t in time_steps, (pm_d, d) in devices
                        container[PSY.get_name(d), t] =
                            PM.var(psi_container.pm, t, pm_v, getfield(pm_d, dir))
                    end
                end
            end
        end
    end
end

function add_pm_con_refs!(
    psi_container::PSIContainer,
    system_formulation::Type{S},
    sys::PSY.System,
) where {S <: PM.AbstractPowerModel}

    time_steps = model_time_steps(psi_container)
    bus_dict = psi_container.pm.ext[:PMmap].bus

    pm_con_names = [
        k
        for
        k in keys(psi_container.pm.con[:nw][1]) if !isempty(PM.con(psi_container.pm, 1, k))
    ]

    pm_con_map = PMconmap(system_formulation)
    for (pm_v, ps_v) in pm_con_map[PSY.Bus]
        if pm_v in pm_con_names
            container = PSI.add_cons_container!(
                psi_container,
                make_constraint_name(ps_v, PSY.Bus),
                [PSY.get_name(b) for b in values(bus_dict)],
                time_steps,
            )
            for t in time_steps, (pm_bus, bus) in bus_dict
                name = PSY.get_name(bus)
                container[name, t] = PM.con(psi_container.pm, t, pm_v)[pm_bus]
            end
        end
    end
end


text:
n_features = (32, 32, 3)
n_targets = (coarse=1, fine=1)

@testset "trainset" begin
    d = CIFAR100()

    @test d.split == :train
    @test extrema(d.features) == (0, 1)
    @test convert2image(d, 1) isa AbstractMatrix{<:RGB}
    @test convert2image(d, 1:2) isa AbstractArray{<:RGB,3}

    test_supervised_array_dataset(d;
        n_features, n_targets, n_obs=50000,
        Tx=Float32, Ty=Int)
end

@testset "testset" begin
    d = CIFAR100(split=:test, Tx=UInt8)

    @test d.split == :test
    @test extrema(d.features) == (0, 255)
    @test convert2image(d, 1) isa AbstractMatrix{<:RGB}
    @test convert2image(d, 1:2) isa AbstractArray{<:RGB,3}

    test_supervised_array_dataset(d;
        n_features, n_targets, n_obs=10000,
        Tx=UInt8, Ty=Int)
end


text:
# Packages used in this chapter
# To check for installed packages in V1.0 we now need to use the PKG API

using Pkg
ipks = Pkg.installed()
pkgs = ["BenchmarkTools","Match","PyPlot","StaticArrays"]
for p in pkgs
    !haskey(ipks,p) && Pkg.add(p)
end

# We can define these here as as required
using BenchmarkTools, Match, PyPlot, StaticArrays

# We need these as some functions have moved from Base to Stdlib
#
using Printf, LinearAlgebra

function hailstone(n)
   k = 1
   a = [n]
   while n > 1
      n = (n % 2 == 0) ? n >> 1 : 3n + 1
      push!(a,n)
      k += 1
   end
   a
end

hailstone(17)

using PyPlot

f(x, p = 0.0) = exp(p*x)*sin(x);
t = range(0.0,stop=8.0*π,length=101);
w = zeros(length(t));
for i = 1:length(w) w[i] = f(t[i], 0.1) end

plot(t, w)

using Match

allodds(x) = @match x begin
    m::Int, if iseven(m) || (m < 3) end        => "Not a valid choice"
    n::Int, if n == 3 || n == 5 || n == 7 end  => "$x is prime"
    _  =>  "By induction all numbers are prime"
end

allodds(3)

allodds(6)

allodds(9)

struct Queen
    x::Int
    y::Int
end

qhorz(qa, qb) = qa.x == qb.x;
qvert(qa, qb) = qa.y == qb.y;
qdiag(qa, qb) = abs(qa.x - qb.x) == abs(qa.y - qb.y);

qhvd(qa, qb) = qhorz(qa, qb) || qvert(qa, qb) || qdiag(qa, qb);
qany(testq, qs) = any(q -> qhvd(testq, q), qs);

function qsolve(nsqsx, nsqsy, nqs, presqs = ())
    nqs == 0 && return presqs
    for xsq in 1:nsqsx
        for ysq in 1:nsqsy
            testq = Queen(xsq, ysq)
            if !qany(testq, presqs)
                tryqs = (presqs..., testq)
                maybe = qsolve(nsqsx, nsqsy, nqs - 1, tryqs)
                maybe !== nothing && return maybe
            end
        end
    end
    return nothing
end

qsolve(nqs) = qsolve(nqs, nqs, nqs)

qsolve(8)

module Vehicles

export Contact, Vehicle, Car, Bike, Yacht, Powerboat, Boat
export Ford, BMW, VW, Scooter, MotorBike, Speedboat

const KNOTS_TO_MPH = 1.151

struct Contact
    name::String
    email::String
    phone::String
end

abstract type Vehicle end

abstract type Car <: Vehicle end
abstract type Bike <: Vehicle  end
abstract type Boat <: Vehicle end

abstract type Powerboat <: Boat end

struct Ford <: Car
    owner::Contact
    make::String
    fuel::String
    color::String
    engine_cc::Int64
    speed_mph::Float64
    function Ford(owner, make, engine_cc,speed_mph)
        new(owner,make,"Petrol","Black",engine_cc,speed_mph)
     end
end

struct BMW <: Car
    owner::Contact
    make::String
    fuel::String
    color::String
    engine_cc::Int64
    speed_mph::Float64
    function BMW(owner,make,engine_cc,speed_mph)
        new(owner,make,"Petrol","Blue",engine_cc,speed_mph)
   end
end

struct VW <: Car
    owner::Contact
    make::String
    fuel::String
    color::String
    engine_cc::Int64
    speed_mph::Float64
end

struct MotorBike <: Bike
    owner::Contact
    make::String
    engine_cc::Int64
    speed_mph::Float64
end

struct Scooter <: Bike
    owner::Contact
    make::String
    engine_cc::Int64
    speed_mph::Float64
end

mutable struct Yacht <: Boat
    owner::Contact
    make::String
    length_m::Float64
end

mutable struct Speedboat <: Powerboat
    owner::Contact
    make::String
    fuel::String
    engine_cc::Int64
    speed_knots::Float64
    length_m::Float64
end

function is_quicker(a::VW, b::BMW)
   if (a.speed_mph == b.speed_mph)
      return nothing
    else
      return(a.speed_mph > b.speed_mph ? a : b)
     end
end

function is_quicker(a::Speedboat, b::Scooter)
   a_mph = KNOTS_TO_MPH * a.speed_knots
   if (a_mph == b.speed_mph)
       return nothing
    else
       return(a_mph > b.speed_mph ? a : b)
    end
end

function is_longer(a::Yacht, b::Speedboat)
   if (a.length_m == b.length_m)
       return nothing
    else
       return(a.length > b.length_m ? a : b)
    end
end

is_quicker(a::BMW, b::VW) = is_quicker(b,a)
is_quicker(a::Scooter, b::Speedboat) = is_quicker(b,a)
is_longer(a::Speedboat, b::Yacht) = is_longer(b,a)

end

using Main.Vehicles;

malcolm = Contact("Malcolm","malcolm@abc.net","07777555999");
myCar = Ford(malcolm, "Model T", 1000, 50.0);
myBike = Scooter(malcolm, "Vespa", 125, 35.0);

james = Contact("James","james@abc.net","07777666888");
jmCar = BMW(james,"Series 500", 3200, 125.0);
jmBoat = Yacht(james,"Oceanis 44",14.6);
jmBike = MotorBike(james, "Harley", 850, 120.0);

david = Contact("David","dave@abc.net","07777222444");
dvCar = VW(david,"Golf", "diesel", "red", 1800, 85.0);
dvBoat = Speedboat(david,"Sealine","petrol", 600, 45.0, 8.2);

myCar.owner = david

dvBoat.owner = malcolm

s = [myCar, jmCar, dvCar]
for c in s
    who = c.owner.name
    model = c.make
    println("$who has a $model car")
end

# Note that Malcolm now owns Dave's boat
#
s =[dvCar, dvBoat]
for c in s
    who = c.owner.name
    model = c.make
    println("$who has a $model vehicle")
end

println("James owns these vehicles: ")
s = [jmCar, jmBike, jmBoat]
for c in s
    println("....\t$(c.make)")
end

# Make this mutable as I might change where I live
#
mutable struct Address
  name::String
  street::String
  city::String
  country::String
  postcode::String
end

postal = Address("Malcolm Sherrington","1 Main Street", "London", "UK", "WC2N 9ZZ");

const Owner = Union{Contact, Address}

struct Skiff <: Boat
   owner::Owner
   make::String
   length_m::Float64
end

y1 = Skiff(malcolm,"Moody 36", 11.02)

y2 = Skiff(postal,"Dufour 44", 13.47)

@assert y1.length_m < y2.length_m

# This module uses Float64 components but could use a parameterised type {T}
#
module V3D

# import Base.+, Base.*, Base./, Base.==, Base.<, Base.>
#
import Base: +, *, /, ==, <, >, zero, one
import LinearAlgebra: norm, dot

export Vec3, norm, dist, zero, one

struct Vec3
    x::Float64
    y::Float64
    z::Float64
end

(+)(a::Vec3, b::Vec3) = Vec3(a.x+b.x, a.y+b.y, a.z+b.z)
(*)(p::Vec3, s::Real) = Vec3(p.x*s, p.y*s, p.z*s)
(*)(s::Real, p::Vec3) = p*s
(/)(p::Vec3, s::Real) = (1.0/s)*p

(==)(a::Vec3, b::Vec3) = (a.x==b.x)&&(a.y==b.y)&&(a.z==b.z) ? true : false;

dot(a::Vec3, b::Vec3) = a.x*b.x + a.y*b.y + a.z*b.z;
norm(a::Vec3) = sqrt(dot(a,a));

zero(a::Vec3) = Vec3(0.0,0.0,0.0);
one(a::Vec3)  = Vec3(1.0,1.0,1.0);

(<)(a::Vec3, b::Vec3) = norm(a) < norm(b) ? true : false;
(>)(a::Vec3, b::Vec3) = norm(a) > norm(b) ? true : false;

dist(a::Vec3, b::Vec3) = sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z))

end

using Main.V3D
using Printf

v1 = Vec3(1.2,3.4,5.6);
v2 = Vec3(2.1,4.3,6.5);
@printf "Distance between vectors is %.3f \n" dist(v1,v2)

vv = [Vec3(rand(),rand(),rand()) for i = 1:1000000];
vs = reshape(vv,1000,1000)

using LinearAlgebra

dv = map(v -> norm(v), vs)
eigvals(dv)

eigvecs(dv)

k = 0
for i in 1:length(vv)
  if norm(vv[i]) < 1.0
    global k +=1
  end
end
@printf "Estimate of PI is %9.5f\n" 6.0*k/length(vv)

#
# Defining the Vec3 parameterically and the rest of the module is unchanged
# i.e.
#
# immutable Vec3P{T}
#    x::T
#    y::T
#    z::T
# end
#
module V3P
#
# import Base.+, Base.*, Base./, Base.norm, Base.==, Base.<, Base.>
#
import Base: +, *, /, ==, <, >
import LinearAlgebra: norm, dot
export Vec3P, norm, distp

struct Vec3P{T}
    x::T
    y::T
    z::T
end

(+)(a::Vec3P, b::Vec3P) = Vec3P(a.x+b.x, a.y+b.y, a.z+b.z)
(*)(p::Vec3P, s::Real) = Vec3P(p.x*s, p.y*s, p.z*s)
(*)(s::Real, p::Vec3P) = p*s
(/)(p::Vec3P, s::Real) = (1.0/s)*p

(==)(a::Vec3P, b::Vec3P) = (a.x==b.x)&&(a.y==b.y)&&(a.z==b.z) ? true : false;

dot(a::Vec3P, b::Vec3P) = a.x*b.x + a.y*b.y + a.z*b.z;
norm(a::Vec3P) = sqrt(dot(a,a));

(<)(a::Vec3P, b::Vec3P) = norm(a) < norm(b) ? true : false;
(>)(a::Vec3P, b::Vec3P) = norm(a) > norm(b) ? true : false;

distp(a::Vec3P, b::Vec3P) = sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z))

end

using Main.V3P

z1 = Vec3P{Complex}(1+2im,2+3im,3+4im);
z2 = Vec3P{Complex}(3-2im,4-3im,5-4im);
z1 + z2

zn = norm(z1 + z2)        # => sqrt(116) but IS a complex number.
@assert zn == sqrt(116)

# Notice that the vaule is returned as a complex.
zn

r1 = Vec3P{Rational}(11//7,13//5,8//17)
r2 = Vec3P{Rational}(17//9,23//15,28//17)
r1 + r2

# Norm is promoted to a Real as it can not be represented as a Rational
#
norm(r1 + r2)

module VNX
using StaticArrays

import Base: +, *, /, ==, <, >
import LinearAlgebra: norm, dot
export VecN, norm, distx

struct VecN
    sv::SVector;
end

sizeof(a::VecN) = length(a.sv)
sOK(a::VecN, b::VecN) =
 (sizeof(a) == sizeof(b)) ? true : throw(BoundsError("Vector of different lengths"));

(+)(a::VecN, b::VecN) = [a.sv[i] + b.sv[i] for i in 1:sizeof(a) if sOK(a,b)]
(*)(x::Real, a::VecN) = [a.sv[i]*x for i in 1:sizeof(a)]
(*)(a::VecN, x::Real) = x*a
(/)(a::VecN, x::Real) = [a.sv[i]/x for i in 1:sizeof(a)]

(==)(a::VecN, b::VecN) = any([(a.sv[i] == b.sv[i]) for i in 1:sizeof(a) if sOK(a,b)])

dot(a::VecN, b::VecN) = sum([a.sv[i]*b.sv[i] for i in 1:sizeof(a) if sOK(a,b)])
norm(a::VecN) = sqrt(dot(a,a));

(<)(a::VecN, b::VecN) = norm(a) < norm(b) ? true : false;
(>)(a::VecN, b::VecN) = norm(a) > norm(b) ? true : false;

distx(a::VecN, b::VecN) = sum(map(x -> x*x,[a.sv[i]-b.sv[i] for i in 1:sizeof(a) if sOK(a,b)]))

end

# Generate K 4-vectors
#
using Main.VNX, StaticArrays

K = 10^5;
vv = Array{VecN}(undef,K);

for j = 1:K
  vv[j] =  VecN(@SVector [rand() for i = 1:4])
end

# Sum up the vectors which lie within the 4-ball
#
s = 0;
for j = 1:K
  if (norm(vv[j]) < 1.0) s += 1 end
end

# Volume of the unit 4-ball is 2*π*π
# The count is 1/16th of the volume
#
mypi = sqrt(32*s/K)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_InstanceViewStatus.jl")
include("model_RunCommandDocument.jl")
include("model_RunCommandDocumentBase.jl")
include("model_RunCommandInput.jl")
include("model_RunCommandInputParameter.jl")
include("model_RunCommandListResult.jl")
include("model_RunCommandParameterDefinition.jl")
include("model_RunCommandResult.jl")


text:
__precompile__()
module MDtables

using StringBuilders, FixedSizeStrings
import Base:print

export SimpleTable, MultilineTable, GridTable, print

abstract type MDtable end

struct SimpleTable{T<:StringBuilder} <: MDtable
    o::T
end

struct MultilineTable{T<:StringBuilder} <: MDtable
    o::T
end

struct GridTable{T<:StringBuilder} <: MDtable
    o::T
end

struct Cell{H, W}
    ss::NTuple{H, FixedSizeString{W}}
end
Cell(a::Vector{<:AbstractString}) = Cell((FixedSizeString.(a)...))
function Cell(x::String)
    a = split(x, '\n')
    h = maximum(length, a)
    a .= rpad.(a, h)
    Cell(a)
end
Cell(x) = Cell(string(x))
geth(x::Cell{H,W}) where {H, W} = H
getw(x::Cell{H,W}) where {H, W} = W
pad(x::Cell{H, W}, h, w) where {H, W} = Cell([rpad.(x.ss, w)..., (" "^w for i in H+1:h)...])

function printline(o, cells, l, delimiter)
    for cell in cells
        append!(o, delimiter)
        append!(o, cell.ss[l])
    end
    append!(o, delimiter)
    append!(o, "\n")
end
printline(o::SimpleTable, cells, l) = printline(o.o, cells, l, " ")
printline(o::MultilineTable, cells, l) = printline(o.o, cells, l, " ")
printline(o::GridTable, cells, l) = printline(o.o, cells, l, "|")
printseparator(o, ws, c, d) = printline(o.o, [Cell((FixedSizeString(c^w), )) for w in ws], 1, d)
printtopbottom(o::SimpleTable, ws) = nothing
printtopbottom(o::MultilineTable, ws) = printseparator(o, ws, "-", "-")
printtopbottom(o::GridTable, ws) = printseparator(o, ws, "-", "+")
printdelimiter(o::SimpleTable, ws) = nothing
printdelimiter(o::MultilineTable, ws) = printseparator(o, ws, " ", " ")
printdelimiter(o::GridTable, ws) = printtopbottom(o, ws)
printunderline(o::SimpleTable, ws) = printseparator(o, ws, "-", " ")
printunderline(o::MultilineTable, ws) = printseparator(o, ws, "-", " ")
printunderline(o::GridTable, ws) = printseparator(o, ws, "=", "+")
function printrow(o, cells)
    for l in 1:geth(cells[1])
        printline(o, cells, l)
    end
end
function getcells(a)
    b = Cell.(a)
    nrow, ncol = size(b)
    hs = [maximum(geth, b[r, :]) for r in 1:nrow]
    ws = [maximum(getw, b[:, c]) for c in 1:ncol]
    for r in 1:nrow, c in 1:ncol
        b[r, c] = pad(b[r, c], hs[r], ws[c])
    end
    return (b, ws, nrow)
end

function print(::Type{T}, a::A) where {T <: MDtable, A <: AbstractArray}
    b, ws, nrow = getcells(a)
    if T <: SimpleTable
        @assert all(geth(bi) == 1 for bi in b) "can't print multiline table as a simple_table"
    end
    o = T(StringBuilder())
    printtopbottom(o, ws)
    printrow(o, b[1, :])
    printunderline(o, ws)
    for r in 2:nrow-1
        printrow(o, b[r, :])
        printdelimiter(o, ws)
    end
    printrow(o, b[nrow, :])
    printtopbottom(o, ws)
    return String(o.o)
end

end


text:
# This function helps to construct `getBlock`, while keeping the data `location_in_blocks` static (local, yet compute-once and persistent).
function getBlock_constructor()
    # Function not in Julia v0.3, but will be added in Julia v0.4.
    cld(x::Int, y::Int) = fld((x - 1), y) + 1

    location_in_blocks = Array(Int, 9, 9)
    for row in 1:9, col in 1:9
        location_in_blocks[row, col] = 3 * (cld(row, 3) - 1) + cld(col, 3)
    end

    return (cell::Cell) -> location_in_blocks[cell.row, cell.col]
end
getBlock = getBlock_constructor()
getRow(cell::Cell) = cell.row
getCol(cell::Cell) = cell.col

Base.getindex(s::GameState, cell::Cell) = s.data[cell.row, cell.col]
Base.setindex!(s::GameState, number::Number, cell::Cell) = (s.data[cell.row, cell.col] = number)

Base.getindex(s::CellsState, cell::Cell) = s.data[cell.row, cell.col]
Base.setindex!(s::CellsState, number::Number, cell::Cell) = (s.data[cell.row, cell.col] = number)

Base.getindex(s::ComponentState, cell::Cell) = s.data[s.findComponentIndex(cell)]
Base.setindex!(s::ComponentState, set::LeftoverNumbers, cell::Cell) = (s.data[s.findComponentIndex(cell)] = set)

Base.length(s::SequenceState) = Base.length(s.data)
Base.getindex(s::SequenceState, i::Int) = s.data[i]
Base.getindex(s::SequenceState, r::UnitRange) = [s.data[i] for i in r]
Base.setindex!(s::SequenceState, cell::Cell, i::Int) = (s.data[i] = cell)

Base.getindex(lc::LevelCount, i::Int) = lc.data[i]
Base.setindex!(lc::LevelCount, value::Int, i::Int) = (lc.data[i] = value)

writeSolution!(solutions::Vector{Matrix{Number}}, s::CellsState) = Base.push!(solutions, deepcopy(s.data))
getTotalCount(lc::LevelCount) = sum(lc.data)

# Set cell to be the given number.
# E.g. if `cell == 1` and `number == 2`,
# then we will set entry of cell (1, 1) with the value of 2.
# Notice that we will also remove 2 from corresponding component so it will not be used for subsequent allocations.
function setCell!(gs::GameState, cell::Cell, number::Number)
    writeNumberToCell!(gs, cell, number)
    removeNumberFromBlock!(gs, cell, number)
    removeNumberFromRow!(gs, cell, number)
    removeNumberFromCol!(gs, cell, number)
end

# Clear cell does the opposite of set cell.
# Clear cell will erase the number allocated to it,
# while reinstating its possibility to be use by other empty cells.
function clearCell!(gs::GameState, cell::Cell)
    reinstateNumberToBlock!(gs, cell)
    reinstateNumberToRow!(gs, cell)
    reinstateNumberToCol!(gs, cell)
    eraseNumberFromCell!(gs, cell)
end

# Write number to cell.
function writeNumberToCell!(gs::GameState, cell::Cell, number::Number)
    gs.cells[cell] = number
end
# Erase number from cell.
function eraseNumberFromCell!(gs::GameState, cell::Cell)
    gs.cells[cell] = BLANK
end

# Remove the number from corresponding (block/row/column) because it has been allocated to a cell within them.
function removeNumberFromBlock!(gs::GameState, cell::Cell, number::Number)
    gs.blocks[cell] = removeNumber(gs.blocks[cell], number)
end
function removeNumberFromRow!(gs::GameState, cell::Cell, number::Number)
    gs.rows[cell] = removeNumber(gs.rows[cell], number)
end
function removeNumberFromCol!(gs::GameState, cell::Cell, number::Number)
    gs.cols[cell] = removeNumber(gs.cols[cell], number)
end

# Reinstate the number to corresponding (block/row/column) because it is freed from a cell within them.
function reinstateNumberToBlock!(gs::GameState, cell::Cell)
    gs.blocks[cell] = includeNumber(gs.blocks[cell], gs.cells[cell])
end
function reinstateNumberToRow!(gs::GameState, cell::Cell)
    gs.rows[cell] = includeNumber(gs.rows[cell], gs.cells[cell])
end
function reinstateNumberToCol!(gs::GameState, cell::Cell)
    gs.cols[cell] = includeNumber(gs.cols[cell], gs.cells[cell])
end

# Get the leftover numbers that can be filled into indicated cell.
getLeftoverNumbers(gs::GameState, cell::Cell) = gs.blocks[cell] ∩ gs.rows[cell] ∩ gs.cols[cell]

(∩)(set₁::LeftoverNumbers, set₂::LeftoverNumbers) = LeftoverNumbers(set₁.val_bits & set₂.val_bits)
function removeNumber(set::LeftoverNumbers, i::Int)
    LeftoverNumbers(set.val_bits & (~(1 << i)))
end
function includeNumber(set::LeftoverNumbers, i::Int)
    LeftoverNumbers(set.val_bits | (1 << i))
end
getSize(set::LeftoverNumbers) = count_ones(set.val_bits)

Base.start(iter::LeftoverNumbers) = iter.val_bits
function Base.next(iter::LeftoverNumbers, state)
    i = state & (-state) #Lowest 1 bit in state.
    state &= ~i
    return (trailing_zeros(i), state)
end
Base.done(iter::LeftoverNumbers, state) = (state == zero(Uint16))

function swapSequenceEntries!(gs::GameState, S1::Int, S2::Int)
    (gs.sequence[S1], gs.sequence[S2]) = (gs.sequence[S2], gs.sequence[S1])
end

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
# This example requires that also the packages Distributions, LabelledArrays,
# CSV and Statistics are installed; install by
# using Pkg; Pkg.add("CSV"); Pkg.add("Distributions"); Pkg.add("LabelledArrays"); Pkg.add("Statistics"); Pkg.add("DataFrames")
using AdaptiveParticleMCMC, LabelledArrays, Distributions, CSV, Statistics, DataFrames

# Define the particle type for the model (here, latent is univariate AR(1))
mutable struct SVParticle
    s::Float64
    SVParticle() = new(0.0) # Void constructor required!
end

# Model parameters:
mutable struct SVParam
    ρ::Float64     # Latent AR(1) coefficient
    σ::Float64     # Latent AR(1) noise sd
    β::Float64     # Latent mean
    σ_s1::Float64  # Latent AR(1) stationary sd
end

# Monthly S&P 500 data (from https://datahub.io/core/s-and-p-500):
data = CSV.read(joinpath(@__DIR__,"sp500post2000.csv"), DataFrame)
sp500_data = diff(log.(data.SP500)) # Monthly log-returns
sp500_data .-= mean(sp500_data)     # Remove trend
# Initialise parameters
sv_par = SVParam(0.9,1.0,0.0,1.0)

# This will be the SequentialMonteCarlo "particle scratch", which
# will contain both model data & parameters, and which will be
# the 'scratch' argument of M_ar1!, lM_ar1, lG_sv
struct SVScratch
    par::SVParam
    y::Vector{Float64}
    SVScratch() = new(sv_par, sp500_data)
end

# Transition *simulator* of stationary zero-mean AR(1) with
# parameters (ρ, σ)
function M_ar1!(x, rng, k, x_prev, scratch)
    if k == 1
        x.s = rand(rng, Normal(0.0, scratch.par.σ_s1))
    else
        x.s = rand(rng, Normal(scratch.par.ρ*x_prev.s, scratch.par.σ))
    end
end

# Log transition *density* of stationary zero-mean AR(1) with parameters (ρ, σ)
# (only used with Particle Gibbs, not with Particle marginal Metropolis-Hastings)
function lM_ar1(k, x_prev, x, scratch)
    if k == 1
        return logpdf(Normal(0.0, scratch.par.σ_s1), x.s)
    else
        return logpdf(Normal(scratch.par.ρ*x_prev.s, scratch.par.σ), x.s)
    end
end

# Potential (stochastic volatility observation log-density)
function lG_sv(k, x, scratch)
    s = exp(.5*(scratch.par.β + x.s))        # Observation sd
    logpdf(Normal(0.0, s), scratch.y[k])   # y[k] ~ N(0, s^2)
end

# We are sampling transformed parameters (logit(ρ), log(σ), β)
# This function sets the model parameters based on sampled (transformed)
# parameters to the 'actual' model parameter values.
inv_logit(x) = 2.0/(1.0+exp(-x)) - 1.0 # (-∞,∞) → (-1,1)
function set_param!(scratch, θ)
    scratch.par.ρ = inv_logit(θ.logit_̢rho)
    scratch.par.σ = exp(θ.log_sigma)
    scratch.par.β = θ.beta
    scratch.par.σ_s1 = scratch.par.σ/sqrt(1.0 - scratch.par.ρ^2) # Stationary variance
end

# Normal prior for the transformed parameters
function prior(theta)
    (logpdf(Normal(6.0,5.0), theta.logit_̢rho)
    +logpdf(Normal(-1.0,5.0), theta.log_sigma)
    +logpdf(Normal(-6.0,5.0), theta.beta))
end

# Create data structures for SequentialMonteCarlo:
##################################################
T = length(sp500_data)
N = 64     # Number of particles
n = 40_000 # Number of PMCMC iterations

# Set up the SMC state
state = SMCState(T, N, SVParticle, SVScratch, set_param!, lG_sv, M_ar1!, lM_ar1)

# Initial (transformed) parameter vector
theta0 = LVector(logit_̢rho=0.0, log_sigma=0.0, beta=0.0)
# Particle marginal Metropolis-Hastings with Adaptive Metropolis
out_pmmh = adaptive_pmmh(theta0, prior, state, n;
  thin=100, show_progress=2, save_paths=true);

# Particle Gibbs with Robust Adaptive Metropolis
out_pg = adaptive_pg(theta0, prior, state, n;
  thin=100, show_progress=2, save_paths=true);

using StatsPlots, Statistics
function quantile_plot!(plt, S, p=0.95; x=1:size(S)[1])
    alpha = (1-p)/2
    qs = [alpha, 0.5, 1-alpha]
    Q = mapslices(x->quantile(x, qs), S, dims=2)
    plot!(plt, x, Q[:,2], ribbon=(Q[:,2]-Q[:,1], Q[:,3]-Q[:,2]), color=:black,
    legend=:false, fillalpha=0.2)
end
function show_out(out; title="")
    labels = [string.(typeof(out.theta0).parameters[4])...]
    p_theta = corrplot(out.Theta', size=(600,600), label=labels,
    title="Parameter posterior$title")
    # The volatilities:
    S = [out.X[j][i].s+out.Theta[3,j] for i=1:length(out.X[1]), j=1:length(out.X)]
    p_paths = plot(xlabel="Time", size=(600,800), ylabel="Log-volatility",
    legend=false, title="Latent posterior median, 50% and 95% credible intervals")
    quantile_plot!(p_paths, S, 0.95; x=data.Date[2:end])
    quantile_plot!(p_paths, S, 0.5; x=data.Date[2:end])
    p_data = plot(data.Date, log.(data.SP500), size=(600,800), ylabel="Log SP500", legend=false)
    plot(p_theta, p_paths, p_data, layout=grid(3,1, heights=[0.7,0.15,0.15]))
end
show_out(out_pg)