text:
# ---------------------------------------------------------------------------- #
#
#   hdgSolveElas.jl
#
#   Solve convection-diffusion equations (n-dimensional)
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    hdgSolve( master::Master, mesh::Mesh, problem::CDR )

Solves convection-diffusion equations for n-dimensional problems.
"""
function hdgSolve( master::Master, mesh::Mesh, problem::CDR; method=:LU)

dim     = mesh.dim
nelem   = size( mesh.nodes, 3 ) # Mumber of elements in mesh
nnodes  = size( mesh.nodes, 1 ) # Number of nodes in one element
nodfac  = master.nodfac # Number of nodes on a face

unkUhat = nodfac                # Total of uhat unknowns per face
nfaces  = dim + 1               # Number of faces per element

# Stability parameter
τ = 5

### Initialize quantities

# Equations of motion (Newton's second law)
A = fill(0.0, nnodes*dim    , nnodes*dim)       # (v_i, q^h_i)_{T^h}
B = fill(0.0, nnodes*dim    , nnodes    )       # <v_i,i,  >_{∂T^h}
C = fill(0.0, nnodes*dim    , nodfac*nfaces)    # <v_i, func(\hat{u}^h_k) >_{∂T^h}

N = fill(0.0, nnodes        , nnodes*dim)       # <v_i, func(\hat{u}^h_k) >_{∂T^h}
D = fill(0.0, nnodes        , nnodes)           # <v_i, func(\hat{u}^h_k) >_{∂T^h}
E = fill(0.0, nnodes        , nodfac*nfaces)    # <v_i, func(\hat{u}^h_k) >_{∂T^h}

K = fill(0.0, nodfac*nfaces , nnodes*dim)       # <v_i, func(\hat{u}^h_k) >_{∂T^h}
L = fill(0.0, nodfac*nfaces , nnodes)           # <v_i, func(\hat{u}^h_k) >_{∂T^h}
M = fill(0.0, nodfac*nfaces , nodfac*nfaces)    # <v_i, func(\hat{u}^h_k) >_{∂T^h}

G = fill(0.0, nodfac*nfaces , 1)                # (v_{i}, b_{i})_{T^h}

# Matrix with unknowns for uhath
H = fill(0.0, nodfac*nfaces , nodfac*nfaces)
R = fill(0.0, nodfac*nfaces , 1            )    # (v_{i}, b_{i})_{T^h}

zm = fill( 0.0, nnodes*dim, 1)

# Matrices that need to be saved in order to recover u and qh
F    = fill(0.0, nnodes        , 1,              nelem)                # (v_{i}, b_{i})_{T^h}
ABND = fill(0.0, nnodes*(dim+1), nnodes*(dim+1), nelem)
CE   = fill(0.0, nnodes*(dim+1), nodfac*nfaces,  nelem)

# Set up column and row indices for sparse matrix
indRow = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )
indCol = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )
indUnk = fill( 0.0, (dim+1)^2 * nelem * unkUhat^2 )

Rfull = fill( 0.0, size(mesh.f,1)*unkUhat ) # RHS vector

rr = 1 # iterator for unknowns in sparsity matrix

∇ϕc = fill(0.0, size(master.∇ϕ))

κ = problem.κ
c = problem.c

# preallocate
jcwd = fill( 0.0, size(master.∇ϕ,2), size(master.∇ϕ,2) )
∂ξ∂x = fill( 0.0, size(master.∇ϕ,2), dim^2 )
∂x∂ξ = fill( 0.0, size(master.∇ϕ,2), dim^2 )
ϕjϕn = [Array{Float64}(nodfac,nodfac) for idx in 1:dim]

for pp in 1:nelem # Loop over all elements

  # zero all matrices
  A .*= 0.0
  B .*= 0.0
  C .*= 0.0
  N .*= 0.0
  D .*= 0.0
  E .*= 0.0
  K .*= 0.0
  L .*= 0.0
  M .*= 0.0
  G .*= 0.0
  H .*= 0.0
  R .*= 0.0

  # Compute Jacobians
  compJacob!( master, mesh.nodes[:,:,pp], ∂ξ∂x, jcwd, ∂x∂ξ )

  pLoc = master.ϕ' * mesh.nodes[:,:,pp]

  ∇ϕc = getderbfel( master, ∂ξ∂x )

  # mass matrices
  ϕj  = master.ϕ * jcwd
  ϕjϕ = master.ϕ * jcwd * master.ϕ'

  # ------------------------- Volume integrals ------------------------------- #
  ## A
  # (v,qₕ)_{Tₕ}
  for ii in 1:dim
    A[1+(ii-1)*nnodes:ii*nnodes,1+(ii-1)*nnodes:ii*nnodes] = ϕjϕ
  end

  ## B
  # (w_{ij}  , ϵ^h_{ij})_{T^h}
  for ii in 1:dim
    B[1+(ii-1)*nnodes:ii*nnodes,:] = ∇ϕc[:,:,ii] * ϕj'
  end

  ## N
  N += -κ * B'

  ## D (second part)
  # (w_{ij,j}, u^h_i)_{T^h} + (w_{ij,i}, u^h_j)_{T^h}
  for ii in 1:dim
    D -= c[ii] * ∇ϕc[:,:,ii] * ϕj'
  end

  ## F
  # (z_{ij}  , σ^h_{ij})_{T^h}
  src       = problem.source(pLoc)
  F[:,:,pp] = ϕj * src

  # -------------------------------------------------------------------------- #

  # ------------------------- Boundary integrals ----------------------------- #
  for qq in 1:nfaces # loop over number of faces

    (ϕdm, pdm, nod, normal, jcwdm) = compJacobFace( mesh, master, pp, qq )

    jcwddm = diagm( jcwdm )

    faceInd = (qq-1) * nodfac + ( 1:nodfac )

    cn = fill(0.0, size(normal,1))
    for ii in 1:dim
      cn += c[ii] * normal[:,ii]
    end

    # mass matrices
    ϕjf  = ϕdm * jcwddm
    ϕjϕf = ϕjf * ϕdm'
    for ii in 1:dim
        ϕjϕn[ii] = ϕjf * ( ϕdm * diagm(normal[:,ii]) )'
    end

    ## C
    # -<v_i, σ^h_{ij} n_j>_{∂T^h}
    for ii in 1:dim
      C[(ii-1)*nnodes + nod,faceInd] -= ϕjϕn[ii]
    end

    ## D
    # <v_i, τ_{ijkl} u^h_k * n_l * n_j >_{∂T^h}
    D[nod,nod] += τ * ϕjϕf

    ## E
    # <v_i, τ_{ijkl} \hat{u}^h_k * n_l * n_j >_{∂T^h}
    E[nod,faceInd] = - τ * ϕjϕf + ϕdm * ( diagm(cn) * jcwddm ) * ϕdm'

    indF = abs( mesh.t2f[pp,qq] )

    if ( mesh.f[ indF, end ] < 0 )
      # At boundary

      bNo = - mesh.f[ indF, end ]

      if problem.bctype[bNo] == 1
        # Dirichlet boundary condition

        ## M
        # (μ_{i} , \hat{u}^h_{i})_{∂Ω_D}
        M[faceInd,faceInd] = ϕjϕf

        # BC RHS
        ## G
        # (μ_{i} , \bar{u}_i)_{∂Ω_D}
        bcout = problem.bcfunc[bNo](pdm)
        G[faceInd, 1] = ϕjf * bcout

      elseif problem.bctype[bNo] == 2
        # Neumann boundary condition

        ## K
        # <μ_{i} , q_{ij}*n_j>_{∂Ω_N}
        for ii in 1:dim
          K[faceInd, (ii-1)*nnodes + nod] = ϕjϕn[ii]
        end

        # BC RHS
        ## G
        # (μ_{i} , \bar{∂u∂x})_{∂Ω_N})
        bcout = problem.bcfunc[bNo](pdm)
        G[faceInd, 1] = ϕjf * bcout
      else
          error("hdgSolveElas:: BC type not recognized. 1 = Dirichlet, 2 = Neumann")
      end

    else
      # In interior

      ## K
      # <μ_{i} , σ^h_{ij}*n_j>_{∂T^h\∂Ω_D}
      for ii in 1:dim
        K[faceInd, (ii-1)*nnodes + nod] = - ϕjϕn[ii]
      end

      ## L
      # <μ_{i} , -τ_{ijkl} u^h_k n_l n_j ) >_{∂T^h\∂Ω_D}

      L[faceInd,nod] = τ * ϕjϕf

      ## M
      # <μ_{i} , τ_{ijkl} \hat{u}^h_k n_l n_j ) >_{∂T^h\∂Ω_D}
      M[faceInd,faceInd] = - τ * ϕjϕf + ϕdm * ( diagm(cn) * jcwddm ) * ϕdm'

    end # boundary if-statement

  end # end loop over element faces
  # -------------------------------------------------------------------------- #

  # ------------------------- Form elemental quantities ---------------------- #
  ABND[:,:,pp] = [ A B;
                   N D ]

  CE[:,:,pp] = [ C; E ]

  H += M - [K  L] * ( ABND[:,:,pp] \ CE[:,:,pp] )
  R += G - [K  L] * ( ABND[:,:,pp] \ [ zm; F[:,:,pp] ] )

  # -------------------------------------------------------------------------- #

  # ------------------------- Fill up complete H and R matrices -------------- #
  for qq = 1:nfaces
    indJ = [i for i=1:nfaces]
    deleteat!(indJ,qq)

    indFall = abs.(mesh.t2f[pp,1:nfaces])
    deleteat!(indFall,qq)

    indF = abs(mesh.t2f[pp,qq])
    Rfull[1+(indF-1)*unkUhat:indF*unkUhat] += R[1+(qq-1)*unkUhat:qq*unkUhat,1]

    indRow[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, unkUhat )
    indCol[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, 1, unkUhat )'[:]
    indUnk[rr:rr-1+unkUhat^2] = H[ 1+(qq-1)*unkUhat:qq*unkUhat, 1+(qq-1)*unkUhat:qq*unkUhat ][:]

    rr = rr + unkUhat^2

    for tt = 1:(nfaces-1)
      indRow[rr:rr-1+unkUhat^2] = repmat(  1+(indF-1)*unkUhat:indF*unkUhat, unkUhat )
      indCol[rr:rr-1+unkUhat^2] = repmat(  1+(indFall[tt]-1)*unkUhat:indFall[tt]*unkUhat, 1, unkUhat )'[:]
      indUnk[rr:rr-1+unkUhat^2] = H[ 1+(qq-1)*unkUhat:qq*unkUhat, 1+(indJ[tt]-1)*unkUhat:indJ[tt]*unkUhat ][:]

      rr = rr + unkUhat^2
    end
  end
  # -------------------------------------------------------------------------- #
end # end element loop

### Compute approximate trace
Hfull = sparse( indRow, indCol, indUnk, size(mesh.f,1) * unkUhat, size(mesh.f,1) * unkUhat )

#fact = ILU.crout_ilu( Hfull, τ = 0.1 )
uhath = [0.]
if method == :LU
    uhath = Hfull \ Rfull
elseif method == :GMRES
    uhath = IterativeSolvers.gmres( Hfull, Rfull, Pl=fact )
end

# ---------------------------------------------------------------------------- #

## Compute approximate scalar value and flux
uhathTri = fill( 0.0, unkUhat*nfaces, 1 )
uh       = fill( 0.0, nnodes,         1,     nelem )
qh       = fill( 0.0, nnodes,         dim,   nelem )

for pp in 1:nelem
    # ----------- Find uhath corresponding to this element ------------------- #
    indF = abs.( mesh.t2f[pp,1:nfaces] )
    for qq in 1:nfaces
      uhathTri[ 1+(qq-1)*unkUhat:qq*unkUhat, 1 ] = uhath[ 1+(indF[qq]-1)*unkUhat:indF[qq]*unkUhat ]
    end
    # ------------------------------------------------------------------------ #
    # ----------- Compute approximate displacement value --------------------- #

    rhsTemp = [ zm; F[:,:,pp] ] - CE[:,:,pp] * uhathTri

    uTemp   = ABND[:,:,pp] \ rhsTemp

    for ii in 1:dim
      qh[:,ii,pp] = uTemp[ 1+(ii-1)*nnodes:ii*nnodes ]
    end
    uh[:,1,pp] = uTemp[ 1+dim*nnodes:(dim+1)*nnodes ]
end

return ( uhath, uh, qh )

end # end hdgSolveCD


text:
import Ecto.MirrorTypes: join_with
using Base.Test

@test [:a, ' ', :b] == join_with([:a,:b], ' ')

ext: jl


text:
module Liu

using ModelingToolkit
using ..NeuronBuilder

import ..get_parameters, ..get_states, ..default_params, ..default_states

include("channels.jl")
include("calc_dynamics.jl")

export get_parameters, get_states, default_params, default_states

end

"""
Deliberately not exporting ion channel names so they don't conflict.

eg Na can only be accessed as Liu.Na from NeuronBuilder

