text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test Base.Cartesian.exprresolve(:(1 + 3)) == 4
ex = Base.Cartesian.exprresolve(:(if 5 > 4; :x; else :y; end))
@test ex.args[2] == QuoteNode(:x)

@test Base.Cartesian.lreplace!("val_col", Base.Cartesian.LReplace{String}(:col, "col", 1)) == "val_1"
@test Base.setindex(CartesianIndex(1,5,4),3,2) == CartesianIndex(1, 3, 4)

@testset "CartesianIndices constructions" begin
    @testset "AbstractUnitRange" begin
        for oinds in [
            (2, 3),
            (UInt8(2), 3),
            (2, UInt8(3)),
            (2, 1:3),
            (Base.OneTo(2), 1:3)
        ]
            R = CartesianIndices(oinds)
            @test size(R) == (2, 3)
            @test axes(R) == (Base.OneTo(2), Base.OneTo(3))
            @test step.(R.indices) == (1, 1)
            @test step(R) == CartesianIndex(1, 1)

            @test R[begin] == CartesianIndex(1, 1)
            @test R[2] == CartesianIndex(2, 1)
            @test R[1, 2] == CartesianIndex(1, 2)
            @test R[end] == CartesianIndex(2, 3)
        end
        @test CartesianIndices((2, 3)) == CartesianIndex(1, 1):CartesianIndex(2, 3)

        R = CartesianIndices((0:5, 0:5))
        @test R[begin] == R[1] == first(R) == CartesianIndex(0, 0)
        @test R[2, 1] == R[2] == CartesianIndex(1, 0)
        @test R[1, 2] == R[7] == CartesianIndex(0, 1)
        @test R[end] == R[length(R)] == last(R) == CartesianIndex(5, 5)

        for oinds in [(2, ), (2, 3), (2, 3, 4)]
            R = CartesianIndices(oinds)
            @test eltype(R) == CartesianIndex{length(oinds)}
            @test ndims(R) == length(oinds)
            @test size(R) == oinds
        end

        # generic iterators doesn't have axes interface
        iter = Iterators.repeated([1 2], 4)
        @test_throws MethodError CartesianIndices(iter)
    end

    @testset "Step Range" begin
        for oinds in [
            (2, 1:2:6),
            (Base.OneTo(2), 1:2:6),
            (UInt8(2), 1:2:6),
            (2, UInt8(1):UInt8(2):UInt8(6))
        ]
            R = CartesianIndices(oinds)
            @test size(R) == (2, 3)
            @test axes(R) == (Base.OneTo(2), Base.OneTo(3))
            @test step.(R.indices) == (1, 2)
            @test step(R) == CartesianIndex(1, 2)

            @test R[begin] == CartesianIndex(1, 1)
            @test R[2] == CartesianIndex(2, 1)
            @test R[1, 2] == CartesianIndex(1, 3)
            @test R[end] == CartesianIndex(2, 5)
        end

        @test CartesianIndices((1:2:5, 1:3:7)) == CartesianIndex(1, 1):CartesianIndex(2,3):CartesianIndex(5,7)

        R = CartesianIndex(0, 0):CartesianIndex(2, 3):CartesianIndex(5, 7)
        @test R[begin] == R[1] == first(R) == CartesianIndex(0, 0)
        @test R[2, 1] == R[2] == CartesianIndex(2, 0)
        @test R[1, 2] == R[4] == CartesianIndex(0, 3)
        @test R[end] == R[length(R)] == last(R) == CartesianIndex(4, 6)

        for oinds in [(1:2:5, ), (1:2:5, 1:3:7), (1:2:5, 1:3:7, 1:4:11)]
            R = CartesianIndices(oinds)
            @test eltype(R) == CartesianIndex{length(oinds)}
            @test ndims(R) == length(oinds)
            @test size(R) == length.(oinds)
        end

        R = CartesianIndices((1:2:5, 7:-3:1))
        @test R == CartesianIndex(1, 7):CartesianIndex(2,-3):CartesianIndex(5, 1)
        @test step.(R.indices) == (2, -3)
        @test R[begin] == R[1] == first(R) == CartesianIndex(1, 7)
        @test R[2, 1] == R[2] == CartesianIndex(3, 7)
        @test R[1, 2] == R[4] == CartesianIndex(1, 4)
        @test R[end] == R[length(R)] == last(R) == CartesianIndex(5, 1)
    end

    @testset "IdentityUnitRange" begin
        function _collect(A)
            rst = eltype(A)[]
            for i in A
                push!(rst, i)
            end
            rst
        end
        function _simd_collect(A)
            rst = eltype(A)[]
            @simd for i in A
                push!(rst, i)
            end
            rst
        end

        for oinds in [
            (Base.IdentityUnitRange(0:1),),
            (Base.IdentityUnitRange(0:1), Base.IdentityUnitRange(0:2)),
            (Base.IdentityUnitRange(0:1), Base.OneTo(3)),
        ]
            R = CartesianIndices(oinds)
            @test axes(R) === oinds
            @test _collect(R) == _simd_collect(R) == vec(collect(R))
        end
        R = CartesianIndices((Base.IdentityUnitRange(0:1), 0:1))
        @test axes(R) == (Base.IdentityUnitRange(0:1), Base.OneTo(2))

    end

    for oinds in [(2, 3), (0:1, 0:2), (0:1:1, 0:1:2), (Base.IdentityUnitRange(0:1), Base.IdentityUnitRange(0:2)) ]
        R = CartesianIndices(oinds)
        @test vec(LinearIndices(R)) == 1:6
    end
    # TODO: non-1 steps are not supported yet, but may change in the future
    @test_throws ArgumentError LinearIndices(CartesianIndices((1:2:5, )))
    @test_throws ArgumentError LinearIndices(CartesianIndices((1:1:5, 1:2:5)))
end

module TestOffsetArray
    isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
    using .Main.OffsetArrays
    using Test

    A = OffsetArray(rand(2, 3), -1, -1)
    R = CartesianIndices(A)
    @test R == CartesianIndices((0:1, 0:2))
    @test axes(R) == (0:1, 0:2)
    for i in eachindex(A)
        @test A[i] == A[R[i]]
    end
    for i in R
        @test A[i] == A[Tuple(i)...]
    end
end

@testset "CartesianIndices getindex" begin
    @testset "AbstractUnitRange" begin
        for oinds in [(2, ), (2, 3), (2, 3, 4)]
            A = rand(1:10, oinds)
            R = CartesianIndices(A)
            @test R == CartesianIndices(oinds)

            @test A[R] == A
            @test axes(A) == axes(R)
            @test all(i->A[i]==A[R[i]], eachindex(A))
            @test all(i->A[i]==A[R[i]], R)
            @test all(i->A[i]==A[R[i]], collect(R))
            @test all(i->i in R, collect(R))
        end
    end

    @testset "StepRange" begin
        for oinds in [(1:2:5, ), (1:2:5, 1:3:7), (1:2:5, 1:3:7, 1:4:11)]
            A = rand(1:10, last.(oinds))
            R = CartesianIndices(A)

            SR = CartesianIndex(first.(oinds)):CartesianIndex(step.(oinds)):CartesianIndex(last.(oinds))
            @test A[oinds...] == A[SR]
            @test A[SR] == A[R[SR]]

            # TODO: A[SR] == A[Linearindices(SR)] should hold for StepRange CartesianIndices
            @test_broken A[SR] == A[LinearIndices(SR)]
        end
    end
end

@testset "range interface" begin
    for (I, i, i_next) in [
        (CartesianIndices((1:2:5, )), CartesianIndex(2, ), CartesianIndex(4, )),
        (1:2:5, 2, 4),
    ]
        # consistent with ranges behavior
        @test !(i in I)
        @test iterate(I, i) == (i_next, i_next)
    end

    # check iteration behavior on boundary
    R = CartesianIndex(1, 1):CartesianIndex(2, 3):CartesianIndex(4, 5)
    @test R.indices == (1:2:3, 1:3:4)
    i = CartesianIndex(4, 1)
    i_next = CartesianIndex(1, 4)
    @test !(i in R) && iterate(R, i) == (i_next, i_next)

    for R in [
        CartesianIndices((1:-1:-1, 1:2:5)),
        CartesianIndices((2, 3)),
        CartesianIndex(1, 2) .- CartesianIndices((1:-1:-1, 1:2:5)),
        CartesianIndex(1, 2) .- CartesianIndices((2, 3)),
    ]
        Rc = collect(R)
        @test all(map(==, R, Rc))
    end
end

@testset "Cartesian simd/broadcasting" begin
    @testset "AbstractUnitRange" begin
        A = rand(-5:5, 64, 64)
        @test abs.(A) == map(abs, A)

        function test_simd(f, @nospecialize(A); init=zero(eltype(A)))
            val_simd = init
            @simd for i in CartesianIndices(A)
                val_simd = f(val_simd, A[i])
            end

            val_iter = init
            for i in CartesianIndices(A)
                val_iter = f(val_iter, A[i])
            end

            @test val_iter == reduce(f, A, init=init)
            @test val_iter ≈ val_simd
        end

        test_simd(+, A)
    end

    R = CartesianIndex(-1, -1):CartesianIndex(6, 7)
    @test R .+ CartesianIndex(1, 2) == CartesianIndex(0, 1):CartesianIndex(7, 9)
    @test R .- CartesianIndex(1, 2) == CartesianIndex(-2, -3):CartesianIndex(5, 5)
    # 37867: collect is needed
    @test collect(CartesianIndex(1, 2) .- R) == CartesianIndex(2, 3):CartesianIndex(-1, -1):CartesianIndex(-5, -5)

    R = CartesianIndex(-1, -1):CartesianIndex(2, 3):CartesianIndex(6, 7)
    @test R .+ CartesianIndex(2, 2) == CartesianIndex(1, 1):CartesianIndex(2, 3):CartesianIndex(8, 9)
    @test R .- CartesianIndex(2, 2) == CartesianIndex(-3, -3):CartesianIndex(2, 3):CartesianIndex(4, 5)
    # 37867: collect is needed
    @test collect(CartesianIndex(1, 1) .- R) == CartesianIndex(2, 2):CartesianIndex(-2, -3):CartesianIndex(-4, -4)
end

@testset "Iterators" begin
    @testset "Reverse" begin
        R = CartesianIndices((0:5, 0:5))
        RR = Iterators.Reverse(R)
        rR = reverse(R)
        @test rR == collect(RR)
        @test rR.indices == (5:-1:0, 5:-1:0)

        @test eltype(RR) == CartesianIndex{2}
        @test size(RR) == size(R)
        @test axes(RR) == axes(R)

        @test first(RR) == last(R) == CartesianIndex(5, 5)
        @test last(RR) == first(R) == CartesianIndex(0, 0)
        RRR = collect(Iterators.Reverse(collect(RR)))
        @test R == RRR
    end

    @testset "collect" begin
        for oinds in [(0:5, ), (2:2:7, ), (2:-1:0, ),
                      (0:5, 2:8), (2:2:7, 3:3:10), (2:-1:0, 2:7),]
            R = CartesianIndices(oinds)
            @test collect(R) == R
        end
    end
end

@testset "set operations" begin
    R1 = CartesianIndices((3, 4, 5))
    R2 = CartesianIndices((-2:2, -3:3, -4:4))
    R = CartesianIndices((2, 3, 4))
    @test intersect(R1, R2) == R
end

# test conversions for CartesianIndex

