text:
########
# Dual #
########

immutable Dual{N,T<:Real} <: Real
    value::T
    partials::Partials{N,T}
end

################
# Constructors #
################

Dual{N,T}(value::T, partials::Partials{N,T}) = Dual{N,T}(value, partials)

function Dual{N,A,B}(value::A, partials::Partials{N,B})
    T = promote_type(A, B)
    return Dual(convert(T, value), convert(Partials{N,T}, partials))
end

Dual(value::Real, partials::Tuple) = Dual(value, Partials(partials))
Dual(value::Real, partials::Tuple{}) = Dual(value, Partials{0,typeof(value)}(partials))
Dual(value::Real, partials::Real...) = Dual(value, partials)

##############################
# Utility/Accessor Functions #
##############################

@inline value(x::Real) = x
@inline value(n::Dual) = n.value

@inline partials(x::Real) = Partials{0,typeof(x)}(tuple())
@inline partials(n::Dual) = n.partials
@inline partials(x::Real, i...) = zero(x)
@inline partials(n::Dual, i) = n.partials[i]
@inline partials(n::Dual, i, j) = partials(n, i).partials[j]
@inline partials(n::Dual, i, j, k...) = partials(partials(n, i, j), k...)

@inline npartials{N}(::Dual{N}) = N
@inline npartials{N,T}(::Type{Dual{N,T}}) = N

@inline degree{T}(::T) = degree(T)
@inline degree{T}(::Type{T}) = 0
degree{N,T}(::Type{Dual{N,T}}) = 1 + degree(T)

@inline numtype(T) = Any
@inline numtype{N,T}(::Dual{N,T}) = T
@inline numtype{N,T}(::Type{Dual{N,T}}) = T

#####################
# Generic Functions #
#####################

macro ambiguous(ex)
    def = ex.head == :macrocall ? ex.args[2] : ex
    f = def.args[1].args[1].args[1]
    return quote
        $(f)(a::Dual, b::Dual) = error("npartials($(typeof(a))) != npartials($(typeof(b)))")
        if !(in($f, (isequal, ==, isless, <, <=, <)))
            $(f)(a::Dual{0}, b::Dual{0}) = Dual($(f)(value(a), value(b)))
            $(f)(a::Dual{0}, b::Dual) = $(f)(value(a), b)
            $(f)(a::Dual, b::Dual{0}) = $(f)(a, value(b))
        end
        $(esc(ex))
    end
end

Base.copy(n::Dual) = n

Base.eps(n::Dual) = eps(value(n))
Base.eps{D<:Dual}(::Type{D}) = eps(numtype(D))

Base.floor{T<:Real}(::Type{T}, n::Dual) = floor(T, value(n))
Base.floor(n::Dual) = floor(value(n))

Base.ceil{T<:Real}(::Type{T}, n::Dual) = ceil(T, value(n))
Base.ceil(n::Dual) = ceil(value(n))

Base.trunc{T<:Real}(::Type{T}, n::Dual) = trunc(T, value(n))
Base.trunc(n::Dual) = trunc(value(n))

Base.round{T<:Real}(::Type{T}, n::Dual) = round(T, value(n))
Base.round(n::Dual) = round(value(n))

Base.hash(n::Dual) = hash(value(n))
Base.hash(n::Dual, hsh::UInt64) = hash(value(n), hsh)

function Base.read{N,T}(io::IO, ::Type{Dual{N,T}})
    value = read(io, T)
    partials = read(io, Partials{N,T})
    return Dual{N,T}(value, partials)
end

function Base.write(io::IO, n::Dual)
    write(io, value(n))
    write(io, partials(n))
end

@inline Base.zero(n::Dual) = zero(typeof(n))
@inline Base.zero{N,T}(::Type{Dual{N,T}}) = Dual(zero(T), zero(Partials{N,T}))

@inline Base.one(n::Dual) = one(typeof(n))
@inline Base.one{N,T}(::Type{Dual{N,T}}) = Dual(one(T), zero(Partials{N,T}))

@inline Base.rand(n::Dual) = rand(typeof(n))
@inline Base.rand{N,T}(::Type{Dual{N,T}}) = Dual(rand(T), zero(Partials{N,T}))
@inline Base.rand(rng::AbstractRNG, n::Dual) = rand(rng, typeof(n))
@inline Base.rand{N,T}(rng::AbstractRNG, ::Type{Dual{N,T}}) = Dual(rand(rng, T), zero(Partials{N,T}))

# Predicates #
#------------#

isconstant(n::Dual) = iszero(partials(n))

@ambiguous Base.isequal{N}(a::Dual{N}, b::Dual{N}) = isequal(value(a), value(b))
@ambiguous @compat(Base.:(==)){N}(a::Dual{N}, b::Dual{N}) = value(a) == value(b)
@ambiguous Base.isless{N}(a::Dual{N}, b::Dual{N}) = value(a) < value(b)
@ambiguous @compat(Base.:<){N}(a::Dual{N}, b::Dual{N}) = isless(a, b)
@ambiguous @compat(Base.:(<=)){N}(a::Dual{N}, b::Dual{N}) = <=(value(a), value(b))

for T in (AbstractFloat, Irrational, Real)
    Base.isequal(n::Dual, x::T) = isequal(value(n), x)
    Base.isequal(x::T, n::Dual) = isequal(n, x)

    @compat(Base.:(==))(n::Dual, x::T) = (value(n) == x)
    @compat(Base.:(==))(x::T, n::Dual) = ==(n, x)

    Base.isless(n::Dual, x::T) = value(n) < x
    Base.isless(x::T, n::Dual) = x < value(n)

    @compat(Base.:<)(n::Dual, x::T) = isless(n, x)
    @compat(Base.:<)(x::T, n::Dual) = isless(x, n)

    @compat(Base.:(<=))(n::Dual, x::T) = <=(value(n), x)
    @compat(Base.:(<=))(x::T, n::Dual) = <=(x, value(n))
end

Base.isnan(n::Dual) = isnan(value(n))
Base.isfinite(n::Dual) = isfinite(value(n))
Base.isinf(n::Dual) = isinf(value(n))
Base.isreal(n::Dual) = isreal(value(n))
Base.isinteger(n::Dual) = isinteger(value(n))
Base.iseven(n::Dual) = iseven(value(n))
Base.isodd(n::Dual) = isodd(value(n))

########################
# Promotion/Conversion #
########################

