text:
#=
A backend for POVRay
=#

export POVRay,
       povray,
       povray_material

# We will use MIME types to encode for POVRay
const MIMEPOVRay = MIME"text/povray"

# First, basic types
show(io::IO, ::MIMEPOVRay, s::String) =
  show(io, s)
show(io::IO, ::MIMEPOVRay, s::Symbol) =
  print(io, s)
show(io::IO, ::MIMEPOVRay, r::Real) =
  show(io, r)
show(io::IO, mime::MIMEPOVRay, v::Vector) =
  begin
    write(io, "<")
    for (i, e) in enumerate(v)
      if i > 1
        write(io, ", ")
      end
      show(io, mime, e)
    end
    write(io, ">")
    nothing
  end
show(io::IO, ::MIMEPOVRay, p::Union{Loc,Vec}) =
  # swap y with z to make it consistent with POVRay coordinate system
  let p = in_world(p)
    print(io, "<$(p.x), $(p.z), $(p.y)>")
  end
show(io::IO, ::MIMEPOVRay, c::RGB) =
  print(io, "color rgb <$(Float64(red(c))), $(Float64(green(c))), $(Float64(blue(c)))>")
show(io::IO, ::MIMEPOVRay, c::RGBA) =
  print(io, "color rgbt <$(Float64(red(c))), $(Float64(green(c))), $(Float64(blue(c))), $(Float64(1-alpha(c)))>")

write_povray_object(f::Function, io::IO, type, material, args...) =
  let mime = MIMEPOVRay()
    write(io, "$(type) {")
    for (i, arg) in enumerate(args)
      if i == 1
        write(io, "\n  ")
      elseif i > 1
        write(io, ", ")
      end
      show(io, mime, arg)
    end
    write(io, '\n')
    if ! isnothing(material)
      write_povray_material(io, material)
    end
    let res = f()
      write(io, "}\n")
      res
    end
  end

write_povray_object(io::IO, type, material, args...) =
  write_povray_object(io, type, material, args...) do
    -1 # Default id
  end

write_povray_value(io::IO, value::Any) =
  show(io, MIMEPOVRay(), value)

write_povray_values(io::IO, values) =
  let first = true
    for value in values
      if !first
        print(io, ", ")
      else
        first = false
      end
      write_povray_value(io, value)
    end
  end

write_povray_param(io::IO, name::String, value::Any) =
  begin
    write(io, "  ", name, " ")
    write_povray_value(io, value)
    write(io, '\n')
  end

write_povray_call(io::IO, name::String, values...) =
  begin
    write(io, "  ", name, "(")
    write_povray_values(io, values)
    write(io, ")\n")
  end

write_povray_material(io::IO, material) =
  begin
    write(io, "  ")
    write_povray_value(io, material)
    write(io, '\n')
  end

write_povray_matrix(io::IO, p::Loc) =
  let t = (translated_cs(p.cs, p.x, p.y, p.z).transform)
    write_povray_param(io, "matrix",
      [t[1,1], t[3,1], t[2,1],
       t[1,2], t[3,2], t[2,2],
       t[1,3], t[3,3], t[2,3],
       t[1,4], t[3,4], t[2,4]])
  end

write_povray_camera(io::IO, camera, target, lens) =
  write_povray_object(io, "camera", nothing) do
    write_povray_param(io, "location", camera)
    write_povray_param(io, "look_at", target)
    write_povray_param(io, "right", Symbol("x*image_width/image_height"))
    write_povray_param(io, "angle", view_angles(lens)[1])
  end

write_povray_pointlight(io::IO, location, color) =
  write_povray_object(io, "light_source", nothing, location, color)

write_povray_polygon(io::IO, mat, vs) =
  write_povray_object(io, "polygon", mat, length(vs)+1, vs..., vs[1])

write_povray_polygons(io::IO, mat, vss) =
  write_povray_object(io, "polygon", mat,
    mapreduce(length, +, vss) + length(vss),
    mapreduce(vs->[vs..., vs[1]], vcat, vss)...)

# mesh2 cannot write materials as mesh did.
write_povray_mesh(io::IO, mat, pts, closed_u, closed_v, smooth_u, smooth_v) =
  let (si, sj) = size(pts),
      vcs = add_z.(pts, 1) .- pts,
      idxs = quad_grid_indexes(si, sj, closed_u, closed_v)
    write_povray_object(io, "mesh2", nothing) do
      write_povray_object(io, "vertex_vectors", nothing, si*sj, reshape(permutedims(pts), :)...)
      # Must understand how to handle smoothness along one direction
      if smooth_u && smooth_v
        write_povray_object(io, "normal_vectors", nothing, si*sj, reshape(permutedims(vcs), :)...)
      end
      write_povray_object(io, "face_indices", nothing, length(idxs), idxs...)
      # Must understand how to handle smoothness along one direction
      #write_povray_object(buf, "normal_indices", nothing, length(idxs), idxs...)
      write_povray_material(io, mat)
    end
  end

write_povray_isosurface(io::IO, mat, func, contained_by, max_gradient) =
  write_povray_object(io, "isosurface", mat) do
    write_povray_object(io, "function", nothing, func)
    write_povray_object(io, "contained_by", nothing, contained_by)
    write_povray_param(io, "max_gradient", max_gradient)
  end

#=
struct POVRayPigment
  color::RGB
end

show(io::IO, mime::MIMEPOVRay, p::POVRayPigment) =
  write_povray_object(io, mime, "pigment") do
    show(io, mime, p.color)
  end

#=
NORMAL:
  normal { [NORMAL_IDENTIFIER] [NORMAL_TYPE] [NORMAL_MODIFIER...] }
NORMAL_TYPE:
  PATTERN_TYPE Amount |
  bump_map { BITMAP_TYPE "bitmap.ext" [BUMP_MAP_MODS...]}
NORMAL_MODIFIER:
  PATTERN_MODIFIER | NORMAL_LIST | normal_map { NORMAL_MAP_BODY } |
  slope_map{ SLOPE_MAP_BODY } | bump_size Amount |
  no_bump_scale Bool | accuracy Float
=#

struct POVRayNormal
  type::String
  amount::Real
end

povray_normal(; bumps) =
  POVRayNormal("bumps", bumps)

show(io::IO, mime::MIMEPOVRay, n::POVRayNormal) =
  write_povray_object(io, mime, "normal") do
    write_povray_param(io, mime, n.type, n.amount)
  end

#=
FINISH:
  finish { [FINISH_IDENTIFIER] [FINISH_ITEMS...] }
FINISH_ITEMS:
  ambient COLOR | diffuse [albedo] Amount [, Amount] | emission COLOR |
  brilliance Amount | phong [albedo] Amount | phong_size Amount | specular [albedo] Amount |
  roughness Amount | metallic [Amount] | reflection COLOR |
  crand Amount | conserve_energy BOOL_ON_OFF |
  reflection { Color_Reflecting_Min [REFLECTION_ITEMS...] } |
  subsurface { translucency COLOR } |
  irid { Irid_Amount [IRID_ITEMS...] }
REFLECTION_ITEMS:
  COLOR_REFLECTION_MAX | fresnel BOOL_ON_OFF |
  falloff FLOAT_FALLOFF | exponent FLOAT_EXPONENT |
  metallic FLOAT_METALLIC
IRID_ITEMS:
  thickness Amount | turbulence Amount
=#

struct POVRayFinish
  ambient::Union{RGB,Nothing}
  diffuse::Union{Real,Nothing}
  emission::Union{RGB,Nothing}
  brilliance::Union{Real,Nothing}
  phong::Union{Real,Nothing}
  phong_size::Union{Real,Nothing}
  specular::Union{Real,Nothing}
  roughness::Union{Real,Nothing}
  metallic::Union{Real,Nothing}
  reflection::Union{RGB,Nothing}
  crand::Union{Real,Nothing}
  conserve_energy::Union{Bool,Nothing}
  subsurface::Any
  irid::Any
end

povray_finish(; ambient::Union{RGB,Nothing}=nothing,
                diffuse::Union{Real,Nothing}=nothing,
                emission::Union{RGB,Nothing}=nothing,
                brilliance::Union{Real,Nothing}=nothing,
                phong::Union{Real,Nothing}=nothing,
                phong_size::Union{Real,Nothing}=nothing,
                specular::Union{Real,Nothing}=nothing,
                roughness::Union{Real,Nothing}=nothing,
                metallic::Union{Real,Nothing}=nothing,
                reflection::Union{RGB,Nothing}=nothing,
                crand::Union{Real,Nothing}=nothing,
                conserve_energy::Union{Bool,Nothing}=nothing,
                subsurface::Any=nothing,
                irid::Any=nothing) =
  POVRayFinish(ambient,
               diffuse,
               emission,
               brilliance,
               phong,
               phong_size,
               specular,
               roughness,
               metallic,
               reflection,
               crand,
               conserve_energy,
               subsurface,
               irid)

show(io::IO, mime::MIMEPOVRay, f::POVRayFinish) =
  write_povray_object(io, mime, "finish") do
    for name in fieldnames(typeof(f))
      let v = getfield(f, name)
        if ! isnothing(v)
          write_povray_param(io, mime, string(name), v)
        end
      end
    end
  end

struct POVRayTexture
  pygment::Union{POVRayPigment,Nothing}
  normal::Union{POVRayNormal,Nothing}
  finish::Union{POVRayFinish,Nothing}
end

struct POVRayLibraryTexture
  name::String
end

show(io::IO, mime::MIMEPOVRay, m::POVRayLibraryTexture) =
  write_povray_object(io, mime, "texture") do
    show(io, mime, m.name)
    #show(io, mime, m.interior)
  end

show(io::IO, mime::MIMEPOVRay, t::POVRayTexture) =
  write_povray_object(io, mime, "texture") do
    show(io, mime, t.pygment)
    show(io, mime, t.normal)
    show(io, mime, t.finish)
  end

struct POVRayInterior
  pygment::Union{POVRayPigment,Nothing}
  normal::Union{POVRayNormal,Nothing}
  finish::Union{POVRayFinish,Nothing}
end

struct POVRayMaterial
  name::String
  texture::POVRayTexture
  interior::Union{POVRayInterior,Nothing}
end
=#

abstract type POVRayMaterial end

struct POVRayDefinition <: POVRayMaterial
  name::String
  kind::String
  description::String
end

convert_to_povray_identifier(name) = replace(name, " " => "_")

write_povray_definition(io::IO, d::POVRayDefinition) =
  write(io, "#declare $(convert_to_povray_identifier(d.name)) =\n  $(d.kind) $(d.description)\n")

show(io::IO, ::MIMEPOVRay, d::POVRayDefinition) =
  write(io, "$(d.kind) { $(convert_to_povray_identifier(d.name)) }")

struct POVRayInclude <: POVRayMaterial
  filename::AbstractString
  kind::String
  name::String
end

write_povray_definition(io::IO, m::POVRayInclude) =
  write(io, "#include \"$(m.filename)\"\n")

show(io::IO, ::MIMEPOVRay, m::POVRayInclude) =
  write(io, "$(m.kind) { $(m.name) }")

export povray_definition, povray_include
const povray_definition = POVRayDefinition
const povray_include = POVRayInclude

const povray_concrete =
  povray_definition("Concrete", "texture", """{
   pigment {
     granite turbulence 1.5 color_map {
       [0  .25 color White color Gray75] [.25  .5 color White color Gray75]
       [.5 .75 color White color Gray75] [.75 1.1 color White color Gray75]}}
  finish {
    ambient 0.2 diffuse 0.3 crand 0.03 reflection 0 }
  normal {
    dents .5 scale .5 }}""")

const povray_neutral =
  povray_definition("Material", "texture",
    "{ pigment { color rgb 0.3 } finish { reflection 0 ambient 0 }}")

povray_material(name::String;
                gray::Real=0.3,
                red::Real=gray, green::Real=gray, blue::Real=gray,
                specularity=0, roughness=0,
                transmissivity=nothing, transmitted_specular=nothing) =
  povray_definition(name, "texture", """{
  pigment { rgb <$(Float64(red)),$(Float64(green)),$(Float64(blue))> }
  finish { specular $(specularity) roughness $(roughness) }
}""")

const povray_grass =
  povray_definition("Grass", "texture",
   "{ pigment{ color <0.2,0.5,0.1>*0.5} normal { bumps 0.5 scale 0.005 }}")

const povray_gray_ground =
  povray_definition("Ground", "texture", "{ pigment { color rgb 0.8 } finish { reflection 0 ambient 0 }}")

#=
IMAGE_MAP:
  pigment {
    image_map {
      [BITMAP_TYPE] "bitmap[.ext]" [gamma GAMMA] [premultiplied BOOL]
      [IMAGE_MAP_MODS...]
      }
  [PIGMENT_MODFIERS...]
  }
 IMAGE_MAP:
  pigment {
   image_map {
     FUNCTION_IMAGE
     }
  [PIGMENT_MODFIERS...]
  }
 BITMAP_TYPE:
   exr | gif | hdr | iff | jpeg | pgm | png | ppm | sys | tga | tiff
 IMAGE_MAP_MODS:
   map_type Type | once | interpolate Type |
   filter Palette, Amount | filter all Amount |
   transmit Palette, Amount | transmit all Amount
 FUNCTION_IMAGE:
   function I_WIDTH, I_HEIGHT { FUNCTION_IMAGE_BODY }
 FUNCTION_IMAGE_BODY:
   PIGMENT | FN_FLOAT | pattern { PATTERN [PATTERN_MODIFIERS] }
=#
povray_image_map_material(name::String; image_map, map_type=0) =
  povray_definition(name, "texture", """{
  pigment { image_map $(image_map) map_type $(map_type) }
  }""")

####################################################
# Sky models

