text:
#=
Write an algorithm to justify text. Given a sequence of words and an integer line length k, return a list of strings which represents each line, fully justified.

More specifically, you should have as many words as possible in each line. There should be at least one space between each word. Pad extra spaces when necessary so that each line has exactly length k. Spaces should be distributed as equally as possible, with the extra spaces, if any, distributed starting from the left.

If you can only fit one word on a line, then you should pad the right-hand side with spaces.

Each word is guaranteed not to be longer than k.

For example, given the list of words ["the", "quick", "brown", "fox", "jumps", "over", "the", "lazy", "dog"] and k = 16, you should return the following:

["the  quick brown", # 1 extra space on the left
"fox  jumps  over", # 2 extra spaces distributed evenly
"the   lazy   dog"] # 4 extra spaces distributed evenly
=#

function justify(text::Array{String,1}, k)
    justified_text = [""]
    for word in text
        justified_text[end], istoolong = addword(justified_text[end], word, k)
        if istoolong
            justified_text[end] = justify(justified_text[end], k)
            push!(justified_text, word)
        end
    end
    justified_text[end] = justify(justified_text[end], k)
    return justified_text
end

function justify(sentence::String, k)
    spaces_to_add = k - length(sentence)
    spaces_in_sentence = findall(isequal(' '), sentence)
    space_pos = spaces_in_sentence[1] - 2
    justified_sentence = sentence
    char_to_skip = 2
    if spaces_in_sentence == 0
        justified_sentence = sentence * repeat(" ", spaces_to_add)
    else
        while spaces_to_add != 0
            space_pos = findnext(isequal(' '), justified_sentence, space_pos + char_to_skip)
            if isnothing(space_pos)
                space_pos = findfirst(isequal(' '), justified_sentence)
                char_to_skip = char_to_skip + 1
            end
            justified_sentence = justified_sentence[1:space_pos] * " " * justified_sentence[space_pos + 1:end]
            spaces_to_add = spaces_to_add - 1
        end
    end
    return justified_sentence
end

function addword(sentence, word, k)
    if sentence == ""
        new_sentence = word
        istoolong = false
    elseif length(sentence) + 1 + length(word) <= k
        new_sentence = sentence * " " * word
        istoolong = false
    else
        new_sentence = sentence
        istoolong = true
    end
    return new_sentence, istoolong
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using PowerSpectra
using Documenter

ENV["PLOTS_DEFAULT_BACKEND"] = "GR"
ENV["GKSwstype"]="nul"
using Plots
using Plots.PlotMeasures: mm

default(
    fontfamily = "Computer Modern", linewidth=1.5,
    titlefontsize=(16+8), guidefontsize=(11+5),
    tickfontsize=(8+4), legendfontsize=(8+4),
    left_margin=5mm, right_margin=5mm)

makedocs(;
    modules=[PowerSpectra],
    authors="Zack Li",
    repo="https://github.com/xzackli/PowerSpectra.jl/blob/{commit}{path}#L{line}",
    sitename="PowerSpectra.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xzackli.github.io/PowerSpectra.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Quickstart" => "quickstart.md",
        "Spectra" => "spectra.md",
        "Covariance" => "covariance.md",
        "Beams" => "beams.md",
        "Utilities" => "util.md",
        "Index" => "module_index.md",
    ],
)

deploydocs(;
    repo="github.com/xzackli/PowerSpectra.jl",
)

ext: jl

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
using PyPlot, Seismic

download("http://seismic.physics.ualberta.ca/data/marmvel.bin","marmvel.bin");
n1 = 751; n2 = 2301; dz =4.0; dx = 4.0;
v = reshape(read(open("marmvel.bin", "r"), Float32, n1*n2), n1, n2)
v = convert(Array{Float64,2}, 1000.*v);
vp = v[400:end, 1:500]
(nz, nx) = size(vp)
# receivers location
ix = collect(1:5:nx); iz = ones(Int64, length(ix)); pos = hcat(iz, ix);
# set up SparseMatrixCSC for computing derivative in spatial
f0 = 30.0; dt = 3e-4;
fidMtx = AcousticSetup(vp, dz, dx, dt, f0, iflag=2);
# source location
isx = 251; isz = 1;
shot = SeisAcousticWave(fidMtx, pos, isz, isx, f0, dt, tmax=1.0, print_flag=true, interval=500)
SeisPlot(shot.d, name="ShotGather.pdf")
SeisPlot(vp, name = "vp.pdf")

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
import DelayDiffEq: MethodOfSteps, Tsit5

@testset "test_delay_system" begin

# For simplicity, equality of DelayLtiSystems are tested over a finite set of frequencies
ω = 0.0:8

@test typeof(promote(delay(0.2), ss(1))[1]) == DelayLtiSystem{Float64,Float64}

@test typeof(promote(delay(0.2), ss(1.0 + im))[1]) == DelayLtiSystem{Complex{Float64}, Float64}

@test sprint(show, ss(1,1,1,1)*delay(1.0)) == "DelayLtiSystem{Float64, Float64}\n\nP: StateSpace{Continuous, Float64}\nA = \n 1.0\nB = \n 0.0  1.0\nC = \n 1.0\n 0.0\nD = \n 0.0  1.0\n 1.0  0.0\n\nContinuous-time state-space model\n\nDelays: [1.0]"

# Extremely basic tests
@test freqresp(delay(1), ω) ≈ reshape(exp.(-im*ω), length(ω), 1, 1) rtol=1e-15
@test freqresp(delay(2.5), ω)[:] ≈ exp.(-2.5im*ω) rtol=1e-15
@test freqresp(3.5*delay(2.5), ω)[:] ≈ 3.5*exp.(-2.5im*ω) rtol=1e-15
@test freqresp(delay(2.5)*1.5, ω)[:] ≈ exp.(-2.5im*ω)*1.5 rtol=1e-15

# Addition of constant
@test evalfr(1 + delay(1.0), 0)[] ≈ 2
@test evalfr(1 - delay(1.0), 0)[] ≈ 0
@test evalfr([2 -delay(1.0)], 0) ≈ [2 -1]

# Stritcly proper system
P1 = DelayLtiSystem(ss(-1.0, 1, 1, 0))
P1_fr = 1 ./ (im*ω .+ 1)
@test freqresp(P1, ω)[:] == P1_fr

# Not stritcly proper system
P2 = DelayLtiSystem(ss(-2.0, -1, 1, 1)) # (s+1)/(s+2)
P2_fr = (im*ω .+ 1) ./ (im*ω .+ 2)
@test freqresp(P2, ω)[:] ≈ P2_fr rtol=1e-15

## Addition
@test freqresp(1 + delay(1), ω)[:] ≈ 1 .+ exp.(-im*ω)
@test freqresp(P1 + delay(1), ω)[:] ≈ P1_fr .+ exp.(-im*ω)

# Substraction
@test freqresp(P1 - delay(1), ω)[:] ≈ P1_fr .- exp.(-im*ω)

## Multiplication by scalar
@test freqresp(2*delay(1), ω)[:] ≈ 2*exp.(-im*ω)
@test freqresp(2*(P1*delay(1)), ω)[:] ≈ 2*P1_fr.*exp.(-im*ω)

## Multiplication
@test freqresp(P1 * delay(1), ω)[:] ≈ P1_fr .* exp.(-im*ω) rtol=1e-15
@test freqresp(delay(1) * P1, ω)[:] ≈ P1_fr .* exp.(-im*ω) rtol=1e-15