@testset "CartesianIndex Conversions" begin
    @test convert(Int, CartesianIndex(42)) === 42
    @test convert(Float64, CartesianIndex(42)) === 42.0
    @test convert(Tuple, CartesianIndex(42, 1)) === (42, 1)

    # can't convert higher-dimensional indices to Int
    @test_throws MethodError convert(Int, CartesianIndex(42, 1))
end

@testset "CartesianIndices overflow" begin
    @testset "incremental steps" begin
        I = CartesianIndices((1:typemax(Int),))
        i = last(I)
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:2:typemax(Int), ))
        i = CartesianIndex(typemax(Int)-1)
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:(typemax(Int)-1),))
        i = CartesianIndex(typemax(Int))
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:2:typemax(Int)-1, ))
        i = CartesianIndex(typemax(Int)-1)
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:typemax(Int), 1:typemax(Int)))
        i = last(I)
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:2:typemax(Int), 1:2:typemax(Int)))
        i = CartesianIndex(typemax(Int)-1, typemax(Int)-1)
        @test iterate(I, i) === nothing

        I = CartesianIndices((1:typemax(Int), 1:typemax(Int)))
        i = CartesianIndex(typemax(Int), 1)
        @test iterate(I, i) === (CartesianIndex(1, 2), CartesianIndex(1,2))

        I = CartesianIndices((1:2:typemax(Int), 1:2:typemax(Int)))
        i = CartesianIndex(typemax(Int)-1, 1)
        @test iterate(I, i) === (CartesianIndex(1, 3), CartesianIndex(1, 3))

        I = CartesianIndices((typemin(Int):(typemin(Int)+3),))
        i = last(I)
        @test iterate(I, i) === nothing

        I = CartesianIndices(((typemin(Int):2:typemin(Int)+3), ))
        i = CartesianIndex(typemin(Int)+2)
        @test iterate(I, i) === nothing
    end

    @testset "decremental steps" begin
        I = Iterators.Reverse(CartesianIndices((typemin(Int):typemin(Int)+10, )))
        i = last(I)
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):2:typemin(Int)+10, )))
        i = last(I)
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):typemin(Int)+10, )))
        i = CartesianIndex(typemin(Int))
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):2:typemin(Int)+10, )))
        i = CartesianIndex(typemin(Int))
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):typemin(Int)+10, typemin(Int):typemin(Int)+10)))
        i = last(I)
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):2:typemin(Int)+10, typemin(Int):2:typemin(Int)+10)))
        i = CartesianIndex(typemin(Int), typemin(Int))
        @test iterate(I, i) === nothing

        I = Iterators.Reverse(CartesianIndices((typemin(Int):typemin(Int)+10, typemin(Int):typemin(Int)+10)))
        i = CartesianIndex(typemin(Int), typemin(Int)+1)
        @test iterate(I, i) === (CartesianIndex(typemin(Int)+10, typemin(Int)), CartesianIndex(typemin(Int)+10, typemin(Int)))

        I = Iterators.Reverse(CartesianIndices((typemin(Int):2:typemin(Int)+10, typemin(Int):2:typemin(Int)+10)))
        i = CartesianIndex(typemin(Int), typemin(Int)+2)
        @test iterate(I, i) === (CartesianIndex(typemin(Int)+10, typemin(Int)), CartesianIndex(typemin(Int)+10, typemin(Int)))

        I = CartesianIndices((typemax(Int):-1:typemax(Int)-10, ))
        i = last(I)
        @test iterate(I, i) === nothing

        I = CartesianIndices((typemax(Int):-2:typemax(Int)-10, ))
        i = last(I)
        @test iterate(I, i) === nothing
    end
end

@testset "CartesianIndices iteration" begin
    I = CartesianIndices((2:4, 0:1, 1:1, 3:5))
    indices = Vector{eltype(I)}()
    for i in I
        push!(indices, i)
    end
    @test length(I) == length(indices)
    @test vec(I) == indices

    empty!(indices)
    I = Iterators.reverse(I)
    for i in I
        push!(indices, i)
    end
    @test length(I) == length(indices)
    @test vec(collect(I)) == indices

    # test invalid state
    I = CartesianIndices((2:4, 3:5))
    @test iterate(I, CartesianIndex(typemax(Int), 3))[1] == CartesianIndex(2,4)
    @test iterate(I, CartesianIndex(typemax(Int), 4))[1] == CartesianIndex(2,5)
    @test iterate(I, CartesianIndex(typemax(Int), 5))    === nothing

    @test iterate(I, CartesianIndex(3, typemax(Int)))[1] == CartesianIndex(4,typemax(Int))
    @test iterate(I, CartesianIndex(4, typemax(Int)))    === nothing
end

@testset "CartesianIndices operations" begin
    I = CartesianIndices((1:3, 4:6))
    J = CartesianIndices((2:4, 3:5))

    @test @inferred(intersect(I, J)) == CartesianIndices((2:3, 4:5))
end

# issue #39705
f39705() = Base.Cartesian.@nany 0 _ -> true
@test f39705() === false

max_stars_repo_name: chenbobby/julia


text:
# functions related to beta distributions

import .RFunctions:
    betapdf,
    betalogpdf,
    betacdf,
    betaccdf,
    betalogcdf,
    betalogccdf,
    betainvcdf,
    betainvccdf,
    betainvlogcdf,
    betainvlogccdf

# pdf for numbers with generic types
betapdf(α::Real, β::Real, x::Number) = x^(α - 1) * (1 - x)^(β - 1) / beta(α, β)

# logpdf for numbers with generic types
betalogpdf(α::Real, β::Real, x::Number) = (α - 1) * log(x) + (β - 1) * log1p(-x) - lbeta(α, β)

ext: jl


text:
module ShirleyRayTracer

using StaticArrays
using LinearAlgebra

const Vec3 = SVector{3, Float64}
const Point3 = SVector{3, Float64}
const Color = SVector{3, Float64}
const Scanline = Vector{Color}

include("PPM.jl")

