text:
using TikzGraphs
using Test

@assert success(`lualatex -v`)
using NBInclude
@nbinclude joinpath(dirname(@__FILE__), "..", "doc", "TikzGraphs.ipynb")


text:
using LinearAlgebra

N = 40
N_t = 300
N_f = 150
#s = tf("s");

""" `toeplitz2(vc, N = length(vc))`
Create toeplitxz matrix by stacking vc as columns, to width N, keeping height at length(vc) """
function toeplitz2(vc, N = size(vc,1))
    M = vc[:,:]
    Ncol = size(vc,2)
    Nrow = size(vc,1)
    for i = 2:N
        M = [M [zeros(i-1, Ncol); vc[1:Nrow-i+1, :]]]
    end
    return M
end

function toeplitz(vc::AbstractVector{T}, vr) where T
    N = length(vc)
    @assert N == length(vr)
    @assert vc[1] == vr[1]
    M = Array{T,2}(undef, N, N)
    # Subdiagonal
    for i = 1:N # Each column
        M[i:N,i] = vc[1:N-i+1]
    end
    for i = 1:N # Each row
        M[i, i:N] = vr[1:N-i+1]
    end
    return M
end

""" Compute A,b for the eqpression
    A*q+b corresponding to S*v=(I-PQ)v """
function Saff_t(TP,v,N)
    A = - T_P*toeplitz2(v,N)
    b = copy(v)
    return A, b
end

""" Compute A,b for the eqpression
    A*q+b corresponding to PS*v=P(I-PQ)v """
function PSaff_t(TP,v,N)
    AQ,bQ = CSaff_t(TP,v,N)
    A = -TP*TP*AQ
    b = TP*(v - TP*bQ)
    return A, b
end

""" Compute A,b for the eqpression
    A*q+b corresponding to CS*v=Qv """
function CSaff_t(TP,v,N)
    A = toeplitz2(v,N)
    b = zero(v)
    return A, b
end

""" Compute A,b for the eqpression
    A*q+b corresponding to CS*v=Qv """
function PCSaff_t(TP,v,N)
    A = TP*toeplitz2(v,N)
    b = zero(v)
    return A, b
end

function Saff_fr(P_fr, Qf)
    # Q_fr = G_q_fr*b
    # S_fr = 1 - P_fr.*Q_fr
    A = -P_fr.*Qf
    b = ones(eltype(Qf), size(P_fr,1))
    return A,b
end

""" T_P = SysP(g)
Given impulse response of system p, generate system (toeplitz) matrix TP """
function SysP(g)
    T_P = toeplitz(g, [g[1]; zeros(length(g)-1)]);
end