povray_realistic_sky_string(turbidity, inner) = """
#version 3.7;
#include "colors.inc"
#include "CIE.inc"
#include "lightsys.inc"
#include "lightsys_constants.inc"
#include "sunpos.inc"
#declare Current_Turbidity = $(turbidity);
//To solve a bug in CIE_Skylight.in
#local fiLuminous=finish{ambient 0 emission 1 diffuse 0 specular 0 phong 0 reflection 0 crand 0 irid{0}}
#include "CIE_Skylight.inc"
global_settings {
  assumed_gamma 1.0
  radiosity {
  }
}
#default {finish {ambient 0 diffuse 1}}
light_source{
  $(inner)
  Light_Color(SunColor,3)
  translate SolarPosition
}
"""
# This one seems to work better. Still 3.5 though.

#=
povray_realistic_sky_string(turbidity, inner) = """
#version 3.5;
#include "colors.inc"
#include "CIE.inc"
#include "lightsys.inc"
#include "lightsys_constants.inc"
#include "sunpos.inc"
global_settings {
  assumed_gamma 1.0
  radiosity {
  }
}
#default {finish {ambient 0 diffuse 1}}
CIE_ColorSystemWhitepoint(Beta_ColSys, Daylight2Whitepoint(Kt_Daylight_Film))
CIE_GamutMapping(off)
#declare Lightsys_Brightness = 1.0;
#declare Lightsys_Filter = <1,1,1>;
#declare Al=38;    // sun altitude
#declare Az=100;   // sun rotation
#declare North=-z;
#declare DomeSize=1e5;
#declare Current_Turbidity = 5.0;
#declare Intensity_Mult = 0.7;
#include "CIE_Skylight.inc"
light_source{ 0
  Light_Color(SunColor,5)
  translate SolarPosition
}
"""
=#

povray_realistic_sky_string(altitude, azimuth, turbidity, withsun) =
  povray_realistic_sky_string(
    turbidity,
    "vrotate(<0,0,1000000000>,<-$(altitude),$(azimuth+180),0>)")

povray_realistic_sky_string(date, latitude, longitude, meridian, turbidity, withsun) =
  povray_realistic_sky_string(
    turbidity,
    """
#local xpto = SunPos($(year(date)), $(month(date)), $(day(date)), $(hour(date)), $(minute(date)), $(meridian), $(latitude), $(longitude));
vrotate(<0,0,1000000000>,<-Al,Az,0>)
""")

####################################################
# Clay models
povray_clay_settings_string() =
"""
#version 3.7;
global_settings {
  assumed_gamma 1.0
  radiosity {
      pretrace_start 64/image_width  //0.04
      pretrace_end 1/image_width     //0.002
      count 1000
      nearest_count 10
      error_bound 0.1
      recursion_limit 1
      low_error_factor 0.2 //0.7
      gray_threshold 0
      minimum_reuse 0.001  //0.01
      brightness 1.2
      adc_bailout 0.01/2
  }
}
sky_sphere {
    pigment {
      color rgb 1
    }
  }
"""

####################################################
abstract type POVRayKey end
const POVRayId = Any
const POVRayRef = NativeRef{POVRayKey, POVRayId}

Base.@kwdef mutable struct POVRayBackend{K,T} <: LazyBackend{K,T}
  shapes::Shapes=Shape[]
  materials::Set{POVRayMaterial}=Set{POVRayMaterial}()
  current_layer::Union{Nothing,AbstractLayer}=nothing
  layers::Dict{AbstractLayer,Vector{Shape}}=Dict{AbstractLayer,Vector{Shape}}()
  sky::String=povray_realistic_sky_string(DateTime(2020, 9, 21, 10, 0, 0), 39, 9, 0, 5, true)
  ground_level::Float64=0.0
  ground_material::POVRayMaterial=povray_gray_ground
  buffer::LazyParameter{IOBuffer}=LazyParameter(IOBuffer, IOBuffer)
  view::View=default_view()
  sun_altitude::Real=90
  sun_azimuth::Real=0
  cached::Bool=false
end

const POVRay = POVRayBackend{POVRayKey, POVRayId}

KhepriBase.backend_name(::POVRay) = "POVRay"

# Traits
# Although we have boolean operation in POVRay at the level of shapes, we do
# not have boolean operations at the level of references
KhepriBase.has_boolean_ops(::Type{POVRay}) = HasBooleanOps{false}()

#=
The POVRay backend cannot realize shapes immediately, only when requested.
=#

save_shape!(b::POVRay, s::Shape) =
  begin
    prepend!(b.shapes, [s])
    if !isnothing(b.current_layer)
      push!(get!(b.layers, b.current_layer, Shape[]), s)
    end
    b.cached = false
    s
  end

save_material(b::POVRay, mat) =
  (push!(b.materials, mat); mat)

used_materials(b::POVRay) = b.materials

KhepriBase.void_ref(b::POVRay) = POVRayRef(-1)

const povray = POVRay()

buffer(b::POVRay) = b.buffer()

KhepriBase.b_trig(b::POVRay, p1, p2, p3, mat) =
  write_povray_object(buffer(b), "triangle", save_material(b, mat), p1, p2, p3)

KhepriBase.b_surface_polygon(b::POVRay, ps, mat) =
  write_povray_object(buffer(b), "polygon", save_material(b, mat), length(ps) + 1, ps..., ps[1])

KhepriBase.b_surface_polygon_with_holes(b::POVRay, ps, qss, mat) =
  let vss = [ps, qss...]
    write_povray_object(buffer(b), "polygon", save_material(b, mat),
      mapreduce(length, +, vss) + length(vss),
      mapreduce(vs->[vs..., vs[1]], vcat, vss)...)
  end

KhepriBase.b_surface_circle(b::POVRay, c, r, mat) =
  	write_povray_object(buffer(b), "disc", save_material(b, mat), c, uvz(c.cs), r)

KhepriBase.b_surface_grid(b::POVRay, ptss, closed_u, closed_v, smooth_u, smooth_v, interpolator, mat) =
  let io = buffer(b),
      mat = save_material(b, mat),
      pts = smooth_u || smooth_v ?
              [location_at(interpolator, u, v)
               for u in division(0, 1, 4*size(ptss, 1)), v in division(0, 1, 4*size(ptss, 2))] :
              ptss,
      (si, sj) = size(pts),
      idxs = quad_grid_indexes(si, sj, closed_u, closed_v)
    write_povray_object(io, "mesh2", nothing) do
      write_povray_object(io, "vertex_vectors", nothing, si*sj, reshape(permutedims(pts), :)...)
      # Must understand how to handle smoothness along one direction
      if smooth_u || smooth_v
        write_povray_object(io, "normal_vectors", nothing, si*sj, reshape(permutedims(add_z.(pts, 1) .- pts), :)...)
      end
      write_povray_object(io, "face_indices", nothing, length(idxs), idxs...)
      #write_povray_object(buf, "normal_indices", nothing, length(idxs), idxs...)
      write_povray_material(io, mat)
    end
  end

KhepriBase.b_cylinder(b::POVRay, cb, r, h, bmat, tmat, smat) =
  let buf = buffer(b)
    write_povray_object(buf, "cylinder", save_material(b, smat), [0,0,0], [0,0,h], r) do
      write_povray_matrix(buf, cb)
    end
    -1
  end

KhepriBase.b_sphere(b::POVRay, c, r, mat) =
  write_povray_object(buffer(b), "sphere", save_material(b, mat), c, r)

KhepriBase.b_torus(b::POVRay, c, ra, rb, mat) =
  let buf = buffer(b)
    write_povray_object(buf, "torus", save_material(b, mat), ra, rb) do
      write_povray_param(buf, "rotate", [90,0,0])
      write_povray_matrix(buf, c)
    end
    -1
  end

KhepriBase.b_box(b::POVRay, c, dx, dy, dz, mat) =
  let buf = buffer(b)
    write_povray_object(buf, "box", save_material(b, mat), [0,0,0], [dx, dy, dz]) do
      write_povray_matrix(buf, c)
    end
    -1
  end

KhepriBase.b_cone(b::POVRay, cb, r, h, bmat, smat) =
  write_povray_object(buffer(b), "cone", save_material(b, smat), cb, r, add_z(cb, h), 0)

KhepriBase.b_cone_frustum(b::POVRay, cb, rb, h, rt, bmat, tmat, smat) =
  write_povray_object(buffer(b), "cone", save_material(b, smat), cb, rb, add_z(cb, h), rt)

# This works (including the smooth bit) BUT
# 1. It does not seem to run faster on POVRay than multiple polygons
# 2. It does not support multiple materials
# 3. The vector must be strictly perpendicular to the profile
# So..., I'm disabling it for now.
# KhepriBase.b_generic_prism(b::POVRay, bs, smooth, v, bmat, tmat, smat) =
#   let buf = buffer(b),
#       cs = cs_from_o_vz(bs[1], -v),
#       ps = [[p.x, p.y] for p in in_cs(bs, cs)]
#     write_povray_object(buf, "prism", save_material(b, smat)) do
#       if smooth
#         write(buf, "linear_sweep cubic_spline ")
#         write_povray_values(buf, [0, norm(v), length(ps) + 3, ps[end], ps..., ps[1], ps[2]])
#       else
#         write(buf, "linear_sweep linear_spline ")
#         write_povray_values(buf, [0, norm(v), length(ps) + 1, ps..., ps[1]])
#       end
#       write_povray_matrix(buf, u0(rotated_x_cs(cs, -pi/2)))
#       -1
#     end
#   end
#

b_isosurface(b::POVRay, frep, bounding_box, mat) =
  write_povray_isosurface(buffer(b), save_material(b, mat), func, bounding_box, max_gradient)

#=
###################################

set_sun(altitude, azimuth, b::POVRay) =
  begin
    b.altitude = altitude
    b.azimuth = azimuth
  end
=#
KhepriBase.b_realistic_sky(b::POVRay, date, latitude, longitude, elevation, meridian, turbidity, withsun) =
  b.sky = povray_realistic_sky_string(date, latitude, longitude, meridian, turbidity, withsun)

# backend_realistic_sky(b::POVRay, altitude, azimuth, turbidity, withsun) =
#   b.sky = povray_realistic_sky_string(altitude, azimuth, turbidity, withsun)

KhepriBase.b_set_ground(b::POVRay, level, mat) =
  (b.ground_level=level; b.ground_material=mat)

KhepriBase.b_delete_all_refs(b::POVRay) =
  begin
    empty!(b.shapes)
    for ss in values(b.layers)
      empty!(ss)
    end
    empty!(b.materials)
    b.cached = false
    nothing
  end

#=
KhepriBase.b_delete_shape(b::POVRay, shape::Shape) =
  let f(s) = s!== shape
    filter!(f, b.shapes)
    for ss in values(b.layers)
      filter!(f, ss)
    end
  end

KhepriBase.b_delete_shapes(b::POVRay, shapes::Shapes) =
  let f(s) = isnothing(findfirst(s1->s1===s, shapes))
    filter!(f, b.shapes)
    for ss in values(b.layers)
      filter!(f, ss)
    end
  end
=#

#=

realize(b::POVRay, s::SweepPath) =
  let vertices = in_world.(path_vertices(s.profile)),
      frames = map_division(identity, s.path, 20),
      buf = buffer(b),
      mat = get_material(b, s)
    write_povray_mesh(
      buf,
      mat,
      [xyz(cx(p), cy(p), cz(p), frame.cs) for p in vertices, frame in frames],
      is_closed_path(s.profile),
      is_closed_path(s.path),
      is_smooth_path(s.profile),
      is_smooth_path(s.path))
    void_ref(b)
  end

# HACK: JUST FOR TESTING
realize(b::POVRay, s::Thicken) =
  realize(b, s.shape)
=#
#=
realize(b::POVRay, s::Move) =
  let buf = buffer(b)
    write_povray_object(buf, "object", nothing) do
      ref(b, s.shape)
      backend_delete_shape(b, s.shape)
      write_povray_param(buf, "translate", s.v)
    end
  end

realize(b::POVRay, s::Scale) =
  let buf = buffer(b),
      trans = s.p-u0()
    write_povray_object(buf, "object", nothing) do
      ref(b, s.shape)
      backend_delete_shape(b, s.shape)
      write_povray_param(buf, "translate", -trans)
      write_povray_param(buf, "scale", s.s)
      write_povray_param(buf, "translate", trans)
    end
  end

realize(b::POVRay, s::Rotate) =
  let buf = buffer(b),
      trans = s.p-u0()
    write_povray_object(buf, "object", nothing) do
      ref(b, s.shape)
      backend_delete_shape(b, s.shape)
      write_povray_param(buf, "translate", -trans)
      write_povray_call(buf, "Axis_Rotate_Trans", s.v, -rad2deg(s.angle))
      write_povray_param(buf, "translate", trans)
    end
  end

realize(b::POVRay, s::UnionShape) =
  let shapes = filter(! is_empty_shape, s.shapes)
    length(shapes) == 1 ?
      (ref(b, shapes[1]); backend_delete_shape(b, shapes[1])) :
      write_povray_object(buffer(b), "union", nothing) do #get_material(b, s)) do
        for ss in shapes
          ref(b, ss)
          backend_delete_shape(b, ss)
        end
      end
    void_ref(b)
  end

realize(b::POVRay, s::IntersectionShape) =
  write_povray_object(buffer(b), "intersection", get_material(b, s)) do
    for ss in s.shapes
      ref(b, ss)
      backend_delete_shape(b, ss)
    end
    void_ref(b)
  end

realize(b::POVRay, s::SubtractionShape3D) =
  write_povray_object(buffer(b), "difference", get_material(b, s)) do
    ref(b, s.shape)
    backend_delete_shape(b, s.shape)
    for ss in s.shapes
      ref(b, ss)
      backend_delete_shape(b, ss)
    end
    void_ref(b)
  end

=#

KhepriBase.b_pointlight(b::POVRay, loc::Loc, color::RGB, range::Real, intensity::Real) =
  write_povray_pointlight(b.buffer(), loc, rgb(red(color)*intensity, green(color)*intensity, blue(color)*intensity))

#=

# BIM

realize_box(b::POVRay, mat, p, dx, dy, dz) =
  let buf = buffer(b),
      bot = p,
      top = add_xyz(p, dx, dy, dz)
    write_povray_object(buf, "box", mat, bot, top)
    void_ref(b)
  end

realize_prism(b::POVRay, top, bot, side, path::PathSet, h::Real) =
  # PathSets require a different approach
  let buf = buffer(b),
      v = planar_path_normal(path)*h,
      bot_vss = map(path_vertices, path.paths),
      top_vss = map(path_vertices, translate(path, v).paths)
    write_povray_polygons(buf, bot, map(reverse, bot_vss))
    write_povray_polygons(buf, top, top_vss)
    for (bot_vs, top_vs) in zip(bot_vss, top_vss)
      for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
        write_povray_polygon(buf, side, vs)
      end
    end
    void_ref(b)
  end

realize_pyramid_frustum(b::POVRay, top, bot, side, bot_vs::Locs, top_vs::Locs, closed=true) =
  let buf = buffer(b)
    if closed
      write_povray_polygon(buf, bot, reverse(bot_vs))
      write_povray_polygon(buf, top, top_vs)
    end
    for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
      write_povray_polygon(buf, side, vs)
    end
    void_ref(b)
  end

#=
POVRay families need to know the different kinds of materials
that go on each surface.
In some cases it might be the same material, but in others, such
as slabs, outside walls, etc, we will have different materials.
=#

const POVRayMaterialFamily = BackendMaterialFamily{POVRayMaterial}
povray_material_family(mat::POVRayMaterial) =
  POVRayMaterialFamily(mat)

const POVRaySlabFamily = BackendSlabFamily{POVRayMaterial}
povray_slab_family(top::POVRayMaterial, bot::POVRayMaterial=top, side::POVRayMaterial=bot) =
  POVRaySlabFamily(top, bot, side)

const POVRayRoofFamily = BackendRoofFamily{POVRayMaterial}
povray_roof_family(top::POVRayMaterial, bot::POVRayMaterial=top, side::POVRayMaterial=bot) =
  POVRayRoofFamily(top, bot, side)

const POVRayWallFamily = BackendWallFamily{POVRayMaterial}
povray_wall_family(right::POVRayMaterial, left::POVRayMaterial=right) =
  POVRayWallFamily(right, left)

=#
export povray_family_materials

povray_family_materials(m1, m2=m1, m3=m2, m4=m3) = (materials=(m1, m2, m3, m4), )

export povray_stone, povray_metal, povray_wood, povray_glass

povray_stone = povray_include("stones2.inc", "texture", "T_Stone35")
povray_metal = povray_include("textures.inc", "texture", "Chrome_Metal")
#povray_wood = povray_include("textures.inc", "texture", "DMFWood1")
povray_wood = povray_include("woods.inc", "texture", "T_Wood10")
povray_glass = povray_include("textures.inc", "material", "M_Glass")

#=

use_family_in_layer(b::POVRay) = true

#=
create_ground_plane(shapes, material=default_povray_ground_material()) =
  if shapes == []
    error("No shapes selected for analysis. Use add-povray-shape!.")
  else
    let (p0, p1) = bounding_box(union(shapes)),
        (center, ratio) = (quad_center(p0, p1, p2, p3),
                  distance(p0, p4)/distance(p0, p2));
     ratio == 0 ?
      error("Couldn"t compute height. Use add-povray-shape!.") :
      let pts = map(p -> intermediate_loc(center, p, ratio*10), [p0, p1, p2, p3]);
         create_surface_layer(pts, 0, ground_layer(), material)
        end
       end
  end
        w = max(floor_extra_factor()*distance(p0, p1), floor_extra_width())
        with(current_layer,floor_layer()) do
          box(xyz(min(p0.x, p1.x)-w, min(p0.y, p1.y)-w, p0.z-1-floor_distance()),
              xyz(max(p0.x, p1.x)+w, max(p0.y, p1.y)+w, p0.z-0-floor_distance()))
        end
      end
    end

=#

realize(b::POVRay, s::Union{Door, Window}) =
  void_ref(b)
=#

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

export_to_povray(b::POVRay, path::String) =
  let buf = b.buffer()
    if ! b.cached
      take!(buf)
      for s in b.shapes
        mark_deleted(b, s)
      end
      i = 1
      while i <= length(b.shapes)
        force_realize(b, b.shapes[i])
        i += 1
      end
      b.cached = true
    end
    open(path, "w") do out
      # write the sky
      write(out, b.sky)
      # write useful include
      write(out, """#include "transforms.inc"\n""")
      # write materials (removing duplicate includes)
      save_material(b, b.ground_material)
      for m in unique(m->m isa POVRayInclude ? m.filename : m, used_materials(b))
        write_povray_definition(out, m)
      end
      # write the ground
      write_povray_object(out, "plane", b.ground_material, :y, b.ground_level)
      # write the objects
      let str = String(take!(buf))
        #save again
        println(buf, str)
        write(out, str)
      end
      # write the view
      write_povray_camera(out, b.view.camera, b.view.target, b.view.lens)
    end
  end

povray_cmd() =
  @static Sys.iswindows() ?
    "C:/Program Files/POV-Ray/v3.7/bin/pvengine64" :
    (@static Sys.isunix() ?
      realpath(normpath(@__DIR__, "..", "bin", "povray")) :
      error("Unsupported operating system!"))

const povray_lib = Parameter(realpath(normpath(@__DIR__, "..", "lib", "include")))
const LightsysIV_lib = Parameter(realpath(normpath(@__DIR__, "..", "lib", "LightsysIV")))

##########################################
export rendered_image_area
rendered_image_area = Parameter{Union{Missing,Tuple{Int,Int,Int,Int}}}(missing)

KhepriBase.b_render_view(b::POVRay, path::String) =
  let povpath = path_replace_suffix(path, ".pov"),
      inipath = path_replace_suffix(path, ".ini"),
      options = @static(Sys.iswindows() ? (film_active() ? ["-D", "/EXIT", "/RENDER"] : ["/RENDER"]) : []),
      cmd = `$(povray_cmd()) $options "$inipath"`
    @info cmd
    export_to_povray(b, povpath)
    open(inipath, "w") do out
      #println(out, "Debug_Console=off")
      println(out, "Library_Path='$(povray_lib())'")
      println(out, "Library_Path='$(LightsysIV_lib())'")
      println(out, "Width=$(render_width())")
      println(out, "Height=$(render_height())")
      println(out, "Antialias=on")
      println(out, "Quality=10")
      println(out, "High_Reproducibility=on")
      println(out, "Input_File_Name='$(povpath)'")
      if ! ismissing(rendered_image_area())
        let (start_col, start_row, end_col, end_row) = rendered_image_area()
          println(out, "Start_Row=", start_row)
          println(out, "Start_Column=", start_col)
          println(out, "End_Row=", end_row)
          println(out, "End_Column=", end_col)
        end
      end
    end
    run(cmd, wait=@static(Sys.iswindows() ? film_active() : true))
  end

export clay_model, realistic_model
clay_model(level::Real=0, b::POVRay=povray) =
  begin
    b.sky = povray_clay_settings_string()
    b.ground_level = level
    b.ground_material = povray_definition("Ground", "texture", "{ pigment { color rgb 3 } finish { reflection 0 ambient 0 }}")
  end

realistic_model(level::Real=0, b::POVRay=povray) =
  begin
    b.sky = povray_realistic_sky_string(DateTime(2020, 9, 21, 10, 0, 0), 39, 9, 0, 5, true)
    b.ground_level = level
    b.ground_material = povray_gray_ground
  end


text:
"""
    getvalue(document::Mongoc.BSON, feature_name::String)

Retrieve the value of a field inside a BSON document. Return `missing` if the field doesn't exist.
"""
function getvalue(document::Mongoc.BSON, field_name::String)
    field_name_split = split(field_name, ".")
    for field_name_level in field_name_split
        if !haskey(document, field_name_level)
            return missing
        elseif typeof(document[field_name_level]) <: Dict
            document = document[field_name_level]
        else
            return document[field_name_level]
        end
    end
end

"""
    getcolnamesdict(collection_fields::Mongoc.Collection; collectiontype::Symbol=:name)

Construct the dict with the format ("features names" => []) of a collection (either instances or decompositions).
See also [`getfeaturesdf`]
"""
function getcolnamesdict(collection_keys::Mongoc.Collection, names::Vector{String}; collectiontype::Symbol=:name)
    collectiontype == :name && (collectiontype = getcollectiontype(collection_keys))
    colnames = Dict(val => Vector{Union{Missing,String}}(undef,0) for val in names)
    colnames["instance_name"] = Vector{Union{Missing,String}}(undef,0)
    if collectiontype == :decompositions
        colnames["nb_added_edges"] = Vector{Union{Missing,String}}(undef,0)
        colnames["origin_oid"] = Vector{Union{Missing,String}}(undef,0)
        colnames["oid"] = Vector{Union{Missing,String}}(undef,0)
        colnames["nb_lc"] = Vector{Union{Missing,String}}(undef,0)
        colnames["dec_type"] = Vector{Union{Missing,String}}(undef,0)
        colnames["merge_treshold"] = Vector{Union{Missing,String}}(undef,0)
        #colnames["merge_type"] = Vector{Union{Missing,String}}(undef,0)
        colnames["origin_nb_added_edges"] = Vector{Union{Missing,String}}(undef,0)
        colnames["origin_dec_type"] = Vector{Union{Missing,String}}(undef,0)
    end
    return colnames
end

"""
    getcollectiontype(collection::Mongoc.Collection)

Retrieve the collection content type of `collection`.
"""
function getcollectiontype(collection::Mongoc.Collection)
    if issubset("instances", collection.name)
        return :instances
    elseif issubset("decompositions", collection.name)
        return :decompositions
    else
        throw(ArgumentError("Could not retrieve if $(collection) corresponds to the instances or decompositions collection."))
    end
end

"""
    fill!(df::DataFrame, documents::Vector{AbstractDict}, names::Vector{String}, collectiontype::Symbol)

Populate `df` with the values in `documents`.
"""
function fill!(df::DataFrame, documents::Vector{N}, names::Vector{String}, collectiontype::Symbol) where N<:AbstractDict
    @info "Filling $(collectiontype) dataframe"
    for (i, document) in enumerate(documents)
        features = Dict{String,Union{String,Missing}}()
        if collectiontype == :instances
            features["instance_name"] = document["_id"]
        else
            !haskey(document["features"], "solver") && continue
            features["instance_name"] = document["_id"]["instance_name"]
            features["oid"] = document["_id"]["oid"]
            features["nb_added_edges"] = string(document["nb_added_edges"])
            features["nb_lc"] = string(getvalue(document, "features.solver.nlc"))
            features["dec_type"] = string(getvalue(document, "dec_type"))
            if document["nb_added_edges"] == -1
                features["origin_oid"] = document["origin_oid"]
                features["merge_treshold"] = string(document["merge_treshold"])
                #features["merge_type"] = string(document["merge_type"])
                features["origin_nb_added_edges"] = string(document["origin_nb_added_edges"])
                if document["origin_nb_added_edges"] == 0
                    features["origin_dec_type"] = string(document["origin_dec_type"])
                else
                    features["origin_dec_type"] = ""
                end
            else
                features["origin_oid"] = ""
                features["merge_treshold"] = ""
                #features["merge_type"] = ""
                features["origin_nb_added_edges"] = ""
                features["origin_dec_type"] = ""
            end
        end
        for feature_name in names
            feature_name == "nb_added_edges" && continue
            feature = getvalue(document, "features.$feature_name")
            features[feature_name] = ismissing(feature) ? feature : string(feature)
        end
        push!(df, features)
    end
end

ext: jl

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
function get_inside_splitters(obj::Characterization)
    return jcall(obj, "getInsideSplitters", NodesSet, ())
end

function get_outside_splitters(obj::Characterization)
    return jcall(obj, "getOutsideSplitters", NodesSet, ())
end

function inside_touching(obj::Characterization)
    return jcall(obj, "insideTouching", SubHyperplane, ())
end

function outside_touching(obj::Characterization)
    return jcall(obj, "outsideTouching", SubHyperplane, ())
end

function touch_inside(obj::Characterization)
    return jcall(obj, "touchInside", jboolean, ())
end

function touch_outside(obj::Characterization)
    return jcall(obj, "touchOutside", jboolean, ())
end


text:
@testset "WignerFunction" begin
    dim1 = 10
    dim2 = 11
    xs = -1:0.1:1
    ps = -1:0.1:1

    # w/o mmap
    for _ in 1:2
        wf = WignerFunction(xs, ps, dim=dim1)
        state = VacuumState(Matrix, dim=dim1)

        w = wf(state)
        ans = real(sum(state .* wf.𝐰, dims=(1, 2)))
        @test all(e == ans[i] for (i, e) in enumerate(w.𝐰_surface))

        @test size(wigner(VacuumState(Matrix, dim=dim2), xs, ps).𝐰_surface) == (length(xs), length(ps))
    end
end

@testset "WignerSurface" begin
    _range = -10:0.1:10
    len = length(_range)
    @test WignerSurface(_range, _range, ones(Float64, len, len)) isa
        WignerSurface{typeof(_range)}
end


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((1, 1, 1, 2, true, true, false), nothing, relaxed_tols),
    ((1, 1, 1, 2, true, false, true),),
    ((1, 1, 1, 2, false, true, false),),
    ((1, 1, 1, 2, false, false, true),),
    ((1, 1, 1, 2, false, false, false),),
    ]
insts["fast"] = [
    ((2, 2, 2, 2, true, true, false),),
    ((2, 2, 2, 2, true, false, true),),
    ((2, 2, 2, 2, false, true, false),),
    ((2, 2, 2, 2, false, false, true),),
    ((2, 2, 2, 2, false, false, false),),
    ((4, 2, 3, 5, true, true, false),),
    ((4, 2, 3, 5, true, false, true),),
    ]