Base.promote_rule{N1,N2,A<:Real,B<:Real}(D1::Type{Dual{N1,A}}, D2::Type{Dual{N2,B}}) = error("can't promote $(D1) and $(D2)")
Base.promote_rule{N,A<:Real,B<:Real}(::Type{Dual{N,A}}, ::Type{Dual{N,B}}) = Dual{N,promote_type(A, B)}
Base.promote_rule{N,T<:Real}(::Type{Dual{N,T}}, ::Type{BigFloat}) = Dual{N,promote_type(T, BigFloat)}
Base.promote_rule{N,T<:Real}(::Type{BigFloat}, ::Type{Dual{N,T}}) = Dual{N,promote_type(BigFloat, T)}
Base.promote_rule{N,T<:Real}(::Type{Dual{N,T}}, ::Type{Bool}) = Dual{N,promote_type(T, Bool)}
Base.promote_rule{N,T<:Real}(::Type{Bool}, ::Type{Dual{N,T}}) = Dual{N,promote_type(Bool, T)}
Base.promote_rule{N,T<:Real,s}(::Type{Dual{N,T}}, ::Type{Irrational{s}}) = Dual{N,promote_type(T, Irrational{s})}
Base.promote_rule{N,s,T<:Real}(::Type{Irrational{s}}, ::Type{Dual{N,T}}) = Dual{N,promote_type(Irrational{s}, T)}
Base.promote_rule{N,A<:Real,B<:Real}(::Type{Dual{N,A}}, ::Type{B}) = Dual{N,promote_type(A, B)}
Base.promote_rule{N,A<:Real,B<:Real}(::Type{A}, ::Type{Dual{N,B}}) = Dual{N,promote_type(A, B)}

Base.convert(::Type{Dual}, n::Dual) = n
Base.convert{N1,N2,T<:Real}(D::Type{Dual{N1,T}}, n::Dual{N2}) = error("can't convert $(typeof(n)) to $(D)")
Base.convert{N,T<:Real}(::Type{Dual{N,T}}, n::Dual{N}) = Dual(convert(T, value(n)), convert(Partials{N,T}, partials(n)))
Base.convert{D<:Dual}(::Type{D}, n::D) = n
Base.convert{N,T<:Real}(::Type{Dual{N,T}}, x::Real) = Dual(convert(T, x), zero(Partials{N,T}))
Base.convert(::Type{Dual}, x::Real) = Dual(x)

Base.promote_array_type{D<:Dual, A<:AbstractFloat}(F, ::Type{D}, ::Type{A}) = promote_type(D, A)
Base.promote_array_type{D<:Dual, A<:AbstractFloat, P}(F, ::Type{D}, ::Type{A}, ::Type{P}) = P
Base.promote_array_type{A<:AbstractFloat, D<:Dual}(F, ::Type{A}, ::Type{D}) = promote_type(D, A)
Base.promote_array_type{A<:AbstractFloat, D<:Dual, P}(F, ::Type{A}, ::Type{D}, ::Type{P}) = P

Base.float{N,T}(n::Dual{N,T}) = Dual{N,promote_type(T, Float16)}(n)

########
# Math #
########

# Addition/Subtraction #
#----------------------#

@ambiguous @inline @compat(Base.:+){N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) + value(n2), partials(n1) + partials(n2))
@inline @compat(Base.:+)(n::Dual, x::Real) = Dual(value(n) + x, partials(n))
@inline @compat(Base.:+)(x::Real, n::Dual) = n + x

@ambiguous @inline @compat(Base.:-){N}(n1::Dual{N}, n2::Dual{N}) = Dual(value(n1) - value(n2), partials(n1) - partials(n2))
@inline @compat(Base.:-)(n::Dual, x::Real) = Dual(value(n) - x, partials(n))
@inline @compat(Base.:-)(x::Real, n::Dual) = Dual(x - value(n), -(partials(n)))
@inline @compat(Base.:-)(n::Dual) = Dual(-(value(n)), -(partials(n)))

# Multiplication #
#----------------#

@inline @compat(Base.:*)(n::Dual, x::Bool) = x ? n : (signbit(value(n))==0 ? zero(n) : -zero(n))
@inline @compat(Base.:*)(x::Bool, n::Dual) = n * x

@ambiguous @inline function @compat(Base.:*){N}(n1::Dual{N}, n2::Dual{N})
    v1, v2 = value(n1), value(n2)
    return Dual(v1 * v2, _mul_partials(partials(n1), partials(n2), v2, v1))
end

@inline @compat(Base.:*)(n::Dual, x::Real) = Dual(value(n) * x, partials(n) * x)
@inline @compat(Base.:*)(x::Real, n::Dual) = n * x

# Division #
#----------#

@ambiguous @inline function @compat(Base.:/){N}(n1::Dual{N}, n2::Dual{N})
    v1, v2 = value(n1), value(n2)
    return Dual(v1 / v2, _div_partials(partials(n1), partials(n2), v1, v2))
end

@inline function @compat(Base.:/)(x::Real, n::Dual)
    v = value(n)
    divv = x / v
    return Dual(divv, -(divv / v) * partials(n))
end

@inline @compat(Base.:/)(n::Dual, x::Real) = Dual(value(n) / x, partials(n) / x)

# Exponentiation #
#----------------#

for f in (macroexpand(:(@compat(Base.:^))), :(NaNMath.pow))
    @eval begin
        @ambiguous @inline function ($f){N}(n1::Dual{N}, n2::Dual{N})
            if isconstant(n2)
                return $(f)(n1, value(n2))
            else
                v1, v2 = value(n1), value(n2)
                expv = ($f)(v1, v2)
                powval = v2 * ($f)(v1, v2 - 1)
                logval = expv * log(v1)
                new_partials = _mul_partials(partials(n1), partials(n2), powval, logval)
                return Dual(expv, new_partials)
            end
        end

        @inline ($f)(::Base.Irrational{:e}, n::Dual) = exp(n)
    end

    for T in (:Integer, :Rational, :Real)
        @eval begin
            @inline function ($f)(n::Dual, x::$(T))
                v = value(n)
                expv = ($f)(v, x)
                deriv = x * ($f)(v, x - 1)
                return Dual(expv, deriv * partials(n))
            end

            @inline function ($f)(x::$(T), n::Dual)
                v = value(n)
                expv = ($f)(x, v)
                deriv = expv*log(x)
                return Dual(expv, deriv * partials(n))
            end
        end
    end
end

# Unary Math Functions #
#--------------------- #

function to_nanmath(x::Expr)
    if x.head == :call
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(x.args[1]))
        return Expr(:call,funsym,[to_nanmath(z) for z in x.args[2:end]]...)
    else
        return Expr(:call,[to_nanmath(z) for z in x.args]...)
    end
end

to_nanmath(x) = x

@inline Base.conj(n::Dual) = n
@inline Base.transpose(n::Dual) = n
@inline Base.ctranspose(n::Dual) = n
@inline Base.abs(n::Dual) = signbit(value(n)) ? -n : n

for fsym in AUTO_DEFINED_UNARY_FUNCS
    v = :v
    deriv = Calculus.differentiate(:($(fsym)($v)), v)

    @eval begin
        @inline function Base.$(fsym)(n::Dual)
            $(v) = value(n)
            return Dual($(fsym)($v), $(deriv) * partials(n))
        end
    end

    # extend corresponding NaNMath methods
    if fsym in NANMATH_FUNCS
        nan_deriv = to_nanmath(deriv)
        @eval begin
            @inline function NaNMath.$(fsym)(n::Dual)
                v = value(n)
                return Dual(NaNMath.$(fsym)($v), $(nan_deriv) * partials(n))
            end
        end
    end
end

#################
# Special Cases #
#################

# Manually Optimized Functions #
#------------------------------#