export Scene, Camera, Point3, Vec3, Color, Scanline
export trace_scanline, render
export magnitude, add!, randf

magnitude(x,y) = sqrt(x^2 + y^2)
magnitude(x,y,z) = sqrt(x^2 + y^2 + z^2)
magnitude(v) = magnitude(v...)

randf(fmin, fmax) = fmin + (fmax-fmin)*rand()
near_zero(v) = v.x < 1e-8 && v.y < 1e-8 && v.z < 1e-8

function random_in_unit_disk()
	x,y = randf(-1, 1), randf(-1, 1)
	while magnitude(x,y)^2 >= 1
     	x,y = randf(-1, 1), randf(-1, 1)
	end
	x,y
end

function random_in_unit_sphere()
	x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	while magnitude(x,y,z)^2 >= 1
		x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	end
	Point3(x,y,z)
end

random_unit_vector() = normalize(random_in_unit_sphere())

function random_in_hemisphere(normal)
    in_unit_sphere = random_in_unit_sphere()
    dot(in_unit_sphere, normal) > 0.0 ? in_unit_sphere : -in_unit_sphere
end

function refract(uv, n, etai_over_etat)
    cos_theta = min(dot(-uv, n), 1.0)
    r_out_perp = etai_over_etat * (uv + cos_theta*n)
    r_out_parallel = -sqrt(abs(1.0 - magnitude(r_out_perp)^2)) * n
    r_out_perp + r_out_parallel
end

struct Ray
	origin::Point3
	direction::Vec3
	udirection::Vec3
	tm::Float64
	Ray(o, d, m) = new(o, d, normalize(d), m)
	Ray(o, d) = Ray(o, d, 0)
	Ray() = new(Vec3(Inf, Inf, Inf), Vec3(Inf, Inf, Inf)) # a type stable sentinel instead of using nothing
end

at(r::Ray, t) = r.origin + t * r.direction

struct Camera
	origin::Point3
	lower_left_corner::Point3
	horizontal::Vec3
	vertical::Vec3
	u::Vec3
	v::Vec3
	w::Vec3
	lens_radius::Float64
	time0::Float64
	time1::Float64
	function Camera(lookfrom, lookat, vup, vfov, aspect_ratio, aperture, focus_dist, _time0=0, _time1=0)
		viewport_height = 2.0 * tan(deg2rad(vfov)/2)
		viewport_width = aspect_ratio * viewport_height

		w = normalize(lookfrom - lookat)
		u = normalize(cross(vup, w))
		v = cross(w, u)

		origin = lookfrom
		horizontal = focus_dist * viewport_width * u
		vertical = focus_dist * viewport_height * v
		lower_left_corner = origin - horizontal/2 - vertical/2 - focus_dist*w

		lens_radius = aperture / 2
		time0 = _time0
		time1 = _time1
		new(origin, lower_left_corner, horizontal, vertical, u, v, w, lens_radius, time0, time1)
	end
	Camera() = Camera(Point3(0,0,-1), Point3(0,0,0), Vec3(0,1,0), 40, 1, 0, 10)
end

abstract type Material end
abstract type Hitable end

include("Hitables.jl")
include("Materials.jl")

struct Scene
	camera::Camera
	hitables::Vector{Hitable}
	Scene(cam) = new(cam, Vector{Hitable}())
end

add!(s::Scene, h::Hitable) = push!(s.hitables, h)

get_ray(scene::Scene, s, t) = get_ray(scene.camera, s, t)

function get_ray(cam::Camera, s, t)
	x,y = cam.lens_radius .* random_in_unit_disk()
	offset = cam.u * x + cam.v *y
	Ray(cam.origin + offset, cam.lower_left_corner + s*cam.horizontal + t*cam.vertical - cam.origin - offset, randf(cam.time0, cam.time1))
end

function trace(scene::Scene, ray::Ray, t_min, t_max)
	closest_t = t_max
	local struck::Hitable

	for hitable in scene.hitables
		t = trace(hitable, ray, t_min, closest_t)
		if t >= 0
			closest_t = t
			struck = hitable
		end
	end

	if @isdefined struck
		Hit(struck, closest_t, ray)
	else
		Hit()
	end
end

function ray_color(scene::Scene, ray::Ray, depth)::Tuple{Float64, Float64, Float64}
	if depth <= 0
        	return 0,0,0
	end

	hit = trace(scene, ray, 0.001, Inf)
	if isdefined(hit, :material)
		s, a = scatter(hit.material, ray, hit)
		if s.origin.x == Inf
			return 0,0,0
		end
		r,g,b = ray_color(scene, s, depth-1)
		return a[1] * r, a[2] * g, a[3] * b
	end

	t::Float64 = 0.5*(ray.udirection.y + 1.0)
	t1m = 1.0 - t
	t1m + 0.5t, t1m + 0.7t, t1m + t
end

function trace_scanline(world, y, nsamples, width, height, max_depth)
	scanline = Scanline(undef, width)
	@simd for x in 1:width
		r=g=b=0.0
		for _ in 1:nsamples
			(r,g,b) = (r,g,b) .+ ray_color(world, get_ray(world, (x + rand()) / width, (y + rand()) / height), max_depth)
		end
		@inbounds scanline[x] = Color(r/nsamples, g/nsamples, b/nsamples)
	end
	scanline
end

function render(scene::Scene, width, height, nsamples=10, max_depth=50)
	scanlines = Vector{Scanline}(undef, height)
	Threads.@threads for y in 1:height
		@inbounds scanlines[y] = trace_scanline(scene, y, nsamples, width, height, max_depth)
    end
	scanlines
end

###
end


text:
## Exercise 4-6
## Write a function called circle that takes a turtle, t, and radius, r, as parameters and that draws an approximate circle by calling polygon with an appropriate length and number of sides. Test your function with a range of values of r.

## TIP
## Figure out the circumference of the circle and make sure that len * n == circumference.
using ThinkJulia

println("Ans: ")

function polygon(turtle::Turtle, nsides::Int, len::Int = 100)
    angle = (360 / nsides)

    for i in 1:nsides
        forward(turtle, len)
        turn(turtle, -angle)
    end
end

function circle(turtle::Turtle, radius::Int)
    circumference = 2 * π * radius
    n = Int(trunc(circumference ÷ 3) + 3)
    len = Int(trunc(circumference / n))
    @show n, len
    polygon(turtle, n, len)
end

@svg begin
    turtle = Turtle()
    circle(turtle, 100)
end

println("End.")


text:
function main(r::Robot)
    steps = 1
    side = Nord
    while ismarker(r)==false
        for _ in 1:2
            special_move!(r,side,steps)
            side=counterclockwise_side(side)
        end
        steps += 1
    end
    println("Marker found")
end
function special_move!(r::Robot,side::HorizonSide,num::Int)
    for _ in 1:num
        if ismarker(r)
            return nothing
        end
        isborder(r, side) ? go_around(r, side) : move!(r,side)
    end
end

function go_around(r, direct_side)
    move_side = clockwise_side(direct_side)
    left_moves = 0
    right_moves = 0
    while isborder(r, direct_side)
        left_moves += 1
        right_moves = left_moves
        while isborder(r,direct_side) && right_moves > 0
            right_moves -= 1
            move!(r,move_side)
        end
        side=inverse(move_side)
    end
    move!(r,direct_side)
    left_moves /= 2
    while left_moves > 0
        left_moves -= 1
        move!(r,move_side)
    end
end

counterclockwise_side(side::HorizonSide) = HorizonSide(mod(Int(side)+1,4))
clockwise_side(side::HorizonSide) = HorizonSide(mod(Int(side)-1,4))
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
"""
    AbstractSet

Abstract supertype for all set types.
"""
abstract type AbstractSet end


text:
__precompile__(true)
module JuML
export importcsv,
       Seq,
       DataFrame,
       Factor,
       Covariate,
       BoolVariate,
       factor,
       covariate,
       getname,
       getlevels,
       isordinal,
       xgblogit,
       cvxgblogit,
       predict,
       getstats,
       getauc,
       getlogloss,
       CovariateStats,
       FactorStats,
       LevelStats

include("const.jl")
include("DataStructures/seq.jl")
include("DataStructures/list.jl")
include("DataStructures/tree.jl")
include("util.jl")
include("DataImport/dataimport.jl")
include("DataFrame/dataframe.jl")

include("Factors/factor.jl")
include("Factors/constfactor.jl")
include("Factors/filefactor.jl")
include("Factors/ordinalfactor.jl")
include("Factors/cachedfactor.jl")
include("Factors/bincovfactor.jl")
include("Factors/boolvarfactor.jl")
include("Factors/maplevelfactor.jl")
include("Factors/permutefactor.jl")
include("Factors/widerfactor.jl")

include("Covariates/covariate.jl")
include("Covariates/constcovariate.jl")
include("Covariates/filecovariate.jl")
include("Covariates/cachedcovariate.jl")
include("Covariates/transcovariate.jl")
include("Covariates/trans2covariate.jl")
include("Covariates/parsefactorcovariate.jl")
include("Covariates/ifelsecovariate.jl")

include("BoolVariates/boolvariate.jl")
include("BoolVariates/transboolvariate.jl")
include("BoolVariates/trans2boolvariate.jl")
include("BoolVariates/transcovboolvariate.jl")
include("BoolVariates/trans2covboolvariate.jl")

include("XGBoost/tree.jl")
include("XGBoost/split.jl")
include("XGBoost/logistic.jl")

end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using DifferentialEquations
using ModelingToolkit
using OrdinaryDiffEq
using StochasticDiffEq
using DiffEqJump
using Random
using Plots

@parameters t β c γ
@variables S(t) I(t) R(t)

N=S+I+R # This is recognized as a derived variable
rxs = [Reaction((β*c)/N, [S,I], [I], [1,1], [2])
       Reaction(γ, [I], [R])]

rs  = ReactionSystem(rxs, t, [S,I,R], [β,c,γ])

tmax = 40.0
tspan = (0.0,tmax);

u0 = [S => 990.0,
      I => 10.0,
      R => 0.0];

p = [β=>0.05,
     c=>10.0,
     γ=>0.25];

Random.seed!(1234);

odesys = convert(ODESystem, rs)
oprob = ODEProblem(odesys, u0, tspan, p)
osol = solve(oprob, Tsit5())
plot(osol)

sdesys = convert(SDESystem, rs)
sprob = SDEProblem(sdesys, u0, tspan, p)
ssol = solve(sprob, LambaEM())
plot(ssol)

jumpsys = convert(JumpSystem, rs)
u0i = [S => 990, I => 10, R => 0]
dprob = DiscreteProblem(jumpsys, u0i, tspan, p)
jprob = JumpProblem(jumpsys, dprob, Direct())
jsol = solve(jprob, SSAStepper())
plot(jsol)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: snowdj/sir-julia

max_issues_repo_name: snowdj/sir-julia

max_forks_repo_name: snowdj/sir-julia


text:
using UnivariateFunctions
using StochasticIntegrals
using DataFrames, Dates
using LinearAlgebra
using Distributions: Normal, MersenneTwister, quantile
using Statistics: var, mean, cov
using Sobol
using Distributions
using Random
const tol = 10*eps()

USD_vol = 0.0001
GBP_vol = 0.00001
GBP_FX_Vol = 0.002

today = Date(2016,12,1)
tommorow = Date(2016,12,2)
date_2020 = Date(2020,12,1)
later_date = Date(2035,12,1)
later_later_date = Date(2080,12,1)

USD_hw_a_curve = PE_Function(USD_vol)
GBP_hw_a_curve = PE_Function(GBP_vol)

brownian_corr_matrix = Hermitian([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
brownian_ids = [:USD_IR, :GBP_IR, :GBP_FX, :BARC]
USD_IR_a_ito  = ItoIntegral(:USD_IR, USD_hw_a_curve)
GBP_IR_a_ito  = ItoIntegral(:GBP_IR, GBP_hw_a_curve)
GBP_FX_ito    = ItoIntegral(:GBP_FX, GBP_FX_Vol)
ito_integrals = Dict([:USD_IR_a, :GBP_IR_a, :GBP_FX] .=> [USD_IR_a_ito, GBP_IR_a_ito, GBP_FX_ito])

ito_set_ = ItoSet(brownian_corr_matrix, brownian_ids, ito_integrals)

for_covar = ForwardCovariance(ito_set_, years_from_global_base(today),
                              years_from_global_base(date_2020))
simp_covar = SimpleCovariance(ito_set_, years_from_global_base(today),
                              years_from_global_base(date_2020))

# Making ItoProcesses

ito_processes = Dict{Symbol,ItoProcess{Float64}}()
update_rates = Dict{Symbol,Exponential}()
for k in keys(ito_integrals)
    ito_processes[k] = ItoProcess( 0.0, 100.0, PE_Function(0.0), ito_integrals[k] )
    update_rates[k] = Exponential(5.0)
end

# Making the series.

for_series_sync = make_ito_process_syncronous_time_series(deepcopy(ito_processes), deepcopy(for_covar), 5.0, 1000)
abs(nrow(for_series_sync) - 3000) < 0.5
simp_series_sync = make_ito_process_syncronous_time_series(deepcopy(ito_processes), deepcopy(simp_covar), 5.0, 1000)
abs(nrow(simp_series_sync) - 3000) < 0.5

for_series_nonsync  = make_ito_process_non_syncronous_time_series(deepcopy(ito_processes), deepcopy(for_covar), update_rates , 1000)
abs(nrow(for_series_nonsync) - 1000) < 0.5
simp_series_nonsync = make_ito_process_non_syncronous_time_series(deepcopy(ito_processes), deepcopy(simp_covar), update_rates , 1000)
abs(nrow(simp_series_nonsync) - 1000) < 0.5

# Testing rolling covariance by date.
rolled_date = ForwardCovariance(for_covar, today, tommorow)
rolled_numbers = ForwardCovariance(for_covar, years_from_global_base(today), years_from_global_base(tommorow))
all(abs.(rolled_date.covariance_ .- rolled_numbers.covariance_) .< 100*eps())


text:
struct RealManifold{N, T <: AbstractInterpolation} <: AbstractManifold{N}
    num_azimuth_angles::Int
    num_elevation_angles::Int
    azimuth_step::Float64
    elevation_step::Float64
    max_elevation::Float64
    expansion_length::Int
    lut::T
end

"""
LUT must be a vector with the length of the number of antenna elements of type
SVector of Matrices.
Each matrix must have the size of #elevation angles x #azimuth angles. The
azimuth angle starts from the x-Axis and increases towards the y-Axis. The
zenith must be at the top of the matrix and the horizont / bottom must be at
the bottom of the matrix. The maximal elevation is given by `max_elevation`.
A value of `π/2` means horizont and `π` until bottom of antenna.
The type of interpolation can be set by the optional parameter `interpolation`.
Possible values are: `Constant`, `Linear` or `Quadratic`. For more information
see regarding the interpolation, see:
https://github.com/JuliaMath/Interpolations.jl.
"""
function RealManifold(
    lut::SVector{N, <: AbstractMatrix{Complex{T}}},
    interpolation::Type{Q} = Constant;
    max_elevation = 1π,
    normalize = true
) where {N, Q <: Interpolations.Degree, T <: Real}
    lut_array = [lut[i][j,k] for i = 1:length(lut), j = 1:size(lut[1],1), k = 1:size(lut[1],2)]
    RealManifold(
        lut_array,
        NumAnts(N),
        Q,
        max_elevation = max_elevation,
        normalize = normalize
    )
end

function RealManifold(
    lut::AbstractArray{Complex{T}, 3},
    num_ants::NumAnts{N},
    interpolation::Type{Q} = Constant;
    max_elevation = 1π,
    normalize = true
) where {N, Q <: Interpolations.Degree, T <: Real}
    test_lut_correctness(lut, max_elevation)
    num_elevation_angles = size(lut, 2)
    num_azimuth_angles = size(lut, 3)
    elevation_step = max_elevation / (num_elevation_angles - 1)
    azimuth_step = 2 * π / num_azimuth_angles
    normalized_lut = normalize ? norm_manifold(lut) : lut
    expansion_length = calc_expansion_length(interpolation)
    expanded_lut = expand(normalized_lut, expansion_length)
    interpolated_lut = interpolate(expanded_lut, (NoInterp(), BSpline(get_elevation_interpolation(Q)), BSpline(get_azimuth_interpolation(Q))))
    RealManifold{N, typeof(interpolated_lut)}(num_azimuth_angles, num_elevation_angles, azimuth_step, elevation_step, max_elevation, expansion_length, interpolated_lut)
end

get_elevation_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC{0} = Q()
get_elevation_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC{1} = Q()
# In Elevation direction the boundary is reflective with a 180° circleshift, see `expand`
get_elevation_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC = Q(Reflect(OnCell()))
get_azimuth_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC{0} = Q()
get_azimuth_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC{1} = Q()
# In Azimuth direction the boundary is periodic so one could use Quadratic(Periodic(OnCell())), but
# this is twice as slow as Quadratic(Reflect(OnCell())), that's why we copy left to the right and right to left
get_azimuth_interpolation(::Type{Q}) where Q <: Interpolations.DegreeBC = Q(Reflect(OnCell()))

function RealManifold(lut::AbstractMatrix{Complex{T}}...; interpolation::Type{Q} = Constant, max_elevation = 1π, normalize = true) where {Q <: Interpolations.Degree, T <: Real}
    RealManifold(SVector(lut), Q, max_elevation = max_elevation, normalize = normalize)
end

function get_steer_vec(manifold::RealManifold{N}, doa::Spherical) where N
    elevation = π / 2 - doa.ϕ # convert to mathematic
    # Rotate azimuth by 180° if elevation is less than zero or greater than π
    azimuth = doa.θ + π * ((elevation < 0) || (elevation > π))
    # Convey the elevation in between 0 <= ϕ <= π
    elevation -= (elevation < 0) * (2 * elevation) +
        (elevation > π && manifold.max_elevation ≈ π) * (2 * (elevation - π))
    elevation > manifold.max_elevation && error("Elevation is above max elevation. Try to increase max elevation.")
    elevation_index = elevation / manifold.elevation_step + 1 + manifold.expansion_length
    azimuth_index = mod(azimuth / manifold.azimuth_step, manifold.num_azimuth_angles) + 1 + manifold.expansion_length
    get_steer_vec(manifold, azimuth_index, elevation_index)
end

function get_steer_vec(manifold::RealManifold{N}, azimuth_index::Real, elevation_index::Real) where N
    steer_vec = MVector{N, ComplexF64}(undef)
    @inbounds for i = 1:length(steer_vec)
        steer_vec[i] = manifold.lut(i, elevation_index, azimuth_index)
    end
    SVector(steer_vec)
end

# See this issue: https://github.com/JuliaGeometry/CoordinateTransformations.jl/issues/25
# if you'd like to know the direction of θ and ϕ
function get_steer_vec(manifold::RealManifold, doa)
    doa_sph = cart2sph(doa) # θ can be between -π and π
    elevation = π / 2 - doa_sph.ϕ # convert to mathematic
    elevation > manifold.max_elevation && error("Elevation is above max elevation. Try to increase max elevation.")
    azimuth = doa_sph.θ + 2π * (doa_sph.θ < 0)
    elevation_index = elevation / manifold.elevation_step + 1 + manifold.expansion_length
    azimuth_index = azimuth / manifold.azimuth_step + 1 + manifold.expansion_length
    get_steer_vec(manifold, azimuth_index, elevation_index)
end

function calc_expansion_length(::Type{<:Interpolations.Degree})
    0
end

function calc_expansion_length(::Type{Quadratic})
    18
end

# Cubic has greater errors than Quadratic when compared with BSpline(Quadratic(Periodic(OnCell())))? Force error
function calc_expansion_length(::Type{Cubic})
    error("Could not find sufficient expansion length. Maybe implement boundary to be periodic?") # 1
end

function expand(lut::AbstractArray{Complex{T}, 3}, num_expand::Int) where {T <: Real}
    num_expand >= 0 || error("Expand number must be greater than zero")
    num_ants = size(lut, 1)
    num_ϕs = size(lut, 2)
    num_θs = size(lut, 3)
    lut_expanded = Array{Complex{T}, 3}(undef, num_ants, num_ϕs + 2 * num_expand, num_θs + 2 * num_expand)
    lut_expanded[:,num_expand + 1:num_ϕs + num_expand,num_expand + 1:num_θs + num_expand] .= lut
    # In Elevation direction the boundary is reflective with a 180° circleshift
    for ant = 1:num_ants
        for i = 1:num_expand
            lut_expanded[ant,num_expand - i + 1,num_expand + 1:num_θs + num_expand] .= circshift(lut_expanded[ant,num_expand + i + 1,num_expand + 1:num_θs + num_expand], floor(Int, num_θs / 2))
            lut_expanded[ant,num_ϕs + num_expand + i,num_expand + 1:num_θs + num_expand] .= circshift(lut_expanded[ant,num_ϕs + num_expand - i,num_expand + 1:num_θs + num_expand], floor(Int, num_θs / 2))
        end
        # In Azimuth direction the boundary is periodic so one could use Quadratic(Periodic(OnCell())), but
        # this is twice as slow as Quadratic(Reflect(OnCell())), that's why we copy left to the right and right to left
        for i = 1:num_expand
            lut_expanded[ant,:,num_θs + num_expand + i] .= lut_expanded[ant,:,num_expand + i]
            lut_expanded[ant,:,num_expand - i + 1] .= lut_expanded[ant,:,num_θs + num_expand - i + 1]
        end
    end
    lut_expanded
end

"""
  norm_manifold(lut::Array{Complex{Float64}, 3})

Normalizes the manifold such that the maximal norm is `sqrt(num_ants)`.
"""
function norm_manifold(lut::AbstractArray{Complex{T}, 3}) where {T <: Real}
    max_norm = mapreduce(norm, max, Slices(lut, True(), False(), False()))
    lut ./ (max_norm / sqrt(size(lut, 1)))
end

function test_lut_correctness(lut::AbstractArray{Complex{T}, 3}, max_el) where T <: Real
    num_els = size(lut, 2)
    num_ants = size(lut, 1)
    near_horizont_index = min(floor(Int, π / 2 / max_el * num_els), num_els)
    lut_reference_one = lut[2:end, :, :] ./ lut[1:1, :, :]
    zenith_angle_var = mean(var(angle.(@view(lut_reference_one[:, 1, :])), dims = 2))
    horizont_angle_var = mean(var(angle.(@view(lut_reference_one[:, near_horizont_index, :])), dims = 2))
    (zenith_angle_var > horizont_angle_var) && error("First row of manifold LUT should be zenith")
    error_az = norm(lut[:,:,1] - lut[:,:,end]) / num_els / num_ants
    (error_az > 2) && error("Azimuth of manifold LUT should vary over columns. First and last column should be similar.")
end


text:
# Goal is to extend the Julia API to be forgiving of mental typos using other language functions

print = println
len = length
type = typeof
gettype = typeof

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
using Test
using WikstromTerelius: PrimeGenerator, validate, order, Enc, Dec, modulus, value, gen_shuffle

q = 11
p = 2*q + 1

@test validate(PrimeGenerator(3, p)) == true
@test validate(PrimeGenerator(11, p)) == false

n = let
    n = 0
    for i in 1:p
        validate(PrimeGenerator(i, p)) && (n+=1)
    end
    n
end
@test n == q - 1

g = PrimeGenerator(3, p)

#q = 17
#g = PrimeGenerator(3, q)

@test g*g^2 == g^3
@test (g^2)^2 == g^4
@test g^(order(g) + 1) == g

h = g^7

@test h*h^2 == h^3
@test (h^2)^2 == h^4
@test h^(order(h) + 1) == h

@test inv(g)*g^2 == g
@test (g^7)^6 == g^(7*6) # This is only true for a cyclic group
@test g*g*g == g^3 # Checking multiplication
@test g^2/g == g

### Let's test ElGammal encryption

sk = 5
pk = g^sk
r = 3
m = g^5
r2 = 4

enc = Enc(pk, g)
dec = Dec(sk)

@test dec(enc(m, r)) == m
@test enc(enc(m, r), r2) == enc(m, r + r2)

### Shuffle generation

sk = 5
pk = g^sk

enc = Enc(pk, g)

m_vec = [g, g^2, g^3]
e_vec = enc.(m_vec, 1) # It is not necessary to randomize encryption for user. It however does make sense to do so for intermidiatery who collects messages from users to not reveal internals.

### The shuffling
r_vec = Int[1, 2, 3]

e_enc = enc.(e_vec, r_vec)
ψ = sortperm(e_enc)
sort!(e_enc)

@test sort(dec.(e_enc)) == sort(m_vec)

m_vec = [g, g^2, g^3]
e_vec = enc.(m_vec, 1)


text:
# Autogenerated wrapper script for HarfBuzz_jll for armv7l-linux-musleabihf
export libharfbuzz, libharfbuzz_gobject, libharfbuzz_subset

using Cairo_jll
using Fontconfig_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using Libffi_jll
JLLWrappers.@generate_wrapper_header("HarfBuzz")
JLLWrappers.@declare_library_product(libharfbuzz, "libharfbuzz.so.0")
JLLWrappers.@declare_library_product(libharfbuzz_gobject, "libharfbuzz-gobject.so.0")
JLLWrappers.@declare_library_product(libharfbuzz_subset, "libharfbuzz-subset.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Fontconfig_jll, FreeType2_jll, Glib_jll, Graphite2_jll, Libffi_jll)
    JLLWrappers.@init_library_product(
        libharfbuzz,
        "lib/libharfbuzz.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_gobject,
        "lib/libharfbuzz-gobject.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_subset,
        "lib/libharfbuzz-subset.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#-----------------------------------------------------------------------# sweepmatrix
"""
`sweepmatrix(x, y, addint = true)`

`sweepmatrix(x, y, wts, addint = true)`

Create the matrix [x y]' * [x y] for sweeping, optionally with weights
"""
function sweepmatrix{T <: LinAlg.BlasFloat}(x::AMat{T}, y::AVec{T}, addint::Bool = true)
    n, p = size(x)
    d = p + 1 + addint
    A = zeros(d, d)

    rng = (1:p) + addint
    # operations: x'x, x'y, y'y
    BLAS.syrk!('U', 'T', 1 / n, x, 0.0, view(A, rng, rng))
    BLAS.gemv!('T', 1 / n, x, y, 0.0, view(A, rng, d))
    A[end, end] = dot(y, y) / n

    if addint
        #operations: 1'1, 1'x, 1'y,
        A[1, 1] = 1.0
        A[1, 2:end - 1] = mean(x, 1)
        A[1, end] = mean(y)
    end
    A
end
function sweepmatrix{T <: LinAlg.BlasFloat}(x::AMat{T}, y::AVec{T}, wts::AVec{T}, addint::Bool = true)
    n, p = size(x)
    d = p + 1 + addint
    wts = wts / sum(wts)
    A = zeros(d, d)
    W = Diagonal(sqrt(wts))
    Wx = W * x
    Wy = W * y

    rng = (1:p) + addint
    # operations: x'Wx, x'Wy, y'Wy
    BLAS.syrk!('U', 'T', 1 / n, Wx, 0.0, view(A, rng, rng))
    BLAS.gemv!('T', 1 / n, Wx, Wy, 0.0, view(A, rng, d))
    A[end, end] = dot(y, Wy) / n

    if addint
        #operations: 1'Wx, 1'Wy,
        A[1, 2:end - 1] = mean(Wx, 1)
        A[1, end] = mean(y, StatsBase.WeightVec(wts))
    end
    A
end

"Add a ridge penalty to the sweepmatrix"
function addridge!{T <: LinAlg.BlasFloat}(A::AMat{T}, λ::T, intercept::Bool)
    rng = (1 + intercept):(size(A, 1) - 1)
    for i in rng
        A[i, i] *= λ
    end
    A
end

# Unpenalized or Ridge-regularized regression
function sweepreg{T <: LinAlg.BlasFloat}(x::AMat{T}, y::AVec{T}, addint::Bool = true; λ::T = zero(T))
    A = sweepmatrix(x, y, addint)
    if λ > zero(T)
        addridge!(A, λ, addint)
    end
    SweepOperator.sweep!(A, 1:size(A, 1) - 1)
    β = A[1:end-1, end]
    SSE = A[end, end]
    minusXTXinv = A[1:end-1, 1:end-1]
    β, SSE, minusXTXinv
end

max_stars_repo_name: joshday/SparseRegressionAlgorithms

max_issues_repo_name: joshday/SparseRegressionAlgorithms

max_forks_repo_name: joshday/SparseRegressionAlgorithms