text:
module test_julia_assign

using Test

(j, j) = (1, 2)
@test j == 2

tup = 1, 2
@test tup == (1, 2)

d, = 1, 2
@test d == 1

end # module test_julia_assign


text:
struct PointSet{Dim,T,V<:AbstractVector{Point{Dim,T}}}
    points::V
end

(==)(x::PointSet, y::PointSet) = x.points == y.points

(transf::Transform)(set::PointSet) = PointSet(map(transf, set.points))

function centroid(set::PointSet)
    points = set.points
    sum(points) / length(points)
end

function boundingelement(set::PointSet{Dim}) where {Dim}
    c = centroid(set)
    coords = map(1:Dim) do i
        maximum(getindex.(set.points, i))
    end
    scale = Scaling(coords - c)
    Translated(Scaled(HyperCube(1.), scale), Translation(c))
end

@generated function PointSet(::Type{<:HyperCube}, P::Type{Point{Dim,T}}) where {Dim,T}
    idxs = CartesianIndices(ntuple(i -> -1:2:1, Dim))
    tuples = getproperty.(idxs, :I)
    :(PointSet(SVector{$(2^Dim),Point{$Dim,$T}}($(tuples...))))
end

PointSet(obj::HyperCube, P) = UniformScaling(obj.radius)(PointSet(typeof(obj), P))
PointSet(obj::Transformed, P) = UniformScaling(obj.radius)(obj.transf(PointSet(obj.obj, P)))

function sort_nearest(set::PointSet, point, dist=HyperSphere(0.))
    f = Translated(dist, Translation(point))
    dists = f.(set.points)
    indices = sortperm(dists)
    set.points[indices]
end

function Base.show(io::IO, set::PointSet{Dim,T}) where {Dim,T}
    print(io, "PointSet{$Dim,$T}(", join(set.points, ", "), ')')
end

ext: jl


text:
using Wkhtmltox
using Base.Test

#### conversion to pdf ####

exfile = joinpath(dirname(@__FILE__),"../examples/example.html")
outfile = tempname() * ".pdf"

pdf_init(0)

ps = PdfSettings("out" => outfile,
                 "orientation" => "Landscape")

@test ps["out"] == outfile
@test ps["orientation"] == "Landscape"

ps["size.pageSize"] = "A4"
@test ps["size.pageSize"] == "A4"

os = PdfObject("page" => exfile)
@test os["page"] == exfile

os2 = PdfObject()
@test os2["page"] == ""
os2["page"] = exfile
@test os2["page"] == exfile

conv = PdfConverter(ps)
push!(conv, os)

@test run(conv) == 1

@test isfile(outfile)

conv = nothing
pdf_deinit()

#### conversion to image ####

exfile = joinpath(dirname(@__FILE__),"../examples/example.html")
outfile = tempname() * ".png"

img_init(1)

is = ImgSettings("in" => exfile,
                 "out" => outfile,
                 "fmt" => "png")

@test is["out"] == outfile
@test is["fmt"] == "png"

conv = ImgConverter(is)

@test run(conv) == 1

@test isfile(outfile)

conv = nothing
img_deinit()


text:
# Autogenerated wrapper script for gmsh_jll for x86_64-linux-musl-cxx03
export gmsh_api, libgmsh

JLLWrappers.@generate_wrapper_header("gmsh")
JLLWrappers.@declare_file_product(gmsh_api)
JLLWrappers.@declare_library_product(libgmsh, "libgmsh.so.4.8")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        gmsh_api,
        "lib/gmsh.jl",
    )

    JLLWrappers.@init_library_product(
        libgmsh,
        "lib/libgmsh.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
LCM Server: an LCM interface to Caesar.jl

=#

using JLD
using Caesar, RoME
using TransformUtils, Rotations, CoordinateTransformations
using Distributions
using PyCall, PyLCM
using LibBSON
using CloudGraphs # for sorryPedro

function gen_bindings()
    @show lcmtpath = joinpath(dirname(@__FILE__),"lcmtypes")
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_point_cloud_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_pose_node_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_pose_pose_nh_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_pose_pose_xyh_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_pose_pose_xyh_nh_t.lcm`)
    run(`lcm-gen -p --ppath $(lcmtpath) $(lcmtpath)/caesar_prior_zpr_t.lcm`)
    println("Adding lcmtypes dir to Python path: $(lcmtpath)")
    unshift!(PyVector(pyimport("sys")["path"]),lcmtpath)
end

println("[Caesar.jl] (re)generating LCM bindings")
gen_bindings()

println("[Caesar.jl] Importing LCM message types")
@pyimport caesar

function initialize!(backend_config,
                    user_config)
    println("[Caesar.jl] Setting up factor graph")
    fg = Caesar.initfg(sessionname=user_config["session"], cloudgraph=backend_config)
    println("[Caesar.jl] Creating SLAM client/object")
    return  SLAMWrapper(fg, nothing, 0)
end

function setRobotParameters!(cg::CloudGraph, session::AbstractString)
  hauvconfig = Dict()
  hauvconfig["robot"] = "hauv"
  hauvconfig["bTc"] = [0.0;0.0;0.0; 1.0; 0.0; 0.0; 0.0]
  hauvconfig["bTc_format"] = "xyzqwqxqyqz"
  # currently unused, but upcoming
  hauvconfig["pointcloud_description_name"] = "BSONpointcloud"
  hauvconfig["pointcloud_color_description_name"] = "BSONcolors"
  # robotdata = json(hauvconfig).data

  # Actually modify the databases
  insertrobotdatafirstpose!(cg, session, hauvconfig)
  nothing
end

"""
Adds pose nodes to graph with a prior on Z, pitch, and roll.
"""
function handle_poses!(slam::SLAMWrapper,
                       message_data)
    message = caesar.pose_node_t[:decode](message_data)
    id = message[:id]
    println("[Caesar.jl] Received pose message for x$(id)")

    mean = message[:mean]
    covar = message[:covar]
    t = [mean[1], mean[2], mean[3]]
    qw = mean[4]
    qxyz = [mean[5], mean[6], mean[7]]
    q = Quaternion(qw,qxyz) # why not a (w,x,y,z) constructor?
    pose = SE3(t,q)
    euler = Euler(q)

    node_label = Symbol("x$(id)")
    xn = addNode!(slam.fg, node_label, labels=["POSE"], dims=6) # this is an incremental inference call
    slam.lastposesym = node_label; # update object

    if id == 0
        println("[Caesar.jl] First pose")
        # this is the first message, and it does not carry odometry, but the prior on the first node.

        # add 6dof prior
        initPosePrior = PriorPose3( MvNormal( veeEuler(pose), diagm([covar...]) ) )
        addFactor!(slam.fg, [xn], initPosePrior)

        # auto init is coming, this code will likely be removed
        initializeNode!(slam.fg, node_label)
    end

    if id == 1
        # set robot parameters in the first pose, this will become a separate node in the future
        # visualization is still assuming parameters in node 1 (not 0), so leaving it here for now
        println("[Caesar.jl] Setting robot parameters")
        setRobotParameters!(slam.fg.cg, slam.fg.sessionname)
    end
end

# handles ZPR priors on poses
function handle_priors!(slam::SLAMWrapper,
                         message_data)

    message = caesar.prior_zpr_t[:decode](message_data)
    id = message[:id]
    println("[Caesar.jl] Adding prior on RPZ to x$(id)")

    node_label = Symbol("x$(id)")
    xn = getVert(slam.fg,node_label)

    z = message[:z]
    pitch = message[:pitch]
    roll = message[:roll]

    var_z = message[:var_z]
    var_pitch = message[:var_pitch]
    var_roll = message[:var_roll]

    rp_dist = MvNormal( [roll;pitch], diagm([var_roll, var_pitch]))
    z_dist = Normal(z, var_z)
    prior_rpz = PartialPriorRollPitchZ(rp_dist, z_dist)
    addFactor!(slam.fg, [xn], prior_rpz)
end

function handle_partials!(slam::SLAMWrapper,
                         message_data)
    message = caesar.pose_pose_xyh_t[:decode](message_data)

    origin_id = message[:node_1_id]
    destination_id = message[:node_2_id]
    origin_label = Symbol("x$(origin_id)")
    destination_label = Symbol("x$(destination_id)")

    println("[Caesar.jl] Adding XYH odometry constraint between(x$(origin_id), x$(destination_id))")

    delta_x = message[:delta_x]
    delta_y = message[:delta_y]
    delta_yaw = message[:delta_yaw]

    var_x = message[:var_x]
    var_y = message[:var_y]
    var_yaw = message[:var_yaw]

    origin_label, destination_label
    xo = getVert(slam.fg,origin_label)
    xd = getVert(slam.fg,destination_label)

    xyh_dist = MvNormal([delta_x, delta_y, delta_yaw], diagm([var_x, var_y, var_yaw]))
    xyh_factor = PartialPose3XYYaw(xyh_dist)
    addFactor!(slam.fg, [xo;xd], xyh_factor)

    initializeNode!(slam.fg, destination_label)
    println()
end

function handle_loops!(slaml::SLAMWrapper,
                       message_data)
    message = caesar.pose_pose_xyh_nh_t[:decode](message_data)

    origin_id = message[:node_1_id]
    destination_id = message[:node_2_id]
    origin_label = Symbol("x$(origin_id)")
    destination_label = Symbol("x$(destination_id)")

    delta_x = message[:delta_x]
    delta_y = message[:delta_y]
    delta_yaw = message[:delta_yaw]

    var_x = message[:var_x]
    var_y = message[:var_y]
    var_yaw = message[:var_yaw]
    confidence = message[:confidence]

    xo = getVert(slaml.fg,origin_label)
    xd = getVert(slaml.fg,destination_label)

    # if (destination_id - origin_id == 1)
    #     warn("Avoiding parallel factor! See: https://github.com/dehann/IncrementalInference.jl/issues/63To ")
    #     return
    # end

    println("[Caesar.jl] Adding XYH-NH loop closure constraint between (x$(origin_id), x$(destination_id))")
    xyh_dist = MvNormal([delta_x, delta_y, delta_yaw], diagm([var_x, var_y, var_yaw]))
    xyh_factor = PartialPose3XYYawNH(xyh_dist ,[1.0-confidence, confidence]) # change to NH
    addFactor!(slaml.fg, [xo;xd], xyh_factor )

    # println("[Caesar.jl] Adding P3P3NH loop closure constraint between (x$(origin_id), x$(destination_id))")

    # # line below fails!
    # lcf = Pose3Pose3NH( MvNormal(veeEuler(rel_pose), diagm(1.0./covar)), [0.5;0.5]) # define 50/50% hypothesis
    # lcf_label = Symbol[origin_label;destination_label]

    # addFactor!(slaml.fg, lcf_label, lcf)
end

# type Cloud
#     points::
#     colors::
# end

# function add_pointcloud!(slaml::SLAMWrapper, nodeID::Symbol, cloud::Cloud )
#     # fetch from database
#     vert = getVert(slaml.fg, nodeID, api=IncrementalInference.dlapi)

#     # add points blob
#     serialized_point_cloud = BSONObject(Dict("pointcloud" => cloud.points))
#     appendvertbigdata!(slaml.fg, vert, "BSONpointcloud", string(serialized_point_cloud).data)

#     # add colors blob
#     serialized_colors = BSONObject(Dict("colors" => cloud.colors))
#     appendvertbigdata!(slaml.fg, vert, "BSONcolors", string(serialized_colors).data)
# end

"""
   handle_clouds(slam::SLAMWrapper, message_data)

   Callback for caesar_point_cloud_t messages. Adds point cloud to SLAM_Client
"""
function handle_clouds!(slam::SLAMWrapper,
                        message_data)
    # TODO: interface here should be as simple as slam_client.add_pointcloud(nodeID, pc::SomeCloudType)

    # TODO: check for empty clouds!

    message = caesar.point_cloud_t[:decode](message_data)

    id = message[:id]

    last_pose = Symbol("x$(id)")
    println("[Caesar.jl] Got cloud $id")

    # 2d arrays of points and colors (from LCM data into arrays{arrays})
    points = [[pt[1], pt[2], pt[3]] for pt in message[:points]]
    colors = [[UInt8(c.data[1]),UInt8(c.data[2]),UInt8(c.data[3])] for c in message[:colors]]

    # TODO: check if vert exists or not (may happen if messages are lost or out of order)
    vert = getVert(slam.fg, last_pose, api=IncrementalInference.dlapi) # fetch from database

    # push to mongo (using BSON as a quick fix)
    # (for deserialization, see src/DirectorVisService.jl:cachepointclouds!)
    serialized_point_cloud = BSONObject(Dict("pointcloud" => points))
    appendvertbigdata!(slam.fg, vert, "BSONpointcloud", string(serialized_point_cloud).data)
    serialized_colors = BSONObject(Dict("colors" => colors))
    appendvertbigdata!(slam.fg, vert, "BSONcolors", string(serialized_colors).data)
end

# this function handles lcm messages
function listener!(slam::SLAMWrapper,
                   lcm_node::LCMCore.LCM)
    # handle traffic
    # TODO: handle termination
    while true
        handle(lcm_node)
    end
end

#println("[Caesar.jl] Prompting user for configuration")
#@load "usercfg.jld"
#user_config["session"] = "SESSHAUVDEV3"
@load "liljon17.jld"
#  include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
# user_config = addrdict
backend_config, user_config = standardcloudgraphsetup(addrdict=user_config)

println("[Caesar.jl] Setting up local solver")
slam_client = initialize!(backend_config,user_config)

# create new handlers to pass in additional data
lcm_pose_handler = (channel, message_data) -> handle_poses!(slam_client, message_data )
lcm_odom_handler = (channel, message_data) -> handle_partials!(slam_client, message_data )
lcm_prior_handler = (channel, message_data) -> handle_priors!(slam_client, message_data )
lcm_cloud_handler = (channel, message_data) -> handle_clouds!(slam_client, message_data )
lcm_loop_handler = (channel, message_data) -> handle_loops!(slam_client, message_data )

# create LCM object and subscribe to messages on the following channels
lcm_node = LCM()
# poses
subscribe(lcm_node, "CAESAR_POSES", lcm_pose_handler)
# factors
subscribe(lcm_node, "CAESAR_PARTIAL_XYH", lcm_odom_handler)
subscribe(lcm_node, "CAESAR_PARTIAL_ZPR", lcm_prior_handler)
subscribe(lcm_node, "CAESAR_PARTIAL_XYH_NH", lcm_loop_handler) # loop closures come in via p3p3nh factors
# sensor data
subscribe(lcm_node, "CAESAR_POINT_CLOUDS", lcm_cloud_handler)

println("[Caesar.jl] Running LCM listener")
listener!(slam_client, lcm_node)

#=
using JLD, Caesar
addrdict = nothing
#to avoid retyping credentials, uncomment the block below

try @load "credentials.jld"
catch addrdict=getcredentials()
@save "credentials.jld" addrdict
end

drawdbdirector(addrdict=addrdict);
=#


text:
using Pkg; Pkg.add("Plots")
using Plots
# IJulia.clear_output(true)

# Define functions and parameters
H=4000.;L=H/10.;A=20.;C=50.;D=2.5;B=8.;E=5.e-4;
ft(x,L)=exp(-x/L)
th1(z)=A*ft.(-z .+ L, L);
th2(z)=C*ft(D*(-z .+ L), L);

# Plot for equally spaced discrete points
npt=1000;
z=-collect(range(0, 1, length=npt)) .* H
phi1=th1.(z)
phi2=th2.(z)
theta=phi1 .- phi2 .+ B .+ E .* z;
plot( theta ,z,label=false, lw=3, color = :black )
plot!(phi1 .+ B ,z,label=false , line=(:dash) )
plot!(phi2 .+ B ,z,label=false , line=(:dash) )

# Plot for ECCO 50-level discrete points
dzECCO50l=[10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01,
    10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82,
    24.85,31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90,
     85.15, 90.18, 93.96, 96.58, 98.25, 99.25,100.01,101.33,
    104.56,111.33,122.83,139.09,158.94,180.83,203.55,226.50,
     249.50,272.50,295.50,318.50,341.50,364.50,387.50,410.50,
    433.50,456.50];
dz=dzECCO50l;
z=-(cumsum(dz).-dz./2);
phi1=th1.(z)
phi2=th2.(z)
theta=phi1 .- phi2 .+ B .+ E .* z;
plot( theta ,z,label=false, lw=3, color = :black )
plot!(phi1 .+ B ,z,label=false , line=(:dash) )
plot!(phi2 .+ B ,z,label=false , line=(:dash) )

max_stars_repo_name: christophernhill/climahacks

max_issues_repo_name: christophernhill/climahacks

max_forks_repo_name: christophernhill/climahacks


text:
"""
A 2D point of type `T`; has fields `x` and `y`, both of type `T`.
A frequent use case is `Point{Int}` for a point in the two-dimensional
Cartesian product of integers.
"""
struct Point{T}
  x::T
  y::T
end

"""
    zipf(N::Int, s::Float64)

Zipf's Law for `N` elements and scaling exponent `s`, i.e. the probability mass function

```math
f(k; s, N) = H_{N,s}^{-1} k^{-s}
```

for ``k = 1, \\dots , N``, where ``H_{N,s}`` is the generalized harmonic
number used for normalization.
"""
function zipf(N::Int, s::Float64)
  out = 1 ./ ((1:N) .^ s)
  out/sum(out)
end

"""
    manhattan(x::Point{Int}, y::Point{Int})

Find the Manhattan distance between two integer points.
"""
function manhattan(x::Point{Int}, y::Point{Int})
  abs(x.x - y.x) + abs(x.y - y.y)
end

"""
    chebyshev(x::Point{Int}, y::Point{Int})

Find the Chebyshev distance between two integer points.
"""
function chebyshev(x::Point{Int}, y::Point{Int})
  max(abs(x.x - y.x), abs(x.y - y.y))
end


text:
# Julia wrapper for header: /usr/include/libavutil/twofish.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_twofish_alloc,
    av_twofish_init,
    av_twofish_crypt

function av_twofish_alloc()
    ccall((:av_twofish_alloc, libavutil), Ptr{AVTWOFISH}, ())
end

function av_twofish_init(ctx, key, key_bits::Integer)
    ccall((:av_twofish_init, libavutil), Cint, (Ptr{AVTWOFISH}, Ptr{UInt8}, Cint), ctx, key, key_bits)
end

function av_twofish_crypt(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_twofish_crypt, libavutil), Cvoid, (Ptr{AVTWOFISH}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end


text:
# Run CountingBag program with with Bag1 or Bag2 as super-class.

function usage(msg = nothing)
    println("Usage: julia run_counting_bag.jl [Bag1|Bag2]")
    msg !== nothing && println(msg)
    exit(1)
end

# Check # of arguments
length(ARGS) != 1 && usage()

# Check argument 1 - must be Bag1 or Bag2
superclass = ARGS[1]
supported = ["Bag1", "Bag2"]
superclass in supported || usage("Please specify either Bag1 or Bag2")

# Compile
run(`cp $(superclass).java Bag.java`)
run(`javac Bag.java CountingBag.java`)
run(`java CountingBag`)


text:
function solve_model(M,rates,dens,cons)
    f = (du,u,p,t) -> M.model(du,u,p,t,dens,cons)
    prob = ODEProblem(f,M.u₀,M.tspan,rates)
    sol = solve(prob,AutoVern7(KenCarp4()), abstol=1e-8,reltol=1e-8)
end

function plot_diss_fit!(plt,M,rates,data;cmap=nothing)

    for i in 1:length(M.cons)
        # Setup the ODE problem, then solve
        sol = solve_model(M,rates,M.dens[i],M.cons[i])

        Σ = sum(Array(sol),dims=1)
        Σ = Im2Re(Σ)
        Σ = normalize_lnp_fit(Σ,data.y[i])
        # h = plot!(sol.t,real.(Σ)',ylims=(-5,0))

        if cmap==nothing
            color = i
        else
            color = cmap
        end
        plot!(plt,sol.t,Σ,
            ylims=(0.0,1.0),
            # xlims=(0,15),
            color=color,
            label=string(data.f[i]))
    end
    if cmap==nothing
        colors = collect(1:length(M.cons))'
    else
        colors = cmap
    end
    plot_dissociations!(plt,data,colors)
    # return plt
end

function Im2Re(S)  # ln of the solutions (gets rid of integration instabilities)
    Σ = Complex.(S)
    Σ = real.(log.(Σ)')
end

function normalize_lnp(lnp)
    normlnpv = lnp ./ log(length(lnp)) .+ 1.0
end

function normalize_lnp_fit(lnp_fit,lnp_data)
    normlnpv_fit = lnp_fit ./ log(length(lnp_data)) .+ 1.0
end

function plot_dissociations!(plt,data,colors)
    normlnp = normalize_lnp.(p2lnp.(data.y))
    ε=0.01
    scatter!(plt,data.t,normlnp,
            xlims = (0.0-ε, maximum(maximum.(data.t))+ε),
            ylims = (0.0-ε, 1.0+ε),
            markercolors = colors, #, markerstrokecolor = colors,
            # markersize = 5, markershape=:square,
            # markeralpha = 0., markerstrokealpha = 1.,
            label=false)
end

function show_diss_fit!(plt,M,rates,data;cmap=nothing)

    plot_diss_fit!(plt,M,rates,data;cmap)
    # return plt
end
# function Im2Re(S)
#     Σ = Complex.(S)
#     Σ = real.(Σ)
#     Σ = log.(Σ')
# end
#

# function plot_diss_fit(rates,dens,cons)
#     for (i,F) in enumerate(cons)
#         # Setup the ODE problem, then solve
#         solve_model(rates,dens,cons)
#
#         Σ = Complex.(sum(Array(sol),dims=1))
#         # h = plot!(sol.t,real.(Σ)',ylims=(-5,0))
#         h = plot!(sol.t,real.(log.(Σ)'),
#             ylims=(-6.,0),
#             xlims=(0,15),
#             color=i,
#             label=string(F),
#             palette = :Dark2_5)
#
#         if generate_data
#             push!(dataset,vec(generate_data(sol;logplot=true)))
#         end
#
#         # display(h)
#     end
#     return h
# end
#


text:
#**************************************************************************************
# Pythical_Model.jl
# =============== part of the GeoEfficiency.jl package.
#
# here is the place where all physical elements is being modeled and created as computer objects.
#
#**************************************************************************************

import Base.show

"""

	Point(Height::Real, Rho::Real)

construct and return a `Point` source that can be a source of itself or an `anchor point` of a source.

`Height` : point height relative to the detector.

`Rho` : point off axis relative to the detector axis of symmetry.
\n*****

	Point(Height::Real)
the same as `Point(Height::Real, Rho::Real)` but return an axial point.

`Height` : point height relative to the detector.
\n*****

	Point()

the same as `Point(Height::Real, Rho::Real)` but ask the user to provide the information from the `console`.
\n*****
`Note please`
\n- each detector interpreted the height in a different way as follow.
\n- for `CylDetector` the point source `height` is consider to be measured from the detector `face surface`.
\n- for `BoreDetector` the point source `height` is consider to be measured from the `detector middle`, +ve value are above the detector center while -ve are below.
\n- for `WellDetector` the point source `height` is considered to be measured from the detector `hole surface`.
"""
type Point
	Height::Float64
	Rho::Float64

	Point(Height::Real, Rho::Real) = new(float(Height), float(Rho))
	function Point(Height::Real)
		Point(Height, 0.0)
	end #function
	function Point()
		print_with_color(:yellow,"\n II- The Radiactive Source Anchor Point:-\n")
		Height = getfloat("\n\t > Height (cm) = ")
		Rho = getfloat("\n\t > Off-axis (cm) = ")
		Point(Height, Rho)
	end #function
end #type
function setHeight!(aPnt::Point, xH::Real)
	aPnt.Height = xH
	return aPnt
end #function
function setHeight(aPnt::Point, xH::Real)
	pnt::Point = deepcopy(aPnt)
	pnt.Height = xH
	return pnt
end #function
function setRho!(aPnt::Point, xRho::Real)
	aPnt.Rho = xRho
	return aPnt
end #function
function setRho(aPnt::Point, xRho::Real)
	pnt::Point = deepcopy(aPnt)
	pnt.Rho = xRho
	return pnt
end #function
id(aPnt::Point) = "Point[Height=$(aPnt.Height), Rho=$(aPnt.Rho)]"
show(pnt::Point) = print(id(pnt))

"""

	source(;isPoint=false)

return a tuple describing the source (`anchorPnt`, `SrcRadius`, `SrcLength`) based on the user input to the `console`.

`aPnt` : the source anchoring point.

`SrcRadius` : source radius.

`SrcLength` : source length.

If `isPoint` is true both `SrcRadius` and `SrcLength` are set to zero.
"""
function source(;isPoint=false)
    anchorPnt = Point()
	isPoint && return (anchorPnt, 0.0, 0.0)

	SrcRadius = getfloat("\n\t > Source Radius (cm) = ")
    if 0.0 != SrcRadius
        SrcLength = getfloat("\n\t > Source Length (cm) = ")
        anchorPnt.Rho = 0.0

	else
        SrcLength = 0.0

	end #if
    return (anchorPnt, SrcRadius, SrcLength)
end #function

"""

abstract base of all the Gamma Detectors.
"""
abstract RadiationDetector
show(io::IO, Detector::RadiationDetector) = print(id(Detector))

"""

	CylDetector(CryRadius::Real, CryLength::Real)

return a `cylindrical` detector.

`CryRadius` : the detector crystal radius.

`CryLength` : the detector crystal length.
\n*****

	CylDetector(CryRadius::Real)
return a cylindrical detector with crystal length 0.0.

`CryRadius` : the detector crystal radius.
\n*****

	CylDetector()

return a cylindrical detector according to the input from the `console`.
"""
immutable CylDetector <: RadiationDetector
	CryRadius::Float64    	#Real
    CryLength::Float64		#Real

	function CylDetector(CryRadius::Real, CryLength::Real)
		@assert CryRadius > 0.0	  	"Crystal Radius: expect +ve number, get $(CryRadius)."
		@assert CryLength >= 0.0  	"Crystal Length: expect +ve number or zero, get) $(CryLength)."
		new(float(CryRadius), float(CryLength))
	end #if
	function CylDetector(CryRadius::Real)
		CylDetector(CryRadius, 0.0)
	end #if
	function CylDetector()
		print_with_color(:yellow," I- The Cylinderical Detetcor physical Dimensions:-\n")
		CryRadius = getfloat("\n\t > Crystal Radius (cm) = ")
		CryLength = getfloat("\n\t > Crystal Length (cm) = ")
		CylDetector(CryRadius, CryLength)
	end #function
end #type
id(Detector::CylDetector) = "CylDetector[CryRadius=$(Detector.CryRadius), CryLength=$(Detector.CryLength)]"

"""

	BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real)

return a `bore-hole` detector.

`CryRadius` : the detector crystal radius.

`CryLength` : the detector crystal length.

`HoleRadius` : the detector hole radius.
\n*****

	BoreDetector()

return a bore-hole detector according to the input from the `console`.
"""
immutable BoreDetector <: RadiationDetector
	CryRadius::Float64    	#Real
    CryLength::Float64    	#Real
	HoleRadius::Float64    	#Real

	function BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real)
		@assert CryRadius > 0.0					"Crystal Radius: expect +ve number, get $(CryRadius)."
		@assert CryLength > 0.0					"Crystal Length: expect +ve number, get $(CryLength)."
		@assert CryRadius > HoleRadius > 0.0	"Hole Radius: expect +ve number Less than 'Crystal Radius=$(CryRadius)', get $(HoleRadius)."
		new(float(CryRadius), float(CryLength), float(HoleRadius))
	end #if
	function BoreDetector()
		print_with_color(:yellow," I- The Bore Hole Detetcor physical Dimensions:-\n")
		CryRadius = getfloat("\n\t > Crystal Radius (cm) = ")
		CryLength = getfloat("\n\t > Crystal Length (cm) = ")
		HoleRadius = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryRadius)
		BoreDetector(CryRadius, CryLength, HoleRadius)
	end #function
end #type
id(Detector::BoreDetector) = "BoreDetector[CryRadius=$(Detector.CryRadius), CryLength=$(Detector.CryLength), HoleRadius=$(Detector.HoleRadius)]"

"""

	WellDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)

return a Well-Type detector.

`CryRadius` : the detector crystal radius.

`CryLength` : the detector crystal length.

`HoleRadius` : the detector hole radius.

`HoleDepth` : the detector hole length.
\n*****

	WellDetector()

return a Well-Type detector according to the input from the `console`.
"""
immutable WellDetector <: RadiationDetector
	CryRadius::Float64
    CryLength::Float64
	HoleRadius::Float64
	HoleDepth::Float64

	function WellDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)
		@assert CryRadius > 0.0					"Crystal Radius: expect +ve number, get $(CryRadius)."
		@assert CryLength > 0.0					"Crystal Length: expect +ve number, get $(CryLength)."
		@assert CryRadius > HoleRadius > 0.0	"Hole Radius: expect +ve number Less than 'Crystal Radius=$(CryRadius)', get $(HoleRadius)."
		@assert CryLength > HoleDepth > 0.0	   	"Hole Depth: expect +ve number Less than 'Crystal Length=$(CryLength)', get $(HoleDepth)."
		new(float(CryRadius), float(CryLength), float(HoleRadius), float(HoleDepth))
	end #if
	function WellDetector()
		print_with_color(:yellow," I- The Well-Type Detetcor physical Dimensions:-\n")
		CryRadius = getfloat("\n\t > Crystal Radius (cm) = ")
		CryLength = getfloat("\n\t > Crystal Length (cm) = ")
		HoleRadius = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryRadius)
		HoleDepth = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryLength)
		WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth)
	end #function
end #type
id(Detector::WellDetector) = "WellDetector[CryRadius=$(Detector.CryRadius), CryLength=$(Detector.CryLength), HoleRadius=$(Detector.HoleRadius), HoleDepth=$(Detector.HoleDepth)]"

"""

	detectorFactory()

construct and return an object of the RadiationDetector type (`CylDetector`, `BoreDetector` or `WellDetector`)
according to the input from the console.
\n*****

it also has methods with argument list.

  `detectorFactory(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)`==> well-type or bore-hole or cylindrical detector
\n*****

  `detectorFactory(CryRadius::Real, CryLength::Real, HoleRadius::Real)` ==> bore-hole or cylindrical detector
\n*****

  `detectorFactory(CryRadius::Real, CryLength::Real)`	==> cylindrical detector
\n*****

  `detectorFactory(CryRadius::Real)` 	==> cylindrical detector
\n*****

`Note please`
\n- if any method with argument(s) take an `invalid` argument it would throw an error.
\n- if the value the last argument is `zero` of a method with `more` than one argument it behave as a missing argument.
"""
function detectorFactory()
	print_with_color(:yellow, "\n I- The Detector physical Dimensions :-\n")
	CryRadius = getfloat("\n\t > Crystal Radius (cm) = ")
	CryLength = getfloat("\n\t > Crystal Length (cm) = ")
	HoleRadius = getfloat("\n(zero for cylindrical detectors) > Hole Radius (cm) = ", 0.0, CryRadius)
	if   0.0 == HoleRadius
		return CylDetector(CryRadius, CryLength)

	else
		HoleDepth = getfloat("\n(zero for Bore-Hole detectors) > Hole Depth (cm) = ", 0.0, CryLength )
		return(0.0 == HoleDepth
				? BoreDetector(CryRadius, CryLength, HoleRadius)
				: WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth)	)
	end #if
end #function
detectorFactory(CryRadius::Real) = CylDetector(CryRadius)
detectorFactory(CryRadius::Real, CryLength::Real) = CylDetector(CryRadius, CryLength)
detectorFactory(CryRadius::Real, CryLength::Real, HoleRadius::Real) = (0.0 == HoleRadius
				?  CylDetector(CryRadius, CryLength)
				: BoreDetector(CryRadius, CryLength, HoleRadius))
detectorFactory(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real) = (0.0 == HoleDepth
				? detectorFactory(CryRadius, CryLength, HoleRadius)
				: WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth))
detectorFactory(Detector::RadiationDetector) = Detector


text:
@inline function simpleEstimate(model::Model, rpio::RPIO{Particle}) where Particle
  potentials = rpio.potentials
  v::Float64 = 0.0
  for p in 1:rpio.n
    _sampleParticles(model, rpio)
    v += log(_evaluatePotentials(model, rpio, p) / rpio.N)
  end
  return v
end


text:
export WaveletClassifier

using ScikitLearnBase
using ScikitLearn.CrossValidation: cross_val_score, StratifiedKFold
import ScikitLearn.@declare_hyperparameters

DIR = @__DIR__

include("$DIR/../tutorial/decoders.jl")
include("$DIR/../tutorial/neurons.jl")
include("$DIR/../tutorial/tutorial_lib.jl")

mutable struct WaveletClassifier <: BaseClassifier
    nsurr::Int
    nscales::Union{Symbol, Int}
    percentile::Int
    maxcoefs::Int
    mincoefs::Int
    model
    chosen_features
    mutual_information
    unbiased_mutual_information
    WaveletClassifier(nsurr, nscales, percentile, maxcoefs, mincoefs) = new(nsurr, nscales, percentile, maxcoefs, mincoefs)
end

function WaveletClassifier(; nsurr=100, nscales=:max, percentile=95, maxcoefs=5, mincoefs=2)
    WaveletClassifier(nsurr, nscales, percentile, maxcoefs, mincoefs)
end

@declare_hyperparameters(WaveletClassifier, Symbol[:nsurr, :nscales, :percentile, :maxcoefs, :mincoefs])

import ScikitLearnBase.is_classifier
is_classifier(clf::WaveletClassifier) = true

function ScikitLearnBase.fit!(lr::WaveletClassifier, X::AbstractArray{XT},
                              y::AbstractArray{yT}) where {XT, yT}
    decomposed_matrix = decompose(X)
    opts = WIopts(lr.nsurr, lr.nscales, lr.percentile, lr.maxcoefs, lr.mincoefs)
    feature_select,(mi,umi)  = WaveletInformation.mi_select(decomposed_matrix,
                                                        y, opts)
    lr.chosen_features = feature_select
    lr.mutual_information = mi
    lr.unbiased_mutual_information = umi

    X = decomposed_matrix[:,feature_select]
    y = y
    n = length(unique(y))
    lr.model = GaussianNB(priors=[1/n for _ in 1:n])
    ScikitLearn.fit!(lr.model, X, y)
    return lr
end

function ScikitLearnBase.predict(lr::WaveletClassifier, X)
    X2 = decompose(X)[:, lr.chosen_features]
    ScikitLearn.predict(lr.model, X2)
end


text:
const FILE_80 = """

 Couenne (/home/bgulcan/.julia/v0.6/AmplNLWriter/.solverdata/tmpMSxFvq.nl Mar  6 2018): Infeasible

Options
3
0
1
0
140
0
110
0
objno 0 220
"""

@testset "Test sol file parsing" begin
    @testset "Issue #80" begin
        io = IOBuffer()
        write(io, FILE_80)
        seekstart(io)
        model = AmplNLWriter.AmplNLMathProgModel("", String[], "")
        model.ncon = 140
        model.nvar = 110
        @test !AmplNLWriter.read_sol(io, model)
        @test model.solve_result_num == 220
        @test all(model.solution .=== NaN)
    end
end


text:
"""
Scale a segmented quantity by its mean and sd
"""
function superimpose(f::Function,q::Quantity,mask::Mask)
    _,qs = segment(q,mask)
    map(f,qs)
end

zscore(x) = (x-mean(x))/std(x)

superimpose(q::Quantity,mask::Mask) = superimpose(zscore,q,mask)

function superimpose(f::Function,q::Quantity,h::Stage,mask::Mask)
    qs = segment(q,mask)
    ts,hs = segment(h,mask)
    ms = map(x->findmax(x)[2],hs)
    maxtimes = [ts[i][ms[i]] for i in eachindex(ms)]
    ts.-maxtimes,superimpose(f,q,mask)
end


text:
module PlotsX

using UnPack, Plots
import FileIO: save

save(fn::String, p::AbstractPlot) = savefig(p, fn)

export save

### Utilities

function autoget_lims(x)
    xlims = [extrema(x[1,:])...]
    dx = xlims[2] - xlims[1]
    xlims += [-0.1dx, +0.1dx]
    ylims = [extrema(x[2,:])...]
    dy = ylims[2] - ylims[1]
    ylims += [-0.1dy, +0.1dy]
    if size(x, 1) == 3
        zlims = [extrema(x[3,:])...]
        dz = zlims[2] - zlims[1]
        zlims += [-0.1dz, +0.1dz]
    else
        zlims = nothing
    end
    return (xlims=xlims, ylims=ylims, zlims=zlims)
end

function autoset_lims!(p, x)
    dim = size(x, 1)
    xlims, ylims, zlims = autoget_lims(x)
    xlims!(p, xlims)
    ylims!(p, ylims)
    !isnothing(zlims) && zlims!(p, zlims)
end
autoset_lims!(x) = autoset_lims!(current(), x)

export autoget_lims, autoset_lims!

### Plots

include("line.jl")
export ErrorBarLines, TwoYAxesLines
include("function.jl")
export OneDimFunction, TwoDimFunction
include("scatter.jl")
export LabelledScatter
include("image.jl")
export ImageGrid, FeatureActivations, make_imggrid, count_leadingzeros
include("path.jl")
export TwoDimPath

end # module


text:
push!(LOAD_PATH,".")
MPRdir=pwd()
using PyPlot
using SIAMFANLEquations
using Printf
using MPResults

max_stars_repo_name: ctkelley/MPResults


text:
using NearestNeighbors
using StatsBase
using Statistics
export VertexNeighborhood,
       KNNNeighborhood,
       create,
       PixelNeighborhood,
       local_scale,
       neighbors,
       RandomNeighborhood,
       CliqueNeighborhood

import SpectralClustering: spatial_position
import Base.ones
"""
```julia
struct RandomKGraph
```
The type RandomKGraph defines the parameters needed to create a random k-graph.
Every vertex it is connected to `k` random neigbors.
# Members
- `number_of_vertices::Integer`. Defines the number of vertices of the graph.
- `k::Integer`. Defines the minimum number of  neighborhood of every vertex.
"""
struct RandomKGraph
    number_of_vertices::Integer
    k::Integer
end
"""
```julia
create(cfg::RandomKGraph)
```
Construct a [`RandomKGraph`](@ref) such that every vertex is connected with other k random vertices.
"""
function create(cfg::RandomKGraph)
    g = Graph(cfg.number_of_vertices)
    for i = 1:cfg.number_of_vertices
        cant = 0
        while cant < cfg.k
            selected = rand(1:cfg.number_of_vertices)
            while selected == i
                selected = rand(1:cfg.number_of_vertices)
            end
            connect!(g, i, selected, rand())
            cant = cant + 1
        end
    end
    return g;
end
"""
```julia
abstract type VertexNeighborhood end
```
The abstract type VertexNeighborhood provides an interface to query for the
neighborhood of a given vertex. Every concrete type that inherit from
VertexNeighborhood must define the function
```julia
neighbors{T<:VertexNeighborhood}(cfg::T, j::Integer, data)
```
which returns the neighbors list of the vertex j for the given data.
"""
abstract type VertexNeighborhood end
"""
```julia
struct PixelNeighborhood  <: VertexNeighborhood
```
`PixelNeighborhood` defines neighborhood for a given pixel based in its spatial location. Given a pixel located at (x,y), returns every pixel inside
\$(x+e,y), (x-e,y)\$ and \$(x,y+e)(x,y-e)\$.

# Members
- e:: Integer. Defines the radius of the neighborhood.

"""
struct PixelNeighborhood  <: VertexNeighborhood
    e::Integer
end
"""
```julia
neighbors(cfg::PixelNeighborhood, j::Integer, img)
```

Returns the neighbors of the pixel j according to the specified in [`PixelNeighborhood`](@ref)
"""
function neighbors(cfg::PixelNeighborhood, j::Integer, img::Matrix{T}) where T <: Colorant
    pos = CartesianIndices(img)[j]
    w_r = max(pos[1] - cfg.e, 1):min(pos[1] + cfg.e, size(img, 1))
    w_c = max(pos[2] - cfg.e, 1):min(pos[2] + cfg.e, size(img, 2))
    return vec(map(x->LinearIndices(img)[x[1],x[2]], CartesianIndices((w_r, w_c))))
end

"""
```julia
struct CliqueNeighborhood <: VertexNeighborhood
```
`CliqueNeighborhood` specifies that the neighborhood for a given vertex \$j\$ in a
graph of \$n\$ vertices are the remaining n-1 vertices.
"""
struct CliqueNeighborhood <: VertexNeighborhood
end
"""
```julia
neighbors(config::CliqueNeighborhood, j::Integer, X)
```
Return every other vertex index different from \$j\$. See [`CliqueNeighborhood`](@ref)
"""
function neighbors(config::CliqueNeighborhood, j::Integer, X)
  return filter!(x->x != j, collect(1:number_of_patterns(X)))
end
"""
```julia
struct KNNNeighborhood <: VertexNeighborhood
    k::Integer
    tree::KDTree
end
```
`KNNNeighborhood` specifies that the neighborhood for a given vertex \$j\$ are the \$k\$ nearest neighborgs. It uses a tree to search the nearest patterns.
# Members
- `k::Integer`. The number of k nearest neighborgs to connect.
- `tree::KDTree`. Internal data structure.
- `f::Function`. Transformation function
"""
struct KNNNeighborhood <: VertexNeighborhood
  k::Integer
  tree::KDTree
  t::Function
end
"""
```julia
KNNNeighborhood(X::Matrix, k::Integer)
```
Create the [`KNNNeighborhood`](@ref) type by building a `k`-nn tre from de data `X`

Return the indexes of the `config.k` nearest neigbors of the data point `j` of the data `X`.
"""
function KNNNeighborhood(X, k::Integer, f::Function = x->x)
   tree = KDTree(hcat([f(get_element(X, j)) for j = 1:number_of_patterns(X)]...))
    return KNNNeighborhood(k, tree, f)
end
neighbors(config::KNNNeighborhood, j::Integer, X) = neighbors(config, get_element(X, j))
function neighbors(config::KNNNeighborhood, data)
    idxs, dists = knn(config.tree, config.t(data), config.k + 1, true)
    return idxs[2:config.k + 1]
end
"""
```@julia
struct RandomNeighborhood <: VertexNeighborhood
    k::Integer
end
```
For a given index `j`return `k` random vertices different from `j`
"""
struct RandomNeighborhood <: VertexNeighborhood
    k::Integer
end
function neighbors(config::RandomNeighborhood, j::Integer, X)
   samples = StatsBase.sample(1:number_of_patterns(X), config.k, replace = false)
    if (in(j, samples))
        filter!(e->e != j, samples)
    end
    while (length(samples) < config.k)
       s  =  StatsBase.sample(1:number_of_patterns(X), 1)[1]
        if (s != j)
            push!(samples, s)
        end
    end
    return samples
end
"""
```julia
weight{T<:DataAccessor}(w::Function,d::T, i::Int,j::Int,X)
```
Invoke the weight function provided to compute the similarity between the pattern `i` and the pattern `j`.
"""
function weight(w::Function, i::Integer, j::Integer, X)
  x_i = get_element(X, i)
  x_j = get_element(X, j)
  return w(i, j, x_i, x_j)
end

"""
```julia
create(w_type::DataType, neighborhood::VertexNeighborhood, oracle::Function,X)
```
Given a [`VertexNeighborhood`](@ref), a simmilarity function `oracle`  construct a simmilarity graph of the patterns in `X`.
"""
function create(w_type::DataType, neighborhood::VertexNeighborhood, oracle::Function, X)
    number_of_vertices = number_of_patterns(X)
    g = Graph(number_of_vertices; weight_type = w_type)
    @Threads.threads for j = 1:number_of_vertices
        neigh = neighbors(neighborhood, j, X)
        x_j = get_element(X, j)
        x_neigh = get_element(X, neigh)
        weights = oracle(j, neigh, x_j, x_neigh)
        connect!(g, j, neigh, weights)
    end
    GC.gc()
    return g
end
"""
```julia
create(neighborhood::VertexNeighborhood, oracle::Function,X)
```
Given a [`VertexNeighborhood`](@ref), a simmilarity function `oracle` construct a simmilarity graph of the patterns in `X`.
"""
function create(neighborhood::VertexNeighborhood, oracle::Function, X)
    create(Float64, neighborhood, oracle, X)
end
"""
```julia
local_scale(neighborhood::KNNNeighborhood, oracle::Function, X; k::Integer = 7)
```
Computes thescale of each pattern according to [Self-Tuning Spectral Clustering](https://papers.nips.cc/paper/2619-self-tuning-spectral-clustering.pdf).
Return a matrix containing for every pattern the local_scale.

# Arguments
    - `neighborhood::KNNNeighborhood`
    - `oracle::Function`
    - `X`
      the data

\"The selection of thescale \$ \\sigma \$ can be done by studying thestatistics of the neighborhoods surrounding points \$ i \$ and \$ j \$ .i \"
Zelnik-Manor and Perona use \$ \\sigma_i = d(s_i, s_K) \$ where \$s_K\$ is the \$ K \$ neighbor of point \$ s_i \$ .
They \"used a single value of \$K=7\$, which gave good results even for high-dimensional data \" .

"""
function local_scale(neighborhood::T, oracle::Function, X; k::Integer = 7, sortdim::Integer=1) where T<:VertexNeighborhood
    sort_data(d::AbstractArray; dims=1) = sort(d)
    sort_data(d::AbstractMatrix; dims=1) = sort(d, dims=dims)
    temp = nothing
    distance_function = nothing
    try
        temp = oracle(get_element(X, 1), get_element(X, [1, 2]))
        distance_function = (a,b,c,d)->oracle(c, d)
    catch e
        temp = oracle(0, [0], get_element(X, 1), get_element(X, [1, 2]))
        distance_function = oracle
    end
    number_of_vertices = number_of_patterns(X)

    scales = zeros(size(temp, 2), number_of_vertices)
    for j = 1:number_of_vertices
        neigh = neighbors(neighborhood, j, X)
        distances = distance_function(j, neigh, get_element(X, j), get_element(X, neigh))

        scales[:, j] .= sort_data(distances, dims=sortdim)[k, :]
    end
    return scales
end

#="""
Given a graph (g) created from a X_prev \in R^{d x n}, updates de graph from
the matrix X \in R^{d x m}, m > n. Adding the correspondent vertices and connecting
them whith the existing ones.
"""
function update!(config::GraphCreationConfig,g::Graph,X)
  number_of_vertices = number_of_patterns(config.da,X)
  old_number_of_vertices = number_of_vertices(g)
  for j=old_number_of_vertices+1:number_of_vertices
    add_vertex!(g)
  end
  for j=old_number_of_vertices+1:number_of_vertices
      neigh = neighbors(config.neighborhood,j,X)
      for i in neigh
          w = weight(config.oracle,i,j,X)
          connect!(g,i,j,w)
      end
  end
end
=#

# Weight functions
constant(k) = (i::Integer, neigh, v, m) = ones(size(m, 2)) * k
ones(i::Integer, neigh, v, m) = ones(size(m, 2))


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using SHTOOLS

makedocs(; sitename="SHTOOLS", format=Documenter.HTML(), modules=[SHTOOLS])

deploydocs(; repo="github.com/eschnett/SHTOOLS.jl.git", devbranch="main",
           push_preview=true)


text:
# Container data structure (structure-of-arrays style) for DG elements
struct ElementContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 4}                   # [variables, i, j, elements]
  u_t::Array{Float64, 4}                 # [variables, i, j, elements]
  u_tmp2::Array{Float64, 4}              # [variables, i, j, elements]
  u_tmp3::Array{Float64, 4}              # [variables, i, j, elements]
  inverse_jacobian::Vector{Float64}      # [elements]
  node_coordinates::Array{Float64, 4}    # [orientation, i, j, elements]
  surface_ids::Matrix{Int}               # [direction, elements]
  surface_flux_values::Array{Float64, 4} # [variables, i, direction, elements]
  cell_ids::Vector{Int}                  # [elements]
end

function ElementContainer2D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, NVARS, n_nodes, n_nodes, capacity)
  u_t = fill(NaN, NVARS, n_nodes, n_nodes, capacity)
  # u_rungakutta is initialized to non-NaN since it is used directly
  u_tmp2 = fill(0.0, NVARS, n_nodes, n_nodes, capacity)
  u_tmp3 = fill(0.0, NVARS, n_nodes, n_nodes, capacity)
  inverse_jacobian = fill(NaN, capacity)
  node_coordinates = fill(NaN, 2, n_nodes, n_nodes, capacity)
  surface_ids = fill(typemin(Int), 2 * 2, capacity)
  surface_flux_values = fill(NaN, NVARS, n_nodes, 2 * 2, capacity)
  cell_ids = fill(typemin(Int), capacity)

  elements = ElementContainer2D{NVARS, POLYDEG}(u, u_t, u_tmp2, u_tmp3, inverse_jacobian, node_coordinates,
                                                surface_ids, surface_flux_values, cell_ids)

  return elements
end

# Return number of elements
nelements(elements::ElementContainer2D) = length(elements.cell_ids)

# Container data structure (structure-of-arrays style) for DG interfaces
struct InterfaceContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 4}      # [leftright, variables, i, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
end

function InterfaceContainer2D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), 2, capacity)
  orientations = fill(typemin(Int), capacity)

  interfaces = InterfaceContainer2D{NVARS, POLYDEG}(u, neighbor_ids, orientations)

  return interfaces
end

# Return number of interfaces
ninterfaces(interfaces::InterfaceContainer2D) = length(interfaces.orientations)

# Container data structure (structure-of-arrays style) for DG boundaries
struct BoundaryContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u::Array{Float64, 4}                # [leftright, variables, i, boundaries]
  neighbor_ids::Vector{Int}           # [boundaries]
  orientations::Vector{Int}           # [boundaries]
  neighbor_sides::Vector{Int}         # [boundaries]
  node_coordinates::Array{Float64, 3} # [orientation, i, elements]
end

function BoundaryContainer2D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u = fill(NaN, 2, NVARS, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)
  neighbor_sides = fill(typemin(Int), capacity)
  node_coordinates = fill(NaN, 2, n_nodes, capacity)

  boundaries = BoundaryContainer2D{NVARS, POLYDEG}(u, neighbor_ids, orientations, neighbor_sides,
                                                   node_coordinates)

  return boundaries
end

# Return number of boundaries
nboundaries(boundaries::BoundaryContainer2D) = length(boundaries.orientations)

# Container data structure (structure-of-arrays style) for DG L2 mortars
# Positions/directions for orientations = 1, large_sides = 2:
# mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
#           |    |
# upper = 2 |    |
#           |    |
#                | 3 = large side
#           |    |
# lower = 1 |    |
#           |    |
struct L2MortarContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u_upper::Array{Float64, 4} # [leftright, variables, i, mortars]
  u_lower::Array{Float64, 4} # [leftright, variables, i, mortars]
  neighbor_ids::Matrix{Int}  # [position, mortars]
  # Large sides: left -> 1, right -> 2
  large_sides::Vector{Int}   # [mortars]
  orientations::Vector{Int}  # [mortars]
end

function L2MortarContainer2D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u_upper = fill(NaN, 2, NVARS, n_nodes, capacity)
  u_lower = fill(NaN, 2, NVARS, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), 3, capacity)
  large_sides = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)

  l2mortars = L2MortarContainer2D{NVARS, POLYDEG}(u_upper, u_lower, neighbor_ids, large_sides, orientations)

  return l2mortars
end

# Return number of L2 mortars
nmortars(l2mortars::L2MortarContainer2D) = length(l2mortars.orientations)

# Allow printing container contents
function Base.show(io::IO, c::L2MortarContainer2D{NVARS, POLYDEG}) where {NVARS, POLYDEG}
  println(io, '*'^20)
  for idx in CartesianIndices(c.u_upper)
    println(io, "c.u_upper[$idx] = $(c.u_upper[idx])")
  end
  for idx in CartesianIndices(c.u_lower)
    println(io, "c.u_lower[$idx] = $(c.u_lower[idx])")
  end
  println(io, "transpose(c.neighbor_ids) = $(transpose(c.neighbor_ids))")
  println(io, "c.large_sides = $(c.large_sides)")
  println(io, "c.orientations = $(c.orientations)")
  println(io, '*'^20)
end

# Container data structure (structure-of-arrays style) for DG Ec mortars
# Positions/directions for large_sides = 1, orientations = 1:
#           |    |
# upper = 2 |    |
#           |    |
#                | 3
#           |    |
# lower = 1 |    |
#           |    |
struct EcMortarContainer2D{NVARS, POLYDEG} <: AbstractContainer
  u_upper::Array{Float64, 3} # [variables, i, mortars]
  u_lower::Array{Float64, 3} # [variables, i, mortars]
  u_large::Array{Float64, 3} # [variables, i, mortars]
  neighbor_ids::Matrix{Int}  # [position, mortars]
  # Large sides: left -> 1, right -> 2
  large_sides::Vector{Int}   # [mortars]
  orientations::Vector{Int}  # [mortars]
end

function EcMortarContainer2D{NVARS, POLYDEG}(capacity::Integer) where {NVARS, POLYDEG}
  # Initialize fields with defaults
  n_nodes = POLYDEG + 1
  u_upper = fill(NaN, NVARS, n_nodes, capacity)
  u_lower = fill(NaN, NVARS, n_nodes, capacity)
  u_large = fill(NaN, NVARS, n_nodes, capacity)
  neighbor_ids = fill(typemin(Int), 3, capacity)
  large_sides = fill(typemin(Int), capacity)
  orientations = fill(typemin(Int), capacity)

  ecmortars = EcMortarContainer2D{NVARS, POLYDEG}(u_upper, u_lower, u_large, neighbor_ids,
                                                  large_sides, orientations)

  return ecmortars
end

# Return number of EC mortars
nmortars(ecmortars::EcMortarContainer2D) = length(ecmortars.orientations)


text:
using LinearAlgebra

mutable struct SR1Hessian
    H::Matrix{Float64}
    x::Vector{Float64}
    g::Vector{Float64}
    utol::Float64
    function SR1Hessian(n, c; utol=1e-8)
        H0 = Matrix(c * I, n, n)
        x0 = zeros(n)
        g0 = zeros(n)
        return new(H0, x0, g0, utol)
    end
end

function update!(sr1::SR1Hessian, x, g)
    s = mul!(sr1.x, x, 1, 1, -1)
    y = mul!(sr1.g, g, 1, 1, -1)

    v = mul!(y, sr1.H, s, -1, 1)
    c = dot(v, s)
    if abs(c) > sr1.utol * norm(s) * norm(v)
        mul!(sr1, v, v', 1 / c, 1)
    end

    copyto!(sr1.x, x)
    copyto!(sr1.g, g)

    return sr1
end

struct CauchyPoint
    xc::Vector{Float64}
    p::Vector{Float64}
    t::Vector{Float64}
    idx::Vector{Int}
    Gp::Vector{Float64}
end

breakpoints(x, p, l, u) = p < 0 ? (x - u) / p : p > 0 ? (x - l) / p : Inf

function _cauchy_point!(c, gx, Hx, ti_m1)
    for i in c.idx
        δt = c.t[i] - ti_m1
        if δt != 0
            mul!(c.Gp, Hx, c.p)
            fc′ = dot(c.p, gx) + dot(c.xc, c.Gp)
            fc′′ = dot(c.p, c.Gp)
            if fc′ > 0
                return nothing
            end
            Δt = -fc′ / fc′′
            if Δt < δt
                mul!(c.xc, c.p, 1, Δt, 1)
                return nothing
            end
            mul!(c.xc, c.p, 1, δt, 1)
            ti_m1 = c.t[i]
        end
        c.p[i] = 0
    end
    return nothing
end

function cauchy_point!(c, gx, Hx, lbound, ubound)
    copyto!(c.p, gx)
    mul!(c.p, Hx, c.xc, 1, 1)
    c.t .= breakpoints.(c.xc, c.p, lbound, ubound)
    idx = sortperm!(c.idx, c.t)
    _cauchy_point!(c, gx, Hx, 0)
    return c
end

function gradient_projection(
    f, x0, lbound, ubound; rtol
)

end

function trust_region_gp(
    f, x0;
    lbound=fill(-Inf, size(x0)), ubound=fill(Inf, size(x0)),
    maxiter=100, rtol=1e-6, Δmax=Inf
)

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# This implements the first attempt at using bitarrays for
# storing and propagating the itemset information at each node

struct Node
    id::Int16
    item_ids::Array{Int16,1}
    transactions::BitArray{1}
    children::Array{Node,1}
    mother::Node
    supp::Int

    function Node(id::Int16, item_ids::Array{Int16,1}, transactions::BitArray{1})
        children = Array{Node,1}(undef, 0)
        nd = new(id, item_ids, transactions, children)
        return nd
    end

    function Node(id::Int16, item_ids::Array{Int16,1}, transactions::BitArray{1}, mother::Node, supp::Int)
        children = Array{Node,1}(undef, 0)
        nd = new(id, item_ids, transactions, children, mother, supp)
        return nd
    end
end

function has_children(nd::Node)
    res = length(nd.children) > 0
    res
end

function older_siblings(nd::Node)
    n_sibs = length(nd.mother.children)

    sib_indcs = map(x -> x.id > nd.id, nd.mother.children)
    return view(nd.mother.children, sib_indcs)
end

# This function is used internally and is the workhorse of the frequent()
# function, which generates a frequent itemset tree. The growtree!() function
# builds up the frequent itemset tree recursively.
function growtree!(nd::Node, minsupp, k, maxdepth)
    sibs = older_siblings(nd)

    for j = 1:length(sibs)
        transacts = nd.transactions .& sibs[j].transactions
        supp = sum(transacts)

        if supp ≥ minsupp
            items = zeros(Int16, k)
            items[1:k-1] = nd.item_ids[1:k-1]
            items[end] = sibs[j].item_ids[end]

            child = Node(Int16(j), items, transacts, nd, supp)
            push!(nd.children, child)
        end
    end
    # Recurse on newly created children
    maxdepth -= 1
    if maxdepth > 1
        for kid in nd.children
            growtree!(kid, minsupp, k+1, maxdepth)
        end
    end
end

function unique_items(transactions::Array{Array{M, 1}, 1}) where {M}
    dict = Dict{M, Bool}()

    for t in transactions
        for i in t
            dict[i] = true
        end
    end
    uniq_items = collect(keys(dict))
    return sort(uniq_items)
end

# This function is used internally by the frequent() function to create the
# initial bitarrays used to represent the first "children" in the itemset tree.
function occurrence(transactions::Array{Array{S, 1}, 1}, uniq_items::Array{S, 1}) where S
    n = length(transactions)
    p = length(uniq_items)

    itm_pos = Dict(zip(uniq_items, 1:p))
    res = falses(n, p)
    for i = 1:n
        for itm in transactions[i]
            j = itm_pos[itm]
            res[i, j] = true
        end
    end
    res
end

"""
    frequent_item_tree(transactions, minsupp, maxdepth)

This function creates a frequent itemset tree from an array of transactions.
The tree is built recursively using calls to the growtree!() function. The
`minsupp` and `maxdepth` parameters control the minimum support needed for an
itemset to be called "frequent", and the max depth of the tree, respectively
"""
function frequent_item_tree(transactions::Array{Array{S, 1}, 1}, uniq_items::Array{S, 1}, minsupp::Int, maxdepth::Int) where S
    occ = occurrence(transactions, uniq_items)

    # Have to initialize `itms` array like this because type inference
    # seems to be broken for this otherwise (using v0.6.0)
    itms = Array{Int16,1}(undef, 1)
    itms[1] = -1
    id = Int16(1)
    transacts = BitArray(undef, 0)
    root = Node(id, itms, transacts)
    n_items = length(uniq_items)

    # This loop creates 1-item nodes (i.e., first children)
    for j = 1:n_items
        supp = sum(occ[:, j])
        if supp ≥ minsupp
            nd = Node(Int16(j), Int16[j], occ[:, j], root, supp)
            push!(root.children, nd)
        end
    end
    n_kids = length(root.children)

    # Grow nodes in breadth-first manner
    for j = 1:n_kids
        growtree!(root.children[j], minsupp, 2, maxdepth)
    end
    root
end

function suppdict_to_dataframe(supp_lkup, item_lkup)
    n_sets = length(supp_lkup)
    df = DataFrame(itemset = fill("", n_sets), supp = zeros(Int, n_sets))
    i = 1

    for (k, v) in supp_lkup
        item_names = map(x -> item_lkup[x], k)
        itemset_string = "{" * join(item_names, " | ") * "}"
        df[i, :itemset] = itemset_string
        df[i, :supp] = v
        i += 1
    end
    df
end

# TODO: Fix the function below so that it handles int `minsupp`

"""
    frequent()
This function just acts as a bit of a convenience function that returns the frequent
item sets and their support count (integer) when given and array of transactions. It
basically just wraps frequent_item_tree() but gives back the plain text of the items,
rather than that Int16 representation.
"""
function frequent(transactions::Array{Array{S, 1}, 1}, minsupp::T, maxdepth) where {T <: Real, S}
    n = length(transactions)
    uniq_items = unique_items(transactions)
    item_lkup = Dict{Int16, S}()
    for (i, itm) in enumerate(uniq_items)
        item_lkup[i] = itm
    end
    if T <: Integer
        supp = minsupp
    elseif T == Float64
        supp = floor(Int, minsupp * n)
        if supp == 0
            warn("Minimum support should not be 0; setting it to 1 now.")
            supp = 1
        end
    end
    freq_tree = frequent_item_tree(transactions, uniq_items, supp, maxdepth)

    supp_lkup = gen_support_dict(freq_tree, n)

    freq = suppdict_to_dataframe(supp_lkup, item_lkup)
    return freq
end

"""
frequent_item_tree(occurrences, minsupp, maxdepth)

This function creates a frequent itemset tree from an occurrence matrix.
The tree is built recursively using calls to the growtree!() function. The
`minsupp` and `maxdepth` parameters control the minimum support needed for an
itemset to be called "frequent", and the max depth of the tree, respectively
"""
function frequent_item_tree(occ::BitArray{2}, minsupp::Int, maxdepth::Int)

    # Have to initialize `itms` array like this because type inference
    # seems to be broken for this otherwise (using v0.6.0)
    itms = Array{Int16,1}(undef, 1)
    itms[1] = -1
    id = Int16(1)
    transacts = BitArray(undef, 0)
    root = Node(id, itms, transacts)
    n_items = size(occ, 2)

    # This loop creates 1-item nodes (i.e., first children)
    for j = 1:n_items
        supp = sum(occ[:, j])
        if supp ≥ minsupp
            nd = Node(Int16(j), Int16[j], occ[:, j], root, supp)
            push!(root.children, nd)
        end
    end
    n_kids = length(root.children)

    # Grow nodes in depth-first manner
    for j = 1:n_kids
        growtree!(root.children[j], minsupp, 2, maxdepth)
    end
    root
end


text:
using Base.Meta, RecipesBase
using Distributions: Chisq, cdf

abstract type Significance end

struct EmpiricalSig{T<:Real} <: Significance
    gt::T
    count::Int
    values::Vector{T}
    p::Float64
    se::Float64

    function EmpiricalSig(gt::T, values::AbstractArray{T,1}) where {T<:Real}
        N = length(values)
        p = count(gt .≤ values) / N
        se = sqrt((p * (1-p))/N)
        new{T}(gt, N, values, p, se)
    end
end

function sig(::Type{EmpiricalSig}, func::Function, args::AbstractArray{<:Number,3}...;
             nperm=1000, parg=1, rng=Random.defautl_rng(), kwargs...)

    gt = func(args...; kwargs...)
    values = zeros(typeof(gt), nperm + 1)
    values[1] = gt
    head, arg, tail = args[1:parg-1], args[parg], args[parg+1:end]
    @views for i in 1:nperm
        perm = randperm(rng, size(arg, 2))
        values[i+1] = func(head..., arg[:,perm,:], tail...; kwargs...)
    end
    EmpiricalSig(gt, values)
end

function sig(::Type{EmpiricalSig}, func::Function, args::AbstractArray{<:Number,2}...;
             nperm=1000, parg=1, rng=Random.default_rng(), kwargs...)

    gt = func(args...; kwargs...)
    values = zeros(typeof(gt), nperm + 1)
    values[1] = gt
    head, arg, tail = args[1:parg-1], args[parg], args[parg+1:end]
    @views for i in 1:nperm
        perm = randperm(rng, size(arg, 2))
        values[i + 1] = func(head..., arg[:,perm], tail...; kwargs...)
    end
    EmpiricalSig(gt, values)
end

function sig(::Type{EmpiricalSig}, func::Function, args::AbstractArray{<:Number,1}...;
             nperm=1000, parg=1, rng=Random.default_rng(), kwargs...)

    gt = func(args...; kwargs...)
    values = zeros(typeof(gt), nperm + 1)
    values[1] = gt
    head, arg, tail = args[1:parg-1], args[parg], args[parg+1:end]
    @views for i in 1:nperm
        perm = randperm(rng, length(arg))
        values[i + 1] = func(head..., arg[perm], tail...; kwargs...)
    end
    EmpiricalSig(gt, values)
end

struct AnalyticSig{T<:Real} <: Significance
    gt::T
    N::Int
    dof::Float64
    dist::Chisq
    p::Float64
    function AnalyticSig(gt::T, N::Int, dof::Integer) where {T<:Real}
        dist = Chisq(dof)
        p = 1 - cdf(dist, 2 * N * gt)
        new{T}(gt, N, dof, dist, p)
    end
end

function sig(::Type{AnalyticSig}, func::Function, args::AbstractArray{<:Integer,3}...;
             nperm=1000, parg=1, rng=Random.default_rng(), kwargs...)
    gt = func(args...; kwargs...)
    dof = prod(map(a -> prod(maximum(a; dims=(2,3))) - 1, args))
    AnalyticSig(gt, size(args[1], 2), dof);
end

function sig(::Type{AnalyticSig}, func::Function, args::AbstractArray{<:Integer,2}...;
             nperm=1000, parg=1, rng=Random.default_rng(), kwargs...)
    gt = func(args...; kwargs...)
    dof = prod(map(a -> prod(maximum(a; dims=2)) - 1, args))
    AnalyticSig(gt, size(args[1], 2), dof);
end

function sig(::Type{AnalyticSig}, func::Function, args::AbstractArray{<:Integer,1}...;
             nperm=1000, parg=1, rng=Random.default_rng(), kwargs...)
    gt = func(args...; kwargs...)
    dof = prod(map(a -> prod(maximum(a)) - 1, args))
    AnalyticSig(gt, length(args[1]), dof);
end

macro sig(method, func, args...)
    if !isexpr(func, :call)
        error("first expression must be a function call")
    end

    kwargs = Dict{Any,Any}(:nperm => 1000, :parg => 1, :rng => :(Random.default_rng()))
    foreach(args) do arg
        if isexpr(arg, :(=))
            kwargs[arg.args[1]] = arg.args[2]
        else
            error("each argument after the first must be a keyword argument")
        end
    end

    kw = if length(func.args) ≥ 2 && isexpr(func.args[2], :parameters)
        append!(func.args[2].args, [Expr(:kw, arg[1], arg[2]) for arg in kwargs])
        func.args[2]
    else
        Expr(:parameters, [Expr(:kw, arg[1], arg[2]) for arg in kwargs]...)
    end

    if length(func.args) == 1
        func.args = [:sig, kw, method, func.args[1]]
    elseif length(func.args) == 2 && isexpr(func.args[2], :parameters)
        func.args = [:sig, kw, method, func.args[1]]
    elseif isexpr(func.args[2], :parameters)
        func.args = [:sig, kw, method, func.args[1], func.args[3:end]...]
    else
        func.args = [:sig, kw, method, func.args...]
    end

    esc(func)
end

@recipe function f(s::EmpiricalSig)
    xguide --> "value"
    yguide --> "PDF"
    @series begin
        seriestype := :density
        label --> nothing
        s.values
    end
    @series begin
        seriestype := :vline
        label --> "ground truth"
        [s.gt]
    end
end

@recipe function f(s::AnalyticSig)
    xguide --> "value"
    yguide --> "PDF"
    @series begin
        label --> nothing
        s.dist
    end
    @series begin
        seriestype := :vline
        label --> "ground truth"
        [s.gt]
    end
end


text:
module Support
  # this module contains supporting functions used by MatlabCompat library

  # Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
  #
  #  MIT license.
  #
  # Permission is hereby granted, free of charge, to any person
  # obtaining a copy of this software and associated documentation files (the
  # "Software"), to deal in the Software without restriction, including without
  # limitation the rights to use, copy, modify, merge, publish, distribute,
  # sublicense, and/or sell copies of the Software, and to permit persons to whom
  # the Software is furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included in all
  # copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  # SOFTWARE.

  using Images
  using ImageView

  export mat2im
  export mat2im
  export rosetta

  function mat2im(array)
    image = grayim(array)
    return image
  end

  function im2mat(image)
    array = reinterpret(Float32, float32(image))
    return convert(Array, array)
  end

  # an m-file parser aimed at converting them as close as possible from MATLAB syntax to Julia syntax
  #using
 function rosetta(filePath...)
    # parse arguments
    if length(filePath) >= 2
      inputMfilePath = filePath[1]
      outputJlFilePath = filePath[2]
    elseif length(filePath) == 1
      inputMfilePath = filePath[1]
    else
      error("not enough input arguments for rosetta()")
    end
    # read the m-file
    mFileContents = open(readlines, inputMfilePath)

    # here we parse the MATLAB/Octave code to be compatible with Julia through MatlabCompat library
    mFileContentsParsed = mFileContents;

    for iLine = 1:size(mFileContents,1)
      # 1. substitute % by # excluding
      if ismatch(r"^(\t*|\s*)%.*", mFileContents[iLine])
        # 1a. match the simplest case when % is in the beginning of the line with any number of white spaces or tabs. If true replace the first occurance of %
        mFileContentsParsed[iLine] = replace(mFileContents[iLine], "%", "#", 1)
      elseif ismatch(r".*\'.*", mFileContents[iLine]) && ismatch(r".*%.*", mFileContents[iLine])
        # 1b. match a complex case where % may be inside of the single quotes - this % shouldn't be replaced
        println("\' and % present");
        fragmentedString = split(mFileContents[iLine], "%")
        numberOfQuotes = 0
        fragmentToComment = 0
        newLine = ""
        firstOccurance = true
        for iFragment = 1:length(fragmentedString)
          # count the quotes
          if ismatch(r"\'", fragmentedString[iFragment])
            numberOfQuotes = numberOfQuotes + length(matchall(r"\'", fragmentedString[iFragment]))
          end
          # if number of quotes is even - they are closed, we can exchange the first occurance of % with # safely
          if (iseven(numberOfQuotes) && numberOfQuotes != 0 && firstOccurance == true)
            newLine = string(newLine, fragmentedString[iFragment], "#");
            firstOccurance = false;
          # is it last fragment?
          elseif (iFragment == length(fragmentedString))
            newLine = string(newLine, fragmentedString[iFragment]);
          else
              newLine = string(newLine, fragmentedString[iFragment], "%");
          end
        end
        mFileContentsParsed[iLine] = newLine
        elseif ~ismatch(r".*\'.*", mFileContents[iLine]) && ismatch(r".*%.*", mFileContents[iLine])
          # 1c. match a case where only % symbols are present
          println("% present");
          mFileContentsParsed[iLine] = replace(mFileContents[iLine], "%", "#")
        else
          println("no comment detected");
          mFileContentsParsed[iLine] = mFileContents[iLine];
        end
        # 2. substitute all single quotes ' by double quotes "
        if ismatch(r".*\'.*", mFileContentsParsed[iLine])
          mFileContentsParsed[iLine] = replace(mFileContentsParsed[iLine], "\'", "\"");
        end
      end
    # 3. append the code array with "using ..."

    extraLines = ["#This Julia file has been generated by rosetta script of MatlabCompat library from an m-file\n\r";
                "#The code generated need further corrections by you. Execute it line by line and correct the errors.\n\r";
                "using MatlabCompat\n\r importall MatlabCompat.ImageTools\n\r importall MatlabCompat.MathTools\n\r"];
    mFileContentsParsed = vcat(extraLines,mFileContentsParsed)
    if length(filePath) >= 2
      # write the jl-file
      jlFileStream = open(outputJlFilePath, "w")
      write(jlFileStream, mFileContentsParsed);
      close(jlFileStream);
      elseif length(filePath) == 1
      return mFileContentsParsed;
    else
      return false
    end
  end
end #End of Support


text:
using PyPlot
using SDELab2
# preferred pdf size
width_inches=4; height_inches=3
#######################################
function get_soln(fcn,quad0,params,phi,extrap_flag)
  println("Computing with dt=",params["MaxStepSize"])
  println(method_code[params["Solver"]])
  tic()
  tspan=[0,T]
  if extrap_flag==true
    quad,stats=SDELab2.sde_gq_solution_extrap(fcn,tspan,
                                      quad0,params)
  else
    quad,stats=SDELab2.sde_gq_solution(fcn,tspan,
                                      quad0,params)

  end
  cpu_time=toq()
  println("Mean number of points ",mean(stats["no_pts"]))
  println(stats["no_pts"])
  println("cpu time ", cpu_time)
  qoi=gauss_quadrature.do_quad(phi,quad)
  return qoi,cpu_time,stats,quad
end
#########################################
function helper(qoi,exact,cpu_time,str)
  println(str)
  println("qoi",qoi)
  error=abs(exact-qoi)
  println("error ", error)
  return error, cpu_time
end
##########################################
function initial_final_cdf_pdf(q0,q1, N, s)
  str2="k-"
  str2b="b:"
  str1="b-"
  clf()
  x,pdf,cdf,x1,cdf1=gauss_quadrature.quad_to_pdf(q1,N[1], s[1])
  q00=gauss_quadrature.QUAD((q0.length+1), [0;q0.X],[0;q0.w])
  x0,pdf0,cdf0,x10,cdf10=gauss_quadrature.quad_to_pdf(q00,N[2],s[2])
  plot(x,pdf,str1)
  plot(x0,pdf0,str2)
  xlabel(utf8("X"))
  ylabel(utf8("pdf"))
  axis([minimum(x),maximum(x),0,maximum([pdf;pdf0]) ])
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(utf8("pdf.pdf"))
  #
  clf()
  plot(x,cdf,str1)
  plot(x1,cdf1,str2b)
  plot(x10,cdf10,str2)
  axis([minimum(x),maximum(x),0,maximum([cdf;cdf0]) ])
  xlabel(utf8("X"))
  ylabel(utf8("cdf"))
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(utf8("cdf.pdf"))
end
#

function initial_final_cdf_pdf(q0,q1)
  str2="k-"
  str2b="b:"
  str1="b-"
  clf()
  x,pdf,cdf,x1,cdf1=gauss_quadrature.quad_to_pdf(q1)#100,0.08)
  q00=gauss_quadrature.QUAD((q0.length+1), [0;q0.X],[0;q0.w])
  x0,pdf0,cdf0,x10,cdf10=gauss_quadrature.quad_to_pdf(q00)#40,0.01)
  plot(x,pdf,str1)
  plot(x0,pdf0,str2)
  xlabel(utf8("X"))
  ylabel(utf8("pdf"))
  axis([minimum(x),maximum(x),0,maximum([pdf;pdf0]) ])
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(utf8("pdf.pdf"))
  #
  clf()
  plot(x,cdf,str1)
  plot(x1,cdf1,str2b)
  plot(x10,cdf10,str2)
  axis([minimum(x),maximum(x),0,maximum([cdf;cdf0]) ])
  xlabel(utf8("X"))
  ylabel(utf8("cdf"))
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(utf8("cdf.pdf"))
end
#
function loglog_safe(a,b,s)
  if length(a)==length(b)
    loglog(a,b,s)
  end
end
#
function time_step_err(dt1,e1,dt2,e2,dt3,e3,s,filename)
  clf()
  str1="b--" # extrap
  str2="g-"
  str3="k-." # 2nd-order
  axis("auto")
  plt[:rcParams]["xtick.labelsize"]="small"
  plt[:rcParams]["ytick.labelsize"]="small"

  loglog_safe(dt1,e1,str1)#dtpts,errors[:,i],str1)
  loglog_safe(dt2,e2,str2)#dtpts,errors_no[:,i])
  loglog_safe(dt3,e3,str3)#dtpts,errors_nd[:,i],str2)
  loglog_safe(dt1,s[1]*e1[1]*(dt1/dt1[1]).^(2),"y-")
  loglog_safe(dt2,s[2]*e2[1]*(dt2/dt2[1]).^(1),"y-")
  xlabel(utf8("time step"))
  ylabel(utf8("error"))
  axis([minimum([dt1;dt2;dt3]), maximum([dt1;dt2;dt3]), minimum([e1;e2;e3]), maximum([e1;e2;e3]) ] )
  axis([5e-4,1e-1,1e-8,1e-2])
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(filename)
end
#
function err_cpu(e1,c1,e2,c2,e3,c3,s,filename)
  clf()
  axis("auto")
  str1="b--" # extrap
  str2="g-"
  str3="k-." # 2nd-order
  plt[:rcParams]["xtick.labelsize"]="small"
  plt[:rcParams]["ytick.labelsize"]="small"
  loglog_safe(e1,c1,str1)#dtpts,errors[:,i],str1)
  loglog_safe(e2,c2,str2)#dtpts,errors_no[:,i])
  loglog_safe(e3,c3,str3)#dtpts,errors_nd[:,i],str2)
  loglog_safe(e1,s[1]*c1[1]*(e1/e1[1]).^(-1/2),"y-")
  loglog_safe(e2,s[2]*c2[1]*(e2/e2[1]).^(-1),"y-")
  xlabel(utf8("error"))
  ylabel(utf8("cpu time (secs)"))

  axis([minimum([e1;e2;e3]), maximum([e1;e2;e3]), minimum([c1;c2;c3]), maximum([c1;c2;c3]) ] )
  gcf()[:set_size_inches](width_inches,height_inches)
  tight_layout()
  savefig(filename)
end


text:
"""
Module for generating and manipulating voxel 3-D images (boxes of voxels).

This module may be helpful when working with CT scan and other medical images
in the context of finite element analysis carried out with the package
FinEtools.jl .
"""
module VoxelBoxModule

import Base.size

mutable struct VoxelBoxVolume{CoordT<:Number,DataT<:Number}
    origin::Array{CoordT,1}
    boxdim::Array{CoordT,1}
    data::Array{DataT,3}
end

function VoxelBoxVolume(::Type{CoordT}, ::Type{DataT}) where {CoordT<:Number,DataT<:Number}
    V = VoxelBoxVolume(zeros(CoordT, 3), zeros(CoordT, 3), fill(zero(DataT), 0, 0, 0));
    return V
end

function VoxelBoxVolume(::Type{DataT}, nvox::Array{Int,1}, boxdim::Array{CoordT,1}) where {CoordT<:Number,DataT<:Number}
    origin = zeros(CoordT,3)
    data = zeros(DataT,nvox...)
    V = VoxelBoxVolume(origin, boxdim, data);
    return V
end

function VoxelBoxVolume(data::Array{DataT,3}, boxdim::Array{CoordT,1})  where {CoordT<:Number,DataT<:Number}
    V = VoxelBoxVolume(CoordT,DataT);
    V.boxdim = deepcopy(boxdim)
    V.data = deepcopy(data)
    return V
end

function voxeldims(V::VoxelBoxVolume)
    nx, ny, nz = size(V.data)
    voxszx=V.boxdim[1]/nx
    voxszy=V.boxdim[2]/ny
    voxszz=V.boxdim[3]/nz
    return (voxszx, voxszy, voxszz)
end

size(V::VoxelBoxVolume, which) = size(V.data, which)
size(V::VoxelBoxVolume) = size(V.data)

struct SolidCF{F<:Function}
    f::F
end

function (o::SolidCF)(x::CoordT, y::CoordT, z::CoordT) where {CoordT<:Number}
    return o.f(x, y, z)
end

function intersectionop(o1::SolidCF, o2::SolidCF)
    function f(x, y, z)
        return o1(x, y, z) && o2(x, y, z)
    end
    return SolidCF(f)
end

function unionop(o1::SolidCF, o2::SolidCF)
    function f(x, y, z)
        return o1(x, y, z) || o2(x, y, z)
    end
    return SolidCF(f)
end

function complementop(o1::SolidCF)
    function f(x, y, z)
        return !o1(x, y, z)
    end
    return SolidCF(f)
end

function differenceop(o1::SolidCF, o2::SolidCF)
    function f(x, y, z)
        return o1(x, y, z) && (!o2(x, y, z))
    end
    return SolidCF(f)
end

"""
    solidsphere(center::Tuple{CoordT, CoordT, CoordT}, r::CoordT) where {CoordT<:Number}

Solid sphere.
"""
function solidsphere(center::Tuple{CoordT, CoordT, CoordT}, r::CoordT) where {CoordT<:Number}
    return SolidCF((x, y, z) -> (r*r - ((x-center[1])^2+(y-center[2])^2+(z-center[3])^2) >= 0.0))
end

"""
    solidhalfspace(center::Tuple{CoordT, CoordT, CoordT},
        normal::Tuple{CoordT, CoordT, CoordT}) where {CoordT<:Number}

Solid halfspace.
"""
function solidhalfspace(center::Tuple{CoordT, CoordT, CoordT},
    normal::Tuple{CoordT, CoordT, CoordT}) where {CoordT<:Number}
    nn = sqrt(normal[1]^2 + normal[2]^2 + normal[3]^2)
    n = (normal[1]/nn, normal[2]/nn, normal[3]/nn)
    function g(x, y, z)
        return (x-center[1])*normal[1] + (y-center[2])*normal[2] + (z-center[3])*normal[3]
    end
    return SolidCF((x, y, z) -> (-g(x, y, z) >= 0.0))
end

"""
    solidbox(corner1::Tuple{CoordT, CoordT, CoordT},
        corner2::Tuple{CoordT, CoordT, CoordT}) where {CoordT<:Number}

Solid box  with faces aligned with the global Cartesian axes.
"""
function solidbox(corner1::Tuple{CoordT, CoordT, CoordT},
    corner2::Tuple{CoordT, CoordT, CoordT}) where {CoordT<:Number}
    cmin = (min(corner1[1], corner2[1]), min(corner1[2], corner2[2]), min(corner1[3], corner2[3]))
    cmax = (max(corner1[1], corner2[1]), max(corner1[2], corner2[2]), max(corner1[3], corner2[3]))
    i1 = intersectionop(solidhalfspace(cmin, (-1.0, 0.0, 0.0)), solidhalfspace(cmax, (+1.0, 0.0, 0.0)))
    i2 = intersectionop(solidhalfspace(cmin, (0.0, -1.0, 0.0)), solidhalfspace(cmax, (0.0, +1.0, 0.0)))
    i3 = intersectionop(solidhalfspace(cmin, (0.0, 0.0, -1.0)), solidhalfspace(cmax, (0.0, 0.0, +1.0)))
    return intersectionop(i1, intersectionop(i2, i3))
end

function solidcylinder(center::Tuple{CoordT, CoordT, CoordT},
    axisdirection::Tuple{CoordT, CoordT, CoordT}, radius::CoordT) where {CoordT<:Number}
    nad = sqrt(axisdirection[1]^2 + axisdirection[2]^2 + axisdirection[3]^2)
    ad = (axisdirection[1]/nad, axisdirection[2]/nad, axisdirection[3]/nad)
    function g(x, y, z)
        xcad = (x-center[1])*ad[1] + (y-center[2])*ad[2] + (z-center[3])*ad[3]
        px = x-center[1] - xcad*ad[1]
        py = y-center[2] - xcad*ad[2]
        pz = z-center[3] - xcad*ad[3]
        return ((px^2+py^2+pz^2) - radius*radius)
    end
    return SolidCF((x, y, z) -> (-g(x, y, z) >= 0.0))
end

"""
    fillsolid!(V::VoxelBoxVolume,
        f::SolidCF, fillvalue::DataT) where {DataT<:Number}

Filled a solid using a solid characteristic function.
"""
function fillsolid!(V::VoxelBoxVolume,
    f::SolidCF, fillvalue::DataT) where {DataT<:Number}
    nx, ny, nz = size(V.data)
    voxszx=V.boxdim[1]/nx
    voxszy=V.boxdim[2]/ny
    voxszz=V.boxdim[3]/nz
    for i = 1:nx
        x = V.origin[1] + voxszx*(i+0.5)
        for j = 1:ny
            y = V.origin[2] + voxszy*(j+0.5)
            for k = 1:nz
                z = V.origin[3] + voxszz*(k+0.5)
                if f(x,y,z)
                    V.data[i,j,k] = fillvalue
                end
            end
        end
    end
    return V
end

"""
    fillvolume!(V::VoxelBoxVolume, fillvalue::DataT) where {DataT<:Number}

Fill the volume with a given value.
"""
function fillvolume!(V::VoxelBoxVolume, fillvalue::DataT) where {DataT<:Number}
    fill!(V.data, fillvalue)
    return V
end

"""
    trim(V::VoxelBoxVolume, emptyvalue)

Trim off pieces of the volume that consist only of the empty value.
"""
function trim(V::VoxelBoxVolume, emptyvalue)
    emptyvalue = convert(eltype(V.data[1]), emptyvalue)
    function sliceisempty(slice, emptyvalue)
        for jx = 1:size(slice, 2)
            for ix = 1:size(slice, 1)
                if slice[ix, jx] != emptyvalue
                    return false
                end
            end
        end
        return true
    end
    xmin = 1
    for i = 1:size(V.data, 1)
        if !sliceisempty(view(V.data, i, :, :), emptyvalue)
            xmin = i; break
        end
    end
    xmax = size(V.data, 1)
    for i = size(V.data, 1):-1:1
        if !sliceisempty(view(V.data, i, :, :), emptyvalue)
            xmax = i; break
        end
    end
    ymin = 1
    for i = 1:size(V.data, 2)
        if !sliceisempty(view(V.data, :, i, :), emptyvalue)
            ymin = i; break
        end
    end
    ymax = size(V.data, 2)
    for i = size(V.data, 2):-1:1
        if !sliceisempty(view(V.data, :, i, :), emptyvalue)
            ymax = i; break
        end
    end
    zmin = 1
    for i = 1:size(V.data, 3)
        if !sliceisempty(view(V.data, :, :, i), emptyvalue)
            zmin = i; break
        end
    end
    zmax = size(V.data, 3)
    for i = size(V.data, 3):-1:1
        if !sliceisempty(view(V.data, :, :, i), emptyvalue)
            zmax = i; break
        end
    end
    if (xmin == 1) && (xmax == size(V.data, 1)) &&
        (ymin == 1) && (ymax == size(V.data, 2))  &&
        (zmin == 1) && (zmax == size(V.data, 3))
        return V
    end
    origin = V.origin .+ vec([xmin, ymin, zmin]) .* voxeldims(V)
    boxdim = vec([(xmax - xmin) (ymax - ymin) (zmax - zmin)]) .* voxeldims(V)
    return VoxelBoxVolume(origin, boxdim, V.data[xmin:xmax, ymin:ymax, zmin:zmax])
end

"""
    pad(V::VoxelBoxVolume, ipad, jpad, kpad, padvalue)

Pad voxel box with a constant value.
"""
function pad(V::VoxelBoxVolume, ipad, jpad, kpad, padvalue)
    # Adjust the number of pixels
    originalvoxeldims = voxeldims(V)
    data = zeros(eltype(V.data[1]), size(V) .+ (sum(ipad), sum(jpad), sum(kpad)))
    fill!(data, padvalue)
    data[ipad[1]+1:ipad[1]+size(V, 1), jpad[1]+1:jpad[1]+size(V, 2), kpad[1]+1:kpad[1]+size(V, 3)] = V.data
    origin = V.origin .- vec([ipad[1], jpad[1], kpad[1]]) .* voxeldims(V)
    boxdim = [x for x in size(data) .* originalvoxeldims]
    return VoxelBoxVolume(origin, boxdim, data)
end

"""
    threshold(V, threshold_value, voxel_below, voxel_above)

Threshold the data.
"""
function threshold(V, threshold_value, voxel_below, voxel_above)
    for k= 1:size(V, 3)
        for j= 1:size(V, 2)
            for i= 1:size(V, 1)
                if V.data[i,j,k] > threshold_value
                    V.data[i,j,k] = voxel_above
                else
                    V.data[i,j,k] = voxel_below
                end
            end
        end
    end
    return V
end

"""
    vtkexport(theFile::String, V::VoxelBoxVolume{CoordT,DataT}) where {CoordT<:Number,DataT<:Number}

Compute.
"""
function vtkexport(theFile::String, V::VoxelBoxVolume{CoordT,DataT}) where {CoordT<:Number,DataT<:Number}
    open(theFile,"w") do fid
        if (fid==-1)
            error(["Could not open " * theFile])
            return nothing
        end
        print(fid,"# vtk DataFile Version 2.0\n");
        print(fid,"Example\n");
        print(fid,"ASCII\n");
        print(fid,"DATASET STRUCTURED_POINTS\n");
        nx, ny, nz = size(V.data)
        print(fid,"DIMENSIONS $(nx) $(ny) $(nz)\n");
        voxszx=V.boxdim[1]/nx
        voxszy=V.boxdim[2]/ny
        voxszz=V.boxdim[3]/nz
        print(fid,"SPACING $(voxszx) $(voxszy) $(voxszz)\n");
        print(fid,"ORIGIN 0 0 0\n");
        print(fid,"POINT_DATA $(nx*ny*nz)\n");
        #     bit , unsigned_char , char , unsigned_short , short , unsigned_int , int ,
        # unsigned_long , long , float , or double
        Types =  Dict{DataType, String}(Int8=>"char" , UInt8=>"unsigned_char",
        UInt16=>"unsigned_short", Int16=>"short",
        UInt32=>"unsigned_int" , Int32=>"int" ,
        UInt64=>"unsigned_long" , Int64=>"long" ,
        Float32=>"float", Float64=> "double")
        typed=Types[DataT]
        print(fid,"SCALARS volume_scalars $typed 1\n");
        print(fid,"LOOKUP_TABLE default\n");
        for j=1:(nx*ny*nz)
            print(fid,"$(V.data[j])\n");
        end
    end
end
export vtkexport

end

# """
#     bar
#
# Compute.
# """
# # Merge the "under" volume into the "over" volume.  The "over" volume voxels have precedence.
# function mergevolume!(Vover::VoxelBoxVolume{CoordT,DataT},
#     Vunder::VoxelBoxVolume{CoordT,DataT},
#     emptyvalue::DataT) where {CoordT<:Number,DataT<:Number}
#     @assert (size(Vover.data) == size(Vunder.data)) "Volumes are incompatible"
#
#     nx, ny, nz = size(Vover.data)
#     voxszx=Vover.dimensions[1]/nx
#     voxszy=Vover.dimensions[2]/ny
#     voxszz=Vover.dimensions[3]/nz
#     for i = 1:nx
#         x = Vover.origin[1] + voxszx*(i+0.5)
#         for j = 1:ny
#             y = Vover.origin[2] + voxszy*(j+0.5)
#             for k = 1:nz
#                 z = Vover.origin[3] + voxszz*(k+0.5)
#                 if (Vunder.data[i,j,k]!=emptyvalue) && (Vover.data[i,j,k]==emptyvalue)
#                     Vover.data[i,j,k] = Vunder.data[i,j,k]
#                 end
#             end
#         end
#     end
#     return Vover
# end


text:
# Spell and einstellgrund encodings in the `vmz_v2_int` file.
#
# Available as a submodule, use as eg `VMZ.sickness`.

module VMZ

import ..@enum_and_Dict

@enum_and_Dict(Einstellgrund,
               einstellgrund_dict,
               "A" => employment_Austria,
               "K" => sickness,
               "M" => nonreporting,
               "P" => retirement,
               "R" => rehab, # NOTE: avoid clash with rehabilitation::Spell
               "S" => other,
               "W" => maternity,
               "X" => precautionary_discontinuation,
               "B" => employment_abroad,
               "V" => forderung_job_search,
               "Z" => forderung_job_offer,
               "F" => forderung_start,
               "-" => unknown,
               "" => unknown, # FIXME should we handle the two types differently?
               "T" => permanent_discontinuation)

@enum_and_Dict(Spell,
               spell_dict,
               "AL" => unemployed,
               "SC" => training,
               "LS" => search_apprentice,
               "SF" => stipend,
               "SR" => rehabilitation,
               "SO" => employment_subsidized,
               "AG" => working_ability_examination,
               "AS" => search_job,
               "AM" => foreign_insurance,
               "AF" => search_early_job,
               "LF" => search_early_training,
               "TA" => partial_foreigner,
               "VM" => preregistered,
               "10" => inability,
               "RE" => rest)

end # module


text:
include("./render.jl")
using .RayTracing
using LinearAlgebra
using Images

const width = 1280
const height = 720
const samples = 100

const logo_img = load("logo.png")

const texture_1 = ConstantTexture(Vec(0.8, 0.8, 0.8))
const texture_2 = ConstantTexture(Vec(0.89, 0.74, 0.25)) #sunshine color
const texture_3 = ImageTexture(logo_img)
const material_1 = Lambertian(texture_1)
const material_2 = Metal(texture_1, 0.1)
const material_3 = Dielectric(1.5)
const material_4 = DiffuseLight(texture_2)
const material_5 = Metal(texture_3, 0.0)

const scene = Scene()
addObject(scene, Object(Sphere(Vec(-0.2, 0.0, -0.2), 0.2), material_1))
addObject(scene, Object(Sphere(Vec(-0.2, 0.0, 0.2), 0.2), material_2))
addObject(scene, Object(Sphere(Vec(0.2, 0.0, 0.2), 0.2), material_3))
addObject(scene, Object(Sphere(Vec(0.2, 0.0, -0.2), 0.2), material_4))
addObject(scene, Object(Rectangle(Vec(1.0, -0.2, 1.0), Vec(-1.0, -0.2, 1.0), Vec(-1.0, -0.2, -1.0)), material_5));

const look_from = Vec(2.0, 3.0, 2.5)
const look_to = Vec(0.0, 0.0, 0.0)
const dist_to_focus = norm(look_from - look_to)
const aperture = 0.1

const camera = Camera(
    look_from,
    look_to,
    Vec(0.0, 1.0, 0.0),
    20.0,
    width / height,
    aperture,
    dist_to_focus)

const image = @time render(scene, camera, (width=width, height=height), samples, Vec(0.5, 0.7, 1.0))

save("output.jpg", image)

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
using JLD
include("src/solver.jl")
k(x, y)=0.6 * sin(2 * pi * x) + 0.2 * sin(4 * pi * y) + 0.1 * sin(10 * pi * x + 8 * pi * y)
N_mesh=512
arr1=[]
arr2=[]
for N_mesh in [256, 512, 1024, 2048, 4096]
    println(N_mesh)
    x=linspace(0,1 - 1 / N_mesh,N_mesh)
    y=linspace(0,1 - 1 / N_mesh,N_mesh)
    sol=zeros(N_mesh,N_mesh)
    for i=1:N_mesh
        for j=1:N_mesh
            sol[i,j]=k(x[i],y[j])
        end
    end
    f=get_f(sol,1.)
    u,err1=solver(f,1.,1e-6)
    push!(arr1,err1)
    u,err2=multigrid_solver(f,1.,1e-6)
    push!(arr2,err2)
end
brr1=[]
brr2=[]

for N_mesh in [256, 512, 1024, 2048, 4096]
    println(N_mesh)
    x=linspace(0,1 - 1 / N_mesh,N_mesh)
    y=linspace(0,1 - 1 / N_mesh,N_mesh)
    sol=zeros(N_mesh,N_mesh)
    for i=1:N_mesh
        for j=1:N_mesh
            sol[i,j]=k(x[i],y[j])
        end
    end
    f=get_f(sol,1.)
    push!(brr1,@elapsed solver(f,1.,1e-6))
    push!(brr2,@elapsed multigrid_solver(f,1.,1e-6))
end

save("plotdata.jld","arr1",arr1,"arr2",arr2,"brr1",brr1,"brr2",brr2)

max_stars_repo_name: JerryLingjieMei/NonlinearPDE

max_issues_repo_name: JerryLingjieMei/NonlinearPDE

max_forks_repo_name: JerryLingjieMei/NonlinearPDE


text:
import java.io.ByteArrayInputStream;
import java.io.ByteArrayOutputStream;
import java.io.IOException;
import java.io.InputStream;
import java.net.URL;
import java.nio.file.Files;
import java.nio.file.OpenOption;
import java.nio.file.Path;
import java.nio.file.StandardOpenOption;
import java.nio.file.attribute.FileAttribute;
import java.nio.file.attribute.PosixFilePermissions;
import java.security.KeyPairGenerator;
import java.security.KeyPair;
import java.security.PublicKey;
import java.security.PrivateKey;

public class URLCache {
    private Path temp;
    private static final FileAttribute[] defaultFileAttributes =
            { PosixFilePermissions.asFileAttribute(PosixFilePermissions.fromString("rwx------")) };
    private static final OpenOption[] defaultOpenOptions =
            { StandardOpenOption.TRUNCATE_EXISTING };
    private final PrivateKey PRIVKEY(H){L} privkey;
    private final PublicKey PUBKEY(H){L} pubkey;

    public URLCache() throws Exception {
	KeyPairGenerator kpg = KeyPairGenerator.getInstance("RSA");
	kpg.initialize(512); // 512 is the keysize
	KeyPair kp = kpg.generateKeyPair();
	privkey = (PRIVKEY(H){L})kp.getPrivate();
	pubkey = (PUBKEY(H){L})kp.getPublic();

        try {
            temp =
                    Files.createTempDirectory("browsercache",
                                              defaultFileAttributes);
        }
        catch (IOException e) {
            e.printStackTrace();
        }
    }

    public InputStream get(URL url) {
        int code = url.toString().hashCode();
        Path file = temp.resolve(Integer.toHexString(code));

        try {
            byte[] bytes = Files.readAllBytes(file);
	    byte[] decrypt = null;
	    try {
		decrypt = CryptoLibrary.decrypt(privkey, bytes);
	    } catch (Exception e) {}
	    return new ByteArrayInputStream(decrypt);
        }
        catch (IOException ioe) {
            return null;
        }
    }

    public void put(URL url, InputStream value) throws IOException {
	byte[]{L} bytes = null;
	try {
	    bytes = CryptoLibrary.encrypt(pubkey, inputStreamToBytes(value));
	} catch (Exception e) {}

	int code = url.toString().hashCode();
	Path file = temp.resolve(Integer.toHexString(code));
	Files.createFile(file, defaultFileAttributes);
	Files.write(file, bytes, defaultOpenOptions);
    }

    private static byte[] inputStreamToBytes(InputStream is) throws IOException {
	ByteArrayOutputStream{L} byteArrayOutputStream =
				     new ByteArrayOutputStream();

	try {
	    int read = is.read();
	    while (read != -1) {
		byteArrayOutputStream.write(read);
		read = is.read();
	    }
	} catch (Exception e) {}

	return byteArrayOutputStream.toByteArray();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using ArgParse
using CSV
using DataFrames
using JLD2
using Sarkses

## =============================================================================
s = ArgParseSettings()
@add_arg_table s begin
    "--sarks"
        help = "jld2 file containing saved trained sarks model"
        required = true
    "--seqs"
        help = "fasta file containing sequences to be scored"
        required = true
    "--k", "-k"
        help = "maximum k-mer length to base smoothed scores on"
        arg_type = Int
        default = 15
end
parsed_args = parse_args(ARGS, s);

sarksFile = parsed_args["sarks"];
seqsFile = parsed_args["seqs"];
k = parsed_args["k"];

## =============================================================================
seqs = Sarkses.readFasta(seqsFile);
@load sarksFile pcSarkses;
mats = [Sarkses.sarksPositionalScoreMat(pcSarks, seqs, k=k)
        for pcSarks in pcSarkses];
CSV.write(stdout,
          DataFrame(hcat(mats...), :auto),
          delim = "\t",
          writeheader = false);

max_stars_repo_name: denniscwylie/deepsarks

max_issues_repo_name: denniscwylie/deepsarks

max_forks_repo_name: denniscwylie/deepsarks


text:
module StructsOfArrays

export StructOfArrays, ScalarRepeat

immutable StructOfArrays{T,N,U<:Tuple} <: AbstractArray{T,N}
    arrays::U
end

type ScalarRepeat{T}
    scalar::T
end
Base.ndims(::ScalarRepeat) = 1
Base.getindex(s::ScalarRepeat, i...) = s.scalar
#should setindex! really be allowed? It will set the index for the whole row...
Base.setindex!{T}(s::ScalarRepeat{T}, value, i...) = (s.scalar = T(value))
Base.eltype{T}(::ScalarRepeat{T}) = T

Base.start(::ScalarRepeat) = 1
Base.next(sr::ScalarRepeat, i) = sr.scalar, i+1
Base.done(sr::ScalarRepeat, i) = false

# since this is used in hot loops, and T.types[.] doesn't play well with compiler
# this needs to be a generated function
@generated function is_tuple_struct{T}(::Type{T})
    is_ts = length(T.types) == 1 && T.types[1] <: Tuple
    :($is_ts)
end
struct_eltypes{T}(struct::T) = struct_eltypes(T)
function struct_eltypes{T}(::Type{T})
    if is_tuple_struct(T) #special case tuple types (E.g. FixedSizeVectors)
        return eltypes = T.types[1].parameters
    else
        return eltypes = T.types
    end
end

make_iterable(x::AbstractArray) = x
make_iterable(x) = ScalarRepeat(x)

@generated function StructOfArrays{T}(::Type{T}, dim1::Integer, rest::Integer...)
    (!isleaftype(T) || T.mutable) && return :(throw(ArgumentError("can only create an StructOfArrays of leaf type immutables")))
    isempty(T.types) && return :(throw(ArgumentError("cannot create an StructOfArrays of an empty or bitstype")))
    dims = (dim1, rest...)
    N = length(dims)
    eltypes  = struct_eltypes(T)
    arrtuple = Tuple{[Array{eltypes[i],N} for i = 1:length(eltypes)]...}

    :(StructOfArrays{T,$N,$arrtuple}(
        ($([:(Array($(eltypes[i]), (dim1, rest...))) for i = 1:length(eltypes)]...),)
    ))
end
StructOfArrays(T::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, dims...)

function StructOfArrays(T::Type, a, rest...)
    arrays = map(make_iterable, (a, rest...))
    N = ndims(arrays[1])
    eltypes = map(eltype, arrays)
    s_eltypes = struct_eltypes(T)
    any(ix->ix[1]!=ix[2], zip(eltypes,s_eltypes)) && throw(ArgumentError(
        "fieldtypes of $T must be equal to eltypes of arrays: $eltypes"
    ))
    any(x->ndims(x)!=N, arrays) && throw(ArgumentError(
        "cannot create an StructOfArrays from arrays with different ndims"
    ))
    arrtuple = Tuple{map(typeof, arrays)...}
    StructOfArrays{T, N, arrtuple}(arrays)
end

Base.linearindexing{T<:StructOfArrays}(::Type{T}) = Base.LinearFast()

@generated function Base.similar{T}(A::StructOfArrays, ::Type{T}, dims::Dims)
    if isbits(T) && length(T.types) > 1
        :(StructOfArrays(T, dims))
    else
        :(Array(T, dims))
    end
end

Base.convert{T,S,N}(::Type{StructOfArrays{T,N}}, A::AbstractArray{S,N}) =
    copy!(StructOfArrays(T, size(A)), A)
Base.convert{T,S,N}(::Type{StructOfArrays{T}}, A::AbstractArray{S,N}) =
    convert(StructOfArrays{T,N}, A)
Base.convert{T,N}(::Type{StructOfArrays}, A::AbstractArray{T,N}) =
    convert(StructOfArrays{T,N}, A)

function Base.size(A::StructOfArrays)
    for elem in A.arrays
        if isa(elem, AbstractArray)
            return size(elem)
        end
    end
    ()
end
Base.size(A::StructOfArrays, d) = size(A)[d]

@generated function Base.getindex{T}(A::StructOfArrays{T}, i::Integer...)
    n = length(struct_eltypes(T))
    Expr(:block, Expr(:meta, :inline),
         :($T($([:(A.arrays[$j][i...]) for j = 1:n]...)))
    )
end

function _getindex{T}(x::T, i)
    is_tuple_struct(T) ? x[i] : getfield(x, i)
end
@generated function Base.setindex!{T}(A::StructOfArrays{T}, x, i::Integer...)
    n = length(struct_eltypes(T))
    quote
        $(Expr(:meta, :inline))
        v = convert(T, x)
        $([:(A.arrays[$j][i...] = _getindex(v, $j)) for j = 1:n]...)
        x
    end
end
end # module