insts["various"] = vcat(insts["fast"], [
    ((4, 2, 3, 5, false, true, false),),
    ((4, 2, 3, 5, false, false, true),),
    ((4, 2, 3, 5, false, false, false),),
    ])
return (PolyNormJuMP, insts)


text:
#just run this script to execute most of the functions

#eg. main(1980,1990,5,0.05)
function main(startYr = 1980, endYr = 1990, windowSize = 5, alpha = 0.05)
    #produces the plots for collusion 1 and 2 way
    include("demo.jl")
    demo(startYr,endYr,windowSize)

    #calls csvDir2AdjList.jl and returns 2 element array with the upper and lower which are the high score pairings and lower score pairings (collusion vs neglect)
    include("countryBiasESC.jl")
    biasESC_uppper_lower = biasesESC(startYr,endYr,windowSize,alpha)

    #generate the plots for the overall associations and trends
    include("biasAnalysis.jl")
    analyzeBiases(biasESC_uppper_lower[1],biasESC_uppper_lower[2])

    #produce the graphs of both analysis
    include("biasGraphView.jl")
    graphAvoid(biasESC_uppper_lower[1],biasESC_uppper_lower[2])

end

max_stars_repo_name: mantzaris/eurovisionBias

max_issues_repo_name: mantzaris/eurovisionBias

max_forks_repo_name: mantzaris/eurovisionBias


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3ecc9360-d570-11eb-0a04-ad57130a3b16
using HypertextLiteral

# ╔═╡ a6033a20-bdb6-41bd-bb31-41e827f756e0
@bind screenWidth @htl("""
	<div>
	<script>
		var div = currentScript.parentElement
		div.value = screen.width
	</script>
	</div>
""")

# ╔═╡ c971f8aa-833f-4ec4-99a6-c3e10437c57a
begin
	cellWidth= min(1000, screenWidth*0.9)
	@htl("""
		<style>
			pluto-notebook {
				margin: auto;
				width: $(cellWidth)px;
			}
		</style>
	""")
end

# ╔═╡ Cell order:
# ╠═3ecc9360-d570-11eb-0a04-ad57130a3b16
# ╠═a6033a20-bdb6-41bd-bb31-41e827f756e0
# ╠═c971f8aa-833f-4ec4-99a6-c3e10437c57a

max_stars_repo_path: 0007/Pluto 0.9screenWidth.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0007/Pluto 0.9screenWidth.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0007/Pluto 0.9screenWidth.jl

max_forks_repo_name: genkuroki/public


text:
module Transducers

import Base: map, filter

function map(f::Function)
    function (rf::Function)
        rf′() = rf()
        rf′(ret) = rf(ret)
        rf′(ret, input) = rf(ret, f(input))
        return rf′
    end
end

function filter(pred::Function)
    function (rf::Function)
        rf′() = rf()
        rf′(ret) = rf(ret)
        rf′(ret, input) = pred(input) ? rf(ret, input) : ret
        return rf′
    end
end

function transduce(xform::Function, rf::Function, init, xs)
    rf = xform(rf)
    return reduce(rf, xs; init=init)
end

transduce(xform::Function, rf::Function, xs) = transduce(xform, rf, rf(), xs)

end


text:
######### Jags line program example  ###########

using Mamba, Jags

ProjDir = joinpath(dirname(@__FILE__))
cd(ProjDir) do

  line = "
  model {
    for (i in 1:n) {
          mu[i] <- alpha + beta*(x[i] - x.bar);
          y[i]   ~ dnorm(mu[i],tau);
    }
    x.bar   <- mean(x[]);
    alpha    ~ dnorm(0.0,1.0E-4);
    beta     ~ dnorm(0.0,1.0E-4);
    tau      ~ dgamma(1.0E-3,1.0E-3);
    sigma   <- 1.0/sqrt(tau);
  }
  "

  monitors = Dict(
    "alpha" => true,
    "beta" => true,
    "tau" => true,
    "sigma" => true
    )

  jagsmodel = Jagsmodel(
    name="line1",
    model=line,
    monitor=monitors,
    deviance=true, dic=true, popt=true,
    pdir=ProjDir
    );

  println("\nJagsmodel that will be used:")
  jagsmodel |> display

  data = Dict(
    "x" => [1, 2, 3, 4, 5],
    "y" => [1, 3, 3, 3, 5],
    "n" => 5
  )

  inits = [
    Dict("alpha" => 0,"beta" => 0,"tau" => 1),
    Dict("alpha" => 1,"beta" => 2,"tau" => 1),
    Dict("alpha" => 3,"beta" => 3,"tau" => 2),
    Dict("alpha" => 5,"beta" => 2,"tau" => 5)
  ]

  println("Input observed data dictionary:")
  data |> display
  println("\nInput initial values dictionary:")
  inits |> display
  println()

  sim = jags(jagsmodel, data, inits, ProjDir)
  describe(sim)
  println()

  ## Brooks, Gelman and Rubin Convergence Diagnostic
  try
    gelmandiag(sim1, mpsrf=true, transform=true) |> display
  catch e
    #println(e)
    gelmandiag(sim, mpsrf=false, transform=true) |> display
  end

  ## Geweke Convergence Diagnostic
  gewekediag(sim) |> display

  ## Highest Posterior Density Intervals
  hpd(sim) |> display

  ## Cross-Correlations
  cor(sim) |> display

  ## Lag-Autocorrelations
  autocor(sim) |> display

  ## Deviance Information Criterion
  #dic(sim) |> display

  ## Plotting
  p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
  draw(p, nrow=4, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:svg)
  draw(p, nrow=4, ncol=4, filename="$(jagsmodel.name)-summaryplot", fmt=:pdf)

  # Below will only work on OSX, please adjust for your environment.
  # JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
  @static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
          for i in 1:4
            isfile("$(jagsmodel.name)-summaryplot-$(i).svg") &&
              run(`open -a $(JULIA_SVG_BROWSER) "$(jagsmodel.name)-summaryplot-$(i).svg"`)
          end
        end : println()

  # Below examples of using other ways to display the simulation results
  (index, chains) = Jags.read_jagsfiles(jagsmodel)

  println()
  #chains[1]["samples"] |> display
  #println()

  if jagsmodel.dic
    (idx0, chain0) = Jags.read_pDfile(jagsmodel)
    #idx0 |> display
    println()
    #chain0[1]["samples"] |> display
  end

  if jagsmodel.dic || jagsmodel.popt
    pDmeanAndpopt = Jags.read_table_file(jagsmodel, data["n"]);
    pDmeanAndpopt["pD.mean"] |> display
    println()
    pDmeanAndpopt["popt"] |> display
  end

end #cd


text:
module Ssfm
using FFTW
using Plots
using Printf
import Plotly
using ProgressMeter
export solve, config

struct Fiber
  T::Float64
  dt::Float64
  L::Float64
  dz::Float64
  alpha::Float64 # dB/km
  beta_2::Float64 # ps^2/km
  gamma::Float64 # km^-1
end

struct Pulse
  P0::Float64
  T0::Float64
  sech_flag::Float64
end

function solve(fiber::Fiber, pulse::Pulse)
  collision = 0;
  alpha = fiber.alpha * 0.23025 / 1e3;
  time_steps = Int(floor(fiber.T/fiber.dt))
  space_steps = Int(floor(fiber.L/fiber.dz))

  LD = pulse.T0^2 / abs(fiber.beta_2)
  LNL = 1 / (fiber.gamma*pulse.P0)

  disp_ratio = fiber.L / LD
  nonl_ratio = fiber.L / LNL
  print("Simulation diagnositics:\n")
  @printf("\tT/T0  = %6.4f  window form factor\n", fiber.T/pulse.T0)
  @printf("\tL/LD  = %6.4f\n", disp_ratio)
  @printf("\tL/LNL = %6.4f\n", nonl_ratio)
  @printf("\tN     = %6.4f \n", sqrt(LD / LNL))
  @printf("time steps : %6d \n", time_steps)
  @printf("space steps: %6d \n", space_steps)

  time = LinRange(-fiber.T/2, fiber.T/2, time_steps)
  space = LinRange(0, fiber.L, space_steps)

  w = 2*pi* LinRange(-1/(2*fiber.dt), 1/(2*fiber.dt), time_steps)
  w = fftshift(w)
  A = zeros(ComplexF64, time_steps, space_steps)
  A_spect = zeros(ComplexF64, time_steps, space_steps)

  if (pulse.sech_flag == 1)
    waveform = sqrt(pulse.P0) .* 2 ./(exp.(time/pulse.T0).+exp.(-time/pulse.T0))
    if collision == 1
      waveform = sqrt(pulse.P0) .* ( 2 ./ (exp.((time.-10*pulse.T0)./pulse.T0).+exp.(-(time.-10*pulse.T0)./pulse.T0)) + 2 ./ (exp.((time.+10*pulse.T0)./pulse.T0).+exp.(-(time.+10*pulse.T0)./pulse.T0)))
    end
  else
    waveform = sqrt(pulse.P0) .* exp.(-(time/pulse.T0) .^ 2 ./ 2)
  end

  A[:, 1] = waveform
#  fig = Plots.plot(time, abs.(A[:, 1]), show=true)
  A_spect[:, 1] = fft(A[:, 1])

  fwd_disp = exp.(fiber.dz * (im * fiber.beta_2/2 * w .^ 2 .- alpha/2))
  @showprogress "Propagating the field... " for n = 1:space_steps-1
    A_spect[:, n] = fft(A[:, n])
    A_spect[:, n+1] = A_spect[:, n] .* fwd_disp
    A[:, n+1] = ifft(A_spect[:, n+1])
    A[:, n+1] = A[:, n+1] .* exp.(im * fiber.dz * fiber.gamma .* (abs.(A[:, n+1])).^2)
  end
  A_spect[:, space_steps] = fft(A[:, space_steps])

  tmargin_l = Int(floor(time_steps/2 - 10 * pulse.T0/fiber.dt))
  tmargin_r = Int(ceil(time_steps/2 + 10 * pulse.T0/fiber.dt))
  t_points = 10000
  z_points = 1000
  t_skip = Int(ceil(20 * pulse.T0/fiber.dt /t_points))
  z_skip = Int(ceil(space_steps/z_points))
  fig2 = Plots.heatmap(abs.(A[tmargin_l:t_skip:tmargin_r, 1:z_skip:space_steps]), show=true)
  # fig3 = Plots.surface(abs.(A[tmargin_l:t_skip:tmargin_r, 1:z_skip:space_steps]), show=true)
end

function mem_estimate(fiber::Fiber)
  return 2* fiber.T/fiber.dt * fiber.L/fiber.dz * 128
end

function config()
  Plots.plotly()
  Plots.default(size=(1300, 800),
                guidefont=("times", 10),
                legendfont=("times", 12),
                tickfont=("times", 10)
                )
end
config()

## Fibers
attenuative = Fiber(20e-9,
               1e-12,
               50000,
               100,
               0.2,
               0e-27,
               0.00e-3)
dispersive = Fiber(2e-9,
               1e-13,
               77500,
               100,
               0.0,
               20e-27,
               0.00e-3)
nonlinear = Fiber(20e-9,
               1e-12,
               27000,
               100,
               0.0,
               0e-27,
               3.00e-3)
solitonic = Fiber(20e-9,
               1e-12,
               27000,
               100,
               0.0,
               -4e-27,
               3.00e-3)
smf_28 =Fiber(10e-9,
              1e-13,
              55500,
              100,
              0.18,
              -21.58e-27,
              0.78e-3)
hnlf = Fiber(1e-9,
            1e-13,
            5300,
            1,
            0.81,
            -7.09e-27,
            10.68e-3)
## Pulses
pulse1 =  Pulse(0.001,
                20e-12,
                0)
pulse2 =  Pulse(0.1,
                20e-12,
                0)
pulse3 =  Pulse(0.46,
                7e-12,
                1)
pulse4 =  Pulse(0.152,
                2.09e-12,
                1)
high_power =  Pulse(0.152*9,
                2.09e-12,
                1)
soliton = Pulse(3.4,
                20e-12,
                1)
## Configuration
configs = [(hnlf, pulse4)      # first order soliton (attenuation limited)
           (hnlf, high_power)] # third order soliton (attenuation limited)

mem_limit = 15000000000 #byte
for (fiber, pulse) in configs
  estimate = mem_estimate(fiber)
  if estimate < mem_limit
    @printf("Estimated memory consumption: %4.1f MiB\n", estimate/(1024^2))
    @time solve(fiber, pulse)
  else
    @printf("Estimated memory consumption (%4.1f MiB) exceed maximum!\n", estimate/(1024^2))
  end
end

end

max_stars_repo_name: lorenzifrancesco/nonlinear-pulse

max_issues_repo_name: lorenzifrancesco/nonlinear-pulse

max_forks_repo_name: lorenzifrancesco/nonlinear-pulse


text:
include("../src/scaler.jl")

bbox = [-1.0 -1.0; 1.0 1.0]
@info "Boundy box = " bbox
@info "x1 = " bbox[1]
@info "y1 = " bbox[2]

scaler = Scaler(bbox)

@testset "Scale some random points" begin
    a = Point2D(-1.0, 0.0)
    b = Point2D(0.0, 1.0)
    c = Point2D(0.5, 0.5)
    esa = Point2D(0.0, 0.5)
    esb = Point2D(0.5, 1.0)
    esc = Point2D(0.75, 0.75)
    sa = scaledpoint2d(a, scaler)
    sb = scaledpoint2d(b, scaler)
    sc = scaledpoint2d(c, scaler)
    @test esa == sa
    @test esb == sb
    @test esc == sc
end

@testset "Proper scale" begin
    expected_scale = 0.5
    @test expected_scale == scaler.scale
end

@testset "Scale then unscale point" begin
    p = Point2D(-1.0, -1.0)
    sp = Point2D(0.0, 0.0)
    spt = scaledpoint2d(p, scaler)
    @test sp == spt
    uspt = unscaledpoint2d(spt, scaler)
    @test p == uspt
end