@test freqresp(P2 * delay(1), ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15
@test freqresp(delay(1) * P2, ω)[:] ≈ P2_fr .* exp.(-im*ω) rtol=1e-15

# Equality
@test P1 == deepcopy(P1)
@test P1 != deepcopy(P2)

# evalfr
s_vec = [0, 1im, 1, 1 + 1im]
@test [evalfr(delay(2), s)[1] for s in s_vec] ≈ [exp(-2*s) for s in s_vec] rtol=1e-16

## Feedback
# The first tests don't include delays, but the linear system is of DelayLtiForm type
# (very simple system so easy to troubleshoot)
@test freqresp(feedback(1.0, P1), ω)[:] ≈ 1 ./ (1 .+  P1_fr) rtol=1e-15
@test freqresp(feedback(P1, 1.0), ω)[:] ≈ P1_fr ./ (1 .+  P1_fr) rtol=1e-15
@test freqresp(feedback(P1, P1), ω)[:] ≈ P1_fr ./ (1 .+ P1_fr .* P1_fr) rtol=1e-15

@test freqresp(feedback(1.0, DelayLtiSystem(ss(0.5))), [0])[:] == [2/3]
@test freqresp(feedback(1.0, P2), ω)[:] ≈ 1 ./ (1 .+ P2_fr)

@test freqresp(feedback(0.5, delay(2.0)), ω) ≈ 0.5 ./ (1 .+ 0.5*exp.(-2im*ω))
@test freqresp(feedback(delay(2.0), 0.5), ω) ≈ exp.(-2im*ω) ./ (1 .+ 0.5*exp.(-2im*ω))

@test freqresp(feedback(P1, delay(1)), ω)[:] ≈ P1_fr ./ (1 .+ exp.(-im*ω) .* P1_fr) rtol=1e-15
@test freqresp(feedback(delay(1), P1), ω)[:] ≈ exp.(-im*ω) ./ (1 .+ exp.(-im*ω) .* P1_fr) rtol=1e-15 #FIXME: Answer is Inf, but should give error.. rtol=1e-15
@test freqresp(feedback(P1*delay(1), 1.0), ω)[:] ≈ P1_fr .* exp.(-im*ω) ./ (1 .+ exp.(-im*ω) .* P1_fr) rtol=1e-15
@test freqresp(feedback(1.0, P1*delay(1)), ω)[:] ≈ 1 ./ (1 .+ exp.(-im*ω) .* P1_fr) rtol=1e-15

@test freqresp(feedback(1.0, P2*0.5*(ss(1.0) + delay(2))), ω)[:] ≈ 1 ./(1 .+ P2_fr .* 0.5.*(1 .+ exp.(-2*im*ω)))

@test freqresp(1.0 + delay(2), ω)[:] ≈ 1 .+ exp.(-2im*ω) rtol=1e-15

G = ss(0.5) + 0.5*delay(2)# + 0.5*delay(3)
G_fr = 0.5 .+ 0.5*exp.(-2*im*ω)# + 0.5*exp.(-3*im*ω)
@test freqresp(G, ω)[:] ≈ G_fr rtol=1e-15

@test freqresp(feedback(1.0, P1*G), ω)[:] ≈ 1 ./(1 .+ P1_fr .* G_fr) rtol=1e-15
@test freqresp(feedback(P1, G), ω)[:] ≈ P1_fr ./(1 .+ P1_fr .* G_fr) rtol=1e-15

@test freqresp(feedback(P2, G), ω)[:] ≈ P2_fr ./(1 .+ P2_fr .* G_fr) rtol=1e-15

s = tf("s")

# Test alternative exp constructor for delays
d = exp(-2*s)
@test freqresp(d, [0, 1, 2]) ≈ [1, exp(-2im), exp(-4im)]

@test_throws ErrorException exp(-s^2 - 2*s)
@test_throws ErrorException exp(-2*s+1) # in principle ok, but not allowed anyway
@test_throws ErrorException exp([-2*s; -s])
@test_throws ErrorException exp(2*s) # Non-causal

Ω = [0.5, 1.0, 1.5, 2.0]
# Test for internal function delayd_ss
@test freqresp(ControlSystems.delayd_ss(1.0, 0.2), Ω)[:] ≈ exp.(-im*Ω) atol=1e-14
@test freqresp(ControlSystems.delayd_ss(3.2, 0.4), Ω)[:] ≈ exp.(-3.2*im*Ω) atol=1e-14
@test_throws ErrorException freqresp(ControlSystems.delayd_ss(3.2, 0.5), Ω)

# Simple tests for c2d of DelayLtiSystems
@test freqresp(c2d(feedback(ss(0,1,1,0), delay(1.5)), 0.5), Ω) ≈ [0.5/((z - 1) + 0.5*z^-3) for z in exp.(im*Ω*0.5)]
@test freqresp(c2d(feedback(delay(1.5), delay(1.0)), 0.5), Ω) ≈ [z^-3/(1 + z^-5) for z in exp.(im*Ω*0.5)]
@test freqresp(c2d(feedback(0.5, delay(1.5)), 0.5), Ω) ≈ [0.5/(1 + 0.5*z^-3) for z in exp.(im*Ω*0.5)]

# Random conversions
sys1 = DelayLtiSystem(1.0/s)
@test sys1.P.A == sys1.P.D == fill(0,1,1)
@test sys1.P.B*sys1.P.C == fill(1,1,1)
@test sys1.Tau == []

## Multiple Delays
G = ss(1.0) + delay(2) + delay(3)
G_fr = 1 .+ exp.(-2*im*ω) .+ exp.(-3*im*ω)
@test freqresp(G, ω)[:] ≈ G_fr rtol=1e-15
@test freqresp(feedback(1.0, G), ω)[:] ≈ 1 ./(1 .+ G_fr) # Somewhat pathological system though

sys = feedback(1.0, P1 * (delay(2) + delay(3)))
expected_sys_fr = 1.0 ./ (P1_fr .* (exp.(-2*im*ω) + exp.(-3*im*ω)) .+ 1)

@test freqresp(sys, ω)[:] ≈ expected_sys_fr rtol=1e-14

@test freqresp(feedback(1.0*P1, G*G*P2*P2), ω)[:] ≈ P1_fr ./(1 .+ (G_fr.^2).*P1_fr.*P2_fr.^2) rtol=1e-15# Somewhat pathological system though

s11 = feedback(ss(1/s), delay(1))
s12 = ss(1/s)
s21 = DelayLtiSystem(ss(1/(s+1)))
s22 = ss(1/(s+10))

s1 = [s11 s12]
s2 = [s21 s22]

f1 = [s1;s2]
f2 = [s11 s12;
      s21 s22]

# Test that different consatenations work
w = 10 .^ (-2:0.1:2)
@test freqresp(f1, w) ≈ freqresp(f2, w) rtol=1e-15

# This used to be a weird bug
@test freqresp(s11, w) ≈ freqresp(f2[1,1], w) rtol=1e-15

@test propertynames(delay(1.0)) == (:P, :Tau)

#FIXME: A lot more tests, including MIMO systems in particular

# Test step
println("Simulating first delay system:")
@time step(delay(1)*tf(1,[1.,1]))
@time step(delay(1)*tf(1,[1,1]))

@time y1, t1, x1 = step([s11;s12], 10)
@time @test y1[2:2,:] ≈ step(s12, t1)[1] rtol = 1e-14

t = 0.0:0.1:10
y2, t2, x2 = step(s1, t)
# TODO Figure out which is inexact here
@test y2[1:1,:,1:1] + y2[1:1,:,2:2] ≈ step(s11, t)[1] + step(s12, t)[1] rtol=1e-14

y3, t3, x3 = step([s11; s12], t)
@test y3[1:1,:,1] ≈ step(s11, t)[1] rtol = 1e-14
@test y3[2:2,:,1] ≈ step(s12, t)[1] rtol = 1e-14

y1, t1, x1 = step(DelayLtiSystem([1.0/s 2/s; 3/s 4/s]), t)
y2, t2, x2 = step([1.0/s 2/s; 3/s 4/s], t)
@test y1 ≈ y2 rtol=1e-14
@test size(x1,2) == length(t)
@test size(x1,3) == 2

##################### Test with known solution
K = 2.0
sys_known = feedback(delay(1)*K/s, 1)

ystep, t, _ = step(sys_known, 3)

function y_expected(t, K)
      if t < 1
            return 0
      elseif t < 2
            return K*(t-1)
      elseif t <= 3
            return K*(t-1)-1/2*K^2*(t-2)^2
      else
            throw(ArgumentError("Test not defined here"))
      end
end

@test ystep' ≈ y_expected.(t, K) atol = 1e-12

function dy_expected(t, K)
      if t < 1
            return 0
      elseif t < 2
            return K
      elseif t <= 3
            return K - K^2*(t-2)
      else
            throw(ArgumentError("Test not defined here"))
      end
end

y_impulse, t, _ = impulse(sys_known, 3)

# TODO Better accuracy for impulse
@test y_impulse' ≈ dy_expected.(t, K) rtol=1e-13
@test maximum(abs, y_impulse' - dy_expected.(t, K)) < 1e-12

y_impulse, t, _ = impulse(sys_known, 3, alg=MethodOfSteps(Tsit5()))
# Two orders of magnitude better with BS3 in this case, which is default for impulse
@test y_impulse' ≈ dy_expected.(t, K) rtol=1e-5
@test maximum(abs, y_impulse' - dy_expected.(t, K)) < 1e-5

## Test delay with D22 term
t = 0:0.01:4

sys = delay(1)

y, t, x = step(sys, t)
@test y[:] ≈ [zeros(100); ones(301)] atol = 1e-12
@test size(x) == (0,401)

sys = delay(1)*delay(1)*1/s

y, t, x = step(sys, t)

y_sol = [zeros(200);0:0.01:2]'

@test maximum(abs,y-y_sol) < 1e-13
@test maximum(abs,x-collect(0:0.01:4)') < 1e-15

# TODO For some reason really bad accuracy here
# Looks like a lag in time
sys = 1/s*delay(1)*delay(1)

y, t, x = step(sys, t)
@test maximum(abs,y-y_sol) < 1e-5
@test maximum(abs,x-y_sol) < 1e-5

t = 0:0.001:0.1
y, t, x = step(sys, t)
@test length(y) == length(t)

##  Test of basic pade functionality

Ω = [0, 0.5, 1, 2, 5]
@test freqresp(pade(1, 1), Ω) == freqresp(tf([-1/2, 1], [1/2, 1]), Ω)
@test freqresp(pade(1, 2), Ω) ≈ freqresp(tf([1/12, -1/2, 1], [1/12, 1/2, 1]), Ω)

for (n, tol)=enumerate([0.05; 1e-3; 1e-5; 1e-7; 1e-11; 1e-14*ones(5)])
    G = pade(0.8, n)

    @test isstable(G)
    @test evalfr(G, 0)[1] ≈ 1
    @test abs(evalfr(G, 2im)[1]) ≈ 1
    @test evalfr(G, 1im)[1] ≈ exp(-0.8im) atol=tol
end

## Test pade applied to DelayLtiSystem

@test freqresp(pade(delay(0.5), 2), Ω) ≈ freqresp(pade(0.5, 2), Ω)

P = delay(1.0) * DemoSystems.lag(T=1)
@test freqresp(pade(feedback(1,P), 2), Ω) == freqresp(feedback(1, pade(P,2)), Ω)

Ω = [0, 0.1, 0.2]
P_wb = DemoSystems.woodberry()

@test freqresp(pade(P_wb, 2), Ω) ≈ freqresp(P_wb, Ω) atol=0.02
@test freqresp(pade(P_wb, 3), Ω) ≈ freqresp(P_wb, Ω) atol=5e-4

@test freqresp(pade(feedback(eye_(2), P_wb), 3), Ω) ≈ freqresp(feedback(eye_(2), P_wb), Ω) atol=1e-4

end


text:
# Put additional units here
@unit lbm        "lbm"       PoundMass   1Unitful.lb            false
@unit mmHg       "mmHg"      MMHg        (1//760)Unitful.atm      false
@unit gal        "gal"       USGallon    231Unitful.inch^3    false

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


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

mutable struct RunCommandDocument <: SwaggerModel
    schema::Any # spec type: Union{ Nothing, String } # spec name: \$schema
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    label::Any # spec type: Union{ Nothing, String } # spec name: label
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    script::Any # spec type: Union{ Nothing, Vector{String} } # spec name: script
    parameters::Any # spec type: Union{ Nothing, Vector{RunCommandParameterDefinition} } # spec name: parameters

    function RunCommandDocument(;schema=nothing, id=nothing, osType=nothing, label=nothing, description=nothing, script=nothing, parameters=nothing)
        o = new()
        validate_property(RunCommandDocument, Symbol("\$schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(RunCommandDocument, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(RunCommandDocument, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(RunCommandDocument, Symbol("label"), label)
        setfield!(o, Symbol("label"), label)
        validate_property(RunCommandDocument, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(RunCommandDocument, Symbol("script"), script)
        setfield!(o, Symbol("script"), script)
        validate_property(RunCommandDocument, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        o
    end
end # type RunCommandDocument

const _property_map_RunCommandDocument = Dict{Symbol,Symbol}(Symbol("\$schema")=>Symbol("schema"), Symbol("id")=>Symbol("id"), Symbol("osType")=>Symbol("osType"), Symbol("label")=>Symbol("label"), Symbol("description")=>Symbol("description"), Symbol("script")=>Symbol("script"), Symbol("parameters")=>Symbol("parameters"))
const _property_types_RunCommandDocument = Dict{Symbol,String}(Symbol("\$schema")=>"String", Symbol("id")=>"String", Symbol("osType")=>"String", Symbol("label")=>"String", Symbol("description")=>"String", Symbol("script")=>"Vector{String}", Symbol("parameters")=>"Vector{RunCommandParameterDefinition}")
Base.propertynames(::Type{ RunCommandDocument }) = collect(keys(_property_map_RunCommandDocument))
Swagger.property_type(::Type{ RunCommandDocument }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_RunCommandDocument[name]))}
Swagger.field_name(::Type{ RunCommandDocument }, property_name::Symbol) =  _property_map_RunCommandDocument[property_name]

const _allowed_RunCommandDocument_osType = ["Windows", "Linux"]

function check_required(o::RunCommandDocument)
    (getproperty(o, Symbol("\$schema")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("label")) === nothing) && (return false)
    (getproperty(o, Symbol("description")) === nothing) && (return false)
    (getproperty(o, Symbol("script")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RunCommandDocument }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "RunCommandDocument", :enum, val, _allowed_RunCommandDocument_osType)
    end
end


text:
using Plots
gr()

pl = plot()

m = BoyerLindquist(M=1.0, a=0.998)
CarterBoyerLindquist.rms(m.M, m.a)
# m = JohannsenAD(M=1.0, a=a)
r_range = R_isco:3.0:30
vϕs = @time find_vϕ_for_orbit_range(m, r_range; upper=0.1)
geodesics = map(i -> geodesic_for(m, r_range[i], vϕs[i]), eachindex(r_range))

for sol in geodesics
    plot!(
        sol,
        vars=(8, 6),
        projection=:polar,
        range=(0.0, last(r_range)),
        legend=false
    )
end

pl


text:
module PackageCompilerX

# TODO: Add good debugging statements
# TODO: sysimage or sysimg...

using Base: active_project
using Libdl: Libdl
using Pkg: Pkg
using UUIDs: UUID
using DocStringExtensions: SIGNATURES, TYPEDEF

export create_sysimage, create_app, audit_app, restore_default_sysimg

include("juliaconfig.jl")

# TODO: Check more carefully how to just use mingw on windows without using cygwin.
function get_compiler()
    if Sys.iswindows()
        return `x86_64-w64-mingw32-gcc`
    else
        if Sys.which("gcc") !== nothing
            return `gcc`
        elseif Sys.which("clang") !== nothing
            return `clang`
        end
        error("could not find a compiler, looked for `gcc` and `clang`")
    end
end

# TODO: Be able to set target for -C?
# TODO: Change to commented default ?
# const DEFAULT_TARGET = "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)"
const DEFAULT_TARGET = "generic"
current_process_sysimage_path() = unsafe_string(Base.JLOptions().image_file)

function get_julia_cmd()
    julia_path = joinpath(Sys.BINDIR, Base.julia_exename())
    cmd = `$julia_path --color=yes --startup-file=no --cpu-target=$DEFAULT_TARGET`
end

all_stdlibs() = readdir(Sys.STDLIB)

function rewrite_sysimg_jl_only_needed_stdlibs(stdlibs::Vector{String})
    sysimg_source_path = Base.find_source_file("sysimg.jl")
    sysimg_content = read(sysimg_source_path, String)
    # replaces the hardcoded list of stdlibs in sysimg.jl with
    # the stdlibs that is given as argument
    return replace(sysimg_content, r"stdlibs = \[(.*?)\]"s => string("stdlibs = [", join(":" .* string.(stdlibs), ",\n"), "]"))
end

function create_fresh_base_sysimage(stdlibs::Vector{String})
    tmp = mktempdir(cleanup=false)
    sysimg_source_path = Base.find_source_file("sysimg.jl")
    base_dir = dirname(sysimg_source_path)
    tmp_corecompiler_ji = joinpath(tmp, "corecompiler.ji")
    tmp_sys_ji = joinpath(tmp, "sys.ji")
    compiler_source_path = joinpath(base_dir, "compiler", "compiler.jl")

    @info "PackageCompilerX: creating base system image (incremental=false), this might take a while..."
    cd(base_dir) do
        # Create corecompiler.ji
        cmd = `$(get_julia_cmd()) --output-ji $tmp_corecompiler_ji -g0 -O0 $compiler_source_path`
        @debug "running $cmd"
        read(cmd)

        # Use that to create sys.ji
        new_sysimage_content = rewrite_sysimg_jl_only_needed_stdlibs(stdlibs)
        new_sysimage_source_path = joinpath(base_dir, "sysimage_packagecompiler_x.jl")
        write(new_sysimage_source_path, new_sysimage_content)
        try
            cmd = `$(get_julia_cmd()) --sysimage=$tmp_corecompiler_ji -g1 -O0 --output-ji=$tmp_sys_ji $new_sysimage_source_path`
            @debug "running $cmd"
            read(cmd)
        finally
            rm(new_sysimage_source_path; force=true)
        end
    end

    return tmp_sys_ji
end

# TODO: Add output file?
function run_precompilation_script(project::String, precompile_file::Union{String, Nothing})
    # TODO: Audit tempname usage
    tracefile = tempname()
    if precompile_file == nothing
        arg = `-e ''`
    else
        arg = `$precompile_file`
    end
    touch(tracefile)
    cmd = `$(get_julia_cmd()) --sysimage=$(current_process_sysimage_path()) --project=$project
            --compile=all --trace-compile=$tracefile $arg`
    @debug "run_precompilation_script: running $cmd"
    run(cmd)
    return tracefile
end

function create_sysimg_object_file(object_file::String, packages::Vector{Symbol};
                            project::String,
                            base_sysimg::String,
                            precompile_execution_file::Union{Vector{String}, Nothing},
                            precompile_statements_file::Union{Vector{String}, Nothing})
    # include all packages into the sysimg
    julia_code = """
        if !isdefined(Base, :uv_eventloop)
            Base.reinit_stdio()
        end
        Base.__init__();
        """
    for package in packages
        julia_code *= "using $package\n"
    end

    # handle precompilation
    precompile_statements = ""
    @debug "running precompilation execution script..."
    tracefiles = String[]
    for file in (precompile_execution_file === nothing ? (nothing,) : precompile_execution_file)
        tracefile = run_precompilation_script(project, file)
        precompile_statements *= "append!(precompile_statements, readlines($(repr(tracefile))))\n"
    end
    if precompile_statements_file != nothing
        for file in precompile_statements_file
            precompile_statements *=
                "append!(precompile_statements, readlines($(repr(file))))\n"
        end
    end

    precompile_code = """
        # This @eval prevents symbols from being put into Main
        @eval Module() begin
            PrecompileStagingArea = Module()
            for (_pkgid, _mod) in Base.loaded_modules
                if !(_pkgid.name in ("Main", "Core", "Base"))
                    eval(PrecompileStagingArea, :(const \$(Symbol(_mod)) = \$_mod))
                end
            end
            precompile_statements = String[]
            $precompile_statements
            for statement in sort(precompile_statements)
                # println(statement)
                try
                    Base.include_string(PrecompileStagingArea, statement)
                catch
                    # See julia issue #28808
                    @error "failed to execute \$statement"
                end
            end
        end # module
        """
    julia_code *= precompile_code

    # finally, make julia output the resulting object file
    @debug "creating object file at $object_file"
    @info "PackageCompilerX: creating system image object file, this might take a while..."

    cmd = `$(get_julia_cmd()) --sysimage=$base_sysimg --project=$project --output-o=$(object_file) -e $julia_code`
    @debug "running $cmd"
    run(cmd)
end

default_sysimg_path() = joinpath(julia_private_libdir(), "sys." * Libdl.dlext)
default_sysimg_name() = basename(default_sysimg_path())
backup_default_sysimg_path() = default_sysimg_path() * ".backup"
backup_default_sysimg_name() = basename(backup_default_sysimg_path())

# TODO: Also check UUIDs for stdlibs, not only names
function gather_stdlibs_project(project::String)
    project_toml_path = abspath(Pkg.Types.projectfile_path(project; strict=true))
    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    @assert ctx.env.manifest !== nothing
    stdlibs = all_stdlibs()
    stdlibs_project = String[]
    for (uuid, pkg) in ctx.env.manifest
        if pkg.name in stdlibs
            push!(stdlibs_project, pkg.name)
        end
    end
    return stdlibs_project
end

"""
    $SIGNATURES
"""
function create_sysimage(packages::Union{Symbol, Vector{Symbol}};
                         sysimage_path::Union{String,Nothing}=nothing,
                         project::String=active_project(),
                         precompile_execution_file::Union{String, Vector{String}, Nothing}=nothing,
                         precompile_statements_file::Union{String, Vector{String}, Nothing}=nothing,
                         incremental::Bool=true,
                         filter_stdlibs=false,
                         replace_defaut::Bool=false)
    if sysimage_path === nothing
        if replace_defaut == false
            error("`sysimage_path` cannot be `nothing` if `replace_defaut` is `false`")
        end
        # We will replace the default sysimage so just put it somewhere for now
        tmp = mktempdir()
        sysimage_path = joinpath(tmp, string("sys.", Libdl.dlext))
    end

    if filter_stdlibs && incremental
        error("must use `incremental=false` to use `filter_stdlibs=true`")
    end

    # Functions lower down handles precompilation file as arrays so convert here
    packages = vcat(packages)
    precompile_execution_file !== nothing && (precompile_execution_file = vcat(precompile_execution_file))
    precompile_statements_file !== nothing && (precompile_statements = vcat(precompile_statements_file))

    if !incremental
        if filter_stdlibs
            stdlibs = gather_stdlibs_project(project)
        else
            stdlibs= all_stdlibs()
        end
        base_sysimg = create_fresh_base_sysimage(stdlibs)
    else
        base_sysimg = current_process_sysimage_path()
    end

    object_file = tempname() * ".o"
    create_sysimg_object_file(object_file, packages;
                              project=project,
                              base_sysimg=base_sysimg,
                              precompile_execution_file=precompile_execution_file,
                              precompile_statements_file=precompile_statements_file)
    create_sysimg_from_object_file(object_file, sysimage_path)
    if replace_defaut
        if !isfile(backup_default_sysimg_path())
            @debug "making a backup of default sysimg"
            cp(default_sysimg_path(), backup_default_sysimg_path())
        end
        @info "PackageCompilerX: default sysimg replaced, restart Julia for the new sysimg to be in effect"
        mv(sysimage_path, default_sysimg_path(); force=true)
    end
    # TODO: Remove object file
end

function create_sysimg_from_object_file(input_object::String, sysimage_path::String)
    julia_libdir = dirname(Libdl.dlpath("libjulia"))

    # Prevent compiler from stripping all symbols from the shared lib.
    # TODO: On clang on windows this is called something else
    if Sys.isapple()
        o_file = `-Wl,-all_load $input_object`
    else
        o_file = `-Wl,--whole-archive $input_object -Wl,--no-whole-archive`
    end
    extra = Sys.iswindows() ? `-Wl,--export-all-symbols` : ``
    run(`$(get_compiler()) -shared -L$(julia_libdir) -o $sysimage_path $o_file -ljulia $extra`)
    return nothing
end

"""
    $SIGNATURES

lalala
"""
function restore_default_sysimg()
    if !isfile(backup_default_sysimg_path())
        error("did not find a backup sysimg")
    end
    cp(backup_default_sysimg_path(), default_sysimg_path(); force=true)
    rm(backup_default_sysimg_path())
    @info "PackageCompilerX: default sysimg restored, restart Julia for the new sysimg to be in effect"
    return nothing
end

const REQUIRES = "Requires" => UUID("ae029012-a4dd-5104-9daa-d747884805df")

# Check for things that might indicate that the app or dependencies
"""
    $SIGNATURES

Check for possible problems with regfards to relocatability at
the project at `app_dir`.

!!! warning
    This cannot guarantee that the project is free of relocatability problems,
    it can only detect some known bad cases and warn about those.
"""
function audit_app(app_dir::String)
    project_toml_path = abspath(Pkg.Types.projectfile_path(app_dir; strict=true))
    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    return audit_app(ctx)
end
function audit_app(ctx::Pkg.Types.Context)
    # Check for Requires.jl usage
    if REQUIRES in ctx.env.project.deps
        @warn "Project has a dependency on Requires.jl, code in `@require` will not be run"
    end
    for (uuid, pkg) in ctx.env.manifest
        if REQUIRES in pkg.deps
            @warn "$(pkg.name) has a dependency on Requires.jl, code in `@require` will not be run"
        end
    end

    # Check for build script usage
    if isfile(joinpath(dirname(ctx.env.project_file), "deps", "build.jl"))
        @warn "Project has a build script, this might indicate that it is not relocatable"
    end
    pkgs = Pkg.Types.PackageSpec[]
    Pkg.Operations.load_all_deps!(ctx, pkgs)
    for pkg in pkgs
        pkg_source = Pkg.Operations.source_path(pkg)
        pkg_source === nothing && continue
        if isfile(joinpath(pkg_source, "deps", "build.jl"))
            @warn "Package $(pkg.name) has a build script, this might indicate that it is not relocatable"
        end
    end
    return
end

"""
    $SIGNATURES
"""
function create_app(package_dir::String,
                    app_dir::String;
                    precompile_execution_file::Union{String, Vector{String}, Nothing}=nothing,
                    precompile_statements_file::Union{String, Vector{String}, Nothing}=nothing,
                    incremental=false,
                    filter_stdlibs=false,
                    audit=true,
                    force=false)
    project_toml_path = abspath(Pkg.Types.projectfile_path(package_dir; strict=true))
    manifest_toml_path = abspath(Pkg.Types.manifestfile_path(package_dir))
    if manifest_toml_path === nothing
        @warn "it is not recommended to create an app without a preexisting manifest"
    end
    project_toml = Pkg.TOML.parsefile(project_toml_path)
    project_path = abspath(package_dir)
    app_name = get(project_toml, "name") do
        error("expected package to have a `name`-entry")
    end
    sysimg_file = app_name * "." * Libdl.dlext

    ctx = Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
    @debug "instantiating project at \"$project_toml_path\""
    Pkg.instantiate(ctx)

    if isdir(app_dir)
        if !force
            error("directory $(repr(app_dir)) already exists, use `force=true` to overwrite (will completely",
                  " remove the directory")
        end
        rm(app_dir; force=true, recursive=true)
    end

    audit && audit_app(ctx)

    mkpath(app_dir)

    bundle_julia_libraries(app_dir)
    bundle_artifacts(ctx, app_dir)

    # TODO: Create in a temp dir and then move it into place?
    binpath = joinpath(app_dir, "bin")
    mkpath(binpath)
    cd(binpath) do
                create_sysimage(Symbol(app_name); sysimage_path=sysimg_file, project=project_path,
                                incremental=incremental, filter_stdlibs=filter_stdlibs,
                                precompile_execution_file = precompile_execution_file,
                                precompile_statements_file = precompile_statements_file)
        create_executable_from_sysimg(; sysimage_path=sysimg_file, executable_path=app_name)
        if Sys.isapple()
            cmd = `install_name_tool -change $sysimg_file @rpath/$sysimg_file $app_name`
            @debug "running $cmd"
            run(cmd)
        end
    end
    return
end

# This requires that the sysimg have been built so that there is a ccallable `julia_main`
# in Main.
function create_executable_from_sysimg(;sysimage_path::String,
                                        executable_path::String)
    flags = join((cflags(), ldflags(), ldlibs()), " ")
    flags = Base.shell_split(flags)
    wrapper = joinpath(@__DIR__, "embedding_wrapper.c")
    if Sys.iswindows()
        rpath = ``
    elseif Sys.isapple()
        rpath = `-Wl,-rpath,'@executable_path' -Wl,-rpath,'@executable_path/../lib'`
    else
        rpath = `-Wl,-rpath,\$ORIGIN:\$ORIGIN/../lib`
    end
    cmd = `$(get_compiler()) -DJULIAC_PROGRAM_LIBNAME=$(repr(sysimage_path)) -o $(executable_path) $(wrapper) $(sysimage_path) -O2 $rpath $flags`
    @debug "running $cmd"
    run(cmd)
    return nothing
end

function bundle_julia_libraries(app_dir)
    app_libdir = joinpath(app_dir, Sys.isunix() ? "lib" : "bin")
    cp(julia_libdir(), app_libdir; force=true)
    # We do not want to bundle the sysimg (nor the backup):
    rm(joinpath(app_libdir, "julia", default_sysimg_name()); force=true)
    rm(joinpath(app_libdir, "julia", backup_default_sysimg_name()); force=true)
    return
end

function bundle_artifacts(ctx, app_dir)
    @debug "bundling artifacts..."

    pkgs = Pkg.Types.PackageSpec[]
    Pkg.Operations.load_all_deps!(ctx, pkgs)

    # Also want artifacts for the project itself
    if ctx.env.pkg !== nothing
        # This is kinda ugly...
        ctx.env.pkg.path = dirname(ctx.env.project_file)
        push!(pkgs, ctx.env.pkg)
    end

    # Collect all artifacts needed for the project
    artifact_paths = String[]
    for pkg in pkgs
        pkg_source_path = Pkg.Operations.source_path(pkg)
        pkg_source_path === nothing && continue
        # Check to see if this package has an (Julia)Artifacts.toml
        for f in Pkg.Artifacts.artifact_names
            artifacts_toml_path = joinpath(pkg_source_path, f)
            if isfile(artifacts_toml_path)
                @debug "bundling artifacts for $(pkg.name)"
                artifact_dict = Pkg.Artifacts.load_artifacts_toml(artifacts_toml_path)
                for name in keys(artifact_dict)
                    @debug "  \"$name\""
                    push!(artifact_paths, Pkg.Artifacts.ensure_artifact_installed(name, artifacts_toml_path))
                end
                break
            end
        end
    end

    # Copy the artifacts needed to the app directory
    artifact_app_path = joinpath(app_dir, "artifacts")
    if !isempty(artifact_paths)
        mkpath(artifact_app_path)
    end
    for artifact_path in artifact_paths
        artifact_name = basename(artifact_path)
        # force=true?
        cp(artifact_path, joinpath(artifact_app_path, artifact_name))
    end
    return
end

end # module


text:
function bisect(f, a, b, tol)
    println("a = ", a, " b = ", b, " tol = ", tol)
    c = 0.0

    while (b - a) / 2 > tol
        c = (b + a) / 2
        if f(c) == 0
            return c
        end

        if f(a) * f(c) < 0
            b = c
        else
            a = c
        end
    end

    c = (b + a) / 2
    return c
end


text:
@testset "Dic constructors" begin
    dict = Dict(2=>3, 3=>4, 4=>2)

    @test @inferred(Dic(2=>3, 3=>4, 4=>2)) isa Dic{Int, Int}
    @test @inferred(Dic(dict)) isa Dic{Int, Int}
    @test @inferred(Dic()) isa Dic{Any, Any}

    @test @inferred(Dic{Int, Int}(2=>3, 3=>4, 4=>2)) isa Dic{Int, Int}
    @test @inferred(Dic{Int, Int}(dict)) isa Dic{Int, Int}
    @test @inferred(Dic{Int, Int}()) isa Dic{Int, Int}
end

@testset "Dic elemtary operations" begin
    d = Dic(1=>2, 2=>3, 3=>4)

    @test length(d) === 3

    @test issetequal(d, [2,3,4])
    @test issetequal(indices(d), [1,2,3])
    @test issetequal(pairs(d), [1=>2, 2=>3, 3=>4])

    @test d[1] === 2
    @test d[2] === 3
    @test d[3] === 4

    @test (d[1] = 10; d[1] === 10)
end

@testset "Dic fancy indexing" begin
    d1 = Dic(1=>2, 2=>3, 3=>1)

    @test d1[d1] isa Dic

    d2 = d1[d1]

    @test length(d2) === 3
    @test d2[1] === 3
    @test d2[2] === 1
    @test d2[3] === 2
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL, LinearAlgebra
GL = ViewerGL

const M44 = convert(GL.Matrix4, Matrix{Float64}(I,4,4))

@testset "Geometry.jl" begin

   points = [0.214359   0.261879  0.317665;
        0.0235449  0.584002  0.477773;
        0.0353606  0.245789  0.0913658;
        0.168101   0.941013  0.120331]
   # function GLHull(points::Array{Float64,2})::GL.GLMesh
   @testset "GLHull" begin
      @test typeof(GL.GLHull(points::Array{Float64,2}))==ViewerGL.GLMesh
      @test ViewerGL.GLMesh.super==Any
      @test supertype(ViewerGL.GLMesh)==Any
      @test GL.GLHull(points::Array{Float64,2}).primitive==4
      @test GL.GLHull(points::Array{Float64,2}).T==
      [1.0  0.0  0.0  0.0;
       0.0  1.0  0.0  0.0;
       0.0  0.0  1.0  0.0;
       0.0  0.0  0.0  1.0];
      @test typeof(GL.GLHull(points))==ViewerGL.GLMesh
   end

   # function GLHull2d(points::Array{Float64,2})::GL.GLMesh # points by row
   @testset "GLHull2d" begin
      points = [0.214359   0.261879;
           0.0235449  0.584002 ;
           0.0353606  0.245789 ;
           0.168101   0.941013 ]
      @test typeof(GL.GLHull2d(points))==ViewerGL.GLMesh
      @test GL.GLHull2d(points::Array{Float64,2}).primitive==1
      @test GL.GLHull2d(points::Array{Float64,2}).T==M44
   end

   # function GLPolygon(V::Lar.Points,copEV::Lar.ChainOp,copFE::Lar.ChainOp)::GL.GLMesh
   @testset "GLPolygon" begin
      V = hcat([[0,0],[1,0],[1,1],[0,1],[.25,.25],[.75,.25],[.75,.75],[.25,.75]]...)
      EV = [[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,5]]
      W = convert(Lar.Points, V')
      cop_EV = Lar.coboundary_0(EV::Lar.Cells)
      cop_EW = convert(Lar.ChainOp, cop_EV)
      V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
      V = GL.two2three(V)
      @test typeof(GL.GLPolygon(V::Lar.Points,copEV::Lar.ChainOp,copFE::Lar.ChainOp))==GL.GLMesh
      @test GL.GLPolygon(V::Lar.Points,copEV::Lar.ChainOp,copFE::Lar.ChainOp).primitive==4
      @test GL.GLPolygon(V::Lar.Points,copEV::Lar.ChainOp,copFE::Lar.ChainOp).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
   end

   # function GLPolygon(V::Lar.Points,EV::Lar.Cells)::GL.GLMesh
   @testset "GLPolygon" begin
      (V, EV) = ([0.43145 0.596771 0.758062 1.0 0.778226 0.919353 0.879033 0.806447 0.778226 0.709677 0.596771 0.262094 0.322578 0.0 0.2379 0.161291 0.467739 0.429435 0.627999 0.627999 0.383062 0.694833 0.653221 0.544027 0.778226 0.848789 0.750707 0.627999 0.694833 0.806447; 0.0 0.233873 0.11694 0.330646 0.625 0.677418 0.810484 0.717742 0.834675 0.743951 0.983871 0.810484 0.625 0.388696 0.439515 0.25403 0.467743 0.625 0.697579 0.506846 0.282258 0.677418 0.439515 0.29787 0.208757 0.439515 0.497366 0.346772 0.267455 0.368691], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 15], [15, 16], [16, 17], [17, 18], [18, 19], [19, 20], [20, 21], [21, 1], [22, 23], [24, 23], [24, 25], [25, 26], [26, 22], [27, 28], [28, 29], [29, 30], [30, 27]])
      @test typeof(GL.GLPolygon(V::Lar.Points,EV::Lar.Cells))==GL.GLMesh
      @test GL.GLPolygon(V::Lar.Points,EV::Lar.Cells).primitive==4
      @test GL.GLPolygon(V::Lar.Points,EV::Lar.Cells).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
   end

   # function GLLar2gl(V::Lar.Points, CV::Lar.Cells)::GL.GLMesh
   @testset "GLLar2gl" begin
      (V, CV) = ([0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0; 0.0 0.0 1.0 1.0 2.0 2.0 0.0 0.0 1.0 1.0 2.0 2.0; 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4, 7, 8, 9, 10], [3, 4, 5, 6, 9, 10, 11, 12]])
      @test typeof(GL.GLLar2gl(V,CV))==ViewerGL.GLMesh
      @test GL.GLLar2gl(V::Lar.Points, CV::Lar.Cells).primitive==4
      @test GL.GLLar2gl(V::Lar.Points, CV::Lar.Cells).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end

   # function GLLines(points::Lar.Points,lines::Lar.Cells)
   @testset "GLLines" begin
      (points, lines) = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0 9.0 5.0 5.0 14.0 13.0 11.0 10.0 10.0 11.0 13.0 14.0 14.0 14.0 15.0 16.0 18.0 19.0 19.0 18.0 16.0 15.0 15.0 16.0 18.0 19.0 20.0 20.0 22.0 24.0 24.0; 6.0 0.0 0.0 1.0 3.0 4.0 4.0 6.0 6.0 0.0 5.0 6.0 6.0 5.0 4.0 3.0 3.0 4.0 6.0 3.0 5.0 6.0 6.0 5.0 4.0 3.0 3.0 2.0 1.0 0.0 0.0 1.0 6.0 0.0 2.0 0.0 6.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [8, 9], [9, 10], [11, 12], [12, 13], [13, 14], [14, 15], [15, 16], [16, 17], [17, 18], [18, 11], [19, 20], [21, 22], [22, 23], [23, 24], [24, 25], [25, 26], [26, 27], [27, 28], [28, 29], [29, 30], [30, 31], [31, 32], [33, 34], [34, 35], [35, 36], [36, 37]])
      @test typeof(GL.GLLines(points::Lar.Points,lines::Lar.Cells))==ViewerGL.GLMesh
      @test GL.GLLines(points::Lar.Points,lines::Lar.Cells).primitive==1
      @test GL.GLLines(points::Lar.Points,lines::Lar.Cells).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end

   # function GLText(string)
   @testset "GLText" begin
      @test typeof(GL.GLText("Plasm"))==ViewerGL.GLMesh
      @test GL.GLText("Plasm").primitive==1
      @test GL.GLText("Plasm").T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end

   # function GLPoints(points::Lar.Points) # points by row
   @testset "GLPoints" begin
      points = rand(50,3)
      @test typeof(GL.GLPoints(points::Lar.Points))==GL.GLMesh
      @test GL.GLPoints(points::Lar.Points).primitive==0
      @test GL.GLPoints(points::Lar.Points).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end

   # function GLPolyhedron(V::Lar.Points, FV::Lar.Cells, T::GL.Matrix4=M44)
   @testset "GLPolyhedron" begin
      (V, FV) = ([0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0; 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4], [5, 6, 7, 8], [1, 2, 5, 6], [3, 4, 7, 8], [1, 3, 5, 7], [2, 4, 6, 8]])
      @test typeof(GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells))==GL.GLMesh
      @test GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells).primitive==4
      @test GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end
