text:
using DataFrames, JSON3, JSONTables, ZipFile, LazyJSON

function get_events(events::String)
    zarchive = ZipFile.Reader("/tmp/events.zip")
    dictio = Dict(zarchive.files[i].name => i for i in eachindex(zarchive.files))
    file_num = dictio[events]
    str = read(zarchive.files[file_num])

    return str
end

events_test = "events_Italy.json"
str = get_events(events_test)
v = LazyJSON.parse(str);
v = convert(Vector{Any},v)
# @time subset_indexes_lazy(v, 2576336)
"""
    function convert_events()
"""
function convert_events(data::DataFrame)

    rename!(data, :id => :event_id,
                    :matchId => :game_id,
                    :matchPeriod => :period_id,
                    :eventSec => :milliseconds,
                    :teamId => :team_id,
                    :playerId => :player_id,
                    :eventId => :type_id,
                    :eventName => :type_name,
                    :subEventId => :subtype_id,
                    :subEventName => :subtype_name)

    data[:, :milliseconds] = data[:,:milliseconds] * 1000
    return data
end

"""
    function subset_index_lazy

Use the vector Any v which is the JSON parsed lazily and the game_id to request data
with matchId correspoding to the game id.
"""
function subset_indexes_lazy(v::Vector{Any}, game_id::Int)
    json_data = Vector{Any}()
    @inbounds for i in eachindex(v)
        if v[i]["matchId"] == game_id
            push!(json_data, v[i])
        end
    end
    # json_data = Vector(Any, 2000)
    # count = 0
    # @inbounds for i in eachindex(v)
    #     if v[i]["matchId"] == game_id
    #         count += 1
    #         json_data[count] = v[i]
    #     end
    # end
    return json_data
end

"""
    function transform_to_df

"""
function transform_to_df(json_subset)
    subset_dic = convert.(Dict, json_subset)

    test_df = DataFrame(:playerId => zeros(length(subset_dic)),
        :matchId => zeros(length(subset_dic)),
        :eventName => "",
        :positions => [Vector{Dict}[] for i in eachindex(subset_dic)],
        :eventId => zeros(length(subset_dic)),
        :subEventName => "",
        :teamId => zeros(length(subset_dic)),
        :id => zeros(length(subset_dic)),
        :matchPeriod => "",
        :subEventId => "",
        :eventSec => zeros(length(subset_dic)),
        :tags => [Vector{Dict}[] for i in eachindex(subset_dic)])

    for i in eachindex(test_df.playerId)
        test_df[i, :playerId] = subset_dic[i]["playerId"]
        test_df[i, :positions] = [convert.(Dict,subset_dic[i]["positions"])]
        test_df[i, :matchId] = subset_dic[i]["matchId"]
        test_df[i, :eventName] = string(subset_dic[i]["eventName"])
        test_df[i, :eventId] = subset_dic[i]["eventId"]
        test_df[i, :subEventName] = string(subset_dic[i]["subEventName"])
        test_df[i, :teamId] = subset_dic[i]["teamId"]
        test_df[i, :id] = subset_dic[i]["id"]
        test_df[i, :matchPeriod] = string(subset_dic[i]["matchPeriod"])
        test_df[i, :subEventId] = string(subset_dic[i]["subEventId"])
        test_df[i, :eventSec] = subset_dic[i]["eventSec"]
        test_df[i, :tags] = [convert.(Dict,subset_dic[i]["tags"])]
    end

    return test_df
end

"""
    function events()
Return a Dataframe with all events in the game
"""
function events(events_Italy, game_id::Int)
    subset_test = subset_indexes_lazy(events_Italy, game_id)
    result = transform_to_df(subset_test)
    data = convert_events(result)
    return data
end

@time events(v, 2576337)


text:
# Setting defaults
settings = settings_simple(
    z_ex = range(0.0, stop = 5.0, length = 100),
    ts = range(0.0, stop = 100, length = 10)
) # everything unmentioned is default

params = parameters_simple(
    υ = 0.1,
    θ = 2.1,
    ζ = 14.5,
    r = 0.05,
    ξ = 1.0
)

params_func = parameters_simple(
    μ = 0.0,
    υ = 0.1,
    θ = 2.1,
    r = t -> 0.05,
    x = t -> 14.5,
    π = (t, z) -> 1)

params_func_varying_1 = merge(params_func, (π = (t, z) -> 1 + 1e-02 * (1-t)/settings.T,))

params_func_varying_2 = merge(params_func, (r = t -> 0.05 - 1e-02 * (1 - t / settings.T),))

params_func_varying_3 = merge(params_func, (r = t -> 0.05 - 1e-02 * (1 - t / settings.T), π = (t, z) -> 1 + 1e-02 * (1-t)/settings.T))

# Solutions.
# Test the stationary residual is close to zero.
    # even by solving with DAE
    residuals, v_ts, g_ts = solve_simple_transition(params_func, settings)
    @test norm(residuals[1]) ≈ 0 atol = 1e-5
    @test norm(residuals[end]) ≈ 0 atol = 1e-5
    @test norm(residuals) ≈ 0 atol = 1e-5

# Solve with time-varying r and π, now with DAE
    # First experiment.
    residuals, v_ts, g_ts = solve_simple_transition(params_func_varying_1, settings)
    @test norm(residuals[1]) ≈ 0 atol = 1e-5
    @test norm(residuals[end]) ≈ 0 atol = 1e-5
    @test norm(residuals) ≈ 0 atol = 1e-5

    # Second experiment.
    residuals, v_ts, g_ts = solve_simple_transition(params_func_varying_2, settings)
    @test norm(residuals[1]) ≈ 0 atol = 1e-5
    @test norm(residuals[end]) ≈ 0 atol = 1e-5
    @test norm(residuals) ≈ 0 atol = 1e-5

    # Third experiment.
    residuals, v_ts, g_ts = solve_simple_transition(params_func_varying_3, settings)
    @test norm(residuals[1]) ≈ 0 atol = 1e-5
    @test norm(residuals[end]) ≈ 0 atol = 1e-5
    @test norm(residuals) ≈ 0 atol = 1e-5

ext: jl


text:
# this file is to run JWAS
using DataFrames,CSV,Random,Distributions,JWAS,DelimitedFiles, LinearAlgebra
Random.seed!(123)

method            = ARGS[1]
rep               = parse(Int, ARGS[2])  #rep1,...,rep10
pct               = parse(Float64, ARGS[3])

@show method,rep,pct

mainpath="/group/qtlchenggrp/tianjing/omics_sample_all/nonlinear/varmi_2/missing_pattern_legarra/"
cd(mainpath*"$pct/$method/rep$rep/")

############ READ DATA ##########
#read missing index
missingIndex=Int.(vec(readdlm(mainpath*"missing_index/missingIndex.pct$pct.rep$rep.txt")))
@show missingIndex[1:3]
@show length(missingIndex)

#read phenotype
phenofile="/group/qtlchenggrp/tianjing/omics_data/nonlinear/varmi_2/data1/y.csv"
phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])

#read omics
omicsfile  = "/group/qtlchenggrp/tianjing/omics_data/nonlinear/varmi_2/data1/omics.csv"
omics  = CSV.read(omicsfile,DataFrame,delim = ',',header=true,missingstrings=["NA"]);
insertcols!(omics,2,:y => phenotypes[:,:y]);
for i in 2:size(omics,2) #skip ID
    omics[!,i]= convert(Vector{Union{Missing,Float64}}, omics[!,i]);
end
omics[956:1055,:y] .= missing  #testing dataset
omics[missingIndex,3:end].=missing; #set missing patterns
@show [count(ismissing,col) for col in eachcol(omics)]

omics_name=["gene$i" for i in 1:1200]

@show size(omics)
@show omics[1:5,1:5];
@show omics[end-4:end,1:5];

#read genotype
if occursin("gblup",method)  #NN-GBLUP
	genofile   = "/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/GRM_geno.csv"
	genotypes  = get_genotypes(genofile;separator=',',method="GBLUP",header=false);
end

model_equation = "y = intercept + genotypes";

model = build_model(model_equation;
                    num_hidden_nodes=1200,
                    latent_traits=omics_name,
                    nonlinear_function="sigmoid"); #"linear" for linear activation function
out   = runMCMC(model,omics,chain_length=5000,printout_model_info=false);

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
using Weave;

list_out_formats()

weave("Operations Research\\01_Introduciton.jmd", doctype=:"md2pdf", out_path=:pwd)
weave("02_SimpleLinearOptimization.jmd", doctype=:"github", out_path=:pwd)
weave("03_Julia-Basics.jmd", doctype=:"github", out_path=:pwd)
weave("03_Gadfly.jmd", doctype=:"github", out_path=:pwd)
weave("04_Numerical-Methods.jmd", doctype=:"github", out_path=:pwd)

max_stars_repo_path: Operations Research/build.jl

max_stars_repo_name: bmoretz/Statistical-Computing

max_issues_repo_path: Operations Research/build.jl

max_issues_repo_name: bmoretz/Statistical-Computing

max_forks_repo_path: Operations Research/build.jl

max_forks_repo_name: bmoretz/Statistical-Computing


text:
## Physics:
##
## We simulate 2d flow of ideal gas through a porous medium
## with the help of Darcy's law.
##
## v⃗ = -μ⁻¹ K̂ ⋅∇P,
##
## where K̂, a second-order tensor in the general case, is
## the specific permeability. It depends only on the
## geometry of the mediumm. Here we assume isotropy
## of space, so K is a scalar. μ is the dynamic viscocity.
##
## Continuity equation:
##
## φ ⋅∂ρ / ∂t + div(ρv⃗) = 0
##
## where φ is the porosity.
##
## So, the equations used are the continuity equation and
## Darcy's equation, which replaces the momentum equation,
## used in our last program.
##
## Another change with respect to the ideal gas program
## is that we will set the boundary conditions using ghost cells.
## This way, we can make use of the central difference scheme for
## boundary conditions which consist of partial derivatives,
## which is second-order.
##
## The boundary conditions consist of walls on the top
## and bottom, and constant pressure on either side.
##
## For now, our velocities will be equal to 0 at the beginning.
##
## ---------
## ->     ->
## ->     ->
## ---------
##
## Initial conditions:
## v⃗ = (u, v) = 0⃗ everywhere
## P = P0 everywhere else
##
## Boundary conditions:
## u = 0 at x = 0, 2 (walls)
## P = Pin at y = 0 and x ∈ [0, 1]
## v = 0 at y = 0 and x ∈ [1, 2]
## P = Pout, dv / dy = 0 at y = 2
## ∂P / ∂y = 0 at y = 0, 2

using DelimitedFiles, PyPlot

function update!(; ρnext, unext, vnext, ρwork, uwork, vwork,
        Δt, Δx, Δy, Pin, Pout, φ, K, μ)
    nx, ny = size(ρnext)
    coeff = 8.314 * 298 / 0.029
    for j = 2 : ny - 1
        for i = 2 : nx - 1
            ## Auxiliary variables to avoid typos when
            ## working with partial derivatives
            ##
            ## NW--N--NE
            ##     |
            ## W---C---E
            ##     |
            ## SW--S--SE

            ρc = ρwork[i, j]
            ρe = ρwork[i + 1, j]
            ρw = ρwork[i - 1, j]
            ρn = ρwork[i, j + 1]
            ρs = ρwork[i, j - 1]

            uc = uwork[i, j]
            ue = uwork[i + 1, j]
            uw = uwork[i - 1, j]

            vc = vwork[i, j]
            vn = vwork[i, j + 1]
            vs = vwork[i, j - 1]

            ## Continuity equation to calculate ρ
            ## φ ⋅∂ρ / ∂t + div(ρv⃗) = 0
            ρnext[i, j] = ρc - Δt / φ *
                (uc * (ρe - ρw) / (2 * Δx) +
                 vc * (ρn - ρs) / (2 * Δy) +
                 ρc * (ue - uw) / (2 * Δx) +
                 ρc * (vn - vs) / (2 * Δy))
        end
    end

    ## Boundary conditions on pressure using ghost cells
    ## P = Pin at y = 0 and x ∈ [0, 1]
    @. @views ρnext[1: nx ÷ 2, 1] = 2 * Pin / coeff - ρnext[1: nx ÷ 2, 2]

    ## P = Pout at y = 2
    @. @views ρnext[:, ny] = 2 * Pout / coeff - ρnext[:, ny - 1]

    ## ∂P / ∂x = 0 at x = 0, 2
    @. @views ρnext[nx, :] = ρnext[nx - 1, :]
    @. @views ρnext[1, :] = ρnext[2, :]

    ## ∂P / ∂y = 0 at y = 0 x ∈ [1, 2]
    @. @views ρnext[nx ÷ 2 + 1: nx, 1] = ρnext[nx ÷ 2 + 1: nx, 2]

    ## Here we use a central difference scheme to discretize the second
    ## order derivatives, so we iterate over the 'inner' grid
    for j = 2 : ny - 1
        for i = 2 : nx - 1
            ## Auxiliar variables
            ## можно через (ρwork + ρnext) / 2
            ρe = ρwork[i + 1, j]
            ρw = ρwork[i - 1, j]
            ρn = ρwork[i, j + 1]
            ρs = ρwork[i, j - 1]

            ## Darcy's law
            ## v⃗ = -μ⁻¹ K̂ ⋅∇P
            unext[i, j] = -K / μ * coeff * (ρe - ρw) / (2 * Δx)
            vnext[i, j] = -K / μ * coeff * (ρn - ρs) / (2 * Δy)
        end
    end

    ## Boundary conditions on velocity, again using ghost cells
    ## u = 0 at x = 0, 2 (walls)
    @. @views unext[1, :] = -unext[2, :]
    @. @views unext[nx, :] = -unext[nx - 1, :]

    ## v = 0 at y = 0 and x ∈ [1, 2]
    @. @views vnext[1, nx ÷ 2 + 1: nx] = -vnext[2, nx ÷ 2 + 1: nx]

    ## dv / dy = 0 at y = 2
    @views vnext[:, ny] = vnext[:, ny - 1]

    ## dv / dy = 0 at y = 0 and x ∈ [0, 1]
    @views vnext[1: nx ÷ 2, 1] =  vnext[1: nx ÷ 2, 2]
end

function filtration!(; ρ, u, v, Δx, Δy, Δt, nsteps, K, φ, μ, Pin, Pout)
    ## Helper arrays to store values from previous iteration
    ρwork, uwork, vwork = ρ, u, v
    unext = copy(u)
    vnext = copy(v)
    ρnext = copy(ρ)

    for t = 1 : nsteps
        update!(; ρnext, unext, vnext, ρwork, uwork, vwork,
                Δt, Δx, Δy, Pin, Pout, φ, K, μ)

        ## Swap next and current and continue iteration
        ρnext, ρwork = ρwork, ρnext
        unext, uwork = uwork, unext
        vnext, vwork = vwork, vnext
    end

    ## If in the end ρwork is not the original ρ,
    ## copy the most recent data into ρ, u, v
    if ρwork !== ρ
        ρ .= ρwork
        u .= uwork
        v .= vwork
    end

    return ρ, u, v
end

function plot(u, v, ρ)
    fig = PyPlot.figure(figsize=(10, 10))
    ax = PyPlot.axes()
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    #legend(loc="upper right",fancybox="true")

    nx, ny= size(u)

    rgx = range(0, 2, length=nx)
    rgy = range(0, 2, length=ny)
    x = [rgx;]'
    y = [rgy;]'
    X = repeat([rgx;]', length(x))
    Y = repeat([rgy;],1, length(y))

    ## Velocity vector field
    ## scale = 100.0
    quiver(x, y, u', v')

    ## Streamplot of velocity vector field
    #streamplot(x, y, u', v')

    ## Density contour map
    pos = ax.contourf(X, Y, ρ', alpha=0.5, cmap=matplotlib.cm.viridis)
    fig.colorbar(pos, ax=ax)
    cp = contour(X, Y, ρ', cmap=matplotlib.cm.viridis)

    PyPlot.title("2d filtration of ideal gas through a porous medium")

    # PyPlot.svg(true)

    savefig("img/porous-medium-darcy.png")
    savefig("img/porous-medium-darcy.svg")
end

let
    duration = 2000
    Δt = 0.001
    nsteps = round(Int64, duration / Δt)
    # nsteps = 1000

    ## Space grid [0, 2] × [0, 2]
    Δx = 0.05
    Δy = Δx
    nx = convert(Int64, 2 / Δx)
    ny = convert(Int64, 2 / Δy)

    ## Porosity
    φ = 0.7

    ## Specific permeability
    K = 1e-12

    ## Dynamic viscocity (for air)
    μ = 18e-6

    ## P = ρRT / μ, (μ - molar mass)
    coeff = 8.314 * 298 / 0.029

    ## Pressures at top and bottom, and initial pressure
    Pin = 1e6
    Pout = 1e5
    P0 = Pout

    u = zeros(nx, ny)
    v = zeros(nx, ny)
    ρ = zeros(nx, ny)

    ρ .= P0 / coeff

    problem = (
               ρ = ρ,
               u = u,
               v = v,
               Δx = Δx,
               Δy = Δy,
               Δt = Δt,
               nsteps = nsteps,
               K = K,
               μ = μ,
               φ = φ,
               Pin = Pin,
               Pout = Pout
              )

    ## Run our solver
    ρ, u, v = filtration!(; problem...)

    # writedlm("density.txt", ρ, ' ')

    ## Plotting
    plot(u, v, ρ)
end

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
function asap3(style = Rectangular(2,1))
    toplevel = TopLevel{2}("asap3")
    toplevel.metadata["style"] = style

    ####################
    # Normal Processor #
    ####################

    # Get a processor tile and instantiate it.
    processor = build_processor_tile(style)
    for r in 0:29, c in 0:31
        add_child(toplevel, processor, CartesianIndex(r,c))
    end
    for r in 30:31, c in 12:19
        add_child(toplevel, processor, CartesianIndex(r,c))
    end

	####################
	# Memory Processor #
	####################

	memory_processor = build_processor_tile(style, include_memory = true)
	for r = 30, c = 0:11
        add_child(toplevel, memory_processor, CartesianIndex(r,c))
	end
	for r = 30, c = 20:31
        add_child(toplevel, memory_processor, CartesianIndex(r,c))
	end

	#################
	# 2 Port Memory #
	#################
    memory_2port = build_memory(2)
	for r = 31, c in (0:2:10)
        add_child(toplevel, memory_2port, CartesianIndex(r,c))
	end
	for r = 31, c in (20:2:30)
        add_child(toplevel, memory_2port, CartesianIndex(r,c))
	end

	#################
	# Input Handler #
	#################
	input_handler = build_input_handler(style)
    add_child(toplevel, input_handler, CartesianIndex(0,-1))

	##################
	# Output Handler #
	##################
	output_handler = build_output_handler(style)
    add_child(toplevel, output_handler, CartesianIndex(0,32))

	connect_processors(toplevel, style)
    connect_io(toplevel, style)
	connect_memories(toplevel, style)
	return toplevel
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Solids (Modia3D/Solids/_module.jl)
#
# Content:
#
#    Mass, centerOfMass, inertiaMatrix of a solid

struct MassProperties <: Modia3D.AbstractMassProperties
   m::Float64                 # mass in [kg]
   rCM::ModiaMath.Vector3D    # center of mass in [m]
   I::SMatrix{3,3,Float64,9}  # inertia matrix in [kg.m^2]

   function MassProperties(mass::Number, centerOfMass::AbstractVector, inertiaMatrix::AbstractMatrix)
      @assert(mass > 0.0)
      new(mass, centerOfMass, inertiaMatrix)
   end
end

MassProperties(;m::Number=0.0, rCM=ModiaMath.ZeroVector3D,
       Ixx::Number=0.0, Iyy::Number=0.0, Izz::Number=0.0,
       Ixy::Number=0.0, Ixz::Number=0.0, Iyz::Number=0.0) = MassProperties(m, rCM, [Ixx Ixy Ixz;
                                                                                Ixy Iyy Iyz;
                                                                                Ixz Iyz Izz])

MassProperties(geo::Geo, mass::Number)                 where Geo <: Modia3D.AbstractSolidGeometry = MassProperties(mass, centroid(geo), inertiaMatrix(geo,mass))
MassProperties(geo::Geo, material::SolidMaterial)      where Geo <: Modia3D.AbstractSolidGeometry = MassProperties(geo, material.density*volume(geo))
MassProperties(geo::Geo, materialName::AbstractString) where Geo <: Modia3D.AbstractSolidGeometry = MassProperties(geo, solidMaterialPalette[materialName])

mutable struct InternalMassProperties <: Modia3D.AbstractMassProperties
   m::Float64                 # mass in [kg]
   rCM::ModiaMath.Vector3D    # center of mass in [m]
   I::SMatrix{3,3,Float64,9}  # inertia matrix in [kg.m^2]
   # Basics.nullMRotation()
   InternalMassProperties() = new(0.0, ModiaMath.ZeroVector3D, SMatrix{3,3,Float64,9}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenLSTO_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenLSTO")
JLLWrappers.@generate_main_file("OpenLSTO", UUID("a318411f-452f-5433-884b-1f6a35676cea"))
end  # module OpenLSTO_jll


text:
using Distributions, CSV, DataFrames, LinearAlgebra, LaTeXStrings

# Load observed data
x_samples_df = CSV.read("gmm_observations.csv")
x_samples = hcat(x_samples_df[!, :X1], x_samples_df[!, :X2])

# Computes πᵢ for all i with respect to the given parameters in θ
# Inputs:
# - θ: A vector containing the current estimates of [α₀, μ₀, Σ₀, μ₁, Σ₁]
# - x: The (n x 2) matrix of observed X values
# Ouputs:
# 𝛑: An n dimensional vector where the ith coordinate is the value of πᵢ
function compute_all_πᵢ(θ, x)
    # Get dataset size
    n = size(x, 1)

    # Get current parameters
    α, μ₀, Σ₀, μ₁, Σ₁ = θ

    # Set distributions
    𝒩₁ = MvNormal(μ₁, Σ₁)
    𝒩₀ = MvNormal(μ₀, Σ₀)

    # Compute all πᵢ
    𝛑 = [α*pdf(𝒩₁, x[i,:])/(α*pdf(𝒩₁, x[i,:]) + (1-α)*pdf(𝒩₀, x[i,:])) for i=1:n]

    return 𝛑
end

# EM algorithm for GMM with observations x and Z ∈ {0,1}
# Inputs:
# - θ₀: The initial estimate of the θ parameter
# - x: The observed X values
# - iterations: The number of iterations for the EM algorithm
# Outputs:
# - θⱼ: The estimate of θ after 1000 iterations
function gmm_EM(θ₀, x, iterations)
    # Get dataset size
    n = size(x, 1)

    # Store the value of θ on jth step
    θⱼ = θ₀

    for j=1:iterations
        # Get current πᵢ values for all i
        𝛑 = compute_all_πᵢ(θⱼ, x)

        # Store the vector containing 1 - 𝛑
        𝛑₁ = 1 .- 𝛑

        # Get next value of α
        α = sum(𝛑)/n

        # Get next value of μ₀
        μ₀ = [dot(𝛑₁, x[:,1]), dot(𝛑₁, x[:,2])]/sum(𝛑₁)

        # Get next value of μ₁
        μ₁ = [dot(𝛑, x[:,1]), dot(𝛑, x[:,2])]/sum(𝛑)

        # Get next value of Σ₀.
        Σ₀ = sum([(x[i,:] - μ₀)*((x[i,:] - μ₀))'*𝛑₁[i] for i=1:n])/sum(𝛑₁)

        # Get next value of Σ₁.
        Σ₁ = sum([(x[i,:] - μ₁)*((x[i,:] - μ₁))'*𝛑[i] for i=1:n])/sum(𝛑)

        # Update θ
        θⱼ = [α, μ₀, Σ₀, μ₁, Σ₁]

    end

    return θⱼ
end

# Initialize EM parameters
α₀ = 0.3
μ₀ = [1.0, 2.0]
μ₁ = [2.0, 3.0]
Σ₀ = 1.0*Matrix(I, 2, 2)
Σ₁ = 1.0*Matrix(I, 2, 2)
θ₀ = [α₀, μ₀, Σ₀, μ₁, Σ₁]

println(gmm_EM(θ₀, x_samples, 20))

## Plots ##
gr(size=(300,300))

function ellipse!(μ,Σ;kw...)
    for i=1:3
        plot!([tuple((i*sqrt(Σ)*[cos(t),sin(t)]+μ)...) for t in range(0,stop=2π,length=100)];
              seriestype=:shape,linealpha=0.5,linecolor=:gray,fillalpha=[0.4,0.2,0.1][i],kw...)
    end
    current()
end

function mixtureplot(X₁,X₂,μ₀,Σ₀,μ₁,Σ₁,Π, plot_title)
    scatter(X₁,X₂;fillalpha=0.5,markerstrokewidth=0.5,
        marker_z=Π,mc=ColorGradient([:lightblue,:orange]),
        colorbar=:false, title=plot_title)
    ellipse!(μ₀,Σ₀,fillcolor=:lightblue)
    ellipse!(μ₁,Σ₁,fillcolor=:orange)
    plot!(;bg=:transparent,xlims=(-4.5,8.5),ylims=(-3.5,10.5),
        leg=false,ticks=:none,ratio=:equal)
end

one_step = gmm_EM(θ₀, x_samples, 1)
one_step_Π = compute_all_πᵢ(one_step, x_samples)
five_steps = gmm_EM(θ₀, x_samples, 5)
five_step_Π = compute_all_πᵢ(ten_steps, x_samples)
twenty_steps = gmm_EM(θ₀, x_samples, 20)
twenty_step_Π = compute_all_πᵢ(twenty_steps, x_samples)

# Zero step plot
mixtureplot(x_samples[:,1], x_samples[:,2], θ₀[2], θ₀[3], θ₀[4], θ₀[5], compute_all_πᵢ(θ₀, x_samples), L"\theta_0")
savefig("../images/gmm_mixtureplot_0.svg")

# One step plot
mixtureplot(x_samples[:,1], x_samples[:,2], one_step[2], one_step[3], one_step[4], one_step[5], one_step_Π,  L"\theta_1")
savefig("../images/gmm_mixtureplot_1.svg")

# Five step plot
mixtureplot(x_samples[:,1], x_samples[:,2], five_steps[2], five_steps[3], five_steps[4], five_steps[5], five_step_Π, L"\theta_5")
savefig("../images/gmm_mixtureplot_5.svg")

# Twenty step plot
mixtureplot(x_samples[:,1], x_samples[:,2], twenty_steps[2], twenty_steps[3], twenty_steps[4], twenty_steps[5], twenty_step_Π,  L"\theta_{20}")
savefig("../images/gmm_mixtureplot_20.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using SplittablesTesting
using Test

@testset "test_ordered" begin
    SplittablesTesting.test_ordered([
        (label = "First Test", data = 1:5),
        (label = "Second Test", data = (a = 1, b = 2, c = 3)),
        zip(1:3, 4:6),
    ])
end

@testset "test_unordered" begin
    SplittablesTesting.test_unordered([
        (label = "First Test", data = 1:5),
        (label = "Second Test", data = (a = 1, b = 2, c = 3)),
        zip(1:3, 4:6),
    ])
end


text:
int main()
{
   int a[13];
   return a[14] + a[-1];
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
using Conda

Conda.add("pymatgen=2022.3.7", channel="conda-forge")


text:
using BenchmarkTools
using CSVReader
using DataFrames
using Dates
using Random
using Test

const iris = "iris.csv"

function testiris(df)
    @test size(df) == (150, 6)
    @test names(df) == ["id", "Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width", "Species"]
    @test unique(df[!, :Species]) == ["setosa", "versicolor", "virginica"]
    @test [sum(df[!, i]) for i in 2:5] ≈ [876.5, 458.6, 563.7, 179.9]
end

function genperf(rows, nfloat, nstr, strsize, nqfloat, nqstr, qstrsize)
    filename = "perf_f$(nfloat)_qf$(nqfloat)_s$(nstr),$(strsize)_qs$(nqstr),$(qstrsize)_$(rows).csv"
    open(filename, "w") do f
        for i in 1:rows
            floats = rand(nfloat)
            qfloats = ['"' * string(rand()) * '"' for _ in 1:nqfloat]
            strs  = [randstring(strsize) for _ in 1:nstr]
            qstrs = ['"' * randstring(qstrsize) * '"' for _ in 1:nqstr]
            println(f, join(vcat(floats, qfloats, strs, qstrs), ","))
        end
    end
    println(Dates.now(), " generated ", filename)
    filename
end

@testset "CSVReader" begin

    # readers available for testing
    readers = [CSVReader.read_csv]

    for reader in readers
        # basic usage
        df = reader("iris.csv")
        testiris(df)

        # specified parsers
        df = reader("iris.csv", parsers"i,f:4,s")
        testiris(df)

        # without headers
        df = reader("iris2.csv"; headers = false)
        rename!(df, Symbol.(["id", "Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width", "Species"]))
        testiris(df)

        # missing data (1st row ok)
        df = reader("missing1.csv")
        @test sum(count(ismissing, df[!,c]) for c in 1:ncol(df)) == 5

        # missing  data (1st row contains missing, causing column to have type String)
        df = reader("missing2.csv")
        @test sum(count(ismissing, df[!,c]) for c in 1:ncol(df)) == 3

        # missing  data (1st row contains missing, causing column to have type String)
        # but, this is corrected via parsers spec
        df = reader("missing2.csv", parsers"i:4")
        @test sum(count(ismissing, df[!,c]) for c in 1:ncol(df)) == 5
    end

    # Performance test - generate test files.
    perf_files = String[]
    try
        for rows in [10, 1000], nfloat = [3, 5], nstr = [3, 5]
            push!(perf_files, genperf(rows, nfloat, nstr, 10, 0, 0, 0))
        end
        for file in perf_files
            for reader in readers
                println(Dates.now(), " [$(reader)] performance testing $file")
                @btime $reader($file)
            end
        end
    finally
        println(Dates.now(), " cleaning up")
        rm.(perf_files)
    end

end


text:
using Adapt
using Statistics
using KernelAbstractions: @kernel, @index, Event
using Oceananigans.Grids

"""
    struct ComputedField{X, Y, Z, A, G, O} <: AbstractField{X, Y, Z, A, G}

Type representing a field computed from an operand.
"""
struct ComputedField{X, Y, Z, S, A, G, O} <: AbstractField{X, Y, Z, A, G}
       data :: A
       grid :: G
    operand :: O
     status :: S

    function ComputedField{X, Y, Z}(data, grid, operand; recompute_safely=true) where {X, Y, Z}
        validate_field_data(X, Y, Z, data, grid)

        # Use FieldStatus if we want to avoid always recomputing
        status = recompute_safely ? nothing : FieldStatus(0.0)

        return new{X, Y, Z, typeof(status), typeof(data),
                   typeof(grid), typeof(operand)}(data, grid, operand, status)
    end

    function ComputedField{X, Y, Z}(data, grid, operand, status) where {X, Y, Z}
        validate_field_data(X, Y, Z, data, grid)

        return new{X, Y, Z, typeof(status), typeof(data),
                   typeof(grid), typeof(operand)}(data, grid, operand, status)
    end
end

"""
    ComputedField(operand; data=nothing)

Returns a field whose data is `computed` from `operand`.
If the keyword argument `data` is not provided, memory is allocated to store
the result. The `arch`itecture of `data` is inferred from `operand`.
"""
function ComputedField(operand; data=nothing, recompute_safely=true)

    loc = location(operand)
    arch = architecture(operand)
    grid = operand.grid

    if isnothing(data)
        data = new_data(arch, grid, loc)
        recompute_safely = false
    end

    return ComputedField{loc[1], loc[2], loc[3]}(data, grid, operand; recompute_safely=recompute_safely)
end

"""
    compute!(comp::ComputedField)

Compute `comp.operand` and store the result in `comp.data`.
"""
function compute!(comp::ComputedField{X, Y, Z}, time=nothing) where {X, Y, Z}
    compute_at!(comp.operand, time) # ensures any 'dependencies' of the computation are computed first

    arch = architecture(comp.data)

    workgroup, worksize = work_layout(comp.grid,
                                      :xyz,
                                      include_right_boundaries=true,
                                      location=(X, Y, Z))

    compute_kernel! = _compute!(device(arch), workgroup, worksize)

    event = compute_kernel!(comp.data, comp.operand; dependencies=Event(device(arch)))

    wait(device(arch), event)

    return nothing
end

compute_at!(field::ComputedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(field, time)

"""Compute an `operand` and store in `data`."""
@kernel function _compute!(data, operand)
    i, j, k = @index(Global, NTuple)
    @inbounds data[i, j, k] = operand[i, j, k]
end

#####
##### Adapt
#####

Adapt.adapt_structure(to, computed_field::ComputedField) = Adapt.adapt(to, computed_field.data)


text:
using DataFrames, Gadfly
set_default_plot_size(6inch, 3inch)

labs = [ "exp", "sqrt", "log", "linear"]
funcs = [ x->60*(1 .- exp.(-0.2*x)), x->sqrt.(x)*10, x->log.(x)*10, x->x*0.6 ]
x = [1.0:30;]
D = vcat([DataFrame(x=x, y=f(x), linev=l) for (f,l) in zip(funcs, labs)]...)

p1 = plot(D, x=:x, y=:y, linestyle=:linev, Geom.line )
p2 = plot(D, x=:x, y=:y, linestyle=:linev, Geom.line,
   Scale.linestyle_discrete(levels=["exp", "log", "sqrt", "linear"]) )
hstack(p1,p2)


text:
# This file includes to blocks used in a digital communication system.

export Detector, isml

"""
    $TYPEDEF

Baseband optimal detector.

!!! note
    If elements of `probs` are equal to each other, than detector detects symbols in the maximum likelihood (ML)
    mode. Otherwise, it detects in the maximum a posterior probabilities (MAP) mode.

# Fields

    $TYPEDFIELDS
"""
struct Detector{T}
    "Basis reference of the detector"
    refs::T
    "A priori probabilities of the message symbols"
    probs::Vector{Float64}
    "2 times the power spectral density of channel noise"
    N0::Float64
end
Detector(refs) = Detector(refs, fill(1 / length(refs), length(refs)), NaN)

"""
    $SIGNATURES

Returns true if `detector` is a maximum likelihood detector, i.e, a priori probabilities of its reference refs are
equal.
"""
isml(detector::Detector) = (probs = detector.probs; all(probs .≈ probs[1]))

(detector::Detector)(rx) = isml(detector) ? mldetect(detector, rx) : mapdetect(detector, rx)

function mapdetect(detector::Detector, rx)
    Es = 1 / 2 * energy.(detector.refs)
    Pm = detector.N0 * log.(detector.probs)
    map(rx) do ri
        argmax(map(s -> real(ri ⋅ s), detector.refs) - Es + Pm)
    end
end

function mldetect(detector::Detector, rx)
    Es = 1 / 2 * energy.(detector.refs)
    map(rx) do ri
        argmax(map(s -> real(ri ⋅ s), detector.refs) - Es)
    end
end


text:
##########
##########
# Create series of networks to run through simulations
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-12-04
##########
##########

# set up =======================================================================

# set packages
using DelimitedFiles

include("../01_species_list_creation.jl")

# create different combinations of networks ====================================

# the different values that have different options are:
#   richness shape
#   biomass shape

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# set parameters
n = 50
richness_shape = "top-heavy"
species_list, network = network_draw(n, richness_shape)
# draw species list
species_list = define_species_list(n, richness_shape)

# test create_interaction_matrix()
interaction_matrix = fill(0, (size(species_list, 1), size(species_list, 1)))
    for i in 1:size(interaction_matrix)[1] #rows
        for j in 1:size(interaction_matrix)[2] #cols
            interaction_matrix[i,j] = determine_interaction(j, i, species_list)
        end
    end
    return interaction_matrix
interaction_matrix = create_interaction_matrix(species_list)
sums = minimum(sum(interaction_matrix, dims = 2) +
            (sum(interaction_matrix, dims = 1)'))
if sums == 0
    # go through while loop to re-draw until connected
    connect_acomplish = false
    connected_counter = 0
    while connect_acomplish == false
        connected_counter += 1
        # draw species list
        species_list = define_species_list(n, richness_shape)
        interaction_matrix = create_interaction_matrix(species_list)

        # turn into graph
        sums = minimum(sum(interaction_matrix, dims = 2) +
            (sum(interaction_matrix, dims = 1)'))
        if sums == 0
            continue
        end

        if sums > 0
            connect_acomplish = true
        end
    end
end
return species_list, interaction_matrix

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in pairs(kwargs) if k in keys)
end

extract_fields(x) = Dict(String(name) => getfield(x, name) for name in fieldnames(typeof(x)))


text:
"""

  DepthwiseConvAllowNinChangeStrategy(newoutputsmax::Integer, multipliersmax::Integer, base, [fallback])
  DepthwiseConvAllowNinChangeStrategy(allowed_new_outgroups::AbstractVector{<:Integer}, allowed_multipliers::AbstractVector{<:Integer}, base, [fallback])

`DecoratingJuMPΔSizeStrategy` which allows both nin and nout of `DepthwiseConv` layers to change independently.

Might cause optimization to take very long time so use with care! Use [`DepthwiseConvSimpleΔSizeStrategy`](@ref)
if `DepthwiseConvAllowNinChangeStrategy` takes too long.

The elements of `allowed_new_outgroups` determine how many extra elements in the output dimension of the weight
shall be tried for each existing output element. For example, for a `DepthwiseConv((k1,k2), nin=>nout))` there
are `nout / nin` elements in the output dimension. With `allowed_new_outgroups = 0:3` it is allowed to insert
0, 1, 2 or 3 new elements in the output dimension between each already existing element (so with `nout / nin`
elements the maximum increase is `3 * nout / nin`).

The elements of `allowed_multipliers` determine the total number of allowed output elements, i.e the allowed
ratios of `nout / nin`.

If `fallback` is not provided, it will be derived from `base`.
"""
struct DepthwiseConvAllowNinChangeStrategy{S,F} <: DecoratingJuMPΔSizeStrategy
  allowed_new_outgroups::Vector{Int}
  allowed_multipliers::Vector{Int}
  base::S
  fallback::F
end
DepthwiseConvAllowNinChangeStrategy(newoutputsmax::Integer, multipliersmax::Integer,base,fb...) = DepthwiseConvAllowNinChangeStrategy(0:newoutputsmax, 1:multipliersmax, base, fb...)

function DepthwiseConvAllowNinChangeStrategy(
  allowed_new_outgroups::AbstractVector{<:Integer},
  allowed_multipliers::AbstractVector{<:Integer},
  base, fb= recurse_fallback(s -> DepthwiseConvAllowNinChangeStrategy(allowed_new_outgroups, allowed_multipliers, s), base))
  return DepthwiseConvAllowNinChangeStrategy(collect(Int, allowed_new_outgroups), collect(Int, allowed_multipliers), base, fb)
end

NaiveNASlib.base(s::DepthwiseConvAllowNinChangeStrategy) = s.base
NaiveNASlib.fallback(s::DepthwiseConvAllowNinChangeStrategy) = s.fallback

NaiveNASlib.add_participants!(s::DepthwiseConvAllowNinChangeStrategy, vs=AbstractVertex[]) = NaiveNASlib.add_participants!(base(s), vs)

"""
  DepthwiseConvSimpleΔSizeStrategy(base, [fallback])

`DecoratingJuMPΔSizeStrategy` which only allows nout of `DepthwiseConv` layers to change.

Use if [`DepthwiseConvAllowNinChangeStrategy`](@ref) takes too long to solve.

The elements of `allowed_multipliers` determine the total number of allowed output elements, i.e the allowed
ratios of `nout / nin`.

If `fallback` is not provided, it will be derived from `base`.
"""
struct DepthwiseConvSimpleΔSizeStrategy{S, F} <: DecoratingJuMPΔSizeStrategy
  allowed_multipliers::Vector{Int}
  base::S
  fallback::F
end

DepthwiseConvSimpleΔSizeStrategy(maxms::Integer, base, fb...) = DepthwiseConvSimpleΔSizeStrategy(1:maxms, base, fb...)
function DepthwiseConvSimpleΔSizeStrategy(ms::AbstractVector{<:Integer}, base, fb=recurse_fallback(s -> DepthwiseConvSimpleΔSizeStrategy(ms, s), base))
  return DepthwiseConvSimpleΔSizeStrategy(collect(Int, ms), base, fb)
end
NaiveNASlib.base(s::DepthwiseConvSimpleΔSizeStrategy) = s.base
NaiveNASlib.fallback(s::DepthwiseConvSimpleΔSizeStrategy) = s.fallback

NaiveNASlib.add_participants!(s::DepthwiseConvSimpleΔSizeStrategy, vs=AbstractVertex[]) = NaiveNASlib.add_participants!(base(s), vs)

recurse_fallback(f, s::AbstractJuMPΔSizeStrategy) = wrap_fallback(f, NaiveNASlib.fallback(s))
recurse_fallback(f, s::NaiveNASlib.DefaultJuMPΔSizeStrategy) = s
recurse_fallback(f, s::NaiveNASlib.ThrowΔSizeFailError) = s
recurse_fallback(f, s::NaiveNASlib.ΔSizeFailNoOp) = s
recurse_fallback(f, s::NaiveNASlib.LogΔSizeExec) = NaiveNASlib.LogΔSizeExec(s.msgfun, s.level, f(s.andthen))

wrap_fallback(f, s) = f(s)
wrap_fallback(f, s::NaiveNASlib.LogΔSizeExec) = NaiveNASlib.LogΔSizeExec(s.msgfun, s.level, f(s.andthen))

function NaiveNASlib.compconstraint!(case, s, ::FluxLayer, data) end
function NaiveNASlib.compconstraint!(case, s::DecoratingJuMPΔSizeStrategy, lt::FluxLayer, data)
   NaiveNASlib.compconstraint!(case, NaiveNASlib.base(s), lt, data)
end
# To avoid ambiguity
function NaiveNASlib.compconstraint!(case::NaiveNASlib.ScalarSize, s::DecoratingJuMPΔSizeStrategy, lt::FluxDepthwiseConv, data)
  NaiveNASlib.compconstraint!(case, NaiveNASlib.base(s), lt, data)
end
function NaiveNASlib.compconstraint!(::NaiveNASlib.ScalarSize, s::AbstractJuMPΔSizeStrategy, ::FluxDepthwiseConv, data, ms=allowed_multipliers(s))

  # Add constraint that nout(l) == n * nin(l) where n is integer
  ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))

   # "data.noutdict[data.vertex] == data.noutdict[ins[i]] * x" where x is an integer variable is not linear
   # Instead we use the good old big-M strategy to set up a series of "or" constraints (exactly one of them must be true).
   # This is combined with the abs value formulation to force equality.
   # multiplier[j] is false if and only if data.noutdict[data.vertex] == ms[j]*data.noutdict[ins[i]]
   # all but one of multiplier must be true
   # Each constraint below is trivially true if multiplier[j] is true since 1e6 is way bigger than the difference between the two variables
  multipliers = @variable(data.model, [1:length(ms)], Bin)
  @constraint(data.model, sum(multipliers) == length(multipliers)-1)
  @constraint(data.model, [i=1:length(ins),j=1:length(ms)], data.noutdict[data.vertex] - ms[j]*data.noutdict[ins[i]] + multipliers[j] * 1e6 >= 0)
  @constraint(data.model, [i=1:length(ins),j=1:length(ms)], data.noutdict[data.vertex] - ms[j]*data.noutdict[ins[i]] - multipliers[j] * 1e6 <= 0)

  # Inputs which does not have a variable, possibly because all_in_Δsize_graph did not consider it to be part of the set of vertices which may change
  # We will constrain data.vertex to have integer multiple of its current size
  fixedins = filter(vin -> vin ∉ ins, inputs(data.vertex))
  if !isempty(fixedins)
    fv_out = @variable(data.model, integer=true)
    @constraint(data.model, [i=1:length(fixedins)], data.noutdict[data.vertex] == nout(fixedins[i]) * fv_out)
  end
end

allowed_multipliers(s::DepthwiseConvAllowNinChangeStrategy) = s.allowed_multipliers
allowed_multipliers(s::DepthwiseConvSimpleΔSizeStrategy) = s.allowed_multipliers
allowed_multipliers(::AbstractJuMPΔSizeStrategy) = 1:10

function NaiveNASlib.compconstraint!(case::NaiveNASlib.NeuronIndices, s::DecoratingJuMPΔSizeStrategy, t::FluxDepthwiseConv, data)
  NaiveNASlib.compconstraint!(case, base(s), t, data)
end
function NaiveNASlib.compconstraint!(case::NaiveNASlib.NeuronIndices, s::AbstractJuMPΔSizeStrategy, t::FluxDepthwiseConv, data)
  # Fallbacks don't matter here since we won't call it from below here, just add default so we don't accidentally crash due to some
  # strategy which hasn't defined a fallback
  if 15 < sum(keys(data.outselectvars)) do v
      layertype(v) isa FluxDepthwiseConv || return 0
      return log2(nout(v)) # Very roughly determined...
  end
    return NaiveNASlib.compconstraint!(case, DepthwiseConvSimpleΔSizeStrategy(10, s, NaiveNASlib.DefaultJuMPΔSizeStrategy()), t, data)
  end
  # The number of allowed multipliers can probably be better tuned, perhaps based on current size.
  return NaiveNASlib.compconstraint!(case, DepthwiseConvAllowNinChangeStrategy(10, 10, s, NaiveNASlib.DefaultJuMPΔSizeStrategy()), t, data)
  #=
  For benchmarking:
    using NaiveNASflux, Flux, NaiveNASlib.Advanced
    function timedwc(ds, ws)
      for w in ws
        for d in ds
          iv = conv2dinputvertex("in", w)
          dv = reduce((v, i) -> fluxvertex("dv$i", DepthwiseConv((1,1), nout(v) => fld1(i, 3) * nout(v)), v), 1:d; init=iv)
          metric = sum(ancestors(dv)) do v
            layer(v) isa DepthwiseConv || return 0
            return log2(nout(v))
          end
          res = @timed Δsize!(ΔNoutRelaxed(outputs(iv)[1] => w; fallback = ΔSizeFailNoOp()))
          println((width=w, depth = d, metric=metric, time=res.time))
          res.time < 5 || break
        end
      end
    end
  =#
end

function NaiveNASlib.compconstraint!(::NaiveNASlib.NeuronIndices, s::DepthwiseConvSimpleΔSizeStrategy, t::FluxDepthwiseConv, data)
  model = data.model
  v = data.vertex
  select = data.outselectvars[v]
  insert = data.outinsertvars[v]

  nin(v)[] == 1 && return # Special case, no restrictions as we only need to be an integer multple of 1

  ngroups = div(nout(v), nin(v)[])
  # Neurons mapped to the same weight are interleaved, i.e layer.weight[:,:,1,:] maps to y[1:ngroups:end] where y = layer(x)
  for group in 1:ngroups
    neurons_in_group = select[group : ngroups : end]
    @constraint(model, neurons_in_group[1] == neurons_in_group[end])
    @constraint(model, [i=2:length(neurons_in_group)], neurons_in_group[i] == neurons_in_group[i-1])

    insert_in_group = insert[group : ngroups : end]
    @constraint(model, insert_in_group[1] == insert_in_group[end])
    @constraint(model, [i=2:length(insert_in_group)], insert_in_group[i] == insert_in_group[i-1])
  end

  NaiveNASlib.compconstraint!(NaiveNASlib.ScalarSize(), s, t, data, allowed_multipliers(s))
end

function NaiveNASlib.compconstraint!(case::NaiveNASlib.NeuronIndices, s::DepthwiseConvAllowNinChangeStrategy, t::FluxDepthwiseConv, data)
  model = data.model
  v = data.vertex
  select = data.outselectvars[v]
  insert = data.outinsertvars[v]

  nin(v)[] == 1 && return # Special case, no restrictions as we only need to be an integer multple of 1?

  # Step 1:
  # Neurons mapped to the same weight are interleaved, i.e layer.weight[:,:,1,:] maps to y[1:ngroups:end] where y = layer(x)
  # where ngroups = nout / nin. For example, nout = 12 and nin = 4 mean size(layer.weight) == (..,3, 4)
  # Thus we have ngroups = 3 and with groupsize of 4.
  # Now, if we were to change nin to 3, we would still have ngroups = 3 but with groupsize 3 and thus nout = 9
  #
  ins = filter(vin -> vin in keys(data.noutdict), inputs(v))
  # If inputs to v are not part of problem we have to keep nin(v) fixed!
  isempty(ins) && return NaiveNASlib.compconstraint!(case, DepthwiseConvSimpleΔSizeStrategy(allowed_multipliers(s), base(s)), t, data)
  # TODO: Check if input is immutable and do simple strat then too?
  inselect = data.outselectvars[ins[]]
  ininsert = data.outinsertvars[ins[]]

  #ngroups = div(nout(v), nin(v)[])
  ningroups = nin(v)[]
  add_depthwise_constraints(model, inselect, ininsert, select, insert, ningroups, s.allowed_new_outgroups, s.allowed_multipliers)
end

function add_depthwise_constraints(model, inselect, ininsert, select, insert, ningroups, allowed_new_outgroups, allowed_multipliers)

  # ningroups is the number of "input elements" in the weight array, last dimension at the time of writing
  # noutgroups is the number of "output elements" in the weight array, second to last dimension at time of writing
  # nin is == ningroups while nout == ningroups * noutgroups

  # When changing nin by Δ, we will get ningroups += Δ which just means that we insert/remove Δ input elements.
  # Inserting one new input element at position i will get us noutgroups new consecutive outputputs at position i
  # Thus nout change by Δ * noutgroups.

  # Example:

  # dc = DepthwiseConv((1,1), 3 => 9; bias=false);

  # size(dc.weight)
  # (1, 1, 3, 3)

  # indata = randn(Float32, 1,1,3, 1);

  # reshape(dc(indata),:)'
  # 1×9 adjoint(::Vector{Float32}) with eltype Float32:
  #  0.315234  0.267166  -0.227757  0.03523  0.0516438  -0.124007  -0.63409  0.0273361  -0.457146

  # Inserting zeros instead of adding/removing just as it is less verbose to do so.
  # dc.weight[:,:,:,2] .= 0;

  # reshape(dc(indata),:)'
  # 1×9 adjoint(::Vector{Float32}) with eltype Float32:
 #  0.315234  0.267166  -0.227757  0.0  0.0  0.0  -0.63409  0.0273361  -0.457146

  # When changing nout by Δ where Δ ==  k * ningroups where k is an integer we can simply just add/remove k output elements
  # Output elements are interleaved in the output of the layer, so if we add/remove output element i, the elements i:ningroups:end
  # will be added/dropped from the layer output.
  # When Δ is not an integer multiple things get pretty hairy and I'm not sure the constraints below handle all causes
  # We then basically need to change nin so that the new nout is an integer multiple of the new nin.

  # Example

  # dc = DepthwiseConv((1,1), 3 => 9; bias=false);

  # reshape(dc(indata),:)'
  # 1×9 adjoint(::Vector{Float32}) with eltype Float32:
  #  0.109549  0.273408  -0.299883  0.0131096  0.00240822  -0.0461159  -0.138605  -0.557134  0.0246111

  # julia> dc.weight[:,:,2,:] .= 0;

  # reshape(dc(indata),:)'
  # 1×9 adjoint(::Vector{Float32}) with eltype Float32:
  # 0.109549  0.0  -0.299883  0.0131096  0.0  -0.0461159  -0.138605  0.0  0.0246111

  # One thing which makes this a bit tricky is that while noutgruops and ningroups can change
  # we are here stuck with the select and insert arrays that correspond to the original sizes
  # As such, keep in mind that noutgroups and ningroups in this function always relate to the
  # lengths of the arrays of JuMP variables that we have right now.

  M = 1e6
  noutgroups = div(length(select), ningroups)

  # select_none_in_group[g] == 1 if we drop group g
  # TODO: Replace by just sum(select_in_group[g,:]) below?
  select_none_in_group = @variable(model, [1:noutgroups], Bin)
  select_in_group = reshape(select, noutgroups, ningroups)
  @constraint(model, sum(select_in_group, dims=2) .<= M .* (1 .- select_none_in_group))

  # Tie selected input indices to selected output indices. If we are to drop any layer output indices, we need to drop the whole group
  # as one whole group is tied to a single element in the output element dimension of the weight array.
  @constraint(model, [g=1:noutgroups, i=1:ningroups], inselect[i] - select_in_group[g, i] + select_none_in_group[g] * M >= 0)
  @constraint(model, [g=1:noutgroups, i=1:ningroups], inselect[i] - select_in_group[g, i] - select_none_in_group[g] * M <= 0)

  # This variable handles the constraint that if we add a new input, we need to add noutgroups new outputs
  # Note that for now, it is only connected to the insert variable through the sum constraint below.
  insert_new_inoutgroups = @variable(model, [[1], 1:ningroups], Int, lower_bound=0, upper_bound=ningroups * maximum(allowed_multipliers))

  insize = @expression(model, sum(inselect) + sum(ininsert))
  outsize = @expression(model, sum(select) + sum(insert))

  # inmultipliers[j] == 1 if nout(v) == allowed_multipliers[j] * nin(v)[]
  inmultipliers = @variable(model, [1:length(allowed_multipliers)], Bin)
  #SOS1 == Only one can be non-zero. Not strictly needed, but it seems like it speeds up the solver
  @constraint(model, inmultipliers in SOS1(1:length(inmultipliers)))
  @constraint(model, sum(inmultipliers) == 1)

  @constraint(model,[i=1:length(ininsert), j=1:length(inmultipliers)], allowed_multipliers[j] * ininsert[i] - insert_new_inoutgroups[1, i] + (1-inmultipliers[j]) * M >= 0)
  @constraint(model,[i=1:length(ininsert), j=1:length(inmultipliers)], allowed_multipliers[j] * ininsert[i] - insert_new_inoutgroups[1, i] - (1-inmultipliers[j]) * M <= 0)

  @constraint(model, [j=1:length(inmultipliers)], allowed_multipliers[j] * insize - outsize + (1-inmultipliers[j]) * M >= 0)
  @constraint(model, [j=1:length(inmultipliers)], allowed_multipliers[j] * insize - outsize - (1-inmultipliers[j]) * M <= 0)

  insert_new_inoutgroups_all_inds = vcat(zeros(noutgroups-1,ningroups), insert_new_inoutgroups)

  # This variable handles the constraint that if we want to increase the output size without changing
  # the input size, we need to insert whole outgroups since all noutgroup outputs are tied to a single
  # output element in the weight array.
  insert_new_outgroups = @variable(model, [1:noutgroups, 1:ningroups], Int, lower_bound = 0, upper_bound=noutgroups * maximum(allowed_new_outgroups))

  # insert_no_outgroup[g] == 1 if we don't insert a new output group after group g
  # TODO: Replace by just sum(insert_new_outgroups[g,:]) below?
  insert_no_outgroup = @variable(model, [1:noutgroups], Bin)
  @constraint(model, sum(insert_new_outgroups, dims=2) .<= M .* (1 .- insert_no_outgroup))

  # When adding a new output group, all inserts must be identical
  # If we don't add any, all inserts are just 0
  @constraint(model, [g=1:noutgroups, i=2:ningroups], insert_new_outgroups[g,i] == insert_new_outgroups[g,i-1])
  @constraint(model, [g=1:noutgroups], insert_new_outgroups[g,1] == insert_new_outgroups[g,end])

  # new_outgroup[g,j] == 1 if we are inserting allowed_new_outgroups[j] new output groups after output group g
  noutmults = 1:length(allowed_new_outgroups)
  new_outgroup = @variable(model, [1:noutgroups, noutmults], Bin)
  #SOS1 == Only one can be non-zero. Not strictly needed, but it seems like it speeds up the solver
  @constraint(model,[g=1:noutgroups], new_outgroup[g,:] in SOS1(1:length(allowed_new_outgroups)))
  @constraint(model,[g=1:noutgroups], sum(new_outgroup[g,:]) == 1)

  groupsum = @expression(model, [g=1:noutgroups], sum(insert_new_outgroups[g,:]) - sum(insert_new_inoutgroups_all_inds[g,:]))
  @constraint(model, [g=1:noutgroups, j=noutmults], groupsum[g] - allowed_new_outgroups[j]*insize + (1-new_outgroup[g,j] + insert_no_outgroup[g]) * M >= 0)
  @constraint(model, [g=1:noutgroups, j=noutmults], groupsum[g] - allowed_new_outgroups[j]*insize - (1-new_outgroup[g,j] + insert_no_outgroup[g]) * M <= 0)

  # Finally, we say what the inserts shall be: the sum of inserts from new inputs and new outputs
  # I think this allows for them to be somewhat independent, but there are still cases when they
  # can't change simultaneously. TODO: Check those cases
  @constraint(model, insert .== reshape(insert_new_outgroups, :) .+ reshape(insert_new_inoutgroups_all_inds,:))

end


text:
@testset "Solvers" begin
  grid = CartesianGrid(10,10)
  vars = (:Z1=>Float64,:Z2=>Float64,:Z3=>Float64)
  prob = SimulationProblem(grid, vars, 1)
  solv = SSolver()
  covars = covariables(prob, solv)
  @test covars[1].names == (:Z1,)
  @test covars[2].names == (:Z2,)
  @test covars[3].names == (:Z3,)
  sol = solve(prob, solv)
  @test sprint(show, sol) == "2D Ensemble"
  @test sprint(show, MIME"text/plain"(), sol) == "2D Ensemble\n  domain: 10×10 CartesianGrid{2,Float64}\n  variables: Z1, Z2 and Z3\n  N° reals:  1"
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + [markdown] toc=true
#
Table of Contents

#

0  曲線座標系での接ベクトルのお話

0.1  概要

0.2  必要な Julia のライブラリを導入

1  曲線座標系の導入

1.1  試しに描画してみる

1.1.1  対話的操作で確認

1.2  曲線による平面の点の特徴付け

1.2.1  確認のため描画

2  曲線座標系での接ベクトル

2.1  一つのパラメータを動かしてみる

2.2  接ベクトルの可視化

2.3  念のため検算

2.3.1  偏微分を計算できるようにする

2.3.2  p, q の型が実数などの数値を取れば ForwardDiffのロジックを使って微分計算をすることになる

2.3.3  SymPy で計算

3  ある方向に対しての接ベクトル

4  接空間が導入できた

5  いやでも，目で見えるじゃん，比較できそうじゃん

6  せやな！それを正当化するのが接続とか共変微分とかの話やで！

# -

# # 曲線座標系での接ベクトルのお話
# ## 概要
#   - 藤原著「情報幾何学の基礎」牧野書店 Chapter1 のモチベーションを Julia を介して理解しようという試みのノートです．
#

# ## 必要な Julia のライブラリを導入

using Plots
using Interact
using ForwardDiff
using SymPy
using LaTeXStrings

# # 曲線座標系の導入
#
# - $p$, $q$ をゼロ以上の実数とし, 次で定まる二次曲線 $C_p$, $C_q$ を考える.
#
# $$
# \begin{equation}
# C_p: y=px^2,\quad C_q: x=qy^2 \label{pqcurve}
# \end{equation}
# $$
#
# - これは簡単な計算によって, 原点及び $P=P(p,q)=(p^{-2/3}q^{-1/3}, p^{-1/3}q^{-2/3})$ という交点を持つ.

# ## 試しに描画してみる
#
# - $(z^1,z^2)$ を標準的な直交座標系とする二次元ユークリッド空間上に曲線 $C_p, C_q$ 及びその交点を描画してみよう．

# +
# P=P(p,q)=(φ¹(p,q),φ²(p,q))
φ¹(p::Real,q::Real) = p^(-2/3)*q^(-1/3)
φ²(p::Real,q::Real) = p^(-1/3)*q^(-2/3)

xlim=(0.,3.)
ylim=(0.,3.)

function vis_pqcurve(p, q)
    x = range(xlim...,length=30) |> collect
    y = range(ylim...,length=30) |> collect
    plot!(x,p*x.^2,color=:blue, label=L"C_p")
    plot!(q*y.^2,y,color=:orange, label=L"C_q")
    scatter!([φ¹(p,q)],[φ²(p,q)], label=L"P")
end

p₀, q₀=(1.,2.4) # 適当に定める.
plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2")
vis_pqcurve(p₀, q₀)
# -
# ### 対話的操作で確認
#
# - せっかくなので $(p,q)$ を動かした時に $C_p$, $C_q$ がどのように振舞うかを確認する.

# +
p_slider = slider(0.001:0.1:3,label="p")
q_slider = slider(0.001:0.1:3,label="q")

cross_pt = map(
    (p,q) -> begin
            plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2")
            vis_pqcurve(p,q)
    end,
    map(observe,[p_slider,q_slider])...,
)

[p_slider,q_slider,cross_pt] .|> display;
# -

# - slider を作成を動かすと曲線が変化する様子がわかる．そして対応する交点も変化する．

# ## 曲線による平面の点の特徴付け
#
# - このようにして二つの二次曲線 $C_p$, $C_q$ を与えると対応する(非自明な交点である) $P$ を与えることができる．
# - 逆に第一象限の点 $P=(z^1, z^2)$, $z^1>0,z^2>0$ が与えられたとき
#
# $$
# \begin{equation}
# p = \frac{z^2}{(z^1)^2} ,\quad q = \frac{z^1}{(z^2)^2} \label{z2p}
# \end{equation}
# $$
#
# で与えられる $p$, $q$ から定まる $C_p$, $C_q$ の交点は $P$ になっている. 実際, $(\ref{z2p})$ は $x=z^1$, $y=z^2$ とおいたときに $(\ref{pqcurve})$ の解になっている.
#
# - さらに, 相異なる $(p,q)$ の組 $(p_1,q_1)$, $(p_2,q_2)$ に対して定まる $P(p_1,q_1)$, $P(p_2,q_2)$は相異なることも容易にわかる．

# ### 確認のため描画
#
# - 下記のコードは直交座標系を持つユークリッド空間の点 $(z^1,z^2)$ を動かした時に対応する $(p,q)$ の位置を描画してくれる.
# - 描画された $(p,q)$ に対応する $C_p$, $C_q$ から定まる $P$ は $(z^1,z^2)$ と一致する

# +
xlim=(0.,3.)
ylim=(0.,3.)

z1_slider = slider(0.001:0.1:3,label="z¹")
z2_slider = slider(0.001:0.1:3,label="z²")

cross_pt = map(
    (z1,z2) -> begin
        p, q = z2/z1^2, z1/z2^2 # (2) の実装
        plt1=plot(xlim=xlim,ylim=ylim,xlabel=L"p", ylabel=L"q",aspect_ratio=:equal)
        scatter!(plt1, [p],[q], label=L"(p,q)")
        plt2=plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2",aspect_ratio=:equal)
        # ここで求めたp, q で定まる二次曲線 C_p, C_q の曲線を描画する．
        vis_pqcurve(p, q)
        # (z^1,z^2) の点を描画. vis_pqcurve 関数で描画する P と重なるはず．
        scatter!([z1],[z2],color=:red, alpha=0.4,label=L"(z^1,z^2)")
        plot(plt1,plt2,layout=(1,2))
    end,
    map(observe,[z1_slider,z2_slider])...,
)

[z1_slider,z2_slider,cross_pt] .|> display;
# -

# - ここからわかったことは二次元ユークリッド空間の第一象限の点は我々が最も馴染みがある直交座標系 $(z^1,z^2)$ による表示以外に別の表示方法が $(\ref{pqcurve})$ の意味で $p$ と $q$ の組み $(p,q)$ で与えられることがわかった. 二次元ユークリッド空間の点を $(\ref{pqcurve})$ の意味で $(p,q)$ と表示する座標系を曲線座標系と呼ぶことにする．

# # 曲線座標系での接ベクトル
#
# - 前章では曲線座標系で通常のユークリッド空間の点を捉えるアプローチを導入した．
# - ここでは曲線座標系での接ベクトルを導入する．曲線座標系での座標軸に沿った接ベクトルは直交座標系の場合と異なり, 注目している点の位置によって依存する. このことを Julia を用いて確認する．

# ## 一つのパラメータを動かしてみる
#
# - 曲線座標系 $(p,q)$ で $q$ を $q=q_0$ と適当な値 $q_0$ と固定し $p$ を変動させてみよう．これを Julia では下記の実装で行うことができる．読者はスライダーを動かすことで対話的にユークリッド空間の点 $P=P(p,q)$ がどのように変動するかを知ることができる.
# - もちろん $p$ を止めて $q$ を動かすこともできる. 適宜 `ui` 変数をいじると良い．

# +
p_slider = slider(0.001:0.1:3,label="p")
q_slider = slider(0.001:0.1:3,label="q")

function vis_interactie_pqcurve(p,q)
    plt1=plot(xlim=xlim,ylim=ylim,xlabel=L"p",ylabel=L"q",aspect_ratio=:equal)
    scatter!([p],[q],label=L"(p,q)")
    plt2=plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2",aspect_ratio=:equal)
    vis_pqcurve(p, q)
    plot(plt1,plt2,layout=(1,2))
end

cross_pt = map(
    (p,q) -> begin
        vis_interactie_pqcurve(p,q)
    end,
    map(observe,[p_slider,q_slider])...,
)

ui=[
    p_slider,
    #q_slider,
    cross_pt
]

ui .|> display;
# -

# - 上記の実装で $p$ を変動することで $P(p,q=q_0)$ がどう振舞うのかが直感的にわかった．
# - $p$ を時間パラメータ $t$ とみなすと $P(p,q=q_0)$ の振る舞いは時刻 $p=p_0=t_0$ にいた粒子 $P(t=t_0):=P(p=t_0,q=q_0)$ が微小時間 $\Delta t$ が経過した時刻 $p=t_0+\Delta t$ での粒子の位置 $P(t=t_0+\Delta t)$ の動きを記述することと解釈することができる．この意味で次の極限
#
# $$
# \begin{equation}
#     \lim_{\Delta t \to 0} \frac{P(t=t_0+\Delta t)-P(t=t_0)}{\Delta t} \label{velocity}
# \end{equation}
# $$
#
# は時刻 $t=t_0$ での速度ベクトルと思える．以下で見るように $(\ref{velocity})$ は $C_{q=q_0}$ の $P(t=t_0)$ での接ベクトルともみなせる.
#
# - 加えて，上記の実装を実行して得られた左の図は $(p,q)$ を直交座標系としてみた $(p,q)$-平面とみなせ，スライド `p_slider ` に対応するスライダーを変化させると $(p,q)$ は $p$ 軸に沿って変動することがわかる．この意味で $(\ref{velocity})$ を $p$ 軸に沿った点 $P=P(t=t_0)$ における接ベクトルと呼ぶことにしよう．というのは $(\ref{velocity})$ は
#
# $$
# \begin{equation}
#     \lim_{\Delta t \to 0} \frac{P(t=t_0+\Delta t)-P(t=t_0)}{\Delta t}
#     =
#     \lim_{\Delta p \to 0} \frac{P(p=p_0+\Delta p,q=q_0)-P(p=p_0,q=q_0)}{\Delta p}
#     \label{lim_along_p_axis}
# \end{equation}
# $$
#
# と変形でき, $(\ref{pqcurve})$ の対応によって, $(p,q)$-平面 で $(p=p_0,q=q_0)$ に向かう $p$ 軸と平行に沿って動かした方向への極限になっているからである．以上の議論から $p$ 軸に沿った点 $P=P(t=t_0)=P(p=p_0,q=q_0)$ における接ベクトルを次で定義する:
#
# $$
# \begin{equation}
# e_p(P)=\left.\frac{\partial P(p,q)}{\partial p}\right|_{(p,q)=(p_0,q_0)} \label{tl_wrt_p}
# \end{equation}
# $$
#
# しばしば $(\ref{tl_wrt_p})$ 右辺はどの点で微分をするかという情報を省いて
#
# $$
# \frac{\partial P}{\partial p}
# $$
#
# のように略記される.
#
#
# - 同様に $q$ 軸に沿った接ベクトル $e_q(P)$ も定義できる.

# ## 接ベクトルの可視化
#
# - $p$ 軸及び $q$ 軸に沿った接ベクトルを Julia で確かめよう.
# - 接ベクトルにおいて数値微分が必要にあるが [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) パッケージで賄っている.
# - 実装の便宜上 $P = P(p,q)=(\varphi^1(p,q),\varphi^2(p,q))$ とおいている．ここで
#
# $$
# \begin{equation}
# \varphi^1(p,q)=p^{-2/3}q^{-1/3},\quad \varphi^2(p,q)=p^{-1/3}q^{-2/3}
# \end{equation}
# $$
#
# である.

# +
φ¹(pq)=φ¹(pq...) # almost same as φ¹(pq[1], pq[2])
φ²(pq)=φ²(pq...)
∇φ¹=pq->ForwardDiff.gradient(φ¹,pq)
∇φ²=pq->ForwardDiff.gradient(φ²,pq)

function vis_tangent_vector(p,q)
    pq=[p,q]
    ∂P∂p = [∇φ¹(pq)[1],∇φ²(pq)[1]]
    ∂P∂q = [∇φ¹(pq)[2],∇φ²(pq)[2]]
    plot!(
        [φ¹(pq),φ¹(pq)+∂P∂p[1]],
        [φ²(pq),φ²(pq)+∂P∂p[2]],
        line=:arrow,
        label=L"e_p(P)"
    )
    plot!(
        [φ¹(pq),φ¹(pq)+∂P∂q[1]],
        [φ²(pq),φ²(pq)+∂P∂q[2]],
        line=:arrow,
        label=L"e_q(P)"
    )
end

p_slider = slider(0.001:0.05:3,label="p")
q_slider = slider(0.001:0.1:3,label="q")

function vis_interactie_pqcurve(p,q)
    plt1=plot(xlim=xlim,ylim=ylim,xlabel=L"p",ylabel=L"q",aspect_ratio=:equal)
    scatter!([p],[q],label=L"(p,q)")
    plt2=plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2",aspect_ratio=:equal)
    vis_pqcurve(p, q)
    vis_tangent_vector(p, q)
    plot(plt1,plt2,layout=(1,2))
end

cross_pt = map(
    (p,q) -> begin
        vis_interactie_pqcurve(p,q)
    end,
    map(observe,[p_slider,q_slider])...,
)

ui=[
    p_slider,
    q_slider,
    cross_pt
]

ui .|> display;
# -

# ## 念のため検算
#
# SymPyを用いて記号計算と組み合わせて確認しておこう
# SymPy.jl に優しい形で `φ¹(p::Sym,q::Sym), φ²(p::Sym,q::Sym)` を定義しておく
#

φ¹(p::Sym,q::Sym) = 1/(∛(p))^2 * 1/∛(q)
φ²(p::Sym,q::Sym) = 1/∛(p) * 1/(∛(q))^2

# ### 偏微分を計算できるようにする

∂φ¹∂p(p::Sym,q::Sym)=diff(φ¹(p,q),p)
∂φ¹∂q(p::Sym,q::Sym)=diff(φ¹(p,q),q)
∂φ²∂p(p::Sym,q::Sym)=diff(φ²(p,q),p)
∂φ²∂q(p::Sym,q::Sym)=diff(φ²(p,q),q)

# ### `p`, `q` の型が実数などの数値を取れば ForwardDiffのロジックを使って微分計算をすることになる

∂φ¹∂p(p::Real,q::Real)=∇φ¹([p,q])[1]
∂φ¹∂q(p::Real,q::Real)=∇φ¹([p,q])[2]
∂φ²∂p(p::Real,q::Real)=∇φ²([p,q])[1]
∂φ²∂q(p::Real,q::Real)=∇φ²([p,q])[2]

# Julia のマルチプルディスパッチの機能を用いて p, q の型から振る舞いを変えることができる．
∂P∂p(p,q) = [∂φ¹∂p(p,q),∂φ²∂p(p,q)]
∂P∂q(p,q) = [∂φ¹∂q(p,q),∂φ²∂q(p,q)]

# ### SymPy で計算
#
# まずは記号計算で偏微分を計算する

∂P∂p(Sym("p"),Sym("q"))

# 次にForwardDiffを呼ぶようにランダムな値を生成しその出力を見る

pq=[rand(), rand()]
∂P∂p(pq...)

# SymPyのオブジェクト p,q に具体的な値を代入して一致しているかを調べる

symp, symq = Sym("p"), Sym("q")
map(x->x(symp=>pq[1],symq=>pq[2]), ∂P∂p(symp,symq))

∂P∂q(Sym("p"),Sym("q"))

# 以下同様に $q$ についての偏微分でも確かめられる

∂P∂q(pq...)

map(x->x(p=>pq[1],q=>pq[2]), ∂P∂q(p,q))

# 良さそうである

# # ある方向に対しての接ベクトル
#
# - 我々は $p$ 軸に沿った接ベクトルを導入した, これは $p$ に関係ないパラメータを(この文脈では $q$ のこと)固定して $p$ を変化させたときの(曲線座標系が導入された) $(z^1,z^2)$ - 平面での点 $P(p,q)$ の動く方向を表していた．$q$ 軸に沿った接ベクトルも同様である．では, $p$, $q$ が同時に動くとどうなるか？その場合の $P(p,q)$ の振る舞いはどう記述できるのかをみてみよう.
# - 時刻 $t$ での パラメータ $p$, $q$ の ($p$, $q$) - 平面でみたときの位置が $(p(t),q(t))$ と書けるとしよう. もう少し形式ばると $p$, $q$ は $t$ についての滑らかな関数となってるとする. ここでの滑らかさは必要に応じて要請される回数ぐらいは微分可能であることを満たすいい感じの条件であるとしておく．
# - 我々の関心はある時刻での $t=t_0$ の微小変化をみたいのでそのある時刻を基準にする, すなわち, $t_0=0$ と仮定して良く. $p(t), q(t)$ の $t$ についての定義域は $t=t_0$ の周りの領域で定義されていれば良い．後々のために $t=0$ での $p(t), q(t)$ の導関数の値 $v_p$, $v_q$を各々次のようにおく:
#
# $$
# \begin{align}
#     v_p := \dot{p}(t=0) := \left. \frac{d p}{d t}(t) \right|_{t=0} \\
#     v_q := \dot{q}(t=0) := \left. \frac{d q}{d t}(t) \right|_{t=0} \\
# \end{align}
# $$

# +
p₀, q₀ = 3., 3. # 時刻 t=0 での初期位置
# v_q のような q の下付きが備わっていないので vp, vq という変数名を採用する
vp, vq = 3*1.5 ,3*1.25 # 右辺は適当に決めた．読者は好きなようにいじって良い

p(t::Float64) = vp*t + p₀
q(t::Float64) = vq*t + q₀
φ¹(t::Float64) = p(t)
φ²(t::Float64) = q(t)
φ(t::Float64) = (p(t),q(t))
# -

# - $t$ を変動させたときの $(p,q)$-平面での $(p(t),q(t))$ の動き，パラメータが定める $(z^1, z^2)$-平面での点 $P:=P(p,q)$ の動きをみておくことにする．

# +
t_slider = slider(-0.5:0.01:0.5,label="t")
# value を指定することで t のスライダーを変化させたときに p,qのスライダーが動くようになる．
# だたしこのアプローチだと，p, q　を動かす操作ができない
p_slider = slider(0.001:0.01:5.,value=map(a->p(a), observe(t_slider)),label="p")
q_slider = slider(0.001:0.01:5.,value=map(a->q(a), observe(t_slider)),label="q")

function vis_interactive_pqcurve(p, q)
    xlim=(0,5)
    ylim=(0,5)
    plt1 = plot(xlim=xlim, ylim=ylim, xlabel=L"p", ylabel=L"q", aspect_ratio=:equal)
    plot!(plt1, [p₀-vp,p₀+vp],[q₀-vq, q₀+vq], color=:green, label=L"(v_p,v_q)")
    scatter!(plt1, [p],[q],label=L"(p,q)")
    plt2 = plot(xlim=xlim, ylim=ylim, xlabel=L"z^1", ylabel=L"z^2", aspect_ratio=:equal)
    vis_pqcurve(p, q)
    vis_tangent_vector(p, q)
    plot(plt1, plt2,layout=(1,2))
end

ui = map(
    t ->begin
        vis_interactive_pqcurve(p(t),q(t))

    end,
    observe(t_slider)
)

vbox([
    hbox(t_slider, p_slider, q_slider),
    ui
])

# -

# - 左の緑色の線は左側の空間上に伸びる $(v_p,v_q)$ を傾きに持つ直線である．$t$ が動くと オレンジ色の $(p,q)$ がその直線状を動くことがわる．
# $P$ もなんとなく$e_p(P)$ または $e_q(P)$ の向く方向へと進んでいる気がする．
# - さて，今一度座標関数系の定義に戻ると $P=P(p,q)=(p^{-2/3}q^{-1/3}, p^{-1/3}q^{-2/3})$ であった，右辺の各々の成分は $p$, $q$ についての関数であるが，我々が今見ている文脈ではさらに $t$ についての関数だとも見ることができる．その意味で $P$ は $t$ 依存するベクトル値関数とみなせる. その意味で得られる $P(t)$ の微分は（各成分に対して）合成関数の連鎖律を用いると
#
# $$
# \begin{equation}
# \frac{d P(t)}{dt} = \dot{p}(t)\left.\frac{\partial P(p,q)}{\partial p}\right|_{(p,q)=(p(t),q(t))} +  \dot{q}(t)\left.\frac{\partial P(p,q)}{\partial q}\right|_{(p,q)=(p(t),q(t))}
# \end{equation}
# $$
#
# となる．これから特に $t=0$ の場合
#
# $$
# \begin{equation}
# \left. \frac{d P(t)}{dt}\right|_{t=0} = v_p e_p(P) +  v_q e_q(P) \label{tangent_vec_along_t}
# \end{equation}
# $$
#
# ということになる．したがって，この値 $\ref{tangent_vec_along_t}$ は $p$,$q$ 軸に沿った接ベクトルのの線形結合でかける．しかも係数は各々の軸に$v_p$,$v_q$が掛かっている．この意味で $\ref{tangent_vec_along_t}$ を $(v_p,v_q)$ 方向の接ベクトルと定義しよう．
# 特に $(v_p,v_q)=(1,0)$ の場合は $p$ 方向の接ベクトルになるが，これは $p$ 軸に沿った接ベクトル $e_p(P)$ そのものである．
#
#

# # 接空間が導入できた

# Remark: 上の実装で確認したように,
#  $(z^1,z^2)$-平面 では $(v_p, v_q)$ の方向に $P=P(p,q)$ は動いていないが $(p,q)$-平面で見ると $(v_p, v_q)$ 方向に伸びる直線上に伸びている様子がわかる．
#
# ところで $e_p(P), e_q(P)$ は $(z^1,z^2)$ の上で $P$ から伸びる一次独立なベクトルの組みになる．これは $P = P(p,q)=(\varphi^1(p,q),\varphi^2(p,q))$ ,ただし,
#
# $$
# \begin{equation}
# \varphi^1(p,q)=p^{-2/3}q^{-1/3},\quad \varphi^2(p,q)=p^{-1/3}q^{-2/3}
# \end{equation}
# $$
#
# で与えられる写像のヤコビ行列の行列式がゼロでないことからわかる．実際，ヤコビ行列は $e_p(P), e_q(P)$ を(縦ベクトルと見てそれらを横に)並べて得られるからである．
#
# 以上のことから次のことがわかる．
#
# - $(z^1,z^2)$-空間上の点 $P$ を曲線座標系 $(p,q)$ の立場から見て解析する際，
# 通常の（直交座標系が導入されている）空間と同様に$P$での接ベクトルの概念が重要になる.
# $P$での接ベクトルは互いに独立な $p$ 方向の接ベクトル $e_p(P)$ と $q$ 方向の接ベクトル $e_p(P)$ の線形結合で表されることがわかった.
# - よって $P$での接ベクトル全体は$e_p(P)$と$e_q(P)$を基底とするベクトル空間をなす.曲線座標系$(p,q)$ が入った空間上の点 $P$ での接空間と呼ぶ. $P$ での接空間を $V_P$ と表記しよう(本当は $T_pM$ のように書いた方がかっこいい)．
#
#
#
#

# - さて, 接ベクトルの可視化の実装を思い出そう. もう一度同じコードを貼り付けている．

# +
φ¹(pq)=φ¹(pq...) # almost same as φ¹(pq[1], pq[2])
φ²(pq)=φ²(pq...)
∇φ¹=pq->ForwardDiff.gradient(φ¹,pq)
∇φ²=pq->ForwardDiff.gradient(φ²,pq)

function vis_tangent_vector(p,q)
    pq=[p,q]
    ∂P∂p = [∇φ¹(pq)[1],∇φ²(pq)[1]]
    ∂P∂q = [∇φ¹(pq)[2],∇φ²(pq)[2]]
    plot!(
        [φ¹(pq),φ¹(pq)+∂P∂p[1]],
        [φ²(pq),φ²(pq)+∂P∂p[2]],
        line=:arrow,
        label=L"e_p(P)"
    )
    plot!(
        [φ¹(pq),φ¹(pq)+∂P∂q[1]],
        [φ²(pq),φ²(pq)+∂P∂q[2]],
        line=:arrow,
        label=L"e_q(P)"
    )
end

p_slider = slider(0.001:0.05:3,label="p")
q_slider = slider(0.001:0.1:3,label="q")

function vis_interactie_pqcurve(p,q)
    plt1=plot(xlim=xlim,ylim=ylim,xlabel=L"p",ylabel=L"q",aspect_ratio=:equal)
    scatter!([p],[q],label=L"(p,q)")
    plt2=plot(xlim=xlim,ylim=ylim,xlabel=L"z^1", ylabel=L"z^2",aspect_ratio=:equal)
    vis_pqcurve(p, q)
    vis_tangent_vector(p, q)
    plot(plt1,plt2,layout=(1,2))
end

cross_pt = map(
    (p,q) -> begin
        vis_interactie_pqcurve(p,q)
    end,
    map(observe,[p_slider,q_slider])...,
)

ui=[
    p_slider,
    q_slider,
    cross_pt
]

ui .|> display;
# -

# - ここで重要なのは 曲線座標系が入った $(z^1,z^2)$-空間では基底であった $e_p(P)$, $e_q(P)$ が変化する. これは互いに異なる点 $P_1$ $P_2$ において $e_p(P_1)$ $e_p(P_2)$ が一般に同じでないことを意味している．何を当たり前のことをいってるんだ？と思われるかもしれない読者はある意味で Julia 実装の世界にのめり込んでしまっている（これは良いこと）．通常の直交座標系が入った通常の空間では接ベクトルの基底は点の位置によらない標準基底がとれてしまうのであまり意識してなかったはずである．
# - Juliaの世界にのめり込んでしまった読者は接空間，接ベクトルという言葉に出会った時「それはどの空間で考えていますか？，どの点での接空間を考えてますか」と質問したくなるだろう. 互いに異なる $P_1$, $P_2$ に対応する接ベクトル $V_{P_1}$, $V_{P_2}$ から一つずつ要素 $v_{P_1}$, $v_{P_2}$ をとってきた時にこれらを直接比較したり足したり，引いたりすることに意味がないとかんじるからである．
#

# # いやでも，目で見えるじゃん，比較できそうじゃん
#
# - $P_1$ と $P_2$ が十分近かったら，近似できそうじゃん？可視化して見えるベクトルを何かしらの方法でズズズーっと平行移動のように引っ張ってきて $P_1$ に持っていってそこで比較すればいいんじゃないの？
#
# # せやな！それを正当化するのが接続とか共変微分とかの話やで！

# - 教科書にある接続係数の話を読んで計算頑張りませう！
# - 糸冬


text:
# Produce benchmarks comparing all functions with each other
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))
using BenchmarkTools, Statistics

# So that we can write a unique loop over all estimators
genentropyf(X, εs::AbstractVector) = genentropy.(Ref(X), εs)
buenooroviof(X, εs) = estimate_r0_buenoorovio(X)
takensf(X, εs) = takens_best_estimate(X, maximum(εs))

Ds = 4:10 #10 # Dimensions
Ns = 10.0 .^ (3.3333:0.3333:5.3333) .|> round .|> Int

methods = [
    "estimate r0",
    "entropy",
    "boxed large ε",
    "boxed small ε",
    "boxed small ε + r0",
    "Takens best",
]

k = 20 # amount of values in ε range
base = MathConstants.e
Ni = 6
Di = 1

function get_times_b(f, X, εs)
    bm = @benchmark f($X, $εs)
    return median(bm).time/1e9
end

function get_times_e(f, X, εs)
    bm = @elapsed f(X, εs)
    return bm
end

get_times = get_times_b

function create_times(X)
    B = zeros(length(methods))
    mi, ma = minmaxima(X)
    εmax = mean(ma - mi)
    r0, ε0 = estimate_r0_buenoorovio(X)

    B[1] = get_times(buenooroviof, X, 0)

    upper_large = log(base, εmax)
    upper_small = log(base, r0)
    lower = log(base, ε0)

    small_range = float(base) .^ range(lower+1, upper_small; length = k)
    if !issorted(small_range)
        println("small range not sorted")
        small_range = float(base) .^ range(lower+1, lower+4; length = k)
    end
    large_range = float(base) .^ range(lower+1, upper_large-2; length = k)
    if !issorted(large_range)
        println("small range not sorted")
        large_range = float(base) .^ range(lower+1, upper_large+4; length = k)
    end

    B[2] = get_times(genentropyf, X, large_range)
    B[3] = get_times(boxed_correlationsum, X, large_range)
    B[4] = get_times(boxed_correlationsum, X, small_range)
    B[5] = B[4] + B[1]
    B[6] = get_times(takensf, X, large_range)
    return B
end

data = Data.lorenz96_chaotic(; D = Ds[Di], N = maximum(Ns))
BN = zeros(length(Ns), length(methods))
for (i, N) in enumerate(Ns)
    @show i, N
    X = data[1:N]
    BN[i, :] = create_times(X)
end

BD = zeros(length(Ds), length(methods))
for (i, D) in enumerate(Ds)
    @show i, D
    X = Data.lorenz96_chaotic(; D, N = Ns[Ni])
    BD[i, :] = create_times(X)
end

# sname = savename("benchmarks", @dict(M, S, T, k), "jld2")

wsave(datadir("benchmarks", "benchmarks.jld2"), (@strdict Ds Ns methods BN BD Di Ni k base))

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
"""
    Continuous(size) <: Block

`Block` for collections of numbers. `obs` is a valid observation
if it's length is `size` and contains `Number`s.
"""
struct Continuous <: Block
    size::Int
end

function checkblock(block::Continuous, x)
    block.size == length(x) && eltype(x) <: Number
end

mockblock(block::Continuous) = rand(block.size)

function blocklossfn(outblock::Continuous, yblock::Continuous)
    outblock.size == yblock.size || error("Sizes of $outblock and $yblock differ!")
    return Flux.Losses.mse
end