@testset "scaledpoint is equal to scaledpoint2d" begin
    p2d = Point2D(-0.3212, 0.6654)
    p = [-0.3212, 0.6654]
    sp2d = scaledpoint2d(p2d, scaler)
    sp = scaledpoint(p, scaler)
    @test getx(sp2d) == sp[1]
    @test gety(sp2d) == sp[2]
end

@testset "unscaledpoint is equal to unscaledpoint2d" begin
    p2d = Point2D(0.3212, 0.6654)
    p = [0.3212, 0.6654]
    up2d = unscaledpoint2d(p2d, scaler)
    up = unscaledpoint(p, scaler)
    @test getx(up2d) == up[1]
    @test gety(up2d) == up[2]
end

@testset "unscale x and y" begin
    p = [-0.3212, 0.6654]
    sp = scaledpoint(p, scaler)
    @test abs(unscalex(sp[1], scaler) - p[1]) < eps(Float64)
    @test abs(unscaley(sp[2], scaler) - p[2]) < eps(Float64)
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
#=
To Do :

Use Static Arrays for better performance of 2D,3D points,vectors and matrices
https://github.com/JuliaArrays/StaticArrays.jl

Performance optimization tips:
https://docs.julialang.org/en/v1/manual/performance-tips/

* Optimization variables
  http://www.juliaopt.org/JuMP.jl/v0.19.0/variables/
   * Lazy broadcasting:
     http://www.cityinthesky.co.uk/2019/01/working-with-broadcasting-in-julia/

=#

#=

Julia Cheatsheet
https://juliadocs.github.io/Julia-Cheat-Sheet/
https://cheatsheets.quantecon.org/
https://learnxinyminutes.com/docs/julia/
https://juliabyexample.helpmanual.io/
https://julia.guide/broadcasting

b = [x for x in 1:0.5:5]'    # Row vector comprehension
a = [2, 3, 4, 5]             # Column vectro comprehension

Outer like operation
floor.(Int, sqrt.(a) .+ sqrt.(b))
@. floor(Int, sqrt(a) + sqrt(b))

ps = [ [ix*0.1 iy*0.1] for ix=1:10, iy=1:10 ]   # 2D comprenesion

in vscode : Open in Bach not Terminal Tools

// https://docs.julialang.org/en/v1/base/file/
pwd()
readdir()
include("pointPlacement.jl")

http://docs.juliaplots.org/latest/attributes/
scatter(y, marker = (:hexagon, 20, 0.6, :green, stroke(3, 0.2, :black, :dot)))
scatter(y, markershape = :hexagon,
    markersize = 20,
    markeralpha = 0.6,
    markercolor = :green,
    markerstrokewidth = 3,
    markerstrokealpha = 0.2,
    markerstrokecolor = :black,
    markerstrokestyle = :dot)

plot(y, xaxis = ("my label", (0,10), 0:0.5:10, :log, :flip, font(20, "Courier")))

plot(y, xlabel = "my label",
    xlims = (0,10),
    xticks = 0:0.5:10,
    xscale = :log,
    xflip = true,
    xtickfont = font(20, "Courier"))

=#

using LinearAlgebra

function hexGrid(n::Int,L::Float64)
    l::Float64  = L/n
    h::Float64  = l*sin(pi/3);
    #np = ( 2*n*n + numerator(n*(n-1)//2) )
    np::Int = ( (2*n+1)*(n+1) + n*n )
    #println( n," ",l," ",h," ",np )
    pts = Array{Float64}(undef,np,2)
    i::Int = 1
    x0::Float64 = -n*l*0.5
    y0::Float64 = n*h
    println( "n, x0, y0 ", n," ", x0," ", y0 )
    y = y0
    for iy::Int in 0:n
        x   = x0
        for ix::Int in 0:(n+iy)
            pts[i,1] = x
            pts[i,2] = y
            x += l
            i += 1
        end
        y  -= h
        x0 -= 0.5*l
    end
    x0 += 0.5*l
    for iy::Int in 0:(n-1)
        x0 += 0.5*l
        x   = x0
        for ix::Int in 0:(2*n-iy-1)
            pts[i,1] = x
            pts[i,2] = y
            x += l
            i += 1
        end
        y  -= h
    end
    #print( "iend = ", i )
    return pts
end

@inline function R2func(r2::Float64)
    r2 = 1-r2
    return r2*r2
end

@inline function R2func_deriv(r2::Float64, invR2::Float64 )
    r2  *= invR2
    drf  = 1-r2
    return drf*drf, -4*drf*invR2
end

R2 = 0.04

function projectRadialToPoints( ps, centers, ys, params )
    np = size(ps,1)
    nc = size(centers,1)
    #println( "np, nc ", np," ", nc )
    for ip in 1:np
        pi = ps[ip,:]
        yi = 0.0::Float64
        for ic in 1:nc
            d  = centers[ic,:] - pi
            r2 = dot(d,d)            ::Float64
            #R2 = params[ic,1]
            #R2 = 0.125
            if r2<R2
                yi += R2func(r2/R2)
            end
        end
        ys[ip] = yi
    end
    return ys
end

function getForces( ps, centers, ys, params, fcs )
    np = size(ps,1)
    nc = size(centers,1)
    fc = zeros(Float64,2)
    dw = -1. /np
    for ic in 1:nc
        pci   = centers[ic,:]
        #param = params[ic,:]
        fill!(fc,0.0)
        for ip in 1:np
            d  = ps[ip,:] - pci
            r2 = dot(d,d)            ::Float64
            #R2 = params[ic,1]
            #R2 = 0.125
            if r2<R2
                ytot  = ys[ip]
                y,dy  = R2func_deriv(r2, 1/R2 )
                fr    = (1-ytot)*dy
                fc   += d*fr
            end
        end
        fcs[ic,:] = fc*dw
        println( fcs )
    end
    return fcs
end

function move( ps, fs, vs; dt=0.1, damp=0.9 )
    #vs[:,:] *= damp ;
    #vs[:,:] += fs*dt   ;
    #ps[:,:] += vs*dt ;
    ps[:,:] += fs*dt ;
end

# ================= Main

nimg = 200
xmax = 1.0
dx   = 2*xmax/nimg

ps = [ [(-xmax+ix*dx) (-xmax+iy*dx)] for iy=1:nimg,ix=1:nimg ]
ps = vcat(ps...)
#xs = ps[:,1]
#ys = ps[:,2]

centers = hexGrid( 2, 0.5 )

#centers = [ [-0.1 0.]; [0.1 0.] ]

println( "centers : ", typeof(centers)," ", centers  )

np = size(ps,1)
nc = size(centers,1)
#;println("ncenters = ", nc)

centers .+= (rand( Float64, nc,2 ) .- 0.5)*0.1

params = rand( Float64, nc,2 )   #;println("params ", params)
ys     = Array{Float64}(undef,np)
fcs    = Array{Float64}(undef,nc,2)
vcs    = zeros(Float64,nc,2)

#println( fcs )

#using Plots

#=
x = 1:0.1:100
y = 1:0.1:100
f(x, y) = begin
        (3x + y ^ 2) * abs(sin(x) + cos(y))
    end
X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))
Z = map(f, X, Y)
heatmap( Z, aspect_ratio=1)
#p1 = contour(x, y, f, fill=true)
#p2 = contour(x, y, Z)
#plot(p1, p2)
=#

Xs = reshape( ps[:,1], (nimg,nimg) )
Ys = reshape( ps[:,2], (nimg,nimg) )
Fs = reshape( ys,      (nimg,nimg) )
#println( Fs )
#p = contour( Xs, Ys, ys, fill=true)
#plot(p)
#@doc heatmap
#heatmap( Fs, aspect_ratio=1)

projectRadialToPoints( ps, centers, ys, params )
getForces( ps, centers, ys, params, fcs )

using PyPlot
imshow( Fs, origin="image", extent=(-1.,1.,-1.,1.) )
colorbar()
#x = range(0; stop=2*pi, length=1000);
#y = sin.(3 * x + 4 * cos.(2 * x));
#plot(x, y, color="red", linewidth=2.0, linestyle="--")
#plot(x, y, ".-" )

#plot(ps[:,1], ps[:,2], ".",  markersize=1.0 )
plot(centers[:,1], centers[:,2], "." )
quiver( centers[:,1], centers[:,2], fcs[:,1], fcs[:,2], )

for iter in 0:100
    close()
    projectRadialToPoints( ps, centers, ys, params )
    getForces( ps, centers, ys, params, fcs )
    #println( fcs )
    #println( centers )
    println( maximum( fcs ) )

    imshow( Fs, origin="image", extent=(-1.,1.,-1.,1.) )
    colorbar()
    plot(centers[:,1], centers[:,2], "." )
    quiver( centers[:,1], centers[:,2], fcs[:,1], fcs[:,2], )
    #axis('equal')
    sleep(0.2)
    move( centers, fcs, vcs, dt=0.5, damp=0.9 )
end

#plot( centers[:,1], centers[:,2] )
#show()

#scatter!( ps[:,1], ps[:,2], markersize = 1, markerstrokewidth = 0, aspect_ratio=:equal )
#scatter!( centers[:,1], centers[:,2], markersize = 2, markerstrokewidth = 0, aspect_ratio=:equal )

#plot( pts[:,1], pts[:,2] )
#scatter( pts[:,1], pts[:,2] )
#scatter( pts[:,1], pts[:,2], markerstrokestyle = :dot )

#pts = hexGrid( 2, 1.0 )
#scatter( pts[:,1], pts[:,2], markersize = 2, markerstrokewidth = 0, aspect_ratio=:equal )
#plot( pts[:,2] )

#show()
#close()

max_stars_repo_name: ProkopHapala/SimpleSimulationEngine

max_issues_repo_name: ProkopHapala/SimpleSimulationEngine

max_forks_repo_name: ProkopHapala/SimpleSimulationEngine


text:
@testset "minimal_inf_test" begin

    @test inf(emptyinterval()) === +Inf

    @test inf(interval(-Inf,+Inf)) === -Inf

    @test inf(interval(1.0,2.0)) === 1.0

    @test inf(interval(-3.0,-2.0)) === -3.0

    @test inf(interval(-Inf,2.0)) === -Inf

    @test inf(interval(-Inf,0.0)) === -Inf

    @test inf(interval(-Inf,-0.0)) === -Inf

    @test inf(interval(-2.0,Inf)) === -2.0

    @test inf(interval(0.0,Inf)) === -0.0

    @test inf(interval(-0.0,Inf)) === -0.0

    @test inf(interval(-0.0,0.0)) === -0.0

    @test inf(interval(0.0,-0.0)) === -0.0

    @test inf(interval(0.0,0.0)) === -0.0

    @test inf(interval(-0.0,-0.0)) === -0.0

end

@testset "minimal_inf_dec_test" begin

    @test isnan(inf(nai()))

    @test inf(DecoratedInterval(emptyinterval(), trv)) === +Inf

    @test inf(DecoratedInterval(interval(-Inf,+Inf), def)) === -Inf

    @test inf(DecoratedInterval(interval(1.0,2.0), com)) === 1.0

    @test inf(DecoratedInterval(interval(-3.0,-2.0), trv)) === -3.0

    @test inf(DecoratedInterval(interval(-Inf,2.0), dac)) === -Inf

    @test inf(DecoratedInterval(interval(-Inf,0.0), def)) === -Inf

    @test inf(DecoratedInterval(interval(-Inf,-0.0), trv)) === -Inf

    @test inf(DecoratedInterval(interval(-2.0,Inf), trv)) === -2.0

    @test inf(DecoratedInterval(interval(0.0,Inf), def)) === -0.0

    @test inf(DecoratedInterval(interval(-0.0,Inf), trv)) === -0.0

    @test inf(DecoratedInterval(interval(-0.0,0.0), dac)) === -0.0

    @test inf(DecoratedInterval(interval(0.0,-0.0), trv)) === -0.0

    @test inf(DecoratedInterval(interval(0.0,0.0), trv)) === -0.0

    @test inf(DecoratedInterval(interval(-0.0,-0.0), trv)) === -0.0

end

@testset "minimal_sup_test" begin

    @test sup(emptyinterval()) === -Inf

    @test sup(interval(-Inf,+Inf)) === +Inf

    @test sup(interval(1.0,2.0)) === 2.0

    @test sup(interval(-3.0,-2.0)) === -2.0

    @test sup(interval(-Inf,2.0)) === 2.0

    @test sup(interval(-Inf,0.0)) === 0.0

    @test sup(interval(-Inf,-0.0)) === 0.0

    @test sup(interval(-2.0,Inf)) === Inf

    @test sup(interval(0.0,Inf)) === Inf

    @test sup(interval(-0.0,Inf)) === Inf

    @test sup(interval(-0.0,0.0)) === 0.0

    @test sup(interval(0.0,-0.0)) === 0.0

    @test sup(interval(0.0,0.0)) === 0.0

    @test sup(interval(-0.0,-0.0)) === 0.0

end

@testset "minimal_sup_dec_test" begin

    @test isnan(sup(nai()))

    @test sup(DecoratedInterval(emptyinterval(), trv)) === -Inf

    @test sup(DecoratedInterval(interval(-Inf,+Inf), def)) === +Inf

    @test sup(DecoratedInterval(interval(1.0,2.0), com)) === 2.0

    @test sup(DecoratedInterval(interval(-3.0,-2.0), trv)) === -2.0

    @test sup(DecoratedInterval(interval(-Inf,2.0), dac)) === 2.0

    @test sup(DecoratedInterval(interval(-Inf,0.0), def)) === 0.0

    @test sup(DecoratedInterval(interval(-Inf,-0.0), trv)) === 0.0

    @test sup(DecoratedInterval(interval(-2.0,Inf), trv)) === Inf

    @test sup(DecoratedInterval(interval(0.0,Inf), def)) === Inf

    @test sup(DecoratedInterval(interval(-0.0,Inf), trv)) === Inf

    @test sup(DecoratedInterval(interval(-0.0,0.0), dac)) === +0.0

    @test sup(DecoratedInterval(interval(0.0,-0.0), trv)) === +0.0

    @test sup(DecoratedInterval(interval(0.0,0.0), trv)) === +0.0

    @test sup(DecoratedInterval(interval(-0.0,-0.0), trv)) === +0.0