#=
   @testset "GLPolyhedron" begin
      (V, FV, EV) = ([0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0; 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4], [5, 6, 7, 8], [1, 2, 5, 6], [3, 4, 7, 8], [1, 3, 5, 7], [2, 4, 6, 8]], Array{Int64,1}[[1, 2], [3, 4], [5, 6
      ], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6], [3, 7], [4, 8]])
      cop_EV = Lar.coboundary_0(EV::Lar.Cells);
      cop_EW = convert(Lar.ChainOp, cop_EV);
      cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
      W = convert(Lar.Points, V');

      V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
      	W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)
      @test typeof(GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells))==GL.GLMesh
      @test GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells).primitive==4
      @test GL.GLPolyhedron(V::Lar.Points, FV::Lar.Cells).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end
=#

   # function GLGrid(V::Lar.Points,CV::Lar.Cells,color=GL.COLORS[1])
   @testset "GLGrid" begin
      V,(VV,EV,FV,CV) = Lar.cuboidGrid([1,2,1],true)
      @test typeof(GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1)))==GL.GLMesh
      @test GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1)).primitive==7
      @test GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1)).T==M44
      @test GL.GLMesh.size==48
      @test GL.GLMesh.isbitstype==false
      @test GL.GLVertexBuffer.isconcretetype==true
   end