"""
    G_q_fr = Qfreq(Ω, N)
    Generate matrix for calculating frequency respose (Q_fr) of Q on grid Ω, such that
    Q_fr = G_q_fr*q
    where q is the impulse response or order N, and Ω∈(0,π)
"""
function Qfreq(Ω, N)
    exp.((-im.*Ω)*(0:N-1)')
end

using ControlSystems
#P = c2d(1/(20*s + 1) * exp(-5*s), 1);
z_d = tf("z", 1)
P_sysd = (1/20)/(z_d-0.95)/z_d^5
P_d(z) = (1/20)/(z-0.95)/z^5
#P(z) = (1/20)/(z-0.95)*z^(-5)

g = impulse(P_sysd, N_t-1)[1][:]

T_P = SysP(g)

Ω = 2*pi*(10 .^range(-3, log10(0.5), length=N_f))

#P_fr = freqresp(P_sysd, Ω)[:]
P_fr = P_d.(exp.(im.*Ω))
G_q_fr = Qfreq(Ω, N)

e1 = [1; zeros(N_t - 1)];
d1 = ones(N_t);

PSA, PSb = PSaff_t(T_P, d1, N)

SA_fr, Sb_fr = Saff_fr(P_fr, G_q_fr)

using Convex

################# Multiple affine and cones in Convex:

b = Variable(N)
y_d = Variable(N_t)
#Q_fr = ComplexVariable(N_f)
S_fr_re = Variable(N_f)
S_fr_im = Variable(N_f)
t1 = Variable()
t2 = Variable()
t3 = Variable(N_f)
t4 = Variable()
#c2 = (Q_fr == G_q_fr*b)

c3_re = (S_fr_re == real(SA_fr)*b + real(Sb_fr))        # Aff
c3_im = (S_fr_im == imag(SA_fr)*b + imag(Sb_fr))        # Aff

c4 = (t2 == t1 - 5^2)                                   # Aff
c5 = sumsquares(b) <= t1;                               # SOC    ∋ (b,t1)
c7 = (t2 <= 0)                                          # NonPos ∋ t2

c8 = (abs2(S_fr_re) + abs2(S_fr_im) <= t3)              # Many SOCs ∋ (S_fr,t3)
c9 = (t3 == 1.6^2)                                      # Affine

#cost = sumsquares(y_d)
c1 = (y_d == PSA*b+PSb)                                 # Affine
c10 = (sumsquares(y_d) <= t4)                           # SOC ∋ (y_d, t4)
c11 = (t4 == 1.1207+0.01)                               # Affine

prob = minimize(0, c1, c3_re, c3_im, c4, c5, c6, c7, c8, c9, c10, c11)

################# As one affine + cone in convex

Astag1 = toeplitz2([1 0 0; zeros(N_f-1, 3)])
Astag2 = toeplitz2([0 1 0; zeros(N_f-1, 3)])
Astag3 = toeplitz2([0 0 1; zeros(N_f-1, 3)])

A1 = [ zeros(N_f,1) -real(SA_fr)    zeros(N_f,1)    Astag2              zeros(N_f,1)    zeros(N_f, N_t) ]
A2 = [ zeros(N_f,1) -imag(SA_fr)    zeros(N_f,1)    Astag3              zeros(N_f,1)    zeros(N_f, N_t) ]
A3 = [   -1         zeros(1,N)      1               zeros(1,N_f*3)      zeros(1, 1)     zeros(1,   N_t) ]
A4 = [ zeros(N_f,1) zeros(N_f,N)    zeros(N_f,1)    Astag1              zeros(N_f,1)    zeros(N_f, N_t) ]
A5 = [ zeros(N_t,1) -PSA            zeros(N_t,1)    zeros(N_t,N_f*3)    zeros(N_t,1)    I               ]
A6 = [   0          zeros(1,N)      0               zeros(1,N_f*3)      1               zeros(1,   N_t) ]

x = [t1;b;t2;
    vcat(([t3[i];S_fr_re[i];S_fr_im[i];] for i in 1:N_f)...);
    t4; y_d]

# For comparison
xc = [sqrt(t1);b;t2;
    vcat(([sqrt(t3[i]);S_fr_re[i];S_fr_im[i];] for i in 1:N_f)...);
    sqrt(t4); y_d]

Aall = [A1;A2;A3;A4;A5;A6]
ball = [real(Sb_fr); imag(Sb_fr); -5^2; fill(1.6^2, N_f); PSb; 1.1207+0.01]

Aff = (Aall*x == ball)

Cone = [sumsquares(b) <= t1; t2 <= 0;
        vcat((abs2(S_fr_re[i]) + abs2(S_fr_im[i]) <= t3[i] for i in 1:N_f)...);
        sumsquares(y_d) <= t4]

prob = minimize(0, Aff, Cone...)
using SCS
@time solve!(prob, SCSSolver())
using Plots
plot(evaluate(b))

############ Using ProximalOperators
using ProximalOperators

ball2 = [real(Sb_fr); imag(Sb_fr); -5; fill(1.6, N_f); PSb; sqrt(1.1207+0.01)]

S1 = IndAffine(Aall, ball2)
fs = (IndSOC(), IndNonpositive(),
        (IndSOC() for i in 1:N_f)...,
        IndSOC())
idxs = (1:(N+1), (N+2):(N+2),
        ((N+3+i-1):(N+3+i+1) for i in 1:3:(N_f*3))...,
        (N+2+3*N_f+1):(N+2+3*N_f+N_t+1))
S2 = SlicedSeparableSum(fs, tuple.(idxs))

y0 = randn(size(Aall,2))
x0 = randn(size(Aall,2))

using FirstOrderSolvers
alg = DR()
problem = FirstOrderSolvers.Feasibility(S1, S2, size(Aall,2))
sol = FirstOrderSolvers.solve!(p, solver, checki=10, eps=1e-9)

using ProximalAlgorithms
solver = ProximalAlgorithms.DouglasRachford(gamma=1.0,maxit=10000,verbose=true)

x0, y0, it = solver(x0, f=S1, g=S2)
@profiler solver(x0, f=S1, g=S2)

tmp1 = similar(x0)
tmp2 = similar(x0)
tmp3 = similar(x0)
@time for i = 1:1000
    prox!(x0, S1, y0)
    tmp1 .= 2. *x0 .- y0
    prox!(tmp2, S2, tmp1)
    tmp3 .= 2. *tmp2 .- tmp1

    i%100 == 0 && (print("Err: "); println(norm(x0-tmp2)))
    i%100 == 0 && (print("Ang: "); println(acos(min(dot(y0-x0, tmp2-tmp1)/norm(x0-y0)/norm(tmp2-tmp1),1))))
    # Shadow in tmp2
    y0 .= (tmp3 .+ y0)./2
end

bsol = tmp2[2:N+1]
plot(bsol)

using Convex

b = Variable(N)

#T_Q = toeplitz2([b; zeros(N_t-N)])
#y_d = T_P*(Matrix{Float64}(I,N_t,N_t) - T_P*T_Q)*d1
#y_d = T_P*d1 - T_P*T_P*toeplitz2(d1,N)*b

y_d = PSA*b+PSb
Q_fr = G_q_fr*b
cost = sumsquares(y_d)
const1 = abs(SA_fr*b+Sb_fr) <= 1.6
const2 = sumsquares(b) <= 5^2
prob = minimize(cost, const1, const2)

using SCS
solve!(prob, SCSSolver())

using Plots
plot(evaluate(b))

plot(Ω./2 ./pi, abs.(1 .- P_fr.*evaluate(Q_fr)), yscale=:log10, xscale=:log10)

C_fr = vec(evaluate(Q_fr) ./ (1 .- P_fr .*evaluate(Q_fr)))

plot(Ω/2/pi, abs.(C_fr), yscale=:log10, xscale=:log10)

ext: jl


text:
using DataFramesMeta, Random, BenchmarkTools, Statistics

N = 1_000_000
K = 10

df = DataFrame(
	id = rand([Symbol("id", i) for i=1:K], N),
	v1 = rand(1:5, N),
	v2 = rand(1:5, N),
	v3 = rand(N)
);

gd = groupby(df, :id);

function complicated_vec(x, y, z)
	mx = mean(x)
	sy = std(y)
	vz = var(z)

	temp = (mx * sy) / vz

	a = @. mx * sy + vz + temp

	ma = mean(a)
	sa = std(a)

	@. (a - ma) / sa
end

function complicated_scalar(x, y, z)
	mx = mean(x)
	sy = std(y)
	vz = var(z)

	temp = (mx * sy) / vz

	a = mx * sy + vz + temp

	a * 1000
end

demean(x) = x .- mean(x)

function DataFrames_timings(df, gd)
	df_res = transform(df,
		:v1 => (t -> t .- mean(t)) => :res1,
		:v2 => demean => :res2,
		[:v1, :v2] => (+) => :res3,
		:id => string => :res4,
		[:v1, :v2, :v3] => complicated_vec => :res5,
		[:v1, :v2, :v3] => ((a, b, c) -> @. a + b * c * c + a) => :res6a,
		[:v1, :v2, :v3] =>
		(
			(a, b, c) -> begin
				d = Vector{Float64}(undef, length(a))
				for i in eachindex(d)
					d[i] = a[i] + b[i] * c[i] * c[i] + a[i]
				end
				d
			end
		) => :res6b
	)

	gd_res = combine(gd,
		:v1 => (t -> mean(t)) => :res7,
		:v2 => mean => :res8,
		[:v1, :v2] => ((t, s) -> std(t) + std(s)) => :res9,
		[:v1, :v2, :v3] => complicated_scalar => :res10,
		[:v1, :v2, :v3] => ((a, b, c)  -> first(a) + mean(b) * std(a) + last(c)) => :res11,
	)

	return(df_res, gd_res)
end

function DataFramesMeta_timings(df, gd)
	df_res = @transform(df,
		res1 = :v1 .- mean(:v1),
		res2 = demean(:v2),
		res3 = :v1 + :v2,
		res4 = string(:id),
		res5 = complicated_vec(:v1, :v2, :v3),
		res6a = @.(:v1 + :v2 + :v3 * :v3 + :v1),
		res6b = begin
			d = Vector{Float64}(undef, length(:v1))
			for i in eachindex(d)
				d[i] = :v1[i] + :v2[i] * :v3[i] * :v3[i] + :v1[i]
			end
			d
		end
	)

	gd_res = @combine(gd,
		res7 = mean(:v1),
		res8 = (t -> mean(t))(:v2),
		res9 = std(:v1) + std(:v2),
		res10 = complicated_scalar(:v1, :v2, :v3),
		res11 = first(:v1) + mean(:v2) * std(:v1) + last(:v3)
	)

	return(df_res, gd_res)
end

println("DataFrames benchmark timings")
@btime DataFrames_timings($df, $gd);
println("DataFramesMeta benchmark, timings")
@btime DataFramesMeta_timings($df, $gd);

N = 10
K = 10

df2 = DataFrame(
	id = rand([Symbol("id", i) for i=1:K], N),
	v1 = rand(1:5, N),
	v2 = rand(1:5, N),
	v3 = rand(N)
);

println("DataFramesMeta raw timing")
@time @select(df2, res1 = :v1 .- mean(:v1));
@time @select(df2, res2 = demean(:v2));
@time @select(df2, res3 = :v1 + :v2);
@time @select(df2, res4 = string(:id));
@time @select(df2, res5 = complicated_vec(:v1, :v2, :v3));
@time @select(df2, res6a = @.(:v1 + :v2 + :v3 * :v3 + :v1));
@time @select(df2, res6b = begin
	d = Vector{Float64}(undef, length(:v1))
	for i in eachindex(d)
		d[i] = :v1[i] + :v2[i] * :v3[i] * :v3[i] + :v1[i]
	end
	d
end);

nothing


text:
using Utils
using ProblemParser

using Chain
using Dictionaries

import IterTools

module Data
using Dictionaries
struct Point3D
    x::Int32
    y::Int32
    z::Int32
end
Point3D() = Point3D(0, 0, 0)
struct Sparse3D
    trues::Indices{Point3D}
end
Sparse3D(points::Vector{Point3D}) = points |> distinct |> Sparse3D
struct AllOrientations
    cube::Sparse3D
end
end
Point3D=Data.Point3D
Sparse3D=Data.Sparse3D
AllOrientations=Data.AllOrientations

Base.:(==)(lhs::Sparse3D, rhs::Sparse3D) = lhs.trues == rhs.trues

rot90(cube::Sparse3D, ::Union{Val{1},Val{:x}})::Sparse3D = map(cube.trues) do p
    Point3D(p.x, -p.z, p.y)
end |> Indices |> Sparse3D

rot90(cube::Sparse3D, ::Union{Val{2},Val{:y}})::Sparse3D = map(cube.trues) do p
    Point3D(p.z, p.y, -p.x)
end |> Indices |> Sparse3D

rot90(cube::Sparse3D, ::Union{Val{3},Val{:z}})::Sparse3D = map(cube.trues) do p
    Point3D(-p.y, p.x, p.z)
end |> Indices |> Sparse3D

mirror(cube::Sparse3D)::Sparse3D = map(cube.trues) do p
    Point3D(-p.x, -p.y, -p.z)
end |> Indices |> Sparse3D

Base.iterate(start::AllOrientations) = start.cube, (start.cube,1)

function Base.iterate(::AllOrientations, (cube, state))
    state > lastindex(rots) && return nothing
    transformed = rots[state](cube)
    transformed, (transformed,state + 1)
end

Base.length(::AllOrientations) = 24
Base.eltype(::AllOrientations) = Sparse3D

# Build all 24 orientations by this series of 24 transformations.
rots = Function[curry2nd(rot90, Val(d)) for d in [:x, :y, :z] for _ in 1:3]
push!(rots, curry2nd(rot90, Val(:x)), curry2nd(rot90, Val(:y)))
push!(rots, mirror)
append!(rots, rots[begin:end-1])
# Make sure it delivers 24 distinct points.
@assert length(Set(c for c in AllOrientations(Sparse3D([Point3D(1,2,3)])))) == 24

using Test
@testset begin
    sample = Sparse3D([Point3D(3,0,0), Point3D(0,5,0), Point3D(0,0,2)])

    @testset "rotate z" begin
        rotations = @chain sample begin
            IterTools.iterated(curry2nd(rot90, Val(:z)), _)
            Iterators.take(5)
            map(s -> s.trues, _)
            collect
        end
        print(rotations)

        @test rotations == [
            sample.trues
            Indices([Point3D(0,3,0), Point3D(-5,0,0), Point3D(0,0,2)])
            Indices([Point3D(-3,0,0), Point3D(0,-5,0), Point3D(0,0,2)])
            Indices([Point3D(0,-3,0), Point3D(5,0,0), Point3D(0,0,2)])
            sample.trues
        ]
    end
end

file = find_input(@__FILE__)#,"example.txt")
load(file) = @chain file begin
    slurp
    parse(Blocks(FirstRest(Lines(), Extract(), Lines(Split(',', Convert(Int32))))), _)
    map(_) do (i, points)
        i, map(tup3 -> Point3D(tup3...), points) |> Sparse3D
    end
end
input = load(file)

manhatten(lhs::Point3D, rhs::Point3D) = abs(lhs.x - rhs.x) + abs(lhs.y - rhs.y) + abs(lhs.z - rhs.z)

fingerprint(cubes::Sparse3D) = Set(manhatten(a, b) for a in cubes.trues for b in cubes.trues if a != b)

function part1(input)
    threshold = 12*11÷2
    base, rest... = input
    basefingerprint = fingerprint(base[2])
    while !isempty(rest)
        found = false
        for idx in eachindex(rest)
            (i, other) = rest[idx]
            println("check $i")
            fp = fingerprint(other)
            intersect!(fp, basefingerprint)
            if length(fp) >= threshold
                println("match for $i")
                deleteat!(rest, idx)
                found = true
                break
            end
        end
        !found && error("no match found in this round")
    end
end
part1(input)

assertequal(part1(input))

part2() = count(Iterators.product(1:input.ur.x, -99:99)) do tup
    any(occursin(input), Point(tup...))
end

assertequal(part2(), 1928)

"""
--- scanner 0 ---
404;-588;-901
528;-643;409
-838;591;734
390;-675;-793
-537;-823;-458
-485;-357;347
-345;-311;381
-661;-816;-575
-876;649;763
-618;-824;-621
553;345;-567
474;580;667
-447;-329;318
-584;868;-557
544;-627;-890
564;392;-477
455;729;728
-892;524;684
-689;845;-530
423;-701;434
7;-33;-71
630;319;-379
443;580;662
-789;900;-551
459;-707;401

--- scanner 1 ---
686;422;578
605;423;415
515;917;-361
-336;658;858
95;138;22
-476;619;847
-340;-569;-846
567;-361;727
-460;603;-452
669;-402;600
729;430;532
-500;-761;534
-322;571;750
-466;-666;-811
-429;-592;574
-355;545;-477
703;-491;-529
-328;-685;520
413;935;-424
-391;539;-444
586;-435;557
-364;-763;-893
807;-499;-711
755;-354;-619
553;889;-390

-618;-824;-621
-537;-823;-458
-447;-329;318
404;-588;-901
544;-627;-890
528;-643;409
-661;-816;-575
390;-675;-793
423;-701;434
-345;-311;381
459;-707;401
-485;-357;347
"""

max_stars_repo_name: CmdQ/AoC

max_issues_repo_name: CmdQ/AoC

max_forks_repo_name: CmdQ/AoC


text:
using Dates
using Random
using Test

using MarketData

using TimeSeries

@testset "readwrite" begin

@testset "readwrite parses IO correctly" begin
    io  = """DateTime,Open
             2010/1/4|9:01:00,7300
             2010/1/4|9:02:00,7316
             2010/1/4|9:03:00,7316
             2010/1/4|9:04:00,7316
             2010/1/4|9:05:00,7316"""
    tm0 = readtimearray(IOBuffer(io), format="yyyy/mm/dd|HH:MM:SS")

    @testset "Input as IOBuffer and specifying DateTime string format for reading" begin
        @test length(tm0)       == 5
        @test size(values(tm0)) == (5,1)
        @test timestamp(tm0)[4] == DateTime(2010,1,4,9,4)
    end
end

@testset "readwrite parses csv correctly" begin
    tm1 = readtimearray(
        joinpath(@__DIR__, "data/datetime3.csv"),
        format="yyyy/mm/dd|HH:MM:SS")
    tm2 = readtimearray(
        joinpath(@__DIR__, "data/datetime3.csv"),
        format="yyyy/mm/dd|HH:MM:SS", meta="foo")
    tm3 = readtimearray(
        joinpath(@__DIR__, "data/read_example_delim.csv"),
        format="dd/mm/yyyy HH:MM", delim=';')
    tm4 = readtimearray(
        joinpath(@__DIR__, "data/headless.csv"),
        format="yyyy-mm-dd HH:MM:SS", header=false)

    @testset "Specifying DateTime string format for reading" begin
        @test length(tm1)        == 5
        @test size(values(tm1))  == (5,1)
        @test timestamp(tm1)[4]  == DateTime(2010,1,4,9,4)
        @test_throws(
            ArgumentError,
            readtimearray(joinpath(dirname(@__FILE__), "data/datetime3.csv")))
    end

    @testset "readtimearray accepts meta field" begin
        @test meta(tm2) == "foo"
    end

    @testset "readtimearray works with arbitrary delimiters" begin
        @test length(tm3)       == 2
        @test size(values(tm3)) == (2,2)
        @test timestamp(tm3)[2] == DateTime(2015,1,1,1,0)
        @test values(tm3)[1,1]  == 10.42
    end

    @testset "headless csv" begin
        @test length(tm4)       == 2
        @test size(values(tm4)) == (2, 4)
        @test values(tm4)       == [1 2 3 4; 5 6 7 8]
    end
end

@testset "readwrite parses MarketData objects correctly" begin
    @testset "1d values array works" begin
        @test typeof(values(cl)) == Array{Float64,1}
    end

    @testset "2d values array works" begin
        @test typeof(values(ohlc)) == Array{Float64,2}
    end

    @testset "timestamp parses to correct type" begin
        @test typeof(timestamp(cl))        == Vector{Date}
        @test typeof(timestamp(datetime1)) == Vector{DateTime}
    end
end

@testset "writetimearray method with default parameters works" begin
    mktemp() do filename, _io
        uohlc = uniformspace(ohlc)
        writetimearray(uohlc, filename)
        readback = readtimearray(filename)

        @test colnames(uohlc)  == colnames(readback)
        @test timestamp(uohlc) == timestamp(readback)
        @test isequal(values(uohlc), values(readback))
    end
end

@testset "writetimearray method with a delimiter works" begin
    mktemp() do filename, _io
        uohlc = uniformspace(ohlc)
        writetimearray(uohlc[1:5], filename, delim=';')
        readback = readtimearray(filename, delim=';')

        @test colnames(uohlc[1:5])  == colnames(readback)
        @test timestamp(uohlc[1:5]) == timestamp(readback)
        @test isequal(values(uohlc[1:5]), values(readback))
    end
end

@testset "writetimearray method with no header works" begin
    mktemp() do filename, _io
        uohlc = uniformspace(ohlc)
        writetimearray(uohlc[1:5], filename, header=false)
        readback = readtimearray(filename, header=false)

        @test timestamp(uohlc[1:5]) == timestamp(readback)
        @test isequal(values(uohlc[1:5]), values(readback))
    end
end

@testset "writetimearray method with a timestamp format works" begin
    mktemp() do filename, _io
        uohlc = uniformspace(ohlc)
        writetimearray(uohlc[1:5], filename, format="yyyy/mm/dd")
        readback = readtimearray(filename, format="yyyy/mm/dd")

        @test timestamp(uohlc[1:5]) == timestamp(readback)
    end
end

end  # @testset "readwrite"


text:
using Test
using Random
using DataValues
# import DataValues: unsafe_getindex_notnull, unsafe_getvalue_notnull

@testset "DataValueArray: Indexing" begin

x = DataValueArray{Int}((5, 2))

for i in eachindex(x)
    x[i] = i
end

for i in eachindex(x)
    y = x[i]
    @test isa(y, DataValue{Int})
    @test get(y) === i
end

for i in eachindex(x)
    x[i] = DataValue{Int}()
end

for i in eachindex(x)
    y = x[i]
    @test isa(y, DataValue{Int})
    @test isna(y)
end

_values = rand(10, 10)
_isna = rand(Bool, 10, 10)
X = DataValueArray(_values, _isna)

# Base.getindex{T, N}(X::DataValueArray{T, N}, I::DataValue{Int}...)
@test_throws DataValueException getindex(X, DataValue{Int}(), DataValue{Int}())
if _isna[1]
    @test isna(getindex(X, DataValue(1)))
else
    @test isequal(getindex(X, DataValue(1)), DataValue(_values[1]))
end

# Scalar getindex
for i = 1:100
    if _isna[i]
        @test isna(X[i])
    else
        @test isequal(X[i], DataValue(_values[i]))
    end
end
for i = 1:10, j = 1:10
    if _isna[i, j]
        @test isna(X[i, j])
    else
        @test isequal(X[i, j], DataValue(_values[i, j]))
    end
end

# getindex with AbstractVectors
rg = 2:9
v = X[rg]
for i = 1:length(rg)
    if _isna[rg[i]]
        @test isna(v[i])
    else
        @test isequal(v[i], DataValue(_values[rg[i]]))
    end
end

v = X[rg, 9]
for i = 1:length(rg)
    if _isna[rg[i], 9]
        @test isna(v[i])
    else
        @test isequal(v[i], DataValue(_values[rg[i], 9]))
    end
end

rg2 = 5:7
v = X[rg, rg2]
for j = 1:length(rg2), i = 1:length(rg)
    if _isna[rg[i], rg2[j]]
        @test isna(v[i, j])
    else
        @test isequal(v[i, j], DataValue(_values[rg[i], rg2[j]]))
    end
end

# range indexing
Z_values = reshape(collect(1:125), (5,5,5))
Z = DataValueArray(Z_values)

@test isequal(Z[1, 1:4, 1], DataValueArray([1, 6, 11, 16]))

# getindex with AbstractVector{Bool}
b = bitrand(10, 10)
rg = (LinearIndices(b))[findall(b)]
v = X[b]
for i = 1:length(rg)
    if _isna[rg[i]]
        @test isna(v[i])
    else
        @test isequal(v[i], DataValue(_values[rg[i]]))
    end
end

# getindex with DataValueVector with null entries throws error
# DA TODO Decide whether I want to keep this
# @test_throws DataValueException X[DataValueArray([1, 2, 3, NA])]

# getindex with DataValueVector and non-null entries
# DA TODO Decide whether I want this or not
# @test isequal(X[DataValueArray([1, 2, 3])], X[[1, 2, 3]])

# indexing with DataValues

n = rand(1:5)
siz = [ rand(2:5) for i in n ]
A = rand(siz...)
M = rand(Bool, siz...)
Z = DataValueArray(A, M)
i = rand(1:length(Z))
@test isequal(Z[DataValue(i)], Z[i])
I = [ rand(1:size(Z,i)) for i in 1:n ]
NI = [ DataValue(i) for i in I ]
@test isequal(Z[NI...], Z[I...])

#----- test setindex! -----#

# setindex! with scalar indices
_values = rand(10, 10)
for i = 1:100
    X[i] = _values[i]
end
@test isequal(X, DataValueArray(_values))

_values = rand(10, 10)
for i = 1:10, j = 1:10
    X[i, j] = _values[i, j]
end
@test isequal(X, DataValueArray(_values))

_values = rand(10, 10)
for i = 1:10, j = 1:10
    X[i, j] = DataValue(_values[i, j])
end
@test isequal(X, DataValueArray(_values))

# ----- test nullify! -----#
_isna = bitrand(10, 10)
for i = 1:100
    if _isna[i]
        X[i] = NA
    end
end

# setindex! with scalar and vector indices
rg = 2:9
_values[rg] .= 1.0
X[rg] .= 1.0
for i = 1:length(rg)
    @test isequal(X[rg[i]], DataValue(1.0))
end

# setindex! with NA and vector indices
rg = 5:13
_isna[rg] .= true
# TODO This should be changed to .=, but currently crashes on 0.7
X[rg] .= Ref(NA)
for i = 1:length(rg)
    @test isna(X[rg[i]])
end

# setindex! with vector and vector indices
rg = 12:67
_values[rg] .= rand(length(rg))
X[rg] .= _values[rg]
for i = 1:length(rg)
    @test isequal(X[rg[i]], DataValue(_values[rg[i]]))
end

#----- test UNSAFE INDEXING -----#

# DA TODO Disabled for now
# X = DataValueArray([1, 2, 3, 4, 5], [true, false, false, false, false])

# @test isequal(unsafe_getindex_notnull(X, 1), DataValue(1))
# @test isequal(unsafe_getindex_notnull(X, 2), DataValue(2))
# @test isequal(unsafe_getvalue_notnull(X, 1), 1)
# @test isequal(unsafe_getvalue_notnull(X, 2), 2)

#----- test Base.checkbounds -----#

X = DataValueArray([1:10...])
b = vcat(false, fill(true, 9))

# Base.checkindex(::Type{Bool}, inds::UnitRange, i::DataValue)
# DA TODO Decide whether I want these
# @test_throws DataValueException checkindex(Bool, 1:1, DataValue{Int}())
# @test checkindex(Bool, 1:10, DataValue(1)) == true
# @test isequal(X[DataValue(1)], DataValue(1))

# Base.checkindex{N}(::Type{Bool}, inds::UnitRange, I::DataValueArray{Bool, N})
# DA TODO Decide whether I want these
# @test checkindex(Bool, 1:5, DataValueArray([true, false, true, false, true]))
@test isequal(X[b], DataValueArray([2:10...]))

# Base.checkindex{T<:Real}(::Type{Bool}, inds::UnitRange, I::DataValueArray{T})
# DA TODO Decide whether I want these
# @test checkindex(Bool, 1:10, DataValueArray([1:10...]))
# @test checkindex(Bool, 1:10, DataValueArray([10, 11])) == false
# @test_throws BoundsError checkbounds(X, DataValueArray([10, 11]))

#---- test Base.to_index -----#

# Base.to_index(X::DataValueArray)
# DA TODO Decide whether I want these
# @test Base.to_index(X) == [1:10...]
# push!(X, DataValue{Int}())
# @test_throws DataValueException Base.to_index(X)

end


text:
const ref_panel_coordinates = [ -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0 ]

function set_panel_vertices_coordinates!( pconn :: Ptr{P4est_wrapper.p4est_connectivity_t},
                                             coarse_discrete_model :: DiscreteModel{2,2},
                                             panel )
  @assert num_cells(coarse_discrete_model) == 6
  @assert panel ≤ num_cells(coarse_discrete_model)
  @assert panel ≥ 1
  trian=Triangulation(coarse_discrete_model)
  cell_vertices=Gridap.Geometry.get_cell_node_ids(trian)
  #println(cell_vertices)
  cell_vertices_panel=cell_vertices[panel]
  conn=pconn[]
  vertices=unsafe_wrap(Array,
                       conn.vertices,
                       length(Gridap.Geometry.get_node_coordinates(coarse_discrete_model))*3)
  for (l,g) in enumerate(cell_vertices_panel)
     # println("XXX $(l) $(g) $(lref)")
     vertices[(g-1)*3+1]=ref_panel_coordinates[(l-1)*2+1]
     vertices[(g-1)*3+2]=ref_panel_coordinates[(l-1)*2+2]
  end
end

function setup_cubed_sphere_coarse_discrete_model()
    # 6 panels (cells), 4 corners (vertices) each panel
    ptr  = [ 1, 5, 9, 13, 17, 21, 25 ]
    data = [ 1,2,3,4, 2,5,4,6, 7,8,5,6, 1,3,7,8, 3,4,8,6, 1,7,2,5  ]
    cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
    node_coordinates = Vector{Point{2,Float64}}(undef,8)
    for i in 1:length(node_coordinates)
      node_coordinates[i]=Point{2,Float64}(0.0,0.0)
    end

    polytope=QUAD
    scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
    cell_types=collect(Fill(1,length(cell_vertex_lids)))
    cell_reffes=[scalar_reffe]
    grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                            cell_vertex_lids,
                                            cell_reffes,
                                            cell_types,
                                            Gridap.Geometry.NonOriented())
    Gridap.Geometry.UnstructuredDiscreteModel(grid)