end

@testset "minimal_mid_test" begin

    @test isnan(mid(emptyinterval()))

    @test mid(interval(-Inf,+Inf)) === 0.0

    @test mid(interval(-0x1.FFFFFFFFFFFFFp1023,+0x1.FFFFFFFFFFFFFp1023)) === 0.0

    @test mid(interval(0.0,2.0)) === 1.0

    @test mid(interval(2.0,2.0)) === 2.0

    @test mid(interval(-2.0,2.0)) === 0.0

    @test mid(interval(0.0,Inf)) === 0x1.FFFFFFFFFFFFFp1023

    @test mid(interval(-Inf,1.2)) === -0x1.FFFFFFFFFFFFFp1023

    @test_broken mid(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022)) === 0.0

    @test mid(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022)) === 0.0

    @test mid(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023)) === 0x1.7FFFFFFFFFFFFP+1023

    @test mid(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022)) === 0x0.0000000000002P-1022

end

@testset "minimal_mid_dec_test" begin

    @test isnan(mid(DecoratedInterval(emptyinterval(), trv)))

    @test isnan(mid(nai()))

    @test mid(DecoratedInterval(interval(-Inf,+Inf), def)) === 0.0

    @test mid(DecoratedInterval(interval(-0x1.FFFFFFFFFFFFFp1023,+0x1.FFFFFFFFFFFFFp1023), trv)) === 0.0

    @test mid(DecoratedInterval(interval(0.0,2.0), com)) === 1.0

    @test mid(DecoratedInterval(interval(2.0,2.0), dac)) === 2.0

    @test mid(DecoratedInterval(interval(-2.0,2.0), trv)) === 0.0

    @test mid(DecoratedInterval(interval(0.0,Inf), trv)) === 0x1.FFFFFFFFFFFFFp1023

    @test mid(DecoratedInterval(interval(-Inf,1.2), trv)) === -0x1.FFFFFFFFFFFFFp1023

    @test_broken mid(DecoratedInterval(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022), trv)) === 0.0

    @test mid(DecoratedInterval(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022), trv)) === 0.0

    @test mid(DecoratedInterval(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023), trv)) === 0x1.7FFFFFFFFFFFFP+1023

    @test mid(DecoratedInterval(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022), trv)) === 0x0.0000000000002P-1022

end

@testset "minimal_rad_test" begin

    @test radius(interval(0.0,2.0)) === 1.0

    @test radius(interval(2.0,2.0)) === 0.0

    @test isnan(radius(emptyinterval()))

    @test radius(interval(-Inf,+Inf)) === Inf

    @test radius(interval(0.0,Inf)) === Inf

    @test radius(interval(-Inf, 1.2)) === Inf

    @test radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022)) === 0x0.0000000000002P-1022

    @test radius(interval(0x0.0000000000001P-1022,0x0.0000000000002P-1022)) === 0x0.0000000000001P-1022

    @test radius(interval(0x1P+0,0x1.0000000000003P+0)) === 0x1P-51

end

@testset "minimal_rad_dec_test" begin

    @test radius(DecoratedInterval(interval(0.0,2.0), trv)) === 1.0

    @test radius(DecoratedInterval(interval(2.0,2.0), com)) === 0.0

    @test isnan(radius(DecoratedInterval(emptyinterval(), trv)))

    @test isnan(radius(nai()))

    @test radius(DecoratedInterval(interval(-Inf,+Inf), trv)) === Inf

    @test radius(DecoratedInterval(interval(0.0,Inf), def)) === Inf

    @test radius(DecoratedInterval(interval(-Inf, 1.2), trv)) === Inf

    @test radius(DecoratedInterval(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022), trv)) === 0x0.0000000000002P-1022

    @test radius(DecoratedInterval(interval(0x0.0000000000001P-1022,0x0.0000000000002P-1022), trv)) === 0x0.0000000000001P-1022

    @test radius(DecoratedInterval(interval(0x1P+0,0x1.0000000000003P+0), trv)) === 0x1P-51

end

@testset "minimal_mid_rad_test" begin

    @test isnan(midpoint_radius(emptyinterval())[1]) && isnan(midpoint_radius(emptyinterval())[2])

    @test midpoint_radius(interval(-Inf,Inf))[1] === 0.0 && midpoint_radius(interval(-Inf,Inf))[2] === Inf

    @test midpoint_radius(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023))[1] === 0.0 && midpoint_radius(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023))[2] === 0x1.FFFFFFFFFFFFFP+1023

    @test midpoint_radius(interval(0.0,2.0))[1] === 1.0 && midpoint_radius(interval(0.0,2.0))[2] === 1.0

    @test midpoint_radius(interval(2.0,2.0))[1] === 2.0 && midpoint_radius(interval(2.0,2.0))[2] === 0.0

    @test midpoint_radius(interval(-2.0,2.0))[1] === 0.0 && midpoint_radius(interval(-2.0,2.0))[2] === 2.0

    @test midpoint_radius(interval(0.0,Inf))[1] === 0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(interval(0.0,Inf))[2] === Inf

    @test midpoint_radius(interval(-Inf, 1.2))[1] === -0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(interval(-Inf, 1.2))[2] === Inf

    @test_broken midpoint_radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022))[1] === 0.0 && midpoint_radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022))[2] === 0x0.0000000000002P-1022

    @test midpoint_radius(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022))[1] === 0.0 && midpoint_radius(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022))[2] === 0x0.0000000000002P-1022

    @test midpoint_radius(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023))[1] === 0x1.7FFFFFFFFFFFFP+1023 && midpoint_radius(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023))[2] === 0x1.0p+1022

    @test midpoint_radius(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022))[1] === 0x0.0000000000002P-1022 && midpoint_radius(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022))[2] === 0x0.0000000000001P-1022

end

@testset "minimal_mid_rad_dec_test" begin

    @test isnan(midpoint_radius(emptyinterval())[1]) && isnan(midpoint_radius(emptyinterval())[2])

    @test isnan(midpoint_radius(nai())[1]) && isnan(midpoint_radius(nai())[2])

    @test midpoint_radius(DecoratedInterval(interval(-Inf,Inf), def))[1] === 0.0 && midpoint_radius(DecoratedInterval(interval(-Inf,Inf), def))[2] === Inf

    @test midpoint_radius(DecoratedInterval(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023), trv))[1] === 0.0 && midpoint_radius(DecoratedInterval(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023), trv))[2] === 0x1.FFFFFFFFFFFFFP+1023

    @test midpoint_radius(DecoratedInterval(interval(0.0,2.0), com))[1] === 1.0 && midpoint_radius(DecoratedInterval(interval(0.0,2.0), com))[2] === 1.0

    @test midpoint_radius(DecoratedInterval(interval(2.0,2.0), dac))[1] === 2.0 && midpoint_radius(DecoratedInterval(interval(2.0,2.0), dac))[2] === 0.0

    @test midpoint_radius(DecoratedInterval(interval(-2.0,2.0), trv))[1] === 0.0 && midpoint_radius(DecoratedInterval(interval(-2.0,2.0), trv))[2] === 2.0

    @test midpoint_radius(DecoratedInterval(interval(0.0,Inf), trv))[1] === 0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(DecoratedInterval(interval(0.0,Inf), trv))[2] === Inf

    @test midpoint_radius(DecoratedInterval(interval(-Inf, 1.2), trv))[1] === -0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(DecoratedInterval(interval(-Inf, 1.2), trv))[2] === Inf

    @test_broken midpoint_radius(DecoratedInterval(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022), trv))[1] === 0.0 && midpoint_radius(DecoratedInterval(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022), trv))[2] === 0x0.0000000000002P-1022

    @test midpoint_radius(DecoratedInterval(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022), trv))[1] === 0.0 && midpoint_radius(DecoratedInterval(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022), trv))[2] === 0x0.0000000000002P-1022

    @test midpoint_radius(DecoratedInterval(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023), trv))[1] === 0x1.7FFFFFFFFFFFFP+1023 && midpoint_radius(DecoratedInterval(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023), trv))[2] === 0x1.0p+1022

    @test midpoint_radius(DecoratedInterval(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022), trv))[1] === 0x0.0000000000002P-1022 && midpoint_radius(DecoratedInterval(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022), trv))[2] === 0x0.0000000000001P-1022

end

@testset "minimal_wid_test" begin

    @test diam(interval(2.0,2.0)) === 0.0

    @test diam(interval(1.0,2.0)) === 1.0

    @test diam(interval(1.0,Inf)) === Inf

    @test diam(interval(-Inf,2.0)) === Inf

    @test diam(interval(-Inf,+Inf)) === Inf

    @test isnan(diam(emptyinterval()))

    @test diam(interval(0x1P+0,0x1.0000000000001P+0)) === 0x1P-52

    @test diam(interval(0x1P-1022,0x1.0000000000001P-1022)) === 0x0.0000000000001P-1022

end

@testset "minimal_wid_dec_test" begin

    @test diam(DecoratedInterval(interval(2.0,2.0), com)) === 0.0

    @test diam(DecoratedInterval(interval(1.0,2.0), trv)) === 1.0

    @test diam(DecoratedInterval(interval(1.0,Inf), trv)) === Inf

    @test diam(DecoratedInterval(interval(-Inf,2.0), def)) === Inf

    @test diam(DecoratedInterval(interval(-Inf,+Inf), trv)) === Inf

    @test isnan(diam(DecoratedInterval(emptyinterval(), trv)))

    @test isnan(diam(nai()))

    @test diam(DecoratedInterval(interval(0x1P+0,0x1.0000000000001P+0), trv)) === 0x1P-52

    @test diam(DecoratedInterval(interval(0x1P-1022,0x1.0000000000001P-1022), trv)) === 0x0.0000000000001P-1022

end

@testset "minimal_mag_test" begin

    @test mag(interval(1.0,2.0)) === 2.0

    @test mag(interval(-4.0,2.0)) === 4.0

    @test mag(interval(-Inf,2.0)) === Inf

    @test mag(interval(1.0,Inf)) === Inf

    @test mag(interval(-Inf,+Inf)) === Inf

    @test isnan(mag(emptyinterval()))

    @test mag(interval(-0.0,0.0)) === 0.0

    @test mag(interval(-0.0,-0.0)) === 0.0

end

@testset "minimal_mag_dec_test" begin

    @test mag(DecoratedInterval(interval(1.0,2.0), com)) === 2.0

    @test mag(DecoratedInterval(interval(-4.0,2.0), trv)) === 4.0

    @test mag(DecoratedInterval(interval(-Inf,2.0), trv)) === Inf

    @test mag(DecoratedInterval(interval(1.0,Inf), def)) === Inf

    @test mag(DecoratedInterval(interval(-Inf,+Inf), trv)) === Inf

    @test isnan(mag(DecoratedInterval(emptyinterval(), trv)))

    @test isnan(mag(nai()))

    @test mag(DecoratedInterval(interval(-0.0,0.0), trv)) === 0.0

    @test mag(DecoratedInterval(interval(-0.0,-0.0), trv)) === 0.0

end

@testset "minimal_mig_test" begin

    @test mig(interval(1.0,2.0)) === 1.0

    @test mig(interval(-4.0,2.0)) === 0.0

    @test mig(interval(-4.0,-2.0)) === 2.0

    @test mig(interval(-Inf,2.0)) === 0.0

    @test mig(interval(-Inf,-2.0)) === 2.0

    @test mig(interval(-1.0,Inf)) === 0.0

    @test mig(interval(1.0,Inf)) === 1.0

    @test mig(interval(-Inf,+Inf)) === 0.0

    @test isnan(mig(emptyinterval()))

    @test mig(interval(-0.0,0.0)) === 0.0

    @test mig(interval(-0.0,-0.0)) === 0.0

end

@testset "minimal_mig_dec_test" begin

    @test mig(DecoratedInterval(interval(1.0,2.0), com)) === 1.0

    @test mig(DecoratedInterval(interval(-4.0,2.0), trv)) === 0.0

    @test mig(DecoratedInterval(interval(-4.0,-2.0), trv)) === 2.0

    @test mig(DecoratedInterval(interval(-Inf,2.0), def)) === 0.0

    @test mig(DecoratedInterval(interval(-Inf,-2.0), trv)) === 2.0

    @test mig(DecoratedInterval(interval(-1.0,Inf), trv)) === 0.0

    @test mig(DecoratedInterval(interval(1.0,Inf), trv)) === 1.0

    @test mig(DecoratedInterval(interval(-Inf,+Inf), trv)) === 0.0

    @test isnan(mig(DecoratedInterval(emptyinterval(), trv)))

    @test isnan(mig(nai()))

    @test mig(DecoratedInterval(interval(-0.0,0.0), trv)) === 0.0

    @test mig(DecoratedInterval(interval(-0.0,-0.0), trv)) === 0.0

end


text:
using PaddleChainRules.Paddle: PaddleModuleWrapper, PaddleFCNet, paddle, to_paddle_tensor, dlpack, pyto_dlpack, pyfrom_dlpack, ispysetup

using Test
using Zygote
import Random
using PyCall
using CUDA
using DLPack
using Functors: fmap

if !ispysetup[]
    return
end

if CUDA.functional()
    @info "using gpu"
else
    @info "using cpu"
end

function compare_grad_wrt_params(model, modelwrap, inputs...)
    pyinputs = fmap(x -> to_paddle_tensor(x,false), inputs)
    paddleoutputs = model(pyinputs...)
    paddlegrad = map(x-> (x.numpy()), paddle.fluid.dygraph.grad(paddleoutputs, model.parameters()))

    grad,  = Zygote.gradient(m->sum(m(inputs...)), modelwrap)
    @test length(paddlegrad) == length(grad.params)
    for i in 1:length(grad.params)
        @test isapprox(sum(paddlegrad[i]), sum(grad.params[i]))
    end
    @test length(grad.params) == length(modelwrap.params)
    @test grad.params[1] !== nothing
    @test grad.params[2] !== nothing
    @test size(grad.params[1]) == size(modelwrap.params[1])
    @test size(grad.params[2]) == size(modelwrap.params[2])