@inline function Base.exp{N}(n::Dual{N})
    expv = exp(value(n))
    return Dual(expv, expv * partials(n))
end

@inline function Base.sqrt{N}(n::Dual{N})
    sqrtv = sqrt(value(n))
    deriv = 0.5 / sqrtv
    return Dual(sqrtv, deriv * partials(n))
end

# Other Functions #
#-----------------#

@inline function calc_atan2(y, x)
    z = y/x
    v= value(z)
    atan2v = atan2(value(y), value(x))
    deriv = inv(one(v) + v*v)
    return Dual(atan2v, deriv * partials(z))
end

@ambiguous @inline Base.atan2{N}(y::Dual{N}, x::Dual{N}) = calc_atan2(y, x)
@inline Base.atan2(y::Real, x::Dual) = calc_atan2(y, x)
@inline Base.atan2(y::Dual, x::Real) = calc_atan2(y, x)

###################
# Pretty Printing #
###################

function Base.show{N}(io::IO, n::Dual{N})
    print(io, "Dual(", value(n))
    for i in 1:N
        print(io, ",", partials(n, i))
    end
    print(io, ")")
end


text:
function parse_header(io)
    seekstart(io)
    rawversion = Array{UInt8}(undef, 6)
    read!(io, rawversion)
    version = String(rawversion)
    if "$version" != "FCS3.0" && version != "FCS3.1"
        warn("$version files are not guaranteed to work")
    end
    seek(io, 10)
    # start, end positions of TEXT, DATA, and ANALYSIS sections
    offsets = Array{Int64}(undef, 6)
    for i in 1:6
        # offsets are encoded as ASCII strings
        raw_str = Array{UInt8}(undef, 8)
        read!(io, raw_str)
        offsets_str = String(raw_str)
        offsets[i] = parse(Int, strip(join(offsets_str)))
    end

    # DATA offsets are larger than 99,999,999bytes
    if offsets[3] == 0 && offsets[4] == 0
        error("Reading of larger FCS not yet implemented, see issues")
    end
    offsets
end

function parse_text(io, start_text::Int, end_text::Int)
    seek(io, start_text)
    # TODO: Check for supplemental TEXT file
    raw_btext = Array{UInt8}(undef, end_text - start_text + 1)
    read!(io, raw_btext)
    raw_text = String(raw_btext)
    delimiter = raw_text[1]

    text_mappings = Dict{String, String}()
    # initialize iterator
    iter_result = iterate(raw_text)
    while iter_result !== nothing
        i, state = iter_result

        # found a new key, value pair
        if i == '$'
            # grab key and ignore escaped delimiters
            key, state = grab_word(raw_text, state, delimiter)
            # grab value and ignore escaped delimiters
            value, state = grab_word(raw_text, state, delimiter)
            # FCS keywords are case insensitive so force them uppercase
            text_mappings["\$"*uppercase(key)] = value
        end
        iter_result = iterate(raw_text, state)
    end
    text_mappings
end

function parse_data(io,
                    start_data::Int,
                    end_data::Int,
                    text_mappings::Dict{String, String})
    seek(io, start_data)
    # Add support for data types other than float
    (text_mappings["\$DATATYPE"] != "F") && error("Non float32 support not implemented yet. Please see github issues for this project.")

    flat_data = Array{Float32}(undef, (end_data - start_data + 1) ÷ 4)
    read!(io, flat_data)
    endian_func = get_endian_func(text_mappings)
    map!(endian_func, flat_data, flat_data)

    n_params = parse(Int, text_mappings["\$PAR"])

    # data should be in multiples of `n_params` for list mode
    (mod(length(flat_data), n_params) != 0) && error("FCS file is corrupt. DATA and TEXT sections don't match.")

    data = Dict{String, Vector{Float32}}()

    for i in 1:n_params
        data[text_mappings["\$P$(i)N"]] = flat_data[i:n_params:end]
    end

    FlowSample(data, text_mappings)
end

ext: jl


text:
using CompScienceMeshes, BEAST
using LinearAlgebra
using Profile
using StaticArrays

ttrc = X->ttrace(X,y)

T= tetmeshsphere(1.0,0.2)
X = nedelecc3d(T)
y = boundary(T)
@show numfunctions(X)

ϵ, μ, ω = 1.0, 1.0, 1.0; κ, η = ω * √(ϵ*μ), √(μ/ϵ)
ϵ_r = 5.0

function tau(x::SVector{U,T}) where {U,T}
    5.0 -1.0
end

χ = tau
K, I, B = VIE.singlelayer2(wavenumber=κ, tau=χ), Identity(), VIE.boundary2(wavenumber=κ, tau=χ)
E = VIE.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

@hilbertspace j
@hilbertspace k
α = ϵ_r
eq = @varform α*I[j,k]-K[j,k]+B[ttrc(j),k] == E[j]

evie = @discretise eq j∈X k∈X
u_n = solve(evie)

#postprocessing
#Φ, Θ = [0.0], range(0,stop=2π,length=100)
#pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ ]
#ffe = potential(VIE.farfield(wavenumber=κ, tau=χ), pts, u_n, X)
#ff = ffe

#Farfield
#plot(xlabel="theta")
#plot!(Θ, norm.(ff), label="far field", title="E-VIE")

using Plots

#Nearfield
Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for z in Z, y in Y]

Enear2 = BEAST.grideval(nfpoints,u_n,X)
Enear2 = reshape(Enear2,100,100)

contour(real.(getindex.(Enear2,1)))
heatmap(Z, Y, real.(getindex.(Enear2,1)))

plot!(Y[2:99],real.(getindex.(Enear2[2:99,50],1)),label="E-VIE (simplex)", linestyle=:dash, linecolor=:darkolivegreen4)


text:
using AbstractPlotting:
    NoConversion,
    convert_arguments,
    conversion_trait

@testset "Conversions" begin

    # NoConversion
    struct NoConversionTestType end
    conversion_trait(::NoConversionTestType) = NoConversion()

    let nctt = NoConversionTestType(),
        ncttt = conversion_trait(nctt)
        @test convert_arguments(ncttt, 1, 2, 3) == (1, 2, 3)
    end

end


text:
# ---
# title: 56. Merge Intervals
# id: problem56
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/merge-intervals/description/>
# hidden: true
# ---
#
# Given an array of `intervals` where `intervals[i] = [starti, endi]`, merge all
# overlapping intervals, and return _an array of the non-overlapping intervals
# that cover all the intervals in the input_.
#
#
#
# **Example 1:**
#
#
#
#     Input: intervals = [[1,3],[2,6],[8,10],[15,18]]
#     Output: [[1,6],[8,10],[15,18]]
#     Explanation: Since intervals [1,3] and [2,6] overlaps, merge them into [1,6].
#
#
# **Example 2:**
#
#
#
#     Input: intervals = [[1,4],[4,5]]
#     Output: [[1,5]]
#     Explanation: Intervals [1,4] and [4,5] are considered overlapping.
#
#
#
#
# **Constraints:**
#
#   * `1 <= intervals.length <= 104`
#   * `intervals[i].length == 2`
#   * `0 <= starti <= endi <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofDensityEstimation
include("runtests.jl")
=#
# TODO: ApproximateTwoSampleKSTest(a, b)

@testset "Test fit" begin
  simdata = CDE.Model.generate_samples(; NC=1000, NT=999, gammaC=.3, gammaT=.3,
                                       etaC=[.5, .5, 0], etaT=[.5, .4, .1],
                                       loc=[-1., 1, 2], scale=[.7, 1.3, 1],
                                       df=[15, 30, 10], skew=[-2, -5, 0])
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module MatrixDecomposition

using BenchmarkTools
using StructuredOptimization
using Images, ImageView

function set_up()
    # dataset from: http://pione.dinf.usherbrooke.ca/static/dataset/clutter/IndianTraffic3.zip
	Frames = 900
	dir = "street"
	img = load("$dir/in000000.jpg")
	n,m = size(img,1),size(img,2)
    frames = collect(1:135:Frames)
	N = length(frames)

	Y = zeros(Float64,n,m,N)
	for f in eachindex(frames)
		a = @sprintf("%6.6i",frames[f])
		img = load("$dir/in$(a).jpg")
		Y[:,:,f] .= convert(Array{Float64},Gray.(img))
		img = []
	end

	Y = reshape(Y,n*m,N)

	F = Variable(n*m,N)
	B = Variable(n*m,N)

    println("frames: $N, number of variables: $(2*N*m*n)")

	R, lambda = 1, 3e-2
	return B, F, Y, R, lambda, n, m, N
end

function run_demo()
	slv = PANOC(tol = 1e-4, verbose = 2)
	setup = set_up()
	solve_problem!(slv,setup...)
	return setup
end

function solve_problem!(slv, B, F, Y, R, lambda, n, m, N)
	it, = @minimize ls(B+F-Y)+lambda*norm(F,1) st rank(B) <= R with slv
	return it
end

function benchmark(;verb = 0, samples = 5, seconds = 100, tol = 1e-4)

	suite = BenchmarkGroup()

	solvers = [
             #  "ZeroFPR",
               "PANOC",
               "PG"
              ]
	slv_opt = [
              # "(verbose = $verb, tol = $tol)",
               "(verbose = $verb, tol = $tol)",
               "(verbose = $verb, tol = $tol)"
              ]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup = set_up()
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable(it = solve_problem!(solver, setup...),
                       setup = (
                                it = 0;
                                setup = deepcopy($setup);
                                solver = deepcopy($solver) ),
                       teardown = (
                                   $its[$solvers[$i]] = it;
                                  ),
                       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("MatrixDecomposition its")
	println(its)
	return results
end

function show_results(B, F, Y, R, lambda, n, m, N)
	F = F.x
	B = B.x

	F[F .!=0] = F[F.!=0]+reshape(B,size(F))[F.!=0]
	F[F.== 0] = 1. #put white in null pixels
	F[F.>1] .= 1; F[F.<0.] .= 0.
	B[B.>1] .= 1; B[B.<0.] .= 0.

	#convert back to images
	Y = reshape(Y,n,m,N)
	B = reshape(B,n,m,N)
	F = reshape(F,n,m,N)
	#
	####videos
	#imshow(Y)
	#imshow(F)
	#imshow(B)
	#
	##images
	idx = round.(Int,linspace(1,N,4))
	imshow(hcat([Y[:,:,i] for i in idx]...))
	imshow(hcat([F[:,:,i] for i in idx]...))
	imshow(hcat([B[:,:,i] for i in idx]...))

    save_stuff = false
    if save_stuff == true
        for (i,idxi) in enumerate(idx)
            save("Y$i.jpg", Y[:,:,idxi])
            save("F$i.jpg", F[:,:,idxi])
            save("B$i.jpg", B[:,:,idxi])
        end
    end

end

end


text:
N = 30
ph_constrained = discrete_phs(N, 0, 1);
set_constraint!(ph_constrained, ph_constrained.B[:,1], ph_constrained.D[1:1,1:1]);
eigval, eigvec = eigen(ph_constrained)
num_freq = frequencies(eigval)
exact_freq = pi*collect(-0:1:N-1);
err = exact_freq - num_freq;

@test abs(err[2]) < 1e-13

ph_c_elim = constraint_elimination(ph_constrained)
eigval, eigvec = eigen(ph_c_elim);
num_freq = frequencies(eigval)
exact_freq = pi*collect(0:1:N-1);
err = exact_freq - num_freq;

@test abs(err[2]) < 1e-13

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
@testset "marketdata_v2" begin
  channel = Channel{Dict}(5)
  println("opening socket to Gemini v2/marketdata...")
  @async marketdata_v2(
    false,
    channel,
    ["candles_1m"],
    ["BTCUSD"]
  )
  for _ in 1:5
    data = take!(channel)
    println("received data...")
  end
  close(channel)
end

max_stars_repo_name: rory-linehan/Gemini

max_issues_repo_name: rory-linehan/Gemini

max_forks_repo_name: rory-linehan/Gemini


text:
function err_fun_generate()
    counter::Integer = 0
    fun = (x::Vector{Float64}) -> begin
        counter += 1
        if counter > 5
            throw(ErrorException("Function can be called only 5 times."))
        end
        5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2
    end
end

@testset "loss error in :CICO_ONE_PASS" begin
    err_fun = err_fun_generate()
    res0 = get_endpoint(
        [3., 4.],
        1,
        err_fun,
        :CICO_ONE_PASS;
        loss_crit = 9.,
        silent = true
    )
    @test res0.value === nothing
    @test res0.profilePoints == ProfilePoint[]
    @test res0.status == :LOSS_ERROR_STOP
    @test res0.counter < 5
    @test typeof(res0.supreme) == Float64
end

@testset "test profile error" begin
    err_func = err_fun_generate()
    prof = LikelihoodProfiler._profile(
        [3., 4.],
        1,
        err_func
    )
    res0 = prof(5.)
    @test res0.ret == :FORCED_STOP
    @test res0.counter == 5
end

@testset "loss error in :LIN_EXTRAPOL" begin
    err_fun = err_fun_generate()
    res0 = get_endpoint(
        [3., 4.],
        1,
        err_fun,
        :LIN_EXTRAPOL;
        loss_crit = 9.,
        silent = true
    )
    @test res0.value === nothing
    @test length(res0.profilePoints) == 1
    @test res0.status == :LOSS_ERROR_STOP
    @test res0.counter < 5
    @test typeof(res0.supreme) == Float64
    pp = res0.profilePoints[1]
    @test pp.ret == :FORCED_STOP
end

@testset "loss error in :QUADR_EXTRAPOL" begin
    err_fun = err_fun_generate()
    res0 = get_endpoint(
        [3., 4.],
        1,
        err_fun,
        :QUADR_EXTRAPOL;
        loss_crit = 9.,
        silent = true
    )
    @test res0.value === nothing
    @test length(res0.profilePoints) == 1
    @test res0.status == :LOSS_ERROR_STOP
    @test res0.counter < 5
    @test typeof(res0.supreme) == Float64
    pp = res0.profilePoints[1]
    @test pp.ret == :FORCED_STOP
end


text:
using Pkg

ENV["PYTHON"] = ""
ENV["CONDA_JL_HOME"] = "D:/Users/Aure/Anaconda3/envs/julia" # path to miniconda3 directory

Pkg.add("PyCall")
Pkg.build("PyCall")

using PyCall

math = pyimport("math")

math.sin(3)

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# Domain and action grounding

include("action.jl")
include("domain.jl")
include("interface.jl")
include("available.jl")
include("execute.jl")
include("transition.jl")
include("utils.jl")


text:
module testConstraints

using Base.Test
using DataFrames
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running constraints tests\n")

##############################
## chkEqualsOne for vectors ##
##############################

## test valid portfolio weights
simVals = rand(8, 1)
simVals = simVals./sum(simVals)

AssetMgmt.chkEqualsOne(simVals)
AssetMgmt.chkEqualsOne(simVals[:])

## test invalid portfolio weights
simVals[2] = simVals[2] - 1e-8

@test_throws AssetMgmt.chkEqualsOne(simVals)
@test_throws AssetMgmt.chkEqualsOne(simVals[:])

##########################
## makeWeights function ##
##########################

## test valid investment weights
simVals = rand(1, 10)
wgts = AssetMgmt.makeWeights(simVals)
AssetMgmt.chkEqualsOne(wgts)

###############################
## chkEqualsOne for matrices ##
###############################

simVals = rand(10, 5)
simVals = AssetMgmt.makeWeights(simVals)

## chkEquals meanwhile defined for matrices
AssetMgmt.chkEqualsOne(simVals)

tm = Timematr(simVals)
AssetMgmt.chkEqualsOne(tm)

## test invalid input
simVals[1, 1] = simVals[1, 1] + 1e-8

tm = Timematr(simVals)
@test_throws AssetMgmt.chkEqualsOne(tm)

end


text:
# ---
# title: 60. Permutation Sequence
# id: problem60
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Backtracking
# link: <https://leetcode.com/problems/permutation-sequence/description/>
# hidden: true
# ---
#
# The set `[1, 2, 3, ..., n]` contains a total of `n!` unique permutations.
#
# By listing and labeling all of the permutations in order, we get the following
# sequence for `n = 3`:
#
#   1. `"123"`
#   2. `"132"`
#   3. `"213"`
#   4. `"231"`
#   5. `"312"`
#   6. `"321"`
#
# Given `n` and `k`, return the `kth` permutation sequence.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 3
#     Output: "213"
#
#
# **Example 2:**
#
#
#
#     Input: n = 4, k = 9
#     Output: "2314"
#
#
# **Example 3:**
#
#
#
#     Input: n = 3, k = 1
#     Output: "123"
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 9`
#   * `1 <= k <= n!`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Rheologies

    using Reexport

    #@reexport using Tensors
    @reexport using JuAFEM
    const J = JuAFEM
    import JuAFEM: QuadratureRule

    #using Parameters
    using InteractiveUtils: subtypes
    using TimerOutputs
    using Printf
    using SparseArrays
    using BlockArrays
    using StaticArrays
    using StructArrays
    using LinearAlgebra

    using ForwardDiff
    using DiffResults
    #import MUMPSjInv
    #import Preconditioners

    #include("convenience_functions.jl")
    include("types/types.jl")
    include("physical_functions.jl")
    include("tangent_operators.jl")
    include("interface.jl")
    include("main.jl")
    include("linear_solve.jl")
    include("nonlinear_iterations.jl")
    include("assemble/assemble_global.jl")
    include("GmshParser.jl")

    # exports #
    # types
    export Model
    export MaterialState, BasicMaterialState, PlasticMaterialState, DamagedPlasticMaterialState
    export AbstractSolver, AbstractLinearSolver, AbstractNonLinearSolver, BackslashSolver, MUMPS, NewtonRaphson, ConjugateGradients,  MINRES, ILU, GMRESIterativeSolvers
    export Rheology, BoundaryConditions, BodyForces, Variables, Clock, PrimitiveVariables
    export Damage
    export Viscosity
    export Elasticity
    export Plasticity, DruckerPrager, VonMises
    export Damage, BRSDamage
    export MixedOutputWriter, VTKOutputWriter, JLD2OutputWriter

    # functions
    export solve, setup_model, create_material_properties

end # module


text:
######################################################################
# Unser erstes Programm
# Twitter: twitter.com/programmonline
# Repl.it: repl.it/@programmonline/RespectfulPastelElectricity#main.jl
######################################################################

println( "Hallo, Welt" )

######################################################################
# (c) by programmierer.online
######################################################################

max_stars_repo_name: programmiereronline/Julia-Kurs

max_issues_repo_name: programmiereronline/Julia-Kurs

max_forks_repo_name: programmiereronline/Julia-Kurs


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V,(VV,EV,FV) = Lar.simplex(2, true)
triangle = (V,EV,FV)
model = Lar.Struct([ triangle, Lar.t(.15,.15), Lar.s(.5,.5), triangle ])
V,EV,FV = Lar.struct2lar(model)
VV = [[k] for k=1:size(V,2)]
GL.VIEW(GL.numbering(.4)((V,[VV, EV, FV]),GL.COLORS[1],0.2));

triangles = Lar.triangulate2d(V, EV)
GL.VIEW([ GL.GLGrid(V,triangles) ])

# generate edges
ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
# remove duplicated edges from triangulation
nodupev = collect(Set(ev))
GL.VIEW(GL.numbering(.4)((V,[VV, nodupev])));


text:
module Rplots

using RCall
const RSOURCE_DIR = "$(pwd())/src/R"
const RPLOT_SOURCE = "$(RSOURCE_DIR)/plots.R"
rplotlib() = R"source($RPLOT_SOURCE, chdir=TRUE)"

end # Rplots

#= TEST
using CytofRepFAM
CytofRepFAM.Rplots.RSOURCE_DIR
CytofRepFAM.Rplots.rplotlib()
=#


text:
# This file simulates an opamp integrator circuit.

using Causal
using Plots

freq = 5e3
T = 1 / freq
r = 10e3
c = 10e-9
τ = r * c

t0, dt, tf = 0, T/1000, 5T

@defmodel model begin
    @nodes begin
        gen = SquarewaveGenerator(high=0.5, low=-0.5, period=T)
        ds = ContinuousLinearSystem(A=fill(0., 1, 1), B=fill(1/τ, 1, 1), C=fill(-1., 1, 1), state=zeros(1))
        writerin = Writer()
        writerout = Writer()
    end
    @branches begin
        gen     =>      ds
        gen     =>      writerin
        ds      =>      writerout
    end
end

sim = simulate!(model, t0, dt, tf)

t, u = read(getnode(model, :writerin).component)
t, y = read(getnode(model, :writerout).component)

p1 = plot(t, u)
    plot!(t, y)
display(p1)


text:
function LazyLoadedCelestialBodies(arg0::DataProvidersManager, arg1::TimeScales, arg2::Frame)
    return LazyLoadedCelestialBodies((DataProvidersManager, TimeScales, Frame), arg0, arg1, arg2)
end

function add_celestial_body_loader(obj::LazyLoadedCelestialBodies, arg0::JString, arg1::CelestialBodyLoader)
    return jcall(obj, "addCelestialBodyLoader", void, (JString, CelestialBodyLoader), arg0, arg1)
end

function add_default_celestial_body_loader(obj::LazyLoadedCelestialBodies, arg0::JString)
    return jcall(obj, "addDefaultCelestialBodyLoader", void, (JString,), arg0)
end

function add_default_celestial_body_loader(obj::LazyLoadedCelestialBodies, arg0::JString, arg1::JString)
    return jcall(obj, "addDefaultCelestialBodyLoader", void, (JString, JString), arg0, arg1)
end

function clear_celestial_body_cache(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "clearCelestialBodyCache", void, ())
end

function clear_celestial_body_cache(obj::LazyLoadedCelestialBodies, arg0::JString)
    return jcall(obj, "clearCelestialBodyCache", void, (JString,), arg0)
end

function clear_celestial_body_loaders(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "clearCelestialBodyLoaders", void, ())
end

function clear_celestial_body_loaders(obj::LazyLoadedCelestialBodies, arg0::JString)
    return jcall(obj, "clearCelestialBodyLoaders", void, (JString,), arg0)
end

function get_body(obj::LazyLoadedCelestialBodies, arg0::JString)
    return jcall(obj, "getBody", CelestialBody, (JString,), arg0)
end

function get_earth(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getEarth", CelestialBody, ())
end

function get_earth_moon_barycenter(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getEarthMoonBarycenter", CelestialBody, ())
end

function get_jupiter(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getJupiter", CelestialBody, ())
end

function get_mars(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getMars", CelestialBody, ())
end

function get_mercury(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getMercury", CelestialBody, ())
end

function get_moon(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getMoon", CelestialBody, ())
end

function get_neptune(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getNeptune", CelestialBody, ())
end

function get_pluto(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getPluto", CelestialBody, ())
end

function get_saturn(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getSaturn", CelestialBody, ())
end

function get_solar_system_barycenter(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getSolarSystemBarycenter", CelestialBody, ())
end

function get_sun(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getSun", CelestialBody, ())
end

function get_uranus(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getUranus", CelestialBody, ())
end

function get_venus(obj::LazyLoadedCelestialBodies)
    return jcall(obj, "getVenus", CelestialBody, ())
end


text:
# This file was generated, do not modify it. # hide
@show y[1:3]
scitype_union(y)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
# Schamber's Vector quant algorithm - It foregoes the a weighted least squares fit for a simple
# homoskedastic fit.  In doing this it can precompute a fitting matrix which requires nothing
# more than a single matrix multiplication to perform the fit.  This makes this mechanism
# extremely quick. This makes it ideal for processing in real-time or HyperSpectrum objects.

struct _VQRefData{T <: AbstractFloat}
    label::ReferenceLabel
    roi::UnitRange{Int}
    charonly::Vector{T}
    sumchar::T
    scale::T
end

function Base.show(io::IO, vqr::_VQRefData)
    print(io,"$(vqr.label)[$(vqr.roi)]")
end

struct VectorQuant{T <: AbstractFloat}
    # Vector(label[1], roi[2], charonly[3], sum(charonly)[4], scale[5])
    references::Vector{_VQRefData{T}}
    vectors::Matrix{T}

    """
        VectorQuant(frefs::Vector{FilteredReference}, filt::TopHatFilter)

    Constructs a structure used to perform accelerated filtered spectrum fits based on the specified
    collection of `FilteredReference`(s), and a `TopHatFilter`.
    """
    function VectorQuant(frefs::Vector{<:FilteredReference}, filt::TopHatFilter{T}) where { T <: AbstractFloat}
        refs = map(frefs) do fref
            _VQRefData(fref.label, fref.roi, fref.charonly, sum(fref.charonly), fref.scale)
        end
        x = zeros(T, (length(filt.filters), length(frefs)))
        for (c, fref) in enumerate(frefs)
            x[fref.ffroi, c] = fref.filtered
        end
        # ((ch × ne)T * (ch × ne))^(-1) * (ch × ne) * (ch × ch) => (ne × ch)
        xTxIxf =
            pinv(transpose(x) * x) *
            transpose(x) *
            filterdata(filt, 1:length(filt.filters))
        return new{T}(refs, xTxIxf)
    end

    VectorQuant(ffrs::FilterFitPacket{S, T}) where { S <: Detector, T <: AbstractFloat } = #
        VectorQuant(ffrs.references, ffrs.filter)
end

NeXLCore.minproperties(::VectorQuant) = (:BeamEnergy, :TakeOffAngle, :)

Base.show(io::IO, vq::VectorQuant) = print(
    io,
    "VectorQuant[\n" * join(map(r -> "\t$r", vq.references), ",\n") * "\n]",
)
"""
    fit_spectrum(
        hs::Spectrum|HyperSpectrum,
        vq::VectorQuant{S <: Detector, T <: AbstractFloat},
        zero = x -> max(Base.zero(T), x)
    )

Fit the spectrum or hyper-spectrum using the vector-quant algorithm. The function `zero` is
applied to the resultant k-ratios before they are returned.  The default simply sets negative
k-ratios to 0.0.  `zero=identity` would leave the negative k-ratios as such.
"""
function fit_spectrum(
    spec::Spectrum,
    vq::VectorQuant{T},
    zero = x -> max(Base.zero(T), x),
) where { T <: AbstractFloat }
    raw = counts(spec, 1:size(vq.vectors, 2), T, true)
    krs = zero.(vq.vectors * raw)
    spsc = T(dose(spec))
    residual = copy(raw)
    for (i, vqr) in enumerate(vq.references)
        residual[vqr.roi] -= krs[i] * vqr.charonly
    end
    peakback = Dict{ReferenceLabel,NTuple{3,T}}()
    dkrs = zeros(T, length(vq.references))
    for (i, vqr) in enumerate(vq.references)
        ii, bb = krs[i] * vqr.sumchar, sum(residual[vqr.roi])
        peakback[vqr.label] = (ii, bb, bb / spsc)
        dkrs[i] = sqrt(max(Base.zero(T), ii + bb)) / vqr.sumchar
    end
    kratios = uvs(
        map(ref -> ref.label, vq.references), #
        map(i -> krs[i] / (vq.references[i].scale * spsc), eachindex(krs)), #
        map(i -> (dkrs[i] / (vq.references[i].scale * spsc))^2, eachindex(krs)),
    )
    return FilterFitResult{T}(
        UnknownLabel(spec),
        kratios,
        1:length(raw),
        raw,
        residual,
        peakback
    )
end
function fit_spectrum(
    hs::HyperSpectrum,
    vq::VectorQuant{T},
    zero = x -> max(Base.zero(T), x),
)::Array{KRatios} where { T <: AbstractFloat }
    krs = zeros(T, length(vq.references), size(hs)...)
    vecs = vq.vectors[:, 1:depth(hs)]
    foreach(i->vecs[i,:]/=vq.references[i].scale, eachindex(vq.references))
    data = T.(counts(hs))  # One large allocation over many smaller????
    # In my testing, ThreadsX produces about a 2.7x speedup for 4 threads on 4 cores
    ThreadsX.foreach(CartesianIndices(hs)) do ci
        krs[:, ci] .= (vecs * view(data, :, ci)) / T(dose(hs, ci))
    end
    # ensure positive...
    map!(zero, krs, krs)
    res = map(filter(ii -> vq.references[ii].label isa CharXRayLabel, eachindex(vq.references))) do i
        lbl = vq.references[i].label
        rprops = properties(spectrum(lbl))
        KRatios(xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :])
    end
    return res
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
###########################################
# From last lecture
###########################################

abstract type Step end

struct GD <: Step
    α::Real
end

optim_step(s::GD, f, g, x) = -s.α*g(x)

function optim(f, g, x, s::Step; max_iter=100)
    for i in 1:max_iter
        x += optim_step(s, f, g, x)
    end
    return x
end

###########################################
# Linear regression
###########################################

# EXTRA

σ(z) = 1/(1+exp(-z))
plot(-10:0.01:10, σ)

# Load data

using Plots
using StatsPlots
using RDatasets
using LinearAlgebra
using Query
using Statistics

iris = dataset("datasets", "iris")

iris[1:5,:]

# Exercise

y = iris.PetalWidth
X = hcat(iris.PetalLength, ones(length(y)))

@df iris scatter(
    :PetalLength,
    :PetalWidth;
    label="",
    xlabel = "Petal length",
    ylabel = "Petal width"
)

# Exercise

w = (X'*X)^(-1)*X'*y
w = (X'*X)^(-1)*(X'*y)
w = X'*X \ X'*y

g(w) = X'*(X*w-y)
w2 = optim([], g, zeros(size(X,2)), GD(1e-4); max_iter=10000)

norm(w-w2)

# EXTRA

w = (X'*X)^(-1)*X'*y
w = (X'*X)^(-1)*(X'*y)
w = X'*X \ X'*y

A = randn(1000,1000)
b = randn(1000)

@time (A*A)*b;
@time A*(A*b);

# Exercise

f_pred(x::Real, w) = ([x 1]*w)[1]
f_pred(x::Real, w) = w[1]*x + w[2]*1

x_lims = extrema(iris.PetalLength) .+ [-0.1, 0.1]

@df iris scatter(
    :PetalLength,
    :PetalWidth;
    xlabel = "Petal length",
    ylabel = "Petal width",
    label = "",
    legend = :topleft,
)

plot!(x_lims, x -> f_pred(x,w); label = "Prediction", line = (:black,3))

###########################################
# Logistic regression
###########################################

# Exercise

iris_reduced = @from i in iris begin
    @where i.Species != "setosa"
    @select {
        i.PetalLength,
        i.PetalWidth,
        intercept = 1,
        i.Species,
        label = i.Species == "virginica",
    }
    @collect DataFrame
end

# EXTRA

iris_reduced2 = iris[iris.Species .!= "setosa", :]
iris_reduced2 = iris_reduced2[:,[3;4;5]]

insertcols!(iris_reduced2, 3, :intercept => 1)
insertcols!(iris_reduced2, 5, :label => iris_reduced2.Species .== "virginica")

isequal(iris_reduced, iris_reduced2)

# Extract data

X = Matrix(iris_reduced[:, 1:3])
y = iris_reduced.label

# Exercise

@df iris_reduced scatter(
    :PetalLength,
    :PetalWidth;
    group = :Species,
    xlabel = "Petal length",
    ylabel = "Petal width",
    legend = :topleft,
)

# EXTRA

ii = iris_reduced.Species .== "virginica"

scatter(iris_reduced.PetalLength[ii,:], iris_reduced.PetalWidth[ii,:];
    xlabel = "Petal length",
    ylabel = "Petal width",
    legend = :topleft,
    label = "virginica"
)

scatter!(iris_reduced.PetalLength[.!ii,:], iris_reduced.PetalWidth[.!ii,:];
    label = "versicolor"
)

# Exercise

function log_reg(X, y, w; max_iter=100, tol=1e-6)
    X_mult = [row*row' for row in eachrow(X)]
    for i in 1:max_iter
        y_hat = σ.(X*w)
        grad = X'*(y_hat.-y) / size(X,1)
        hess = y_hat.*(1 .-y_hat).*X_mult |> mean
        w -= hess \ grad
    end
    return w
end

w = log_reg(X, y, zeros(size(X,2)))

# EXTRA

function g(w, X, y)
    n = size(X,1)
    y_hat = σ.(X*w)
    return X'*(y_hat - y) / n
end

optim([], w -> g(w, X, y), [0;0;0], GD(5e-1); max_iter=1000000)

# Plot the separating hyperplane

separ(x, w) = (-w[3]-w[1]*x)/w[2]

xlims = extrema(iris_reduced.PetalLength) .+ [-0.1, 0.1]
ylims = extrema(iris_reduced.PetalWidth) .+ [-0.1, 0.1]

@df iris_reduced scatter(
    :PetalLength,
    :PetalWidth;
    group = :Species,
    xlabel = "Petal length",
    ylabel = "Petal width",
    legend = :topleft,
    xlims,
    ylims,
)

plot!(xlims, x -> separ(x,w); label = "Separation", line = (:black,3))

# Check optimality

y_hat = σ.(X*w)
grad = X'*(y_hat.-y) / size(X,1)
norm(grad)

# Exercise

pred = y_hat .>= 0.5
"Correct number of predictions: " * string(sum(pred .== y))
"Wrong   number of predictions: " * string(sum(pred .!= y))

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
using Oceananigans.Grids: metrics_precomputed, on_architecture, pop_flat_elements
import Oceananigans.Grids: architecture, size, new_data, halo_size
import Oceananigans.Grids: with_halo, on_architecture

struct MultiRegionGrid{FT, TX, TY, TZ, P, G, D, Arch} <: AbstractMultiRegionGrid{FT, TX, TY, TZ, Arch}
    architecture :: Arch
    partition :: P
    region_grids :: G
    devices :: D

    function MultiRegionGrid{FT, TX, TY, TZ}(arch::A, partition::P,
                                             region_grids::G,
                                             devices::D) where {FT, TX, TY, TZ, P, G, D, A}

        return new{FT, TX, TY, TZ, P, G, D, A}(arch, partition, region_grids, devices)
    end
end

@inline isregional(mrg::MultiRegionGrid)        = true
@inline getdevice(mrg::MultiRegionGrid, i)      = getdevice(mrg.region_grids, i)
@inline switch_device!(mrg::MultiRegionGrid, i) = switch_device!(getdevice(mrg, i))
@inline devices(mrg::MultiRegionGrid)           = devices(mrg.region_grids)
@inline sync_all_devices!(mrg::MultiRegionGrid) = sync_all_devices!(devices(mrg))

@inline  getregion(mrg::MultiRegionGrid, r) = _getregion(mrg.region_grids, r)
@inline _getregion(mrg::MultiRegionGrid, r) =  getregion(mrg.region_grids, r)

@inline Base.length(mrg::MultiRegionGrid)   = Base.length(mrg.region_grids)

const ImmersedMultiRegionGrid = MultiRegionGrid{FT, TX, TY, TZ, P, <:MultiRegionObject{<:Tuple{Vararg{<:ImmersedBoundaryGrid}}}} where {FT, TX, TY, TZ, P}

function MultiRegionGrid(global_grid; partition = XPartition(2), devices = nothing)

    if length(partition) == 1
        return global_grid
    end

    @warn "MultiRegion functionalities are experimental: help the development by reporting bugs or non-implemented features!"

    arch = architecture(global_grid)

    devices = validate_devices(partition, arch, devices)
    devices = assign_devices(partition, devices)

    global_grid  = on_architecture(CPU(), global_grid)
    local_size   = MultiRegionObject(partition_size(partition, global_grid), devices)
    local_extent = MultiRegionObject(partition_extent(partition, global_grid), devices)
    local_topo   = MultiRegionObject(partition_topology(partition, global_grid), devices)

    global_topo  = topology(global_grid)

    FT   = eltype(global_grid)

    args = (Reference(global_grid),
            Reference(arch),
            local_topo,
            local_size,
            local_extent,
            Reference(partition),
            Iterate(1:length(partition)))

    region_grids = construct_regionally(construct_grid, args...)

    return MultiRegionGrid{FT, global_topo[1], global_topo[2], global_topo[3]}(arch, partition, region_grids, devices)
end

function construct_grid(grid::RectilinearGrid, child_arch, topo, size, extent, args...)
    halo = halo_size(grid)
    size = pop_flat_elements(size, topo)
    halo = pop_flat_elements(halo, topo)
    FT   = eltype(grid)
    return RectilinearGrid(child_arch, FT; size = size, halo = halo, topology = topo, extent...)
end

function construct_grid(grid::LatitudeLongitudeGrid, child_arch, topo, size, extent, args...)
    halo = halo_size(grid)
    FT   = eltype(grid)
    lon, lat, z = extent
    return LatitudeLongitudeGrid(child_arch, FT;
                                 size = size, halo = halo, radius = grid.radius,
                                 latitude = lat, longitude = lon, z = z, topology = topo,
                                 precompute_metrics = metrics_precomputed(grid))
end

function construct_grid(ibg::ImmersedBoundaryGrid, child_arch, topo, local_size, extent, partition, region)
    boundary = partition_immersed_boundary(ibg.immersed_boundary, partition, local_size, region, child_arch)
    return ImmersedBoundaryGrid(construct_grid(ibg.underlying_grid, child_arch, topo, local_size, extent), boundary)
end

partition_immersed_boundary(b, args...) =
    getname(b)(partition_global_array(getproperty(b, propertynames(b)[1]), args...))

function reconstruct_global_grid(mrg)
    size    = reconstruct_size(mrg, mrg.partition)
    extent  = reconstruct_extent(mrg, mrg.partition)
    topo    = topology(mrg)
    switch_device!(mrg.devices[1])
    return construct_grid(mrg.region_grids[1], architecture(mrg), topo, size, extent)
end

function reconstruct_global_grid(mrg::ImmersedMultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    underlying_mrg = MultiRegionGrid{FT, TX, TY, TZ}(architecture(mrg),
                                                     mrg.partition,
                                                     construct_regionally(getproperty, mrg, :underlying_grid),
                                                     mrg.devices)

    global_grid     = on_architecture(CPU(), reconstruct_global_grid(underlying_mrg))
    cpu_mrg         = on_architecture(CPU(), mrg)
    local_boundary  = construct_regionally(getproperty, cpu_mrg, :immersed_boundary)
    local_array     = construct_regionally(getproperty, local_boundary, propertynames(local_boundary[1])[1])
    local_array     = construct_regionally(getinterior, local_array, mrg)
    global_boundary = getname(local_boundary[1])(reconstruct_global_array(local_array, mrg.partition, architecture(mrg)))
    return on_architecture(architecture(mrg), ImmersedBoundaryGrid(global_grid, global_boundary))
end

getinterior(array::AbstractArray{T, 2}, grid) where T = array[1:grid.Nx, 1:grid.Ny]
getinterior(array::AbstractArray{T, 3}, grid) where T = array[1:grid.Nx, 1:grid.Ny, 1:grid.Nz]
getinterior(func::Function, grid) = func

multi_region_object_from_array(a::AbstractArray, grid) = arch_array(architecture(grid), a)

function multi_region_object_from_array(a::AbstractArray, mrg::MultiRegionGrid)
    local_size = construct_regionally(size, mrg)
    arch = architecture(mrg)
    a    = arch_array(CPU(), a)
    ma   = construct_regionally(partition_global_array, a, mrg.partition, local_size, Iterate(1:length(mrg)), arch)
    return ma
end

####
#### Utilitites for MultiRegionGrid
####

new_data(FT::DataType, mrg::MultiRegionGrid, args...) = construct_regionally(new_data, FT, mrg, args...)

function with_halo(new_halo, mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    new_grids = construct_regionally(with_halo, new_halo, mrg)
    return MultiRegionGrid{FT, TX, TY, TZ}(mrg.architecture, mrg.partition, new_grids, mrg.devices)
end

function on_architecture(::CPU, mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    new_grids = construct_regionally(on_architecture, CPU(), mrg)
    devices   = Tuple(CPU() for i in 1:length(mrg))
    return MultiRegionGrid{FT, TX, TY, TZ}(CPU(), mrg.partition, new_grids, devices)
end

Base.summary(mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    "MultiRegionGrid{$FT, $TX, $TY, $TZ} with $(summary(mrg.partition)) on $(string(typeof(mrg.region_grids[1]).name.wrapper))"

Base.show(io::IO, mrg::MultiRegionGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    print(io, "MultiRegionGrid{$FT, $TX, $TY, $TZ} partitioned on $(architecture(mrg)): \n",
              "├── grids: $(summary(mrg.region_grids[1])) \n",
              "├── partitioning: $(summary(mrg.partition)) \n",
              "└── devices: $(devices(mrg))")

function Base.:(==)(mrg1::MultiRegionGrid, mrg2::MultiRegionGrid)
    #check if grids are of the same type
    vals = construct_regionally(Base.:(==), mrg1, mrg2)
    return all(vals.regions)
end

# These are not used in the code, should we remove them?

grids(mrg::MultiRegionGrid) = mrg.region_grids

getmultiproperty(mrg::MultiRegionGrid, x::Symbol) = construct_regionally(Base.getproperty, grids(mrg), x)

const MRG = MultiRegionGrid

@inline Base.getproperty(ibg::MRG, property::Symbol) = get_multi_property(ibg, Val(property))
@inline get_multi_property(ibg::MRG, ::Val{property}) where property = getproperty(getindex(getfield(ibg, :region_grids), 1), property)
@inline get_multi_property(ibg::MRG, ::Val{:partition}) = getfield(ibg, :partition)
@inline get_multi_property(ibg::MRG, ::Val{:region_grids}) = getfield(ibg, :region_grids)
@inline get_multi_property(ibg::MRG, ::Val{:devices}) = getfield(ibg, :devices)