end


text:
@testset "4.2.2.3 (g cos)^p (a+b cos)^m (c+d cos)^n" begin
    (a, c, e, f, x, ) = @variables a c e f x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*problems*of*the*form*(g*cos(e+f*x))^p*(a+b*cos(e+f*x))^m*(c+d*cos(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*cos(e+f*x))^p*(a+a*cos(e+f*x))^m*(c-c*cos(e+f*x))^n=#

    @test_int [((a + a*cos(e + f*x))^2*sec(e + f*x)^2)/(-c + c*cos(e + f*x)), x, 6, -((3*a^2*atanh(sin(e + f*x)))/(c*f)) + (4*a^2*sin(e + f*x))/(c*f*(1 - cos(e + f*x))) - (a^2*tan(e + f*x))/(c*f)]
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.2.3 (g cos)^p (a+b cos)^m (c+d cos)^n.jl


text:
module Schauder
export level, levelK, vectoroflevels, number, finger_pm, finger_permute, permutationmatrix, pickup!, pickup, drop!, drop, dropB1, dropB2, hat, fe_transf, fe_transfB2, fe_transfB1

#%  .. currentmodule:: Schauder
#%
#%  .. _modschauder:
#%
#%  Module Schauder
#%  ---------------
#%
#%  Introduction
#%  ~~~~~~~~~~~~
#%
#%  In the following ``hat(x)`` is the piecewise linear function taking values
#%  values ``(0,0), (0.5,1), (1,0)`` on the interval ``[0,1]`` and ``0`` elsewhere.
#%
#%  The Schauder basis of level ``L > 0`` in the interval ``[a,b]`` can be defined recursively
#%  from ``n = 2^L-1`` classical finite elements :math:`\psi_i(x)` on the grid
#%
#%  	``a + (1:n)/(n+1)*(b-a)``.
#%
#%
#%  Assume that f is expressed as linear combination
#%
#%  	:math:`f(x) = \sum_{i =1}^n c_i \psi_i(x)`
#%
#%  with
#%
#%  	:math:`\psi_{2j-1}(x) = hat(nx-j + 1)` 	for :math:`j = 1 \dots 2^{L-1}`
#%
#%  and
#%
#%  	:math:`\psi_{2j}(x) = hat(nx-j + 1/2)` 	for :math:`j = 1 \dots 2^{L-1}-1`
#%
#%  Note that these coefficients are easy to find for the finite element basis, just
#%
#%  .. code-block:: matlab
#%
#%  	function fe_transf(f, a,b, L)
#%  	    n = 2^L-1
#%  	    return map(f, a + (1:n)/(n+1)*(b-a))
#%  	end
#%

#%
#%  Then the coefficients of the same function with respect to the Schauder basis
#%
#%  	:math:`f(x) = \sum_{i = 1}^n c_i \phi_i(x)`
#%
#%  where for ``L = 2``
#%
#%  	:math:`\phi_2(x) = 2 hat(x)`
#%
#%  	:math:`\phi_1(x) = \psi_1(x) = hat(2x)`
#%
#%  	:math:`\phi_3(x) = \psi_3(x) = hat(2x-1)`
#%
#%  can be computed directly, but also using the recursion
#%
#%  	:math:`\phi_2(x) = \psi_1(x) + 2\psi_2(x) + \psi_2(x)`
#%
#%  This can be implemented inplace (see ``pickup()``), and is used throughout.
#%

#%
#% .. code-block:: matlab
#%
#%  	for l in 1:L-1
#%  		b = sub(c, 2^l:2^l:n)
#%  		b[:] *= 0.5
#%  		a = sub(c, 2^(l-1):2^l:n-2^(l-1))
#%  		a[:] -= b[:]
#%  		a = sub(c, 2^l+2^(l-1):2^l:n)
#%  		a[:] -= b[1:length(a)]
#%  	end
#%

#%  Reference
#%  ~~~~~~~~~
#%
#%
#%  .. function:: pickup!(x)
#%
#%  	Inplace computation of 2^L-1 Schauder-Faber coefficients from
#%  	``2^L-1`` overlapping finite-element coefficients ``x``.
#%
#%  	-- inverse of ``Schauder.drop``
#%
#%  	-- L = level(xj)
#%
#%  .. function:: drop!(x)
#%
#%  	Inplace computation of 2^L-1 finite element coefficients from
#%  	2^L-1 Faber schauder coefficients ``x``.
#%
#%  	-- inverse of ``Schauder.pickup``
#%
#%  .. function:: finger_permute(x)
#%
#%  	Reorders vector ``x`` or matrix ``A`` according to the reordering
#%  	of the elements of a Faber-Schauder-basis from
#%  	left to right, from bottom to top.
#%
#%  .. function:: finger_pm(L, K)
#%
#%  	Returns the permuation used in ``finger_permute``.
#%  	Memoized reordering of faber schauder elements from low level to high level. The last K elements/rows are left untouched.
#%

order(a,b) = a > b ? (b,a): (a,b)

function ilogbi(n::Integer)
	Base.exponent(float(n))
end

#%  .. function:: level(x)
#%
#%  	Gives the no. of levels of the biggest Schauder basis with less then length(x) elements.
#%  		level(x) = ilogb(size(x,1)+1)
#%
level(x) = ilogbi(size(x,1)+1)

#%  .. function:: level(x, K)
#%
#%  	Gives the no. of levels ``l`` of the biggest Schauder basis with less then length(x) elements
#%  	and the number of additional elements ``n-2^l+1``.
#%
function levelK(x)
	n = size(x,1)
	l = ilogbi(n+1)
	(l, n-2^l+1)
end

#%  .. function:: vectoroflevels(L, K)
#%
#%  	Gives a vector with the level of the hierarchical elements.
#%

function vectoroflevels(L, K)
	n = 2^L
	o = zeros(Int32, n-1+K)
	i = 0
	for l in 0:L-1
	   for j in 2^l:2^l:n-1
		o[j] = l
	   end
	end
	# additional elements (levels) at the end
	for l in n^2:n^2-1+K
		o[l] = l
	end
	o
end

function number(L, K)
	n = 2^L
	o = zeros(Int32, n-1+K)
	i = 0
	for l in 0:L-1
	   for j in 2^l:2^(l+1):n-1
		i+=1
		o[j] = i
	   end
	end
	for l in n:n-1+K
		i+=1
		o[l] = i
	end
	o
end
order_tbl = nothing

#memoized reordering of faber schauder elements from low level to high level
function finger_pm(L, K)
	global order_tbl
	if order_tbl == nothing
		order_tbl = Dict()
	end

	if haskey(order_tbl, (L,K))
               pm = order_tbl[(L,K)]
	else
		pm = sortperm(number(L,K))
		assign(order_tbl,pm, (L,K))

	end
	pm
end

# permutationmatrix(sortperm(number(L))
#  pm * pickup_Sigma!(fe_sigma(rand(100), 0, 1, 3),3) * pm'
function permutationmatrix(p)
	d = length(p)
	P = zeros(eltype(p), d,d)
	for i in 1:d
		P[i, p[i]] = 1
	end
	P
end

function finger_permute(x::Vector)
	(L,K) = levelK(x)
	x[finger_pm(L,K)] #undo with ipermute!(x, pm)
end

function finger_permute(Sigma::Matrix)
	(L,K) = levelK(Sigma)
	pm = permutationmatrix(finger_pm(L,K))
	pm * Sigma * pm'
end

#leaves K elements untouched
function pickup!(xj)
	L = level(xj)
	n = 2^L-1

	for l in 1:L-1
		b = sub(xj, 2^l:2^l:n)
		b[:] *= 0.5
		a = sub(xj, 2^(l-1):2^l:n-2^(l-1))
		a[:] -= b[:]
		a = sub(xj, 2^l+2^(l-1):2^l:n)
		a[:] -= b[1:length(a)]
	end
	xj
end

#leaves K elements untouched
function drop!(xj)
	L = level(xj)
	n = 2^L - 1

	for l in L-1:-1:1
		b = sub(xj, 2^l:2^l:n)
		a = sub(xj, 2^(l-1):2^l:n-2^(l-1))
		a[:] += b[:]
		a = sub(xj, 2^l+2^(l-1):2^l:n)
		a[:] += b[1:length(a)]
		b[:] *= 2

	end
	xj
end

compose(f, g) = x -> f(g(x))
drop = compose(drop!,copy)
pickup = compose(pickup!,copy)

function dropB1(yj)
	xj = copy(yj)
	(L, K) = levelK(xj)
	assert(K==1)
	n = 2^L-1
	xj[1:n] .+= xj[end]
	xj
end

function dropB2(yj)
	xj = copy(yj)
	(L, K) = levelK(xj)
	assert(K==2)
	n = 2^L-1

	for i in 1:n
		xj[i] += xj[end]*i/(n+1) + xj[end-1]*(n+1-i)/(n+1)
	end
	xj = xj[[n+1, 1:n, n+2]] #permute
end

#%  .. function:: hat(x)
#%
#%  	Hat function. Piecewise linear functions with values (-inf,0), (0,0),(0.5,1), (1,0), (inf,0).
#%  	-- ``x`` vector or number
function hat(x)
	max(1. .- abs(2.*x .- 1.),0.)
end

function fe_transf(f, a,b, L)
	n::Float64 = 2^L-1
	return map!(f, a .+ [1.:n]/(n+1.)*(b-a))
end

# versions for different basis elements

# phi1(x) = x
# phi2(x) = 1-x
function fe_transfB2(f, a,b, L)
	n::Float64 = 2^L-1
		return	[map!(f, a .+ [1.:n]/(n+1.)*(b-a)) .- [1.:n]/(n+1.)*f(a) .- [n:-1.:1.]/(n+1.)*f(b) ,f(a) , f(b)]
end
# phi1(x) = c
# periodic boundary condition
function fe_transfB1(f, a,b, L)
	n::Float64 = 2^L-1
	return	[map!(f, a .+ [1:n]/(n+1)*(b-a)) .- 0.5(f(a) + f(b)), 0.5(f(a) + f(b))]
end

include("npbayes.jl")

end


text:
"""Constant for base url and url suffix."""
const BASE_URL = "https://api.twitter.com/1.1/statuses/"
const URL_SUFFIX = ".json"

"""
    get_tweets(a::Authentictor, p::ParamsGetTweet)

Try to collect tweets based on given object of type Autenticator and ParamsGetTweet. The `get_tweets` is a wrapper around the `oauth_request_resource` method from `OAuth` package which is responsible for sending the `get` request to the `statuses/user_timeline.json`.

# Example
```
# Here Twitter OAuth token and secret is stored as environment varialbe whith name
# CKEY, CSEC, OTOK, OSEC
julia> creds = Authentictor(ENV["CKEY"], ENV["CSEC"], ENV["OTOK"], ENV["OSEC"])
julia> params = ParamsGetTweet("realDonaldTrump", count=100)
julia> tweets = get_tweets(creds, params);
```
That should return 100 latest tweets from the profile @realDonaldTrump.
"""
function get_tweets(a::Authentictor, p::ParamsGetTweets)
    url = joinpath(BASE_URL, "user_timeline.json")
    # paramter convertion to dictionary
    params = params_to_dict(p)
    # Response Object
    resp = oauth_request_resource(
        url,
        "GET",
        params,
        a.consumer_key,
        a.consumer_secret,
        a.oauth_token,
        a.oauth_secret
    )
    if resp.status == 200
        newdata = JSON.parse(String(resp.body))
        return newdata
    else
        error("Twitter API returned $(r.status) status")
    end
end

"""
    post_tweet(a::Authentictor, p::ParamsPostTweet)

Try to collect tweets based on given object of type Autenticator and ParamsPostTweet. The `post_tweet` is a wrapper around the `oauth_request_resource` method from `OAuth` package which is responsible for sending the `get` request to the `statuses/update.json`.

# Example
```
# Here Twitter OAuth token and secret is stored as environment varialbe whith name
# CKEY, CSEC, OTOK, OSEC
julia> creds = Authentictor(ENV["CKEY"], ENV["CSEC"], ENV["OTOK"], ENV["OSEC"])
julia> resource_params = ParamsPostTweet("asdf asdf asdf Test tweet")
julia> tweets = post_tweet(creds, resource_params);
```
That should return 100 latest tweets from the profile @realDonaldTrump.
"""
function post_tweet(a::Authentictor, p::ParamsPostTweet)
    url = joinpath(BASE_URL, "update.json")
    # paramter convertion to dictionary
    params = params_to_dict(p)
    # Response Object
    resp = oauth_request_resource(
        url,
        "POST",
        params,
        a.consumer_key,
        a.consumer_secret,
        a.oauth_token,
        a.oauth_secret
    )
    if resp.status == 200
        newdata = JSON.parse(String(resp.body))
        return newdata
    else
        error("Twitter API returned $(r.status) status")
    end
end


text:
# GROUP BY clause.

mutable struct GroupClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    by::Vector{SQLClause}

    GroupClause(;
                over = nothing,
                by = SQLClause[]) =
        new(over, by)
end

GroupClause(by...; over = nothing) =
    GroupClause(over = over, by = SQLClause[by...])

"""
    GROUP(; over = nothing, by = [])
    GROUP(by...; over = nothing)

A `GROUP BY` clause.

# Examples

```jldoctest
julia> c = FROM(:person) |>
           GROUP(:year_of_birth) |>
           SELECT(:year_of_birth, AGG("COUNT", OP("*")));

julia> print(render(c))
SELECT
  "year_of_birth",
  COUNT(*)
FROM "person"
GROUP BY "year_of_birth"
```
"""
GROUP(args...; kws...) =
    GroupClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(GROUP), pats::Vector{Any}) =
    dissect(scr, GroupClause, pats)

function PrettyPrinting.quoteof(c::GroupClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(GROUP))
    append!(ex.args, quoteof(c.by, ctx))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::GroupClause, c′) =
    GroupClause(over = rebase(c.over, c′), by = c.by)


text:
using cme257hw6
using Base.Test

function R(z::Complex)
    return z
end

@test cme257hw6.jset(R, collect(10:100), collect(0:100), 100, 10.) == zeros(91,101)
#this should be the zero matrix; everything immediately falls outside of the convergence theshhold

@test cme257hw6.jset(R, collect(10:100), collect(0:100), 100, 1000.) == 101*ones(91,101)
#everything begins in the convergence area, so everything stays convergent.


text:
module TestDataMigration
using Test

using Catlab, Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra
using Catlab.Graphs.BasicGraphs: TheoryGraph, TheoryWeightedGraph
using Catlab.Graphs.BipartiteGraphs: TheoryUndirectedBipartiteGraph

# Pullback data migration
#########################

@present TheoryDDS(FreeSchema) begin
  X::Ob
  Φ::Hom(X,X)
end

@abstract_acset_type AbstractDDS
@acset_type DDS(TheoryDDS, index=[:Φ]) <: AbstractDDS

h = Graph(3)
add_parts!(h, :E, 3, src = [1,2,3], tgt = [2,3,1])

# Identity data migration.
@test h == Graph(h, Dict(:V => :V, :E => :E),
                    Dict(:src => :src, :tgt => :tgt))

# Migrate DDS → Graph.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[2,3,1])
X = TheoryDDS[:X]
@test h == Graph(dds, Dict(:V => :X, :E => :X),
                 Dict(:src => id(X), :tgt => :Φ))

h2 = copy(h)
migrate!(h2, dds, Dict(:V => :X, :E => :X),
                  Dict(:src => id(X), :tgt => :Φ))
@test h2 == ob(coproduct(h, h))

# Migrate DDS → DDS by advancing four steps.
@test dds == DDS(dds, Dict(:X => :X),
                      Dict(:Φ => [:Φ, :Φ, :Φ, :Φ]))

@present TheoryLabeledDDS <: TheoryDDS begin
  Label::AttrType
  label::Attr(X, Label)
end
@acset_type LabeledDDS(TheoryLabeledDDS, index=[:Φ, :label])

S, ϕ, Label, label = generators(TheoryLabeledDDS)
V, E, s, t, Weight, weight = generators(TheoryWeightedGraph)

ldds = LabeledDDS{Int}()
add_parts!(ldds, :X, 4, Φ=[2,3,4,1], label=[100, 101, 102, 103])

wg = WeightedGraph{Int}(4)
add_parts!(wg, :E, 4, src=[1,2,3,4], tgt=[2,3,4,1], weight=[101, 102, 103, 100])

@test wg == WeightedGraph{Int}(ldds,
  Dict(:V => :X, :E => :X),
  Dict(:src => id(S), :tgt => :Φ, :weight => [:Φ, :label]))

@test Presentation(Graph(1)) == TheoryGraph
@test Presentation(ldds) == TheoryLabeledDDS

F = FinFunctor(
  Dict(V => S, E => S),
  Dict(s => id(S), t => ϕ, weight => compose(ϕ, label)),
  TheoryWeightedGraph, TheoryLabeledDDS
)
ΔF = Delta(F, LabeledDDS{Int}, WeightedGraph{Int})
@test wg == ΔF(ldds)

idF = FinFunctor(
  Dict(X => X, Label => Label),
  Dict(ϕ => ϕ, label => label),
  TheoryLabeledDDS, TheoryLabeledDDS
)
@test ldds == LabeledDDS{Int}(ldds, idF)

# Left pushforward data migration
#################################

Σ = Sigma

V1B, V2B, EB = generators(TheoryUndirectedBipartiteGraph, :Ob)
srcB, tgtB = generators(TheoryUndirectedBipartiteGraph, :Hom)

VG, EG = generators(TheoryGraph, :Ob)
srcG, tgtG = generators(TheoryGraph, :Hom)

F = FinFunctor(
  Dict(V1B => VG, V2B => VG, EB => EG),
  Dict(srcB => srcG, tgtB => tgtG),
  TheoryUndirectedBipartiteGraph, TheoryGraph
)

idF = FinFunctor(
  Dict(VG => VG, EG => EG),
  Dict(srcG => srcG, tgtG => tgtG),
  TheoryGraph, TheoryGraph
)

ΣF = Σ(F, UndirectedBipartiteGraph, Graph)
X = UndirectedBipartiteGraph()

Y = ΣF(X)
@test nparts(Y, :V) == 0
@test nparts(Y, :E) == 0

X = @acset UndirectedBipartiteGraph begin
  V₁ = 4
  V₂ = 3
  E = 4
  src = [1,2,2,3]
  tgt = [1,1,2,3]
end

Y = ΣF(X)
@test nparts(Y, :V) == 7
@test nparts(Y, :E) == 4
@test length(Y[:src] ∩ Y[:tgt]) == 0

@test Σ(idF, Graph, Graph)(Y) == Y

@present ThSpan(FreeSchema) begin
  (L1, L2, A)::Ob
  l1::Hom(A, L1)
  l2::Hom(A, L2)
end
@acset_type Span(ThSpan, index=[:l1, :l2])

X = @acset Span begin
  L1 = 3
  L2 = 4
  A = 3
  l1 = [1,1,2]
  l2 = [1,2,3]
end

@present ThInitial(FreeSchema) begin
  I::Ob
end
@acset_type Initial(ThInitial)

L1, L2, A = generators(ThSpan, :Ob)
l1, l2 = generators(ThSpan, :Hom)
I = generators(ThInitial)[1]

bang = FinFunctor(
  Dict(L1 => I, L2 => I, A => I),
  Dict(l1 => id(I), l2 => id(I)),
  ThSpan, ThInitial
)

Σbang = Σ(bang, Span, Initial)
Y = Σbang(X)

@test nparts(Y, :I) == 4

vertex = FinFunctor(Dict(I => VG), Dict(), ThInitial, TheoryGraph)
edge = FinFunctor(Dict(I => EG), Dict(), ThInitial, TheoryGraph)

Z = Σ(vertex, Initial, Graph)(Y)
@test nparts(Z, :V) == 4
@test nparts(Z, :E) == 0

Z = Σ(edge, Initial, Graph)(Y)
@test nparts(Z, :V) == 8
@test nparts(Z, :E) == 4
@test Z[:src] ∪ Z[:tgt] == 1:8

end


text:
# INSTRUCTIONS:

# Running this file is fairly straightforward:
# From the terminal call julia with this file as the first argument, as well as one of three other arguments:
# 1. A number, corresponding to the ID of a SuiteSparse: Matrix Collection matrix.
# 2. The path to an .mtx file.
# 3. The path to a file containing either of the above options on each line.

# This would look something like
# >julia benchmarks.jl 1375
# or
# >julia benchmarks.jl ~/mymtx.mtx

# CHANGING THE SHARED LIBRARY.
# Further instructions on changing the shared library programmatically can be found in the docs.
# However, simply changing the LocalPreferences.toml file will suffice for this benchmark script.

# Some options can be found further down under SETTINGS
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using SuiteSparseMatrixCollection
using MatrixMarket
using SuiteSparseGraphBLAS
using BenchmarkTools
using SparseArrays
using LinearAlgebra

#OPTIONS SET 1:
# Maximum number of samples taken for each benchmark
BenchmarkTools.DEFAULT_PARAMETERS.samples = 10
# Total amount of time allowed for each benchmark, minimum of 1 sample taken.
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 60

# Comment or uncomment this line to disable or enable MKLSparse respectively.
# This will only work for SpMM and SpMV and only operates on CSC.
using MKLSparse

# Change this to change the size of the dense RHS of csrtimesfull and csctimesfull
const sizefullrhs = [1,2,4]

const suite = BenchmarkGroup()
const ssmc = ssmc_db()

function AxB_allbycol(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: CSC = CSC * Full\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :bycol)
    gbset(G, :format, :bycol) #set to CSC

    C = GBMatrix{eltype(G)}(size(G, 1), size(m2, 2))
    gbset(C, :format, :bycol)

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)

    printstyled("\nGBMatrix:\n", bold=true)
    for n ∈ nthreads
        printstyled("\nC = S * F with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        #print burble for checking
        gbset(:burble, true)
        mul!(C, G, m2)
        gbset(:burble, false)
        B = @benchmark mul!($C, $G, $m2)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
    gbset(G, :format, :byrow) #set back to CSR
end

function AxB_allbyrow(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: CSR = CSR * Full\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :byrow)
    gbset(G, :format, :byrow)

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)

    printstyled("\nGBMatrix:\n", bold=true)
    for n ∈ nthreads
        printstyled("\nC' = S' * F' with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        #print burble for checking
        gbset(:burble, true)
        mul(G, m2)
        gbset(:burble, false)
        B = @benchmark mul($G, $m2)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
end

function AxB_ColxRow(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: ByRow = CSC * Full_byrow\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :byrow)
    gbset(G, :format, :bycol)

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)

    printstyled("\nGBMatrix:\n", bold=true)

    for n ∈ nthreads
        printstyled("\nC' = S * F' with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        #print burble for checking
        gbset(:burble, true)
        mul(G, m2)
        gbset(:burble, false)
        B = @benchmark mul($G, $m2)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
    gbset(G, :format, :byrow) #set back to CSR
end

function CaccumAxB_allbycol(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: Full += CSC * Full\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :bycol)
    gbset(G, :format, :bycol) #set to CSC

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)

    C = GBMatrix(size(G, 1), size(m2, 2), 0.0)
    gbset(C, :sparsity_control, :full)
    gbset(C, :format, :bycol)

    printstyled("\nGBMatrix:\n", bold=true)
    #print burble for checking

    for n ∈ nthreads
        printstyled("\nF += S * F with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        gbset(:burble, true)
        mul!(C, G, m2; accum=+)
        gbset(:burble, false)
        B = @benchmark mul!($C, $G, $m2; accum=+)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
    gbset(G, :format, :byrow) #set back to CSR
end

function CaccumAxB_allbyrow(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: Full_byrow += CSR * Full_byrow\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :byrow)
    gbset(G, :format, :byrow)

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)
    gbset(:burble, true)
    C = GBMatrix(size(G, 1), size(m2, 2), 0.0)
    gbset(C, :sparsity_control, :full)
    gbset(C, :format, :byrow)

    printstyled("\nGBMatrix:\n", bold=true)
    #print burble for checking
    gbset(:burble, false)
    for n ∈ nthreads
        printstyled("\nF' += S' * F' with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        gbset(:burble, true)
        mul!(C, G, m2; accum=+)
        gbset(:burble, false)
        B = @benchmark mul!($C, $G, $m2; accum=+)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
end

function CaccumAxB_CRC(S, G, nthreads, sizerhs)
    printstyled("\nMETHOD: Full_bycol += CSR * Full_bycol\n", color=:green)
    GC.gc()
    m = rand(size(S, 2), sizerhs)
    m2 = GBMatrix(m)
    println("Size of dense matrix is $(size(m))")

    gbset(m2, :format, :bycol)
    gbset(G, :format, :byrow)

    printstyled("\nSparseMatrixCSC:\n", bold=true)
    A = @benchmark $S * $m
    show(stdout, MIME("text/plain"), A)
    gbset(:burble, true)
    C = GBMatrix(size(G, 1), size(m2, 2), 0.0)
    gbset(C, :sparsity_control, :full)
    gbset(C, :format, :bycol)

    printstyled("\nGBMatrix:\n", bold=true)
    #print burble for checking
    gbset(:burble, false)
    for n ∈ nthreads
        printstyled("\nF += S' * F with $n threads: \n", bold=true, color=:green)
        gbset(:nthreads, n)
        gbset(:burble, true)
        mul!(C, G, m2; accum=+)
        gbset(:burble, false)
        B = @benchmark mul!($C, $G, $m2; accum=+)
        show(stdout, MIME("text/plain"), B)
        tratio = ratio(median(A), median(B))
        color = tratio.time >= 1.0 ? :green : :red
        printstyled("\nMedian speedup over SparseArrays using $n threads is: $(string(tratio))\n"; bold=true, color)
    end
end

# OPTIONS SET 2:
# run these functions for benchmarking:
const functorun = [AxB_allbycol, AxB_ColxRow, CaccumAxB_allbycol, CaccumAxB_allbyrow, CaccumAxB_CRC]
#= The choices are:
AxB_allbycol - S * F
AxB_ColxRow - S' * F
CaccumAxB_allbycol - F += S * F
CaccumAxB_allbyrow - F' += S' * F'
CaccumAxB_CRC - F += S' * F
Please open an issue or message me for further functions to add here.
=#

# run with these nthread settings, add or remove to/from vector.
const threadlist = [1, 4, 8, Sys.CPU_THREADS ÷ 2, Sys.CPU_THREADS]

function singlebench(pathornum)
    x = tryparse(Int64, pathornum)
    if x !== nothing
        ssmc[x, :real] == true || throw(ArgumentError("SSMC ID must be for a matrix with real values"))
        path = joinpath(fetch_ssmc(ssmc[x, :group], ssmc[x, :name]), "$(ssmc[x, :name]).mtx")
    elseif isfile(pathornum)
        path = pathornum
    else
        throw(ErrorException("Argument is not a path or SuiteSparseMatrixCollection ID number"))
    end
    name = basename(path)
    mmpath = @time MatrixMarket.mmread(path)
    S = convert(SparseMatrixCSC{Float64}, mmpath)
    G = GBMatrix(S)
    gbset(G, :format, :byrow)
    diag(G)
    printstyled("\n#################################################################################\n"; bold=true, color=:green)
    printstyled("Benchmarking $name:\n"; bold=true, color=:green)
    printstyled("#################################################################################\n"; bold=true, color=:green)
    for i ∈ sizefullrhs
        printstyled("\n================================================================================= $name : Using a size $i B matrix"; bold=true, color=:green)
        for f ∈ functorun
            f(S, G, threadlist, i)
        end
    end
end

if length(ARGS) != 0
    if isfile(ARGS[1])
        if splitext(ARGS[1])[2] == ".mtx"
            singlebench(ARGS[1])
        else
            lines = readlines(ARGS[1])
            filter!((x) -> !occursin("#", x), lines)
            singlebench.(lines)
        end
    elseif tryparse(Int64, ARGS[1]) !== nothing
        singlebench(ARGS[1])
    else
        throw(ArgumentError("The first argument must a file with a list of SuiteSparse ID numbers or paths to MatrixMarket files"))
    end
end


text:
function rewire!(network, nodes, ID, options; remove=false)
    isempty(options) && return

    target = rand(options)

    # number of values
    val_range = axes(nodes[ID].values, 1)
    val_norm = 1 / (val_range[end] * 10)
    p = val_norm * sum(
        i->abs(nodes[ID].values[i]-nodes[target].values[i]),
        val_range)
    if !remove
        p = 1 - p
    end

    # p=1-(((abs(nodes[ID].values[1]-nodes[target].values[1]))
    #     +(abs(nodes[ID].values[2]-nodes[target].values[2]))
    #     +(abs(nodes[ID].values[3]-nodes[target].values[3]))
    #     +(abs(nodes[ID].values[4]-nodes[target].values[4]))
    #     +(abs(nodes[ID].values[5]-nodes[target].values[5])))
    #     /50)
    if rand() ≤ p
        #Add pr remove edge
        if remove
            rem_edge!(network, target, ID)
        else
            add_edge!(network, target, ID)
        end

        #Adjust energy
        if nodes[target].vote[] == nodes[ID].vote[]
            j = remove ? J : -J
            nodes[ID].energy[] += j
        end
    end
end

function procedure2(nodes, network, data; ID, N, T)

    #Add new row to data frame
    push!(data, zeros(10+1))

    #Select unconnected node, connect with probability w
    options = Int[]
    for i in 1:N
        if has_edge(network, i, ID) == false
            push!(options, i)
        end
    end

    rewire!(network, nodes, ID, options, remove=false)

    #Select edge, disconnect with probabilty w^-1
    options = inneighbors(network, ID)

    rewire!(network, nodes, ID, options, remove=true)

    # if options != []
    #     target = options[rand(1:end)] #Choose a random in-neighbor
    #     p=(((abs(nodes[ID].values[1]-nodes[target].values[1]))
    #     +(abs(nodes[ID].values[2]-nodes[target].values[2]))
    #     +(abs(nodes[ID].values[3]-nodes[target].values[3]))
    #     +(abs(nodes[ID].values[4]-nodes[target].values[4]))
    #     +(abs(nodes[ID].values[5]-nodes[target].values[5])))
    #     /50)
    #     if rand(0.0:1.0)<=p
    #         #Remove edge
    #         rem_edge!(Network,target,ID)

    #         #Adjust energy
    #         if nodes[target].vote == nodes[ID].vote
    #             nodes[ID].energy += J
    #         end
    #     end
    # end

    #Select new random preference & track preference change
    oldVal = nodes[ID].vote[]
    newVal = rand(1:10)
    nodes[ID].vote[] = newVal

    #Substract current node energy from former node energy
    ΔE = dE(nodes, ID, network) - nodes[ID].energy[]

    p = -(ΔE)/T

    #If ΔE<0 apply it:
    if ΔE < 0
        nodes[ID].energy[] += ΔE
        for i in outneighbors(network, ID)
            nodes[i].energy[] = dE(nodes, i, oldVal, newVal)
        end
    #If ΔE>0 apply it with following probability:
    elseif rand() < exp(p)
        nodes[ID].energy[] += ΔE
        for i in outneighbors(network, ID)
            nodes[i].energy[] = dE(nodes, i, oldVal, newVal)
        end
    else
        nodes[ID].vote[] = oldVal
    end

    #Log preferences
    trackPreference!(data, oldVal, nodes[ID].vote[])

    #Compute system energy
    computeEnergy!(data, nodes)

    #Optional: change core opinions
end

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
module TermoQuimica

using Unitful
using NLsolve

export EOS, general, vanlaar, ideal, margules, wilson, wilsonmod, nrtl, uniquac, unifac, unifacmod, saft
include("EOS.jl")
include("general.jl")
include("ideal.jl")
include("margules.jl")
include("vanlaar.jl")
include("wilson.jl")
include("wilsonmod.jl")
include("nrtl.jl")
include("uniquac.jl")
include("unifac.jl")
include("unifacmod.jl")
include("saft.jl")

end # module

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


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

using Markdown

# ╔═╡ 8b27c417-eb11-4bb2-b11d-3835de736c35
md"""
### Notebook 1
#### 1.1
Buscar `convert` y `parse` en la documentación.
"""

# ╔═╡ 88aa3f49-0f95-43aa-bd7a-8c1d79f2d341
?convert;
# ?parse

# ╔═╡ 2ea62f54-efd4-420c-909e-0003784d8dc1
md"""
#### 1.2
Asignar `365` a una variable llamada `days`. Convierte `days` a flotante.
"""

# ╔═╡ 82b41bf8-5766-4fd2-807f-7e0dba424633
days = 365
convert(Float64, days)

# ╔═╡ de189aeb-361d-4f47-96ae-2a670056a590
md"""
#### 1.3
Ver que sucede ejectuando

```julia
convert(Int64, '1')
```
and

```julia
parse(Int64, '1')
```
"""

# ╔═╡ 73f6a473-6aa8-4f08-9aea-e74ce01e1b7b
# Esto regresa el código ascii (un entero) asociado con el caracter '1'
convert(Int64, '1')

# ╔═╡ 2a6dceb7-f609-499a-b736-117b72380ae2
# Esto regresa un entero encapsulado entre comillas
parse(Int64, '1')

# ╔═╡ 70b7b5fe-f290-4087-927e-21c80320548b
md"""
### Notebook 2
#### 2.1
Crea una cadena que dice "hola" 1000 veces
"""

# ╔═╡ 6e768f73-8239-4d14-ae4c-adfe7c9c77f7
"hola"^1000

# ╔═╡ a176f8f7-afc4-43d3-a960-3e3a6ea7ff79
md"""
#### 2.2
Agrega dos números dentro de una cadena
"""

# ╔═╡ 09643980-8c7b-40ce-bb73-9565408555bf
m, n = 1, 1
"$m + $n = $(m + n)"

# ╔═╡ b533f6f7-63c1-415d-acb9-0ec8790976ac
md"""
### Notebook 3

#### 3.1
Crea un arreglo, `arreglo`, que es un arreglo en 1D de 2 elementos, cada uno conteniendo el número 0.
Agrega a `arreglo` para agregar un segundo número, `1`, a cada arreglo.
"""

# ╔═╡ da8e8e66-ed67-4064-9d66-ba829c350380
a_ray = [[0], [0]]
push!(a_ray[1], 1)
push!(a_ray[2], 1)
a_ray

# ╔═╡ 3b35ffd4-15a4-46f5-8cdf-191411ec46db
md"""
### 3.2
Trata de agregar "Emergencia" a `miagenda` con el valor `911`. Trata de agregar `911` como un entero y no como cadena. ¿Porqué no funciona?
"""

# ╔═╡ 67cc87a9-abe7-42eb-bf27-b0f1aa7fc4ed
miagenda = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ╔═╡ 0438c676-2e85-45bb-8a90-0e692ddb4f71
miagenda["Emergency"] = 911
#=

Julia infiere que "miagenda" toma ambos llaves
y valores del tipo "String". Podemos ver que miagenda
es un Dict{String,String} con 2 entradas. Esto significa que
Julia no va a aceptar enteros como valores en miagenda.

=#

# ╔═╡ 26508993-2fea-4de0-8b3f-366d1ef77eb7
md"""
#### 3.3
Crea un nuevo diccionario que se llame `agenda_flexible` que tenga el número de Jenny guardado como cadena y el de los Cazafantasmas como entero.
"""

# ╔═╡ d8351db6-dc61-45c2-a771-0a353cb6ceb7
agenda_flexible = Dict("Jenny" => "867-5309", "Ghostbusters" => 5552368)

# ╔═╡ 5560c6a3-a7e1-4963-bc16-4999b1eb1157
md"""
#### 3.4
Add the key "Emergency" with the value `911` (an integer) to `flexible_phonebook`.
"""

# ╔═╡ fc29f9aa-b3b0-4fd2-b0bc-96301f950600
flexible_phonebook["Emergency"] = 911

# ╔═╡ d77e981c-9440-4882-abd5-ca391a054621
md"""
##### 3.5
3.5 ¿Porqué podemos agregar un entero como valor a `agenda_flexible` pero no a  `miagenda`? ¿Cómo pudimos haber inicializado `miagenda` para que aceptara enteros como valores?
"""

# ╔═╡ b1830c1b-3f00-42f9-81bf-adecfa05c161
#=

Julia infiere que miagenda_flexible toma valores del tipo
Any. A diferencia de miagenda, miagenda_flexible es un
Dict{String,Any} con 2 entradas.

Para evitar esto, podemos inicializar miagenda a un
diccionario vacío y agregamos entradas después. O podemos
decirle a Julia explícitamente que queremos un diccionario
que acepte objectos del tipo Any como valores.

¡Ve los ejemplos!

=#

miagenda = Dict()

# ╔═╡ 03aac330-27fc-459e-8261-0fdac3721197
miagenda = Dict{String, Any}("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ╔═╡ af8072fd-3175-49c6-9248-41f48588d5a2
md"""
### Notebook 4

#### 4.1

4.1 Crea un diccionario `squares`, que tiene llaves de valores de 1 a 100. El valor asociado a cada llave es el cuadrado de la llave. Guarda los valores asociados a las llaves pares como enteros y las impares como cadenas. Por ejemplo,

```julia
squares[10] == 100
squares[11] == "121"
```

(¡No necesitas condicionales para esto!)
"""

# ╔═╡ 72506456-8104-45b8-8243-722bd3bbce1f
squares = Dict()
iterable = range(1, 2, 50)
#  Otra opción
# iterable = 1:2:99
for key in iterable
    squares[key] = "$(key^2)"
    squares[key + 1] = (key + 1)^2
end

@show squares[10]
@show squares[11]

# ╔═╡ ace49f1a-2149-4d68-81ea-965d4d0e9197
md"""
#### 4.2

4.2 Usa `fill` para crea una matriz de `10x10` de solo `0`'s. Pobla las primeras 10 entradas con el índice de esa entrada. ¿Julia usa el orden de primero columna o primero renglón? (O sea, ¿el "segundo" elemento es el de la primera columna en el primer renglón, ó es el de el primer renglón en la segunda columna?)
"""

# ╔═╡ 6c9e8943-0997-4f01-bc16-68bf62f72d21
A = fill(0, (10, 10))
for i in 1:10
    A[i] = i
end
A
# ¡Julia usa order columnar! Los primeros 10 elementos van a poblar la primera columna.

# ╔═╡ a544197a-a657-453c-80d5-b08e0d2dcdf6
md"""
### Notebook 5

#### 5.1

5.1 Reescribe FizzBuzz sin usar `elseif`.
"""

# ╔═╡ 2e786f7a-e0a4-4e7c-90b6-7081864b6fa2
N = 16
if (N % 3 == 0) & (N % 5 == 0)
    println("FizzBuzz")
else
    if (N % 3 == 0)
        println("Fizz")
    else
        if (N % 5 == 0)
            println("Buzz")
        else
            println(N)
        end
    end
end

# ╔═╡ 6b91710e-5b70-434e-917c-d8c1e5f9809f
md"""
#### 5.2

Reescribe FizzBuzz usando el operador ternario.
"""

# ╔═╡ 3add75d7-88e6-4eb3-ad52-ace52755cf58
((N % 3 == 0) & (N % 5 == 0)) ? println("FizzBuzz") : ((N % 3 == 0) ? println("Fizz") : ((N % 5 == 0) ? println("Buzz") : println(N)))

# ╔═╡ 71cd7876-82ca-4b09-b5b9-d79160b8061b
md"""
### Notebook 6

#### 6.1

6.1 En vez de broadcastear `f` sobre `v`,  pudimos haber hecho `v .^ 2`.

Sin declarar una nueva funcion, agrega 1 a cada elemento de una matriz de `3x3` llena de `0`'s.
"""

# ╔═╡ 4594cabf-46b0-42dc-b027-92c42de34d15
fill(0, (3, 3)) .+ 1

# ╔═╡ d36ccfa0-6a44-4644-aa69-a46c09668c4d
md"""
#### 6.2

6.2 En vez de broadcastear `f` sobre el vector `v` con la sintaxis de punto, aplica `f` a todos los elementos de`v` usando `map` como función.
"""

# ╔═╡ c0ab3d4b-93f1-4588-b7fa-c007039889a4
f(x) = x^2
v = [1, 2, 3]
map(f, v)

# ╔═╡ 06339cb3-9e6b-4214-9276-86574cc5d6bb
md"""
#### 6.3

Una cifra de César recorre cada letra un número determinado de plazas más adelante en el abecedario. Un corrimiento, o shift, de 1 manda "A" a "B". Escribe una función llamada `cesar` que toma una cadena como input y un corrimiento y regresa una cadena desencriptada tal que obtengas

```julia
cesar("abc", 1)
"bcd"

cesar("hello", 4)
"lipps"
```
"""

# ╔═╡ bb0149cd-6e2d-4bf1-af68-71b3ddae9a61
caesar(input_string, shift) = map(x -> x + shift, input_string)

# ╔═╡ 83177d5d-28a4-4d52-8df4-4b316e5e3080
md"""
### Notebook 7

#### 7.1

7.1 Usa el paquete de  (código fuente en https://github.com/JuliaMath/Primes.jl) para encontrar el número primer más grande menor a 1,000,000
"""

# ╔═╡ d09fbb3f-9cd8-4c8a-b2f8-2258579a71b5
#Pkg.add("Primes")
using Primes
maximum(primes(1000000))

# ╔═╡ e5de2132-0a97-4802-a10a-cd0e770c8985
md"""
### Notebook 8

#### 8.1

8.1 Grafica y vs x para `y = x^2` usando el backend de PyPlot.
"""

# ╔═╡ 77929858-3d6a-4832-ac1a-726ef652aea9
using Plots
pyplot()
x = 1:10
y = x .^ 2
plot(x, y)

# ╔═╡ bb14c31f-4d17-4fe6-86a0-7c2799c078f4
md"""
### Notebook 9

#### 9.1

Agrega un método para `+` que aplique un cifrado de César a una cadena (cómo en el notebook 6) tal que

```julia
"hello" + 4 == "lipps"
```
"""

# ╔═╡ 085df56e-9938-45bb-944b-21173b3f6c76
import Base: +
+(x::String, y::Int) = map(x -> x + y, x)

# ╔═╡ a515c3be-a27b-4cbf-9417-7b8446399ad2
md"""
#### 9.2

Checa que has extendido propiamente `+` recorriendo la próxima cadena para atrás por 3 letras:

"Gr#qrw#phggoh#lq#wkh#diidluv#ri#gudjrqv#iru#|rx#duh#fuxqfk|#dqg#wdvwh#jrrg#zlwk#nhwfkxs1"
"""

# ╔═╡ 05f18b69-67ae-40fc-aefd-622af5f46a96
"Gr#qrw#phggoh#lq#wkh#diidluv#ri#gudjrqv#iru#|rx#duh#fuxqfk|#dqg#wdvwh#jrrg#zlwk#nhwfkxs1" + -3

# ╔═╡ c59e1737-d5fa-4cfe-9bd8-f0c8eb677231
md"""
### Notebook 10

#### 10.1

10.1 Usa `circshift` para obtener una matriz con las columnas de A cíclicamente recorridas a la derecha por 3 columnas.

Empezando con

```
A = [
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 ]
```

Quieres obtener

```
A = [
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 ]
```
"""

# ╔═╡ 242397e3-e3e7-4f76-817c-2bccfb69009a
circshift(A, (0, 4))

# ╔═╡ 79802f7a-50f2-4ae4-8991-046cd8fff9c6
md"""
#### 10.2

10.2 Toma el producto de un vector `v` con sí mismo.
"""

# ╔═╡ fd8a44d2-ca0b-469c-a33e-ecf02dec6955
v = [1, 2, 3]
v * v'

# ╔═╡ 0e6ddaa6-3d8c-4db7-904a-e8eb76acb083
md"""
#### 10.3
10.3 Toma el producto de un vector `v` con sí mismo.
"""

# ╔═╡ ec5e8a05-43f4-4f8a-a6e0-cdc7213a2d82
v' * v

# ╔═╡ 595844c9-cc40-4aaf-862a-e8becf9b4f0c
md"""
### Notebook 11

#### 11.1

 ¿Cuáles son los eigenvalores de la Matriz A

```
A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]
```
"""

# ╔═╡ f4e78715-446b-4434-b10e-a39ed3cac57c
A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]

# ╔═╡ bc2e20b3-7f7d-418e-a5ac-0bb90bfcecbc
eigdec = eigfact(A)
eigdec[:values]

# ╔═╡ f66596d8-a5c2-4025-995e-a0a1a3ebbade
md"""
#### 11.2

Crea una matriz diagonal de los eigenvalores de A
"""

# ╔═╡ 3881814e-0956-4588-831c-a6afc18c4290
Diagonal(eigdec[:values])

# ╔═╡ 985f658f-81ad-46d9-8158-673a10c77628
md"""
#### 11.3

Realiza un factorización de Hessenberg sobre la matriz A. Verifica que `A = QHQ'`.
"""

# ╔═╡ dd6cf48b-a756-4d1f-95a9-89649892bab7
F = hessfact(A)

# ╔═╡ 6cf155cd-07e5-466d-9e98-903450b969fc
isapprox(A, F[:Q] * F[:H] * F[:Q]')

# ╔═╡ 3426c47a-b63c-4d3c-89d3-24bbd751300c

# ╔═╡ Cell order:
# ╟─8b27c417-eb11-4bb2-b11d-3835de736c35
# ╠═88aa3f49-0f95-43aa-bd7a-8c1d79f2d341
# ╟─2ea62f54-efd4-420c-909e-0003784d8dc1
# ╠═82b41bf8-5766-4fd2-807f-7e0dba424633
# ╟─de189aeb-361d-4f47-96ae-2a670056a590
# ╠═73f6a473-6aa8-4f08-9aea-e74ce01e1b7b
# ╠═2a6dceb7-f609-499a-b736-117b72380ae2
# ╟─70b7b5fe-f290-4087-927e-21c80320548b
# ╠═6e768f73-8239-4d14-ae4c-adfe7c9c77f7
# ╟─a176f8f7-afc4-43d3-a960-3e3a6ea7ff79
# ╠═09643980-8c7b-40ce-bb73-9565408555bf
# ╟─b533f6f7-63c1-415d-acb9-0ec8790976ac
# ╠═da8e8e66-ed67-4064-9d66-ba829c350380
# ╟─3b35ffd4-15a4-46f5-8cdf-191411ec46db
# ╠═67cc87a9-abe7-42eb-bf27-b0f1aa7fc4ed
# ╠═0438c676-2e85-45bb-8a90-0e692ddb4f71
# ╟─26508993-2fea-4de0-8b3f-366d1ef77eb7
# ╠═d8351db6-dc61-45c2-a771-0a353cb6ceb7
# ╟─5560c6a3-a7e1-4963-bc16-4999b1eb1157
# ╠═fc29f9aa-b3b0-4fd2-b0bc-96301f950600
# ╟─d77e981c-9440-4882-abd5-ca391a054621
# ╠═b1830c1b-3f00-42f9-81bf-adecfa05c161
# ╠═03aac330-27fc-459e-8261-0fdac3721197
# ╟─af8072fd-3175-49c6-9248-41f48588d5a2
# ╠═72506456-8104-45b8-8243-722bd3bbce1f
# ╟─ace49f1a-2149-4d68-81ea-965d4d0e9197
# ╠═6c9e8943-0997-4f01-bc16-68bf62f72d21
# ╟─a544197a-a657-453c-80d5-b08e0d2dcdf6
# ╠═2e786f7a-e0a4-4e7c-90b6-7081864b6fa2
# ╟─6b91710e-5b70-434e-917c-d8c1e5f9809f
# ╠═3add75d7-88e6-4eb3-ad52-ace52755cf58
# ╟─71cd7876-82ca-4b09-b5b9-d79160b8061b
# ╠═4594cabf-46b0-42dc-b027-92c42de34d15
# ╟─d36ccfa0-6a44-4644-aa69-a46c09668c4d
# ╠═c0ab3d4b-93f1-4588-b7fa-c007039889a4
# ╟─06339cb3-9e6b-4214-9276-86574cc5d6bb
# ╠═bb0149cd-6e2d-4bf1-af68-71b3ddae9a61
# ╟─83177d5d-28a4-4d52-8df4-4b316e5e3080
# ╠═d09fbb3f-9cd8-4c8a-b2f8-2258579a71b5
# ╟─e5de2132-0a97-4802-a10a-cd0e770c8985
# ╠═77929858-3d6a-4832-ac1a-726ef652aea9
# ╟─bb14c31f-4d17-4fe6-86a0-7c2799c078f4
# ╠═085df56e-9938-45bb-944b-21173b3f6c76
# ╟─a515c3be-a27b-4cbf-9417-7b8446399ad2
# ╠═05f18b69-67ae-40fc-aefd-622af5f46a96
# ╟─c59e1737-d5fa-4cfe-9bd8-f0c8eb677231
# ╠═242397e3-e3e7-4f76-817c-2bccfb69009a
# ╟─79802f7a-50f2-4ae4-8991-046cd8fff9c6
# ╠═fd8a44d2-ca0b-469c-a33e-ecf02dec6955
# ╟─0e6ddaa6-3d8c-4db7-904a-e8eb76acb083
# ╠═ec5e8a05-43f4-4f8a-a6e0-cdc7213a2d82
# ╟─595844c9-cc40-4aaf-862a-e8becf9b4f0c
# ╠═f4e78715-446b-4434-b10e-a39ed3cac57c
# ╠═bc2e20b3-7f7d-418e-a5ac-0bb90bfcecbc
# ╟─f66596d8-a5c2-4025-995e-a0a1a3ebbade
# ╠═3881814e-0956-4588-831c-a6afc18c4290
# ╟─985f658f-81ad-46d9-8158-673a10c77628
# ╠═dd6cf48b-a756-4d1f-95a9-89649892bab7
# ╠═6cf155cd-07e5-466d-9e98-903450b969fc
# ╠═3426c47a-b63c-4d3c-89d3-24bbd751300c

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using Cairn

using Base.Test
import SimpleFileServer
using Logging
@Logging.configure(level=DEBUG)

# write your own tests here

workers = addprocs(10)
@info "Added 10 Procs"
mktempdir() do tmpd
    upload_name = "TEST_FILE"
    name, handle = mktemp(tmpd)
    test_data = Mmap.mmap(handle, Array{UInt8, 1}, 1000)
    rand!(test_data)
    Mmap.sync!(test_data)
    # Launch Consensus Nodes
    peers = [("localhost", x) for x = 8001:8003 ]
    @everywhere using Cairn
    @sync for (_, port ) = peers
        @spawnat pop!(workers) begin
            @async mktempdir() do f
                Cairn.Server.start(Cairn.Server.make(port, port - 8000, peers, f))
            end
        end
    end
    @info "Cairn has started"

    # Launch Some FS Nodes
    fs_nodes = [("localhost", x) for x = 8004:8010]
    @sync for (_, port) = fs_nodes
        @spawnat pop!(workers) let port = port
            cli = Cairn.Client.make("localhost", 8001)
            add_comm = Cairn.RPC.AddStorageNode("localhost", port)
            setup = () -> Cairn.Client.command(cli, add_comm)
            @async mktempdir() do f
                SimpleFileServer.start(port, f, setup)
            end

        end
    end
    @info "FS Nodes started"

    # Connect to the cluster
    c = Cairn.Client.make("localhost", 8001)
    while true
        try
            n = Cairn.Client.upload(c, upload_name, test_data, 3)
            @debug "$n Replicas have file"
            if n > 0
                break
            end
        catch err
            @debug err
        end
        wait(Timer(1))
    end
    tmp = ASCIIString(mktemp()[1])
    while true
        try
            d  = Cairn.Client.download(c, upload_name, tmp)
            @test d == test_data
            break
        catch err
            @debug err
        end
        wait(Timer(1))
    end

end


text:
using Flux
using Flux: mse, throttle
using Base.Iterators: repeated

# input
x = [0 0 1 1
     0 1 0 1
     1 1 1 1]
# output
y = [0 0 1 1]

# model
m = Dense(3,1,σ)    # only one layer (excl. input layer) with 3 inputs, 1 output and a sigmoid function

# loss function
loss(x,y) = mse(m(x), y)           # calculates the mean squared error

dataset = repeated((x, y), 10000)  # Train the model 10000 times with the same dataset
opt = SGD(params(m), 0.25)         # Gradient descent with learning rate 0.25
evalcb = () -> @show(loss(x, y))   # output the progress

Flux.train!(loss, dataset, opt, cb = throttle(evalcb, 0.5))

# output final error
@show(loss(x, y))

# testing the results
@show m([1,0,0])    # should be close to 1
@show m([0,1,0])    # should be close to 0


text:
using Documenter, Archimed

makedocs(;
    modules=[Archimed],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/VEZY/Archimed.jl/blob/{commit}{path}#L{line}",
    sitename="Archimed.jl",
    authors="remi.vezy <VEZY@users.noreply.github.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/VEZY/Archimed.jl",
)


text:
f = open("jawiki-country.json", "r")

for line in readlines(f)
  if ismatch(r"\"title\": \"イギリス\"}", line)
    line = replace(line,r"^\{\"text\":\s", "")
    line = replace(line,r",\s\"title\":\s\".+\"\}$", "")
    println(replace(line,"\\n", "\n"))
  end
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

function set_terminal(term::AbstractString)
    global gaston_config

    if validate_terminal(term)
        gaston_config.terminal = term
    else
        error(string("Terminal type ", term, " not supported."))
    end
    return term
end

# For terminals that support file output, set the file name
function set_filename(s::AbstractString)
    global gaston_config

    gaston_config.outputfile = s
end

# Set default values for CurveConf and AxesConf. Return set value upon
# success, error out otherwise.
function set_default_legend(s::AbstractString)
    global gaston_config

    gaston_config.legend = s
end

function set_default_plotstyle(s::AbstractString)
    global gaston_config

    if validate_2d_plotstyle(s) || validate_3d_plotstyle(s) ||
        validate_image_plotstyle(s)
        gaston_config.plotstyle = s
    else
        error(string("Plotstyle ", s, " not supported."))
    end
end

function set_default_color(s::AbstractString)
    global gaston_config

    # We can't validate colors yet
    gaston_config.color = s
end

function set_default_marker(s::AbstractString)
    global gaston_config

    if validate_marker(s)
        gaston_config.marker = s
    else
        error(string("Marker ", s, " not supported."))
    end
end

function set_default_linewidth(i::Real)
    global gaston_config

    gaston_config.linewidth = i
end

function set_default_pointsize(i::Real)
    global gaston_config

    gaston_config.pointsize = i
end

function set_default_title(s::AbstractString)
    global gaston_config

    gaston_config.title = s
end

function set_default_xlabel(s::AbstractString)
    global gaston_config

    gaston_config.xlabel = s
end

function set_default_ylabel(s::AbstractString)
    global gaston_config

    gaston_config.ylabel = s
end

function set_default_zlabel(s::AbstractString)
    global gaston_config

    gaston_config.zlabel = s
end

function set_default_fill(s::AbstractString)
	global gaston_config

	if validate_fillstyle(s)
		gaston_config.fill = s
	else
		error(string("Fill style ", s, " not supported."))
	end
end

function set_default_grid(s::AbstractString)
	global gaston_config

	if validate_grid(s)
		gaston_config.grid = s
	else
		error(string("Grid style ", s, " not supported."))
	end
end

function set_default_box(s::AbstractString)
    global gaston_config

    # We can't validate box yet
    gaston_config.box = s
end

function set_default_axis(s::AbstractString)
    global gaston_config

    if validate_axis(s)
        gaston_config.marker = s
    else
        error(string("Axis ", s, " not supported."))
    end
end

function set_default_xrange(s::AbstractString)
    global gaston_config

    if validate_range(s)
        gaston_config.xrange = s
    else
        error(string("Range ", s, " not supported."))
    end
end

function set_default_yrange(s::AbstractString)
    global gaston_config

    if validate_range(s)
        gaston_config.yrange = s
    else
        error(string("Range ", s, " not supported."))
    end
end

function set_default_zrange(s::AbstractString)
    global gaston_config

    if validate_range(s)
        gaston_config.zrange = s
    else
        error(string("Range ", s, " not supported."))
    end
end

# functions to set print parameters
function set_print_color(s::AbstractString)
    global gaston_config

    gaston_config.print_color = b
end

function set_print_fontface(s::AbstractString)
    global gaston_config

    gaston_config.print_fontface = s
end

function set_print_fontsize(s::Real)
    global gaston_config

    gaston_config.print_fontsize = s
end

function set_print_fontscale(r::Real)
    global gaston_config

    gaston_config.print_fontscale = r
end

function set_print_linewidth(r::Real)
    global gaston_config

    gaston_config.print_linewidth = r
end

function set_print_size(s::AbstractString)
    global gaston_config

    gaston_config.print_size = s
end


text:
module PublicSuffix

using Compat

export Domain, tld_exists, public_suffix, PublicSuffixList
# TODO: puny code support
#export puny_encode, puny_decode
#include("punycode.jl")

type DPart
    name::AbstractString
    is_leaf::Bool
    is_exception::Bool
    children::Dict{AbstractString,DPart}
end

type PublicSuffixList
    _dtree::DPart

    function PublicSuffixList(list_file::AbstractString="")
        if isempty(list_file)
            list_file = joinpath(dirname(@__FILE__), "..", "data", "effective_tld_names.dat")
        end

        root = DPart("", false, false, Dict{AbstractString,DPart}())
        open(list_file, "r") do f
            for line in readlines(f)
                (isempty(line) || @compat(startswith(line, "//")) || isspace(line[1])) && continue
                line = split(line)[1]
                add_node(line, root)
            end
        end
        new(root)
    end

    function add_node(fullname, root::DPart)
        is_exception = false
        if @compat(startswith(fullname, '!'))
            is_exception = true
            fullname = fullname[2:end]
        end

        # parse fullname, hook into _dtree
        parts = split(fullname, '.')
        nparts = length(parts)
        node = root
        for idx in 1:nparts
            idx = nparts+1-idx
            part = parts[idx]
            if haskey(node.children, part)
                node = node.children[part]
            elseif !is_exception && haskey(node.children, "*")
                node = node.children["*"]
            else
                new_node = DPart(part, false, false, Dict{AbstractString,DPart}())
                node.children[part] = new_node
                node = new_node
            end
        end
        node.is_exception = is_exception
        node.is_leaf = true
        nothing
    end
end

type Domain
    full::AbstractString
    sub_domain::AbstractString
    public_suffix::AbstractString
    top_domain::AbstractString
    typ::Type    # valid types: IPv6, IPv4, Domain

    function Domain(domain::AbstractString, list::PublicSuffixList=_def_list)
        try
            ip = parseip(domain)
            # this is an ip
            return new(domain, "", "", "", typeof(ip))
        end
        domain = lowercase(domain)
        (@compat(startswith(domain, '.')) || endswith(domain, '.')) && error("invalid domain name $domain")

        parts = split(domain, '.')
        nparts = length(parts)
        node = list._dtree
        pub_nodes = AbstractString[]
        sub_nodes = AbstractString[]
        top_domain = AbstractString[]
        tld_mode = true
        is_exception = false

        for idx in 1:nparts
            idx = nparts+1-idx
            part = parts[idx]
            isempty(part) && error("invalid domain name $domain")

            if tld_mode && haskey(node.children, part)
                node = node.children[part]
                if node.is_exception
                    insert!(sub_nodes, 1, part)
                    tld_mode = false
                    is_exception = true
                else
                    insert!(pub_nodes, 1, part)
                end
                isempty(top_domain) && insert!(top_domain, 1, part)
            elseif tld_mode && haskey(node.children, "*")
                node = node.children["*"]
                insert!(pub_nodes, 1, part)
                tld_mode = false
            else
                insert!(sub_nodes, 1, part)
            end
        end
        isempty(sub_nodes) && error("incomplete domain name $domain")
        if isempty(top_domain)
            if length(sub_nodes) > 1
                insert!(pub_nodes, 1, pop!(sub_nodes))
                insert!(pub_nodes, 1, pop!(sub_nodes))
            end
        else
            insert!(pub_nodes, 1, is_exception ? sub_nodes[end] : pop!(sub_nodes))
        end
        new(domain, join(sub_nodes, '.'), join(pub_nodes, '.'), isempty(top_domain) ? "" : top_domain[1], Domain)
    end
end

_def_list = PublicSuffixList()

function tld_exists(tld::AbstractString; list::PublicSuffixList=_def_list)
    parts = split(tld, '.')
    nparts = length(parts)
    node = list._dtree
    for idx in 1:nparts
        idx = nparts+1-idx
        part = parts[idx]
        if haskey(node.children, part)
            node = node.children[part]
        elseif haskey(node.children, "*")
            node = node.children["*"]
        else
            return false
        end
    end
    return node.is_leaf
end

public_suffix(domain::AbstractString; list::PublicSuffixList=_def_list) = return Domain(domain, list).public_suffix

end # module