end

function compare_grad_wrt_inputs(model, modelwrap, inputs)
    pyinputs = to_paddle_tensor(inputs,false)
    paddleoutputs = model(pyinputs)
    paddlegrad = map(x-> (x.numpy()), paddle.fluid.dygraph.grad(paddleoutputs, pyinputs))[1]
    grad, = Zygote.gradient(z->sum(modelwrap(z)), inputs)
    @test size(grad) == size(inputs)
    @test length(paddlegrad) == length(grad)
    @test isapprox(sum(paddlegrad), sum(grad))
end

batchsize = 1
indim = 3
outdim = 2
hiddendim = 4

@testset "mlp" begin
    mlp = paddle.nn.Sequential(
        paddle.nn.Linear(indim, hiddendim),
        paddle.nn.Sigmoid(),
        paddle.nn.Linear(hiddendim, outdim)
    )
    mlpwrap = PaddleModuleWrapper(mlp)
    if CUDA.functional()
        mlpwrap = fmap(CUDA.cu, mlpwrap)
    end
    x = randn(Float32, indim, batchsize)
    if CUDA.functional()
        x = CUDA.cu(x)
    end
    y = mlpwrap(x)
    @test size(y) == (outdim, batchsize)
    compare_grad_wrt_params(mlp, mlpwrap, x)
    compare_grad_wrt_inputs(mlp, mlpwrap, x)
end

@testset "General Net" begin
    NN = paddle.nn.Linear(indim, outdim)

    NNwrap = PaddleModuleWrapper(NN)
    if CUDA.functional()
        NNwrap = fmap(CUDA.cu, NNwrap)
    end
    x = randn(Float32, indim, batchsize)
    if CUDA.functional()
        x = CUDA.cu(x)
    end
    y = NNwrap(x)
    @test size(y) == (outdim, batchsize)
    compare_grad_wrt_params(NN, NNwrap, x)
    compare_grad_wrt_inputs(NN, NNwrap, x)
end

@testset "FC Net" begin
    NN = paddle.nn.Sequential(
        paddle.nn.Linear(indim, hiddendim),
        paddle.nn.Sigmoid(),
        paddle.nn.Linear(hiddendim, outdim)
    )

    NNwrap = PaddleFCNet(indim, outdim, 2, hiddendim; dtype="float32", activation="sigmoid")

    # set same parameters
    pyparams = NN.parameters()
    jlparams = map(x->DLPack.wrap(x, pyto_dlpack), pyparams)
    map((p1,p2)->copy!(p1, p2), NNwrap.params, jlparams)

    if CUDA.functional()
        NNwrap = fmap(CUDA.cu, NNwrap)
    end
    x = randn(Float32, indim, batchsize)
    if CUDA.functional()
        x = CUDA.cu(x)
    end
    y = NNwrap(x)
    @test size(y) == (outdim, batchsize)
    compare_grad_wrt_params(NN, NNwrap, x)
    compare_grad_wrt_inputs(NN, NNwrap, x)

end

@testset "Error throw" begin
    NN = paddle.nn.Sequential(
        paddle.nn.Conv2D(1,2,5)
    )
    @test_throws ErrorException("Unsupported layer type") PaddleModuleWrapper(NN)
    @test_throws ErrorException("Unsupported activation type") PaddleFCNet(indim, outdim, 2, hiddendim; dtype="float32", activation="celu")
    @test_throws ErrorException("Unsupported data type, use \"float16\", \"float32\" or \"float64\"") PaddleFCNet(indim, outdim, 2, hiddendim; dtype="int8", activation="sigmoid")
end

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
#################### Dependent ####################

const depfxargs = [(:model, MCPhylo.Model)]

#################### Base Methods ####################

dims(d::AbstractDependent) = size(d)

function names(d::AbstractDependent)
  names(d, d.symbol)
end
"""
    setmonitor!(d::AbstractDependent, monitor::Bool)

Specify node elements to be included in monitored MCMC sampler output.

Returns d with its monitor field updated to reflect the specified monitoring.

* `d` : node whose elements contain sampled MCMC values.

* `monitor` : boolean indicating whether all elements are monitored.

"""
function setmonitor!(d::AbstractDependent, monitor::Bool)
  value = monitor ? Int[0] : Int[]
  setmonitor!(d, value)
end
"""
    setmonitor!(d::AbstractDependent, monitor::Vector{Int})

Specify node elements to be included in monitored MCMC sampler output.

Returns d with its monitor field updated to reflect the specified monitoring.

* `d` : node whose elements contain sampled MCMC values.

* `monitor` : vector of element-wise indices of elements to monitor.
"""
function setmonitor!(d::AbstractDependent, monitor::Vector{Int})
  values = monitor
  if !isempty(monitor)
    n = isa(d, AbstractTreeStochastic) ? length(unlist_tree(d)) : length(unlist(d))
    if n > 0
      if monitor[1] == 0
        values = collect(1:n)
      elseif minimum(monitor) < 1 || maximum(monitor) > n
        throw(BoundsError())
      end
    end
  end
  d.monitor = values
  d
end
#
"""
    setmonitor!(d::AbstractTreeStochastic, monitor::Bool)

Specify node elements to be included in monitored MCMC sampler output.

Returns `d` with its `monitor` field updated to reflect the specified monitoring.

* `d` : node whose elements contain sampled MCMC values.

* `monitor` : boolean indicating whether all elements are monitored.
"""
function setmonitor!(d::AbstractTreeStochastic, monitor::Bool)
  value = monitor ? Int[0] : Int[]
  setmonitor!(d, value)
end

"""
    setmonitor!(d::AbstractTreeStochastic, monitor::Vector{Int})

Specify node elements to be included in monitored MCMC sampler output.

Returns `d` with its `monitor` field updated to reflect the specified monitoring.

* `d` : node whose elements contain sampled MCMC values.

* `monitor` : vector of element-wise indices of elements to monitor.
"""
function setmonitor!(d::AbstractTreeStochastic, monitor::Vector{Int})
  values = monitor
  if !isempty(monitor)
    n = length(unlist_tree(d))
    if n > 0
      if monitor[1] == 0
        values = collect(1:n)
      elseif minimum(monitor) < 1 || maximum(monitor) > n
        throw(BoundsError())
      end
    end
  end
  d.monitor = values
  d
end

#################### Distribution Fallbacks ####################

unlist(d::AbstractDependent, transform::Bool=false) =
  unlist(d, d.value, transform)

unlist(d::AbstractDependent, x::Real, transform::Bool=false) = [x]

unlist(d::AbstractDependent, x::AbstractArray, transform::Bool=false) = vec(x)

relist(d::AbstractDependent, x::AbstractArray, transform::Bool=false) =
  relistlength(d, x, transform)[1]

logpdf(d::AbstractDependent, transform::Bool=false) = 0.0

logpdf(d::AbstractDependent, x, transform::Bool=false) = 0.0

gradlogpdf(d::AbstractDependent, x, transform::Bool=false) = 0.0

gradlogpdf(d::AbstractDependent, transform::Bool=false) = 0.0

#################### Logical ####################

@promote_scalarvariate ScalarLogical

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

"""
    Logical(f::Function, monitor::Union{Bool, Vector{Int}}=true)

Constructor for a Logical model node. This function assumes a the output of the
logical operation to be scalar.

* `f`: Function specifying the deterministic operation performed on its arguments. These arguments are other nodes of the model.

* `monitor` : Indicates whether the results should be monitored, i.e.saved.
"""
function Logical(f::Function, monitor::Union{Bool, Vector{Int}}=true)
  value = Float64(NaN)
  fx, src = modelfxsrc(depfxargs, f)
  l = ScalarLogical(value, :nothing, Int[], fx, src, Symbol[])
  setmonitor!(l, monitor)
end

Logical(f::Function, d::Integer, args...) = Logical(d, f, args...)

"""
    Logical(d::Integer, f::Function, monitor::Union{Bool, Vector{Int}}=true)

Constructor for a Logical model node.

* `d` : Specifies the dimension of the output.

* `f` : Specifies the deterministic operation performed on its arguments. These arguments are other nodes of the model.

* `monitor` : Indicates whether the results should be monitored, i.e.saved.
"""
function Logical(d::Integer, f::Function,
                 monitor::Union{Bool, Vector{Int}}=true)
  value = Array{Float64}(undef, fill(0, d)...)
  fx, src = modelfxsrc(depfxargs, f)
  l = ArrayLogical(value, :nothing, Int[], fx, src, Symbol[])
  setmonitor!(l, monitor)
end

Logical(f::Function, d::T, args...)  where T<:GeneralNode = Logical(d, f, args...)

"""
    Logical(d::T, f::Function, monitor::Union{Bool, Vector{Int}}=true) where T<:GeneralNode

Constructor for a Logical model node, which can hold a Node structure, i.e. a tree.

* `f` is a function specifying the deterministic operation performed on its arguments. These arguments are other nodes of the model.

* `monitor` indicates whether the results should be monitored, i.e. saved.
"""
function Logical(d::T, f::Function,
                 monitor::Union{Bool, Vector{Int}}=true) where T<:GeneralNode
  value = T()
  fx, src = modelfxsrc(depfxargs, f)
  l = TreeLogical(value, :nothing, Int[], fx, src, Symbol[])
  setmonitor!(l, monitor)
end

ScalarLogical(x::T) where T <: Real = x

#################### Updating ####################
"""
    setinits!(l::AbstractLogical, m::Model, ::Any=nothing)

Set initial values for a logical node.

Returns the result of a call to setmonitor!(l, l.monitor) or setmonitor!(d, d.monitor).

* `l` : logical node to which to assign initial values.

* `m` : model containing the node.
"""
function setinits!(l::AbstractLogical, m::Model, ::Any=nothing)
  l.value = l.eval(m)
  setmonitor!(l, l.monitor)
end
"""
    setinits!(d::TreeLogical, m::Model, x::T) where {T<:Node}

Set initial values for a logical node.

Returns the result of a call to `setmonitor!(l, l.monitor)` or `setmonitor!(d, d.monitor)`.

* `l`  : logical node to which to assign initial values.

* `m`  : model containing the node.
"""
function setinits!(d::TreeLogical, m::Model, x::T) where {T<:Node}
    d.value = d.eval(m)
    setmonitor!(d, d.monitor)
end # function

"""
    update!(l::AbstractLogical, m::Model)

Update the values of a logical node according to its relationship with others in a model.

Returns the node with its values updated.

* `l` : logical node to update.

* `m` : model containing the node.
"""
function update!(l::AbstractLogical, m::Model)
  l.value = l.eval(m)
  l
end

#################### Distribution Methods ####################

relistlength(d::ScalarLogical, x::AbstractArray, transform::Bool=false) =
  (x[1], 1)

function relistlength(d::ArrayLogical, x::AbstractArray, transform::Bool=false)
  n = length(d)
  value = reshape(x[1:n], size(d))
  (value, n)
end

#################### Stochastic ####################

#################### Base Methods ####################

@promote_scalarvariate ScalarStochastic

function showall(io::IO, s::AbstractStochastic)
  show(io, s)
  print(io, "\n\nDistribution:\n")
  show(io, s.distr)
  print(io, "\nFunction:\n")
  show(io, "text/plain", first(code_typed(s.eval)))
  print(io, "\n\nSource Nodes:\n")
  show(io, s.sources)
  print(io, "\n\nTarget Nodes:\n")
  show(io, s.targets)
end

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

"""
    Stochastic(f::Function, monitor::Union{Bool, Vector{Int}}=true)

Constructor for a Stochastic model node. This function assumes a the output of the
logical operation to be scalar.

* `f` : Specifies the distributional relationship between the arguments and the node. These arguments are other nodes of the model.

* `monitor` : Indicates whether the results should be monitored, i.e. saved.
"""
function Stochastic(f::Function, monitor::Union{Bool, Vector{Int}}=true)
  value = Float64(NaN)
  fx, src = modelfxsrc(depfxargs, f)
  s = ScalarStochastic(value, :nothing, Int[], fx, src, Symbol[],
                       NullUnivariateDistribution())
  setmonitor!(s, monitor)
end

Stochastic(f::Function, d::Integer, args...) = Stochastic(d, f, args...)

function Stochastic_cu(d::Integer, f::Function,
                    monitor::Union{Bool, Vector{Int}}=true)

  value = CuArray{Float64}(undef, fill(0, d)...)
  fx, src = modelfxsrc(depfxargs, f)
  s = ArrayStochastic(value, :nothing, Int[], fx, src, Symbol[],
                      NullUnivariateDistribution())
  setmonitor!(s, monitor)
end

function Stochastic_ncu(d::Integer, f::Function,
                    monitor::Union{Bool, Vector{Int}}=true)

  value = Array{Float64}(undef, fill(0, d)...)

  fx, src = modelfxsrc(depfxargs, f)
  s = ArrayStochastic(value, :nothing, Int[], fx, src, Symbol[],
                      NullUnivariateDistribution())
  setmonitor!(s, monitor)
end

"""
    Stochastic(d::Integer, f::Function, monitor::Union{Bool,
        Vector{Int}}=true, cuda::Bool=false)

Constructor for a Stochastic model node.

* `d` : Specifies dimensions of the output.

* `f` : Specifies the distributional relationship between the arguments and the node. These arguments are other nodes of the model.

* `monitor` : indicates whether the results should be monitored, i.e. saved.

* `cuda` : indicates whether the function supports cuda functionality and the data is
in the respective format. THIS FEATURE IS NOT FULLY SUPPORTED.
"""
function Stochastic(d::Integer, f::Function,
                    monitor::Union{Bool, Vector{Int}}=true, cuda::Bool=false)
  if cuda
     return Stochastic_cu(d, f, monitor)
  else
     return Stochastic_ncu(d, f, monitor)
  end