Deliberately exporting other structs that act on components. because it is one struct, not defined separaetly in different modules, that acts on module-specific structs
"""


text:
# @test_throw, with additional testing for the exception message
macro test_throws_message(f, typ, ex...)
    quote
        msg = ""
        @test_throws $(esc(typ)) try
            $(esc(ex...))
        catch err
            msg = sprint(showerror, err)
            rethrow()
        end

        if !$(esc(f))(msg)
            # @test should return its result, but doesn't
            errmsg = "Failed to validate error message\n" * msg
            @error errmsg
        end
        @test $(esc(f))(msg)
    end
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::T) where T <: Union{Int32,UInt32} =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", T, Tuple{T}, i)
@inline sink(i::T) where T <: Union{Int64,UInt64} =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", T, Tuple{T}, i)

# the GPU runtime library
module TestRuntime
    # dummy methods
    signal_exception() = return
    # HACK: if malloc returns 0 or traps, all calling functions (like jl_box_*)
    #       get reduced to a trap, which really messes with our test suite.
    malloc(sz) = Ptr{Cvoid}(Int(0xDEADBEEF))
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end

struct TestCompilerParams <: AbstractCompilerParams end
GPUCompiler.runtime_module(::CompilerJob{<:Any,TestCompilerParams}) = TestRuntime


text:
@testset "Dinic" begin
    # Construct DiGraph
    flow_graph = DiGraph(8)

    # Load custom dataset
    flow_edges = [
        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
    ]

    capacity_matrix = zeros(Int, nv(flow_graph), nv(flow_graph))

    for e in flow_edges
        u, v, f = e
        add_edge!(flow_graph, u, v)
        capacity_matrix[u, v] = f
    end

    # Construct the residual graph
    for fg in (flow_graph, DiGraph{UInt8}(flow_graph), DiGraph{Int16}(flow_graph))
      residual_graph = @inferred(LightGraphs.residual(fg))

      # Test with default distances
      @test @inferred(LightGraphs.dinic_impl(residual_graph, 1, 8, LightGraphs.DefaultCapacity(residual_graph)))[1] == 3

      # Test with capacity matrix
      @test @inferred(LightGraphs.dinic_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28

      # Test on disconnected graphs
      function test_blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)
          #disconnect source
          h = copy(residual_graph)
          for dst in collect(neighbors(residual_graph, source))
              rem_edge!(h, source, dst)
          end
          @test @inferred(LightGraphs.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0

          #disconnect target and add unreachable vertex
          h = copy(residual_graph)
          for src in collect(in_neighbors(residual_graph, target))
              rem_edge!(h, src, target)
          end
          @test @inferred(LightGraphs.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0

          # unreachable vertex (covers the case where a vertex isn't reachable from the source)
          h = copy(residual_graph)
          add_vertex!(h)
          add_edge!(h, nv(residual_graph) + 1, target)
          capacity_matrix_ = vcat(hcat(capacity_matrix, zeros(Int, nv(residual_graph))), zeros(Int, 1, nv(residual_graph) + 1))
          flow_graph_  = vcat(hcat(flow_matrix, zeros(Int, nv(residual_graph))), zeros(Int, 1, nv(residual_graph) + 1))

          @test @inferred(LightGraphs.blocking_flow(h, source, target, capacity_matrix_, flow_graph_)) > 0

          #test with connected graph
          @test @inferred(LightGraphs.blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)) > 0
      end

      flow_matrix = zeros(Int, nv(residual_graph), nv(residual_graph))
      test_blocking_flow(residual_graph, 1, 8, capacity_matrix, flow_matrix)
    end
end


text:
"""
    embedding(small::AbstractMesh, big::AbstractMesh) -> S::AbstractSparseArray
"""
function embedding(small::AbstractMesh, big::AbstractMesh)

    @assert dimension(small) == 1
    @assert dimension(big)   == 1

    T = coordtype(small)
    tol = sqrt(eps(T))

    # store the edges of the big mesh in an octree
    ctrs_big = [cartesian(center(chart(big,v))) for v in cells(big)]
    tree_big = Octree(ctrs_big, zeros(size(ctrs_big)))

    rows = 1:numcells(small)
    cols = zeros(numcells(small))
    sgns = zeros(numcells(small))

    cells_big = collect(cells(big))
    for (i,c) in enumerate(cells(small))
        # V = vertices(small,c)
        V = chart(small, c).vertices
        ch = chart(small,c)
        x = cartesian(center(ch))
        pred(c,s) = fitsinbox(x, 0.0, c, s+tol)
        found = false
        for box in boxes(tree_big, pred)
            for j in box
                c′ = cells_big[j]
                # V′ = vertices(big,c′)
                V′ = chart(big,c′).vertices
                ch′ = chart(big, c′)
                if norm(V[1]-V′[1]) + norm(V[2]-V′[2]) < tol
                    cols[i] = j
                    sgns[i] = +1
                    found = true
                    break
                elseif norm(V[1]-V′[2]) + norm(V[2]-V′[1]) < tol
                    cols[i] = j
                    sgns[i] = -1
                    found = true
                    break
                end
            end
            found && break
        end
    end

    @assert !any(cols .== 0)
    @assert !any(sgns .== 0)
    sparse(rows, cols, sgns, numcells(small), numcells(big))
end


text:
module SPSSTools

using Tables

export value_labels

"""
    value_labels(data, value::Symbol, label::Symbol; varname = "varname", out = nothing, charq = :auto)

VALUE LABELS command with value-label pairs from data table.

* `charq` - character quoting, `:auto` or Tuple{Bool, Bool}

Example:

 `value_labels(data, :value, :label; varname = "name",  charq = (false, true))`

Pattern:

VALUE LABELS
varname
val1 'label1'
val2 'label2'
val3 'label3'.
EXECUTE.
"""
function value_labels(data, value::Symbol, label::Symbol; varname = "varname", out = nothing, charq = :auto)
    cols = Tables.columns(data)
    res  = Tuple(Tables.getcolumn(cols, y) for y in (value, label))
    str = "VALUE LABELS \n$(varname)"
    if isa(charq, Symbol) && charq == :auto
        if eltype(res[1]) <: Number
            valstr = ""
        else
            valstr = "'"
        end
        if eltype(res[2]) <: Number
            labstr = ""
        else
            labstr = "'"
        end
    else
        if charq[1] == true valstr = "'" else valstr = "" end
        if charq[2] == true labstr = "'" else labstr = "" end
    end
    for i = 1:length(res[1])
        str = str * "\n$valstr$(res[1][i])$valstr $labstr$(res[2][i])$labstr"
    end
    str = str * ".\nEXECUTE."
    str
end

end # module


text:
using TensorToolbox
using LinearAlgebra
using Arpack
#using TensorDecompositions

"""
hyper parameters max_iter and tol are decieded by sklearn
"""
function lranmf_mu(X, r ; max_iter=200, tol = 1.0E-3, verbose = false)
    # lranmf_mu is effective when r << min(n,m)
    # input is a onnegative matrix
    # proposed by Guoxu Zhou in 2012
    # https://ieeexplore.ieee.org/document/6166354

    n, m = size(X)
    epsilon = 0.000001

    # See step1 in section II-B in the paper
    if r == min(n, m)
        svd_X = svd(X)
    else
        svd_X = svds(X; nsv=r, ritzvec=true)[1]
    end
    Achil = svd_X.U * diagm(svd_X.S)
    Bchil = svd_X.V

    # See step2 in section II-B in the paper
    A = rand(n, r)
    B = rand(m, r)
    cost_at_init = norm(X - A*B')
    previous_cost = cost_at_init
    for iter = 1:max_iter
        B .= B .* ( max.( Bchil*(Achil' * A), epsilon ) ) ./ ( B*(A'*A) )
        A .= A .* ( max.( Achil*(Bchil' * B), epsilon ) ) ./ ( A*(B'*B) )

        if tol > 0 && iter % 10 == 0
            cost = norm(X - A*B')
            if verbose
                println("iter: $iter cost: $cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end

    # A * B' is rank-r matrix
    return A, B'
end

function lraSNTD(Y, reqrank)
    # Sequential nonnegative Tucker based on lraNMF
    # input is a nonnegative tensor
    # proposed by Guoxu Zhou in 2012
    # https://ieeexplore.ieee.org/document/6166354

    N = ndims(Y)
    input_tensor_shape = size(Y)

    # See section IV-2) in the paper
    # get non-negative factors
    A = []
    for n=1:N
        # Yn is Matricization of tensor Y by mode n
        Yn = tenmat(Y, n)
        An, _ = lranmf_mu(Yn, reqrank[n])
        push!(A, An)
        # ttm(Y,An,n) is mode-n product of tensor X and matrix An
        Y = ttm(Y, pinv(An), n)
    end

    # reproduce
    Y = ttm(Y, A[1], 1)
    for n=2:N
        Y = ttm(Y, A[n], n)
    end

    return Y
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
"""
    ProjectedPointDistribution(M::Manifold, d, proj!, p)

Generates a random point in ambient space of `M` and projects it to `M`
using function `proj!`. Generated arrays are of type `TResult`, which can be
specified by providing the `p` argument.
"""
struct ProjectedPointDistribution{TResult,TM<:Manifold,TD<:Distribution,TProj} <:
       MPointDistribution{TM}
    manifold::TM
    distribution::TD
    proj!::TProj
end

function ProjectedPointDistribution(
    M::Manifold,
    d::Distribution,
    proj!,
    ::TResult,
) where {TResult}
    return ProjectedPointDistribution{TResult,typeof(M),typeof(d),typeof(proj!)}(
        M,
        d,
        proj!,
    )
end

"""
    projected_distribution(M::Manifold, d, [p=rand(d)])

Wrap the standard distribution `d` into a manifold-valued distribution. Generated
points will be of similar type to `p`. By default, the type is not changed.
"""
function projected_distribution(M::Manifold, d, p=rand(d))
    return ProjectedPointDistribution(M, d, project!, p)
end

function Random.rand(
    rng::AbstractRNG,
    d::ProjectedPointDistribution{TResult},
) where {TResult}
    p = convert(TResult, rand(rng, d.distribution))
    return d.proj!(d.manifold, p, p)
end

function Distributions._rand!(
    rng::AbstractRNG,
    d::ProjectedPointDistribution,
    p::AbstractArray{<:Number},
)
    Distributions._rand!(rng, d.distribution, p)
    return d.proj!(d.manifold, p, p)
end

Distributions.support(d::ProjectedPointDistribution) = MPointSupport(d.manifold)

"""
    ProjectedFVectorDistribution(type::VectorBundleFibers, p, d, project!)

Generates a random vector from ambient space of manifold `type.manifold`
at point `p` and projects it to vector space of type `type` using function
`project!`, see [`project`](@ref) for documentation.
Generated arrays are of type `TResult`.
"""
struct ProjectedFVectorDistribution{
    TResult,
    TSpace<:VectorBundleFibers,
    ManifoldPoint,
    TD<:Distribution,
    TProj,
} <: FVectorDistribution{TSpace,ManifoldPoint}
    type::TSpace
    point::ManifoldPoint
    distribution::TD
    project!::TProj
end

function ProjectedFVectorDistribution(
    type::VectorBundleFibers,
    p,
    d::Distribution,
    project!,
    ::TResult,
) where {TResult}
    return ProjectedFVectorDistribution{
        TResult,
        typeof(type),
        typeof(p),
        typeof(d),
        typeof(project!),
    }(
        type,
        p,
        d,
        project!,
    )
end

function Random.rand(
    rng::AbstractRNG,
    d::ProjectedFVectorDistribution{TResult},
) where {TResult}
    X = convert(TResult, reshape(rand(rng, d.distribution), size(d.point)))
    return d.project!(d.type, X, d.point, X)
end

function Distributions._rand!(
    rng::AbstractRNG,
    d::ProjectedFVectorDistribution,
    X::AbstractArray{<:Number},
)
    # calling _rand!(rng, d.d, v) doesn't work for all arrays types
    return copyto!(X, rand(rng, d))
end

"""
    normal_tvector_distribution(M::Euclidean, p, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `p`.
"""
function normal_tvector_distribution(M::Manifold, p, σ)
    d = Distributions.MvNormal(zero(vec(p)), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(M), p, d, project!, p)
end

function Distributions.support(tvd::ProjectedFVectorDistribution)
    return FVectorSupport(tvd.type, tvd.point)
end


text:
export DefiniteIntegral,DefiniteLineIntegral

abstract type CalculusFunctional{S,T} <: Operator{T} end

@functional CalculusFunctional

##TODO: Add ConcreteOp

macro calculus_functional(Op)
    ConcOp=Meta.parse("Concrete"*string(Op))
    WrappOp=Meta.parse(string(Op)*"Wrapper")
    return esc(quote
        abstract type $Op{SSS,TTT} <: CalculusFunctional{SSS,TTT} end
        struct $ConcOp{S,T} <: $Op{S,T}
            domainspace::S
        end
        struct $WrappOp{BT<:Operator,S<:Space,T} <: $Op{S,T}
            op::BT
        end

        @wrapper $WrappOp

        # We expect the operator to be real/complex if the basis is real/complex
        $ConcOp(dsp::Space) = $ConcOp{typeof(dsp),prectype(dsp)}(dsp)

        $Op() = $Op(UnsetSpace())
        $Op(dsp) = $ConcOp(dsp)
        $Op(d::Domain) = $Op(Space(d))

        ApproxFunBase.promotedomainspace(::$Op,sp::Space) = $Op(sp)

        Base.convert(::Type{Operator{T}},Σ::$ConcOp) where {T} =
            (T==eltype(Σ) ? Σ : $ConcOp{typeof(Σ.domainspace),T}(Σ.domainspace))::Operator{T}

        ApproxFunBase.domain(Σ::$ConcOp) = domain(Σ.domainspace)
        ApproxFunBase.domainspace(Σ::$ConcOp) = Σ.domainspace

        Base.getindex(::$ConcOp{UnsetSpace},kr::AbstractRange) =
            error("Spaces cannot be inferred for operator")

        $WrappOp(op::Operator) =
            $WrappOp{typeof(op),typeof(domainspace(op)),eltype(op)}(op)

        Base.convert(::Type{Operator{T}},Σ::$WrappOp) where {T} =
            (T==eltype(Σ) ? Σ : $WrappOp(convert(Operator{T},Σ.op)))::Operator{T}
    end)
end

@calculus_functional(DefiniteIntegral)
@calculus_functional(DefiniteLineIntegral)

#default implementation

DefiniteIntegral(sp::UnsetSpace) = ConcreteDefiniteIntegral(sp)
DefiniteLineIntegral(sp::UnsetSpace) = ConcreteDefiniteLineIntegral(sp)

function DefiniteIntegral(sp::Space)
    if typeof(canonicaldomain(sp)) == typeof(domain(sp))
        # try using `Integral`
        Q = Integral(sp)
        rsp = rangespace(Q)
        DefiniteIntegralWrapper((Evaluation(rsp,rightendpoint)-Evaluation(rsp,leftendpoint))*Q)
    else
        # try mapping to canonical domain
        M = Multiplication(fromcanonicalD(sp),setcanonicaldomain(sp))
        Op = DefiniteIntegral(rangespace(M))*M
        DefiniteIntegralWrapper(SpaceOperator(Op,sp,rangespace(Op)))
    end
end

function DefiniteLineIntegral(sp::Space)
    if typeof(canonicaldomain(sp)) == typeof(domain(sp))
        error("Override DefiniteLineIntegral for $sp")
    end

    M = Multiplication(abs(fromcanonicalD(sp)),setcanonicaldomain(sp))
    Op = DefiniteLineIntegral(rangespace(M))*M
    DefiniteLineIntegralWrapper(SpaceOperator(Op,sp,rangespace(Op)))
end

#TODO: Remove SPECIALOPS reimplement
# *{T,D<:DefiniteIntegral,M<:Multiplication}(A::TimesFunctional{T,D,M},b::Fun) = bilinearform(A.op.f,b)
# *{T,D<:DefiniteLineIntegral,M<:Multiplication}(A::TimesFunctional{T,D,M},b::Fun) = linebilinearform(A.op.f,b)
# *{T,D<:Union{DefiniteIntegral,DefiniteLineIntegral},
#   M<:Multiplication,V}(A::FunctionalOperator{TimesFunctional{T,D,M},V},b::Fun) =
#     Fun(A.op*b)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

x0 = [0.0, 0.0]
widths = [1.0, 1.0]
nelements = [2, 2]
polyorder = 2
numqp = required_quadrature_order(polyorder) + 2

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)
quad = tensor_product_quadrature(2,4)
refpoints = interpolation_points(elasticitybasis)
dim,nf = size(interpolation_points(levelsetbasis))

dgmesh = CutCellDG.DGMesh(x0,widths,nelements,refpoints)
CutCellDG.make_vertical_periodic!(dgmesh)

testcellconnectivity = [2  1  4  3
                        3  4  0  0
                        2  1  4  3
                        0  0  1  2]
@test allequal(testcellconnectivity,CutCellDG.cell_connectivity(dgmesh,:,:))


text:
"""
@fastfunc functionname polynomial

Generates a function for fast numerical evaluation. The number of input
arguments of the generated function will match the number of variables in the
input polynomial.

Note that the polynomial is transformed using `@fastmath`, which calls
functions that may violate strict IEEE semantics.

If you only plan a few evaluations of polynomial use the standard method (see
e.g. TypedPolynomials for details). Only for a large number of evaluations the
one can break even with the function generated by this method due to the just
in time compilation overhead.

# Examples
```
julia> using SphericalHarmonics

julia> @polyvar x y z
(x, y, z)

julia> p = 15.0*x*y^2+7.5*x*z^13
7.5xz¹³ + 15.0xy²

julia> foo = @fastfunc p

julia> foo(1.0,2.0,3.0)
1.19574825e7
```

# Non-global scope

Usage from within local scope requires special care to avoid issue #4. So
instead of `foo(1.0,2.0,3.0)` we need `Base.invokelatest(foo, 1.0,2.0,3.0)`.

```
julia> using SphericalHarmonics

julia> @polyvar x y z
(x, y, z)

julia> function useInsideFunctionScope()
           p = 15.0*x*y^2+7.5*x*z^13
           foo = @fastfunc p
           Base.invokelatest(foo, 1.0,2.0,3.0)
       end
useInsideFunctionScope (generic function with 1 method)

julia> useInsideFunctionScope()
1.19574825e7
```
"""
macro fastfunc(polynomial)
	return quote
		local polystr = string($(esc(polynomial)))
		local vars = string(tuple(variables($(esc(polynomial)))...))
        # create expression for function definition
		eval(Meta.parse(vars*" -> @fastmath "*polystr))
	end
end

"""
fastfunc(polynomial::Polynomial)

Generates a function for fast numerical evaluation. The number of input
arguments of the generated function will match the number of variables in the
input polynomial.

Note that the polynomial is transformed using `@fastmath`, which calls
functions that may violate strict IEEE semantics.

If you only plan a few evaluations of polynomial use the standard method (see
e.g. TypedPolynomials for details). Only for a large number of evaluations the
one can break even with the function generated by this method due to the just
in time compilation overhead.

# Examples

```
julia> using SphericalHarmonics

julia> @polyvar x y z
(x, y, z)

julia> p = 15.0*x*y^2+7.5*x*z^13
7.5xz¹³ + 15.0xy²

julia> foo = fastfunc(p)

julia> foo(1.0,2.0,3.0)
1.19574825e7
```
"""
function fastfunc(polynomial)
    expr = Meta.parse(string(variables(polynomial))*" -> "*string(polynomial))
    # use fastmath
    expr = Base.FastMath.make_fastmath(expr)
    return mk_function(expr)
end


text:
@inline now_localzone() = TimeZones.now(TimeZones.localzone())

function utc_to_string(zdt::TimeZones.ZonedDateTime)
    zdt_as_utc = TimeZones.astimezone(zdt, TimeZones.tz"UTC")
    return Dates.format(zdt_as_utc, "yyyy-mm-dd-HH-MM-SS-sss")
end


text:
#adapt FFTViews code
using Base: tail, unsafe_length, @propagate_inbounds
# using FFTW

# A custom rangetype that will be used for indices and never throws a
# boundserror because the domain is actually periodic.
using CustomUnitRanges
include(CustomUnitRanges.filename_for_urange)

@static if isdefined(Base, :IdentityUnitRange)
    const indextypes = (URange, Base.IdentityUnitRange{<:URange})
    const CircVRange{T} = Union{URange{T}, Base.Slice{URange{T}}, Base.IdentityUnitRange{URange{T}}}
    indrange(i) = Base.IdentityUnitRange(URange(first(i), last(i)))
else
    const indextypes = (URange, Base.Slice{<:URange})
    const CircVRange{T} = Union{URange{T}, Base.Slice{URange{T}}}
    indrange(i) = Base.Slice(URange(first(i), last(i)))
end

for T in indextypes
    @eval begin
        Base.checkindex(::Type{Bool}, ::$T, ::Base.Slice) = true
        Base.checkindex(::Type{Bool}, ::$T, ::Base.LogicalIndex) = true
        Base.checkindex(::Type{Bool}, ::$T, ::Real) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractRange) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractVector{Bool}) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractArray{Bool}) = true
        Base.checkindex(::Type{Bool}, ::$T, ::AbstractArray) = true
    end
    if isdefined(Base, :IdentityUnitRange)
        @eval Base.checkindex(::Type{Bool}, ::$T, ::Base.IdentityUnitRange) = true
    end
end

export CircView

abstract type AbstractCircView{T,N} <: AbstractArray{T,N} end

struct CircView{T,N,A<:AbstractArray} <: AbstractCircView{T,N}
    parent::A

    function CircView{T,N,A}(parent::A) where {T,N,A}
        new{T,N,A}(parent)
    end
end

CircView(parent::AbstractArray{T,N}) where {T,N} = CircView{T,N,typeof(parent)}(parent)
CircView{T,N}(dims::Dims{N}) where {T,N} = CircView(Array{T,N}(undef, dims))
CircView{T}(dims::Dims{N}) where {T,N} = CircView(Array{T,N}(undef, dims))

# Note: there are no bounds checks because it's all periodic
@inline @propagate_inbounds function Base.getindex(F::CircView{T,N}, I::Vararg{Int,N}) where {T,N}
    P = parent(F)
    @inbounds ret = P[reindex(CircView, axes(P), I)...]
    ret
end

@inline @propagate_inbounds function Base.setindex!(F::CircView{T,N}, val, I::Vararg{Int,N}) where {T,N}
    P = parent(F)
    @inbounds P[reindex(CircView, axes(P), I)...] = val
end

Base.parent(F::AbstractCircView) = F.parent
Base.axes(F::AbstractCircView) = map(indrange, axes(parent(F)))
Base.size(F::AbstractCircView) = size(parent(F))

function Base.similar(A::AbstractArray, T::Type, shape::Tuple{CircVRange,Vararg{CircVRange}})
    all(x->first(x)==1, shape) || throw(BoundsError("cannot allocate CircView with the first element of the range non-zero"))
    CircView(similar(A, T, map(length, shape)))
end

function Base.similar(f::Union{Function,Type}, shape::Tuple{CircVRange,Vararg{CircVRange}})
    all(x->first(x)==1, shape) || throw(BoundsError("cannot allocate CircView with the first element of the range non-zero"))
    CircView(similar(f, map(length, shape)))
end

Base.reshape(F::CircView{_,N}, ::Type{Val{N}}) where {_,N}   = F
Base.reshape(F::CircView{_,M}, ::Type{Val{N}}) where {_,M,N} = CircView(reshape(parent(F), Val(N)))

# FFTW.fft(F::CircView; kwargs...) = fft(parent(F); kwargs...)
# FFTW.rfft(F::CircView; kwargs...) = rfft(parent(F); kwargs...)
# FFTW.fft(F::CircView, dims; kwargs...) = fft(parent(F), dims; kwargs...)
# FFTW.rfft(F::CircView, dims; kwargs...) = rfft(parent(F), dims; kwargs...)

@inline reindex(::Type{V}, inds, I) where {V} = (_reindex(V, inds[1], I[1]), reindex(V, tail(inds), tail(I))...)
reindex(::Type{V}, ::Tuple{}, ::Tuple{}) where {V} = ()
_reindex(::Type{CircView}, ind, i) = modrange(i, ind)

modrange(i, rng::AbstractUnitRange) = mod(i-first(rng), unsafe_length(rng))+first(rng)

## test
using BenchmarkTools
x = LinRange(-5,5,1000)

@btime x[200:800]

xc = CircView(x)
@btime xc[-800:-200]

## test FFTViews separately
using FFTViews
x = LinRange(-5,5,1000)

@btime x[200:800]

xc = FFTView(x)
@btime @view xc[200:800]


text:
# Fourier demo with gtk-server
# #
# Note: requires package Gaston for plotting, or adjust the code to work with your
# plotting front-end.

using GtkServer, Gaston, FFTW

# plot
function doplots(t, A, f₀, ϕ, fs)
    s = A.*sin.(2π*f₀*t .+ ϕ)
    S = fftshift(fft(s))./length(s)
    F = fftshift(fftfreq(length(s), fs))

    p1 = plot(t, s, handle = 1,
              Axes(title = "'Sinusoidal in time domain'", xlabel = "'Time (s)'",
                   ylabel = "'Amplitude'", yrange = (-5.5, 5.5), ytics = -6:2:6))
    p2 = plot(F, abs.(S), handle = 2,
              Axes(title = "'Magnitude spectrum'", xlabel = "'Frequency (Hz)'",
                   yrange = (0, 3)))
    p3 = plot(F, angle.(S), handle = 3,
              Axes(title = "'Phase spectrum'", xlabel = "'Frequency (Hz)'",
                   yrange = (-3, 3)))
    P = plot([p1, p2, p3], handle = 4)
    display(P)
end

function demo()
    # Set initial values for amplitude, frequency and phase, and display
    # the sine function.
    fs = 100
    t = 0.0:1.0/fs:1.0
    A = 1.0
    f₀ = 1.0
    ϕ = 0.0

    # Do initial plot with default values
    doplots(t, A, f₀, ϕ, fs)

    # Main program
    # This program interactively plots a sine function in time, and its
    # magnitude and phase spectra. The amplitude, frequency and phase of
    # sine function are selected through three sliders.
    gtk_init()
    # 'win' is a handle (ID number) to the main GUI window.
    gtk_send("gtk_init NULL NULL")
    win = gtk_send("gtk_window_new 0")
    # Set window's title, size and position.
    # Note that the message sent to gtk-server is always a string.
    gtk_send("gtk_window_set_title $win 'Fourier Series Demo'")
    gtk_send("gtk_window_set_default_size $win 300 200")
    gtk_send("gtk_window_set_position $win 1")
    # We'll set our widgets in a table, one widget per row-column combination.
    tbl = gtk_send("gtk_table_new 8 10 0")
    gtk_send("gtk_container_add $win $tbl")
    # Create three labels...
    lblA = gtk_send("gtk_label_new Amplitude:")
    lblF = gtk_send("gtk_label_new Frequency:")
    lblP = gtk_send("gtk_label_new Phase:")
    # and put them in specific position on our table.
    gtk_send("gtk_table_attach_defaults $tbl $lblA 1 3 1 2")
    gtk_send("gtk_table_attach_defaults $tbl $lblF 1 3 3 4")
    gtk_send("gtk_table_attach_defaults $tbl $lblP 1 3 5 6")
    # horizontal separator
    hsep = gtk_send("gtk_hseparator_new")
    gtk_send("gtk_table_attach_defaults $tbl $hsep 1 9 6 7")
    # Some Julia advertisement
    adv = gtk_send("gtk_label_new 'Created with Julia!'")
    gtk_send("gtk_table_attach_defaults $tbl $adv 1 6 7 8")
    # Create (and put on the table) the quit button.
    but = gtk_send("gtk_button_new_from_stock 'gtk-quit'")
    gtk_send("gtk_table_attach_defaults $tbl $but 7 9 7 8")
    # Create three sliders...
    sliA = gtk_send("gtk_hscale_new_with_range 1 5 0.5")
    sliF = gtk_send("gtk_hscale_new_with_range 1 40 1")
    sliP = gtk_send("gtk_hscale_new_with_range -1.6 1.6 0.1")
    # and put them on the table.
    gtk_send("gtk_table_attach_defaults $tbl $sliA 4 9 1 2")
    gtk_send("gtk_table_attach_defaults $tbl $sliF 4 9 3 4")
    gtk_send("gtk_table_attach_defaults $tbl $sliP 4 9 5 6")
    # Now we need to actually display all the widgets we've created.
    gtk_send("gtk_widget_show_all $win")

    # endless loop; we only break out of it when the quit button is clicked
    while true
        # Up to this point, we have sent commands to gtk-server in order
        # to set up the GUI and display it. However, once the GUI is
        # created, gtk-server will send messages (events) to us whenever
        # the user clicks something on the GUI. In order to read these
        # events, we send gtk-server the message "gtk_server_callback wait".
        # So, if the user has clicked on a widget on the GUI, then
        # 'event' contains the widget's handle (ID number).
        event = gtk_send("gtk_server_callback wait")
        # The user clicked on the windows 'close' button.
        if event == win
            break
        end
        # The user clicked on the quit button.
        if event == but
            break
        end
        # The user clicked on the 'amplitude' slider. Now we need
        # to go to gtk-server again to ask for the new slider's value.
        if event == sliA
            # here we get the new value
            A_g = gtk_send("gtk_range_get_value $sliA")
            A = parse(Float32, A_g)
            # plot again with the new sine values
            P = doplots(t, A, f₀, ϕ, fs)
        end
        # Frequency slider.
        if event == sliF
            f₀_g = gtk_send("gtk_range_get_value $sliF")
            f₀ = parse(Float32, f₀_g)
            P = doplots(t, A, f₀, ϕ, fs)
        end
        # ϕ slider.
        if event == sliP
            ϕ_g = gtk_send("gtk_range_get_value $sliP")
            ϕ = parse(Float32, ϕ_g)
            P = doplots(t, A, f₀, ϕ, fs)
        end
    end

    # Clean up -- close all figures and tell gtk-server to quit.
    closeall()
    gtk_exit()
end


text:
module QualityControl

using Statistics, SparseArrays

export selectg,
       scstats,
       emptyrate

function ismtgene(genelist)
    startswith.(uppercase.(genelist),"MT-")
end

function selectg(X,genelist)
    ng=size(X,1);
    i=vec(sum(!iszero,X,dims=2)./ng.>0.05)
    X=X[i,:];
    genelist=genelist[i];
    return X,genelist
end

function scstats(X)
    logmean=log10.(mean(X,dims=2));
    logvar=log10.(var(X,dims=2));
    dropoutrate=mean(X.==0,dims=2);
    return logmean,logvar,dropoutrate
end

function emptyrate(X)
    # nnz(sparse(X))
    count(!iszero,X)./count(isreal,X)
end

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ cab5e5da-8609-11ec-2b2c-2d98fe66bb99
begin
    import Pkg
	ENV["JULIA_MARGO_LOAD_PYPLOT"] = "no thank you"
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="ClimateMARGO", version=v"0.3.2"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="HypertextLiteral", version="0.9"),
    ])

	using Plots
	using Plots.Colors
	using ClimateMARGO
	using ClimateMARGO.Models
	using ClimateMARGO.Optimization
	using ClimateMARGO.Diagnostics
	using PlutoUI
    using HypertextLiteral

	Plots.default(linewidth=5)
end;

# ╔═╡ cc8614a3-dd69-42d4-9954-44885bac93eb
md"## Exploring uncertainty in climate forcing, feedbacks, and ocean heat uptake.

*by Henri F. Drake*
"

# ╔═╡ 84ff8ba3-fd8a-4bc9-ab4e-7f763f98d34d
md"##### Anthropogenic forcing"

# ╔═╡ 8037ac4e-5c70-4aa7-9982-7442baaa83a0
md"##### The climate response"

# ╔═╡ 4575f6fd-d1ec-4fd1-bd56-268fab59a533

# ╔═╡ 7ea69497-8b83-45f9-93d3-bbce4c3c10b2
reveal_box = @bind reveal CheckBox(default=false);

# ╔═╡ 71975b7f-e8b7-4097-b573-8e383bb83868
if reveal
md"The analytical scalings are:

$\text{TCR}\propto \frac{a}{\lambda + \kappa}, \quad\quad\quad \text{ECS} \propto \frac{a}{\lambda}, \quad\quad\quad \tau = \frac{C_{D}}{\lambda}\,\frac{\lambda + \kappa}{\kappa}.$"
else
md"Reveal solution: $(reveal_box)"
end

# ╔═╡ 8b17a78f-f3d8-4777-becd-01de58974de3

# ╔═╡ 56e32b2b-8830-4f2c-97ea-7bc5820dca1b
ρ_slider = @bind ρ Slider(0.:0.5:3, default=3, show_value=true);

# ╔═╡ 7e63eaf5-6841-4f2d-9388-427d9905fd2e
do_optimize_box = @bind do_optimize CheckBox(default=false);

# ╔═╡ a442349e-f6ba-4841-93a6-d7ea7235371a
if reveal
md"#### Exercise 2: policy solutions
Explore the sensivity of ''optimal'' policy solutions. Check the box: $(do_optimize_box)
"
end

# ╔═╡ 053db221-7f7f-4358-915f-bcfb9f6b6783
if do_optimize
	md"*Note the emissions reductions (and carbon dioxide removal) and corresponding warming reductions above!*"
end

# ╔═╡ 2f4a79c4-5455-4b8d-9a85-1aac1cdad3f7
if do_optimize
	md"#### Exercise 3: Intergenerational (in)equity

You may have guessed the ''optimal'' level of warming depends strongly on economic preferences, e.g. how much the model values the welfare of future generations.

How does the ''optimal'' policy response change as we decrease the ''discount rate'' towards intergenerational equity?

 $(ρ_slider) % per year
	"
end

# ╔═╡ 42a58da0-bffc-4fd1-ab7b-0a90b0664e4e
begin
	if do_optimize
		last_year = 2800
	else
		last_year = 4000
	end
	end_year_slider = @bind end_year Slider(2100:50:last_year, default=2200, show_value=true);
end;

# ╔═╡ 97e66e74-7f38-4811-8d11-f69df3e5ce49
md"
#### Exercise 1: transient and equilibrium global warming
Lengthen the simulation's time horizon and graphically explore the scaling of the following emergent properties with the four key parameters above:

A) the Transient Climate Response (TCR), i.e. the warming in 2150 when forcing stabilizes

B) the Equilibrium Climate Sensivity (ECS), i.e. the warming when temperatures have equilibrated with

C) the timescale $\tau_{D}$ of the equilibration

$(end_year_slider)
"

# ╔═╡ 56b88523-aa8e-42d7-8637-17b0a0c6fd22

# ╔═╡ b72a5d5f-1806-4ae0-95e1-d3b1e4b845ce

# ╔═╡ 4fee292d-9cf1-40ec-b96d-0eacb0135cae
md"# Appendix"

# ╔═╡ 1a3b7318-d0a8-424a-96db-1deb03da716b
cmip5_std = Dict("a"=>0.9*log(2.), "λ"=>0.31, "κ"=>0.18, "Cd"=>62.)

# ╔═╡ 8742ea1c-54af-4a7d-8f41-1566a71c4d94
blob(el, color = "blue") = @htl("""
$(el)
""");

# ╔═╡ 24d098eb-a548-4c80-935b-0e9bec229105
function default_parameters(years)::ClimateModelParameters
	result = deepcopy(ClimateMARGO.IO.included_configurations["default"])
	result.domain = years isa Domain ? years : Domain(step(years), first(years), last(years))
	result.economics.baseline_emissions = ramp_emissions(result.domain)
    result.economics.extra_CO₂ = zeros(size(result.economics.baseline_emissions))
	return result
end;

# ╔═╡ 226ae32f-36d4-4563-8af4-1aa9d16702c6
begin
	# Create default parameter values
	years = 2020:6.0:end_year;
	default_params = default_parameters(years);
	default_params.physics.a = 6.9*log(2.);
	phys = default_params.physics
end

# ╔═╡ c834b860-866d-4f21-b51a-4fa301d17348
a_slider = @bind a Slider(
	round(phys.a-2*cmip5_std["a"], digits=2)
	:round(4*cmip5_std["a"]/100., digits=2)
	:round(phys.a+2*cmip5_std["a"], digits=2),
	default=phys.a, show_value=true
);

# ╔═╡ 27c8bca8-55f9-4289-a839-cebbbf696169
λ_slider = @bind λ Slider(
	round(phys.B-2*cmip5_std["λ"], digits=2)
	:round(4*cmip5_std["λ"]/100., digits=2)
	:round(phys.B+2*cmip5_std["λ"], digits=2),
	default=phys.B, show_value=true
);

# ╔═╡ 672287f1-f1bd-4899-9d7e-29a6a5646f10
κ_slider = @bind κ Slider(
	round(phys.κ-2*cmip5_std["κ"], digits=2)
	:round(4*cmip5_std["κ"]/100., digits=2)
	:round(phys.κ+2*cmip5_std["κ"], digits=2),
	default=phys.κ, show_value=true
);

# ╔═╡ 58441503-2998-472e-b466-5d4c832f9042
Cd_slider = @bind Cd Slider(
	max(Int(round(phys.Cd-2*cmip5_std["Cd"], digits=0)), 20)
	:Int(round(4*cmip5_std["Cd"]/100., digits=0))
	:Int(round(phys.Cd+2*cmip5_std["Cd"], digits=0)),
	default=phys.Cd, show_value=true
);

# ╔═╡ 9c7d7e6c-d2d9-4754-b2e3-72c3b9e0ec06
let
	blob(
		@htl("""

Energy Balance Model Parameters

Defaults are CMIP5 multi-model mean; ranges are ±2σ.

		Value

		CO2 forcing parameter, a [W m⁻²]
		$(a_slider)

		Feedback parameter, λ [W m⁻² K⁻¹]
		$(λ_slider)

		 Heat uptake rate, κ [W m⁻² K⁻¹]
		$(κ_slider)

		 Deep ocean heat capacity, C₀ [W yr m⁻² K⁻¹]
		$(Cd_slider)

		"""),
		"#c0ecff33"
	)
end

# ╔═╡ 8c1ddae9-a9dd-4556-9027-4f701a064604
begin
	# Overwrite defaults with slider values
	params = deepcopy(default_params)
	params.economics.ρ = ρ/100. + params.economics.γ;
	params.economics.β = 0.003;
	params.physics.T0 = 1.0;
	params.physics.a = a;
	params.physics.B = λ;
	params.physics.κ = κ;
	params.physics.Cd = Cd;
	m = ClimateModel(params);
	max_deployment = Dict("mitigate"=>1., "remove"=>0.5, "geoeng"=>0., "adapt"=>0.)
	if do_optimize
		opt = optimize_controls!(
			m, obj_option="net_benefit", max_deployment=max_deployment
		);
	end
	m;
end

# ╔═╡ 1285f653-a4a8-4c5c-b395-2f3ffd24811d
begin
	p_emit = plot(t(m), effective_emissions(m), color=:grey, linewidth=3, alpha=0.5, label="baseline fossil growth")
	if do_optimize
		plot!(p_emit, t(m), effective_emissions(m, M=true, R=true), linewidth=5, color=:deepskyblue2, label="''optimal'' policy")
	end
	plot!(p_emit, t(m), 0. *t(m), linestyle=:dash, color=:black, alpha=0.3, linewidth = 2.5, label="")
	plot!(p_emit, ylim=(-5., 15.), xlim=(t(m)[1], t(m)[end]))
	plot!(p_emit, ylabel="emissions [ppm/year]", xlabel="year")
	plot!(p_emit, size=(700, 225), margins=2.75Plots.Measures.mm)
end

# ╔═╡ 13b0c090-ff52-4eca-98a5-1306f859c58f
begin
	p = plot(t(m), T(m), color=:grey, linewidth=3, alpha=0.5, label="baseline fossil growth")
	if do_optimize
		plot!(t(m), T(m, M=true, R=true), linewidth=5, color=:deepskyblue2, label="''optimal'' policy")
	end
	if minimum(T(m, M=true, R=true)) < 0.
		plot!(ylim=(-2, max(8, maximum(T(m, M=true, R=true)))), xlim=(t(m)[1], t(m)[end]))
	else
		plot!(ylim=(0, max(8, maximum(T(m, M=true, R=true)))), xlim=(t(m)[1], t(m)[end]))
	end
	plot!(ylabel="warming (relative to P-I)", xlabel="year")
	plot!(yticks=[-1., 0., 1., 1.5, 2., 2.5, 3., 4., 5., 6., 7., 8., 9, 10.])
	plot!(size=(700, 250), margins=2.75Plots.Measures.mm)
	plot!(legend=:topleft)
end

# ╔═╡ 10befe5b-9e7c-40a6-9079-1ab5a380888b
if do_optimize
	p_discount = plot(
		t(m),
		(1 .- (m.economics.ρ - m.economics.γ)).^((t(m) .- t(m)[1])) * 100., label=""
	)
	plot!(xlabel="year", ylabel="discount factor [%]", ylim = (0, 100))
end

# ╔═╡ Cell order:
# ╟─cc8614a3-dd69-42d4-9954-44885bac93eb
# ╟─84ff8ba3-fd8a-4bc9-ab4e-7f763f98d34d
# ╟─1285f653-a4a8-4c5c-b395-2f3ffd24811d
# ╟─8037ac4e-5c70-4aa7-9982-7442baaa83a0
# ╟─13b0c090-ff52-4eca-98a5-1306f859c58f
# ╟─9c7d7e6c-d2d9-4754-b2e3-72c3b9e0ec06
# ╟─97e66e74-7f38-4811-8d11-f69df3e5ce49
# ╟─4575f6fd-d1ec-4fd1-bd56-268fab59a533
# ╟─71975b7f-e8b7-4097-b573-8e383bb83868
# ╟─7ea69497-8b83-45f9-93d3-bbce4c3c10b2
# ╟─a442349e-f6ba-4841-93a6-d7ea7235371a
# ╟─053db221-7f7f-4358-915f-bcfb9f6b6783
# ╟─8b17a78f-f3d8-4777-becd-01de58974de3
# ╟─2f4a79c4-5455-4b8d-9a85-1aac1cdad3f7
# ╟─10befe5b-9e7c-40a6-9079-1ab5a380888b
# ╟─42a58da0-bffc-4fd1-ab7b-0a90b0664e4e
# ╟─56e32b2b-8830-4f2c-97ea-7bc5820dca1b
# ╟─7e63eaf5-6841-4f2d-9388-427d9905fd2e
# ╟─56b88523-aa8e-42d7-8637-17b0a0c6fd22
# ╟─b72a5d5f-1806-4ae0-95e1-d3b1e4b845ce
# ╟─4fee292d-9cf1-40ec-b96d-0eacb0135cae
# ╠═cab5e5da-8609-11ec-2b2c-2d98fe66bb99
# ╠═226ae32f-36d4-4563-8af4-1aa9d16702c6
# ╠═8c1ddae9-a9dd-4556-9027-4f701a064604
# ╠═1a3b7318-d0a8-424a-96db-1deb03da716b
# ╟─8742ea1c-54af-4a7d-8f41-1566a71c4d94
# ╟─c834b860-866d-4f21-b51a-4fa301d17348
# ╟─27c8bca8-55f9-4289-a839-cebbbf696169
# ╟─672287f1-f1bd-4899-9d7e-29a6a5646f10
# ╟─58441503-2998-472e-b466-5d4c832f9042
# ╟─24d098eb-a548-4c80-935b-0e9bec229105

max_stars_repo_name: ClimateMARGO/interactive

max_issues_repo_name: ClimateMARGO/interactive

max_forks_repo_name: ClimateMARGO/interactive


text:
module ProjectEuler

include("utils.jl")
using .Utils

function problem001()
    #=
    Problem 1:
        If we list all the natural numbers below 10 that are multiples of 3 or 5,
        we get 3, 5, 6 and 9. The sum of these multiples is 23.
        Find the sum of all the multiples of 3 or 5 below 1000.
    =#
    return sum(i for i in 1:999 if i%3==0 || i%5==0) # => 234168
end

function problem002()
    #=
    Problem 2:
        Each new term in the Fibonacci sequence is generated by
        adding the previous two terms. By starting with 1 and 2,
        the first 10 terms will be:
            1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
        By considering the terms in the Fibonacci sequence whose
        values do not exceed four million, find the sum of the
        even-valued terms.
    =#
    sum = 0
    f1,f2 = 1,1
    while f1 <= 4_000_000
        f1, f2 = f1+f2, f1
        if f1 % 2 == 0
            sum += f1
        end
    end
    return sum
end

function problem003()
    #=
    Problem 3:
        The prime factors of 13195 are 5, 7, 13 and 29.
        What is the largest prime factor of the number 600851475143 ?
    =#
    factors = Utils.factorize(600851475143)
    return maximum(i -> i[1],factors) # largest factor
end

function problem004()
    #=
    Problem 4:
        A palindromic number reads the same both ways. The
        largest palindrome made from the product of two 2-digit
        numbers is 9009 = 91 × 99.

        Find the largest palindrome made from the product of two
        3-digit numbers.
    =#
    palindrome_products = (i*j for i=100:999, j=100:999 if string(i*j)==reverse(string(i*j)))
    return maximum(palindrome_products)
end

function problem005()
    #=
    Problem 5:
        2520 is the smallest number that can be divided by each of
        the numbers from 1 to 10 without any remainder.

        What is the smallest positive number that is evenly
        divisible by all of the numbers from 1 to 20?
        =#

        # The desired number will be the product of the factors of
        # each number from 1 to 20, excluding duplicates. For example,
        # 16 = 2^4, so that also covers 2,4,and 8
    return (2^4 * # 2,4,8,16
        3^2 * # 3, 9
        5 *
        7 *
        11 *
        13 *
        17 *
        19
    )
end

function problem006()
    #=
    Problem 6:
        The sum of the squares of the first ten natural numbers is,
            1^2 + 2^2 + ... + 10^2 = 385
        The square of the sum of the first ten natural numbers is,
            (1 + 2 + ... + 10)^2 = 55^2 = 3025
        Hence the difference between the sum of the squares of the
        first ten natural numbers and the square of the sum is
        3025 − 385 = 2640.
        Find the difference between the sum of the squares of the first
        one hundred natural numbers and the square of the sum.
    =#
    return sum(1:100)^2 - sum(i^2 for i=1:100)
end

function problem007()
    #=
    Problem 7:
    By listing the first six prime numbers:
        2, 3, 5, 7, 11, and 13,
    we can see that the 6th prime is 13.
    What is the 10,001'st prime number?
    =#
    for (i,p) in enumerate(Utils.Primes())
        if i == 10001
            return p
        end
    end
end

function problem008()
    #=
    Problem 8
        The four adjacent digits in the 1000-digit number that
        have the greatest product are 9 × 9 × 8 × 9 = 5832.

        73167176531330624919225119674426574742355349194934
        96983520312774506326239578318016984801869478851843
        85861560789112949495459501737958331952853208805511
        12540698747158523863050715693290963295227443043557
        66896648950445244523161731856403098711121722383113
        62229893423380308135336276614282806444486645238749
        30358907296290491560440772390713810515859307960866
        70172427121883998797908792274921901699720888093776
        65727333001053367881220235421809751254540594752243
        52584907711670556013604839586446706324415722155397
        53697817977846174064955149290862569321978468622482
        83972241375657056057490261407972968652414535100474
        82166370484403199890008895243450658541227588666881
        16427171479924442928230863465674813919123162824586
        17866458359124566529476545682848912883142607690042
        24219022671055626321111109370544217506941658960408
        07198403850962455444362981230987879927244284909188
        84580156166097919133875499200524063689912560717606
        05886116467109405077541002256983155200055935729725
        71636269561882670428252483600823257530420752963450

        Find the thirteen adjacent digits in the 1000-digit number
        that have the greatest product. What is the value
        of this product?
    =#

    number_str = """
    73167176531330624919225119674426574742355349194934
    96983520312774506326239578318016984801869478851843
    85861560789112949495459501737958331952853208805511
    12540698747158523863050715693290963295227443043557
    66896648950445244523161731856403098711121722383113
    62229893423380308135336276614282806444486645238749
    30358907296290491560440772390713810515859307960866
    70172427121883998797908792274921901699720888093776
    65727333001053367881220235421809751254540594752243
    52584907711670556013604839586446706324415722155397
    53697817977846174064955149290862569321978468622482
    83972241375657056057490261407972968652414535100474
    82166370484403199890008895243450658541227588666881
    16427171479924442928230863465674813919123162824586
    17866458359124566529476545682848912883142607690042
    24219022671055626321111109370544217506941658960408
    07198403850962455444362981230987879927244284909188
    84580156166097919133875499200524063689912560717606
    05886116467109405077541002256983155200055935729725
    71636269561882670428252483600823257530420752963450
    """
    number = reshape(readdlm(IOBuffer(strip(replace(number_str,""," "))),Int)',1,:) # 1x1000 array
    products = (prod(number[i:(i+12)]) for i in 1:(length(number)-12)) # get product of every rolling 13-digit subarray
    return maximum(products)
end

function problem009()
    #=
    Problem 9:
    A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,
    a^2 + b^2 = c^2
    For example, 3^2 + 4^2 = 9 + 16 = 25 = 5^2.

    There exists exactly one Pythagorean triplet for which a + b + c = 1000.
    Find the product abc.
    =#
end

end # module


text:
module TestWith

using FGenerators
using FGenerators: @with
using Test

callonce(f) = f()

@fgenerator function generate123()
    @yield 1
    ans = @with(callonce()) do
        @yield 2
        :returning
    end
    @test ans == :returning
    @yield 3
end

@testset begin
    @test collect(generate123()) == 1:3
end

end  # module


text:
# This file was generated, do not modify it. # hide
@show size(auto)
@show nrow(auto)
@show ncol(auto)


text:
# degree 2 bezier curve

function animate_bezier(px,py;steps=100)
    n = length(px)-1
    # saving all in between de Casteljau points
    bs = [zeros(2,r) for r=n+1:-1:0]
    # base case Bézier points
    bs[1][1,:] = px
    bs[1][2,:] = py
    points = zeros(steps,2)
    colors = [:green,:orange,:red,:yellow]
    c = 1
    anim = @animate for t in LinRange(0,1,steps)
        plot(;size=(700,500), axisratio=:equal, legend=false)
        plot!(px, py, linetype=:scatter)
        for i=1:n
            for j=0:n-i
                # linear interpolation between two points
                new_b = (1-t)*bs[i][:,j+1]+t*bs[i][:,j+2]
                bs[i+1][:,j+1] = new_b
                # drawing the line for linear interpolation as well as the specific point for t
                plot!(bs[i][1,j+1:j+2],bs[i][2,j+1:j+2], legend=false, color=colors[i])
                plot!([new_b[1]],[new_b[2]], linetype=:scatter, legend=false, color=colors[i])
                if i == n
                    points[c,:] = [new_b[1],new_b[2]]
                    # draw the curve until point t = LinRange(0,1,steps)[c]
                    plot!(points[1:c,1],points[1:c,2], color=:black)
                    c += 1
                end
            end
        end
    end
    gif(anim, "bezier_d2.gif", fps=30)
end

px = [0, 3, 7]
py = [2, 9, 3]
animate_bezier(px, py)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# # Cellular Automaton
#
# ## Abstract
# We have an arrangement of atoms/cells. Each cell belongs to one in a finite
# set of states. With the progression of a unit time cells change their state
# in a parallel fashion determined by the rules of the system.
#
# Depending on the initial state and thr rules of the system we get a variety
# of behaviours. Some systems lead to homogenous states, some exhibit periodic
# oscillations between states, some are chaotic non-deterministic behaviours,
# ,.etc . The study of these discrete and abstract computational systems is
# called cellular automaton.
# ![https://plato.stanford.edu/entries/cellular-automata/](https://plato.stanford.edu/entries/cellular-automata/)
# is a great resource to understand and explore these systems.
#
# # Game of Life
#
# Game of Life is a popular cellular automaton devised by British mathematician
# John Conway. It consists of a two dimensional grid of cells with a start state
# with each cell being either dead or live at any point in time.
#
# The rules of the system are as follows :-
# * A cell with less than two neighbours dies due to underpopulation.
# * A cell with two or three neighbours lives.
# * A cell with greater than three neighbours dies by overpopulation.
# * A dead cell with three live neighbours becomes live.
#
# We will implement this game of life in minecraft to generate some interesting
# patterns.

using PiCraft

# First let's write our function which will transform from one state to another.

function nextState(A::Array{Bool, 2})
    B = similar(A)
    dims = size(A)
    for i = 1:dims[1], j = 1:dims[2]
        liveNeighbours = 0
        for p = -1:1, q = -1:1
            if checkbounds(Bool, A, i + p, j + q) && (p == 0 && q == 0)
                liveNeighbours += 1
            end
        end
        if A[i, j] == false && liveNeighbours == 3
            B[i, j] = true
        elseif liveNeighbours < 2 || liveNeighbours > 3
            B[i, j] = false
        end
    end
    return B
end

# Now to render an arbitrary 2-D board in Minecraft
# Let live cells be white wool and dead be black wool

function displayState(p::Tuple{Int, Int, Int}, A::Array{Bool, 2})
    dims = size(A)
    for i = 1:dims[1], j = 1:dims[2]
        if A[i, j]
            setBlock(p .+ (i, j, 0), Block(35))
        else
            setBlock(p .+ (i, j, 0), Block(35, 15))
        end
    end
end

p = getTile() .+ (3, 4, 3)
A = rand(Bool, 10, 10)

for i = 1:20
    displayState(A)
    A = nextState(A)
    sleep(0.5)
end

# A variety of interesting patterns can be found ![here](http://web.stanford.edu/~cdebs/GameOfLife/)
# It is encouraged that you experiment with various patterns.

# # 3-D implementation
# If we want to consider a 3-D implementation of the Game of Life then we must
# realise that in 3-D a cell will have 26 Neighbours (3^3 - 1).
#
# We generally express a 3-D system as a `Life wxyz` system.
# The rules of such a system are :
#
# * The number of live neighbours must be between `w` and `x` to sustain life.
# * The number of live neighbours must be between `y` and `z` to bring a dead
#   back to life.
#
# Conway's Game of Life is hence a `Life 2333` system.
# Since a large number of combinations are possible they have been studied
# extensively for their properties. Most systems quickly collapse or explode.
# Life 4555 and 5766 are the only valid systems which exhibit life like
# properties.
#
# Let's implement Life 4555. For that we will rewrite the `nextState` and `displayState`
# funcitons.

function nextState(A::Array{Bool, 3})
    B = similar(A)
    dims = size(A)
    for i = 1:dims[1], j = 1:dims[2], k = 1:dims[3]
        liveNeighbours = 0
        for p = -1:1, q = -1:1, r = -1:1
            if checkbounds(Bool, A, i + p, j + q, k + r) && (p == 0 && q == 0 && r == 0)
                liveNeighbours += 1
            end
        end
        if A[i, j, k] == false && liveNeighbours == 5
            B[i, j, k] = true
        elseif liveNeighbours < 4 || liveNeighbours > 5
            B[i, j, k] = false
        end
    end
    return B
end

function displayState(p::Tuple{Int, Int, Int}, A::Array{Bool, 3})
    dims = size(A)
    for i = 1:dims[1], j = 1:dims[2], k = 1:dims[3]
        if A[i, j, k]
            setBlock(p .+ (i, j, k), Block(35))
        else
            setBlock(p .+ (i, j, k), Block(35, 15))
        end
    end
end

p = getTile() .+ (3, 4, 3)
A = rand(Bool, 7, 7, 7)

for i = 1:20
    displayState(A)
    A = nextState(A)
    sleep(0.5)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7ef4b0b3-9d0a-4356-bd16-14487f14e0ab
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 0cd92d78-1da6-435b-ad8e-60af29502e1d
begin
	# Necessary packages
	# If a package is missing, follow the instructions.
	using LightGraphs
	using GraphPlot
	using SparseArrays
	using LinearAlgebra
	using Plots
	using Distances
	using Arpack
	using Images
	using Random
end

# ╔═╡ cb61f761-b2ee-4818-9149-f85d95b76a1e
md"""
# Spectral Graph Bipartitioning

Many data clustering problems can be interpreted as clustering of vertices of graphs. __Graph bipartitioning problem__ is to partition vertices into subsets such that the connections within subsets are stronger than the connections between different subsets.

Partition of the vertices into two subsetts is done according to signs of the eigenvectors of the second smallest eigenvalue of the Laplacian matrix.

__Prerequisites__

The reader should be familiar with the basic graph theory, linear algebra and, in particular,  eigenvalues and eigenvectors.

__Competences__

The reader should be able to apply graph spectral bipartitioning and recursive bipartitioning to data clustering problems.

__Credits.__ The notebook was initially derived from M.Sc. Thesis of Ivančica Mirošević.
"""

# ╔═╡ 5926c79a-6038-4dbc-b49a-5b05d02250a0
md"""
# Graphs

For more details, see [W. H. Haemers, Matrices and Graphs, in L. Hogben, Ed., 'Handbook of Linear Algebra', pp. 39.1-39.14, CRC Press, Boca Raton, 2014.](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and [S. Butler and F. Chung, Spectral Graph Theory, ibid., pp. 47.1-47.6](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

A __weighted graph__ is an ordered triplet $G=(V,E,\omega)$, where $V=\{1,2,3,...,n\}$ is the set of __vertices__ , $E=\{(i,j)\}$ is a set of __edges__ connecting vertices, and $\omega$ is a set of __weights__ of edges. We assume $G$ is undirected.
"""

# ╔═╡ a66a6e80-1f80-11eb-15f8-71ffaa5050e7
md"
__Adjacency matrix__ of graph $G$ is the matrix $A$ defined as

$$A_{ij}=\begin{cases} 1 \quad \textrm{if}\ (i,j)\in E, \\
0\quad  \textrm{otherwise} \end{cases}.$$

__Weight matrix__ of graph $G$ is the matrix $W$ defined as

$W_{ij}=\begin{cases} \omega(e) \quad \textrm{if}\ e=(i,j)\in E, \\
0\quad  \textrm{otherwise} \end{cases}.$

__Laplacian matrix__ of graph $G$ is the matrix

$L=D-W,$

where
$D=\mathop{\mathrm{diag}}(d_1,d_2,\ldots,d_n)$ with $d_i=\sum_{k=1}^n W_{ik}$ for $i=1,\ldots,n$.

__Normalized Laplacian matrix__ is the matrix

$L_n=D^{-1/2} L D^{-1/2}\equiv D^{-1/2} (D-W) D^{-1/2}$

(__diagonally scaled $L$__).

__Incidence matrix__ of graph $G$ is the $|V|\times |E|$ matrix $I_G$. Each row of $I_G$ corresponds to a vertex of $G$ and each column corresponds to an edge of $G$.
In the column corresponding to en edge $e=(i,j)$, all elements are zero except the ones in the $i$-th and $j$-th row, which are equal to $\sqrt{\omega(e)}$ and $-\sqrt{\omega(e)}$, respectively.
"

# ╔═╡ d4cf12ce-12fb-4562-bbea-e9b85cba94c2
md"""
### Examples

Graph types and algorithms are implemented in the package [LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl).

Plotting graphs is done by the packages
[GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl).

As a small inconvenience, we can only plot unweighted graphs and plot weights as edge labels.
"""

# ╔═╡ 87cca1d3-b944-4d74-bf4c-1562f34cf91a
begin
	# Sources, targets and weights
	n=7
	sn=[1,1,1,2,2,3,3,3,5,5,6]
	tn=[2,3,4,4,5,4,6,7,6,7,7]
	wn=[2,3,4,7,1,3,2,1,7,3,5]
	[sn tn wn]
end

# ╔═╡ 3d7b4a01-9a7d-4111-a45d-1d31f1ef6278
begin
	# Create the graph
	G=Graph(n)
	for i=1:length(sn)
	    add_edge!(G,sn[i],tn[i])
	end
	G
end

# ╔═╡ 3b5ed1fa-75b0-43fb-a918-4c2e8f52d039
# What is the optimal bipartition?
gplot(G, nodelabel=1:n, edgelabel=wn)

# ╔═╡ 8f7d1f6c-ee4c-407f-80d9-d2abc2af949e
begin
	# We define some functions
	function WeightMatrix(src::Array,dst::Array,weights::Array)
	    n=nv(G)
	    sparse([src;dst],[dst;src],[weights;weights],n,n)
	end

	Laplacian(W::AbstractMatrix)=spdiagm(0=>vec(sum(W,dims=2)))-W

	function NormalizedLaplacian(L::AbstractMatrix)
	    D=1.0./sqrt.(diag(L))
		Diagonal(D)*L*Diagonal(D)
	end
end

# ╔═╡ b7f127d8-e8cd-4771-8091-c88f9ffe6b08
W=WeightMatrix(sn,tn,wn)

# ╔═╡ 3f5447bb-95cd-4520-b0b9-e85f33b323bd
Matrix(W)

# ╔═╡ ecb2a46e-f9dc-4e57-812f-0a2788afb202
begin
	L=Laplacian(W)
	Matrix(L)
end

# ╔═╡ 86e5f00f-13ee-4867-b0cd-187f2c782a46
Lₙ=NormalizedLaplacian(L)

# ╔═╡ cf26dbe9-3c4f-4c95-88ef-f197f360d759
Matrix(Lₙ)

# ╔═╡ 8fcd2440-4013-4b37-a747-2c0c7f6d4aeb
issymmetric(Lₙ)

# ╔═╡ 0c8d3572-3436-446f-952a-b19d0370ab38
# Let us compute the incidence matrix
function IncidenceMatrix(G::Graph, weights::Array)
    A=zeros(nv(G),ne(G))
    k=1
    for a in edges(G)
        A[a.dst,k]=sqrt.(weights[k])
        A[a.src,k]=-sqrt(weights[k])
        k+=1
    end
    A
end

# ╔═╡ 1437aadc-1c26-4eff-bb22-9cf9c15e1f5f
Iᵧ=IncidenceMatrix(G,wn)

# ╔═╡ 6529324b-49f7-4f01-b5cd-33668979b56b
md"""
## Facts

1.  $L=I_{G}I_{G}^{T}$.

2.  $L$ is symmetric PSD matrix.

3.  $L\mathbf{1}=0$ for $\mathbf{1}=[1,...,1]^{T}$, thus $0$ is an eigenvalue of $L$  and $\mathbf{1}$ is the corresponding eigenvector.

4. If $G$ has $c$ connected components, then $L$ has $c$ eigenvalues equal to $0$.

5. For every $x\in \mathbb{R}^{n}$, it holds
$x^{T}L x=\sum\limits_{i<j}W_{ij}(x_{i}-x_{j})^{2}$.

6. For every $x\in\mathbb{R}^{n}$ and $\alpha,\beta\in\mathbb{R}$, it holds
$(\alpha x+\beta \mathbf{1})^{T} L (\alpha x+\beta \mathbf{1})
=\alpha^{2} x^{T}L x$.

7. Assume that the eigenvalues of $L$ are increasingly ordered. Then,

$$
0=\lambda_1(L)\leq \lambda_2(L)\leq \cdots \leq\lambda_{n}(L)\leq
2\max\limits_{i=1,\cdots ,n}d_{i}.$$

8.  $\sigma(L_n) \subseteq [0,2]$.
"""

# ╔═╡ a9ba64a6-bd17-4b86-acb4-e7c7d8f87dc1
md"""
### Examples
"""

# ╔═╡ 6a6f5db5-7dd4-4a26-ac1a-ef3ca3a5d6be
# Fact 1
norm(L-Iᵧ*Iᵧ')

# ╔═╡ e11bc23d-03e2-4780-9488-93f0adfc2fcb
# Facts 2 and 7
issymmetric(L), eigs(L)[1], 2*maximum(diag(L))

# ╔═╡ 7192f7d0-28b9-11eb-28dd-8db0019504ca
eigen(Matrix(L))

# ╔═╡ e5890dca-368f-437e-bc45-6a395714dae2
# Fact 3
L*ones(n)

# ╔═╡ 61a978a7-b2d2-43e7-b3f7-95bf4c4f54b9
begin
	# Fact 5
	x=rand(n)
	x'*L*x, sum([W[i,j]*(x[i]-x[j])^2 for i=1:n, j=1:n])/2
end

# ╔═╡ da7f85ed-73b0-44b9-a97f-4a8ed5e9a1d1
begin
	# Fact 6
	α,β=rand(),rand()
	(α*x+β*ones(n))'*L*(α*x+β*ones(n)), α^2*x'*L*x
end

# ╔═╡ fed7f592-be03-4baf-8ace-ff9ee78fee71
# Fact 8
eigvals(Matrix(Lₙ))

# ╔═╡ 9f31928e-c933-43bb-9a5c-bf55f29737cd
md"""
# Bipartitioning

## Definitions

Let $\pi=\{V_{1},V_{2}\}$ be a partition of $V$ with $V_1,V_2\neq \emptyset$.

__Cut__ of partition $\pi$ is the sum of weights of all
edges between $V_1$ and $V_2$,

$$\mathop{\mathrm{cut}}(\pi)\equiv \mathop{\mathrm{cut}}(V_1,V_2)=\sum\limits_{{\displaystyle i\in V_{1} \atop \displaystyle j\in V_{2}}}W_{ij}.$$

__Weight__ of vertex $i\in V$ is the sum of the weights of all egdges emanating from $i$,
$\omega(i)=\sum\limits_{j=1}^{n}W_{ij}$.

__Weight__ of a subset $\bar V\subset V$ is the sum of the weights of all vertices in $\bar V$,

$\omega(\bar V)=\sum\limits_{\displaystyle i\in\bar V} \omega(i)$.

__Proportional cut__ of partition $\pi$ is

$$
\mathop{\mathrm{pcut}}(\pi)=\displaystyle\frac{\mathop{\mathrm{cut}}(\pi)}{|V_{1}|}+\frac{\mathop{\mathrm{cut}}(\pi)}{|V_{2}|}.$$

__Normalized cut__ of partition $\pi$ is

$$
\mathop{\mathrm{ncut}}(\pi)=\displaystyle\frac{\mathop{\mathrm{cut}}(\pi)}{\omega(V_{1})}+\frac{\mathop{\mathrm{cut}}(\pi)}{\omega(V_{2})}.$$
"""

# ╔═╡ 1e6d4c72-f25d-4bf5-99d6-ae0446be25dd
md"""
### Example

Consider the following partitions (all edges have unit weights):

 $(load(\"./files/cut2.png\"))

Left partition is $\pi$, right partition is $\pi'$.

|     Cut \ Partition  | $\pi$            |  $\pi'$     |
| ------- | ---------------- | ------------|
| $\mathop{\mathrm{cut}}$  |  $2$        |     $3$        |
| $\mathop{\mathrm{pcut}}$ | $\frac{2}{1}+\frac{2}{11}=2.18$|    $\frac{3}{6}+\frac{3}{6}=1$      |
| $\mathop{\mathrm{ncut}}$ | $\frac{2}{2}+\frac{2}{50}=1.04$ |  $\frac{3}{27}+\frac{3}{25}=0.23$|
"""

# ╔═╡ 02035eed-6d49-4253-a755-68a999e7e90e
md"""
## Facts

1. The informal description of the bipartitioning problem can be formulated as two problems,

$$
\mathop{\textrm{arg min}}\limits_{\pi} \mathop{\mathrm{pcut}}(\pi) \quad \textrm{or} \quad
\mathop{\textrm{arg min}}\limits_{\pi} \mathop{\mathrm{ncut}}(\pi).$$

The first problem favors partitions into subsets with similar numbers of vertices, while the second problem favors partitions into subsets with similar weights.

2. Both problems are NP-hard.

3. __Approximate solutions can be computed by suitable relaxations in $O(n^2)$ operations.__

4. The partition $\pi$ is defined by the vector $y$ such that

$$
y_{i}=
\begin{cases}
\frac{1}{2} & \text{for } i\in V_1 \\
-\frac{1}{2} & \text{for } i\in V_2
\end{cases}$$

The proportional cut problem can be formulated as the  __discrete proportional cut__ problem

$$
\underset{\displaystyle \big|\mathbf{y}^{T}\mathbf{1} \big|\leq \beta}
{\min\limits_{\displaystyle y_{i}\in \{-\frac{1}{2},\frac{1}{2}\}}}
\frac{1}{2}\sum_{i,j}(y_{i}-y_{j})^{2}W_{ij}.$$

Parameter $\beta$ controls the number of vertices in each subset.

5. The normalized cut problem can be formulated as the __discrete normalized cut__ problem

$$
\underset{\displaystyle \big|y^{T}D\mathbf{1} \big|\leq \beta}
{\min\limits_{\displaystyle y_{i}\in \{-\frac{1}{2},\frac{1}{2}\}}}
\frac{1}{2}\sum_{i,j}(y_{i}-y_{j})^{2}W_{ij}.$$

Parameter $\beta$ controls the weights of each subset.

6. Using the Fact 5 above, the discrete proportional cut problem can be formulated as the __relaxed proportional cut__ problem

$$
\underset{\displaystyle y^{T}y=1}{\underset{\displaystyle \big| y^{T}\mathbf{1} \big|
\leq 2\frac{\beta}{\sqrt{n}}}
{\min\limits_{\displaystyle y\in \mathbb{R}^{n}}}} y^{T}L y.$$

Similarly, the discrete normalized cut problem can be formulated as the __relaxed normalized cut__ problem

$$
\underset{\displaystyle y^{T}Dy=1}{\underset{\displaystyle \big| y^{T}D\mathbf{1}\big|
\leq \displaystyle \frac{\beta}{\sqrt{\theta n}}}{\min\limits_{\displaystyle y\in
\mathbb{R}^{n}}}}y^{T}L_n y.$$

7. __The Main Theorem.__ Let $A\in \mathbb{R}^{n\times n}$ be a symmetric matrix with eigenvalues $\lambda _{1}<\lambda _{2}<\lambda_{3}\leq \cdots \leq \lambda _{n}$ and let $v^{[1]},v^{[2]},\ldots,v^{[n]}$ be the corresponding eigenvectors. For the fixed $0\leq \alpha <1$, the solution of the problem

$$
\underset{\displaystyle y^{T}y=1}{\underset{\displaystyle \left|y^{T}v^{[1]}\right|\leq \alpha}
{\min\limits_{\displaystyle y\in \mathbb{R}^{n}}}} y^{T}Ay$$

is $y=\pm \alpha v^{[1]}\pm \sqrt{1-\alpha^{2}}v^{[2]}$.

For the proof see [D. J. Higham and M. Kibble, A Unified View of Spectral Clustering, Theorem 3.1, p. 7](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.115.1591&rep=rep1&type=pdf).

8. For $0\leq \beta <\frac{n}{2}$, the solution of the relaxed proportional cut problem is

$$
y=\pm \frac{2\beta}{n\sqrt{n}}\mathbf{1}\pm
\sqrt{1-4\frac{\beta ^{2}}{n^{2}}}v^{[2]},$$

where $v^{[2]}$ is an eigenvector corresponding to $\lambda_2(L)$. $v^{[2]}$ the __Fiedler vector__. Since the first summand carries no information, $V$ is partitioned according to the signs of the components of $v^{[2]}$:

$$
V_{1}=\{i:v^{[2]}_i <0\}, \quad V_{2}=\{i:v^{[2]}_i \geq 0\}.$$

_Notice that the value of $\beta$ is irrelevant for the solution._

9. For $0\leq \beta <\sqrt{\theta n}\left\Vert D^{\frac{1}{2}}\mathbf{1} \right\Vert _{2},$ the solution of the relaxed normalized cut problem is

$$
y=\pm \frac{\beta }{\sqrt{\theta n}\left\Vert
D^{\frac{1}{2}} \mathbf{1}\right\Vert _{2}^{2}}\mathbf{1}\pm
\sqrt{1-\frac{\beta ^{2}}{ \theta n\left\Vert
D^{\frac{1}{2}}\mathbf{1}\right\Vert _{2}^{2}}}D^{-\frac{1
}{2}} v_n^{[2]},$$

where $v_n^{[2]}$ is an eigenvector corresponding to $\lambda_2(L_n)$. $V_n$ is partitioned according to the signs of the components of $v_n^{[2]}$, as above.

10. Neither of the relaxed algorithms is guaranteed to solve exactly the true proportional / normalized cut problem. However, the computed solutions are in the right direction. Whether to use proportional or normalized cut formulation, depends upon the specific problem.
"""

# ╔═╡ c188b7cf-03e7-40b3-bf90-773978c7b410
# Voila!
eigs(L,nev=2,which=:SM, v0=ones(n))

# ╔═╡ f60b51b6-79c3-4a3f-8bfb-8f198cdb92a8
# For the normalized cut
eigs(Lₙ,nev=2,which=:SM, v0=ones(n))

# ╔═╡ bfac7a55-0d82-4549-b999-7db3272eaa20
md"""
### Concentric rings

A __complete graph__ has edges connecting each pair of vertices.

To a set of points $X=\{x_{1},x_{2},\cdots ,x_{m}\}$ , where $x_{i}\in\mathbb{R}^{n}$, we assign a weighted complete graph $G=(V,E)$ with $m$ vertices, where the vertex $j\in V$ corresponds to the point $x_j\in X$.

The main idea is to assign weight of an edge $e=(i,j)$ which reflects the distance between $x_i$ and $x_j$, something like $\omega(e)=\displaystyle\frac{1}{\mathop{\mathrm{dist}}(x_i,x_j)}$.

However, this has to be implemented with care. For example, using simple Euclidean distance yield the same results as the function `kmeans()`. In this example we use Gaussian kernel, that is

$$
\omega(e)=e^{\displaystyle -\|x_i-x_j\|_2^2/\sigma^2},$$

where the choice of $\sigma$ is based on experience.

The computation of various distances is implemented in the package [Distances.jl](https://github.com/JuliaStats/Distances.jl).

We will construct the Laplace matrix directly.
"""

# ╔═╡ 80ed91ad-6f9a-4f12-a1ee-bc8dc6b58168
begin
	# Two concentric circles
	k=2
	# Center
	Random.seed!(541)
	# center=[rand(-5:5),rand(-5:5)]
	center=[0,0]
	# Radii
	radii=randperm(10)[1:k]
	# Number of points in circles
	sizes=rand(1000:2000,k)
	center,radii,sizes
end

# ╔═╡ 2d8689c9-1b9e-4360-865c-23c2c8c5dcb9
begin
	# Generate points
	X=Array{Float64}(undef,2,sum(sizes))
	csizes=cumsum(sizes)
	# Random angles
	ϕ=2*π*rand(sum(sizes))
	for i=1:csizes[1]
		X[:,i]=center+radii[1]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
	end
	for j=2:k
		for i=csizes[j-1]+1:csizes[j]
			X[:,i]=center+radii[j]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
		end
	end
	scatter(X[1,:],X[2,:],title="Concentric rings", aspect_ratio=1,label="Points")
end

# ╔═╡ fa4ae713-a0b2-46ce-b5f7-558b26d82019
# Weight matrix
W₁=1 ./pairwise(SqEuclidean(),X)

# ╔═╡ 760e8ca2-596f-435f-befb-140b3d204f7f
begin
	# Laplacian matrix
	m=csizes[end]
	for i=1:m
	    W₁[i,i]=0
	end
	L₁=Diagonal(vec(sum(W₁,dims=2)))-W₁
	# Check Fact 3
	norm(L₁*ones(m))
end

# ╔═╡ 75956488-76ab-4262-a32f-0c53e1b3dc17
# Notice λ₁=0
E=eigs(L₁,nev=2,which=:SM, v0=ones(m))

# ╔═╡ a55260dd-f05d-4a65-9fdf-cb3fd4638857
begin
	# Define clusters
	C=ones(Int64,m)
	C[findall(E[2][:,2].>0)].=2
	C
end

# ╔═╡ ba164ff3-5ef4-4480-8af2-9dd9a99b1e8d
# Yet another plotting function
function plotKpartresult(C::Vector,X::Array)
	scatter(aspect_ratio=1)
    k=maximum(C)
    for j=1:k
        scatter!(X[1,findall(C.==j)],X[2,findall(C.==j)],label="Cluster $j")
    end
	scatter!(aspect_ratio=1)
end

# ╔═╡ 46de28bc-73ed-495c-a2e0-f565c1ce5651
plotKpartresult(C,X)

# ╔═╡ e004e4bf-ca04-4b41-947b-9f1ecd058abb
md"""
This is the same partitioning as obtained by `kmeans()`. Let us try Gaussian kernel. A rule of thumb is: if rings are close, use $\sigma<1$, if rings are apart, use $\sigma>1$.
"""

# ╔═╡ bb9d4d98-3aff-439e-b8a3-8347cf345839
begin
	σ=0.7 # 0.1
	W₂=exp.(-pairwise(SqEuclidean(),X)/σ^2)-I
	L₂=Diagonal(vec(sum(W₂,dims=2)))-W₂
	E₂=eigs(L₂,nev=2,which=:SM, v0=ones(m))
	C₂=ones(Int64,m)
	C₂[findall(E₂[2][:,2].>0)].=2
	plotKpartresult(C₂,X)
end

# ╔═╡ 012fa54a-5171-4980-a0ba-474a22c25981
md"""
# Recursive bipartitioning

## Definitions

Let $G=(V,E)$ be a weighted graph with weights $\omega$.

Let $\pi_k =\{V_{1},V_{2},...,V_{k}\}$ be a $k$-partition of $V$, with $V_i\neq \emptyset$ for $i=1,\ldots,k$.

The previous definition of $cut(\pi)\equiv cut(\pi_2)$ extends naturally to $k$-partition.

A __cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{cut}}(\pi_k)=\sum\limits_{\displaystyle i<j} \mathop{\mathrm{cut}}(V_{i},V_{j}),$$

where $\mathop{\mathrm{cut}}(V_{i},V_{j})$ is interpreted as a cut of the bipartition of the subgraph of $G$ with vertices $V_1\cup V_2$.

__Proportional cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{pcut}}(\pi_k)=\underset{i<j}{\sum\limits_{i,j=1}^{k}} \left(
\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{|V_{i}|}+\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{|V_{j}|}\right) =
\sum_{i=1}^{k}\frac{\mathop{\mathrm{cut}}(V_{i},V\backslash V_{i})}{|V_{i}|}.$$

__Normalized cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{ncut}}(\pi_k)=\underset{i<j}{\sum\limits_{i,j=1}^{k}} \left(
\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{\omega(V_{i})}+\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{\omega(V_{j})}\right) =
\sum_{i=1}^{k}\frac{\mathop{\mathrm{cut}}(V_{i},V\backslash V_{i})}{ \omega(V_{i})}.$$

## Facts

If we want to cluster vertices of graph $G=(V,E)$ into $k$ clusters, we can apply the following recursive algorithm:

1. __Initialization.__ Compute the bipartition $\pi=\{V_{1},V_{2}\}$ of $V$. Set the counter $c=2$.

2. __Recursion.__ While $c<k$ repeat:

    1. Compute the bipartition of each subset of $V$.

    2. Among all $(c+1)$-partitions, choose the one with the smallest $\mathop{\mathrm{pcut}}(\pi_{c+1})$ or $\mathop{\mathrm{ncut}}(\pi_{c+1})$, respectively.

    3. Set $c=c+1$.

3. __Stop.__

There is no guarantee for optimality of this algorithm. Clearly, the optimal $k$-partiton may be a subpartition of one of the discarded partitions.
"""

# ╔═╡ f3eb8ad4-55a8-4e76-a581-07f9a7cc75f7

# ╔═╡ Cell order:
# ╟─7ef4b0b3-9d0a-4356-bd16-14487f14e0ab
# ╟─cb61f761-b2ee-4818-9149-f85d95b76a1e
# ╟─5926c79a-6038-4dbc-b49a-5b05d02250a0
# ╟─a66a6e80-1f80-11eb-15f8-71ffaa5050e7
# ╟─d4cf12ce-12fb-4562-bbea-e9b85cba94c2
# ╠═0cd92d78-1da6-435b-ad8e-60af29502e1d
# ╠═87cca1d3-b944-4d74-bf4c-1562f34cf91a
# ╠═3d7b4a01-9a7d-4111-a45d-1d31f1ef6278
# ╠═3b5ed1fa-75b0-43fb-a918-4c2e8f52d039
# ╠═8f7d1f6c-ee4c-407f-80d9-d2abc2af949e
# ╠═b7f127d8-e8cd-4771-8091-c88f9ffe6b08
# ╠═3f5447bb-95cd-4520-b0b9-e85f33b323bd
# ╠═ecb2a46e-f9dc-4e57-812f-0a2788afb202
# ╠═86e5f00f-13ee-4867-b0cd-187f2c782a46
# ╠═cf26dbe9-3c4f-4c95-88ef-f197f360d759
# ╠═8fcd2440-4013-4b37-a747-2c0c7f6d4aeb
# ╠═0c8d3572-3436-446f-952a-b19d0370ab38
# ╠═1437aadc-1c26-4eff-bb22-9cf9c15e1f5f
# ╟─6529324b-49f7-4f01-b5cd-33668979b56b
# ╟─a9ba64a6-bd17-4b86-acb4-e7c7d8f87dc1
# ╠═6a6f5db5-7dd4-4a26-ac1a-ef3ca3a5d6be
# ╠═e11bc23d-03e2-4780-9488-93f0adfc2fcb
# ╠═7192f7d0-28b9-11eb-28dd-8db0019504ca
# ╠═e5890dca-368f-437e-bc45-6a395714dae2
# ╠═61a978a7-b2d2-43e7-b3f7-95bf4c4f54b9
# ╠═da7f85ed-73b0-44b9-a97f-4a8ed5e9a1d1
# ╠═fed7f592-be03-4baf-8ace-ff9ee78fee71
# ╟─9f31928e-c933-43bb-9a5c-bf55f29737cd
# ╟─1e6d4c72-f25d-4bf5-99d6-ae0446be25dd
# ╟─02035eed-6d49-4253-a755-68a999e7e90e
# ╠═c188b7cf-03e7-40b3-bf90-773978c7b410
# ╠═f60b51b6-79c3-4a3f-8bfb-8f198cdb92a8
# ╟─bfac7a55-0d82-4549-b999-7db3272eaa20
# ╠═80ed91ad-6f9a-4f12-a1ee-bc8dc6b58168
# ╠═2d8689c9-1b9e-4360-865c-23c2c8c5dcb9
# ╠═fa4ae713-a0b2-46ce-b5f7-558b26d82019
# ╠═760e8ca2-596f-435f-befb-140b3d204f7f
# ╠═75956488-76ab-4262-a32f-0c53e1b3dc17
# ╠═a55260dd-f05d-4a65-9fdf-cb3fd4638857
# ╠═ba164ff3-5ef4-4480-8af2-9dd9a99b1e8d
# ╠═46de28bc-73ed-495c-a2e0-f565c1ce5651
# ╟─e004e4bf-ca04-4b41-947b-9f1ecd058abb
# ╠═bb9d4d98-3aff-439e-b8a3-8347cf345839
# ╟─012fa54a-5171-4980-a0ba-474a22c25981
# ╠═f3eb8ad4-55a8-4e76-a581-07f9a7cc75f7

max_stars_repo_path: src/Module C - Applications/L10 Spectral Graph Bipartitioning.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using Makie, GLMakie, ModernGL
function rotate(p::Array{N}, angle) where {N <: Number}

    matrix = [cos(angle) -sin(angle) 0; sin(angle) cos(angle) 0; 0 0 1]

    output = matrix*vcat(p, 1)
    return output[1:2]
end

function plot_makie()
    points = [Point{3,Float64}(vcat(rotate([1,0],2*pi*i/100), 1)) for i = 1:101]

    for i = 1:100
        scene = Scene()
        angle = 2*pi*i/100
        point = Point{3,Float64}(vcat(rotate([1,0], angle), 1))
        scatter!(scene, point, color = :red, markersize = 50)
        lines!(scene, points, color = :black)
        cam3d!(scene)
        zlims!(scene, (0,2))

        save("scene"*lpad(i, 4, "0")*".png", scene; resolution = (600,400))
    end
end

max_stars_repo_name: leios/simuleios

max_issues_repo_name: leios/simuleios

max_forks_repo_name: leios/simuleios


text:
export Polynomial, print

"""
    Polynomial(poly::SortedDict{Exponent, Number}, degree::Degree)

Define a mathematical polynomial, that is a linear combinaison of product of
`Variable` and conjugated `Variable`.

### Attributes
- `poly` : a dictionary associating an `Exponent` to a a number.
- `degree`: a `Degree` object indicating the global degree of the polynomial :
the maximum of the explicit exponent degrees and of the conjugated exponent
degrees.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > p = a*b*conj(b) + (2+3im) * b^3
julia > p.poly = SortedDict(a*b*conj(b)=>1, b^3=>2+3im)
julia > p.degree = (3,1)
```
"""
struct Polynomial <: AbstractPolynomial
    poly::SortedDict{Exponent, Number}
    degree::Degree

    function Polynomial(poly::SortedDict{Exponent, Number})
        degexpl, degconj = 0, 0
        for (expo, λ) in poly
            if λ!=0
                degexpl = max(degexpl, expo.degree.explvar)
                degconj = max(degconj, expo.degree.conjvar)
            else
                delete!(poly, expo)
            end
        end
        return new(poly, Degree(degexpl, degconj))
    end
end

#############################
## Constructors
#############################
Polynomial() = Polynomial(SortedDict{Exponent, Number}())

#############################
## Equality, hash
#############################
function ==(p1::Polynomial, p2::Polynomial)
  if (length(p1) != length(p2)) || (p1.degree != p2.degree)
    return false
  end

  state1 = start(p1)
  state2 = start(p2)
  while !done(p1, state1) && !done(p2, state2)
    (i1, state1) = next(p1, state1)
    (i2, state2) = next(p2, state2)

    if !(i1 == i2)
      return false
    end
  end
  return true
end
!=(pol1::Polynomial, pol2::Polynomial) = !(pol1 == pol2)

#############################
## Print
#############################
function Base.print(io::IO, P::Polynomial)
  poly = P.poly
  i = length(poly)
  un = Exponent()
  for (expo, λ) in poly

    if expo == un
      print(io, "$λ")
    else
      if λ != 1
        print(io, "(", λ, ")*")
      end
      print(io, expo)
    end
    if i > 1
      print(io, " + ")
    end
    i -= 1
  end
end


text:
# Cost functions

mse(ŷ, y) = sum((ŷ .- y).^2)/length(y)

crossentropy(ŷ::AbstractVecOrMat, y::AbstractVecOrMat) =
  -sum(y .* log.(ŷ)) / size(y, 2)

@deprecate logloss(x, y) crossentropy(x, y)

function logitcrossentropy(logŷ::AbstractVecOrMat, y::AbstractVecOrMat)
  logŷ = logŷ .- maximum(logŷ, 1)
  ypred = logŷ .- log.(sum(exp.(logŷ), 1))
  -sum(y .* ypred) / size(y, 2)
end


text:
# Tom Anderson
# vdiv_filter.jl
# MIT License
# 10/30/2016
# Demo program to show the frequency response of an
# RC lowpass filter and an RLC lowpass filter

include("/Users/toma/tools/julia/rc_graph.jl")

fstart = 1
fstop = 10meg
npts = 2810
freq = logsweep(fstart, fstop, npts)

# RC lowpass filter
# R=100 Ohms, C= 0.01uF
filt1 = db(vdiv(zr(100,freq), zc(.01u,freq)))

# RLC lowpass filter with peaking.
# R=100 Ohms, L=1mH, c=0.01uF
filt2 = db(vdiv(zr(100,freq) + zl(1m,freq), zc(.01u,freq)))

plt1 = FramedPlot()
setattr(plt1, "xlog", true)
setattr(plt1, "xrange", (fstart, fstop))
setattr(plt1.frame1, "draw_grid", true)
setattr(plt1.x1, "label", "Frequency, Hz")
setattr(plt1.y1, "label", "LPF out, dB")

a = Curve(freq, filt1, "color", "blue")
b = Curve(freq, filt2, "color", "black")
add(plt1, a)
add(plt1, b)

max_stars_repo_name: tomacorp/PyPost

max_issues_repo_name: tomacorp/PyPost

max_forks_repo_name: tomacorp/PyPost


text:
using DynamicGrids, Dispersal, Test, Colors, Flatten, LossFunctions, Optim
using Dispersal: stepfromframe, SingleCoreReplicates, ThreadedReplicates
using DynamicGrids: grid2image

@testset "step from frame" begin
    @test stepfromframe(1, 1, 1) == 1
    @test stepfromframe(1, 5, 1) == 5
    @test stepfromframe(2, 1, 1) == 1
    @test stepfromframe(2, 2, 1) == 1
    @test stepfromframe(2, 3, 1) == 2
    @test stepfromframe(2, 4, 1) == 2
    @test stepfromframe(2, 5, 1) == 3
    @test stepfromframe(3, 4, 1) == 2
    @test stepfromframe(3, 5, 1) == 2
end

@testset "SimpleObjective" begin

    init = [1.0 0.0 1.0
            0.0 0.0 0.0
            1.0 1.0 0.0]

    target = [8.0 0.0 8.0
              0.0 0.0 0.0
              8.0 8.0 0.0]

    # Test the last simulation frame matches the target when the correct rate is used
    tspan = 1:4
    ngroups = 1
    groupsize = 1
    rate = log(2.0)
    ruleset = Ruleset(ExactExponentialGrowth(intrinsicrate=rate))
    output = ArrayOutput(init; tspan=tspan)
    sim!(output, ruleset)

    @test output[end] == target

    # Test optim finds the correct intrinsic growth rate
    output = ArrayOutput(init; tspan=tspan)
    objective = SimpleObjective(target)
    loss = LogitDistLoss()

    p = Parametriser(ruleset, output, objective, identity, loss, ngroups, groupsize);
    @test DynamicGrids.ruleset(p::Parametriser) === ruleset
    @test Dispersal.output(p::Parametriser) === output
    @test Dispersal.transform(p::Parametriser) == p.transform
    @test Dispersal.objective(p::Parametriser) === objective
    @test Dispersal.loss(p::Parametriser) === loss
    @test Dispersal.ngroups(p::Parametriser) == 1
    @test Dispersal.groupsize(p::Parametriser) == 1
    @test Dispersal.threading(p::Parametriser) == SingleCoreReplicates()
    @test DynamicGrids.tspan(p::Parametriser) == 1:4
    @test Dispersal.data(p::Parametriser) == p.data
    @test Dispersal.targetbuffer(p::Parametriser) == target
    @test Dispersal.targetbuffer(p::Parametriser) !== target
    @test Dispersal.predictionbuffer(p::Parametriser) == target
    @test Dispersal.predictionbuffer(p::Parametriser) !== target

    res = Optim.optimize(p, [log(1.8)], [log(2.2)], [log(1.85)], SAMIN(), Optim.Options(iterations=1000))
    @test res.minimizer[1] ≈ rate atol=4

end

@testset "RegionObjective" begin
    init =  [1.0  1.0  0.5
             0.0  0.0  0.0
             0.0  0.0  0.0]

    regionlookup = [1  2  3
                    1  2  3
                    2  2  0]

    occurance = Bool[0  1  0  # 1
                     1  0  1  # 2
                     1  1  0] # 3

    ngroups = 1
    groupsize = 1
    framesperstep = 2
    steps = 3
    tspan = 1:framesperstep * steps
    detectionthreshold = 0.1
    transform = x -> 2x - 1

    objective = RegionObjective(detectionthreshold, regionlookup, occurance, framesperstep, 1)
    output = RegionOutput(init; tspan=tspan, objective=objective)
    ruleset = Ruleset(ExactExponentialGrowth(intrinsicrate = log(2.0)))
    loss = ZeroOneLoss()

    @testset "SingleCoreReplicates" begin
        threading = SingleCoreReplicates()
        parametriser = Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize, threading)
        parametriser(flatten(ruleset))
        res = Optim.optimize(parametriser, [log(1.8)], [log(2.2)], [log(1.85)], SAMIN(), Optim.Options(iterations=1000))
    end

    @testset "ThreadedReplicates" begin
        threading = ThreadedReplicates()
        parametriser = Parametriser(ruleset, output, objective, transform, loss, ngroups, groupsize, threading)
        parametriser(flatten(ruleset))
        res = Optim.optimize(parametriser, [log(1.8)], [log(2.2)], [log(1.85)], SAMIN(), Optim.Options(iterations=1000))
    end

    # TODO test the result with a real dispersal function and predict parameters

end


text:
module ThermoModels

using ThermoState
using ThermoState.Types
using ThermoState.QuickStates
using ThermoState.StatePoints

import ThermoState: pressure,temperature,mass,moles,molar_mass
import ThermoState: mass_volume, mol_volume, total_volume
import ThermoState: mass_density, mol_density
import ThermoState: mass_enthalpy, mol_enthalpy,total_enthalpy
import ThermoState: mass_gibbs, mol_gibbs, total_gibbs
import ThermoState: mass_helmholtz, mol_helmholtz, total_helmholtz
import ThermoState: mass_internal_energy, mol_internal_energy, total_internal_energy
import ThermoState: mass_entropy, mol_entropy, total_entropy
import ThermoState: mass_fraction, mol_fraction
import ThermoState: mass_number, mol_number
import ThermoState: options, phase, quality
import ThermoState: mol_cp,mol_cv,mass_cv,mass_cp
import ThermoState: sound_speed
import ThermoState: molecular_weight

using StaticArrays
using DiffResults
using ForwardDiff

#units
using Unitful
using Unitful: @u_str

#for cubic solvers
import Polynomials

import Roots
using Roots: find_zero

#for pt-flash
import NLSolversBase
using NLSolversBase: only_fg!,OnceDifferentiable

import Optim
using Optim: optimize

import LinearAlgebra
using LinearAlgebra: dot,norm

#to support static arrays correctly
import BangBang
using BangBang: @!

# Write your package code here.
include("base/utils.jl")
include("base/suva.jl")
include("base/types.jl")
include("base/mixing.jl")
include("base/rr_flash.jl")

#this is for using saturation models on estimation of helmholtz models
include("base/saturation.jl")
include("models/saturation.jl")

include("base/helmholtz/helmholtz.jl")

include("base/helmholtz/cubic.jl")
include("base/helmholtz/volume_solver.jl")

include("base/helmholtz/equilibria_single.jl")
include("base/helmholtz/equilibria_pt.jl")

#each model file exports its own model
include("models/iapws95.jl")
include("models/gerg2008.jl")
include("models/cubic/cubics.jl")

#include("models/lennon2000air.jl")
module Saturation
    import ..ThermoModels: Antoine, LeeKesler,Edalat,Sanjari,WaterSat
    import ..ThermoModels: VdWSatApprox, RKSatApprox, PRSatApprox

    export Antoine, LeeKesler,Edalat,Sanjari,WaterSat
    export VdWSatApprox,RKSatApprox,PRSatApprox
end

export VanDerWaals, RedlichKwong
export IAPWS95, GERG2008
export Saturation
export equilibria,kvalues
export compressibility_factor, acentric_factor
export fugacity_coeff_impl,kvalues,flash_impl,normalizefrac!!

end


text:
## Empirical evaluation of our model (comparison with 6 largest food webs archived on Mangal)

# read metadata of networks archieved on Mangal
mangal_networks = CSV.read(joinpath("data", "raw", "mangal_networks.csv"))

# filter for food webs (predation or herbivory links) archieved on Mangal
mangal_fw = mangal_networks[(mangal_networks[:,:P] .> 0) .| (mangal_networks[:,:H] .> 0), :]

# sort by species richness and filter for 6 largest webs
sort!(mangal_fw, :S, rev=true)
fw_large = first(mangal_fw, 6)

# convert these large food webs to unipartite networks
fw_large_UN = convert.(UnipartiteNetwork, network.(fw_large.id))

# number of species and links in the networks
fw_large_S = richness.(fw_large_UN)
fw_large_L = links.(fw_large_UN)

# empirical degree distributions
fw_large_dd = values.(degree.(fw_large_UN))

# maxent distributions
dd_maxent = maxent_dd.(fw_large_S, fw_large_L)

# empirical degree distributions of 6 largest food webs
dd_emp = empirical_dd.(fw_large_dd)

## plots

# plot empirical and maxent degree distributions
function plot_dd(dd_emp, dd_maxent)
    scatter(dd_emp[:,1], dd_emp[:,2], color=:grey, alpha=0.4, label="", markersize=2,
    framestyle=:box, dpi=200)
    plot!(dd_maxent, color=pal_median, linewidth=2, linestyle=:solid, label="")
    xaxis!(xlabel="k", xlims=(0, maximum(dd_emp[:,1])+5))
    yaxis!(ylabel="p(k)", ylims=(0,maximum(dd_emp[:,2])*1.1))
end

plot_dd1 = plot_dd(dd_emp[1], dd_maxent[1])
plot_dd2 = plot_dd(dd_emp[2], dd_maxent[2])
plot_dd3 = plot_dd(dd_emp[3], dd_maxent[3])
plot_dd4 = plot_dd(dd_emp[4], dd_maxent[4])
plot_dd5 = plot_dd(dd_emp[5], dd_maxent[5])
plot_dd6 = plot_dd(dd_emp[6], dd_maxent[6])

plot(plot_dd1, plot_dd2, plot_dd3, plot_dd4, plot_dd5, plot_dd6, layout=(2,3),
    legend=false,
    title=["S=$(fw_large_S[i]), L=$(fw_large_L[i])" for j in 1:1, i in 1:6], titleloc=:right, titlefont=font(6))

savefig(joinpath("fig", "emp_dd"))

max_stars_repo_name: FrancisBanville/Trophic-METE

max_issues_repo_name: FrancisBanville/Trophic-METE

max_forks_repo_name: FrancisBanville/Trophic-METE


text:
module LARX

include("larx_node_description.jl")
include("larx_rules_prototypes.jl")
include("larx_update_rules.jl")
include("util.jl")

end

max_stars_repo_name: wmkouw/LARX

max_issues_repo_name: wmkouw/LARX

max_forks_repo_name: wmkouw/LARX


text:
using ModernGraphs
using Test

@testset "ModernGraphs.jl" begin
    # Write your tests here.
end


text:
using Oceananigans.BoundaryConditions: BC, Periodic

struct PressureSolver{T, A, W, S, R, C}
          type :: T
  architecture :: A
   wavenumbers :: W
       storage :: S
    transforms :: R
     constants :: C
end

abstract type PressureSolverType end

struct TriplyPeriodic       <: PressureSolverType end
struct HorizontallyPeriodic <: PressureSolverType end
struct Channel              <: PressureSolverType end
struct Box                  <: PressureSolverType end

poisson_bc_symbol(::BC) = :N
poisson_bc_symbol(::BC{<:Periodic}) = :P
poisson_bc_symbol(::Nothing) = :P

#####
##### Special forms for IncompressibleModel constructor
#####

PressureSolver(::Nothing, arch, grid, pressure_bcs) = PressureSolver(arch, grid, pressure_bcs)

# In principle we should check that the grid used to construct pressure_solver
# and the argument `grid` are identical. We don't this right now.
PressureSolver(pressure_solver::PressureSolver, arch, grid, pressure_bcs) = pressure_solver

function PressureSolver(arch, grid, pressure_bcs, planner_flag=FFTW.PATIENT)
    x = poisson_bc_symbol(pressure_bcs.x.left)
    y = poisson_bc_symbol(pressure_bcs.y.left)
    z = poisson_bc_symbol(pressure_bcs.z.left)
    bc_symbol = Symbol(x, y, z)

    if bc_symbol == :PPP
        return TriplyPeriodicPressureSolver(arch, grid, pressure_bcs, planner_flag)
    elseif bc_symbol == :PPN
        return HorizontallyPeriodicPressureSolver(arch, grid, pressure_bcs, planner_flag)
    elseif bc_symbol == :PNN
        return ChannelPressureSolver(arch, grid, pressure_bcs, planner_flag)
    elseif bc_symbol == :NNN
        return BoxPressureSolver(arch, grid, pressure_bcs, planner_flag)
    else
        throw(ArgumentError("Unsupported pressure boundary conditions: $bc_symbol"))
    end
end


text:
transD_GP.getchi2forall(opt, fsize=6)
ax = gcf().axes;
ax[2].set_ylim(100, 200)
χ² = sum(2*csem.ndata)
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
savefig("csem_conv_ns_1.png", dpi=300)
transD_GP.getchi2forall(optlog10λ, fsize=6)
ax = gcf().axes;
ax[2].set_ylim(100, 200)
χ² = sum(2*csem.ndata)
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
savefig("csem_conv_ns_2.png", dpi=300)
opt.xall[:] .= zall
transD_GP.plot_posterior(csem, opt, optlog10λ, burninfrac=0.5,
    figsize=(7.8,4), cmappdf="inferno", qp1=0.1, qp2=0.9, nbins=100)
ax = gcf().axes
ax[1].step(log10.(ρ[2:end]), z[2:end], alpha=0.4, color="w")
ax[1].set_xlim(-1.0959789629712384, 2.5166552528237873)
savefig("csem_post_ns.png", dpi=300)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA