text:
using Colors
using Javis

function ground(args...)
    background("black")
    sethue("white")
    return nothing
end

function title_object(ptitle)
    fontsize(20)
    text("Little Prince 44 sunsets in a day", ptitle, valign = :middle, halign = :center)
    return nothing
end

function intro_text_object(pt_text, width_text)
    fontsize(20)
    #=     textwrap("Here comes some introductory story about Little Prince"
            * "and his little planet, and blah, blah, blah", width_text,
            pt_text) =#
    textbox(
        [
            "Here is the story",
            "of this little fella",
            "known as the Little Prince",
            "and the day he saw",
            "44 sunsets,",
            "on a single day,",
            "walking around",
            "his little planet.",
        ],
        pt_text,
    )
    return nothing
end

function sun_object(c, r, d, θ)
    sun_color = RGB(1.0, 0.9, 0.2)
    sethue(sun_color)
    circle(c + d * Point(cos(θ), sin(θ)), r, :fill)
    return nothing
end

function planet_object(c, r, θ)
    planet_shadow_color = RGB(0.01, 0.04, 0.12)
    planet_sunny_color = RGB(0.4, 0.8, 1.0)
    sethue(planet_shadow_color)
    # circle(c, r, :fill)
    for ζ = 0:0.02:1
        α = π * ζ / 2
        p = 2.5
        sethue((1 - ζ^p) * planet_sunny_color + ζ^p * planet_shadow_color)
        P1 = c + r * Point(cos(θ + α), sin(θ + α))
        P2 = c + r * Point(cos(θ - α), sin(θ - α))
        arc2r(c, P1, P2, :fill)
    end
    stem_color = 0 < θ < π ? sethue("darkgreen") : sethue("lightgreen")
    flower_color = 0 < θ < π ? sethue("darkred") : sethue("red")
    sethue(stem_color)
    rect(c - r * Point(0, 1.0) - Point(r / 40.0, r / 10.0), r / 20.0, r / 5.0, :fill)
    sethue(flower_color)
    circle(c - r * Point(0.0, 1.0) - Point(0.0, r / 6.0), r / 12, :fill)
    return nothing
end

function make_scene()
    nframes = 400

    width, height = 854, 480 # approximatey 16:9 aspect ratio
    origin = Point(div(height, 2) - div(width, 2), 0)
    ptitle = origin - Point(0, 0.9 * div(height, 2))

    pt_intro_text = origin + Point(div(height, 2), -div(height, 5))
    width_intro_text = div(height, 2)

    vid = Video(width, height)

    fadein = 20
    fadeout = 10

    # initial frame, last frame and number of frames
    title_n0, title_n1 = 11, 100; title_nf = title_n1 - title_n0 + 1
    intro_n0, intro_n1 = 121, 300; intro_nf = intro_n1 - intro_n0 + 1

    Background(1:nframes, ground)

    title = Object(title_n0:title_n1, (args...) -> title_object(ptitle))
    act!(
        title,
        [
            Action(1:fadein, appear(:fade)),
            Action(title_nf-fadeout:title_nf, disappear(:fade)),
        ],
    )

    intro_text = Object(
        intro_n0:intro_n1,
        (args...) -> intro_text_object(pt_intro_text, width_intro_text),
    )
    act!(
        intro_text,
        [
            Action(1:fadein, appear(:fade)),
            Action(intro_nf-fadeout:intro_nf, disappear(:fade)),
        ],
    )

    sun = Object(1:nframes, (args...; θ) -> sun_object(origin, 20, 220, θ))
    act!(sun, Action(1:nframes, change(:θ, 0 => 2π)))

    planet = Object(1:nframes, (args...; θ) -> planet_object(origin, 80, θ))
    act!(planet, Action(1:nframes, change(:θ, 0 => 2π)))

    render(vid; pathname = joinpath("..", "images", "44sunsets.gif"), framerate = 15)
end

make_scene()


text:
module Statistikk

using ConjugatePriors
using Distributions

import Base: rand
import ConjugatePriors: NormalInverseChisq
import Distributions: invsqrt2π, log2π, logpdf, pdf, zval
import Distributions: MixtureModel, Normal
import InteractiveUtils: @which

export UnnormalizedCategorical

include("conjugate.jl")
include("distributions.jl")
include("missings.jl")
include("vector.jl")
include("dict.jl")

end

ext: jl

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk


text:
struct Manifold{T}
    penetration::T
    normal::SA.SVector{2, T}
    contact::SA.SVector{2, T}
end

get_penetration(manifold::Manifold) = manifold.penetration
get_normal(manifold::Manifold) = manifold.normal
get_tangent(manifold::Manifold) = rotate_minus_90(get_normal(manifold))
get_contact(manifold::Manifold) = manifold.contact

#####
# StdCircle vs. StdCircle
#####

function Manifold(a::StdCircle{T}, b::StdCircle{T}, pos_ba::SA.SVector{2, T}) where {T}
    r_a = get_radius(a)
    r_b = get_radius(b)
    normal = LA.normalize(pos_ba)
    if all(isfinite.(normal))
        penetration = r_a + r_b - LA.norm(pos_ba)
        contact = (r_a - penetration / 2) * normal
        return Manifold(penetration, normal, contact)
    else
        penetration = r_a + r_b
        normal = SA.SVector(one(T), zero(T))
        contact = (r_a - penetration / 2) * normal
        return Manifold(penetration, normal, contact)
    end
end

Manifold(a::StdCircle{T}, b::StdCircle{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T} = Manifold(a, b, pos_ba)

#####
# HyperRectangle vs. HyperSphere
#####

function get_closest_edge_from_inside(rect::StdRect{T}, pos::SA.SVector{2, T}) where {T}
    half_width = get_half_width(rect)
    half_height = get_half_height(rect)

    x = get_x(pos)
    y = get_y(pos)

    penetration, edge_id = findmin((y + half_height, half_width - x, half_height - y, x + half_width))

    return (penetration, edge_id)
end

function Manifold(a::StdRect{T}, b::StdCircle{T}, pos_ba::SA.SVector{2, T}) where {T}
    r_b = get_radius(b)
    projection = get_projection(a, pos_ba)
    vec = pos_ba - projection
    normal = LA.normalize(vec)
    if all(isfinite.(normal))
        penetration = r_b - LA.norm(vec)
        contact = pos_ba - (r_b - penetration / 2) * normal
        return Manifold(penetration, normal, contact)
    else
        d, edge_id = get_closest_edge_from_inside(a, pos_ba)
        penetration = r_b + d
        normal = get_normals(a)[edge_id]
        contact = pos_ba - (r_b - penetration / 2) * normal
        return Manifold(penetration, normal, contact)
    end
end

function Manifold(a::StdCircle{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}) where {T}
    pos_ab = -pos_ba
    manifold_ab = Manifold(b, a, pos_ab)

    penetration = get_penetration(manifold_ab)
    normal = manifold_ab |> get_normal |> rotate_180
    contact = pos_ba + get_contact(manifold_ab)
    return Manifold(penetration, normal, contact)
end

Manifold(a::StdRect{T}, b::StdCircle{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T} = Manifold(a, b, pos_ba)

function Manifold(a::StdCircle{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T}
    dir_ab = invert_relative_direction(dir_ba)
    pos_ab = -rotate(pos_ba, dir_ab)
    manifold_ab = Manifold(b, a, pos_ab, dir_ab)

    penetration = get_penetration(manifold_ab)
    normal = get_relative_direction(dir_ab, rotate_180(get_normal(manifold_ab)))
    contact = rotate(get_contact(manifold_ab) - pos_ab, dir_ba)
    return Manifold(penetration, normal, contact)
end

#####
# StdRect vs. StdRect
#####

function Manifold(a::StdRect{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}) where {T}
    half_width_a = get_half_width(a)
    half_height_a = get_half_height(a)

    bottom_left = max.(get_bottom_left(a), get_bottom_left(b, pos_ba))
    top_right = min.(get_top_right(a), get_top_right(b, pos_ba))

    contact = (top_right + bottom_left) / 2
    x_contact = get_x(contact)
    y_contact = get_y(contact)

    half_widths_intersection = (top_right - bottom_left) / 2
    half_width_intersection = half_widths_intersection[1]
    half_height_intersection = half_widths_intersection[2]

    penetration, edge_id = findmin((half_height_a + y_contact + half_height_intersection,
                                    half_width_a - x_contact + half_width_intersection,
                                    half_height_a - y_contact + half_height_intersection,
                                    half_width_a + x_contact + half_width_intersection,
                                   ))

    normal = get_normals(a)[edge_id]

    return Manifold(penetration, normal, contact)
end

function get_clipped_vertices(a::StdRect{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T}
    half_width_a = get_half_width(a)
    half_height_a = get_half_height(a)

    vertices_ba = get_vertices(b, pos_ba, dir_ba)
    VecType = typeof(pos_ba)

    initial_vertices = (vertices_ba..., vertices_ba[1])
    final_vertices = VecType[]

    # e1
    for i in 1:length(initial_vertices) - 1
        v1 = initial_vertices[i]
        v2 = initial_vertices[i + 1]

        x1 = v1[1]
        y1 = v1[2]

        x2 = v2[1]
        y2 = v2[2]

        # in-in
        if (y1 >= -half_height_a) && (y2 >= -half_height_a)
            push!(final_vertices, v2)
        # in-out
        elseif (y1 >= -half_height_a) && (y2 < -half_height_a)
            x_intersection = x1 + (x2 - x1) * (y1 + half_height_a) / (y1 - y2)
            if isfinite(x_intersection)
                push!(final_vertices, VecType(x_intersection, -half_height_a))
            else
                push!(final_vertices, (v1 + v2) / 2)
            end
        # out-in
        elseif (y1 < -half_height_a) && (y2 >= -half_height_a)
            x_intersection = x1 + (x2 - x1) * (y1 + half_height_a) / (y1 - y2)
            if isfinite(x_intersection)
                push!(final_vertices, VecType(x_intersection, -half_height_a))
                push!(final_vertices, v2)
            else
                push!(final_vertices, (v1 + v2) / 2)
                push!(final_vertices, v2)
            end
        end
    end

    # e2
    initial_vertices = (final_vertices..., final_vertices[1])
    final_vertices = VecType[]

    for i in 1:length(initial_vertices) - 1
        v1 = initial_vertices[i]
        v2 = initial_vertices[i + 1]

        x1 = v1[1]
        y1 = v1[2]

        x2 = v2[1]
        y2 = v2[2]

        # in-in
        if (x1 <= half_width_a) && (x2 <= half_width_a)
            push!(final_vertices, v2)
        # in-out
        elseif (x1 <= half_width_a) && (x2 > half_width_a)
            y_intersection = y1 + (y2 - y1) * (half_width_a - x1) / (x2 - x1)
            if isfinite(y_intersection)
                push!(final_vertices, VecType(half_width_a, y_intersection))
            else
                push!(final_vertices, (v1 + v2) / 2)
            end
        # out-in
        elseif (x1 > half_width_a) && (x2 <= half_width_a)
            y_intersection = y1 + (y2 - y1) * (half_width_a - x1) / (x2 - x1)
            if isfinite(y_intersection)
                push!(final_vertices, VecType(half_width_a, y_intersection))
                push!(final_vertices, v2)
            else
                push!(final_vertices, (v1 + v2) / 2)
                push!(final_vertices, v2)
            end
        end
    end

    initial_vertices = (final_vertices..., final_vertices[1])
    final_vertices = VecType[]

    # e3
    for i in 1:length(initial_vertices) - 1
        v1 = initial_vertices[i]
        v2 = initial_vertices[i + 1]

        x1 = v1[1]
        y1 = v1[2]

        x2 = v2[1]
        y2 = v2[2]

        # in-in
        if (y1 <= half_height_a) && (y2 <= half_height_a)
            push!(final_vertices, v2)
        # in-out
        elseif (y1 <= half_height_a) && (y2 > half_height_a)
            x_intersection = x1 + (x2 - x1) * (half_height_a - y1) / (y2 - y1)
            if isfinite(x_intersection)
                push!(final_vertices, VecType(x_intersection, half_height_a))
            else
                push!(final_vertices, (v1 + v2) / 2)
            end
        # out-in
        elseif (y1 > half_height_a) && (y2 <= half_height_a)
            x_intersection = x1 + (x2 - x1) * (half_height_a - y1) / (y2 - y1)
            if isfinite(x_intersection)
                push!(final_vertices, VecType(x_intersection, half_height_a))
                push!(final_vertices, v2)
            else
                push!(final_vertices, (v1 + v2) / 2)
                push!(final_vertices, v2)
            end
        end
    end

    # e4
    initial_vertices = (final_vertices..., final_vertices[1])
    final_vertices = VecType[]

    for i in 1:length(initial_vertices) - 1
        v1 = initial_vertices[i]
        v2 = initial_vertices[i + 1]

        x1 = v1[1]
        y1 = v1[2]

        x2 = v2[1]
        y2 = v2[2]

        # in-in
        if (x1 >= -half_width_a) && (x2 >= -half_width_a)
            push!(final_vertices, v2)
        # in-out
        elseif (x1 >= -half_width_a) && (x2 < -half_width_a)
            y_intersection = y1 + (y2 - y1) * (x1 + half_width_a) / (x1 - x2)
            if isfinite(y_intersection)
                push!(final_vertices, VecType(-half_width_a, y_intersection))
            else
                push!(final_vertices, (v1 + v2) / 2)
            end
        # out-in
        elseif (x1 < -half_width_a) && (x2 >= -half_width_a)
            y_intersection = y1 + (y2 - y1) * (x1 + half_width_a) / (x1 - x2)
            if isfinite(y_intersection)
                push!(final_vertices, VecType(-half_width_a, y_intersection))
                push!(final_vertices, v2)
            else
                push!(final_vertices, (v1 + v2) / 2)
                push!(final_vertices, v2)
            end
        end
    end

    return final_vertices
end

function get_centroid(vertices::Vararg{SA.SVector{2, T}}) where {T}
    VecType = typeof(vertices[1])
    area = zero(T)
    c_x_num = zero(T)
    c_y_num = zero(T)
    den = zero(T)

    augmented_vertices = (vertices..., vertices[1])

    for i in 1:length(augmented_vertices) - 1
        v_i = augmented_vertices[i]
        v_i_plus_1 = augmented_vertices[i + 1]

        x_i = v_i[1]
        y_i = v_i[2]

        x_i_plus_1 = v_i_plus_1[1]
        y_i_plus_1 = v_i_plus_1[2]

        cross = x_i * y_i_plus_1 - x_i_plus_1 * y_i
        c_x_num += (x_i + x_i_plus_1) * cross
        c_y_num += (y_i + y_i_plus_1) * cross
        den += cross
    end

    c_x = c_x_num / (3 * den)
    c_y = c_y_num / (3 * den)
    if isfinite(c_x) && isfinite(c_y)
        return VecType(c_x, c_y)
    else
        centroid = zero(VecType)
        for vertex in vertices
            centroid = centroid + vertex
        end
        return centroid / length(vertices)
    end
end

function get_contact(a::StdRect{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T}
    clipped_vertices = get_clipped_vertices(a, b, pos_ba, dir_ba)
    return get_centroid(clipped_vertices...)
end

function get_candidate_support(a::StdRect{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T}
    half_width_a = get_half_width(a)
    half_height_a = get_half_height(a)

    vertices_ba = get_vertices(b, pos_ba, dir_ba)

    value_1, _, vertex_id_1 = findmax(vertex -> vertex[2], vertices_ba)
    max_penetration_1 = half_height_a + value_1

    value_2, _, vertex_id_2 = findmin(vertex -> vertex[1], vertices_ba)
    max_penetration_2 = half_width_a - value_2

    value_3, _, vertex_id_3 = findmin(vertex -> vertex[2], vertices_ba)
    max_penetration_3 = half_height_a - value_3

    value_4, _, vertex_id_4 = findmax(vertex -> vertex[1], vertices_ba)
    max_penetration_4 = half_width_a + value_4

    penetration, (_, vertex_id), edge_id = findmin(x -> x[1], ((max_penetration_1, vertex_id_1),
                                          (max_penetration_2, vertex_id_2),
                                          (max_penetration_3, vertex_id_3),
                                          (max_penetration_4, vertex_id_4)))

    return penetration, vertex_id, edge_id
end

function Manifold(a::StdRect{T}, b::StdRect{T}, pos_ba::SA.SVector{2, T}, dir_ba::SA.SVector{2, T}) where {T}
    penetration_ba, vertex_id_b, edge_id_a = get_candidate_support(a, b, pos_ba, dir_ba)

    pos_ab, dir_ab = invert(pos_ba, dir_ba)
    penetration_ab, vertex_id_a, edge_id_b = get_candidate_support(b, a, pos_ab, dir_ab)

    contact = get_contact(a, b, pos_ba, dir_ba)

    if penetration_ba <= penetration_ab
        penetration = penetration_ba
        normal = get_normals(a)[edge_id_a]
        return Manifold(penetration, normal, contact)
    else
        penetration = penetration_ab
        normal = get_normals(b, dir_ba)[edge_id_b]
        return Manifold(penetration, normal, contact)
    end
end


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

#This file contains the functions hecke_operator_on_basis and hecke_operator_on_qexp.

"""
    hecke_operator_on_basis(B, n, k)

Compute the matrix of the Hecke operator T_`n` of weight `k`
relative to the given basis `B` of q-expansions for a space of
modular forms.

# Arguments
- `B::Array`: array of q-expansions
- `n::Integer`: integer >=1
- `k::Integer`: weight

# Examples still missing
"""
function hecke_operator_on_basis(B::Array{fmpz_rel_series,1}, n::Int, k::Int)

   #error handling
   if n < 1
      error("n (=$n) must be a positive integer")
   end
   if isa(B, Array) == false
      error("B (=$B) must be an array")
   end

   #construct dxd matrix
   ring = base_ring(B[1])
   d = length(B)			#check if this is always correct
   S = MatrixSpace(ring, d, d)
   matrix = S()

   #compute Tn(f) to precision d+1 for each element f of B
   for j in 1:d
      f = B[j]
      T_f = hecke_operator_on_qexp(f, n, k, d+1)
      #check if cusp form
      if d == dim_Sk(k)
         for i in 1:d
            #Tn for the jth element of B corresponds to the jth row
            matrix[j,i] = coeff(T_f,i)
         end
      else
         for i in 0:d-1
            #Tn for the jth element of B corresponds to the jth row
            matrix[j,i+1] = coeff(T_f,i)
         end
      end
   end

   return matrix

end

"""
    divisors(n)

Return an array consisting of all divisors of `n` if n is
nonzero, raise an error otherwise.
"""
function divisors(n::Int)

   #error handling
   if n == 0
      error("n must be nonzero")
   end

   n = abs(n)
   array = [div for div in 1:n if n%div == 0]
   return array

end

"""
    hecke_operator_on_qexp(f, n, k, prec=nothing)

Compute the image of the q-expansion `f` of a modular form under
the Hecke operator T_`n` of weight `k`. Return a power series to
precision `prec`.

# Arguments
- `f::RelSeriesElem`: q-expansion
- `n::Integer`: integer >=1
- `k::Integer`: weight
- `prec::Integer=nothing`: precision of the output

# Examples still missing
"""
function hecke_operator_on_qexp(f::fmpz_rel_series, n::Int, k::Int, prec::Union{Int, Nothing}=nothing)

   max_prec = Int(ceil(f.prec / Int(n)))	#check if works also for n nonprime
   if prec == nothing
      prec = max_prec
   elseif prec > max_prec
      error("desired precision is too high given precision of f")
   end

   R, q = PowerSeriesRing(base_ring(f), prec, "q")
   T_p = R(0)

   l = k-1
   for m in 0:prec-1		#start with 0 to deal with all modular forms
      sum = 0
      array = divisors(gcd(n,m)) 	#check if there is a quicker method
      for i in 1:length(array)
         d = array[i]
         if (m*n)%(d*d) == 0
            index = (m*n)//(d*d)
            sum += d^l * coeff(f, Int(index))
         end
      end
      T_p += sum*q^m
   end

   return T_p

end


text:
module Vizinanigans

using GLMakie, Statistics, Printf

include("utils.jl")
include("./Cartesian/volumeslice.jl")
include("./Cartesian/visualize.jl")

end # module


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using FactCheck

function pooling_forward(input, mask, output, kernel, stride, pad, _type)
  width, height, channels, num = size(input)
  pooled_width, pooled_height, _, _ = size(output)
  kernel_size = kernel * kernel

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*stride - pad + 1
          wstart = (pw-1)*stride - pad + 1
          hend = min(hstart + kernel - 1, height)
          wend = min(wstart + kernel - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          region = sub(input, wstart:wend, hstart:hend, c, n)
          if _type == :max
              index = indmax(region)
              mask[pw, ph, c, n] = index # note this is local index in region
              output[pw, ph, c, n] = region[index]
          elseif _type == :mean
              output[pw, ph, c, n] = sum(region) / kernel_size
          end
        end
      end
    end
  end
end

function convolution_forward(filter::Array, bias::Array, input::Array,
                             output::Array, stride, pad, kernel)
    width, height, channel, num = size(input)
    width_out, height_out, n_filter, _ = size(output)

    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out,
        k = 1:channel, p = 1:kernel, q = 1:kernel

        in_y = (y-1) * stride - pad + p
        in_x = (x-1) * stride - pad + q
        if (in_y >= 1 && in_y <= height && in_x >= 1 && in_x <= width)
            filter_idx = ((k - 1) * kernel + (p - 1)) * kernel + q
            output[x, y, o, n] += input[in_x, in_y, k, n] *
                                  filter[filter_idx, o]
        end
    end

    # add bias
    for n = 1:num, o = 1:n_filter, y = 1:height_out, x = 1:width_out
        output[x, y, o, n] += bias[o]
    end

    return output
end

facts("Testing Conv-Relu-Pool fusion") do
    net = Net(8)
    data,  data_value  = MemoryDataLayer(net, :data, (224, 224, 3))
    label, label_value = MemoryDataLayer(net, :label, (1,))
    data_value[:]  = rand(Float32, size(data_value)...) * 10
    label_value[:] = map(floor, rand(Float32, size(label_value)...) * 10)
    conv1        = ConvolutionLayer(:conv1, net, data, 10, 3, 1, 1)
    relu1        = ReLULayer(:relu1, net, conv1)
    pool1        = MaxPoolingLayer(:pool1, net, relu1, 2, 2, 0)

    conv2_1      = ConvolutionLayer(:conv2_1, net, pool1, 10, 3, 1, 1)
    relu2_1      = ReLULayer(:relu2_1, net, conv2_1)
    conv2_2      = ConvolutionLayer(:conv2_2, net, relu2_1, 10, 3, 1, 1)
    relu2_2      = ReLULayer(:relu2_2, net, conv2_2)
    pool2        = MaxPoolingLayer(:pool2, net, relu2_2, 2, 2, 0)

    init(net)

    input     = get_buffer(net, :datavalue)
    filters1  = get_buffer(net, :conv1weights)
    bias1     = get_buffer(net, :conv1bias)
    rand!(bias1)
    conv1_out = zeros(get_buffer(net, :conv1value))

    mask1     = get_buffer(net, :pool1maxidx)
    mask1_out = zeros(mask1)
    pool1_out = zeros(get_buffer(net, :pool1value))

    filters2_1  = get_buffer(net, :conv2_1weights)
    bias2_1     = get_buffer(net, :conv2_1bias)
    rand!(bias2_1)
    conv2_1_out = zeros(get_buffer(net, :conv2_1value))

    filters2_2  = get_buffer(net, :conv2_2weights)
    bias2_2     = get_buffer(net, :conv2_2bias)
    rand!(bias2_2)
    conv2_2_out = zeros(get_buffer(net, :conv2_2value))

    mask2     = get_buffer(net, :pool2maxidx)
    mask2_out = zeros(mask2)
    pool2_out = zeros(get_buffer(net, :pool2value))
    forward(net; phase=Latte.Test)
    convolution_forward(filters1, bias1, input, conv1_out, 1, 1, 3)
    map!((x) -> x > 0.0f0 ? x : 0.0f0, conv1_out)
    @fact conv1_out --> roughly(get_buffer(net, :conv1value))
    pooling_forward(conv1_out, mask1_out, pool1_out, 2, 2, 0, :max)
    @fact pool1_out --> roughly(get_buffer(net, :pool1value))

    # Sometimes wrong because the c version chooses the wrong index on very
    # similar numbers, fast math comparison?
    @pending mask1_out --> mask1

    convolution_forward(filters2_1, bias2_1, pool1_out, conv2_1_out, 1, 1, 3)
    map!((x) -> x > 0.0f0 ? x : 0.0f0, conv2_1_out)
    @fact conv2_1_out --> roughly(get_buffer(net, :conv2_1value))
    convolution_forward(filters2_2, bias2_2, conv2_1_out, conv2_2_out, 1, 1, 3)
    map!((x) -> x > 0.0f0 ? x : 0.0f0, conv2_2_out)
    @fact conv2_2_out --> roughly(get_buffer(net, :conv2_2value))
    pooling_forward(conv2_2_out, mask2_out, pool2_out, 2, 2, 0, :max)
    @fact pool2_out --> roughly(get_buffer(net, :pool2value))

    # Sometimes wrong because the c version chooses the wrong index on very
    # similar numbers, fast math comparison?
    @pending mask2_out --> roughly(mask2)

    #=
    for i in CartesianRange(size(mask1_out))
        if mask1_out[i] != mask1[i]
            println("Error at at index $(i)")
            println("    mask1_out = $(mask1_out[i])")
            println("    mask1     = $(mask1[i])")
            pw, ph, c, n = i.I
            stride = 2
            pad = 0
            kernel = 2
            width, height, channels, batch_size = size(conv1_out)
            hstart = (ph-1)*stride - pad + 1
            wstart = (pw-1)*stride - pad + 1
            hend = min(hstart + kernel - 1, height)
            wend = min(wstart + kernel - 1, width)
            hstart = max(1, hstart)
            wstart = max(1, wstart)
            region = conv1_out[wstart:wend, hstart:hend, c, n]
            println((wstart:wend, hstart:hend, c, n))
            println(region)
            println(indmax(region))
            exit(1)
        end
    end
    for i in CartesianRange(size(mask2_out))
        if mask2_out[i] != mask2[i]
            println("Error at at index $(i)")
            println("    mask2_out = $(mask2_out[i])")
            println("    mask2     = $(mask2[i])")
            width, height, channels, batch_size = size(conv2_2_out)
            pw, ph, c, n = i.I
            stride = 2
            pad = 0
            kernel = 2
            hstart = (ph-1)*stride - pad + 1
            wstart = (pw-1)*stride - pad + 1
            hend = min(hstart + kernel - 1, height)
            wend = min(wstart + kernel - 1, width)
            hstart = max(1, hstart)
            wstart = max(1, wstart)
            println(size(conv2_2_out))
            region = conv2_2_out[wstart:wend, hstart:hend, c, n]
            println((wstart:wend, hstart:hend, c, n))
            println(region)
            println(indmax(region))
            println(indmax(sub(conv2_2_out, wstart:wend, hstart:hend, c, n)))
            exit(1)
        end
    end
    =#
end


text:
# Parameter Learning

module ParameterLearning

using SPNetworks
import SPNetworks:
    Node, SumNode, ProductNode, LeafNode, CategoricalDistribution, IndicatorFunction, GaussianDistribution,
    isleaf, isprod, issum,
    logpdf!, plogpdf!,
    vardims, layers
import Random

"""
    Parameter Learning Algorithm
"""
abstract type ParameterLearner end

"""
    Computes the partial derivatives of the network output w.r.t. to each node at the current evaluation. Assumes network has been evaluated first in original domain.

    Returns a vector of derivatives.
"""
# function backpropagaten!(spn::SumProductNetwork)::Vector{Float64}
#     # Assumes network has been evaluted at some assignment
#     # Backpropagate derivatives
#     diff = zeros(Float64,length(spn._nodes))
#     diff[1] = 1.0
#     for (i,node) in enumerate(spn._nodes)
#         if isa(node, SumNode)
#             for j in spn._backward[i]
#                 diff[j] += spn._weights[j,i]*diff[i]
#             end
#         elseif isa(node, ProductNode)
#             for j in spn._backward[i]
#                 diff[j] += diff[i]*spn._values[i]/spn._values[j]
#             end
#         end
#     end
#     return diff
# end

# """
#  Computes derivatives for values propagated in log domain.
# Returns vector of derivatives.
# """
# function backpropagate(spn::SumProductNetwork)::Vector{Float64}
#     # Assumes network has been evaluted at some assignment using logpdf!
#     # Backpropagate derivatives
#     diff = zeros(Float64,length(spn))
#     values = similar(diff)
#     logpdf(spn,values)
#     backpropagate(spn,values,diff)
#     # diff[1] = 1.0
#     # for i = 1:length(spn)
#     #     if isa(spn[i], SumNode)
#     #         for j in children(spn,i)
#     #             diff[j] += spn._weights[j,i]*diff[i]
#     #         end
#     #     elseif isa(spn[i], ProductNode)
#     #         for j in children(spn,i)
#     #             diff[j] += diff[i]*exp(spn._values[i]-spn._values[j])
#     #         end
#     #     end
#     # end
#     return diff
# end

"""
    backpropagate!(diff,spn,values)

 Computes derivatives for values propagated in log domain and stores results in given vector diff.
"""
function backpropagate!(diff::Vector{Float64},spn::SumProductNetwork,values::Vector{Float64})
    # Assumes network has been evaluted at some assignment using logpdf!
    # Backpropagate derivatives
    @assert length(diff) == length(spn) == length(values)
    fill!(diff,0.0)
    @inbounds diff[1] = 1.0
    for i = 1:length(spn)
        @inbounds node = spn[i]
        if issum(node)
            @inbounds for (k,j) in enumerate(node.children)
               diff[j] += node.weights[k]*diff[i]
            end
        elseif isprod(node)
            @inbounds for j in node.children
                if isfinite(values[j])
                    # @assert isfinite(exp(values[i]-values[j]))  "contribution to derivative of ($i,$j) is not finite: $(values[i]), $(values[j]), $(exp(values[i]-values[j]))"
                    diff[j] += diff[i]*exp(values[i]-values[j])
                else
                    δ = exp(sum(values[k] for k in node.children if k ≠ j))
                    # @assert isfinite(δ)  "contribution to derivative of ($i,$j) is not finite: $(values[i]), $(values[j]), $(δ)"
                    diff[j] += diff[i]*δ
                end
            end
        end
    end
end

# compute log derivatives (not working!)
function logbackpropagate(spn::SumProductNetwork,values::Vector{Float64},diff::Vector{Float64})
    # Assumes network has been evaluted at some assignment using logpdf!
    # Backpropagate derivatives
    @assert length(diff) == length(spn) == length(values)
    # create storage for each computed value (message from parent to child)
    from = []
    to = []
    for i = 1:length(spn)
        if isa(spn[i], SumNode) || isa(spn[i], ProductNode)
            for j in children(spn,i)
                push!(from,i)
                push!(to,j)
            end
        end
    end
    cache = sparse(to,from,ones(Float64,length(to)))
    #fill!(diff,0.0)
    logdiff = zeros(Float64,length(diff))
    for i = 1:length(spn)
        if i == 1
            diff[i] == 1.0
        else
            cache_vals = nonzeros(cache[i,:]) # incoming arc values
            offset = maximum(cache_vals)
            logdiff[i] = offset + log(sum(exp.(cache_vals.-offset)))
            diff[i] = isfinite(logdiff[i]) ? exp(logdiff[i]) : 0.0
        end
        if isa(spn[i], SumNode)
            for j in children(spn,i)
                #@inbounds diff[j] += getweight(spn,i,j)*diff[i]
                cache[j,i] = logweight(spn,i,j) + logdiff[i]
            end
        elseif isa(spn[i], ProductNode)
            for j in children(spn,i)
                #@inbounds diff[j] += diff[i]*exp(values[i]-values[j])
                cache[j,i] = logdiff[i] + values[i]-values[j]
            end
        end
    end
end

"""
 Computes derivatives for given vector of values propagated in log domain.

Returns vector of derivatives.
"""
function backpropagate(spn::SumProductNetwork,values::Vector{Float64})::Vector{Float64}
    diff = Array{Float64}(undef,length(spn))
    backpropagate!(diff,spn,values)
    return diff
end

"""
Random initialization of weights
"""
function initialize(learner::ParameterLearner) #spn::SumProductNetwork)
    spn = learner.spn
    sumnodes = filter(i -> isa(spn[i], SumNode), 1:length(spn))
    for i in sumnodes
        #@inbounds ch = spn[i].children  # children(spn,i)
        @inbounds Random.rand!(spn[i].weights)
        @inbounds spn[i].weights ./= sum(spn[i].weights)
        @assert sum(spn[i].weights) ≈ 1.0 "Unnormalized weight vector at node $i: $(sum(spn[i].weights)) | $(spn[i].weights)"
    end
    gaussiannodes = filter(i -> isa(spn[i],GaussianDistribution), 1:length(spn))
    for i in gaussiannodes
        @inbounds spn[i].mean = rand()
        @inbounds spn[i].variance = 1.0
    end
end

include("ExpectationMaximization.jl")
include("SquaredExpectationMaximization.jl")

end


text:
@testset "Objective terms" begin
    local_opt = ECOS.Optimizer # Support for SOCP.

    g = MetaDiGraph()
    add_vertex!(g, :name, "A")
    add_vertex!(g, :name, "B")
    add_vertex!(g, :name, "C")
    add_vertex!(g, :name, "D")
    add_edge!(g, 1, 2, :capacity, 2.0)
    add_edge!(g, 1, 3, :capacity, 2.0)
    add_edge!(g, 2, 4, :capacity, 2.0)
    add_edge!(g, 3, 4, :capacity, 2.0)
    add_edge!(g, 4, 1, :capacity, 0.001) # inedge(src), outedge(dst)

    k = SimpleDiGraph(4)
    add_edge!(k, 1, 4)
    d = Edge(1, 4)
    dm = Dict{Edge{Int}, Float64}(d => 4.0)

    mt = ModelType(Load(), MinimumTotal(), FlowFormulation(), false, CuttingPlane(), ObliviousUncertainty(), UncertainDemand())
    rd = RoutingData(g, k, local_opt, mt)
    m = Model(local_opt)
    set_silent(m)
    rm = basic_routing_model_unitary(m, rd, FlowFormulation())
    dm = Dict{Edge{Int}, Float64}(d => 4.0)

    @testset "Nonexistent function" begin
        @test_throws ErrorException objective_edge_expression(rm, TestEdgeWise(), Edge(1, 2))
    end

    @testset "Using a local traffic matrix" begin
        rd_dm = RoutingData(g, k, local_opt, mt, traffic_matrix=dm)
        m_dm = Model(local_opt)
        set_silent(m_dm)
        rm_dm = basic_routing_model_unitary(m_dm, rd_dm, FlowFormulation())
        @test objective_edge_expression(rm_dm, Load(), Edge(1, 2), dm) == objective_edge_expression(rm_dm, Load(), Edge(1, 2))
    end

    @testset "Load" begin
        @test_throws ErrorException objective_edge_expression(rm, Load(), Edge(1, 2))
        @test objective_edge_expression(rm, Load(), Edge(1, 2), dm) == total_flow_in_edge(rm, Edge(1, 2), dm) / 2.0
        @test objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 2 * rm.routing[d, Edge(1, 2)]
    end

    @testset "Kleinrock" begin
        # Needs a SOCP solver.
        cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 0.75)
        obj = objective_edge_expression(rm, KleinrockLoad(), Edge(1, 2), dm)
        @objective(rm.model, Min, obj)
        optimize!(rm.model)
        @test termination_status(rm.model) == MOI.OPTIMAL
        @test value(obj) ≈ 0.75 / 0.25
        delete(rm.model, cref)
    end

    @testset "Fortz-Thorup" begin
        # - Below 1/3.
        cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 0.3)
        obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        @objective(rm.model, Min, obj)
        optimize!(rm.model)
        @test termination_status(rm.model) == MOI.OPTIMAL
        @test value(obj) ≈ 0.3
        delete(rm.model, cref)

        # - Between 1/3 and 2/3.
        cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 0.6)
        obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        @objective(rm.model, Min, obj)
        optimize!(rm.model)
        @test termination_status(rm.model) == MOI.OPTIMAL
        @test value(obj) ≈ 1.13333333 atol=1.e-5
        delete(rm.model, cref)

        # - Between 2/3 and 9/10.
        cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 0.8)
        obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        @objective(rm.model, Min, obj)
        optimize!(rm.model)
        @test termination_status(rm.model) == MOI.OPTIMAL
        @test value(obj) ≈ 2.66666667 atol=1.e-5
        delete(rm.model, cref)

        # - Between 9/10 and 1.
        cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 0.95)
        obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        @objective(rm.model, Min, obj)
        optimize!(rm.model)
        @test termination_status(rm.model) == MOI.OPTIMAL
        @test value(obj) ≈ 7.16666667 atol=1.e-5
        delete(rm.model, cref)

        # TODO: the other parts are disabled, it seems like the values are too large for the solver to handle (!?).

        # # - Between 1 and 11/10.
        # cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 1.05)
        # obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        # @objective(rm.model, Min, obj)
        # optimize!(rm.model)
        # @test termination_status(rm.model) == MOI.OPTIMAL
        # @test value(obj) ≈ 35.6666667 atol=1.e-5
        # delete(rm.model, cref)

        # # - Over 11/10.
        # cref = @constraint(rm.model, objective_edge_expression(rm, Load(), Edge(1, 2), dm) == 1.2)
        # obj = objective_edge_expression(rm, FortzThorupLoad(), Edge(1, 2), dm)
        # @objective(rm.model, Min, obj)
        # optimize!(rm.model)
        # @test termination_status(rm.model) == MOI.OPTIMAL
        # @test value(obj) ≈ 560.666667 atol=1.e-5
        # delete(rm.model, cref)
    end
end


text:
function newton(r::Function, x; epsilon = 1e-8, max_num_iter = 20)
    num_iter = 0
    r_x = r(x)
    while num_iter < max_num_iter && abs(r_x) > epsilon
        x += r_x
        r_x = r(x)
        num_iter += 1
    end

    if abs(r_x) <= epsilon
        return x
    else
        return NaN
    end
end

newton(x -> (cos(x)-x)/(sin(x)+1), 0.5)

function polyval(P,x)
    dQ = 0
    Q = P[1]
    for i in 2:length(P)
        dQ = dQ*x + Q
        Q = Q*x + P[i]
    end
    return Q, dQ
end

function r(P, x) # what is it?
    y, dy = polyval(P, x)
    return -y/dy
end

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

using LinearAlgebra

function reverse_gauss_1(A::AbstractMatrix{T}, b::AbstractVector{T}) where T
    x = similar(b)
    N = size(A, 1)
    for k in 0:N-1
        x[N-k] = (b[N-k] - sum(A[N-k,N-k+1:end] .* x[N-k+1:end])) / A[N-k,N-k]
    end
    return x
end

function reverse_gauss_2(A::AbstractMatrix{T}, b::AbstractVector{T}) where T
    x = similar(b)
    N = size(A, 1)
    for k in 0:N-1
        x[N-k] = (b[N-k] - sum(@view(A[N-k,N-k+1:end]) .* @view(x[N-k+1:end]))) / A[N-k,N-k]
    end
    return x
end

function reverse_gauss_3(A::AbstractMatrix{T}, b::AbstractVector{T}) where T
    x = similar(b)
    N = size(A, 1)
    for k in 0:N-1
        x[N-k] = (b[N-k] - dot(@view(A[N-k,N-k+1:end]), @view(x[N-k+1:end]))) / A[N-k,N-k]
    end
    return x
end

function reverse_gauss_4(A,b)
    #A - транспонированная матрица
    N = size(A, 2)
    x = Vector{eltype(Ab)}(undef,N)

    for k in 0:N-1
        x[N-k] = (b[N-k] - dot(@view(A[N-k+1:N, N-k]), @view(x[N-k+1:N])))/A[N-k,N-k]
    end
    return x
end

function reverse_gauss_5(Ab)
    #Ab - транспонированная расширенная матрица
    N = size(A, 2)
    x = Vector{eltype(Ab)}(undef,N)

    for k in 0:N-1
        x[N-k] = (Ab[N-k, N] - dot(@view(Ab[N-k+1:N, N-k]), @view(x[N-k+1:N])))/Ab[N-k,N-k]
    end
    return x
end

function upper_triangular(N::Integer)
    A = randn(N,N)
    _, A = lu(A)
    return A
end

N=10_000
A = upper_triangular(N)
b = ones(N)
At=A'
Abt=[A b]'

println("1:")
@time reverse_gauss_1(A,b);
println("2:")
@time reverse_gauss_2(A,b);
println("3:")
@time reverse_gauss_3(A,b);
println("4:")
@time reverse_gauss_3(At,b); # - это лучший вариант

println("5:")
@time reverse_gauss_5(Abt);
println("...................")

# #  function merge(A, B)
# #  	C = promote_type(eltype(A),eltype(B))[]
# #     j = 1
# #     k = 1
# #     for i = 1:(size(A,1)+size(B,1))
# #        if j>size(A,1)
# # 		   push!(C, B[k])
# # 		   k+=1
# # 		elseif k>size(B,1)
# # 			push!(C, A[j])
# # 			j+=1
# # 		elseif A[j]>B[k]
# #         	push!(C, B[k])
# #         	k+=1
# #     	else
# #     		push!(C, A[j])
# #         	j+=1
# #     	end
# # 	end
# # 	return C
# # end
#
# function merg(A, B)
# 	C = Vector{promote_type(eltype(A),eltype(B))}(undef, length([A ; B]))
# 	j = 1
#     k = 1
#     for i = 1:length(C)
#        if j>length(A)
# 		   C[i] = B[k]
# 		   k+=1
# 	   elseif k>length(B)
# 			C[i] = A[j]
# 			j+=1
# 		elseif A[j]>B[k]
#         	C[i] = B[k]
#         	k+=1
#     	else
#     		C[i] = A[j]
#         	j+=1
#     	end
# 	end
# 	return C
# end
#
#
# function merg!(A, B, C)
# 	if (length(C)!=length([A ; B]))
# 		C = Vector{promote_type(eltype(A),eltype(B))}(undef, length([A ; B]))
# 	end
# 	j = 1
#     k = 1
#     for i = 1:length(C)
#        if j>length(A)
# 		   C[i] = B[k]
# 		   k+=1
# 	   elseif k>length(B)
# 			C[i] = A[j]
# 			j+=1
# 		elseif A[j]>B[k]
#         	C[i] = B[k]
#         	k+=1
#     	else
#     		C[i] = A[j]
#         	j+=1
#     	end
# 	end
# end
#
#
# function part_sort(A, t)
# 	a = []
# 	b = []
# 	c = []
# 	for i = 1:length(A)
# 		if A[i]<t
# 			push!(a, A[i])
# 		elseif A[i]==t
# 			push!(b, A[i])
# 		else
# 			push!(c, A[i])
# 		end
# 	end
# 	[a;b;c]
# end
#
# function part_sort(A, t)
# 	a = []
# 	b = []
# 	c = []
# 	for i = 1:length(A)
# 		if A[i]<t
# 			push!(a, A[i])
# 		elseif A[i]==t
# 			push!(b, A[i])
# 		else
# 			push!(c, A[i])
# 		end
# 	end
# 	[a;b;c]
# end

struct Polynomial{T}
    coeff::Vector{T}

    function Polynomial{T}(coeff) where T
        n = 0
        for c in reverse(coeff)
            if c == 0
                n+=1
            end
        end
        new(coeff[1:end-n])
    end
end

deg(p::Polynomial{T}) where{T} = length(p.coeff) - 1

function eyler(n)
    # n - заданное число
    s = 0.0
    a = 1.0
    for k in 1:n+1
    # границы значений индекса k по отношению к прежнему варианту теперь смещены на 1 вправо
        s += a
        a /= k
    end
end

function eyler()
    s = 0
    a = 1
    k = 0
    while s + a != s
        k += 1
        s += a
        a /= k
    end
    return s
end # Объяснение состоит в том, что
 # при вычислении s+a должны будут складываться
 # соответствующие 53-х битные мантиссы (здесь речь идет
 #  о формате типа Float64). Но сначала порядок меньшего
 #   слагаемого должен быть приведен к порядку большего
 #   слагаемого путем "сдвига" его мантиссы вправо с заменой
 #    "освобождающихся" битов нулями. Таким образом, если сдвиг
 #     составит 53 бита или больше, то мантисса меньшего числа
 #      окажется полностью замененной нулями.

 function eyler_sin(x)
     a = x
     s = 0
     e = 0.000000001
     k = 1
     while abs(a) > e
         s += a
         k += 1
         a *= -x*x/(2*k-1)/(2*k-2)
     end
     s
 end

# sin function
#  function Base.sin(x,e)
#     xx=x^2
#     a=x
#     k=1
#     s=typeof(x)(0) # - преобразование к 0 нужного типа, что обеспечит стабильность типа переменной s
#     while abs(a)>e
#         s+=a
#         a=-a*xx/2k/(2k+1)
#         k+=1
#     end
#     #УТВ: |sin(x)-s|<= ε
#     return s
# end

function Base.sin(x,e)
    xx=x^2
    a=x
    m=2
    s=typeof(x)(0)
    while abs(a)>e
        s+=a
        a=-a*xx/m/(m+1)
        m+=2
    end
    #УТВ: |sin(x)-s|<= ε
    return s
end

# the best sin function
function sin_(x)
    xx=x^2
    a=x
    m=2
    s=typeof(x)(0)
    while s+a != s
        s+=a
        a=-a*xx/m/(m+1)
        m+=2
    end
    #УТВ: |sin(x)-s| <= 1е-16
    return s
end

function exp_(x, eps)
    sum = 1;
    f = x
    k = 2
    while f > eps
        sum += f
        f *= x/k
        k += 1
    end
    sum
end

function harmonic_sum()
    s=0.0
    k=1
    a=1.0
    while s+a != s
        a=1/k
        s+=a
        k+=1
    end
    return s
end

function eyler_cos(x, n)
    mx2 = -x*x
    part = 1
    sum = 1
    for i = 2:2:n
        part *= mx2/((i-1)*i)
        sum += part
    end
    sum
end

function cos_(x, e = 0.0000000000001)
    s = 0
    an = 1
    n = 0
    while abs(an) > e
        s += an
        n += 1
        an *= -x*x/(2*n-1)/(2*n)
    end
    s
end

using Plots

x = 0:0.1:10
p = plot()
plot!(x, eyler_cos.(x, 2))
plot!(x, eyler_cos.(x, 4))
plot!(x, eyler_cos.(x, 8))
plot!(x, eyler_cos.(x, 16))
display(p)

using SpecialFun

function bessel(nu::Real, x::AbstractFloat)
    if isinteger(nu)
        if typemin(Cint) <= nu <= typemax(Cint)
            return bessel(Cint(nu), x)
        end
    elseif x < 0
        throw(DomainError())
    end
    real(bessel(float(nu), complex(x)))
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
# This file includes intregrator construction

function construct_integrator(deproblem, input, righthandside, state, t, modelargs=(), solverargs=();
    alg=nothing, stateder=state, modelkwargs=NamedTuple(), solverkwargs=NamedTuple(), numtaps=3)
    # If needed, construct interpolant for input.
    interpolant = input === nothing ? nothing : Interpolant(numtaps, length(input))

    # Construct the problem
    if deproblem == SDEProblem
        problem = deproblem(righthandside[1], righthandside[2], state, (t, Inf), interpolant, modelargs...;
        modelkwargs...)
    elseif deproblem == DDEProblem
        problem = deproblem(righthandside[1], state, righthandside[2], (t, Inf), interpolant, modelargs...;
         modelkwargs...)
    elseif deproblem == DAEProblem
        problem = deproblem(righthandside, stateder, state, (t, Inf), interpolant, modelargs...;
        modelkwargs...)
    else
        problem = deproblem(righthandside, state, (t, Inf), interpolant, modelargs...;
        modelkwargs...)
    end

    # Initialize the integrator
    init(problem, alg, solverargs...; save_everystep=false, dense=true, solverkwargs...)
end

#= This function checks whether the syntax is of the form

    @my_macro_to_define_new_dynamical_system mutable struct NewSystem{T, S} <: SuperTypeName
        # fields
    end

    where @my_macro_to_define_new_dynamical_system is any macro used to define new dynamical system such as @def_ode_system, @def_sde_system, etc.
=#
function checksyntax(ex::Expr, supertypename::Symbol)
    ex.head == :struct && ex.args[1] ||
        error("Invalid usage. The expression should start with `mutable struct`.\n$ex")
    ex.args[2].head == :(<:) && ex.args[2].args[2] == supertypename ||
        error("Invalid usage. The type should be a subtype of $supertypename.\n$ex")
end


text:
"""
    inclusion_exclusion(P)

Return ℙ(∪ᵢAᵢ) where the `i`th entry in `P` stores ℙ(Aᵢ).

Assumes that A₁, A₂, …, Aₙ are mutually independent events.
"""
function inclusion_exclusion(P::Vector{Float64})
    @assert length(P) != 0

    if length(P) == 1
        return P[1]
    end

    m = length(P) ÷ 2
    first_half = P[1:m]
    second_half = P[m+1:end]

    p = inclusion_exclusion(first_half)
    q = inclusion_exclusion(second_half)
    return p + q - p * q
end

inclusion_exclusion([0.1, 0.25, 0.5])


text:
import Base.==

"""

`intersect(dim, arrs...)`

Take intersection of arrays of type `LabeledArray`/`DictArray`/`AbstractArrayWrapper`. The order is preserved.

##### Arguments

* `dim` : direction along which to intersect.
* `arrs...` : arrays to intersect.

##### Examples

```julia
julia> intersect(1, nalift([1 2 3;4 5 6]), nalift([1 2 3;5 5 6]))
1x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)

julia> intersect(2, nalift([1 2 3;4 5 6]), nalift([1 2 3;5 5 6]))
2x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(2)  Nullable(3)
 Nullable(5)  Nullable(6)

julia> intersect(1, darr(a=[:x,:y,:z]), darr(a=[:x,:x,:y]), darr(a=[:y,:y,:y]))
1 DictArray

a
--
y

julia> intersect(1, larr(a=[1 2 3;4 5 6], axis1=[:X,:Y]), larr(a=[1 2 3;4 3 2], axis1=[:X,:Y]))
1 x 3 LabeledArray

  |1 |2 |3
--+--+--+--
  |a |a |a
--+--+--+--
X |1 |2 |3

julia> intersect(1, larr(a=[1 2 3;4 5 6], axis1=[:X,:Y]), larr(a=[1 2 3;4 3 2], axis1=[:Z,:Y]))
0 x 3 LabeledArray

 |1 |2 |3
-+--+--+--
 |a |a |a
```

"""
Base.intersect(dim::Integer, arr0::LabeledArray, arr_rest::LabeledArray...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might need to be optimized later.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  elems = Dict{Any,Int}()
  isdefaxis = isa(permuted_arrs[1].axes[ndimsarr], DefaultAxis)
  permuted_arrs1 = permuted_arrs[1]
  lastsize = size(permuted_arrs[1], ndimsarr)

  # certainly, this logic is suboptimal...
  if isdefaxis
    for i in 1:lastsize
      full_coords = (remainingcoords..., i)
      # using getindexvalue would have caused the NullException() in the intersect call below.
      # TODO: need to make this intersect method type stable.
      oneslice = getindex(permuted_arrs1.data, full_coords...)
      #push!(elems, oneslice => i)
      if !haskey(elems, oneslice)
        elems[oneslice] = i
      end
    end
  else
    for i in 1:lastsize
      full_coords = (remainingcoords..., i)
      oneslice = getindex(permuted_arrs1.data, full_coords...)
      labels = permuted_arrs1.axes[ndimsarr][i]
      #push!(elems, (oneslice,dummy_nullable_wrap_if_necessary(labels)) => i)
      k = (oneslice,dummy_nullable_wrap_if_necessary(labels))
      if !haskey(elems, k)
        elems[k] = i
      end
    end
  end

  result_arrs = if isdefaxis
    map(permuted_arrs) do arr
      map(1:size(arr, ndimsarr)) do i
        full_coords = (remainingcoords..., i)
        getindex(arr.data, full_coords...)
      end
    end
  else
    map(permuted_arrs) do arr
      map(1:size(arr, ndimsarr)) do i
        full_coords = (remainingcoords..., i)
        oneslice = getindex(arr.data, full_coords...)
        labels = arr.axes[ndimsarr][i]
        (oneslice, dummy_nullable_wrap_if_necessary(labels))
      end
    end
  end
  intersect_slices = intersect(result_arrs...)
  coords_toshow = [elems[slice] for slice in intersect_slices]
  result = permuted_arrs[1][remainingcoords..., coords_toshow]
  permutedims_if_necessary(result, invperm(newdirs))
end

# this is to avoid null exception in doing intersect.
immutable DummyNullableWrapper{T<:Nullable}
  x::T
end

(==)(x::DummyNullableWrapper, y::DummyNullableWrapper) = isnull(x.x) && isnull(y.x) || ignabool(naop_eq(x.x, y.x))

dummy_nullable_wrap_if_necessary(x::Nullable) = DummyNullableWrapper(x)
dummy_nullable_wrap_if_necessary(x) = x

Base.intersect(dim::Int, arr0::DictArray, arr_rest::DictArray...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  result_arrs = map(permuted_arrs) do arr
    map(1:size(arr, ndimsarr)) do i
      full_coords = (remainingcoords..., i)
      getindex(arr, full_coords...)
    end
  end
  intersect_slices = intersect(result_arrs...)
  result = similar(permuted_arrs[1],
                   (size(permuted_arrs[1])[1:end-1]...,
                    length(intersect_slices)))
  for (i,arr) in enumerate(intersect_slices)
    result[remainingcoords...,i] = arr
  end
  permutedims_if_necessary(result, invperm(newdirs))
end

Base.intersect{T,U}(dim::Int,
                    arr0::AbstractArrayWrapper{Nullable{T}},
                    arr_rest::AbstractArrayWrapper{Nullable{U}}...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  result_arrs = map(permuted_arrs) do arr
    map(1:size(arr, ndimsarr)) do i
      full_coords = (remainingcoords..., i)
      getindex(arr, full_coords...)
    end
  end
  intersect_slices = intersect(result_arrs...)
  result = similar(permuted_arrs[1],
                   (size(permuted_arrs[1])[1:end-1]...,
                    length(intersect_slices)))
  #                 (size(result_arrs[1])..., length(intersect_slices)))
  for (i,arr) in enumerate(intersect_slices)
    result[remainingcoords...,i] = arr
  end
  permutedims_if_necessary(result, invperm(newdirs))
end