end

function generate_ptr(ncells)
  nvertices=4
  ptr  = Vector{Int}(undef,ncells+1)
  ptr[1]=1
  for i=1:ncells
    ptr[i+1]=ptr[i]+nvertices
  end
  ptr
end

function generate_cell_coordinates_and_panels(parts,
                                   coarse_discrete_model,
                                   ptr_pXest_connectivity,
                                   ptr_pXest,
                                   ptr_pXest_ghost)

  Dc=2
  PXEST_CORNERS=4
  pXest_ghost = ptr_pXest_ghost[]
  pXest = ptr_pXest[]

  # Obtain ghost quadrants
  ptr_ghost_quadrants = Ptr{P4est_wrapper.p4est_quadrant_t}(pXest_ghost.ghosts.array)

  tree_offsets = unsafe_wrap(Array, pXest_ghost.tree_offsets, pXest_ghost.num_trees+1)
  dcell_coordinates_and_panels=map_parts(parts) do part
     ncells=pXest.local_num_quadrants+pXest_ghost.ghosts.elem_count
     panels = Vector{Int}(undef,ncells)
     data = Vector{Point{Dc,Float64}}(undef,ncells*PXEST_CORNERS)
     ptr  = generate_ptr(ncells)
     current=1
     current_cell=1
     vxy=Vector{Cdouble}(undef,Dc)
     pvxy=pointer(vxy,1)
     for itree=1:pXest_ghost.num_trees
       tree = p4est_tree_array_index(pXest.trees, itree-1)[]
       if tree.quadrants.elem_count > 0
          set_panel_vertices_coordinates!( ptr_pXest_connectivity, coarse_discrete_model, itree)
       end
       for cell=1:tree.quadrants.elem_count
          panels[current_cell]=itree
          quadrant=p4est_quadrant_array_index(tree.quadrants, cell-1)[]
          for vertex=1:PXEST_CORNERS
            GridapP4est.p4est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                  p4est_topidx_t(itree-1),
                                                  quadrant.x,
                                                  quadrant.y,
                                                  quadrant.level,
                                                  Cint(vertex-1),
                                                  pvxy)
            data[current]=Point{Dc,Float64}(vxy...)
            current=current+1
          end
          current_cell=current_cell+1
       end
     end

     # Go over ghost cells
     for i=1:pXest_ghost.num_trees
      if tree_offsets[i+1]-tree_offsets[i] > 0
        set_panel_vertices_coordinates!( ptr_pXest_connectivity, coarse_discrete_model, i)
      end
      for j=tree_offsets[i]:tree_offsets[i+1]-1
          panels[current_cell]=i
          quadrant = ptr_ghost_quadrants[j+1]
          for vertex=1:PXEST_CORNERS
            GridapP4est.p4est_get_quadrant_vertex_coordinates(ptr_pXest_connectivity,
                                                     p4est_topidx_t(i-1),
                                                     quadrant.x,
                                                     quadrant.y,
                                                     quadrant.level,
                                                     Cint(vertex-1),
                                                     pvxy)

          #  if (MPI.Comm_rank(comm.comm)==0)
          #     println(vxy)
          #  end
          data[current]=Point{Dc,Float64}(vxy...)
          current=current+1
         end
         current_cell=current_cell+1
       end
     end
     Gridap.Arrays.Table(data,ptr), panels
  end
end

function generate_cube_grid_geo(cell_coordinates_and_panels)
  function map_panel_xy_2_xyz(xy,panel)
    a,b=xy
    if panel==1
      x=Point(1.0,a,b)
    elseif panel==2
      x=Point(-a,1.0,b)
    elseif panel==3
      x=Point(-1.0,b,a)
    elseif panel==4
      x=Point(-b,-1.0,a)
    elseif panel==5
      x=Point(-b,a,1.0)
    elseif panel==6
      x=Point(-a,b,-1.0)
    end
    x
  end

  map_parts(cell_coordinates_and_panels) do (cell_coordinates,panels)
     ptr  = generate_ptr(length(cell_coordinates))
     data = collect(1:length(cell_coordinates)*4)
     cell_vertex_lids=Gridap.Arrays.Table(data,ptr)
     node_coordinates=Vector{Point{3,Float64}}(undef,length(cell_coordinates)*4)

     cache=array_cache(cell_coordinates)
     current=1
     for cell=1:length(cell_coordinates)
        current_cell_coordinates=getindex!(cache,cell_coordinates,cell)
        for i=1:length(current_cell_coordinates)
          node_coordinates[current]=map_panel_xy_2_xyz(current_cell_coordinates[i],panels[cell])
          current=current+1
        end
     end

     polytope=QUAD
     scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
     cell_types=collect(Fill(1,length(cell_vertex_lids)))
     cell_reffes=[scalar_reffe]
     grid=Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                      cell_vertex_lids,
                                      cell_reffes,
                                      cell_types,
                                      Gridap.Geometry.NonOriented())
     grid
  end
end

function generate_cube_grid_top(cell_vertex_lids_nlvertices)
  map_parts(cell_vertex_lids_nlvertices) do (cell_vertex_lids,nlvector)
     node_coordinates=Vector{Point{2,Float64}}(undef,nlvector)
     polytope=QUAD
     scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
     cell_types=collect(Fill(1,length(cell_vertex_lids)))
     cell_reffes=[scalar_reffe]
    #  if (part==2)
    #   println(cell_vertex_lids)
    #  end
     cell_vertex_lids_gridap=Gridap.Arrays.Table(cell_vertex_lids.data,cell_vertex_lids.ptrs)
     grid=Gridap.Geometry.UnstructuredGrid(node_coordinates,
                                      cell_vertex_lids_gridap,
                                      cell_reffes,
                                      cell_types,
                                      Gridap.Geometry.NonOriented())
     grid
  end
end

function CubedSphereDiscreteModel(
  parts::MPIData,
  num_uniform_refinements::Int;
  radius=1,
  p4est_verbosity_level=P4est_wrapper.SC_LP_DEFAULT)

  comm = parts.comm

  sc_init(parts.comm, Cint(true), Cint(true), C_NULL, p4est_verbosity_level)
  p4est_init(C_NULL, p4est_verbosity_level)

  Dc=2

  coarse_discrete_model=setup_cubed_sphere_coarse_discrete_model()

  ptr_pXest_connectivity=GridapP4est.setup_pXest_connectivity(coarse_discrete_model)

  # Create a new forest
  ptr_pXest = GridapP4est.setup_pXest(Val{Dc},comm,ptr_pXest_connectivity,num_uniform_refinements)

  # Build the ghost layer
  ptr_pXest_ghost=GridapP4est.setup_pXest_ghost(Val{Dc},ptr_pXest)

  cellindices = GridapP4est.setup_cell_prange(Val{Dc},parts,ptr_pXest,ptr_pXest_ghost)

  ptr_pXest_lnodes=GridapP4est.setup_pXest_lnodes(Val{Dc}, ptr_pXest, ptr_pXest_ghost)

  cell_vertex_gids=GridapP4est.generate_cell_vertex_gids(ptr_pXest_lnodes,cellindices)

  cell_vertex_lids_nlvertices=GridapP4est.generate_cell_vertex_lids_nlvertices(cell_vertex_gids)

  cell_coordinates_and_panels=generate_cell_coordinates_and_panels(parts,
                                             coarse_discrete_model,
                                             ptr_pXest_connectivity,
                                             ptr_pXest,
                                             ptr_pXest_ghost)

  cube_grid_geo=generate_cube_grid_geo(cell_coordinates_and_panels)
  cube_grid_top=generate_cube_grid_top(cell_vertex_lids_nlvertices)

  # # do_on_parts(comm, cell_coordinates) do part, cell_coords
  # #   if part==1
  # #     println(cell_coords[9:12])
  # #   end
  # # end

  ddiscretemodel=
  map_parts(cube_grid_geo,cube_grid_top) do cube_grid_geo, cube_grid_top
    D2toD3AnalyticalMapCubedSphereDiscreteModel(cube_grid_geo, cube_grid_top, radius=radius)
  end

  # Write forest to VTK file
  #p4est_vtk_write_file(unitsquare_forest, C_NULL, "my_step")

  # Destroy lnodes
  GridapP4est.p4est_lnodes_destroy(ptr_pXest_lnodes)
  GridapP4est.p4est_ghost_destroy(ptr_pXest_ghost)
  # Destroy the forest
  GridapP4est.p4est_destroy(ptr_pXest)
  # Destroy the connectivity
  GridapP4est.p4est_connectivity_destroy(ptr_pXest_connectivity)

  sc_finalize()

  GridapDistributed.DistributedDiscreteModel(ddiscretemodel,cellindices)
end

struct D2toD3AnalyticalMapCubedSphereTriangulation{M} <: Triangulation{2,3}
  model::M
end

function Gridap.CellData.get_normal_vector(trian::D2toD3AnalyticalMapCubedSphereTriangulation)
  cell_normal = Gridap.Geometry.get_facet_normal(trian)
  Gridap.CellData.GenericCellField(cell_normal,trian,ReferenceDomain())
end

function _unit_outward_normal(v::Gridap.Fields.MultiValue{Tuple{2,3}})
  n1 = v[1,2]*v[2,3] - v[1,3]*v[2,2]
  n2 = v[1,3]*v[2,1] - v[1,1]*v[2,3]
  n3 = v[1,1]*v[2,2] - v[1,2]*v[2,1]
  n = VectorValue(n1,n2,n3)
  n/norm(n)
end

function Gridap.Geometry.get_facet_normal(trian::D2toD3AnalyticalMapCubedSphereTriangulation)
  # Get the Jacobian of the cubed sphere mesh
  map   = get_cell_map(trian)
  Jt    = lazy_map(∇,map)
  p=lazy_map(Operation(_unit_outward_normal),Jt)
  p
end

# Triangulation API

# Delegating to the underlying face Triangulation

Gridap.Geometry.get_cell_coordinates(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_coordinates(trian.model.cube_grid_geo)

Gridap.Geometry.get_reffes(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_reffes(trian.model.cube_grid_geo)

Gridap.Geometry.get_cell_type(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_type(trian.model.cube_grid_geo)

Gridap.Geometry.get_node_coordinates(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_node_coordinates(trian.model.cube_grid_geo)

Gridap.Geometry.get_cell_node_ids(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = Gridap.Geometry.get_cell_node_ids(trian.model.cube_grid_geo)

Gridap.Geometry.get_cell_map(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = trian.model.cell_map

Gridap.Geometry.get_background_model(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = trian.model

function Gridap.Geometry.get_glue(a::D2toD3AnalyticalMapCubedSphereTriangulation,D::Val{2})
  nc=num_cells(a.model.cube_model_top)
  tface_to_mface=Gridap.Fields.IdentityVector(nc)
  tface_to_mface_map=Fill(Gridap.Fields.GenericField(identity),nc)
  mface_to_tface=tface_to_mface
  Gridap.Geometry.FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

Gridap.Geometry.get_grid(trian::D2toD3AnalyticalMapCubedSphereTriangulation) = trian.model.cube_grid_geo

struct D2toD3AnalyticalMapCubedSphereDiscreteModel{T,B,C} <: Gridap.Geometry.DiscreteModel{2,3}
  cell_map::T
  cube_model_top::B
  cube_grid_geo::C
  function D2toD3AnalyticalMapCubedSphereDiscreteModel(
      cube_grid_geo::Gridap.Geometry.UnstructuredGrid{2,3},
      cube_grid_top::Gridap.Geometry.UnstructuredGrid{2,2};
      radius=1)
    m1=Fill(Gridap.Fields.GenericField(MapCubeToSphere(radius)),num_cells(cube_grid_geo))
    m2=get_cell_map(cube_grid_geo)
    m=lazy_map(∘,m1,m2)

    cube_model_top=Gridap.Geometry.UnstructuredDiscreteModel(cube_grid_top)

    # Build output object
    T=typeof(m)
    B=typeof(cube_model_top)
    C=typeof(cube_grid_geo)
    GC.gc()
    new{T,B,C}(m,cube_model_top,cube_grid_geo)
  end
end

# IMPORTANT NOTE: this method is needed as its default definition in Gridap
#   num_point_dims(model::DiscreteModel) = num_point_dims(get_grid_topology(model))
# returns Dp=2 as the topological dimension of the model is 2
Gridap.Geometry.num_point_dims(::D2toD3AnalyticalMapCubedSphereDiscreteModel) = 3
Gridap.Geometry.get_cell_map(model::D2toD3AnalyticalMapCubedSphereDiscreteModel) = model.cell_map
Gridap.Geometry.get_grid(model::D2toD3AnalyticalMapCubedSphereDiscreteModel) = model.cube_grid_geo
Gridap.Geometry.get_grid_topology(model::D2toD3AnalyticalMapCubedSphereDiscreteModel) = Gridap.Geometry.get_grid_topology(model.cube_model_top)
Gridap.Geometry.get_face_labeling(model::D2toD3AnalyticalMapCubedSphereDiscreteModel) = Gridap.Geometry.get_face_labeling(model.cube_model_top)
function Gridap.Geometry.Triangulation(a::D2toD3AnalyticalMapCubedSphereDiscreteModel)
  D2toD3AnalyticalMapCubedSphereTriangulation(a)
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# L1EpsilonInsLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(L1EpsilonInsLoss(1.5), c=:olive, linewidth=2, annotation=(0, 1.1, L"L\left( r \right) = max\{0, |r| - \epsilon \}", 15))
Plots.plot!(L1EpsilonInsLoss(1), c=:orangered4, linewidth=2.5)
Plots.plot!(L1EpsilonInsLoss(.5), c=:springgreen4, linewidth=2.5)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
function a(b::Int, c::Int)
    e = b*1.9
    f = c/0.5

    return e/f
end


text:
using DiffEqBase
using Test

mutable struct TestDiffEqOperator{T} <: DiffEqBase.AbstractDiffEqLinearOperator{T}
    m::Int
    n::Int
end

TestDiffEqOperator(A::AbstractMatrix{T}) where {T} =
    TestDiffEqOperator{T}(size(A)...)

Base.size(A::TestDiffEqOperator) = (A.m, A.n)

A = TestDiffEqOperator([0 0; 0 1])
B = TestDiffEqOperator([0 0 0; 0 1 0; 0 0 2])

@test_throws ErrorException AffineDiffEqOperator{Int64}((A,B),())


text:
"""
Package for simulating quantum error correction using stabilizer codes.
"""
module Qecsim

# imports
using Reexport

# exports
export QecsimError
include("error.jl")

# include core sub-modules (and reexport)
include("PauliTools.jl")
include("Model.jl")
include("App.jl")
@reexport using .PauliTools
@reexport using .Model
@reexport using .App

# include implementation sub-modules (not reexported)
include("models/BasicModels.jl")
include("models/GenericModels.jl")

end


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

mutable struct VirtualMachineScaleSetNetworkConfiguration <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, VirtualMachineScaleSetNetworkConfigurationProperties } # spec name: properties

    function VirtualMachineScaleSetNetworkConfiguration(;id=nothing, name=nothing, properties=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetNetworkConfiguration, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualMachineScaleSetNetworkConfiguration, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineScaleSetNetworkConfiguration, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type VirtualMachineScaleSetNetworkConfiguration

const _property_map_VirtualMachineScaleSetNetworkConfiguration = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_VirtualMachineScaleSetNetworkConfiguration = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("name")=>"String", Symbol("properties")=>"VirtualMachineScaleSetNetworkConfigurationProperties")
Base.propertynames(::Type{ VirtualMachineScaleSetNetworkConfiguration }) = collect(keys(_property_map_VirtualMachineScaleSetNetworkConfiguration))
Swagger.property_type(::Type{ VirtualMachineScaleSetNetworkConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetNetworkConfiguration[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetNetworkConfiguration }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetNetworkConfiguration[property_name]

function check_required(o::VirtualMachineScaleSetNetworkConfiguration)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetNetworkConfiguration }, name::Symbol, val)
end


text:
## This is until we have a good replacement for `info(io, msg)` on 0.7+
function info(io::IO, args...)
    printstyled(io, "Info: "; color=Base.info_color())
    printstyled(io, args...; color=Base.info_color())
    println(io)
end
function warn(io::IO, args...)
    printstyled(io, "Warning: "; color=Base.warn_color())
    printstyled(io, args...; color=Base.warn_color())
    println(io)
end


text:
function plot_poincare_error(t::Vector, I::Vector, filename; plot_title=nothing)
    if plot_title == nothing
        plot_title=L"\Delta I"
    end

    xf, yf = get_default_formatter(t)

    I_error = (I - I[1]) / I[1]

    fig = figure(figsize=(6,5))
    subplots_adjust(left=0.15, right=0.93, top=0.94, bottom=0.12)
    plot(t, I_error, ".", markersize=2)
    xticks(linspace(t[1], t[end], 6))
    xlim(t[1], t[end])
    xlabel(L"$t$")
    ylabel(plot_title)
    ax = gca()
    ax[:xaxis][:set_major_formatter](xf)
    ax[:yaxis][:set_major_formatter](yf)
    ax[:yaxis][:set_label_coords](-0.1,0.5)
    savefig(filename)
    close(fig)
end

function plot_poincare_error(t::Vector, I::Vector, K::Vector, Δt, filename; plot_title=nothing)
    if plot_title == nothing
        plot_title=L"\Delta I"
    end

    xf, yf = get_default_formatter(t)

    I_error  = zeros(I)
    I_error .= (I .- I[1]) .- Δt^2 .* (K .- K[1])
    I_error ./= (I[1] - Δt^2 * K[1])

    fig = figure(figsize=(6,5))
    subplots_adjust(left=0.15, right=0.93, top=0.94, bottom=0.12)
    plot(t, I_error, ".", markersize=2)
    xticks(linspace(t[1], t[end], 6))
    xlim(t[1], t[end])
    xlabel(L"$t$")
    ylabel(plot_title)
    ax = gca()
    ax[:xaxis][:set_major_formatter](xf)
    ax[:yaxis][:set_major_formatter](yf)
    ax[:yaxis][:set_label_coords](-0.1,0.5)
    savefig(filename)
    close(fig)
end

function plot_poincare_trajectories(sol::Solution, nplot, filename, dpi=100; xmin=-0.6, xmax=+0.6, ymin=-0.6, ymax=+0.6)
    if nplot > 0
        fig = figure(figsize=(10,8))
        subplots_adjust(left=0.0, right=0.7, top=1.0, bottom=0.0)
        ax  = gca(projection="3d")

        plot3D(sol.q.d[1,1,:], sol.q.d[2,1,:], sol.q.d[3,1,:])

        for i in 1:sol.ni
            if mod(i, div(sol.ni, nplot)) == 0
                plot3D(sol.q.d[1,:,i], sol.q.d[2,:,i], sol.q.d[3,:,i])
            end
        end

        xlabel(L"$x$", fontsize=18, labelpad=10)
        ylabel(L"$y$", fontsize=18, labelpad=10)
        zlabel(L"$z$", fontsize=18, labelpad=10)

        xlim(xmin, xmax)
        ylim(ymin, ymax)

        savefig(filename, dpi=dpi)
        close(fig)
    end
end

function plot_poincare_loop(sol::Solution, nplot, filename, dpi=100; xmin=-0.6, xmax=+0.6, ymin=-0.6, ymax=+0.6)
    if nplot > 0
        fig = figure(figsize=(10,8))
        subplots_adjust(left=0.0, right=0.7, top=1.0, bottom=0.0)
        ax  = gca(projection="3d")

        for i in 1:sol.nt+1
            if mod(i-1, div(sol.nt, nplot)) == 0
                # tpower      = log10(sol.t.t[end])
                # exponent    = @sprintf("%2d", tpower)
                # coefficient = @sprintf("%2.2f", sol.t.t[i] / 10^tpower)
                # tstr = "$coefficient \\times 10^\{ $exponent \}"
                # plot3D(sol.q.d[1,i,:], sol.q.d[2,i,:], sol.q.d[3,i,:], label="\$ t = $tstr \$")
                # coefficient = @sprintf("%4i", sol.t.t[i])
                # tstr = "$coefficient"
                tstr = string(sol.t.t[i])
                plot3D(sol.q.d[1,i,:], sol.q.d[2,i,:], sol.q.d[3,i,:], label="\$ t = $tstr \$")
            end
        end

        xlabel(L"$x$", fontsize=18, labelpad=10)
        ylabel(L"$y$", fontsize=18, labelpad=10)
        zlabel(L"$z$", fontsize=18, labelpad=10)

        xlim(xmin, xmax)
        ylim(ymin, ymax)

        legend(loc=2, bbox_to_anchor=(1.05, 0.8))
        savefig(filename, dpi=dpi)
        close(fig)
    end
end

function plot_poincare_surface(sol::Solution, nplot, filename, dpi=100; xmin=-0.15, xmax=+0.15, ymin=-0.15, ymax=+0.15)
    if nplot > 0
        fig = figure(figsize=(10,8))
        subplots_adjust(left=0.0, right=0.7, top=1.0, bottom=0.0)
        ax  = gca(projection="3d")

        nplot > sol.nt ? nplot = sol.nt : nothing

        for i in 0:sol.nt
            if mod(i, div(sol.nt, nplot)) == 0
                tpower      = log10(sol.t[end])
                exponent    = @sprintf("%2d", floor(Int64, tpower))
                coefficient = @sprintf("%2.2f", sol.t[i] / 10^floor(tpower))
                tstr = "$coefficient \\times 10^{ $exponent }"
                plot3D(sol.q.d[1,i+1,:], sol.q.d[2,i+1,:], sol.q.d[3,i+1,:], ".", label="\$ t = $tstr \$")
           end
        end

        xlabel(L"$x$", fontsize=18, labelpad=10)
        ylabel(L"$y$", fontsize=18, labelpad=10)
        zlabel(L"$z$", fontsize=18, labelpad=10)

        xlim(xmin, xmax)
        ylim(ymin, ymax)

        legend(loc=2, bbox_to_anchor=(1.05, 0.8))
        savefig(filename, dpi=dpi)
        close(fig)
    end
end


text:
using Jecco
using Plots
gr()

ts_multi = OpenPMDTimeSeries("./data")

it = 48

phi_1, grid_1 = get_field(ts_multi, it=it, field="phi c=1")
phi_2, grid_2 = get_field(ts_multi, it=it, field="phi c=2")
phi_3, grid_3 = get_field(ts_multi, it=it, field="phi c=3")
phi_4, grid_4 = get_field(ts_multi, it=it, field="phi c=4")

u_1, x, y = grid_1[:]
u_2, x, y = grid_2[:]
u_3, x, y = grid_3[:]
u_4, x, y = grid_4[:]

# heatmap(x, y, phi_1[1,:,:])
heatmap(x, y, phi_2[end,:,:])

# plot!(u_1, phi_1[:,10,10], marker=:circle)
# plot!(u_2, phi_2[:,10,10], marker=:circle)
# plot!(u_3, phi_3[:,10,10], marker=:circle)
# plot!(u_4, phi_4[:,10,10], marker=:circle)


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots, Plots.Measures
using HDF5
using DelimitedFiles
using Distributions
using StatsPlots

#Initialise variables
GeneIsoformRelationships = "data/gencode_human_gene_isoform_relationships.txt"
NumSimulations = 100
NumIsoformsToSimulate = 1
filteringThreshold = 4
alpha = nothing
beta = nothing
scnorm_counts = ["H1_24_scnorm_counts.txt",
"H1_96_scnorm_counts.txt",
"H9_24_scnorm_counts.txt",
"H9_96_scnorm_counts.txt"]
pFP = 0.01
pFN = 0.04
isoChoiceModels = ["Bernoulli", "Normal","Uniform"]

for countsMatrix in scnorm_counts

    output = "figures/figure6/" * countsMatrix

    for model in 1:length(isoChoiceModels)

        output = "figures/figure6/" * isoChoiceModels[model] * "_" * countsMatrix

        out = globalPredict("data/" * countsMatrix, GeneIsoformRelationships,
        NumSimulations, output, NumIsoformsToSimulate, filteringThreshold,
        false, alpha, beta, pFN, pFP, "none", isoChoiceModels[model],
        "figures/figure6_data/" * isoChoiceModels[model] * "_" * countsMatrix)
    end
end

###########################################################
# Make figures
function make_hist(path, titleArg)
    Arr = readdlm(path)

    if titleArg == "3 Isoforms" && path == "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100)

    elseif titleArg == "Real"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        xlabel = "Mean No. Isoforms\nper Gene per Cell")
    end
    return p1
end

function make_drop_plots(path, output)

    test = load(path, "data")

    arr1 = []
    arr2 = []
    arr3 = []
    arr4 = []

    for (key, value) in test

        new_val = sort(value)
        push!(arr1, new_val[1])
        push!(arr2, new_val[2])
        push!(arr3, new_val[3])
        push!(arr4, new_val[4])
    end

    p1 = histogram(arr1, normalize = :pdf, title = "Rank 1", xlim = (0,1), ylim = (0,10) , ylabel="Probability Density" )
    p2 = histogram(arr2, normalize = :pdf, title = "Rank 2", xlim = (0,1), ylim = (0,10) )
    p3 = histogram(arr3, normalize = :pdf, title = "Rank 3", xlim = (0,1), ylim = (0,10) , ylabel="Probability Density", xlabel = "p(Dropout)")
    p4 = histogram(arr4, normalize = :pdf, title = "Rank 4", xlim = (0,1), ylim = (0,10), xlabel = "p(Dropout)" )

    pSave = plot(p1,p2,p3,p4, layout = (2,2), legend = :false,
    size = (500 * 8,500 *8), axis = (font(120)), titlefontsize=120, margin = 40mm)
    savefig(pSave, output)
    return pSave
end

function make_expr_plots(path, output)

    test = load(path, "data")

    arr1 = []
    arr2 = []
    arr3 = []
    arr4 = []

    for (key, value) in test

        new_val = sort(value, rev = true)
        push!(arr1, new_val[1])
        push!(arr2, new_val[2])
        push!(arr3, new_val[3])
        push!(arr4, new_val[4])
    end

    p1 = histogram(log.(arr1 .+ 1), normalize = :pdf, title = "Rank 1", xlim = (0,10), ylim = (0,2), ylabel="Probability Density" )
    p2 = histogram(log.(arr2 .+ 1), normalize = :pdf, title = "Rank 2", xlim = (0,10), ylim = (0,2) )
    p3 = histogram(log.(arr3 .+ 1), normalize = :pdf, title = "Rank 3", xlim = (0,10), ylim = (0,2), xlabel="log(counts + 1)", ylabel="Probability Density" )
    p4 = histogram(log.(arr4 .+ 1), normalize = :pdf, title = "Rank 4", xlim = (0,10), ylim = (0,2), xlabel="log(counts + 1)")

    pSave = plot(p1,p2,p3,p4, layout = (2,2), legend = :false,
    size = (500 * 8,500 *8), axis = (font(120)), titlefontsize=120, margin = 40mm)
    savefig(pSave, output)
    return pSave
end

normal_hist = make_hist("figures/figure6_data/Normal_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")
normal_hist = plot(normal_hist, ylabel = "Probability Density")
bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")
uniform_hist = make_hist("figures/figure6_data/Uniform_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")

bernoulli_dist = bar(Bernoulli(0.25), bar_width = 0.1, ylim =(0,1), legend = false, title = "Bernoulli", xticks = 0:0.5:1, xlabel = "p(choice)")
normal_dist = plot(TruncatedNormal(0.25, 0.06, 0.0, 1.0), legend = false, title = "Normal", linestyle=:solid, linewidth=2*8, xticks = 0:0.5:1, ylabel = "Probability Density", xlabel = "p(choice)")
uniform_dist = histogram([ 0.25], ylim = (0,1), legend = false, title = "p=0.25", xlim = (0,1), bins = 0:0.1:1, xticks = 0:0.5:1, xlabel = "p(choice)")

p_A = plot(normal_dist, bernoulli_dist, uniform_dist,
normal_hist, bernoulli_hist, uniform_hist,
layout = (2,3), legend = false,
size = (1000 * 8,500 *8), axis = (font(120)), titlefontsize=120, margin = 80mm)

savefig(p_A, "figures/figure6/different_dists.png")

expr_H1_24 = make_expr_plots("figures/pdrop_scnorm_H1_24.pngexpr_dict.jld", "figures/figure6/H1_24_expr.png")
expr_H1_96 = make_expr_plots("figures/pdrop_scnorm_H1_96.pngexpr_dict.jld", "figures/figure6/H1_96_expr.png")
expr_H9_24 = make_expr_plots("figures/pdrop_scnorm_H9_24.pngexpr_dict.jld", "figures/figure6/H9_24_expr.png")
expr_H9_96 = make_expr_plots("figures/pdrop_scnorm_H9_96.pngexpr_dict.jld", "figures/figure6/H9_96_expr.png")

drop_H1_24 = make_drop_plots("figures/pdrop_scnorm_H1_24.png_dict.jld", "figures/figure6/H1_24_drop.png")
drop_H1_96 = make_drop_plots("figures/pdrop_scnorm_H1_96.png_dict.jld", "figures/figure6/H1_96_drop.png")
drop_H9_24 = make_drop_plots("figures/pdrop_scnorm_H9_24.png_dict.jld", "figures/figure6/H9_24_drop.png")
drop_H9_96 = make_drop_plots("figures/pdrop_scnorm_H9_96.png_dict.jld", "figures/figure6/H9_96_drop.png")

#make a supplementary fig
H1_24_normal_hist = make_hist("figures/figure6_data/Normal_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "")
H1_24_normal_hist = plot(H1_24_normal_hist, ylabel = "Probability Density")
H1_24_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "")
H1_24_uniform_hist = make_hist("figures/figure6_data/Uniform_H1_24_scnorm_counts.txt_1_isoforms_array.txt", "")

H9_24_normal_hist = make_hist("figures/figure6_data/Normal_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "")
H9_24_normal_hist = plot(H9_24_normal_hist, ylabel = "Probability Density")
H9_24_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "")
H9_24_uniform_hist = make_hist("figures/figure6_data/Uniform_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "")

H9_96_normal_hist = make_hist("figures/figure6_data/Normal_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "")
H9_96_normal_hist = plot(H9_96_normal_hist, ylabel = "Probability Density")
H9_96_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "")
H9_96_uniform_hist = make_hist("figures/figure6_data/Uniform_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "")

H1_96_normal_hist = make_hist("figures/figure6_data/Normal_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")
H1_96_normal_hist = plot(H1_96_normal_hist, ylabel = "Probability Density")
H1_96_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")
H1_96_uniform_hist = make_hist("figures/figure6_data/Uniform_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "")

p_supp = plot(normal_dist, bernoulli_dist, uniform_dist,
H1_96_normal_hist, H1_96_bernoulli_hist, H1_96_uniform_hist,
H1_24_normal_hist, H1_24_bernoulli_hist, H1_24_uniform_hist,
H9_96_normal_hist, H9_96_bernoulli_hist, H9_96_uniform_hist,
H9_24_normal_hist, H9_24_bernoulli_hist, H9_24_uniform_hist,
layout = (5,3), legend = false,
size = (1000 * 8,1200 *8), axis = (font(120)), titlefontsize=120, margin = 80mm)

savefig(p_supp, "figures/supplementary_figs/different_dists.png")

###############################################################
using JLD2
using DelimitedFiles
using HypothesisTests
using DataFrames

df = DataFrame(IsoformsSimulated=Int64[], Comparison = String[], pValue = Float64[])

normal_H1_96 = readdlm("figures/figure6_data/Normal_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
bernoulli_H1_96 = readdlm("figures/figure6_data/Bernoulli_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
uniform_hist_H1_96 = readdlm("figures/figure6_data/Uniform_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["H1_96"], pValue=[pvalue(KSampleADTest(vec(normal_H1_96), vec(bernoulli_H1_96), vec(uniform_hist_H1_96)))]))

normal_H1_24 = readdlm("figures/figure6_data/Normal_H1_24_scnorm_counts.txt_1_isoforms_array.txt")
bernoulli_H1_24 = readdlm("figures/figure6_data/Bernoulli_H1_24_scnorm_counts.txt_1_isoforms_array.txt")
uniform_hist_H1_24 = readdlm("figures/figure6_data/Uniform_H1_24_scnorm_counts.txt_1_isoforms_array.txt")
append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["H1_24"], pValue=[pvalue(KSampleADTest(vec(normal_H1_24), vec(bernoulli_H1_24), vec(uniform_hist_H1_24)))]))

normal_H9_96 = readdlm("figures/figure6_data/Normal_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
bernoulli_H9_96 = readdlm("figures/figure6_data/Bernoulli_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
uniform_hist_H9_96 = readdlm("figures/figure6_data/Uniform_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["H9_96"], pValue=[pvalue(KSampleADTest(vec(normal_H9_96), vec(bernoulli_H9_96), vec(uniform_hist_H9_96)))]))

normal_H9_24 = readdlm("figures/figure6_data/Normal_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
bernoulli_H9_24 = readdlm("figures/figure6_data/Bernoulli_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
uniform_hist_H9_24 = readdlm("figures/figure6_data/Uniform_H9_24_scnorm_counts.txt_1_isoforms_array.txt")
append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["H9_24"], pValue=[pvalue(KSampleADTest(vec(normal_H9_24), vec(bernoulli_H9_24), vec(uniform_hist_H9_24)))]))

jldopen("figures/figure6_data/Supfig_AD_test_results.jld2", "w") do file
    file["df"] = df
end

#############################################
# Overlap
#############################################

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if titleArg == "3 Isoforms" && path == "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1),
        ylims = (0, 10),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 50,
        color = palette(:default)[2],)

    elseif titleArg == "Real"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1),
        ylims = (0, 10),
        normalize=:pdf,
        xlabel = "Overlap With Ground Truth",
        bins = 50,
        color = palette(:default)[2],)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1),
        ylims = (0, 10),
        normalize=:pdf,
        bins = 50,
        color = palette(:default)[2],
        xlabel = "Overlap With Ground Truth")
    end
    return p1
end

#make a supplementary fig
H1_24_normal_hist = make_hist("figures/figure6_data/Normal_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H1_24_normal_hist = plot(H1_24_normal_hist, ylabel = "Probability Density")
H1_24_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H1_24_uniform_hist = make_hist("figures/figure6_data/Uniform_H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")

H9_24_normal_hist = make_hist("figures/figure6_data/Normal_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H9_24_normal_hist = plot(H9_24_normal_hist, ylabel = "Probability Density")
H9_24_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H9_24_uniform_hist = make_hist("figures/figure6_data/Uniform_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")

H9_96_normal_hist = make_hist("figures/figure6_data/Normal_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H9_96_normal_hist = plot(H9_96_normal_hist, ylabel = "Probability Density")
H9_96_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H9_96_uniform_hist = make_hist("figures/figure6_data/Uniform_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")

H1_96_normal_hist = make_hist("figures/figure6_data/Normal_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H1_96_normal_hist = plot(H1_96_normal_hist, ylabel = "Probability Density")
H1_96_bernoulli_hist = make_hist("figures/figure6_data/Bernoulli_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")
H1_96_uniform_hist = make_hist("figures/figure6_data/Uniform_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "")

p_supp = plot(normal_dist, bernoulli_dist, uniform_dist,
H1_96_normal_hist, H1_96_bernoulli_hist, H1_96_uniform_hist,
H1_24_normal_hist, H1_24_bernoulli_hist, H1_24_uniform_hist,
H9_96_normal_hist, H9_96_bernoulli_hist, H9_96_uniform_hist,
H9_24_normal_hist, H9_24_bernoulli_hist, H9_24_uniform_hist,
layout = (5,3), legend = false,
size = (1000 * 8,1200 *8), axis = (font(120)), titlefontsize=120, margin = 80mm)

savefig(p_supp, "figures/supplementary_figs/different_dists_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Geometry{Dim,T}

A geometry embedded in a `Dim`-dimensional space with coordinates of type `T`.
"""
abstract type Geometry{Dim,T} end

"""
    embeddim(geometry)

Return the number of dimensions of the space where the `geometry` is embedded.
"""
embeddim(::Type{<:Geometry{Dim,T}}) where {Dim,T} = Dim
embeddim(g::Geometry) = embeddim(typeof(g))

"""
    paramdim(geometry)

Return the number of parametric dimensions of the `geometry`. For example, a
sphere embedded in 3D has 2 parametric dimension (polar and azimuthal angles).
"""
paramdim(g::Geometry) = paramdim(typeof(g))

"""
    coordtype(geometry)

Return the machine type of each coordinate used to describe the `geometry`.
"""
coordtype(::Type{<:Geometry{Dim,T}}) where {Dim,T} = T
coordtype(g::Geometry) = coordtype(typeof(g))

"""
    measure(geometry)

Return the measure of the `geometry`, i.e. the length, area, or volume.
"""
function measure end

"""
    boundary(geometry)

Return the boundary of the `geometry`.
"""
function boundary end

# ----------------
# IMPLEMENTATIONS
# ----------------
include("primitives.jl")
include("polytopes.jl")

"""
    GeometrySet(geometries)

A set of `geometries` seen as a single geometry.

## Example

In a geographic map, countries can be described with
multiple polygons (a.k.a. MultiPolygon).
"""
struct GeometrySet{Dim,T,G<:Geometry{Dim,T}} <: Geometry{Dim,T}
  geometries::Vector{G}
end


text:
using Base: AbstractFloat
########################################################################################################
## Details
# This file will create the type structures for the UQ inference problems.
#
#########################################################################################################

abstract type ArbitraryDistribution end

mutable struct PotentialDistribution <: ArbitraryDistribution
    potential::Potentials.ArbitraryPotential          # potential defined in Potential.jl
    prior :: ContinuousMultivariateDistribution
    distribution::Function
    x::NamedTuple
    t::TransformVariables.AbstractTransform
    samples::Vector
    function PotentialDistribution(potential::Potentials.ArbitraryPotential, x::NamedTuple, prior::ContinuousMultivariateDistribution, likelihood::Function)
        t = namedtp_to_vec(x)
        distribution(x) = likelihood(x) + prior.logpdf(x)
        new(potential, prior, distribution, x, t, zeros(0))
    end
end

function (potential_dist::PotentialDistribution)(θ)
    x = inverse(potential_dist.t, θ)
    potential_dist.distribution(x) + potential_dist.prior.logpdf(x)
end

#########################################################################################################
########################################### Prior #######################################################
#########################################################################################################

# Define an ArbitraryPrior type
abstract type ArbitraryPrior <: ContinuousMultivariateDistribution end

function Base.length(d::ArbitraryPrior)
    return length(d.x)
end

Distributions.rand(rng::AbstractRNG, d::ArbitraryPrior) = d.sampler(rng, d)
function Distributions._rand!(rng::AbstractRNG, d::ArbitraryPrior, x::AbstractVector)
    x = d.sampler(rng, d)
end

Distributions._logpdf(d::ArbitraryPrior, x::AbstractArray) = d.logpdf(x)
Bijectors.bijector(d::ArbitraryPrior) = Bijectors.Identity{1}()

########################################### LJ #######################################################
struct LJ_Prior_Distribution <: ArbitraryPrior
    y_train :: Vector{Potentials.Configuration}
    x       :: Vector{Float64}
    t       :: TransformVariables.AbstractTransform
    sampler :: Function
    logpdf  :: Function
end

Bijectors.bijector(d::LJ_Prior_Distribution) = Bijectors.Logit{1}(0.5, 1.5)

########################################### SNAP  #######################################################
struct SNAP_Prior_Distribution <: ArbitraryPrior
    A :: Array
    b :: Vector
    x :: Vector
    t :: TransformVariables.AbstractTransform
    sampler :: Function
    logpdf :: Function
end
function SNAP_Prior_Distribution(A::Array, b::Vector, v::NamedTuple)
    t = namedtp_to_vec(v)
    x = inverse(t, v)
    sampler(rng::AbstractRNG, d::ArbitraryPrior) = randn(length(d))
    function logpdf(θ::Vector)
        n = length(θ)
        loglikelihood( MvNormal( zeros(n), I(n) ), θ )
    end
    SNAP_Prior_Distribution(A, b, x, t, sampler, logpdf)
end

#####################################################################################################


text:
export trouve_intervalleASecA
function trouve_intervalleASecA(h :: LineModel,
                                t₀ :: Float64,
                                tₘ :: Float64;
                                ϵ :: Float64=1e-10,
                                verbose :: Bool=false)

        h₀=obj(h,t₀)
        dh₀=grad(h,t₀)

        γ=0.8
        t=tₘ
        tp=t₀
        tqnp=t₀
        hk=0
        hkm1=0
        gkm1=0
        hp=0
        i=0

        hk=obj(h,t)
        gk=grad(h,t)
        hkm1=obj(h,tqnp)
        gkm1=grad(h,tqnp)

        verbose && @printf("iter tp        gkm1        hkm1         t        gk        hk\n")
        verbose && @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n", i, tp,gkm1,hkm1,t,gk,hk)

        while  i<50
          if (hk>h₀+0.01*(t-t₀)*dh₀) ||((hk>hkm1) & (i>1))
            (topt,iter)=zoom_SecA(h,tp,t,verbose=false)
            return (topt,iter)
          end

          if (abs(gk)<=ϵ)
            iter=i
            topt=t
            return (topt,iter)
          end

          if (gk>=0)
            (topt,iter)=zoom_SecA(h,t,tp,verbose=false)
            return (topt,iter)
         end

          s = t-tqnp
          y = gk-gkm1

          Γ=3*(gk+gkm1)*s-6*(hk-hkm1)
          if y*s+Γ < eps(Float64)*(s^2)
            yt=y
          else
            yt=y+Γ/s
          end

          dS=-gk*s/yt

          if ((tp-t)*dS>0) & (dS/(tp-t)<γ)
            tplus = t + dS
            hplus = obj(h, tplus)
            gplus = grad(h,tplus)
            verbose && println("N")
          else
            tplus = (t+tp)/2
            hplus = obj(h, tplus)
            gplus = grad(h,tplus)
            verbose && println("B")
          end

          if t>tp
            if gplus<0
              tp=t
              tqnp=t
              t=tplus
            else
              tqnp=t
              t=tplus
            end
          else
            if gplus>0
              tp=t
              tqnp=t
              t=tplus
            else
              tqnp=t
              t=tplus
            end
          end

          #mise à jour des valeurs
          hkm1=hk
          gkm1=gk
          hk=hplus
          gk=gplus
          i=i+1
          verbose && @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n", i, tim1,gkm1,hkm1,ti,gk,hk)
        end
end


text:
#importall CUTEst
#using SparseArrays
#import CUTEst

export get_original_x, get_y

mutable struct Class_bounds
    l_i::Array{Int64,1}
    u_i::Array{Int64,1}
    l::Array{Float64,1}
    u::Array{Float64,1}

    function Class_bounds(lb_vec::Array{Float64,1}, ub_vec::Array{Float64,1})
        @assert(length(lb_vec) == length(ub_vec))
        this = new([],[],[],[]);
        for i = 1:length(lb_vec)
            if lb_vec[i] > -Inf
               push!(this.l_i, i)
               push!(this.l, lb_vec[i])
            end

            if ub_vec[i] < Inf
               push!(this.u_i, i)
               push!(this.u, ub_vec[i])
            end
        end

        return this;
    end
end

function _i_not_fixed(m::NLPModels.AbstractNLPModel)
    return (1:m.meta.nvar)[m.meta.lvar .!= m.meta.uvar]
end

mutable struct Class_CUTEst <: abstract_nlp
    nlp::NLPModels.AbstractNLPModel

    bcon::Class_bounds
    bvar::Class_bounds

    function Class_CUTEst(nlp::NLPModels.AbstractNLPModel)
        ind = _i_not_fixed(nlp)
        return new(nlp, Class_bounds(nlp.meta.lcon[:], nlp.meta.ucon[:]), Class_bounds(nlp.meta.lvar[ind], nlp.meta.uvar[ind]))
    end
end

#=
function con_info(m::Class_CUTEst, index::Int64)
    con = false
    if index in cons_indicies(m)
      con = true
      if index < length(nlp.bcon.l_i)
        #[lb(a, m.bcon); ub(a, m.bcon); lb(x, m.bvar); ub(x, m.bvar)]
      end
    end

    return;
end=#

function suggested_starting_point(m::Class_CUTEst)
    ind = _i_not_fixed(m.nlp)
    return deepcopy(m.nlp.meta.x0[ind])
end

function ncon(m::Class_CUTEst)
    return nbounds_orginal(m) + ncons_orginal(m)
end

function lb(x::Array{Float64,1}, bd::Class_bounds)
    return x[bd.l_i] - bd.l
end

function ub(x::Array{Float64,1}, bd::Class_bounds)
    return bd.u - x[bd.u_i]
end

function linear_cons(m::Class_CUTEst)
    is_linear = zeros(m.nlp.meta.ncon)
    is_linear[m.nlp.meta.lin] .= 1.0
    vec = [is_linear[m.bcon.l_i]; is_linear[m.bcon.u_i]; ones(nbounds_orginal(m))]

    return 1.0 .== vec
end

function ineq_cons(m::Class_CUTEst)
    is_ineq = zeros(m.nlp.meta.ncon)
    is_ineq[m.nlp.meta.lcon .== m.nlp.meta.ucon] .= 1.0
    vec = [is_ineq[m.bcon.l_i]; is_ineq[m.bcon.u_i]; ones(nbounds_orginal(m))]

    return 1.0 .== vec
end

function nbounds_orginal(nlp::Class_CUTEst)
    return length(nlp.bvar.l_i) + length(nlp.bvar.u_i)
end

function ncons_orginal(nlp::Class_CUTEst)
    return length(nlp.bcon.l_i) + length(nlp.bcon.u_i)
end

function cons_indicies(nlp::Class_CUTEst)
    m = ncons_orginal(nlp)
    if m > 0
      return 1:m
    else
      return []
    end
end

#function linear_indicies(nlp::Class_CUTEst)
#    #nlp.nlp.meta.lin
#end

function bound_indicies(nlp::Class_CUTEst)
    m = ncons_orginal(nlp)
    r = nbounds_orginal(nlp)
    if r > 0
      return (m + 1):(m + r)
    else
      return [];
    end
end

function y_l_con(y::Array{Float64,1}, m::Class_CUTEst)
    return y[1:length(m.bcon.l)]
end

function y_u_con(y::Array{Float64,1}, m::Class_CUTEst)
    n_lcon = length(m.bcon.l)
    return y[(n_lcon + 1):(n_lcon + length(m.bcon.u))]
end

function eval_f(m::Class_CUTEst, x::Array{Float64,1})
    return obj(m.nlp, _cute_x(m, x) )
end

function eval_a(m::Class_CUTEst, x::Array{Float64,1})
    a = cons(m.nlp, _cute_x(m, x) )
    if a == Float64[]
      a = zeros(1)
    end
    return [lb(a, m.bcon); ub(a, m.bcon); lb(x, m.bvar); ub(x, m.bvar)];
end

function _cute_x(m::Class_CUTEst, x::Array{Float64,1})
    ind = _i_not_fixed(m.nlp)
    #@show length(ind)
    if(length(x) != length(ind))
      error("$(length(x)) = length(x) != length(i_not_fixed) = $(length(ind))")
    end

    cute_x = deepcopy(m.nlp.meta.lvar) # get correct values of fixed variables
    cute_x[ind] = x

    return cute_x
end

function eval_jac(m::Class_CUTEst, x::Array{Float64,1})
    cute_x = _cute_x(m, x)
    J_full_T = jac(m.nlp, cute_x)'
    J_T = J_full_T[_i_not_fixed(m.nlp),:];
    #my_eye = SparseArrays.speye(length(x))
    my_eye = SparseMatrixCSC{Float64}(LinearAlgebra.I, length(x), length(x))
    Q_T = [J_T[:,m.bcon.l_i] -J_T[:,m.bcon.u_i] my_eye[:,m.bvar.l_i] -my_eye[:,m.bvar.u_i]];
    return SparseArrays.sparse(Q_T')

    #return @time [J[m.bcon.l_i,:]; -J[m.bcon.u_i,:]; my_eye[m.bvar.l_i,:]; -my_eye[m.bvar.u_i,:]];
end

function eval_grad_f(m::Class_CUTEst, x::Array{Float64,1})
    return grad(m.nlp, _cute_x(m, x))[_i_not_fixed(m.nlp)]
end

function get_constrduals(m::Class_CUTEst, y::Array{Float64,1})
     constrdual = zeros(m.nlp.meta.ncon)
     st_u = length(m.bcon.l) #+ ncons_orginal(m)
     constrdual[m.bcon.l_i] += y[1:st_u]
     constrdual[m.bcon.u_i] -= y[(st_u+1):(length(m.bcon.u) + st_u)]
     return constrdual
end

function get_reducedcosts(m::Class_CUTEst, y::Array{Float64,1})
    rc = zeros(m.nlp.meta.nvar)
    st_l = ncons_orginal(m)
    st_u = length(m.bvar.l) + ncons_orginal(m)
    rc[m.bvar.l_i] += y[(st_l+1):st_u]
    rc[m.bvar.u_i] -= y[(st_u+1):(length(m.bvar.u) + st_u)]
    return rc
end

function y_cons_net(m::Class_CUTEst, y::Array{Float64,1})
    return;
end

function eval_jtprod(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1})
    return jtprod(nlp_raw, x, y)
end

function eval_schur_and_jtprod(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1}, s::Array{Float64,1})

end

#=function ncon(m::Class_CUTEst)
    return length(m.bcon.l) + length(m.bcon.u)
end=#

#=function nvar(m::Class_CUTEst)
    return m.nlp.meta.nvar
end=#
#=
function make_symmetric(M::SparseMatrixCSC{Float64,Int32})
    n = size(M,1)
    new_M = spzeros(n,n)
    rows = rowvals(M)
    vals = nonzeros(M)
    for col = 1:n
      for j in nzrange(M, col)
         row = rows[j]
         val = vals[j]
         # perform sparse wizardry...
         new_M[col,row] = vals[j]
         new_M[row,col] = vals[j]
      end
    end

    return new_M
end=#

function eval_lag_hess(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1}, w::Float64)
    y_cons = zeros(m.nlp.meta.ncon)
    y_cons[m.bcon.l_i] -= y_l_con(y, m)
    y_cons[m.bcon.u_i] += y_u_con(y, m)

    #H = hess(m.nlp, _cute_x(m, x), obj_weight=w, y=y_cons);
    #H = hess(m.nlp, _cute_x(m, x), y=y_cons; obj_weight=w);
    H = hess(m.nlp, _cute_x(m, x), y_cons; obj_weight=w);
    #H = hess(m.nlp, _cute_x(m, x), w, y_cons);
    #@show typeof(H)

    ind = _i_not_fixed(m.nlp)

    H_not_fixed = H[ind,ind]
    #H_true = SparseArrays.sparse(Symmetric(H_not_fixed,:L) + spzeros(length(x),length(x)))
    #H_true = H_not_fixed + H_not_fixed' - spdiagm(diag(H_not_fixed))
    #convert()
    #@time H_true = make_symmetric(H_not_fixed)
    #@show LinearAlgebra.norm(H_true - H2,Inf)

    return H_not_fixed
end

function eval_Jt_prod(m::Class_CUTEst, x::Array{Float64,1}, y::Array{Float64,1})
    return jtprod(m.nlp, x, y)
end


text:
# 单独保存一个.jl文件（头文件），和刚才的notebook放在同一路径下
module FancyStructs
  struct Coordinate{T}
    x::T
    y::T
    z::T
  end
  function CoordinateAdd(obj1::FancyStructs.Coordinate, obj2::FancyStructs.Coordinate)
    println("(", obj1.x, ",", obj1.y, "," , obj1.z,")",  " + ", "(", obj2.x, ",", obj2.y, "," , obj2.z,")",
        " = ", "(",obj1.x + obj2.x, ",", obj1.y + obj2.y, "," , obj1.z + obj2.z,")"  )
  end
end


text:
import Base.print

#=
print : print Assoc in tabular form.
=#
function print(A::Assoc)
    dump(full(A))
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
isdefined(Base, :__precompile__) && __precompile__()

module TexExtensions

using Compat

@compat import Base.show

@Base.Multimedia.textmime "text/mathtex+latex"

const tm = MIME("text/mathtex+latex")
const T = MIME"text/mathtex+latex"

@compat function show(io::IO, ::T, x::AbstractString)
    write(io,"\\text{")
    write(io,x)
    write(io,"}")
end

@compat function show(io::IO, ::T, x::Rational)
    write(io,"\\frac{")
    show(io,tm,x.num)
    write(io,"}{")
    show(io,tm,x.den)
    write(io,"}")
end

# This code is largely copied from grisu.jl, we really should have to copy this mutch

_writefloat(io::IO, x::AbstractFloat, mode::Int32, n::Int, t) =
    _writefloat(io, x, mode, n, t, "NaN", "Inf")
_writefloat(io::IO, x::Float32, mode::Int32, n::Int, t) =
    _writefloat(io, x, mode, n, t, "NaN32", "Inf32")
_writefloat(io::IO, x::Float16, mode::Int32, n::Int, t) =
    _writefloat(io, x, mode, n, t, "NaN16", "Inf16")

function _writefloat(io::IO, x::AbstractFloat, mode::Int32, n::Int, typed, nanstr, infstr)
    if isnan(x) return write(io, "\\text{", typed ? nanstr : "NaN","}"); end
    if isinf(x)
        if x < 0 write(io,"\{\}-") end
        write(io, "\\infty")
        return
    end
    if typed && isa(x,Float16) write(io, "\\text{float16("); end
    len, pt, neg, pdigits = grisu_fmt(x, mode, n)
    if mode == Base.Grisu.PRECISION
        while len > 1 && Base.Grisu.DIGITS[len] == '0'
            len -= 1
        end
    end
    # Make sure latex understands this is negations
    if neg write(io,"\{\}-") end
    if pt <= -4 || pt > 6 # .00001 to 100000.
        # => #.#######e###
        write(io, pdigits, 1)
        write(io, '.')
        if len > 1
            write(io, pdigits+1, len-1)
        else
            write(io, '0')
        end
        write(io, "\\times 10^{",dec(pt-1),"}")
        if typed && isa(x,Float16)
            write(io, ")\\}")
        end
        return
    elseif pt <= 0
        # => 0.00########
        write(io, "0.")
        while pt < 0
            write(io, '0')
            pt += 1
        end
        write(io, pdigits, len)
    elseif pt >= len
        # => ########00.0
        write(io, pdigits, len)
        while pt > len
            write(io, '0')
            len += 1
        end
        write(io, ".0")
    else # => ####.####
        write(io, pdigits, pt)
        write(io, '.')
        write(io, pdigits+pt, len-pt)
    end
    if typed && isa(x,Float16) write(io, ")\\}"); end
    nothing
end

if VERSION < v"0.4-dev"
    eval(quote
        function grisu_fmt(x, mode, n)
            @Base.Grisu.grisu_ccall x mode n
            pdigits = pointer(Base.Grisu.DIGITS)
            neg = Base.Grisu.NEG[1]
            len = int(Base.Grisu.LEN[1])
            pt  = int(Base.Grisu.POINT[1])
            len, pt, neg, pdigits
        end
    end)
else
    grisu_fmt(x, mode, n) = Base.Grisu.grisu(x, mode, n)
end

@compat function show(io::IO, ::T, x::BigFloat)
    e = Array(Culong,1)
    str = ccall((:mpfr_get_str,:libmpfr),Ptr{UInt8},(Ptr{UInt8},Ptr{Culong},Cint,Csize_t,Ptr{Void},Int32),C_NULL,e,
        10,0,&x,Base.MPFR.ROUNDING_MODE[end])
    ret = bytestring(str)
    ccall((:mpfr_free_str,:libmpfr),Void,(Ptr{UInt8},),str)
    if ret[1] == '-'
        radix = 2
    else
        radix = 1
    end
    write(io,ret[1:radix])
    write(io,'.')
    write(io,ret[(radix+1):end])
    write(io, "\\times 10^{",dec(e[1]),"}")
end

@compat show(io::IO, ::T, x::Float64) = _writefloat(io, x, Base.Grisu.SHORTEST, 0, true)
@compat show(io::IO, ::T, x::Float32) = _writefloat(io, x, Base.Grisu.SHORTEST_SINGLE, 0, true)
@compat show(io::IO, ::T, x::Float16) = _writefloat(io, x, Base.Grisu.PRECISION, 4, true)

@compat show(io::IO, ::T, x::Integer) = show(io,x)

end # module


text:
function defaultLPsolverfor{N,T}(p::Rep{N,T})
  if vrepiscomputed(p)
    SimpleVRepPolyhedraModel{N,T}()
  else
    MathProgBase.defaultLPsolver
  end
end


text:
module Benchmarks

    using Compat

    export @benchmark

    include("01_clock_resolution.jl")
    include("02_environment.jl")
    include("03_samples.jl")
    include("04_results.jl")
    include("benchmarkable.jl")
    include("ols.jl")
    include("execute.jl")
    include("benchmark.jl")
    include("wrapper_types.jl")
end


text:
#=
These functions are needed for specific hardware implementations
=#

#=
Audio
=#

set_audio(fp::Array{FPGA,1},h::Gui_Handles,r::RHD2000)=set_audio_fpga(fp,h.sc.spike,r.refs[h.sc.spike])

function set_audio(fpga::DArray{FPGA,1,Array{FPGA,1}},han::Gui_Handles,rhd::RHD2000)

    remotecall_wait(((x,h,ii)->set_audio_fpga(localpart(x),h,ii)),2,fpga,han.sc.spike,rhd.refs[han.sc.spike])

end

function set_audio_fpga(fpga::Array,ii,refs)

    if fpga[1].amps[1] != 255
        selectDacDataStream(fpga[1],0,div(ii-1,32))
        selectDacDataChannel(fpga[1],0,rem(ii-1,32))

        if refs>0
            enableDac(fpga[1],1,true)
            selectDacDataStream(fpga[1],1,div(refs-1,32))
            selectDacDataChannel(fpga[1],1,rem(refs-1,32))
        else
            enableDac(fpga[1],1,false)
        end
    end
    nothing
end

#=
Bandwidth Adjustment
=#

change_bandwidth(fpgas::Array{FPGA,1},lower,upper,dsp_lower)=change_bandwidth_fpga(fpgas,lower,upper,dsp_lower)

change_bandwidth(fpgas::DArray{FPGA,1,Array{FPGA,1}},lower,upper,dsp_lower)=change_bandwidth_fpga(fpgas,lower,upper,dsp_lower)

function change_bandwidth_fpga(fpgas::Array,lower,upper,dsp_lower)

    for fpga in fpgas
        setLowerBandwidth(lower,fpga.r)
        setUpperBandwidth(upper,fpga.r)
        setDspCutoffFreq(dsp_lower,fpga.r)
        commandList=createCommandListRegisterConfig(zeros(Int32,1),true,fpga.r)
        uploadCommandList(fpga,commandList, "AuxCmd3", 1)
    end
    nothing
end

function change_bandwidth_fpga(fpgas::DArray,lower,upper,dsp_lower)

    @sync begin
        for p in procs(fpgas)
            @async remotecall_wait((ff)->change_bandwidth(localpart(ff),lower,upper,dsp_lower),p,fpgas)
        end
    end
    nothing
end

#=
Threshold
=#

function send_thres_to_ic(han::Gui_Handles,fpga::Array{FPGA,1})
    nothing
end

function send_thres_to_ic(han::Gui_Handles,fpga::DArray{FPGA,1,Array{FPGA,1}})
    nothing
end

#=
Select New Channel
=#

new_single_channel(han::Gui_Handles,rhd::RHD2000,s,fpga::Array{FPGA,1})=new_single_channel_fpga(han,rhd,s,fpga)

new_single_channel(han::Gui_Handles,rhd::RHD2000,s,fpga::DArray{FPGA,1,Array{FPGA,1}})=new_single_channel_fpga(han,rhd,s,fpga)

function new_single_channel_fpga(han::Gui_Handles,rhd::RHD2000,s,fpga)

    han.total_clus[han.sc.spike] = han.sc.total_clus

    han.sc.spike=han.chan_per_display*han.num16-han.chan_per_display+han.num

    han.sc.total_clus = han.total_clus[han.sc.spike]

    SpikeSorting.clear_c2(han.sc.c2,han.sc.spike)
    han.sc.ctx2=Gtk.getgc(han.sc.c2)
    han.sc.ctx2s=copy(han.sc.ctx2)

    #Audio output
    set_audio(fpga,han,rhd)

    #Display Gain
    setproperty!(han.sc.gain_widgets.gainbox,:value,round(Int,han.scale[han.sc.spike,1]*-1000))

    #Display Threshold
    get_thres(han,s)

    han.buf.ind=1
    han.buf.count=1

    #Get Cluster
    get_cluster(han,s)

    #Update treeview
    SpikeSorting.update_treeview(han.sc)

    #update selected cluster
    SpikeSorting.select_unit(han.sc)

    #Sort Button
    if han.sc.sort_cb
        SpikeSorting.draw_templates(han.sc)
    end

    han.spike_changed=false

    nothing
end


text:
import ForwardDiff: valtype
const ∇ = Tensors.gradient
const Δ = Tensors.hessian

function Ψ(C, μ, Kb)
    T = ForwardDiff.valtype(eltype(C))
    detC = det(C)
    J = sqrt(detC)
    Ĉ = detC^(convert(T, -1 / 3)) * C
    return (convert(T, μ) * (tr(Ĉ) - 3) / 2 + convert(T, Kb) * (J - 1)^2)
end

function S(C, μ, Kb)
    T = ForwardDiff.valtype(eltype(C))
    I = one(C)
    J = sqrt(det(C))
    invC = inv(C)
    return  convert(T, μ) * det(C)^(convert(T, -1/3)) * (I - tr(C) * invC / 3) + convert(T, Kb) * (J - 1) * J * invC
end

const μ = 1e10;
const Kb = 1.66e11;

Ψ(C) = Ψ(C, μ, Kb)
S(C) = S(C, μ, Kb)

SUITE["gradient"] = BenchmarkGroup(["ad"])
SUITE["hessian"] = BenchmarkGroup(["ad"])

for dim in (ALL_DIMENSIONS ? (1,2,3) : (3,))
    for T in (Float64, Float32)
        V2 = tensor_dict[(dim, 2, T)]
        V2 = V2' ⋅ V2
        V2sym = symtensor_dict[(dim, 2, T)]
        V2sym = V2sym' ⋅ V2sym

        @assert eltype(Ψ(V2)) == T
        @assert eltype(S(V2)) == T
        @assert eltype(∇(Ψ, V2)) == T
        @assert eltype(Δ(Ψ, V2)) == T

        SUITE["gradient"]["∇Ψ(SymmetricTensor{2, $dim, $T})"]       = @benchmarkable ∇(Ψ, $V2sym)
        SUITE["gradient"]["∇Ψ(Tensor{2, $dim, $T})"]                = @benchmarkable ∇(Ψ, $V2)
        SUITE["gradient"]["∇Ψ(SymmetricTensor{2, $dim, $T}) - ana"] = @benchmarkable S($V2sym)
        SUITE["gradient"]["∇Ψ(Tensor{2, $dim, $T}) - ana"]          = @benchmarkable S($V2)
        SUITE["gradient"]["∇S(SymmetricTensor{2, $dim, $T})"]       = @benchmarkable ∇(S, $V2sym)
        SUITE["gradient"]["∇S(Tensor{2, $dim, $T})"]                = @benchmarkable ∇(S, $V2)

        SUITE["hessian"]["ΔΨ(SymmetricTensor{2, $dim, $T})"] = @benchmarkable Δ(Ψ, $V2sym)
        SUITE["hessian"]["ΔΨ(Tensor{2, $dim, $T})"]          = @benchmarkable Δ(Ψ, $V2)
    end
end


text:
# test on partial products with SpecialEuclidean(2)

using Manifolds
using ApproxManifoldProducts
using Test
# using Random
using FileIO

##

@testset "partial product with a SpecialEuclidean(2)" begin
##

datafile = joinpath(@__DIR__, "testdata", "partialtest.jld2")
dict = load(datafile)
pts1 = dict["pts1"]
pts2 = dict["pts2"]

randU = Float64[]
randN = Float64[]

len = length(pts1)

# define test manifold
M = SpecialEuclidean(2)
e0 = identity_element(M)

# p1 full SpecialEuclidean(2)
p1 = manikde!(M, pts1)
p1_ = marginal(p1,[1;2])

# p2 only Translation(2) part
_p2 = manikde!(M, pts2)
p2 = marginal(_p2, [1;2])

# product of full and marginal
# p12 = p1*p2

## =======================FIRST PRODUCT================================================

selectedLabels=Vector{Vector{Int}}()
p12 = manifoldProduct([p1; p2]; addEntropy=false,
                                recordLabels=true,
                                selectedLabels=selectedLabels,
                                _randU=randU,
                                _randN=randN  )
#
selectedLabels
#

@test !isPartial(p12)

p12_ = marginal(p12, [1;2])
@test isPartial(p12_)

_p12_ = manikde!(TranslationGroup(2), getPoints(p12_))

## intermediate test, check product of selected kernels match what is in the marginal

for sidx = 1:len
  bw1 = getBW(p1)[:,1] .^2
  bw2 = getBW(p2, false)[:,1] .^2

  u1 = pts1[selectedLabels[sidx][1]]
  u2 = pts2[selectedLabels[sidx][2]]

  u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);

  @test isapprox( u12.parts[1], getPoints(p12)[sidx].parts[1] )
end

## now check the marginal dimensions only

# now do submanifold dimensions separately as reference test -- should get a similar result
pts1_ = getPoints(p1_)
pts2_ = getPoints(p2)

## Do the translation part separate

for sidx = 1:len
  bw1 = getBW(p1_)[:,1] .^2
  bw2 = getBW(p2)[:,1] .^2

  u1 = pts1_[selectedLabels[sidx][1]]
  u2 = pts2_[selectedLabels[sidx][2]]

  u12, = calcProductGaussians(TranslationGroup(2), [u1,u2], [bw1,bw2])

  @test isapprox( u12, getPoints(p12)[sidx].parts[1] )
end

## =======================SECOND PRODUCT================================================
## NEW PRODUCT ON ONLY THE TRANSLATION PART, COMPARE WITH ABOVE RESULT

p1__ = manikde!(TranslationGroup(2), pts1_)
p2__ = manikde!(TranslationGroup(2), pts2_)

# p12__ = p1__*p2__
selectedLabels__=Vector{Vector{Int}}()
p12__ = manifoldProduct([p1__; p2__]; addEntropy=false,
                                      recordLabels=true,
                                      selectedLabels=selectedLabels__,
                                      _randU=randU,
                                      _randN=randN  )
#
selectedLabels__

## compare calcProduct of full*partial with selection from partial*partial

sidx = 1
for sidx = 1:len

bw1 = getBW(p1)[:,1] .^2
bw2 = getBW(p2, false)[:,1] .^2

# full-partial points, but selected from partial-partial product
u1 = pts1[selectedLabels__[sidx][1]]
u2 = pts2[selectedLabels__[sidx][2]]

u12, = calcProductGaussians(M, [u1,u2], [bw1,bw2]);
u12_, = calcProductGaussians(TranslationGroup(2), [u1.parts[1],u2.parts[1]], [bw1[1:2],bw2[1:2]]);

@test isapprox( u12.parts[1], u12_ )
@test isapprox( getPoints(p12__)[sidx], u12_ )

end

##

@test mmd(_p12_, p12__) < 0.001

##
end

## DEBUG PLOTTING ==========================================================================

## Plots showing the problem, p12 was wrong!!

# using Cairo, RoMEPlotting
# Gadfly.set_default_plot_size(35cm,20cm)

# n=10; plotKDE([p1_;p2; p12], levels=3, selectedPoints=selectedLabels[n:n])
# n=3; plotKDE([p1__;p2__; p12__], levels=3, selectedPoints=selectedLabels__[n:n])
# plotKDE([p1__; p2__; p12__])

# plotKDE([_p12_; p12__])

##


text:
## convolution

function conv2d!{T}(y::CuArray{T}, x::CuArray{T}, w::CuArray{T};
                    handle=cudnnhandle(), algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0,
                    alpha=1, beta=0, o...) # padding=0, stride=1, upscale=1, mode=0
    @cuda(cudnn, cudnnConvolutionForward,
          (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
          handle,Ref(T(alpha)),TD(x),x.ptr,FD(w),w.ptr,CD(w,x;o...),algo,workSpace,
          workSpaceSizeInBytes,Ref(T(beta)),TD(y),y.ptr)
    return y
end

function conv2d{T}(x::CuArray{T}, w::CuArray{T};
                  handle=cudnnhandle(), algo=0, workSpace=C_NULL, workSpaceSizeInBytes=0, alpha=1,
                  o...) # padding=0, stride=1, upscale=1, mode=0
    y = similar(x, cdims(w,x;o...))
    conv2d!(y, x, w; handle=handle, algo=algo, workSpace=workSpace,
            workSpaceSizeInBytes=workSpaceSizeInBytes, alpha=1, o...)
    return y
end

function conv2_grad_x!{T}(dx::CuArray{T}, x::CuArray{T}, w::CuArray{T}, dy::CuArray{T};
                          handle=cudnnhandle(), alpha=1, algo=0, beta=0, workSpace=C_NULL,
                          workSpaceSizeInBytes=0, o...) # padding=0, stride=1, upscale=1, mode=0
    cudnnVersion = cudnn_version()
    if cudnnVersion >= 4000
        @cuda(cudnn,cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w.ptr,TD(dy),dy.ptr,CD(w,x;o...),algo,workSpace,
              workSpaceSizeInBytes,Ref(T(beta)),TD(dx),dx.ptr)
    elseif cudnnVersion >= 3000
        @cuda(cudnn,cudnnConvolutionBackwardData_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w.ptr,TD(dy),dy.ptr,CD(w,x;o...),algo,workSpace,
              workSpaceSizeInBytes,Ref(T(beta)),TD(dx),dx.ptr)
    else
        @cuda(cudnn,cudnnConvolutionBackwardData,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),FD(w),w.ptr,TD(dy),dy.ptr,CD(w,x;o...),
              Ref(T(beta)),TD(dx),dx.ptr)
    end
    return dx
end

function conv2d_grad_x{T}(x::CuArray{T}, w::CuArray{T}, dy::CuArray{T};
                          handle=cudnnhandle(), alpha=1, algo=0, workSpace=C_NULL,
                          workSpaceSizeInBytes=0, o...) # padding=0, stride=1, upscale=1, mode=0
    dx = similar(x)
    conv2_grad_x!(dx, x, w, dy; handle=handle, alpha=alpha, algo=algo, workSpace=workSpace,
                   workSpaceSizeInBytes=workSpaceSizeInBytes, o...)
    return dx
end

function conv2d_grad_w!{T}(dw::CuArray{T}, x::CuArray{T}, w::CuArray{T}, dy::CuArray{T};
                           handle=cudnnhandle(), alpha=1, beta=0, algo=0, workSpace=C_NULL,
                           workSpaceSizeInBytes=0,
                           o...) # padding=0, stride=1, upscale=1, mode=0
    cudnnVersion = cudnn_version()
    if cudnnVersion >= 4000
        @cuda(cudnn,cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x.ptr,TD(dy),dy.ptr,CD(w,x;o...),algo,workSpace,
              workSpaceSizeInBytes,Ref(T(beta)),FD(dw),dw.ptr)
    elseif cudnnVersion >= 3000
        @cuda(cudnn,cudnnConvolutionBackwardFilter_v3,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,UInt32,Cptr,Csize_t,Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x.ptr,TD(dy),dy.ptr,CD(w,x;o...),algo,workSpace,
              workSpaceSizeInBytes,Ref(T(beta)),FD(dw),dw.ptr)
    else
        @cuda(cudnn,cudnnConvolutionBackwardFilter,
              (Cptr,Ptr{T},Cptr,Ptr{T},Cptr,Ptr{T},Cptr,       Ptr{T},Cptr,Ptr{T}),
              handle,Ref(T(alpha)),TD(x),x.ptr,TD(dy),dy.ptr,CD(w,x;o...),Ref(T(beta)),FD(dw),dw.ptr)
    end
    return dw
end

function conv2d_grad_w{T}(x::CuArray{T}, w::CuArray{T}, dy::CuArray{T};
                          handle=cudnnhandle(), alpha=1, algo=0, workSpace=C_NULL,
                          workSpaceSizeInBytes=0,
                   o...) # padding=0, stride=1, upscale=1, mode=0
    dw = similar(w)
    conv2d_grad_w!(dw, x, w, dy; handle=handle, alpha=alpha, algo=algo,  workSpace=workSpace,
                   workSpaceSizeInBytes=workSpaceSizeInBytes, o...)
    return dw
end

## deconvolution

function deconv2d(x, w; o...)
    y = similar(x, dcdims(w, x; o...))
    return conv2_grad_x(y, w, x; o...)
end

function deconv2_grad_w(x, w, dy; o...)
    return conv2_grad_w(dy, w, x; o...)
end

function deconv2_grad_x(x, w, dy; o...)
    return conv2d(dy, w; o...)
end


text:
#module DDForDP
using Statistics
using UUIDs

include("./utilities/utilitiesMain.jl")
include("./intervalTrees/treesMain.jl")
include("./model/modelMain.jl")

@enum RRBSetting sop
include("./solvers/solversMain.jl")
#end #module

include("./io/ioMain.jl")
include("./examples/examples.jl")

function testRelaxedBound()
    model = testModelMain()

    (bound, path) = calculateRelaxedBound(model, 2)

    println(bound)
    for (key,value) in path
        println("   $key => $value")
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PicoSAT"
version = v"965.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://fmv.jku.at/picosat/picosat-965.tar.gz", "15169b4f28ba8f628f353f6f75a100845cdef4a2244f101a02b6e5a26e46a754"),
    DirectorySource("./bundled")
    ]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/picosat*
cp ../{Makefile,config.h} .
if [[ ${target} == *musl* ]]; then
    sed -i 's!sys/unistd.h!unistd.h!g' picosat.c
fi
make -j${nproc}
make install
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("picosat", :picosat),
    ExecutableProduct("picomus", :picomus),
    ExecutableProduct("picomcs", :picomcs),
    ExecutableProduct("picogcnf", :picogcnf),
    LibraryProduct("libpicosat",:libpicosat)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module TestFilters

using ..Fixtures
using Log4jl
import Log4jl: FILTER, MARKER, filter, isfiltered, MarkerFilter, ThresholdFilter, FilterResult, state, start, stop
using FactCheck

evnt = Fixtures.TestEvent() # Level = INFO

facts("Filters") do
    context("need to implement certain methods") do
        @fact filter(Fixtures.TESTFILTER, evnt) --> FilterResult.NEUTRAL
        @fact filter(Fixtures.TESTFILTER, Fixtures.TESTLEVEL, MARKER(Fixtures.TESTMARKER), Fixtures.TESTMSG) --> FilterResult.NEUTRAL
        @fact isfiltered(FILTER(Fixtures.TESTFILTER), evnt) --> false
    end
    context("can filter on event marker") do
        flt = MarkerFilter(Fixtures.TESTMARKER)
        @fact isfiltered(FILTER(flt), evnt) --> false
        @fact filter(flt, evnt) --> FilterResult.NEUTRAL
        flt = MarkerFilter(Fixtures.TESTMARKER, FilterResult.DENY, FilterResult.ACCEPT)
        @fact isfiltered(FILTER(flt), evnt) --> true
        @fact filter(flt, evnt) --> FilterResult.DENY
        @fact filter(flt, Fixtures.TESTLEVEL, MARKER(Fixtures.TESTMARKER), Fixtures.TESTMSG) --> FilterResult.DENY
        @fact state(flt) --> Log4jl.LifeCycle.INITIALIZED
        start(FILTER(flt))
        @fact state(flt) --> Log4jl.LifeCycle.STARTED
        stop(FILTER(flt))
        @fact state(flt) --> Log4jl.LifeCycle.STOPPED
    end

    context("can filter on event level") do
        flt = ThresholdFilter() # Level = ERROR
        @fact isfiltered(FILTER(flt), evnt) --> true
        @fact filter(flt, evnt) --> FilterResult.DENY
        flt = ThresholdFilter(Log4jl.Level.DEBUG)
        @fact isfiltered(FILTER(flt), evnt) --> false
        @fact filter(flt, evnt) --> FilterResult.NEUTRAL
        flt = ThresholdFilter(Log4jl.Level.WARN, FilterResult.ACCEPT, FilterResult.NEUTRAL)
        @fact isfiltered(FILTER(flt), evnt) --> false
        @fact filter(flt, evnt) --> FilterResult.NEUTRAL
        @fact filter(flt, Fixtures.TESTLEVEL, MARKER(Fixtures.TESTMARKER), Fixtures.TESTMSG) --> FilterResult.NEUTRAL
        @fact state(flt) --> Log4jl.LifeCycle.INITIALIZED
        start(flt)
        @fact state(flt) --> Log4jl.LifeCycle.STARTED
        stop(flt)
        @fact state(flt) --> Log4jl.LifeCycle.STOPPED
    end
end

end