end

Stochastic(f::Function, d::T, args...)  where T<:GeneralNode = Stochastic(d, f, args...)

"""

    Stochastic(d::N, f::Function, monitor::Union{Bool, Vector{Int}}=true) where N<:GeneralNode

Constructor for a Stochastic model node, which can hold a Node structure, i.e. a tree.

* `f` : Specifies the distributional reslationship between the arguments and the node.
These arguments are other nodes of the model.

* `monitor` : Indicates whether the results should be monitored, i.e. saved.
"""
function Stochastic(d::N, f::Function, monitor::Union{Bool, Vector{Int}}=true) where N<:GeneralNode
    value = Node()
    fx, src = modelfxsrc(depfxargs, f)
    s = TreeStochastic(value, :nothing, Int[], fx, src, Symbol[],
                      NullUnivariateDistribution())
    setmonitor!(s, monitor)
end

ScalarStochastic(x::T) where T <: Real = x

#################### Updating ####################
"""
    setinits!(s::ScalarStochastic, m::Model, x::Real)

Set initial values for a stochastic node.

Returns the node with its assigned initial values.

* `s` : stochastic node to which to assign initial values.

* `m` : model containing the node.

* `x` : values to assign to the node.
"""
function setinits!(s::ScalarStochastic, m::Model, x::Real)
  s.value = convert(Float64, x)
  s.distr = s.eval(m)
  setmonitor!(s, s.monitor)
end
"""
    setinits!(s::ArrayStochastic, m::Model, x::DenseArray)

Set initial values for a stochastic node.

Returns the node with its assigned initial values.

* `s` : stochastic node to which to assign initial values.

* `m` : model containing the node.

* `x` : values to assign to the node.
"""
function setinits!(s::ArrayStochastic, m::Model, x::DenseArray)
  s.value = convert(typeof(s.value), copy(x))
  s.distr = s.eval(m)
  if isa(s.distr, PhylogeneticDistribution)
    distrdims = dims(s.distr)
    for (ind, di) in enumerate(dims(s))
      if ind != 2
        di != distrdims[ind] && throw(DimensionMismatch("incompatible distribution for stochastic node"))
      end
    end
  elseif !isa(s.distr, UnivariateDistribution) && dims(s) != dims(s.distr)
    throw(DimensionMismatch("incompatible distribution for stochastic node"))
  end
  setmonitor!(s, s.monitor)
end

function setinits!(s::AbstractStochastic, m::Model, x)
  throw(ArgumentError("incompatible initial value for node : $(s.symbol)"))
end
"""
    setinits!(d::TreeStochastic, m::Model, x::T) where {T<:GeneralNode}

Set initial values for a stochastic node.

Returns the node with its assigned initial values.

* `s` or `d` : stochastic node to which to assign initial values.

* `m` : model containing the node.

* `x` : values to assign to the node.
"""
function setinits!(d::TreeStochastic, m::Model, x::T) where {T<:GeneralNode}
    d.value = deepcopy(x)
    d.distr = d.eval(m)
    insupport(d.distr, x) || throw(ArgumentError("The supplied tree does not match the topological tree constraints."))
    setmonitor!(d, d.monitor)
end # function

"""
    update!(s::AbstractStochastic, m::Model)

Update the values of a stochastic node according to its relationship with others in a model.

Returns the node with its values updated.

* `s` : stochastic node to update.

* `m` : model containing the node.
"""
function update!(s::AbstractStochastic, m::Model)
  s.distr = s.eval(m)
  s
end

#################### Distribution Methods ####################

function unlist(s::AbstractStochastic, transform::Bool=false)
  unlist(s, s.value, transform)
end

function unlist(s::AbstractStochastic, x::Real, transform::Bool=false)
  unlist(s, [x], transform)
end

function unlist(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  transform ? unlist_sub(s.distr, link_sub(s.distr, x)) :
              unlist_sub(s.distr, x)
end

function relist(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  relistlength(s, x, transform)[1]
end

function relistlength(s::AbstractStochastic, x::AbstractArray,
                      transform::Bool=false)
  value, n = relistlength_sub(s.distr, s, x)
  (transform ? invlink_sub(s.distr, value) : value, n)
end

function relistlength(s::AbstractTreeStochastic, x::AbstractArray,
                      transform::Bool=false)
  value, n = relistlength_sub(s.distr, s, x)

  (transform ? invlink_sub(s.distr, value) : value, n)
end

function relistlength(s::AbstractTreeStochastic, x::N,
                      transform::Bool=false) where N<:GeneralNode
  value, n = relistlength_sub(s.distr, s, x)
  (transform ? invlink_sub(s.distr, value) : value, n)
end

function logpdf(s::AbstractTreeStochastic, transform::Bool=false)
  logpdf(s, s.value, transform)
end

function logpdf(s::AbstractStochastic, transform::Bool=false)
  logpdf(s, s.value, transform)
end

function conditional_likelihood(s::AbstractStochastic, args...)
  conditional_likelihood(s, s.value, args...)
end

function conditional_likelihood(s::AbstractStochastic, x::AbstractArray, args...)
  logcond(s.distr, x, args...)
end

function pseudologpdf(s::AbstractStochastic, x::Real, transform::Bool=false)
  logpdf(s, x, transform)
end

function pseudologpdf(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  logpdf(s, x, transform)
end

function rand(s::AbstractStochastic, x::Int64)
  rand(s.distr, x)
end

function gradlogpdf(s::AbstractStochastic)
  gradlogpdf(s, s.value)
end

function gradlogpdf(s::AbstractLogical)
  gradlogpdf(s, s.value)
end

function gradlogpdf(s::AbstractTreeStochastic, x::N, transform::Bool=false) where N<:GeneralNode
  gradlogpdf(s.distr, x)
end

function gradlogpdf(s::AbstractStochastic, x::AbstractArray)
  gradlogpdf_sub(s.distr, x)
end

function logpdf(s::AbstractStochastic, x::AbstractArray, transform::Bool=false)
  logpdf_sub(s.distr, x, transform)
end

function logpdf(s::AbstractStochastic, x::Real, transform::Bool=false)
  logpdf_sub(s.distr, x, transform)
end

function logpdf(s::AbstractTreeStochastic, x::N, transform::Bool=false) where N<:GeneralNode
  logpdf_sub(s.distr, x, transform)
end

rand(s::AbstractStochastic) = rand_sub(s.distr, s.value)


text:
# This file contains code that helps to query from the registry to determine the format

"""
`unknown(f)` returns true if the format of `f` is unknown.
"""
unknown(::Type{format"UNKNOWN"}) = true
unknown(::Type{DataFormat{sym}}) where {sym} = false

unknown(::File{F}) where {F} = unknown(F)
unknown(::Stream{F}) where {F} = unknown(F)

const unknown_df = DataFormat{:UNKNOWN}

"""
`info(fmt)` returns the magic bytes/extension information for
`DataFormat` `fmt`.
"""
info(::Type{DataFormat{sym}}) where {sym} = sym2info[sym]

"`magic(fmt)` returns the magic bytes of format `fmt`"
magic(fmt::Type{<:DataFormat})= UInt8[info(fmt)[1]...]

"""
`skipmagic(s)` sets the position of `Stream` `s` to be just after the magic bytes.
For a plain IO object, you can use `skipmagic(io, fmt)`.
"""
skipmagic(s::Stream{F}) where {F} = (skipmagic(stream(s), F); s)
function skipmagic(io, fmt::Type{DataFormat{sym}}) where sym
    magic, _ = sym2info[sym]
    skipmagic(io, magic)
    nothing
end
skipmagic(io, magic::Function) = nothing
skipmagic(io, magic::NTuple{N,UInt8}) where {N} = seek(io, length(magic))
function skipmagic(io, magic::Tuple)
    lengths = map(length, magic)
    all(x-> lengths[1] == x, lengths) && return seek(io, lengths[1]) # it doesn't matter what magic bytes get skipped as they all have the same length
    magic = [magic...]
    sort!(magic, lt = (a,b)-> length(a) >= length(b)) # start with longest first, to avoid overlapping magic bytes
    seekend(io)
    len = position(io)
    seekstart(io)
    filter!(x-> length(x) <= len, magic) # throw out magic bytes that are longer than IO
    tmp = read(io, length(first(magic))) # now, first is both the longest and guaranteed to fit into io, so we can just read the bytes
    for m in magic
        if magic_equal(m, tmp)
            seek(io, length(m))
            return nothing
        end
    end
    error("tried to skip magic bytes of an IO that does not contain the magic bytes of the format. IO: $io")
end
function magic_equal(magic, buffer)
    for (i,elem) in enumerate(magic)
        buffer[i] != elem && return false
    end
    true
end

"""
`query(filename)` returns a `File` object with information about the
format inferred from the file's extension and/or magic bytes.
"""
function query(filename::AbstractString)
    _, ext = splitext(filename)
    if haskey(ext2sym, ext)
        sym = ext2sym[ext]
        no_magic = !hasmagic(sym)
        if lensym(sym) == 1 && (no_magic || !isfile(filename)) # we only found one candidate and there is no magic bytes, or no file, trust the extension
            return File{DataFormat{sym}}(filename)
        elseif !isfile(filename) && lensym(sym) > 1
            return File{DataFormat{sym[1]}}(filename)
        end
        if no_magic && !hasfunction(sym)
            error("Some formats with extension ", ext, " have no magic bytes; use `File{format\"FMT\"}(filename)` to resolve the ambiguity.")
        end
    end
    !isfile(filename) && return File{unknown_df}(filename) # (no extension || no magic byte) && no file
    # Otherwise, check the magic bytes
    file!(query(open(filename), abspath(filename)))
end

lensym(s::Symbol) = 1
lensym(v::Vector) = length(v)

hasmagic(s::Symbol) = hasmagic(sym2info[s][1])
hasmagic(v::Vector) = any(hasmagic, v)

hasmagic(t::Tuple) = !isempty(t)
hasmagic(::Any) = false   # for when magic is a function

hasfunction(s::Symbol) = hasfunction(sym2info[s][1])
hasfunction(v::Vector) = any(hasfunction, v)
hasfunction(s::Any) = true #has function
hasfunction(s::Tuple) = false #has magic

"""
`query(io, [filename])` returns a `Stream` object with information about the
format inferred from the magic bytes.
"""
query(io::IO, filename) = query(io, String(filename))

function query(io::IO, filename::Union{Nothing, String} = nothing)
    magic = Vector{UInt8}()
    pos = position(io)
    for p in magic_list
        m = first(p)
        length(m) == 0 && continue
        while length(m) > length(magic)
            if eof(io)
                seek(io, pos)
                return Stream{unknown_df, typeof(io)}(io, filename)
            end
            push!(magic, read(io, UInt8))
        end
        if iter_eq(magic, m)
            seek(io, pos)
            return Stream{DataFormat{last(p)},typeof(io)}(io, filename)
        end
    end
    if seekable(io)
        for p in magic_func
            seek(io, pos)
            f = first(p)
            try
                if f(io)
                    return Stream{DataFormat{last(p)},typeof(io)}(seek(io, pos), filename)
                end
            catch e
                println("There was an error in magick function $f")
                println("Please open an issue at FileIO.jl. Error:")
                println(e)
            end
        end
        seek(io, pos)
    end
    Stream{unknown_df,typeof(io)}(io, filename)
end

seekable(io::IOBuffer) = io.seekable
seekable(::IOStream) = true
seekable(::Any) = false

function iter_eq(A, B)
    length(A) == length(B) || return false
    i,j = 1,1
    for _=1:length(A)
        a=A[i]; b=B[j]
        a == b && (i+=1; j+=1; continue)
        a == UInt32('\r') && (i+=1; continue) # this seems like the shadiest solution to deal with windows \r\n
        b == UInt32('\r') && (j+=1; continue)
        return false #now both must be unequal, and no \r windows excemption any more
    end
    true
end


text:
module OmegaRejectionSample
import Base.Threads

import ..OmegaCore
using ..Util, ..Sample, ..Conditioning, ..TrackError, ..RNG
const OC = OmegaCore
# import ..OmegaCore: randsample

export RejectionSample

"Rejection Sampling  Algorithm"
struct RejectionSampleAlg end
const RejectionSample = RejectionSampleAlg()

OmegaCore.defrandalg(args...) = RejectionSample

function omegarandsample1(rng,
                          ΩT::Type{OT},
                          y,
                          alg::RejectionSampleAlg) where OT
  @label restart
  ω = ΩT()
  ω_ = tagrng(ω, rng)
  !y(ω_) && @goto restart
  ω
end

"`n` samples from ω::ΩT such that `y(ω)` is true"
function OC.omegarandsample(rng,
                       ΩT::Type{OT},
                       y,
                       n,
                       alg::RejectionSampleAlg) where OT
  hack = OT()
  ωsamples = Vector{typeof(hack)}(undef, n)  # Fixme
  accepted = 0
  i = 1
  # rngs = OC.duplicaterng(rng, Threads.nthreads())
  # Threads.@threads for i = 1:n
  #   rng = rngs[Threads.threadid()]
  #   @inbounds ωsamples[i] = omegarandsample1(rng, ΩT, y, n, alg)
  # end
  for i = 1:n
    # rng = rngs[Threads.threadid()]
    @inbounds ωsamples[i] = omegarandsample1(rng, ΩT, y, alg)
  end

  ωsamples
end

function OC.randsample(rng,
                       ΩT::Type{OT},
                       x,
                       n,
                       alg::RejectionSampleAlg) where {OT}
  # introduce conditions
  # y = OC.mem(OC.indomain(x))
  y = condvar(x, Bool)
  ωsamples = OC.omegarandsample(rng, ΩT, y, n, alg)
  # map(OC.mem(x), ωsamples)
  map(x, ωsamples)
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect